How should standard errors for mixed effects model estimates be calculated?

Question

In particular, how should the standard errors of the fixed effects in a linear mixed effects model be calculated (in a frequentist sense)?

I have been lead to believe that the typical estimates ($Var(\hat\beta)=(X'VX)^{-1}$), such as those presented in Laird and Ware [1982] will give SE's that are underestimated in size because the estimated variance components are treated as though they are the true values.

I have noticed that the SE's produced by the lme and summary functions in the nlme package for R are not simply equal to the square root of the diagonals of the variance-covariance matrix given above. How are they calculated?

I am also under the impression, that Bayesians use inverse gamma priors for the estimation of variance components? Do these give the same results (in the right setting) as lme?

Answer 1

My initial thought was that, for ordinary linear regression, we just plug in our estimate of the residual variance, $\sigma^2$, as if it were the truth.

However, take a look at McCulloch and Searle (2001) Generalized, linear and mixed models, 1st edition, Section 6.4b, "Sampling variance". They indicate that you can't just plug in the estimates of the variance components:

Instead of dealing with the variance (matrix) of a vector $X \hat{\beta}$ we consider the simpler case of the scalar $l^' \hat{\beta}$ for estimable $l^'\beta$ (i.e., $l^' = t^'X$ for some $t^'$).

For known $V$, we have from (6.21) that $\text{var}(l^' \beta^0) = l^'(X^'V^{-1}X)^- l$. A replacement for this when $V$ is not known is to use $l^'(X^'\hat{V}^{-1}X)^- l$, which is an estimate of $\text{var}(l^'\beta^0) = \text{var}[l^' (X^' V^{-1} X)^- X^' V^{-1} y]$. But it is not an estimate of $\text{var}(l^'\hat{\beta}) = \text{var}[l^' (X^' \hat{V}^{-1} X)^- X^' \hat{V}^{-1} y]$. The latter requires taking account of the variability of $\hat{V}$ as well as that in $y$. To deal with this, Kackar and Harville (1984, p. 854) observe that (in our notation) $l^' \hat{\beta} - l^' \beta$ can be expressed as the sum of two independent parts, $l^' \hat{\beta} - l^' \beta^0$ and $l^'\beta^0 - l^' \beta$. This leads to $\text{var}(l^' \hat{\beta})$ being expressed as a sum of two variances which we write as

$$\text{var}(l^'\hat{\beta}) = ... \approx l^'(X^' V^{-1} X)l + l^' T \; l$$

They go on to explain $T$.

So this answers the first part of your question and indicates that your intuition was correct (and mine was wrong).