Skip to content

lme4 to get a vcov for fixed-effects variance

2 messages · wilson 1998, Ben Bolker

Original
wilson 1998
# 
Following the previous email to make it just text only questions: I?ll restate my questions:

Hi Ben Bolker,

I would like to know about some mechanism about vcov output from glmer model from here<https://cran.r-project.org/web/packages/lme4/vignettes/lmer.pdf>

So I am running a model using glmer under poisson model
fit1 = glmer(count~spp + (1|site),data=Salamanders, family=poisson)
using Salamanders data<https://www.rdocumentation.org/packages/glmmTMB/versions/1.0.2.1/topics/Salamanders> from glmmTMB package in R Studio

Then it could produce the variance covariance matrix about the fixed effects given below:, also since full=T it also gives the extra column and extra row to account for the random effect as well.

vcov(fit1,full=T)


I noticed in page 26 in section 5.1 in lme4 documentation (equation 54 and equation 55) from this link: https://cran.r-project.org/web/packages/lme4/vignettes/lmer.pdf#page=14&zoom=100,568,614


However, How do I reproduce the output for vcov manually? I know that under equation (55) it says that they could reproduce using chol2inv(RX) to represent R_X^-1 (R^T_X)^-1

I don?t get it how to produce the R_X, L(theta) ,R_ZX as well? It is stated that there could be Cholesky decomposition, but what if I would like to get this manually for each component?
The main point is that I would like to know how could the lme4 package get the output for the vcov(fit1,full=T). And in this case, I also would like to know how to extract R_X and R_ZX from the fit1?

Thankyou.

Regards,
Wilson



Below is the previous email with the comment

Your request to the R-sig-mixed-models mailing list

    Posting of your message titled "lme4 to get a vcov for
fixed-effects variance"

has been rejected by the list moderator.  The moderator gave the
following reason for rejecting your request:

"Could you please resend without the embedded text from the JSS paper?
(The mailing list is very old-fashioned and only likes plain-text
messages.) You can point to the specific section in the JSS paper
(e.g. "in section 5.1 'Theory underlying the output module', in the
first subsection 'Covariance matrix of the fixed-effect
coefficients').

If you could also clarify your question that would be useful: do you
want to know how RX is computed during the estimation procedure (see
eq 51 of the JSS paper), or how to extract it from the fit [which is
stated in the text you quoted], or ... ?

  sincerely
   Ben Bolker"

Any questions or comments should be directed to the list administrator
at:

    r-sig-mixed-models-owner at r-project.org
Ben Bolker
# 
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])
On 2/24/21 8:30 PM, wilson 1998 wrote: