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.
Usage
newpickle(alphaceesigma, fixed, random, obj, y, origin, zwz, deriv = 0)
Arguments
alphaceesigma
the parameter value where the function is evaluated,
a numeric vector, see details.
fixed
the model matrix for fixed effects. The number of rows
is nrow(obj$data). The number of columns is the number of fixed
effects.
random
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.
obj
aster model object, the result of a call to aster.
y
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
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.
zwz
A possible value of t(Z) W Z, where Z is the
model matrix for all random effects and W is the variance matrix of
the response. May be missing, in which case it is calculated from
alphaceesigma. See details.
deriv
Number of derivatives wanted, either zero or one.
Must be zero if zwz is missing.
Details
Define
p(alpha, c, sigma) = m(a + M alpha + Z A c) + t(c) c / 2 + log det[A t(Z) W(a + M alpha + Z A c) Z A + I]
where m is minus the log likelihood function of a saturated aster model,
where W is the Hessian matrix of m,
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),
Z is a known matrix, the model matrix for random effects
(either the argument random of this functions if it is a matrix or
Reduce(cbind, random) if random is a list of matrices),
where A is a diagonal matrix whose diagonal is the vector
rep(sigma, times = nrand)
where nrand is sapply(random, ncol)
when random is a list of
matrices and ncol(random) when random is a matrix,
and where I is the identity matrix.
This function evaluates p(alpha, c, sigma)
when zwz is missing.
Otherwise it evaluates the same thing except that
t(Z) W(a + M alpha + Z A c) Z
is replaced by zwz.
Note that A is a function of sigma although the
notation does not explicitly indicate this.
Value
a list with components value and gradient,
the latter missing if deriv == 0.
Note
Not intended for use by naive users. Use reaster,
which calls them.
Examples
data(radish)
pred <- c(0,1,2)
fam <- c(1,3,2)
### need object of type aster to supply to penmlogl and pickle
aout <- aster(resp ~ varb + fit : (Site * Region + Block + Pop),
pred, fam, varb, id, root, data = radish)
### model matrices for fixed and random effects
modmat.fix <- model.matrix(resp ~ varb + fit : (Site * Region),
data = radish)
modmat.blk <- model.matrix(resp ~ 0 + fit:Block, data = radish)
modmat.pop <- model.matrix(resp ~ 0 + fit:Pop, data = radish)
rownames(modmat.fix) <- NULL
rownames(modmat.blk) <- NULL
rownames(modmat.pop) <- NULL
idrop <- match(aout$dropped, colnames(modmat.fix))
idrop <- idrop[! is.na(idrop)]
modmat.fix <- modmat.fix[ , - idrop]
nfix <- ncol(modmat.fix)
nblk <- ncol(modmat.blk)
npop <- ncol(modmat.pop)
alpha.start <- aout$coefficients[match(colnames(modmat.fix),
names(aout$coefficients))]
cee.start <- rep(0, nblk + npop)
sigma.start <- rep(1, 2)
alphaceesigma.start <- c(alpha.start, cee.start, sigma.start)
foo <- newpickle(alphaceesigma.start, fixed = modmat.fix,
random = list(modmat.blk, modmat.pop), obj = aout)