Deschutes Redband, 2016
This analysis is an update to the 2015 Deschutes River Redband report. Data from 2016 are included here. Future analyses could incorporate previous years of data in order to assess change through time.
The Data
cd('\\Kalawatseti\home\falcym\docs\Projects\People\Seals')
[num,txt]=xlsread('Redband2016','Redband2016','A1:E411');
num2=num(num(:,4)<10,:);%remove 98 year old fish
section=txt([false;num(:,4)<10],1);%remove header and remove 98 year old fish
T=table(section, num2(:,2), num2(:,3), num2(:,4),'VariableNames',{'Location' 'FL_mm' 'Wt_g' 'Age'})
gscatter(T.Age,T.FL_mm,section,'brgc','o'), ylim([0 500]), xlabel('Age'), ylabel('Fork Length (mm)')
A Single von Bertalanffy Model
This is the model I will fit:
where is asymptotic average maximum body size, K is a growth rate coefficient that determines how quickly the max is attained, and is the hypothetical age at which the species has zero length. ϵ is normally distributed error. are data (T.FL_mm) andt are data (T.Age)
Maximize the Likelihood
f=@(params) VB(params,T.Age,T.FL_mm)
[x, fval,exitflag,output,grad,hessian]=fminunc(f,[400,0.5,0.5,100]);%I want that Hessian
exitflag
output
fval
%[x, fval,exitflag,output]=fminsearch(f,[400,0.5,0.5,100])% same result?
%%
AICc=2*fval + 2*length(x)*(length(T.FL_mm)/(length(T.FL_mm)-length(x)-1))
Tabular Results
varcov=inv(hessian);
stders=sqrt(diag(varcov));
c=stders*sqrt(chi2inv(0.95,1));
out=[x', stders, x'-c, x'+c];
table(x', stders, x'-c, x'+c,...
'RowNames',{'L_inf';'K';'t_0';'SD'},'VariableNames',{'Est' 'SE' 'L95CI' 'U95CI'})
Graphical Result
gscatter(T.Age,T.FL_mm,section,'brgc','o'), ylim([0 500]), xlabel('Age'), ylabel('Fork Length (mm)')
hold on
xi=1:0.01:max(T.Age);
yi=x(1).*(1-exp(-x(2).*(xi-x(3))));
plot(xi,yi,'k')
Multiple von Bertalanffy Models
Now I fit von Bertalanffy models to fish from each stream section. These models will share the error term, $\epsilon$, but all other paramters are estimated separately. Thus we have,
%
% for the p = 1,2,3,4 different stream sections.
%
% <include>VB2.m</include>
Maximize Likelihood
f2=@(params) VB2(params,T.Age,T.FL_mm,num2(:,1))
options=optimoptions(@fminunc,'MaxFunctionEvaluations',5000,'FunctionTolerance',1e-12); %This is 13 dimensional problem now.
[x, fval,exitflag,output,grad,hessian]=fminunc(f2,[reshape(repmat([400,0.5,0.5],[1,4]),1,[]),100],options);%I want that Hessian
exitflag
output
fval
AICc=2*fval + 2*length(x)*(length(T.FL_mm)/(length(T.FL_mm)-length(x)-1))
Tabular Results
varcov=inv(hessian);
stders=sqrt(diag(varcov));
c=stders*sqrt(chi2inv(0.95,1));
out=[x', stders, x'-c, x'+c];
names={'Nena','Jones','N.Junction','Trout Cr'};
%for i=1:4;
%names(i)
%table(x(i*3-2:i*3)', stders(i*3-2:i*3), x(i*3-2:i*3)'-c(i*3-2:i*3), x(i*3-2:i*3)'+c(i*3-2:i*3),...
% 'RowNames',{'L_inf';'K';'t_0'},'VariableNames',{'Est' 'SE' 'L95CI' 'U95CI'})
%end
i=1
names(i)
table(x(i*3-2:i*3)', stders(i*3-2:i*3), x(i*3-2:i*3)'-c(i*3-2:i*3), x(i*3-2:i*3)'+c(i*3-2:i*3),...
'RowNames',{'L_inf';'K';'t_0'},'VariableNames',{'Est' 'SE' 'L95CI' 'U95CI'})
i=2
names(i)
table(x(i*3-2:i*3)', stders(i*3-2:i*3), x(i*3-2:i*3)'-c(i*3-2:i*3), x(i*3-2:i*3)'+c(i*3-2:i*3),...
'RowNames',{'L_inf';'K';'t_0'},'VariableNames',{'Est' 'SE' 'L95CI' 'U95CI'})
i=3
names(i)
table(x(i*3-2:i*3)', stders(i*3-2:i*3), x(i*3-2:i*3)'-c(i*3-2:i*3), x(i*3-2:i*3)'+c(i*3-2:i*3),...
'RowNames',{'L_inf';'K';'t_0'},'VariableNames',{'Est' 'SE' 'L95CI' 'U95CI'})
i=4
names(i)
table(x(i*3-2:i*3)', stders(i*3-2:i*3), x(i*3-2:i*3)'-c(i*3-2:i*3), x(i*3-2:i*3)'+c(i*3-2:i*3),...
'RowNames',{'L_inf';'K';'t_0'},'VariableNames',{'Est' 'SE' 'L95CI' 'U95CI'})
'Shared SD'
table(x(13), stders(13), x(13)'-c(13), x(13)'+c(13),...
'RowNames',{'SD'},'VariableNames',{'Est' 'SE' 'L95CI' 'U95CI'})
Graphical Results
figure
gscatter(T.Age,T.FL_mm,section,'brgc','o'), ylim([0 500]), xlabel('Age'), ylabel('Fork Length (mm)')
hold on
xi=1:0.01:max(T.Age);
yi1=x(1).*(1-exp(-x(2).*(xi-x(3))));
yi2=x(4).*(1-exp(-x(5).*(xi-x(6))));
yi3=x(7).*(1-exp(-x(8).*(xi-x(9))));
yi4=x(10).*(1-exp(-x(11).*(xi-x(12))));
plot(xi,yi1,'b')
plot(xi,yi2,'r')
plot(xi,yi3,'g')
plot(xi,yi4,'c')
Comentary on von Bertalanffy growth differences between river sections
The relative AICc values (3992.7 < 4052.2) indicate that the model with separate growth curves for each river section is "better" than the model that ignores river section. This assessment does not try to determine which componets of growth are most likely to vary between river sections (i.e. ). It is possible that some growth parameters can be shared across river sections. This could be explored further.
A Single Length-weight Model
A power function can describe the relationship between length (L) and weight (W)of the ith fish: where α and β are fitted parameters, and error, ϵ is lognormally distributed. Taking logs we obtain: and ϵ is now normally distributed.
%
% <include>LW.m</include>
%
Maximize Likelihood
f=@(params) LW(params,T.FL_mm,T.Wt_g)
options=optimoptions(@fmincon,'MaxFunctionEvaluations',5000,'MaxIteration',5000);
%Going to add consraints to the optimization so that hessian does't get
%whacky
A=[];
b=[];
Aeq=[];
beq=[];
lb=[0 0 0];
ub=[];
[x, fval,exitflag,output,lambda,grad,hessian]=fmincon(f,[0.5,2,0.2],A,b,Aeq,beq,lb,ub);%I want that Hessian
exitflag
output
fval
%%
% AICc
AICc=2*fval + 2*length(x)*(length(T.FL_mm)/(length(T.FL_mm)-length(x)-1))
Tabular Results
varcov=inv(hessian);
stders=sqrt(diag(varcov));
c=stders*sqrt(chi2inv(0.95,1));
out=[x', stders, x'-c, x'+c];
table(x', stders, x'-c, x'+c,...
'RowNames',{'alpha';'beta';'SD'},'VariableNames',{'Est' 'SE' 'L95CI' 'U95CI'})
Graphical Results
figure
gscatter(T.FL_mm,T.Wt_g,section,'brgc','o')
xlim([0 500])
ylim([0 1400])
hold on
Li=0.01:500;
Wi=x(1).*Li.^x(2);
plot(Li,Wi,'-k')
xlabel('Fork Length (mm)')
ylabel('Weight')
Multiple Length-weight Models
I now fit a unique length-weight relationship to each river segment. As with the multiple von Bertalanffy model, I will force each river segment to have teh same error variance. where α and β are fitted parameters, and error, ϵ is lognormally distributed. Taking logs we obtain: and ϵ is again normally distributed and shared among all sections (no 'p' subscript).
%
% <include>LW2.m</include>
%
Maximize Likelihood
f2=@(params) LW2(params, T.FL_mm, T.Wt_g, num2(:,1))
%Going to add consraints to the optimization so that hessian does't get
%whacky
A=[];
b=[];
Aeq=[];
beq=[];
lb=[1e-10 0.1 1e-10 0.1 1e-10 0.1 1e-10 0.1 0.05];
ub=[];
options=optimoptions(@fmincon,'MaxFunctionEvaluations',5000,'MaxIteration',...
5000,'StepTolerance',1e-20,'ConstraintTolerance',1e-20);
[x, fval,exitflag,output,lambda,grad,hessian]=fmincon(f2,...
[reshape(repmat([7e-05,2.7],[1,4]),1,[]),0.25],A,b,Aeq,beq,lb,ub,[],options);
%repmat([7e-05,2.7,0.25],[1,4]),A,b,Aeq,beq,lb,ub,[],options)
%{
options=optimset('TolFun',1e-12,'MaxFunEvals',50000);
[x, fval,exitflag,output]=fminsearch(f2,repmat([7e-05,2.7,0.25],[1,4]),options)% same result
%[x, fval,exitflag,output]=fminsearch(f2,[reshape(repmat([7e-05,2.6],[1,4]),1,[]),0.25],options)% same result
options=optimoptions('fminunc','TolFun',1e-20,'MaxFunEvals',10000, 'FunctionTolerance',1e-50);
[x, fval,exitflag,output]=fminunc(f2,[reshape(repmat([9e-05,2.9],[1,4]),1,[]),0.28],options)% same result
%}
exitflag
output
fval
%%
AICc=2*fval + 2*length(x)*(length(T.FL_mm)/(length(T.FL_mm)-length(x)-1))
Tabular Results
varcov=inv(hessian);
stders=sqrt(diag(varcov));
c=stders*sqrt(chi2inv(0.95,1));
out=[x', stders, x'-c, x'+c];
names={'Nena','Jones','N.Junction','Trout Cr'};
%for i=1:4;
%names(i)
%table(x(i*2-1:i*2)', stders(i*2-1:i*2), x(i*2-1:i*2)'-c(i*2-1:i*2), x(i*2-1:i*2)'+c(i*2-1:i*2),...
% 'RowNames',{'alpha';'beta'},'VariableNames',{'Est' 'SE' 'L95CI' 'U95CI'})
%end
i=1;
names(i)
table(x(i*2-1:i*2)', stders(i*2-1:i*2), x(i*2-1:i*2)'-c(i*2-1:i*2), x(i*2-1:i*2)'+c(i*2-1:i*2),...
'RowNames',{'alpha';'beta'},'VariableNames',{'Est' 'SE' 'L95CI' 'U95CI'})
i=2;
names(i)
table(x(i*2-1:i*2)', stders(i*2-1:i*2), x(i*2-1:i*2)'-c(i*2-1:i*2), x(i*2-1:i*2)'+c(i*2-1:i*2),...
'RowNames',{'alpha';'beta'},'VariableNames',{'Est' 'SE' 'L95CI' 'U95CI'})
i=3;
names(i)
table(x(i*2-1:i*2)', stders(i*2-1:i*2), x(i*2-1:i*2)'-c(i*2-1:i*2), x(i*2-1:i*2)'+c(i*2-1:i*2),...
'RowNames',{'alpha';'beta'},'VariableNames',{'Est' 'SE' 'L95CI' 'U95CI'})
i=4;
names(i)
table(x(i*2-1:i*2)', stders(i*2-1:i*2), x(i*2-1:i*2)'-c(i*2-1:i*2), x(i*2-1:i*2)'+c(i*2-1:i*2),...
'RowNames',{'alpha';'beta'},'VariableNames',{'Est' 'SE' 'L95CI' 'U95CI'})
'Shared SD'
table(x(9), stders(9), x(9)'-c(9), x(9)'+c(9),...
'RowNames',{'SD'},'VariableNames',{'Est' 'SE' 'L95CI' 'U95CI'})
Graphical Results
figure
gscatter(T.FL_mm,T.Wt_g,section,'brgc','o')
xlim([0 500])
ylim([0 1400])
hold on
Li=0.01:500;
Wi1=x(1).*Li.^x(2);
plot(Li,Wi1,'-b')
Wi2=x(3).*Li.^x(4);
plot(Li,Wi2,'-r')
Wi3=x(5).*Li.^x(6);
plot(Li,Wi3,'-g')
Wi4=x(7).*Li.^x(8);
plot(Li,Wi4,'-c')
xlabel('Fork Length (mm)')
ylabel('Weight (g)')
Comentary on length-weight differences between river sections
The relative AICc values (-3.3036 < 37.443) provides strong evidence for differences between river sections. As with the multiple von Bertalanffy assessment, this analysis does not try to determine which parameters (α or β) differ between sections. A better model could involve sharing parameters between sets of sections or relaxing the assumption that that magnitude of errors is identical for all sections.
Catch-at-age Analysis
Catch at age analysis assumes no emigration or immigration. Fish grow until they are "recruited into the sampling gear." Numbers at age then decline exponentially with a rate of decline equal to the instantaneous total mortality rate. The model is:
where is the number of fish caught in year y, is the initial recruitment into the cohort, and Z is the instantaneous rate of total mortality.The mortality parameter Z is estimated with linear regerssion after a log transformation. It is customary to exclude fish that are not yet recruited into the sampling gear from the regression.
x=tabulate(T.Age);
x1=tabulate(T.Age(num2(:,1)==1));
x2=tabulate(T.Age(num2(:,1)==2));
x3=tabulate(T.Age(num2(:,1)==3));
x4=tabulate(T.Age(num2(:,1)==4));
figure
plot(1:5,x1(:,3),'--.b','Markersize',24)
hold on
plot(1:6,x2(:,3),'--.r','Markersize',24)
plot(1:6,x3(:,3),'--.g','Markersize',24)
plot(1:5,x4(:,3),'--.c','Markersize',24)
plot(1:6,x(:,3),'--.k','Markersize',24)
xlabel('Age')
ylabel('Relative Frequency')
legend('Nena','Jones','N.Junction','Trout Cr','Combined')
figure
plot(1:5,log(x1(:,2)),'--.b','Markersize',24)
hold on
plot(1:6,log(x2(:,3)),'--.r','Markersize',24)
plot(1:6,log(x3(:,3)),'--.g','Markersize',24)
plot(1:5,log(x4(:,3)),'--.c','Markersize',24)
plot(1:6,log(x(:,3)),'--.k','Markersize',24)
xlabel('Age')
ylabel('Log(catch)')
legend('Nena','Jones','N.Junction','Trout Cr','Combined','Location','SW')
The multiple von Bertalanffy model revealed that Trout Cr fish grow the fastest. Here, we see that Trout Cr fish are recruited into the fishing gear earlier than the others. The absense of a linear decline in log(catch) after the time of recruitment (age 3) indicates that the assumed exponential decline with age may not be appropriate. It may work for Jones, however.
[b, bint, r, rint,stats]=regress(log(x2(3:6,3)),[ones(4,1),(3:6)']);
Pvalue=stats(3)
Rsquared=stats(1)
slope=b(2) %This is instantaneous mortality rate, Z
The annual mortality rate, A, is
A=1-exp(b(2))
%%
% Note that the annual mortality rate is only computed for Jones, and that a
% regression was fit to only 4 data points.
Condition Factor
The condition factor is computed for "small" fish (fork length <= 300 mm) and "large" fish (fork length > 300 mm). Within each size category, c, compute the mean .
for p=1:4
K_hat(p,1)= nanmean((T.Wt_g(T.FL_mm<=300 & num2(:,1)==p)./T.FL_mm(T.FL_mm<=300 & num2(:,1)==p).^3).*100000);
K_hat(p,2)= nanmean((T.Wt_g(T.FL_mm>300 & num2(:,1)==p)./T.FL_mm(T.FL_mm>300 & num2(:,1)==p).^3).*100000);
end
K_hat
Now I'll do some emprical (nonparametric) bootsrapping to generate confidence intervals.
The difference between the empirical sample mean condition factor and the true mean μ is . I use an emprical bootstrap to to approximate the distribution . This is done by randomly sampling individual W and L data, with replacement, and then computing
detla=NaN([9999, 2,4]); %preallocate for speed
for rep=1:9999
for p=1:4
delta(rep,1,p)=nanmean(datasample((T.Wt_g(T.FL_mm<=300 & num2(:,1)==p)./T.FL_mm(T.FL_mm<=300 & num2(:,1)==p).^3).*100000,length(T.FL_mm(T.FL_mm<=300 & num2(:,1)==p))))-K_hat(p,1);
delta(rep,2,p)=nanmean(datasample((T.Wt_g(T.FL_mm>300 & num2(:,1)==p)./T.FL_mm(T.FL_mm>300 & num2(:,1)==p).^3).*100000,length(T.FL_mm(T.FL_mm<=300 & num2(:,1)==p))))-K_hat(p,2);
end
end
figure('Position',[1300,200,500,600])
subplot(4,2,1),hist(delta(:,1,1),50),title('Nena, Small'),xlabel('\delta^*'),ylabel('Frequency'),ylim([0 800]),xlim([-1 1])
subplot(4,2,2),hist(delta(:,2,1),50),title('Nena, Large'),xlabel('\delta^*'),ylabel('Frequency'),ylim([0 800]),xlim([-0.5 0.5])
subplot(4,2,3),hist(delta(:,1,2),50),title('Jones, Small'),xlabel('\delta^*'),ylabel('Frequency'),ylim([0 800]),xlim([-0.5 0.5])
subplot(4,2,4),hist(delta(:,2,2),50),title('Jones, Large'),xlabel('\delta^*'),ylabel('Frequency'),ylim([0 800]),xlim([-0.5 0.5])
subplot(4,2,5),hist(delta(:,1,3),50),title('N.Junction, Small'),xlabel('\delta^*'),ylabel('Frequency'),ylim([0 800]),xlim([-0.5 0.5])
subplot(4,2,6),hist(delta(:,2,3),50),title('N.Junction, Large'),xlabel('\delta^*'),ylabel('Frequency'),ylim([0 800]),xlim([-0.5 0.5])
subplot(4,2,7),hist(delta(:,1,4),50),title('Trout Cr, Small'),xlabel('\delta^*'),ylabel('Frequency'),ylim([0 800]),xlim([-0.5 0.5])
subplot(4,2,8),hist(delta(:,2,4),50),title('Trout Cr, Large'),xlabel('\delta^*'),ylabel('Frequency'),ylim([0 800]),xlim([-0.5 0.5])
hold on
Then find values of at the 5th and 95th percetiles.
count=0;
for p=1:4
for c=1:2
count=count+1;
pt(p,:,c)=prctile(delta(:,c,p),[5 95]);
subplot(4,2,count),line([pt(p,1,c),pt(p,1,c)],[0 800])
subplot(4,2,count),line([pt(p,2,c),pt(p,2,c)],[0 800])
end
end
The boostrap confidence interval forμ is .
%get sample size
for p=1:4
N(p,1)=length(T.FL_mm(T.FL_mm<=300 & num2(:,1)==p));
N(p,2)=length(T.FL_mm(T.FL_mm>300 & num2(:,1)==p));
end
Small=table(K_hat(:,1),K_hat(:,1)+pt(:,1,1),K_hat(:,1)+pt(:,2,1),N(:,1),...
'RowNames',{'Nena';'Jones';'N.Junction';'Trout Cr'},'VariableNames',{'Est' 'L95CI' 'U95CI' 'N'})
Large=table(K_hat(:,2),K_hat(:,2)+pt(:,1,2),K_hat(:,2)+pt(:,2,2),N(:,2),...
'RowNames',{'Nena';'Jones';'N.Junction';'Trout Cr'},'VariableNames',{'Est' 'L95CI' 'U95CI' 'N'})
figure('position',[1300,300, 800,300])
subplot(1,2,1),plot(1:4,K_hat(:,1),'.k','MarkerSize',20),xlim([0.5 4.5]),title('Small')
xticklabels({'Nena','Jones','N. Junct','Trout Cr.'})
ylabel('Condition Factor, K, and 95% CI')
hold on
for p=1:4
line([p,p],[table2array(Small(p,2)) table2array(Small(p,3))])
end
subplot(1,2,2),plot(1:4,K_hat(:,2),'.k','MarkerSize',20),xlim([0.5 4.5]),title('Large')
xticklabels({'Nena','Jones','N. Junct','Trout Cr.'})
ylabel('Condition Factor, K, and 95% CI')
hold on
for p=1:4
line([p,p],[table2array(Large(p,2)) table2array(Large(p,3))])
end
Commentary on condition factors
The data plotted in the figures in the section on the length-weight relastionship reveal two "outlier" small fish in Nena. These outliers caused the strange resampling distribution for small Nena fish in this section. To a lesser extent, the exact same thing occurs with N. Junction.
Since we have no reason to discard the "outliers" as a data entry errors, the results presented here suggest that small fish in Jones and Trout Cr have signifigantly lower condition factors than small fish in Nena and N. Junction. For large fish, Jones has signifigantly lower condition factor than Nena and Trout cr. These visual assessments of the confidence intervals do not attempt to control for multiple tests (e.g. Bonferonni), and Tukey HSD or similar method should be used test for differences between specific pairs of sections.