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