Define a model family

Constructs a family object for use with frm(). The log-density function must be vectorized and AD-compatible: it is evaluated on RTMB 'advector' objects during taping, so it must use RTMB-overloaded operations (RTMB and RTMBdist d* functions, plain arithmetic) and must not branch on parameter values.

Usage

frmtmb_family(
  family,
  dpars,
  links,
  lpdf,
  valid_y = NULL,
  init_dpars = list(),
  type = "continuous",
  post = list(),
  sim = NULL,
  sim_ctx = NULL,
  sim_refusal = NULL,
  primary_dpars = "mu",
  lcdf = NULL,
  lccdf = NULL,
  cdf_log = NULL,
  required_aterms = character(0),
  accepts_aterms = NULL,
  se_dpar = NULL,
  exclusive_aterms = list(),
  family_finalize = NULL,
  extra_pars = NULL,
  drop_intercept = FALSE,
  structure = NULL,
  aterm_data = NULL,
  simplex_pars = NULL,
  mi_chart = NULL
)

custom_family(
  family,
  dpars,
  links,
  lpdf,
  valid_y = NULL,
  init_dpars = list(),
  type = "continuous",
  post = list(),
  sim = NULL,
  sim_ctx = NULL,
  sim_refusal = NULL,
  primary_dpars = "mu",
  lcdf = NULL,
  lccdf = NULL,
  cdf_log = NULL,
  required_aterms = character(0),
  accepts_aterms = NULL,
  se_dpar = NULL,
  exclusive_aterms = list(),
  family_finalize = NULL,
  extra_pars = NULL,
  drop_intercept = FALSE,
  structure = NULL,
  aterm_data = NULL,
  simplex_pars = NULL,
  mi_chart = NULL
)

Arguments

family

Character name of the family.

dpars

Character vector of distributional parameter names. The first entry must be "mu".

links

Named list mapping each dpar to a link name or a link object. frmtmb-links lists the names and says what a link object has to carry.

lpdf

Function (y, dpars, aterms) returning the vectorized log-density. dpars is a named list of advector vectors; aterms is a named list of numeric addition-term values (for example trials). Read dpars by name, never by position: it also carries reserved entries that are not distributional parameters. Take log(mu), log(1 - mu) and 1 - mu from frmtmb-robust-dpars rather than writing them out, because an inverse link saturates and the plain arithmetic returns NaN for the value and the gradient alike in the tail.

valid_y

Optional function (y, aterms) that signals an error for invalid responses. Called once at assembly time.

init_dpars

Optional named list of functions (y, aterms) giving a response-scale starting value per dpar (applied to the intercept through the link).

A function that declares a THIRD argument, (y, aterms, resid), is given the response with its mu predictor taken out: any offset() subtracted, then least squares on that predictor's own design. A start that depends on the SHAPE of the residual rather than on the response cannot be read off y: a covariate can give the raw response the opposite skew, whether it enters as a column or as an offset (see skew_normal()). When mu is an intercept alone and carries no offset there is nothing to take out and resid is the response itself, so such a model starts exactly where it did. A two-argument function is called with two arguments, so nothing already written has to change.

The value an initializer returns names the PREDICTOR, and it is placed in whichever coefficients that design uses to carry one: an intercept when there is one, and otherwise the least-squares solution, which puts the value in every cell of ~ 0 + g. A design that was skipped for having no intercept used to start every coefficient at zero.

type

One of "continuous", "discrete", "ordinal", "categorical". The last two say the modelled response is a distribution over 1..K categories rather than a number, so fitted() returns an n x K probability matrix; "ordinal" shares one latent predictor across the categories and "categorical" gives each non-reference category its own.

post

Named list of numeric helper functions used by fitted(), predict() and residuals(): mean_fn(dpars, aterms) (the response mean), var_fn(dpars, aterms) (for pearson residuals) and dev_fn(y, dpars, aterms) (the unit deviance 2 * (loglik of the saturated fit - loglik at the fitted value), for residuals(type = "deviance")). A family that omits one is refused by the method that needs it.

post$fit_check(fit, resp) is different in kind: it is run once, when a fit FINISHES, and its return value is discarded. It is where a family says something about where the optimizer landed, which nothing else can: logLik() reads the optimizer's own value, so a family whose likelihood is floored or degenerate in some region had no way to report it. Warn from it rather than stopping; a hook that throws is caught, reported as a warning naming the family, and the fit is returned regardless.

post$curve_check(fit, resp) is the refusal that fit_check is not. frm_curve(), frm_curve_deriv() and frm_curve_feature() run it before they read a curve off a response of this family, and discard its return value. Stop from it when the fit has no curve to report, for example because its likelihood is a floor artifact at the fitted parameters: the curve of such a fit is an artifact too. The error reaches the caller unchanged. frmtmb.survival's royston_parmar() uses it to refuse a fit whose cumulative hazard is not monotone.

post$stationary is a declaration rather than a function: a named list, one entry per dpar, each list(at =, tol =, from =). It says that the likelihood has a stationary point where that dpar's PREDICTOR equals at on the LINK scale at every row. A fit whose fitted predictor is within tol of that everywhere is refitted from each value in from, with the coefficients solved so the predictor takes that value, and whichever optimum is best is kept.

Reading the predictor rather than the coefficients is what makes this work for a design without an intercept: alpha ~ 0 + g sits on the same point as alpha ~ g with no coefficient named "(Intercept)" to read or to write a restart into. Declaring a stationary point is ALSO what earns a dpar a starting value on an intercept-less design at all: a dpar that declares none keeps the plain rule, an intercept or nothing, because a least-squares start is right on the predictor scale and can be ruinous on the parameter scale. Non-finite entries in from are dropped, so from = c(Inf, -2) restarts from one side only. Declare one only for a point the likelihood has at EVERY sample, such as alpha = 0 for the skew normal, where the information is singular and no convergence test can separate a stall from a maximum. The extra fits only ever run when the declared point is reached.

mixture() carries a component's declaration through under the dpar name the mixture gives it (alpha1, alpha2, ...).

sim

Optional numeric simulator (dpars, aterms, n) returning n response draws; used by simulate(), posterior_predict() and frm_simulate(). It is stateless and rowwise: it sees the distributional parameters and nothing else. A family whose extra parameters (extra_pars) enter the draw declares a fourth argument extra instead.

sim_ctx

Optional structured simulator (ctx) for a family whose draws are not rowwise (see Structured simulators). It takes precedence over sim.

sim_refusal

Optional one-sentence reason why the family has no simulator, appended to the refusal each entry point raises. Use it when the omission is a decision rather than a gap.

primary_dpars

Which dpars receive the main model formula (default "mu"). Families with several location predictors (for example multinomial's per-category mu2, mu3, ...) list them all; these live in the beta parameter vector and are integrated out under REML.

lcdf

Optional vectorized AD CDF (q, dpars, aterms) returning probabilities; enables cens() and trunc() addition terms. It is the only thing either term asks of a family, whatever the family's type. A discrete family is read on the integer lattice (see Censoring a discrete response), so trunc(lb = ) calls this at one below the bound.

lccdf

Optional vectorized AD LOG SURVIVOR function (q, dpars, aterms) returning log(1 - F(q)) directly. A family that declares it scores the upper tail from it instead of from log(1 - F), which cannot be accurate once F rounds to one (see Censoring and the representable tail). It is optional and independent of lcdf: a family that supplies only lccdf accepts right censoring and refuses left censoring, interval censoring and trunc(), each by name.

cdf_log

Optional vectorized AD function (q, dpars, aterms) returning list(lF = log F(q), lS = log S(q)), both tails on the log scale. A family that declares it has every censored row and every truncation window read from it, and lcdf and lccdf are derived from it when they are not given. It is how a family whose two tails are both exact passes that to core: a probability-scale lcdf ends at 1e-308.

required_aterms

The addition-term values the density cannot do without, named as they reach aterms: "vint1", "vreal2", "trials". A character vector names the terms it needs ALL of. A LIST adds the alternative: an element of length one is a term that must be there, and an element of length more than one is a set of spellings, ANY one of which will do, so required_aterms = list(c("dec", "vint1"), "vreal1") reads as one of dec or vint1, and vreal1. A family that reads the same datum from either of two terms declares it that way instead of hand-rolling the refusal.

Frame assembly refuses a model that leaves a requirement unmet, naming the family, the terms and the spelling that supplies them; for a set of alternatives it names all of them and writes the first into the example formula. Without the declaration an absent term reaches the density as NULL, arithmetic on it gives numeric(0), and the log-likelihood becomes a sum over nothing: a fit that returns, with a log-likelihood of zero. Declare every per-row datum the density indexes.

accepts_aterms

The addition terms this family reads or lets the core act on, named as a formula writes them and without parentheses: c("weights", "trials", "cens"). Frame assembly refuses any other term on the response, by name, and lists the ones the family takes. character(0) declares a family that takes none. NULL, the default, accepts every registered term, which is what a family written before this argument existed keeps.

required_aterms is a conjunction of what the density cannot do without; this is the complementary allow-list, and the two are read together, so a required term need not be repeated here. Without a declaration an unread term is parsed, stored on the fit and silently ignored: wiener() accepted a vint() it cannot use, and lba() accepted a dec(), both giving a fit bit-identical to the one without the term.

Naming "se" here is more than an allow-list entry: it is how a family OPTS IN to se(). That term is the one whose entire effect is inside the density - the core hands over aterms[["se"]] and maps out the residual scale it replaces, and does nothing else with it - so a family that does not declare it is refused the term rather than given one it would ignore. NULL, which accepts every other term, declares nothing and so does not open this one. Read aterms[["se"]] as the known standard deviation, and aterms[["se_sigma"]] to honor se(x, sigma = TRUE), which asks for the known and estimated scales in quadrature.

se_dpar

The dpar that a known standard error replaces, named so that the core can map it out: se_dpar = "tau" maps out tau exactly as the convention maps out sigma. NA declares that se() replaces NO dpar, which is the shape of a family whose whole scale IS the known standard error. NULL, the default, reads the convention: the dpar named sigma, if the family has one.

Only a family that declares se() reaches this. se() without sigma = TRUE says the residual scale is known, so the dpar it replaces has to stop being estimated. A dpar the density never reads is a flat direction and a NaN standard error, so a declaring family with no sigma and another free dpar is refused: the core cannot tell a second SCALE from a genuine SHAPE. This argument is how the family says which one it has. se_dpar = NA is a promise that every remaining dpar is read alongside the known standard error, the way a skew or a tail index is.

exclusive_aterms

Sets of addition-term values that say the SAME thing to the density, so that at most one of each set may be supplied. Named as required_aterms names them, in values rather than in terms: exclusive_aterms = list(c("dec", "vint1")). A character vector is one such set; a list is several. The FIRST name of a set is the spelling the density reads, and the refusal tells the user to keep it.

This is what an allow-list cannot express, because both spellings are legitimately on it. wiener() reads its boundary indicator from dec() and falls back to vint1, so rt | dec(u) + vint(1 - u) passed every guard and fitted with a log-likelihood bit-identical to the dec()-only model: the user said one thing twice, and inconsistently, and nothing complained. Declare the alternatives of an any-of required_aterms group here when the density reads only one of them; leave them undeclared when it reads both, as gddm() does, where dec() and vint1 carry different data.

Read WITH required_aterms, not instead of it: an any-of group in one and the same set in the other together mean "exactly one". Declaring a set that holds two values required_aterms demands TOGETHER is refused here, because no model could then be fitted.

family_finalize

Optional function (fam, y, aterms) returning a family. It runs once at frame assembly, after the response is coerced and validated and before any link is used, and whatever it returns is the family the rest of the fit sees. It is how a family derives a link bound, a default, or an extra slot FROM the data instead of asking the user for a quantity the framework already holds (see Deriving a family from the data).

extra_pars

Optional function (y, aterms) returning a named list of numeric starting vectors for family-level parameters outside the dpar system (for example ordinal thresholds). They join the parameter template under their own names and reach lpdf as its fourth argument. Keep each on a scale near 1: frm(importance =) differentiates by steps of 1e-4 on it.

drop_intercept

If TRUE, the intercept column is removed from the main formula's design matrix (ordinal families: thresholds take its place).

structure

Optional frmtmb_structure() for a family whose likelihood does not factorize over rows (a group-level mixture(), a hidden Markov chain). It carries the non-rowwise log-likelihood, the frame block that likelihood reads, and the capability flags that say which post-fit methods the family can answer. lpdf stays required even then, for the rowwise contract, and may be a stub that refuses.

aterm_data

Optional function (y, aterms) returning a named list of DATA the densities read, built once from the observed response at frame assembly (step 3 of "Slot call order") and added to the addition terms every slot receives. It is for a quantity that is a function of the observed response and not a parameter, such as a spline basis over the observed times.

simplex_pars

Optional named character vector that declares which extra parameters (extra_pars) are simplexes. A name is the simplex's brms name, which is its prior class and the stem of its draw names. A value is the extra parameter that holds the simplex's D - 1 free coordinates. c(sbhaz = "sbhaz_raw") is frmtmb.survival's cox() baseline.

The densities read the weights as mo_simplex(extra[[par]], aterms[[paste0(name, "_chart")]]). The second argument is the chart of the coordinates: NULL for the chart of the fit, and "softmax" on frmtmb.sample's sampler. A declared simplex gets what a mo() simplex gets (see "Monotonic effects" in frm()): the search off a pole and off a face that the likelihood rises from (frmtmb_control(mo_search = )), the softmax chart and brms's dirichlet(1) default on the sampler, and the prior set_prior("dirichlet(...)", class = name). A name must not be one of frmtmb's own prior classes other than "sbhaz", which is brms's class for a baseline hazard.

mi_chart

The chart of a missing value of this family's response under mi(), for a family with "mi" in accepts_aterms: "log" for support on (0, Inf), "logit" for support on (0, 1), "circle" for support on (-pi, pi), and "identity" (the default) for a density on the whole line. A family with a bounded support must set one: on the identity chart a latent value leaves the support, and the density there is not a number. The latent value is held on this chart with the map's log-Jacobian in the objective, so that it is unconstrained, and its density has an interior mode whatever the family's parameters: an exponential density has none in the value itself. "no_mode" is the identity chart for a density with no curvature at its mode, such as asym_laplace(), whose log density is piecewise linear: frm() refuses a value another density reads, because a Laplace approximation needs that curvature, and frmtmb.sample::frm_sample() samples it. A value no density reads needs no chart on frm(), where its row leaves the objective.

Value

An object of class frmtmb_family. $ on it also answers brms's link fields, computed from links when read: link, the name of the link for the mean (for an ordinal family, the distribution function its thresholds are read through; "logit" for a categorical or multinomial family, whose categories share one joint link; "identity" for a family with no mu, such as a mixture()), linkfun and linkinv for that link, and link_<dpar>, the link name of every other distributional parameter, as in beta_binomial()$link_phi. Because they are read from links each time, a family_finalize() that replaces a link, or any other edit, cannot leave them stale.

The family records the package that built it: the namespace of the function that called frmtmb_family(), read with frm_family_package(), and "frmtmb" for a family built outside any package. A saved fit keeps that record. When the fit is read back in a session where that package is not installed, frmtmb stops at the first method that needs the family's functions (summary(), vcov(), predict(), fitted(), simulate(), update(), frm_curve() and the rest) with a message that names the package, rather than with the "could not find function" of one of its internals. A mixture() records the packages of its components the same way.

$ does not partial-match on a family object. fam$lpd is an error naming the field it would have matched, rather than silently returning it. Use the full name.

Structured simulators

Some families cannot draw a response one row at a time: a group-level mixture() draws one class per group, a multivariate normal mixture (frmtmb.latent's mixture_mvn()) needs the class covariances, which are family-level extras rather than dpars, and a hidden Markov family walks a chain per sequence. Those families supply sim_ctx(ctx) instead of sim(dpars, aterms, n).

ctx is a list with fit (any object carrying spec, frame and estimates - a fitted model, one posterior draw, or the de novo shim), family, rspec, resp, dpars (the evaluated numeric distributional parameters), aterms, n, extra (the family-level extra parameters) and the frame structures autocor and block (the structured family's own data; see frmtmb_structure()). Read its fields with [[ ]].

The same sim_ctx() serves [simulate()], posterior_predict() and [frm_simulate()]. Because a structured draw covers whole sequences or groups, trunc() rejection and newdata cannot apply to it and are refused.

[ ]: R:%20 [simulate()]: R:simulate() [frm_simulate()]: R:frm_simulate()

Slot call order

The order the slots run in is part of the contract, because a family that derives anything from the data depends on it. Measured on an instrumented family, one fit of one response:

valid_y(y, aterms), once, at frame assembly, with the response coerced and the addition terms evaluated.
family_finalize(fam, y, aterms), once, immediately after, and still before any link function is called.
aterm_data() and extra_pars(y, aterms), once each, at assembly.
init_dpars[[dpar]](y, aterms), once per dpar, when the starting values are built; each value goes straight through that dpar's linkfun.
linkinv and lpdf, on the tape, from then on.

post$mean_fn, post$var_fn, post$dev_fn and sim are never called by frm(). They run only when a post-fit method asks for them.

So a link may read anything valid_y or family_finalize computed, and neither of those may read anything a link produced.

Deriving a family from the data

Some families are not fully determined until the response is in hand. A shifted family whose density is zero below a non-decision time wants a link bounded above by min(y), so that the constraint is structural rather than left to the optimizer. The bound is a property of the data, and the family object is built before frm() sees any.

family_finalize closes that gap. It receives the family, the coerced response and the evaluated addition terms, and returns the family the fit will use:

shifted <- frmtmb_family(
  "shifted", dpars = c("mu", "ndt"),
  links = list(mu = "log", ndt = "log"),
  lpdf = function(y, dpars, aterms) { ... },
  family_finalize = function(fam, y, aterms) {
    # a logit onto (0, min(y)): ndt cannot reach the smallest
    # observation, whatever the optimizer tries
    ub <- min(y)
    fam$links$ndt <- list(
      name    = paste0("scaled_logit(0, ", signif(ub, 4), ")"),
      linkfun = function(mu) log(mu / (ub - mu)),
      linkinv = function(eta) ub / (1 + exp(-eta)),
      mu_eta  = function(eta) {
        p <- 1 / (1 + exp(-eta))
        ub * p * (1 - p)
      }
    )
    fam
  }
)

Replacing a link this way replaces it everywhere: the starting values, the tape, and every post-fit method read the finalized family. The alternative an extension reached for before this slot existed was to have valid_y write the bound into an environment the link closures read at call time, which works only for as long as the call order happens to hold and leaves the family object lying about what it is.

Censoring and the representable tail

For a family that declares cdf_log, every censored row and every truncation window is the mass of a window (lo, hi] (row_lpdf()), and core forms it on the log scale from both tails of the family's CDF: from lS = log S where S(lo) >= 1/2, as lS(lo) + log(1 - exp(lS(hi) - lS(lo))), and from lF = log F below that. A log near zero keeps every digit of its complement, so each form keeps the digits the other loses. A window open above is lS(lo) and one open below is lF(hi), so a right-censored row and trunc(lb = ) read the log survivor itself.

The probability scale cannot do this. A double cannot represent the complement of a probability that has rounded to one, so log(1 - F) past that point is not merely inaccurate: it is CONSTANT, and its gradient is exactly zero. An optimizer then prices such a row the same however far it moves. That failure is silent: the fit converges, and logLik() and AIC() report the floored number. Measured on a standard normal tail, one process:

z	log(1 - pnorm(z))	pnorm(z, lower.tail = FALSE, log.p = TRUE)
8.0	-35.013	-35.013
8.3	-Inf	-37.494
37	-Inf	-689.031
500	-Inf	-125007.13

Every built-in family with a CDF declares cdf_log except exponential(), weibull(), frechet(), gen_extreme_value() and discrete_weibull(). A family without it is read through lcdf and lccdf as before: a right-censored row is lccdf (or log(1 - lcdf) without it), a left-censored row log(F(y) - F(lb)), an interval log(F(y2) - F(y)), the window log(F(ub) - F(lb)), and trunc(lb = ) with no upper bound lccdf(lb). Each slot is called where it always was, because a family may read the censoring codes to know which rows a call is for, as the frmtmb.eam Wiener families do, and may supply an lcdf that is not the complement of its lccdf. The probability scale ends at 1e-308, below which a left-censored row is log(0), and an interval in an upper tail loses its digits; a family that declares cdf_log has neither limit.

The built-in cdf_log pairs compute the tail on the far side of q directly, as an integral or a series whose terms are all positive: the incomplete beta for student(), Beta(), binomial() and negbinomial(), the incomplete gamma for Gamma() and for poisson() past a count of 30 (a smaller count is a finite sum of its masses), Owen's T as an angle integral for skew_normal(), the density's own angle integral for von_mises(), and pnorm()'s two log tails for gaussian() and lognormal(). inverse.gaussian() and exgaussian() are log-scale differences. A zero-inflated or hurdle family is built from its base family's pair. Two are complements in part of the range: brmsfamily("com_poisson"), whose support has no top, takes its survivor as the complement of the summed CDF until that is below about 1e-3 and as an upper sum beyond; and a hurdle family takes F as the complement of S where S < 1/2.

RTMB::ppois(lower.tail = FALSE, log.p = TRUE) and RTMB::pgamma()'s upper tail do not tape (inside MakeADFun they reach stats::ppois() and stats::pgamma() and stop at "Non-numeric argument to mathematical function"), which is why poisson had no lccdf through frmtmb 0.72.0; its survivor is exact now (-2773.28 at q = 700, lambda = 5, where log(1 - F) is -Inf).

Censoring a discrete response

A discrete response's bounds are read as brms 2.23.0 reads them:

term	means	scored as
`cens()` `0` "none"	`Y == y`	`f(y)`
`cens()` `-1` "left"	`Y <= y`	`F(y)`
`cens()` `1` "right"	`Y > y`	`S(y)`
`cens()` `2` "interval"	`y < Y <= y2`	`F(y2) - F(y)`
`trunc(lb, ub)`	`lb <= Y <= ub`	`F(ub) - F(lb - 1)`
`trunc(lb)`	`Y >= lb`	`S(lb - 1)`

A censoring bound is the edge of a continuous-style window (lo, hi], and only the truncation lower bound is inclusive. trunc(lb = 0), or any lb <= 0 on a count, is no truncation and is accepted, as it is in brms. Every bound must be a whole number: a fractional one would be read as its floor. An interval whose ends agree (y2 == y) is an event of probability zero and is refused, as on a continuous response; an exact observation is code 0. To score "6 or more", right censor at 5.

Under trunc() a censored row is the event inside the window, on any response: right censoring at y is P(y < Y <= ub) and left censoring P(lb <= Y <= y), each over the window's mass. brms 2.23.0 divides P(Y > y) or P(Y <= y), the event over the whole line, by the window's mass, which can exceed one, so the two differ where a censored event reaches past a bound: right censoring with an upper bound, left censoring with a lower one. With no bound on that side, and for every interval, they agree.

simulate(censored = TRUE) records a draw above a right-censoring point k as k, and a draw at k is the uncensored value k, so the simulated values are the same under either reading of a bound.

residuals(type = "osa") is refused on a censored discrete fit: an uncensored row's support is [lo + 1, hi] rather than the [lo, hi] the one-step window is built on.

Tape-safe scope

lpdf and lcdf run with RTMB's tape-safe c(), [<- and diag<- in scope automatically (the "c" <- RTMB::ADoverload("c") boilerplate is spliced in unless the function already binds it), so base spellings keep the automatic-differentiation class. A helper the density CALLS still needs its own bindings: lexical scope does not travel into other functions.

Examples

# a custom family is a plain R log-density over taped parameters
dd <- data.frame(y = rbinom(100, 5, 0.4),
                 size = 5, x = rnorm(100))
fam <- custom_family(
  "vbinom", dpars = "mu", links = list(mu = "logit"),
  lpdf = function(y, dpars, aterms) {
    RTMB::dbinom(y, aterms$vint1, dpars$mu, log = TRUE)
  },
  # the density indexes vint1, so a model without it is refused
  # rather than fitted against a zero-length log-likelihood
  required_aterms = "vint1",
  type = "discrete"
)
fit <- frm(bf(y | vint(size) ~ x) + fam, data = dd)
fixef(fit)
#>              Estimate  Est.Error       Q2.5      Q97.5
#> Intercept -0.58378094 0.09329314 -0.7666321 -0.4009297
#> x         -0.02406052 0.08624634 -0.1931002  0.1449792