https://github.com/cran/aster
Tip revision: aa47935123bfca8a22cbc8345d658d0c1713a289 authored by Charles J. Geyer on 14 December 2023, 15:20:02 UTC
version 1.1-3
version 1.1-3
Tip revision: aa47935
pickle.R
### minus the PQL criterion for estimating the variance components is
###
### pout$value + log(det(I + H)) / 2
###
### where pout is the object returned by the penmlogl function, and where
###
### H = D^(1 / 2) t(Z) W Z D^(1 / 2)
###
### (see equation (27) of inst/doc/re.tex) and I is the identity matrix.
### H is the random effects block of
###
### pout$mlogl.hessian
###
### If lambda is the vector of eigenvalues of H, then sum(log1p(lambda))
### calculates log(det(H + I))
checkargs <- function(sigma, parm, fixed, random, obj, y, origin, cache, ...)
{
stopifnot(inherits(obj, "aster"))
if (! missing(y)) {
stopifnot(is.matrix(y))
stopifnot(is.numeric(y))
stopifnot(is.finite(y))
stopifnot(dim(y) == dim(obj$x))
}
if (! missing(origin)) {
stopifnot(is.matrix(origin))
stopifnot(is.numeric(origin))
stopifnot(is.finite(origin))
stopifnot(dim(origin) == dim(obj$origin))
}
if (! missing(cache))
stopifnot(is.environment(cache))
stopifnot(is.matrix(fixed))
stopifnot(is.numeric(fixed))
stopifnot(is.finite(fixed))
nfix <- ncol(fixed)
stopifnot(is.matrix(random) | is.list(random))
if (! is.list(random))
random <- list(random)
for (i in seq(along = random)) {
foo <- random[[i]]
if (! is.matrix(foo))
stop("random not matrix or list of matrices")
if (nrow(foo) != nrow(fixed))
stop("fixed and random effects model matrices with different nrow")
if (! all(is.finite(foo)))
stop("some random effects model matrix not all finite")
}
nrand <- sapply(random, ncol)
stopifnot(is.numeric(sigma))
stopifnot(is.finite(sigma))
if (length(sigma) != length(nrand))
stop("length(sigma) != number of random effects model matrices")
if ((! missing(cache)) && exists("parm", envir = cache, inherits = FALSE)) {
stopifnot(is.numeric(cache$parm))
stopifnot(is.finite(cache$parm))
if (length(cache$parm) != nfix + sum(nrand))
stop("length(cache$parm) != total number of effects (fixed and random)")
} else {
stopifnot(is.numeric(parm))
stopifnot(is.finite(parm))
if (length(parm) != nfix + sum(nrand))
stop("length(parm) != total number of effects (fixed and random)")
}
}
pickle <- function(sigma, parm, fixed, random, obj, y, origin, cache, ...)
{
checkargs(sigma, parm, fixed, random, obj, y, origin, cache, ...)
if (missing(y)) y <- obj$x
if ((! missing(cache)) && exists("parm", envir = cache, inherits = FALSE))
parm <- cache$parm
trustargs <- list(objfun = penmlogl, parinit = parm, sigma = sigma,
fixed = fixed, random = random, obj = obj, y = y, ...)
if (is.null(trustargs$rinit)) trustargs$rinit <- 1
if (is.null(trustargs$rmax)) trustargs$rmax <- 10
if (! missing(origin)) trustargs$origin <- origin
trustargs$cache <- NULL
tout <- do.call(trust, trustargs)
stopifnot(tout$converged)
if ((! missing(cache))) cache$parm <- tout$argument
nfix <- ncol(fixed)
bigh <- tout$hessian
bigh <- bigh[1:nrow(bigh) > nfix, , drop = FALSE]
bigh <- bigh[ , 1:ncol(bigh) > nfix, drop = FALSE]
bigh <- bigh + diag(ncol(bigh))
baz <- diag(chol(bigh))
return(tout$value + sum(log(baz)))
}
### like above except now t(Z) W Z is fixed matrix supplied in argument zwz
### should be symmetric and positive semidefinite
### suitable zwz can be obtained from makezwz
makezwz <- function(sigma, parm, fixed, random, obj, y, origin)
{
checkargs(sigma, parm, fixed, random, obj, y, origin)
nfix <- ncol(fixed)
if (! is.list(random)) random <- list(random)
nrand <- sapply(random, ncol)
modmat <- fixed
for (i in seq(along = random))
modmat <- cbind(modmat, random[[i]], deparse.level = 0)
scalevec <- rep(c(1, sigma), times = c(nfix, nrand))
if (missing(y)) y <- obj$x
bar <- mlogl(scalevec * parm, obj$pred, obj$fam, y,
obj$root, modmat, deriv = 2, famlist = obj$famlist, origin = origin)
zwz <- bar$hessian
zwz <- zwz[1:nrow(zwz) > nfix, , drop = FALSE]
zwz <- zwz[ , 1:ncol(zwz) > nfix, drop = FALSE]
return(zwz)
}
pickle1 <- function(sigma, parm, fixed, random, obj, y, origin,
cache, zwz, deriv = 0, ...)
{
checkargs3part1(fixed, random, obj, y, origin)
nfix <- ncol(fixed)
if (! is.list(random)) random <- list(random)
nrand <- sapply(random, ncol)
if (! missing(cache)) {
stopifnot(is.environment(cache))
if (exists("parm", envir = cache, inherits = FALSE))
parm <- cache$parm
}
stopifnot(is.vector(parm))
if(length(parm) != nfix + sum(nrand))
stop("length(parm) != number of fixed effects + number of random effects")
idx <- seq(along = parm)
alpha <- parm[idx <= nfix]
cee <- parm[idx > nfix]
checkargs3part2(alpha, cee, sigma, nfix, nrand)
stopifnot(is.matrix(zwz))
stopifnot(is.numeric(zwz))
stopifnot(is.finite(zwz))
if (any(dim(zwz) != sum(nrand)))
stop("zwz not square matrix with dimension = number of random effects")
stopifnot(length(deriv) == 1)
stopifnot(deriv %in% c(0, 1))
if (missing(y)) y <- obj$x
trustargs <- list(objfun = penmlogl, parinit = parm, sigma = sigma,
fixed = fixed, random = random, obj = obj, y = y, ...)
if (is.null(trustargs$rinit)) trustargs$rinit <- 1
if (is.null(trustargs$rmax)) trustargs$rmax <- 10
if (! missing(origin)) trustargs$origin <- origin
trustargs$cache <- NULL
tout <- do.call(trust, trustargs)
stopifnot(tout$converged)
if ((! missing(cache))) cache$parm <- tout$argument
alpha <- tout$argument[idx <= nfix]
cee <- tout$argument[idx > nfix]
pout <- pickleHelper(alpha, cee, sigma, fixed, random, obj, y,
origin, zwz, deriv)
if (deriv == 0)
return(pout)
return(list(value = pout$value, gradient = pout$gradient$pt))
}
# note: we have implementation of deriv = 2, but it doesn't work because
# of inexactness of computer arithmetic.
# thus we only allow deriv = 0 or deriv = 1
pickle2 <- function(alphasigma, parm, fixed, random, obj, y, origin,
cache, zwz, deriv = 0, ...)
{
checkargs3part1(fixed, random, obj, y, origin)
nfix <- ncol(fixed)
if (! is.list(random)) random <- list(random)
nrand <- sapply(random, ncol)
if (! missing(cache)) {
stopifnot(is.environment(cache))
if (exists("parm", envir = cache, inherits = FALSE))
parm <- cache$parm
}
stopifnot(is.vector(alphasigma))
if(length(alphasigma) != nfix + length(nrand))
stop("length(alphasigma) != number of fixed effects + number of variance components")
idx <- seq(along = alphasigma)
alpha <- alphasigma[idx <= nfix]
cee <- parm
sigma <- alphasigma[idx > nfix]
checkargs3part2(alpha, cee, sigma, nfix, nrand)
stopifnot(is.matrix(zwz))
stopifnot(is.numeric(zwz))
stopifnot(is.finite(zwz))
if (any(dim(zwz) != sum(nrand)))
stop("zwz not square matrix with dimension = number of random effects")
stopifnot(length(deriv) == 1)
stopifnot(deriv %in% c(0, 1))
if (missing(y)) y <- obj$x
trustargs <- list(objfun = penmlogl2, parinit = parm, alpha = alpha,
sigma = sigma, fixed = fixed, random = random, obj = obj, y = y, ...)
if (is.null(trustargs$rinit)) trustargs$rinit <- 1
if (is.null(trustargs$rmax)) trustargs$rmax <- 10
if (! missing(origin)) trustargs$origin <- origin
trustargs$cache <- NULL
tout <- do.call(trust, trustargs)
stopifnot(tout$converged)
if ((! missing(cache))) cache$parm <- tout$argument
cee <- tout$argument
pout <- pickleHelper(alpha, cee, sigma, fixed, random, obj, y,
origin, zwz, deriv)
if (deriv == 0)
return(pout)
grad <- c(pout$gradient$pa, pout$gradient$pt)
if (deriv == 1)
return(list(value = pout$value, gradient = grad))
paa <- pout$hessian$paa
pac <- pout$hessian$pac
pat <- pout$hessian$pat
pcc <- pout$hessian$pcc
pct <- pout$hessian$pct
ptt <- pout$hessian$ptt
pcc.inv <- chol2inv(chol(pcc))
qaa <- paa - pac %*% pcc.inv %*% t(pac)
qat <- pat - pac %*% pcc.inv %*% pct
qtt <- ptt - t(pct) %*% pcc.inv %*% pct
hess <- rbind(cbind(qaa, qat), cbind(t(qat), qtt))
return(list(value = pout$value, gradient = grad, hessian = hess))
}