From 9b9f098475307c480088bcd41584281b1672c741 Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Thu, 13 Feb 2025 17:25:39 -0800 Subject: [PATCH 01/14] Make various minor roxygen/package structure/vignette suggestions. --- DESCRIPTION | 6 +-- LICENSE | 3 ++ R/FORLOOP.R | 2 +- R/LINPRED.R | 59 +++++++++++++++-------------- R/LM.R | 33 +++++++++-------- R/setPriors.R | 33 +++++++++-------- README.md | 2 +- vignettes/nimbleMacros.Rmd | 76 +++++++++++++++++++++----------------- 8 files changed, 117 insertions(+), 97 deletions(-) create mode 100644 LICENSE diff --git a/DESCRIPTION b/DESCRIPTION index b375270..be38edf 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: nimbleMacros Version: 0.1.0 Date: 2025-01-28 -Title: Macros Generating 'nimble' Code for Linear Models +Title: Macros Generating 'nimble' Code for Linear Models. Authors@R: c(person("Ken", "Kellner", email="contact@kenkellner.com", role=c("cre","aut")), person("Perry", "de Valpine", role = "aut"), @@ -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/R/FORLOOP.R b/R/FORLOOP.R index 28d6a8f..9a438c9 100644 --- a/R/FORLOOP.R +++ b/R/FORLOOP.R @@ -17,7 +17,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..eff97e9 100644 --- a/R/LINPRED.R +++ b/R/LINPRED.R @@ -10,20 +10,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. 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 priorSpecs Prior specifications, generated with \code{setPrior()} +#' @param modMatNames 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 noncenter Logical indicating whether to use noncentered parameterization. +#' Default is \code{FALSE}. +#' @param centerVar Grouping variable (covariate) to 'center' the random effects on. By +#' default all random effects have mean 0 as with \code{lme4}. #' #' @author Ken Kellner #' @@ -33,7 +34,7 @@ #' mu[1:3] <- LINPRED(~x + x2) #' }) #' -#' mod <- nimbleModel(code, constants=constants) +#' mod <- nimbleModel(code, constants = constants) #' mod$getCode() NULL @@ -92,29 +93,31 @@ 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 priorSpecs 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 modMatNames 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 noncenter Logical indicating whether to use noncentered parameterization. +#' Default is \code{FALSE}. +#' @param centerVar Grouping variable (covariate) to 'center' the random effects on. By +#' default all random effects have mean 0 as with \code{lme4}. #' +#' #' @author Ken Kellner #' #' @examples @@ -1378,12 +1381,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..eebe27f 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,19 @@ #' @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 priorSpecs List of prior specifications, generated using +#' \code{setPriors()}. +#' @param modMatNames Logical indicating if parameters be named so they match the +#' names one would get from R's \code{model.matrix}. #' #' @author Ken Kellner #' @@ -35,7 +36,7 @@ #' LM(y ~ x + x2) #' }) #' -#' mod <- nimbleModel(code, constants=constants) +#' mod <- nimbleModel(code, constants = constants) #' mod$getCode() NULL @@ -81,7 +82,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]])))) @@ -145,10 +146,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 +168,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/setPriors.R b/R/setPriors.R index 49d68e1..6e51fe4 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 @@ -28,6 +21,16 @@ #' #' @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, + eta = 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 @@ -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 @@ -116,7 +119,7 @@ removeBracket <- function(node){ #' @param ... Character strings that categorize the parameter and match #' the names of elements in \code{priorSpecs}. The order is important: #' the first match found is used. -#' @param priorSpecs A named list of prior settings, e.g. as generated by +#' @param priorSpecs A named list of prior settings, e.g., as generated by #' \code{setPriors} #' #' @author Ken Kellner @@ -161,10 +164,10 @@ matchPrior <- function(parName, ..., priorSpecs){ return(priorSpecs[[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/vignettes/nimbleMacros.Rmd b/vignettes/nimbleMacros.Rmd index 89d0ad8..c8f19df 100644 --- a/vignettes/nimbleMacros.Rmd +++ b/vignettes/nimbleMacros.Rmd @@ -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 @@ -96,7 +96,7 @@ 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 @@ -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 based on 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) ``` @@ -372,7 +372,7 @@ Notice that unlike with `LM`, we explicitly define the dimensions of `mu` and `T ``` r code <- nimbleCode({ - mu[1:N] <- LINPRED(~Time[1:N], priorSpecs=NULL) + mu[1:N] <- LINPRED(~Time[1:N], priorSpecs = 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, priorSpecs = 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], priorSpecs = NULL) }) mod <- nimbleModel(code, chick_const2) mod$getCode() @@ -512,7 +512,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], priorSpecs = NULL) }) mod <- nimbleModel(code, chick_const) mod$getCode() @@ -642,6 +642,8 @@ 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. @@ -736,7 +738,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 `|`. @@ -884,7 +886,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 +929,14 @@ 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 sclar 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 +1035,15 @@ 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. +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 +1055,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. @@ -1085,7 +1095,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"`. +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 +1122,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 +1134,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 +1159,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 +1176,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) ``` @@ -1190,11 +1200,11 @@ 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 + # Modify the constants. modelInfo$constants$new <- 1 list(code = output_code, modelInfo = modelInfo) @@ -1205,7 +1215,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) ``` From 75d2791d96f573c0579397ce280ec7a7fe9b6732 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Fri, 14 Feb 2025 11:44:16 -0500 Subject: [PATCH 02/14] Update tests to reflect new eta value --- tests/testthat/test_LINPRED.R | 8 ++++---- tests/testthat/test_LM.R | 2 +- tests/testthat/test_priors.R | 12 ++++++------ 3 files changed, 11 insertions(+), 11 deletions(-) diff --git a/tests/testthat/test_LINPRED.R b/tests/testthat/test_LINPRED.R index cc3ebb7..9a38aa8 100644 --- a/tests/testthat/test_LINPRED.R +++ b/tests/testthat/test_LINPRED.R @@ -1031,7 +1031,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 @@ -1106,7 +1106,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 +1161,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 +1198,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 diff --git a/tests/testthat/test_LM.R b/tests/testthat/test_LM.R index 9d6eeae..343eeda 100644 --- a/tests/testthat/test_LM.R +++ b/tests/testthat/test_LM.R @@ -379,7 +379,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 diff --git a/tests/testthat/test_priors.R b/tests/testthat/test_priors.R index 50ae8d6..13b0308 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)), eta=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)), eta=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)), eta=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)), eta=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)), eta=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)), eta=1) ) # Prior not quoted should error From dd09cc585819483655066d984212e7402d69dcc3 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Fri, 14 Feb 2025 11:53:36 -0500 Subject: [PATCH 03/14] Apply roxygen changes and fix a couple typos --- DESCRIPTION | 2 +- man/FORLOOP.Rd | 2 +- man/LINPRED.Rd | 25 +++++++++++++------------ man/LINPRED_PRIORS.Rd | 25 +++++++++++++------------ man/LM.Rd | 25 +++++++++++++------------ man/matchPrior.Rd | 4 ++-- man/setPriors.Rd | 6 ++++-- man/uppertri_mult_diag.Rd | 8 ++++---- vignettes/nimbleMacros.Rmd | 2 +- 9 files changed, 52 insertions(+), 47 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index be38edf..63b9a41 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: nimbleMacros Version: 0.1.0 Date: 2025-01-28 -Title: Macros Generating 'nimble' Code for Linear Models. +Title: Macros Generating 'nimble' Code for Linear Models Authors@R: c(person("Ken", "Kellner", email="contact@kenkellner.com", role=c("cre","aut")), person("Perry", "de Valpine", role = "aut"), diff --git a/man/FORLOOP.Rd b/man/FORLOOP.Rd index 4fb2f1d..e24a18e 100644 --- a/man/FORLOOP.Rd +++ b/man/FORLOOP.Rd @@ -18,7 +18,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..0bbde71 100644 --- a/man/LINPRED.Rd +++ b/man/LINPRED.Rd @@ -7,27 +7,28 @@ \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{priorSpecs}{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{modMatNames}{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{noncenter}{Logical indicating whether to use noncentered parameterization. +Default is \code{FALSE}.} -\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 all random effects have mean 0 as with \code{lme4}.} } \description{ Converts an R formula into corresponding code for a linear predictor in NIMBLE model code. @@ -40,7 +41,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..4ce4f6d 100644 --- a/man/LINPRED_PRIORS.Rd +++ b/man/LINPRED_PRIORS.Rd @@ -2,32 +2,33 @@ % 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{priorSpecs}{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{modMatNames}{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{noncenter}{Logical indicating whether to use noncentered parameterization. +Default is \code{FALSE}.} -\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 all random effects have mean 0 as with \code{lme4}.} } \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..39ce30c 100644 --- a/man/LM.Rd +++ b/man/LM.Rd @@ -7,28 +7,29 @@ \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{priorSpecs}{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{modMatNames}{Logical indicating if parameters be named so they match the +names one would get from R's \code{model.matrix}.} } \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 +41,7 @@ code <- nimbleCode({ LM(y ~ x + x2) }) -mod <- nimbleModel(code, constants=constants) +mod <- nimbleModel(code, constants = constants) mod$getCode() } \author{ diff --git a/man/matchPrior.Rd b/man/matchPrior.Rd index 7e3d73b..ce0436e 100644 --- a/man/matchPrior.Rd +++ b/man/matchPrior.Rd @@ -14,12 +14,12 @@ including brackets/indices} the names of elements in \code{priorSpecs}. The order is important: the first match found is used.} -\item{priorSpecs}{A named list of prior settings, e.g. as generated by +\item{priorSpecs}{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 diff --git a/man/setPriors.Rd b/man/setPriors.Rd index e7670f3..0995bee 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, + eta = 1, ... ) } @@ -34,9 +34,11 @@ correlation matrix in correlated random slope and intercept models} } \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/vignettes/nimbleMacros.Rmd b/vignettes/nimbleMacros.Rmd index c8f19df..b16d9de 100644 --- a/vignettes/nimbleMacros.Rmd +++ b/vignettes/nimbleMacros.Rmd @@ -935,7 +935,7 @@ 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 sclar 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). +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. From 3c255a91003556df059c95792ebff6c116f19ee4 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Sat, 15 Feb 2025 15:03:17 -0500 Subject: [PATCH 04/14] Don't export formula handlers --- NAMESPACE | 4 - R/LINPRED.R | 12 +- R/{formulaFunction.R => formulaHandler.R} | 170 +++++++++++----------- man/formulaHandler_I.Rd | 37 ----- man/formulaHandler_log.Rd | 44 ------ man/formulaHandler_offset.Rd | 44 ------ man/formulaHandler_scale.Rd | 46 ------ tests/testthat/test_formulaHandlers.R | 52 +++++++ 8 files changed, 145 insertions(+), 264 deletions(-) rename R/{formulaFunction.R => formulaHandler.R} (64%) delete mode 100644 man/formulaHandler_I.Rd delete mode 100644 man/formulaHandler_log.Rd delete mode 100644 man/formulaHandler_offset.Rd delete mode 100644 man/formulaHandler_scale.Rd 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/LINPRED.R b/R/LINPRED.R index eff97e9..d1f399d 100644 --- a/R/LINPRED.R +++ b/R/LINPRED.R @@ -366,9 +366,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, ...) 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/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/tests/testthat/test_formulaHandlers.R b/tests/testthat/test_formulaHandlers.R index 7561d00..19fc1bd 100644 --- a/tests/testthat/test_formulaHandlers.R +++ b/tests/testthat/test_formulaHandlers.R @@ -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)) + ) + +}) From 6835e7297b929f53358f7f7aa7800c04f3e98f6a Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Sat, 15 Feb 2025 15:08:59 -0500 Subject: [PATCH 05/14] Fix test after updating default eta prior --- tests/testthat/test_LM.R | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/testthat/test_LM.R b/tests/testthat/test_LM.R index 343eeda..7614344 100644 --- a/tests/testthat/test_LM.R +++ b/tests/testthat/test_LM.R @@ -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 @@ -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) }) From 9d0cf9a877537ee11a42c9078e35a45d529881d4 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Sun, 16 Feb 2025 11:10:11 -0500 Subject: [PATCH 06/14] Error-trap formulas with a LHS in LINPRED --- R/LINPRED.R | 15 ++++++++++++--- tests/testthat/test_LINPRED.R | 16 ++++++++++++++++ 2 files changed, 28 insertions(+), 3 deletions(-) diff --git a/R/LINPRED.R b/R/LINPRED.R index d1f399d..dc846d6 100644 --- a/R/LINPRED.R +++ b/R/LINPRED.R @@ -45,7 +45,7 @@ function(stoch, LHS, formula, link=NULL, coefPrefix=quote(beta_), noncenter = 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) @@ -136,8 +136,7 @@ function(formula, coefPrefix=quote(beta_), sdPrefix=NULL, priorSpecs=setPriors() modMatNames=FALSE, noncenter = 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) @@ -171,6 +170,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) diff --git a/tests/testthat/test_LINPRED.R b/tests/testthat/test_LINPRED.R index 9a38aa8..d56f523 100644 --- a/tests/testthat/test_LINPRED.R +++ b/tests/testthat/test_LINPRED.R @@ -386,6 +386,22 @@ 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 + + # 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) From 368f7d20070e849eb039a9f83444e3117f9895a2 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Sun, 16 Feb 2025 16:47:42 -0500 Subject: [PATCH 07/14] Move some edits to the vignette source --- vignettes/nimbleMacros.Rmd | 46 +++++++++++++++++----- vignettes/nimbleMacros.Rmd.orig | 70 ++++++++++++++++++++------------- 2 files changed, 79 insertions(+), 37 deletions(-) diff --git a/vignettes/nimbleMacros.Rmd b/vignettes/nimbleMacros.Rmd index b16d9de..0c656e5 100644 --- a/vignettes/nimbleMacros.Rmd +++ b/vignettes/nimbleMacros.Rmd @@ -214,7 +214,7 @@ summary(rmod) ## Time -0.422 ``` -To fit this model with NIMBLE based on using `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 @@ -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, priorSpecs=NULL) }) mod <- nimbleModel(code, chick_const) mod$getCode() @@ -512,7 +512,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], priorSpecs=NULL) }) mod <- nimbleModel(code, chick_const) mod$getCode() @@ -642,7 +642,8 @@ 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. +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 @@ -724,7 +725,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 @@ -935,7 +936,8 @@ 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). +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. @@ -1043,7 +1045,8 @@ fun <- function(input){ } ``` -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.) +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 @@ -1095,7 +1098,8 @@ 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. +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 @@ -1185,6 +1189,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() ``` @@ -1204,7 +1220,7 @@ fun4 <- function(stoch, lhs, sdev, modelInfo, .env){ output_code <- substitute(LHS ~ dnorm(0, sd = SDEV), list(LHS = lhs, SDEV = sdev)) - # Modify the constants. + # Modify the constants modelInfo$constants$new <- 1 list(code = output_code, modelInfo = modelInfo) @@ -1224,6 +1240,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..857417f 100644 --- a/vignettes/nimbleMacros.Rmd.orig +++ b/vignettes/nimbleMacros.Rmd.orig @@ -55,7 +55,7 @@ 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")] @@ -75,7 +75,7 @@ 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({ @@ -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) ``` @@ -203,7 +203,7 @@ Notice that unlike with `LM`, we explicitly define the dimensions of `mu` and `T ```{r} code <- nimbleCode({ - mu[1:N] <- LINPRED(~Time[1:N], priorSpecs=NULL) + mu[1:N] <- LINPRED(~Time[1:N], priorSpecs = NULL) }) ``` @@ -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], priorSpecs = NULL) }) mod <- nimbleModel(code, chick_const2) mod$getCode() @@ -333,6 +333,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 +376,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 `|`. @@ -446,7 +449,7 @@ 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 +481,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 +553,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} @@ -585,6 +598,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} MYDIST <- list(process = fun) @@ -603,8 +617,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 +628,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 +645,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 +661,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 +677,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 +691,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} From 1f79c73b5a9e796dda5b2341d163898c9351a126 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Sun, 16 Feb 2025 16:55:59 -0500 Subject: [PATCH 08/14] Clarify use of factors and LINPRED_PRIORS in vignette --- vignettes/nimbleMacros.Rmd | 4 ++++ vignettes/nimbleMacros.Rmd.orig | 4 ++++ 2 files changed, 8 insertions(+) diff --git a/vignettes/nimbleMacros.Rmd b/vignettes/nimbleMacros.Rmd index 0c656e5..4262799 100644 --- a/vignettes/nimbleMacros.Rmd +++ b/vignettes/nimbleMacros.Rmd @@ -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: @@ -805,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 `priorSpecs` 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. diff --git a/vignettes/nimbleMacros.Rmd.orig b/vignettes/nimbleMacros.Rmd.orig index 857417f..82ebc0c 100644 --- a/vignettes/nimbleMacros.Rmd.orig +++ b/vignettes/nimbleMacros.Rmd.orig @@ -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: @@ -407,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 `priorSpecs` 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. From 61808dfdb842017d43d15e464ad60d6ea029898b Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Mon, 17 Feb 2025 14:51:40 -0500 Subject: [PATCH 09/14] Add test to use later when use3pieces error trapping in nimble is improved --- tests/testthat/test_LINPRED.R | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/tests/testthat/test_LINPRED.R b/tests/testthat/test_LINPRED.R index d56f523..e9787f0 100644 --- a/tests/testthat/test_LINPRED.R +++ b/tests/testthat/test_LINPRED.R @@ -392,6 +392,13 @@ test_that("LINPRED error traps LHS in formula", { 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)) From 9ff8ee00faf35efd04c323773e589c19d6e9e045 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Thu, 20 Feb 2025 08:59:19 -0500 Subject: [PATCH 10/14] Add roxygen for modelInfo and .env arguments --- R/FORLOOP.R | 2 ++ R/LINPRED.R | 5 ++++- R/LM.R | 2 ++ man/FORLOOP.Rd | 4 ++++ man/LINPRED.Rd | 4 ++++ man/LINPRED_PRIORS.Rd | 4 ++++ man/LM.Rd | 4 ++++ 7 files changed, 24 insertions(+), 1 deletion(-) diff --git a/R/FORLOOP.R b/R/FORLOOP.R index 9a438c9..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({ diff --git a/R/LINPRED.R b/R/LINPRED.R index dc846d6..3692a91 100644 --- a/R/LINPRED.R +++ b/R/LINPRED.R @@ -25,6 +25,8 @@ #' Default is \code{FALSE}. #' @param centerVar Grouping variable (covariate) to 'center' the random effects on. By #' default all random effects have mean 0 as with \code{lme4}. +#' @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 #' @@ -116,7 +118,8 @@ unpackArgs=TRUE #' Default is \code{FALSE}. #' @param centerVar Grouping variable (covariate) to 'center' the random effects on. By #' default all random effects have mean 0 as with \code{lme4}. -#' +#' @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 #' diff --git a/R/LM.R b/R/LM.R index eebe27f..3cd2dc6 100644 --- a/R/LM.R +++ b/R/LM.R @@ -24,6 +24,8 @@ #' \code{setPriors()}. #' @param modMatNames 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 #' diff --git a/man/FORLOOP.Rd b/man/FORLOOP.Rd index e24a18e..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 diff --git a/man/LINPRED.Rd b/man/LINPRED.Rd index 0bbde71..8ef949e 100644 --- a/man/LINPRED.Rd +++ b/man/LINPRED.Rd @@ -29,6 +29,10 @@ Default is \code{FALSE}.} \item{centerVar}{Grouping variable (covariate) to 'center' the random effects on. By default all random effects have mean 0 as with \code{lme4}.} + +\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. diff --git a/man/LINPRED_PRIORS.Rd b/man/LINPRED_PRIORS.Rd index 4ce4f6d..a33c4f9 100644 --- a/man/LINPRED_PRIORS.Rd +++ b/man/LINPRED_PRIORS.Rd @@ -24,6 +24,10 @@ Default is \code{FALSE}.} \item{centerVar}{Grouping variable (covariate) to 'center' the random effects on. By default all random effects have mean 0 as with \code{lme4}.} + +\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 diff --git a/man/LM.Rd b/man/LM.Rd index 39ce30c..f5a0163 100644 --- a/man/LM.Rd +++ b/man/LM.Rd @@ -25,6 +25,10 @@ default is no prefix.} \item{modMatNames}{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 From f6ca74fbe1ea216336c7eb47ec0764cee9b67171 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Thu, 20 Feb 2025 09:34:37 -0500 Subject: [PATCH 11/14] Change noncenter argument to noncentered; improve docs for centerVar and noncentered --- R/LINPRED.R | 50 ++++++++++++++++++++++++++--------- man/LINPRED.Rd | 20 +++++++++++--- man/LINPRED_PRIORS.Rd | 22 +++++++++++---- tests/testthat/test_LINPRED.R | 8 +++--- 4 files changed, 74 insertions(+), 26 deletions(-) diff --git a/R/LINPRED.R b/R/LINPRED.R index 3692a91..9be36fd 100644 --- a/R/LINPRED.R +++ b/R/LINPRED.R @@ -21,10 +21,22 @@ #' @param priorSpecs Prior specifications, generated with \code{setPrior()} #' @param modMatNames 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 noncenter Logical indicating whether to use noncentered parameterization. -#' Default is \code{FALSE}. -#' @param centerVar Grouping variable (covariate) to 'center' the random effects on. By -#' default all random effects have mean 0 as with \code{lme4}. +#' @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 #' @@ -44,7 +56,7 @@ NULL LINPRED <- nimble::buildMacro( function(stoch, LHS, formula, link=NULL, coefPrefix=quote(beta_), sdPrefix=NULL, priorSpecs=setPriors(), modMatNames = FALSE, - noncenter = FALSE, centerVar=NULL, modelInfo, .env){ + noncentered = FALSE, centerVar=NULL, modelInfo, .env){ # Make sure formula is in correct format formula <- check_formula(formula) @@ -80,10 +92,10 @@ function(stoch, LHS, formula, link=NULL, coefPrefix=quote(beta_), if(!is.null(priorSpecs)){ priorCode <- substitute(nimbleMacros::LINPRED_PRIORS(FORMULA, coefPrefix=COEFPREFIX, sdPrefix=SDPREFIX, priorSpecs=PRIORSET, modMatNames=MODMAT, - noncenter = UNCENTER, centerVar=CENTERVAR), + noncentered = UNCENTER, centerVar=CENTERVAR), list(COEFPREFIX=coefPrefix, FORMULA=formula, SDPREFIX=sdPrefix, PRIORSET=priorSpecs, MODMAT=modMatNames, - UNCENTER = noncenter, CENTERVAR=centerVar)) + UNCENTER = noncentered, CENTERVAR=centerVar)) code <- embedLinesInCurlyBrackets(list(code, priorCode)) } @@ -114,10 +126,22 @@ unpackArgs=TRUE #' setPriors() #' @param modMatNames 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 noncenter Logical indicating whether to use noncentered parameterization. -#' Default is \code{FALSE}. -#' @param centerVar Grouping variable (covariate) to 'center' the random effects on. By -#' default all random effects have mean 0 as with \code{lme4}. +#' @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 #' @@ -136,7 +160,7 @@ NULL #' @export LINPRED_PRIORS <- nimble::buildMacro( function(formula, coefPrefix=quote(beta_), sdPrefix=NULL, priorSpecs=setPriors(), - modMatNames=FALSE, noncenter = FALSE, centerVar=NULL, modelInfo, .env){ + modMatNames=FALSE, noncentered = FALSE, centerVar=NULL, modelInfo, .env){ # Make sure formula is in correct format formula <- check_formula(formula) @@ -155,7 +179,7 @@ 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) + modMatNames = modMatNames, noncenter=noncentered) # Get complete prior code code <- getPriors(components) diff --git a/man/LINPRED.Rd b/man/LINPRED.Rd index 8ef949e..032c079 100644 --- a/man/LINPRED.Rd +++ b/man/LINPRED.Rd @@ -24,11 +24,23 @@ Default is no prefix.} \item{modMatNames}{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 indicating whether to use noncentered parameterization. -Default is \code{FALSE}.} +\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 variable (covariate) to 'center' the random effects on. By -default all random effects have mean 0 as with \code{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} diff --git a/man/LINPRED_PRIORS.Rd b/man/LINPRED_PRIORS.Rd index a33c4f9..0d3189f 100644 --- a/man/LINPRED_PRIORS.Rd +++ b/man/LINPRED_PRIORS.Rd @@ -19,11 +19,23 @@ setPriors()} \item{modMatNames}{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 indicating whether to use noncentered parameterization. -Default is \code{FALSE}.} - -\item{centerVar}{Grouping variable (covariate) to 'center' the random effects on. By -default all random effects have mean 0 as with \code{lme4}.} +\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 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} diff --git a/tests/testthat/test_LINPRED.R b/tests/testthat/test_LINPRED.R index e9787f0..15b268d 100644 --- a/tests/testthat/test_LINPRED.R +++ b/tests/testthat/test_LINPRED.R @@ -903,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) @@ -941,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) @@ -979,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) @@ -1251,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) From fcaa3662ab5072952d51d132720df8dc1219c323 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Thu, 20 Feb 2025 12:20:48 -0500 Subject: [PATCH 12/14] Change eta to lkjPrior --- R/LINPRED.R | 6 +++--- R/setPriors.R | 6 +++--- man/setPriors.Rd | 4 ++-- tests/testthat/test_LINPRED.R | 4 ++-- tests/testthat/test_priors.R | 12 ++++++------ 5 files changed, 16 insertions(+), 16 deletions(-) diff --git a/R/LINPRED.R b/R/LINPRED.R index 9be36fd..5486e5f 100644 --- a/R/LINPRED.R +++ b/R/LINPRED.R @@ -1355,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 diff --git a/R/setPriors.R b/R/setPriors.R index 6e51fe4..f8998ef 100644 --- a/R/setPriors.R +++ b/R/setPriors.R @@ -15,7 +15,7 @@ #' 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 #' @@ -49,7 +49,7 @@ setPriors <- function(intercept = quote(dnorm(0, sd = 1000)), sd = quote(dunif(0, 100)), factor = NULL, continuous = NULL, - eta = 1, + lkjShape = 1, ...){ # Get specific prior names extra <- list(...) @@ -75,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 } diff --git a/man/setPriors.Rd b/man/setPriors.Rd index 0995bee..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, + lkjShape = 1, ... ) } @@ -27,7 +27,7 @@ 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} diff --git a/tests/testthat/test_LINPRED.R b/tests/testthat/test_LINPRED.R index 15b268d..e476af2 100644 --- a/tests/testthat/test_LINPRED.R +++ b/tests/testthat/test_LINPRED.R @@ -1073,8 +1073,8 @@ 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) }) diff --git a/tests/testthat/test_priors.R b/tests/testthat/test_priors.R index 13b0308..e493e2e 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) + 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) + 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) + 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) + "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) + "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) + sd=quote(dnorm(0, sd = 3)), lkjShape=1) ) # Prior not quoted should error From 35ee4d2b67fa26b6ab6ec1bf12dcc735f7232cab Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Thu, 20 Feb 2025 13:04:56 -0500 Subject: [PATCH 13/14] modMatNames -> modelMatrixNames --- R/LINPRED.R | 14 +++++++------- R/LM.R | 6 +++--- man/LINPRED.Rd | 2 +- man/LINPRED_PRIORS.Rd | 2 +- man/LM.Rd | 2 +- tests/testthat/test_LINPRED.R | 16 ++++++++-------- tests/testthat/test_LM.R | 12 ++++++------ tests/testthat/test_priors.R | 2 +- vignettes/nimbleMacros.Rmd | 4 ++-- vignettes/nimbleMacros.Rmd.orig | 4 ++-- 10 files changed, 32 insertions(+), 32 deletions(-) diff --git a/R/LINPRED.R b/R/LINPRED.R index 5486e5f..1da6a17 100644 --- a/R/LINPRED.R +++ b/R/LINPRED.R @@ -19,7 +19,7 @@ #' @param sdPrefix All dispersion parameters will begin with this prefix. #' Default is no prefix. #' @param priorSpecs Prior specifications, generated with \code{setPrior()} -#' @param modMatNames Logical indicating if parameters should be named so they match the +#' @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, @@ -55,7 +55,7 @@ NULL #' @export LINPRED <- nimble::buildMacro( function(stoch, LHS, formula, link=NULL, coefPrefix=quote(beta_), - sdPrefix=NULL, priorSpecs=setPriors(), modMatNames = FALSE, + sdPrefix=NULL, priorSpecs=setPriors(), modelMatrixNames = FALSE, noncentered = FALSE, centerVar=NULL, modelInfo, .env){ # Make sure formula is in correct format @@ -91,10 +91,10 @@ function(stoch, LHS, formula, link=NULL, coefPrefix=quote(beta_), # Add code for priors to output if needed if(!is.null(priorSpecs)){ priorCode <- substitute(nimbleMacros::LINPRED_PRIORS(FORMULA, coefPrefix=COEFPREFIX, sdPrefix=SDPREFIX, - priorSpecs=PRIORSET, modMatNames=MODMAT, + priorSpecs=PRIORSET, modelMatrixNames=MODMAT, noncentered = UNCENTER, centerVar=CENTERVAR), list(COEFPREFIX=coefPrefix, FORMULA=formula, SDPREFIX=sdPrefix, - PRIORSET=priorSpecs, MODMAT=modMatNames, + PRIORSET=priorSpecs, MODMAT=modelMatrixNames, UNCENTER = noncentered, CENTERVAR=centerVar)) code <- embedLinesInCurlyBrackets(list(code, priorCode)) } @@ -124,7 +124,7 @@ unpackArgs=TRUE #' Default is no prefix. #' @param priorSpecs List of prior specifications, generated using \code{setPriors}. #' setPriors() -#' @param modMatNames Logical indicating if parameters should be named so they match the +#' @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, @@ -160,7 +160,7 @@ NULL #' @export LINPRED_PRIORS <- nimble::buildMacro( function(formula, coefPrefix=quote(beta_), sdPrefix=NULL, priorSpecs=setPriors(), - modMatNames=FALSE, noncentered = FALSE, centerVar=NULL, modelInfo, .env){ + modelMatrixNames=FALSE, noncentered = FALSE, centerVar=NULL, modelInfo, .env){ # Make sure formula is in correct format formula <- check_formula(formula) @@ -179,7 +179,7 @@ function(formula, coefPrefix=quote(beta_), sdPrefix=NULL, priorSpecs=setPriors() components <- buildPriors(components, coefPrefix=safeDeparse(coefPrefix), sdPrefix=sdPrefix, modelInfo = modelInfo, priorSpecs=priorSpecs, - modMatNames = modMatNames, noncenter=noncentered) + modMatNames = modelMatrixNames, noncenter=noncentered) # Get complete prior code code <- getPriors(components) diff --git a/R/LM.R b/R/LM.R index 3cd2dc6..ebda757 100644 --- a/R/LM.R +++ b/R/LM.R @@ -22,7 +22,7 @@ #' default is no prefix. #' @param priorSpecs List of prior specifications, generated using #' \code{setPriors()}. -#' @param modMatNames Logical indicating if parameters be named so they match the +#' @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 @@ -63,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 @@ -114,7 +114,7 @@ LM <- list(process = function(code, modelInfo, .env){ # 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), + modelMatrixNames=MODMAT), list(IDX=idx, FORM=form, COEFPREFIX=coefPrefix, SDPREFIX=sdPrefix, PRIORS=priorSpecs, LINK=link, MODMAT=modMatNames)) pars_added <- list(quote(mu_)) diff --git a/man/LINPRED.Rd b/man/LINPRED.Rd index 032c079..dcac7f5 100644 --- a/man/LINPRED.Rd +++ b/man/LINPRED.Rd @@ -21,7 +21,7 @@ Default is no prefix.} \item{priorSpecs}{Prior specifications, generated with \code{setPrior()}} -\item{modMatNames}{Logical indicating if parameters should be named so they match the +\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{noncentered}{Logical indicating whether to use noncentered parameterization diff --git a/man/LINPRED_PRIORS.Rd b/man/LINPRED_PRIORS.Rd index 0d3189f..95d04d9 100644 --- a/man/LINPRED_PRIORS.Rd +++ b/man/LINPRED_PRIORS.Rd @@ -16,7 +16,7 @@ Default is no prefix.} \item{priorSpecs}{List of prior specifications, generated using \code{setPriors}. setPriors()} -\item{modMatNames}{Logical indicating if parameters should be named so they match the +\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{noncentered}{Logical indicating whether to use noncentered parameterization diff --git a/man/LM.Rd b/man/LM.Rd index f5a0163..c6120e9 100644 --- a/man/LM.Rd +++ b/man/LM.Rd @@ -23,7 +23,7 @@ default is no prefix.} \item{priorSpecs}{List of prior specifications, generated using \code{setPriors()}.} -\item{modMatNames}{Logical indicating if parameters be named so they match the +\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} diff --git a/tests/testthat/test_LINPRED.R b/tests/testthat/test_LINPRED.R index e476af2..1d38ae7 100644 --- a/tests/testthat/test_LINPRED.R +++ b/tests/testthat/test_LINPRED.R @@ -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() diff --git a/tests/testthat/test_LM.R b/tests/testthat/test_LM.R index 7614344..d3fa0d1 100644 --- a/tests/testthat/test_LM.R +++ b/tests/testthat/test_LM.R @@ -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, priorSpecs=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, priorSpecs=pr) }) mod.D93 <- nimbleModel(code.D93, constants=constants) @@ -428,7 +428,7 @@ test_that("Run lme4::glmer() example", { 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) + family=binomial, modelMatrixNames=TRUE) }) mod_gm1 <- nimbleModel(code_gm1, constants=constants) diff --git a/tests/testthat/test_priors.R b/tests/testthat/test_priors.R index e493e2e..7e521db 100644 --- a/tests/testthat/test_priors.R +++ b/tests/testthat/test_priors.R @@ -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 4262799..ed682e9 100644 --- a/vignettes/nimbleMacros.Rmd +++ b/vignettes/nimbleMacros.Rmd @@ -866,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, priorSpecs = pr, modelMatrixNames = TRUE) }) mod <- nimbleModel(code, chick_const) mod$getCode() diff --git a/vignettes/nimbleMacros.Rmd.orig b/vignettes/nimbleMacros.Rmd.orig index 82ebc0c..665fdb8 100644 --- a/vignettes/nimbleMacros.Rmd.orig +++ b/vignettes/nimbleMacros.Rmd.orig @@ -443,11 +443,11 @@ 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, priorSpecs = pr, modelMatrixNames = TRUE) }) mod <- nimbleModel(code, chick_const) mod$getCode() From 1c8c91c505a3b35253ef250f89b36ff344be6b18 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Wed, 5 Mar 2025 11:14:23 -0500 Subject: [PATCH 14/14] Change argument priorSpecs to priors --- R/LINPRED.R | 26 +++++++++++------------ R/LM.R | 12 +++++------ R/setPriors.R | 28 ++++++++++++------------- man/LINPRED.Rd | 2 +- man/LINPRED_PRIORS.Rd | 2 +- man/LM.Rd | 2 +- man/matchPrior.Rd | 12 +++++------ tests/testthat/test_LINPRED.R | 10 ++++----- tests/testthat/test_LM.R | 8 +++---- tests/testthat/test_formulaHandlers.R | 6 +++--- tests/testthat/test_priors.R | 18 ++++++++-------- vignettes/nimbleMacros.Rmd | 30 +++++++++++++-------------- vignettes/nimbleMacros.Rmd.orig | 30 +++++++++++++-------------- 13 files changed, 93 insertions(+), 93 deletions(-) diff --git a/R/LINPRED.R b/R/LINPRED.R index 1da6a17..d16f6e9 100644 --- a/R/LINPRED.R +++ b/R/LINPRED.R @@ -18,7 +18,7 @@ #' 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, generated with \code{setPrior()} +#' @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 @@ -55,7 +55,7 @@ NULL #' @export LINPRED <- nimble::buildMacro( function(stoch, LHS, formula, link=NULL, coefPrefix=quote(beta_), - sdPrefix=NULL, priorSpecs=setPriors(), modelMatrixNames = FALSE, + sdPrefix=NULL, priors=setPriors(), modelMatrixNames = FALSE, noncentered = FALSE, centerVar=NULL, modelInfo, .env){ # Make sure formula is in correct format @@ -89,12 +89,12 @@ 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, modelMatrixNames=MODMAT, + priors=PRIORSET, modelMatrixNames=MODMAT, noncentered = UNCENTER, centerVar=CENTERVAR), list(COEFPREFIX=coefPrefix, FORMULA=formula, SDPREFIX=sdPrefix, - PRIORSET=priorSpecs, MODMAT=modelMatrixNames, + PRIORSET=priors, MODMAT=modelMatrixNames, UNCENTER = noncentered, CENTERVAR=centerVar)) code <- embedLinesInCurlyBrackets(list(code, priorCode)) } @@ -122,7 +122,7 @@ unpackArgs=TRUE #' 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, generated using \code{setPriors}. +#' @param priors List of prior specifications, generated using \code{setPriors}. #' setPriors() #' @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}. @@ -159,14 +159,14 @@ NULL #' @export LINPRED_PRIORS <- nimble::buildMacro( -function(formula, coefPrefix=quote(beta_), sdPrefix=NULL, priorSpecs=setPriors(), +function(formula, coefPrefix=quote(beta_), sdPrefix=NULL, priors=setPriors(), modelMatrixNames=FALSE, noncentered = FALSE, centerVar=NULL, modelInfo, .env){ # Make sure formula is in correct format 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 @@ -178,7 +178,7 @@ function(formula, coefPrefix=quote(beta_), sdPrefix=NULL, priorSpecs=setPriors() components <- buildPriors(components, coefPrefix=safeDeparse(coefPrefix), sdPrefix=sdPrefix, modelInfo = modelInfo, - priorSpecs=priorSpecs, + priorSpecs=priors, modMatNames = modelMatrixNames, noncenter=noncentered) # Get complete prior code @@ -1093,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 @@ -1151,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 @@ -1160,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 @@ -1201,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)) }) diff --git a/R/LM.R b/R/LM.R index ebda757..1eee246 100644 --- a/R/LM.R +++ b/R/LM.R @@ -20,7 +20,7 @@ #' 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, generated using +#' @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}. @@ -99,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"){ @@ -113,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, + 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)) diff --git a/R/setPriors.R b/R/setPriors.R index f8998ef..91c767e 100644 --- a/R/setPriors.R +++ b/R/setPriors.R @@ -117,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 @@ -127,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) @@ -148,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) @@ -158,10 +158,10 @@ 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("matchPrior: Unable to set prior for parameter: ", parName, call.=FALSE) diff --git a/man/LINPRED.Rd b/man/LINPRED.Rd index dcac7f5..6fe8621 100644 --- a/man/LINPRED.Rd +++ b/man/LINPRED.Rd @@ -19,7 +19,7 @@ default is \code{"beta_"} (so 'x' becomes 'beta_x', etc.)} \item{sdPrefix}{All dispersion parameters will begin with this prefix. Default is no prefix.} -\item{priorSpecs}{Prior specifications, generated with \code{setPrior()}} +\item{priors}{Prior specifications, generated with \code{setPrior()}} \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}.} diff --git a/man/LINPRED_PRIORS.Rd b/man/LINPRED_PRIORS.Rd index 95d04d9..4a2d21d 100644 --- a/man/LINPRED_PRIORS.Rd +++ b/man/LINPRED_PRIORS.Rd @@ -13,7 +13,7 @@ default is \code{"beta_"} (so 'x' becomes 'beta_x', etc.)} \item{sdPrefix}{All dispersion parameters will begin with this prefix. Default is no prefix.} -\item{priorSpecs}{List of prior specifications, generated using \code{setPriors}. +\item{priors}{List of prior specifications, generated using \code{setPriors}. setPriors()} \item{modelMatrixNames}{Logical indicating if parameters should be named so they match the diff --git a/man/LM.Rd b/man/LM.Rd index c6120e9..bb56abc 100644 --- a/man/LM.Rd +++ b/man/LM.Rd @@ -20,7 +20,7 @@ default is \code{"beta_"} (so 'x' becomes 'beta_x', etc.)} \item{sdPrefix}{Character. All dispersion parameters will begin with this prefix. default is no prefix.} -\item{priorSpecs}{List of prior specifications, generated using +\item{priors}{List of prior specifications, generated using \code{setPriors()}.} \item{modelMatrixNames}{Logical indicating if parameters be named so they match the diff --git a/man/matchPrior.Rd b/man/matchPrior.Rd index ce0436e..334225e 100644 --- a/man/matchPrior.Rd +++ b/man/matchPrior.Rd @@ -4,17 +4,17 @@ \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{ @@ -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/tests/testthat/test_LINPRED.R b/tests/testthat/test_LINPRED.R index 1d38ae7..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) @@ -473,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) @@ -1076,7 +1076,7 @@ test_that("correlated random effects", { # 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) @@ -1402,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 d3fa0d1..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,7 +36,7 @@ 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)"), modelMatrixNames=TRUE)) + code <- quote(LM(y ~ x + x2, priors=setPriors(sd="dunif(0, 5)"), modelMatrixNames=TRUE)) out <- LM$process(code, modelInfo, environment()) code <- nimbleCode({ @@ -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(), - modelMatrixNames=TRUE, priorSpecs=pr) + modelMatrixNames=TRUE, priors=pr) }) mod.D93 <- nimbleModel(code.D93, constants=constants) @@ -427,7 +427,7 @@ 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, + 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 19fc1bd..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") }) diff --git a/tests/testthat/test_priors.R b/tests/testthat/test_priors.R index 7e521db..bf352d6 100644 --- a/tests/testthat/test_priors.R +++ b/tests/testthat/test_priors.R @@ -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", { diff --git a/vignettes/nimbleMacros.Rmd b/vignettes/nimbleMacros.Rmd index ed682e9..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) ``` @@ -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))) ``` @@ -101,7 +101,7 @@ Two things to note: ``` r code <- nimbleCode({ - LM(am ~ scale(mpg), family = binomial, priorSpecs=pr) + LM(am ~ scale(mpg), family = binomial, priors = pr) }) ``` @@ -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() @@ -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() @@ -513,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() @@ -545,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() @@ -572,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() @@ -608,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). @@ -806,7 +806,7 @@ 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 `priorSpecs` argument above). +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). @@ -837,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() @@ -871,7 +871,7 @@ We can replicate this structure in NIMBLE model code by setting the `LINPRED_PRI ``` r code <- nimbleCode({ - LINPRED_PRIORS(~Time + Diet, priorSpecs = pr, modelMatrixNames = TRUE) + LINPRED_PRIORS(~Time + Diet, priors = pr, modelMatrixNames = TRUE) }) mod <- nimbleModel(code, chick_const) mod$getCode() @@ -1088,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()) ``` ``` diff --git a/vignettes/nimbleMacros.Rmd.orig b/vignettes/nimbleMacros.Rmd.orig index 665fdb8..e3e6be8 100644 --- a/vignettes/nimbleMacros.Rmd.orig +++ b/vignettes/nimbleMacros.Rmd.orig @@ -51,7 +51,7 @@ 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) ``` @@ -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))) ``` @@ -79,7 +79,7 @@ Two things to note: ```{r} code <- nimbleCode({ - LM(am ~ scale(mpg), family = binomial, priorSpecs=pr) + LM(am ~ scale(mpg), family = binomial, priors = pr) }) ``` @@ -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() @@ -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() @@ -268,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() @@ -287,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() @@ -301,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() @@ -320,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} @@ -408,7 +408,7 @@ 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 `priorSpecs` argument above). +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). @@ -430,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() @@ -447,7 +447,7 @@ We can replicate this structure in NIMBLE model code by setting the `LINPRED_PRI ```{r} code <- nimbleCode({ - LINPRED_PRIORS(~Time + Diet, priorSpecs = pr, modelMatrixNames = TRUE) + LINPRED_PRIORS(~Time + Diet, priors = pr, modelMatrixNames = TRUE) }) mod <- nimbleModel(code, chick_const) mod$getCode() @@ -597,7 +597,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()) ``` Finally, `nimble` needs a way to know that this function is intended to be a macro, and that it is called `MYDIST`.