https://github.com/cran/aster
Raw File
Tip revision: aa47935123bfca8a22cbc8345d658d0c1713a289 authored by Charles J. Geyer on 14 December 2023, 15:20:02 UTC
version 1.1-3
Tip revision: aa47935
quickle.Rd
\name{quickle}
\alias{quickle}
\title{Penalized Quasi-Likelihood for Aster Models}
\usage{
quickle(alphanu, bee, fixed, random, obj, y, origin, zwz, deriv = 0)
}
\note{
  Not intended for use by naive users.  Use \code{\link{summary.reaster}},
  which calls it.
}
\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 \dQuote{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.
}
\arguments{
  \item{alphanu}{the parameter vector value where the function is evaluated,
    a numeric vector, see details.}
  \item{bee}{the random effects vector that is used as the starting point
    for the inner optimization, which maximizes the penalized log likelihood
    to find the optimal random effects vector matching \code{alphanu}.}
  \item{fixed}{the model matrix for fixed effects.  The number of rows
    is \code{nrow(obj$data)}.  The number of columns is the number of fixed
    effects.}
  \item{random}{the model matrix or matrices for random effects.
    The number of rows is \code{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.}
  \item{obj}{aster model object, the result of a call to \code{\link{aster}}.}
  \item{y}{response vector.  May be omitted, in which case \code{obj$x}
    is used.  If supplied, must be a matrix of the same dimensions as
    \code{obj$x}.}
  \item{origin}{origin of aster model.  May be omitted, in which case
    default origin (see \code{\link{aster}}) is used.  If supplied, must be
    a matrix of the same dimensions \code{obj$x}.}
  \item{zwz}{A possible value of \eqn{Z^T W Z}{t(Z) W Z}, where \eqn{Z} is the
    model matrix for all random effects and \eqn{W} is the variance matrix of
    the response.  See details.  Typically constructed by
    the function \code{\link{makezwz}}.}
  \item{deriv}{Number of derivatives wanted, zero, one, or two.}
}
\details{
Define
\deqn{p(\alpha, b, \nu) = m(a + M \alpha + Z b) + {\textstyle \frac{1}{2}} b^T D^{- 1} b + {\textstyle \frac{1}{2}} \log \det[Z^T W Z D + I]}{p(alpha, b, nu) = m(a + M alpha + Z b) + t(b) D^(- 1) b / 2 + log det[t(Z) W Z D + I] / 2}
where \eqn{m} is minus the log likelihood function of a saturated aster model,
where \eqn{a} is a known vector (the \emph{offset vector} in the terminology
of \code{\link{glm}} but the \emph{origin} in the terminology
of \code{\link{aster}}),
where \eqn{M} is a known matrix, the model matrix for fixed effects
(the argument \code{fixed} of this function),
where \eqn{Z} is a known matrix, the model matrix for random effects
(either the argument \code{random} of this function if it is a matrix or
\code{Reduce(cbind, random)} if \code{random} is a list of matrices),
where \eqn{D} is a diagonal matrix whose diagonal is the vector
\code{rep(nu, times = nrand)}
where \code{nrand} is \code{sapply(random, ncol)}
when \code{random} is a list of
matrices and \code{ncol(random)} when \code{random} is a matrix,
where \eqn{W} is an arbitrary symmetric positive semidefinite matrix
(\eqn{Z^T W Z}{t(Z) W Z} is the argument \code{zwz} of this function),
and where \eqn{I} is the identity matrix.
Note that \eqn{D} is a function of \eqn{\nu}{nu}
although the notation does not explicitly indicate this.

The argument \code{alphanu} of this function is the concatenation
of the parameter vectors \eqn{\alpha}{alpha} and \eqn{\nu}{\nu}.
The argument \code{bee} of this function is a possible value of \eqn{b}.
The length of \eqn{\alpha}{alpha} is the column dimension of \eqn{M}.
The length of \eqn{b} is the column dimension of \eqn{Z}.
The length of \eqn{\nu} is the length of the argument \code{random}
of this function if it is a list and is one otherwise.

Let \eqn{b^*}{bstar} denote the minimizer
of \eqn{p(\alpha, b, \nu)}{p(alpha, b, nu)} considered as a function of
\eqn{b} for fixed \eqn{\alpha}{alpha} and \eqn{\nu}{nu}, so \eqn{b^*}{bstar}
is a function of \eqn{\alpha}{alpha} and \eqn{\nu}{nu}.
This function evaluates
\deqn{q(\alpha, \nu) = p(\alpha, b^*, \nu)}{q(alpha, nu) = p(alpha, bstar, nu)}
and its gradient vector and Hessian matrix (if requested).
Note that \eqn{b^*}{bstar} is a function of \eqn{\alpha}{alpha}
and \eqn{\nu}{nu} although the notation does not explicitly indicate this.
}
\value{
  a list with some of the following components: \code{value}, \code{gradient},
  \code{hessian}, \code{alpha}, \code{bee}, \code{nu}.  The first three are
  the requested derivatives.  The second three are the corresponding parameter
  values: \code{alpha} and \code{nu} are the corresponding parts of the
  argument \code{alphanu}, the value of \code{bee} is the result of the inner
  optimization (\eqn{b^*}{bstar} in the notation in details),
  not the argument \code{bee} of this function.
}
\examples{
data(radish)

pred <- c(0,1,2)
fam <- c(1,3,2)

rout <- reaster(resp ~ varb + fit : (Site * Region),
    list(block = ~ 0 + fit : Block, pop = ~ 0 + fit : Pop),
    pred, fam, varb, id, root, data = radish)

alpha.mle <- rout$alpha
bee.mle <- rout$b
nu.mle <- rout$sigma^2
zwz.mle <- rout$zwz
obj <- rout$obj
fixed <- rout$fixed
random <- rout$random
alphanu.mle <- c(alpha.mle, nu.mle)

qout <- quickle(alphanu.mle, bee.mle, fixed, random, obj,
    zwz = zwz.mle, deriv = 2)
}
\keyword{misc}
back to top