In this lab, we applied Dynamic Linear Models to examine some of the time-varying properties of the spawner-recruit relationship for Pacific salmon. We downloaded sockeye data from the Kvichak river system in Bristol Bay Alaska using the MARSS package. This included spawner and recruit data spanning from 1956 to 2005.
Environmental covariates
We tested the data support for the inclusion of annual summer and winter PDO values in our models. We also tested the support for sea surface temperature (SST). These data were also downloaded from the MARSS package in R.
To examine our model assumptions, we inspected acf plots of our model residuals.
Data
The data are a dataframe with columns for brood year (brood_year), number of spawners (spawners), number of recruits (recruits) and PDO at year \(t-2\) in summer (pdo_summer_t2) and in winter (pdo_winter_t2).
Load the data
Code
## library(devtools)## Windows users will likely need to set this#Sys.setenv("R_REMOTES_NO_ERRORS_FROM_WARNINGS" = "true")#devtools::install_github("nwfsc-timeseries/atsalibrary")library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.3 ✔ readr 2.1.4
✔ forcats 1.0.0 ✔ stringr 1.5.0
✔ ggplot2 3.4.2 ✔ tibble 3.2.1
✔ lubridate 1.9.2 ✔ tidyr 1.3.0
✔ purrr 1.0.1
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ purrr::%||%() masks base::%||%()
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Code
library(MARSS)## get datadata(KvichakSockeye, package="atsalibrary")SR_data <- KvichakSockeyehead(SR_data)
#create log ratio for ytSR_data$y <-log(SR_data$recruits/SR_data$spawners)#date to columnsdat <-t(SR_data)colnames(dat) <- SR_data$brood_yeardat <- dat[-1,]dat <- dat[c(5,1,3,4,2),]## plot response variabledat.ts <-ts(log(SR_data$recruits/SR_data$spawners), start =1952, end =2005, frequency =1)plot.ts(dat.ts)
Results
Task 1: Reduced Ricker
We first fit a reduced form of Ricker’s model that includes only a time-varying level (\(\alpha_t\)) and observation error (\(v_t\)). That is,
This model assumes no density-dependent survival in that the number of recruits is an ascending function of spawners.
Code
# NO density dependence / stochastic level mod1 <-list(Z ="identity",U="zero",R=matrix("r",nrow=1),B="identity",A="zero",Q=matrix("q",nrow=1))m1 <-MARSS(dat[1,],model = mod1)
Success! abstol and log-log tests passed at 29 iterations.
Alert: conv.test.slope.tol is 0.5.
Test with smaller values (<0.1) to ensure convergence.
MARSS fit is
Estimation method: kem
Convergence test: conv.test.slope.tol = 0.5, abstol = 0.001
Estimation converged in 29 iterations.
Log-likelihood: -52.53672
AIC: 111.0734 AICc: 111.6888
Estimate
R.r 0.228
Q.q 0.281
x0.x0 0.952
Initial states (x0) defined at t=0
Standard errors have not been calculated.
Use MARSSparamCIs to compute CIs and bias estimates.
This simplified model show some signs of autocorrelation in residuals.
Task 2: Full Ricker Model
We then fit the full Ricker model with a time-varying intercept and a static effect of spawners.
Code
#we model the underlying states of alpha (brood-year productivity) and beta#beta is prevented from changing.spawn <-matrix(dat[2,-c(1:4)],nrow=1) #not sure if we should scale this or notratio <- SR_data$y[-c(1:4)] #remove responses that corresponded to Spawner = NAm <-2#number of parameters in process TT <-length(spawn) #number of data pointsB <-diag(m) #"identity"U <-matrix(0, nrow = m, ncol =1) #"zero"Q <-matrix(list("q.alpha", 0, 0, 0), nrow =2) # Setting the effect of spawners to staticZ <-array(NA, c(1, m, TT)) ## NxMxT; empty for nowZ[1,1,] <-rep(1, TT) ## Nx1; 1's for interceptZ[1,2,] <- spawn ## Nx1; predictor variable A =matrix(0) #"zero"R <-matrix("r")inits_list <-list(x0 =matrix(c(0, 0), nrow = m)) mod_list <-list(B = B, U = U, Q = Q, Z = Z, A = A, R = R)## fit the model with MARSSm2 <-MARSS(ratio, mod_list, inits = inits_list)
Success! abstol and log-log tests passed at 29 iterations.
Alert: conv.test.slope.tol is 0.5.
Test with smaller values (<0.1) to ensure convergence.
MARSS fit is
Estimation method: kem
Convergence test: conv.test.slope.tol = 0.5, abstol = 0.001
Estimation converged in 29 iterations.
Log-likelihood: -51.81022
AIC: 111.6204 AICc: 112.6731
Estimate
R.r 2.19e-01
Q.q.alpha 2.92e-01
x0.X1 9.35e-01
x0.X2 7.02e-06
Initial states (x0) defined at t=0
Standard errors have not been calculated.
Use MARSSparamCIs to compute CIs and bias estimates.
\(\alpha_t\) plot and the AICc.
The time series of \(\alpha_t\) for this model is an estimate of the stock productivity in the absence of density-dependent effects. Here we plot the ts of \(\alpha_t\) and note the AICc for this model.
Code
alpha <-as.numeric(m2$states[1,])beta <-as.numeric(m2$states[2,]) # Here, beta does not time vary (this is the static effect of spawners)alpha.se <-as.numeric(m2$states.se[1,])beta.se <-as.numeric(m2$states.se[2,])#from first modelplot(alpha1~SR_data$brood_year,type='l', ylim=c(-3,3), xlim=c(1950, 2005))lines(alpha1+2*alpha.se1~SR_data$brood_year, lty="dashed")lines(alpha1-2*alpha.se1~SR_data$brood_year, lty="dashed")#second modellines(alpha~SR_data$brood_year[-c(1:4)],type='l',col="red")lines(alpha+2*alpha.se~SR_data$brood_year[-c(1:4)], lty="dashed", col="red")lines(alpha-2*alpha.se~SR_data$brood_year[-c(1:4)], lty="dashed", col="red")
Code
#AIC valuem1$AICc
[1] 111.6888
Code
m2$AICc #AICc is within a single point.
[1] 112.6731
How do these estimates of productivity compare to those from the previous question?
We were expecting the alphas (productivity) from model 2 to be smaller than model 1, however the estimated values of alpha are very similar in each model.
The model AICcs are also very similar and only vary by a single point.
Model diagnostics
Code
#plot model diagnostics. #autoplot.marssMLE(m2)plot(m2, plot.type ="qqplot.std.model.resids.ytt1")
plot type = qqplot.std.model.resids.ytt1
Code
plot(m2, plot.type ="acf.std.model.resids.ytt1")
plot type = acf.std.model.resids.ytt1
Similar to our first model, the ACF plot still shows a lag at 5 and 10.
Task 3: Expanded Ricker with summer PDO
Next we fit the expanded Ricker model including the summer PDO index as the covariate (pdo_summer_t2).
Code
#first model we assume there is no density dependence#we model the underlying states of alpha (brood-year productivity) and beta#beta is prevented from changing.pdo_s <-matrix(dat[3,-c(1:4)], nrow=1)m <-3#number of parameters in process TT <-length(spawn) #number of data pointsB <-diag(m) #"identity"U <-matrix(0, nrow = m, ncol =1) #"zero"Q <-matrix(list("q.alpha", 0,0, #we no longer want alpha to vary0, 0,0, #Do we want beta to vary now?0,0,"q.pdo_s"), nrow =3, byrow =TRUE) # to have characters and numbers in same matrix use a listZ <-array(NA, c(1, m, TT)) ## NxMxT; empty for nowZ[1,1,] <-rep(1, TT) ## Nx1; 1's for interceptZ[1,2,] <- spawn ## Nx1; predictor variable spawnersZ[1,3,] <- pdo_s ## Nx1; predictor PDO summerA =matrix(0) #"zero"R <-matrix("r")inits_list2 <-list(x0 =matrix(c(0,0,0), nrow = m)) mod_list2 <-list(B = B, U = U, Q = Q, Z = Z, A = A, R = R)cont_list <-list(maxit=10000)# run modelm3 <-MARSS(ratio, mod_list2, inits = inits_list2, method ="BFGS")
Success! Converged in 21 iterations.
Function MARSSkfas used for likelihood calculation.
MARSS fit is
Estimation method: BFGS
Estimation converged in 21 iterations.
Log-likelihood: -51.78324
AIC: 115.5665 AICc: 117.8998
Estimate
R.r 2.00e-01
Q.q.alpha 3.17e-01
Q.q.pdo_s 1.73e-14
x0.X1 1.01e+00
x0.X2 5.49e-06
x0.X3 2.85e-02
Initial states (x0) defined at t=0
Standard errors have not been calculated.
Use MARSSparamCIs to compute CIs and bias estimates.
Plot the ts of \(\delta_t\) and note the AICc for this model.
The plot for state = X3 is the time series of \(delta_t\). The model estimates X3 to be 0.03 and static.
Code
autoplot(m3, plot.type ="xtT")
Model diagnostics
Code
plot(m3)
plot type = fitted.ytT Observations with fitted values
Hit <Return> to see next plot (q to exit):
plot type = xtT Estimated states
Hit <Return> to see next plot (q to exit):
plot type = model.resids.ytt1
Hit <Return> to see next plot (q to exit):
plot type = std.model.resids.ytT
Hit <Return> to see next plot (q to exit):
plot type = std.state.resids.xtT
Hit <Return> to see next plot (q to exit):
plot type = qqplot.std.model.resids.ytt1
Hit <Return> to see next plot (q to exit):
plot type = qqplot.std.state.resids.xtT
Hit <Return> to see next plot (q to exit):
plot type = acf.std.model.resids.ytt1
Task 4: Expanded Ricker with winter PDO
Finally, we fit the expanded model including the winter PDO index as the covariate (pdo_winter_t2).
Code
pdo_w <-matrix(dat[4,-c(1:4)], nrow=1)m <-3#number of parameters in process TT <-length(spawn) #number of data pointsB <-diag(m) #"identity"U <-matrix(0, nrow = m, ncol =1) #"zero"Q <-matrix(list("q.alpha", 0,0, #we no longer want alpha to vary0, 0,0, #Do we want beta to vary now?0,0,"q.pdo_w"), nrow =3, byrow =TRUE) # to have characters and numbers in same matrix use a listZ <-array(NA, c(1, m, TT)) ## NxMxT; empty for nowZ[1,1,] <-rep(1, TT) ## Nx1; 1's for interceptZ[1,2,] <- spawn ## Nx1; predictor variable spawnersZ[1,3,] <- pdo_w ## Nx1; predictor PDO summer#Z[1,4,] <- pdo_w ## Nx1; predictor PDO winterA =matrix(0) #"zero"R <-matrix("r")inits_list4 <-list(x0 =matrix(c(0,0,0), nrow = m)) mod_list4 <-list(B = B, U = U, Q = Q, Z = Z, A = A, R = R)cont_list <-list(maxit=10000)# run modelm4 <-MARSS(ratio, mod_list4, inits = inits_list4, control = cont_list, method ="BFGS")
Success! Converged in 18 iterations.
Function MARSSkfas used for likelihood calculation.
MARSS fit is
Estimation method: BFGS
Estimation converged in 18 iterations.
Log-likelihood: -51.38573
AIC: 114.7715 AICc: 117.1048
Estimate
R.r 1.96e-01
Q.q.alpha 3.11e-01
Q.q.pdo_w 4.17e-12
x0.X1 1.13e+00
x0.X2 1.32e-06
x0.X3 1.54e-01
Initial states (x0) defined at t=0
Standard errors have not been calculated.
Use MARSSparamCIs to compute CIs and bias estimates.
\(\alpha_t\), \(\delta_t\) plots and the AICc.
What is the mean level of productivity? Plot the ts of \(\delta_t\) and note the AICc for this model. The mean productivity level is 0.28.
plot type = fitted.ytT Observations with fitted values
Hit <Return> to see next plot (q to exit):
plot type = xtT Estimated states
Hit <Return> to see next plot (q to exit):
plot type = model.resids.ytt1
Hit <Return> to see next plot (q to exit):
plot type = std.model.resids.ytT
Hit <Return> to see next plot (q to exit):
plot type = std.state.resids.xtT
Hit <Return> to see next plot (q to exit):
plot type = qqplot.std.model.resids.ytt1
Hit <Return> to see next plot (q to exit):
plot type = qqplot.std.state.resids.xtT
Hit <Return> to see next plot (q to exit):
plot type = acf.std.model.resids.ytt1
Task 5: Best Model Forecasts
Model comparison
Code
#AIC valuem1$AICc
[1] 111.6888
Code
m2$AICc #AICc is within a single point.
[1] 112.6731
Code
m3$AICc # getting worse
[1] 117.8998
Code
m4$AICc
[1] 117.1048
Based on AICc, which of the models above is the most parsimonious? Is it well behaved (i.e., are the model assumptions met)? Plot the model forecasts for the best model. Is this a good forecast model?
The best model is our first model with just a time varying intercept. However, the model with a static effect of spawners had almost equal support. The model is well behaved and all model assumptions are met.
Model forecasts
Code
fr1 <-predict(m1, type="ytT", n.ahead=1)plot(fr1)
Task 6: Other Environmental Covariates
As our final task, we considered the effect of sea surface tmeperature on the time-varying productivity of Kvichak salmon.
Code
bb_data <-readRDS(here::here("Lab-1", "Data_Images", "bristol_bay_data_plus_covariates.rds"))sst <- bb_data %>%filter(system =="Kvichak", brood_yr <=2005)%>%select(brood_yr, env_sst)%>%group_by(brood_yr)%>%summarize(sst =mean(env_sst))spawn1957 <- spawn[-1]ratio1957 <- ratio[-1]sst1957 <-t(sst$sst)m <-3#number of parameters in process TT <-length(spawn1957) #number of data pointsB <-diag(m) #"identity"U <-matrix(0, nrow = m, ncol =1) #"zero"Q <-matrix(list("q.alpha", 0,0, 0, 0,0, #Do we want beta to vary now?0,0,"q.sst"), nrow =3, byrow =TRUE) # to have characters and numbers in same matrix use a listZ <-array(NA, c(1, m, TT)) ## NxMxT; empty for nowZ[1,1,] <-rep(1, TT) ## Nx1; 1's for interceptZ[1,2,] <- spawn1957 ## Nx1; predictor variable spawnersZ[1,3,] <- sst1957 ## Nx1; predictor PDO summerA =matrix(0) #"zero"R <-matrix("r")inits_list2 <-list(x0 =matrix(c(0,0,0), nrow = m)) mod_list2 <-list(B = B, U = U, Q = Q, Z = Z, A = A, R = R)# run modelm5 <-MARSS(ratio1957, mod_list2, inits = inits_list2)
Success! abstol and log-log tests passed at 144 iterations.
Alert: conv.test.slope.tol is 0.5.
Test with smaller values (<0.1) to ensure convergence.
MARSS fit is
Estimation method: kem
Convergence test: conv.test.slope.tol = 0.5, abstol = 0.001
Estimation converged in 144 iterations.
Log-likelihood: -48.97244
AIC: 109.9449 AICc: 112.3449
Estimate
R.r 1.60e-01
Q.q.alpha 3.73e-02
Q.q.sst 7.65e-03
x0.X1 2.48e+00
x0.X2 -5.17e-06
x0.X3 -3.75e-01
Initial states (x0) defined at t=0
Standard errors have not been calculated.
Use MARSSparamCIs to compute CIs and bias estimates.
Code
m5$AICc
[1] 112.3449
\(\alpha_t\), \(\delta_t\) plots and the AICc.
What is the mean level of productivity? Plot the ts of \(\delta_t\) and note the AICc for this model.
The mean level of productivity in this model was 2.48. We found that this model did not test higher than the reduced Ricker model.
Model Diagnostics
Also plot appropriate model diagnostics.
Code
plot(m5)
MARSSresiduals.tt1 reported warnings. See msg element of returned residuals object.
plot type = fitted.ytT Observations with fitted values
Hit <Return> to see next plot (q to exit):
plot type = xtT Estimated states
Hit <Return> to see next plot (q to exit):
plot type = model.resids.ytt1
Hit <Return> to see next plot (q to exit):
plot type = std.model.resids.ytT
Hit <Return> to see next plot (q to exit):
plot type = std.state.resids.xtT
Hit <Return> to see next plot (q to exit):
plot type = qqplot.std.model.resids.ytt1
Hit <Return> to see next plot (q to exit):
plot type = qqplot.std.state.resids.xtT
Hit <Return> to see next plot (q to exit):
plot type = acf.std.model.resids.ytt1
Discussion
We found that the Kvichak sockeye spawner recruit time series was best explained by the simplest model with just a time varying intercept. Of the three covariates (summer pdo, winter pdo and SST) we tested, the model with SST had the lowest AIC but not lower than the simplest model. We found that the effect of summer and winter pdo were not time varying. This might be because a time varying intercept explained most of the variance.
Team contributions
Dylan wrote the original code for reduced ricker and full ricker models. Maria wrote the code for the PDO models. Liz wrote the code for the SST and coded the models in matrix form. All team members contributed to writing and organization of the lab document.
Source Code
---title: "Team 1 - Lab 5"subtitle: "Dynamic Linear Models (DLMs)"author: "Dylan, Liz, Maria"date: May 11, 2023output: html_document: code-folding: true toc: true toc_float: true---```{r setup, include = FALSE}options(dplyr.summarise.inform = FALSE)```***# Methods## DataIn this lab, we applied Dynamic Linear Models to examine some of the time-varying properties of the spawner-recruit relationship for Pacific salmon. We downloaded sockeye data from the Kvichak river system in Bristol Bay Alaska using the MARSS package. This included spawner and recruit data spanning from 1956 to 2005.## Environmental covariatesWe tested the data support for the inclusion of annual summer and winter PDO values in our models. We also tested the support for sea surface temperature (SST). These data were also downloaded from the MARSS package in R. ## Models in matrix formWe fit the following four models:Model 1: Reduced Ricker model$$\begin{equation}y_t = \begin{bmatrix}1\end{bmatrix}\begin{bmatrix} \alpha \\ \end{bmatrix}_t + v_t \\\begin{bmatrix} \alpha \\ \end{bmatrix}_t \ = \begin{bmatrix} \alpha \\ \end{bmatrix}_{t-1} + \begin{bmatrix} w_{\alpha} \\ \end{bmatrix}_t\end{equation}$$$$w_t \sim N(0, \begin{bmatrix} q_\alpha \end{bmatrix})$$Model 2: Full Ricker model$$\begin{equation}y_t = \begin{bmatrix} 1 & S_t \end{bmatrix}\begin{bmatrix} \alpha \\ \beta \end{bmatrix}_t + v_t \\ \begin{bmatrix} \alpha \\ \beta \end{bmatrix}_t = \begin{bmatrix} \alpha \\ \beta \end{bmatrix}_{t-1} + \begin{bmatrix} w_{\alpha} \\ w_{\beta} \end{bmatrix}_t \end{equation}$$$$w_t \sim N(0, \begin{bmatrix} q_\alpha & 0 \\ 0 & 0 \end{bmatrix})$$Model 3: Full Ricker model w/ summer pdo as a covariate$$\begin{equation}y_t = \begin{bmatrix} 1 & S_t & PDOsummer_t \end{bmatrix}\begin{bmatrix} \alpha \\ \beta \\ \delta\\ \end{bmatrix}_t + v_t \\ \begin{bmatrix} \alpha \\ \beta \\ \delta \\ \end{bmatrix}_t = \begin{bmatrix} \alpha \\ \beta \\ \delta \\ \end{bmatrix}_{t-1} + \begin{bmatrix} w_{\alpha} \\ w_{\beta} \\ w_{\delta} \\ \end{bmatrix}_t \end{equation}$$$$w_t \sim N(0, \begin{bmatrix} q_\alpha & 0 & 0\\ 0 & 0 & 0 \\ 0 & 0 & q_\delta \end{bmatrix})$$Model 4: Full Ricker model w/ winter pdo as a covariate$$\begin{equation}y_t = \begin{bmatrix} 1 & S_t & PDOwinter_t \end{bmatrix}\begin{bmatrix} \alpha \\ \beta \\ \delta\\ \end{bmatrix}_t + v_t \\ \begin{bmatrix} \alpha \\ \beta \\ \delta \\ \end{bmatrix}_t = \begin{bmatrix} \alpha \\ \beta \\ \delta \\ \end{bmatrix}_{t-1} + \begin{bmatrix} w_{\alpha} \\ w_{\beta} \\ w_{\delta} \\ \end{bmatrix}_t \end{equation}$$$$w_t \sim N(0, \begin{bmatrix} q_\alpha & 0 & 0\\ 0 & 0 & 0 \\ 0 & 0 & q_\delta \end{bmatrix})$$## Model diagnositicsTo examine our model assumptions, we inspected acf plots of our model residuals. # DataThe data are a dataframe with columns for brood year (`brood_year`), number of spawners (`spawners`), number of recruits (`recruits`) and PDO at year $t-2$ in summer (`pdo_summer_t2`) and in winter (`pdo_winter_t2`).## Load the data```{r dlm-load-atsa}## library(devtools)## Windows users will likely need to set this#Sys.setenv("R_REMOTES_NO_ERRORS_FROM_WARNINGS" = "true")#devtools::install_github("nwfsc-timeseries/atsalibrary")library(tidyverse)library(MARSS)## get datadata(KvichakSockeye, package="atsalibrary")SR_data <- KvichakSockeyehead(SR_data)```## Wrangle the data```{r wrangle_data}#create log ratio for ytSR_data$y <- log(SR_data$recruits/SR_data$spawners)#date to columnsdat <- t(SR_data)colnames(dat) <- SR_data$brood_yeardat <- dat[-1,]dat <- dat[c(5,1,3,4,2),]## plot response variabledat.ts <- ts(log(SR_data$recruits/SR_data$spawners), start = 1952, end = 2005, frequency =1)plot.ts(dat.ts)```# Results## Task 1: Reduced Ricker1. We first fit a reduced form of Ricker's model that includes only a time-varying level ($\alpha_t$) and observation error ($v_t$). That is,$$\begin{align*}\text{log}(R_t) &= \alpha_t + \text{log}(S_t) + v_t \\\text{log}(R_t/S_t) &= \alpha_t + v_t\end{align*}$$This model assumes no density-dependent survival in that the number of recruits is an ascending function of spawners. ```{r reduced Ricker}# NO density dependence / stochastic level mod1 <- list( Z = "identity", U="zero", R=matrix("r",nrow=1), B="identity", A="zero", Q=matrix("q",nrow=1))m1 <- MARSS(dat[1,],model = mod1)```### $\alpha_t$ plot and the AICc. ```{r alpha values}#plot alpha valuesalpha1 <- as.numeric(m1$states)alpha.se1 <- as.numeric(m1$states.se)plot(alpha1~SR_data$brood_year,type='l', ylim=c(-3,3))lines(alpha1+2*alpha.se1~SR_data$brood_year, lty="dashed")lines(alpha1-2*alpha.se1~SR_data$brood_year, lty="dashed")#AICc valuem1$AICc```The AICc of our first model is 111.688. ### Model diagnostics```{r}#plot model residuals diagnostics. res <-residuals(m1)[,7]tmp <-which(!is.na(res))res <- res[tmp]acf(res)#or #autoplot.marssMLE(m1)plot(m1, plot.type ="qqplot.std.model.resids.ytt1")plot(m1, plot.type ="acf.std.model.resids.ytt1")```This simplified model show some signs of autocorrelation in residuals.## Task 2: Full Ricker ModelWe then fit the full Ricker model with a time-varying intercept and a static effect of spawners. ```{r model_1, cache = TRUE}#we model the underlying states of alpha (brood-year productivity) and beta#beta is prevented from changing.spawn <- matrix(dat[2,-c(1:4)],nrow=1) #not sure if we should scale this or notratio <- SR_data$y[-c(1:4)] #remove responses that corresponded to Spawner = NAm <- 2 #number of parameters in process TT <- length(spawn) #number of data pointsB <- diag(m) #"identity"U <- matrix(0, nrow = m, ncol = 1) #"zero"Q <- matrix(list("q.alpha", 0, 0, 0), nrow = 2) # Setting the effect of spawners to staticZ <- array(NA, c(1, m, TT)) ## NxMxT; empty for nowZ[1,1,] <- rep(1, TT) ## Nx1; 1's for interceptZ[1,2,] <- spawn ## Nx1; predictor variable A = matrix(0) #"zero"R <- matrix("r")inits_list <- list(x0 = matrix(c(0, 0), nrow = m)) mod_list <- list(B = B, U = U, Q = Q, Z = Z, A = A, R = R)## fit the model with MARSSm2 <- MARSS(ratio, mod_list, inits = inits_list)```### $\alpha_t$ plot and the AICc. The time series of $\alpha_t$ for this model is an estimate of the stock productivity in the absence of density-dependent effects. Here we plot the ts of $\alpha_t$ and note the AICc for this model. ```{r}alpha <-as.numeric(m2$states[1,])beta <-as.numeric(m2$states[2,]) # Here, beta does not time vary (this is the static effect of spawners)alpha.se <-as.numeric(m2$states.se[1,])beta.se <-as.numeric(m2$states.se[2,])#from first modelplot(alpha1~SR_data$brood_year,type='l', ylim=c(-3,3), xlim=c(1950, 2005))lines(alpha1+2*alpha.se1~SR_data$brood_year, lty="dashed")lines(alpha1-2*alpha.se1~SR_data$brood_year, lty="dashed")#second modellines(alpha~SR_data$brood_year[-c(1:4)],type='l',col="red")lines(alpha+2*alpha.se~SR_data$brood_year[-c(1:4)], lty="dashed", col="red")lines(alpha-2*alpha.se~SR_data$brood_year[-c(1:4)], lty="dashed", col="red")#AIC valuem1$AICcm2$AICc #AICc is within a single point.```How do these estimates of productivity compare to those from the previous question? We were expecting the alphas (productivity) from model 2 to be smaller than model 1, however the estimated values of alpha are very similar in each model. The model AICcs are also very similar and only vary by a single point.### Model diagnostics ```{r}#plot model diagnostics. #autoplot.marssMLE(m2)plot(m2, plot.type ="qqplot.std.model.resids.ytt1")plot(m2, plot.type ="acf.std.model.resids.ytt1")```Similar to our first model, the ACF plot still shows a lag at 5 and 10.## Task 3: Expanded Ricker with summer PDONext we fit the expanded Ricker model including the summer PDO index as the covariate (`pdo_summer_t2`). ```{r model_3, cache = TRUE}#first model we assume there is no density dependence#we model the underlying states of alpha (brood-year productivity) and beta#beta is prevented from changing.pdo_s <- matrix(dat[3,-c(1:4)], nrow=1)m <- 3 #number of parameters in process TT <- length(spawn) #number of data pointsB <- diag(m) #"identity"U <- matrix(0, nrow = m, ncol = 1) #"zero"Q <- matrix(list("q.alpha", 0,0, #we no longer want alpha to vary 0, 0,0, #Do we want beta to vary now? 0,0,"q.pdo_s"), nrow = 3, byrow = TRUE) # to have characters and numbers in same matrix use a listZ <- array(NA, c(1, m, TT)) ## NxMxT; empty for nowZ[1,1,] <- rep(1, TT) ## Nx1; 1's for interceptZ[1,2,] <- spawn ## Nx1; predictor variable spawnersZ[1,3,] <- pdo_s ## Nx1; predictor PDO summerA = matrix(0) #"zero"R <- matrix("r")inits_list2 <- list(x0 = matrix(c(0,0,0), nrow = m)) mod_list2 <- list(B = B, U = U, Q = Q, Z = Z, A = A, R = R)cont_list <- list(maxit=10000)# run modelm3 <- MARSS(ratio, mod_list2, inits = inits_list2, method = "BFGS")```### $\alpha_t$, $\delta_t$ plots and the AICc.What is the mean level of productivity? The mean level of productivity is 0.26.```{r}alpha3 <-as.numeric(m3$states[1,]) #static numberbeta3 <-as.numeric(m3$states[2,])pdo_s <-as.numeric(m3$states[3,])alpha.se3 <-as.numeric(m3$states.se[1,])beta.se3 <-as.numeric(m3$states.se[2,])pdo_s.se3 <-as.numeric(m3$states.se[3,])mean_produc_m3 <-mean(alpha3)#from first modelplot(alpha1~SR_data$brood_year,type='l', ylim=c(-3,3), xlim=c(1950, 2005))lines(alpha1+2*alpha.se1~SR_data$brood_year, lty="dashed")lines(alpha1-2*alpha.se1~SR_data$brood_year, lty="dashed")#second modellines(alpha~SR_data$brood_year[-c(1:4)],type='l',col="red")lines(alpha+2*alpha.se~SR_data$brood_year[-c(1:4)], lty="dashed", col="red")lines(alpha-2*alpha.se~SR_data$brood_year[-c(1:4)], lty="dashed", col="red")#third modellines(alpha3~SR_data$brood_year[-c(1:4)],type='l', col="blue")lines(alpha3+2*alpha.se3~SR_data$brood_year[-c(1:4)], lty="dashed", col="blue")lines(alpha3-2*alpha.se3~SR_data$brood_year[-c(1:4)], lty="dashed", col="blue")m1$states[-c(1:4)] #alpha from model 1 trimmed to same length as model 2 for comparisonm2$states[1,] #alpha from model 2m3$states[1,] #alpha from model 3m2$states[2,] # beta from model 2m3$states[2,] #beta from model 3m3$states[3,] #estimate of pdo_s coeficient was allowed to vary in Q but was estimated as static#AIC valuem1$AICcm2$AICc #AICc is within a single point.m3$AICc # getting worseprint(mean_produc_m3)```Plot the ts of $\delta_t$ and note the AICc for this model.The plot for state = X3 is the time series of $delta_t$. The model estimates X3 to be 0.03 and static.```{r}autoplot(m3, plot.type ="xtT")```### Model diagnostics```{r}plot(m3)```## Task 4: Expanded Ricker with winter PDOFinally, we fit the expanded model including the winter PDO index as the covariate (`pdo_winter_t2`). ```{r}pdo_w <-matrix(dat[4,-c(1:4)], nrow=1)m <-3#number of parameters in process TT <-length(spawn) #number of data pointsB <-diag(m) #"identity"U <-matrix(0, nrow = m, ncol =1) #"zero"Q <-matrix(list("q.alpha", 0,0, #we no longer want alpha to vary0, 0,0, #Do we want beta to vary now?0,0,"q.pdo_w"), nrow =3, byrow =TRUE) # to have characters and numbers in same matrix use a listZ <-array(NA, c(1, m, TT)) ## NxMxT; empty for nowZ[1,1,] <-rep(1, TT) ## Nx1; 1's for interceptZ[1,2,] <- spawn ## Nx1; predictor variable spawnersZ[1,3,] <- pdo_w ## Nx1; predictor PDO summer#Z[1,4,] <- pdo_w ## Nx1; predictor PDO winterA =matrix(0) #"zero"R <-matrix("r")inits_list4 <-list(x0 =matrix(c(0,0,0), nrow = m)) mod_list4 <-list(B = B, U = U, Q = Q, Z = Z, A = A, R = R)cont_list <-list(maxit=10000)# run modelm4 <-MARSS(ratio, mod_list4, inits = inits_list4, control = cont_list, method ="BFGS")```### $\alpha_t$, $\delta_t$ plots and the AICc.What is the mean level of productivity? Plot the ts of $\delta_t$ and note the AICc for this model. The mean productivity level is 0.28.```{r}alpha4 <-as.numeric(m4$states[1,]) #static numberbeta4 <-as.numeric(m4$states[2,])pdo_w4 <-as.numeric(m4$states[3,])alpha.se4 <-as.numeric(m4$states.se[1,])beta.se4 <-as.numeric(m4$states.se[2,])pdo_s.se4 <-as.numeric(m4$states.se[3,])mean_produc_m4 <-mean(alpha4)print(mean_produc_m4)#from first modelplot(alpha1~SR_data$brood_year,type='l', ylim=c(-3,3), xlim=c(1950, 2005))lines(alpha1+2*alpha.se1~SR_data$brood_year, lty="dashed")lines(alpha1-2*alpha.se1~SR_data$brood_year, lty="dashed")#second modellines(alpha~SR_data$brood_year[-c(1:4)],type='l',col="red")lines(alpha+2*alpha.se~SR_data$brood_year[-c(1:4)], lty="dashed", col="red")lines(alpha-2*alpha.se~SR_data$brood_year[-c(1:4)], lty="dashed", col="red")#third modellines(alpha3~SR_data$brood_year[-c(1:4)],type='l', col="blue")lines(alpha3+2*alpha.se3~SR_data$brood_year[-c(1:4)], lty="dashed", col="blue")lines(alpha3-2*alpha.se3~SR_data$brood_year[-c(1:4)], lty="dashed", col="blue")#fourth modellines(alpha4~SR_data$brood_year[-c(1:4)],type='l', col="green")lines(alpha4+2*alpha.se4~SR_data$brood_year[-c(1:4)], lty="dashed", col="green")lines(alpha4-2*alpha.se4~SR_data$brood_year[-c(1:4)], lty="dashed", col="green")``````{r}autoplot(m4, plot.type ="xtT")```### Model DiagnosticsAlso plot appropriate model diagnostics.```{r}plot(m4)```# Task 5: Best Model Forecasts## Model comparison```{r}#AIC valuem1$AICcm2$AICc #AICc is within a single point.m3$AICc # getting worsem4$AICc```5. Based on AICc, which of the models above is the most parsimonious? Is it well behaved (i.e., are the model assumptions met)? Plot the model forecasts for the best model. Is this a good forecast model?The best model is our first model with just a time varying intercept. However, the model with a static effect of spawners had almost equal support. The model is well behaved and all model assumptions are met. ## Model forecasts```{r}fr1 <-predict(m1, type="ytT", n.ahead=1)plot(fr1)```# Task 6: Other Environmental CovariatesAs our final task, we considered the effect of sea surface tmeperature on the time-varying productivity of Kvichak salmon. ```{r}bb_data <-readRDS(here::here("Lab-1", "Data_Images", "bristol_bay_data_plus_covariates.rds"))sst <- bb_data %>%filter(system =="Kvichak", brood_yr <=2005)%>%select(brood_yr, env_sst)%>%group_by(brood_yr)%>%summarize(sst =mean(env_sst))spawn1957 <- spawn[-1]ratio1957 <- ratio[-1]sst1957 <-t(sst$sst)m <-3#number of parameters in process TT <-length(spawn1957) #number of data pointsB <-diag(m) #"identity"U <-matrix(0, nrow = m, ncol =1) #"zero"Q <-matrix(list("q.alpha", 0,0, 0, 0,0, #Do we want beta to vary now?0,0,"q.sst"), nrow =3, byrow =TRUE) # to have characters and numbers in same matrix use a listZ <-array(NA, c(1, m, TT)) ## NxMxT; empty for nowZ[1,1,] <-rep(1, TT) ## Nx1; 1's for interceptZ[1,2,] <- spawn1957 ## Nx1; predictor variable spawnersZ[1,3,] <- sst1957 ## Nx1; predictor PDO summerA =matrix(0) #"zero"R <-matrix("r")inits_list2 <-list(x0 =matrix(c(0,0,0), nrow = m)) mod_list2 <-list(B = B, U = U, Q = Q, Z = Z, A = A, R = R)# run modelm5 <-MARSS(ratio1957, mod_list2, inits = inits_list2)m5$AICc```### $\alpha_t$, $\delta_t$ plots and the AICc.What is the mean level of productivity? Plot the ts of $\delta_t$ and note the AICc for this model. ```{r}alpha_m5 <-as.numeric(m5$states[1,]) alpha_m5_se <-as.numeric(m5$states.se[1,])sst_m5 <-as.numeric(m5$states[3,])sst_m5_se <-as.numeric(m5$states.se[3,])mean_produc_m5 <-mean(alpha_m5)# Alpha plotplot(alpha_m5~sst$brood_yr,type='l', ylim =c(0,5),xlim=c(1950, 2005))lines((alpha_m5+2*alpha_m5_se)~sst$brood_yr, lty="dashed", col="red")lines(alpha_m5-2*alpha_m5_se~sst$brood_yr, lty="dashed", col="red")text(1960, 0.4, "mean producivity = 2.48")# Delta plotplot(sst_m5~sst$brood_yr,type='l',ylim =c(-1.1, 0.3), xlim=c(1950, 2005))lines((sst_m5+2*sst_m5_se)~sst$brood_yr, lty="dashed", col="red")lines(sst_m5-2*sst_m5_se~sst$brood_yr, lty="dashed", col="red")#AIC valuem5$AICc```The mean level of productivity in this model was 2.48. We found that this model did not test higher than the reduced Ricker model. ### Model DiagnosticsAlso plot appropriate model diagnostics.```{r}plot(m5)```# DiscussionWe found that the Kvichak sockeye spawner recruit time series was best explained by the simplest model with just a time varying intercept. Of the three covariates (summer pdo, winter pdo and SST) we tested, the model with SST had the lowest AIC but not lower than the simplest model. We found that the effect of summer and winter pdo were not time varying. This might be because a time varying intercept explained most of the variance. # Team contributionsDylan wrote the original code for reduced ricker and full ricker models. Maria wrote the code for the PDO models. Liz wrote the code for the SST and coded the models in matrix form. All team members contributed to writing and organization of the lab document.
