Standard error of the AR(1) correlation parameter in glmmTMB

Question

Using AR(1) model in glmmTMB, with code based on the vignette (btw, haven't found any other reference documentation, that would precisely define the ar1() functionality and the corresponding summary output etc.; is there something else than this vignette?):

library(glmmTMB)

n <- 25                                              ## Number of time points
x <- MASS::mvrnorm(mu = rep(0,n),
             Sigma = 0.7 ^ as.matrix(dist(1:n)) )    ## Simulate the process using the MASS package

times <- factor(1:n, levels=1:n)
head(levels(times))

group <- factor(rep(1,n))
dat0 <- data.frame(x, times, group)

g2 <- glmmTMB(x ~ ar1(times + 0 | group), data=dat0) 
# the AR(1) correlation should be ~ 0.7 for large enough samples

summary(g2)

# compare with the raw correlation coefficient:
cor.test(head(x, -1), tail(x, -1))

It seems to be fine, but I can't seem to figure out how do I get the standard error of the Corr parameter? (The AR(1) correlation)? I suppose all the optimized parameters should always have the standard error available... without any need to use some bootstrapping etc...

PS: I tried things like confint(g2) or confint(g2, parm = "theta_") but it seems to return something else...

PS: In the end, I rather need the CI. The standard error was just means to get it :-)

EDIT: Ben's current hack/solution doesn't work in general, when I add more random effects:

n <- 10                                              ## Number of time points
time.reff <- rnorm(rep(0,n), 0, 0.2)
x <- MASS::mvrnorm(1000, mu = time.reff,
             Sigma = .7 ^ as.matrix(dist(1:n)) )    ## Simulate the process using the MASS package

library(reshape2)

d <- melt(x)
d$time <- factor(d$Var2)
d$group <- d$Var1

g3 <- glmmTMB(value ~ (1|time) + ar1(time + 0 | group), data = d)

Now, when trying to extract it from the glmmTMB object, it all has just similar confusing names:

vcov(g3, full = TRUE)
#                          (Intercept) disp~(Intercept)  theta_1|time.1 theta_time+0|group.1 theta_time+0|group.2
# (Intercept)           0.0044675603195  0.0000006258999 -0.000037166774      0.0000001358308      0.0000003024984
# disp~(Intercept)      0.0000006258999  0.0763324452714  0.000014868716     -0.0014493749352      0.0063889358336
# theta_1|time.1       -0.0000371667736  0.0000148687157  0.051209915518      0.0000014854011      0.0000051939167
# theta_time+0|group.1  0.0000001358308 -0.0014493749352  0.000001485401      0.0001635499153      0.0000657916127
# theta_time+0|group.2  0.0000003024984  0.0063889358336  0.000005193917      0.0000657916127      0.0009709511260

Answer 1

I didn't see anything in my quick perusal of the glmmTMB help file. However, it would be easy enough to make something clean that would do that.

Data

library(glmmTMB)
set.seed(519)
  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

times <- factor(1:n, levels=1:n)
head(levels(times))
#> [1] "1" "2" "3" "4" "5" "6"

group <- factor(rep(1,n))
dat0 <- data.frame(x, times, group)

Model

g2 <- glmmTMB(x ~ ar1(times + 0 | group), data=dat0)

Extract AR parameter

Here's a function called get_AR() that will pull the AR parameter out in its theta form. It creates a normal theory CI on theta and then does an end-point transformation to get the CI on phi the transformation of theta into the AR parameter.

theta2phi <- function (theta) theta/sqrt(1+theta^2) # AR1 link function, see https://cran.r-project.org/web/packages/glmmTMB/vignettes/covstruct.html

get_AR <- function(x){
    pars <- x$sdr$par.fixed
    sds <- sqrt(diag(x$sdr$cov.fixed))
    ind <- grep("theta", names(pars))
    t2 <- pars[ind[2]]
    s2 <- sds[ind[2]]
    list(
        est = theta2phi(unname(t2)),
        ci = theta2phi(qnorm(c(0.025, 0.975), t2, s2))
    )
}
get_AR(g2)
# $est
# [1] 0.6632351 
#
# $ci
# [1] 0.1793698 0.8465077 

Compare the result above with the model output below to confirm that the transformation to phi was successful. The AR1 parameter (Corr) in the summary below is the same (approx 0.66).

summary(g2)
#>  Family: gaussian  ( identity )
#> Formula:          x ~ ar1(times + 0 | group)
#> Data: dat0
#> 
#>      AIC      BIC   logLik deviance df.resid 
#>     68.7     73.6    -30.4     60.7       21 
#> 
#> Random effects:
#> 
#> Conditional model:
#>  Groups   Name   Variance  Std.Dev.  Corr      
#>  group    times1 1.158e+00 1.076e+00 0.66 (ar1)
#>  Residual        1.191e-09 3.451e-05           
#> Number of obs: 25, groups:  group, 1
#> 
#> Dispersion estimate for gaussian family (sigma^2): 1.19e-09 
#> 
#> Conditional model:
#>             Estimate Std. Error z value Pr(>|z|)
#> (Intercept)  0.06542    0.45390   0.144    0.885

^{Created on 2024-12-09 with reprex v2.1.0}