• Time series can be useful for identifying structure, improving precision, and accuracy of forecasts

• Modeling multivariate time series

  • e.g. MARSS() function, with each observed time series mapped to a single discrete state

• Using DFA

  • Structure determined by factor loadings

• Making inference about population as a whole involves jointly modeling both time series

• Last week, in talking about Gaussian process models, we showed the number of covariance matrix parameters = \(m*(m+1)/2\)

  • problematic for > 5 or so time series

• MARSS solutions: diagonal matrices or ‘equalvarcov’

• DFA runs into same issues with estimating unconstrained \(R\) matrix

• Sites separated by large distances may be grouped together

• Sites close to one another may be found to have very different dynamics

• Puget Sound Chinook salmon
  • 21 populations generally cluster into 2-3 groups based on genetics
  • Historically large hatchery programs
• Hood canal harbor seals
  • Visited by killer whales

• Adjacent sites can be allowed to covary

• Estimated parameters greatly reduced to 2-5

Point referenced data (aka geostatistical)

• Typically 2-D, but could be 1-D or 3-D (depth, altitude)
• May be fixed station or random (e.g. trawl surveys)

Point pattern data

• Spatially referenced based on outcomes (e.g. presence)
• Inference focused on describing clustering (or not)

Areal data

• Locations occur in blocks
• counties, management zones, etc.

CAR (conditionally autoregressive models)

• Better suited for Bayesian methods
• Goal of both is to write the distribution of a single prediction analytically in terms of the joint (y1, y2)

SAR (simultaneous autregressive models)

• Better suited for ML methods
• Simultaneously model distribution of predicted values

‘Autoregressive’ in the sense of spatial dependency / correlation between locations

\[{ Y }_{ i }=\textbf{b}{ X }_{ i }+{ \phi }_{ i }+{ \varepsilon }_{ i }\]

\({ X }_{ i }\) are predictors (regression) \({ \phi }_{ i }i\) spatial component, (aka markov random field) \({ \varepsilon }_{ i }\) residual error term

• Create spatial adjacency matrix W, based on neighbors, e.g. 
• W(i,j) = 1 if neighbors, 0 otherwise
• W often row-normalized (rows sum to 1)
• Diagonal elements W(i,i) are 0

In matrix form,

\[{ \phi\sim N\left( 0,{ \left( I-\rho W \right) }^{ -1 }\widetilde { D } \right) \\ { \widetilde { D } }_{ ii }={ \sigma }_{ i } }\] * Implemented in ‘spdep’, ‘CARBayes’, etc

• Each element has conditional distribution dependent on others,

\(\phi_i \mid \phi_j, \sim \mathrm{N} \left( \sum_{j = 1}^n W_{ij} \phi_j, {\sigma}^2 \right)\) for \(j \neq i\)

• Simultaneous autoregressive model

\[{ \phi \sim N\left( 0,{ \left( I-\rho W \right) }^{ -1 }\widetilde { D } { \left( I-\rho W \right) }^{ -1 } \right) \\ { \widetilde { D } }_{ ii }={ \sigma }_{ i } }\]

• Remember that the CAR was

\[{ \phi \sim N\left( 0,{ \left( I-\rho W \right) }^{ -1 }\widetilde { D } \right) \\ { \widetilde { D } }_{ ii }={ \sigma }_{ i } }\]

• \(I-\rho W\)

• Example \(W\) matrix,

##      [,1] [,2] [,3]
## [1,]    0    1    0
## [2,]    1    0    1
## [3,]    0    1    0

• \((I-0.3 W)\), \(\rho=0.3\)

##      [,1] [,2] [,3]
## [1,]  1.0 -0.3  0.0
## [2,] -0.3  1.0 -0.3
## [3,]  0.0 -0.3  1.0

• \((I-0.3 W)\widetilde { D }\)
• Let’s assume same variance ~ 0.1
• D = diag(0.1,3)

##       [,1]  [,2]  [,3]
## [1,]  0.10 -0.03  0.00
## [2,] -0.03  0.10 -0.03
## [3,]  0.00 -0.03  0.10

• Adjacency matrix W can also instead be modified to include distance

• Models spatial dependency as a function of single parameter \(\rho\)

• Models don’t include time dimension in spatial field

  • One field estimated for all time steps
  • How could we do this?

?brms::sar
?brms::car

Data from COVID tracking project link

• Suppose we wanted to model log difference of positive % cases

d = dplyr::filter(d, state %in% rownames(W)) %>% 
  dplyr::mutate( 
    pct = positive/totalTestResults) %>% 
  dplyr::filter(yday %in% c(122,366)) %>% 
  dplyr::group_by(state) %>% 
  dplyr::arrange(yday) %>% 
  dplyr::summarize(diff = log(pct[2])-log(pct[1]))

# fit a CAR model
fit <- brm(diff ~ 1 + car(W), 
           data = d, data2 = list(W = W), chains=1) 
summary(fit)

Wall (2004) “A close look at the spatial structure implied by the CAR and SAR models”.

• Note: not a 1-to-1 mapping of distance and correlation

• model elements of Q as functions
• Create matrix D, as pairwise distances
• This can be 1-D, or any dimension
  • We’ll use Euclidian distances for 2-D

Should spatial dependency be included? If so, how to model it?

• Constant
• Time varying
• Autoregressive
• Random walk
• Independent variation

1. Generalized least squares

2. Bayesian methods in spBayes, glmmfields

3. INLA models

4. Spatial GAMs

5. TMB (VAST, sdmTMB)

• Generalized least squares function
  • similar syntax to lm, glm, etc.
• Flexible correlation structures
  • corExp()
  • corGaus()
  • corLin()
  • corSpher()
• Allows irregularly spaced data / NAs
  • unlike Arima(), auto.arima(), etc.

We’ll use Snow Water Equivalent (SWE) data in Washington state

70 SNOTEL sites

• we’ll focus only on Cascades

1981-2013

Initially start using just the February SWE data

1518 data points (only 29 missing values!)

mod.exp = gls(Feb ~ elev, 
  correlation = corExp(form=~lat+lon,nugget=T), 
  data = y[which(is.na(y$Feb)==F & y$Water.Year==2013),])
AIC(mod.exp) = 431.097

mod.gaus = gls(Feb ~ elev, 
  correlation = corGaus(form=~lat+lon,nugget=T), 
  data = y[which(is.na(y$Feb)==F & y$Water.Year==2013),])
AIC(mod.gaus) = 433.485

var.exp <- Variogram(mod.exp, form =~ lat+lon)
plot(var.exp,main="Exponential",ylim=c(0,1))

var.gaus <- Variogram(mod.gaus, form =~ lat+lon)
plot(var.gaus,main="Gaussian",ylim=c(0,1))

Semivariance = \(0.5 \quad Var({x}_{1},{x}_{2})\)
Semivariance = Sill - \(Cov({x}_{1},{x}_{2})\)

corExp and corGaus spatial structure useful for wide variety of models / R packages

Linear/non-linear mixed effect models

• lme() / nlme() in nlme package

Generalized linear mixed models

• glmmPQL() in MASS package

Generalized additive mixed models

• gamm() in mgcv package

Previous approaches modeled errors (or random effects) as correlated

• Spatial GAMs generally model mean as spatially correlated

First a simple GAM, with latitude and longtide.

• Note we’re not transforming to UTM, but probably should
• Note we’re not including tensor smooths te(), but probably should

d = dplyr::filter(d, !is.na(Water.Year), !is.na(Feb)) 
mod = gam(Feb ~ s(Water.Year) + 
    s(Longitude, Latitude), data=d)
d$resid = resid(mod)

First, residuals through time

• We’ll use a subset of years here

• What do these interactions mean in the context of our SNOTEL data?

1. First, spatial field may vary by year

mod2 = gam(Feb ~ s(Water.Year) + 
    s(Longitude, Latitude, by=as.factor(Water.Year)), data=d)
d$resid = resid(mod2)

1. Spatial effect may be continuous and we can include space-time interaction

• Why ti() instead of te() ? ti() only is the interaction where te() includes the main effects

mod2 = gam(Feb ~ s(Longitude, Latitude) + s(Water.Year) + 
    ti(Longitude, Latitude,Water.Year, d=c(2,1)), data=d)
d$resid = resid(mod2)

Quickly evolving features for spatio-temporal models

• Random effects gamm
• Interface with inla ginla
• Bayesian estimation bam

Lots of detailed examples and resources

*e.g. ESA workshop by Simpson/Pederson/Ross/Miller

https://noamross.github.io/mgcv-esa-workshop/

• Hard to specify covariance structure explicitly

• Hard to add constraints, like making spatial field be an AR(1) process

• Motivates more custom approaches

• Package on CRAN, spBayes

• Follows general form of CAR / SAR models, in doing spatial regression

\[\textbf{ Y }=\textbf{b}\textbf{ X }+{ \omega }+{ \varepsilon }\]
* \(\omega\) is some spatial field / process of interest
* \(\epsilon\) is the residual error

• spLM / spGLM

• spMvLM / spMvGLM

• spDynLM: space continuous, time discrete

• Slightly more complicated syntax than functions we’ve been working with:

Specify:

• Priors on parameters
• Tuning parameters for Metropolis sampling (jumping variance)
• Starting / initial values
• Covariance structure (“exponential”, “gaussian”, “matern”, etc)
• Number of MCMC iterations / burn-in, etc.

# This syntax is dependent on model parameters. See vignette
priors <- list("beta.Norm"=list(rep(0,p), 
  diag(1000,p)), "phi.Unif"=c(3/1, 3/0.1), 
  "sigma.sq.IG"=c(2, 2), "tau.sq.IG"=c(2, 0.1))

# Phi is spatial scale parameter, sigma.sq is spatial variance, 
# tau.sq = residual
starting <- list("phi"=3/0.5, "sigma.sq"=50, "tau.sq"=1)

# variance of normal proposals for Metropolis algorithm
tuning <- list("phi"=0.1, "sigma.sq"=0.1, "tau.sq"=0.1)

m.1 <- spLM(y~X-1, coords=cords, 
  n.samples=10000, 
  cov.model = "exponential", priors=priors, 
  tuning=tuning, starting=starting)

##recover beta and spatial random effects
burn.in <- 5000
m.1 <- spRecover(m.1, start=burn.in, verbose=FALSE)

Output is of class ‘mcmc’

Nelson et al. 2018

Joseph 2016

Smith & Edwards 2020

Finley & Banerjee 2020