4  Climate Variables

Author

E Holmes

Published

April 22, 2025

Code 
library(tidyverse)
options(dplyr.summarise.inform = FALSE)

Climate Indices

See also https://psl.noaa.gov/data/atmoswrit/corr/timeseries/

Index Region Influences Notes Data Source
PDO (Pacific Decadal Osc.) North Pacific Ocean temps, productivity, stratification Warm phase → poor survival in Alaska, good in CA (and vice versa) PDO CSV
ONI (Oceanic Niño Index) Equatorial Pacific (global effects) Upwelling, SST, marine food web El Niño often harmful for juvenile salmon on West Coast ONI TXT
NPGO (North Pacific Gyre Osc) North Pacific Salinity, nutrient flux, productivity Positively correlated with salmon survival (esp. CA Current) NPGO TXT
Cold Pool Index Bering Sea shelf Bottom temperature, habitat compression Cold pool extent influences migration, prey availability, and predation coldpool r package
NEP SST (NE Pacific SST) Coastal Alaska, BC, WA, OR, CA Marine heatwaves, thermal stress associated with reduced survival
Upwelling Index (Bakun) California/Oregon coast Nutrient delivery, primary productivity Drives ocean growth conditions for smolts NOAA Upwelling
PNA (Pacific-North American) North America/western U.S. Winter temps, river flow, snowpack Influences freshwater survival, stream temperatures info PNA TXT
MEI / MEIv2 (Multivariate ENSO Index) ENSO-linked Broad atmospheric/ocean composite More integrated than ONI; includes wind, pressure, SST MEIv2 TXT
West Coast Stream Flow Daily West Coast Stream Flow

Loading data

This can be painful. Share your code in the class repo so others can use your work! Post here

Code 
pdo_url <- "https://www.ncei.noaa.gov/pub/data/cmb/ersst/v5/index/ersst.v5.pdo.dat"
pdo <- read.table(pdo_url, skip = 2, header = FALSE, check.names = FALSE)
pdo[pdo==99.99] <- NA
pdo_vec <- as.numeric(t(pdo[, 2:13]))
pdots <- ts(pdo_vec, freq=12, start=c(1854, 1))
plot(pdots)

Prepping your covariates for MARSS()

You need a matrix:

  • Columns = time and number of columns must match your data
  • Rows = 1 row for each covariate

Example for PDO

Let’s say we want to use the April-June average PDO.

Code 
# col 1 is year
years <- which(pdo[,1] %in% 2000:2025)
spring_pdo <- rowMeans(pdo[years, 5:7], na.rm = TRUE)
cmat <- matrix(spring_pdo, nrow=1)
cmat
           [,1]      [,2]      [,3] [,4]      [,5] [,6]       [,7]       [,8]
[1,] -0.6933333 -1.173333 -1.236667 0.05 0.1666667 0.83 -0.2333333 -0.2733333
          [,9]     [,10]      [,11] [,12]     [,13] [,14] [,15] [,16]    [,17]
[1,] -1.753333 -1.426667 -0.3366667 -0.92 -1.586667 -0.77  0.45  0.68 1.403333
         [,18]      [,19]      [,20]      [,21] [,22]     [,23]     [,24]
[1,] 0.4466667 -0.6433333 0.09333333 -0.8666667  -1.9 -1.883333 -2.673333
         [,25] [,26]
[1,] -2.746667   NaN

This is now \(\mathbf{c}_t\) in \[ \mathbf{x}_t = \mathbf{x}_{t-1} + \mathbf{u} + \mathbf{C}\mathbf{c}_t + \mathbf{w}_t \]

Use PDO to predict steelhead returns

Idea is that pdo (\(p_t\)) is correlated with the fluctuations in log spawner returns (\(x_t\)).

\[ x_t = x_{t-1} + u + (\alpha p_t + w_t) \] But we have 4 steelhead time series for this example. We can write this as: \[ \begin{bmatrix}x_{1,t}\\x_{2,t}\\x_{3,t}\\x_{4,t}\end{bmatrix}= \begin{bmatrix}x_{1,t-1}\\x_{2,t-1}\\x_{3,t-1}\\x_{4,t-1}\end{bmatrix} + \begin{bmatrix}u_1\\u_2\\u_3\\u_4\end{bmatrix}+ \begin{bmatrix}\alpha_1\\\alpha_2\\\alpha_3\\\alpha_4\end{bmatrix} p_t + \begin{bmatrix}w_1\\w_2\\w_3\\w_4\end{bmatrix} \]

What about the \(w_t\)’s? How are they related? Let’s say they are independent with unequal variances (the \(x_t\) are diagonal and unequal). What about our observations? Let’s say each \(y_t\) measures one \(x_t\) with independent and equal error (R is diagonal and equal). Let’s fit this model with MARSS.

Wrangle the raw spawner data to get the \(x_t\).

Code 
library(tidyverse)
load("Data_Images/esa-salmon.rda")
esu <- unique(esa.salmon$esu_dps)
esunum <- which(esu == "Steelhead (Upper Columbia River DPS)")
esuname <- esu[esunum]

dat <- esa.salmon %>% 
  subset(esu_dps == esuname) %>% # get only this ESU
  mutate(log.spawner = log(value)) %>% # create a column called log.spawner
  select(esapopname, spawningyear, log.spawner) %>% # get just the columns that I need
  pivot_wider(names_from = "esapopname", values_from = "log.spawner") %>% 
  column_to_rownames(var = "spawningyear") %>% # make the years rownames
  as.matrix() %>% # turn into a matrix with year down the rows
  t() # make time across the columns
# MARSS complains if I don't do this
dat[is.na(dat)] <- NA
# Clean up rownames
tmp <- rownames(dat)
tmp <- stringr::str_replace(tmp, "Steelhead [(]Upper Columbia River DPS[)]", "")
tmp <- stringr::str_replace(tmp, "River - summer", "")
tmp <- stringr::str_trim(tmp)
rownames(dat) <- tmp

Get the PDO for the years in our data:

Code 
# col 1 is year
years <- which(pdo[,1] %in% colnames(dat))
spring_pdo <- rowMeans(pdo[years, 5:7], na.rm = TRUE)
cmat <- matrix(spring_pdo, nrow=1)

Check that cmat and dat have the same number of columns.

Code 
dim(cmat)
[1]  1 63
Code 
dim(dat)
[1]  4 63

Specify a model

Code 
# pdo
mod.list1 <- list(
  U = "unequal",
  R = "diagonal and equal",
  Q = "diagonal and unequal",
  c = cmat,
  C = "unequal"
)

# no pdo
mod.list2 <- list(
  U = "unequal",
  R = "diagonal and equal",
  Q = "diagonal and unequal")
Code 
library(MARSS)
Warning: package 'MARSS' was built under R version 4.4.3
Code 
fit1 <- MARSS(dat, model=mod.list1)
Success! abstol and log-log tests passed at 229 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 229 iterations. 
Log-likelihood: -152.4383 
AIC: 338.8765   AICc: 343.0398   
 
                            Estimate
R.diag                        0.1005
U.X.Entiat                   -0.0144
U.X.Methow                   -0.0122
U.X.Okanogan                 -0.0462
U.X.Wenatchee                -0.0140
Q.(X.Entiat,X.Entiat)         0.1178
Q.(X.Methow,X.Methow)         0.1200
Q.(X.Okanogan,X.Okanogan)     0.1045
Q.(X.Wenatchee,X.Wenatchee)   0.3614
x0.X.Entiat                   5.1194
x0.X.Methow                   6.9854
x0.X.Okanogan                 7.0155
x0.X.Wenatchee                8.1083
C.X.Entiat                   -0.0789
C.X.Methow                   -0.0516
C.X.Okanogan                 -0.0973
C.X.Wenatchee                 0.0377
Initial states (x0) defined at t=0

Standard errors have not been calculated. 
Use MARSSparamCIs to compute CIs and bias estimates.
Code 
fit2 <- MARSS(dat, model=mod.list2)
Success! abstol and log-log tests passed at 304 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 304 iterations. 
Log-likelihood: -154.5953 
AIC: 335.1906   AICc: 337.6012   
 
                            Estimate
R.diag                       0.08492
U.X.Entiat                   0.00213
U.X.Methow                  -0.00557
U.X.Okanogan                -0.00920
U.X.Wenatchee               -0.02427
Q.(X.Entiat,X.Entiat)        0.14095
Q.(X.Methow,X.Methow)        0.14800
Q.(X.Okanogan,X.Okanogan)    0.14483
Q.(X.Wenatchee,X.Wenatchee)  0.38490
x0.X.Entiat                  5.57099
x0.X.Methow                  7.40887
x0.X.Okanogan                7.10616
x0.X.Wenatchee               8.06274
Initial states (x0) defined at t=0

Standard errors have not been calculated. 
Use MARSSparamCIs to compute CIs and bias estimates.

Which one fits better based on AIC? The model without PDO is better (lower AICc).

Code 
fit1$AICc
[1] 343.0398
Code 
fit2$AICc
[1] 337.6012

Resources

Chapter 7 MARSS models. ATSA Lab Book. https://atsa-es.github.io/atsa-labs/chap-mss.html

Chapter 8 MARSS models with covariates. ATSA Lab Book. https://atsa-es.github.io/atsa-labs/chap-msscov.html

Chapter 16 Modeling cyclic sockeye https://atsa-es.github.io/atsa-labs/chap-cyclic-sockeye.html