Nothing
## backward compat (copied from lme4)
if((Rv <- getRversion()) < "3.2.1") {
lengths <- function (x, use.names = TRUE) vapply(x, length, 1L, USE.NAMES = use.names)
}
rm(Rv)
## generate a list with names equal to values
## See also: \code{tibble::lst}, \code{Hmisc::llist}
namedList <- function (...) {
L <- list(...)
snm <- sapply(substitute(list(...)), deparse)[-1]
if (is.null(nm <- names(L)))
nm <- snm
if (any(nonames <- nm == ""))
nm[nonames] <- snm[nonames]
setNames(L, nm)
}
## Sparse Schur complement (Marginal of precision matrix)
##' @importFrom Matrix Cholesky solve
GMRFmarginal <- function(Q, i, ...) {
ind <- seq_len(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, as.matrix(Q[i0, i1, drop = FALSE]))
ans
}
parallel_default <- function(parallel=c("no","multicore","snow"),ncpus=1) {
## boilerplate parallel-handling stuff, copied from lme4
if (missing(parallel)) parallel <- getOption("profile.parallel", "no")
parallel <- match.arg(parallel)
do_parallel <- (parallel != "no" && ncpus > 1L)
if (do_parallel && parallel == "multicore" &&
.Platform$OS.type == "windows") {
warning("no multicore on Windows, falling back to non-parallel")
parallel <- "no"
}
return(list(parallel=parallel,do_parallel=do_parallel))
}
## from length of (strictly) lower triangle, compute dimension of matrix
get_matdim <- function(ntri) {
as.integer(round(0.5 * (1 + sqrt(1 + 8 * ntri))))
}
##' transform correlation parameters to and from glmmTMB parameterization
##' @param theta vector of internal correlation parameters (elements of scaled Cholesky factor, in \emph{row-major} order)
##' @param return_val return a vector of correlation values from the lower triangle ("vec"), or the full correlation matrix ("mat")?
##' @return a vector of correlation values (\code{get_cor}) or glmmTMB scaled-correlation parameters (\code{put_cor})
##' @details
##' \code{\link{get_cor}} transforms from the glmmTMB parameterization (components of a \code{theta} parameter vector) to correlations;
##' \code{\link{put_cor}} does the inverse transformations, from correlations to \code{theta} values.
##'
##' These functions follow the definition at \url{http://kaskr.github.io/adcomp/classdensity_1_1UNSTRUCTURED__CORR__t.html}:
##' if \eqn{L} is the lower-triangular matrix with 1 on the diagonal and the correlation parameters in the lower triangle, then the correlation matrix is defined as \eqn{\Sigma = D^{-1/2} L L^\top D^{-1/2}}{Sigma = sqrt(D) L L' sqrt(D)}, where \eqn{D = \textrm{diag}(L L^\top)}{D = diag(L L')}. For a single correlation parameter \eqn{\theta_0}{theta0} (i.e. the correlation in a 2x2 correlation matrix), this works out to \eqn{\rho = \theta_0/\sqrt{1+\theta_0^2}}{rho = theta0/sqrt(1+theta0^2)}. The \code{get_cor} function returns the elements of the lower triangle of the correlation matrix, in column-major order.
##'
##' These functions also work for AR1 correlation parameters.
##' @examples
##' th0 <- 0.5
##' stopifnot(all.equal(get_cor(th0), th0/sqrt(1+th0^2)))
##' set.seed(101)
##' ## pick 6 values for a random 4x4 correlation matrix
##' print(C <- get_cor(rnorm(6), return_val = "mat"), digits = 3)
##' ## transform a correlation matrix to a theta vector
##' cor_mat <- matrix(c(1,0.3,0.1,
##' 0.3,1,0.2,
##' 0.1,0.2,1), ncol = 3)
##' put_cor(cor_mat, "mat")
##' put_cor(cor_mat[lower.tri(cor_mat)], "vec")
##' ## test: round-trip
##' stopifnot(all.equal(get_cor(put_cor(C), return_val = "mat"), C))
##' @export
get_cor <- function(theta, return_val = c("vec", "mat")) {
return_val <- match.arg(return_val)
n <- get_matdim(length(theta))
R <- diag(n)
R[upper.tri(R)] <- theta
R[] <- crossprod(R) # R <- t(R) %*% R
scale <- 1 / sqrt(diag(R))
R[] <- scale * R * rep(scale, each = n) # R <- cov2cor(R)
if (return_val == "mat") return(R)
return(R[lower.tri(R)])
}
##' @rdname get_cor
##' @param C a correlation matrix
##' @param input_val input a vector of correlation values from the lower triangle ("vec"), or the full correlation matrix ("mat")?
##' @export
put_cor <- function(C, input_val = c("mat", "vec")) {
input_val <- match.arg(input_val)
if (input_val == "vec") {
## construct matrix
M <- diag(get_matdim(length(C)))
M[lower.tri(M)] <- C
M[upper.tri(M)] <- t(M)[upper.tri(M)]
C <- M
}
cc2 <- chol(C)
scale <- diag(cc2)
cc2 <- cc2 %*% diag(1/scale)
cc2[upper.tri(cc2)]
}
hasRandom <- function(x) {
pl <- getParList(x)
return(length(unlist(pl[grep("^theta",names(pl))]))>0)
}
##' retrieve parameters by name or index
##' @param parm parameter specifier
##' @param object fitted glmmTMB object
##' @param full include simple dispersion parameter?
##' @param include_nonest include mapped parameter indices for non-estimated (mapped or rank-deficient/dropped) parameters?
##' @noRd
getParms <- function(parm=NULL, object, full=FALSE, include_nonest = FALSE) {
vv <- vcov(object, full=TRUE, include_nonest = include_nonest)
sds <- sqrt(diag(vv))
pnames <- names(sds) <- rownames(vv) ## parameter names (user-facing)
ee <- object$obj$env
## don't use object$obj$env$random; we want to keep "beta" vals, which may be
## counted as "random" if using REML
drop_rand <- function(x) x[!x %in% c("b", "bzi", "bdisp")]
if (!include_nonest) {
intnames <- drop_rand(names(ee$last.par))
} else {
pl <- ee$parList()
intnames <- drop_rand(rep(names(pl), lengths(pl)))
}
if (length(pnames) != length(sds)) { ## shouldn't happen ...
stop("length mismatch between internal and external parameter names")
}
if (is.null(parm)) {
if (!full && trivialDisp(object)) {
parm <- grep("betadisp", intnames, invert=TRUE)
} else {
parm <- seq_along(sds)
}
}
if (is.character(parm)) {
if (identical(parm,"psi_")) {
parm <- grep("^psi",intnames)
} else if (identical(parm,"theta_")) {
parm <- grep("^theta",intnames)
} else if (identical(parm,"beta_")) {
if (trivialDisp(object)) {
## include conditional and zi params
## but not dispersion params
parm <- grep("^beta(zi)?$",intnames)
} else {
parm <- grep("beta",intnames)
}
} else if (identical(parm, "disp_") ||
identical(parm, "sigma")) {
parm <- grep("^betadisp", intnames)
} else { ## generic parameter vector
nparm <- match(parm,pnames)
if (any(is.na(nparm))) {
stop("unrecognized parameter names: ",
parm[is.na(nparm)])
}
parm <- nparm
}
}
return(parm)
}
isREML <- function(x) {
if (is.null(REML <- x$modelInfo$REML)) {
## let vcov work with old (pre-REML option) stored objects
REML <- FALSE
}
return(REML)
}
## action: message, warning, stop
check_dots <- function(..., .ignore = NULL, .action="stop") {
L <- list(...)
if (length(.ignore)>0) {
L <- L[!names(L) %in% .ignore]
}
if (length(L)>0) {
FUN <- get(.action)
FUN("unknown arguments: ",
paste(names(L), collapse=","))
}
return(NULL)
}
if (getRversion()<"4.0.0") {
deparse1 <- function (expr, collapse = " ", width.cutoff = 500L, ...) {
paste(deparse(expr, width.cutoff, ...), collapse = collapse)
}
}
## in case these are useful, we can document and export them later ...
#' @importFrom stats rnbinom qnbinom dnbinom pnbinom
rnbinom1 <- function(n, mu, phi) {
## var = mu*(1+phi) = mu*(1+mu/k) -> k = mu/phi
rnbinom(n, mu=mu, size=mu/phi)
}
dnbinom1 <- function(x, mu, phi, ...) {
dnbinom(x, mu=mu, size=mu/phi, ...)
}
pnbinom1 <- function(q, mu, phi, ...) {
pnbinom(q, mu=mu, size=mu/phi, ...)
}
qnbinom1 <- function(p, mu, phi, ...) {
qnbinom(p, mu=mu, size=mu/phi, ...)
}
nullSparseMatrix <- function() {
argList <- list(
dims=c(0,0),
i=integer(0),
j=integer(0),
x=numeric(0))
if (utils::packageVersion("Matrix")<"1.3.0") {
do.call(Matrix::sparseMatrix, c(argList, list(giveCsparse=FALSE)))
} else {
do.call(Matrix::sparseMatrix, c(argList, list(repr="T")))
}
}
## simplified version of glmmTMB package checking
##' @param this_pkg downstream package being tested
##' @param dep_pkg upstream package on which \code{this_pkg} depends
##' @param dep_type "ABI" or "package"
##' @param built_version a \code{numeric_version} object indicating what version of \code{dep_pkg} was used to build \code{this_pkg}
##' @param warn (logical) warn if condition not met?
##' @importFrom utils packageVersion
##' @noRd
check_dep_version <- function(this_pkg = "glmmTMB", dep_pkg = "TMB", dep_type = "package",
built_version = .TMB.build.version,
warn = TRUE) {
## FIXME: replace by TMB.Version() when available ?
cur_version <- packageVersion(dep_pkg)
result_ok <- cur_version == built_version
if (!result_ok) {
warning(
sprintf("%s version mismatch: \n", dep_type),
sprintf("%s was built with %s %s version %s\n",
this_pkg, dep_pkg, dep_type, built_version),
sprintf("Current %s %s version is %s\n",
dep_pkg, dep_type, cur_version),
sprintf("Please re-install %s from source ", this_pkg),
"or restore original ",
sQuote(dep_pkg), " package",
" (see '?reinstalling' for more information)"
)
}
return(result_ok)
}
#' @name reinstalling
#' @rdname reinstalling
#' @title Reinstalling binary dependencies
#'
#' @description The \code{glmmTMB} package depends on several upstream packages, which it
#' uses in a way that depends heavily on their internal (binary) structure.
#' Sometimes, therefore, installing an update to one of these packages will
#' require that you re-install a \emph{binary-compatible} version of \code{glmmTMB},
#' i.e. a version that has been compiled with the updated version of the upstream
#' package.
#' \itemize{
#' \item If you have development tools (compilers etc.) installed, you
#' should be able to re-install a binary-compatible version of the package by running
#' \code{install.packages("glmmTMB", type="source")}. If you want to install
#' the development version of \code{glmmTMB} instead, you can use
#' \code{remotes::install_github("glmmTMB/glmmTMB/glmmTMB")}.
#' (On Windows, you can install development tools following the instructions at
#' \url{https://cran.r-project.org/bin/windows/Rtools/}; on MacOS, see
#' \url{https://mac.r-project.org/tools/}.)
#'
#' \item If you do \emph{not} have development tools and can't/don't want to
#' install them (and so can't install packages with compiled code from source),
#' you have two choices:
#' \itemize{
#' \item revert the upstream package(s) to their previous binary version. For example, using the
#' \code{checkpoint} package:
#' \preformatted{
#' ## load (installing if necessary) the checkpoint package
#' while (!require("checkpoint")) install.packages("checkpoint")
#' ## retrieve build date of installed version of glmmTMB
#' bd <- as.character(asDateBuilt(
#' packageDescription("glmmTMB",fields="Built")))
#' oldrepo <- getOption("repos")
#' use_mran_snapshot(bd) ## was setSnapshot() pre-checkpoint v1.0.0
#' install.packages("TMB")
#' options(repos=oldrepo) ## restore original repo
#' }
#' A similar recipe (substituting \code{Matrix} for \code{TMB} and \code{TMB} for \code{glmmTMB})
#' can be used if you get warnings about an incompatibility between \code{TMB} and \code{Matrix}.
#' \item hope that the glmmTMB maintainers have posted a binary
#' version of the package that works with your system; try installing it via
#' \code{install.packages("glmmTMB",repos="https://glmmTMB.github.io/glmmTMB/repos",type="binary")}
#' If this doesn't work, please file an issue (with full details about your
#' operating system and R version) asking the maintainers to build and
#' post an appropriate binary version of the package.
#' }
#' }
NULL
#' Check OpenMP status
##'
##' Checks whether OpenMP has been successfully enabled for this
##' installation of the package. (Use the \code{parallel} argument
##' to \code{\link{glmmTMBControl}}, or set \code{options(glmmTMB.cores=[value])},
##' to specify that computations should be done in parallel.) To further
##' trace OpenMP settings, use \code{options(glmmTMB_openmp_debug = TRUE)}.
##' @seealso \code{\link[TMB]{benchmark}}, \code{\link{glmmTMBControl}}
##' @return \code{TRUE} or \code{FALSE} depending on availability of OpenMP,
##' @aliases openmp
##' @export
omp_check <- function() {
.Call("omp_check", PACKAGE="glmmTMB")
}
get_pars <- function(object, unlist=TRUE) {
ee <- object$obj$env
x <- ee$last.par.best
## work around built-in default to parList, which
## is bad if no random component
if (length(ee$random)>0) x <- x[-ee$random]
p <- ee$parList(x=x)
if (!unlist) return(p)
p <- unlist(p[names(p)!="b"]) ## drop primary RE
names(p) <- gsub("[0-9]+$","",names(p)) ## remove disambiguators
return(p)
}
## replacement for (unexported) TMB:::isNullPointer
isNullPointer <- function(x) {
attributes(x) <- NULL
identical(x, new("externalptr"))
}
#' conditionally update glmmTMB object fitted with an old TMB version
#'
#' @rdname gt_load
#' @param oldfit a fitted glmmTMB object
#' @param update_gauss_disp update \code{betadisp} from variance to SD parameterization?
#' @export
up2date <- function(oldfit, update_gauss_disp = FALSE) {
openmp(1) ## non-parallel/make sure NOT grabbing all the threads!
obj <- oldfit$obj
ee <- obj$env
if (isNullPointer(oldfit$obj$env$ADFun$ptr)) {
pars <- c(grep("last\\.par", names(ee), value = TRUE), "par",
"parfull")
## using ee$parList() rather than ee$parameters should help
## with mapped parameter s... ??
params <- ee$parList()
if (length(ee$map) > 0) {
for (n in names(ee$map)) {
ee$parameters[[n]] <- params[[n]]
}
}
## change name of thetaf to psi
if ("thetaf" %in% names(params)) {
ee$parameters$psi <- params$thetaf
ee$parameters$thetaf <- NULL
pars <- c(grep("last\\.par", names(ee), value = TRUE),
"par")
for (p in pars) {
if (!is.null(nm <- names(ee[[p]]))) {
names(ee[[p]])[nm == "thetaf"] <- "psi"
}
}
}
if ("betad" %in% names(ee$parameters)) { #FIXME: DRY
ee$parameters$betadisp <- params$betad
ee$parameters$betad <- NULL
pars <- c(grep("last\\.par", names(ee), value = TRUE),
"par")
for (p in pars) {
if (!is.null(nm <- names(ee[[p]]))) {
names(ee[[p]])[nm == "betad"] <- "betadisp"
}
}
ee$data$Xdisp <- ee$data$Xd
ee$data$Xd <- NULL
ee$data$dispoffset <- ee$data$doffset
ee$data$doffset <- NULL
}
if(!"Zdisp" %in% names(ee$data)) {
ee$data$Zdisp <- new("dgTMatrix",Dim=c(as.integer(nrow(ee$data$Xdisp)),0L)) ## matrix(0, ncol=0, nrow=nobs)
ee$parameters$bdisp <- rep(0, ncol(ee$data$Zdisp))
ee$parameters$thetadisp <- numeric(0)
}
if (!"aggregate" %in% names(ee$data)) {
ee$data[["aggregate"]] <- numeric(0)
}
ee2 <- oldfit$sdr$env
if ("thetaf" %in% names(ee2$parameters)) {
ee2$parameters$psi <- ee2$parameters$thetaf
ee2$parameters$thetaf <- NULL
}
for (i in seq_along(ee$data$terms)) {
ee$data$terms[[i]]$simCode <- .valid_simcode[["random"]]
}
for (i in seq_along(ee$data$termszi)) {
ee$data$termszi[[i]]$simCode <- .valid_simcode[["random"]]
}
if ("betad" %in% names(ee2$parameters)) { #FIXME: DRY
ee2$parameters$betadisp <- ee2$parameters$betad
ee2$parameters$betad <- NULL
}
## prior_ivars, prior_fvars are defined in priors.R
if (!"prior_distrib" %in% names(ee$data)) {
## these are DATA_IVECTOR but apparently after processing
## TMB turns these into numeric ... ??
for (v in prior_ivars) ee$data[[v]] <- numeric(0)
for (v in prior_fvars) ee$data[[v]] <- numeric(0)
}
## switch from variance to SD parameterization
if (update_gauss_disp &&
family(oldfit)$family == "gaussian") {
ee$parameters$betadisp <- params$betadisp/2
for (p in pars) {
if (!is.null(nm <- names(ee[[p]]))) {
ee[[p]][nm == "betadisp"] <- ee[[p]][nm == "betadisp"]/2
}
if (!is.null(nm <- names(oldfit$fit[[p]]))) {
oldfit$fit[[p]][nm == "betadisp"] <- oldfit$fit[[p]][nm == "betadisp"]/2
}
}
}
oldfit$obj <- with(ee,
TMB::MakeADFun(data,
parameters,
map = map,
random = random,
silent = silent,
DLL = "glmmTMB"))
oldfit$obj$env$last.par.best <- ee$last.par.best
##
}
for (t in c("condReStruc", "ziRestruc", "dispRestruc")) {
for (i in seq_along(oldfit$modelInfo$reStruc[[t]])) {
oldfit$modelInfo$reStruc[[t]][[i]]$simCode <- .valid_simcode[["random"]]
}
}
## changed format of 'parallel' control to add autopar info
if (length(p <- oldfit$modelInfo$parallel) <= 1) {
if (!(is.null(p) || is.numeric(p))) {
stop("oldfit$modelInfo$parallel has an unexpected value")
}
oldfit$modelInfo$parallel <- list(n = p, autopar = NULL)
}
## dispersion was NULL rather than 1 in old R versions ...
omf <- oldfit$modelInfo$family
if (getRversion() >= "4.3.0" &&
!("dispersion" %in% names(omf))) {
## don't append() or c(), don't want to lose class info
oldfit$modelInfo$family$dispersion <- 1
}
if (!"priors" %in% names(oldfit$modelInfo)) {
## https://stackoverflow.com/questions/7944809/assigning-null-to-a-list-element-in-r
## n.b. can't use ...$priors <- NULL
oldfit$modelInfo["priors"] <- list(NULL)
}
if (!"aggregate" %in% names(ee$data)) {
ee$data[["aggregate"]] <- numeric(0)
}
if ("Xd" %in% names(ee$data)) {
ee$data[["Xdisp"]] <- ee$data[["Xd"]]
ee$data[["Xd"]] <- NULL
}
if ("XdS" %in% names(ee$data)) {
ee$data[["XdispS"]] <- ee$data[["XdS"]]
ee$data[["XdS"]] <- NULL
}
return(oldfit)
}
#' Load data from system file, updating glmmTMB objects
#'
#' @param fn partial path to system file (e.g. test_data/foo.rda)
#' @param verbose print names of updated objects?
#' @param mustWork fail if file not found?
#' @param \dots values passed through to \code{up2date}
#' @export
gt_load <- function(fn, verbose=FALSE, mustWork = FALSE, ...) {
sf <- system.file(fn, package = "glmmTMB")
found_file <- file.exists(sf)
if (mustWork && !found_file) {
stop("couldn't find system file ", sf)
}
L <- load(sf)
for (m in L) {
if (inherits(get(m), "glmmTMB")) {
if (verbose) cat(m,"\n")
assign(m, up2date(get(m), ...))
}
assign(m, get(m), parent.env(), envir = parent.frame())
}
return(found_file)
}
#' truncated distributions
#'
#' Probability functions for k-truncated Poisson and negative binomial distributions.
#' @param x value
#' @param size number of trials/overdispersion parameter
#' @param mu mean parameter
#' @param k truncation parameter
#' @param log (logical) return log-probability?
#' @export
dtruncated_nbinom2 <- function(x, size, mu, k=0, log=FALSE) {
y <- ifelse(x<=k,-Inf,
dnbinom(x, mu=mu, size=size, log=TRUE) -
pnbinom(k, mu=mu, size=size, lower.tail=FALSE,
log.p=TRUE))
if (log) return(y) else return(exp(y))
}
#' @rdname dtruncated_nbinom2
#' @param lambda mean parameter
#' @importFrom stats dpois
#' @export
dtruncated_poisson <- function(x,lambda,k=0,log=FALSE) {
y <- ifelse(x<=k,-Inf,
dpois(x,lambda,log=TRUE) -
ppois(k, lambda=lambda, lower.tail=FALSE,
log.p=TRUE))
if (log) return(y) else return(exp(y))
}
#' @rdname dtruncated_nbinom2
#' @param phi overdispersion parameter
#' @export
dtruncated_nbinom1 <- function(x, phi, mu, k=0, log=FALSE) {
## V=mu*(1+phi) = mu*(1+mu/k) -> k=mu/phi
size <- mu/phi
y <- ifelse(x<=k,-Inf,
dnbinom(x,mu=mu, size=size,log=TRUE) -
pnbinom(k, mu=mu, size=size, lower.tail=FALSE,
log.p=TRUE))
if (log) return(y) else return(exp(y))
}
## utilities for constructing lists of parameter names
## for matching map names vs nameList components ...
par_components <- c("beta","betazi","betadisp","theta","thetazi","psi")
## all parameters, including both mapped and rank-dropped
getParnames <- function(object, full, include_dropped = TRUE, include_mapped = TRUE,
include_mapequal = FALSE) {
mkNames <- function(tag="") {
X <- getME(object,paste0("X",tag))
if (is.null(X) && tag == "disp") stop("are you using a stored model from an earlier glmmTMB version? try running up2date()")
dropped <- attr(X, "col.dropped") %||% numeric(0)
ntot <- ncol(X) + length(dropped)
## identical instead of ==; ncol(X) may be NULL for older models
if (ntot == ncol(X) || !include_dropped) {
nn <- colnames(X)
} else {
nn <- character(ntot)
nn[-dropped] <- colnames(X)
nn[ dropped] <- names(dropped)
}
if (trivialFixef(nn, tag)
&& !(tag == "disp" && (!trivialDisp(object) && !zeroDisp(object)))
## if 'full', keep disp even if trivial, if used by family
&& !(full && tag =="disp" &&
(usesDispersion(family(object)$family) && !zeroDisp(object)))) {
return(character(0))
}
if (tag == "") return(nn)
return(paste(tag,nn,sep="~"))
}
nameList <- setNames(Map(mkNames, c("", "zi", "disp")),
names(cNames))
if(full) {
## FIXME: haven't really decided if we should drop the
## trivial variance-covariance dispersion parameter ??
## if (trivialDisp(object))
## res <- covF[-nrow(covF),-nrow(covF)]
reNames <- function(tag) {
re <- object$modelInfo$reStruc[[paste0(tag,"ReStruc")]]
num_theta <- vapply(re,"[[","blockNumTheta", FUN.VALUE = numeric(1))
nn <- mapply(function(n,L) paste(n, seq(L), sep="."),
names(re), num_theta)
if (length(nn) == 0) return(nn)
return(paste("theta",gsub(" ", "", unlist(nn)), sep="_"))
}
## nameList for estimated variables;
nameList <- c(nameList,
list(theta = reNames("cond"),
thetazi = reNames("zi"),
thetadisp = reNames("disp")))
##
if (length(fp <- family_params(object)) > 0) {
nameList <- c(nameList, list(psi = names(fp)))
}
}
if (!include_mapped) {
map <- object$obj$env$map
if (length(map)>0) {
for (m in seq_along(map)) {
if (!include_mapequal) {
missmap <- which(is.na(map[[m]]) | duplicated(map[[m]]))
} else {
missmap <- which(is.na(map[[m]]))
}
if (length(missmap)>0) {
w <- match(names(map)[m], par_components) ##
if (length(nameList)>=w) { ## may not exist if !full
nameList[[w]] <- nameList[[w]][-missmap]
}
} ## if NA/duplicated values
} ## for (m in seq_along(map))
} ## if (length(map) > 0)
}
return(nameList)
}
## OBSOLETE (?)
## reassign predvars to have term vars in the right order,
## but with 'predvars' values inserted where appropriate
fix_predvars <- function(pv,tt) {
if (length(tt)==3) {
## convert two-sided to one-sided formula
tt <- RHSForm(tt, as.form=TRUE)
}
## ugh, deparsing again ...
tt_vars <- vapply(attr(tt, "variables"), deparse1, character(1))[-1]
## remove terminal paren - e.g. match term poly(x, 2) to
## predvar poly(x, 2, <stuff>)
## beginning of string, including open-paren, colon
## but not *first* comma nor arg ...
## could possibly try init_regexp <- "^([^,]+).*" ?
init_regexp <- "^([(^:_.[:alnum:]]+).*"
tt_vars_short <- gsub(init_regexp,"\\1",tt_vars)
if (is.null(pv) || length(tt_vars)==0) return(NULL)
new_pv <- quote(list())
## maybe multiple variables per pv term ... [-1] ignores head
## FIXME: test for really long predvar strings ????
pv_strings <- vapply(pv,deparse1,FUN.VALUE=character(1))[-1]
pv_strings <- gsub(init_regexp,"\\1",pv_strings)
for (i in seq_along(tt_vars)) {
w <- match(tt_vars_short[[i]],pv_strings)
if (!is.na(w)) {
new_pv[[i+1]] <- pv[[w+1]]
} else {
## insert symbol from term vars
new_pv[[i+1]] <- as.symbol(tt_vars[[i]])
}
}
return(new_pv)
}
collapse_list <- function(pList) {
## workaround to get non-unique names ...
pList <- mapply(function(x, n) { setNames(x, rep(n, length(x))) },
pList, names(pList))
pvec <- unlist(unname(pList))
}
make_pars <- function(pList, ..., include_extra = TRUE) {
## FIXME: check for name matches, length matches etc.
## (useful errors)
## better to split by name first??
L <- list(...)
unmatched <- setdiff(names(L), names(pList))
if (length(unmatched) > 0) {
warning(sprintf("unmatched parameter names: %s",
paste(unmatched, collapse =", ")))
}
if (!include_extra) L <- L[intersect(names(L), names(pList))]
for (nm in names(L)) {
if ((len1 <- length(pList[[nm]])) == (len2 <- length(L[[nm]]))) {
pList[[nm]] <- L[[nm]]
} else {
## skip cases with different length (== partially-mapped vectors)
plen <- len1 + length(attr(pList[[nm]], "map"))
## if L[[nm]] is a list, that's because we're
## passing a partial b vector ...
if (plen != len2 && !is.list(L[[nm]])) {
warning(sprintf("length mismatch in component %s (%d != %d); not setting",
nm, len1, len2))
}
}
}
return(collapse_list(pList))
}
##' helper function to modify simulation settings for random effects
##'
##' This modifies the TMB object \emph{in place} (beware!)
##' Ultimately this will allow \code{terms} to be a vector of term names,
##' with a matching \code{val} vector to specify the behaviour for each term
##'
##' @param g a TMB object
##' @param val a legal setting for sim codes ("zero", "random", or "fix")
##' @param terms which terms to apply this to
##' @export
##'
set_simcodes <- function(g, val = "zero", terms = "ALL") {
ee <- g$env
if (terms != "ALL") stop("termwise setting of simcodes not implemented yet")
if (terms == "ALL") {
for (i in seq_along(ee$data$terms)) {
ee$data$terms[[i]]$simCode <- .valid_simcode[[val]]
}
}
}
##' Simulate from covariate/metadata in the absence of a real data set (EXPERIMENTAL)
##'
##' See \code{vignette("sim", package = "glmmTMB")} for more details and examples,
##' and \code{vignette("covstruct", package = "glmmTMB")}
##' for more information on the parameterization of different covariance structures.
##'
##' @inheritParams glmmTMB
##' @param object a \emph{one-sided} model formula (e.g. \code{~ a + b + c}
##' (peculiar naming is for consistency with the generic function, which typically
##' takes a fitted model object)
##' @param nsim number of simulations
##' @param seed random-number seed
##' @param newdata a data frame containing all variables listed in the formula,
##' \emph{including} the response variable (which needs to fall within
##' the domain of the conditional distribution, and should probably not
##' be all zeros, but whose value is otherwise irrelevant)
##' @param newparams a list of parameters containing sub-vectors
##' (\code{beta}, \code{betazi}, \code{betadisp}, \code{theta}, etc.) to
##' be used in the model. If \code{b} is specified in this list, then the conditional modes/BLUPs
##' will be set to these values; otherwise they will be drawn from the appropriate Normal distribution.
##' See \code{vignette("covstruct", package = "glmmTMB")} for details on the parameterizations used
##' for various random-effects models (i.e., \code{theta}).
##' @param ... other arguments to \code{glmmTMB} (e.g. \code{family})
##' @param return_val what information to return: "sim" (the default) returns a list of vectors of simulated outcomes; "pars" returns the default parameter vector (this variant does not require \code{newparams} to be specified, and is useful for figuring out the appropriate dimensions of the different parameter vectors); "object" returns a fake \code{glmmTMB} object (useful, e.g., for retrieving the Z matrix (\code{getME(simulate_new(...), "Z")}) or covariance matrices (\code{VarCorr(simulate_new(...))}) implied by a particular set of input data and parameter values)
##' @details Use the \code{weights} argument to set the size/number of trials per observation for binomial-type models; the default is 1 for every observation (i.e., Bernoulli trials)
##' @seealso \code{\link{glmmTMB}}, \code{\link{family_glmmTMB}} (for conditional distribution parameterizations [\code{betadisp}]), \code{\link{put_cor}} (for correlation matrix parameterizations)
##' @examples
##' ## use Salamanders data for observational design and covariate values
##' ## parameters used here are sensible, but do not fit the original data
##' params <- list(beta = c(2, 1),
##' betazi = c(-0.5, 0.5), ## logit-linear model for zi
##' betadisp = log(2), ## log(NB dispersion)
##' theta = log(1)) ## log(among-site SD)
##' sim_count <- simulate_new(~ mined + (1|site),
##' newdata = Salamanders,
##' zi = ~ mined,
##' family = nbinom2,
##' seed = 101,
##' newparams = params
##' )
##' ## simulate_new with return="sim" always returns a list of response vectors
##' Salamanders$sim_count <- sim_count[[1]]
##' summary(glmmTMB(sim_count ~ mined + (1|site), data=Salamanders, ziformula=~mined, family=nbinom2))
##' ## return a glmmTMB object
##' sim_obj <- simulate_new(~ mined + (1|site),
##' return_val = "object",
##' newdata = Salamanders,
##' zi = ~ mined,
##' family = nbinom2,
##' newparams = params)
##' ## simulate Gaussian data, multivariate random effect
##' data("sleepstudy", package = "lme4")
##' sim_obj <- simulate_new(~ 1 + (1|Subject) + ar1(0 + factor(Days)|Subject),
##' return_val = "pars",
##' newdata = sleepstudy,
##' family = gaussian,
##' newparams = list(beta = c(280, 1),
##' betad = log(2), ## log(residual std err)
##' theta = c(log(2), ## log(SD(subject))
##' log(2), ## log(SD(slope))
##' ## AR1 correlation = 0.2
##' put_cor(0.2, input_val = "vec"))
##' )
##' )
##'
##' @export
simulate_new <- function(object,
nsim = 1,
seed = NULL,
family = gaussian,
newdata, newparams, ...,
return_val = c("sim", "pars", "object")) {
return_val <- match.arg(return_val)
family <- get_family(family, deparse(substitute(family)))
## truncate
if (length(object) == 3) stop("simulate_new should take a one-sided formula")
newparams0 <- newparams
## store original params
## (in case we need both complete-b and unmapped-b versions)
## fill in fake LHS
form <- object
form[[3]] <- form[[2]]
form[[2]] <- quote(..y)
## insert a legal value: 1.0 is OK as long as family != "beta"
## (note the family *function* is 'beta_family' but the internal
## $family value is 'beta')
newdata[["..y"]] <- if (family$family == "beta") 0.5 else 1.0
r1 <- glmmTMB(form,
data = newdata,
family = family,
## make sure optim doesn't actually do anything
## (if return_val is "object")
control = glmmTMBControl(optCtrl = list(iter.max = 0)),
...,
doFit = FALSE)
## sort out components of b (if necessary)
if ("b" %in% names(newparams)) {
components <- c("cond", "zi")
cnames <- paste0(components, "ReStruc")
tnames <- c("terms", "termszi")
if (!is.list(newparams$b)) {
b_inds <- seq_along(newparams$b)
for (i in seq_along(cnames)) { ## loop over RE categories in object ('cond' and 'zi')
for (j in seq_along(r1[[cnames[i]]])) { ## loop over terms within the category
r1[[cnames[i]]][[j]]$simCode <- .valid_simcode[["fix"]]
r1$data.tmb[[tnames[i]]][[j]]$simCode <- .valid_simcode[["fix"]]
}
}
} else {
restrucs <- r1[cnames]
b_inds <- get_b_inds(restrucs, names(newparams$b))
b_terms <- get_b_inds(restrucs, names(newparams$b),
ret_val = "terms")
## b_terms gives indices in overall sequence, not
## indices per term -- we have to split these by
## term. Should probably be handled upstream?
for (i in seq_along(cnames)) {
nre <- length(r1[[i]])
cur_b <- which(b_terms < nre)
if (length(cur_b) > 0) {
for (j in b_terms[cur_b]) {
## ugh, have to set in two places: do this upstream
## of glmmTMB() call?
r1[[cnames[i]]][[j]]$simCode <- .valid_simcode[["fix"]]
r1$data.tmb[[tnames[i]]][[j]]$simCode <- .valid_simcode[["fix"]]
}
b_terms <- b_terms[-cur_b]
}
b_terms <- b_terms - nre
}
}
n_b <- length(r1$parameters$b)
b_fac <- seq(n_b)
b_fac[unlist(b_inds)] <- NA
## TMB complains if number of levels doesn't match values:
## make into a factor *after* setting NA values
b_fac <- factor(b_fac)
new_b <- rep(0, n_b)
new_b[unlist(b_inds)] <- unlist(newparams$b)
r1$parameters$b <- new_b
r1$map <- r1$mapArg <- list(b = b_fac)
}
## construct TMB object, but don't fit it
## (for cnms etc., simulations, etc.)
r2 <- fitTMB(r1, doOptim = FALSE)
set_b <- function(x, b, map = NULL) {
if (!is.null(map)) {
b <- b[!is.na(map)]
}
pList <- split(x, names(x))
pList$b <- b
return(collapse_list(pList))
}
pars <- do.call("make_pars",
c(list(r2$env$parameters), newparams,
list(include_extra = FALSE)))
if (!is.null(seed)) set.seed(seed)
if (return_val %in% c("pars", "object")) {
b_vals <- r2$simulate(par = pars)$b
if (return_val == "pars") {
pars <- do.call("make_pars",
c(list(r2$env$parameters), newparams,
list(include_extra = TRUE)))
return(set_b(pars, b_vals))
}
## FIXME: need to set start= values as well when simulating
## object??
## what do we do when b is specified/partially specified?
for (nm in names(newparams)) {
r1$parameters[[nm]] <- newparams[[nm]]
}
r3 <- suppressWarnings(fitTMB(r1, doOptim = TRUE))
r3$fit$parfull <- set_b(r3$fit$parfull, b_vals, r1$map$b)
return(r3)
}
replicate(nsim, r2$simulate(par = pars)$yobs, simplify = FALSE)
}
set_class <- function(x, cls, prepend = TRUE) {
if (is.null(x)) return(NULL)
if (!prepend) class(x) <- cls
else class(x) <- c(cls, class(x))
x
}
## convert from parameter name to component name or vice versa
## first name shoudl be em
compsyn <- c(cond = "", zi = "zi", disp = "d")
match_names <- function(x, to_parvec = FALSE, prefix = "beta") {
if (to_parvec) {
## "cond" -> "theta" etc.
return(paste0(prefix, compsyn[x]))
} else {
## "beta" -> "cond" etc.
x <- gsub(prefix, "", x)
return(names(compsyn)[match(x, compsyn)])
}
}
get_family <- function(family, f_name) {
if (is.character(family)) {
if (family=="beta") {
family <- "beta_family"
warning("please use ",sQuote("beta_family()")," rather than ",
sQuote("\"beta\"")," to specify a Beta-distributed response")
}
family <- get(family, mode = "function", envir = parent.frame(2))
}
if (is.function(family)) {
## call family with no arguments
family <- family()
}
## FIXME: what is this doing? call to a function that's not really
## a family creation function?
if (isS4(family) || is.null(family$family)) {
stop(sprintf("after evaluation, 'family' (specified as '%s') must be a list containing a '$family' element", f_name))
}
return(family)
}
## ugh, better way to do this?
get_re_names <- function(re) {
names(unlist(sapply(re, function(x) sapply(x, function(y) NA))))
}
#' @param reStrucs a list containing conditional and z-i RE structures
#' @param b_names vector of names matching RE terms
#' @examples
#' data("sleepstudy", package = "lme4")
#' fm1 <- glmmTMB(Reaction ~ 1 + (1|Subject) + ar1(0+factor(Days)|Subject), sleepstudy)
#' re <- fm1$modelInfo$reStruc
#' get_b_inds(re, "1|Subject")
#' @noRd
get_b_terms <- function(nms, inms) {
squash_ws <- function(x) gsub(" ", "", x)
nms <- squash_ws(nms)
inms <- squash_ws(inms)
w <- match(nms, inms)
unmatched <- which(is.na(w))
if (length(unmatched)>0) {
w[unmatched] <- match(nms[unmatched],
gsub("^(cond|zi)ReStruc\\.", "", inms))
}
if (any(is.na(w))) {
stop("unmatched RE terms")
}
return(w)
}
get_b_inds <- function(reStrucs, b_names, ret_val = c("indices", "terms")) {
ret_val <- match.arg(ret_val)
bfun <- function(x) {
if (length(x) == 0) return(numeric(0))
with(x, blockSize * blockReps)
}
retrms <- unlist(sapply(reStrucs, function(x) sapply(x, bfun)))
## set up indices ...
if (ret_val == "terms") {
return(get_b_terms(b_names, get_re_names(reStrucs)))
}
w <- get_b_terms(b_names, get_re_names(reStrucs))
inds <- cumsum(c("start" = 0, retrms))
## first try to match full name (component + term)
## set specified values
res <- lapply(w, function(i) seq(inds[i]+1, inds[i+1]))
names(res) <- b_names
res
}
## add negative-value check to binomial initialization method
our_binom_initialize <- function(family) {
newtest <- substitute(
## added test for glmmTMB
if (any(y<0)) {
stop(sprintf('negative values not allowed for the %s family', FAMILY))
}
, list(FAMILY=family))
b0 <- binomial()$initialize
b0[[length(b0)+1]] <- newtest
return(b0)
}
#' retrieve current value of TMB autopar setting
get_autopar <- function() {
attr(TMB::openmp(DLL = "glmmTMB"), "autopar")
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.