Coho Parr per Female Spawner
Matt Falcy
November 8, 2017
The issue
Spawner-recruit data are notoriously “noisy.” But not if the recruits are Oregon coastal coho parr. Here we look at adult female-to-parr data. It has been customary to present data amalgamated at the scale of most of the oregon coast. But we can break the data down into a finer spatial scale that more closely reflects how density-dependence is operating. Here are the data at the scale of the “North Coast”, “Mid Coast”, “Umpqua” and “Mid South” strata. An even finer scale would be desireable but alas not possible.
The data
setwd("//Kalawatseti/home/falcym/docs/Projects/People/Constable")
D<-read.csv('FemalesParrStrata.csv')
D
## Brood.Year Females Parr PopID Pop.name
## 1 1998 1170.500 175942.6 1 NC
## 2 1999 4047.000 1475425.1 1 NC
## 3 2000 9109.000 1870350.1 1 NC
## 4 2001 16343.000 2092193.5 1 NC
## 5 2002 24922.729 2805006.0 1 NC
## 6 2003 27614.892 2420595.1 1 NC
## 7 2004 14835.255 2451858.8 1 NC
## 8 2005 6488.461 4041800.3 1 NC
## 9 2006 10819.138 2925199.0 1 NC
## 10 2007 9181.857 1065508.9 1 NC
## 11 2008 12602.850 2384639.6 1 NC
## 12 2009 21563.662 2227115.3 1 NC
## 13 2010 24787.210 2134810.2 1 NC
## 14 2011 23231.812 1658095.2 1 NC
## 15 2012 4075.000 1319598.7 1 NC
## 16 2013 5478.000 968781.0 1 NC
## 17 1998 1220.500 579785.3 2 MC
## 18 1999 4809.000 1829197.0 2 MC
## 19 2000 7781.000 2307964.0 2 MC
## 20 2001 10537.500 2062590.6 2 MC
## 21 2002 45372.910 2509646.5 2 MC
## 22 2003 37792.079 1932978.8 2 MC
## 23 2004 21432.770 1753236.8 2 MC
## 24 2005 27650.735 3413790.1 2 MC
## 25 2006 9727.667 2464333.9 2 MC
## 26 2007 6358.905 2313408.6 2 MC
## 27 2008 36587.505 3866857.8 2 MC
## 28 2009 49153.970 2397812.2 2 MC
## 29 2010 31231.665 2305825.2 2 MC
## 30 2011 64226.162 2901726.5 2 MC
## 31 2012 17972.000 3211346.2 2 MC
## 32 2013 22046.000 2948208.5 2 MC
## 33 1998 4875.500 758353.9 3 Ump
## 34 1999 4288.000 2722720.2 3 Ump
## 35 2000 7297.000 1483618.6 3 Ump
## 36 2001 17542.000 2449398.9 3 Ump
## 37 2002 21588.099 2406163.4 3 Ump
## 38 2003 17537.647 2654087.1 3 Ump
## 39 2004 13983.590 1460830.6 3 Ump
## 40 2005 21776.384 1479648.2 3 Ump
## 41 2006 7592.179 2758898.7 3 Ump
## 42 2007 5343.453 1704212.9 3 Ump
## 43 2008 20272.768 3089470.9 3 Ump
## 44 2009 27531.971 2161491.5 3 Ump
## 45 2010 36917.496 3205401.3 3 Ump
## 46 2011 51107.652 2057585.9 3 Ump
## 47 2012 10977.000 1915522.9 3 Ump
## 48 2013 16858.000 1775416.6 3 Ump
## 49 1998 8971.000 1028818.7 4 MS
## 50 1999 10377.000 2194129.8 4 MS
## 51 2000 12658.500 2869687.1 4 MS
## 52 2001 36395.500 3038375.3 4 MS
## 53 2002 33123.685 2718526.1 4 MS
## 54 2003 34119.665 2530359.0 4 MS
## 55 2004 41421.516 3203429.6 4 MS
## 56 2005 22490.968 3379361.5 4 MS
## 57 2006 34439.407 3693253.4 4 MS
## 58 2007 14081.874 3819113.7 4 MS
## 59 2008 28098.812 4859646.7 4 MS
## 60 2009 38759.052 3280938.6 4 MS
## 61 2010 61651.656 4980190.2 4 MS
## 62 2011 58754.739 4583892.2 4 MS
## 63 2012 21067.000 6301492.8 4 MS
## 64 2013 25272.000 2767419.5 4 MS
I’ll maximize the likelihood of the parameters (\(\alpha\), \(\beta\), \(\epsilon\)) of a Beverton-Holt model given coho female spawner and parr recruit data:
Model 1: A universal curve
\[ \LARGE R_t=\frac{\alpha S_t}{1+\frac{\alpha}{\beta} S_t}e^{\epsilon_t}\]
#Ignore strata
LogLik<-function(D,params){
LL<-vector()
for (i in 1:dim(D)[1]){
LL[i]<-dnorm(log(D[i,3]),log(params[1])+log(D[i,2])-log(1+params[1]/params[2]*D[i,2]),params[3],log=TRUE)
}
-1*sum(LL)
}
out<-optim(par=c(500,3000000,10),fn=LogLik,D=D,control=list(reltol=1e-16,maxit=5000))
#out<-optim(par=c(450,3500000,0.5),fn=LogLik,D=D,method="L-BFGS-B",lower=c(20,1000,0.01),upper=c(1000,10000000,500),control=list(factr=1e-16))
plot(D[,2],D[,3],xlab='Females',ylab='Parr')
xi<-seq(1,max(D[,2]),by=1000)
yi<-(out$par[1]*xi)/(1+out$par[1]/out$par[2]*xi)
lines(xi,yi)
That’s nice. But perhaps we should allow each management area to have a uniqe asymptote, \(\beta\)
Model 2: Universal \(\alpha\), strata-specific \(\beta\)
Expand the Beverton-Holt model to accomodate sptatial sturcture of the data. For strata \(i\),
\[ \LARGE R_{i,t}=\frac{\alpha S_{i,t}}{1+\frac{\alpha}{\beta_i}S_{i,t}}e^{\epsilon_{i,t}}\]
LogLik2<-function(D,params){
LL2<-vector()
for (i in 1:dim(D)[1]){
LL2[i]<-dnorm(log(D[i,3]),log(params[1])+log(D[i,2])-log(1+params[1]/params[1+D[i,4]]*D[i,2]),params[6],log=TRUE)
}
-1*sum(LL2)
}
#out2<-optim(par=c(450,3500000,3500000,3500000,4000000,0.5),fn=LogLik2,D=D,method="L-BFGS-B",lower=c(20,1000,1000,1000,1000,0.01),upper=c(1000,10000000,10000000,10000000,10000000,500),control=list(factr=1e-16))
out2<-optim(par=c(450,3500000,3500000,3500000,4000000,0.3),fn=LogLik2,D=D,control=list(factr=1e-16,maxit=5000),hessian=TRUE)
plot(D[1:16,2],D[1:16,3],xlab='Females',ylab='Parr',col=2,xlim=c(0,max(D[,2])),ylim=c(0,max(D[,3])),pch=16)
lines(D[17:32,2],D[17:32,3],col=3,type='p',pch=16)
lines(D[33:48,2],D[33:48,3],col=4,type='p',pch=16)
lines(D[49:64,2],D[49:64,3],col=5,type='p',pch=16)
xi<-seq(1,max(D[,2]),by=1000)
yi1<-(out2$par[1]*xi)/(1+out2$par[1]/out2$par[2]*xi)
yi2<-(out2$par[1]*xi)/(1+out2$par[1]/out2$par[3]*xi)
yi3<-(out2$par[1]*xi)/(1+out2$par[1]/out2$par[4]*xi)
yi4<-(out2$par[1]*xi)/(1+out2$par[1]/out2$par[5]*xi)
lines(xi,yi1,col=2)
lines(xi,yi2,col=3)
lines(xi,yi3,col=4)
lines(xi,yi4,col=5)
legend(x=40000,y=6e+06,c("NC","MC","Ump","MS"),col=c(2,3,4,5),pch=c(16,16,16,16))
Compare AICs of the two models
#Model1
2*out$value + 2*length(out$par)
## [1] 56.85294
#Model2
2*out2$value + 2*length(out2$par)
## [1] 57.21791
Relatively lower AIC values indicate better models. However, these two AIC values are just 0.365 units part, suggesting that there is very weak evidence that Model 1 is indeed better than Model 2.
Plot residuals and look for trend
Now let’s plot residuals from Model 2 and look for a trend.
resid<-vector()
for (i in 1:16){
resid[i]<-D[i,3]-(out2$par[1]*D[i,2])/(1+out2$par[1]/out2$par[2]*D[i,2])
}
for (i in 17:32){
resid[i]<-D[i,3]-(out2$par[1]*D[i,2])/(1+out2$par[1]/out2$par[3]*D[i,2])
}
for (i in 33:48){
resid[i]<-D[i,3]-(out2$par[1]*D[i,2])/(1+out2$par[1]/out2$par[4]*D[i,2])
}
for (i in 49:64){
resid[i]<-D[i,3]-(out2$par[1]*D[i,2])/(1+out2$par[1]/out2$par[5]*D[i,2])
}
plot(seq(1998,2013),resid[1:16],xlab='Year',ylab='Residuals',type='b',col=2,ylim=c(-1.5e6,3e6))
lines(seq(1998,2013),resid[17:32],type='b',col=3)
lines(seq(1998,2013),resid[33:48],type='b',col=4)
lines(seq(1998,2013),resid[49:64],type='b',col=5)
legend(x=2000,y=2500000,c("NC","MC","Ump","MS"),col=c(2,3,4,5),pch=c(16,16,16,16))
abline(0,0)
Not much, but the Mid-South might be on an upward trend. Let’s assess signifigance. The best way to do this would be to include time in the S/R model rather than doing two separate analysis. Let’s just look a simple linear regression of residuals on time for now:
summary(lm(resid[49:64]~seq(1998, 2013)))
##
## Call:
## lm(formula = resid[49:64] ~ seq(1998, 2013))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1764534 -449762 -100121 539790 2067685
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -279241496 101080348 -2.763 0.0153 *
## seq(1998, 2013) 139318 50401 2.764 0.0152 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 929400 on 14 degrees of freedom
## Multiple R-squared: 0.3531, Adjusted R-squared: 0.3069
## F-statistic: 7.641 on 1 and 14 DF, p-value: 0.01522
Yes, the MS residuals have statistically signifigant positive trend (P=0.015). Note that this means that he assumption of stationarity in the S/R relationship is not valid. Further AIC analyses could be done to determine if the best model of the data should include time as a covariate.