AR1 in glmmTMB vs nlme

[mattansb]

I'm trying to fix a growth model such that:

1. Every person (group) has their own intercept and slope for time
2. There is an autocorrelation between the residuals

library(glmmTMB)
library(nlme)

# Data ------------------------------------

set.seed(1234)
n <- 25                                              ## Number of time points
x <- MASS::mvrnorm(mu = rep(0,n),
                   Sigma = .7 ^ as.matrix(dist(1:n)) )    ## Simulate the process using the MASS package
y <- x + rnorm(n)      
times <- 1:n
group <- factor(rep(1,n))


dat0 <- data.frame(y, times, group)

in nlme::lme() I would do it like this:

# With nlme ---------------------------------------------------------------

ma1 <- lme(y ~ times, data = dat0,
           random = ~ times | group)

mb1 <- lme(y ~ times, data = dat0,
           random = ~ times | group,
           correlation = corAR1(form = ~ times | group))

AIC(ma1, mb1) # df increased by 1
#>     df      AIC
#> ma1  6 78.62858
#> mb1  7 80.10863


# AR1 parameter:
mb1[["modelStruct"]][["corStruct"]]
#> Correlation structure of class corAR1 representing
#>       Phi 
#> 0.1781138

I expect to do this in glmmTMB like the following, but I run into the following issues:

1. convergence problem
2. the number of parameters has increased by 2 and not by 1
3. The ar parameter is difference from nlme.

# glmmTMB -----------------------------------------------------------------

ma2 <- glmmTMB(y ~ times + (1 + times | group), 
               data = dat0)

mb2 <- glmmTMB(y ~ times + (1 + times | group) + 
                 ar1(factor(times) + 0 | group), 
               data = dat0)
#> Warning message:
#> In finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old) :
#>   Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')


AIC(ma2, mb2) # df increased by 2 - Why?
#>     df      AIC
#> ma2  6 71.04545
#> mb2  8       NA

VarCorr(mb2)
#> 
#> Conditional model:
#>  Groups   Name           Std.Dev.   Corr                   
#>  group    (Intercept)    9.1176e-06                        
#>           times          6.7320e-07 0.976                  
#>  group.1  factor(times)1 3.8456e-01 0.616 (ar1) 0.616 (ar1)
#>  Residual                6.9058e-01  

I also tried without level 1 variance (dispformula), but still get different results.

mb3 <- glmmTMB(y ~ times + (1 + times | group) + 
                 ar1(factor(times) + 0 | group), 
               dispformula = ~0,
               data = dat0)
#> Warning message:
#> In finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old) :
#>   Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')


VarCorr(mb3)
#> 
#> Conditional model:
#>  Groups  Name           Std.Dev.   Corr                   
#>  group   (Intercept)    1.2204e-05                        
#>          times          3.1057e-09 0.996                  
#>  group.1 factor(times)1 7.8840e-01 0.060 (ar1) 0.060 (ar1)

What am I doing wrong?

[bbolker]

For the convergence problem, this is just glmmTMB being pickier than lme. lme is notorious (in my experience) for hiding poorly estimated random effects components, e.g.:

intervals(ma1)

Error in intervals.lme(ma1) :
cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance
Consider 'which = "fixed"'

As for the df increasing by 2: you need to set dispformula = ~0, otherwise you are estimating both an AR1 and a 'nugget' variance. You also need REML=TRUE to match the lme default.

mc2 <- update(mb2, REML = TRUE, dispformula = ~0)

This gives the same AR1 correlation estimate (print(VarCorr(mc2), digits = 8)). The glmmTMB object also gets the same (restricted) log-likelihood, if you know how to ask for it: mc2$obj$fn().

[mattansb]

Okay, so REML = TRUE, dispformula = ~0 gives the same results as lme() - thanks!

But you're also saying there's an issue here that lme() isn't flagging?
Is this related to the dataset, of the model? Because from what I can see random-growth model with autocorrelation is quite common...

[bbolker]

In this example you're trying to fit a random slopes model (the random = ~ times | group component, or + (times|group) in glmmTMB) to a data set with a single group. I don't see how that's ever going to work sensibly ... (I'm guessing this is just a "thinko" on your part?)

[mattansb]

Thanks @bbolker, that was very silly of me 🤦‍♂️.

Here is a more complete example with the sleepstudy dataset that also gives a warning and differs from lme():

library(glmmTMB)
#> Warning: package 'glmmTMB' was built under R version 4.3.2
library(nlme)

# Data ------------------------------------

data(sleepstudy, package = "lme4")

mb1 <- lme(Reaction ~ Days, 
           random = ~ Days | Subject,
           correlation = corAR1(form = ~ Days | Subject),
           data = sleepstudy,
           method = "REML", 
           control = lmeControl(opt = "optim"))

mb1[["modelStruct"]][["corStruct"]]
#> Correlation structure of class corAR1 representing
#>       Phi 
#> 0.4870368
VarCorr(mb1)
#> Subject = pdLogChol(Days) 
#>             Variance  StdDev    Corr  
#> (Intercept) 221.39047 14.879196 (Intr)
#> Days         22.65628  4.759861 0.897 
#> Residual    930.29232 30.500694


mb2 <- glmmTMB(Reaction ~ Days + (1 + Days | Subject) + 
                 ar1(factor(Days) + 0 | Subject),
               dispformula = ~0,
               REML = TRUE,
               data = sleepstudy,
               control = glmmTMBControl(optimizer = optim))
#> Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
#> problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')

VarCorr(mb2)
#> 
#> Conditional model:
#>  Groups    Name          Std.Dev.  Corr                   
#>  Subject   (Intercept)    2.072765                        
#>            Days           0.090318 -0.927                 
#>  Subject.1 factor(Days)0 49.355403 0.799 (ar1) 0.799 (ar1)

^{Created on 2024-02-03 with reprex v2.0.2}

The fixed effects are virtually the same, but the SE for the intercept is substatially different:

parameters::compare_parameters(mb1, mb2, select = "{estimate}, ({se})")
#> Parameter    |            mb1 |             mb2
#> -----------------------------------------------
#> (Intercept)  | 252.24, (6.85) | 253.73, (11.22)
#> Days         |  10.47, (1.53) |  10.47, ( 1.69)
#> -----------------------------------------------
#> Observations |            180 |             180

(I am trying to "convert" some users from SPSS to glmmTMB, but can only find SPSS -> nlme guides...)

[bbolker]

I don't know exactly. I did a little bit of digging, with the following conclusions:

• this is still a very unstable fit (mainly, I think, because a random effect for variation in slopes across subjects and an AR1 correlation term within subjects with a strong positive correlation are measuring approximately the same thing). intervals(mb1) reports that the 95% CI on the correlation between slopes and intercepts goes from -0.998 to 0.999 ... (the interval for the SD of the intercept is also large, 5-40, although that's not as obviously pathological ...) (mb1 reports the point estimate as 0.89, mb2 reports -0.93 ...) (I initially called it "singular", but I no longer think that's quite right)
• in this case mb2 does come up with a much worse overall fit mb2$obj$fn() is 869.6, -1*logLik(mb1) is 862.1 ... it would be interesting to investigate the likelihood surface (slightly tedious with 5 parameters, but not impossible), or to see how glmmTMB does from multiple starting values/if started from the lme estimates (converting lme parameter estimates into the glmmTMB parameterization is, again, slightly tedious but certainly do-able ...)
• the broom.mixed::tidy method for glmmTMB objects with AR1 components needs some love ...
• The covariance matrix is computed quite differently for glmmTMB and nlme. I initially thought it might have to do with which parameters are conditioned on (e.g., is the uncertainty due to uncertainty in the RE parameters included in the computation?), but a little bit of experimentation makes me less sure about that ...

More thoughts

• What would SPSS estimate/report for this example?
• It would be interesting to know which answer is closer to correct, i.e. in simulations of this type, what are the coverages of confidence intervals computed by nlme vs glmmTMB?
• This example illustrates some of the potential issues with singular and/or "non-positive-definite Hessian" fits. In a perfect world there's nothing wrong with a singular fit (it represents a valid optimization of the ML/REML criterion, it just happens to lie on a boundary of the feasible space), but ...
  • it may suggest that the fit is likely to be numerically unstable
  • it may suggest that there's lots of uncertainty in the RE parameters, in turn suggesting that conditioning on those parameter estimates when computing other quantities (such as the covariance matrix of the fixed parameters) might be dicey
  • more generally, it may suggest that (RE)MLE-based approaches, which do more conditioning and less averaging across uncertainty, may be less reliable
  • computations that depend on Wald approximations (i.e. relying on the curvature of the REML/likelihood surface at the estimated parameters) may break in interesting ways

[mattansb]

This is very interesting... Thanks for taking the time to look into / think about this.

I will run the analysis in SPSS and report back.

[mattansb]

While I work on the SPSS results, just as a point of reference, it seems like {brms} gives results that match nlme::lme() more closely than {glmmTMB} (though the HDI are quite wide):

library(brms)

# Data ------------------------------------

data(sleepstudy, package = "lme4")

mb3 <- brm(Reaction ~ Days + (Days | Subject) + ar(time = Days, gr = Subject), 
           data = sleepstudy, 
           backend = "cmdstanr")

variables(mb3)

pars <- as.matrix(mb3, variable = c("b_Intercept", "b_Days", 
                                    "sd_Subject__Intercept", "sd_Subject__Days", "cor_Subject__Intercept__Days",
                                    "sigma", "ar[1]")) |> 
  posterior::rvar(nchains = 4)

bayestestR::describe_posterior(pars, centrality = "map", ci_method = "hdi", test = NULL)
#> Summary of Posterior Distribution
#> 
#> Parameter                       |    MAP |           95% CI
#> -----------------------------------------------------------
#> x[b_Intercept]                  | 254.22 | [239.45, 268.61]
#> x[b_Days]                       |  10.24 | [  7.13,  13.98]
#> x[sd_Subject__Intercept]        |  18.96 | [  1.04,  34.70]
#> x[sd_Subject__Days]             |   5.07 | [  1.51,   9.52]
#> x[cor_Subject__Intercept__Days] |   0.82 | [ -0.32,   1.00]
#> x[sigma]                        |  27.16 | [ 24.00,  31.09]
#> x[ar[1]]                        |   0.52 | [  0.29,   0.81]

[bbolker]

Any further comments on this?

[mattansb]

I'm having a hard time tracking down someone who is proficient enough in SPSS :/
Still working on it...

[mattansb]

Alright, we go it! (SPSS syntax and output: MIXED.docx)

We can now compare:

Paramaeter	nlme	glmmTMB	brms	SPSS
Fixed
Intercept	252.24	253.73	254.22	252.14
Days	10.46	10.47	10.24	10.46
Random
Var(Intercept)	221.39	4.29	359.48	212.91
Var(Days)	22.65	0.008	25.70	22.06
Rho(Intercept, Days)	0.89	-0.92	0.82	0.99
Residual structure
Sigma^2	930.29	2435.95	737.66	875.89
Rho AR1	0.48	0.79	0.52	0.45

[bbolker]

OK, glmmTMB is clearly the odd one out (but also, pretty clear that this is a weakly identified problem. Not surprising that brms would be somewhat different (priors, marginal posterior medians (means?) rather than (RE)MLEs, etc.), but even var(intercept) and rho(intercept, days) are quite different between nlme and SPSS (maybe? REML vs ML?). My first guess/exploration will be considering different starting values ...

[mattansb]

SPSS is using REML.

Yes, this is a very "tight" model, but it is used in the growth models literature (might work better with more time points?).

Happy to help in any way that I can.

[strengejacke]

SPSS has a lot of quirks that aren't immediately clearly visible. E.g., SPSS uses PQL approximation in GLMMs by default, therefore, coefficients differ from what you get from R packages. There are too many hidden default settings in SPSS that make it almost useless to try to get exact results between R and SPSS (at least for everything beyond simple linear models).