R: Calculates the variances of a set of predicted effects from a...

mat.Vpredicts {dae}

R Documentation

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

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.

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.

Author(s)

Chris Brien

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

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

[Package dae version 3.2.28 Index]