Skip to content

extracting gradient and hessian matrix from lme4

3 messages · Joshua Wiley, Ben Bolker

#
Hi,

I am trying to use a multivariate mixed effects linear model to
examine mediation.  This works fine.  The final step is to compute the
indirect effect and its standard error.  The indirect effect is easy
(product of coefficients plus their covariance).  For the standard
error, I need the gradient (D) and the hessian (H):
the variance is then:

D'\Sigma(\theta)D + \frac{1}{2}Tr{(\mathbf{H}\boldsymbol{\Sigma}(\theta))^{2}}

This is all given in the Appendix of
http://www.unc.edu/~dbauer/manuscripts/bauer-preacher-gil-PM-2006.pdf

Is there a way to get this out of a mer class object?  Looking at
class?mer, the appropriate bits of vcov give me $\Sigma(\theta)$.  @V
seems like it would give me the gradient but is null for a basic lmer
model.

Thanks,

Josh
1 day later
#
Joshua Wiley <jwiley.psych at ...> writes:
If you're willing to try out the development version (i.e., lme4
from r-forge), I think you can do this as follows:

fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
fm1Fun <- update(fm1,devFunOnly=TRUE)
library(numDeriv)
fm1_thpar <- getME(fm1,"theta")
h <- hessian(fm1Fun,fm1_thpar)

  and similarly for the gradient.

  Let me know how it goes.

  Ben Bolker
11 days later
#
Hi Ben,

Many thanks for the help.  I tried your suggestion out and it seemed
to work (and I learned a bit about lme4 in the process :)

library(lme4)
fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
fm1Fun <- update(fm1,devFunOnly=TRUE)
library(numDeriv)
fm1_thpar <- getME(fm1,"theta")
h <- hessian(fm1Fun, fm1_thpar)
g <- grad(fm1Fun, fm1_thpar)

which I can use (I think) to get standard errors of the variance parameters.

library(MASS)
sqrt(diag(ginv(h)))

which plug into a longer formula that attempts to test the
significance of indirect effects.  Given that the variance parameters
are not normally distributed, my hunch is that even though both fixed
and random effects (and their variances/standard errors) are being
built into the mediation test, it is probably not well-behaved either,
but it is nice to be able to try to replicate models in the article.
Even if they are not perfectly accurate, I am hoping I can use them as
a sanity check for when I play with some mcmc and bootstrapping.

Thanks again!

Josh

Session info below just for the record if anyone else is trying to try this.

R Under development (unstable) (2012-03-29 r58868)
Platform: x86_64-pc-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] MASS_7.3-17        reshape2_1.2.1     numDeriv_2012.3-1  lme4_0.999902344-0
[5] Matrix_1.0-6       lattice_0.20-6

loaded via a namespace (and not attached):
[1] grid_2.16.0    minqa_1.2.0    nlme_3.1-103   plyr_1.7.1     splines_2.16.0
[6] stringr_0.6    tools_2.16.0
On Mon, Mar 19, 2012 at 12:15 PM, Ben Bolker <bbolker at gmail.com> wrote: