Set up a model formula

Specifies a model with brms-compatible syntax. Distributional parameters (dpars) can get their own formulas with the full predictor grammar, or be fixed to constants: bf(y ~ x + (1 | g), sigma ~ z + (1 | g)), bf(y ~ x, sigma = 1). Nonlinear formulas (nl = TRUE) and multivariate models (see mvbf()) use the same grammar.

Usage

bf(
  formula,
  ...,
  family = NULL,
  nl = NULL,
  center = NULL,
  cmc = NULL,
  autocor = NULL
)

Arguments

formula: The model formula for mu, or a formula bf() already built. Given one of those and nothing else, bf() returns it unchanged, as brms's bf() does; further arguments add to it.
...: Two-sided formulas for other dpars (the left-hand side names the dpar, e.g. sigma ~ z, or several sharing one right-hand side, e.g. b1 + b2 ~ 1), one-sided formulas named by their dpar (sigma = ~ z, the same formula), named scalars fixing a dpar to a constant on the response scale (e.g. sigma = 1), or a dpar equated to another of the same class, as sigma1 = "sigma2" in a mixture. See the section on equated parameters.
family: Optional family; can also be attached with + or passed to frm(), which uses gaussian() when nothing names one.
nl: Nonlinear-formula flag: the main formula becomes a nonlinear expression of named parameters, each given its own ... formula with the full predictor grammar. NULL, the default as in brms, means FALSE for a new formula and keeps the setting of a formula bf() already built. The body is R code on the tape, and it can call the functions of the Stan language that brms's own bodies call: inv_logit(), logit(), inv_cloglog(), Phi(), square(), cbrt(), exp2(), pow(), inv(), inv_sqrt(), inv_square(), hypot(), log1m(), step(), log1p_exp(), log1m_exp(), log_diff_exp(), log_sum_exp(), log_inv_logit() and log1m_inv_logit(), with Stan's values and derivatives. A function of the same name found from the formula's environment comes first, which includes the packages on the search path, except for stats::step(). x^b and pow(x, b) with x = 0 on a row have Stan's gradient. On a family with several location parameters (categorical(), multinomial(), dirichlet() and the other simplex families, a mixture()), the one body is every location's, with one set of nonlinear parameters, as in brms. On an ordinal family the body is the latent predictor, and the thresholds stay.
center: Whether the location formula's intercept is brms's class "Intercept", a density on the intercept at the means of the predictors. NULL, the default, means TRUE for a new formula and keeps the setting of a formula bf() already built. FALSE makes the intercept an ordinary coefficient, class "b" with coef "Intercept", as 0 + Intercept in the formula does; see the section on brms's reserved Intercept below. It changes priors only: the likelihood and the maximum likelihood fit are the same. It applies to the location formula alone, as in brms; give a parameter formula its own with lf().
cmc: Cell-mean coding, as in brms. NULL, the default, means TRUE for a new formula and keeps the setting of a formula bf() already built. With TRUE, R's rule applies: a formula without an intercept, y ~ 0 + g, gives a factor one column per level. With FALSE, the factor keeps its contrasts and only the intercept column goes, so y ~ 0 + g has the treatment columns gb, gc, and not ga; an ordered factor keeps its polynomial contrasts. The group-level terms of the formula follow the same rule: (0 + g | h) has the coefficients gb, gc, and so does each member of (0 + g | mm(h1, h2)). A formula with an intercept does not change. It applies to the location formula alone, as in brms; give a parameter formula its own with lf(). cmc = FALSE is refused on ar1(), hetar1(), cs(), homcs(), toep(), homtoep() and the distance structures ou(), exp(), gau(), mat(), which read one coefficient per level of their factor.
autocor: A one-sided formula of autocorrelation terms, as in brms, e.g. ~ ar(week, subj). The terms are added to the location formula, so bf(y ~ x, autocor = ~ ar(week, subj)) is bf(y ~ x + ar(week, subj)). With nl = TRUE this is the way to give a nonlinear model a term, since a body is evaluated and not read as terms. See autocor-terms.

Value

An object of class frmtmb_formula.

Details

The left-hand side accepts addition terms after |: y | weights(w) ~ ... and y | trials(n) ~ .... As in brms, an addition term's argument is evaluated on the data rows, so it can be an expression, as in weights(w * 2), trials(n + 1) or trunc(lb = min(y) - 1), and a single value is used for every row. Every variable an addition term reads must be a column of the data, as brms requires: trials(k) or weights(wt * k) with k defined only outside the data is refused by name. Predictions on newdata and refits read the term again, and a value from outside the data may be missing or different by then, most harmfully in trials(), where the change would go unseen. Function calls such as log() are allowed, and predictor formulas and offset() keep R's usual rule of looking in the formula environment too. An expression that gives NA on a row is refused, as brms refuses it. weights(w, scale = TRUE) scales the weights to a mean of one, as brms does. Weights multiply each row's log-likelihood. When every row has a random effect of its own in mu (one level per row, or ar1(t + 0 | g) with one row per cell) and the weights sum to more than the number of rows, a family whose dispersion formula has an intercept has no maximum likelihood, and frm() refuses the model by name, unless a prior on the dispersion provably bounds it. Every linear predictor accepts lme4-style random effects (1 | g), (1 + x | g), (x || g), and explicit covariance-structure wrappers us(x | g) and diag(x | g).

Attach a family with +, for example bf(y ~ x) + gaussian(), or pass one to frm(). A model that names no family is gaussian.

Equated parameters

A named character string equates a distributional parameter to another of the same class, as in brms: in bf(y ~ x, sigma1 = "sigma2") + mixture(gaussian(), gaussian()), the two components share one residual standard deviation. The model has one parameter, sigma2, and sigma1 takes its value. sigma2 is the name that fixef(), summary() and get_prior() report and that set_prior() addresses; variables() lists sigma1 as well, with the value of sigma2, as brms lists it as a transformed parameter. The fit is the same as that of a hand-written model with one shared parameter. A target that is reported as a coefficient on its link scale, such as a transition cell of frmtmb.latent::hmm(), has no natural value to repeat, and the equated name is then not listed. coef() has no block for the equated parameter.

brms's rules apply, with brms's messages. The two names must be of the same class (sigma1 and sigma2, not sigma1 and nu2). A location parameter (mu1) cannot be equated. The parameter on the right cannot itself be fixed, equated or predicted by a formula. A mixture's weights (theta1) cannot be equated: brms stops with an internal error on this, and frmtmb refuses it. lf() takes the same spelling, bf(y ~ x) + lf(sigma1 = "sigma2").

brms's reserved `Intercept`

In a formula without an intercept, Intercept is a reserved name, as in brms: y ~ 0 + Intercept + x is the model y ~ 1 + x, with the intercept as an ordinary population-level coefficient. Factors get treatment contrasts, as they do beside an intercept. The likelihood is the same, so a maximum likelihood fit is the same; what changes is the prior. brms places a class "Intercept" prior at the means of the predictors, and a class "b" prior on this coefficient at zero, and so does frmtmb: set_prior(..., class = "b") reaches it, class = "Intercept" does not. bf(y ~ x, center = FALSE) gives the same model. The spelling works in every linear formula: the location, a distributional parameter (sigma ~ 0 + Intercept + z) and a nonlinear parameter (a ~ 0 + Intercept + x, where it changes nothing, because brms never centers a nonlinear parameter).

Intercept must be a term of its own; Intercept:x is refused. An ordinal family refuses 0 + Intercept, as brms does, because its thresholds take the intercept's place. A data column named Intercept must hold only ones in such a model, as brms requires; in a formula with an intercept, Intercept is an ordinary variable. Prediction needs no Intercept column in newdata.

Missing values under `mi()`

bf(y | mi() ~ x) takes a response with missing values, for the 14 continuous families brms takes it for. What frm() does with a value depends on whether another density reads it:

No density in the objective reads it: It integrates to 1, so it leaves the objective. The fit is the complete-case one, exactly. A reader whose own density has left does not read: z | mi() ~ mi(y) with z missing where y is takes both rows out, and a reader's weight of 0 does not read. The row stays in the frame and predicts; fit$estimates$miss has no entry for it.
A density reads it: (a predictor's mi(y), mi(sd), a cov = FALSE autocorrelation term.) It stays latent and is integrated out by a Laplace approximation, on a chart that holds it in the family's support. asym_laplace() and a mixture() have no mode for it to expand about and are refused here.

The approximation is exact only for a gaussian value that gaussian means read. Elsewhere it can move an estimate. The largest shift from an exact maximum, in its standard errors, on 200 rows: 2.74 exgaussian with 60 read by a slope of 0.05, 2.14 student with the same, 0.98 beta and 0.96 Gamma with 100 read by a slope of 0.7, 0.14 exgaussian with 60 read by a slope of 0.7. A reader that carries little leaves the row nearly unread, with an unread row's error. Per value, against integrate(), on 80 rows with 8 missing and a gaussian reader: 0 gaussian, -0.114 Gamma, -0.098 exponential, -0.078 weibull, +0.008 frechet, -0.125 beta, -0.079 exgaussian, -0.076 von_mises, -0.037 lognormal, -0.023 skew_normal, +0.021 inverse.gaussian, +0.008 gen_extreme_value. von_mises loses -1 to -4.2 per value at a mean near the end of its window, where the value has two modes, and estimates there can be off by up to 17 standard errors.

After the fit, frm() warns, with class frmtmb_mi_laplace_error, when the approximation moves an estimate. Three inner Hessians over every read value give Tierney-Kadane's term, a gate that costs 0.1 percent of a fit. Past it, an integrate() over 10 of the values at the estimate and at one standard error above it, per parameter of the read response, tells how far a parameter moves: a move of 0.5 or more warns, and fit$mi_laplace holds the moves. That stage costs about 37 percent of the fit. For a von_mises value, 10 percent or more of read rows with pi - abs(mu) < 2 / sqrt(kappa) warn with no integral. check_mi in frmtmb_control() sets what happens. The shift from an exact maximum, in standard errors, over 144 fits: a shift of 1 or more warned on 75 of 76, and a shift below 0.5 on 3 of 39, and the largest silent shift was 1.01.

frmtmb.sample::frm_sample() holds every value and has no such error. Its draws name a value as brms does, Ymi_<resp>[k], on the response's scale. If the reader carries little about the value, deleting those rows gives nearly the same model without the error.

Examples

# brms-style model formulas: attach a family with `+`
bf(y ~ x + (1 | g)) + gaussian()
#> y ~ x + (1 | g)
#> Family: gaussian 
# distributional parameters get their own formulas or constants
bf(y ~ x, sigma ~ x)
#> y ~ x
#> sigma ~ x
bf(y ~ x, shape = 2) + Gamma()
#> y ~ x
#> shape = 2 
#> Family: Gamma 
# nonlinear models declare parameter formulas and nl = TRUE
bf(y ~ a * exp(-b * x), a ~ 1, b ~ 1 + (1 | g), nl = TRUE)
#> y ~ a * exp(-b * x) (nonlinear)
#> a ~ 1
#> b ~ 1 + (1 | g)
# an intercept that a class "b" prior reaches; the two are one model
bf(y ~ 0 + Intercept + x)
#> y ~ 0 + Intercept + x
bf(y ~ x, center = FALSE)
#> y ~ x
# treatment contrasts without the intercept column
bf(y ~ 0 + g, cmc = FALSE)
#> y ~ 0 + g (cmc = FALSE)
# a mixture whose components share one standard deviation
bf(y ~ x, sigma1 = "sigma2") + mixture(gaussian(), gaussian())
#> y ~ x
#> sigma1 = sigma2 
#> Family: mixture(gaussian, gaussian)