spring_chinook_status_assessment.pdf |
Spring Chinook Status Assessment: McKenzie, Clackamas, and Sandy River Populations
Pinnipeds Mortaltiy Rates
cd('N:\docs\Projects\CHS_McKenzieClackSandy\McKenzie')
clear M
M=[2014 1703 496 712 780; 2015 4149 899 172 557;2016 2252 650 768 915; 2017 1824 399 181 270];
M=array2table(M, 'VariableNames',{'Year','H_ChS','W_ChS','StS','StW' });
Total morts in 2016 and 2017 are for Falls stratum only. Expand for River stratum by factor described in Wright et al. 2017 Table 5:
M{3,2:5}=M{3,2:5}/(1-2870/7455);
M{4,2:5}=M{4,2:5}/(1-1615/4288)
N=[2014 6412; 2015 8948; 2016 6631; 2017 5905];
N=array2table(N, 'VariableNames',{'Year','W_ChS'})
Mortality Rate 2014 to 2017
MR1417=M{:,3}./(M{:,3} + N{:,2})
Plot Functional Relationship
scatter((M{:,3} + N{:,2}),MR1417,'k','filled')
xlabel('Spring Chinook Abundance at WF')
ylabel('Mortality Rate')
title('Predation Functional Relationship')
ylim([0 0.22])
xlim([6000, 10000])
for i=1:4
text((M{i,3} + N{i,2})+30,MR1417(i),num2str(M{i,1}))
end
Assume linear decrease in harvest rate from 2014 back to 2002. Wright et al. (2014): "predation losses were generally a few hundred or less" during late 1990s-2003 How does 150 salmonid morts scale to these WCHS?
% Mean proportion of all mort that is WCHS
mean_prop=mean(M{:,3}./sum(M{:,2:5},2));
MR_pre2004=mean(mean_prop*150./(mean_prop*150+M{:,3}));
% Make linear increase between 2003 and 2014;
MRstep=(MR1417(1)-MR_pre2004)/(2014-2004);
MR=zeros(2016-2002+1,2);
MR(:,1)=(2002:2016)';
MR(1:2,2)=MR_pre2004;
for i=3:(length(MR)-4)
MR(i,2)=MR(i-1,2)+MRstep;
end
MR((end-3):end,2)=MR1417;
figure
plot(MR(1:11,1),MR(1:11,2),'o-k','markerfacecolor','white','markersize',6)
hold on
plot(MR(12:end,1),MR(12:end,2),'.-k','markersize',14)
xlabel('Year')
ylabel('Mortality Rate from CSL at WF')
title('McKenzie CHS')
line([MR(11,1),MR(12,1)],[MR(11,2),MR(12,2)],'col','k')
legend('Inferred','Observed')
Load McKenzie NOS, HOS, and age data
cd('N:\docs\Projects\CHS_McKenzieClackSandy\McKenzie')
Age36=xlsread('Streamnet_CoordinatedAssessmentsFINAL.xlsx','McKenzie','D11:G25')
NOS=xlsread('Streamnet_CoordinatedAssessmentsFINAL.xlsx','McKenzie','J50:J64')
H1=xlsread('Streamnet_CoordinatedAssessmentsFINAL.xlsx','McKenzie','E50:E64')
H2=xlsread('Streamnet_CoordinatedAssessmentsFINAL.xlsx','McKenzie','F50:F64')
[num,txt,raw] =xlsread('Streamnet_CoordinatedAssessmentsFINAL.xlsx','McKenzie','G50:G64');
H3=[raw{:}]';
pHOS=nansum([H1,H2,H3],2)./(NOS+nansum([H1,H2,H3],2))
%plot(2002:2016,NOS,'.-k')
%xlabel('Year')
%ylabel('Natural Origin Spawners')
%title('McKenzie')
Do Trend Analysis on NOS
cd('N:\docs\Projects\CHS_McKenzieClackSandy')
NOSpred=vertcat(NOS,repmat(NaN,50,1));
clear init initStructs
fixedyear=floor(length(NOS)/2);
for chain=1:4
init(chain) = struct('alpha', normrnd(log(NOS(1)),1,[1, 1]),'beta', normrnd(0,0.1,[1, 1]),...
'tauepsilon', normrnd(1,0.1,[1, 1]));
end
initStructs=init;
Data = struct('N',NOSpred, 'Nyears', length(NOSpred),'fixedyear',fixedyear);
doparallel = 1; % do not use parallelization
fprintf( 'Running JAGS...\n' );
tic
[samples, stats, structArray] = matjagsParallel( ...
Data, ... % Observed data
fullfile(pwd, 'Trend.txt'), ... % File that contains model definition
init, ... % Initial values for latent variables
'doparallel' , doparallel, ... % Parallelization flag
'nchains', 4,... % Number of MCMC chains
'nburnin', 300000,... % Number of burnin steps
'nsamples', 10000, ... % Number of samples to extract
'thin', 41, ... % Thinning parameter
'dic', 1, ... % Do the DIC?
'monitorparams', {'alpha','beta','sd_epsilon','B','fitted'}, ... % List of latent variables to monitor
'savejagsoutput' , 1 , ... % Save command line output produced by JAGS?
'verbosity' , 0 , ... % 0=do not produce any output; 1=minimal text output; 2=maximum text output
'cleanup' , 0 ); % clean up of temporary files?
toc
delete(gcp)
Prob_decline=sum(sum(samples.B(:,:)<0))/(4*10000)
stats.Rhat
%{
subplot(2,1,1),plot(samples.alpha(1,:))
hold on
subplot(2,1,1),plot(samples.alpha(2,:))
subplot(2,1,1),plot(samples.alpha(3,:))
subplot(2,1,1),plot(samples.alpha(4,:))
ylabel('\alpha')
subplot(2,1,2),plot(samples.beta(1,:))
hold on
subplot(2,1,2),plot(samples.beta(2,:))
subplot(2,1,2),plot(samples.beta(3,:))
subplot(2,1,2),plot(samples.beta(4,:))
ylabel('\beta')
%}
figure
plot(samples.alpha(1,:))
hold on
plot(samples.alpha(2,:))
plot(samples.alpha(3,:))
plot(samples.alpha(4,:))
ylabel('\alpha')
figure
plot(samples.beta(1,:))
hold on
plot(samples.beta(2,:))
plot(samples.beta(3,:))
plot(samples.beta(4,:))
ylabel('\beta')
figure
[xi,f]=ksdensity(reshape(samples.B,[],1));
%subplot(3,1,1),
figure, plot(f,xi)
CIofB=hpdi(reshape(samples.B,[],1),95)
xlabel('Geometric Mean Rate of Interannual Change, %')
ylabel({'Probability Density (blue)','95% interval (red)'})
title('McKenzie')
%subplot(3,1,1),
line([CIofB(1),CIofB(1)],[0 0.25],'col',[1 0 0])
%subplot(3,1,1),
line([CIofB(2),CIofB(2)],[0 0.25],'col',[1 0 0])
horizon=15;
for t=1:(15+horizon)
CI(t,1:2)=hpdi(reshape(samples.fitted(:,:,t),[],1),95);
end
figure
xlim([2002 2015+16])
%subplot(3,1,2),
fill([2002:(2016+horizon) (horizon+2016):-1:2002],[CI(:,1); flipud(CI(:,2))],'b')
alpha 0.2
hold on
%subplot(3,1,2),
plot(2002:2016,NOS,'.-k','markersize',14)
xlabel('Year')
ylabel('Natural Origin Spawners')
title('McKenzie ChS')
Compute Recruits
Apply the Mortality Rates to McKenzie Look at three different assumptions of Relative Reproductive Success (RRS) of hatchery fish.
% Harvest 2002-2017. Table 13 of "Fisheries Managmenet and Evalutiaon for 2017 Willamette River Spring Chinook."
L.ColCom=[0.024, 0.011, 0.042, 0.020, 0.080, 0.027, 0.005, 0.015, 0.060, 0.047, 0.025, 0.046, 0.030, 0.034, 0.016, 0.014];
L.ColRec=[0.011, 0.012, 0.010, 0.007, 0.008, 0.008, 0.002, 0.008, 0.018, 0.006, 0.010, 0.005, 0.009, 0.007, 0.003, 0.001];
L.Will= [0.030, 0.024, 0.027, 0.030, 0.043, 0.042, 0.048, 0.036, 0.069, 0.062, 0.045, 0.028, 0.039, 0.040, 0.022, 0.024];
L.Clack= [0.049, 0.008, 0.003, 0.005, 0.003, 0.001, 0.002, 0.002, 0.002, 0.002, 0.003, 0.003, 0.001, 0.001, 0.001, 0.003];
U.Will= [0.003, 0.003, 0.000, 0.001, 0.001, 0.001, 0.001, 0.001, 0.002, 0.001, 0.002, 0.001, 0.001, 0.001, 0.001, 0.001];
McKen= [0.000, 0.054, 0.019, 0.007, 0.014, 0.013, 0.001, 0.013, 0.010, 0.011, 0.013, 0.008, 0.005, 0.004, 0.005, 0.002];
HR=1-(1-L.ColCom).*(1-L.ColRec).*(1-L.Will).*(1-U.Will).*(1-McKen)
MR
RRS=[1 0.5 0];
for i=1:1
RRSi=repmat(RRS(i),[length(H1),1]);
S_McK=nansum([NOS,(RRSi.*[H1,H2,H3])],2);
NA36=NOS.*Age36;
NA36=NA36./( (1-(MR(:,2))').*(1-L.ColCom(1:15)).*(1-L.ColRec(1:15)).*(1-L.Will(1:15)).*(1-U.Will(1:15)).*(1-McKen(1:15)))';%McKenzie
for t=1:9
%R(t)=sum(diag(NA36((t+3:t+6),:)))./(1-MR(i+6,2));
R_McK(t)=sum(diag(NA36((t+3:t+6),:)));
end
%{
subplot(1,3,i),scatter(S(1:9),R)
xlabel({'Spawners', 'natural and RRS-adjusted hatchery fish'})
ylabel('Natural-origin Recruits')
title(['RRS= ' num2str(RRS(i)) '' ])
xlim([0 6500])
ylim([0 6500])
hold on
line([0 6500],[0 6500],'col',[0 0 0])
%}
end
%R/S bump
%Use i=2, RRS=0.5 for this bump exercise.
median(R_McK'./S_McK(1:9))
std(R_McK'./S_McK(1:9))
Load SR data from Clack and Sandy
cd('N:\docs\Projects\CHS_McKenzieClackSandy\Clackamas')
%S_Clack=xlsread('Clackamas run size 2002-2017.xlsx','Recruits are Preharvest','V3:V12') %RRS=0.5
S_Clack=xlsread('Clackamas run size 2002-2017.xlsx','Recruits are Preharvest','T3:T12') %RRS=1
R_Clack=xlsread('Clackamas run size 2002-2017.xlsx','Recruits are Preharvest','S3:S12')
cd('N:\docs\Projects\CHS_McKenzieClackSandy\Sandy')
%S_Sandy=xlsread('SandyCHS.xlsx','SandyChS','Q6:Q15') %RRS=0.5
S_Sandy=xlsread('SandyCHS.xlsx','SandyChS','O6:O15') %RRS=1
R_Sandy=xlsread('SandyCHS.xlsx','SandyChS','N6:N15')
%Format McK for consistency with Sandy and Clack
S_McK(length(R_McK)+1:end)=[];
R_McK=R_McK';
Preliminaries for JAGS
S=vertcat(S_McK,S_Clack,S_Sandy);
R=vertcat(R_McK,R_Clack,R_Sandy);
p=vertcat(repmat(1,[9,1]),repmat(2,[10,1]),repmat(3,[10,1]));
Data = struct('S',S,'logR',log(R),'p',p,'Tmax',length(S));
nchains = 4; % How Many Chains?
nburnin = 35000; % How Many Burn-in Samples?
nsamples = 10000; % How Many Recorded Samples?
% To run a different recrutitment model, change the last number in the
% lines in the loop below to match the needed number of paramters.
clear init I
for i=1:nchains
I.a= 5+2*rand([1,3]);
%I.b= 1000+2000*rand([1,4]);%BH
I.b= 0.01*rand([1,3]);%RK
I.sig_tau = 0+2*rand([1,1]); %Log scale
init(i) = I; % init is a structure array that has the initial values for all latent variables for each chain
end
Calling JAGS to fit recruitment model
%
% Here are the JAGS files for the three Ricker recruitment models
%
% Model 1:
% <https://falcy.weebly.com/uploads/1/1/4/9/114983925/model1.txt>
% Model 2:
% <https://falcy.weebly.com/uploads/1/1/4/9/114983925/model2.txt>
cd('N:\docs\Projects\CHS_McKenzieClackSandy')
doparallel =1;
%fprintf( 'Running JAGS...\n' );
tic
[samples, stats, structArray] = matjagsParallel( ...
Data, ... % Observed data
fullfile(pwd, 'Model2.txt'), ... % File that contains model definition
init, ... % Initial values for latent variables
'doparallel' , doparallel, ... % Parallelization flag
'nchains', nchains,... % Number of MCMC chains
'nburnin', nburnin,... % Number of burnin steps
'nsamples', nsamples, ... % Number of samples to extract
'thin', 13, ... % Thinning parameter
'dic', 1, ... % Do the DIC?
'monitorparams', {'a','b','sig','density'}, ... % List of latent variables to monitor
'savejagsoutput' , 1 , ... % Save command line output produced by JAGS?
'verbosity' , 0 , ... % 0=do not produce any output; 1=minimal text output; 2=maximum text output
'cleanup' , 0 ); % clean up of temporary files?
delete(gcp)
Compute WAIC
%lppd = log pointwise predictive density
%lppd=sum(sum(log(1/(nsamples*nchains).*sum(sum(samples.density,2),4)),1),3);
%The above lppd is idential to the one below
reshape1=permute(samples.density,[2,4,3,1]);
reshape2=reshape(reshape1,size(reshape1,1)*size(reshape1,2),[]);%rows are mcmc, cols are data
lppd=sum(log((1/(nsamples*nchains))*sum(reshape2)));
pwaic=sum(var(log(reshape2),0,1));
WAIC=-2*(lppd-pwaic)
Inspect model fit
stats.Rhat.a
stats.Rhat.b
figure
p=3;
plot(samples.a(1,:,p))
hold on
plot(samples.a(2,:,p))
plot(samples.a(3,:,p))
plot(samples.a(4,:,p))
Plot the fitted model
figure('Position',[2100 300 560*0.8 420*2])
names={'McKenzie','Clackamas','Sandy'};
mx=[max(S_McK),max(S_Clack),max(S_Sandy)];
Sx=horzcat([S_McK;NaN],S_Clack,S_Sandy);
Rx=horzcat([R_McK;NaN],R_Clack,R_Sandy);
for pop=1:3;
clear Ri Si
Si=1:max(Sx(:,pop));
%Ri=((stats.mean.a.*Si).*exp(-stats.mean.b(pop).*Si))';
%subplot(2,2,pop),plot(Si,Ri)
%hold on
%subplot(2,2,pop),scatter(N{1:end-6,1+pop},R2(:,pop))
bi=(reshape(samples.b(:,:,pop),[],1));
ai=(reshape(samples.a(:,:),[],1));
for rep=1:100
rani=randi(40000);
for i=1:length(Si)
Ri(i,rep)=(Si(i).*ai(rani)).*exp(-bi(rani).*Si(i));
end
end
for i=1:100
subplot(3,1,pop),plot(Si,Ri(:,i),'col',[0.1,0.1,0.1,0.1])
hold on
end
subplot(3,1,pop),scatter(Sx(:,pop),Rx(:,pop),'or')
xlim([0 max(max(Sx(:,pop)),max(Rx(:,pop)))])
ylim([0 max(max(Sx(:,pop)),max(Rx(:,pop)))])
line([0 max(max(Sx(:,pop)),max(Rx(:,pop)))],[0,max(max(Sx(:,1),max(Rx(:,pop))))],'Color',[0 0 1])
Ri=((stats.mean.a(pop).*Si).*exp(-stats.mean.b(pop).*Si))';
subplot(3,1,pop),plot(Si,Ri,'g')
xlabel('Spawners, RRS=0.5')
ylabel('Wild Recruits')
title(names(pop))
end
Do PVA with MCMC samples
This code is set up to run a PVA on the N.Santiam. Change "pop" immediately below to 2, or 3, or 4 for S.santiam, Calapooia, and Molallarespectively. Scenarios 1 through 4 simulate different effects of sea lions. Note that this PVA is capable of incorporating ocean harvest. There is a lot of useless dividing by 1 in this PVA because there it does not attempt to incorporate ocean harvest by age.
clear X
pop=1;
scenario=5;
mean_age=mean([zeros(length(find(1-isnan(Rx(:,1)))),1),Age36(1:length(find(1-isnan(Rx(:,1)))),:),zeros(length(find(1-isnan(Rx(:,1)))),1)]);
tmax=160;% This allows a burn-in period in the PVA
reps=100;%per rand param
randparams=1000;% rand param
out1=[];
switch scenario
case 1
%Scenario 1: McKenzie: random sample of history
terminal_h =1-((1-MR(:,2))'.*(1-L.ColCom(1:15)).*(1-L.ColRec(1:15)).*(1-L.Will(1:15)).*(1-U.Will(1:15)).*(1-McKen(1:15)));
case 2% McKenzie: Max CSL, recent harvest
terminal_h =1-((1-max(MR(:,2)))'.*(1-L.ColCom(end)).*(1-L.ColRec(end)).*(1-L.Will(end)).*(1-U.Will(end)).*(1-McKen(end)));
case 3 %Clackamas
terminal_h =1-((1-L.ColCom(end)).*(1-L.ColRec(end)).*(1-L.Will(end)).*(1-L.Clack(end)));
case 4%Sandy
terminal_h =1-((1-L.ColCom(end)).*(1-L.ColRec(end)).*(1-L.Clack(end))); %Use Clack sport fishery as proxy
case 5% McKenzie: No CSL, recent harvest
terminal_h =1-((1-L.ColCom(end)).*(1-L.ColRec(end)).*(1-L.Will(end)).*(1-U.Will(end)).*(1-McKen(end)));
case 6% McKenzie: Historical CSL, recent harvest
terminal_h =1-((1-L.ColCom(1:15)).*(1-L.ColRec(end)).*(1-L.Will(end)).*(1-U.Will(end)).*(1-McKen(end)));
end
Coshak=ones(9,6);%No expansion for ocean harvest. Holdover from ChF PVA
QET=100;
RFT=QET;
TermR=Rx(find(1-isnan(Rx(:,1))),1).*[zeros(length(find(1-isnan(Rx(:,1)))),1),Age36(1:length(find(1-isnan(Rx(:,1)))),:),zeros(length(find(1-isnan(Rx(:,1)))),1)];
CK=Coshak;
TotR=Coshak.*TermR;
for i=1:length(Coshak)
TermR_AComp(i,:)=TermR(i,1:6)./sum(TermR(i,1:6),2);%Terminal Recruit Age Comp
TotR_AComp(i,:)=TotR(i,1:6)./sum(TotR(i,1:6),2);%Terminal Recruit Age Comp
end
tic
%Sample param estimates
for randparam = 1:randparams
j=[randi([1,4]),randi(length(samples.a))];
a=samples.a(j(1),j(2),pop);
b=samples.b(j(1),j(2),pop);
for i=1:length(Rx(find(1-isnan(Rx(:,1)))))
z5(i)=log(Rx(i,pop))-(log(a)+log(Sx(i,pop))-b*Sx(i,pop));%Ricker
end
rho=corr(z5(1:(length(z5)-1))',z5(2:length(z5))');
%QET1=0;
%QET2=0;
%QET3=0;
QET4=0;
for rep=1:reps
%qet1=0;
%qet2=0;
%qet3=0;
qet4=0;
X=zeros(tmax,29);
%BEGIN simulate normally distributed deviates and Coshak Expansion
sig=std(z5);
X(1,9)=0;
for t=2:tmax
X(t,9)=rho*X(t-1,9)+sig*sqrt(1-rho^2)*randn;
i=randperm(length(Coshak),1);
X(t,11:16)=CK(i,:);
X(t,17:22)=TotR_AComp(i,:);
X(t,29)=terminal_h(randperm(length(terminal_h),1));%terminal harvest
end
%END simulate normally distributed deviates and Coshak Expansion
N0=nanmean(Sx(:,pop));
Ni=N0.*mean_age;
X(1,1)=Ni(1);
X(1:2,2)=Ni(2);
X(1:3,3)=Ni(3);
X(1:4,4)=Ni(4);
X(1:5,5)=Ni(5);
X(1:6,6)=Ni(6);
for t =2:(tmax-7)
St=sum(X(t,1:7));
X(t,8)=St;
lnRt=log(a)+log(St)-b*St+X(t,9);%RK
%Rt=log(a)+log(St)-log(1+a/b*St)+X(t,8);%BH
if St>RFT
X(t,10)=exp(lnRt);
else
X(t,10)=0; %Implement RFT
end
X(t,23:28)=X(t,10).*X(t,17:22)./X(t,11:16) ;%Terminal Recruits by age
for i=2:6
X(t+i,i)=nansum([X(t+i,1),X(t,22+i)*(1-X(t+i,29))]);% Make Spawners
end
%{
if isempty(Ph)==1
X(t+1:t+7,1:7)=X(t+1:t+7,1:7)+diag(floor(X(t,11).*A));%No hatchery fish
else
k=randi(length(Ph));
H=Ph(k)*X(t,11)/(1-Ph(k));%Ph=H/(H+W). Rearranging gives H=(Ph*W)/(1-Ph)
X(t+1:t+7,1:7)=X(t+1:t+7,1:7)+diag(floor((H+X(t,11)).*A));%hatchery fish
end
%}
end%time t
u=X(51:150,8)<QET;
%{
%if min(X(21:120,9))<QET
% ext=ext+1;
%end
%for i = 1:length(u)
% qet1=qet1+u(i);
%end
%for i = 1:length(u)-1
% qet2=qet2+u(i)*u(i+1);
%end
%}
%for i = 1:length(u)-2
% qet3=qet3+u(i)*u(i+1)*u(i+2);
%end
for i = 1:length(u)-3
qet4=qet4+u(i)*u(i+1)*u(i+2)*u(i+3);
end
%QET1=QET1+(qet1>0);
%QET2=QET2+(qet2>0);
%QET3=QET3+(qet3>0);
QET4=QET4+(qet4>0);
end%rep
out1=[out1; a b QET4/reps];
end% randperm- samples of params
delete(gcp)
mins=toc/60
mean(out1(:,3))
Here is a single simulation with the actual timeseries superimposed. This intends to demonstrate that the PVA centers abundance appropriately, captures stochasticisity and autocorrlation. In short, the simulated timeseries should "look like" the obseved time series.
figure
plot(1:100,X(51:150,8),'.-b')
%ylim([0 4000])
hold on
plot(1:9,S_McK,'.-r')
xlabel('Year','Fontsize',14)
ylabel('Spawner Abundance','Fontsize',14)
%title(names(pop),'Fontsize',14)
legend('Simulated','Empirical')