diff --git a/DESCRIPTION b/DESCRIPTION index b375270..63b9a41 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -16,8 +16,8 @@ Suggests: knitr, markdown, testthat -Description: Macros to generate 'nimble' code from a concise syntax. Included are macros for generating linear modeling code using a formula-based syntax, and for 'for' loops. (de Valpine et al. (2015) ). -License: GPL-3 +Description: Macros to generate 'nimble' code from a concise syntax. Included are macros for generating linear modeling code using a formula-based syntax and for 'for' loops. +License: BSD_3_clause + file LICENSE | GPL (>= 2) NeedsCompilation: no Encoding: UTF-8 VignetteBuilder: knitr diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..de82bb8 --- /dev/null +++ b/LICENSE @@ -0,0 +1,3 @@ +YEAR: 2025 +COPYRIGHT HOLDER: Ken Kellner, Perry de Valpine, Christopher Paciorek, Daniel Turek +ORGANIZATION: University of California diff --git a/NAMESPACE b/NAMESPACE index 296be55..715aafa 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,10 +4,6 @@ export(FORLOOP) export(LINPRED) export(LINPRED_PRIORS) export(LM) -export(formulaHandler_I) -export(formulaHandler_log) -export(formulaHandler_offset) -export(formulaHandler_scale) export(matchPrior) export(setPriors) export(uppertri_mult_diag) diff --git a/R/FORLOOP.R b/R/FORLOOP.R index 28d6a8f..d4093ae 100644 --- a/R/FORLOOP.R +++ b/R/FORLOOP.R @@ -10,6 +10,8 @@ #' @author Ken Kellner and Perry de Valpine #' #' @param code The right-hand side of a parameter declaration +#' @param modelInfo Used internally by nimbleMacros; a list of model information such as constants and dimensions +#' @param .env Used internally by nimbleMacros; the environment where the model was created #' #' @examples #' code <- nimbleCode({ @@ -17,7 +19,7 @@ #' mu[1:n] <- FORLOOP(beta[1] + beta[2]*x[1:n]) #' }) #' -#' mod <- nimbleModel(code, constants=list(n=10)) +#' mod <- nimbleModel(code, constants = list(n=10)) #' mod$getCode() NULL diff --git a/R/LINPRED.R b/R/LINPRED.R index f260008..d16f6e9 100644 --- a/R/LINPRED.R +++ b/R/LINPRED.R @@ -10,20 +10,35 @@ #' @param formula An R formula, possibly with the parameters followed by #' brackets containing indices. If there are no indices, the macro attempts #' to guess the correct indices from the context. The formula must be -#' right-hand side only (e.g. ~x). This must always be the first argument supplied -#' to LINPRED. -#' @param link A link function which will be applied to the +#' right-hand side only (e.g. \code{~x}). This must always be the first argument supplied +#' to \code{LINPRED}. +#' @param link A link function that will be applied to the #' left-hand-side (the response) in the final linear predictor. Default is none. #' @param coefPrefix All model coefficient names will begin with this prefix. -#' default is beta_ (so x becomes beta_x, etc.) +#' default is \code{"beta_"} (so 'x' becomes 'beta_x', etc.) #' @param sdPrefix All dispersion parameters will begin with this prefix. -#' default is no prefix. -#' @param priorSpecs Prior specifications, should be generated with setPrior() -#' @param modMatNames Logical, should parameters be named so they match the -#' names you would get from R's model.matrix function? -#' @param noncenter Logical; use noncentered parameterization? -#' @param centerVar Grouping covariate to 'center' on in parameterization. By -#' default all random effects have mean 0 as with lme4. +#' Default is no prefix. +#' @param priors Prior specifications, generated with \code{setPrior()} +#' @param modelMatrixNames Logical indicating if parameters should be named so they match the +#' names one would get from R's \code{model.matrix}. Default is \code{FALSE}. +#' @param noncentered Logical indicating whether to use noncentered parameterization +#' for random effects. Default is \code{FALSE}. Under the noncentered parameterization, +#' random effects have a standard normal prior (\code{beta_x_raw ~ dnorm(0, sd = 1)}) +#' and are then scaled by the hyperparameters (mean and SD), i.e., +#' \code{beta_x = mu_beta + sd_beta * beta_x_raw}. This parameterization can improve +#' mixing in some cases. +#' @param centerVar Grouping variable (covariate) to 'center' the random effects on. +#' By default (when NULL), random effects come from normal distributions with mean 0 +#' as with \code{lme4}. For example, for random intercepts by grouping variable \code{x}, +#' the linear predictor would be \code{beta_Intercept + beta_x[x[i]]} and the +#' prior for the random effects would be \code{beta_x ~ dnorm(0, sd_x)}. When +#' \code{centerVar = x}, the linear predictor would be \code{beta_x[x[i]]} +#' and the random effect prior would be \code{beta_x ~ dnorm(beta_Intercept, sd = sd_x)}. +#' That is, the mean of the random effects is now \code{beta_Intercept}. +#' These two formulations should yield the same results. Note that this option +#' is unrelated to the \code{noncentered} argument despite the similar names. +#' @param modelInfo Used internally by nimbleMacros; a list of model information such as constants and dimensions +#' @param .env Used internally by nimbleMacros; the environment where the model was created #' #' @author Ken Kellner #' @@ -33,18 +48,18 @@ #' mu[1:3] <- LINPRED(~x + x2) #' }) #' -#' mod <- nimbleModel(code, constants=constants) +#' mod <- nimbleModel(code, constants = constants) #' mod$getCode() NULL #' @export LINPRED <- nimble::buildMacro( function(stoch, LHS, formula, link=NULL, coefPrefix=quote(beta_), - sdPrefix=NULL, priorSpecs=setPriors(), modMatNames = FALSE, - noncenter = FALSE, centerVar=NULL, modelInfo, .env){ + sdPrefix=NULL, priors=setPriors(), modelMatrixNames = FALSE, + noncentered = FALSE, centerVar=NULL, modelInfo, .env){ # Make sure formula is in correct format - formula <- stats::as.formula(formula) + formula <- check_formula(formula) # Get index range on LHS to use if the RHS formulas do not specify them LHS_ind <- extractBracket(LHS) @@ -74,13 +89,13 @@ function(stoch, LHS, formula, link=NULL, coefPrefix=quote(beta_), code <- substitute(LHS <- nimbleMacros::FORLOOP(RHS), list(LHS = LHS, RHS = RHS)) # Add code for priors to output if needed - if(!is.null(priorSpecs)){ + if(!is.null(priors)){ priorCode <- substitute(nimbleMacros::LINPRED_PRIORS(FORMULA, coefPrefix=COEFPREFIX, sdPrefix=SDPREFIX, - priorSpecs=PRIORSET, modMatNames=MODMAT, - noncenter = UNCENTER, centerVar=CENTERVAR), + priors=PRIORSET, modelMatrixNames=MODMAT, + noncentered = UNCENTER, centerVar=CENTERVAR), list(COEFPREFIX=coefPrefix, FORMULA=formula, SDPREFIX=sdPrefix, - PRIORSET=priorSpecs, MODMAT=modMatNames, - UNCENTER = noncenter, CENTERVAR=centerVar)) + PRIORSET=priors, MODMAT=modelMatrixNames, + UNCENTER = noncentered, CENTERVAR=centerVar)) code <- embedLinesInCurlyBrackets(list(code, priorCode)) } @@ -92,29 +107,44 @@ unpackArgs=TRUE ) -#' Macro to build code for priors on a linear predictor from R formula +#' Macro to build code for priors on a linear predictor from an R formula #' #' Generates appropriate priors for a linear predictor derived from an #' R formula. As such it makes the most sense to use this macro together with -#' the LINPRED macro which takes similar arguments. +#' the LINPRED macro, which takes similar arguments. #' #' @name LINPRED_PRIORS #' @author Ken Kellner #' -#' @param formula An R formula The formula must be right-hand side only (e.g. ~x). -#' This must always be the first argument supplied to LINPRED_PRIORS +#' @param formula An R formula The formula must be right-hand side only (e.g., \code{~x}). +#' This must always be the first argument supplied to \code{LINPRED_PRIORS}. #' @param coefPrefix All model coefficient names will begin with this prefix. -#' default is beta_ (so x becomes beta_x, etc.) +#' default is \code{"beta_"} (so 'x' becomes 'beta_x', etc.) #' @param sdPrefix All dispersion parameters will begin with this prefix. -#' default is no prefix. -#' @param priorSpecs List of prior specifications, should be generated using +#' Default is no prefix. +#' @param priors List of prior specifications, generated using \code{setPriors}. #' setPriors() -#' @param modMatNames Logical, should parameters be named so they match the -#' names you would get from R's model.matrix function? -#' @param noncenter Logical, use noncentered parameterization? -#' @param centerVar Grouping covariate to 'center' on in parameterization. By -#' default all random effects have mean 0 as with lme4. -#' +#' @param modelMatrixNames Logical indicating if parameters should be named so they match the +#' names one would get from R's \code{model.matrix}. Default is \code{FALSE}. +#' @param noncentered Logical indicating whether to use noncentered parameterization +#' for random effects. Default is \code{FALSE}. Under the noncentered parameterization, +#' random effects have a standard normal prior (\code{beta_x_raw ~ dnorm(0, sd = 1)}) +#' and are then scaled by the hyperparameters (mean and SD), i.e., +#' \code{beta_x = mu_beta + sd_beta * beta_x_raw}. This parameterization can improve +#' mixing in some cases. +#' @param centerVar Grouping variable (covariate) to 'center' the random effects on. +#' By default (when NULL), random effects come from normal distributions with mean 0 +#' as with \code{lme4}. For example, for random intercepts by grouping variable \code{x}, +#' the linear predictor would be \code{beta_Intercept + beta_x[x[i]]} and the +#' prior for the random effects would be \code{beta_x ~ dnorm(0, sd_x)}. When +#' \code{centerVar = x}, the linear predictor would be \code{beta_x[x[i]]} +#' and the random effect prior would be \code{beta_x ~ dnorm(beta_Intercept, sd = sd_x)}. +#' That is, the mean of the random effects is now \code{beta_Intercept}. +#' These two formulations should yield the same results. Note that this option +#' is unrelated to the \code{noncentered} argument despite the similar names. +#' @param modelInfo Used internally by nimbleMacros; a list of model information such as constants and dimensions +#' @param .env Used internally by nimbleMacros; the environment where the model was created +#' #' @author Ken Kellner #' #' @examples @@ -129,15 +159,14 @@ NULL #' @export LINPRED_PRIORS <- nimble::buildMacro( -function(formula, coefPrefix=quote(beta_), sdPrefix=NULL, priorSpecs=setPriors(), - modMatNames=FALSE, noncenter = FALSE, centerVar=NULL, modelInfo, .env){ +function(formula, coefPrefix=quote(beta_), sdPrefix=NULL, priors=setPriors(), + modelMatrixNames=FALSE, noncentered = FALSE, centerVar=NULL, modelInfo, .env){ # Make sure formula is in correct format - if(formula[[1]] != quote(`~`)) formula <- c(quote(`~`),formula) - formula <- stats::as.formula(formula) + formula <- check_formula(formula) # Evaluate prior settings - priorSpecs <- eval(priorSpecs, envir=.env) + priors <- eval(priors, envir=.env) # Split formula into components and process the components components <- buildLP(formula, defaultBracket = "[]", # default bracket not used below @@ -149,8 +178,8 @@ function(formula, coefPrefix=quote(beta_), sdPrefix=NULL, priorSpecs=setPriors() components <- buildPriors(components, coefPrefix=safeDeparse(coefPrefix), sdPrefix=sdPrefix, modelInfo = modelInfo, - priorSpecs=priorSpecs, - modMatNames = modMatNames, noncenter=noncenter) + priorSpecs=priors, + modMatNames = modelMatrixNames, noncenter=noncentered) # Get complete prior code code <- getPriors(components) @@ -168,6 +197,16 @@ use3pieces=FALSE, unpackArgs=TRUE ) +# Check validity of formula +check_formula <- function(formula){ + if(formula[[1]] != quote(`~`)) formula <- c(quote(`~`),formula) + formula <- stats::as.formula(formula) + if(length(formula) > 2){ + stop("check_formula: Formula should be RHS-only, such as ~1 or ~x, but not y~x", + call.=FALSE) + } + formula +} # This function starts with a formula plus options # Splits the formula into separate components (terms) @@ -363,9 +402,17 @@ processFormulaHandler <- function(x, defaultBracket, coefPrefix="beta_", cand <- paste0("formulaHandler_", func) + # Look for formula handler + # First look in user's environment, then look in nimbleMacros processor_available <- FALSE - if(exists(cand, envir = env)){ - processor <- get(cand, envir = env) + in_nimbleMacros <- exists(cand, envir = environment()) + in_user_env <- exists(cand, envir = env) + if(in_nimbleMacros | in_user_env){ + if(in_user_env){ + processor <- get(cand, envir = env) + } else { + processor <- get(cand, envir = environment()) # fallback to nimbleMacros + } if(inherits(processor, "nimbleFormulaHandler")){ processor_available <- TRUE out <- processor(x, defaultBracket, coefPrefix, sdPrefix, modelInfo, env, ...) @@ -1046,7 +1093,7 @@ addPriorsCode.formulaComponentIntercept <- function(x, priorSpecs, ...){ # Create prior code par <- str2lang(x$parameter) prior <- nimbleMacros::matchPrior(par, "intercept", - priorSpecs = priorSpecs) + priors = priorSpecs) code <- substitute(PAR ~ PRIOR, list(PAR=par, PRIOR = prior)) x$priorCode <- code x @@ -1104,7 +1151,7 @@ addPriorsCode.formulaComponentFixed <- function(x, priorSpecs, coefPrefix, modMa # specifically we need to add a line equating the default and model matrix names (NODE<-MODPAR) if(modMatNames & !node_modpar_match){ prior <- nimbleMacros::matchPrior(modmat_par, type, "coefficient", - priorSpecs = priorSpecs) + priors = priorSpecs) code_part <- substitute({ NODE <- MODPAR MODPAR ~ PRIOR @@ -1113,7 +1160,7 @@ addPriorsCode.formulaComponentFixed <- function(x, priorSpecs, coefPrefix, modMa # Otherwise just use the original name and we can skip one line of BUGS # code equating to the model matrix and default names prior <- nimbleMacros::matchPrior(node, type, "coefficient", - priorSpecs = priorSpecs) + priors = priorSpecs) code_part <- substitute(NODE ~ PRIOR, list(NODE = node, PRIOR=prior)) } code_part @@ -1154,7 +1201,7 @@ addPriorsCode.formulaComponentRandom <- function(x, priorSpecs, sdPrefix = NULL, } code1 <- lapply(sd_names, function(z){ - prior <- nimbleMacros::matchPrior(str2lang(z), "sd", priorSpecs = priorSpecs) + prior <- nimbleMacros::matchPrior(str2lang(z), "sd", priors = priorSpecs) substitute(SD ~ PRIOR, list(SD = str2lang(z), PRIOR = prior)) }) @@ -1308,12 +1355,12 @@ correlatedRandomPrior <- function(x, priorSpecs, sdPrefix, sd_name, modelInfo, c Ustar_name <- as.name(paste0("Ustar_", rfact)) U_name <- as.name(paste0("U_", rfact)) # LKJ distribution shape parameter - eta <- priorSpecs$eta - if(is.null(eta)) eta <- 1.3 + lkj_shape <- priorSpecs$lkjShape + if(is.null(lkj_shape)) lkj_shape <- 1 u <- substitute({ USTAR[1:NP, 1:NP] ~ dlkj_corr_cholesky(ETA, NP) U[1:NP, 1:NP] <- uppertri_mult_diag(USTAR[1:NP, 1:NP], SDS[1:NP]) - }, list(USTAR=Ustar_name, U=U_name, NP=as.numeric(np), SDS=str2lang(sd_vec), ETA=eta) + }, list(USTAR=Ustar_name, U=U_name, NP=as.numeric(np), SDS=str2lang(sd_vec), ETA=lkj_shape) ) # Get means for each random effect, default to 0 @@ -1378,12 +1425,12 @@ correlatedRandomPrior <- function(x, priorSpecs, sdPrefix, sd_name, modelInfo, c #' uppertri_mult_diag #' #' nimbleFunction needed when fitting correlated random effects. -#' Generates upper triangular Cholesky factor of covariance matrix (U in code) -#' from upper tri Cholesky factor of correlation matrix (Ustar in code) +#' Generates upper triangular Cholesky factor of covariance matrix ("U" in code) +#' from upper triangular Cholesky factor of correlation matrix ("Ustar" in code) #' and vector of standard deviations. Taken from the NIMBLE manual, -#' section 5.2.4.1.2 LKJ distribution for correlation matrices. +#' section 5.2.4.1.2 (LKJ Distribution for Correlation Matrices). #' -#' @param mat upper triangular Cholesky factor of correlation matrix (Ustar) +#' @param mat upper triangular Cholesky factor of correlation matrix ("Ustar") #' @param vec vector of standard deviations for individual random effects #' #' @name uppertri_mult_diag diff --git a/R/LM.R b/R/LM.R index b333331..1eee246 100644 --- a/R/LM.R +++ b/R/LM.R @@ -1,7 +1,8 @@ #' Macro for fitting linear models, GLMs, and GLMMs #' #' This macro generates code for LMs, GLMs, and GLMMs using formula notation -#' and arguments similar to R functions such as lm(), glm(), and lmer()/glmer(). +#' and arguments similar to R functions such as \code{lm()}, \code{glm()}, +#' and \code{lmer()}/\code{glmer()}. #' Currently only normal, Poisson, and binomial models are supported. #' #' @name LM @@ -10,19 +11,21 @@ #' @param formula An R formula, possibly with the parameters followed by #' brackets containing indices. If there are no indices, the macro attempts #' to guess the correct indices from the context. Formulas can include -#' random effects via lme4-style notation (e.g. ~ x + (1|group)) +#' random effects via lme4-style notation (e.g., \code{~ x + (1|group)}) #' @param family A description of the error distribution and link function to #' be used in the model. This can be a character string naming a family function, -#' a family function or the result of a call to a family function. See ?family. -#' Supported families are gaussian (default), binomial, and poisson. -#' @param coefPrefix All model coefficient names will begin with this prefix. -#' default is beta_ (so x becomes beta_x, etc.) -#' @param sdPrefix All dispersion parameters will begin with this prefix. +#' a family function or the result of a call to a family function. See \code{?family}. +#' Supported families are \code{gaussian} (default), \code{binomial}, and \code{poisson}. +#' @param coefPrefix Character. All model coefficient names will begin with this prefix. +#' default is \code{"beta_"} (so 'x' becomes 'beta_x', etc.) +#' @param sdPrefix Character. All dispersion parameters will begin with this prefix. #' default is no prefix. -#' @param priorSpecs List of prior specifications, should be generated using -#' setPriors() -#' @param modMatNames Logical, should parameters be named so they match the -#' names you would get from R's model.matrix function? +#' @param priors List of prior specifications, generated using +#' \code{setPriors()}. +#' @param modelMatrixNames Logical indicating if parameters be named so they match the +#' names one would get from R's \code{model.matrix}. +#' @param modelInfo Used internally by nimbleMacros; a list of model information such as constants and dimensions +#' @param .env Used internally by nimbleMacros; the environment where the model was created #' #' @author Ken Kellner #' @@ -35,7 +38,7 @@ #' LM(y ~ x + x2) #' }) #' -#' mod <- nimbleModel(code, constants=constants) +#' mod <- nimbleModel(code, constants = constants) #' mod$getCode() NULL @@ -60,7 +63,7 @@ LM <- list(process = function(code, modelInfo, .env){ link <- if(family$link == "identity") NULL else as.name(family$link) # Get modMatNames - modMatNames <- RHS$modMatNames + modMatNames <- RHS$modelMatrixNames if(is.null(modMatNames)) modMatNames <- FALSE # Create binomial sample size @@ -81,7 +84,7 @@ LM <- list(process = function(code, modelInfo, .env){ # Check that LHS is in the constants, otherwise this won't work LHS_name <- safeDeparse(LHS) if(! LHS_name %in% names(modelInfo$constants)){ - stop(paste0("If you don't provide dimensions for the data ", LHS_name, + stop(paste0("LM: If you don't provide dimensions for the data ", LHS_name, ", you must include ", LHS_name, " in the constants instead."), call.=FALSE) } idx <- substitute(1:LEN, list(LEN=as.numeric(length(modelInfo$constants[[LHS_name]])))) @@ -96,7 +99,7 @@ LM <- list(process = function(code, modelInfo, .env){ # RHS formula form <- RHS[[2]] - priorSpecs <- if(is.null(RHS$priorSpecs)) quote(setPriors()) else RHS$priorSpecs + priors <- if(is.null(RHS$priors)) quote(setPriors()) else RHS$priors coefPrefix <- if(is.null(RHS$coefPrefix)) quote(beta_) else RHS$coefPrefix sdPrefix <- RHS$sdPrefix if(family$family == "gaussian"){ @@ -110,16 +113,16 @@ LM <- list(process = function(code, modelInfo, .env){ dataDec <- getDataDistCode(family$family, LHS, idx, par2) # FIXME: LHS par should not be fixed at mu LP <- substitute(mu_[IDX] <- LINPRED(FORM, link=LINK, coefPrefix=COEFPREFIX, - sdPrefix=SDPREFIX, priorSpecs=PRIORS, - modMatNames=MODMAT), + sdPrefix=SDPREFIX, priors=PRIORS, + modelMatrixNames=MODMAT), list(IDX=idx, FORM=form, COEFPREFIX=coefPrefix, SDPREFIX=sdPrefix, - PRIORS=priorSpecs, LINK=link, MODMAT=modMatNames)) + PRIORS=priors, LINK=link, MODMAT=modMatNames)) pars_added <- list(quote(mu_)) out <- list(dataDec, LP) if (family$family == "gaussian"){ - priorSpecs <- eval(priorSpecs, envir=.env) - sdPrior <- matchPrior(par2, "sd", priorSpecs = priorSpecs) + priors <- eval(priors, envir=.env) + sdPrior <- matchPrior(par2, "sd", priors = priors) sigprior <- substitute(SDRES ~ SDPRIOR, list(SDPRIOR=sdPrior, SDRES=par2)) pars_added <- c(pars_added, list(par2)) out <- c(out, list(sigprior)) @@ -145,10 +148,10 @@ processFamily <- function(fam){ } if(is.null(fam$family)){ - stop("'family' not recognized") + stop("processFamily: 'family' not recognized: ", fam, call. = FALSE) } if(!fam$link %in% c("log", "identity","logit","probit","cloglog")){ - stop("Unsupported link function", call.=FALSE) + stop("processFamily: Unsupported link function: ", fam$link, call.=FALSE) } fam } @@ -167,7 +170,7 @@ getDataDistCode <- function(family, response, idx, dispParName){ out <- substitute(DAT ~ FORLOOP(DIST(mu_[IDX], size=N)), list(DAT=response, DIST=quote(dbinom), N=dispParName, IDX=idx)) } else { - stop("Family not supported", call.=FALSE) + stop("getDataDistCode: Family not supported: ", family, call.=FALSE) } out } diff --git a/R/formulaFunction.R b/R/formulaHandler.R similarity index 64% rename from R/formulaFunction.R rename to R/formulaHandler.R index d606d6d..c585be6 100644 --- a/R/formulaFunction.R +++ b/R/formulaHandler.R @@ -1,24 +1,23 @@ -#' Function to handle offset() in LINPRED -#' -#' Translates offset() in an R formula passed to LINPRED into corresponding -#' NIMBLE code for an offset in the linear predictor. This is used internally -#' by \code{LINPRED} and should not be called directly. New formula functions -#' should have the same arguments, naming structure, class (\code{nimbleFormulaHandler}) -#' and return object class (\code{formulaComponent}). -#' -#' @param x A \code{formulaComponentFunction} object created from an offset() term -#' @param defaultBracket The bracket from the LHS of LINPRED -#' @param coefPrefix The prefix to use for any new linear predictor parameters created -#' @param sdPrefix The prefix to use for any new standard deviation parameters created -#' @param modelInfo Named list containing model information including constants -#' @param env Environment in which the LINPRED macro was called -#' @param ... Not currently used -#' -#' @return An object of class \code{formulaComponent}. -#' -#' @author Ken Kellner -#' -#' @export +# Function to handle offset() in LINPRED +# +# Translates offset() in an R formula passed to LINPRED into corresponding +# NIMBLE code for an offset in the linear predictor. This is used internally +# by \code{LINPRED} and should not be called directly. New formula functions +# should have the same arguments, naming structure, class (\code{nimbleFormulaHandler}) +# and return object class (\code{formulaComponent}). +# +# @param x A \code{formulaComponentFunction} object created from an offset() term +# @param defaultBracket The bracket from the LHS of LINPRED +# @param coefPrefix The prefix to use for any new linear predictor parameters created +# @param sdPrefix The prefix to use for any new standard deviation parameters created +# @param modelInfo Named list containing model information including constants +# @param env Environment in which the LINPRED macro was called +# @param ... Not currently used +# +# @return An object of class \code{formulaComponent}. +# +# @author Ken Kellner +# formulaHandler_offset <- function(x, defaultBracket, coefPrefix, sdPrefix, modelInfo, env, ...){ # Get the code inside offset() @@ -75,29 +74,28 @@ get_all_names_recursive <- function(code){ NULL } -#' Function to handle scale() in LINPRED -#' -#' Translates scale() in an R formula passed to LINPRED into corresponding -#' NIMBLE code (and new constant) for a scaled covariate. This is used internally -#' by \code{LINPRED} and should not be called directly. New formula functions -#' should have the same arguments, naming structure, class (\code{nimbleFormulaHandler}) -#' and return object class (\code{formulaComponent}). Note: when applied to a -#' matrix or array covariate, scale() will calculate mean/SD relative to the entire -#' matrix/array, NOT column-wise as is the case if you use scale() in base R. -#' -#' @param x A formulaComponentFunction object created from a scale() term -#' @param defaultBracket The bracket from the LHS of LINPRED -#' @param coefPrefix The prefix to use for any new linear predictor parameters created -#' @param sdPrefix The prefix to use for any new standard deviation parameters created -#' @param modelInfo Named list containing model information including constants -#' @param env Environment in which the LINPRED macro was called -#' @param ... Not currently used -#' -#' @return An object of class \code{formulaComponentFixed}. -#' -#' @author Ken Kellner -#' -#' @export +# Function to handle scale() in LINPRED +# +# Translates scale() in an R formula passed to LINPRED into corresponding +# NIMBLE code (and new constant) for a scaled covariate. This is used internally +# by \code{LINPRED} and should not be called directly. New formula functions +# should have the same arguments, naming structure, class (\code{nimbleFormulaHandler}) +# and return object class (\code{formulaComponent}). Note: when applied to a +# matrix or array covariate, scale() will calculate mean/SD relative to the entire +# matrix/array, NOT column-wise as is the case if you use scale() in base R. +# +# @param x A formulaComponentFunction object created from a scale() term +# @param defaultBracket The bracket from the LHS of LINPRED +# @param coefPrefix The prefix to use for any new linear predictor parameters created +# @param sdPrefix The prefix to use for any new standard deviation parameters created +# @param modelInfo Named list containing model information including constants +# @param env Environment in which the LINPRED macro was called +# @param ... Not currently used +# +# @return An object of class \code{formulaComponentFixed}. +# +# @author Ken Kellner +# formulaHandler_scale <- function(x, defaultBracket, coefPrefix, sdPrefix, modelInfo, env, ...){ # Identify which interaction terms involve scale @@ -166,28 +164,27 @@ formulaHandler_scale <- function(x, defaultBracket, coefPrefix, sdPrefix, modelI class(formulaHandler_scale) <- c(class(formulaHandler_scale), "nimbleFormulaHandler") -#' Function to handle I() in LINPRED -#' -#' Translates I() in an R formula passed to LINPRED into corresponding -#' NIMBLE code (and new constant) for a scaled covariate. -#' Only allows for expressions involving one covariate (not functions of covariates). -#' This is used internally by \code{LINPRED} and should not be called directly. -#' New formula functions should have the same arguments, naming structure, class -#' (\code{nimbleFormulaHandler}) and return object class (\code{formulaComponent}). -#' -#' @param x A formulaComponentFunction object created from an I() term -#' @param defaultBracket The bracket from the LHS of LINPRED -#' @param coefPrefix The prefix to use for any new linear predictor parameters created -#' @param sdPrefix The prefix to use for any new standard deviation parameters created -#' @param modelInfo Named list containing model information including constants -#' @param env Environment in which the LINPRED macro was called -#' @param ... Not currently used -#' -#' @return An object of class \code{formulaComponentFixed}. -#' -#' @author Ken Kellner -#' -#' @export +# Function to handle I() in LINPRED +# +# Translates I() in an R formula passed to LINPRED into corresponding +# NIMBLE code (and new constant) for a scaled covariate. +# Only allows for expressions involving one covariate (not functions of covariates). +# This is used internally by \code{LINPRED} and should not be called directly. +# New formula functions should have the same arguments, naming structure, class +# (\code{nimbleFormulaHandler}) and return object class (\code{formulaComponent}). +# +# @param x A formulaComponentFunction object created from an I() term +# @param defaultBracket The bracket from the LHS of LINPRED +# @param coefPrefix The prefix to use for any new linear predictor parameters created +# @param sdPrefix The prefix to use for any new standard deviation parameters created +# @param modelInfo Named list containing model information including constants +# @param env Environment in which the LINPRED macro was called +# @param ... Not currently used +# +# @return An object of class \code{formulaComponentFixed}. +# +# @author Ken Kellner +# formulaHandler_I <- function(x, defaultBracket, coefPrefix, sdPrefix, modelInfo, env, ...){ # Identify which interaction terms involve scale @@ -250,27 +247,26 @@ formulaHandler_I <- function(x, defaultBracket, coefPrefix, sdPrefix, modelInfo, class(formulaHandler_I) <- c(class(formulaHandler_I), "nimbleFormulaHandler") -#' Function to handle log() in LINPRED -#' -#' Translates log() in an R formula passed to LINPRED into corresponding -#' NIMBLE code (and new constant) for a log-transformed covariate. This is used internally -#' by \code{LINPRED} and should not be called directly. New formula functions -#' should have the same arguments, naming structure, class (\code{nimbleFormulaHandler}) -#' and return object class (\code{formulaComponent}). -#' -#' @param x A formulaComponentFunction object created from a log() term -#' @param defaultBracket The bracket from the LHS of LINPRED -#' @param coefPrefix The prefix to use for any new linear predictor parameters created -#' @param sdPrefix The prefix to use for any new standard deviation parameters created -#' @param modelInfo Named list containing model information including constants -#' @param env Environment in which the LINPRED macro was called -#' @param ... Not currently used -#' -#' @return An object of class \code{formulaComponentFixed}. -#' -#' @author Ken Kellner -#' -#' @export +# Function to handle log() in LINPRED +# +# Translates log() in an R formula passed to LINPRED into corresponding +# NIMBLE code (and new constant) for a log-transformed covariate. This is used internally +# by \code{LINPRED} and should not be called directly. New formula functions +# should have the same arguments, naming structure, class (\code{nimbleFormulaHandler}) +# and return object class (\code{formulaComponent}). +# +# @param x A formulaComponentFunction object created from a log() term +# @param defaultBracket The bracket from the LHS of LINPRED +# @param coefPrefix The prefix to use for any new linear predictor parameters created +# @param sdPrefix The prefix to use for any new standard deviation parameters created +# @param modelInfo Named list containing model information including constants +# @param env Environment in which the LINPRED macro was called +# @param ... Not currently used +# +# @return An object of class \code{formulaComponentFixed}. +# +# @author Ken Kellner +# formulaHandler_log <- function(x, defaultBracket, coefPrefix, sdPrefix, modelInfo, env, ...){ # Identify which interaction terms involve scale diff --git a/R/setPriors.R b/R/setPriors.R index 49d68e1..91c767e 100644 --- a/R/setPriors.R +++ b/R/setPriors.R @@ -1,16 +1,9 @@ #' Set up prior values for different categories of nodes #' #' Generates a list of custom specifications of priors for parameters in the model. -#' It is possible to set priors for a category of parameters (e.g. intercept, +#' It is possible to set priors for a category of parameters (e.g., intercept, #' coefficient, sd, factor, continuous) or to set a prior for a specific #' parameter name (optionally including brackets with indices). -#' Exact name matches including brackets/indices are used first, followed by -#' name matches without indices, followed by data type (factor/continuous) -#' followed by parameter type (intercept/coefficient/sd). -#' Arguments can be supplied as quoted code, a character string, or -#' as a list of prior components. For example, to specify the prior -#' \code{dnorm(0, sd = 10)} you could specify \code{quote(dnorm(0, sd = 10))}, -#' or \code{"dnorm(0, sd = 10)"}, or \code{list("dnorm", 0, sd = 10)}. #' #' @name setPriors #' @author Ken Kellner @@ -22,12 +15,22 @@ #' to factor data #' @param continuous Prior specifications for slope coefficients corresponding #' to continuous data -#' @param eta Value of shape parameter for LKJ distribution prior, used for +#' @param lkjShape Value of shape parameter for LKJ distribution prior, used for #' correlation matrix in correlated random slope and intercept models #' @param ... Specific parameters, optionally with brackets/indices #' #' @seealso [nimble::dlkj_corr_cholesky] for more on the LKJ distribution #' +#' @details +#' +#' Exact name matches including brackets/indices are used first, followed by +#' name matches without indices, followed by data type (factor/continuous) +#' followed by parameter type (intercept/coefficient/sd). +#' Arguments can be supplied as quoted code, a character string, or +#' as a list of prior components. For example, to specify the prior +#' \code{dnorm(0, sd = 10)} you could specify \code{quote(dnorm(0, sd = 10))}, +#' or \code{"dnorm(0, sd = 10)"}, or \code{list("dnorm", 0, sd = 10)}. +#' #' @author Ken Kellner #' #' @examples @@ -46,7 +49,7 @@ setPriors <- function(intercept = quote(dnorm(0, sd = 1000)), sd = quote(dunif(0, 100)), factor = NULL, continuous = NULL, - eta = 1.3, + lkjShape = 1, ...){ # Get specific prior names extra <- list(...) @@ -57,7 +60,7 @@ setPriors <- function(intercept = quote(dnorm(0, sd = 1000)), # Check if any outputs are numeric (indicating they were passed without quote()) lapply(1:length(out), function(i) if(is.numeric(out[[i]])){ - stop("Prior code for ", names(out)[i], " should be wrapped in quote()", call.=FALSE) + stop("setPriors: Prior code for ", names(out)[i], " should be wrapped in quote()", call.=FALSE) }) # Remove any null values @@ -72,7 +75,7 @@ setPriors <- function(intercept = quote(dnorm(0, sd = 1000)), }) # Add eta/shape parameter (adding here because we don't want to check it above) - out$eta <- eta + out$lkjShape <- lkjShape out } @@ -82,7 +85,7 @@ setPriors <- function(intercept = quote(dnorm(0, sd = 1000)), convertListToPrior <- function(input){ if(!is.list(input)) return(input) if(is.function(input[[1]])){ - stop("Function name must be a character string or wrapped in quote()", call.=FALSE) + stop("convertListToPrior: Function name must be a character string or wrapped in quote()", call.=FALSE) } out <- input @@ -101,7 +104,7 @@ removeBracket <- function(node){ #' Match a prior from a list of prior settings #' #' Attempts to determine which prior to put on a parameter based on a list of settings, -#' such as the output from setPriors(). The function follows the following +#' such as the output from \code{setPriors()}. The function follows the following #' search pattern: (1) looks for an exact match to the parameter name including brackets; #' (2) a match to the parameter name without brackets; (3) goes through each value #' supplied to \code{...} in order and looks for a match in the names of the @@ -114,9 +117,9 @@ removeBracket <- function(node){ #' @param parName Parameter to get a prior for, as quoted code/name, possibly #' including brackets/indices #' @param ... Character strings that categorize the parameter and match -#' the names of elements in \code{priorSpecs}. The order is important: +#' the names of elements in \code{priors}. The order is important: #' the first match found is used. -#' @param priorSpecs A named list of prior settings, e.g. as generated by +#' @param priors A named list of prior settings, e.g., as generated by #' \code{setPriors} #' #' @author Ken Kellner @@ -124,20 +127,20 @@ removeBracket <- function(node){ #' @examples #' pr <- setPriors(intercept = quote(dunif(-3, 3)), 'alpha' = quote(dunif(0,1)), #' 'alpha[2]' = "dnorm(0, 3)") -#' matchPrior(quote(alpha), priorSpecs=pr) -#' matchPrior(quote(alpha[2]), priorSpecs=pr) -#' matchPrior(quote(intercept), priorSpecs=pr) +#' matchPrior(quote(alpha), priors=pr) +#' matchPrior(quote(alpha[2]), priors=pr) +#' matchPrior(quote(intercept), priors=pr) #' #' @export -matchPrior <- function(parName, ..., priorSpecs){ +matchPrior <- function(parName, ..., priors){ - if(is.call(priorSpecs) | is.name(priorSpecs)){ - priorSpecs <- eval(priorSpecs) + if(is.call(priors) | is.name(priors)){ + priors <- eval(priors) } # 1. If exact prior name is specified in priors par_char <- safeDeparse(parName) - possible_prior <- priorSpecs[[par_char]] + possible_prior <- priors[[par_char]] if(!is.null(possible_prior)){ checkValidPrior(possible_prior) return(possible_prior) @@ -145,7 +148,7 @@ matchPrior <- function(parName, ..., priorSpecs){ # 2. If prior name without brackets is specified par_nobracks <- safeDeparse(removeBracket(parName)) - possible_prior <- priorSpecs[[par_nobracks]] + possible_prior <- priors[[par_nobracks]] if(!is.null(possible_prior)){ checkValidPrior(possible_prior) return(possible_prior) @@ -155,16 +158,16 @@ matchPrior <- function(parName, ..., priorSpecs){ par_types <- list(...) for (i in par_types){ if(is.null(i) || is.na(i)) next - if(i %in% names(priorSpecs)){ - if(is.null(priorSpecs[[i]])) next - checkValidPrior(priorSpecs[[i]]) - return(priorSpecs[[i]]) + if(i %in% names(priors)){ + if(is.null(priors[[i]])) next + checkValidPrior(priors[[i]]) + return(priors[[i]]) } } - stop("Unable to set prior for parameter ", parName, call.=FALSE) + stop("matchPrior: Unable to set prior for parameter: ", parName, call.=FALSE) } checkValidPrior <- function(prior){ - if(!(is.call(prior) | is.name(prior))) stop("Invalid prior ", prior, call.=FALSE) + if(!(is.call(prior) | is.name(prior))) stop("Invalid prior: ", prior, call.=FALSE) invisible() } diff --git a/README.md b/README.md index 4581a42..1eddf98 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,7 @@ status](https://github.com/nimble-dev/nimbleMacros/workflows/R-CMD-check/badge.svg)](https://github.com/nimble-dev/nimbleMacros/actions) `nimbleMacros` is an R package that extends [NIMBLE](https://r-nimble.org/), an R package for hierarchical statistical modeling. -The package provides a set of model macros, which are special code chunks which NIMBLE automatically expands into larger blocks of valid NIMBLE model code. +The package provides a set of model macros, which are special code chunks that NIMBLE automatically expands into larger blocks of valid NIMBLE model code. For example, the following NIMBLE model code contains a macro called `LINPRED`, and the objective is to create a linear predictor containing an intercept and a slope for one covariate, `x`. diff --git a/man/FORLOOP.Rd b/man/FORLOOP.Rd index 4fb2f1d..9f88294 100644 --- a/man/FORLOOP.Rd +++ b/man/FORLOOP.Rd @@ -5,6 +5,10 @@ \title{Macro to build for loop(s) from code with index ranges in brackets} \arguments{ \item{code}{The right-hand side of a parameter declaration} + +\item{modelInfo}{Used internally by nimbleMacros; a list of model information such as constants and dimensions} + +\item{.env}{Used internally by nimbleMacros; the environment where the model was created} } \description{ This macro takes a line of NIMBLE model code with index ranges inside brackets @@ -18,7 +22,7 @@ code <- nimbleCode({ mu[1:n] <- FORLOOP(beta[1] + beta[2]*x[1:n]) }) -mod <- nimbleModel(code, constants=list(n=10)) +mod <- nimbleModel(code, constants = list(n=10)) mod$getCode() } \author{ diff --git a/man/LINPRED.Rd b/man/LINPRED.Rd index 8bd5b5d..6fe8621 100644 --- a/man/LINPRED.Rd +++ b/man/LINPRED.Rd @@ -7,27 +7,44 @@ \item{formula}{An R formula, possibly with the parameters followed by brackets containing indices. If there are no indices, the macro attempts to guess the correct indices from the context. The formula must be -right-hand side only (e.g. ~x). This must always be the first argument supplied -to LINPRED.} +right-hand side only (e.g. \code{~x}). This must always be the first argument supplied +to \code{LINPRED}.} -\item{link}{A link function which will be applied to the +\item{link}{A link function that will be applied to the left-hand-side (the response) in the final linear predictor. Default is none.} \item{coefPrefix}{All model coefficient names will begin with this prefix. -default is beta_ (so x becomes beta_x, etc.)} +default is \code{"beta_"} (so 'x' becomes 'beta_x', etc.)} \item{sdPrefix}{All dispersion parameters will begin with this prefix. -default is no prefix.} +Default is no prefix.} -\item{priorSpecs}{Prior specifications, should be generated with setPrior()} +\item{priors}{Prior specifications, generated with \code{setPrior()}} -\item{modMatNames}{Logical, should parameters be named so they match the -names you would get from R's model.matrix function?} +\item{modelMatrixNames}{Logical indicating if parameters should be named so they match the +names one would get from R's \code{model.matrix}. Default is \code{FALSE}.} -\item{noncenter}{Logical; use noncentered parameterization?} +\item{noncentered}{Logical indicating whether to use noncentered parameterization +for random effects. Default is \code{FALSE}. Under the noncentered parameterization, +random effects have a standard normal prior (\code{beta_x_raw ~ dnorm(0, sd = 1)}) +and are then scaled by the hyperparameters (mean and SD), i.e., +\code{beta_x = mu_beta + sd_beta * beta_x_raw}. This parameterization can improve +mixing in some cases.} -\item{centerVar}{Grouping covariate to 'center' on in parameterization. By -default all random effects have mean 0 as with lme4.} +\item{centerVar}{Grouping variable (covariate) to 'center' the random effects on. +By default (when NULL), random effects come from normal distributions with mean 0 +as with \code{lme4}. For example, for random intercepts by grouping variable \code{x}, +the linear predictor would be \code{beta_Intercept + beta_x[x[i]]} and the +prior for the random effects would be \code{beta_x ~ dnorm(0, sd_x)}. When +\code{centerVar = x}, the linear predictor would be \code{beta_x[x[i]]} +and the random effect prior would be \code{beta_x ~ dnorm(beta_Intercept, sd = sd_x)}. +That is, the mean of the random effects is now \code{beta_Intercept}. +These two formulations should yield the same results. Note that this option +is unrelated to the \code{noncentered} argument despite the similar names.} + +\item{modelInfo}{Used internally by nimbleMacros; a list of model information such as constants and dimensions} + +\item{.env}{Used internally by nimbleMacros; the environment where the model was created} } \description{ Converts an R formula into corresponding code for a linear predictor in NIMBLE model code. @@ -40,7 +57,7 @@ code <- nimbleCode({ mu[1:3] <- LINPRED(~x + x2) }) -mod <- nimbleModel(code, constants=constants) +mod <- nimbleModel(code, constants = constants) mod$getCode() } \author{ diff --git a/man/LINPRED_PRIORS.Rd b/man/LINPRED_PRIORS.Rd index d95ef51..4a2d21d 100644 --- a/man/LINPRED_PRIORS.Rd +++ b/man/LINPRED_PRIORS.Rd @@ -2,32 +2,49 @@ % Please edit documentation in R/LINPRED.R \name{LINPRED_PRIORS} \alias{LINPRED_PRIORS} -\title{Macro to build code for priors on a linear predictor from R formula} +\title{Macro to build code for priors on a linear predictor from an R formula} \arguments{ -\item{formula}{An R formula The formula must be right-hand side only (e.g. ~x). -This must always be the first argument supplied to LINPRED_PRIORS} +\item{formula}{An R formula The formula must be right-hand side only (e.g., \code{~x}). +This must always be the first argument supplied to \code{LINPRED_PRIORS}.} \item{coefPrefix}{All model coefficient names will begin with this prefix. -default is beta_ (so x becomes beta_x, etc.)} +default is \code{"beta_"} (so 'x' becomes 'beta_x', etc.)} \item{sdPrefix}{All dispersion parameters will begin with this prefix. -default is no prefix.} +Default is no prefix.} -\item{priorSpecs}{List of prior specifications, should be generated using +\item{priors}{List of prior specifications, generated using \code{setPriors}. setPriors()} -\item{modMatNames}{Logical, should parameters be named so they match the -names you would get from R's model.matrix function?} +\item{modelMatrixNames}{Logical indicating if parameters should be named so they match the +names one would get from R's \code{model.matrix}. Default is \code{FALSE}.} -\item{noncenter}{Logical, use noncentered parameterization?} +\item{noncentered}{Logical indicating whether to use noncentered parameterization +for random effects. Default is \code{FALSE}. Under the noncentered parameterization, +random effects have a standard normal prior (\code{beta_x_raw ~ dnorm(0, sd = 1)}) +and are then scaled by the hyperparameters (mean and SD), i.e., +\code{beta_x = mu_beta + sd_beta * beta_x_raw}. This parameterization can improve +mixing in some cases.} -\item{centerVar}{Grouping covariate to 'center' on in parameterization. By -default all random effects have mean 0 as with lme4.} +\item{centerVar}{Grouping variable (covariate) to 'center' the random effects on. +By default (when NULL), random effects come from normal distributions with mean 0 +as with \code{lme4}. For example, for random intercepts by grouping variable \code{x}, +the linear predictor would be \code{beta_Intercept + beta_x[x[i]]} and the +prior for the random effects would be \code{beta_x ~ dnorm(0, sd_x)}. When +\code{centerVar = x}, the linear predictor would be \code{beta_x[x[i]]} +and the random effect prior would be \code{beta_x ~ dnorm(beta_Intercept, sd = sd_x)}. +That is, the mean of the random effects is now \code{beta_Intercept}. +These two formulations should yield the same results. Note that this option +is unrelated to the \code{noncentered} argument despite the similar names.} + +\item{modelInfo}{Used internally by nimbleMacros; a list of model information such as constants and dimensions} + +\item{.env}{Used internally by nimbleMacros; the environment where the model was created} } \description{ Generates appropriate priors for a linear predictor derived from an R formula. As such it makes the most sense to use this macro together with -the LINPRED macro which takes similar arguments. +the LINPRED macro, which takes similar arguments. } \examples{ constants <- list(x = rnorm(3), x2 = factor(letters[1:3])) diff --git a/man/LM.Rd b/man/LM.Rd index f3f463f..bb56abc 100644 --- a/man/LM.Rd +++ b/man/LM.Rd @@ -7,28 +7,33 @@ \item{formula}{An R formula, possibly with the parameters followed by brackets containing indices. If there are no indices, the macro attempts to guess the correct indices from the context. Formulas can include -random effects via lme4-style notation (e.g. ~ x + (1|group))} +random effects via lme4-style notation (e.g., \code{~ x + (1|group)})} \item{family}{A description of the error distribution and link function to be used in the model. This can be a character string naming a family function, -a family function or the result of a call to a family function. See ?family. -Supported families are gaussian (default), binomial, and poisson.} +a family function or the result of a call to a family function. See \code{?family}. +Supported families are \code{gaussian} (default), \code{binomial}, and \code{poisson}.} -\item{coefPrefix}{All model coefficient names will begin with this prefix. -default is beta_ (so x becomes beta_x, etc.)} +\item{coefPrefix}{Character. All model coefficient names will begin with this prefix. +default is \code{"beta_"} (so 'x' becomes 'beta_x', etc.)} -\item{sdPrefix}{All dispersion parameters will begin with this prefix. +\item{sdPrefix}{Character. All dispersion parameters will begin with this prefix. default is no prefix.} -\item{priorSpecs}{List of prior specifications, should be generated using -setPriors()} +\item{priors}{List of prior specifications, generated using +\code{setPriors()}.} -\item{modMatNames}{Logical, should parameters be named so they match the -names you would get from R's model.matrix function?} +\item{modelMatrixNames}{Logical indicating if parameters be named so they match the +names one would get from R's \code{model.matrix}.} + +\item{modelInfo}{Used internally by nimbleMacros; a list of model information such as constants and dimensions} + +\item{.env}{Used internally by nimbleMacros; the environment where the model was created} } \description{ This macro generates code for LMs, GLMs, and GLMMs using formula notation -and arguments similar to R functions such as lm(), glm(), and lmer()/glmer(). +and arguments similar to R functions such as \code{lm()}, \code{glm()}, +and \code{lmer()}/\code{glmer()}. Currently only normal, Poisson, and binomial models are supported. } \examples{ @@ -40,7 +45,7 @@ code <- nimbleCode({ LM(y ~ x + x2) }) -mod <- nimbleModel(code, constants=constants) +mod <- nimbleModel(code, constants = constants) mod$getCode() } \author{ diff --git a/man/formulaHandler_I.Rd b/man/formulaHandler_I.Rd deleted file mode 100644 index 9a5c9b3..0000000 --- a/man/formulaHandler_I.Rd +++ /dev/null @@ -1,37 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/formulaFunction.R -\name{formulaHandler_I} -\alias{formulaHandler_I} -\title{Function to handle I() in LINPRED} -\usage{ -formulaHandler_I(x, defaultBracket, coefPrefix, sdPrefix, modelInfo, env, ...) -} -\arguments{ -\item{x}{A formulaComponentFunction object created from an I() term} - -\item{defaultBracket}{The bracket from the LHS of LINPRED} - -\item{coefPrefix}{The prefix to use for any new linear predictor parameters created} - -\item{sdPrefix}{The prefix to use for any new standard deviation parameters created} - -\item{modelInfo}{Named list containing model information including constants} - -\item{env}{Environment in which the LINPRED macro was called} - -\item{...}{Not currently used} -} -\value{ -An object of class \code{formulaComponentFixed}. -} -\description{ -Translates I() in an R formula passed to LINPRED into corresponding -NIMBLE code (and new constant) for a scaled covariate. -Only allows for expressions involving one covariate (not functions of covariates). -This is used internally by \code{LINPRED} and should not be called directly. -New formula functions should have the same arguments, naming structure, class -(\code{nimbleFormulaHandler}) and return object class (\code{formulaComponent}). -} -\author{ -Ken Kellner -} diff --git a/man/formulaHandler_log.Rd b/man/formulaHandler_log.Rd deleted file mode 100644 index 73e02f7..0000000 --- a/man/formulaHandler_log.Rd +++ /dev/null @@ -1,44 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/formulaFunction.R -\name{formulaHandler_log} -\alias{formulaHandler_log} -\title{Function to handle log() in LINPRED} -\usage{ -formulaHandler_log( - x, - defaultBracket, - coefPrefix, - sdPrefix, - modelInfo, - env, - ... -) -} -\arguments{ -\item{x}{A formulaComponentFunction object created from a log() term} - -\item{defaultBracket}{The bracket from the LHS of LINPRED} - -\item{coefPrefix}{The prefix to use for any new linear predictor parameters created} - -\item{sdPrefix}{The prefix to use for any new standard deviation parameters created} - -\item{modelInfo}{Named list containing model information including constants} - -\item{env}{Environment in which the LINPRED macro was called} - -\item{...}{Not currently used} -} -\value{ -An object of class \code{formulaComponentFixed}. -} -\description{ -Translates log() in an R formula passed to LINPRED into corresponding -NIMBLE code (and new constant) for a log-transformed covariate. This is used internally -by \code{LINPRED} and should not be called directly. New formula functions -should have the same arguments, naming structure, class (\code{nimbleFormulaHandler}) -and return object class (\code{formulaComponent}). -} -\author{ -Ken Kellner -} diff --git a/man/formulaHandler_offset.Rd b/man/formulaHandler_offset.Rd deleted file mode 100644 index b362c54..0000000 --- a/man/formulaHandler_offset.Rd +++ /dev/null @@ -1,44 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/formulaFunction.R -\name{formulaHandler_offset} -\alias{formulaHandler_offset} -\title{Function to handle offset() in LINPRED} -\usage{ -formulaHandler_offset( - x, - defaultBracket, - coefPrefix, - sdPrefix, - modelInfo, - env, - ... -) -} -\arguments{ -\item{x}{A \code{formulaComponentFunction} object created from an offset() term} - -\item{defaultBracket}{The bracket from the LHS of LINPRED} - -\item{coefPrefix}{The prefix to use for any new linear predictor parameters created} - -\item{sdPrefix}{The prefix to use for any new standard deviation parameters created} - -\item{modelInfo}{Named list containing model information including constants} - -\item{env}{Environment in which the LINPRED macro was called} - -\item{...}{Not currently used} -} -\value{ -An object of class \code{formulaComponent}. -} -\description{ -Translates offset() in an R formula passed to LINPRED into corresponding -NIMBLE code for an offset in the linear predictor. This is used internally -by \code{LINPRED} and should not be called directly. New formula functions -should have the same arguments, naming structure, class (\code{nimbleFormulaHandler}) -and return object class (\code{formulaComponent}). -} -\author{ -Ken Kellner -} diff --git a/man/formulaHandler_scale.Rd b/man/formulaHandler_scale.Rd deleted file mode 100644 index 6840f2d..0000000 --- a/man/formulaHandler_scale.Rd +++ /dev/null @@ -1,46 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/formulaFunction.R -\name{formulaHandler_scale} -\alias{formulaHandler_scale} -\title{Function to handle scale() in LINPRED} -\usage{ -formulaHandler_scale( - x, - defaultBracket, - coefPrefix, - sdPrefix, - modelInfo, - env, - ... -) -} -\arguments{ -\item{x}{A formulaComponentFunction object created from a scale() term} - -\item{defaultBracket}{The bracket from the LHS of LINPRED} - -\item{coefPrefix}{The prefix to use for any new linear predictor parameters created} - -\item{sdPrefix}{The prefix to use for any new standard deviation parameters created} - -\item{modelInfo}{Named list containing model information including constants} - -\item{env}{Environment in which the LINPRED macro was called} - -\item{...}{Not currently used} -} -\value{ -An object of class \code{formulaComponentFixed}. -} -\description{ -Translates scale() in an R formula passed to LINPRED into corresponding -NIMBLE code (and new constant) for a scaled covariate. This is used internally -by \code{LINPRED} and should not be called directly. New formula functions -should have the same arguments, naming structure, class (\code{nimbleFormulaHandler}) -and return object class (\code{formulaComponent}). Note: when applied to a -matrix or array covariate, scale() will calculate mean/SD relative to the entire -matrix/array, NOT column-wise as is the case if you use scale() in base R. -} -\author{ -Ken Kellner -} diff --git a/man/matchPrior.Rd b/man/matchPrior.Rd index 7e3d73b..334225e 100644 --- a/man/matchPrior.Rd +++ b/man/matchPrior.Rd @@ -4,22 +4,22 @@ \alias{matchPrior} \title{Match a prior from a list of prior settings} \usage{ -matchPrior(parName, ..., priorSpecs) +matchPrior(parName, ..., priors) } \arguments{ \item{parName}{Parameter to get a prior for, as quoted code/name, possibly including brackets/indices} \item{...}{Character strings that categorize the parameter and match -the names of elements in \code{priorSpecs}. The order is important: +the names of elements in \code{priors}. The order is important: the first match found is used.} -\item{priorSpecs}{A named list of prior settings, e.g. as generated by +\item{priors}{A named list of prior settings, e.g., as generated by \code{setPriors}} } \description{ Attempts to determine which prior to put on a parameter based on a list of settings, -such as the output from setPriors(). The function follows the following +such as the output from \code{setPriors()}. The function follows the following search pattern: (1) looks for an exact match to the parameter name including brackets; (2) a match to the parameter name without brackets; (3) goes through each value supplied to \code{...} in order and looks for a match in the names of the @@ -29,9 +29,9 @@ prior value. \examples{ pr <- setPriors(intercept = quote(dunif(-3, 3)), 'alpha' = quote(dunif(0,1)), 'alpha[2]' = "dnorm(0, 3)") -matchPrior(quote(alpha), priorSpecs=pr) -matchPrior(quote(alpha[2]), priorSpecs=pr) -matchPrior(quote(intercept), priorSpecs=pr) +matchPrior(quote(alpha), priors=pr) +matchPrior(quote(alpha[2]), priors=pr) +matchPrior(quote(intercept), priors=pr) } \author{ diff --git a/man/setPriors.Rd b/man/setPriors.Rd index e7670f3..299765d 100644 --- a/man/setPriors.Rd +++ b/man/setPriors.Rd @@ -10,7 +10,7 @@ setPriors( sd = quote(dunif(0, 100)), factor = NULL, continuous = NULL, - eta = 1.3, + lkjShape = 1, ... ) } @@ -27,16 +27,18 @@ to factor data} \item{continuous}{Prior specifications for slope coefficients corresponding to continuous data} -\item{eta}{Value of shape parameter for LKJ distribution prior, used for +\item{lkjShape}{Value of shape parameter for LKJ distribution prior, used for correlation matrix in correlated random slope and intercept models} \item{...}{Specific parameters, optionally with brackets/indices} } \description{ Generates a list of custom specifications of priors for parameters in the model. -It is possible to set priors for a category of parameters (e.g. intercept, +It is possible to set priors for a category of parameters (e.g., intercept, coefficient, sd, factor, continuous) or to set a prior for a specific parameter name (optionally including brackets with indices). +} +\details{ Exact name matches including brackets/indices are used first, followed by name matches without indices, followed by data type (factor/continuous) followed by parameter type (intercept/coefficient/sd). diff --git a/man/uppertri_mult_diag.Rd b/man/uppertri_mult_diag.Rd index 268d7aa..1da2a5d 100644 --- a/man/uppertri_mult_diag.Rd +++ b/man/uppertri_mult_diag.Rd @@ -7,14 +7,14 @@ uppertri_mult_diag(mat, vec) } \arguments{ -\item{mat}{upper triangular Cholesky factor of correlation matrix (Ustar)} +\item{mat}{upper triangular Cholesky factor of correlation matrix ("Ustar")} \item{vec}{vector of standard deviations for individual random effects} } \description{ nimbleFunction needed when fitting correlated random effects. -Generates upper triangular Cholesky factor of covariance matrix (U in code) -from upper tri Cholesky factor of correlation matrix (Ustar in code) +Generates upper triangular Cholesky factor of covariance matrix ("U" in code) +from upper triangular Cholesky factor of correlation matrix ("Ustar" in code) and vector of standard deviations. Taken from the NIMBLE manual, -section 5.2.4.1.2 LKJ distribution for correlation matrices. +section 5.2.4.1.2 (LKJ Distribution for Correlation Matrices). } diff --git a/tests/testthat/test_LINPRED.R b/tests/testthat/test_LINPRED.R index cc3ebb7..04f101c 100644 --- a/tests/testthat/test_LINPRED.R +++ b/tests/testthat/test_LINPRED.R @@ -118,7 +118,7 @@ test_that("LINPRED basic fixed effects models", { # Set custom priors priors <- setPriors(intercept="dnorm(0, sd=1)", coefficient="dnorm(0, sd=2)") code <- nimbleCode({ - mu[1:n] <- LINPRED(~x3, priorSpecs=priors) + mu[1:n] <- LINPRED(~x3, priors=priors) }) mod <- nimbleModel(code, constants = modInfo$constants) @@ -135,7 +135,7 @@ test_that("LINPRED basic fixed effects models", { # No priors code <- nimbleCode({ - mu[1:n] <- LINPRED(~x3, priorSpecs=NULL) + mu[1:n] <- LINPRED(~x3, priors=NULL) }) mod <- nimbleModel(code, constants = modInfo$constants) @@ -150,7 +150,7 @@ test_that("LINPRED basic fixed effects models", { # Modmatnames = TRUE code <- nimbleCode({ - mu[1:n] <- LINPRED(~x + x3, modMatNames = TRUE) + mu[1:n] <- LINPRED(~x + x3, modelMatrixNames = TRUE) }) mod <- nimbleModel(code, constants = modInfo$constants) @@ -223,9 +223,9 @@ test_that("LINPRED basic fixed effects models", { beta_x3 = 0, beta_x_x3 = c(0, 0, 0)) ) - # continuous-factor interaction with modMatNames = TRUE + # continuous-factor interaction with modelMatrixNames = TRUE code <- nimbleCode({ - mu[1:n] <- LINPRED(~x*x3, modMatNames = TRUE) + mu[1:n] <- LINPRED(~x*x3, modelMatrixNames = TRUE) }) mod <- nimbleModel(code, constants = modInfo$constants) @@ -274,9 +274,9 @@ test_that("LINPRED basic fixed effects models", { list(beta_x3 = 0, beta_x_x3 = c(0, 0, 0)) ) - # continuous-factor interaction with modMatNames = TRUE + # continuous-factor interaction with modelMatrixNames = TRUE code <- nimbleCode({ - mu[1:n] <- LINPRED(~x*x3, modMatNames = TRUE) + mu[1:n] <- LINPRED(~x*x3, modelMatrixNames = TRUE) }) mod <- nimbleModel(code, constants = modInfo$constants) @@ -340,9 +340,9 @@ test_that("LINPRED basic fixed effects models", { 0, 0, 0, 0, 0), dim = 3:2)) ) - # Factor-factor interaction with modMatNames = TRUE + # Factor-factor interaction with modelMatrixNames = TRUE code <- nimbleCode({ - mu[1:n] <- LINPRED(~x*x2, modMatNames = TRUE) + mu[1:n] <- LINPRED(~x*x2, modelMatrixNames = TRUE) }) mod <- nimbleModel(code, constants = modInfo$constants) @@ -378,7 +378,7 @@ test_that("LINPRED basic fixed effects models", { # Covariate not in constants works (but is assumed to be continuous) code <- nimbleCode({ - mu[1:n] <- LINPRED(~x5, modMatNames = TRUE) + mu[1:n] <- LINPRED(~x5, modelMatrixNames = TRUE) }) mod <- nimbleModel(code, constants = modInfo$constants) mod$getCode() @@ -386,6 +386,29 @@ test_that("LINPRED basic fixed effects models", { nimbleOptions(enableMacroComments = TRUE) }) +test_that("LINPRED error traps LHS in formula", { + nimbleOptions(enableMacroComments = FALSE) + set.seed(123) + modInfo <- list(constants=list(y = rnorm(10), x=round(rnorm(10), 3), n = 10)) + + # Missing LHS error trapping must be handled by buildMacro() in nimble + # Enable the test below when the use3pieces error trap PR in nimble is merged + #code <- quote(LINPRED(~1)) + + #expect_error( + # LINPRED$process(code, modelInfo=modInfo, environment()), + # "This macro must be used as part of an assignment" + #) + + # LHS in formula + code <- quote(mu[1:n] <- LINPRED(y~1)) + + expect_error( + LINPRED$process(code, modelInfo=modInfo, environment()), + "Formula should be RHS-only" + ) + +}) test_that("LINPRED with uncorrelated random effects", { nimbleOptions(enableMacroComments = FALSE) @@ -450,7 +473,7 @@ test_that("LINPRED with uncorrelated random effects", { # Set custom priors pr <- setPriors(intercept="dnorm(0, sd=1)", sd="dunif(0, 1)") code <- nimbleCode({ - mu[1:n] <- LINPRED(~x + (1|group), priorSpecs=pr) + mu[1:n] <- LINPRED(~x + (1|group), priors=pr) }) mod <- nimbleModel(code, constants = modInfo$constants) @@ -880,7 +903,7 @@ test_that("noncentered parameterization with uncorrelated random effects", { x4 = factor(sample(letters[6:8], 10, replace=T)), n=10)) code <- nimbleCode({ - mu[1:n] <- LINPRED(~x + (x||group), noncenter=TRUE) + mu[1:n] <- LINPRED(~x + (x||group), noncentered=TRUE) }) mod <- nimbleModel(code, constants=modInfo$constants) @@ -918,7 +941,7 @@ test_that("noncentered parameterization with uncorrelated random effects", { # With centering variable also code <- nimbleCode({ - mu[1:n] <- LINPRED(~x + (x||group), noncenter=TRUE, centerVar=group) + mu[1:n] <- LINPRED(~x + (x||group), noncentered=TRUE, centerVar=group) }) mod <- nimbleModel(code, constants=modInfo$constants) @@ -956,7 +979,7 @@ test_that("noncentered parameterization with uncorrelated random effects", { # Factor slope code <- nimbleCode({ - mu[1:n] <- LINPRED(~x + (x2||group), noncenter=TRUE) + mu[1:n] <- LINPRED(~x + (x2||group), noncentered=TRUE) }) mod <- nimbleModel(code, constants=modInfo$constants) @@ -1031,7 +1054,7 @@ test_that("correlated random effects", { sd_x_group ~ dunif(0, 100) re_sds_group[1] <- sd_group re_sds_group[2] <- sd_x_group - Ustar_group[1:2, 1:2] ~ dlkj_corr_cholesky(1.3, 2) + Ustar_group[1:2, 1:2] ~ dlkj_corr_cholesky(1, 2) U_group[1:2, 1:2] <- uppertri_mult_diag(Ustar_group[1:2, 1:2], re_sds_group[1:2]) re_means_group[1] <- 0 @@ -1050,10 +1073,10 @@ test_that("correlated random effects", { beta_x_group = c(0, 0, 0), sd_group = 1, sd_x_group = 1) ) - # Set eta value - pr <- setPriors(eta=3) + # Set LKJ shape value + pr <- setPriors(lkjShape=3) code <- nimbleCode({ - mu[1:n] <- LINPRED(~x + (x|group), priorSpecs=pr) + mu[1:n] <- LINPRED(~x + (x|group), priors=pr) }) mod <- nimbleModel(code, constants=modInfo$constants) @@ -1106,7 +1129,7 @@ test_that("correlated random effects", { re_sds_group[1] <- sd_group re_sds_group[2] <- sd_x_group re_sds_group[3] <- sd_x3_group - Ustar_group[1:3, 1:3] ~ dlkj_corr_cholesky(1.3, 3) + Ustar_group[1:3, 1:3] ~ dlkj_corr_cholesky(1, 3) U_group[1:3, 1:3] <- uppertri_mult_diag(Ustar_group[1:3, 1:3], re_sds_group[1:3]) re_means_group[1] <- 0 @@ -1161,7 +1184,7 @@ test_that("Centering with correlated random effects", { sd_x_group ~ dunif(0, 100) re_sds_group[1] <- sd_group re_sds_group[2] <- sd_x_group - Ustar_group[1:2, 1:2] ~ dlkj_corr_cholesky(1.3, 2) + Ustar_group[1:2, 1:2] ~ dlkj_corr_cholesky(1, 2) U_group[1:2, 1:2] <- uppertri_mult_diag(Ustar_group[1:2, 1:2], re_sds_group[1:2]) re_means_group[1] <- beta_Intercept @@ -1198,7 +1221,7 @@ test_that("Centering with correlated random effects", { sd_x_group ~ dunif(0, 100) re_sds_group[1] <- sd_group re_sds_group[2] <- sd_x_group - Ustar_group[1:2, 1:2] ~ dlkj_corr_cholesky(1.3, 2) + Ustar_group[1:2, 1:2] ~ dlkj_corr_cholesky(1, 2) U_group[1:2, 1:2] <- uppertri_mult_diag(Ustar_group[1:2, 1:2], re_sds_group[1:2]) re_means_group[1] <- 0 @@ -1228,7 +1251,7 @@ test_that("Noncentered parameterization doesn't work with correlated random effe x=round(rnorm(10),3), x3=round(rnorm(10), 3), x4 = factor(sample(letters[6:8], 10, replace=T)), n=10)) # Noncentered doesn't work with correlated random effects - code <- quote(LINPRED_PRIORS(~x + (x|group), noncenter=TRUE)) + code <- quote(LINPRED_PRIORS(~x + (x|group), noncentered=TRUE)) expect_error(LINPRED_PRIORS$process(code, modInfo, NULL)$code, "Noncentered") nimbleOptions(enableMacroComments = TRUE) @@ -1379,7 +1402,7 @@ test_that("Nested random effects", { # w:x notation # Linear predictor - code <- quote(mu[1:n] <- LINPRED(~x3 + (x3||x:w), priorSpecs=NULL)) + code <- quote(mu[1:n] <- LINPRED(~x3 + (x3||x:w), priors=NULL)) out <- LINPRED$process(code, modelInfo=modInfo, NULL) # Make sure new combined levels constant is added expect_equal( diff --git a/tests/testthat/test_LM.R b/tests/testthat/test_LM.R index 9d6eeae..9c2341e 100644 --- a/tests/testthat/test_LM.R +++ b/tests/testthat/test_LM.R @@ -8,7 +8,7 @@ test_that("LM for linear regression", { modelInfo <- list(constants=dat) code <- nimbleCode({ - LM(y ~ x + x2, priorSpecs=setPriors(sd="dunif(0, 5)")) + LM(y ~ x + x2, priors=setPriors(sd="dunif(0, 5)")) }) mod <- nimbleModel(code, constants=dat) @@ -36,11 +36,11 @@ test_that("LM for linear regression", { expect_equal(mod$getConstants()$x2, c(1,2,3)) # With model matrix names - code <- quote(LM(y ~ x + x2, priorSpecs=setPriors(sd="dunif(0, 5)"), modMatNames=TRUE)) + code <- quote(LM(y ~ x + x2, priors=setPriors(sd="dunif(0, 5)"), modelMatrixNames=TRUE)) out <- LM$process(code, modelInfo, environment()) code <- nimbleCode({ - LM(y ~ x + x2, modMatNames = TRUE) + LM(y ~ x + x2, modelMatrixNames = TRUE) }) mod <- nimbleModel(code, constants=dat) @@ -234,7 +234,7 @@ test_that("Run lm example", { lm.D9 <- stats::lm(weight ~ group) code.D9 <- nimbleCode({ - LM(weight ~ group, modMatNames = TRUE) + LM(weight ~ group, modelMatrixNames = TRUE) }) constants=data.frame(weight=weight, group=group) mod.D9 <- nimbleModel(code.D9, constants=constants) @@ -264,7 +264,7 @@ test_that("Run lm example", { lm.D90 <- stats::lm(weight ~ group - 1) # omitting intercept code.D90 <- nimbleCode({ - LM(weight ~ group - 1, modMatNames = TRUE) + LM(weight ~ group - 1, modelMatrixNames = TRUE) }) mod.D90 <- nimbleModel(code.D90, constants=constants) expect_equal( @@ -309,7 +309,7 @@ test_that("Run glm example", { pr <- setPriors(intercept="dnorm(0, sd=10)", coefficient="dnorm(0, sd = 10)") code.D93 <- nimbleCode({ LM(counts ~ outcome + treatment, family = poisson(), - modMatNames=TRUE, priorSpecs=pr) + modelMatrixNames=TRUE, priors=pr) }) mod.D93 <- nimbleModel(code.D93, constants=constants) @@ -349,7 +349,6 @@ test_that("Run glm example", { }) test_that("Run lme4::lmer() example", { - skip_on_ci() skip_on_cran() nimbleOptions(enableMacroComments = FALSE) sleepstudy <- readRDS('sleepstudy.Rds') # from lme4 package @@ -379,7 +378,7 @@ test_that("Run lme4::lmer() example", { sd_Days_Subject ~ dunif(0, 100) re_sds_Subject[1] <- sd_Subject re_sds_Subject[2] <- sd_Days_Subject - Ustar_Subject[1:2, 1:2] ~ dlkj_corr_cholesky(1.3, 2) + Ustar_Subject[1:2, 1:2] ~ dlkj_corr_cholesky(1, 2) U_Subject[1:2, 1:2] <- uppertri_mult_diag(Ustar_Subject[1:2, 1:2], re_sds_Subject[1:2]) re_means_Subject[1] <- 0 @@ -409,7 +408,7 @@ test_that("Run lme4::lmer() example", { `Days|Subject SD` = 6.766, Cor = -0.255, sigma = 25.526) comp <- cbind(lmer=lmer_est, nimble=nim_est) # comparison - expect_equivalent(nim_est, c(244.733, 11.268, 33.393, 7.366, -0.165, 25.966)) + expect_equivalent(nim_est, c(247.552, 11.408, 33.200, 7.267, -0.214, 25.807)) nimbleOptions(enableMacroComments = TRUE) }) @@ -428,8 +427,8 @@ test_that("Run lme4::glmer() example", { constants$noncase <- constants$size - constants$incidence pr <- setPriors(intercept="dnorm(0, sd=10)", coefficient="dnorm(0, sd=10)") code_gm1 <- nimbleCode({ - LM(cbind(incidence, noncase) ~ period + (1 | herd), priorSpecs=pr, - family=binomial, modMatNames=TRUE) + LM(cbind(incidence, noncase) ~ period + (1 | herd), priors=pr, + family=binomial, modelMatrixNames=TRUE) }) mod_gm1 <- nimbleModel(code_gm1, constants=constants) diff --git a/tests/testthat/test_formulaHandlers.R b/tests/testthat/test_formulaHandlers.R index 7561d00..fa0151f 100644 --- a/tests/testthat/test_formulaHandlers.R +++ b/tests/testthat/test_formulaHandlers.R @@ -8,13 +8,13 @@ test_that("Error when function in formula is unsupported", { x2=factor(sample(letters[4:5], 10, replace=T)), x3=round(rnorm(10),3))) - code <- quote(y[1:n] <- LINPRED(~test(x), priorSpecs=NULL)) + code <- quote(y[1:n] <- LINPRED(~test(x), priors=NULL)) expect_error(LINPRED$process(code, modelInfo=modInfo, .env=environment()), "No formula handler") - code <- quote(y[1:n] <- LINPRED(~test(x) + (1|x2), priorSpecs=NULL)) + code <- quote(y[1:n] <- LINPRED(~test(x) + (1|x2), priors=NULL)) expect_error(LINPRED$process(code, modelInfo=modInfo, .env=environment()), "No formula handler") - code <- quote(y[1:n] <- LINPRED(~x3 + test(x[1:10]), priorSpecs=NULL)) + code <- quote(y[1:n] <- LINPRED(~x3 + test(x[1:10]), priors=NULL)) expect_error(LINPRED$process(code, modelInfo=modInfo, .env=environment()), "No formula handler") }) @@ -621,3 +621,55 @@ test_that("log formula handler works", { nimbleOptions(enableMacroComments=TRUE) }) + +test_that("Can overwrite existing formulaHandler", { + + # Create new version of scale() that just returns vector of 1s + formulaHandler_scale <- function(x, defaultBracket, coefPrefix, + sdPrefix, modelInfo, env, ...){ + + const <- modelInfo$constants[[deparse(x$lang[[2]])]] + const <- rep(1, length(const)) + add_const <- list(x_test = const) + # Add the new constant to the object + if(is.null(x$constants)) x$constants <- list() + x$constants <- utils::modifyList(x$constants, add_const) + x$lang <- quote(x_test) + # Update class + class(x)[1] <- "formulaComponentFixed" + + x + } + class(formulaHandler_scale) <- c(class(formulaHandler_scale), "nimbleFormulaHandler") + + + nimbleOptions(enableMacroComments=FALSE) + set.seed(123) + constants <- list(y=rnorm(3), x = runif(3, 5, 10), z = rnorm(3), n=3, + x2=matrix(runif(3,5,10), 3, 1), z2=matrix(rnorm(3), 3, 1), + x3 = matrix(runif(9,5,10), 3, 3)) + + # Make sure the new scale is being used instead of the default one + code <- nimbleCode({ + mu[1:n] <- LINPRED(~scale(x) + z) + }) + mod <- nimbleModel(code, constants=constants) + + expect_equal( + mod$getCode(), + quote({ + for (i_1 in 1:n) { + mu[i_1] <- beta_Intercept + beta_x_test * x_test[i_1] + + beta_z * z[i_1] + } + beta_Intercept ~ dnorm(0, sd = 1000) + beta_x_test ~ dnorm(0, sd = 1000) + beta_z ~ dnorm(0, sd = 1000) + }) + ) + expect_equivalent( + mod$getConstants()$x_test, + rep(1, length(constants$x)) + ) + +}) diff --git a/tests/testthat/test_priors.R b/tests/testthat/test_priors.R index 50ae8d6..bf352d6 100644 --- a/tests/testthat/test_priors.R +++ b/tests/testthat/test_priors.R @@ -8,7 +8,7 @@ test_that("setPriors",{ expect_equal(setPriors(), list(intercept=quote(dnorm(0, sd = 1000)), coefficient=quote(dnorm(0, sd = 1000)), - sd=quote(dunif(0, 100)), eta=1.3) + sd=quote(dunif(0, 100)), lkjShape=1) ) expect_equal( @@ -16,7 +16,7 @@ test_that("setPriors",{ list(intercept=quote(dnorm(0, sd = 1000)), coefficient=quote(dnorm(0, sd = 1000)), sd=quote(dunif(0, 100)), - factor=quote(dnorm(0, 1)), eta=1.3) + factor=quote(dnorm(0, 1)), lkjShape=1) ) expect_equal( @@ -24,7 +24,7 @@ test_that("setPriors",{ list(intercept=quote(dnorm(0, sd = 1000)), coefficient=quote(dnorm(0, sd = 1000)), sd=quote(dunif(0, 100)), - continuous=quote(dnorm(0, 1)), eta=1.3) + continuous=quote(dnorm(0, 1)), lkjShape=1) ) expect_equal( @@ -32,20 +32,20 @@ test_that("setPriors",{ list(intercept=quote(dnorm(0, sd = 1000)), coefficient=quote(dnorm(0, sd = 1000)), sd=quote(dunif(0, 100)), - "alpha[1]"=quote(dnorm(0, 1)), eta=1.3) + "alpha[1]"=quote(dnorm(0, 1)), lkjShape=1) ) expect_equal( setPriors("alpha[1]" = "dnorm(0, 1)"), list(intercept=quote(dnorm(0, sd = 1000)), coefficient=quote(dnorm(0, sd = 1000)), sd=quote(dunif(0, 100)), - "alpha[1]"=quote(dnorm(0, 1)), eta=1.3) + "alpha[1]"=quote(dnorm(0, 1)), lkjShape=1) ) expect_equal( setPriors(sd = list("dnorm", 0, sd = 3)), list(intercept=quote(dnorm(0, sd = 1000)), coefficient=quote(dnorm(0, sd = 1000)), - sd=quote(dnorm(0, sd = 3)), eta=1.3) + sd=quote(dnorm(0, sd = 3)), lkjShape=1) ) # Prior not quoted should error @@ -89,32 +89,32 @@ test_that("matchPrior", { alpha = quote(dnorm(0, 5)), "alpha[1]" = quote(dnorm(0, 6))) - expect_equal(matchPrior(quote(beta), "intercept", priorSpecs=priors), + expect_equal(matchPrior(quote(beta), "intercept", priors=priors), quote(dnorm(0,1))) - expect_equal(matchPrior(quote(beta), "coefficient", priorSpecs=priors), + expect_equal(matchPrior(quote(beta), "coefficient", priors=priors), quote(dnorm(0,2))) - expect_equal(matchPrior(quote(beta), "continuous", "coefficient", priorSpecs=priors), + expect_equal(matchPrior(quote(beta), "continuous", "coefficient", priors=priors), quote(dnorm(0,3))) - expect_equal(matchPrior(quote(beta), "factor", "coefficient", priorSpecs=priors), + expect_equal(matchPrior(quote(beta), "factor", "coefficient", priors=priors), quote(dnorm(0,4))) - expect_equal(matchPrior(quote(alpha), "coefficient", priorSpecs=priors), + expect_equal(matchPrior(quote(alpha), "coefficient", priors=priors), quote(dnorm(0,5))) - expect_equal(matchPrior(quote(alpha[2]), "coefficient", priorSpecs=priors), + expect_equal(matchPrior(quote(alpha[2]), "coefficient", priors=priors), quote(dnorm(0,5))) - expect_equal(matchPrior(quote(alpha[1]), "coefficient", priorSpecs=priors), + expect_equal(matchPrior(quote(alpha[1]), "coefficient", priors=priors), quote(dnorm(0,6))) # Possible errors - expect_error(matchPrior(quote(beta), "fake", priorSpecs=priors)) + expect_error(matchPrior(quote(beta), "fake", priors=priors)) priors$bad <- "dnorm(0, 1)" - expect_error(matchPrior(quote(beta), "bad", priorSpecs=priors)) + expect_error(matchPrior(quote(beta), "bad", priors=priors)) }) test_that("spaces in factor levels are handled", { @@ -125,7 +125,7 @@ test_that("spaces in factor levels are handled", { dat <- list(y=y, x=x, z=z, n=3) code <- nimbleCode({ - LINPRED_PRIORS(~x + z, modMatNames=TRUE) + LINPRED_PRIORS(~x + z, modelMatrixNames=TRUE) }) nimbleOptions(enableMacroComments=FALSE) diff --git a/vignettes/nimbleMacros.Rmd b/vignettes/nimbleMacros.Rmd index 89d0ad8..1c1e116 100644 --- a/vignettes/nimbleMacros.Rmd +++ b/vignettes/nimbleMacros.Rmd @@ -49,7 +49,7 @@ We'll model transmission type (`am`, where `am = 1` is manual transmission) as a ``` r -rmod <- glm(am ~ scale(mpg), data = mtcars, family=binomial) +rmod <- glm(am ~ scale(mpg), data = mtcars, family = binomial) summary(rmod) ``` @@ -74,7 +74,7 @@ summary(rmod) ## Number of Fisher Scoring iterations: 5 ``` -To fit this model with the `LM` macro, we'll start by organizing the data and constants into a list. +To create this model with the `LM` macro, we'll start by organizing the data and constants into a list. ``` r @@ -88,7 +88,7 @@ Note that these distribution definitions are wrapped in `quote`; you can also su ``` r -pr <- setPriors(intercept=quote(dunif(-10,10)), +pr <- setPriors(intercept = quote(dunif(-10,10)), coefficient = quote(dnorm(0, 0.1))) ``` @@ -96,12 +96,12 @@ Finally, we specify the NIMBLE code, which is a single call to the `LM` macro th Two things to note: 1. We include a call to `scale` in the formula as with `glm()`; `LM` will do the scaling automatically. -2. You don't need to specify dimensions of the variables with `LM`, as this is handled automatically, as `LM` assumes all data and constants are vectors of the same length. In more complex situations you can use the `LINPRED` macro; see the next section. +2. You don't need to specify dimensions of the variables with `LM`, as this is handled automatically, with `LM` assuming all data and constants are vectors of the same length. In more complex situations you can use the `LINPRED` macro; see the next section. ``` r code <- nimbleCode({ - LM(am ~ scale(mpg), family = binomial, priorSpecs=pr) + LM(am ~ scale(mpg), family = binomial, priors = pr) }) ``` @@ -137,11 +137,11 @@ mod$getCode() ## } ``` -Finally, we'll fit the model which yields similar results to `glm`. +Finally, we'll fit the model using MCMC, which yields similar results to `glm`. ``` r -samples <- nimbleMCMC(mod, nchain=1, niter=3000, nburnin=2000, +samples <- nimbleMCMC(mod, nchain = 1, niter = 3000, nburnin = 2000, samplesAsCodaMCMC = TRUE) ``` @@ -178,7 +178,7 @@ summary(samples) ## Example linear mixed model using `ChickWeights` The `LM` macro can also generate code for random effects models. -To illustrate this we will use the built-in `ChickWeight` dataset which contains repeated measurements of chick weights over time. +To illustrate this we will use the built-in `ChickWeight` dataset, which contains repeated measurements of chick weights over time. First, fit the model in R using the `lme4` package: ```r @@ -214,14 +214,14 @@ summary(rmod) ## Time -0.422 ``` -To fit this model with `LM`, we'll first organize the `ChickWeight` data into a list of constants and data for NIMBLE. +To fit this model with NIMBLE using `LM`, we'll first organize the `ChickWeight` data into a list of constants and data for NIMBLE. ``` r chick_const <- list(weight = ChickWeight$weight, Time = ChickWeight$Time, Chick = ChickWeight$Chick) ``` -We can fit this model with `LM` using the same formula notation as `lme4` with the `(1|Chick)` part of the formula indicating random intercepts by chick. +We can create the model with `LM` using the same formula notation as `lme4` with the `(1|Chick)` part of the formula indicating random intercepts by chick. ``` r @@ -273,7 +273,7 @@ Fit the model: ``` r -samples <- nimbleMCMC(mod, nchain=1, niter=3000, nburnin=2000, +samples <- nimbleMCMC(mod, nchain = 1, niter = 3000, nburnin = 2000, samplesAsCodaMCMC = TRUE) ``` @@ -366,13 +366,13 @@ chick_const$N <- nrow(ChickWeight) ``` In the model code, we specify the parameter name for the linear predictor (`mu`) and calculate it with the `LINPRED` macro and the same formula as above. -Here we also set an argument `priorSpecs = NULL`; we'll come back to this later. +Here we also set an argument `priors = NULL`; we'll come back to this later. Notice that unlike with `LM`, we explicitly define the dimensions of `mu` and `Time`. ``` r code <- nimbleCode({ - mu[1:N] <- LINPRED(~Time[1:N], priorSpecs=NULL) + mu[1:N] <- LINPRED(~Time[1:N], priors = NULL) }) ``` @@ -423,7 +423,7 @@ When all the data and constants used by `LINPRED` have the same dimensions, you ``` r code <- nimbleCode({ - mu[1:N] <- LINPRED(~Time, priorSpecs=NULL) + mu[1:N] <- LINPRED(~Time, priors = NULL) }) mod <- nimbleModel(code, chick_const) mod$getCode() @@ -447,7 +447,7 @@ For example, we'll convert `Time` into a one-column matrix, but keep `mu` as a v ``` r chick_const2 <- chick_const -chick_const2$Time <- matrix(chick_const2$Time, ncol=1) +chick_const2$Time <- matrix(chick_const2$Time, ncol = 1) str(chick_const2) ``` @@ -464,7 +464,7 @@ To handle this we just need to specify different dimensions for each in the call ``` r code <- nimbleCode({ - mu[1:N] <- LINPRED(~Time[1:N,1], priorSpecs=NULL) + mu[1:N] <- LINPRED(~Time[1:N,1], priors = NULL) }) mod <- nimbleModel(code, chick_const2) mod$getCode() @@ -485,6 +485,7 @@ mod$getCode() ## Factor/categorical covariates The `LINPRED` macro can automatically handle categorical covariates (what R calls factors). +A categorical covariate can be provided in the constants as either an R factor, or as a character vector/matrix/array (which will be converted automatically to a factor). For example, suppose we included a covariate for diet type (`Diet`) in the model. The `Diet` covariate is already coded as an R factor: @@ -512,7 +513,7 @@ chick_const$Diet <- ChickWeight$Diet ``` r code <- nimbleCode({ - mu[1:N] <- LINPRED(~Time[1:N] + Diet[1:N], priorSpecs=NULL) + mu[1:N] <- LINPRED(~Time[1:N] + Diet[1:N], priors = NULL) }) mod <- nimbleModel(code, chick_const) mod$getCode() @@ -544,7 +545,7 @@ Suppose we wanted to force `mu` to be positive; we could use a log link function ``` r code <- nimbleCode({ - mu[1:N] <- LINPRED(~Time[1:N], priorSpecs = NULL, link = log) + mu[1:N] <- LINPRED(~Time[1:N], priors = NULL, link = log) }) mod <- nimbleModel(code, chick_const) mod$getCode() @@ -571,7 +572,7 @@ For example, to build a linear predictor with the `Time` covariate scaled to hav ``` r code <- nimbleCode({ - mu[1:N] <- LINPRED(~scale(Time[1:N]), priorSpecs = NULL, link = log) + mu[1:N] <- LINPRED(~scale(Time[1:N]), priors = NULL, link = log) }) mod <- nimbleModel(code, chick_const) mod$getCode() @@ -607,7 +608,7 @@ See the source code of the package, specifically the file `formulaFunction.R`, f ## Get default priors -If we do not manually set the `priorSpecs` argument, `LINPRED` will automatically generate a default set of priors for the two parameters it created (`beta_Intercept` and `beta_Time`). +If we do not manually set the `priors` argument, `LINPRED` will automatically generate a default set of priors for the two parameters it created (`beta_Intercept` and `beta_Time`). Internally, `LINPRED` is calling the `LINPRED_PRIORS` macro (described below). @@ -642,6 +643,9 @@ mod$getCode() We now have priors for each parameter. Note that as mentioned previously, the first value of the `beta_Diet` vector is fixed at 0, serving as the reference level. +Warning: one should check the default priors carefully to make sure that they make sense for your problem. +In particular, if the values in your data are of magnitudes such that parameters of the linear predictor could be large in magnitude, the default priors could have a strong influence on your results. + ## Suppress the macro comments We can also suppress the added macro comments to make things a little more concise. @@ -722,7 +726,7 @@ mod$getCode() ## sd_Time_Chick ~ dunif(0, 100) ## re_sds_Chick[1] <- sd_Chick ## re_sds_Chick[2] <- sd_Time_Chick -## Ustar_Chick[1:2, 1:2] ~ dlkj_corr_cholesky(1.3, 2) +## Ustar_Chick[1:2, 1:2] ~ dlkj_corr_cholesky(1, 2) ## U_Chick[1:2, 1:2] <- uppertri_mult_diag(Ustar_Chick[1:2, ## 1:2], re_sds_Chick[1:2]) ## re_means_Chick[1] <- 0 @@ -736,7 +740,7 @@ mod$getCode() ## } ``` -The correlated random effects model here uses an LKJ distribution prior; read more about this prior [here](https://r-nimble.org/html_manual/cha-writing-models.html#lkj-distribution-for-correlation-matrices). +The correlated random effects model here uses an [LKJ distribution prior](https://r-nimble.org/html_manual/cha-writing-models.html#lkj-distribution-for-correlation-matrices). As with `lme4` we can also specify uncorrelated random slopes and intercepts using `||` instead of `|`. @@ -802,6 +806,9 @@ mod$getCode() # `LINPRED_PRIORS` This macro generates a set of priors for parameters of a linear predictor, also based on the R formula. +In most cases, we won't use this macro directly; rather it will be called automatically when we use `LINPRED` (see discussion of `priors` argument above). +However, if for some reason we want to generate code for only the priors, and not the corresponding linear predictor, we can use `LINPRED_PRIORS` separately as we show below. + The formula `~Time` implies two parameters as we described above, and `LINPRED_PRIORS` will generate a prior for each depending on the type of the parameter (e.g. intercept vs. slope). It is not necessary to specify the dimensions in the formula when using `LINPRED_PRIORS` directly. @@ -830,7 +837,7 @@ pr <- setPriors(intercept = quote(dnorm(-5, 5)), coefficient = quote(dnorm(0, sd = 5))) code <- nimbleCode({ - LINPRED_PRIORS(~Time + Diet, priorSpecs = pr) + LINPRED_PRIORS(~Time + Diet, priors = pr) }) mod <- nimbleModel(code, chick_const) mod$getCode() @@ -859,12 +866,12 @@ colnames(model.matrix(~Time + Diet, ChickWeight)) ``` Here `model.matrix` has created several separate parameters `Diet2`, `Diet3`, and `Diet4` rather than a vector. -We can replicate this structure in NIMBLE model code by setting the `LINPRED_PRIORS` option `modMatNames = TRUE`. +We can replicate this structure in NIMBLE model code by setting the `LINPRED_PRIORS` option `modelMatrixNames = TRUE`. ``` r code <- nimbleCode({ - LINPRED_PRIORS(~Time + Diet, priorSpecs = pr, modMatNames = TRUE) + LINPRED_PRIORS(~Time + Diet, priors = pr, modelMatrixNames = TRUE) }) mod <- nimbleModel(code, chick_const) mod$getCode() @@ -884,7 +891,7 @@ mod$getCode() ## } ``` -With this option set a few extra lines of NIMBLE model code are added to yield parameter names that match `model.matrix`. +When this option is set, a few extra lines of NIMBLE model code are added to yield parameter names that match `model.matrix`. # `FORLOOP` @@ -927,6 +934,15 @@ mod$getCode() ## } ``` +Note that one can use vectorization in NIMBLE model code directly, such as + +``` +mu[1:N] <- alpha + beta * X[1:N] +``` + +The distinction is that that creates a single multivariate node, `mu[1:N]`, while the loop creates `N` separate scalar nodes. +Which is best will depend on the model structure and the algorithm used with the model [as discussed in the NIMBLE user manual](https://r-nimble.org/html_manual/cha-writing-models.html#subsec:vectorized-versus-scalar-declarations). + It may seem a bit trivial to make a macro for such a simple chunk of code. However it can be nice to use this when you have nested for loops, which are a bit more complicated. @@ -1025,15 +1041,16 @@ For example: ``` r fun <- function(input){ - LHS <- input[[2]] + lhs <- input[[2]] sdev <- input[[3]]$sdev output <- substitute(LHS ~ dnorm(0, sd = SDEV), - list(LHS = LHS, SDEV = sdev)) + list(LHS = lhs, SDEV = sdev)) output } ``` Note that we have to do a bit of clunky code-processing to extract the left-hand-side of the assignment and the argument value. +(The use of the `[[` extraction with a code object extracts information from the abstract syntax tree representation of an R expression.) ``` r @@ -1045,12 +1062,12 @@ fun(input) ## y ~ dnorm(0, sd = 2) ``` -In practice, `nimble` expects this function to have two additional arguments: `modelInfo` and `.env`. +In practice, `nimble` expects the macro function converter (`fun` in the example above) to have two additional arguments: `modelInfo` and `.env`. When processing a macro, `nimble` will provide a named list to the `modelInfo` argument which contains some additional model information. Most importantly, this list includes `constants`, which can be used and/or changed by the macro. -As you can probably guess, `nimble` provides the R environment in which `nimbleModel` was called to the `.env` argument, to help with scoping issues. +And, as you can probably guess, `nimble` provides the R environment in which `nimbleModel` was called to the `.env` argument, to help with scoping issues. -In addition to these two arguments, the function also needs a bit more structure in the output. +In addition to these two arguments, the function also needs to return a particular structure in the output. The output should be a named list with two elements: `code`, which contains the output code (as in the original function above), and `modelInfo` which is the list of model information, which we could modify or leave as-is. @@ -1071,7 +1088,7 @@ fun <- function(input_code, modelInfo, .env){ ``` r input <- quote(y ~ MYDIST(sdev = 2)) -fun(input, modelInfo = list(constants=NULL), .env = environment()) +fun(input, modelInfo = list(constants = NULL), .env = environment()) ``` ``` @@ -1086,6 +1103,7 @@ fun(input, modelInfo = list(constants=NULL), .env = environment()) Finally, `nimble` needs a way to know that this function is intended to be a macro, and that it is called `MYDIST`. The final step then is to create a wrapper `list` with our function inserted as an element called `process`, and assign the class `"model_macro"`. +(In the next section we'll see that the helper function `buildMacro` will create this wrapper list for us.) ``` r @@ -1112,8 +1130,8 @@ mod$getCode() ## Using `buildMacro` to create the macro -Manually creating the macro is fine, but gets complicated if we want to extract a lot of separate pieces of information from the code line, e.g. if our macro has many arguments. -As an alternative, we can use the `buildMacro` function included with `nimble` which handles some of this code processing for us. +Manually creating the macro is fine, but it gets complicated if we want to extract a lot of separate pieces of information from the code line, e.g., if our macro has many arguments. +As an alternative, we can use the `buildMacro` function included with `nimble`, which handles some of this code processing for us. The first argument to `buildMacro`, is a function structured like the one we wrote above `fun`. The other two arguments `use3pieces` and `unpackArgs` determine how the arguments of `fun` should be structured. @@ -1124,7 +1142,7 @@ In this case, all `buildMacro` will do is create the required macro list structu ``` r -MYDIST2 <- buildMacro(fun, use3pieces=FALSE, unpackArgs = FALSE) +MYDIST2 <- buildMacro(fun, use3pieces = FALSE, unpackArgs = FALSE) ``` @@ -1149,10 +1167,10 @@ With these settings, our function's first three arguments need to be `stoch`, `L ``` r -fun3 <- function(stoch, LHS, RHS, modelInfo, .env){ - sdev <- RHS$sdev +fun3 <- function(stoch, lhs, rhs, modelInfo, .env){ + sdev <- rhs$sdev output_code <- substitute(LHS ~ dnorm(0, sd = SDEV), - list(LHS = LHS, SDEV = sdev)) + list(LHS = lhs, SDEV = sdev)) # Modify the constants modelInfo$constants$new <- 1 @@ -1166,7 +1184,7 @@ Note that in this simple example we aren't using `stoch` because we want our out ``` r -MYDIST3 <- buildMacro(fun3, use3pieces=TRUE, unpackArgs = FALSE) +MYDIST3 <- buildMacro(fun3, use3pieces = TRUE, unpackArgs = FALSE) ``` @@ -1175,6 +1193,18 @@ code <- nimbleCode({ y ~ MYDIST3(sdev = 2) }) mod <- nimbleModel(code = code, data = list(y = 1)) +``` + +``` +## Error in (function (stoch, lhs, rhs, modelInfo, .env) : +## unused arguments (LHS = quote(y), RHS = quote(MYDIST3(sdev = 2))) +``` + +``` +## Error: Model macro MYDIST3 failed. +``` + +``` r mod$getCode() ``` @@ -1190,9 +1220,9 @@ Thus, we need to include all these arguments explicitly as arguments in our func ``` r -fun4 <- function(stoch, LHS, sdev, modelInfo, .env){ +fun4 <- function(stoch, lhs, sdev, modelInfo, .env){ output_code <- substitute(LHS ~ dnorm(0, sd = SDEV), - list(LHS = LHS, SDEV = sdev)) + list(LHS = lhs, SDEV = sdev)) # Modify the constants modelInfo$constants$new <- 1 @@ -1205,7 +1235,7 @@ Even less code processing is required in this case, since `buildMacro` will spli ``` r -MYDIST4 <- buildMacro(fun4, use3pieces=TRUE, unpackArgs = TRUE) +MYDIST4 <- buildMacro(fun4, use3pieces = TRUE, unpackArgs = TRUE) ``` @@ -1214,6 +1244,18 @@ code <- nimbleCode({ y ~ MYDIST4(sdev = 2) }) mod <- nimbleModel(code = code, data = list(y = 1)) +``` + +``` +## Error in (function (stoch, lhs, sdev, modelInfo, .env) : +## unused argument (LHS = quote(y)) +``` + +``` +## Error: Model macro MYDIST4 failed. +``` + +``` r mod$getCode() ``` diff --git a/vignettes/nimbleMacros.Rmd.orig b/vignettes/nimbleMacros.Rmd.orig index 81acf09..e3e6be8 100644 --- a/vignettes/nimbleMacros.Rmd.orig +++ b/vignettes/nimbleMacros.Rmd.orig @@ -51,11 +51,11 @@ Here's a logistic regression in R using the built-in `mtcars` dataset. We'll model transmission type (`am`, where `am = 1` is manual transmission) as a function of standardized miles per gallon. ```{r} -rmod <- glm(am ~ scale(mpg), data = mtcars, family=binomial) +rmod <- glm(am ~ scale(mpg), data = mtcars, family = binomial) summary(rmod) ``` -To fit this model with the `LM` macro, we'll start by organizing the data and constants into a list. +To create this model with the `LM` macro, we'll start by organizing the data and constants into a list. ```{r} const <- mtcars[c("am", "mpg")] @@ -67,7 +67,7 @@ We'll specify a uniform(-10, 10) distribution for the intercept and a normal dis Note that these distribution definitions are wrapped in `quote`; you can also supply them as strings. ```{r} -pr <- setPriors(intercept=quote(dunif(-10,10)), +pr <- setPriors(intercept = quote(dunif(-10,10)), coefficient = quote(dnorm(0, 0.1))) ``` @@ -75,11 +75,11 @@ Finally, we specify the NIMBLE code, which is a single call to the `LM` macro th Two things to note: 1. We include a call to `scale` in the formula as with `glm()`; `LM` will do the scaling automatically. -2. You don't need to specify dimensions of the variables with `LM`, as this is handled automatically, as `LM` assumes all data and constants are vectors of the same length. In more complex situations you can use the `LINPRED` macro; see the next section. +2. You don't need to specify dimensions of the variables with `LM`, as this is handled automatically, with `LM` assuming all data and constants are vectors of the same length. In more complex situations you can use the `LINPRED` macro; see the next section. ```{r} code <- nimbleCode({ - LM(am ~ scale(mpg), family = binomial, priorSpecs=pr) + LM(am ~ scale(mpg), family = binomial, priors = pr) }) ``` @@ -90,10 +90,10 @@ mod <- nimbleModel(code = code, constants = const) mod$getCode() ``` -Finally, we'll fit the model which yields similar results to `glm`. +Finally, we'll fit the model using MCMC, which yields similar results to `glm`. ```{r} -samples <- nimbleMCMC(mod, nchain=1, niter=3000, nburnin=2000, +samples <- nimbleMCMC(mod, nchain = 1, niter = 3000, nburnin = 2000, samplesAsCodaMCMC = TRUE) summary(samples) ``` @@ -101,7 +101,7 @@ summary(samples) ## Example linear mixed model using `ChickWeights` The `LM` macro can also generate code for random effects models. -To illustrate this we will use the built-in `ChickWeight` dataset which contains repeated measurements of chick weights over time. +To illustrate this we will use the built-in `ChickWeight` dataset, which contains repeated measurements of chick weights over time. First, fit the model in R using the `lme4` package: ```r @@ -137,13 +137,13 @@ summary(rmod) ## Time -0.422 ``` -To fit this model with `LM`, we'll first organize the `ChickWeight` data into a list of constants and data for NIMBLE. +To fit this model with NIMBLE using `LM`, we'll first organize the `ChickWeight` data into a list of constants and data for NIMBLE. ```{r} chick_const <- list(weight = ChickWeight$weight, Time = ChickWeight$Time, Chick = ChickWeight$Chick) ``` -We can fit this model with `LM` using the same formula notation as `lme4` with the `(1|Chick)` part of the formula indicating random intercepts by chick. +We can create the model with `LM` using the same formula notation as `lme4` with the `(1|Chick)` part of the formula indicating random intercepts by chick. ```{r} code <- nimbleCode({ @@ -161,7 +161,7 @@ mod$getCode() Fit the model: ```{r} -samples <- nimbleMCMC(mod, nchain=1, niter=3000, nburnin=2000, +samples <- nimbleMCMC(mod, nchain = 1, niter = 3000, nburnin = 2000, samplesAsCodaMCMC = TRUE) summary(samples) ``` @@ -198,12 +198,12 @@ chick_const$N <- nrow(ChickWeight) ``` In the model code, we specify the parameter name for the linear predictor (`mu`) and calculate it with the `LINPRED` macro and the same formula as above. -Here we also set an argument `priorSpecs = NULL`; we'll come back to this later. +Here we also set an argument `priors = NULL`; we'll come back to this later. Notice that unlike with `LM`, we explicitly define the dimensions of `mu` and `Time`. ```{r} code <- nimbleCode({ - mu[1:N] <- LINPRED(~Time[1:N], priorSpecs=NULL) + mu[1:N] <- LINPRED(~Time[1:N], priors = NULL) }) ``` @@ -228,7 +228,7 @@ When all the data and constants used by `LINPRED` have the same dimensions, you ```{r} code <- nimbleCode({ - mu[1:N] <- LINPRED(~Time, priorSpecs=NULL) + mu[1:N] <- LINPRED(~Time, priors = NULL) }) mod <- nimbleModel(code, chick_const) mod$getCode() @@ -239,7 +239,7 @@ For example, we'll convert `Time` into a one-column matrix, but keep `mu` as a v ```{r} chick_const2 <- chick_const -chick_const2$Time <- matrix(chick_const2$Time, ncol=1) +chick_const2$Time <- matrix(chick_const2$Time, ncol = 1) str(chick_const2) ``` @@ -247,7 +247,7 @@ To handle this we just need to specify different dimensions for each in the call ```{r} code <- nimbleCode({ - mu[1:N] <- LINPRED(~Time[1:N,1], priorSpecs=NULL) + mu[1:N] <- LINPRED(~Time[1:N,1], priors = NULL) }) mod <- nimbleModel(code, chick_const2) mod$getCode() @@ -256,6 +256,7 @@ mod$getCode() ## Factor/categorical covariates The `LINPRED` macro can automatically handle categorical covariates (what R calls factors). +A categorical covariate can be provided in the constants as either an R factor, or as a character vector/matrix/array (which will be converted automatically to a factor). For example, suppose we included a covariate for diet type (`Diet`) in the model. The `Diet` covariate is already coded as an R factor: @@ -267,7 +268,7 @@ chick_const$Diet <- ChickWeight$Diet ```{r} code <- nimbleCode({ - mu[1:N] <- LINPRED(~Time[1:N] + Diet[1:N], priorSpecs=NULL) + mu[1:N] <- LINPRED(~Time[1:N] + Diet[1:N], priors = NULL) }) mod <- nimbleModel(code, chick_const) mod$getCode() @@ -286,7 +287,7 @@ Suppose we wanted to force `mu` to be positive; we could use a log link function ```{r} code <- nimbleCode({ - mu[1:N] <- LINPRED(~Time[1:N], priorSpecs = NULL, link = log) + mu[1:N] <- LINPRED(~Time[1:N], priors = NULL, link = log) }) mod <- nimbleModel(code, chick_const) mod$getCode() @@ -300,7 +301,7 @@ For example, to build a linear predictor with the `Time` covariate scaled to hav ```{r} code <- nimbleCode({ - mu[1:N] <- LINPRED(~scale(Time[1:N]), priorSpecs = NULL, link = log) + mu[1:N] <- LINPRED(~scale(Time[1:N]), priors = NULL, link = log) }) mod <- nimbleModel(code, chick_const) mod$getCode() @@ -319,7 +320,7 @@ See the source code of the package, specifically the file `formulaFunction.R`, f ## Get default priors -If we do not manually set the `priorSpecs` argument, `LINPRED` will automatically generate a default set of priors for the two parameters it created (`beta_Intercept` and `beta_Time`). +If we do not manually set the `priors` argument, `LINPRED` will automatically generate a default set of priors for the two parameters it created (`beta_Intercept` and `beta_Time`). Internally, `LINPRED` is calling the `LINPRED_PRIORS` macro (described below). ```{r} @@ -333,6 +334,9 @@ mod$getCode() We now have priors for each parameter. Note that as mentioned previously, the first value of the `beta_Diet` vector is fixed at 0, serving as the reference level. +Warning: one should check the default priors carefully to make sure that they make sense for your problem. +In particular, if the values in your data are of magnitudes such that parameters of the linear predictor could be large in magnitude, the default priors could have a strong influence on your results. + ## Suppress the macro comments We can also suppress the added macro comments to make things a little more concise. @@ -373,7 +377,7 @@ mod <- nimbleModel(code, chick_const) mod$getCode() ``` -The correlated random effects model here uses an LKJ distribution prior; read more about this prior [here](https://r-nimble.org/html_manual/cha-writing-models.html#lkj-distribution-for-correlation-matrices). +The correlated random effects model here uses an [LKJ distribution prior](https://r-nimble.org/html_manual/cha-writing-models.html#lkj-distribution-for-correlation-matrices). As with `lme4` we can also specify uncorrelated random slopes and intercepts using `||` instead of `|`. @@ -404,6 +408,9 @@ mod$getCode() # `LINPRED_PRIORS` This macro generates a set of priors for parameters of a linear predictor, also based on the R formula. +In most cases, we won't use this macro directly; rather it will be called automatically when we use `LINPRED` (see discussion of `priors` argument above). +However, if for some reason we want to generate code for only the priors, and not the corresponding linear predictor, we can use `LINPRED_PRIORS` separately as we show below. + The formula `~Time` implies two parameters as we described above, and `LINPRED_PRIORS` will generate a prior for each depending on the type of the parameter (e.g. intercept vs. slope). It is not necessary to specify the dimensions in the formula when using `LINPRED_PRIORS` directly. @@ -423,7 +430,7 @@ pr <- setPriors(intercept = quote(dnorm(-5, 5)), coefficient = quote(dnorm(0, sd = 5))) code <- nimbleCode({ - LINPRED_PRIORS(~Time + Diet, priorSpecs = pr) + LINPRED_PRIORS(~Time + Diet, priors = pr) }) mod <- nimbleModel(code, chick_const) mod$getCode() @@ -436,17 +443,17 @@ colnames(model.matrix(~Time + Diet, ChickWeight)) ``` Here `model.matrix` has created several separate parameters `Diet2`, `Diet3`, and `Diet4` rather than a vector. -We can replicate this structure in NIMBLE model code by setting the `LINPRED_PRIORS` option `modMatNames = TRUE`. +We can replicate this structure in NIMBLE model code by setting the `LINPRED_PRIORS` option `modelMatrixNames = TRUE`. ```{r} code <- nimbleCode({ - LINPRED_PRIORS(~Time + Diet, priorSpecs = pr, modMatNames = TRUE) + LINPRED_PRIORS(~Time + Diet, priors = pr, modelMatrixNames = TRUE) }) mod <- nimbleModel(code, chick_const) mod$getCode() ``` -With this option set a few extra lines of NIMBLE model code are added to yield parameter names that match `model.matrix`. +When this option is set, a few extra lines of NIMBLE model code are added to yield parameter names that match `model.matrix`. # `FORLOOP` @@ -478,6 +485,15 @@ mod <- nimbleModel(code, constants) mod$getCode() ``` +Note that one can use vectorization in NIMBLE model code directly, such as + +``` +mu[1:N] <- alpha + beta * X[1:N] +``` + +The distinction is that that creates a single multivariate node, `mu[1:N]`, while the loop creates `N` separate scalar nodes. +Which is best will depend on the model structure and the algorithm used with the model [as discussed in the NIMBLE user manual](https://r-nimble.org/html_manual/cha-writing-models.html#subsec:vectorized-versus-scalar-declarations). + It may seem a bit trivial to make a macro for such a simple chunk of code. However it can be nice to use this when you have nested for loops, which are a bit more complicated. @@ -541,27 +557,28 @@ For example: ```{r} fun <- function(input){ - LHS <- input[[2]] + lhs <- input[[2]] sdev <- input[[3]]$sdev output <- substitute(LHS ~ dnorm(0, sd = SDEV), - list(LHS = LHS, SDEV = sdev)) + list(LHS = lhs, SDEV = sdev)) output } ``` Note that we have to do a bit of clunky code-processing to extract the left-hand-side of the assignment and the argument value. +(The use of the `[[` extraction with a code object extracts information from the abstract syntax tree representation of an R expression.) ```{r} input <- quote(y ~ MYDIST(sdev = 2)) fun(input) ``` -In practice, `nimble` expects this function to have two additional arguments: `modelInfo` and `.env`. +In practice, `nimble` expects the macro function converter (`fun` in the example above) to have two additional arguments: `modelInfo` and `.env`. When processing a macro, `nimble` will provide a named list to the `modelInfo` argument which contains some additional model information. Most importantly, this list includes `constants`, which can be used and/or changed by the macro. -As you can probably guess, `nimble` provides the R environment in which `nimbleModel` was called to the `.env` argument, to help with scoping issues. +And, as you can probably guess, `nimble` provides the R environment in which `nimbleModel` was called to the `.env` argument, to help with scoping issues. -In addition to these two arguments, the function also needs a bit more structure in the output. +In addition to these two arguments, the function also needs to return a particular structure in the output. The output should be a named list with two elements: `code`, which contains the output code (as in the original function above), and `modelInfo` which is the list of model information, which we could modify or leave as-is. ```{r} @@ -580,11 +597,12 @@ fun <- function(input_code, modelInfo, .env){ ```{r} input <- quote(y ~ MYDIST(sdev = 2)) -fun(input, modelInfo = list(constants=NULL), .env = environment()) +fun(input, modelInfo = list(constants = NULL), .env = environment()) ``` Finally, `nimble` needs a way to know that this function is intended to be a macro, and that it is called `MYDIST`. The final step then is to create a wrapper `list` with our function inserted as an element called `process`, and assign the class `"model_macro"`. +(In the next section we'll see that the helper function `buildMacro` will create this wrapper list for us.) ```{r} MYDIST <- list(process = fun) @@ -603,8 +621,8 @@ mod$getCode() ## Using `buildMacro` to create the macro -Manually creating the macro is fine, but gets complicated if we want to extract a lot of separate pieces of information from the code line, e.g. if our macro has many arguments. -As an alternative, we can use the `buildMacro` function included with `nimble` which handles some of this code processing for us. +Manually creating the macro is fine, but it gets complicated if we want to extract a lot of separate pieces of information from the code line, e.g., if our macro has many arguments. +As an alternative, we can use the `buildMacro` function included with `nimble`, which handles some of this code processing for us. The first argument to `buildMacro`, is a function structured like the one we wrote above `fun`. The other two arguments `use3pieces` and `unpackArgs` determine how the arguments of `fun` should be structured. @@ -614,7 +632,7 @@ This is identical to how we wrote the function above. In this case, all `buildMacro` will do is create the required macro list structure and assign the class for us. ```{r} -MYDIST2 <- buildMacro(fun, use3pieces=FALSE, unpackArgs = FALSE) +MYDIST2 <- buildMacro(fun, use3pieces = FALSE, unpackArgs = FALSE) ``` ```{r} @@ -631,10 +649,10 @@ Specifically, it will identify if the assignment is stochastic (`stoch`, `TRUE` With these settings, our function's first three arguments need to be `stoch`, `LHS`, and `RHS`, in that order (plus we keep the always-required `modelInfo` and `.env`). ```{r} -fun3 <- function(stoch, LHS, RHS, modelInfo, .env){ - sdev <- RHS$sdev +fun3 <- function(stoch, lhs, rhs, modelInfo, .env){ + sdev <- rhs$sdev output_code <- substitute(LHS ~ dnorm(0, sd = SDEV), - list(LHS = LHS, SDEV = sdev)) + list(LHS = lhs, SDEV = sdev)) # Modify the constants modelInfo$constants$new <- 1 @@ -647,7 +665,7 @@ We do a little less code processing this time since `buildMacro` has already sep Note that in this simple example we aren't using `stoch` because we want our output code to always be stochastic. ```{r} -MYDIST3 <- buildMacro(fun3, use3pieces=TRUE, unpackArgs = FALSE) +MYDIST3 <- buildMacro(fun3, use3pieces = TRUE, unpackArgs = FALSE) ``` ```{r} @@ -663,9 +681,9 @@ Now, `buildMacro` will dig into the right-hand-side of the macro and extract the Thus, we need to include all these arguments explicitly as arguments in our function. ```{r} -fun4 <- function(stoch, LHS, sdev, modelInfo, .env){ +fun4 <- function(stoch, lhs, sdev, modelInfo, .env){ output_code <- substitute(LHS ~ dnorm(0, sd = SDEV), - list(LHS = LHS, SDEV = sdev)) + list(LHS = lhs, SDEV = sdev)) # Modify the constants modelInfo$constants$new <- 1 @@ -677,7 +695,7 @@ fun4 <- function(stoch, LHS, sdev, modelInfo, .env){ Even less code processing is required in this case, since `buildMacro` will split out the `sdev` argument for us. ```{r} -MYDIST4 <- buildMacro(fun4, use3pieces=TRUE, unpackArgs = TRUE) +MYDIST4 <- buildMacro(fun4, use3pieces = TRUE, unpackArgs = TRUE) ``` ```{r}