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

[John Maindonald]

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
>