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
Independent random effects with US and AR1 structures using glmmTMB
1 message · Leigh Ann Starcevich