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