ar_stuff.R
John Flournoy
2026-10-12
library(simsem)
library(nlme)
library(lmerTest)
library(lme4)
library(semPlot)
library(dplyr)
library(tidyr)
reps_per_n = 1e4
rerun = FALSE
N_seq = c(seq(30, 150, 40), seq(300, 1050, 150))
cores = 10
# install.packages(c('simsem', 'nlme', 'lmerTest', 'lme4', 'semPlot', 'dplyr', 'tidyr'))We want to sample data from a null model where x and y are correlated over time (due to some shared influences), but are not causally related. This is the null model we’re implicitly testing against, so we need to know that our procedure’s error rate performance is as expected (that is, that if we set \(\alpha=.05\) we’re actually controlling or error rate at 5%). I’ve created a model that has x and y observed at 12 occasions, with x and y having high stability (.7) and being somewhat correlated (.5).
someData <- simulateData(model=generatingModel, sample.nobs=250, empirical=T,
auto.fix.first = FALSE, auto.var = FALSE, auto.fix.single = FALSE,
auto.cov.lv.x = FALSE, auto.cov.y = FALSE)
fit.DGM <- lavaan(make_clpm_lavaan(generate = F), someData, fixed.x=F)
semPaths(fit.DGM, what='est', rotation=2, exoCov=T, exoVar=T, layout = 'tree2',
residuals = F, structural = F, sizeLat = 3, sizeMan = 3, sizeMan2 = 2,
intercepts = F, layoutSplit = F, bifactor = c('omega', 'kappa', paste0('q', 1:12)))
# summary(fit.DGM, standardize=T)
someData_long <- someData %>%
mutate(id = 1:n()) %>%
gather(key, value, -id) %>%
extract(key, c('var', 'wave'), '(x|y)(\\d+)') %>%
spread(var, value) %>%
mutate(wave = as.numeric(wave)) %>%
arrange(id, wave)
someData_long$gmean_x <- mean(someData_long$x)
someData_long$gmean_y <- mean(someData_long$y)
someData_long <- someData_long %>%
group_by(id) %>%
mutate(x_lag = lag(x),
y_lag = lag(y),
id_mean_x = mean(x),
wcen_x = x - id_mean_x,
gcen_x = id_mean_x - gmean_x,
wcen_x_lag = lag(wcen_x),
id_mean_y = mean(y),
wcen_y = y - id_mean_y,
gcen_y = id_mean_y - gmean_y,
wcen_y_lag = lag(wcen_y))Now we can see what happens if we predict \(y_{t}\) from \(x_{t-1}\) without accounting for previous levels of y.
lag1.misspec <- lme(y ~ 1 + x_lag, random = ~ 1 | id, data = someData_long,
na.action = na.omit)
summary(lag1.misspec)## Linear mixed-effects model fit by REML
## Data: someData_long
## AIC BIC logLik
## 8727.472 8751.147 -4359.736
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 0.6900746 1.092353
##
## Fixed effects: y ~ 1 + x_lag
## Value Std.Error DF t-value p-value
## (Intercept) 0.0000000 0.04836027 2499 0.00000 1
## x_lag 0.2798822 0.01872749 2499 14.94499 0
## Correlation:
## (Intr)
## x_lag 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.367431244 -0.637850873 0.006606042 0.633221867 3.104646784
##
## Number of Observations: 2750
## Number of Groups: 250We see a very robust effect of lagged x on y, even though no causal effect exists, by construction. This is due to the fact the x_lag is correlated with y_lag, and y_lag is correlated with y.
lag1.ylag <- lme(y ~ 1 + y_lag + x_lag, random = ~ 1 | id, data = someData_long,
na.action = na.omit)
summary(lag1.ylag)## Linear mixed-effects model fit by REML
## Data: someData_long
## AIC BIC logLik
## 7833.269 7862.86 -3911.634
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 0.03921928 0.9997851
##
## Fixed effects: y ~ 1 + y_lag + x_lag
## Value Std.Error DF t-value p-value
## (Intercept) 0.0000000 0.01922583 2498 0.00000 1
## y_lag 0.6980061 0.01650849 2498 42.28165 0
## x_lag 0.0000000 0.01650849 2498 0.00000 1
## Correlation:
## (Intr) y_lag
## y_lag 0.0
## x_lag 0.0 -0.5
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.19824100 -0.67240624 -0.01534057 0.64763590 3.48659284
##
## Number of Observations: 2750
## Number of Groups: 250Once we control for lagged y, which is the true causal effect, we see that lagged x does not have an influence on y
Can we use the residual autocorrelation structure in place of y_lag?
lag1.yAR <-lme(y ~ 1 + x_lag, random = ~ 1 | id, data = someData_long,
correlation = corAR1(form = ~ wave | id),
na.action = na.omit)
summary(lag1.yAR) ## Linear mixed-effects model fit by REML
## Data: someData_long
## AIC BIC logLik
## 7929.139 7958.732 -3959.57
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 0.1451315 1.343548
##
## Correlation Structure: AR(1)
## Formula: ~wave | id
## Parameter estimate(s):
## Phi
## 0.6758963
## Fixed effects: y ~ 1 + x_lag
## Value Std.Error DF t-value p-value
## (Intercept) 0.00000000 0.05045094 2499 0.000000 1.0000
## x_lag 0.02840258 0.01930589 2499 1.471188 0.1414
## Correlation:
## (Intr)
## x_lag 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.78583351 -0.66899662 0.00104045 0.64184993 3.11438368
##
## Number of Observations: 2750
## Number of Groups: 250Including the residual AR stucture helps, but we still see a positive estimate for the effect of x_lag. Since we’re working with data that reflects exaclty the data generating model, we should get exactly the same numbers that we put in – so this positive estimate is cause for concern. How much does this inflate our false positive rate if we sample randomly from the data generating model (as we do when we run a study)?
If we’re just estimating contemporaneous effects without allowing for possible lagged causes what do we get?
cont1.y <-lme(y ~ 1 + x, random = ~ 1 | id, data = someData_long,
na.action = na.omit)
summary(cont1.y) ## Linear mixed-effects model fit by REML
## Data: someData_long
## AIC BIC logLik
## 8825.798 8849.821 -4408.899
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 0.6377892 0.973903
##
## Fixed effects: y ~ 1 + x
## Value Std.Error DF t-value p-value
## (Intercept) 0.0 0.04408245 2749 0.00000 1
## x 0.5 0.01581666 2749 31.61223 0
## Correlation:
## (Intr)
## x 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.713824188 -0.651624345 -0.004876129 0.649875333 3.634580842
##
## Number of Observations: 3000
## Number of Groups: 250cont1.yAR <-lme(y ~ 1 + x, random = ~ 1 | id, data = someData_long,
correlation = corAR1(form = ~ wave | id),
na.action = na.omit)
summary(cont1.yAR) ## Linear mixed-effects model fit by REML
## Data: someData_long
## AIC BIC logLik
## 7710.68 7740.708 -3850.34
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 0.002047129 1.151257
##
## Correlation Structure: AR(1)
## Formula: ~wave | id
## Parameter estimate(s):
## Phi
## 0.6742302
## Fixed effects: y ~ 1 + x
## Value Std.Error DF t-value p-value
## (Intercept) 0.0 0.04108791 2749 0.00000 1
## x 0.5 0.01581666 2749 31.61223 0
## Correlation:
## (Intr)
## x 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.391260089 -0.666467292 0.004466756 0.658487307 3.519422128
##
## Number of Observations: 3000
## Number of Groups: 250This looks right, ’cause this is very close to the generating model with no possible confounds.
Numerous simulations
generate_data <- function(n = 250, gen_mod, wcen = FALSE){
someData <- simulateData(model=gen_mod, sample.nobs=n, empirical=F)
someData_long <- someData %>%
mutate(id = 1:n()) %>%
gather(key, value, -id) %>%
extract(key, c('var', 'wave'), '(x|y)(\\d+)') %>%
spread(var, value) %>%
mutate(wave = as.numeric(wave)) %>%
arrange(id, wave) %>%
group_by(id) %>%
mutate(x_lag = lag(x),
y_lag = lag(y))
if(wcen){
someData_long$gmean_x <- mean(someData_long$x)
someData_long$gmean_y <- mean(someData_long$y)
someData_long <- someData_long %>%
group_by(id) %>%
mutate(id_mean_x = mean(x),
wcen_x = x - id_mean_x,
gcen_x = id_mean_x - gmean_x,
wcen_x_lag = lag(wcen_x),
id_mean_y = mean(y),
wcen_y = y - id_mean_y,
gcen_y = id_mean_y - gmean_y,
wcen_y_lag = lag(wcen_y))
}
return(someData_long)
}
run_model <- function(n = 250, gen_mod, wcen = FALSE, use_lag_y = FALSE, use_lme4 = FALSE, wcen_re = FALSE){
someData_long <- generate_data(n = n, gen_mod = gen_mod, wcen = wcen)
if(wcen){
if(use_lag_y){
if(use_lme4){
if(wcen_re){
lag1.yAR <-lmer(y ~ 1 + wcen_y_lag + wcen_x_lag + gcen_x + (1 + wcen_x_lag || id), data = someData_long,
REML = TRUE,
na.action = na.omit)
} else {
lag1.yAR <-lmer(y ~ 1 + wcen_y_lag + wcen_x_lag + gcen_x + (1 | id), data = someData_long,
REML = TRUE,
na.action = na.omit)
}
} else {
lag1.yAR <-lme(y ~ 1 + wcen_y_lag + wcen_x_lag + gcen_x, random = ~ 1 | id, data = someData_long,
na.action = na.omit,
control=lmeControl(opt = "optim"))
}
} else {
lag1.yAR <-lme(y ~ 1 + wcen_x_lag + gcen_x, random = ~ 1 | id, data = someData_long,
correlation = corAR1(form = ~ wave | id),
na.action = na.omit,
control=lmeControl(opt = "optim"))
}
if(use_lme4){
asum <- coef(summary(lag1.yAR, ddf = 'Satterthwaite'))
x_lag_p <- asum['wcen_x_lag', 'Pr(>|t|)']
x_lag_bias <- asum['wcen_x_lag', 'Estimate']
} else {
asum <- summary(lag1.yAR)
x_lag_p <- asum$tTable['wcen_x_lag', 'p-value']
x_lag_bias <- asum$tTable['wcen_x_lag', 'Value']
}
} else {
if(use_lag_y){
lag1.yAR <-lme(y ~ 1 + y_lag + x_lag, random = ~ 1 | id, data = someData_long,
na.action = na.omit)
} else {
lag1.yAR <-lme(y ~ 1 + x_lag, random = ~ 1 | id, data = someData_long,
correlation = corAR1(form = ~ wave | id),
na.action = na.omit)
}
asum <- summary(lag1.yAR)
x_lag_p <- asum$tTable['x_lag', 'p-value']
x_lag_bias <- asum$tTable['x_lag', 'Value']
}
return(data.frame(n = n, x_lag_p = x_lag_p, x_lag_bias = x_lag_bias))
}
if(rerun){
library(parallel)
system.time(no_wcen_ar_reps <- mclapply(N_seq,
function(n) {
somereps <- dplyr::bind_rows(
replicate(reps_per_n,
run_model(n, generatingModel, wcen = F),
simplify = F))
return(somereps)
}, mc.cores = cores))
saveRDS(no_wcen_ar_reps, 'no_wcen_ar.RDS')
system.time(wcen_ar_reps <- mclapply(N_seq,
function(n) {
somereps <- dplyr::bind_rows(
replicate(reps_per_n,
run_model(n, generatingModel, wcen = T),
simplify = F))
return(somereps)
}, mc.cores = cores))
saveRDS(wcen_ar_reps, 'wcen_ar.RDS')
system.time(no_wcen_ylag_reps <- mclapply(N_seq,
function(n) {
somereps <- dplyr::bind_rows(
replicate(reps_per_n,
run_model(n, generatingModel, wcen = F, use_lag_y = T),
simplify = F))
return(somereps)
}, mc.cores = cores))
saveRDS(no_wcen_ylag_reps, 'no_wcen_ylag_reps.RDS')
system.time(wcen_ylag_reps <- mclapply(N_seq,
function(n) {
somereps <- dplyr::bind_rows(
replicate(reps_per_n,
run_model(n, generatingModel, wcen = T, use_lag_y = T),
simplify = F))
return(somereps)
}, mc.cores = cores))
saveRDS(wcen_ylag_reps, 'wcen_ylag_reps.RDS')
system.time(wcen_ylag_l4_reps <- mclapply(N_seq,
function(n) {
somereps <- dplyr::bind_rows(
replicate(reps_per_n,
run_model(n, generatingModel, wcen = T, use_lag_y = T,
use_lme4 = T, wcen_re = F),
simplify = F))
return(somereps)
}, mc.cores = cores))
saveRDS(wcen_ylag_l4_reps, 'wcen_ylag_l4_reps.RDS')
system.time(wcen_ylag_l4_wcenre_reps <- mclapply(N_seq,
function(n) {
somereps <- dplyr::bind_rows(
replicate(reps_per_n,
run_model(n, generatingModel, wcen = T, use_lag_y = T,
use_lme4 = T, wcen_re = T),
simplify = F))
return(somereps)
}, mc.cores = cores))
saveRDS(wcen_ylag_l4_wcenre_reps, 'wcen_ylag_l4_wcenre_reps.RDS')
} else {
no_wcen_ar_reps <- readRDS('no_wcen_ar.RDS')
wcen_ar_reps <- readRDS('wcen_ar.RDS')
no_wcen_ylag_reps <- readRDS('no_wcen_ylag_reps.RDS')
wcen_ylag_reps <- readRDS('wcen_ylag_reps.RDS')
wcen_ylag_l4_reps <- readRDS('wcen_ylag_l4_reps.RDS')
wcen_ylag_l4_wcenre_reps <- readRDS('wcen_ylag_l4_wcenre_reps.RDS')
}Using AR residual structure to account for stability
lag1.yAR <-lme(y ~ 1 + x_lag, random = ~ 1 | id, data = someData_long,
correlation = corAR1(form = ~ wave | id),
na.action = na.omit)plot_results(no_wcen_ar_reps, max_fp = 1)## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
What we see here is a consistent bias in the estimates of the effect of lagged x on y. As our sample size increases, we are better and better powered, which means we erroneously reject the null more often due to this bias. In other words, if we are powered to detect an effect of lagged x on y, we are powered to mistake this bias for a causal effect, even using the AR structure to adjust for the stability of y.
Using lagged DV to account for stability
lag1.yAR <-lme(y ~ 1 + y_lag + x_lag, random = ~ 1 | id, data = someData_long,
na.action = na.omit)plot_results(no_wcen_ylag_reps)
Here we see no bias in the estimate of the effect of lagged x on y, and correspondingly our error rate is controlled appropriately.
AR with within-person-centering
We might want to create separate IVs for the within-person versus the between-person effect. In this example, I’ve created a variable that is just each person’s mean across all waves, centered at the group mean (gcen) and a within-person variable that is each person’s score at a wave minus their mean.
This is in line with the advice in:
Wang, L. (Peggy), & Maxwell, S. E. (2015). On disaggregating between-person and within-person effects with longitudinal data using multilevel models. Psychological Methods, 20(1), 63–83. https://doi.org/10.1037/met0000030
lag1.yAR <-lme(y ~ 1 + wcen_x_lag + gcen_x, random = ~ 1 | id, data = someData_long,
correlation = corAR1(form = ~ wave | id),
na.action = na.omit)plot_results(wcen_ar_reps)
Interstingly, when we within-person center the lagged x variable, we see a negative bias. This is also picked up in the p-values as our sample size increases.
Lagged y with within-person-centering
Note that when we within-person center the IV, x, we also within-person center the lagged DV, y, and include that (the between-person component of y is redundant with the random intercept. If you don’t do this, wow, the false positive error rate goes way up. I’m not quite sure why this is (simulations not shown).
lag1.yAR <-lme(y ~ 1 + wcen_y_lag + wcen_x_lag + gcen_x, random = ~ 1 | id, data = someData_long,
na.action = na.omit) plot_results(wcen_ylag_reps)
The error rate in this case is still higher than the nominal rate, but since bias is 0, this may be due to a problem with the standard errors, or the DF being used. Perhaps another possibility is that the random effects are mispecified. I’m not sure why this is, but we can test whether using the Satterthwaite correction, or including a random effect of the within-person predictor helps.
Lagged y with within-person-centering, Satterthwaite correction
Perhaps using the Satterthwaite correction will appropriately adjust the p-values.
lag1.yAR <-lmer(y ~ 1 + wcen_y_lag + wcen_x_lag + gcen_x + (1 | id), data = someData_long,
REML = TRUE,
na.action = na.omit)
summary(lag1.yAR, ddf = 'Satterthwaite')plot_results(wcen_ylag_l4_reps)
Bummer. Still a higher than expected error rate.
Lagged y with within-person-centering, Satterthwaite correction + WCEN RE
The last thing I’ll try is to add the within-person centered lagged predictor as random effect.
lag1.yAR <-lmer(y ~ 1 + wcen_y_lag + wcen_x_lag + gcen_x + (1 + wcen_x_lag || id), data = someData_long,
REML = TRUE,
na.action = na.omit)
summary(lag1.yAR, ddf = 'Satterthwaite')plot_results(wcen_ylag_l4_wcenre_reps)
That looks like it helped a little bit, but it also seems to induce a small negative bias.
What’s going on?
The problem with the naive model is an ordinary omitted-variable problem, just one that is easy to miss because the omitted variable is the outcome itself, one wave earlier. In the generating model \(y_{t-1}\) causes \(y_{t}\) (stability of .7), and \(x_{t-1}\) is correlated with \(y_{t-1}\) (they share correlated residuals, at .5), so there is an open path from \(x_{t-1}\) to \(y_{t}\) that runs through \(y_{t-1}\) and has nothing to do with \(x\) causing \(y\). Starting from wave 1, where both variances are 1 and the correlation is .5, the covariance of \(y_{2}\) with \(x_{1}\) is .7 × .5 = .35, and a regression of \(y_{t}\) on \(x_{t-1}\) that does not condition on \(y_{t-1}\) reads that covariance as an effect. It is also worth noticing what the random intercept is doing in that model. There are no stable between-person differences in this generating model (it is a CLPM, with no trait component), and yet the naive model above estimates a random-intercept SD of .69 (against a residual SD of 1.09), which falls to .04 once \(y_{t-1}\) is in the model. Twelve waves of a process with stability .7 make people look as though they differ in their typical level, because someone who happens to start high stays high for a while; the random intercept absorbs some of that apparent stability, \(x_{t-1}\) picks up a good deal of the rest, and neither one is a substitute for conditioning on where \(y\) actually was last wave.
The AR(1) residual structure gets surprisingly close (in the single data set above it takes the estimate from .28 to .03), but it is a model of the errors rather than a predictor, and it does not quite close the path. My best guess at why (supported by a quick side simulation, not shown) is each person’s first observation. Fitting AR(1) errors amounts to quasi-differencing, i.e., regressing \(y_{t} - \phi y_{t-1}\) on \(x_{t-1} - \phi x_{t-2}\), and with \(\phi\) near .7 that removes the \(y_{t-1}\) pathway from every row except the first one for each person, which has nothing earlier to difference against and so enters as the naive regression of \(y_{2}\) on \(x_{1}\), confound and all. On top of that, \(\phi\) is estimated rather than known, and it comes in a little under .7 (.68 in the example above), which leaves a sliver of \(y_{t-1}\) in the other rows too. The important thing about this bias is that it is a property of each person’s series (one contaminated row in eleven), not of the sample size, so it does not shrink as \(N\) grows. That is what the simulations show: the average estimate sits at about .03 from \(N\) = 30 to \(N\) = 750, and the false-positive rate climbs from .11 to .72 over the same range (see the tables below), because the standard error keeps shrinking around an estimate that is not centered on zero.
That pattern turns out to be a fairly useful diagnostic for the other specifications. A constant bias produces a false-positive rate that rises with \(N\); a false-positive rate that is inflated but flat across \(N\), with an average estimate at zero, points instead to standard errors (or degrees of freedom) that are too optimistic. With 10,000 replications per cell, the Monte Carlo standard error of a false-positive rate near .05 is about .002, so differences of a percentage point or so are real. Reading the tables with that in mind:
- Including lagged \(y\) as a predictor, with no centering, is the only specification that holds the error rate at the nominal level: .049 to .055 for \(N\) from 30 to 600, with average estimates within .001 of zero.
- Within-person centering with AR(1) residuals swaps the positive bias for a small negative one (about −.009, steady across \(N\)), and the false-positive rate rises accordingly, from .06 at \(N\) = 30 to .155 at \(N\) = 900.
- Within-person centering with lagged \(y\) leaves the average estimate within .001 of zero, but the false-positive rate sits between about .074 and .083 at every \(N\), with or without the Satterthwaite correction. Because it does not move between \(N\) = 30 and \(N\) = 1050, where the degrees of freedom number in the thousands, a small-sample degrees-of-freedom problem seems unlikely to be the explanation, which is presumably why Satterthwaite does not help.
- Adding a random slope for the centered lagged \(x\) brings the false-positive rate down at small \(N\) (.054 at \(N\) = 30), but it introduces a small negative bias (about −.004), and the rate creeps up with \(N\) (to about .075 to .078 at the largest samples), which is the signature of bias rather than of bad standard errors.
I do not have a settled account of why person-mean centering misbehaves here, so what follows is speculative. The person means are computed from all 12 waves, including \(y_{t}\) and everything after it, so the centered lagged predictors are correlated with the current residual; this is the short-panel problem econometricians call Nickell bias, and with only 12 waves it would not be expected to wash out. Another, perhaps contributing, possibility is that the generating model is not stationary: every residual variance (including wave 1) is fixed at 1, so the variances of \(x\) and \(y\) grow from 1 toward about 1.96 (i.e., 1/(1 − .7²)) over the first several waves, and none of the fitted models allow the variance to change across waves. Neither of these has been tested here.
Overall, in these simulations the safest thing to do with a lagged (Granger-type) question is to put the prior level of the outcome in the model as a predictor, and to not count on the AR(1) residual structure to do that job, since it leaves a bias that more data only makes easier to mistake for an effect. Note well that the lagged-\(y\) model is also the true model for these data, so these results are the best case for it. In real data, where people plausibly do differ in their typical levels of \(x\) and \(y\), a lagged outcome alongside a random intercept has a well-known bias of its own (the lagged outcome is correlated with the random intercept), and the principled move is to separate the stable between-person part from the within-person dynamics explicitly, e.g., with a random-intercept cross-lagged panel model (RI-CLPM) or latent person-mean centering. Observed person-mean centering, as tried here, did not get the error rate under control.
Some limits on all of this:
- There is one generating model: a CLPM with 12 waves, a stability of .7 for both variables, a contemporaneous correlation of .5, and a cross-lagged effect of exactly zero. None of these were varied, so I cannot say how the bias scales with stability, correlation, or the number of waves (though the first-row account above predicts it should shrink with more waves).
- Because the generating model has no stable trait component, these
results say nothing directly about the case that motivates separating
within- from between-person variance in the first place.
make_riclpm_lavaan()is defined above but not used to generate any data; running the same specifications on RI-CLPM data is the obvious next simulation. - Only the null was simulated, so there is nothing here about power, or about how well any of these models recover a nonzero cross-lagged effect.
- Some sample sizes failed to finish for some specifications and are missing from the plots above and from the tables below (listed under the tables).
| n | y ~ x_lag, AR(1) | y ~ y_lag + x_lag | centered x_lag, AR(1) | centered y_lag + x_lag (nlme) | centered y_lag + x_lag (lme4, Satterthwaite) | centered y_lag + x_lag (lme4, Satterthwaite, random slope) |
|---|---|---|---|---|---|---|
| 30 | 0.115 | 0.052 | 0.059 | 0.082 | 0.075 | 0.054 |
| 70 | 0.165 | 0.055 | 0.062 | 0.077 | 0.083 | 0.061 |
| 110 | 0.200 | 0.051 | 0.067 | 0.076 | 0.081 | 0.063 |
| 150 | 0.245 | 0.051 | 0.071 | 0.076 | 0.080 | 0.069 |
| 300 | 0.394 | 0.050 | 0.089 | 0.078 | 0.074 | 0.063 |
| 450 | 0.527 | 0.049 | 0.102 | 0.079 | 0.081 | 0.069 |
| 600 | 0.637 | 0.049 | 0.119 | 0.078 | 0.079 | 0.069 |
| 750 | 0.721 | failed | 0.135 | 0.078 | 0.080 | 0.077 |
| 900 | failed | failed | 0.155 | 0.078 | 0.082 | 0.078 |
| 1050 | failed | failed | failed | 0.075 | 0.081 | 0.075 |
| n | y ~ x_lag, AR(1) | y ~ y_lag + x_lag | centered x_lag, AR(1) | centered y_lag + x_lag (nlme) | centered y_lag + x_lag (lme4, Satterthwaite) | centered y_lag + x_lag (lme4, Satterthwaite, random slope) |
|---|---|---|---|---|---|---|
| 30 | 0.0322 | 0.0008 | -0.0067 | -0.0007 | -0.0008 | -0.0038 |
| 70 | 0.0315 | -0.0000 | -0.0089 | -0.0012 | -0.0003 | -0.0046 |
| 110 | 0.0299 | 0.0000 | -0.0087 | -0.0008 | -0.0006 | -0.0039 |
| 150 | 0.0300 | 0.0006 | -0.0092 | -0.0008 | -0.0011 | -0.0050 |
| 300 | 0.0294 | -0.0001 | -0.0096 | -0.0011 | -0.0008 | -0.0041 |
| 450 | 0.0293 | -0.0001 | -0.0092 | -0.0008 | -0.0008 | -0.0039 |
| 600 | 0.0292 | -0.0001 | -0.0093 | -0.0007 | -0.0010 | -0.0041 |
| 750 | 0.0290 | failed | -0.0093 | -0.0008 | -0.0006 | -0.0043 |
| 900 | failed | failed | -0.0096 | -0.0007 | -0.0008 | -0.0041 |
| 1050 | failed | failed | failed | -0.0008 | -0.0008 | -0.0040 |
Cells where the simulation run failed (no results saved):
- y ~ x_lag, AR(1), N = 900
- y ~ x_lag, AR(1), N = 1050
- y ~ y_lag + x_lag, N = 750
- y ~ y_lag + x_lag, N = 900
- y ~ y_lag + x_lag, N = 1050
- centered x_lag, AR(1), N = 1050