Hey all,
I simulated longitudinal data with an AR correlation structure (with the help of the function of Ben Bolker). I expanded the simulation in a multi trial context. I set phi of the correlation matrix to 0.8. I fitted the model with glmmTMB and nlme and got 0.8 as estimate. But fitting the model with brms, I got an estimate of about 0.87. Did I something wrong? I do not understand, why I observed this difference between brms and the other packages. Thanks in advance for any help. In the following is my R code and the results:.
Thanks in advance for any suggestion.
#### Simulation of the data
library(tidyverse)
library(brms)
library(glmmTMB)
library(nlme)
library(rlist)
simCor1 <- function(phi = 0.8,
sdgrp = 2,
sdres = 1,
npergrp = 20,
ngrp = 20,
seed = NULL,
## set linkinv/simfun for GLMM sims
linkinv = identity,
simfun = identity)
{
if (!is.null(seed)) set.seed(seed)
cmat <- sdres*phi^abs(outer(0:(npergrp-1),0:(npergrp-1),"-"))
errs <- MASS::mvrnorm(ngrp,mu=rep(0,npergrp),Sigma=cmat)
ranef <- rnorm(ngrp,mean=0,sd=sdgrp)
d <- data.frame(f=rep(1:ngrp,each=npergrp))
eta <- ranef[as.numeric(d$f)] + c(t(errs)) ## unpack errors by row
mu <- linkinv(eta)
d$y <- simfun(mu)
d$tt <- factor(rep(1:npergrp,ngrp))
return(d)
}
d <- NULL
for(i in 1:10)
{
d_temp <- simCor1(phi = 0.8,
sdgrp = 2,
sdres = 1,
seed = (i+100))
d_temp$study <- paste0("study_",i)
d_temp$y <- d_temp$y + rnorm(1,1,0.5)
d_temp$treat <- rep(c(rep(1,10),
rep(0,10)),20)
d_temp[which(d_temp$treat == 1),"y"] <- d_temp[which(d_temp$treat == 1),"y"] + rnorm(1,1,0.1)
d[[i]] <- d_temp
rm(d_temp)
}
d <- list.rbind(d)
d$f_study <- paste0(d$study,"_",d$f)
##### Fitting of the models
m1 <- glmmTMB(y~1 + treat + (1|study) + (1|f_study) + ar1(tt - 1|f_study),
data=d,
family=gaussian)
m2 <- lme(y ~ treat,
random = list(study = ~ 1, f_study = ~ 1),
correlation = corAR1(),
data = d)
m3 <- brm(y ~ 1 + treat + (1|study) + (1|f_study) + ar(p = 1, gr = study:f_study),
data = d,
cores = 4,
iter = 2000)
###### Summary glmmTMB model:
summary(m1)
Family: gaussian ( identity )
Formula: y ~ 1 + treat + (1 | study) + (1 | f_study) + ar1(tt - 1 | f_study)
Data: d
AIC BIC logLik deviance df.resid
7909.8 7953.8 -3947.9 7895.8 3993
Random effects:
Conditional model:
Groups Name Variance Std.Dev. Corr
study (Intercept) 4.312e-01 0.6566454
f_study (Intercept) 4.763e+00 2.1823505
f_study.1 tt1 9.824e-01 0.9911571 0.80 (ar1)
Residual 1.971e-08 0.0001404
Number of obs: 4000, groups: study, 10; f_study, 200
Dispersion estimate for gaussian family (sigma^2): 1.97e-08
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.09107 0.26256 4.155 3.25e-05 ***
treat 1.04571 0.04021 26.004 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##### Summary nlme model:
summary(m2)
Linear mixed-effects model fit by REML
Data: d
AIC BIC logLik
7913.153 7950.915 -3950.577
Random effects:
Formula: ~1 | study
(Intercept)
StdDev: 0.7122897
Formula: ~1 | f_study %in% study
(Intercept) Residual
StdDev: 2.18227 0.9915965
Correlation Structure: AR(1)
Formula: ~1 | study/f_study
Parameter estimate(s):
Phi
0.8039477
Fixed effects: y ~ treat
Value Std.Error DF t-value p-value
(Intercept) 1.091068 0.27668685 3799 3.943332 1e-04
treat 1.045712 0.04021853 3799 26.000754 0e+00
Correlation:
(Intr)
treat -0.073
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.038103847 -0.541361092 0.005007161 0.562445221 2.690507449
Number of Observations: 4000
Number of Groups:
study f_study %in% study
10 200
##### Summary of brms model :
summary(m3)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: y ~ 1 + treat + (1 | study) + (1 | f_study) + ar(p = 1, gr = study:f_study)
Data: d (Number of observations: 4000)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Correlation Structures:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
ar[1] 0.87 0.01 0.85 0.90 1.00 1874 2254
Group-Level Effects:
~f_study (Number of levels: 200)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.27 0.12 2.05 2.53 1.01 541 798
~study (Number of levels: 10)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.80 0.34 0.15 1.54 1.02 142 71
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 1.10 0.31 0.48 1.73 1.00 953 1427
treat 1.05 0.04 0.97 1.12 1.00 3486 2712
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.61 0.01 0.59 0.62 1.00 3419 2688
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).