dae (version 3.1-23)

mat.Vpredicts: Calculates the variances of a set of predicted effects from a mixed model, based on supplied matrices or formulae.

Description

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 formulae 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 formulae.

Usage

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

The n x r incidence matrix for the r nuisance randm effects, or a formula from which the matrix can be generated. The default is 0, which implies that there are no random effects.

G

The r x r variance matrix of the r nuisance random effects. 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.

Value

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 \(\bold{Y} = \bold{X\beta} + \bold{Ww} + \bold{Zu} + \bold{e}\), where \(\bold{W}\) is the incidence matrix for the target predicted effects \(\bold{w}\), \(\bold{X}\) is the is incidence matrix for the fixed nuisance effects \(\bold{\beta}\), \(\bold{Z}\) is the is incidence matrix for the random nuisance effects \(\bold{u}\), \(\bold{e}\) are the residuals; the \(\bold{u}\) are assumed to have variance matrix \(\bold{G}\) so that their contribution to the variance matrix for \(\bold{Y}\) is \(\bold{Vu} = \bold{ZGZ}^T\) and \(\bold{e}\) is assumed to have variance matrix \(\bold{R}\). If the target effects are random then the variance matrix for \(\bold{w}\) is \(\bold{G}_t\) so that their contribution to the variance matrix for \(\bold{Y}\) is \(\bold{WG}_t\bold{W}^T\).

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.

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.

See Also

designAmeasures, mat.Vpred.

Examples

# NOT RUN {
## 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)

## 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)
# }