[R-sig-ME] lme4 to get a vcov for fixed-effects variance

[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 using 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:
> 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 using r-project.org
> 
> 
> 	[[alternative HTML version deleted]]
> 
> 
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>