diff --git a/DESCRIPTION b/DESCRIPTION
index a8f9645..26e1489 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: nimbleMacros
-Version: 0.0.12
-Date: 2025-01-14
+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"))
diff --git a/NAMESPACE b/NAMESPACE
index 7bf36ec..34c7bda 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,10 +1,14 @@
 # Generated by roxygen2: do not edit by hand
 
 export(FORLOOP)
+export(IFormulaFunction)
 export(LINPRED)
 export(LINPRED_PRIORS)
 export(LM)
+export(logFormulaFunction)
 export(matchPrior)
+export(offsetFormulaFunction)
+export(scaleFormulaFunction)
 export(setPriors)
 export(uppertri_mult_diag)
 importFrom(nimble,nimMatrix)
diff --git a/R/LINPRED.R b/R/LINPRED.R
index ff603ed..514222d 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,59 @@ 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)
+
+  # 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
+  if(exists(cand, envir = env)){
+    processor <- get(cand, envir = env)
+    if(inherits(processor, "formulaFunction")){
+      processor_available <- TRUE
+      out <- processor(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
+}
+
+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
 # Returns the updated component
diff --git a/R/formulaFunction.R b/R/formulaFunction.R
new file mode 100644
index 0000000..2885856
--- /dev/null
+++ b/R/formulaFunction.R
@@ -0,0 +1,326 @@
+#' 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{formulaFunction}) 
+#' 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}.
+#'
+#' @export
+offsetFormulaFunction <- 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 <- utils::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) <- c(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
+}
+
+#' 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{formulaFunction}) 
+#' 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}.
+#'
+#' @export
+scaleFormulaFunction <- function(x, defaultBracket, coefPrefix, sdPrefix, modelInfo, env, ...){
+
+  # Identify which interaction terms involve scale
+  trms <- splitInteractionTerms(x$lang)
+  has_scale <- !sapply(trms, is.name)
+
+  bracks <- extractAllBrackets(x$lang)
+
+  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
+      # 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))
+      }
+
+      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(scaleFormulaFunction) <- c(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).
+#' 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}.
+#'
+#' @export
+IFormulaFunction <- 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, 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 <- 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(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/IFormulaFunction.Rd b/man/IFormulaFunction.Rd
new file mode 100644
index 0000000..cf6e387
--- /dev/null
+++ b/man/IFormulaFunction.Rd
@@ -0,0 +1,34 @@
+% 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}
+\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/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/man/offsetFormulaFunction.Rd b/man/offsetFormulaFunction.Rd
new file mode 100644
index 0000000..bf9bed4
--- /dev/null
+++ b/man/offsetFormulaFunction.Rd
@@ -0,0 +1,41 @@
+% 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}
+\usage{
+offsetFormulaFunction(
+  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{formulaFunction}) 
+and return object class (\code{formulaComponent}).
+}
diff --git a/man/scaleFormulaFunction.Rd b/man/scaleFormulaFunction.Rd
new file mode 100644
index 0000000..1081b2c
--- /dev/null
+++ b/man/scaleFormulaFunction.Rd
@@ -0,0 +1,43 @@
+% 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}
+\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. 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}). 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_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..6682162
--- /dev/null
+++ b/tests/testthat/test_formulaFunctions.R
@@ -0,0 +1,621 @@
+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)
+})
+
+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),
+                    x3 = matrix(runif(9,5,10), 3, 3))
+
+  # 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
+  )
+
+  # 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) 
+  })
+  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)
+})
+
+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)
+})
+
+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 a50db13..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`
@@ -490,6 +491,49 @@ mod$getCode()
 ## }
 ```
 
+
+## Functions in formulas
+
+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
+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..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)
@@ -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()`, `log()`, 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`).