Rogue Spring Chinook Spawner-Recruit Analysis
Matt Falcy
July 11, 2017
This analysis is a spawner-Recruit analysis of Rogue River Spring Chinook. The purpose of this document is describe the steps taken to obtain results and facilitate reproducibility.
First, look at the data received from the Upper Rogue District.
(SR<-read.csv('ChS for Matt_SR.csv',header=T))
## Brood Wild Hatchery Spawners Age2 Age3 Age4 Age5 Age6 Recruits
## 1 1981 9545 176 9721 670 2595 2505 3487 464 9966
## 2 1982 14208 296 14504 3814 2467 7312 764 147 12489
## 3 1983 6941 118 7059 504 1487 3878 1190 0 33680
## 4 1984 5080 109 5189 1154 573 3188 257 16 45945
## 5 1985 17478 474 17952 3655 2943 7775 3579 0 22126
## 6 1986 25645 1509 27154 6387 12126 7391 1113 135 11276
## 7 1987 24330 1863 26193 2890 7723 14580 1000 0 3992
## 8 1988 42624 1070 43694 2595 9970 27810 3319 0 3909
## 9 1989 12526 2021 14547 777 1705 8185 3880 0 5767
## 10 1990 7491 730 8221 454 806 5736 1081 145 7735
## 11 1991 3435 402 3837 85 630 1882 1241 0 14478
## 12 1992 1391 190 1581 206 95 835 445 0 11610
## 13 1993 6254 674 6928 362 1191 3488 1805 83 10474
## 14 1994 4197 410 4607 29 333 2108 1952 186 6374
## 15 1995 18945 2842 21786 1934 3838 11707 4159 148 6156
## 16 1996 9295 1195 10490 931 1848 5637 2002 71 5976
## 17 1997 9599 1288 10887 967 1918 5850 2078 74 8693
## 18 1998 3684 491 4175 371 735 2243 797 28 10794
## 19 1999 5952 601 6553 582 1154 3521 1251 45 16044
## 20 2000 3443 1073 4516 401 796 2427 862 31 13236
## 21 2001 9339 957 10296 914 1814 5533 1965 70 8689
## 22 2002 6987 1632 8619 765 1518 4631 1645 59 6109
## 23 2003 19270 903 20173 1791 3554 10840 3851 137 5296
## 24 2004 13255 1040 14295 1269 2518 7681 2729 97 5400
## 25 2005 5803 491 6294 559 1109 3382 1201 43 5282
## 26 2006 4763 278 5041 448 888 2709 962 34 11252
## 27 2007 3465 1695 5160 458 909 2773 985 35 9868
## 28 2008 3970 1458 5428 482 956 2917 1036 37 13736
## 29 2009 5234 666 5900 524 1039 3171 1126 40 15690
## 30 2010 9556 925 10481 139 1972 7730 641 0 6927
## 31 2011 9940 390 10330 290 2397 5525 2001 56 15685
## 32 2012 14400 819 15219 122 3889 9439 1771 0 0
## 33 2013 12147 572 12719 76 868 9955 1761 76 0
## 34 2014 5593 771 6364 181 1324 3304 1556 0 0
## 35 2015 15320 245 15565 163 1781 10987 2633 0 0
## 36 2016 9573 66 9639 319 1554 4470 3298 0 0
It’s not yet clear how Recruits were computed. First, I’ll check if hatchery fish were included in the age break down
sum(SR[1,5:9])
## [1] 9721
Yes, hatchery fish are in the age break down. I know this because summing across the ages (9721) equals Wild + Hatchery in year 1 (1981). Now I’ll see if the Recruits reported here are just the diagonal sum of the abudnance-by-age data.
t<-1
sum(diag(as.matrix(SR[(t+2):(t+6),5:9])))
## [1] 9965
Yes, within rounding error these “Recruits” include hatchery fish. I want recruits of naturally spawning fish (exclude hatchery-origin fish) so I’m going to recompute recruits myself. First thing I need to do is get the age composition of spawners from the these data.
Age<-SR[,5:9]/rowSums(as.matrix(SR[,5:9]))
Now get the harvest rate data because I’m going to define recruits at a pre-freshwater fishing stage.
(HR<-read.csv('Harvest.csv'))
## Year LowerHR UpperHR
## 1 1961 0.138967979 0.11956522
## 2 1962 0.113953545 0.11117899
## 3 1963 0.145538820 0.09354527
## 4 1964 0.098098307 0.14174194
## 5 1965 0.133319707 0.16674788
## 6 1966 0.123046087 0.23510469
## 7 1967 0.189965244 0.21978953
## 8 1968 0.161824786 0.35655296
## 9 1969 0.147498322 0.24266202
## 10 1970 0.127372066 0.16969831
## 11 1971 0.218967563 0.15347342
## 12 1972 0.189110939 0.17270414
## 13 1973 0.087844629 0.17028629
## 14 1974 0.255348469 0.19283321
## 15 1975 0.163040588 0.25003296
## 16 1976 0.153161052 0.25061089
## 17 1977 0.113101464 0.31559872
## 18 1978 0.079071381 0.18998952
## 19 1979 0.125727029 0.12829906
## 20 1980 0.082025129 0.22574848
## 21 1981 0.120170300 0.22032976
## 22 1982 0.125870930 0.17140627
## 23 1983 0.124931513 0.15126707
## 24 1984 0.067145396 0.11077460
## 25 1985 0.059759073 0.10573412
## 26 1986 0.065713873 0.11518547
## 27 1987 0.036998977 0.13836083
## 28 1988 0.081598885 0.15000000
## 29 1989 0.115931834 0.15000000
## 30 1990 0.147560182 0.15000000
## 31 1991 0.157773924 0.15000000
## 32 1992 0.130314216 0.15000000
## 33 1993 0.244009880 0.04000000
## 34 1994 0.109188069 0.08000000
## 35 1995 0.151376579 0.02000000
## 36 1996 0.259009392 0.00000000
## 37 1997 0.131821625 0.00000000
## 38 1998 0.196431648 0.00000000
## 39 1999 0.147017124 0.00000000
## 40 2000 0.131197582 0.00000000
## 41 2001 0.089373417 0.00000000
## 42 2002 0.098559164 0.00000000
## 43 2003 0.185327805 0.00000000
## 44 2004 0.089727872 0.00000000
## 45 2005 0.203719728 0.00000000
## 46 2006 0.135423192 0.00000000
## 47 2007 0.102088856 0.00000000
## 48 2008 0.000985653 0.00000000
## 49 2009 0.011006977 0.00000000
## 50 2010 0.059942993 0.00000000
## 51 2011 0.040681838 0.00000000
## 52 2012 0.043183391 0.00000000
## 53 2013 0.032690326 0.00000000
## 54 2014 0.025960151 0.00000000
## 55 2015 0.081114707 0.00000000
## 56 2016 NA 0.00000000
HR[56,2]<-mean(HR[53:55,2])#Fill-in missing harvest rate data with recent observations
Get wild, preharvest abudnance at age through time.
Wentry<-SR$Wild/((1-HR[21:56,2])*(1-HR[21:56,3]))
WA<-Wentry*Age #abundance of wild fish at age
Now I can compute Recruits associated with each spawning year
Recruits<-vector()
for (t in 1:31){
Recruits[t]<-sum(diag(as.matrix(WA[(t+2):(t+6),])),na.rm=T)#Assume 0 6-year olds 1n 2017 and take Recruits to 2011
}
Plot the difference between my new Recruit data (X-axis) an the existing Recruit data.
plot(Recruits,SR$Recruits[1:31])# small difference wrt orignial computation of "Recruits."
Plot the spawner-recruit data using my new Recruit data
plot(SR$Spawners[1:31],Recruits,xlab='Spawners',xlim=c(0,50000),ylim=c(0,50000))
abline(0,1)
Here are the new Recruit data
(Recruits)
## [1] 11617.998 14333.836 38555.598 55323.341 26432.203 13960.519 4889.520
## [8] 4803.797 6881.319 8502.772 15530.804 12927.850 11096.681 6884.187
## [15] 6364.622 5775.954 8469.501 10938.029 17815.910 14107.790 9658.572
## [22] 6301.092 4485.865 4314.992 4564.922 10800.361 9783.259 13622.618
## [29] 15357.098 6777.612 16473.543
Get environmental covariates and then standardize them to a z-score.
(Env<-read.csv('Enviros.csv',header=T))
## SpawnYear SnowMaySept NPGO PDO Logerwell
## 1 1981 NA -4.2475335 0.45 109
## 2 1982 NA -0.4714075 10.43 126
## 3 1983 NA 4.1503529 1.94 112
## 4 1984 NA -4.0124191 2.50 48
## 5 1985 NA -4.1887797 3.87 89
## 6 1986 NA 1.9730959 9.86 81
## 7 1987 NA 8.9062497 2.40 68
## 8 1988 NA 3.2485294 1.80 97
## 9 1989 345 -1.6725678 1.64 81
## 10 1990 420 -3.6522981 -1.07 99
## 11 1991 265 -7.2568509 6.97 123
## 12 1992 605 -12.2753825 11.07 161
## 13 1993 310 -6.0855296 -0.40 87
## 14 1994 570 -5.2425061 5.81 95
## 15 1995 570 -5.7675932 4.15 120
## 16 1996 570 -3.3982884 11.92 146
## 17 1997 485 1.6454235 -0.37 105
## 18 1998 810 8.5242468 -5.13 91
## 19 1999 425 10.2024408 -3.58 72
## 20 2000 270 10.5842462 -4.22 61
## 21 2001 395 7.1067948 -0.26 80
## 22 2002 415 4.8042687 3.42 112
## 23 2003 440 1.8509556 2.96 110
## 24 2004 355 -8.2001023 3.48 145
## 25 2005 645 -3.3779722 0.28 112
## 26 2006 390 6.8078528 0.91 74
## 27 2007 620 10.8123831 -7.63 89
## 28 2008 500 4.8410953 -1.11 85
## 29 2009 465 7.0594872 -3.53 100
## 30 2010 700 6.0914540 -6.45 100
## 31 2011 505 10.1559965 -7.79 121
## 32 2012 380 1.1754769 -3.47 100
Env<-Env[-32,]
#Z-score the Enviros
mu<-colMeans(Env[,2:5],na.rm='T')
sig<-apply(Env[,2:5],2,sd,na.rm='T')
Zenv<-matrix(NA,31,4)
for (i in 1:4){
Zenv[,i]<-(Env[,i+1]-mu[i])/sig[i]
}
That’s all the data I need to fit stock-recruitment models. Now I’m going to write three JAGS files that specifiy the stock-recruitment models I want to fit. This first model (Mod1) is a Ricker model that won’t accept any covariates.
Mod1<-function(){
for (t in 1:T){
logR[t]~dnorm(mu[t],tau)
mu[t]<-log(alpha)+log(S[t])-(beta*S[t])#This is Ricker model
density[t]<-dnorm(logR[t],mu[t],tau)# Want this for WAIC calculation
}
#Derived parameters of interest
Smsy<-log(alpha)/beta*(0.5-0.07*log(alpha))
Hmsy<-0.5*log(alpha)-0.07*(log(alpha))^2
#priors
alpha~dunif(1,100)
beta~dunif(0,0.01)
tau<-pow(sigma,-2)
sigma~dunif(0,3)
sig<-sqrt(1/tau)
}
This is an alternative Ricker model that accpets a single environmental covariate. Subsequently, I’ll call this model three different times to test three different covariates (NPGO, PDO, Logerwell).
Mod2<-function(){
for (t in 1:T){
logR[t]~dnorm(mu[t],tau)
mu[t]<-log(alpha)+log(S[t])-(beta*S[t])+gamma1*env[t]#This is Ricker model
density[t]<-dnorm(logR[t],mu[t],tau)# Want this for WAIC calculation
}
#Derived parameters of interest
Smsy<-log(alpha)/beta*(0.5-0.07*log(alpha))
Hmsy<-0.5*log(alpha)-0.07*(log(alpha))^2
#priors
alpha~dunif(1,100)
beta~dunif(0,0.01)
gamma1~dnorm(0,100)
tau<-pow(sigma,-2)
sigma~dunif(0,3)
sig<-sqrt(1/tau)
}
Finally, this third model is also Ricker that accepts a single covariate. However, this is different than the model above because it will take a covariate with missing values (the snowpack covariate is missing its first eight values).
Mod3<-function(){
for (t in 1:8){
env[t]~dnorm(0,1)
}
for (t in 1:T){
logR[t]~dnorm(mu[t],tau)
mu[t]<-log(alpha)+log(S[t])-(beta*S[t])+gamma1*env[t]#This is Ricker model
density[t]<-dnorm(logR[t],mu[t],tau)# Want this for WAIC calculation
}
#Derived parameters of interest
Smsy<-log(alpha)/beta*(0.5-0.07*log(alpha))
Hmsy<-0.5*log(alpha)-0.07*(log(alpha))^2
#priors
alpha~dunif(1,100)
beta~dunif(0,0.01)
gamma1~dnorm(0,100)
tau<-pow(sigma,-2)
sigma~dunif(0,3)
sig<-sqrt(1/tau)
}
Set up preliminaries for JAGS. “gamma1” is the coefficient on the covariate. So if Mod1 (no covariate) is desired, then “gamma1” should be deleted from the “Init.fun”" and “params” below. Further, the last number of the line at the bottom of this chunk specifies which environmental covariate to use. Delete this if using Mod1.
nc<-4 #Number of chains
#Initialize MCMC chains
Init.fun<-function(){
list("alpha"=runif(1,1,5),"beta"=runif(1,0.0001,0.001),"sigma"=runif(1,0,3),"gamma1"=rnorm(1,0,0.3))
}
# List the parameters to save as output
params=as.character(c("alpha","beta","gamma1","sig","Smsy","Hmsy","density"))
#Repackage empirical data
D<-list(logR=log(Recruits),S=SR$Spawners[1:31],T=31,env=Zenv[,3])
Call up JAGS to run the MCMC. This model fits very easy and fast, but other applications may need to increase burnin, iterations, or thining.
#Call JAGS into a parallel process
library(R2jags)
## Loading required package: rjags
## Loading required package: coda
## Linked to JAGS 4.2.0
## Loaded modules: basemod,bugs
##
## Attaching package: 'R2jags'
## The following object is masked from 'package:coda':
##
## traceplot
jags.out<-jags.parallel(data=D, inits=Init.fun, parameters.to.save=params, model.file = Mod2,
n.chains = nc,n.iter = 100000, n.burnin = 5000,
n.thin = 13,
n.cluster= nc, DIC = TRUE,
working.directory = NULL, jags.seed = 123, digits=5,
jags.module = c("dic"))
Now inspect the output from JAGS to ensure that the model converged nicely.
## Mean SD Naive SE Time-series SE
## alpha 3.776710e+00 8.188375e-01 4.789261e-03 4.932776e-03
## beta 1.004777e-04 1.482411e-05 8.670405e-08 8.865780e-08
## density[1] 5.517362e-01 7.636725e-02 4.466609e-04 4.466820e-04
## density[10] 4.633546e-01 7.086382e-02 4.144721e-04 4.149406e-04
## density[11] 4.501376e-01 8.718779e-02 5.099486e-04 5.115484e-04
## density[12] 2.493572e-01 1.060263e-01 6.201321e-04 6.252289e-04
## density[13] 5.524269e-01 7.754695e-02 4.535607e-04 4.535798e-04
## density[14] 4.476035e-01 8.326615e-02 4.870115e-04 4.923147e-04
## density[15] 4.840569e-01 8.556692e-02 5.004684e-04 4.977737e-04
## density[16] 2.498342e-01 9.026829e-02 5.279660e-04 5.279741e-04
## density[17] 4.551768e-01 6.716757e-02 3.928533e-04 3.924122e-04
## density[18] 5.470387e-01 7.803259e-02 4.564012e-04 4.564145e-04
## density[19] 4.815634e-01 7.425690e-02 4.343177e-04 4.319640e-04
## density[2] 5.486665e-01 7.885436e-02 4.612076e-04 4.788128e-04
## density[20] 5.012721e-01 7.820343e-02 4.574004e-04 4.618324e-04
## density[21] 5.046802e-01 7.135129e-02 4.173233e-04 4.161294e-04
## density[22] 3.098411e-01 6.505094e-02 3.804734e-04 3.823635e-04
## density[23] 2.973521e-01 8.193761e-02 4.792411e-04 4.792154e-04
## density[24] 1.746536e-01 5.737865e-02 3.355993e-04 3.382806e-04
## density[25] 2.096069e-01 6.339788e-02 3.708049e-04 3.752028e-04
## density[26] 5.591354e-01 7.812996e-02 4.569707e-04 4.569868e-04
## density[27] 5.456314e-01 8.099389e-02 4.737214e-04 4.707317e-04
## density[28] 5.422750e-01 7.562548e-02 4.423224e-04 4.489121e-04
## density[29] 5.176079e-01 7.532539e-02 4.405672e-04 4.451956e-04
## density[3] 1.672762e-01 5.959723e-02 3.485754e-04 3.510978e-04
## density[30] 3.643349e-01 8.338514e-02 4.877075e-04 4.823037e-04
## density[31] 5.199118e-01 7.945331e-02 4.647108e-04 4.686333e-04
## density[4] 5.240140e-02 3.641495e-02 2.129857e-04 2.154216e-04
## density[5] 2.594158e-01 7.243839e-02 4.236815e-04 4.211826e-04
## density[6] 3.328471e-01 1.201858e-01 7.029493e-04 7.067479e-04
## density[7] 4.775981e-01 9.857697e-02 5.765622e-04 5.736542e-04
## density[8] 2.723875e-01 1.679758e-01 9.824659e-04 9.894184e-04
## density[9] 3.908767e-01 6.545150e-02 3.828162e-04 3.820820e-04
## deviance 6.504787e+01 2.812219e+00 1.644826e-02 1.625565e-02
## gamma1 2.530222e-02 8.005542e-02 4.682324e-04 4.647415e-04
## Hmsy 5.305210e-01 6.734180e-02 3.938723e-04 4.029503e-04
## sig 7.089442e-01 9.943788e-02 5.815976e-04 5.879011e-04
## Smsy 5.318183e+03 5.061969e+02 2.960671e+00 2.960618e+00
## [[1]]
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## alpha beta density[1] density[10] density[11] density[12]
## 0.25134 0.33341 0.99351 0.64293 0.20359 0.25953
## density[13] density[14] density[15] density[16] density[17] density[18]
## 0.85198 -0.07700 0.73742 -0.57088 0.72360 0.42897
## density[19] density[2] density[20] density[21] density[22] density[23]
## 0.03565 0.63871 0.10638 0.90425 -0.33124 0.08427
## density[24] density[25] density[26] density[27] density[28] density[29]
## -0.56463 -0.33967 0.71411 0.59544 0.53941 0.28085
## density[3] density[30] density[31] density[4] density[5] density[6]
## -0.49299 0.48106 0.15035 -0.09050 -0.71739 -0.38965
## density[7] density[8] density[9] deviance gamma1 Hmsy
## 0.72612 -0.62869 0.40563 1.21082 0.79637 0.03987
## sig Smsy
## -0.44225 -0.31690
##
##
## [[2]]
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## alpha beta density[1] density[10] density[11] density[12]
## -0.08826 -0.37309 0.78980 0.20738 -0.72797 -1.18222
## density[13] density[14] density[15] density[16] density[17] density[18]
## 0.50006 0.85710 1.44157 1.18924 0.46720 0.53236
## density[19] density[2] density[20] density[21] density[22] density[23]
## 0.81045 0.42043 0.75479 0.63031 0.58838 0.36104
## density[24] density[25] density[26] density[27] density[28] density[29]
## -0.09195 -0.09726 0.44942 -0.00751 0.50285 0.75234
## density[3] density[30] density[31] density[4] density[5] density[6]
## -0.62988 -0.71029 1.23573 -0.56993 -0.80159 -0.85934
## density[7] density[8] density[9] deviance gamma1 Hmsy
## 1.10722 -0.03707 0.76922 -0.19376 -1.36815 -0.51273
## sig Smsy
## -0.85620 -0.68439
##
##
## [[3]]
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## alpha beta density[1] density[10] density[11] density[12]
## -0.66680 -1.20263 -0.86472 -0.24268 -1.76376 -1.64421
## density[13] density[14] density[15] density[16] density[17] density[18]
## -0.73200 1.04589 -1.25764 1.05462 -0.52725 -0.79451
## density[19] density[2] density[20] density[21] density[22] density[23]
## -0.55371 -0.71531 -0.73201 -0.67263 0.92090 -0.83537
## density[24] density[25] density[26] density[27] density[28] density[29]
## 0.21412 0.89781 -0.82748 -0.66052 -1.12683 -0.79294
## density[3] density[30] density[31] density[4] density[5] density[6]
## 0.05222 -0.81977 0.09274 0.30179 1.02577 0.63396
## density[7] density[8] density[9] deviance gamma1 Hmsy
## -1.53652 1.56106 -0.58163 0.45298 -1.05700 -0.80060
## sig Smsy
## 0.35551 0.84487
##
##
## [[4]]
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## alpha beta density[1] density[10] density[11] density[12]
## -0.84905 -0.10038 0.57565 1.20940 0.13055 0.06537
## density[13] density[14] density[15] density[16] density[17] density[18]
## 0.52776 0.35262 0.51638 -0.73244 1.22931 -0.32791
## density[19] density[2] density[20] density[21] density[22] density[23]
## -1.23587 0.52986 -1.01749 1.01420 0.59097 0.48155
## density[24] density[25] density[26] density[27] density[28] density[29]
## 0.33808 0.79497 0.28368 0.62701 -0.37753 -0.95570
## density[3] density[30] density[31] density[4] density[5] density[6]
## -0.94561 1.76911 -0.92377 -0.62533 -0.83841 0.06448
## density[7] density[8] density[9] deviance gamma1 Hmsy
## 0.62275 -0.43608 0.98780 -0.05254 1.25453 -0.64439
## sig Smsy
## 0.04977 -0.74602
## [[1]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## Burn-in Total Lower bound Dependence
## (M) (N) (Nmin) factor (I)
## alpha 26 47944 3746 12.8
## beta 26 47944 3746 12.8
## density[1] 26 48490 3746 12.9
## density[10] 26 50713 3746 13.5
## density[11] 26 49036 3746 13.1
## density[12] 26 47411 3746 12.7
## density[13] 26 48490 3746 12.9
## density[14] 26 50713 3746 13.5
## density[15] 26 47944 3746 12.8
## density[16] 26 49036 3746 13.1
## density[17] 26 49582 3746 13.2
## density[18] 26 50713 3746 13.5
## density[19] 26 47944 3746 12.8
## density[2] 26 49582 3746 13.2
## density[20] 26 48490 3746 12.9
## density[21] 26 50713 3746 13.5
## density[22] 26 50141 3746 13.4
## density[23] 26 47944 3746 12.8
## density[24] 26 49036 3746 13.1
## density[25] 26 49036 3746 13.1
## density[26] 26 50141 3746 13.4
## density[27] 26 47944 3746 12.8
## density[28] 26 51285 3746 13.7
## density[29] 26 48490 3746 12.9
## density[3] 26 51857 3746 13.8
## density[30] 26 47944 3746 12.8
## density[31] 26 49582 3746 13.2
## density[4] 39 52455 3746 14.0
## density[5] 26 48490 3746 12.9
## density[6] 26 49036 3746 13.1
## density[7] 26 47411 3746 12.7
## density[8] 26 49036 3746 13.1
## density[9] 26 47411 3746 12.7
## deviance 26 47944 3746 12.8
## gamma1 26 48490 3746 12.9
## Hmsy 26 47944 3746 12.8
## sig 26 51285 3746 13.7
## Smsy 26 50141 3746 13.4
##
##
## [[2]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## Burn-in Total Lower bound Dependence
## (M) (N) (Nmin) factor (I)
## alpha 26 49582 3746 13.2
## beta 26 46891 3746 12.5
## density[1] 26 47411 3746 12.7
## density[10] 26 51285 3746 13.7
## density[11] 26 50141 3746 13.4
## density[12] 26 49582 3746 13.2
## density[13] 26 46891 3746 12.5
## density[14] 39 53053 3746 14.2
## density[15] 26 47944 3746 12.8
## density[16] 26 49582 3746 13.2
## density[17] 26 49582 3746 13.2
## density[18] 26 47411 3746 12.7
## density[19] 26 48490 3746 12.9
## density[2] 26 48490 3746 12.9
## density[20] 26 48490 3746 12.9
## density[21] 26 49582 3746 13.2
## density[22] 26 50713 3746 13.5
## density[23] 26 46891 3746 12.5
## density[24] 26 46358 3746 12.4
## density[25] 26 50141 3746 13.4
## density[26] 26 48490 3746 12.9
## density[27] 26 47944 3746 12.8
## density[28] 26 49036 3746 13.1
## density[29] 26 49036 3746 13.1
## density[3] 26 49036 3746 13.1
## density[30] 26 48490 3746 12.9
## density[31] 26 49036 3746 13.1
## density[4] 26 49582 3746 13.2
## density[5] 26 48490 3746 12.9
## density[6] 26 49036 3746 13.1
## density[7] 26 47944 3746 12.8
## density[8] 26 48490 3746 12.9
## density[9] 26 48490 3746 12.9
## deviance 26 50141 3746 13.4
## gamma1 26 49036 3746 13.1
## Hmsy 26 49582 3746 13.2
## sig 26 49582 3746 13.2
## Smsy 26 49582 3746 13.2
##
##
## [[3]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## Burn-in Total Lower bound Dependence
## (M) (N) (Nmin) factor (I)
## alpha 26 49582 3746 13.2
## beta 26 48490 3746 12.9
## density[1] 26 48490 3746 12.9
## density[10] 26 50141 3746 13.4
## density[11] 26 50141 3746 13.4
## density[12] 26 48490 3746 12.9
## density[13] 26 50713 3746 13.5
## density[14] 26 49582 3746 13.2
## density[15] 26 50141 3746 13.4
## density[16] 26 49608 3746 13.2
## density[17] 26 50713 3746 13.5
## density[18] 26 48490 3746 12.9
## density[19] 26 49582 3746 13.2
## density[2] 26 49036 3746 13.1
## density[20] 26 49582 3746 13.2
## density[21] 26 49036 3746 13.1
## density[22] 26 46891 3746 12.5
## density[23] 26 47944 3746 12.8
## density[24] 26 47424 3746 12.7
## density[25] 26 48490 3746 12.9
## density[26] 26 49582 3746 13.2
## density[27] 26 47411 3746 12.7
## density[28] 26 48490 3746 12.9
## density[29] 26 48490 3746 12.9
## density[3] 26 49036 3746 13.1
## density[30] 26 48490 3746 12.9
## density[31] 26 47944 3746 12.8
## density[4] 26 50141 3746 13.4
## density[5] 26 48490 3746 12.9
## density[6] 26 49582 3746 13.2
## density[7] 26 49036 3746 13.1
## density[8] 26 49036 3746 13.1
## density[9] 26 49582 3746 13.2
## deviance 26 47411 3746 12.7
## gamma1 26 47944 3746 12.8
## Hmsy 26 49582 3746 13.2
## sig 26 49036 3746 13.1
## Smsy 26 49036 3746 13.1
##
##
## [[4]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## Burn-in Total Lower bound Dependence
## (M) (N) (Nmin) factor (I)
## alpha 26 47944 3746 12.8
## beta 26 50713 3746 13.5
## density[1] 26 49036 3746 13.1
## density[10] 26 49036 3746 13.1
## density[11] 26 49582 3746 13.2
## density[12] 26 49036 3746 13.1
## density[13] 26 48490 3746 12.9
## density[14] 26 50141 3746 13.4
## density[15] 26 48789 3746 13.0
## density[16] 26 47944 3746 12.8
## density[17] 26 49036 3746 13.1
## density[18] 26 47944 3746 12.8
## density[19] 26 47411 3746 12.7
## density[2] 26 48490 3746 12.9
## density[20] 26 47411 3746 12.7
## density[21] 26 50713 3746 13.5
## density[22] 26 47411 3746 12.7
## density[23] 26 48490 3746 12.9
## density[24] 26 47411 3746 12.7
## density[25] 26 47944 3746 12.8
## density[26] 26 47944 3746 12.8
## density[27] 26 49036 3746 13.1
## density[28] 26 47411 3746 12.7
## density[29] 26 46358 3746 12.4
## density[3] 26 51857 3746 13.8
## density[30] 26 47411 3746 12.7
## density[31] 26 47944 3746 12.8
## density[4] 26 49582 3746 13.2
## density[5] 26 48490 3746 12.9
## density[6] 26 48490 3746 12.9
## density[7] 26 49036 3746 13.1
## density[8] 26 49036 3746 13.1
## density[9] 26 46891 3746 12.5
## deviance 26 48490 3746 12.9
## gamma1 26 46891 3746 12.5
## Hmsy 26 47944 3746 12.8
## sig 26 48490 3746 12.9
## Smsy 26 47944 3746 12.8
Yes, everyting looks great and all the models converge really fast. Now compute summary statitics and look at the posterior probability density for the covariate coefficient, “gamma1”. Replace “gamma1” with “Smsy” or “Hmsy” or “alpha” to look at different parameters.
jags.out$BUGSoutput$summary
## mean sd 2.5% 25%
## Hmsy 5.305210e-01 6.734180e-02 3.913585e-01 4.874152e-01
## Smsy 5.318183e+03 5.061969e+02 4.434719e+03 4.982304e+03
## alpha 3.776710e+00 8.188375e-01 2.446880e+00 3.205171e+00
## beta 1.004777e-04 1.482411e-05 7.146359e-05 9.064020e-05
## density[1] 5.517362e-01 7.636725e-02 4.076684e-01 4.993601e-01
## density[2] 5.486665e-01 7.885436e-02 3.990020e-01 4.947310e-01
## density[3] 1.672762e-01 5.959723e-02 6.368158e-02 1.236118e-01
## density[4] 5.240140e-02 3.641495e-02 6.945556e-03 2.506342e-02
## density[5] 2.594158e-01 7.243839e-02 1.248887e-01 2.077911e-01
## density[6] 3.328471e-01 1.201858e-01 1.105322e-01 2.439983e-01
## density[7] 4.775981e-01 9.857697e-02 2.706533e-01 4.148525e-01
## density[8] 2.723875e-01 1.679758e-01 1.766994e-02 1.284791e-01
## density[9] 3.908767e-01 6.545150e-02 2.636817e-01 3.460420e-01
## density[10] 4.633546e-01 7.086382e-02 3.284737e-01 4.144689e-01
## density[11] 4.501376e-01 8.718779e-02 2.767951e-01 3.913561e-01
## density[12] 2.493572e-01 1.060263e-01 6.763030e-02 1.700368e-01
## density[13] 5.524269e-01 7.754695e-02 4.050406e-01 4.995100e-01
## density[14] 4.476035e-01 8.326615e-02 2.843631e-01 3.915740e-01
## density[15] 4.840569e-01 8.556692e-02 3.132572e-01 4.271239e-01
## density[16] 2.498342e-01 9.026829e-02 8.898999e-02 1.841465e-01
## density[17] 4.551768e-01 6.716757e-02 3.272644e-01 4.086916e-01
## density[18] 5.470387e-01 7.803259e-02 3.982156e-01 4.930695e-01
## density[19] 4.815634e-01 7.425690e-02 3.388960e-01 4.308268e-01
## density[20] 5.012721e-01 7.820343e-02 3.499619e-01 4.478744e-01
## density[21] 5.046802e-01 7.135129e-02 3.702489e-01 4.555087e-01
## density[22] 3.098411e-01 6.505094e-02 1.857432e-01 2.646634e-01
## density[23] 2.973521e-01 8.193761e-02 1.436987e-01 2.396102e-01
## density[24] 1.746536e-01 5.737865e-02 7.315092e-02 1.331000e-01
## density[25] 2.096069e-01 6.339788e-02 9.639677e-02 1.638414e-01
## density[26] 5.591354e-01 7.812996e-02 4.115633e-01 5.057931e-01
## density[27] 5.456314e-01 8.099389e-02 3.894988e-01 4.908225e-01
## density[28] 5.422750e-01 7.562548e-02 3.997135e-01 4.900126e-01
## density[29] 5.176079e-01 7.532539e-02 3.745015e-01 4.655252e-01
## density[30] 3.643349e-01 8.338514e-02 2.044508e-01 3.079647e-01
## density[31] 5.199118e-01 7.945331e-02 3.659403e-01 4.662486e-01
## deviance 6.504787e+01 2.812219e+00 6.164461e+01 6.301744e+01
## gamma1 2.530222e-02 8.005542e-02 -1.335310e-01 -2.845592e-02
## sig 7.089442e-01 9.943788e-02 5.457647e-01 6.393382e-01
## 50% 75% 97.5% Rhat n.eff
## Hmsy 5.337880e-01 5.763006e-01 6.537735e-01 1.000969 29000
## Smsy 5.280439e+03 5.609372e+03 6.423706e+03 1.001051 21000
## alpha 3.693504e+00 4.241316e+00 5.603138e+00 1.000989 29000
## beta 1.003083e-04 1.102293e-04 1.296948e-04 1.001044 22000
## density[1] 5.497116e-01 6.017558e-01 7.088358e-01 1.000937 29000
## density[2] 5.470448e-01 6.006390e-01 7.100769e-01 1.000977 29000
## density[3] 1.633785e-01 2.071054e-01 2.923231e-01 1.000944 29000
## density[4] 4.401267e-02 7.107780e-02 1.444059e-01 1.000953 29000
## density[5] 2.576054e-01 3.086526e-01 4.067500e-01 1.001019 28000
## density[6] 3.326727e-01 4.192499e-01 5.617411e-01 1.001047 21000
## density[7] 4.825107e-01 5.451951e-01 6.606371e-01 1.001041 23000
## density[8] 2.576605e-01 4.071668e-01 5.936255e-01 1.001078 17000
## density[9] 3.900714e-01 4.348779e-01 5.205958e-01 1.001043 22000
## density[10] 4.623262e-01 5.104815e-01 6.062370e-01 1.000970 29000
## density[11] 4.508389e-01 5.098467e-01 6.179095e-01 1.000948 29000
## density[12] 2.420729e-01 3.214062e-01 4.698901e-01 1.000995 29000
## density[13] 5.505904e-01 6.033661e-01 7.107692e-01 1.000945 29000
## density[14] 4.473300e-01 5.033869e-01 6.111022e-01 1.001005 29000
## density[15] 4.847785e-01 5.420716e-01 6.497303e-01 1.001002 29000
## density[16] 2.455579e-01 3.114419e-01 4.365841e-01 1.000995 29000
## density[17] 4.543187e-01 5.000943e-01 5.904237e-01 1.000979 29000
## density[18] 5.456183e-01 5.988122e-01 7.049640e-01 1.000961 29000
## density[19] 4.808712e-01 5.310219e-01 6.281404e-01 1.001001 29000
## density[20] 5.009676e-01 5.534047e-01 6.550823e-01 1.000985 29000
## density[21] 5.028418e-01 5.521368e-01 6.492011e-01 1.000948 29000
## density[22] 3.090072e-01 3.537623e-01 4.390465e-01 1.000976 29000
## density[23] 2.957029e-01 3.527917e-01 4.614104e-01 1.001119 13000
## density[24] 1.716892e-01 2.129832e-01 2.944164e-01 1.001051 21000
## density[25] 2.064564e-01 2.517516e-01 3.402885e-01 1.000991 29000
## density[26] 5.574156e-01 6.106411e-01 7.178131e-01 1.000953 29000
## density[27] 5.439407e-01 5.997139e-01 7.091782e-01 1.000947 29000
## density[28] 5.406258e-01 5.920958e-01 6.966011e-01 1.000953 29000
## density[29] 5.165132e-01 5.678465e-01 6.693831e-01 1.000980 29000
## density[30] 3.635165e-01 4.198762e-01 5.314725e-01 1.001085 16000
## density[31] 5.194583e-01 5.729368e-01 6.785948e-01 1.001036 24000
## deviance 6.440476e+01 6.633834e+01 7.224261e+01 1.000982 29000
## gamma1 2.544000e-02 7.965241e-02 1.823404e-01 1.001075 17000
## sig 6.973710e-01 7.666768e-01 9.357520e-01 1.000942 29000
jags.out$BUGSoutput$mean
## $Hmsy
## [1] 0.530521
##
## $Smsy
## [1] 5318.183
##
## $alpha
## [1] 3.77671
##
## $beta
## [1] 0.0001004777
##
## $density
## [1] 0.5517362 0.5486665 0.1672762 0.0524014 0.2594158 0.3328471 0.4775981
## [8] 0.2723875 0.3908767 0.4633546 0.4501376 0.2493572 0.5524269 0.4476035
## [15] 0.4840569 0.2498342 0.4551768 0.5470387 0.4815634 0.5012721 0.5046802
## [22] 0.3098411 0.2973521 0.1746536 0.2096069 0.5591354 0.5456314 0.5422750
## [29] 0.5176079 0.3643349 0.5199118
##
## $deviance
## [1] 65.04787
##
## $gamma1
## [1] 0.02530222
##
## $sig
## [1] 0.7089442
out2<-as.mcmc(jags.out)
summary(mcmc.list(out2))$statistics
## Mean SD Naive SE Time-series SE
## alpha 3.776710e+00 8.188375e-01 4.789261e-03 4.932776e-03
## beta 1.004777e-04 1.482411e-05 8.670405e-08 8.865780e-08
## density[1] 5.517362e-01 7.636725e-02 4.466609e-04 4.466820e-04
## density[10] 4.633546e-01 7.086382e-02 4.144721e-04 4.149406e-04
## density[11] 4.501376e-01 8.718779e-02 5.099486e-04 5.115484e-04
## density[12] 2.493572e-01 1.060263e-01 6.201321e-04 6.252289e-04
## density[13] 5.524269e-01 7.754695e-02 4.535607e-04 4.535798e-04
## density[14] 4.476035e-01 8.326615e-02 4.870115e-04 4.923147e-04
## density[15] 4.840569e-01 8.556692e-02 5.004684e-04 4.977737e-04
## density[16] 2.498342e-01 9.026829e-02 5.279660e-04 5.279741e-04
## density[17] 4.551768e-01 6.716757e-02 3.928533e-04 3.924122e-04
## density[18] 5.470387e-01 7.803259e-02 4.564012e-04 4.564145e-04
## density[19] 4.815634e-01 7.425690e-02 4.343177e-04 4.319640e-04
## density[2] 5.486665e-01 7.885436e-02 4.612076e-04 4.788128e-04
## density[20] 5.012721e-01 7.820343e-02 4.574004e-04 4.618324e-04
## density[21] 5.046802e-01 7.135129e-02 4.173233e-04 4.161294e-04
## density[22] 3.098411e-01 6.505094e-02 3.804734e-04 3.823635e-04
## density[23] 2.973521e-01 8.193761e-02 4.792411e-04 4.792154e-04
## density[24] 1.746536e-01 5.737865e-02 3.355993e-04 3.382806e-04
## density[25] 2.096069e-01 6.339788e-02 3.708049e-04 3.752028e-04
## density[26] 5.591354e-01 7.812996e-02 4.569707e-04 4.569868e-04
## density[27] 5.456314e-01 8.099389e-02 4.737214e-04 4.707317e-04
## density[28] 5.422750e-01 7.562548e-02 4.423224e-04 4.489121e-04
## density[29] 5.176079e-01 7.532539e-02 4.405672e-04 4.451956e-04
## density[3] 1.672762e-01 5.959723e-02 3.485754e-04 3.510978e-04
## density[30] 3.643349e-01 8.338514e-02 4.877075e-04 4.823037e-04
## density[31] 5.199118e-01 7.945331e-02 4.647108e-04 4.686333e-04
## density[4] 5.240140e-02 3.641495e-02 2.129857e-04 2.154216e-04
## density[5] 2.594158e-01 7.243839e-02 4.236815e-04 4.211826e-04
## density[6] 3.328471e-01 1.201858e-01 7.029493e-04 7.067479e-04
## density[7] 4.775981e-01 9.857697e-02 5.765622e-04 5.736542e-04
## density[8] 2.723875e-01 1.679758e-01 9.824659e-04 9.894184e-04
## density[9] 3.908767e-01 6.545150e-02 3.828162e-04 3.820820e-04
## deviance 6.504787e+01 2.812219e+00 1.644826e-02 1.625565e-02
## gamma1 2.530222e-02 8.005542e-02 4.682324e-04 4.647415e-04
## Hmsy 5.305210e-01 6.734180e-02 3.938723e-04 4.029503e-04
## sig 7.089442e-01 9.943788e-02 5.815976e-04 5.879011e-04
## Smsy 5.318183e+03 5.061969e+02 2.960671e+00 2.960618e+00
out3<-do.call(rbind.data.frame,out2)
d<-density(out3$gamma1)
plot(d,xlab='gamma1',main='')
polygon(d,col='skyblue',border='black',lwd=1)
Compute WAIC for the model. Do this for Mod1, Mod2 with covs=c(NPGO, PDO, Logerwell) and Mod3 with cov=SnowMaySept. That’s five models to compare. It turns out that WAIC is very similar for all models. Must specify “Den” to match the 31 desnities. Other models may stick these between out3[5:35]
names(out3)
## [1] "alpha" "beta" "density[1]" "density[10]" "density[11]"
## [6] "density[12]" "density[13]" "density[14]" "density[15]" "density[16]"
## [11] "density[17]" "density[18]" "density[19]" "density[2]" "density[20]"
## [16] "density[21]" "density[22]" "density[23]" "density[24]" "density[25]"
## [21] "density[26]" "density[27]" "density[28]" "density[29]" "density[3]"
## [26] "density[30]" "density[31]" "density[4]" "density[5]" "density[6]"
## [31] "density[7]" "density[8]" "density[9]" "deviance" "gamma1"
## [36] "Hmsy" "sig" "Smsy"
Den<-out3[3:33]
lppd<-sum(log((1/dim(Den)[1])*colSums(Den)))
pwaic<-sum(diag(var(log(Den))))
(WAIC<--2*(lppd-pwaic))
## [1] 68.84242
Compute the 95% higheste probability density interval (HPDI).
#HPDinterval(mcmc(jags.out$BUGSoutput$sims.matrix[,37]),prop=0.95) #isolate a variable
HPDinterval(as.mcmc(out3),prop=0.95)# Do everything
## lower upper
## alpha 2.357866e+00 5.456534e+00
## beta 7.177892e-05 1.299963e-04
## density[1] 4.043279e-01 7.042776e-01
## density[10] 3.274733e-01 6.046680e-01
## density[11] 2.805214e-01 6.204926e-01
## density[12] 5.323347e-02 4.496747e-01
## density[13] 3.968165e-01 7.010752e-01
## density[14] 2.763533e-01 6.029875e-01
## density[15] 3.181956e-01 6.533969e-01
## density[16] 7.754263e-02 4.212680e-01
## density[17] 3.218260e-01 5.844038e-01
## density[18] 3.929280e-01 6.981015e-01
## density[19] 3.415654e-01 6.300355e-01
## density[2] 3.951378e-01 7.055952e-01
## density[20] 3.505128e-01 6.553741e-01
## density[21] 3.647591e-01 6.426802e-01
## density[22] 1.844754e-01 4.369881e-01
## density[23] 1.350425e-01 4.516779e-01
## density[24] 6.797977e-02 2.869581e-01
## density[25] 9.036855e-02 3.321078e-01
## density[26] 4.081404e-01 7.137888e-01
## density[27] 3.883465e-01 7.070864e-01
## density[28] 4.017059e-01 6.983661e-01
## density[29] 3.736820e-01 6.681505e-01
## density[3] 5.524163e-02 2.814785e-01
## density[30] 2.016935e-01 5.276779e-01
## density[31] 3.646583e-01 6.770468e-01
## density[4] 1.290203e-03 1.241640e-01
## density[5] 1.217654e-01 4.015567e-01
## density[6] 1.063051e-01 5.560265e-01
## density[7] 2.729310e-01 6.622921e-01
## density[8] 7.918456e-04 5.572591e-01
## density[9] 2.608393e-01 5.164452e-01
## deviance 6.125583e+01 7.049977e+01
## gamma1 -1.394176e-01 1.759585e-01
## Hmsy 3.912283e-01 6.535649e-01
## sig 5.224635e-01 9.015117e-01
## Smsy 4.404970e+03 6.381787e+03
## attr(,"Probability")
## [1] 0.9499863
Fancy interactive 3D plot of the covariance between alpha and beta. Click and drag the figure.
library(rgl)
library(MASS)
a<-out3$alpha
b<-out3$beta
biv<- kde2d(a, b, n = 50)
col <- heat.colors(length(biv$z))[rank(biv$z)]
persp3d(x=biv,col=col,xlab="Alpha",ylab="Beta",zlab='Frequency')
Plot the spawner-recruit model. In light of the covariance between alpha and beta we have a problem in describing the 95%credible interval for the expected recruitment function. A numerical solution is to get the HPDI of Recuits after sampling the MCMC samples thousands of times. This will take a few minutes.
plot(SR$Spawners[1:31],Recruits,xlim=c(0,max(c(Recruits,SR$Spawners[1:31]))),ylim=c(0,max(c(Recruits,SR$Spawners[1:31]))),xlab='Spawners',ylab='Recruits')
Sx<-seq(1,max(SR$Spawners[1:31]))
ax<-summary(mcmc.list(out2))$statistics[1]
bx<-summary(mcmc.list(out2))$statistics[2]
Rx<-(ax*Sx)*exp(-bx*Sx)
lines(Sx,Rx,type='l')
abline(0,1,col='2')
HPDI<-matrix(NA,max(SR$Spawners[1:31]),2)
for (t in 1:max(SR$Spawners[1:31])){
i<-sample(1:29232,10000,replace=T)
ai<-out3[i,1]
bi<-out3[i,2]
Ri<-ai*t*exp(-bi*t)
HPDI[t,1]<-HPDinterval(as.mcmc(Ri),prop=0.95)[1]
HPDI[t,2]<-HPDinterval(as.mcmc(Ri),prop=0.95)[2]
}
lines(1:max(SR$Spawners[1:31]),HPDI[,1],type='l',col='cyan')
lines(1:max(SR$Spawners[1:31]),HPDI[,2],type='l',col='cyan')