newpickle {aster} | R Documentation |
Penalized Quasi-Likelihood for Aster Models
Description
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 |
random |
the model matrix or matrices for random effects.
The number of rows is |
obj |
aster model object, the result of a call to |
y |
response vector. May be omitted, in which case |
origin |
origin of aster model. May be omitted, in which case
default origin (see |
zwz |
A possible value of |
deriv |
Number of derivatives wanted, either zero or one.
Must be zero if |
Details
Define
where is minus the log likelihood function of a saturated aster model,
where
is the Hessian matrix of
,
where
is a known vector (the offset vector in the terminology
of
glm
but the origin in the terminology
of aster
), where is a known matrix, the model matrix for
fixed effects (the argument
fixed
of this function),
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 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 is the identity matrix.
This function evaluates
when
zwz
is missing.
Otherwise it evaluates the same thing except that
is replaced by zwz
.
Note that is a function of
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
.
Actually no longer used by other functions in this package.
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)