Mixed model with autoregressive autocorrelation (first order) yield different estimates for phi between brms and glmmTMB/nlme

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


If TMB/nlme are doing max marginal fits like lme4, then you shouldn’t expect to get the same answer as in a full Bayesian posterior

To test your own code without using these packages, you can do is simulate data with known parameters then fit and see if you recover them in their posterior intervals at the nominal coverage. The formal version of this is called “simulation-based calibration”. If you can recover this way, your model’s doing the right thing.