extracting gradient and hessian matrix from lme4

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.

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:

Joshua Wiley <jwiley.psych at ...> writes:

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.

?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

--
Joshua Wiley
Ph.D. Student, Health Psychology
Programmer Analyst II, Statistical Consulting Group
University of California, Los Angeles
https://joshuawiley.com/