[R-sig-ME] reference for standard error of estimated variance components

[Ben Bolker]

  Another option for diagnosing problems with the fit is to
compute the Hessian of the variance-covariance parameters, which
are fitted with a Cholesky-factor parameterization ... this
doesn't actually get you the uncertainty of the variance-covariance
matrices, but it gives something that is arguably more useful.
https://github.com/lme4/lme4/issues/120 has some useful machinery
written by Rune Haubo which we eventually hope to integrate into
the package.

  This should be a lot cheaper than likelihood profiling too
(takes only O(p^2) function evaluations where p is the number
of model parameters)

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

Various diagnostics:

rcond(h)
diag(solve(h))

On 13-11-15 08:10 PM, John Maindonald wrote:
> See help(mcmcsamp), which tells you that mcmcsamp() is no longer
> available, but gives a number of other recourses. For example:
> 
>> library(DAAG) science1.lmer <- lmer(like ~ sex + PrivPub +
>> (1:school:class), data=science, na.action=na.exclude) 
>> confint(science1.lmer)
> Computing profile confidence intervals ...
>> round(confint(science1.lmer), 2)
> Computing profile confidence intervals ... 2.5 % 97.5 % .sig01
> 0.42   0.71 .sigma         1.68   1.82 (Intercept)    4.40   5.04 
> sexm          -0.01   0.37 PrivPubpublic  0.05   0.78
> 
> These are in the same ballpark as the estimates we gave in 'Data
> Analysis & Graphics Using R', CUP, 3rd edn 2010, p.316, using
> mcmcsamp():
> 
> Compare
>> round(confint(science1.lmer)[1:2,]^2, 2)  ## Convert to variances
> Computing profile confidence intervals ... 2.5 % 97.5 % .sig01  0.18
> 0.51 .sigma  2.83   3.30
> 
> with mcmcsamp(), from a pre-1.0 version of lme4 0.148 0.43 2.84
> 3.28
> 
> John Maindonald             email: john.maindonald at anu.edu.au phone :
> +61 2 (6125)3473    fax  : +61 2(6125)5549 Centre for Mathematics &
> Its Applications, Room 1194, John Dedman Mathematical Sciences
> Building (Building 27) Australian National University, Canberra ACT
> 0200.
> 
> 
> On 16/11/2013, at 10:49 AM, Jake Westfall <jake987722 at hotmail.com>
> wrote:
> 
>> Hi all,
>> 
>> Wondering if anyone knew off the top of their head, or could very
>> quickly locate, a reference describing how to compute the
>> variance/standard error of a variance component as, e.g., estimated
>> by lmer(). I am fully aware that this might not be a very sensible
>> thing to do because we strongly expect for the distribution of
>> these estimates to be asymmetric and etc. Let me explain my
>> motivation.
>> 
>> I am interested in these standard errors for purely for
>> troubleshooting purposes; specifically, helping diagnose situations
>> where we have a bad model specification given the data at hand. I
>> find that in many situations, particularly those involving crossed
>> random effects, lmer() happily provides estimates for models that
>> include certain variance components that we know are in-principle
>> inestimable given the structure of the dataset.
>> 
>> For example, random factors A and B are crossed, but there is only
>> 1 observation at each level of the crossing, so the A*B variance is
>> confounded with the residual variance, yet lmer() lists distinct
>> estimates for A*B variance and residual variance. Another example:
>> random A and random B are crossed, but B is nested in fixed X, yet
>> somehow lmer() will provide estimates for both the random A*X
>> interaction (which can be estimated fine) and also a random B*X
>> interaction (which cannot be).
>> 
>> Furthermore, there is no error message indicating that anything
>> fishy is going on in these cases, and the estimated variances even
>> look rather sensible (I suspect it is doing something like
>> double-counting random terms that are confounded and thus providing
>> similar but non-identical estimates for the two confounded
>> variances, but I don't know). So a less experienced user may never
>> get any clue that the model does not make sense given the data.
>> (FYI, I report this behavior using lme4 version < 1.0, and have not
>> yet verified on version > 1.0, but my understanding is that if
>> anything the newer lme4 functions tend to provide *fewer* warnings
>> and errors rather than more, so I would be pretty surprised if
>> everything blows up the way it really ought to with these models in
>> the newest versions.)
>> 
>> My suspicion is that if one were to compute the standard errors of
>> the confounded variance components, they would be unreasonably
>> large, and this might serve as the only obvious clue to many users
>> that all is not well. So I want to try this to see if it works as a
>> reasonably informative indicator of model misspecification.
>> 
>> I was sure I would be able to find the variance formula I needed
>> from Littell et al. "SAS System for Mixed Models," but no luck so
>> far. Googling has also not turned it up yet, although some
>> searching remains to be done. So I thought I would ask here to see
>> if any of you had leads.
>> 
>> Jake [[alternative HTML version deleted]]
>> 
>> _______________________________________________ 
>> R-sig-mixed-models at r-project.org mailing list 
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>> 
> 
> _______________________________________________ 
> R-sig-mixed-models at r-project.org mailing list 
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>