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 `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`. 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),
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)
```

[Package aster version 1.1-2 Index]