Post-fit Sampling Covariance of BLUPs

Computes the plug-in sampling variance or covariance of fitted random-effect predictors, $\hat u$, using Henderson's mixed-model coefficient system. This is distinct from the prior variance of the random effects and their prediction error variance (PEV).

Usage

postVarU(object, mode=1L)

Arguments

object: A Gaussian mmes Henderson fit with stored C and Cscale. Direct-engine, PQL and matrix-free solveOnly fits are not supported.
mode: 1 (default) returns diagonal BLUP sampling variances in uVarList; 2 additionally returns the joint sampling covariance matrix VarU; 0 clears only VarU-derived outputs.

Details

For $y=X\beta+Zu+e$, with $\operatorname{Var}(u)=G$ and $\operatorname{Var}(e)=R$, the unscaled Henderson matrix is $$C_0=\left[\begin{array}{cc} X'R^{-1}X & X'R^{-1}Z \\ Z'R^{-1}X & Z'R^{-1}Z+G^{-1} \end{array}\right].$$ Its random-effect inverse block is $Q=(C_0^{-1})_{uu}=\operatorname{Var}(u-\hat u)$. Consequently, $$\operatorname{Var}(\hat u)=G-Q.$$ For known covariance parameters this identity is exact. At estimated parameters the result is a plug-in covariance, omitting covariance-parameter estimation uncertainty. The same identity is used for ML and REML fits.

The implementation does not construct or invert the observation covariance $V=ZGZ'+R$. Error covariance is computed by solving the stored Henderson system and multiplying by Cscale to recover the original response scale. For each independent random term, $G_k=\Sigma_k\otimes A_k$, where $\Sigma_k$ is the fitted covariance-coordinate matrix and $A_k$ the relationship covariance. Products with $A_k$ are evaluated by solving its retained sparse precision, not by explicitly forming its inverse or the full joint $G$.

mode=1 processes identity contrasts in batches of at most 128 rows, so neither a full coefficient inverse nor a full random-effect covariance matrix is stored. It requires solves for all requested effect diagonals and can be expensive for large models. mode=2 requires quadratic storage in the number of random effects. Although independent terms have zero prior cross-covariance, their BLUPs can have nonzero sampling cross-covariance; the joint VarU includes these cross-term entries.

Rotation and factor-score augmentation are mapped to the public uList coordinates before computing error covariance. New fits retain randomPrecision in original level order, so external Gu objects are not needed after fitting. Older ordinary fits attempt to rebuild relationship precisions from saved model inputs, requiring the original formula objects; dimension checks cannot detect all changes to those objects. Older rotated or augmented fits must be refitted. No clipping is performed: small negative variances from numerical cancellation remain visible.

Existing Ci and uPevList are left unchanged. The variance here is not the empirical variance across BLUP values and is not a posterior variance of $u$. It is the repeated-sampling covariance of $\hat u$.

Value

The original fitted object with uVarList, a named list matching the dimensions and labels of uList, and uVarMode. With mode=2, VarU is the joint covariance in term order, with each uList matrix stacked column by column. Its labels contain term, covariance-coordinate and main-effect labels. Modes 0 and 1 remove any previous full VarU; mode 0 also removes uVarList.

Examples

data(DT_example, package="enhancer")
model <- mmes(Yield~Env, random=~Name, data=DT_example,
              computeCi=0, verbose=FALSE)
model <- postVarU(model)
head(model$uVarList[[1]])
full <- postVarU(model, mode=2)
head(diag(full$VarU))

prediction <- predict(model, D="Name", PEV=TRUE, VarU=TRUE)
head(prediction$pvals)
prediction$PEV[1:3, 1:3]
prediction$VarU[1:3, 1:3]
prediction$sampling.vcov[1:3, 1:3]