[R-sig-ME] Independent random effects with US and AR1 structures using glmmTMB

[Leigh Ann Starcevich]

I am modeling counts collected at several Visits at several Sites.  The
data structure is similar to the "dp" data set in this document, where f =
Site and tt = Visit:

http://bbolker.github.io/mixedmodels-misc/notes/corr_braindump.html

To reproduce, use the following function from
https://github.com/bbolker/mixedmodels-misc/blob/master/notes/lme4ord_poiscorsim.R
:

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

library(glmmTMB)
set.seed(101)


# Simulate data
dp <- simCor1(phi=0.8,sdgrp=2,sdres=1,seed=101,
              linkinv=exp,simfun=function(x) rpois(length(x),lambda=x))

This model is used:

glmmTMB_pois_fit <- glmmTMB(y~1 + (1|f) + ar1(tt-1|f), data=dp, family
= poisson)

and the following output is obtained:

#  Groups Name        Variance Std.Dev. Corr

#  f      (Intercept) 5.3065   2.3036
#  f.1    tt1         0.9996   0.9998   0.77 (ar1)
# Number of obs: 400, groups:  f, 20

# Conditional model:
#             Estimate Std. Error z value Pr(>|z|)
# (Intercept)  -0.8319     0.5626  -1.479    0.139

summary(glmmTMB_pois_fit)$varcor$cond

# only one correlation matrix provided in output

Because we are modeling the AR1 structure by each level of the grouping
factor, f, I expected to obtain a different AR1 correlation estimate for
each level of f.  However, the model summary indicates that only one AR1
correlation is estimated.  Is this a mean across levels of f?

If the AR1 correlation structure is modeled across all Sites for a
single-level group variable, ff, how are these parameterizations
different?  The results are VERY different:

# Model one AR1 correlation across levels of f using single-level grouping
factor ff
dp$ff = as.factor(rep(1,nrow(dp)))
glmmTMB_pois_fit2 <- glmmTMB(y~1 + (1|f) + ar1(tt-1|ff), data=dp,
 family=poisson)
summary(glmmTMB_pois_fit2)

# Conditional model:
#  Groups Name        Variance Std.Dev. Corr
#  f      (Intercept) 5.5511   2.3561
#  ff     tt1         0.1049   0.3238   0.04 (ar1)
# Number of obs: 400, groups:  f, 20; ff, 1
#
# Conditional model:
#             Estimate Std. Error z value Pr(>|z|)
# (Intercept)  -0.5406     0.5616  -0.963    0.336

summary(glmmTMB_pois_fit)$varcor$cond

# only one correlation matrix provided in output

Thanks for your help,

--Leigh Ann Starcevich

	[[alternative HTML version deleted]]