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

[Douglas Bates]

On Fri, Mar 30, 2012 at 3:44 PM, Joshua Wiley <jwiley.psych at gmail.com> wrote:
> 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.

I have been travelling and haven't tracked messages on the list
closely so I might have missed something here.  What do you mean by
"the variance parameters"?  The theta parameters aren't variances and
covariances.  They are values from the relative covariance factor.
For the model you fit the first element of theta is the standard
deviation of the intercept random effects divided by the standard
deviation of the per-observation noise.  The second and third elements
only make sense in terms of the Cholesky factor of the relative
variance-covariance matrix.

This choice is not arbitrary.  Although we are accustomed to thinking
in terms of variances and covariances or in terms of standard
deviations and correlations, these are difficult scales on which to
optimize because the constraints on these values are complicated
nonlinear constraints.

So as long as you recognize that the "variance parameters" aren't
variances or covariances you should be okay.

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