These functions are the interface frmtmb.sample is written
against. They were internal until the sampling surface moved to its
own package, and they are documented here as one contract rather
than one page each, because what has to be reviewed is the set.
They are lower-level than the rest of frmtmb: they take assembled
frames, resolved prior entries and single random-effect blocks, not
formulas. Most users want frm(), set_prior() and
frmtmb.sample::frm_sample() instead. Nothing here validates its
arguments the way the user-facing surface does.
The objective seam
build_objective(frame) returns the bare negative log-likelihood
closure of an assembled frame, ready for RTMB::MakeADFun(). It
honors three fields a caller may set on a COPY of the frame:
frame$map; frame$ncp_blocks, the indices of the
random-effect blocks to build in their non-centered form
(b = L(theta) z); and frame$prior_only. With prior_only = TRUE
the closure leaves out every response density and keeps the rest:
the densities of the random-effect blocks, of the me() latent
values and their measurements, and the measurement density of a
mi(sdy) response, observed values around their latent ones. A
missing value of an mi() term has no density left, as in brms.
These are the terms brms keeps outside
if (!prior_only) under brm(sample_prior = "only"), so this closure
plus neg_log_prior_fn()'s is the prior of the whole model. A
prior-augmented objective is this closure plus
neg_log_prior_fn()'s.
gp_tape_scope(frame, expr) evaluates expr, the call that builds a
tape of such a closure, with the tape settings a frame with an exact
gp() block needs: the field's inner Hessian is dense, and without them
its tape passes 13 GB at 500 positions. They change no value. Its parameter
list goes through gp_tape_order(frame, pars), which puts b first: the
settings give a wrong gradient otherwise. With hessian = FALSE,
gp_tape_scope() leaves out the inner Hessian, for a tape that never reads
it.
gp_nugget_message(fit) is the warning of a fit whose noise has reached an
exact gp()'s nugget, or NULL.
row_lpdf(fam, y, yobs, dpars, aterms, extra) is the per-row
log-density composition the objective itself runs, with cens() and
trunc() folded in; it runs on numeric dpar values as readily as on
the tape, which is what makes a pointwise log_lik() reproduce the
fitted density exactly instead of approximating it.
with_cs_offsets(fit, rspec, dpv) takes the WHOLE list
eval_dpars() returns, one element per response (rspec is
ignored), and gives it back with the category-specific (cs())
offsets of the fitted rows added as each response's .cs; on a model
without cs() terms it returns dpv unchanged, so it is safe to
call unconditionally before row_lpdf(). In frmtmb 0.62.0 and
earlier this paragraph said it took one response's list; a caller
that passed one got no offsets, because they were written a level
deeper than the simulator reads.
ce_by_response(resps, one, effects) joins the conditional-effects panel
sets
one(r, effects) gives each response, in brms's order. ce_na_vars(x, pop_level) names the grid columns a display holds at NA.
autocor_for_newdata(fit, ac, rspec, newdata) is a residual correlation term
over the groups and times of newdata. resp_refuse_single() is brms's
refusal of a method that reads one response of a multivariate model.
update_merge_prior(old, new) is brms's merge of update(prior = ): new
rows take their slots and the other old rows stay.
cs_offsets_add(fit, resp, newdata, dpv) takes ONE response's
dpar-value list, however it was built, and adds that response's
.cs evaluated at newdata (the fitted rows when NULL); it is the
one to use when drawing at new data. us_chol_cor(theta, K) is the
unstructured correlation matrix of a thetar segment, which a
set_rescor(TRUE) model's joint row density needs.
rescor_row_loglik(fit, dpv) is that joint density itself, one value
per row, at the dpar values dpv (the whole list eval_dpars()
returns): the multivariate normal of a gaussian model or the
multivariate t of a Student-t one, with its one shared nu. A
pointwise log_lik() of a rescor model reads it rather than
assuming the normal.
Its resp argument (default every response) names the responses
whose joint density is wanted, as brms's log_lik(resp =) asks.
predict_rescor_draw(fs, rspecs, dps, ok) is one JOINT predictive
draw of the responses rspecs (a named subset of
fs$spec$responses) at the parameters of the fit fs, on the
rows where ok is TRUE: dps holds each response's dpar values,
with mu and sigma, paired with rspecs by name when it is
named (the whole list eval_dpars() returns serves any subset) and
by position otherwise. It is brms's multivariate normal draw, or its
multivariate t with the one shared nu, and it returns a list of
one vector per response, or NULL when no row is ok. Both
predict() and posterior_predict() of a rescor model draw through
it, so the two cannot drift apart.
autocor_row_loglik(fit, dpv, resp) is the pointwise log-density of
a response with a residual covariance term (the covariance form of
ar(), ma() and arma(), cosy(), unstr(), fcor(), a term of
frmtmb_register_autocor()), at dpv (the whole list eval_dpars()
returns) and the fit's thetaac: brms's leave-one-out conditional
density of each row given the other rows of its group, one value per
row (frmtmb-autocor, "The pointwise log-density"). The values do
not sum to the log-likelihood. On a maximum-likelihood fit it is the
value at the estimates.
frame_ac_terms(frame) lists a frame's autocorrelation terms by
response: the residual terms of frame$autocor, and each latent residual
(frmtmb-autocor) with theta_idx into theta past its log sderr.
latres_new_spec(fit, rspec, newdata) builds what a draw of the latent
residuals of newdata's rows needs, once per call, or NULL for a response
with no such term, and refuses new rows under cov = FALSE in brms's words.
latres_shifted(fit, lat, resp) gives the fit back with one draw of them
in the predictor of mu. autocor_cond_strip_latent(fit) marks a fit whose
predictions read a grid: its latent residuals are not read there, as brms's
incl_autocor = FALSE drops them.
arma_cond_resp(fit) names the responses carrying brms's default
cov = FALSE form of ar(), ma() or arma(), and
arma_cond_dpars(fit, dpv) gives the whole eval_dpars() list back
with each such response's mu moved to brms's one-step conditional
mean, the regression on the earlier OBSERVED residuals of the row's
group. The two answer one question between them: which pointwise
log-density a row gets, and what mu it is evaluated at. Under
cov = FALSE every row keeps the family's own density; a cov = TRUE
covariance block has none, and its row gets autocor_row_loglik()
instead. Call arma_cond_dpars() before row_lpdf() (and before
rescor_row_loglik()), which is the order the taped objective uses
and which makes the row density read the objective's own mu instead
of a second one. It is safe to call unconditionally: with no such
block it returns dpv unchanged.
arma_cond_fill_dpars(fit, rspec, newdata, dpars_fn) gives one
response's distributional parameters on newdata for a predictive
draw, where dpars_fn(fit) computes them for a fit. Under
cov = FALSE it takes brms's treatment of a missing response: a row
of newdata whose response is NA (every row, when the column is
absent) is filled with a draw from the family at its one-step mean,
in each group's time order, before the rows after it read its
residual. It uses the random-number stream, and with no missing
response it is dpars_fn(fit).
arma_cond_fill_epred(fit, rspec, newdata, re_formula, dpar) is the
expected response of one draw on newdata with brms's fill: each
missing cov = FALSE response is a draw at its shifted mean, so
posterior_epred() carries the spread of the unobserved past, as
brms's does. It returns NULL when nothing needs filling, and the
caller keeps its own route.
response_codes_newdata(rspec, y, what) codes a response read from
newdata as the fit coded its own: a bernoulli response takes the 0
and 1 of the fit's two values (brms's level order), and is refused,
naming what, when it holds a third value. Any other response is
returned unchanged.
subset_resp_check(fit, resp, what) and
subset_newdata(fit, resp, newdata) carry brms's rule for a
multivariate model whose responses use subset(): each response has
its own rows, so a method is asked for one response at a time, and
on newdata it answers for the rows where that response's subset is
TRUE. The first refuses anything but a single resp on such a model
and is silent on every other; the second returns newdata unchanged
when the response has no subset().
A mo() simplex is held in coordinates whose chart suits maximum
likelihood and not a sampler (?frm, "Monotonic effects").
mo_simplex(z, chart) maps coordinates to the simplex and
mo_coords(w, chart) back, in the fit's chart (chart = NULL) or
the softmax (chart = "softmax", log ratios to the first weight).
mo_chart_frame(frame, chart) is a copy of the frame whose objective
reads every simplex in chart, and mo_frame_terms(fit) lists each
term's simplex component, its size D (0 for one category, 1 for two,
a weight fixed at 1), brms's names for its weights,
simo_<label>1[k], and the label; simplex_frame_terms(fit) adds
each simplex a family declares (frmtmb_family(simplex_pars =))
under its brms name, sbhaz[k] for frmtmb.survival's cox() baseline,
with kind the
family's name; the chart frame reads those the same way.
summary_mo_frame(fit, prob) is the weights'
table summary() prints: each weight with its delta-method standard
error and interval, NA for a weight at 0 or 1 (named by the
face attribute).
The prior seam
The prior VOCABULARY - set_prior(), prior_normal() and its
relatives, get_prior(), prior_summary() - is ordinary exported
API and is not part of this page. What is here is the RESOLUTION
machinery underneath it: as_priorlist() coerces the accepted
spellings (including a brms brmsprior) to one,
resolve_prior_input() turns a priorlist plus a model into
list(entries, lower, upper) on the internal parameter scale,
neg_log_prior_fn() turns resolved entries into a tapeable closure,
resolve_bounds() turns user-spelled bounds into internal-scale
vectors over the outer parameters, spec_target() names the slot
one specification addresses, and spec_spelling() gives back the
class and dpar it was WRITTEN with, which for a density on a
distributional parameter itself is that parameter's own class rather
than the intercept slot the resolver assigns to.
vcov_estimated(fit) is the covariance of EVERY estimated
coefficient under the internal names, which is what pairs with
estimated_coef_names(); vcov() itself covers brms's
population-level block and leaves out an intercept-only
distributional parameter. It does NOT cover an ordinal fit's
thresholds or its cs() coefficients, which live outside beta
and betad: the interop seams pair with a longer vector that has
them, and vcov(fit, full = TRUE) is where they are under
confint()'s names.
fam_is_category_valued(fam) says whether a response is a set of
categories, so that a predictive summary is a set of proportions
rather than a mean, and predict_category_props(fit, rspec, draws)
builds brms's P(Y = k) matrix from a matrix of simulated or
sampled responses. brms_summary_matrix(), brms_summary_array(),
brms_summarize_draws() and brms_prob_cols() build brms's
Estimate / Est.Error / Q columns from a point estimate with a
standard error, or from draws. brms_fixef_rows(fit) gives brms's
population-level rows in brms's order: which coefficient-table rows
they are, and the ordinal thresholds and cs() coefficients that are
not, each with the map from the internal vector to the reported value,
so a draws object can report the same rows fixef() does.
ord_delta_info(fit) lists brms's delta parameters of an ordinal
fit with equidistant thresholds: the name, where the internal value
sits, and the map from it to delta. ord_thres_linpred(fit, newdata, resp, re_formula, allow_new_levels) is brms's linear
predictor with the thresholds included at the fit's estimates, the
n x K1 matrix one draw of posterior_linpred(incl_thres = TRUE)
holds.
ord_lp_blocks(frame, spec) lists an ordinal model's threshold
blocks with the linear predictor and the dpar that class
"Intercept" addresses each by, which is where brms puts its
default on the thresholds.
autocor_flat_priors(fit) lists brms's flat default on each
residual-correlation class held through a transform (ar, ma,
cosy, a registry term's lagsar): one row per class and response,
with the ends of a uniform density on the natural parameter that
spans the whole image of the term's map, so that it is flat
wherever the parameter can go, and the narrower bound brms declares
where it has one (lb = 0 on cosy). A cov = FALSE term has no
row, because its coefficients are the internal parameters.
autocor_bound_chart(fit, lower, upper) gives, for each bounded
first-order residual-correlation parameter, a chart that samples it
through one logistic onto its natural interval, with the map to the
internal parameter and its log Jacobian.
frmtmb_register_prior_defaults() is the other direction: it lets a
package tell get_prior() what defaults it would apply.
get_prior() reads the registry under route = "sample" ONLY, and
refuses that route when nothing is registered. Its default
route = "fit" reports what frm() applies and reads nothing, so a
registration cannot change the answer of a call that did not ask
about the sampling route.
The covstruct block readers
Each takes ONE element of frame$re_blocks and answers one
structural question about it. ncp_eligible() says whether the
block has a linear Cholesky factor and only standard-deviation and
correlation parameters, so that b = L(theta) z is a bijection;
ncp_scale_b() and ncp_unscale_b() are that map and its inverse.
covstruct_has_chol() is the first half of eligibility on its own
(a factor is registered at all), block_sd_idx() gives the block's
standard-deviation positions within its theta segment,
block_n_cor() counts its correlation positions, block_cor_prior()
says whether an LKJ density fits them ("none", "lkj" or
"unsupported"), and is_student_block() reports a Student-t
latent, which is a scale mixture and so has no linear factor at all.
The registry these read is not exported as a list. A package adds a
structure to it with frmtmb_register_covstruct(), whose arguments
are the API, and asks about a block through these eight questions.
The simulator seam
sim_can() and sim_note() say whether a family can be simulated
from and what is missing when it cannot; sim_context() builds the
context one draw is taken from and sim_draw() takes it.
sim_is_structured() reports whether the draw walks a fitted
structure - a hidden state sequence, a group-level latent class, a
correlated residual - which is what makes newdata and re_formula
refusable rather than merely unimplemented on a predictive method.
brms's names and brms's VarCorr
brms_coef_names(fit) gives brms's name for each estimated
coefficient, in estimated_coef_names() order (b_Intercept,
b_sigma_Intercept, or sigma for a sigma with no formula); brms_re_rnames(fit, bk) gives brms's
coefficient names for one random-effect block (Intercept,
sigma_Intercept), and brms_block_has_r(bk) says whether the
block's sampled values are brms's r_ coefficients.
brms_coef_table(fit) adds, per coefficient, whether brms reports it
as the distributional parameter itself on the natural scale (sigma)
with its linkinv and linkfun; brms_stan_name() is brms's
response spelling (y_a is ya), brms_group_name(bk) its group
spelling, brms_levels(bk, for_r) its level spelling, and
brms_re_parts(fit, bk) a block's prefixes and renamed coefficients;
brms_par_labels(fit, include_random) labels a
sampled parameter vector in template order, in brms's spelling where
brms has a parameter with the same content, and
gp_brms_natural(fit) maps the theta entries it labels sdgp_
and lscale_ to brms's scales and back. mi_sampling(fit) gives a
fit back the mi() rows maximum likelihood left out of its objective,
which a sampler keeps; mi_sampling(fit, "laplace") gives an unfitted
object that design, for a sampler that integrates the values out;
mi_sampling(fit, "natural") maps a missing
value from the chart the tape holds it on to the value brms reports
under Ymi_<resp>[k], and back. varcorr_matrices(fit, theta) is the covariance matrix of every block at theta,
varcorr_layout(fit) is brms's VarCorr() grouping for the fit,
and varcorr_values(fit, vals, comp, layout) gives the standard
deviations, correlations and covariances of that layout at one
parameter vector in the hyp_vals_only() layout.
Parameter labeling and fitted quantities
par_name_bare() is the draws-side spelling of a parameter name,
with parentheses dropped (Intercept, not (Intercept));
outer_par_names() names the outer parameter vector in
confint() row order; estimated_coef_names() names the estimated
coefficient vector with mapped betad entries excluded;
log_sd_theta_index() says which theta positions are log standard
deviations, so that a variance at its boundary can be told from an
AR(1) coefficient at its own. sdr_of() is the cached
TMB::sdreport(), and require_fitted() is the refusal a
maximum-likelihood-only method owes an object that was assembled but
never optimized.
The hypothesis expression engine
hyp_parse_all() turns hypothesis strings into expressions and
per-expression directions, once. hyp_vals_only(fit) reads the
fit's current estimates into a flat named value vector plus a
parallel component vector (which template component each value came
from); hyp_env_vals(fit, vals, comp) takes exactly that pair and
builds the evaluation environment, and hyp_eval() evaluates one
expression in it (hyp_eval_in(ex, values) over any named list keyed
by parameter name, and hyp_expr_vars(ex) the names it reads) - so a
caller with many parameter vectors parses
once, then per vector swaps the estimates in and calls
hyp_vals_only() and hyp_env_vals() again.
hyp_tail_p() is the tail probability of a directional claim.
hyp_class_prefix(class, group) is brms's name prefix for its
class and group arguments. hyp_labels() writes brms's
Hypothesis label for each string, hyp_samples_frame(m, k) lays a
draws matrix out as brms's samples frame (H1, H2, ...), and
hyp_brms_result() assembles the brms-shaped list every
hypothesis() method returns, so no method builds that shape on its
own.
The conditional-effects engine
ce_grids_build() builds the prediction grids, effect list,
condition sets and base values; ce_boot_one() evaluates one grid
at one parameter vector and flattens it; ce_frame() lays out one
grid's rows in the column order brms's frame carries (the varied
predictor(s), every other model variable at its held value,
cond__, cats__, effect1__, effect2__), and ce_finalize()
assembles the per-effect data frames into the returned object with
the attributes plot() reads. ce_cats_display() says whether the
display is per-category and ce_display_kind() which of the three
displays a call asks for, ce_structure_check() is the refusal a
structured likelihood owes a grid, and ce_re_formula() resolves
the random-effect argument of such a call.
ce_pred_dpar() says what to PREDICT, which is not always the dpar
the display is labeled with. It returns NULL - the expected
response - whenever no dpar was named and the family's mean is not
the inverse link of mu, or the response carries trunc() bounds.
Call it rather than passing a resolved dpar straight to
ce_boot_one(): on a zero-inflated fit the mu predictor alone is
1 / (1 - zi) times the mean, on a mixture it is one component's
mean, and on a response truncated below at zero it can be negative.
The predicate underneath stays private, because what an extension
needs is the decision and not the test that makes it.
Three serve brms's display options on both surfaces.
ce_spaghetti(grid, m, effects) lays out brms's "spaghetti"
attribute from a draws-by-rows matrix m over the rows of grid;
ce_spaghetti_check() raises brms's refusal of spaghetti together
with surface; ce_check_distance() validates select_points and
too_far. ce_grids_build() takes surface and too_far, and
ce_finalize() takes select_points, pred_vars (the grid's own
pred_vars) and draws_by_eff, a list by effect of the draws
matrices whose columns are that effect's rows.
conditional_smooths() has its engine here too. cs_build() returns
one grid per smooth term, each with the design A of the term over
its coefficients; cs_coef(fit, tm) is those coefficients at the
parameters a fit object carries, so a sampler takes one product
A %*% cs_coef() per draw; cs_frame() lays one term's summary out
in brms's columns and cs_finalize() makes the returned list.
cs_probs() resolves brms's prob and deprecated probs.
ce_dots() is the argument surface itself: it pulls
allow_new_levels out of a call's dots, returns whether it was set,
and REFUSES whatever is left, naming it. An extension that
hand-rolls the same check accepts a different set of arguments from
the fit method, which is how allow_new_levels came to work on a
fit and warn on draws.
Three more serve one purpose between them: brms's rule for which
group a grid row reads. ce_group_vars() names the grouping
variables a re_formula = NULL call blanks where the grid does not
set them. ce_level_plan(fit, grids, base, re_formula) then decides,
per row and per group-level block, whether the row reads a level the
fit saw or a NEW one (unset, unseen, or a combination a nested term
never saw), and places each new one at a placeholder level of its
block. A row reads only the block of its own level of a gr(g, by = f) term, and each new member value of an mm() term is its own new
level. ce_plan_eval(f, plan, gi, cache, ...) evaluates grid gi at
the parameters of f, drawing every new level from the covariance
f's own theta implies, once per distinct level: a sampler calls
it per POSTERIOR DRAW, with a fresh cache environment per draw,
where the fit method calls it per bootstrap replicate, so the
construction is the same and the frames stay comparable.
ce_plan_has_new() says whether a plan draws anything. They replace
ce_new_level_spec(), ce_boot_grids() and ce_draw_new_levels(),
which placed one new group over every row and are removed.
method = "predict" has three more. ce_method() resolves brms's
method names ("posterior_predict" is "predict"),
ce_predict_family_check(rspec) refuses a family that cannot draw a
response row by row, and ce_sim_one(fit, nd, resp, re_form, allow_new_levels) draws one response per row of grid nd at the
parameters of fit, so that one call per posterior draw gives the
posterior predictive draws. ce_plan_eval() takes it through its
one argument, a function of the fit and a grid, so that a
simulated response reads the new-level draws of the plan.
pp_check_resp(fit, resp) is the response pp_check() checks, with
brms's refusals of a missing, plural or unknown resp on a
multivariate model, so the fit and draws methods say the same.
The two-dialect argument seam
frmtmb answers to two argument dialects: a brms-named function takes
re_formula, frmtmb's own fit surface takes lme4's re_formula, and
the brms-named ones accept both. arg_unset() is the "not supplied"
marker a formal defaults to when NULL and NA are both real
settings and neither can double as unset; re_form_arg() resolves
the pair, refusing rather than guessing when both are given.
The borrowed-generic seam
An extension that DEFINES a generic under a name another package
owns breaks that package. UseMethod() reads the method table of
the namespace where the generic it reached was defined, so the
extension's generic, reached first on the search path, sends the
owner's objects into a table with no method for them. Delayed
registration with S3method(pkg::generic, class) does not prevent
this: it puts the extension's method in the owner's table, and it
does nothing about the extension's rival generic.
frm_install_generics(pkgname, owners) is the repair frmtmb applies
to itself. owners is a named list with one entry per generic name,
each a character vector of the packages that own that name, in the
order to prefer when more than one is loaded and none is attached.
Each name becomes an active binding that returns the owner's generic
while the owner is loaded, and the package's own generic while it is
not. Call it FIRST in the extension's own .onLoad(), with the
extension's pkgname: the binding can only be installed while the
namespace is unsealed, which is during that call and at no later
time.
Three rules come with it. The extension's own generic must be a bare
UseMethod(), because whenever the owner is loaded that generic is
not the one that runs. Its formals must match the owner's. And every
method must be registered twice, with S3method(generic, class) for
the package's own table and S3method(owner::generic, class) for
each owner's, or it is unreachable while that owner is loaded.
See also
frmtmb-extension-api for the family-level accessors a
structured family uses, frmtmb_structure() for the protocol
those serve, and frmtmb_register_compat() for the compatibility
registry, which was documented here while it had no page of its
own.
Examples
set.seed(1)
dd <- data.frame(x = rnorm(60), g = factor(rep(1:6, 10)))
dd$y <- rnorm(60, 1 + 0.5 * dd$x + rnorm(6, 0, 0.5)[dd$g], 1)
fit <- frm(bf(y ~ x + (1 | g)) + gaussian(), data = dd)
#> Boundary (singular) fit: sd_g__Intercept is at zero. The likelihood does not change as it moves further that way, so it has no standard error (NaN in summary(), vcov() and confint()); the other standard errors are kept. lme4 reports such a fit as singular; diagnose() lists the components at the boundary
# the objective seam: the bare likelihood closure of the frame
nll <- build_objective(fit$frame)
nll(fit$estimates)
#> [1] 19.59269
# a block reader: this one has a diagonal factor and one sd, so it
# can be sampled non-centered
bk <- fit$frame$re_blocks[[1]]
c(eligible = ncp_eligible(bk), n_cor = block_n_cor(bk),
cor_prior = block_cor_prior(bk))
#> eligible n_cor cor_prior
#> "TRUE" "0" "none"
# parameter labeling, in the two spellings
head(outer_par_names(fit))
#> [1] "(Intercept)" "x" "sigma_(Intercept)"
#> [4] "theta_1"
head(par_name_bare(outer_par_names(fit)))
#> [1] "Intercept" "x" "sigma_Intercept" "theta_1"