Solved – Estimates of the variance of the variance component of a mixed effects model

anovamaximum likelihoodmixed modelr

Say $y=X\beta+ Zu +\epsilon$ is our mixed effects model where $u=(u_1,..,u_r)$ and $u_{j} \stackrel{i.i.d.}{\sim} N(0, \sigma^2_{a})$ for $j=1,…,r$ and $\epsilon=(\epsilon_1,…,\epsilon_n)$ are i.i.d. $N(0, \sigma^2_{b})$, furthermore $\epsilon_j$ and $u_i$ are also assumed to be independent for all $j$'s and all $i$'s.

I am interested in Var($\widehat{\sigma_{a}^2}$) i.e., the variance of the estimate of $\sigma^2_{a}$. In the R package ${\bf lme4}$ they do provide the M.L.Es but they do not provide the estimate of the Var($\widehat{\sigma_{a}^2}$). I don't think there is any closed form expression for the estimate of Var($\widehat{\sigma_{a}^2}$) and I was wondering if anyone knew of a R implementation (or any easily implementable algorithm) of how to calculate this quantity.

Best Answer

Here is the analysis with R-package VCA V1.2:

> library(VCA)
> data(sleepstudy)
> fit <- anovaMM(Reaction~Days*(Subject), sleepstudy)
> inf <- VCAinference(fit, VarVC=TRUE)
> print(inf, what="VCA")



Inference from Mixed Model Fit
------------------------------

> VCA Result:
-------------

        [Fixed Effects]

     int     Days 
251.4051  10.4673 


        [Variance Components]

  Name         DF    SS          MS         VC        %Total  SD      CV[%]   Var(VC)   
1 total        11.21                        1388.5416 100     37.2631 12.4831           
2 Subject      17    250618.1083 14742.2417 698.5289  50.3067 26.4297 8.8539  94751.0064
3 Days:Subject 17    60322.0013  3548.353   35.0717   2.5258  5.9221  1.9839  204.4845  
4 error        144   94311.5079  654.941    654.941   47.1675 25.5918 8.5732  5914.196  

Mean: 298.5079 (N = 180) 

Experimental Design: unbalanced

Fixed effects are equal and variance components of the ANOVA Type1-estimators are, except for Subject which is a bit larger (conservatively estimated), also equal to REML-estimators. Column "Var(VC)" contains variances of variance components according to Giebrecht and Burns (1985). The complete covariance matrix for variance components can also be extracted:

> vcovVC(fit)
               Subject Days:Subject       error
Subject      94751.006   -128.55799 -1523.85985
Days:Subject  -128.558    204.48451   -47.53872
error        -1523.860    -47.53872  5914.19600
attr(,"method")
[1] "gb"
Related Question