Creates a JAGS formula syntax, prepares data input, and
returns modified prior list for further processing in the JAGS_fit
function.
JAGS_formula(
formula,
parameter,
data,
prior_list,
formula_scale = NULL,
prior_random = NULL,
random_effects_compile = NULL
)formula specifying the right hand side of the assignment (the
left hand side is ignored), or a BayesTools_random_effects object returned
by random_effects_formula(). If the formula suppresses the intercept with
0 + or -1, it is converted to include an intercept with a
spike(0) prior.
The formula can also have a "log(intercept)" attribute set to TRUE
to generate syntax of the form log(intercept) + sum(beta_i * x_i), which
is useful for parameters that must be positive (e.g., standard deviation).
In that case, the intercept prior must have strictly positive support.
valid unindexed JAGS node name of the parameter to be created with the formula
data.frame containing predictors included in the formula. A
logical predictor is coded as a two-level indicator; an interaction of a
logical predictor with a factor in which model.matrix() would code it
by indicators (e.g., flag:g without g) is unavailable, in a
fixed formula and in a random-effect term alike,
unless formula_scale standardizes the predictor, which makes it
numeric. Otherwise convert the predictor to a factor or to numeric 0/1.
named list of prior distribution of parameters specified within
the formula. The name "intercept" is reserved for the intercept:
no other predictor or prior name may contain the substring "intercept".
When suppressing the intercept in the formula, an
"intercept" prior can be explicitly specified; otherwise,
prior("spike", list(0)) is
automatically added, or prior("spike", list(1)) when the formula uses
the "log(intercept)" attribute. The list can also include two special entries:
"__default_continuous"A prior to use for any continuous predictors (including the intercept) that are not explicitly specified in the prior list. This can also be a zero-argument function returning a prior object; in that case it is evaluated only if a missing continuous term needs the default.
"__default_factor"A prior to use for any factor predictors (including interactions involving factors) that are not explicitly specified in the prior list. This can also be a zero-argument function returning a prior object; in that case it is evaluated only if a missing factor term needs the default.
These default priors allow for more concise specification when many predictors
share the same prior distribution. For continuous formula terms,
prior_none() is canonicalized to a point prior at zero.
The prior of a formula term can carry a "multiply_by" attribute (a
number or the name of another parameter) that multiplies the term's
contribution to the linear predictor; the intercept prior cannot be scaled
and a "multiply_by" attribute on it is an error. For a mixture or
spike-and-slab prior, set the attribute on that prior itself; a component
carrying one is an error.
named list specifying whether to standardize continuous predictors.
If NULL (default), no standardization is applied. If a named list is provided,
continuous predictors with TRUE values will be standardized (mean-centered and
scaled by standard deviation). The intercept is never standardized. Scaled
predictors are centered in every fixed and random term, including terms
without a free intercept, so the model is specified on the standardized
scale; see Details.
optional prior_random() object defining random-effect
standard-deviation, covariance, monitoring, and prediction policies. Required
when formula contains random effects.
optional random_effects_compile() object
specifying which resolved random-effect blocks should be compiled as sampled
random effects and which should be compiled as structural marginalized
blocks.
JAGS_formula returns a list containing the formula JAGS syntax,
JAGS data object, modified prior_list, and (if standardization was applied) a
formula_scale list with standardization information for back-transformation.
When a formula suppresses the intercept with 0 + or
-1, the function adds an intercept back to the compiled formula and
includes a point prior that contributes zero on the formula scale: spike(0)
ordinarily and spike(1) for a log-transformed intercept.
Scaling with formula_scale replaces a continuous predictor \(x\)
by \((x - m) / s\) in every term, where \(m\) and \(s\) are its mean
and standard deviation in the fitting data. Centering also applies to terms
without a free intercept, which is intended: the model is specified on the
standardized scale. For example, ~ 0 + x with a scaled x fits
\(\mu = b (x - m) / s\), a line through \((m, 0)\) rather than through
the origin; on the original scale its slope is \(b / s\) and its implied
intercept is \(-b m / s\). Original-scale summaries
(transform_scaled = TRUE) report that implied intercept although the
standardized intercept is fixed at zero. Likewise, a random-effect block
without a random intercept, such as (0 + x | g), fits group
contributions \(u (x - m) / s\). It implies an original-scale random
intercept \(-u m / s\) (SD \(\tau |m| / s\) for a scaled-slope SD
\(\tau\)), perfectly correlated with the original-scale slope
\(u / s\); random-effect summaries report the block's own coefficients,
such as sd(x) \(= \tau / s\), and not this implied intercept.
Leave a predictor unscaled when a line or group slope through the origin of
\(x\) is intended.
Independence of random-effect coefficients is likewise specified on the
centered scale. An independent block such as (1 + x || g) (or a
diag() block) with a scaled x has independent
group coefficients \(u_0\) and \(u_1\) with SDs \(\tau_0\) and
\(\tau_1\). On the original scale, the group intercept
\(u_0 - u_1 m / s\) and slope \(u_1 / s\) are correlated with
$$\rho = -\frac{\tau_1 m / s}{\sqrt{\tau_0^2 + (\tau_1 m / s)^2}},$$
for example \(\rho \approx -0.93\) for equal SDs and
\(m / s = 2.5\). The original-scale summaries of such a block report its
SDs, \(\sqrt{\tau_0^2 + (\tau_1 m / s)^2}\)
and \(\tau_1 / s\), but no
correlation, because the block has no correlation parameter. Use a
correlated block, such as (1 + x | g), to estimate the correlation,
or leave x unscaled for independence on the original scale.
A homogeneous-SD block (id(), random(covariance = "id"), or
diag(..., hom = TRUE)) shares one SD among its columns, which cannot
carry the different original-scale SDs of the columns. A homogeneous block
with more than one column that contains a column of a scaled predictor is
therefore refused by JAGS_formula(), whether its shared SD is sampled
or fixed by a point prior: use a heterogeneous block, or leave the predictor
unscaled. Single-column homogeneous blocks, such as
id(0 + x | g), are not affected.
Factor contrasts are owned by prior_factor() and stored formula metadata,
not inferred from whether an intercept is present. Thus ~ 0 + group with
an independent factor prior means a structural zero intercept plus one
coefficient per group level. With treatment or mean-difference contrasts,
the same no-intercept expression preserves that selected basis.
An interaction without one of its lower-order terms, such as g:x
in ~ g + g:x, codes that factor by level indicators, as
stats::model.matrix() does, and has one coefficient per level
(e.g., one slope per level of g). Such a term does not use the
contrast of that factor, so its prior may use a different contrast than the
factor's terms that code it by its contrast (e.g., a mean-difference
g with an independent g:x). When no term codes the factor by
its contrast (e.g., ~ g:x + g:z), its indicator-coded terms must
agree on the contrast instead. Mean-difference and orthonormal priors,
other than a point mass at zero, are unavailable for such a term; use a
treatment or independent prior for per-level coefficients, or add the
lower-order term (here ~ g * x) to keep the contrast.
When using default priors ("__default_continuous" or "__default_factor"),
explicitly specified priors for individual terms take precedence over the defaults.
The defaults are only applied to terms that are not already in the prior list.
Formula random effects require prior_random. Random-effect SD priors in
prior_list using "term|group" names are no longer supported.
Continuous fixed-effect terms must expand to one design-matrix column.
Matrix-valued continuous predictors are not currently supported.
Fixed formulas support literal data-column names and standard formula
operators, including interaction expansion such as (x + z)^2 with
literal whole degrees from 2 through .Machine$integer.max. This is
formula expansion, not arithmetic powers; I(x^2) remains unsupported.
Degree 0 or 1 and symbolic degrees retain the base-R formula refusal.
Dot expansion, offset(), inline transformations, and
arbitrary calls are rejected; create explicit data columns for transformed
predictors. The BayesTools expression(...) facility accepts finite
numeric constants, formula or model data, JAGS-style i indexing,
sampled scalar or one-dimensional indexed parameters, arithmetic operators,
and abs(), exp(), log(), or sqrt(). Expressions
retain literal JAGS indexing: users must write x or x[i] as
required by the JAGS data shape. Their parsed syntax and dependencies are
persisted and replayed draw by draw during prediction and marginal-
likelihood reconstruction. Parameter dependencies must be owned by a prior
or declared through add_parameters; opaque deterministic nodes and
formula-output dependencies are rejected by JAGS_fit() because their
defining JAGS graph is unavailable during replay.
Random-effect predictors likewise support literal data-column names and
formula operators, but not inline transformations or arbitrary calls.
Create transformed random slopes as explicit data columns. Grouping terms
support variables, : interactions, and / nesting.
Random-effect structures have two left-side grammars. id(), diag(), and
us() / un() use a coefficient formula, so 1, 0, and -1 control the
random intercept and slopes or interactions generate coefficient columns.
Plain bars default to us() and double bars to diag(). A factor's concrete
random basis is inherited from resolved contrast metadata by default and can
be set for one block with random_block(contrasts = ...).
cs() / hcs(), ar1() / ar() / har(), and car() instead use a
structure-owned index specification. They reject intercept controls and
block contrast overrides. See random_effects_formula() and
prior_random() for their index and covariance contracts.
Categorical predictor and grouping levels must not contain BayesTools'
reserved internal tokens, such as __xXx__.
The predictor name intercept is reserved for the formula intercept.
Refusals that callers may need to tell apart are errors of class
BayesTools_formula_transform_unavailable (parents error and
condition); their reason field names the cause:
"single_level_factor": a factor with a single level in the data
has no contrasts to code (the field predictor names it). The index
of a cs(), hcs(), ar1(), or har() block is coded
by one indicator per level and may have a single level.
"logical_factor_interaction": a logical predictor in an
interaction with a factor that model.matrix() would code by
indicators, in a fixed formula or a random-effect term (see the
data argument; the fields predictors and terms name
them).
"unsupported_random_factor_contrast": a factor used only in a
random-effect block that carries a contrasts matrix or an unsupported
contrast name (the field predictor names it; set the coding with
random_block(contrasts = )).
"homogeneous_scaled_block": a homogeneous-SD random block with
a standardized slope (see Details; the fields parameter and
coordinates name the formula parameter and the shared SD).
"basis_dependent_factor_prior": the prior of an interaction
term of two or more factors, of which at least one is coded by a
mean-difference or orthonormal contrast, depends on the basis of that
contrast and would not treat the factor levels exchangeably. Only
mean-difference and orthonormal priors, and treatment or independent
priors built from zero-mean untruncated normal and zero spike priors,
are accepted on such a term; when every coded factor has two levels,
treatment and independent priors symmetric about zero are accepted as
well. Other treatment and independent priors and every ordered prior are
refused (the fields terms and predictors name the term and
the coded factors; a refused ordered prior also has class
BayesTools_ordered_unavailable).
Other invalid inputs stop with ordinary errors.
# simulate data
set.seed(1)
df <- data.frame(
y = rnorm(60),
x_cont = rnorm(60),
x_bin = rbinom(60, 1, .5),
x_fac3 = factor(rep(c("A", "B", "C"), 20), levels = c("A", "B", "C")),
x_fac4 = factor(rep(c("A", "B", "C", "D"), 15), levels = c("A", "B", "C", "D"))
)
# specify priors with intercept
prior_list <- list(
"intercept" = prior("normal", list(0, 1)),
"x_cont" = prior("normal", list(0, .5)),
"x_fac3" = prior_factor("normal", list(0, 1), contrast = "treatment"),
"x_fac4" = prior_factor("mnormal", list(0, 1), contrast = "orthonormal"),
"x_fac3:x_fac4" = prior_factor("mnormal", list(0, .5), contrast = "orthonormal")
)
# create the formula object
formula_obj <- JAGS_formula(
formula = ~ x_cont + x_fac3 * x_fac4,
parameter = "mu", data = df, prior_list = prior_list)
# using -1 notation (automatically adds spike(0) intercept)
prior_list_no_intercept <- list(
"x_fac3" = prior_factor("normal", list(0, 1), contrast = "treatment")
)
formula_no_intercept <- JAGS_formula(
formula = ~ x_fac3 - 1,
parameter = "mu", data = df, prior_list = prior_list_no_intercept)
# Equivalent to specifying intercept = prior("spike", list(0))
# using default priors for simpler specification
prior_list_defaults <- list(
"__default_continuous" = prior("normal", list(0, 1)),
"__default_factor" = prior_factor("normal", list(0, 0.5), contrast = "treatment")
)
formula_defaults <- JAGS_formula(
formula = ~ x_cont + x_fac3,
parameter = "mu", data = df, prior_list = prior_list_defaults)
# intercept and x_cont get the default continuous prior
# x_fac3 gets the default factor prior