From b558fd98668dbba8c41762429dd77da92725ea38 Mon Sep 17 00:00:00 2001
From: Ken Kellner <ken@kenkellner.com>
Date: Wed, 15 Jan 2025 12:07:07 -0500
Subject: [PATCH 1/9] Add formula function processing and support for offset()

---
 NAMESPACE                              |   1 +
 R/LINPRED.R                            |  60 ++++++++++++--
 R/formulaFunction.R                    |  69 ++++++++++++++++
 man/offsetFormulaFunction.Rd           |  12 +++
 tests/testthat/test_LINPRED.R          |  18 -----
 tests/testthat/test_formulaFunctions.R | 108 +++++++++++++++++++++++++
 6 files changed, 243 insertions(+), 25 deletions(-)
 create mode 100644 R/formulaFunction.R
 create mode 100644 man/offsetFormulaFunction.Rd
 create mode 100644 tests/testthat/test_formulaFunctions.R

diff --git a/NAMESPACE b/NAMESPACE
index 7bf36ec..4c9e164 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -5,6 +5,7 @@ export(LINPRED)
 export(LINPRED_PRIORS)
 export(LM)
 export(matchPrior)
+export(offsetFormulaFunction)
 export(setPriors)
 export(uppertri_mult_diag)
 importFrom(nimble,nimMatrix)
diff --git a/R/LINPRED.R b/R/LINPRED.R
index ff603ed..9b22aba 100644
--- a/R/LINPRED.R
+++ b/R/LINPRED.R
@@ -50,7 +50,7 @@ function(stoch, LHS, formula, link=NULL, coefPrefix=quote(beta_),
     # Split formula into components and process the components
     components <- buildLP(formula, defaultBracket = LHS_ind, 
                     coefPrefix = safeDeparse(coefPrefix), sdPrefix=sdPrefix,
-                    modelInfo = modelInfo, centerVar = centerVar)
+                    modelInfo = modelInfo, centerVar = centerVar, env=.env)
     
     # Update modelInfo
     modelInfo <- updateModelInfo(modelInfo, components)
@@ -138,7 +138,7 @@ function(formula, coefPrefix=quote(beta_), sdPrefix=NULL, priorSpecs=setPriors()
   # Split formula into components and process the components
   components <- buildLP(formula, defaultBracket = "[]", # default bracket not used below 
                     coefPrefix = safeDeparse(coefPrefix), sdPrefix=sdPrefix,
-                    modelInfo = modelInfo, centerVar = centerVar)
+                    modelInfo = modelInfo, centerVar = centerVar, env=.env)
     
   # Update constants in modelInfo
   modelInfo <- updateModelInfo(modelInfo, components)
@@ -169,11 +169,11 @@ unpackArgs=TRUE
 # Splits the formula into separate components (terms)
 # Then iteratively adds information to each component until eventually
 # the linear predictor code for the components can be added
-buildLP <- function(formula, defaultBracket, coefPrefix="beta_", sdPrefix=NULL, modelInfo, centerVar=NULL){
+buildLP <- function(formula, defaultBracket, coefPrefix="beta_", sdPrefix=NULL, modelInfo, centerVar=NULL, env){
   comps <- separateFormulaComponents(formula)
-  # Functions will be handled here; for now an error
-  is_function <- sapply(comps, function(x) inherits(x, "formulaComponentFunction"))
-  if(any(is_function)) stop("Functions in formulas not yet supported", call.=FALSE)
+  # Process formula functions, erroring if an unsupported one is found
+  comps <- lapply(comps, processFormulaFunction, defaultBracket = defaultBracket,
+                  coefPrefix=coefPrefix, sdPrefix=sdPrefix, modelInfo = modelInfo, env = env)
   comps <- lapply(comps, addTermsAndBrackets, defaultBracket = defaultBracket, constants = modelInfo$constants)
   # Update constants in modelInfo with any new constants before moving on
   # constants may have been created by addTermsAndBrackets if there were nested random effects
@@ -294,7 +294,16 @@ createInterceptComponent <- function(formula){
 # presence of (), e.g. scale(), I()
 createFixedComponents <- function(formula){
   fixed <- reformulas::nobars(formula) # remove random terms first
-  fixed_terms <- attr(stats::terms(fixed), "term.labels")
+  trms <- stats::terms(fixed)
+  fixed_terms <- attr(trms, "term.labels")
+  
+  # Handle offset special case
+  off <- attr(trms, "offset")
+  if(!is.null(off)){
+    off <- rownames(attr(trms, "factors"))[off]
+    fixed_terms <- c(fixed_terms, off)
+  }
+
   if(length(fixed_terms) == 0) return(NULL)
   components <- lapply(fixed_terms, function(x){
     if(grepl("(", x, fixed=TRUE)){ # is this a function component?
@@ -322,6 +331,43 @@ createRandomComponents <- function(formula){
 }
 
 
+# processFormulaFunctions------------------------------------------------------
+# For formulaComponentFunction objects
+# If the formula component is FUN(), then this function looks for
+# a corresponding formulaFunction class object in the environment called 
+# "FUNFormulaFunction". This function is then run on the component, which
+# should return a component with updated linear predictor and priors slots
+# to be used in final code compilation.
+# It is also possible that the processing function could do some work and
+# then change the class of the output component to formulaComponentFixed so
+# that it will be further processed in later steps.
+processFormulaFunction <- function(x, defaultBracket, coefPrefix="beta_", 
+                                    sdPrefix=NULL, modelInfo, env, ...){
+
+  if(!inherits(x, "formulaComponentFunction")) return(x)
+
+  func <- safeDeparse(x$lang[[1]])
+  cand <- paste0(func, "FormulaFunction")
+  
+  processor_available <- FALSE
+  if(exists(cand, envir = env)){
+    processor <- get(cand, envir = env)
+    if(inherits(processor, "formulaFunction")){
+      processor_available <- TRUE
+      out <- processor$process(x, defaultBracket, coefPrefix, sdPrefix, modelInfo, env, ...)
+      if(!inherits(out, "formulaComponent")){
+        stop("Processing function doesn't return formulaComponent object", call.=FALSE)
+      }
+    }
+  }
+  if(!processor_available){
+    stop("No processing function for ", func, "() available", call.=FALSE)
+  }
+
+  out
+}
+
+
 # addTermsAndBrackets----------------------------------------------------------
 # Takes formula component, adds formula terms to term slot and splits out brackets into bracket slot
 # Returns the updated component
diff --git a/R/formulaFunction.R b/R/formulaFunction.R
new file mode 100644
index 0000000..5fd8b32
--- /dev/null
+++ b/R/formulaFunction.R
@@ -0,0 +1,69 @@
+#' 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.
+#'
+#' @name offsetFormulaFunction
+#'
+#' @param x A formulaComponentFunction object created from an offset() term
+#'
+#'
+NULL
+
+#' @export
+offsetFormulaFunction <- list(process = 
+  function(x, defaultBracket, coefPrefix, sdPrefix, modelInfo, env, ...){
+ 
+  # Get the code inside offset()
+  interior <- x$lang[[2]]
+
+  # Extract any brackets
+  brack <- extractAllBrackets(interior)
+  if(is.null(brack)) brack <- defaultBracket
+  if(length(brack) > 1) stop("Can't handle multiple brackets inside offset()", call.=FALSE)
+
+  # Remove brackets from code inside offset
+  interior <- removeSquareBrackets(interior)
+
+  # Handle any additional functions inside offset (such as log)
+  if(length(interior) > 1){ # If >1, there must be another function call
+    # Get the name of the function
+    func <- interior[[1]]
+    # Get all names inside the call and create the name of the new constant
+    all_names <- get_all_names(interior)
+    const_name <- sapply(all_names, safeDeparse)
+    const_name <- paste(const_name, collapse="_")
+    # Try to evaluate the code
+    new_const <- try(eval(interior, modelInfo$constants))
+    if(inherits(new_const, "try-error")){
+      stop("Problem evaluating function inside offset()", call.=FALSE)
+    }
+    # Add the new constant to the output
+    add_const <- list(new_const)
+    names(add_const) <- const_name
+    if(is.null(x$constants)) x$constants <- list()
+    x$constants <- modifyList(x$constants, add_const)
+    # Update the name of the term to put in the linear predictor
+    interior <- str2lang(const_name)
+  }
+  
+  # Combine the constant name and bracket to create the linear predictor term
+  code <- paste0(safeDeparse(interior), brack)  
+  x$linPredCode <- code
+  x$priorCode <- NULL   # no prior code for an offset
+  x
+})
+class(offsetFormulaFunction) <- "formulaFunction"
+
+# Extract all 'names' from a call
+# So qnorm(log(x)) would yield qnorm, log, x
+get_all_names <- function(code){
+  out <- get_all_names_recursive(code)
+  unlist(out)
+}
+
+get_all_names_recursive <- function(code){
+  if(is.name(code)) return(code)
+  if(is.call(code)) return(lapply(code, get_all_names))
+  NULL
+}
diff --git a/man/offsetFormulaFunction.Rd b/man/offsetFormulaFunction.Rd
new file mode 100644
index 0000000..36055a6
--- /dev/null
+++ b/man/offsetFormulaFunction.Rd
@@ -0,0 +1,12 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/formulaFunction.R
+\name{offsetFormulaFunction}
+\alias{offsetFormulaFunction}
+\title{Function to handle offset() in LINPRED}
+\arguments{
+\item{x}{A formulaComponentFunction object created from an offset() term}
+}
+\description{
+Translates offset() in an R formula passed to LINPRED into corresponding
+NIMBLE code for an offset in the linear predictor.
+}
diff --git a/tests/testthat/test_LINPRED.R b/tests/testthat/test_LINPRED.R
index 0341944..36334fc 100644
--- a/tests/testthat/test_LINPRED.R
+++ b/tests/testthat/test_LINPRED.R
@@ -1366,24 +1366,6 @@ test_that("LINPRED with factor array covariate", {
 })
 
 
-test_that("LINPRED errors when there are functions in the formula", {
-  set.seed(123)
-  modInfo <- list(constants=list(y = rnorm(10), x=factor(sample(letters[1:3], 10, replace=T)),
-                    x2=factor(sample(letters[4:5], 10, replace=T)),
-                    x3=round(rnorm(10),3)))
-
-  code <- quote(y[1:n] <- LINPRED(~scale(x), priorSpecs=NULL))
-  expect_error(LINPRED$process(code, modelInfo=modInfo, .env=NULL), "Functions")
-
-  code <- quote(y[1:n] <- LINPRED(~scale(x) + (1|x2), priorSpecs=NULL))
-  expect_error(LINPRED$process(code, modelInfo=modInfo, .env=NULL), "Functions")
-
-  code <- quote(y[1:n] <- LINPRED(~x3 + I(x[1:10]), priorSpecs=NULL))
-  expect_error(LINPRED$process(code, modelInfo=modInfo, .env=NULL), "Functions")
-
-})
-
-
 test_that("Nested random effects", {
   nimbleOptions(enableMacroComments = FALSE)
   set.seed(123)
diff --git a/tests/testthat/test_formulaFunctions.R b/tests/testthat/test_formulaFunctions.R
new file mode 100644
index 0000000..35f90bb
--- /dev/null
+++ b/tests/testthat/test_formulaFunctions.R
@@ -0,0 +1,108 @@
+context("formulaFunctions")
+
+test_that("Error when function formula is unsupported", {
+  set.seed(123)
+  modInfo <- list(constants=list(y = rnorm(10), x=factor(sample(letters[1:3], 10, replace=T)),
+                    x2=factor(sample(letters[4:5], 10, replace=T)),
+                    x3=round(rnorm(10),3)))
+
+  code <- quote(y[1:n] <- LINPRED(~test(x), priorSpecs=NULL))
+  expect_error(LINPRED$process(code, modelInfo=modInfo, .env=environment()), "No processing")
+
+  code <- quote(y[1:n] <- LINPRED(~test(x) + (1|x2), priorSpecs=NULL))
+  expect_error(LINPRED$process(code, modelInfo=modInfo, .env=environment()), "No processing")
+
+  code <- quote(y[1:n] <- LINPRED(~x3 + test(x[1:10]), priorSpecs=NULL))
+  expect_error(LINPRED$process(code, modelInfo=modInfo, .env=environment()), "No processing")
+
+})
+
+test_that("offset formula function works", {
+  nimbleOptions(enableMacroComments=FALSE)
+  set.seed(123)
+  constants <- list(y = rnorm(5), x = rnorm(5), z = runif(5, 0, 1), 
+                  z2 = matrix(runif(5, 0, 1), 5, 1), n = 5)
+
+  # Basic offset
+  code <- nimbleCode({
+    mu[1:n] <- LINPRED(~x + offset(z))
+  })
+  mod <- nimbleModel(code, constants=constants)
+  
+  expect_equal(
+    mod$getCode(),
+    quote({
+    for (i_1 in 1:n) {
+        mu[i_1] <- beta_Intercept + beta_x * x[i_1] + z[i_1]
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_x ~ dnorm(0, sd = 1000)
+    })
+  )
+
+  # Bracket
+  code <- nimbleCode({
+    mu[1:n] <- LINPRED(~x + offset(z2[1:n,1]))
+  })
+  mod <- nimbleModel(code, constants=constants)
+  
+  expect_equal(
+    mod$getCode(),
+    quote({
+    for (i_1 in 1:n) {
+        mu[i_1] <- beta_Intercept + beta_x * x[i_1] + z2[i_1, 1]
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_x ~ dnorm(0, sd = 1000)
+    })
+  )
+
+  # Nested functions inside offset
+  code <- nimbleCode({
+    mu[1:n] <- LINPRED(~x + offset(pnorm(log(z), 0)))
+  })
+  mod <- nimbleModel(code, constants=constants)
+  
+  expect_equal(
+    mod$getCode(),
+    quote({
+    for (i_1 in 1:n) {
+        mu[i_1] <- beta_Intercept + beta_x * x[i_1] + pnorm_log_z[i_1]
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_x ~ dnorm(0, sd = 1000)
+    })
+  )
+  # This should create a new constant
+  expect_equal(mod$getConstants()$pnorm_log_z,
+               pnorm(log(constants$z), 0))
+ 
+  # Nested functions inside offset with bracket
+  code <- nimbleCode({
+    mu[1:n] <- LINPRED(~x + offset(pnorm(log(z2[1:n,1]), 0)))
+  })
+  mod <- nimbleModel(code, constants=constants)
+  
+  expect_equal(
+    mod$getCode(),
+    quote({
+    for (i_1 in 1:n) {
+        mu[i_1] <- beta_Intercept + beta_x * x[i_1] + pnorm_log_z2[i_1,1]
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_x ~ dnorm(0, sd = 1000)
+    })
+  )
+  # This should create a new constant
+  expect_equal(mod$getConstants()$pnorm_log_z2,
+               pnorm(log(constants$z2), 0))
+
+  # Error if function eval doesn't work
+  code <- quote(mu[1:n] <- LINPRED(~x + offset(test(z))))
+  options(show.error.messages=FALSE)
+  expect_error(LINPRED$process(code, modelInfo=list(constants=constants), environment()),
+               "Problem evaluating")
+  options(show.error.messages=TRUE)
+
+  nimbleOptions(enableMacroComments=TRUE)
+})

From 2f0ec60c207500dad88a06d61eba602a8f054c40 Mon Sep 17 00:00:00 2001
From: Ken Kellner <ken@kenkellner.com>
Date: Thu, 16 Jan 2025 10:40:27 -0500
Subject: [PATCH 2/9] Handle interactions with formula functions

---
 R/LINPRED.R         | 18 +++++++++++++-
 R/formulaFunction.R | 58 +++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 75 insertions(+), 1 deletion(-)

diff --git a/R/LINPRED.R b/R/LINPRED.R
index 9b22aba..70a77cc 100644
--- a/R/LINPRED.R
+++ b/R/LINPRED.R
@@ -346,7 +346,17 @@ processFormulaFunction <- function(x, defaultBracket, coefPrefix="beta_",
 
   if(!inherits(x, "formulaComponentFunction")) return(x)
 
-  func <- safeDeparse(x$lang[[1]])
+  # Identify the first function that appears in the term (in case it's an interaction)
+  trms <- splitInteractionTerms(x$lang)
+  trms <- trms[!sapply(trms, is.name)]
+  funcs <- sapply(trms, function(x) x[[1]])
+
+  if(any(funcs[[1]] != funcs)){
+    stop("Interactions with multiple different formula functions not supported", call.=FALSE)
+  }
+
+  func <- safeDeparse(funcs[[1]])
+
   cand <- paste0(func, "FormulaFunction")
   
   processor_available <- FALSE
@@ -367,6 +377,12 @@ processFormulaFunction <- function(x, defaultBracket, coefPrefix="beta_",
   out
 }
 
+splitInteractionTerms <- function(code){
+  out <- removeSquareBrackets(code)
+  out <- safeDeparse(out)
+  out <- strsplit(out, ":")[[1]]
+  lapply(out, str2lang)
+}
 
 # addTermsAndBrackets----------------------------------------------------------
 # Takes formula component, adds formula terms to term slot and splits out brackets into bracket slot
diff --git a/R/formulaFunction.R b/R/formulaFunction.R
index 5fd8b32..3009a94 100644
--- a/R/formulaFunction.R
+++ b/R/formulaFunction.R
@@ -67,3 +67,61 @@ get_all_names_recursive <- function(code){
   if(is.call(code)) return(lapply(code, get_all_names))
   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.
+#'
+#' @name scaleFormulaFunction
+#'
+#' @param x A formulaComponentFunction object created from an scale() term
+#'
+#'
+NULL
+
+#' @export
+offsetFormulaFunction <- list(process = 
+  function(x, defaultBracket, coefPrefix, sdPrefix, modelInfo, env, ...){
+ 
+  # Get the code inside offset()
+  interior <- x$lang[[2]]
+
+  # Extract any brackets
+  brack <- extractAllBrackets(interior)
+  if(is.null(brack)) brack <- defaultBracket
+  if(length(brack) > 1) stop("Can't handle multiple brackets inside offset()", call.=FALSE)
+
+  # Remove brackets from code inside offset
+  interior <- removeSquareBrackets(interior)
+
+  # Handle any additional functions inside offset (such as log)
+  if(length(interior) > 1){ # If >1, there must be another function call
+    # Get the name of the function
+    func <- interior[[1]]
+    # Get all names inside the call and create the name of the new constant
+    all_names <- get_all_names(interior)
+    const_name <- sapply(all_names, safeDeparse)
+    const_name <- paste(const_name, collapse="_")
+    # Try to evaluate the code
+    new_const <- try(eval(interior, modelInfo$constants))
+    if(inherits(new_const, "try-error")){
+      stop("Problem evaluating function inside offset()", call.=FALSE)
+    }
+    # Add the new constant to the output
+    add_const <- list(new_const)
+    names(add_const) <- const_name
+    if(is.null(x$constants)) x$constants <- list()
+    x$constants <- modifyList(x$constants, add_const)
+    # Update the name of the term to put in the linear predictor
+    interior <- str2lang(const_name)
+  }
+  
+  # Combine the constant name and bracket to create the linear predictor term
+  code <- paste0(safeDeparse(interior), brack)  
+  x$linPredCode <- code
+  x$priorCode <- NULL   # no prior code for an offset
+  x
+})
+class(offsetFormulaFunction) <- "formulaFunction"
+

From 6087c9e3da91bf340b51cd7c7312a3c2c29b3c3a Mon Sep 17 00:00:00 2001
From: Ken Kellner <ken@kenkellner.com>
Date: Thu, 16 Jan 2025 11:41:50 -0500
Subject: [PATCH 3/9] Add scale() support to LINPRED

---
 NAMESPACE                              |   1 +
 R/formulaFunction.R                    |  94 ++++++++-----
 man/scaleFormulaFunction.Rd            |  12 ++
 tests/testthat/test_formulaFunctions.R | 186 +++++++++++++++++++++++++
 4 files changed, 257 insertions(+), 36 deletions(-)
 create mode 100644 man/scaleFormulaFunction.Rd

diff --git a/NAMESPACE b/NAMESPACE
index 4c9e164..af1261a 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -6,6 +6,7 @@ export(LINPRED_PRIORS)
 export(LM)
 export(matchPrior)
 export(offsetFormulaFunction)
+export(scaleFormulaFunction)
 export(setPriors)
 export(uppertri_mult_diag)
 importFrom(nimble,nimMatrix)
diff --git a/R/formulaFunction.R b/R/formulaFunction.R
index 3009a94..79ee43f 100644
--- a/R/formulaFunction.R
+++ b/R/formulaFunction.R
@@ -75,53 +75,75 @@ get_all_names_recursive <- function(code){
 #'
 #' @name scaleFormulaFunction
 #'
-#' @param x A formulaComponentFunction object created from an scale() term
+#' @param x A formulaComponentFunction object created from a scale() term
 #'
 #'
 NULL
 
 #' @export
-offsetFormulaFunction <- list(process = 
+scaleFormulaFunction <- list(process = 
   function(x, defaultBracket, coefPrefix, sdPrefix, modelInfo, env, ...){
- 
-  # Get the code inside offset()
-  interior <- x$lang[[2]]
 
-  # Extract any brackets
-  brack <- extractAllBrackets(interior)
-  if(is.null(brack)) brack <- defaultBracket
-  if(length(brack) > 1) stop("Can't handle multiple brackets inside offset()", call.=FALSE)
+  # Identify which interaction terms involve scale
+  trms <- splitInteractionTerms(x$lang)
+  has_scale <- !sapply(trms, is.name)
 
-  # Remove brackets from code inside offset
-  interior <- removeSquareBrackets(interior)
+  bracks <- extractAllBrackets(x$lang)
 
-  # Handle any additional functions inside offset (such as log)
-  if(length(interior) > 1){ # If >1, there must be another function call
-    # Get the name of the function
-    func <- interior[[1]]
-    # Get all names inside the call and create the name of the new constant
-    all_names <- get_all_names(interior)
-    const_name <- sapply(all_names, safeDeparse)
-    const_name <- paste(const_name, collapse="_")
-    # Try to evaluate the code
-    new_const <- try(eval(interior, modelInfo$constants))
-    if(inherits(new_const, "try-error")){
-      stop("Problem evaluating function inside offset()", call.=FALSE)
+  for (i in 1:length(trms)){
+    if(!has_scale[i]){
+      trm_idx <- safeDeparse(trms[[i]])
+      if(trm_idx %in% names(bracks)){
+        brack <- bracks[[trm_idx]]
+      } else {
+        brack <- defaultBracket
+      }
+      trms[[i]] <- str2lang(paste0(safeDeparse(trms[[i]]), brack))
+    } else {
+      interior <- trms[[i]][[2]]
+      if(!is.name(interior)) stop("Can't handle expression inside scale()", call.=FALSE)
+
+      # Extract any brackets
+      trm_idx <- safeDeparse(interior)
+      if(trm_idx %in% names(bracks)){
+        brack <- bracks[[trm_idx]]
+      } else {
+        brack <- defaultBracket
+      }
+
+      # Remove brackets from code
+      interior <- removeSquareBrackets(interior)
+  
+      # New term
+      new_term <- paste0(safeDeparse(interior), "_scaled")
+
+      # Get constant
+      const <- modelInfo$constants[[safeDeparse(interior)]]
+      if(is.null(const)) stop("Covariate inside scale() is missing from constants", call.=FALSE)
+      if(!is.numeric(const)) stop("Covariate inside scale() must be numeric", call.=FALSE)
+
+      # Make sure constant retains dimensions after being scaled
+      out <- as.numeric(scale(const))
+      if(!is.null(dim(const))){
+        out <- array(out, dim=dim(const))
+      }
+
+      add_const <- list(out)
+      names(add_const) <- new_term 
+
+      # Add the new constant to the object
+      if(is.null(x$constants)) x$constants <- list()
+      x$constants <- modifyList(x$constants, add_const)
+      trms[[i]] <- str2lang(paste0(new_term, brack))
     }
-    # Add the new constant to the output
-    add_const <- list(new_const)
-    names(add_const) <- const_name
-    if(is.null(x$constants)) x$constants <- list()
-    x$constants <- modifyList(x$constants, add_const)
-    # Update the name of the term to put in the linear predictor
-    interior <- str2lang(const_name)
   }
+
+  # Updated interaction
+  x$lang <- str2lang(paste(sapply(trms, safeDeparse), collapse=":"))
+
+  # Update class
+  class(x)[1] <- "formulaComponentFixed"
   
-  # Combine the constant name and bracket to create the linear predictor term
-  code <- paste0(safeDeparse(interior), brack)  
-  x$linPredCode <- code
-  x$priorCode <- NULL   # no prior code for an offset
   x
 })
-class(offsetFormulaFunction) <- "formulaFunction"
-
+class(scaleFormulaFunction) <- "formulaFunction"
diff --git a/man/scaleFormulaFunction.Rd b/man/scaleFormulaFunction.Rd
new file mode 100644
index 0000000..b1db1ea
--- /dev/null
+++ b/man/scaleFormulaFunction.Rd
@@ -0,0 +1,12 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/formulaFunction.R
+\name{scaleFormulaFunction}
+\alias{scaleFormulaFunction}
+\title{Function to handle scale() in LINPRED}
+\arguments{
+\item{x}{A formulaComponentFunction object created from a scale() term}
+}
+\description{
+Translates scale() in an R formula passed to LINPRED into corresponding
+NIMBLE code (and new constant) for a scaled covariate.
+}
diff --git a/tests/testthat/test_formulaFunctions.R b/tests/testthat/test_formulaFunctions.R
index 35f90bb..3f03d38 100644
--- a/tests/testthat/test_formulaFunctions.R
+++ b/tests/testthat/test_formulaFunctions.R
@@ -106,3 +106,189 @@ test_that("offset formula function works", {
 
   nimbleOptions(enableMacroComments=TRUE)
 })
+
+test_that("scale formula function works", {
+  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))
+
+  # Basic scale
+  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_scaled * x_scaled[i_1] + 
+            beta_z * z[i_1]
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_x_scaled ~ dnorm(0, sd = 1000)
+    beta_z ~ dnorm(0, sd = 1000)
+    })
+  )
+  expect_equivalent(
+    mod$getConstants()$x_scaled,
+    scale(constants$x)
+  )
+  # x is not removed from constants
+  expect_equal(
+    mod$getConstants()$x,
+    constants$x
+  )
+
+  # Interaction with non-scaled term
+  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_scaled_z * x_scaled[i_1] * 
+            z[i_1]
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_x_scaled_z ~ dnorm(0, sd = 1000)
+    })
+  )
+  expect_equivalent(
+    mod$getConstants()$x_scaled,
+    scale(constants$x)
+  )
+
+  # Interaction and linear term
+  code <- nimbleCode({
+    mu[1:n] <- LINPRED(~scale(x) + 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_scaled * x_scaled[i_1] + 
+            beta_x_scaled_z * x_scaled[i_1] * z[i_1]
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_x_scaled ~ dnorm(0, sd = 1000)
+    beta_x_scaled_z ~ dnorm(0, sd = 1000)
+    })
+  )
+  expect_equivalent(
+    mod$getConstants()$x_scaled,
+    scale(constants$x)
+  )
+
+  # Non-scaled term comes first in interaction
+  code <- nimbleCode({
+    mu[1:n] <- LINPRED(~z:scale(x)) 
+  })
+  mod <- nimbleModel(code, constants=constants)
+  expect_equal(
+    mod$getCode(),
+    quote({
+    for (i_1 in 1:n) {
+        mu[i_1] <- beta_Intercept + beta_z_x_scaled * z[i_1] * 
+            x_scaled[i_1]
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_z_x_scaled ~ dnorm(0, sd = 1000)
+    })
+  )
+  expect_equivalent(
+    mod$getConstants()$x_scaled,
+    scale(constants$x)
+  )
+
+  # Interaction of two scaled terms
+  code <- nimbleCode({
+    mu[1:n] <- LINPRED(~scale(x):scale(z)) 
+  })
+  mod <- nimbleModel(code, constants=constants)
+  expect_equal(
+    mod$getCode(),
+    quote({
+    for (i_1 in 1:n) {
+        mu[i_1] <- beta_Intercept + beta_x_scaled_z_scaled * 
+            x_scaled[i_1] * z_scaled[i_1]
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_x_scaled_z_scaled ~ dnorm(0, sd = 1000)
+    })
+  )
+  expect_equivalent(
+    mod$getConstants()$x_scaled,
+    scale(constants$x)
+  )
+  expect_equivalent(
+    mod$getConstants()$z_scaled,
+    scale(constants$z)
+  )
+
+  # Scale with brackets
+  code <- nimbleCode({
+    mu[1:n] <- LINPRED(~scale(x2[1:n,1]):z)
+  })
+  mod <- nimbleModel(code, constants=constants)
+  expect_equal(
+    mod$getCode(),
+    quote({
+    for (i_1 in 1:n) {
+        mu[i_1] <- beta_Intercept + beta_x2_scaled_z * x2_scaled[i_1, 
+            1] * z[i_1]
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_x2_scaled_z ~ dnorm(0, sd = 1000)
+  })
+  )
+  expect_equivalent(
+    mod$getConstants()$x2_scaled,
+    scale(constants$x2)
+  )
+  expect_equal(dim(mod$getConstants()$x2_scaled), dim(constants$x2))
+
+  # Scale with brackets on the other interaction term
+  code <- nimbleCode({
+    mu[1:n] <- LINPRED(~scale(x2[1:n,1]):z2[1:n,1])
+  })
+  mod <- nimbleModel(code, constants=constants)
+  expect_equal(
+    mod$getCode(),
+    quote({
+    for (i_1 in 1:n) {
+        mu[i_1] <- beta_Intercept + beta_x2_scaled_z2 * x2_scaled[i_1, 
+            1] * z2[i_1, 1]
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_x2_scaled_z2 ~ dnorm(0, sd = 1000)
+  })
+  )
+
+  # Error if covariate is not in constants
+  code <- quote(mu[1:n] <- LINPRED(~scale(test)))
+  expect_error(
+    LINPRED$process(code, list(constants=constants), environment()),
+    "Covariate inside"
+  )
+
+  # Error if expression inside scale
+  code <- quote(mu[1:n] <- LINPRED(~scale(x*x)))
+  expect_error(
+    LINPRED$process(code, list(constants=constants), environment()),
+    "expression"
+  )
+
+  # Error if multiple types of functions in an interaction
+  code <- quote(mu[1:n] <- LINPRED(~scale(x):test(z)))
+  expect_error(
+    LINPRED$process(code, list(constants=constants), environment()),
+    "multiple different formula functions"
+  )
+
+  nimbleOptions(enableMacroComments=TRUE)
+})

From ed3d93eed510c30183ad66a344d90990bb849848 Mon Sep 17 00:00:00 2001
From: Ken Kellner <ken@kenkellner.com>
Date: Thu, 16 Jan 2025 14:02:29 -0500
Subject: [PATCH 4/9] Add I() support to LINPRED

---
 DESCRIPTION                            |   2 +-
 NAMESPACE                              |   1 +
 R/formulaFunction.R                    |  77 +++++++++++++++
 man/IFormulaFunction.Rd                |  13 +++
 tests/testthat/test_formulaFunctions.R | 124 +++++++++++++++++++++++++
 5 files changed, 216 insertions(+), 1 deletion(-)
 create mode 100644 man/IFormulaFunction.Rd

diff --git a/DESCRIPTION b/DESCRIPTION
index a8f9645..6ccb48a 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: nimbleMacros
 Version: 0.0.12
-Date: 2025-01-14
+Date: 2025-01-16
 Title: Macros for 'nimble' Code
 Authors@R: person("Ken", "Kellner", email="contact@kenkellner.com", 
                   role=c("cre","aut"))
diff --git a/NAMESPACE b/NAMESPACE
index af1261a..f85a273 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 export(FORLOOP)
+export(IFormulaFunction)
 export(LINPRED)
 export(LINPRED_PRIORS)
 export(LM)
diff --git a/R/formulaFunction.R b/R/formulaFunction.R
index 79ee43f..a5d284e 100644
--- a/R/formulaFunction.R
+++ b/R/formulaFunction.R
@@ -147,3 +147,80 @@ scaleFormulaFunction <- list(process =
   x
 })
 class(scaleFormulaFunction) <- "formulaFunction"
+
+
+#' 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).
+#'
+#' @name IFormulaFunction
+#'
+#' @param x A formulaComponentFunction object created from an I() term
+#'
+#'
+NULL
+
+#' @export
+IFormulaFunction <- list(process = 
+  function(x, defaultBracket, coefPrefix, sdPrefix, modelInfo, env, ...){
+
+  # Identify which interaction terms involve scale
+  trms <- splitInteractionTerms(x$lang)
+  has_I <- !sapply(trms, is.name)
+
+  bracks <- extractAllBrackets(x$lang)
+
+  for (i in 1:length(trms)){
+    if(!has_I[i]){
+      trm_idx <- safeDeparse(trms[[i]])
+      if(trm_idx %in% names(bracks)){
+        brack <- bracks[[trm_idx]]
+      } else {
+        brack <- defaultBracket
+      }
+      trms[[i]] <- str2lang(paste0(safeDeparse(trms[[i]]), brack))
+    } else {
+      interior <- trms[[i]][[2]]
+
+      # Extract any brackets
+      trm_idx <- all.vars(interior)
+      if(length(trm_idx) > 1) stop("Cannot handle more than one variable in I()", call.=FALSE)
+      if(trm_idx %in% names(bracks)){
+        brack <- bracks[[trm_idx]]
+      } else {
+        brack <- defaultBracket
+      }
+  
+      # New term
+      new_term <- safeDeparse(interior)
+      new_term <- gsub(" ", "", new_term)
+      new_term <- gsub("\\^|\\*|\\+|\\-|\\/", "_", new_term)
+
+      # Get constant
+      const <- modelInfo$constants[[trm_idx]]
+      if(is.null(const)) stop("Covariate inside I() is missing from constants", call.=FALSE)
+      if(!is.numeric(const)) stop("Covariate inside I() must be numeric", call.=FALSE)
+
+      # Evaluate expression
+      out <- eval(interior, env = modelInfo$constants)
+      add_const <- list(out)
+      names(add_const) <- new_term 
+
+      # Add the new constant to the object
+      if(is.null(x$constants)) x$constants <- list()
+      x$constants <- modifyList(x$constants, add_const)
+      trms[[i]] <- str2lang(paste0(new_term, brack))
+    }
+  }
+
+  # Updated interaction
+  x$lang <- str2lang(paste(sapply(trms, safeDeparse), collapse=":"))
+
+  # Update class
+  class(x)[1] <- "formulaComponentFixed"
+  
+  x
+})
+class(IFormulaFunction) <- "formulaFunction"
diff --git a/man/IFormulaFunction.Rd b/man/IFormulaFunction.Rd
new file mode 100644
index 0000000..367886d
--- /dev/null
+++ b/man/IFormulaFunction.Rd
@@ -0,0 +1,13 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/formulaFunction.R
+\name{IFormulaFunction}
+\alias{IFormulaFunction}
+\title{Function to handle I() in LINPRED}
+\arguments{
+\item{x}{A formulaComponentFunction object created from an I() term}
+}
+\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).
+}
diff --git a/tests/testthat/test_formulaFunctions.R b/tests/testthat/test_formulaFunctions.R
index 3f03d38..6e20f74 100644
--- a/tests/testthat/test_formulaFunctions.R
+++ b/tests/testthat/test_formulaFunctions.R
@@ -292,3 +292,127 @@ test_that("scale formula function works", {
 
   nimbleOptions(enableMacroComments=TRUE)
 })
+
+test_that("I() formula function works", {
+  nimbleOptions(enableMacroComments=FALSE)
+  set.seed(123)
+  constants <- list(x=rnorm(3), z=rnorm(3), a=matrix(rnorm(3), 3, 1), n=3)
+
+  code <- nimbleCode({
+    mu[1:n] <- LINPRED(~x + I(z^2))
+  })
+  mod <- nimbleModel(code, constants=constants)
+  
+  # Basic example
+  expect_equal(
+    mod$getCode(),
+    quote({
+    for (i_1 in 1:n) {
+        mu[i_1] <- beta_Intercept + beta_x * x[i_1] + beta_z_2 * 
+            z_2[i_1]
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_x ~ dnorm(0, sd = 1000)
+    beta_z_2 ~ dnorm(0, sd = 1000)
+    })
+  )
+  expect_equal(
+    mod$getConstants()$z_2,
+    constants$z^2
+  )
+
+  # Involved in interaction
+  code <- nimbleCode({
+    mu[1:n] <- LINPRED(~x:I(z^2))
+  })
+  mod <- nimbleModel(code, constants=constants)
+  expect_equal(
+    mod$getCode(),
+    quote({
+    for (i_1 in 1:n) {
+        mu[i_1] <- beta_Intercept + beta_x_z_2 * x[i_1] * z_2[i_1]
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_x_z_2 ~ dnorm(0, sd = 1000)
+    })
+  )
+  expect_equal(
+    mod$getConstants()$z_2,
+    constants$z^2
+  )
+
+  # Subtraction operation
+  code <- nimbleCode({
+    mu[1:n] <- LINPRED(~I(x-1))
+  })
+  mod <- nimbleModel(code, constants=constants)
+  
+  expect_equal(
+    mod$getCode(),
+    quote({
+    for (i_1 in 1:n) {
+        mu[i_1] <- beta_Intercept + beta_x_1 * x_1[i_1]
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_x_1 ~ dnorm(0, sd = 1000)
+    })
+  )
+  expect_equal(
+    mod$getConstants()$x_1,
+    constants$x - 1
+  )
+
+  # Indices
+  code <- nimbleCode({
+    mu[1:n] <- LINPRED(~x + I(a[1:n,1]^2))
+  })
+  mod <- nimbleModel(code, constants=constants)
+  
+  expect_equal(
+    mod$getCode(),
+    quote({
+    for (i_1 in 1:n) {
+        mu[i_1] <- beta_Intercept + beta_x * x[i_1] + beta_a_2 * 
+            a_2[i_1, 1]
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_x ~ dnorm(0, sd = 1000)
+    beta_a_2 ~ dnorm(0, sd = 1000)
+    })
+  )
+  expect_equal(
+    mod$getConstants()$a_2,
+    constants$a^2
+  )
+
+  code <- nimbleCode({
+    mu[1:n] <- LINPRED(~x + z:I(a[1:n,1]^2))
+  })
+  mod <- nimbleModel(code, constants=constants)
+  
+  expect_equal(
+    mod$getCode(),
+    quote({
+    for (i_1 in 1:n) {
+        mu[i_1] <- beta_Intercept + beta_x * x[i_1] + beta_z_a_2 * 
+            z[i_1] * a_2[i_1, 1]
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_x ~ dnorm(0, sd = 1000)
+    beta_z_a_2 ~ dnorm(0, sd = 1000)
+    })
+  )
+  expect_equal(
+    mod$getConstants()$a_2,
+    constants$a^2
+  )
+
+  # Two variables in I() not allowed
+  code <- quote(mu[1:n] <- LINPRED(~I(x+z)))
+  expect_error(
+    LINPRED$process(code, list(constants=constants), environment()),
+    "more than one"
+  )
+
+  nimbleOptions(enableMacroComments=TRUE)
+})

From 9885b8d304b1f293b268f47f043c7b9a81c5c0de Mon Sep 17 00:00:00 2001
From: Ken Kellner <ken@kenkellner.com>
Date: Thu, 16 Jan 2025 14:20:06 -0500
Subject: [PATCH 5/9] Add vignette section

---
 vignettes/nimbleMacros.Rmd      | 43 +++++++++++++++++++++++++++++++++
 vignettes/nimbleMacros.Rmd.orig | 25 +++++++++++++++++++
 2 files changed, 68 insertions(+)

diff --git a/vignettes/nimbleMacros.Rmd b/vignettes/nimbleMacros.Rmd
index a50db13..b8a890c 100644
--- a/vignettes/nimbleMacros.Rmd
+++ b/vignettes/nimbleMacros.Rmd
@@ -490,6 +490,49 @@ mod$getCode()
 ## }
 ```
 
+
+## Functions in formulas
+
+The `LINPRED` macro supports a handful of functions commonly included in formulas: `offset()`, `scale()`, and `I()`.
+For example, to build a linear predictor with the `Time` covariate scaled to have mean 0 and SD 1:
+
+
+``` r
+code <- nimbleCode({
+  mu[1:N] <- LINPRED(~scale(Time[1:N]), priorSpecs = NULL, link = log)
+})
+mod <- nimbleModel(code, chick_const)
+mod$getCode()
+```
+
+```
+## {
+##     "# LINPRED"
+##     "  ## nimbleMacros::FORLOOP"
+##     for (i_1 in 1:N) {
+##         log(mu[i_1]) <- beta_Intercept + beta_Time_scaled * Time_scaled[i_1]
+##     }
+##     "  ## ----"
+##     "# ----"
+## }
+```
+
+This adds a new constant `Time_scaled` to the constants:
+
+
+``` r
+head(mod$getConstants()$Time_scaled)
+```
+
+```
+## [1] -1.5858774 -1.2899493 -0.9940213 -0.6980932 -0.4021652 -0.1062371
+```
+
+`LINPRED` uses a modular system to handle formula functions.
+You can add support for additional functions in `LINPRED` formulas by writing a function with a specific name, class, and output.
+See the source code of the package, specifically the file `formulaFunction.R`, for examples of how to create these functions.
+
+
 ## 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`).
diff --git a/vignettes/nimbleMacros.Rmd.orig b/vignettes/nimbleMacros.Rmd.orig
index 34f2a22..ef3adc6 100644
--- a/vignettes/nimbleMacros.Rmd.orig
+++ b/vignettes/nimbleMacros.Rmd.orig
@@ -278,6 +278,31 @@ mod <- nimbleModel(code, chick_const)
 mod$getCode()
 ```
 
+
+## Functions in formulas
+
+The `LINPRED` macro supports a handful of functions commonly included in formulas: `offset()`, `scale()`, and `I()`.
+For example, to build a linear predictor with the `Time` covariate scaled to have mean 0 and SD 1:
+
+```{r}
+code <- nimbleCode({
+  mu[1:N] <- LINPRED(~scale(Time[1:N]), priorSpecs = NULL, link = log)
+})
+mod <- nimbleModel(code, chick_const)
+mod$getCode()
+```
+
+This adds a new constant `Time_scaled` to the constants:
+
+```{r}
+head(mod$getConstants()$Time_scaled)
+```
+
+`LINPRED` uses a modular system to handle formula functions.
+You can add support for additional functions in `LINPRED` formulas by writing a function with a specific name, class, and output.
+See the source code of the package, specifically the file `formulaFunction.R`, for examples of how to create these functions.
+
+
 ## 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`).

From 0ef8891adeb8820f392349ff9da06325b28561d2 Mon Sep 17 00:00:00 2001
From: Ken Kellner <ken@kenkellner.com>
Date: Sat, 18 Jan 2025 19:43:37 -0500
Subject: [PATCH 6/9] formula functions are now functions instead of lists

---
 R/LINPRED.R                  |  2 +-
 R/formulaFunction.R          | 77 +++++++++++++++++++++---------------
 man/IFormulaFunction.Rd      | 21 ++++++++++
 man/offsetFormulaFunction.Rd | 33 +++++++++++++++-
 man/scaleFormulaFunction.Rd  | 31 ++++++++++++++-
 5 files changed, 129 insertions(+), 35 deletions(-)

diff --git a/R/LINPRED.R b/R/LINPRED.R
index 70a77cc..514222d 100644
--- a/R/LINPRED.R
+++ b/R/LINPRED.R
@@ -364,7 +364,7 @@ processFormulaFunction <- function(x, defaultBracket, coefPrefix="beta_",
     processor <- get(cand, envir = env)
     if(inherits(processor, "formulaFunction")){
       processor_available <- TRUE
-      out <- processor$process(x, defaultBracket, coefPrefix, sdPrefix, modelInfo, env, ...)
+      out <- processor(x, defaultBracket, coefPrefix, sdPrefix, modelInfo, env, ...)
       if(!inherits(out, "formulaComponent")){
         stop("Processing function doesn't return formulaComponent object", call.=FALSE)
       }
diff --git a/R/formulaFunction.R b/R/formulaFunction.R
index a5d284e..3ceb2f9 100644
--- a/R/formulaFunction.R
+++ b/R/formulaFunction.R
@@ -1,18 +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.
+#' 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{formulaFunction}) 
+#' and return object class (\code{formulaComponent}).
 #'
-#' @name offsetFormulaFunction
+#' @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
 #'
-#' @param x A formulaComponentFunction object created from an offset() term
+#' @return An object of class \code{formulaComponent}.
 #'
-#'
-NULL
-
 #' @export
-offsetFormulaFunction <- list(process = 
-  function(x, defaultBracket, coefPrefix, sdPrefix, modelInfo, env, ...){
+offsetFormulaFunction <- function(x, defaultBracket, coefPrefix, sdPrefix, modelInfo, env, ...){
  
   # Get the code inside offset()
   interior <- x$lang[[2]]
@@ -42,7 +47,7 @@ offsetFormulaFunction <- list(process =
     add_const <- list(new_const)
     names(add_const) <- const_name
     if(is.null(x$constants)) x$constants <- list()
-    x$constants <- modifyList(x$constants, add_const)
+    x$constants <- utils::modifyList(x$constants, add_const)
     # Update the name of the term to put in the linear predictor
     interior <- str2lang(const_name)
   }
@@ -52,8 +57,8 @@ offsetFormulaFunction <- list(process =
   x$linPredCode <- code
   x$priorCode <- NULL   # no prior code for an offset
   x
-})
-class(offsetFormulaFunction) <- "formulaFunction"
+}
+class(offsetFormulaFunction) <- c(class(offsetFormulaFunction), "formulaFunction")
 
 # Extract all 'names' from a call
 # So qnorm(log(x)) would yield qnorm, log, x
@@ -71,18 +76,23 @@ get_all_names_recursive <- function(code){
 #' 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.
-#'
-#' @name scaleFormulaFunction
+#' 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{formulaFunction}) 
+#' and return object class (\code{formulaComponent}).
 #'
 #' @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}.
 #'
-NULL
-
 #' @export
-scaleFormulaFunction <- list(process = 
-  function(x, defaultBracket, coefPrefix, sdPrefix, modelInfo, env, ...){
+scaleFormulaFunction <- function(x, defaultBracket, coefPrefix, sdPrefix, modelInfo, env, ...){
 
   # Identify which interaction terms involve scale
   trms <- splitInteractionTerms(x$lang)
@@ -133,7 +143,7 @@ scaleFormulaFunction <- list(process =
 
       # Add the new constant to the object
       if(is.null(x$constants)) x$constants <- list()
-      x$constants <- modifyList(x$constants, add_const)
+      x$constants <- utils::modifyList(x$constants, add_const)
       trms[[i]] <- str2lang(paste0(new_term, brack))
     }
   }
@@ -145,8 +155,8 @@ scaleFormulaFunction <- list(process =
   class(x)[1] <- "formulaComponentFixed"
   
   x
-})
-class(scaleFormulaFunction) <- "formulaFunction"
+}
+class(scaleFormulaFunction) <- c(class(scaleFormulaFunction), "formulaFunction")
 
 
 #' Function to handle I() in LINPRED
@@ -154,17 +164,22 @@ class(scaleFormulaFunction) <- "formulaFunction"
 #' 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).
-#'
-#' @name IFormulaFunction
+#' 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{formulaFunction}) 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}.
 #'
-NULL
-
 #' @export
-IFormulaFunction <- list(process = 
-  function(x, defaultBracket, coefPrefix, sdPrefix, modelInfo, env, ...){
+IFormulaFunction <- function(x, defaultBracket, coefPrefix, sdPrefix, modelInfo, env, ...){
 
   # Identify which interaction terms involve scale
   trms <- splitInteractionTerms(x$lang)
@@ -204,13 +219,13 @@ IFormulaFunction <- list(process =
       if(!is.numeric(const)) stop("Covariate inside I() must be numeric", call.=FALSE)
 
       # Evaluate expression
-      out <- eval(interior, env = modelInfo$constants)
+      out <- eval(interior, envir = modelInfo$constants)
       add_const <- list(out)
       names(add_const) <- new_term 
 
       # Add the new constant to the object
       if(is.null(x$constants)) x$constants <- list()
-      x$constants <- modifyList(x$constants, add_const)
+      x$constants <- utils::modifyList(x$constants, add_const)
       trms[[i]] <- str2lang(paste0(new_term, brack))
     }
   }
@@ -222,5 +237,5 @@ IFormulaFunction <- list(process =
   class(x)[1] <- "formulaComponentFixed"
   
   x
-})
-class(IFormulaFunction) <- "formulaFunction"
+}
+class(IFormulaFunction) <- c(class(IFormulaFunction), "formulaFunction")
diff --git a/man/IFormulaFunction.Rd b/man/IFormulaFunction.Rd
index 367886d..cf6e387 100644
--- a/man/IFormulaFunction.Rd
+++ b/man/IFormulaFunction.Rd
@@ -3,11 +3,32 @@
 \name{IFormulaFunction}
 \alias{IFormulaFunction}
 \title{Function to handle I() in LINPRED}
+\usage{
+IFormulaFunction(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{formulaFunction}) and return object class (\code{formulaComponent}).
 }
diff --git a/man/offsetFormulaFunction.Rd b/man/offsetFormulaFunction.Rd
index 36055a6..bf9bed4 100644
--- a/man/offsetFormulaFunction.Rd
+++ b/man/offsetFormulaFunction.Rd
@@ -3,10 +3,39 @@
 \name{offsetFormulaFunction}
 \alias{offsetFormulaFunction}
 \title{Function to handle offset() in LINPRED}
+\usage{
+offsetFormulaFunction(
+  x,
+  defaultBracket,
+  coefPrefix,
+  sdPrefix,
+  modelInfo,
+  env,
+  ...
+)
+}
 \arguments{
-\item{x}{A formulaComponentFunction object created from an offset() term}
+\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.
+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{formulaFunction}) 
+and return object class (\code{formulaComponent}).
 }
diff --git a/man/scaleFormulaFunction.Rd b/man/scaleFormulaFunction.Rd
index b1db1ea..21ad705 100644
--- a/man/scaleFormulaFunction.Rd
+++ b/man/scaleFormulaFunction.Rd
@@ -3,10 +3,39 @@
 \name{scaleFormulaFunction}
 \alias{scaleFormulaFunction}
 \title{Function to handle scale() in LINPRED}
+\usage{
+scaleFormulaFunction(
+  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.
+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{formulaFunction}) 
+and return object class (\code{formulaComponent}).
 }

From 350d2d7149fae18373014f9136fcff1068f4bdc8 Mon Sep 17 00:00:00 2001
From: Ken Kellner <ken@kenkellner.com>
Date: Fri, 24 Jan 2025 14:45:09 -0500
Subject: [PATCH 7/9] Fix scale() on non-vectors

---
 R/formulaFunction.R                    |  7 +++++--
 man/scaleFormulaFunction.Rd            |  4 +++-
 tests/testthat/test_formulaFunctions.R | 26 +++++++++++++++++++++++++-
 3 files changed, 33 insertions(+), 4 deletions(-)

diff --git a/R/formulaFunction.R b/R/formulaFunction.R
index 3ceb2f9..a379213 100644
--- a/R/formulaFunction.R
+++ b/R/formulaFunction.R
@@ -79,7 +79,9 @@ get_all_names_recursive <- function(code){
 #' 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{formulaFunction}) 
-#' and return object class (\code{formulaComponent}).
+#' 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
@@ -133,7 +135,8 @@ scaleFormulaFunction <- function(x, defaultBracket, coefPrefix, sdPrefix, modelI
       if(!is.numeric(const)) stop("Covariate inside scale() must be numeric", call.=FALSE)
 
       # Make sure constant retains dimensions after being scaled
-      out <- as.numeric(scale(const))
+      # Note: this scales across the entire constant, not by column!
+      out <- as.numeric(scale(as.numeric(const)))
       if(!is.null(dim(const))){
         out <- array(out, dim=dim(const))
       }
diff --git a/man/scaleFormulaFunction.Rd b/man/scaleFormulaFunction.Rd
index 21ad705..1081b2c 100644
--- a/man/scaleFormulaFunction.Rd
+++ b/man/scaleFormulaFunction.Rd
@@ -37,5 +37,7 @@ 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{formulaFunction}) 
-and return object class (\code{formulaComponent}).
+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.
 }
diff --git a/tests/testthat/test_formulaFunctions.R b/tests/testthat/test_formulaFunctions.R
index 6e20f74..72d6913 100644
--- a/tests/testthat/test_formulaFunctions.R
+++ b/tests/testthat/test_formulaFunctions.R
@@ -111,7 +111,8 @@ test_that("scale formula function works", {
   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))
+                    x2=matrix(runif(3,5,10), 3, 1), z2=matrix(rnorm(3), 3, 1),
+                    x3 = matrix(runif(9,5,10), 3, 3))
 
   # Basic scale
   code <- nimbleCode({
@@ -141,6 +142,29 @@ test_that("scale formula function works", {
     constants$x
   )
 
+  # Scale with non-vectors
+  code <- nimbleCode({
+    mu[1:n,1:3] <- LINPRED(~scale(x3[1:n,1:3])) 
+  })
+  mod <- nimbleModel(code, constants=constants)
+  expect_equal(
+    mod$getCode(),
+    quote({
+    for (i_1 in 1:n) {
+        for (i_2 in 1:3) {
+            mu[i_1, i_2] <- beta_Intercept + beta_x3_scaled * 
+                x3_scaled[i_1, i_2]
+        }
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_x3_scaled ~ dnorm(0, sd = 1000)
+    })
+  )
+  x3_scale <- mod$getConstants()$x3_scale
+  expect_equal(dim(x3_scale), dim(constants$x3))
+  expect_equal(mean(x3_scale), 0, tol=1e-6)
+  expect_equal(sd(x3_scale), 1, tol=1e-6)
+
   # Interaction with non-scaled term
   code <- nimbleCode({
     mu[1:n] <- LINPRED(~scale(x):z) 

From 7692423588ffa3365fec7a86a0aa78a2427e87a7 Mon Sep 17 00:00:00 2001
From: Ken Kellner <ken@kenkellner.com>
Date: Fri, 24 Jan 2025 15:21:35 -0500
Subject: [PATCH 8/9] Add support for log() in formulas

---
 NAMESPACE                              |   1 +
 R/formulaFunction.R                    |  82 +++++++++++
 man/logFormulaFunction.Rd              |  41 ++++++
 tests/testthat/test_formulaFunctions.R | 179 +++++++++++++++++++++++++
 vignettes/nimbleMacros.Rmd             |  21 +--
 vignettes/nimbleMacros.Rmd.orig        |  10 +-
 6 files changed, 319 insertions(+), 15 deletions(-)
 create mode 100644 man/logFormulaFunction.Rd

diff --git a/NAMESPACE b/NAMESPACE
index f85a273..34c7bda 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -5,6 +5,7 @@ export(IFormulaFunction)
 export(LINPRED)
 export(LINPRED_PRIORS)
 export(LM)
+export(logFormulaFunction)
 export(matchPrior)
 export(offsetFormulaFunction)
 export(scaleFormulaFunction)
diff --git a/R/formulaFunction.R b/R/formulaFunction.R
index a379213..2885856 100644
--- a/R/formulaFunction.R
+++ b/R/formulaFunction.R
@@ -242,3 +242,85 @@ IFormulaFunction <- function(x, defaultBracket, coefPrefix, sdPrefix, modelInfo,
   x
 }
 class(IFormulaFunction) <- c(class(IFormulaFunction), "formulaFunction")
+
+
+#' 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{formulaFunction}) 
+#' 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}.
+#'
+#' @export
+logFormulaFunction <- function(x, defaultBracket, coefPrefix, sdPrefix, modelInfo, env, ...){
+
+  # Identify which interaction terms involve scale
+  trms <- splitInteractionTerms(x$lang)
+  has_log <- !sapply(trms, is.name)
+
+  bracks <- extractAllBrackets(x$lang)
+
+  for (i in 1:length(trms)){
+    if(!has_log[i]){
+      trm_idx <- safeDeparse(trms[[i]])
+      if(trm_idx %in% names(bracks)){
+        brack <- bracks[[trm_idx]]
+      } else {
+        brack <- defaultBracket
+      }
+      trms[[i]] <- str2lang(paste0(safeDeparse(trms[[i]]), brack))
+    } else {
+      interior <- trms[[i]][[2]]
+      if(!is.name(interior)) stop("Can't handle expression inside log()", call.=FALSE)
+
+      # Extract any brackets
+      trm_idx <- safeDeparse(interior)
+      if(trm_idx %in% names(bracks)){
+        brack <- bracks[[trm_idx]]
+      } else {
+        brack <- defaultBracket
+      }
+
+      # Remove brackets from code
+      interior <- removeSquareBrackets(interior)
+  
+      # New term
+      new_term <- paste0(safeDeparse(interior), "_log")
+
+      # Get constant
+      const <- modelInfo$constants[[safeDeparse(interior)]]
+      if(is.null(const)) stop("Covariate inside log() is missing from constants", call.=FALSE)
+      if(!is.numeric(const)) stop("Covariate inside log() must be numeric", call.=FALSE)
+
+      out <- log(const)
+
+      add_const <- list(out)
+      names(add_const) <- new_term 
+
+      # Add the new constant to the object
+      if(is.null(x$constants)) x$constants <- list()
+      x$constants <- utils::modifyList(x$constants, add_const)
+      trms[[i]] <- str2lang(paste0(new_term, brack))
+    }
+  }
+
+  # Updated interaction
+  x$lang <- str2lang(paste(sapply(trms, safeDeparse), collapse=":"))
+
+  # Update class
+  class(x)[1] <- "formulaComponentFixed"
+  
+  x
+}
+class(logFormulaFunction) <- c(class(logFormulaFunction), "formulaFunction")
diff --git a/man/logFormulaFunction.Rd b/man/logFormulaFunction.Rd
new file mode 100644
index 0000000..fca4426
--- /dev/null
+++ b/man/logFormulaFunction.Rd
@@ -0,0 +1,41 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/formulaFunction.R
+\name{logFormulaFunction}
+\alias{logFormulaFunction}
+\title{Function to handle log() in LINPRED}
+\usage{
+logFormulaFunction(
+  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{formulaFunction}) 
+and return object class (\code{formulaComponent}).
+}
diff --git a/tests/testthat/test_formulaFunctions.R b/tests/testthat/test_formulaFunctions.R
index 72d6913..6682162 100644
--- a/tests/testthat/test_formulaFunctions.R
+++ b/tests/testthat/test_formulaFunctions.R
@@ -440,3 +440,182 @@ test_that("I() formula function works", {
 
   nimbleOptions(enableMacroComments=TRUE)
 })
+
+test_that("log formula function works", {
+  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))
+
+  # Basic log
+  code <- nimbleCode({
+    mu[1:n] <- LINPRED(~log(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_log * x_log[i_1] + 
+            beta_z * z[i_1]
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_x_log ~ dnorm(0, sd = 1000)
+    beta_z ~ dnorm(0, sd = 1000)
+    })
+  )
+  expect_equivalent(
+    mod$getConstants()$x_log,
+    log(constants$x)
+  )
+
+  # log with non-vectors
+  code <- nimbleCode({
+    mu[1:n,1:3] <- LINPRED(~log(x3[1:n,1:3])) 
+  })
+  mod <- nimbleModel(code, constants=constants)
+  expect_equal(
+    mod$getCode(),
+    quote({
+    for (i_1 in 1:n) {
+        for (i_2 in 1:3) {
+            mu[i_1, i_2] <- beta_Intercept + beta_x3_log * 
+                x3_log[i_1, i_2]
+        }
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_x3_log ~ dnorm(0, sd = 1000)
+    })
+  )
+  expect_equal(mod$getConstants()$x3_log, log(constants$x3))
+
+  # Interaction with non-scaled term
+  code <- nimbleCode({
+    mu[1:n] <- LINPRED(~log(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_log_z * x_log[i_1] * 
+            z[i_1]
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_x_log_z ~ dnorm(0, sd = 1000)
+    })
+  )
+  expect_equivalent(
+    mod$getConstants()$x_log,
+    log(constants$x)
+  )
+
+  # Non-log term comes first in interaction
+  code <- nimbleCode({
+    mu[1:n] <- LINPRED(~z:log(x)) 
+  })
+  mod <- nimbleModel(code, constants=constants)
+  expect_equal(
+    mod$getCode(),
+    quote({
+    for (i_1 in 1:n) {
+        mu[i_1] <- beta_Intercept + beta_z_x_log * z[i_1] * 
+            x_log[i_1]
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_z_x_log ~ dnorm(0, sd = 1000)
+    })
+  )
+  expect_equivalent(
+    mod$getConstants()$x_log,
+    log(constants$x)
+  )
+
+  # Interaction of two logged terms
+  code <- nimbleCode({
+    mu[1:n] <- LINPRED(~log(x):log(z)) 
+  })
+  expect_warning(mod <- nimbleModel(code, constants=constants), "NaN")
+  expect_equal(
+    mod$getCode(),
+    quote({
+    for (i_1 in 1:n) {
+        mu[i_1] <- beta_Intercept + beta_x_log_z_log * 
+            x_log[i_1] * z_log[i_1]
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_x_log_z_log ~ dnorm(0, sd = 1000)
+    })
+  )
+  expect_equivalent(
+    mod$getConstants()$x_log,
+    log(constants$x)
+  )
+  expect_equivalent(
+    mod$getConstants()$z_log,
+    expect_warning(log(constants$z))
+  )
+
+  # log with brackets
+  code <- nimbleCode({
+    mu[1:n] <- LINPRED(~log(x2[1:n,1]):z)
+  })
+  mod <- nimbleModel(code, constants=constants)
+  expect_equal(
+    mod$getCode(),
+    quote({
+    for (i_1 in 1:n) {
+        mu[i_1] <- beta_Intercept + beta_x2_log_z * x2_log[i_1, 
+            1] * z[i_1]
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_x2_log_z ~ dnorm(0, sd = 1000)
+  })
+  )
+  expect_equivalent(
+    mod$getConstants()$x2_log,
+    log(constants$x2)
+  )
+
+  # Scale with brackets on the other interaction term
+  code <- nimbleCode({
+    mu[1:n] <- LINPRED(~log(x2[1:n,1]):z2[1:n,1])
+  })
+  mod <- nimbleModel(code, constants=constants)
+  expect_equal(
+    mod$getCode(),
+    quote({
+    for (i_1 in 1:n) {
+        mu[i_1] <- beta_Intercept + beta_x2_log_z2 * x2_log[i_1, 
+            1] * z2[i_1, 1]
+    }
+    beta_Intercept ~ dnorm(0, sd = 1000)
+    beta_x2_log_z2 ~ dnorm(0, sd = 1000)
+  })
+  )
+
+  # Error if covariate is not in constants
+  code <- quote(mu[1:n] <- LINPRED(~log(test)))
+  expect_error(
+    LINPRED$process(code, list(constants=constants), environment()),
+    "Covariate inside"
+  )
+
+  # Error if expression inside scale
+  code <- quote(mu[1:n] <- LINPRED(~log(x*x)))
+  expect_error(
+    LINPRED$process(code, list(constants=constants), environment()),
+    "expression"
+  )
+
+  # Error if multiple types of functions in an interaction
+  code <- quote(mu[1:n] <- LINPRED(~log(x):scale(z)))
+  expect_error(
+    LINPRED$process(code, list(constants=constants), environment()),
+    "multiple different formula functions"
+  )
+
+  nimbleOptions(enableMacroComments=TRUE)
+})
diff --git a/vignettes/nimbleMacros.Rmd b/vignettes/nimbleMacros.Rmd
index b8a890c..f9efd24 100644
--- a/vignettes/nimbleMacros.Rmd
+++ b/vignettes/nimbleMacros.Rmd
@@ -78,7 +78,7 @@ To fit this model with the `LM` macro, we'll start by organizing the data and co
 
 
 ``` r
-const <- list(mpg = as.numeric(scale(mtcars$mpg)), am = mtcars$am)
+const <- mtcars[c("am", "mpg")]
 ```
 
 Next we'll specify the prior distributions we want to use for the intercept and the slope coefficient(s).
@@ -93,14 +93,15 @@ pr <- setPriors(intercept=quote(dunif(-10,10)),
 ```
 
 Finally, we specify the NIMBLE code, which is a single call to the `LM` macro that should look familiar to our call to `glm` above.
+
 Two things to note:
-1. We don't include a call to `scale` in the formula (`LM` does not support this). Instead we standardized the `mpg` values when we put them into our list `const`.
+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.
 
 
 ``` r
 code <- nimbleCode({
-  LM(am ~ mpg, family = binomial, priorSpecs=pr)
+  LM(am ~ scale(mpg), family = binomial, priorSpecs=pr)
 })
 ```
 
@@ -131,15 +132,15 @@ summary(samples)
 ## 1. Empirical mean and standard deviation for each variable,
 ##    plus standard error of the mean:
 ## 
-##                   Mean     SD Naive SE Time-series SE
-## beta_Intercept -0.4651 0.4693  0.01484        0.03274
-## beta_mpg        2.0391 0.7063  0.02234        0.04881
+##                    Mean     SD Naive SE Time-series SE
+## beta_Intercept  -0.4651 0.4693  0.01484        0.03274
+## beta_mpg_scaled  2.0391 0.7063  0.02234        0.04881
 ## 
 ## 2. Quantiles for each variable:
 ## 
-##                   2.5%     25%     50%    75%  97.5%
-## beta_Intercept -1.4296 -0.7595 -0.4548 -0.134 0.4101
-## beta_mpg        0.7872  1.5560  2.0039  2.464 3.5055
+##                    2.5%     25%     50%    75%  97.5%
+## beta_Intercept  -1.4296 -0.7595 -0.4548 -0.134 0.4101
+## beta_mpg_scaled  0.7872  1.5560  2.0039  2.464 3.5055
 ```
 
 ## Example linear mixed model using `ChickWeights`
@@ -493,7 +494,7 @@ mod$getCode()
 
 ## Functions in formulas
 
-The `LINPRED` macro supports a handful of functions commonly included in formulas: `offset()`, `scale()`, and `I()`.
+The `LINPRED` macro supports a handful of functions commonly included in formulas: `offset()`, `scale()`, `log()`, and `I()`.
 For example, to build a linear predictor with the `Time` covariate scaled to have mean 0 and SD 1:
 
 
diff --git a/vignettes/nimbleMacros.Rmd.orig b/vignettes/nimbleMacros.Rmd.orig
index ef3adc6..1b8728e 100644
--- a/vignettes/nimbleMacros.Rmd.orig
+++ b/vignettes/nimbleMacros.Rmd.orig
@@ -58,7 +58,7 @@ summary(rmod)
 To fit this model with the `LM` macro, we'll start by organizing the data and constants into a list.
 
 ```{r}
-const <- list(mpg = as.numeric(scale(mtcars$mpg)), am = mtcars$am)
+const <- mtcars[c("am", "mpg")]
 ```
 
 Next we'll specify the prior distributions we want to use for the intercept and the slope coefficient(s).
@@ -72,20 +72,20 @@ pr <- setPriors(intercept=quote(dunif(-10,10)),
 ```
 
 Finally, we specify the NIMBLE code, which is a single call to the `LM` macro that should look familiar to our call to `glm` above.
+
 Two things to note:
-1. We don't include a call to `scale` in the formula (`LM` does not support this). Instead we standardized the `mpg` values when we put them into our list `const`.
+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.
 
 ```{r}
 code <- nimbleCode({
-  LM(am ~ mpg, family = binomial, priorSpecs=pr)
+  LM(am ~ scale(mpg), family = binomial, priorSpecs=pr)
 })
 ```
 
 Finally, we'll fit the model which yields similar results to `glm`.
 
 ```{r}
-
 samples <- nimbleMCMC(code, constants=const, nchain=1, niter=3000, nburnin=2000,
                       samplesAsCodaMCMC = TRUE)
 summary(samples)
@@ -281,7 +281,7 @@ mod$getCode()
 
 ## Functions in formulas
 
-The `LINPRED` macro supports a handful of functions commonly included in formulas: `offset()`, `scale()`, and `I()`.
+The `LINPRED` macro supports a handful of functions commonly included in formulas: `offset()`, `scale()`, `log()`, and `I()`.
 For example, to build a linear predictor with the `Time` covariate scaled to have mean 0 and SD 1:
 
 ```{r}

From 112739612f45499c8978ec61c6961eae5188387e Mon Sep 17 00:00:00 2001
From: Ken Kellner <ken@kenkellner.com>
Date: Fri, 24 Jan 2025 15:22:08 -0500
Subject: [PATCH 9/9] Bump version

---
 DESCRIPTION | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 6ccb48a..26e1489 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: nimbleMacros
-Version: 0.0.12
-Date: 2025-01-16
+Version: 0.0.13
+Date: 2025-01-24
 Title: Macros for 'nimble' Code
 Authors@R: person("Ken", "Kellner", email="contact@kenkellner.com", 
                   role=c("cre","aut"))