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)

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

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)

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))

ggplot() + geom_point(data = NULL, aes(y = resid, x = data.temporalCor$x))

ggplot(data = NULL, aes(y = resid, x = data.temporalCor$year)) +
    geom_point() + geom_line() + geom_hline(yintercept = 0, linetype = "dashed")

plot(acf(resid, lag = 40))

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)

stan_ac(data.temporalCor.rstan)

stan_rhat(data.temporalCor.rstan)

stan_ess(data.temporalCor.rstan)

Conclusions - all diagnostics seem fine.

Perform model validation.

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))

ggplot() + geom_point(data = NULL, aes(y = resid, x = data.temporalCor$x))

ggplot(data = NULL, aes(y = resid, x = data.temporalCor$year)) +
    geom_point() + geom_line() + geom_hline(yintercept = 0, linetype = "dashed")

plot(acf(resid, lag = 40))

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)

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)

stan_ac(data.temporalCor.brms$fit)

stan_rhat(data.temporalCor.brms$fit)

stan_ess(data.temporalCor.brms$fit)

Conclusions - all diagnostics seem fine.

Perform model validation.

resid = resid(data.temporalCor.brms, type = "pearson")[, "Estimate"]
fit = fitted(data.temporalCor.brms)[, "Estimate"]
ggplot() + geom_point(data = NULL, aes(y = resid, x = fit))

ggplot() + geom_point(data = NULL, aes(y = resid, x = data.temporalCor$x))

ggplot(data = NULL, aes(y = resid, x = data.temporalCor$year)) +
    geom_point() + geom_line() + geom_hline(yintercept = 0, linetype = "dashed")

plot(acf(resid, lag = 40))

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)

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

General
JAGS
STAN
BRMS

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)

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.

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))

ggplot() + geom_point(data = NULL, aes(y = resid, x = data.temporalCor$x))

ggplot(data = NULL, aes(y = resid, x = data.temporalCor$year)) +
    geom_point() + geom_line() + geom_hline(yintercept = 0, linetype = "dashed")

plot(acf(resid, lag = 40))

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)

stan_ac(data.temporalCor.rstan)

stan_rhat(data.temporalCor.rstan)

stan_ess(data.temporalCor.rstan)

Conclusions - all diagnostics seem fine.

Perform model validation.

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))

ggplot() + geom_point(data = NULL, aes(y = resid, x = data.temporalCor$x))

ggplot(data = NULL, aes(y = resid, x = data.temporalCor$year)) +
    geom_point() + geom_line() + geom_hline(yintercept = 0, linetype = "dashed")

plot(acf(resid, lag = 40))

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)

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)

stan_ac(data.temporalCor.brms$fit)

stan_rhat(data.temporalCor.brms$fit)

stan_ess(data.temporalCor.brms$fit)

Conclusions - all diagnostics seem fine.

Perform model validation.

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))

ggplot() + geom_point(data = NULL, aes(y = resid, x = data.temporalCor$x))

ggplot(data = NULL, aes(y = resid, x = data.temporalCor$year)) +
    geom_point() + geom_line() + geom_hline(yintercept = 0, linetype = "dashed")

plot(acf(resid, lag = 40))

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)

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