Loading Web-Font TeX/Math/Italic

Jump to main navigation


Tutorial 8.3b - Dealing with temporal autocorrelation

05 Feb 2018

Preparation

The following packages will be required.
library(tidyverse)
library(coda)
library(R2jags)
library(gridExtra)
library(broom)

Important opening note

Dealing with temporal autocorrelation and analysing temporal trends are not the same thing. The former attempts to counter the lack of independence associated with temporal data whereas the later attempts to model the influence of temporal patterns. This tutorial will focus only on temporal autocorrelation, time series analyses will be the focus of another tutorial.

Overview

The following representation of the standard linear model highlights the various assumptions that are imposed upon the underlying data. The current tutorial will focus directly on issues related to the nature of residual variability.

Please make sure you are familiar with the general concepts and frequentist approaches to dealing with temporal autocorrelation as outlined in Tutorial 8.3a.

Recall from the assumptions outlined in Tutorial 8.3a, in order that the estimated parameters represent the underlying populations in an unbiased manner, the residuals (and thus each each observation) must be independent.

However, what if we were sampling a population over time and we were interested in investigating how changes in a response relate to changes in a predictor (such as rainfall). For any response that does not 'reset' itself on a regular basis, the state of the population (the value of its response) at a given time is likely to be at least partly dependent on the state of the population at the sampling time before.

We can further generalized the above into: \begin{array}{rcl} y_{i}&\sim&\mathcal{Dist}(\mu_{i})\\ \mu_{i}&=&\mathbf{X}\boldsymbol{\beta} + \mathbf{Z}\boldsymbol{\gamma}\\ \end{array}

where \mathbf{X} and \boldsymbol{\beta} represent the 'fixed' data structure (model matrix) and 'fixed' effects respectively and \mathbf{Z} and \boldsymbol{\gamma} represent the 'varying' (=random in frequentist) data structure and 'varying' effects (including residuals) respectively.

In simple regression, there are no 'varying' effects, and thus: \begin{array}{rcl} \boldsymbol{\gamma}&\sim&\mathcal{MVN}(0,\Sigma)\\ \end{array}

where \mathcal{MVN} is a multi-variate normal and \Sigma represents a variance-covariance matrix of the form: \Sigma =\frac{\sigma^2}{1-\rho^2} \begin{pmatrix} 1&\rho^{\phi_{1,2}}&\cdots&\rho^{\phi_{1,n}}\\ \rho^{\phi_{2,1}}&1&\cdots&\vdots\\ \vdots&\cdots&1&\vdots\\ \rho^{\phi_{n,1}}&\cdots&\cdots&1\\ \end{pmatrix}

Notice that this introduces a very large number of additional parameters that require estimating: \sigma^2 (error variance), \rho (base autocorrelation) and each of the individual covariances (\rho^{\phi_{n,n}}). Hence, there are always going to be more parameters to estimate than there are date avaiable to use to estimate these paramters.

We typically make one of a number of alternative assumptions so as to make this task more manageable.

  • When we assume that all residuals are independent (regular regression), ie \rho=0, \Sigma essential is equal to \sigma^2 I and we simply use: \boldsymbol{\gamma} \sim{} \mathcal{N}(0,\sigma^2)
  • We could assume there is a reasonably simple pattern of correlation that declines over time. The simplest of these is a first order autoregressive (AR1) structure in which exponent on the correlation declines linearly according to the time lag (|t-s|). \Sigma =\frac{\sigma^2}{1-\rho^2} \begin{pmatrix} 1&\rho&\cdots&\rho^{|t-s|}\\ \rho&1&\cdots&\vdots\\ \vdots&\cdots&1&\vdots\\ \rho^{|t-s|}&\cdots&\cdots&1\\ \end{pmatrix}
    Note, in making this assumption, we are also assuming that the degree of correlation is dependent only on the lag and not on when the lag occurs (stationarity). That is all lag 1 residual pairs will have the same degree of correlation, all the lag 2 pairs will have the same correlation and so on.

First order autocorrelation

Consider an example, in which the number of individuals at time 2 will be partly dependent on the number of individuals present at time 1. Clearly then, the observations (and thus residuals) are not fully independent - there is an auto-regressive correlation dependency structure. We could accommodate this lack of independence by fitting a model that incorporates a AR1 variance-covariance structure. Alternatively, we fit the following model: \begin{array}{rcl} y_{it}&\sim&\mathcal{Dist}(\mu_{it})\\ \mu_{it}&=&\mathbf{X}\beta + \rho\times\varepsilon_{i,t-1}+\gamma_{it}\\ \gamma&\sim&N(0,\sigma^2) \end{array}

In this version of the model, we are stating that the expected value of an observation is equal to the regular linear predictor plus the autocorrelation parameter (\rho) multipled by the residual associated with the previous observation plus the regular independently distributed noise (\sigma^2). Such a model is substantially faster to fit, although along with stationarity assumes in estimating the autocorrelation parameter, only the smallest lags are used.

To see this in action, we will first generate some temporally auto-correlated data.

set.seed(16)
n = 50
a <- 20  #intercept
b <- 0.2  #slope
x <- round(runif(n, 1, n), 1)  #values of the year covariate
year <- 1:n
sigma <- 20
rho <- 0.8

## define a constructor for a first-order
## correlation structure
ar1 <- corAR1(form = ~year, value = rho)
## initialize this constructor against our data
AR1 <- Initialize(ar1, data = data.frame(year))
## generate a correlation matrix
V <- corMatrix(AR1)
## Cholesky factorization of V
Cv <- chol(V)
## simulate AR1 errors
e <- t(Cv) %*% rnorm(n, 0, sigma)  # cov(e) = V * sig^2
## generate response
y <- a + b * x + e
data.temporalCor = data.frame(y = y, x = x, year = year)
write.table(data.temporalCor, file = "../downloads/data/data.temporalCor.csv",
    sep = ",", quote = F, row.names = FALSE)
pairs(data.temporalCor)
plot of chunk tut8.3bS5.1a

We will now proceed to analyse these data via both of the above techniques for each of JAGS, RSTAN and BRMS:

  • incorporating AR1 residual autocorrelation structure
  • incorporating lagged residuals into the model

Incorporating lagged residuals

  • Define the model
    modelString = " 
      model {
      #Likelihood
      for (i in 1:n) {
      fit[i] <- inprod(beta[],X[i,])
      y[i] ~ dnorm(mu[i],tau.cor)
      }
      e[1] <- (y[1] - fit[1])
      mu[1] <- fit[1]
      for (i in 2:n) {
      e[i] <- (y[i] - fit[i]) #- phi*e[i-1]
      mu[i] <- fit[i] + phi * e[i-1]
      }
      #Priors
      phi ~ dunif(-1,1)
      for (i in 1:nX) {
      beta[i] ~ dnorm(0,1.0E-6)
      }
      sigma <- z/sqrt(chSq)    # prior for sigma; cauchy = normal/sqrt(chi^2)
      z ~ dnorm(0, 0.04)I(0,)
      chSq ~ dgamma(0.5, 0.5)  # chi^2 with 1 d.f.
      tau <- pow(sigma, -2)
      tau.cor <- tau #* (1- phi*phi)
      }
      "
    
  • Arrange the data
    Xmat = model.matrix(~x, data.temporalCor)
    data.temporalCor.list <- with(data.temporalCor, list(y = y, X = Xmat,
        n = nrow(data.temporalCor), nX = ncol(Xmat)))
    
  • Fit the model
    library(R2jags)
    
    data.temporalCor.r2jags <- jags(data = data.temporalCor.list,
        inits = NULL, parameters.to.save = c("beta", "sigma", "phi"),
        model.file = textConnection(modelString), n.chains = 3, n.iter = 10000,
        n.burnin = 5000, n.thin = 10)
    
    Compiling model graph
       Resolving undeclared variables
       Allocating nodes
    Graph information:
       Observed stochastic nodes: 50
       Unobserved stochastic nodes: 5
       Total graph size: 421
    
    Initializing model
    
  • Explore MCMC mixing diagnostics
    library(coda)
    data.mcmc = as.mcmc(data.temporalCor.r2jags)
    plot(data.mcmc)
    
    plot of chunk tut8.3bbFitJAGS.AR.4
    plot of chunk tut8.3bbFitJAGS.AR.4
    autocorr.diag(data.mcmc)
    
                 beta[1]      beta[2]     deviance
    Lag 0    1.000000000  1.000000000  1.000000000
    Lag 10   0.001092035 -0.018372250 -0.005516129
    Lag 50   0.002004868  0.003609002 -0.020406834
    Lag 100  0.018638496 -0.027413117  0.004258020
    Lag 500 -0.003754928  0.004363078  0.004921508
                     phi        sigma
    Lag 0    1.000000000  1.000000000
    Lag 10  -0.012457395  0.005843107
    Lag 50  -0.001470570  0.014372526
    Lag 100 -0.014432385  0.023423489
    Lag 500  0.006227912 -0.002678010
    
    Conclusions - all diagnostics seem fine.
  • Perform model validation.

    Whenever we fit a model that incorporates changes to the variance-covariance structures, we need to explore modified standardized residuals. In this case, the raw residuals should be updated to reflect the autocorrelation (subtract residual from previous time weighted by the autocorrelation parameter) before standardizing by sigma \begin{align} Res_i =& Y_{i} - \mu_{i}\\ Res_{i+1} =& Res_{i+1} - \rho Res_{i}\\ Res_i =& \frac{Res_i}{\sigma}\\ \end{align}

    mcmc = data.temporalCor.r2jags$BUGSoutput$sims.matrix
    # generate a model matrix
    newdata = data.temporalCor
    Xmat = model.matrix(~x, newdata)
    ## get median parameter estimates
    wch = grep("beta", colnames(mcmc))
    coefs = mcmc[, wch]
    fit = coefs %*% t(Xmat)
    resid = -1 * sweep(fit, 2, data.temporalCor$y, "-")
    n = ncol(resid)
    resid[, -1] = resid[, -1] - (resid[, -n] * mcmc[, "phi"])
    resid = apply(resid, 2, median)/median(mcmc[, "sigma"])
    fit = apply(fit, 2, median)
    ggplot() + geom_point(data = NULL, aes(y = resid, x = fit))
    
    plot of chunk tut8.3bbFitJAGS.AR.5
    ggplot() + geom_point(data = NULL, aes(y = resid, x = data.temporalCor$x))
    
    plot of chunk tut8.3bbFitJAGS.AR.5
    ggplot(data = NULL, aes(y = resid, x = data.temporalCor$year)) +
        geom_point() + geom_line() + geom_hline(yintercept = 0, linetype = "dashed")
    
    plot of chunk tut8.3bbFitJAGS.AR.5
    plot(acf(resid, lag = 40))
    
    plot of chunk tut8.3bbFitJAGS.AR.5
    Conclusions - no obvious autocorrelation or other issues with residuals remaining
  • Explore parameter estimates
    tidyMCMC(as.mcmc(data.temporalCor.r2jags), conf.int = TRUE, conf.method = "HPDinterval")
    
          term    estimate   std.error    conf.low
    1  beta[1]  34.7359826 11.98372345  12.4617763
    2  beta[2]   0.3630489  0.08829582   0.1864085
    3 deviance 387.1079201  2.57015258 383.6776550
    4      phi   0.8923788  0.07090857   0.7570717
    5    sigma  11.6040465  1.19635680   9.4197453
        conf.high
    1  57.9110641
    2   0.5271743
    3 392.4009332
    4   0.9993796
    5  14.1018961
    
  • Define the model
    stanString = "
      data {
      int<lower=1> n;
      vector [n] y;
      int<lower=1> nX;
      matrix[n,nX] X;
      }
      transformed data {
      }
      parameters {
      vector[nX] beta;
      real<lower=0> sigma;
      real<lower=-1,upper=1> phi;
      }
      transformed parameters {
      vector[n] mu;
      vector[n] epsilon;
      mu = X*beta;
      epsilon[1] = y[1] - mu[1];
      for (i in 2:n) {
      epsilon[i] = (y[i] - mu[i]);
      mu[i] = mu[i] + phi*epsilon[i-1];
      }
      }
      model {
      phi ~ uniform(-1,1);
      beta ~ normal(0,100);
      sigma ~ cauchy(0,5);
      y ~ normal(mu, sigma);
      }
      generated quantities {
      }
      "
    
  • Arrange the data
    Xmat = model.matrix(~x, data.temporalCor)
    data.temporalCor.list <- with(data.temporalCor, list(y = y, X = Xmat,
        n = nrow(data.temporalCor), nX = ncol(Xmat)))
    
  • Fit the model
    library(rstan)
    
    data.temporalCor.rstan <- stan(data = data.temporalCor.list,
        model_code = stanString, chains = 3, iter = 2000, warmup = 500,
        thin = 3)
    
    In file included from /usr/local/lib/R/site-library/BH/include/boost/config.hpp:39:0,
                     from /usr/local/lib/R/site-library/BH/include/boost/math/tools/config.hpp:13,
                     from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core/var.hpp:7,
                     from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
                     from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core.hpp:12,
                     from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/mat.hpp:4,
                     from /usr/local/lib/R/site-library/StanHeaders/include/stan/math.hpp:4,
                     from /usr/local/lib/R/site-library/StanHeaders/include/src/stan/model/model_header.hpp:4,
                     from file68ea474eec5a.cpp:8:
    /usr/local/lib/R/site-library/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
     #  define BOOST_NO_CXX11_RVALUE_REFERENCES
     ^
    <command-line>:0:0: note: this is the location of the previous definition
    
    SAMPLING FOR MODEL '4ca8d7c936cc998c61c008bee11475a5' NOW (CHAIN 1).
    
    Gradient evaluation took 3e-05 seconds
    1000 transitions using 10 leapfrog steps per transition would take 0.3 seconds.
    Adjust your expectations accordingly!
    
    
    Iteration:    1 / 2000 [  0%]  (Warmup)
    Iteration:  200 / 2000 [ 10%]  (Warmup)
    Iteration:  400 / 2000 [ 20%]  (Warmup)
    Iteration:  501 / 2000 [ 25%]  (Sampling)
    Iteration:  700 / 2000 [ 35%]  (Sampling)
    Iteration:  900 / 2000 [ 45%]  (Sampling)
    Iteration: 1100 / 2000 [ 55%]  (Sampling)
    Iteration: 1300 / 2000 [ 65%]  (Sampling)
    Iteration: 1500 / 2000 [ 75%]  (Sampling)
    Iteration: 1700 / 2000 [ 85%]  (Sampling)
    Iteration: 1900 / 2000 [ 95%]  (Sampling)
    Iteration: 2000 / 2000 [100%]  (Sampling)
    
     Elapsed Time: 0.072464 seconds (Warm-up)
                   0.088587 seconds (Sampling)
                   0.161051 seconds (Total)
    
    
    SAMPLING FOR MODEL '4ca8d7c936cc998c61c008bee11475a5' NOW (CHAIN 2).
    
    Gradient evaluation took 1.4e-05 seconds
    1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
    Adjust your expectations accordingly!
    
    
    Iteration:    1 / 2000 [  0%]  (Warmup)
    Iteration:  200 / 2000 [ 10%]  (Warmup)
    Iteration:  400 / 2000 [ 20%]  (Warmup)
    Iteration:  501 / 2000 [ 25%]  (Sampling)
    Iteration:  700 / 2000 [ 35%]  (Sampling)
    Iteration:  900 / 2000 [ 45%]  (Sampling)
    Iteration: 1100 / 2000 [ 55%]  (Sampling)
    Iteration: 1300 / 2000 [ 65%]  (Sampling)
    Iteration: 1500 / 2000 [ 75%]  (Sampling)
    Iteration: 1700 / 2000 [ 85%]  (Sampling)
    Iteration: 1900 / 2000 [ 95%]  (Sampling)
    Iteration: 2000 / 2000 [100%]  (Sampling)
    
     Elapsed Time: 0.111481 seconds (Warm-up)
                   0.149398 seconds (Sampling)
                   0.260879 seconds (Total)
    
    
    SAMPLING FOR MODEL '4ca8d7c936cc998c61c008bee11475a5' NOW (CHAIN 3).
    
    Gradient evaluation took 1.5e-05 seconds
    1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
    Adjust your expectations accordingly!
    
    
    Iteration:    1 / 2000 [  0%]  (Warmup)
    Iteration:  200 / 2000 [ 10%]  (Warmup)
    Iteration:  400 / 2000 [ 20%]  (Warmup)
    Iteration:  501 / 2000 [ 25%]  (Sampling)
    Iteration:  700 / 2000 [ 35%]  (Sampling)
    Iteration:  900 / 2000 [ 45%]  (Sampling)
    Iteration: 1100 / 2000 [ 55%]  (Sampling)
    Iteration: 1300 / 2000 [ 65%]  (Sampling)
    Iteration: 1500 / 2000 [ 75%]  (Sampling)
    Iteration: 1700 / 2000 [ 85%]  (Sampling)
    Iteration: 1900 / 2000 [ 95%]  (Sampling)
    Iteration: 2000 / 2000 [100%]  (Sampling)
    
     Elapsed Time: 0.088959 seconds (Warm-up)
                   0.086188 seconds (Sampling)
                   0.175147 seconds (Total)
    
  • Explore MCMC mixing diagnostics
    library(coda)
    stan_trace(data.temporalCor.rstan)
    
    plot of chunk tut8.3bbFitRSTAN.AR.4
    stan_ac(data.temporalCor.rstan)
    
    plot of chunk tut8.3bbFitRSTAN.AR.4
    stan_rhat(data.temporalCor.rstan)
    
    plot of chunk tut8.3bbFitRSTAN.AR.4
    stan_ess(data.temporalCor.rstan)
    
    plot of chunk tut8.3bbFitRSTAN.AR.4
    Conclusions - all diagnostics seem fine.
  • Perform model validation.

    Whenever we fit a model that incorporates changes to the variance-covariance structures, we need to explore modified standardized residuals. In this case, the raw residuals should be updated to reflect the autocorrelation (subtract residual from previous time weighted by the autocorrelation parameter) before standardizing by sigma \begin{align} Res_i =& Y_{i} - \mu_{i}\\ Res_{i+1} =& Res_{i+1} - \rho Res_{i}\\ Res_i =& \frac{Res_i}{\sigma}\\ \end{align}

    mcmc = as.matrix(data.temporalCor.rstan)
    wch = grep("beta", colnames(mcmc))
    # generate a model matrix
    newdata = data.frame(x = data.temporalCor$x)
    Xmat = model.matrix(~x, newdata)
    ## get median parameter estimates
    coefs = mcmc[, wch]
    fit = coefs %*% t(Xmat)
    resid = -1 * sweep(fit, 2, data.temporalCor$y, "-")
    n = ncol(resid)
    resid[, -1] = resid[, -1] - (resid[, -n] * mcmc[, "phi"])
    resid = apply(resid, 2, median)/median(mcmc[, "sigma"])
    fit = apply(fit, 2, median)
    ggplot() + geom_point(data = NULL, aes(y = resid, x = fit))
    
    plot of chunk tut8.3bbFitRSTAN.AR.5
    ggplot() + geom_point(data = NULL, aes(y = resid, x = data.temporalCor$x))
    
    plot of chunk tut8.3bbFitRSTAN.AR.5
    ggplot(data = NULL, aes(y = resid, x = data.temporalCor$year)) +
        geom_point() + geom_line() + geom_hline(yintercept = 0, linetype = "dashed")
    
    plot of chunk tut8.3bbFitRSTAN.AR.5
    plot(acf(resid, lag = 40))
    
    plot of chunk tut8.3bbFitRSTAN.AR.5
    fit = coefs %*% t(Xmat)
    ## draw samples from this model
    yRep = sapply(1:nrow(mcmc), function(i) rnorm(nrow(data.temporalCor),
        fit[i, ], mcmc[i, "sigma"]))
    ggplot() + geom_density(data = NULL, aes(x = as.vector(yRep),
        fill = "Model"), alpha = 0.5) + geom_density(data = data.temporalCor,
        aes(x = y, fill = "Obs"), alpha = 0.5)
    
    plot of chunk tut8.3bbFitRSTAN.AR.5
    Conclusions - no obvious autocorrelation or other issues with residuals remaining
  • Explore parameter estimates
    tidyMCMC(data.temporalCor.rstan, par = c("beta", "phi", "sigma"),
        conf.int = TRUE, conf.method = "HPDinterval", rhat = TRUE,
        ess = TRUE)
    
         term   estimate   std.error   conf.low
    1 beta[1] 34.0763301 11.72124133 12.7029023
    2 beta[2]  0.3661828  0.08890794  0.1907748
    3     phi  0.8895282  0.06989449  0.7659398
    4   sigma 11.6031493  1.19587520  9.2675954
       conf.high     rhat  ess
    1 57.6730924 1.000200 1336
    2  0.5442838 1.000230 1470
    3  0.9998038 1.001021 1312
    4 13.7677156 1.000743 1394
    
  • Fit the model
    library(brms)
    
    data.temporalCor.brms = brm(y ~ x, data = data.temporalCor,
        autocor = cor_ar(form = ~year, p = 1), iter = 2000,
        warmup = 200, chains = 3, thin = 2, refresh = 0,
        prior = c(prior(normal(0, 100), class = "Intercept"),
            prior(normal(0, 10), class = "b"), prior(cauchy(0,
                5), class = "sigma")))
    
    Gradient evaluation took 5.1e-05 seconds
    1000 transitions using 10 leapfrog steps per transition would take 0.51 seconds.
    Adjust your expectations accordingly!
    
    
    
     Elapsed Time: 0.195714 seconds (Warm-up)
                   0.320546 seconds (Sampling)
                   0.51626 seconds (Total)
    
    
    Gradient evaluation took 2.5e-05 seconds
    1000 transitions using 10 leapfrog steps per transition would take 0.25 seconds.
    Adjust your expectations accordingly!
    
    
    
     Elapsed Time: 0.164758 seconds (Warm-up)
                   0.478901 seconds (Sampling)
                   0.643659 seconds (Total)
    
    
    Gradient evaluation took 2.8e-05 seconds
    1000 transitions using 10 leapfrog steps per transition would take 0.28 seconds.
    Adjust your expectations accordingly!
    
    
    
     Elapsed Time: 0.202119 seconds (Warm-up)
                   0.378927 seconds (Sampling)
                   0.581046 seconds (Total)
    
  • Explore MCMC mixing diagnostics
    library(coda)
    stan_trace(data.temporalCor.brms$fit)
    
    plot of chunk tut8.3bbFitBRMS.AR.4
    stan_ac(data.temporalCor.brms$fit)
    
    plot of chunk tut8.3bbFitBRMS.AR.4
    stan_rhat(data.temporalCor.brms$fit)
    
    plot of chunk tut8.3bbFitBRMS.AR.4
    stan_ess(data.temporalCor.brms$fit)
    
    plot of chunk tut8.3bbFitBRMS.AR.4
    Conclusions - all diagnostics seem fine.
  • Perform model validation.

    Whenever we fit a model that incorporates changes to the variance-covariance structures, we need to explore modified standardized residuals. In this case, the raw residuals should be updated to reflect the autocorrelation (subtract residual from previous time weighted by the autocorrelation parameter) before standardizing by sigma \begin{align} Res_i =& Y_{i} - \mu_{i}\\ Res_{i+1} =& Res_{i+1} - \rho Res_{i}\\ Res_i =& \frac{Res_i}{\sigma}\\ \end{align}

    resid = resid(data.temporalCor.brms, type = "pearson")[, "Estimate"]
    fit = fitted(data.temporalCor.brms)[, "Estimate"]
    ggplot() + geom_point(data = NULL, aes(y = resid, x = fit))
    
    plot of chunk tut8.3bbFitBRMS.AR.5
    ggplot() + geom_point(data = NULL, aes(y = resid, x = data.temporalCor$x))
    
    plot of chunk tut8.3bbFitBRMS.AR.5
    ggplot(data = NULL, aes(y = resid, x = data.temporalCor$year)) +
        geom_point() + geom_line() + geom_hline(yintercept = 0, linetype = "dashed")
    
    plot of chunk tut8.3bbFitBRMS.AR.5
    plot(acf(resid, lag = 40))
    
    plot of chunk tut8.3bbFitBRMS.AR.5
    y_pred = posterior_predict(data.temporalCor.brms)
    newdata = data.temporalCor %>% cbind(t(y_pred)) %>% gather(key = "Rep",
        value = "Value", -y, -x)
    ggplot(newdata, aes(Value, x = x)) + geom_violin(color = "blue",
        fill = "blue", alpha = 0.5) + geom_violin(data = data.temporalCor,
        aes(y = y, x = x), fill = "red", color = "red", alpha = 0.5)
    
    plot of chunk tut8.3bbFitBRMS.AR.5
    Conclusions - no obvious autocorrelation or other issues with residuals remaining
  • Explore parameter estimates
    tidyMCMC(data.temporalCor.brms$fit, par = c("b_Intercept", "b_x",
        "ar", "sigma"), conf.int = TRUE, conf.method = "HPDinterval",
        rhat = TRUE, ess = TRUE)
    
             term   estimate   std.error   conf.low
    1 b_Intercept 33.6180282 11.42842969 13.5440475
    2         b_x  0.3672319  0.08984579  0.1909742
    3       ar[1]  0.8870894  0.07165605  0.7501506
    4       sigma 11.6119833  1.18136401  9.4317751
       conf.high      rhat  ess
    1 57.1223686 1.0003986 1849
    2  0.5455241 0.9994711 2301
    3  0.9999494 1.0003625 1809
    4 13.9504086 0.9995145 2203
    

incorporating AR1 residual autocorrelation structure

  • Define the model.
    modelString = "
      model {
      #Likelihood
      for (i in 1:n) {
      mu[i] <- inprod(beta[],X[i,])
      }
      y[1:n] ~ dmnorm(mu[1:n],Omega)
      for (i in 1:n) {
      for (j in 1:n) {
      Sigma[i,j] <- sigma2*(equals(i,j) + (1-equals(i,j))*pow(phi,abs(i-j))) 
      }
      }
      Omega <- inverse(Sigma)
      
      #Priors
      phi ~ dunif(-1,1)
      for (i in 1:nX) {
      beta[i] ~ dnorm(0,1.0E-6)
      }
      sigma <- z/sqrt(chSq)    # prior for sigma; cauchy = normal/sqrt(chi^2)
      z ~ dnorm(0, 0.04)I(0,)
      chSq ~ dgamma(0.5, 0.5)  # chi^2 with 1 d.f.
      sigma2 = pow(sigma,2)
      #tau.cor <- tau #* (1- phi*phi)
      }
      "
    
  • Arrange the data
    Xmat = model.matrix(~x, data.temporalCor)
    data.temporalCor.list <- with(data.temporalCor, list(y = y, X = Xmat,
        n = nrow(data.temporalCor), nX = ncol(Xmat)))
    
  • Fit the model. This will require a much higher number of samples, and will take a long time to run! You have been warned.
    library(R2jags)
    
    data.temporalCor.r2jags <- jags(data = data.temporalCor.list,
        inits = NULL, parameters.to.save = c("beta", "sigma", "phi"),
        model.file = textConnection(modelString), n.chains = 3, n.iter = 1e+05,
        n.burnin = 5000, n.thin = 100)
    
    Compiling model graph
       Resolving undeclared variables
       Allocating nodes
    Graph information:
       Observed stochastic nodes: 1
       Unobserved stochastic nodes: 5
       Total graph size: 40526
    
    Initializing model
    
  • Explore MCMC mixing diagnostics
    library(coda)
    data.mcmc = as.mcmc(data.temporalCor.r2jags)
    plot(data.mcmc)
    
    plot of chunk tut8.3bbFitJAGS.AR1.4
    plot of chunk tut8.3bbFitJAGS.AR1.4
    autocorr.diag(data.mcmc)
    
                  beta[1]      beta[2]     deviance
    Lag 0     1.000000000  1.000000000  1.000000000
    Lag 100  -0.020652417 -0.024430350 -0.012402940
    Lag 500   0.011247465 -0.004110036 -0.003832072
    Lag 1000  0.002106766  0.018811984  0.025926215
    Lag 5000  0.014759734  0.002515957 -0.010421871
                     phi        sigma
    Lag 0     1.00000000 1.0000000000
    Lag 100   0.04459154 0.0652759005
    Lag 500   0.03102583 0.0049869372
    Lag 1000 -0.03004530 0.0001472537
    Lag 5000  0.02039965 0.0082833443
    
    Conclusions - all diagnostics seem fine.
  • Perform model validation.

    Whenever we fit a model that incorporates changes to the variance-covariance structures, we need to explore modified standardized residuals. In this case, the raw residuals should be updated to reflect the autocorrelation (subtract residual from previous time weighted by the autocorrelation parameter) before standardizing by sigma \begin{align} Res_i =& Y_{i} - \mu_{i}\\ Res_{i+1} =& Res_{i+1} - \rho Res_{i}\\ Res_i =& \frac{Res_i}{\sigma}\\ \end{align}

    mcmc = data.temporalCor.r2jags$BUGSoutput$sims.matrix
    # generate a model matrix
    newdata = data.temporalCor
    Xmat = model.matrix(~x, newdata)
    ## get median parameter estimates
    wch = grep("beta", colnames(mcmc))
    coefs = mcmc[, wch]
    fit = coefs %*% t(Xmat)
    resid = -1 * sweep(fit, 2, data.temporalCor$y, "-")
    n = ncol(resid)
    resid[, -1] = resid[, -1] - (resid[, -n] * mcmc[, "phi"])
    resid = apply(resid, 2, median)/median(mcmc[, "sigma"])
    fit = apply(fit, 2, median)
    ggplot() + geom_point(data = NULL, aes(y = resid, x = fit))
    
    plot of chunk tut8.3bbFitJAGS.AR1.5
    ggplot() + geom_point(data = NULL, aes(y = resid, x = data.temporalCor$x))
    
    plot of chunk tut8.3bbFitJAGS.AR1.5
    ggplot(data = NULL, aes(y = resid, x = data.temporalCor$year)) +
        geom_point() + geom_line() + geom_hline(yintercept = 0, linetype = "dashed")
    
    plot of chunk tut8.3bbFitJAGS.AR1.5
    plot(acf(resid, lag = 40))
    
    plot of chunk tut8.3bbFitJAGS.AR1.5
    Conclusions - no obvious autocorrelation or other issues with residuals remaining
  • Explore parameter estimates
    tidyMCMC(as.mcmc(data.temporalCor.r2jags), conf.int = TRUE, conf.method = "HPDinterval")
    
          term    estimate   std.error    conf.low
    1  beta[1]  21.7994706 15.20123414  -6.1813173
    2  beta[2]   0.3673528  0.09013621   0.1837226
    3 deviance 387.0301283  2.71846642 383.0946706
    4      phi   0.8488231  0.06921159   0.7188739
    5    sigma  24.2744187 10.52484850  12.9694607
        conf.high
    1  51.6173567
    2   0.5305505
    3 392.1218115
    4   0.9769154
    5  40.1024073
    
  • Define the model
    						  functions { 
    						    matrix cov_matrix_ar1(real ar, real sigma, int nrows) { 
    						    matrix[nrows, nrows] mat; 
    						    vector[nrows - 1] gamma; 
    						    mat = diag_matrix(rep_vector(1, nrows)); 
    						    for (i in 2:nrows) { 
    						      gamma[i - 1] = pow(ar, i - 1); 
    						      for (j in 1:(i - 1)) { 
    						        mat[i, j] = gamma[i - j]; 
    						        mat[j, i] = gamma[i - j]; 
    						      } 
    						    } 
    						    return sigma^2 / (1 - ar^2) * mat; 
    						    }
    						  } 
    
    						  data { 
    						    int<lower=1> n;  // total number of observations 
    						    vector[n] y;  // response variable
    						    int<lower=1> nX;
    						    matrix[n,nX] X;
    						  } 
    						  transformed data {
    						  vector[n] se2 = rep_vector(0, n); 
    						  } 
    						  parameters { 
    						    vector[nX] beta;
    						    real<lower=0> sigma;  // residual SD 
    						    real <lower=-1,upper=1> phi;  // autoregressive effects 
    						  } 
    						  transformed parameters { 
    						  } 
    						  model {
    						    matrix[n, n] res_cov_matrix;
    						    matrix[n, n] Sigma; 
    						    vector[n] mu = X*beta;
    						    res_cov_matrix = cov_matrix_ar1(phi, sigma, n);
    						    Sigma = res_cov_matrix + diag_matrix(se2);
    						    Sigma = cholesky_decompose(Sigma); 
    
    						    // priors including all constants
    						    beta ~ student_t(3,30,30);
    						    sigma ~ cauchy(0,5);
    						    y ~ multi_normal_cholesky(mu,Sigma);
    						  } 
    						  generated quantities { 
    						  }
    
  • Arrange the data
    Xmat = model.matrix(~x, data.temporalCor)
    data.temporalCor.list <- with(data.temporalCor, list(y = y, X = Xmat,
        n = nrow(data.temporalCor), nX = ncol(Xmat)))
    
  • Fit the model
    library(rstan)
    
    data.temporalCor.rstan <- stan(data = data.temporalCor.list,
        model_code = stanString, chains = 3, iter = 2000, warmup = 500,
        thin = 3)
    
    In file included from /usr/local/lib/R/site-library/BH/include/boost/config.hpp:39:0,
                     from /usr/local/lib/R/site-library/BH/include/boost/math/tools/config.hpp:13,
                     from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core/var.hpp:7,
                     from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
                     from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core.hpp:12,
                     from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/mat.hpp:4,
                     from /usr/local/lib/R/site-library/StanHeaders/include/stan/math.hpp:4,
                     from /usr/local/lib/R/site-library/StanHeaders/include/src/stan/model/model_header.hpp:4,
                     from file7fdca3a04fa.cpp:8:
    /usr/local/lib/R/site-library/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
     #  define BOOST_NO_CXX11_RVALUE_REFERENCES
     ^
    <command-line>:0:0: note: this is the location of the previous definition
    
    SAMPLING FOR MODEL 'b9378eb4f1c2ff31d7e84cefd3f15961' NOW (CHAIN 1).
    
    Gradient evaluation took 0.000364 seconds
    1000 transitions using 10 leapfrog steps per transition would take 3.64 seconds.
    Adjust your expectations accordingly!
    
    
    Iteration:    1 / 2000 [  0%]  (Warmup)
    Iteration:  200 / 2000 [ 10%]  (Warmup)
    Iteration:  400 / 2000 [ 20%]  (Warmup)
    Iteration:  501 / 2000 [ 25%]  (Sampling)
    Iteration:  700 / 2000 [ 35%]  (Sampling)
    Iteration:  900 / 2000 [ 45%]  (Sampling)
    Iteration: 1100 / 2000 [ 55%]  (Sampling)
    Iteration: 1300 / 2000 [ 65%]  (Sampling)
    Iteration: 1500 / 2000 [ 75%]  (Sampling)
    Iteration: 1700 / 2000 [ 85%]  (Sampling)
    Iteration: 1900 / 2000 [ 95%]  (Sampling)
    Iteration: 2000 / 2000 [100%]  (Sampling)
    
     Elapsed Time: 2.15712 seconds (Warm-up)
                   2.92353 seconds (Sampling)
                   5.08065 seconds (Total)
    
    
    SAMPLING FOR MODEL 'b9378eb4f1c2ff31d7e84cefd3f15961' NOW (CHAIN 2).
    
    Gradient evaluation took 0.000218 seconds
    1000 transitions using 10 leapfrog steps per transition would take 2.18 seconds.
    Adjust your expectations accordingly!
    
    
    Iteration:    1 / 2000 [  0%]  (Warmup)
    Iteration:  200 / 2000 [ 10%]  (Warmup)
    Iteration:  400 / 2000 [ 20%]  (Warmup)
    Iteration:  501 / 2000 [ 25%]  (Sampling)
    Iteration:  700 / 2000 [ 35%]  (Sampling)
    Iteration:  900 / 2000 [ 45%]  (Sampling)
    Iteration: 1100 / 2000 [ 55%]  (Sampling)
    Iteration: 1300 / 2000 [ 65%]  (Sampling)
    Iteration: 1500 / 2000 [ 75%]  (Sampling)
    Iteration: 1700 / 2000 [ 85%]  (Sampling)
    Iteration: 1900 / 2000 [ 95%]  (Sampling)
    Iteration: 2000 / 2000 [100%]  (Sampling)
    
     Elapsed Time: 2.83276 seconds (Warm-up)
                   2.55861 seconds (Sampling)
                   5.39137 seconds (Total)
    
    
    SAMPLING FOR MODEL 'b9378eb4f1c2ff31d7e84cefd3f15961' NOW (CHAIN 3).
    
    Gradient evaluation took 0.000219 seconds
    1000 transitions using 10 leapfrog steps per transition would take 2.19 seconds.
    Adjust your expectations accordingly!
    
    
    Iteration:    1 / 2000 [  0%]  (Warmup)
    Iteration:  200 / 2000 [ 10%]  (Warmup)
    Iteration:  400 / 2000 [ 20%]  (Warmup)
    Iteration:  501 / 2000 [ 25%]  (Sampling)
    Iteration:  700 / 2000 [ 35%]  (Sampling)
    Iteration:  900 / 2000 [ 45%]  (Sampling)
    Iteration: 1100 / 2000 [ 55%]  (Sampling)
    Iteration: 1300 / 2000 [ 65%]  (Sampling)
    Iteration: 1500 / 2000 [ 75%]  (Sampling)
    Iteration: 1700 / 2000 [ 85%]  (Sampling)
    Iteration: 1900 / 2000 [ 95%]  (Sampling)
    Iteration: 2000 / 2000 [100%]  (Sampling)
    
     Elapsed Time: 2.64989 seconds (Warm-up)
                   2.32676 seconds (Sampling)
                   4.97665 seconds (Total)
    
  • Explore MCMC mixing diagnostics
    library(coda)
    stan_trace(data.temporalCor.rstan)
    
    plot of chunk tut8.3bbFitRSTAN.AR1.4
    stan_ac(data.temporalCor.rstan)
    
    plot of chunk tut8.3bbFitRSTAN.AR1.4
    stan_rhat(data.temporalCor.rstan)
    
    plot of chunk tut8.3bbFitRSTAN.AR1.4
    stan_ess(data.temporalCor.rstan)
    
    plot of chunk tut8.3bbFitRSTAN.AR1.4
    Conclusions - all diagnostics seem fine.
  • Perform model validation.

    Whenever we fit a model that incorporates changes to the variance-covariance structures, we need to explore modified standardized residuals. In this case, the raw residuals should be updated to reflect the autocorrelation (subtract residual from previous time weighted by the autocorrelation parameter) before standardizing by sigma \begin{align} Res_i =& Y_{i} - \mu_{i}\\ Res_{i+1} =& Res_{i+1} - \rho Res_{i}\\ Res_i =& \frac{Res_i}{\sigma}\\ \end{align}

    mcmc = as.matrix(data.temporalCor.rstan)
    wch = grep("beta", colnames(mcmc))
    # generate a model matrix
    newdata = data.frame(x = data.temporalCor$x)
    Xmat = model.matrix(~x, newdata)
    ## get median parameter estimates
    coefs = mcmc[, wch]
    fit = coefs %*% t(Xmat)
    resid = -1 * sweep(fit, 2, data.temporalCor$y, "-")
    n = ncol(resid)
    resid[, -1] = resid[, -1] - (resid[, -n] * mcmc[, "phi"])
    resid = apply(resid, 2, median)/median(mcmc[, "sigma"])
    fit = apply(fit, 2, median)
    ggplot() + geom_point(data = NULL, aes(y = resid, x = fit))
    
    plot of chunk tut8.3bbFitRSTAN.AR1.5
    ggplot() + geom_point(data = NULL, aes(y = resid, x = data.temporalCor$x))
    
    plot of chunk tut8.3bbFitRSTAN.AR1.5
    ggplot(data = NULL, aes(y = resid, x = data.temporalCor$year)) +
        geom_point() + geom_line() + geom_hline(yintercept = 0, linetype = "dashed")
    
    plot of chunk tut8.3bbFitRSTAN.AR1.5
    plot(acf(resid, lag = 40))
    
    plot of chunk tut8.3bbFitRSTAN.AR1.5
    fit = coefs %*% t(Xmat)
    ## draw samples from this model
    yRep = sapply(1:nrow(mcmc), function(i) rnorm(nrow(data.temporalCor),
        fit[i, ], mcmc[i, "sigma"]))
    ggplot() + geom_density(data = NULL, aes(x = as.vector(yRep),
        fill = "Model"), alpha = 0.5) + geom_density(data = data.temporalCor,
        aes(x = y, fill = "Obs"), alpha = 0.5)
    
    plot of chunk tut8.3bbFitRSTAN.AR1.5
    Conclusions - no obvious autocorrelation or other issues with residuals remaining
  • Explore parameter estimates
    tidyMCMC(data.temporalCor.rstan, par = c("beta", "phi", "sigma"),
        conf.int = TRUE, conf.method = "HPDinterval", rhat = TRUE,
        ess = TRUE)
    
         term   estimate   std.error   conf.low
    1 beta[1] 23.2545995 12.95018459 -0.5973877
    2 beta[2]  0.3671429  0.08696398  0.1963584
    3     phi  0.8626731  0.06820634  0.7318934
    4   sigma 11.4014846  1.18633413  9.2937517
       conf.high      rhat  ess
    1 51.0486563 1.0013252 1034
    2  0.5295868 1.0001775 1500
    3  0.9868331 1.0054995 1237
    4 13.7749828 0.9989293 1500
    
  • Fit the model. Note, I have increased the thinnning rate.
    library(brms)
    
    data.temporalCor.brms = brm(y ~ x, data = data.temporalCor,
        autocor = cor_ar(form = ~year, p = 1, cov = TRUE),
        iter = 2000, warmup = 200, chains = 3, thin = 10,
        refresh = 0, prior = c(prior(normal(0, 100), class = "Intercept"),
            prior(normal(0, 10), class = "b"), prior(cauchy(0,
                5), class = "sigma")))
    
    Gradient evaluation took 0.000423 seconds
    1000 transitions using 10 leapfrog steps per transition would take 4.23 seconds.
    Adjust your expectations accordingly!
    
    
    
     Elapsed Time: 2.61491 seconds (Warm-up)
                   4.80638 seconds (Sampling)
                   7.42129 seconds (Total)
    
    
    Gradient evaluation took 0.000311 seconds
    1000 transitions using 10 leapfrog steps per transition would take 3.11 seconds.
    Adjust your expectations accordingly!
    
    
    
     Elapsed Time: 1.75944 seconds (Warm-up)
                   4.23725 seconds (Sampling)
                   5.99669 seconds (Total)
    
    
    Gradient evaluation took 0.000261 seconds
    1000 transitions using 10 leapfrog steps per transition would take 2.61 seconds.
    Adjust your expectations accordingly!
    
    
    
     Elapsed Time: 2.13783 seconds (Warm-up)
                   4.77644 seconds (Sampling)
                   6.91427 seconds (Total)
    
  • Explore MCMC mixing diagnostics
    library(coda)
    stan_trace(data.temporalCor.brms$fit)
    
    plot of chunk tut8.3bbFitBRMS.AR1.4
    stan_ac(data.temporalCor.brms$fit)
    
    plot of chunk tut8.3bbFitBRMS.AR1.4
    stan_rhat(data.temporalCor.brms$fit)
    
    plot of chunk tut8.3bbFitBRMS.AR1.4
    stan_ess(data.temporalCor.brms$fit)
    
    plot of chunk tut8.3bbFitBRMS.AR1.4
    Conclusions - all diagnostics seem fine.
  • Perform model validation.

    Whenever we fit a model that incorporates changes to the variance-covariance structures, we need to explore modified standardized residuals. In this case, the raw residuals should be updated to reflect the autocorrelation (subtract residual from previous time weighted by the autocorrelation parameter) before standardizing by sigma \begin{align} Res_i =& Y_{i} - \mu_{i}\\ Res_{i+1} =& Res_{i+1} - \rho Res_{i}\\ Res_i =& \frac{Res_i}{\sigma}\\ \end{align}

    resid = resid(data.temporalCor.brms, type = "pearson")[, "Estimate"]
    n = length(resid)
    resid[-1] = resid[-1] - (resid[-n] * mean(as.matrix(data.temporalCor.brms)[,
        "ar[1]"]))
    fit = fitted(data.temporalCor.brms)[, "Estimate"]
    ggplot() + geom_point(data = NULL, aes(y = resid, x = fit))
    
    plot of chunk tut8.3bbFitBRMS.AR1.5
    ggplot() + geom_point(data = NULL, aes(y = resid, x = data.temporalCor$x))
    
    plot of chunk tut8.3bbFitBRMS.AR1.5
    ggplot(data = NULL, aes(y = resid, x = data.temporalCor$year)) +
        geom_point() + geom_line() + geom_hline(yintercept = 0, linetype = "dashed")
    
    plot of chunk tut8.3bbFitBRMS.AR1.5
    plot(acf(resid, lag = 40))
    
    plot of chunk tut8.3bbFitBRMS.AR1.5
    y_pred = posterior_predict(data.temporalCor.brms)
    newdata = data.temporalCor %>% cbind(t(y_pred)) %>% gather(key = "Rep",
        value = "Value", -y, -x, -year)
    ggplot(newdata, aes(Value, x = x)) + geom_violin(color = "blue",
        fill = "blue", alpha = 0.5) + geom_violin(data = data.temporalCor,
        aes(y = y, x = x), fill = "red", color = "red", alpha = 0.5)
    
    plot of chunk tut8.3bbFitBRMS.AR1.5
    Conclusions - no obvious autocorrelation or other issues with residuals remaining
  • Explore parameter estimates
    tidyMCMC(data.temporalCor.brms$fit, par = c("b_Intercept", "b_x",
        "ar", "sigma"), conf.int = TRUE, conf.method = "HPDinterval",
        rhat = TRUE, ess = TRUE)
    
             term   estimate   std.error    conf.low
    1 b_Intercept 21.6943704 21.71927159 -23.5374801
    2         b_x  0.3693860  0.08928101   0.1839314
    3       ar[1]  0.8687691  0.07187332   0.7461086
    4       sigma 11.3925161  1.16831533   9.4208507
       conf.high      rhat ess
    1 76.4179267 0.9989062 468
    2  0.5433781 1.0017556 540
    3  0.9973380 0.9957717 540
    4 13.8989327 0.9965453 540