##' Calculate one-step-ahead (OSA) residuals for a latent variable
##' model. (\emph{Modified from TMB version to allow easier parallel
##' computations})
##'
##' @title Calculate one-step-ahead (OSA) residuals for a latent variable model.
##' @param obj Output from \code{MakeADFun}.
##' @param observation.name Character naming the observation in the template.
##' @param data.term.indicator Character naming an indicator data variable in the template (not required by all methods - see details).
##' @param method Method to calculate OSA (see details).
##' @param subset Index vector of observations that will be added one by one during OSA. By default \code{1:length(observations)} (with \code{conditional} subtracted).
##' @param conditional Index vector of observations that are fixed during OSA. By default the empty set.
##' @param discrete Are observations discrete? (assumed FALSE by default)
##' @param discreteSupport Possible outcomes of discrete distribution (\code{method="oneStepGeneric"} only).
##' @param range Possible range of the observations.
##' @param seed Randomization seed (discrete case only). If \code{NULL} the RNG seed is untouched by this routine.
##' @param parallel Run in parallel using the \code{parallel} package?
##' @param ncores number of threads to run in parallel. Ignored if parallel = FALSE
##' @param trace Trace progress?
##' @param reverse Do calculations in opposite order to improve stability ? (currently enabled by default for \code{oneStepGaussianOffMode} method only)
##' @param ... Control parameters for OSA method
##' @return \code{data.frame} with OSA \emph{standardized} residuals
##' in column \code{residual}. Depending on the method the output may
##' also include OSA expected observation in column \code{mean}.
##'
##' @importFrom parallel detectCores mclapply
##' @importFrom stats splinefun runif qnorm optimHess
##' @importFrom graphics abline
##' @keywords internal
TMBoneStepPredict <- function(obj,
## Names of data objects (not all are optional)
observation.name = NULL,
data.term.indicator = NULL,
method=c(
"oneStepGaussianOffMode",
"fullGaussian",
"oneStepGeneric",
"oneStepGaussian",
"cdf"),
subset = NULL,
conditional = NULL,
discrete = NULL,
discreteSupport = NULL,
range = c(-Inf, Inf),
seed = 123,
parallel = FALSE,
ncores = 2,
trace = TRUE,
reverse = (method == "oneStepGaussianOffMode"),
...
){
if (missing(observation.name))
stop("'observation.name' must define a data component")
if (!(observation.name %in% names(obj$env$data)))
stop("'observation.name' must be in data component")
method <- match.arg(method)
if (is.null(data.term.indicator)){
if(method != "fullGaussian"){
stop(paste0("method='",method,"' requires a 'data.term.indicator'"))
}
}
if (!missing(discreteSupport) && !missing(range))
stop("Cannot specify both 'discreteSupport' and 'range'")
obs <- as.vector(obj$env$data[[observation.name]])
if(is.null(discrete)){
ndup <- sum(duplicated(obs))
if(ndup > 0){
warning("Observations do not look continuous. Number of duplicates = ", ndup)
stop("Argument 'discrete' (TRUE/FALSE) must be specified.")
}
discrete <- FALSE
} else {
stopifnot(is.logical(discrete))
}
## Using wrong method for discrete data ?
if (discrete){
if (! (method %in% c("oneStepGeneric", "cdf")) ){
stop(paste0("method='",method,"' is not for discrete observations."))
}
}
## Default subset/permutation:
if(is.null(subset)){
subset <- 1:length(obs)
subset <- setdiff(subset, conditional)
}
## Check
if(!is.null(conditional)){
if(length(intersect(subset, conditional)) > 0){
stop("'subset' and 'conditional' have non-empty intersection")
}
}
unconditional <- setdiff(1:length(obs), union(subset, conditional))
## Args to construct copy of 'obj'
args <- as.list(obj$env)[intersect(names(formals(MakeADFun)), ls(obj$env))]
## Use the best encountered parameter for new object
if(length(obj$env$random))
args$parameters <- obj$env$parList(par = obj$env$last.par.best)
else
args$parameters <- obj$env$parList(obj$env$last.par.best)
## Fix all non-random components of parameter list
names.random <- unique(names(obj$env$par[obj$env$random]))
names.all <- names(args$parameters)
fix <- setdiff(names.all, names.random)
map <- lapply(args$parameters[fix], function(x)factor(x*NA))
ran.in.map <- names.random[names.random %in% names(args$map)]
if(length(ran.in.map)) map <- c(map, args$map[ran.in.map]) # don't overwrite random effects mapping
args$map <- map ## Overwrite map
## Find randomeffects character
args$random <- names.random
args$regexp <- FALSE
## Move data$name to parameter$name
args$parameters[observation.name] <- args$data[observation.name]
args$data[observation.name] <- NULL
## Make data.term.indicator in parameter list
if(!is.null(data.term.indicator)){
one <- rep(1, length(obs))
zero <- rep(0, length(obs))
if(method=="cdf"){
args$parameters[[data.term.indicator]] <- cbind(one, zero, zero)
} else {
args$parameters[[data.term.indicator]] <- cbind(one)
}
}
## Pretend these are *not observed*:
if(length(unconditional)>0){
if(is.null(data.term.indicator))
stop("Failed to disable some data terms (because 'data.term.indicator' missing)")
args$parameters[[data.term.indicator]][unconditional, 1] <- 0
}
## Pretend these are *observed*:
if(length(conditional)>0){
if(is.null(data.term.indicator))
stop("Failed to enable some data terms (because 'data.term.indicator' missing)")
args$parameters[[data.term.indicator]][conditional, 1] <- 1
}
## Make map for observations and indicator variables:
makeFac <- function(x){
fac <- as.matrix(x)
fac[] <- 1:length(x)
fac[conditional, ] <- NA
fac[unconditional, ] <- NA
fac[subset, ] <- 1:(length(subset)*ncol(fac)) ## Permutation
factor(fac)
}
map <- list()
map[[observation.name]] <- makeFac(obs)
if(!is.null(data.term.indicator)){
map[[data.term.indicator]] <- makeFac(args$parameters[[data.term.indicator]])
}
args$map <- c(args$map, map)
## New object be silent
args$silent <- TRUE
## Create new object
newobj <- do.call("MakeADFun", args)
## Helper function to loop through observations:
nm <- names(newobj$par)
obs.pointer <- which(nm == observation.name)
if(method=="cdf"){
tmp <- matrix(which(nm == data.term.indicator), ncol=3)
data.term.pointer <- tmp[,1]
lower.cdf.pointer <- tmp[,2]
upper.cdf.pointer <- tmp[,3]
} else {
data.term.pointer <- which(nm == data.term.indicator)
lower.cdf.pointer <- NULL
upper.cdf.pointer <- NULL
}
observation <- local({
obs.local <- newobj$par
i <- 1:length(subset)
function(k, y=NULL, lower.cdf=FALSE, upper.cdf=FALSE){
## Disable all observations later than k:
obs.local[data.term.pointer[k<i]] <- 0
## On request, overwrite k'th observation:
if(!is.null(y)) obs.local[obs.pointer[k]] <- y
## On request, get tail probs rather than point probs:
if(lower.cdf | upper.cdf){
obs.local[data.term.pointer[k]] <- 0
if(lower.cdf) obs.local[lower.cdf.pointer[k]] <- 1
if(upper.cdf) obs.local[upper.cdf.pointer[k]] <- 1
}
obs.local
}
})
## Parallel case: overload lapply
if(parallel){
## mclapply uses fork => must set nthreads=1
## modified from TMB to ensure correct DLL selected (in case glmmTMB DLL is loaded)
nthreads.restore <- TMB::openmp(DLL = "mpmm")
on.exit( TMB::openmp( nthreads.restore, DLL = "mpmm" ), add=TRUE)
# TMB::openmp(ncores)
# requireNamespace("parallel") # was library(parallel)
# lapply <- parallel::mclapply
}
## Trace one-step functions
tracefun <- function(k)if(trace)print(k)
## Apply a one-step method and generate common output assuming
## the method generates at least:
## * nll
## * nlcdf.lower
## * nlcdf.upper
applyMethod <- function(oneStepMethod){
ord <- seq_along(subset)
if (reverse) ord <- rev(ord)
if(parallel) pred <- do.call("rbind", mclapply(ord, oneStepMethod, mc.cores = ncores))
else pred <- do.call("rbind", lapply(ord, oneStepMethod))
pred <- as.data.frame(pred)[ord, ]
pred$Fx <- 1 / ( 1 + exp(pred$nlcdf.lower - pred$nlcdf.upper) )
pred$px <- 1 / ( exp(-pred$nlcdf.lower + pred$nll) +
exp(-pred$nlcdf.upper + pred$nll) )
if(discrete){
if(!is.null(seed)){
## Restore RNG on exit:
Random.seed <- .GlobalEnv$.Random.seed
on.exit(.GlobalEnv$.Random.seed <- Random.seed)
set.seed(seed)
}
U <- runif(nrow(pred))
} else {
U <- 0
}
pred$residual <- qnorm(pred$Fx - U * pred$px)
pred
}
## ######################### CASE: oneStepGaussian
if(method == "oneStepGaussian"){
p <- newobj$par
newobj$fn(p) ## Test eval
oneStepGaussian <- function(k){
tracefun(k)
index <- subset[k]
f <- function(y){
newobj$fn(observation(k, y))
}
g <- function(y){
newobj$gr(observation(k, y))[obs.pointer[k]]
}
opt <- nlminb(obs[index], f, g)
H <- optimHess(opt$par, f, g)
c(observation=obs[index], mean=opt$par, sd=sqrt(1/H))
}
ord <- seq_along(subset)
if (reverse) ord <- rev(ord)
if(parallel) pred <- do.call("rbind", mclapply(ord, oneStepGaussian, mc.cores = ncores))
else pred <- do.call("rbind", lapply(ord, oneStepGaussian))
pred <- as.data.frame(pred)[ord, ]
pred$residual <- (pred$observation-pred$mean)/pred$sd
}
## ######################### CASE: oneStepGaussianOffMode
if(method == "oneStepGaussianOffMode"){
p <- newobj$par
newobj$fn(p) ## Test eval
newobj$env$random.start <- expression({last.par[random]})
oneStepGaussian <- function(k){
tracefun(k)
index <- subset[k]
f <- function(y){
newobj$fn(observation(k, y))
}
g <- function(y){
newobj$gr(observation(k, y))[obs.pointer[k]]
}
c(observation=obs[index], nll = f(obs[index]), grad = g(obs[index]))
}
ord <- seq_along(subset)
if (reverse) ord <- rev(ord)
if(parallel) pred <- do.call("rbind", mclapply(ord, oneStepGaussian, mc.cores = ncores))
else pred <- do.call("rbind", lapply(ord, oneStepGaussian))
pred <- as.data.frame(pred)[ord, ]
################### Convert value and gradient to residual
## Need Lambert W function: x = W(x) * exp( W(x) ) , x > 0
## Vectorized in x and tested on extreme cases W(.Machine$double.xmin)
## and W(.Machine$double.xmax).
W <- function(x){
## Newton: f(y) = y * exp(y) - x
## f'(y) = y * exp(y) + exp(y)
rel.tol <- sqrt(.Machine$double.eps)
logx <- log(x)
fdivg <- function(y)(y - exp(logx - y)) / (1 + y)
y <- pmax(logx, 0)
while( any( abs( logx - log(y) - y) > rel.tol, na.rm=TRUE) ) {
y <- y - fdivg(y)
}
y
}
getResid <- function(value, grad){
Rabs <- sqrt( W( exp( 2*(value - log(sqrt(2*pi)) + log(abs(grad))) ) ) )
R <- sign(grad) * Rabs
R
}
nll0 <- newobj$fn(observation(0))
R <- getResid( diff( c(nll0, pred$nll) ), pred$grad )
M <- pred$observation - ifelse(pred$grad != 0, R * (R / pred$grad), 0)
pred$mean <- M
pred$residual <- R
}
## ######################### CASE: oneStepGeneric
if((method == "oneStepGeneric") && missing(discreteSupport)){
p <- newobj$par
newobj$fn(p) ## Test eval
newobj$env$value.best <- -Inf ## <-- Never overwrite last.par.best
nan2zero <- function(x)if(!is.finite(x)) 0 else x
## Set default configuration for this method (modify with '...'):
formals(tmbprofile)$ytol <- 10 ## Tail tolerance (increase => more tail)
formals(tmbprofile)$ystep <- .5 ## Grid spacing (decrease => more accuracy)
## Handle discrete case
if(discrete){
formals(tmbprofile)$h <- 1
integrate <- function(f, lower, upper, ...){
grid <- ceiling(lower):floor(upper)
list( value = sum( f(grid) ) )
}
}
oneStepGeneric <- function(k){
tracefun(k)
ans <- try({
index <- subset[k]
f <- function(y){
newobj$fn(observation(k, y))
}
nll <- f(obs[index]) ## Marginal negative log-likelihood
newobj$env$last.par.best <- newobj$env$last.par ## <-- used by tmbprofile
slice <- tmbprofile(newobj, k, slice=TRUE,
parm.range = range,...)
spline <- splinefun(slice[[1]], slice[[2]])
spline.range <- range(slice[[1]])
if(trace >= 2){
plotfun <- function(slice, spline){
plot(slice, type="p", level=NULL)
plot(spline, spline.range[1], spline.range[2], add=TRUE)
abline(v=obs[index], lty="dashed")
}
if(trace >= 3){
slice$value <- exp( -(slice$value - nll) )
plotfun(slice, function(x)exp(-(spline(x) - nll)))
}
else
plotfun(slice, spline)
}
F1 <- integrate(function(x)exp(-(spline(x) - nll)),
spline.range[1],
obs[index])$value
F2 <- integrate(function(x)exp(-(spline(x) - nll)),
obs[index] + discrete,
spline.range[2])$value
mean <- integrate(function(x)exp(-(spline(x) - nll)) * x,
spline.range[1],
spline.range[2])$value / (F1 + F2)
## Was:
## F1 <- integrate(Vectorize( function(x)nan2zero( exp(-(f(x) - nll)) ) ), -Inf, obs[index])$value
## F2 <- integrate(Vectorize( function(x)nan2zero( exp(-(f(x) - nll)) ) ), obs[index], Inf)$value
nlcdf.lower = nll - log(F1)
nlcdf.upper = nll - log(F2)
c(nll=nll, nlcdf.lower=nlcdf.lower, nlcdf.upper=nlcdf.upper, mean=mean)
})
if(is(ans, "try-error")) ans <- NaN
ans
}
pred <- applyMethod(oneStepGeneric)
}
## ######################### CASE: oneStepDiscrete
if((method == "oneStepGeneric") && !missing(discreteSupport)){
p <- newobj$par
newobj$fn(p) ## Test eval
obs <- as.integer(round(obs))
if(is.null(discreteSupport)){
warning("Setting 'discreteSupport' to ",min(obs),":",max(obs))
discreteSupport <- min(obs):max(obs)
}
oneStepDiscrete <- function(k){
tracefun(k)
ans <- try({
index <- subset[k]
f <- function(y){
newobj$fn(observation(k, y))
}
nll <- f(obs[index]) ## Marginal negative log-likelihood
F <- Vectorize(function(x)exp(-(f(x) - nll))) (discreteSupport)
F1 <- sum( F[discreteSupport <= obs[index]] )
F2 <- sum( F[discreteSupport > obs[index]] )
nlcdf.lower = nll - log(F1)
nlcdf.upper = nll - log(F2)
c(nll=nll, nlcdf.lower=nlcdf.lower, nlcdf.upper=nlcdf.upper)
})
if(is(ans, "try-error")) ans <- NaN
ans
}
pred <- applyMethod(oneStepDiscrete)
}
## ######################### CASE: fullGaussian
if(method == "fullGaussian"){
## Same object with y random:
args2 <- args
args2$random <- c(args2$random, observation.name)
## Change map: Fix everything except observations
fix <- data.term.indicator
args2$map[fix] <- lapply(args2$map[fix],
function(x)factor(NA*unclass(x)))
newobj2 <- do.call("MakeADFun", args2)
newobj2$fn() ## Test-eval to find mode
mode <- newobj2$env$last.par
GMRFmarginal <- function (Q, i, ...) {
ind <- 1:nrow(Q)
i1 <- (ind)[i]
i0 <- setdiff(ind, i1)
if (length(i0) == 0)
return(Q)
Q0 <- as(Q[i0, i0, drop = FALSE], "symmetricMatrix")
L0 <- Cholesky(Q0, ...)
ans <- Q[i1, i1, drop = FALSE] - Q[i1, i0, drop = FALSE] %*%
solve(Q0, Q[i0, i1, drop = FALSE])
ans
}
h <- newobj2$env$spHess(mode, random=TRUE)
i <- which(names(newobj2$env$par[newobj2$env$random]) == observation.name)
Sigma <- solve( as.matrix( GMRFmarginal(h, i) ) )
res <- obs[subset] - mode[i]
L <- t(chol(Sigma))
pred <- data.frame(residual = as.vector(solve(L, res)))
}
## ######################### CASE: cdf
if(method == "cdf"){
p <- newobj$par
newobj$fn(p) ## Test eval
cdf <- function(k){
tracefun(k)
nll <- newobj$fn(observation(k))
nlcdf.lower <- newobj$fn(observation(k, lower.cdf = TRUE))
nlcdf.upper <- newobj$fn(observation(k, upper.cdf = TRUE))
c(nll=nll, nlcdf.lower=nlcdf.lower, nlcdf.upper=nlcdf.upper)
}
pred <- applyMethod(cdf)
}
pred
}
## Goodness of fit residuals based on an approximate posterior
## sample. (\emph{Beta version; may change without notice})
##
## Denote by \eqn{(u, x)} the pair of the true un-observed random effect
## and the data. Let a model specification be given in terms of the
## estimated parameter vector \eqn{\theta} and let \eqn{u^*} be a
## sample from the conditional distribution of \eqn{u} given
## \eqn{x}. If the model specification is correct, it follows that the
## distribution of the pair \eqn{(u^*, x)} is the same as the distribution
## of \eqn{(u, x)}. Goodness-of-fit can thus be assessed by proceeding as
## if the random effect vector were observed, i.e check that \eqn{u^*}
## is consistent with prior model of the random effect and that \eqn{x}
## given \eqn{u^*} agrees with the observation model.
##
## This function can carry out the above procedure for many TMB models
## under the assumption that the true posterior is well approximated by a
## Gaussian distribution.
##
## First a draw from the Gaussian posterior distribution \eqn{u^*} is
## obtained based on the mode and Hessian of the random effects given the
## data.
## This sample uses sparsity of the Hessian and will thus work for large systems.
##
## An automatic standardization of the sample can be carried out \emph{if
## the observation model is Gaussian} (\code{fullGaussian=TRUE}). In this
## case the prior model is obtained by disabling the data term and
## calculating mode and Hessian. A \code{data.term.indicator} must be
## given in order for this to work. Standardization is performed using
## the sparse Cholesky of the prior precision.
## By default, this step does not use a fill reduction permutation \code{perm=FALSE}.
## This is often superior wrt. to interpretation of the.
## the natural order of the parameter vector is used \code{perm=FALSE}
## which may be superior wrt. to interpretation. Otherwise
## \code{perm=TRUE} a fill-reducing permutation is used while
## standardizing.
## @references Waagepetersen, R. (2006). A Simulation-based Goodness-of-fit Test for Random Effects in Generalized Linear Mixed Models. Scandinavian journal of statistics, 33(4), 721-731.
## @param obj TMB model object from \code{MakeADFun}.
## @param observation.name Character naming the observation in the template.
## @param data.term.indicator Character naming an indicator data variable in the template. Only used if \code{standardize=TRUE}.
## @param standardize Logical; Standardize sample with the prior covariance ? Assumes all latent variables are Gaussian.
## @param as.list Output posterior sample, and the corresponding standardized residual, as a parameter list ?
## @param perm Logical; Use a fill-reducing ordering when standardizing ?
## @param fullGaussian Logical; Flag to signify that the joint distribution of random effects and data is Gaussian.
## @return List with components \code{sample} and \code{residual}.
oneSamplePosterior <- function(obj,
observation.name = NULL,
data.term.indicator = NULL,
standardize = TRUE,
as.list = TRUE,
perm = FALSE,
fullGaussian = FALSE){
## Draw Gaussian posterior sample
tmp <- obj$env$MC(n=1, keep=TRUE, antithetic=FALSE)
samp <- as.vector( attr(tmp, "samples") )
## If standardize
resid <- NULL
if (standardize) {
## Args to construct copy of 'obj'
args <- as.list(obj$env)[intersect(names(formals(MakeADFun)), ls(obj$env))]
## Use the best encountered parameter for new object
args$parameters <- obj$env$parList(par = obj$env$last.par.best)
## Make data.term.indicator in parameter list
obs <- obj$env$data[[observation.name]]
nobs <- length(obs)
zero <- rep(0, nobs)
if ( ! fullGaussian )
args$parameters[[data.term.indicator]] <- zero
## Fix all non-random components of parameter list
names.random <- unique(names(obj$env$par[obj$env$random]))
names.all <- names(args$parameters)
fix <- setdiff(names.all, names.random)
map <- lapply(args$parameters[fix], function(x)factor(x*NA))
args$map <- map ## Overwrite map
## If 'fullGaussian == TRUE' turn 'obs' into a random effect
if (fullGaussian) {
names.random <- c(names.random, observation.name)
args$parameters[[observation.name]] <- obs
}
## Find randomeffects character
args$random <- names.random
args$regexp <- FALSE
## New object be silent
args$silent <- TRUE
## Create new object
newobj <- do.call("MakeADFun", args)
## Construct Hessian and Cholesky
newobj$fn()
## Get Cholesky and prior mean
## FIXME: We are using the mode as mean. Consider skewness
## correction similar to 'bias.correct' in 'sdreport'.
L <- newobj$env$L.created.by.newton
mu <- newobj$env$last.par
## If perm == FALSE redo Cholesky with natural ordering
if ( ! perm ) {
Q <- newobj$env$spHess(mu, random=TRUE)
L <- Matrix::Cholesky(Q, super=TRUE, perm=FALSE)
}
## If 'fullGaussian == TRUE' add 'obs' to the sample
if (fullGaussian) {
tmp <- newobj$env$par * NA
tmp[names(tmp) == observation.name] <- obs
tmp[names(tmp) != observation.name] <- samp
samp <- tmp
}
## Standardize ( P * Q * P^T = L * L^T )
r <- samp - mu
rp <- r[L@perm + 1]
Lt <- Matrix::t(
as(L, "sparseMatrix")
)
resid <- as.vector( Lt %*% rp )
}
if (as.list) {
if (standardize) obj <- newobj
par <- obj$env$last.par.best
asList <- function(samp) {
par[obj$env$random] <- samp
samp <- obj$env$parList(par=par)
nm <- unique(names(obj$env$par[obj$env$random]))
samp[nm]
}
samp <- asList(samp)
if (!is.null(resid))
resid <- asList(resid)
}
ans <- list()
ans$sample <- samp
ans$residual <- resid
ans
}
if(FALSE) {
library(TMB)
runExample("MVRandomWalkValidation", exfolder="../../tmb_examples/validation")
set.seed(1)
system.time( qw <- TMB:::oneSamplePosterior(obj, "obs", "keep") )
qqnorm(as.vector(qw$residual$u)); abline(0,1)
runExample("rickervalidation", exfolder="../../tmb_examples/validation")
set.seed(1)
system.time( qw <- TMB:::oneSamplePosterior(obj, "Y", "keep") )
qqnorm(as.vector(qw$residual$X)); abline(0,1)
runExample("ar1xar1")
set.seed(1)
system.time( qw <- TMB:::oneSamplePosterior(obj, "N", "keep") )
qqnorm(as.vector(qw$residual$eta)); abline(0,1)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.