Middle Fork John Day Spring Chinook Life Cycle Model
Matt Falcy
September 6, 2019
1. The Data
These data received from Ian Tattam, ODFW.
setwd('//Kalawatseti/home/falcym/docs/Projects/People/Tattam')
Z<-read.csv('ExportFishData.csv',header=T,stringsAsFactors = FALSE)
Z
## BroodYear Redds Spawners_Redd TotalSpawners Age3 Age4 Age5
## 1 2002 389 2.89 1124 0.025 0.960 0.015
## 2 2003 236 2.92 689 0.056 0.810 0.135
## 3 2004 319 2.27 724 0.009 0.950 0.037
## 4 2005 178 2.10 374 0.013 0.960 0.026
## 5 2006 199 2.42 482 0.032 0.960 0.006
## 6 2007 85 2.95 251 0.031 0.860 0.108
## 7 2008 169 2.16 365 0.058 0.870 0.073
## 8 2009 251 3.24 813 0.051 0.880 0.067
## 9 2010 197 2.71 534 0.065 0.850 0.083
## 10 2011 505 3.89 1964 0.036 0.920 0.040
## 11 2012 493 3.01 1484 0.000 0.873 0.127
## 12 2013 113 4.68 529 0.149 0.634 0.217
## 13 2014 518 2.87 1487 0.026 0.896 0.077
## 14 2015 442 2.35 1039 0.009 0.786 0.205
## 15 2016 160 2.91 466 0.040 0.640 0.320
## 16 2017 29 2.76 80 0.125 0.750 0.125
## FallMigrants SpringMigrants TotalMigrants FallSurvival SpringSurvival
## 1 723 19425 20149 0.28 0.64
## 2 818 19355 20173 0.28 0.70
## 3 NA 18465 18465 0.28 0.68
## 4 127 13454 13581 0.28 0.68
## 5 NA 7382 7382 0.28 0.74
## 6 8187 30332 38519 0.28 0.70
## 7 1203 34262 35465 0.28 0.69
## 8 2746 12457 15203 0.28 0.90
## 9 18655 29214 47869 0.28 0.80
## 10 16870 32271 49142 0.28 0.68
## 11 556 18347 18903 0.28 0.54
## 12 58 6249 6307 0.28 0.27
## 13 10368 15990 26358 0.28 0.62
## 14 5615 26795 32410 0.28 0.69
## 15 3930 19436 23366 0.28 0.57
## 16 NA NA NA NA NA
## SpringMigrants_JDDam FallMigrants_JDDam Smolts_JDDam SAR
## 1 12508 203 12711 0.033
## 2 13529 229 13758 0.015
## 3 12588 NA 12588 0.021
## 4 9107 36 9143 0.051
## 5 5492 NA 5492 0.062
## 6 21315 2295 23609 0.075
## 7 23679 337 24016 0.051
## 8 11214 770 11983 0.008
## 9 23325 5229 28554 0.036
## 10 21790 4729 26518 0.059
## 11 9950 156 10105 0.045
## 12 1700 16 1716 0.051
## 13 9850 2906 12756 0.026
## 14 18381 1574 19955 NA
## 15 11079 1102 12180 NA
## 16 NA NA NA NA
2. Does the proportion migrating in Fall depend on the number of redds?
It would be interesting to find density-dependence here.
pMF<-Z$FallMigrants/(Z$FallMigrants+Z$SpringMigrants)
plot(Z$Redds,pMF,pch=21,cex=1.2,bg='red',xlab='Redds',ylab='Proportion Migrating in Fall')
cor(Z$Redds,pMF,use="complete.obs")
## [1] 0.314736
summary(lm(Z$Redds~pMF))
##
## Call:
## lm(formula = Z$Redds ~ pMF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -222.12 -123.97 -45.11 148.43 249.05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 233.80 65.43 3.573 0.00437 **
## pMF 344.97 313.68 1.100 0.29493
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 157.5 on 11 degrees of freedom
## (3 observations deleted due to missingness)
## Multiple R-squared: 0.09906, Adjusted R-squared: 0.01715
## F-statistic: 1.209 on 1 and 11 DF, p-value: 0.2949
Answer: No. THe plot shows nothing. Correlation coefficient (0.31) isnโt too bad, but the regression as a P-value of 0.3 and an R^2 of 0.017.
3. Is there evidence of density-dependence from TotalSpawners to TotalMigrants?
plot(Z$TotalSpawners,Z$TotalMigrants,pch=21,cex=1.2,bg='red',xlab='Total Spawners',ylab='Total Migrants',xlim=c(0, 2100),ylim=c(0, 50000))
text(Z$TotalMigrants~Z$TotalSpawners,labels=Z$BroodYear,cex=0.9,font=1,pos=4)
Answer: No relationship. Fit a BH model anyway just for curiousity.
4. Maximum likelihood estimates of multiple Beverton-Holt models
4.1 Model BH0: Beverton-Holt, no covariates
BH0<-function(data,params){
-sum(dnorm(log(data[,1]),log(params[1])+log(data[,2])-log(1+params[2]*data[,2]),params[3],log=TRUE))
}
BH0fit<-optim(par=c(50,1,3),fn=BH0,data=cbind(Z$TotalMigrants[1:15],Z$TotalSpawners[1:15]),control=list(reltol=1e-16,maxit=2000))
Si<-seq(0,max(Z$TotalMigrants[1:14]))
Ri<-(BH0fit$par[1]*Si)/(1+BH0fit$par[2]*Si)
plot(Z$TotalSpawners,Z$TotalMigrants,pch=21,cex=1.2,bg='red',xlab='Total Spawners',ylab='Total Migrants',xlim=c(0, 2100),ylim=c(0, 50000))
text(Z$TotalMigrants~Z$TotalSpawners,labels=Z$BroodYear,cex=0.9,font=1,pos=4)
lines(Si,Ri)
This is really questionable. The optimizer generated NaNs, and the point estimate of intrinsit productivity isโฆ
BH0fit$par[1]
## [1] 1243.589
โฆmigrants per spawner.
4.2 Model BH1: Beverton-Holt with July mean flow.
covs<-read.csv('ExportCovs.csv',header=T)
covs
## Brood.Year MeanQ JulyT95 JulyMaxT PDO
## 1 2002 54.0 19 26.8 2.96
## 2 2003 80.0 9 na 3.48
## 3 2004 51.4 9 na 0.28
## 4 2005 57.0 17 28.2 0.91
## 5 2006 31.5 17 27.9 -7.63
## 6 2007 107.7 5 13.7 -1.11
## 7 2008 68.0 0 24.9 -3.53
## 8 2009 96.1 3 15.8 -6.45
## 9 2010 262.4 3 21.7 -7.79
## 10 2011 60.7 4 25.1 -3.47
## 11 2012 39.8 13 27.1 5.07
## 12 2013 58.9 8 25.5 8.08
## 13 2014 31.2 4 27.2 6.60
## 14 2015 36.9 2 26.1 2.18
## 15 2016 50.5 7 25.2 NA
zmeanQ<-(covs$MeanQ[1:14]-mean(covs$MeanQ[1:15]))/sd(covs$MeanQ[1:15])
BH1<-function(data,params){
-sum(dnorm(log(data[,1]),log(params[1])+log(data[,2])-log(1+params[2]*data[,2])+params[3]*data[,3],params[4],log=TRUE))
}
BH1fit<-optim(par=c(100,0.01,0,0.5),fn=BH1,data=cbind(Z$TotalMigrants[1:15],Z$TotalSpawners[1:15],zmeanQ),control=list(reltol=1e-16,maxit=2000),hessian=TRUE)
Si<-seq(0,max(Z$TotalMigrants[1:15]))
Ri<-(BH1fit$par[1]*Si)/(1+BH1fit$par[2]*Si)
plot(Z$TotalSpawners,Z$TotalMigrants,pch=21,cex=1.2,bg='red',xlab='Total Spawners',ylab='Total Migrants',xlim=c(0, 2100),ylim=c(0, 55000))
text(Z$TotalMigrants~Z$TotalSpawners,labels=Z$BroodYear,cex=0.9,font=1,pos=4)
lines(Si,Ri)
Sx<-Z$TotalSpawners[1:15]
Ri<-(BH1fit$par[1]*Sx)/(1+BH1fit$par[2]*Sx)
Rx<-(BH1fit$par[1]*Sx)/(1+BH1fit$par[2]*Sx)*exp(BH1fit$par[3]*zmeanQ)
for (t in 1:15){
lines(c(Sx[t],Sx[t]),c(Ri[t],Rx[t]))
}
This seems to be an improvement. Slope at the origin is much lower. Look at the output from optim()
BH1fit
## $par
## [1] 1.420521e+02 5.043384e-03 2.918441e-01 5.227222e-01
##
## $value
## [1] 11.5535
##
## $counts
## function gradient
## 1155 NA
##
## $convergence
## [1] 0
##
## $message
## NULL
##
## $hessian
## [,1] [,2] [,3] [,4]
## [1,] 2.720531e-03 -58.99243 1.581856e-03 8.881784e-09
## [2,] -5.899243e+01 1396083.83428 2.501670e+02 3.859125e+02
## [3,] 1.581856e-03 250.16696 5.107837e+01 2.620570e-06
## [4,] 8.881784e-09 385.91246 2.620570e-06 1.097978e+02
Indeed, the point estimate of intrinsic productivity is now much lower (115.35). The numbers in โ$parโ (short for parameters) above are just point estimates. What about standard errors or confidence intervals of the paramters? In the chunk of code above, I asked the optimizer to return the Hessian. Iโll use this to get the SEs and CIs. A Hessian is a square matrix of second-order partial derivates evaluted at the maximum likelihood estimates of each paramter:
\[H = \begin{pmatrix} \frac{\partial^2f}{\partial x_1^2} & \frac{\partial^2f}{\partial x_1 \partial x_2} & \cdots & \frac{\partial^2f}{\partial x_1 \partial x_n} \\ \frac{\partial^2f}{\partial x_2^2} & \frac{\partial^2f}{\partial x_2^2} & \cdots & \frac{\partial^2f}{\partial x_2 \partial x_n} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^2f}{\partial x_n \partial x_1} & \frac{\partial^2f}{\partial x_n \partial x_2} & \cdots & \frac{\partial^2f}{\partial x_n^2} \end{pmatrix}\]
This matrix gives the magnitue of the curvature in all directions of the multidimenesional likelihood surface. A more โspikedโ surface means that a parameter can only take on a narrow range of values, whereas a flatter surface means that the parameter estimates are more uncertain.
We can get standard errors by calculating the inverse of the square root of the diagonal elements of the observed Fisher Information matrix \(\mathbf I(\mathbf \theta)\), where \(\mathbf \theta\) is a vector of all the paramters. The Fisher Information matrix is found by inverting the Hessian matrix.
\(\mathbf I(\mathbf \theta) = \mathbf H^{-1}\)
It is important to note that the Hessian pertains to a negative log likelihood, becuase thatโs what we minimized in with the optimizer. If we had maximized the log likelihood without minimizing the negative log likelihood, then weโd have to take the negtaive of the Hessian before inverting it.
And now:
\(s_\mathbf \theta = \frac{1}{\sqrt{I(\theta)}}\).
From here the 95% confidence intervals are easy. Hereโs how to implement all this in Rโฆ
fisher.info<-solve(BH1fit$hessian)
sig<-sqrt(diag(fisher.info))
upper<-BH1fit$par+1.96*sig
lower<-BH1fit$par-1.96*sig
interval<-data.frame(BH1fit$par,lower=lower,upper=upper)
interval
## BH1fit.par lower upper
## 1 1.420521e+02 10.6370436776 273.46724501
## 2 5.043384e-03 -0.0007631154 0.01084988
## 3 2.918441e-01 0.0156973562 0.56799092
## 4 5.227222e-01 0.3345615223 0.71088297
The effect of mean July flow (parameter #3) is statistically significant because its 95% CIs do not overlap zero.
4.3 Model BH2: Beverton-Holt with July high air temps above 95F
zT95<-(covs$JulyT95[1:15]-mean(covs$JulyT95[1:15]))/sd(covs$JulyT95[1:15])
BH2<-function(data,params){
-sum(dnorm(log(data[,1]),log(params[1])+log(data[,2])-log(1+params[2]*data[,2])+params[3]*data[,3],params[4],log=TRUE))
}
BH2fit<-optim(par=c(10,0.01,0,0.5),fn=BH2,data=cbind(Z$TotalMigrants[1:15],Z$TotalSpawners[1:15],zT95),control=list(reltol=1e-16,maxit=2000),hessian=TRUE)
Si<-seq(0,max(Z$TotalMigrants[1:14]))
Ri<-(BH2fit$par[1]*Si)/(1+BH2fit$par[2]*Si)
plot(Z$TotalSpawners,Z$TotalMigrants,pch=21,cex=1.2,bg='red',xlab='Total Spawners',ylab='Total Migrants',xlim=c(0, 2100),ylim=c(0, 55000))
text(Z$TotalMigrants~Z$TotalSpawners,labels=Z$BroodYear,cex=0.9,font=1,pos=4)
lines(Si,Ri)
Sx<-Z$TotalSpawners[1:15]
Ri<-(BH2fit$par[1]*Sx)/(1+BH2fit$par[2]*Sx)
Rx<-(BH2fit$par[1]*Sx)/(1+BH2fit$par[2]*Sx)*exp(BH2fit$par[3]*zT95)
for (t in 1:15){
lines(c(Sx[t],Sx[t]),c(Ri[t],Rx[t]))
}
Good. And now the 95% confidence intervals:
fisher.info<-solve(BH2fit$hessian)
sig<-sqrt(diag(fisher.info))
upper<-BH2fit$par+1.96*sig
lower<-BH2fit$par-1.96*sig
interval<-data.frame(BH2fit$par,lower=lower,upper=upper)
interval
## BH2fit.par lower upper
## 1 564.66000501 -1.315616e+03 2444.9362270
## 2 0.02478292 -6.295471e-02 0.1125205
## 3 -0.35943219 -6.069804e-01 -0.1118840
## 4 0.47250771 3.033510e-01 0.6416644
Notice that intrinsic prodcutivity is much higher than in the previous model, and that temperature has a significant effect.
4.4 Model BH3: Beverton-Holt that includes both flow and temperature
BH3<-function(data,params){
-sum(dnorm(log(data[,1]),log(params[1])+log(data[,2])-log(1+params[2]*data[,2])+params[3]*data[,3]+params[4]*data[,4],params[5],log=TRUE))
}
BH3fit<-optim(par=c(100,0.01,0,0,0.5),fn=BH3,data=cbind(Z$TotalMigrants[1:15],Z$TotalSpawners[1:15],zmeanQ,zT95),control=list(reltol=1e-16,maxit=2000),hessian=TRUE)
Si<-seq(0,max(Z$TotalMigrants[1:15]))
Ri<-(BH3fit$par[1]*Si)/(1+BH3fit$par[2]*Si)
plot(Z$TotalSpawners,Z$TotalMigrants,pch=21,cex=1.2,bg='red',xlab='Total Spawners',ylab='Total Migrants',xlim=c(0, 2100),ylim=c(0, 55000))
text(Z$TotalMigrants~Z$TotalSpawners,labels=Z$BroodYear,cex=0.9,font=1,pos=4)
lines(Si,Ri)
Sx<-Z$TotalSpawners[1:15]
Ri<-(BH3fit$par[1]*Sx)/(1+BH3fit$par[2]*Sx)
Rx<-(BH3fit$par[1]*Sx)/(1+BH3fit$par[2]*Sx)*exp(BH3fit$par[3]*zmeanQ)*exp(BH3fit$par[4]*zT95)
for (t in 1:14){
lines(c(Sx[t],Sx[t]),c(Ri[t],Rx[t]))
}
Good. And now the 95% confidence intervals:
fisher.info<-solve(BH3fit$hessian)
sig<-sqrt(diag(fisher.info))
upper<-BH3fit$par+1.96*sig
lower<-BH3fit$par-1.96*sig
interval<-data.frame(BH3fit$par,lower=lower,upper=upper)
interval
## BH3fit.par lower upper
## 1 108.522828072 35.668181452 1.813775e+02
## 2 0.003463763 0.000431818 6.495709e-03
## 3 0.279340153 -0.019476051 5.781564e-01
## 4 -0.088205363 -0.438768903 2.623582e-01
## 5 0.514425598 0.264602561 7.642486e-01
When both flow and temperature are used as covariates, flow becomes insignificant.
4.5 AICc
From here we need AICc to ajudicate between these four models. \(AICc = AIC+\frac{2k^{2}+2k}{n-k-1}=-2ln\left(L\right)+2k+\frac{2k^{2}+2k}{n-k-1}\). We already have the negative log liklihood \(-ln\left(L\right)\) because that is what the optimizer maximized (it actually minimized the negative). This is recorded in fit$value.
AICcBH0<- 2*BH0fit$value +2*3 + (2*3^2+2*3)/(14-3-1)
AICcBH1<- 2*BH1fit$value +2*4 + (2*4^2+2*4)/(14-4-1)
AICcBH2<- 2*BH2fit$value +2*4 + (2*4^2+2*4)/(14-4-1)
AICcBH3<- 2*BH3fit$value +2*5 + (2*5^2+2*5)/(14-5-1)
data.frame(AICcBH0,AICcBH1,AICcBH2,AICcBH3)
## AICcBH0 AICcBH1 AICcBH2 AICcBH3
## 1 34.94471 35.55145 32.52156 38.58562
The Beverton-Holt model with temperature is the best model.
5. Is there density-dependence in survival of Spring migrants to John Day Dam?
par(mfrow=c(1,2))
plot(Z$SpringMigrants,Z$SpringSurvival,pch=21,cex=1.2,bg='red',xlab='SpringMigrants',ylab='SpringSurvival')
text(Z$SpringMigrants~Z$SpringSurvival,labels=Z$BroodYear,cex=0.9,font=1,pos=4)
plot(Z$TotalMigrants,Z$SpringSurvival,pch=21,cex=1.2,bg='red',xlab='TotalMigrants',ylab='SpringSurvival')
text(Z$TotalMigrants~Z$SpringSurvival,labels=Z$BroodYear,cex=0.9,font=1,pos=4)
Answer: No
6. Can we predict observed spawners using (i) smolt abundance at JD, (ii)CWT-based SARs, and (iii) spawner age?
First get recruit age composition
RAge<-matrix(NA,11,3)
NAge<-Z$TotalSpawners*Z[,5:7]
NAge<-Z$TotalSpawners*Z[,5:7]
for (t in 1:11){
RAge[t,1]<-NAge[t+3,1]/(NAge[t+3,1]+NAge[t+4,2]+NAge[t+5,3])
RAge[t,2]<-NAge[t+4,2]/(NAge[t+3,1]+NAge[t+4,2]+NAge[t+5,3])
RAge[t,3]<-NAge[t+5,3]/(NAge[t+3,1]+NAge[t+4,2]+NAge[t+5,3])
}
Now do the arithmetic to make spawenrs-at-age
pred<-matrix(NA,16,4)
pred[1:16,1]<-seq(2002,2017)
for (t in 1:(16-5)){
pred[t+3,2]<-Z$Smolts_JDDam[t]*Z$SAR[t]*RAge[t,1]
pred[t+4,3]<-Z$Smolts_JDDam[t]*Z$SAR[t]*RAge[t,2]
pred[t+5,4]<-Z$Smolts_JDDam[t]*Z$SAR[t]*RAge[t,3]
}
Now plot the predictive performance
Spawner.hat<-cbind(seq(2007,2015), rowSums(pred[6:14,2:4]), Z$TotalSpawners[6:14])
plot(Spawner.hat[,3],Spawner.hat[,2],pch=21,cex=1.2,bg='red',xlab='Observed Spawner Abudnance',ylab='Predicted Spawner Abudnance',xlim=c(0,2000),ylim=c(0,2000))
abline(0,1)
Slightly biased toward underpredicting. Create a correction factor.
cf<-lm(Spawner.hat[,2]~ -1 + Spawner.hat[,3])
Now look at the predictive ability using the correction factor
plot(Spawner.hat[,3],Spawner.hat[,2]/cf$coefficients,pch=21,cex=1.2,bg='red',xlab='Observed Spawner Abudnance',ylab='Predicted Spawner Abudnance',xlim=c(0,2000),ylim=c(0,2000))
abline(0,1)
This is pretty good.
7. Build preliminaries for a basic life-cycle simulation
Here are the stages of the life-cycle model:
- Spawners to TotalMigrants: Beverton-Holt model with covaraite. Include random residual error.
- TotalMigrants to Smolts at JD Dam: Random proportion migrating in Fall vs.ย Spring; random survival of Spring migrants.
Smolts at JD Dam to Sapwners: Random SAR, random age composition, fixed correction factor.
Spawner to TotalMigrants
TM<-function(S){
exp(rnorm(1,log(BH2fit$par[1]) + log(S) -log(1+BH2fit$par[2]*S) +BH2fit$par[3]*rnorm(1,0,1),BH2fit$par[4]))
}
2a. Get random proportion of migrating in Fall. Estimate theoretical distribution so proportion can be more extrement than observations.
beta<-function(data,params){
-sum(dbeta(data,params[1],params[2],log=TRUE))
}
betafit<-optim(par=c(1,1),fn=beta,data=pMF[!is.na(pMF)])
rpMP<-function(){
rbeta(1,betafit$par[1],betafit$par[2])
}
2b. Get Smolts at JD Dam. Estimate theoretical distributuion so Spring survivals can be more extreme than observations.
Ssurv<-function(data,params){
-sum(dbeta(data,params[1],params[2],log=TRUE))
}
Ssurfit<-optim(par=c(1,1),fn=Ssurv,data=Z$SpringSurvival[1:15])
Smolts<-function(TM){
rpMPi<-rpMP()
rpMPi*TM*0.28 + (1-rpMPi)*TM*rbeta(1,Ssurfit$par[1],Ssurfit$par[2])
}
3a. Estimate theoretical distribution so SAR can be more extreme than observations.
dSAR<-function(data,params){
-sum(dbeta(data,params[1],params[2],log=TRUE))
}
SARfit<-optim(par=c(1,1),fn=dSAR,data=Z$SAR[!is.na(Z$SAR)])
## Warning in dbeta(data, params[1], params[2], log = TRUE): NaNs produced
SAR<-function(){
rbeta(1,SARfit$par[1],SARfit$par[2])
}
8. Conduct basic life-cycle simulation
tAge<-RAge[sample(1:11,100,replace=TRUE),]#randomize existing emprical recruit age composition
N<-matrix(NA,105,3)
TotMig<-vector()
NJD<-vector()
N[1:5,1:3]<-as.matrix(Z$TotalSpawners[12:16]*Z[12:16,5:7])#Initialize
for (t in 1:100){
S<-sum(N[t,])
TotMig[t]<-TM(S)
NJD[t]<-Smolts(TotMig[t])
R<-NJD[t]*SAR()*tAge[t,]/cf$coefficients
#R<-Smolts(TM(S))*SAR()*tAge[t,]/cf$coefficients
N[t+3,1]<-R[1]
N[t+4,2]<-R[2]
N[t+5,3]<-R[3]
}
plot simulated Spawner abundnaces at age
plot(1:100,N[1:100,2],xlab='Years',ylab='Abundance',col=2,type='l')
lines(1:100,N[1:100,1],col=3)
lines(1:100,N[1:100,3],col=4)
legend('topleft',legend=c('Age4','Age3','Age5'),lty=c(1,1,1),col=c(2,3,4))
Compare simulated abudnance of TotalMigrants with empricial:
plot(1:100,TotMig,xlab='Years',ylab='Total Migrants',col=1,type='l')
lines(40:55,Z$TotalMigrants,col=2)
legend('topleft',legend=c('Simulated','Observed'),lty=c(1,1),col=c(1,2))
In the figure above, we hope the mean and varaince of the red line โlooks likeโ the mean and variance of the black line. We expect to see more extreme simulated values (black) than emprical (red) simply because there are more realizations of the simulated abudnances. This looks good.
Now compare simulated smolt abudnance at JD Dam with empirical:
plot(1:100,NJD,xlab='Years',ylab='Smolts at JD Dam',col=1,type='l')
lines(40:55,Z$Smolts_JDDam,col=2)
legend('topleft',legend=c('Simulated','Observed'),lty=c(1,1),col=c(1,2))
This alos looks good.
Finally, compare the simulated TotalSpawners with emprical:
plot(1:100,rowSums(N[1:100,]),xlab='Years',ylab='Total Spawners',col=1,type='l')
lines(40:55,Z$TotalSpawners,col=2)
legend('topleft',legend=c('Simulated','Observed'),lty=c(1,1),col=c(1,2))
Excellent. The life-cycle model is nicely capturing the central location and randomness in all three life stages.
9. Population Viability Analysis (PVA)
9.a No change in temperature
A critical uncertainty in PVA is the Reproductive Failure Threshold (RFT). This is the abundance of spawners below wich no production occurs. This may arise from Allee effects (mate finding). Another cricical uncertainty in PVA is the Quasi-Extinction Threshold (QET). This is the abundance below which we declare extinction. This number is greater than zero because the the PVA does not inclde genetic drift and inbreeding depression. Think of it as the event horizon of the extinction vortex. Because RFT and QET are unknown, I will vary these across multiple levels and repeat the PVA for all combinations of RFT and QET.
The settings on the PVA below (reps and levels of QET and RFT) will take about 1 hour and 20 minutes to complete.
PVAout<-matrix(NA,16,16)
Reps<-10000
ptm<-proc.time()
for (QETi in 1:16){
for (RFTi in 1:16){
QET<-40+QETi*10
RFT<-40+RFTi*10
QET3<-0
for (rep in 1:Reps){
for (t in 1:100){
qet3<-0
S<-sum(N[t,])
if (S<RFT){
S=0
}
R<-Smolts(TM(S))*SAR()*tAge[t,]/cf$coefficients
N[t+3,1]<-R[1]
N[t+4,2]<-R[2]
N[t+5,3]<-R[3]
}#t
# 3 consecutive years below QET
u<-as.numeric(rowSums(N[1:100,])<QET)
for (i in 1:(length(u)-4)){
qet3<-qet3+u[i]*u[i+1]*u[i+2]
}
QET3<-QET3+(qet3>0)
}#rep
PVAout[QETi,RFTi]<-QET3/Reps
}#RFTi
}#QETi
proc.time() - ptm
## user system elapsed
## 4255.64 0.58 4262.80
x<-as.vector(t(replicate(16,seq(50,200,by=10))))
y<-as.vector(replicate(16,seq(50,200,by=10)))
z<-as.vector(PVAout)
library(akima)
zz<-interp(x,y,z,nx=50,ny=50)
#filled.contour(zz,nlevels=50,color.palette=colorRampPalette(c('white','blue','yellow','red','darkred')),xlab="RFT",ylab="QET",main="Pr[extirpation] over 100 years of Coho at Huntley")
#contour(zz,add=T,axes=T,frame.plot=T)
filled.contour(zz,nlevels=50,frame.plot=T,key.title=title(main="Prob"),color.palette=colorRampPalette(c('white','blue','yellow','red','darkred')),xlab="RFT",ylab="QET",main="Viabilty of Middle Fork John Day Spring Chinook",
plot.axes={contour(zz,add=T,levels=c(0.01,0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.7),drawlabels=T,lwd=2,labcex=1);
axis(1);
axis(2)} )
The โProbโ plotted in the figure is the probablity of spawenr abundance dropping below QET for any 3 consecutive years over a 100year period.
9.b Increased Temperature
This PVA will assume that temperatures have increased by 1 standard deviation relative to the empricial data.
Redefine the following function defining Spawner to TotalMigrants
TM1<-function(S){
exp(rnorm(1,log(BH2fit$par[1]) + log(S) -log(1+BH2fit$par[2]*S) +BH2fit$par[3]*rnorm(1,1,1),BH2fit$par[4]))
}
Then just rerun the PVA by substituting TM1() for TM()
PVAout<-matrix(NA,16,16)
Reps<-10000
ptm<-proc.time()
for (QETi in 1:16){
for (RFTi in 1:16){
QET<-40+QETi*10
RFT<-40+RFTi*10
QET3<-0
for (rep in 1:Reps){
for (t in 1:100){
qet3<-0
S<-sum(N[t,])
if (S<RFT){
S=0
}
R<-Smolts(TM1(S))*SAR()*tAge[t,]/cf$coefficients
N[t+3,1]<-R[1]
N[t+4,2]<-R[2]
N[t+5,3]<-R[3]
}#t
# 3 consecutive years below QET
u<-as.numeric(rowSums(N[1:100,])<QET)
for (i in 1:(length(u)-4)){
qet3<-qet3+u[i]*u[i+1]*u[i+2]
}
QET3<-QET3+(qet3>0)
}#rep
PVAout[QETi,RFTi]<-QET3/Reps
}#RFTi
}#QETi
proc.time() - ptm
## user system elapsed
## 4292.75 0.44 4301.58
x<-as.vector(t(replicate(16,seq(50,200,by=10))))
y<-as.vector(replicate(16,seq(50,200,by=10)))
z<-as.vector(PVAout)
library(akima)
zz<-interp(x,y,z,nx=50,ny=50)
#filled.contour(zz,nlevels=50,color.palette=colorRampPalette(c('white','blue','yellow','red','darkred')),xlab="RFT",ylab="QET",main="Pr[extirpation] over 100 years of Coho at Huntley")
#contour(zz,add=T,axes=T,frame.plot=T)
filled.contour(zz,nlevels=50,frame.plot=T,key.title=title(main="Prob"),color.palette=colorRampPalette(c('white','blue','yellow','red','darkred')),xlab="RFT",ylab="QET",main="Viabilty of Middle Fork John Day Spring Chinook",
plot.axes={contour(zz,add=T,levels=c(0.01,0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.7),drawlabels=T,lwd=2,labcex=1);
axis(1);
axis(2)} )
Compare this PVA output with the first PVA output to see the effect of warming.