newpickle {aster} | R Documentation |

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.

newpickle(alphaceesigma, fixed, random, obj, y, origin, zwz, deriv = 0)

`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 |

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.

a list with components `value`

and `gradient`

,
the latter missing if `deriv == 0`

.

Not intended for use by naive users. Use `reaster`

.
Actually no longer used by other functions in this package.

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)

[Package *aster* version 1.1-2 Index]