# Calculates the variances of a set of predicted effects from a mixed model, based on supplied matrices or formulae. For `n` observations, `w` effects to be predicted, `f` nuiscance fixed effects, `r` nuisance random effects and `n` residuals, the variances of a set of predicted effects is calculated using the incidence matrix for the effects to be predicted and, optionally, a variance matrix of these effects, an incidence matrix for the nuisance fixed factors and covariates, the variance matrix of the nuisance random effects and the residual variance matrix. The `matrices` can be supplied directly or using `formula`e and a `matrix` specifying the variances of the nuisance random effects. The difference between `mat.Vpredicts` and `mat.Vpred` is that the former has different names for equivalent arguments and the latter does not allow for the use of `formula`e. ```r mat.Vpredicts(target, Gt = 0, fixed = ~ 1, random, G, R, design, eliminate, keep.order = TRUE, result = "variance.matrix") ``` ## Arguments - `target`: The `n x w` incidence `matrix` for the `w` effects targetted for prediction, or a `formula` from which the `matrix` can be generated. - `Gt`: The value of the variance component for the targetted effects or the `w x w` variance `matrix` of the `w` targetted effects. If the targetted effects are fixed, set `Gt` to 0. - `fixed`: The `n x f` incidence `matrix` for the `f` nuisance fixed effects and covariates, or a `formula` from which the `matrix` can be generated. The default is a `formula` for an intercept-only model. - `random`: A `formula` or a `matrix`. If a `formula`, it specifies the random effects from which the `matrix` for the contribution of the random effects to the variance matrix can be generated. If it is a matrix, it must be an `n x n` `matrix` and will be passed on to form the variance matrix of the observations. The default is 0, which implies that there are no random effects. - `G`: This term only needs to be set if `random` is a `formula`. Then it is set to a `list`, in which each component is either a single value or a `matrix`; there needs to be a component for each term in the expanded `formula`, with the order of the terms and components matching. If it is a single value, a diagonal matrix of dimension equal to the product of the numbers of levels of the factors in its term. If it is a matrix, its dimension must be equal to the product of the numbers of levels of the factors in its term. - `R`: The `n x n` residual variance `matrix`. If `R` is not set in the call, then it is set to the identity `matrix`. - `design`: A `data.frame` containing the design to be used in an experiment from which predictions are to be obtained. It is not required when the only `formula` specified is an intercept-only `formula`. - `eliminate`: The `n x n` `projector` onto the subspace corresponding to the effects to be eliminated from the information matrix prior to inverting it to form the variance `matrix` of the predicted effects. It is only appropriate to use this option when the effects to be predicted are fixed. - `keep.order`: A `logical` indicating whether the terms should keep their position in the expanded `formula` projector, or reordered so that main effects precede two-factor interactions, which precede three-factor interactions and so on. - `result`: A `character` indicating which matrix is to be returned: `variance.matrix` or `information.matrix`. ## Returns A `w x w` `matrix` containing the variances and covariances of the predicted effects or the information matrix for the effects, depending on the setting of `result`. The matrix has its rank as an attribute. ## Details The mixed model for which the predictions are to be obtained is of the form $Y = X\beta + Ww + Zu + e$, where $W$ is the incidence matrix for the `target` predicted effects $w$, $X$ is the is incidence matrix for the `fixed` nuisance effects $\beta$, $Z$ is the is incidence matrix for the `random` nuisance effects $u$, $e$ are the residuals; the $u$ are assumed to have variance matrix $G$ so that their contribution to the variance matrix for $Y$ is $Vu = Z %*% G %*% t(Z)$ and $e$ is assumed to have variance matrix $R$. If the `target` effects are random then the variance matrix for $w$ is $Gt$ so that their contribution to the variance matrix for $Y$ is $W %*% Gt %*% t(W)$. As described in Hooks et al. (2009, Equation 19), the information matrix is calculated as `A <- t(W) %*% Vinv %*% W + ginv(Gg) - A%*%ginv(t(X)%*%Vinv%*%X)%*%t(A)`, where `Vinv <- ginv(Vu + R)`, `A = t(W) %*% Vinv %*% X` and ginv(B) is the unique Moore-Penrose inverse of B formed using the eigendecomposition of B. Then, if `eliminate` is set and the effects to be predicted are fixed then the reduced information matrix is calculated as `A <- (I - eliminate) Vinv (I - eliminate)`. Finally, if `result` is set to `variance.matrix`, the variance of the predicted effects is calculated: `Vpred <- ginv(A)` and returned; otherwise the information matrix A is returned. The rank of the matrix to be returned is obtain via a singular value decomposition of the information matrix, it being the number of nonzero eigenvalues. An eigenvalue is regarded as zero if it is less than `daeTolerance`, which is initially set to.Machine$double.eps ^ 0.5 (about 1.5E-08). The function `set.daeTolerance` can be used to change `daeTolerance`. ## Author(s) Chris Brien ## See Also `designAmeasures`, `mat.random`, `mat.Vpred`. ## References Hooks, T., Marx, D., Kachman, S., and Pedersen, J. (2009). Optimality criteria for models with random effects. **Revista Colombiana de Estadistica**, 32 , 17-31. Smith, A. B., D. G. Butler, C. R. Cavanagh and B. R. Cullis (2015). Multi-phase variety trials using both composite and individual replicate samples: a model-based design approach. **Journal of Agricultural Science**, 153 , 1017-1029. ## Examples ```r ## Reduced example from Smith et al. (2015) ## Generate two-phase design mill.fac <- fac.gen(list(Mrep = 2, Mday = 2, Mord = 3)) field.lay <- fac.gen(list(Frep = 2, Fplot = 4)) field.lay$Variety <- factor(c("D","E","Y","W","G","D","E","M"), levels = c("Y","W","G","M","D","E")) start.design <- cbind(mill.fac, field.lay[c(3,4,5,8,1,7,3,4,5,8,6,2),]) rownames(start.design) <- NULL ## Set gammas terms <- c("Variety", "Frep", "Frep:Fplot", "Mrep", "Mrep:Mday", "Mrep:Mday:Mord") gammas <- c(1, 0.1, 0.2, 0.3, 0.2, 1) names(gammas) <- terms ## Specify matrices to calculate the variance matrix of the predicted fixed Variety effects W <- model.matrix(~ -1 + Variety, start.design) Vu <- with(start.design, fac.vcmat(Mrep, gammas["Mrep"]) + fac.vcmat(fac.combine(list(Mrep,Mday)), gammas["Mrep:Mday"]) + fac.vcmat(Frep, gammas["Frep"]) + fac.vcmat(fac.combine(list(Frep,Fplot)), gammas["Frep:Fplot"])) R <- diag(1, nrow(start.design)) ## Calculate variance matrix Vp <- mat.Vpredicts(target = W, random=Vu, R=R, design = start.design) ## Calculate the variance matrix of the predicted random Variety effects using formulae Vp <- mat.Vpredicts(target = ~ -1 + Variety, Gt = 1, fixed = ~ 1, random = ~ -1 + Mrep/Mday + Frep/Fplot, G = as.list(gammas[c(4,5,2,3)]), R = R, design = start.design) designAmeasures(Vp) ## Calculate the variance matrix of the predicted fixed Variety effects, ## elminating the grand mean n <- nrow(start.design) Vp.reduc <- mat.Vpredicts(target = ~ -1 + Variety, random = ~ -1 + Mrep/Mday + Frep/Fplot, G = as.list(gammas[c(4,5,2,3)]), eliminate = projector(matrix(1, nrow = n, ncol = n)/n), design = start.design) designAmeasures(Vp.reduc) ```

mat.Vpredicts function