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 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.
Usage
mat.Vpredicts(target, Gt = 0, fixed = ~ 1, random, G, R, design,
eliminate, keep.order = TRUE, result = "variance.matrix")
Arguments
target 
The 
Gt 
The value of the variance component for the targetted effects or the

fixed 
The 
random 
A 
G 
This term only needs to be set if 
R 
The 
design 
A 
eliminate 
The 
keep.order 
A 
result 
A 
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 MoorePenrose 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.5E08). 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, 1731.
Smith, A. B., D. G. Butler, C. R. Cavanagh and B. R. Cullis (2015). Multiphase variety trials using both composite and individual replicate samples: a modelbased design approach. Journal of Agricultural Science, 153, 10171029.
See Also
designAmeasures
, mat.random
, mat.Vpred
.
Examples
## Reduced example from Smith et al. (2015)
## Generate twophase 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)