Evaluates the objective function for approximate maximum likelihood for an aster model with random effects. Uses Laplace approximation to integrate out the random effects analytically. The “quasi” in the title is a misnomer in the context of aster models but the acronym PQL for this procedure is well-established in the generalized linear mixed models literature.
quickle(alphanu, bee, fixed, random, obj, y, origin, zwz, deriv = 0)a list with some of the following components: value, gradient,
hessian, alpha, bee, nu. The first three are
the requested derivatives. The second three are the corresponding parameter
values: alpha and nu are the corresponding parts of the
argument alphanu, the value of bee is the result of the inner
optimization (\(b^*\) in the notation in details),
not the argument bee of this function.
the parameter vector value where the function is evaluated, a numeric vector, see details.
the random effects vector that is used as the starting point
for the inner optimization, which maximizes the penalized log likelihood
to find the optimal random effects vector matching alphanu.
the model matrix for fixed effects. The number of rows
is nrow(obj$data). The number of columns is the number of fixed
effects.
the model matrix or matrices for random effects.
The number of rows is nrow(obj$data). The number of columns
is the number of random effects in a group. Either a matrix or a list
each element of which is a matrix.
aster model object, the result of a call to aster.
response vector. May be omitted, in which case obj$x
is used. If supplied, must be a matrix of the same dimensions as
obj$x.
origin of aster model. May be omitted, in which case
default origin (see aster) is used. If supplied, must be
a matrix of the same dimensions obj$x.
A possible value of \(Z^T W Z\), where \(Z\) is the
model matrix for all random effects and \(W\) is the variance matrix of
the response. See details. Typically constructed by
the function makezwz.
Number of derivatives wanted, zero, one, or two.
Define
$$p(\alpha, b, \nu) = m(a + M \alpha + Z b) + {\textstyle \frac{1}{2}} b^T D^{- 1} b + {\textstyle \frac{1}{2}} \log \det[Z^T W Z D + I]$$
where \(m\) is minus the log likelihood function of a saturated aster model,
where \(a\) is a known vector (the offset vector in the terminology
of glm but the origin in the terminology
of aster),
where \(M\) is a known matrix, the model matrix for fixed effects
(the argument fixed of this function),
where \(Z\) is a known matrix, the model matrix for random effects
(either the argument random of this function if it is a matrix or
Reduce(cbind, random) if random is a list of matrices),
where \(D\) is a diagonal matrix whose diagonal is the vector
rep(nu, times = nrand)
where nrand is sapply(random, ncol)
when random is a list of
matrices and ncol(random) when random is a matrix,
where \(W\) is an arbitrary symmetric positive semidefinite matrix
(\(Z^T W Z\) is the argument zwz of this function),
and where \(I\) is the identity matrix.
Note that \(D\) is a function of \(\nu\)
although the notation does not explicitly indicate this.
The argument alphanu of this function is the concatenation
of the parameter vectors \(\alpha\) and \(\nu\).
The argument bee of this function is a possible value of \(b\).
The length of \(\alpha\) is the column dimension of \(M\).
The length of \(b\) is the column dimension of \(Z\).
The length of \(\nu\) is the length of the argument random
of this function if it is a list and is one otherwise.
Let \(b^*\) denote the minimizer of \(p(\alpha, b, \nu)\) considered as a function of \(b\) for fixed \(\alpha\) and \(\nu\), so \(b^*\) is a function of \(\alpha\) and \(\nu\). This function evaluates $$q(\alpha, \nu) = p(\alpha, b^*, \nu)$$ and its gradient vector and Hessian matrix (if requested). Note that \(b^*\) is a function of \(\alpha\) and \(\nu\) although the notation does not explicitly indicate this.
data(radish)
pred <- c(0,1,2)
fam <- c(1,3,2)
rout <- reaster(resp ~ varb + fit : (Site * Region),
list(block = ~ 0 + fit : Block, pop = ~ 0 + fit : Pop),
pred, fam, varb, id, root, data = radish)
alpha.mle <- rout$alpha
bee.mle <- rout$b
nu.mle <- rout$sigma^2
zwz.mle <- rout$zwz
obj <- rout$obj
fixed <- rout$fixed
random <- rout$random
alphanu.mle <- c(alpha.mle, nu.mle)
qout <- quickle(alphanu.mle, bee.mle, fixed, random, obj,
zwz = zwz.mle, deriv = 2)
Run the code above in your browser using DataLab