[R-sig-ME] extracting gradient and hessian matrix from lme4

[Joshua Wiley]

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:
> 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
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models



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