newpickle {aster}

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 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 Z^T 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) + c^T c / 2 + \log \det[A Z^T 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

Z^T 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. 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)