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

[Ben Bolker]

On 12-04-01 03:50 AM, Douglas Bates wrote:
> 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.

  The warning is worthwhile.  I haven't taken enough time to figure out
the context (Appendix of

http://www.unc.edu/~dbauer/manuscripts/bauer-preacher-gil-PM-2006.pdf

as stated), but to me it actually seems that the context there is
actually using the RE variances as plug-in estimates (i.e. ignoring
their uncertainty) ... ??


> 
>> 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