Give a predictor a varying coefficient

Marks a predictor inside a bartisan() formula as one whose coefficient is a function of the other predictors rather than a constant, so that the coefficient gets a forest and a prior of its own. vc() is a marker rather than a function: it is read out of the formula it appears in and throws an error if it is called on its own.

Usage

vc(x, modifiers = NULL, center = "auto")

Arguments

x: the predictor whose coefficient varies, named as a bare variable rather than as an expression. A numeric variable gets one forest; a factor gets one per level, coded symmetrically (i.e., with no level held out as a reference) as multinomial() codes its predictors.
modifiers: a one-sided formula naming the variables this coefficient's forest may split on. Default is NULL to allow every predictor in the model except x itself. A . in this formula stands for those same predictors, so ~ . is the default written out and ~ . + w or ~ . - z adds to or removes from it; see Details. A variable the model formula's fixed part does not carry may be named here, and it then modifies the coefficient without entering the control function. Naming something that is not a column of data at all is an error rather than a silent restriction.
center: the value of x at which the control function is read, given as either a string or a number. Allowable options include "auto" (the default), "mean", "zero", "mid", and "estimate". "mean" centers x, "zero" leaves it alone, "mid" uses the midpoint of its range, a number uses that number, and "auto" picks between "zero" and "mean" by the covariate. For a factor, center is "mean" or the name of a level to report against. "estimate" draws the coding rather than fixing it, and needs a covariate with between two and twenty distinct values. See Details.

Value

Nothing; vc() is never evaluated, and calling it directly is an error.

Details

The model is

$$g(\mu_i) = f_0(Z_i) + \sum_j (X_{ij} - c_j) f_j(Z_i)$$

with a forest for the control function $f_0$ and one for each varying coefficient $f_j$. Every forest is fitted at once, so the coefficient has a prior of its own rather than being whatever difference a single forest with x among its predictors happens to produce.

Setting `modifiers`

By default a coefficient may vary with every predictor in the model except x itself, and the modifiers argument narrows that.

Within modifiers, . stands for the predictors of the control function rather than, as in a model formula, for every column of the data. That makes vc(z, ~ .) the same model as vc(z), and vc(z, ~ . - x1) the default with one modifier taken away, which is the short way to leave a coefficient free to vary with everything but one thing. vc(z, ~ . + w) adds one, and is how a variable the fixed part leaves out (below) is given to a coefficient without listing the rest. This is the sense . has in update.formula(), not the one it has in lm().

It can also widen it. A variable named here that the fixed part leaves out is given a column of its own and reaches this coefficient's forest and no other, so the coefficient varies with it while the control function stays blind to it. That asymmetry is the point rather than a side effect: it is how an effect is allowed to depend on something a prognostic function must not see. The case that needs it is a modifier from which the treatment could be reconstructed: a control function allowed to see such a variable can represent the treatment itself, so the effect is not identified, even though the effect may legitimately vary with it. Naming it as a modifier alone keeps it out of the control function. Only an explicit modifiers formula reaches such a variable: a bare vc(z) means the predictors of the model, which these are not.

A covariate with a varying coefficient is kept out of the control function, since $f_0(Z) + z f_1(Z)$ is not identified when z is among $f_0$'s predictors: any function of z moves between the two. Reached through . it is dropped quietly, because . did not name it; named outright, the model is fitted as asked and a warning says why that is a choice.

A numeric covariate may modify its own coefficient, and that is how the effect stops being linear in it. vc(z, ~ z + x1) fits $z f_1(z, x_1)$, so the slope moves across z's range and the dose response is a curve rather than a line, which is worth reaching for whenever the effect of a continuous predictor might not be proportional to it. A categorical covariate is removed from its own forests instead: a level's indicator is nonzero only on the rows where that level holds, and the variable is constant on exactly those rows, so such a split separates rows that contribute from rows that contribute nothing.

At the other extreme, ~ 1 names nothing at all, which leaves the coefficient's forest no predictor to split on: every tree in it is a stump, so the coefficient is one drawn number and x enters as a linear term while the rest of the model stays nonparametric. It is drawn under the leaf prior rather than a prior written for a regression coefficient, so it is shrunk toward zero. Comparing vc(z, ~ 1) against vc(z) with loo() is a test of whether the effect of z varies at all; vignette("varying") works one through.

Setting `center`

Centering is a reparameterization of $f_0$ alone: every coefficient and every estimand is identical under any choice, and the choice changes only what the control function means. "auto" picks by the covariate. For a 0/1 covariate it uses zero, so $f_0$ is the surface among the untreated, which is a quantity with its own meaning. For any other numeric covariate it uses the mean, because zero may be nowhere near the data and the control function there would be an extrapolation.

A factor is always fitted mean-centered and gets one forest per level, coded symmetrically rather than as contrasts against a level that happened to sort first. That coding leaves the reference a reporting choice: the center argument names the level coef() reports against, and no refit is needed to change it.

Drawing the Coding (`center = "estimate"`)

This is different in kind from the choices above. Rather than subtract a number from x, it gives each of x's values a coefficient of its own and draws it, so that every contrast carries the same prior whatever the number of values and no value is a reference. coef() returns the identified contrasts.

At two values it restricts nothing, and bcf() uses it: with it, the answer stops depending on which level was written as 1. Above two values every contrast becomes one shared shape times a scalar, where the symmetric coding gives each level its own, so it is the parsimonious model against a general one. Reach for it when the levels plausibly differ in degree rather than in kind.

Families with Several Additive Predictors

gaussian_ls(), zi_poisson() and the rest fit one forest per distributional parameter, and each parameter's formula carries its own vc() terms. The forests are then two-dimensional (a control function and its coefficients, for each parameter) and named accordingly, mean, mean:z, log_sd and log_sd:z for y ~ x1 + x2 + vc(z) under gaussian_ls(), the names per-forest settings are keyed by. One formula reaches every parameter, and a list of formulas gives each its own; vignette("varying") shows both.

So the same covariate may have a coefficient on more than one parameter: $z$ shifting the mean and widening the spread are different questions, and both are answered at once. coef() returns one column per coefficient, named for its forest.

Some consequences follow from the parameters being different. A group intercept from (1 | g) reaches every control function and no coefficient, since a group-varying coefficient is a random slope. And center = "estimate" is judged per parameter: the drawn coding needs a leaf target that is quadratic in the predictor it feeds, which gaussian_ls() is in the mean and is not in the log standard deviation, so the same request is accepted on one and refused on the other.

The two multinomial families are the exception and refuse vc(). Their forests are the levels of one parameter rather than separate parameters, identified only up to a function they all share, and reporting removes it; a coefficient forest per level would add one such direction per coefficient and the reporting does not carry them.

vignette("varying") is the worked account of all of this: each argument on the rhc data, the constant-against-varying comparison, reading the fit with coef() and estimate_effect(), and how bcf() writes its model in these terms.

Examples

data("rhc")

set.seed(123)

# The effect of right heart catheterization on death,
# free to vary with every other covariate. `rhc` is left
# out of the control function, which keeps the two
# forests identified
model <- death ~ age + sex + race + edu + aps + meanbp +
  resp + hema + pafi + paco2 + crea + surv2m + card +
  vc(rhc)

fit <- bartisan(model, data = rhc, family = binomial(),
                num_trees = 10, num_burn = 50,
                num_draws = 50)

# One coefficient per patient: the coefficient function
# at each observation
head(coef(fit))
#>            rhc
#> [1,] 0.5098281
#> [2,] 0.2377270
#> [3,] 0.2123576
#> [4,] 0.2790761
#> [5,] 0.3253521
#> [6,] 0.4182815

# The same effect, free to vary with severity of illness
# alone
model2 <- death ~ age + sex + race + edu + aps + meanbp +
  resp + hema + pafi + paco2 + crea + surv2m + card +
  vc(rhc, ~ aps)

fit2 <- bartisan(model2, data = rhc, family = binomial(),
                 num_trees = 10, num_burn = 50,
                 num_draws = 50)

head(coef(fit2))
#>            rhc
#> [1,] 0.4040734
#> [2,] 0.2496033
#> [3,] 0.1571161
#> [4,] 0.4555922
#> [5,] 0.3282443
#> [6,] 0.1985482