willamettesthd.pdf |
Viability of Willamette River Winter Steelhead: An assessment of the effects of sea lions at Willamette Falls
This website is a companion to an eponymous report . The code and data provided here are intended to facilitate technical review and generate reproducible results.
North Santiam
[num1,txt1,raw1]=xlsread('North Santiam Basin STW summary.xlsx','1985-2017','A2:D143');
T=cell2table(raw1,'VariableNames',{'Creek','Date','Miles','Redd'})
T2=tabulate(txt1(:,1))
X=[];
for i=1:length(T2)
X=[X;repmat(i, [cell2mat(T2(i,2)) 1]),num1(strcmp(txt1(:,1),T2(i)),1:3)];
end
NS=NaN(length(1985:2017),length(T2),2);
for i=1:length(X)
NS(X(i,2)-1984,X(i,1),:)=X(i,3:4);
end
NS=cat(2,reshape(repmat(1985:2017,[1,2])',[],1,2),NS)
%But "Mainstem" is actually 6 different surveys (next),cedar and snake less than 4 years
NS(:,8:10,:)=[]; %delete
%now add back individual mainstem surveys
[num1,txt1,raw1]=xlsread('North Santiam Basin STW summary.xlsx','Mainstem','C15:N21');
num1=flip(num1,2);%Time is backwards
NS(:,8:13,:)=NaN;
for i =1:6
NS(18:27,7+i,2)=num1(i+1,1:10);
end
%Add mainstem miles
NS(18:27,8,1)=num1(2,12);
NS(18:27,9,1)=num1(3,12);
NS(18:27,10,1)=num1(4,12);
NS(18:27,11,1)=num1(5,12);
NS(18:27,12,1)=num1(6,12);
NS(18:27,13,1)=num1(7,12);
NS_redd_mi=array2table([(1985:2017)', NS(:,2:end,2)./NS(:,2:end,1)],...
'VariableNames',{'Year',cell2mat(T2(1)), cell2mat(T2(2)), cell2mat(T2(3)), cell2mat(T2(4)), cell2mat(T2(5)),cell2mat(T2(6)), 'Mint_FB', 'FB_Meh', 'Meh_StayIs','StayIs_Stay','Stay_Green','Green_Mouth' })
%Compare mean of ratio to ratio of mean
MoR= nansum(NS(:,2:end,2)./NS(:,2:end,1),2)./sum(~isnan(NS(:,2:end,2)./NS(:,2:end,1)),2 );
%MoR(all(isnan(NS(:,6:11,2))&isnan(NS(:,6:11,1)),2))=NaN
RoM=nansum(NS(:,2:end,2),2)./nansum(NS(:,2:end,1),2);
%RoM(all(isnan(NS(:,6:11,2))&isnan(NS(:,6:11,1)),2))=NaN
NS_ratios=array2table([(1985:2017)' MoR RoM],'VariableNames',{'Year','MoR','RoM'})
South Santiam
clearvars -except NS NS_redd_mi NS_ratios
[num1,txt1,raw1]=xlsread('SSantiam.xlsx','CrabTree','A3:I35');
[num2,txt2,raw2]=xlsread('SSantiam.xlsx','Thomas','A3:I35');
[num3,txt3,raw3]=xlsread('SSantiam.xlsx','Wiley','A6:I37');%ignore 1966 and 1980
SS(:,1,1)=num1(:,2);
SS(:,1,2)=num1(:,3);
SS(1:(end-1),2,1)=num2(:,2);
SS(1:(end-1),2,2)=num2(:,3);
SS(:,3,1)=num3(:,2);
SS(:,3,2)=num3(:,3);
SS=cat(2,reshape(repmat(1985:2016,[1,2])',[],1,2),SS);
SS_redd_mi=array2table([(1985:2016)',SS(:,2:end,2)./SS(:,2:end,1)],...
'VariableNames',{'Year','CrabTree','Thomas','Wiley' })
%Compare mean of ratio to ratio of mean
MoR= nansum(SS(:,2:end,2)./SS(:,2:end,1),2)./sum(~isnan(SS(:,2:end,2)./SS(:,2:end,1)),2 );
RoM=nansum(SS(:,2:end,2),2)./nansum(SS(:,2:end,1),2);
SS_ratios=array2table([(1985:2016)' MoR RoM],'VariableNames',{'Year','MoR','RoM'})
Calapooia
clearvars -except SS SS_redd_mi SS_ratios NS NS_redd_mi NS_ratios
[num1,txt1,raw1]=xlsread('Calapooia.xlsx','Mainstem','A8:C39');
[num2,txt2,raw2]=xlsread('Calapooia.xlsx','NFk','A4:C35');
[num3,txt3,raw3]=xlsread('Calapooia.xlsx','Potts','A2:C33');%ignore 1966 and 1980
Ca(:,1,1)=num1(:,3);
Ca(:,1,2)=num1(:,2);
Ca(:,2,1)=num2(:,3);
Ca(:,2,2)=num2(:,2);
Ca(:,3,1)=num3(:,3);
Ca(:,3,2)=num3(:,2);
Ca=cat(2,reshape(repmat(1985:2016,[1,2])',[],1,2),Ca);
%The sheet titled 'Calapooia StW spawning Summary.xls' notes several years
%when the count is biased low becuase of survey conditions or timing.
%Remove these values and impute them subsequently...
%remove=[1988; 1990; 1993; 1995; 1996; 1999; 2000; 2005];
%[c,ia,ib]=intersect(Ca(:,1,1),remove(:));
%Ca(ia,2,1)=NaN;
%Ca(ia,2,2)=NaN;
Ca_redd_mi=array2table([(1985:2016)',Ca(:,2:end,2)./Ca(:,2:end,1)],...
'VariableNames',{'Year','CrabTree','Thomas','Wiley' })
Ca_redd_mi=array2table([(1985:2016)',Ca(:,2:end,2)./Ca(:,2:end,1)],...
'VariableNames',{'Year','Mainstem','NFk','Potts' })
%Compare mean of ratio to ratio of mean
MoR= nansum(Ca(:,2:end,2)./Ca(:,2:end,1),2)./sum(~isnan(Ca(:,2:end,2)./Ca(:,2:end,1)),2 );
RoM=nansum(Ca(:,2:end,2),2)./nansum(Ca(:,2:end,1),2);
Ca_ratios=array2table([(1985:2016)' MoR RoM],'VariableNames',{'Year','MoR','RoM'})
Mollala
clearvars -except SS SS_redd_mi SS_ratios NS NS_redd_mi NS_ratios Ca Ca_redd_mi Ca_ratios
num1=xlsread('Molalla.xlsx','Butte 12.6 to 13.2','A2:C33');
num2=xlsread('Molalla.xlsx','Butte 15.1 to 15.6','A2:C33');
num3=xlsread('Molalla.xlsx','Mill','A2:C33');
num4=xlsread('Molalla.xlsx','Camp','A2:C33');
num5=xlsread('Molalla.xlsx','Mid Fk','A2:C33');
num6=xlsread('Molalla.xlsx','Lukens','A2:C33');
num7=xlsread('Molalla.xlsx','Dead Horse','A2:C33');
num8=xlsread('Molalla.xlsx','North Fk','A2:C33');
num9=xlsread('Molalla.xlsx','Cooper','A2:C33');
num10=xlsread('Molalla.xlsx','Mainstem 45.8 to 46.2','A2:C33');
num11=xlsread('Molalla.xlsx','Mainstem 44.9 to 45.8','A2:C33');
num12=xlsread('Molalla.xlsx','Mainstem 42.1 to 43.4','A2:C33');
Mo(:,1,1)=num1(:,3);
Mo(:,1,2)=num1(:,2);
Mo(:,2,1)=num2(:,3);
Mo(:,2,2)=num2(:,2);
Mo(:,3,1)=num3(:,3);
Mo(:,3,2)=num3(:,2);
Mo(:,4,1)=num4(:,3);
Mo(:,4,2)=num4(:,2);
Mo(:,5,1)=num5(:,3);
Mo(:,5,2)=num5(:,2);
Mo(:,6,1)=num6(:,3);
Mo(:,6,2)=num6(:,2);
Mo(:,7,1)=num7(:,3);
Mo(:,7,2)=num7(:,2);
Mo(:,8,1)=num8(:,3);
Mo(:,8,2)=num8(:,2);
Mo(:,9,1)=num9(:,3);
Mo(:,9,2)=num9(:,2);
Mo(:,10,1)=num10(:,3);
Mo(:,10,2)=num10(:,2);
Mo(:,11,1)=num11(:,3);
Mo(:,11,2)=num11(:,2);
Mo(:,12,1)=num12(:,3);
Mo(:,12,2)=num12(:,2);
%plotmatrix(Mo(:,:,2))
Mo=cat(2,reshape(repmat(1985:2016,[1,2])',[],1,2),Mo);
Mo_redd_mi=array2table([(1985:2016)',Mo(:,2:end,2)./Mo(:,2:end,1)],...
'VariableNames',{'Year','Butte_12_6','Butte_15_1',...
'Mill','Camp','Mid_Fk','Lukens','Dead_Horse','North_Fk','Cooper',...
'Mainstem_45_8','Mainstem_44_9','Mainstem_42_1'})
%Compare mean of ratio to ratio of mean
MoR= nansum(Mo(:,2:end,2)./Mo(:,2:end,1),2)./sum(~isnan(Mo(:,2:end,2)./Mo(:,2:end,1)),2 );
RoM=nansum(Mo(:,2:end,2),2)./nansum(Mo(:,2:end,1),2);
Mo_ratios=array2table([(1985:2016)' MoR RoM],'VariableNames',{'Year','MoR','RoM'})
clearvars -except SS SS_redd_mi SS_ratios NS NS_redd_mi NS_ratios Ca Ca_redd_mi Ca_ratios Mo Mo_redd_mi Mo_ratios
Get dam counts
WF=xlsread('Steelhead_WLC - Mapes Edits.xls','AllData&Method','C21:C52');
[Mint1, Mint2, Mint3]=xlsread('Steelhead_WLC - Mapes Edits.xls','AllData&Method','G21:G52');
Mint=cell2mat(Mint3);
[Fost1, Fost2, Fost3]=xlsread('Steelhead_WLC - Mapes Edits.xls','AllData&Method','H21:H52');
Fost=cell2mat(Fost3);
Impute all the missing redd/mi
% <include>Impute3.m</include>
%
X=horzcat(NS(1:(end-1),2:end,2)./NS(1:(end-1),2:end,1), SS(:,2:end,2)./SS(:,2:end,1),...
Ca(:,2:end,2)./Ca(:,2:end,1),Mo(:,2:end,2)./Mo(:,2:end,1),WF,Mint,Fost);
Z=Impute3(X)
%% Get mean of ratio, MoR=mean(count/mi) and relative density
%Can't back calculate redd counts without a length to get RoM
%C=Z.*horzcat(NS(1:(end-1),2:end,1), SS(:,2:end,1),Ca(:,2:end,1),Mo(:,2:end,1));
%how many sites, i, in each population
i_NS=size(NS,2)-1; %subtract 1 for the year column
i_SS=size(SS,2)-1;
i_Ca=size(Ca,2)-1;
i_Mo=size(Mo,2)-1;
NS_MoR= mean(Z(:,1:i_NS),2);
SS_MoR= mean(Z(:,i_NS+1:(i_NS+i_SS)),2);
Ca_MoR= mean(Z(:,(i_NS+i_SS+1):(i_NS+i_SS+i_Ca)),2);
Mo_MoR= mean(Z(:,(i_NS+i_SS+i_Ca+1):(i_NS+i_SS+i_Ca+i_Mo)),2);
MoR=array2table([(1985:2016)',NS_MoR,SS_MoR,Ca_MoR,Mo_MoR],...
'VariableNames',{'Year','N_Santiam','S_Santiam','Calapooia','Molalla' })
RelDensity=array2table([(1985:2016)',MoR{:,2:end}./sum(MoR{:,2:end},2)],...
'VariableNames',{'Year','N_Santiam','S_Santiam','Calapooia','Molalla' })
Apportion Willamette Falls counts to the populations
WF1=WF.*(39+29+7+31)/170; %Jepson et al. 2014 radio telemetry of fish at SS NS Ca Mo
WF2=WF1-Z(:,32)-Fost; %subtract Minto(with imputed years) and Foster
miles=[98 130 55 154]; %spawning miles [NS(below Mint), SS(below Foster), Ca, Mo]
%These mileages computed by Erin Gilbert
N=MoR{:,2:5}.*miles./sum(MoR{:,2:5}.*miles,2).*(WF2+[Z(:,32) Fost zeros(length(WF2),2)]); %add Minto and Foster back
N=array2table([(1985:2016)',N], 'VariableNames',{'Year','N_Santiam','S_Santiam','Calapooia','Molalla' })
figure
plot(1985:2016,N{:,2:5},'.-','MarkerSize',12),xlim([1985,2016]),title('Winter Steelhead')
legend('N. Santiam','S. Santiam','Calapooia','Molalla')
clear WF1 WF2 miles NS_MoR SS_MoR Ca_MoR Mo_MoR
Get age
age1=xlsread('Age Compilation.xlsx','StW Age Data Compilation','B4:D52');
age2=xlsread('Age Compilation.xlsx','StW Age Data Compilation','F4:H14');
age3=xlsread('Age Compilation.xlsx','StW Age Data Compilation','K4:M18');
A2012=xlsread('StW age structure.xlsx','2012','E2:E97');
A2013=xlsread('StW age structure.xlsx','2013','E2:E118');
A2014=xlsread('StW age structure.xlsx','2013','E2:E112');
A2012=[repmat(2012,[size(A2012),1]),A2012, ones(length(A2012),1)];
A2013=[repmat(2013,[size(A2013),1]),A2013, ones(length(A2013),1)];
A2014=[repmat(2014,[size(A2014),1]),A2014, ones(length(A2014),1)];
age=vertcat(age1,age2,age3,A2012,A2013,A2014);
A=zeros(length(unique(age(:,1))),7);
A(:,1)=unique(age(:,1));
for i =1:size(A,1)
Ai=age(age(:,1)==A(i,1),2:3);
for j=1:size(Ai,1)
A(i,1+Ai(j,1))=A(i,1+Ai(j,1))+Ai(j,2);
end
end
%A=array2table(A, 'VariableNames',{'Year','age1','age2','age3','age4','age5','age6' });
clear Ai age age1 age2 age3 A2012 A2013 A2014
Age=NaN(length(1985:2016),7);
Age(:,1)=1985:2016;
Age((A(:,1)-1984),2:7)=A(:,2:7);
Age=array2table(Age, 'VariableNames',{'Year','age1','age2','age3','age4','age5','age6' });
Age
California sea lion mortality
%Basic pinniped mort data from Wright et al.(2016) report
M=[2014 1703 496 712 780; 2015 4149 899 172 557;2016 2252 650 768 915];
M=array2table(M, 'VariableNames',{'Year','H_ChS','W_ChS','StS','StW' })
%Figure 2 (page14) of Wright et al. (2016) shows the spatial and...
%temporal coverage of the mortality estimates.
%In 2016, only falls stratum was estimated. Wright estimates the "river"
%stratum proportion is 0.385. Expand 2016 to include "river"
M{3,5}=M{3,5}/(1-0.385);
%Mean proportion of Pinn mort that is StW
mean_prop=mean(M{:,5}./sum(M{:,2:5},2));
%Remove harvest of StW that aren't in the 4 focal populations
M{:,5}= M{:,5}.*(39+29+7+31)/170 %Jepson et al.
%Start date of the mortality estimate misses a lot of the StW run.
%Get daily counts of StW at Willamette Falls and compare with CSL
%timeseries. Note that Nov1 is start of "winter" run.
figure('Position',[1050 400 900 800])
N14=xlsread('FallsCounts.xlsx','2013_2014','C2:E213');
%yyaxis left
t=datetime(2013, 11, 1)+ caldays(1:212);
subplot(3,1,1),plot(t,N14(:,2),'.-'),ylim([0 max(N14(:,2))]),title('Willamette Falls'),ylabel('StW Daily Count')
hold on
subplot(3,1,1),plot([t(123),t(123)],[0 max(N14(:,2))],'col',[1 0 0])
%yyaxis right
%plot(t,N(:,3)),ylim([0 max(N(:,3))]),title('Willamette Falls'),ylabel('CSL Max Count')
sum(N14(1:123,2))/sum(N14(:,2))% 23% percent of run not included in mortality study
N15=xlsread('FallsCounts.xlsx','2014_2015','C2:E213');
%yyaxis left
t=datetime(2014, 11, 1)+ caldays(1:212);
subplot(3,1,2),plot(t,N15(:,2),'.-'),ylim([0 max(N15(:,2))]),title('Willamette Falls'),ylabel('StW Daily Count')
hold on
subplot(3,1,2),plot([t(101),t(101)],[0 max(N15(:,2))],'col',[1 0 0])
%yyaxis right
%plot(t,N(:,3)),ylim([0 max(N(:,3))]),title('Willamette Falls'),ylabel('CSL Max Count')
sum(N15(1:101,2))/sum(N15(:,2))% 30% of run not included in mortaltiy study
N16=xlsread('FallsCounts.xlsx','2015_2016','C2:E214');
t=datetime(2015, 11, 1)+ caldays(1:213);
subplot(3,1,3)
yyaxis left
plot(t,N16(:,2),'.-'),ylim([0 max(N16(:,2))]),title('Willamette Falls'),ylabel('StW Daily Count')
yyaxis right
subplot(3,1,3),plot(t,N16(:,3),'.-'),ylim([0 max(N16(:,3))]),title('Willamette Falls'),ylabel('CSL Max Count')
hold on
subplot(3,1,3),plot([t(93),t(93)],[0 max(N16(:,3))],'col',[1 0 0])
sum(N16(1:93,2))/sum(N16(:,2))% 22% of run not included in mortaltiy study
%smooth CSL data
CSL=smooth(93:209,N16(93:209,3),0.4,'loess');
figure('Position',[1050 600 700 300])
plot(t(93:209),N16(93:209,3),'.-r')
hold on
plot(t(93:209),CSL,'-g')
plot(t(1:93),repmat(CSL(1),93,1),'-k')
ylabel('CSL max count')
legend('observed','loess smooth','Assumed','Location','northwest')
CSL=[repmat(CSL(1),92,1); CSL; [0 0 0 0]'];
% For 2016, compute "overlap of the StW and CSL
i=sum(CSL(93:end).*N16(93:end,2));%interaction index for mort estimate
i_star=sum(CSL.*N16(:,2));%complete interaction index
EF16=i_star/i
% For 2015, compute "overlap of the StW and CSL
N15=xlsread('FallsCounts.xlsx','2014_2015','D2:D213');
CSL(121)=[];%delete leap year
%mort study begins Feb9, which is StW day 101
i=sum(CSL(101:end).*N15(101:end));%interaction index for mort estimate
i_star=sum(CSL.*N15);%complete interaction index
EF15=i_star/i
% For 2014, compute "overlap of the StW and CSL
N14=xlsread('FallsCounts.xlsx','2013_2014','D2:D213');
%mort study begins March3, which is StW day 123
i=sum(CSL(123:end).*N14(123:end));%interaction index for mort estimate
i_star=sum(CSL.*N14);%complete interaction index
EF14=i_star/i
M{:,5}= M{:,5}.*[EF14;EF15;EF16]% expand for proportion of run not surveyed for mort
%Assume no mortality prior to 1995
M1985_1994=repmat(0,length(1985:1994),1);
%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 StW?
%start+RiverStratum*ProportionStW*FocalPops*RunTimeAdjustment
M1995_2003(1:length(1995:2003))=150/(1-0.385)*mean_prop*(39+29+7+31)/170*mean([EF14,EF15,EF16]);
%linear increas from 2004 to 2013
Nsteps=2013-2004+2;
step=(M{1,5}-M1995_2003(end))/Nsteps;
M2004_2013(1)=M1995_2003(end)+step;
for i=2:length(2004:2013)
M2004_2013(i)= M2004_2013(i-1)+step;
end
Mall=horzcat((1985:2016)',vertcat(M1985_1994,M1995_2003',M2004_2013',M{:,5}));
figure
plot(Mall(:,1),Mall(:,2),'.-')
ylabel('CSL Predation on Focal Winter Steelhead')
RelAbund=N{:,2:5}./sum(N{:,2:5},2);
Mort=array2table([Mall(:,1:2),Mall(:,2).*RelAbund],'VariableNames',{'Year','Total_Mort','N_Santiam','S_Santiam','Calapooia','Molalla' })
HR=array2table([Mort{:,1},Mort{:,3:6}./( Mort{:,3:6}+N{:,2:5})],'VariableNames',{'Year','N_Santiam','S_Santiam','Calapooia','Molalla' })
clear A Ca CSL EF14 EF15 Fost Fost1 Fost2 Fost3 i i_CA i_Mo i_NS i_SS
clear i_star j mean_prop Mint Mint1 dMint2 Mint3 Mo N14 N15 N16 NS SS t WF X Z
clear EF16 i_Ca Mint2
Get proportions of hatchery-origin spawners from previous planning documents
Ph_NS=xlsread('Steelhead_WLC - Mapes Edits.xls','Nsant','B14:B37') ;
Ph_NS=vertcat(Ph_NS,zeros(8,1)); %no hatchery fish 2009-2016;
Ph_SS=xlsread('Steelhead_WLC - Mapes Edits.xls','SSantTotal','N25:N48') ;
Ph_SS=vertcat(Ph_SS,zeros(8,1)); %no hatchery fish 2009-2016;
Ph_Ca=zeros(32,1); %never hatchery fish in Calapooia
Ph_Mo=xlsread('Steelhead_WLC - Mapes Edits.xls','Molalla','D15:D38') ;
Ph_Mo=vertcat(Ph_Mo,zeros(8,1)); %no hatchery fish 2009-2016;
Ph=array2table([(1985:2016)',Ph_NS, Ph_SS, Ph_Ca, Ph_Mo],'VariableNames',{'Year','N_Santiam','S_Santiam','Calapooia','Molalla' } )
clear Ph_NS Ph_SS Ph_Ca Ph_Mo
cd('\\Kalawatseti\home\falcym\docs\Projects\WillamettePinnipedSteelhead')
Calculate recruits
% Use same age comp data for all 4 pops on years when data exist. Use
% mean composition when data do not exist.
Age2=Age{:,2:7}./sum(Age{:,2:7},2);
mean_age=nanmean(Age2);
Age2(isnan(Age2(:,1)),:)=repmat(mean_age,[sum(isnan(Age2(:,1))),1]);
%matrix of Wild spawner abundance at age
for p=1:4
Wvec(:,:,p)= (1-Ph{:,p+1}).*N{:,p+1}.*Age2;
end
%Post-mortality Recruits
for p=1:4
for t=1:(size(Wvec,1)-6)
R1(t,p)=sum(diag(Wvec(t+(1:6),:,p)));
end
end
%Pre-mortality Recruits
%Use Conservation Plan's harvest rates of 0.23 through 1992 and 0.07 after.
Harv=vertcat(repmat(0.23,[length(1985:1992),4]),repmat(0.07,[length(1993:2016),4]));
R2=R1./(1-Harv(7:end,:))+Mort{7:end,3:6};
names={'N.Santiam','S.Santiam','Calapooia','Molalla' };
for p=1:4
subplot(2,2,p),scatter(N{1:end-6,1+p},R2(:,p))
title(names(p)),xlabel('Spawners'),ylabel('Recruits')
end
Preliminaries for JAGS
Data = struct('S',N{1:end-6,2:5},'logR',log(R2));
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,1]);
%I.b= 1000+2000*rand([1,4]);%BH
I.b= 0.01*rand([1,4]);%RK
I.sig_tau = 0+2*rand([1,4]); %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:
% <http://people.oregonstate.edu/~falcym/Model1.txt>
% Model 2:
% <http://people.oregonstate.edu/~falcym/Model2.txt>
% Model 3:
% <http://people.oregonstate.edu/~falcym/Model3.txt>
doparallel =1;
%fprintf( 'Running JAGS...\n' );
tic
[samples, stats, structArray] = matjagsParallel( ...
Data, ... % Observed data
fullfile(pwd, 'Model3.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=1;
plot(samples.sig(1,:,p))
hold on
plot(samples.sig(2,:,p))
plot(samples.sig(3,:,p))
plot(samples.sig(4,:,p))
%figure
%{
for p=1:4
subplot(4,2,p*2-1),plot(samples.a(1,:,p));
hold on
subplot(4,2,p*2-1),plot(samples.a(2,:,p));
subplot(4,2,p*2-1),plot(samples.a(3,:,p));
subplot(4,2,p*2-1),plot(samples.a(4,:,p));
subplot(4,2,p*2),plot(samples.b(1,:,p));
hold on
subplot(4,2,p*2),plot(samples.b(2,:,p));
subplot(4,2,p*2),plot(samples.b(3,:,p));
subplot(4,2,p*2),plot(samples.b(4,:,p));
end
%}
figure
subplot(5,1,1),plot(samples.a(1,:));
hold on
subplot(5,1,1),plot(samples.a(2,:));
subplot(5,1,1),plot(samples.a(3,:));
subplot(5,1,1),plot(samples.a(4,:));
for p=1:4
subplot(5,1,1+p),plot(samples.b(1,:,p));
hold on
subplot(5,1,1+p),plot(samples.b(2,:,p));
subplot(5,1,1+p),plot(samples.b(3,:,p));
subplot(5,1,1+p),plot(samples.b(4,:,p));
end
figure
for pop=1:4;
clear Ri Si
Si=1:max(N{1:end-6,1+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(2,2,pop),plot(Si,Ri(:,i),'col',[0.1,0.1,0.1,0.1])
hold on
end
subplot(2,2,pop),scatter(N{1:end-6,1+pop},R2(:,pop),'or')
xlim([0 max(max(N{1:end-6,1+pop}),max(R2(:,pop)))])
ylim([0 max(max(N{1:end-6,1+pop}),max(R2(:,pop)))])
line([0 max(max(N{1:end-6,1+pop}),max(R2(:,pop)))],[0,max(max(N{1:end-6,1+pop},max(R2(:,pop))))],'Color',[0 0 1])
Ri=((stats.mean.a.*Si).*exp(-stats.mean.b(pop).*Si))';
subplot(2,2,pop),plot(Si,Ri,'g')
xlabel('Spawners')
ylabel('Wild Recruits')
title(names(pop))
end
% plot BH
%{
pop=1;
Si=1:max(N{1:end-6,1+pop});
Ri=(stats.mean.a.*Si)./(1+stats.mean.a/stats.mean.b(1).*Si);
plot(Si,Ri)
hold on
scatter(N{1:end-6,1+pop},R2(:,pop))
Si=1:max(N{1:end-6,1+pop});
for rep=1:100
rani=randi(40000);
for i=1:length(Si)
Ri(i,rep)=(Si(i)*ai(rani))/(1+ai(rani)/bi(rani)*Si(i));
end
end
figure
hold on
for i=1:100
plot(Si,Ri(:,i),'col',[0.1,0.1,0.1,0.1])
end
scatter(N{1:end-6,1+pop},R(:,pop),'or')
xlim([0 max(max(S(1:31),max(R)))])
ylim([0 max(max(S(1:31),max(R)))])
line([0 max(max(S(1:31),max(R)))],[0,max(max(S(1:31),max(R)))],'Color',[0 0 1])
xlabel('Spawners')
ylabel('Wild Recruits')
title(names(pop))
%}
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 Molalla
% respectively. 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
% not attempt to incorporate ocean harvest by age.
pop=1;
scenario=2;
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: No pinniped mortality. Add 0.07 for incidental fishery
terminal_h=zeros([26,1])+0.07;
case 2
%Scenario 2: Lowest empiricical pinniped mortality (2015), 13% Add 0.07 for incidental
terminal_h =repmat(HR{(end-1),1+pop}+0.07,[26,1]);
case 3
%Scenario 3: 2016 mortality repeats in the future. Add 0.07 for incidental
terminal_h=repmat(HR{end,1+pop}+0.07,[26,1]);
case 4
%Scenario 4: 2017 mortality (33%) Add 0.07 for incidental fishery
terminal_h= repmat(0.33+0.07,[26,1]);
end
Coshak=ones(26,6);%No expansion for ocean harvest. Holdover from ChF PVA
QET=100;
RFT=QET;
TermR=R2(:,pop).*Age2(1:26,:);
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));
b=samples.b(j(1),j(2),pop);
for i=1:length(R2)
z5(i)=log(R2(i))-(log(a)+log(N{i,1+pop})-b*N{i,1+p});%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=mean(N{:,1+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
mins=toc/60
mean(out1(:,3))
Here is the a single simulation with the actual timeseries superimposed. This is 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
subplot(2,2,1),plot(X(51:150,8),'.-b')
%ylim([0 4000])
hold on
plot(33:64,N{:,pop+1},'.-r')
%plot(41:66,xlsread('ChF_Revision',population,'B4:B29'),'.-r')%Elk River.
xlabel('Year','Fontsize',14)
ylabel('Spawner Abundance','Fontsize',14)
title(names(pop),'Fontsize',14)
legend('Simulated','Empirical')