lme4 to get a vcov for fixed-effects variance

There are still some aspects of your question that aren't clear to me, 
but I might be able to shed some light on some things.

(1) when fitting a GLMM (rather than a LMM), lme4 *doesn't* use the RX 
component as described in the JSS paper, except as a fallback solution; 
it instead directly uses the inverse of the Hessian of the objective 
(deviance) function (H^(-1)), unless there appear to be numeric problems 
with H (e.g. it is non-positive-definite). (See code below.)

(2) the full covariance matrix estimated this way is *approximately* the 
same as you get for glmmTMB, except that glmmTMB is estimating 
log(theta) rather than theta.  After adjusting for this, the estimated 
variances disagree by about 10% (too bad, but not shocking for this kind 
of calculation).

(3) As for the actual process of extracting RX (if you wanted to compute 
that way) - it's spelled out in the paper. I'm not at all sure what you 
mean by "get this manually for each component".  To extract R_X and 
R_ZX, use getME() (see ?getME)

library(glmmTMB)
library(lme4)
## use only three species, to make cov matrices smaller/easier to examine
ss <- subset(Salamanders,spp %in% levels(spp)[1:3])
fit1 <- glmer(count~spp + (1|site),data=ss, family=poisson)
fit2 <- glmmTMB(count~spp + (1|site),data=ss, family=poisson)

## glmmTMB "full" covariance matrix
V1 <- vcov(fit2,full=TRUE)

## default glmer covariance matrix, converted back to
## a regular base-R matrix (from a Matrix).
V2 <- as.matrix(vcov(fit1))

## glmer covariance matrix using RX as in JSS paper.
V3 <- as.matrix(vcov(fit1,use.hessian=FALSE))
## warning message

##
neword <- c(2:4,1) ## put theta param last to match glmmTMB
V4 <- 2*solve(fit1 at optinfo$derivs$Hessian)[neword,neword]

## RX-based computation from scratch
sigma2 <- 1  ## sigma==1 for Poisson model
RX <- getME(fit1, "RX")
V5 <- sigma2 * chol2inv(RX)

all.equal(unname(V3),V5,tol=1e-12)

th <- getME(fit1,"theta")
## var(log(theta)) = var(theta)*d(log(theta))/d theta = var(theta)/theta^2
all.equal(V4[4,4]/unname(th)^2,V1[4,4])

lme4 to get a vcov for fixed-effects variance

Thread (2 messages)