Runs tmbstan::tmbstan() on the model's objective and returns the
draws with frmtmb parameter names. Given a frmtmb::frm() fit it samples the
fitted objective, initialized at the ML estimates, which shortens
warmup considerably. Given a formula and data it assembles the same
objective without optimizing anything first (see Sampling from a
formula).
Arguments
- fit
A
frmtmb_fit, or abf()/formula to assemble and sample directly (thendatais required), or afrmtmb::frm_multiple()result (see Multiply imputed data).- data
A data frame of model data, when
fitis a formula.- family
Family, when
fitis a plain formula that does not carry one (frm_sample(bf(y ~ x), data = dd, family = poisson()); the+spellingbf(y ~ x) + poisson()works too).- ...
Passed to
tmbstan::tmbstan()(chains,iter,laplace(see Laplace draws),cores, ...).coresparallelizes over chains on every platform. On Windows the chains run on socket workers, each of which rebuilds the tape from the serialized objective closure (tmbstan retapes on the worker; the closures are self-contained), giving draws identical to a sequential run at the same seed. The per-worker startup (a new R process, package load, retape) is a fixed several seconds, so short chains gain nothing; long chains approach a per-chain speedup. A core count inherited fromoptions(mc.cores), which rstan reads whencoresis not given, behaves the same way.seedfixes Stan's stream, the init jitter (seeinit_jitter) and the prior draws ofsample_prior = "yes";NA, brms's default, leaves all three to the session.save_warmupis rstan's and brms's:TRUE, the default, keeps each chain's warmup iterations beside the draws, whereinc_warmup = TRUEreads them (see Warmup).- prior
Priors: a
frmtmb::set_prior()specification, abrmsprioror a prior table, or the string"flat". A specificationset_prior("", class = ...)(or a table row whose prior is(flat)and that carries no bound) is brms's(flat)on that slot, in place of its default. A slot with no default of its own, such as agrouporcoefrow under a class row, inherits the class row's density then, as brms's empty row does, so afrmtmb::default_prior()table ofroute = "sample"with its class row edited applies the edit to every group. A bound on one coefficient (coefwithlborub) is refused with brms's message: brms, and this function, take bounds per class.set_prior(class = "Intercept", nlpar = "a")is brms's class rowclass = "b", nlpar = "a"wherea ~ 1, its bound included; with a bound, it is refused whereahas more than an intercept. The objective is re-taped with the prior terms added; a fitted model itself is unchanged. On both paths the brms default priors apply to whatever this argument, and a MAP fit's own prior, leave alone (see Default priors), andprior = "flat"replaces them with brms's(flat)in every slot, flat on brms's own parameters. Abrmspriorobject built by brms's ownprior()is translated row by row. The argument takes brms's spelling,prior. On a fit that carries its own prior, a specification here for a slot the fit's prior also names (the same class, coef, group, resp, dpar and nlpar) REPLACES the stored one, as brms'supdate(prior = )replaces a row; it is not refused as a duplicate, and specifications for other slots stack with the stored ones, the more specific applying. Two specifications for one slot within this argument are refused. Thepriorsof releases before 0.43 is gone rather than aliased, and because this function's...would otherwise swallow it, the old name is refused by name.- sample_prior
brms's
sample_prior:"no"(the default),"yes"or"only", orTRUEfor"yes"andFALSEfor"no"."yes"adds draws from the priors to the draws, under brms's names (prior_b,prior_Intercept,prior_sigma, ...), whichprior_draws()returns andfrmtmb::hypothesis()reads for the evidence ratio of a point hypothesis."only"samples the prior of the whole model with the likelihood left out, so thatposterior_predict(),pp_check()and conditional effects give the prior predictive distribution. See Prior draws.- init
Initialization. On a fit the default (
"last.par.best") starts chain 1 exactly at the ML mode and every further chain at the mode plus a normal perturbation of sdinit_jitteron the unconstrained scale. The mode anchor keeps warmup short; the jitter keeps the chains overdispersed enough for Rhat to retain power against multimodality (the standard objection to identical mode starts)."random"requests Stan's own overdispersed initialization, and is the default from a formula, where there is no mode to start at.- init_jitter
Per-chain perturbation sd for the default init;
0starts every chain exactly at the mode. Withseedgiven (in...), the jitter is drawn from that seed rather than from the session's random stream, soseedreproduces every chain, asbrm(seed = )does. Without one it is drawn from the session's stream, soset.seed()decides it. Either way rstan's sampler still advances the session's stream (it draws each chain's permutation of its draws), as it does underbrm().- reparameterize
Sample the qualifying random-effect blocks in their non-centered form (
b = L(theta) z), which removes the funnel of the centered joint posterior.TRUEby default. Not every block qualifies, and the call says which did not and why; see Reparameterization, which also explains why nothing downstream can tell the two routes apart. The mode-anchoredinitstill starts at the ML mode: it is mapped throughz0 = L(theta_hat)^-1 b_hat. Seeded draws differ between the two routes, soreparameterize = FALSEis also the way to reproduce a run made before this default existed.- data2, start, fit_control, na.action, REML
As in
frmtmb::frm()(fit_controlis that function'scontrol). They belong to the formula path, where they assemble the model this function samples, and on a fitted object they are refused by name rather than accepted and discarded: the model is already assembled, so nothing would read them.fit_controlis not the sampler's control list;controlis.- control
The SAMPLER's control list, brms's spelling and brms's meaning:
control = list(adapt_delta = 0.99, max_treedepth = 12)reachestmbstan::tmbstan()and then rstan unchanged, so abrm()call ports across with this argument included.NULL, the default, leaves rstan's own defaults in place, and so doeslist(). The names are rstan's (adapt_delta,max_treedepth,stepsize,metric, ...) and an unrecognized one is refused by name, because rstan answers an unknown option by returning an empty fit rather than by raising.frmtmb_control()fields arriving here are refused by name too and point atfit_control: earlier releases spelled the fit-time optionscontrol, and a silent reinterpretation would handgrad_tolto the sampler. So is an ABBREVIATION of the name: onlycontrolbinds to this argument, and a shorter spelling would otherwise travel through...and partial-match rstan's owncontrol, reaching the sampler unchecked.- .diagnostic
Internal, and named with a dot because there is no reason to set it by hand.
TRUEturns the default priors off and silences the disclosure, and leaves the fit's own objective alone: a MAP fit's penalty is taped into it, so no prior is resolved again and no tape is rebuilt.check_laplace()sets it, because it measures the Laplace and Wald approximations against the density the fit maximized and a default prior would change what is being measured. Writeprior = "flat"for the user-facing opt-out, which goes further and rebuilds a MAP fit's objective without its penalty.
Value
An object of class frmtmb_draws: list with the stanfit,
a draws matrix with named columns (as.matrix() method), the
originating fit, and, when any block was non-centered, a
reparam note saying which. The stanfit holds the parameters as
Stan sampled them, so on a non-centered run its random-effect
columns are z while the draws matrix holds b. The draws matrix
also carries brms's covariance summaries under brms's names,
sd_<group>__<coef>, cor_<group>__<coef1>__<coef2>, a smooth's
sds_, rescor__ and the residual structures' ar[1], cosy and
cortime__, which frmtmb::variables() lists for the fit. They
are derived from each draw's internal covariance parameters
(theta_1, ...), which stay in the matrix because the model reads
them, so as_draws_df(draws)$sd_g__Intercept works as it does on a
brms fit. With sample_prior = "yes" the matrix ends with the
prior draws. The object also holds what converts the saved warmup
in the stanfit to the draws' columns (and the prior draws of its
rows), the sample_prior setting and the arguments of the call,
which update() reads.
Details
Both routes sample a posterior under the same default priors (see Default priors), brms 2.23's own weakly-informative ones. What the fit adds is a starting point, not a different density: the ML mode anchors the chains and shortens warmup, and the model is the one the fit already carries.
The two differ only in what else is in the stack. A fit made with
frm(prior = ) is a MAP fit, and it carries its own prior into
sampling: the priors that apply are then this call's prior = first,
the fit's own next, and the defaults under both, most explicit
winning per slot it addresses. prior = "flat" opts out of all of
it and puts brms's (flat) in every slot, on either route (see
Default priors).
check_laplace() is the one caller that stays unpenalized by
default. It measures the Laplace and Wald approximations against the
shape of the very objective the fit maximized, and a default prior
would change the thing being measured, so it asks for the flat
density explicitly rather than inheriting whatever this function
defaults to.
Sampling from a formula
frm_sample(bf(y ~ x + (1 | g)), data = dd, family = gaussian())
parses, assembles and tapes exactly as frmtmb::frm() does, stops before the
optimizer, and hands the objective to Stan. Every pre-optimizer
refusal still applies (REML, quadrature, the mixture and frmtmb.latent::hmm()
guards). There is no mode, so the default init is "random":
Stan's own overdispersed initialization on the unconstrained scale,
inside any bounds the prior carries.
The returned object supports the whole draws surface -
summary(), frmtmb::fixef(), frmtmb::VarCorr(), frmtmb::ranef(), frmtmb::hypothesis(),
posterior_epred(), posterior_predict(), posterior_linpred(),
frmtmb::pp_check(), frmtmb::as_draws() - because those read the model frame and
one draw at a time. The methods that report a maximum-likelihood
quantity refuse instead of inventing one: check_laplace() (it
compares NUTS against a mode that does not exist here), and on the
embedded object reachable as x$fit, summary(), vcov(),
confint(), logLik(), frmtmb::fixef(), frmtmb::ranef(), frmtmb::VarCorr(),
predict(), fitted(), residuals() and simulate().
Reparameterization
A random-effect block has a FUNNEL in its centered joint posterior:
the width of the prior on b is a standard deviation that is being
sampled at the same time, so the region NUTS must explore narrows as
that standard deviation shrinks, and one step size cannot fit both
ends of it.
reparameterize = TRUE (the default) samples z ~ N(0, I) instead
and computes b = L(theta) z on the tape, with L the block's own
Cholesky factor, which is brms's construction. Each draw is mapped back
through ITS OWN theta, so the random-effect columns of the draws
matrix
hold the same quantity in the same order under the same names as
reparameterize = FALSE gives, and every method downstream
(posterior_epred(), frmtmb::ranef(), frmtmb::log_lik(),
frmtmb::loo(), frmtmb::conditional_effects(),
frmtmb::hypothesis()) reads them without knowing
which route produced them. Only the stanfit inside the object
carries z.
Priors are unaffected. They are declared on the OUTER parameters
(beta, theta, betad), and theta means the same thing on both
routes; nothing prior-able is reparameterized. What priors DO decide
is which blocks are eligible, below.
Which blocks, and why not all of them. A block is non-centered when
two things hold, and they are the same thing twice: every parameter it
has is a standard deviation or a correlation with a Cholesky factor
registered for its structure, and every one of those parameters
carries a PRIOR. Both are about not handing the sampler a direction
it can run away in.
Non-centering gives NUTS the run of theta's whole range, and the
centered funnel is what was keeping a chain out of the parts of that
range where the posterior is flat or improper. Removing the funnel
without closing those off first trades a slow chain for a wrong one.
In practice that means a default-prior run, on EITHER route, whose
priors cover every standard deviation and every correlation,
non-centers
(1 | g) and any block written one term at a time, diag() and
homdiag() blocks, mgcv::s() smooths, equalto(), gr(cov = ),
and the CORRELATED blocks (x | g), cs(), ar1() and hetar1().
A k = Hilbert-space gp() block non-centers too, since its
factor is diagonal and the defaults prior its sd and its length
scales (brms's classes "sdgp" and "lscale").
rr() and an exact gp() are already non-centered by
construction: the coefficients of rr() are its standard normal
factors, and those of an exact gp() are brms's zgp, with the field
sdgp * L zgp, L the Cholesky factor of the kernel over a farthest-first
order of the field's distinct positions, taken on the tape at a cost of
n^3 / 3 per gradient for n positions, at a nugget of 1e-12 * sdgp^2.
Drawing that field does not converge at sigma near 1e-3 and below on a
small n: at n 15, Rhat on sdgp read 1.04 to 2.37 at 1e-3, 1.36 to 1.95 at
1e-4 and 1.56 to 2.47 at 1e-6, and sigma's median up to 394 times the
truth at 1e-6. brms's program does not mix there either.
laplace = TRUE does, Rhat at most 1.02.
What stays centered is a run whose variance parameters carry
no density at all: check_laplace(), which asks for that on purpose.
Under prior = "flat" a standard deviation is flat on its own scale,
as brms's (flat) is, which is a density on theta that vanishes
at sd = 0, so such a run non-centers too.
The call message()s every block it left centered, with the reason.
A model made only of those samples centered throughout and says so
rather than failing. The reasons:
NO DENSITY on the block's log standard deviation. Send
sddown withzwhere it is:b = sd zgoes to zero, the likelihood settles on the model without that random effect, and the density stops changing: a flat tail with nothing to stop a chain in it. Measured on a six-group random-intercept model, one chain of 2000, three seeds: under flat priors the non-centered chain walksthetato -1e15 at a bulk-ESS of 1; under the defaultstudent_t(3, 0, s), which makes that tail integrable, 174 to 284 against 3 to 48 centered.NO DENSITY on the block's CORRELATION, which is the same argument: flat on frmtmb's unbounded correlation parameter is
(1 - rho^2)^-3/2, improper, with all its mass at|rho| = 1. Before 0.39 that kept every correlated block centered, because no proper correlation prior existed to close the tail;lkj(eta)now does, and it is a default, so(Days | Subject),cs(),ar1()andhetar1()non-center like any other block.A Student-t latent (
gr(dist = "student")): a scale mixture, not a linear factor.car(),spde()andgr(prec = ): sparse precisions whose factor is dense.The spatial covariances (
ou,exp,gau,mat): a full factorization per gradient evaluation.toep(): not positive definite everywhere in its parameterization, sob = L zis not a bijection there.
What it is worth. One chain of 2000 iterations, three seeds
(dev/benchmarks.md). Where each group's own data say little (the
regime the funnel lives in) it is decisive: 80 groups of 2 binary
observations run at a min-ESS of 236 against 5 centered, 55 times the
effective draws per second. Where the groups are informative it is a
wash: the epilepsy GLMM and an uncorrelated sleepstudy are within
noise of the centered chain either way. It is not a blanket speed-up,
and the blocks it declines to touch are the ones where it would have
done harm.
laplace = TRUE integrates the random effects out, which removes the
funnel by itself, so that route ignores reparameterize entirely.
The ML fit is untouched either way. The Laplace approximation
integrates b out, and that integral is invariant under a linear
change of the integrated variable, so frm() has nothing to change
and no fitted object ever carries a non-centered tape.
Laplace draws
laplace = TRUE passes to tmbstan::tmbstan(), which then samples
the outer parameters alone (the coefficients, the distributional
parameters and the covariance parameters) and integrates the rest
out with the Laplace approximation at every step: the group-level
coefficients, the coefficients of a smooth, and the values of an
mi() term. The draws hold no columns for those. brms has no such
route, because its sampler always visits the group-level
coefficients.
What the draws answer. Every summary of the outer parameters:
summary(), fixef(), VarCorr(), posterior_summary(), the
as_draws_*() converters, hypothesis() on the coefficients and on
sd_ and cor_ quantities, and mcmc_plot(). Also every predictive
quantity that reads no integrated value: posterior_epred(),
posterior_linpred() and posterior_predict() with
re_formula = NA on a model whose only integrated values are
group-level coefficients, a distributional parameter with no
group-level term, and the default curves of conditional_effects().
These are computed from the draws as they are and are exact.
What they refuse, by name. A quantity that reads an integrated
value: a prediction at the observed group levels (the default
re_formula, with or without newdata), any prediction on a model
with a smooth or an mi() term that reads them, ranef(),
coef(), log_lik() and everything built on it (loo(), waic(),
psis(), the loo_*() functions, pp_check()'s loo_* types), and
bayes_R2(). The
refusal names the function you called and says what to do instead;
it suggests re_formula = NA only where that function takes it and
it computes on the model. Whether a call reads an integrated value is
found out by evaluating it at one draw with those values set to NA
and to 0: a result that changes read them. A draw whose result has
non-finite values where the first draw's has none is refused the same
way, so a value that vanishes at the first draw only cannot slip
through as NaN.
Why they refuse rather than fill the values in. laplace = TRUE
samples an approximate posterior, and its use is to check that
approximation against draws that sample the group effects. For
predictions, and for anything else that conditions on the group
effects, sample without laplace = TRUE. Filling the values in would
also be wrong or random: a prediction at each draw's conditional
modes leaves out the uncertainty of the group effects given the
outer parameters, and understated the spread of the in-sample
posterior_epred() by up to a factor of four on a gaussian model
with 12 groups of 10 rows; drawing them from their Laplace
conditional matched the full draws there, but it would make
posterior_epred() depend on the random seed and disagree with
log_lik() and ranef() about what a draw contains.
Two edge cases. On a model with nothing to integrate out,
laplace = TRUE is the model itself, so the call says so and samples
it on the full route, with draws identical to a call without the
argument. On a REML fit whose objective is sampled as it stands
(check_laplace()), the objective integrates the coefficients out
too, and that combination is refused: the draws would hold no
coefficient.
Default priors
Both routes default to brms 2.23's own weakly-informative
priors, read off brms::default_prior() on matched models, family by
family. Write y* for the response transformed by a link and
s = max(2.5, round(mad(y*), 1)):
| class | default | scale |
b (slopes) | flat | - |
Intercept | student_t(3, round(median(y*), 1) - m, s) | link |
Intercept (ordinal thresholds) | student_t(3, 0, 2.5) | thresholds |
delta (equidistant thresholds) | (flat), lb = 0 where ordered | natural, log-Jacobian applied |
sd | student_t(3, 0, s) | natural sd, log-Jacobian applied |
meanme (latent mean of a me() term) | (flat) | natural |
sdme (latent sd of a me() term) | (flat), lb = 0 | natural sd, log-Jacobian applied |
cor, rescor, lncor, cortime, corme | lkj(1) | correlation matrix |
ar, ma (cov = TRUE), cosy, lagsar | flat, see below | natural |
simo, sbhaz | dirichlet(1) | the simplex |
theta (mixture weights) | dirichlet(1) | the simplex |
sigma | student_t(3, 0, s) | natural |
shape, phi | gamma(0.01, 0.01) | natural |
shape of the negative binomials | inv_gamma(0.4, 0.3) | natural |
sigma of brmsfamily("negbinomial2"), 1 / shape | gamma(0.4, 0.3) | natural |
nu | gamma(2, 0.1) | natural |
kappa | gamma(2, 0.01); gamma(0.01, 0.01) for xbeta() | natural |
beta (exgaussian) | gamma(1, 0.1) | natural |
zi, hu, zoi, coi, bias, quantile | beta(1, 1) | natural |
bs | gamma(1, 1) | natural |
ndt | uniform(0, min(y)) | natural |
alpha | normal(0, 4) | natural |
xi | normal(0, 2.5) on tmp_xi | see below |
The rows marked natural are for a parameter without a formula: the
density is on the parameter itself, with the Jacobian of its link,
which is what brms's bounds do on its constrained parameters. A
parameter with a formula, shape ~ x or shape ~ 1, has its default
on its intercept instead, on the link scale, as in brms:
student_t(3, 0, s) for sigma (with y* on sigma's link),
student_t(3, 0, 2.5) for shape and phi, normal(2.7, 0.8) for
nu, normal(5.0, 0.8) for kappa, normal(1.7, 1.3) for beta,
logistic(0, 1) for zi, hu, zoi, coi, bias, quantile and
a mixture weight, normal(-0.6, 1.3) for bs, normal(0, 4) for
xi and alpha, and normal(0, 1) for an ordinal disc
(lognormal(0, 1) under its identity link). ndt with a formula
stays flat, as in brms. A family brms does not have takes the
default of the dpar's name, which is what brms gives a custom family.
For sigma without a formula y* is the response itself (logged for
a log-scale family), whatever the mu link: brms reads sigma's
own link there, which is the identity. com_poisson()'s dispersion,
frmtmb's nu, takes brms's shape default, because it is the same
quantity under brms's name. A mixture component's dpars take the
class's default (sigma2 is class sigma) and no family's own: brms
asks the mixture, which states none.
The support of each natural default is brms's bound on the class:
frmtmb's link maps onto it (log onto (0, Inf), logm1 onto
(1, Inf) for nu, logit onto (0, 1)). Two classes differ, and
the default is still the faithful one. ndt of
frmtmb::shifted_lognormal() is on a log link, and its upper bound,
brms's min_Y, is the density's: the uniform is NaN past it, where
the likelihood has no support either. ndt of frmtmb.eam's
wiener() is on a link bounded by the fastest response, so the
uniform is flat over exactly brms's interval. A scalar xi of
frmtmb::gen_extreme_value() is sampled through brms's tmp_xi (see
the message), and is above -1, where brms's range goes on.
A residual-correlation coefficient that frmtmb holds through a
transform gets brms's (flat) on the coefficient itself, written as
a uniform density over the whole interval the transform maps onto and
carried inward with its Jacobian: uniform(-1, 1) on a first-order
cov = TRUE ar or ma, a uniform on frmtmb.spatial's sar()
interval for lagsar and errorsar, and for cosy a uniform on its
positive-definite window with lb = 0, which is brms's [0, 1]. A
cov = TRUE term of order two or more, which brms does not have,
takes a uniform on a box that holds its stationary (or invertible)
region, which is flat on that region. A prior of your own on one of
these classes keeps the default's bound on each side it leaves
unset, and a bounds-only one keeps the flat density, as brms keeps
both. The cov = FALSE coefficients are unconstrained and are the
sampled parameters, so brms's flat needs no density there. A bounded
one of these parameters is sampled through one logistic onto its
interval, as Stan's bound is in brms.
Under a cov = TRUE ar(), ma() or arma(), sigma is the
innovation sd, as in brms (frmtmb::frmtmb-autocor), in these draws
as in frmtmb::frm(). So the default student_t(3, 0, s) on class
"sigma", a prior of your own with its lb and ub, and under
sigma ~ x the default and your prior on the intercept of sigma
(the log innovation sd), are brms's priors on the same parameter,
and the draws of sigma are brms's.
Every default with a bound (nu above 1 on an identity link, sigma
above 0 on one, ndt below the smallest response, cosy above 0)
keeps that bound under a prior of your own on the same slot, on each
side the prior leaves unset, as brms keeps the bound it declares.
A mixture's weights without a formula are a simplex, and
dirichlet(1) on it is the flat density there, with the Jacobian of
the log ratios the sampler works on. No frmtmb::set_prior() class
names them (theta is frmtmb's covariance vector), so the default
is disclosed in the message, and prior = "flat" keeps it: it is
brms's (flat) on the weights.
A mo() simplex is sampled in softmax coordinates under
dirichlet(1), brms's default and the flat density on the simplex,
which prior = "flat" therefore keeps; the fit's own coordinates
(see "Monotonic effects" in frmtmb::frm()) have folds at the faces
that wall a density off into separate sheets. The draws carry brms's
names, simo_<term>1[k], one column per weight. A simo row in
prior =, set_prior("dirichlet(c(5, 1, 1))", class = "simo", coef = "mox1"), replaces dirichlet(1) on that simplex, as in brms,
with the softmax's Jacobian applied.
The disc default matters more than most: the likelihood cannot
tell an intercept in disc apart from the scale of the thresholds,
so this prior is what places it.
m is the mean of the predictor's offset(), 0 without one. brms
takes it off because an offset moves the intercept by that much, and
only under the links below that transform the response; under any
other link the location stays 0.
The sd default is written once per response, distributional
parameter and nonlinear parameter that has a random-effect standard
deviation, as brms writes one row for each, with y* on the mu
link. A class "sd" specification reaches only its own response,
distributional parameter and nonlinear parameter (see
frmtmb::set_prior()), so a set_prior(class = "sd") of your own
replaces the location parameter's default and leaves, for example, a
phi block's default in place, as in brms.
A multivariate model reads y* and s off each response for that
response's rows, which are written with resp, as brms's are.
The link is transformed only for identity, log, inverse,
sqrt and 1/mu^2 - brms's own list - with a log-scale family
(frmtmb::lognormal() and its relatives, whose mu link is spelled
identity but whose response is logged inside the density) counted
as log. Under any other link (logit, cloglog, ...), for a
response that was a factor, and for a family whose parameters run
over categories (ordinal, categorical, multinomial), the location is
0 and the scale 2.5 - which is what brms reports for a bernoulli
model. Under a log, inverse or 1/mu^2 link the ZEROS of the
response are replaced by 0.1 before the transform, element by
element, exactly as brms does; the non-zero values are left alone,
so a count response with median 1 keeps log(1) = 0 as its
location.
The call message()s the priors it chose, one compact line per
class, AND one line per slot it deliberately left flat, so a model
that gets few defaults - or none - says so rather than looking flat
by accident. Wrap it in suppressMessages() to silence that.
prior_summary() on the returned draws reproduces the chosen priors
exactly, except the mixture weights' Dirichlet, which no class names.
What is deliberately NOT matched. Each of these is named in the message whenever the model has one.
A
toep()correlation stays flat. Its banded parameterization is not positive definite everywhere, so there is no correlation matrix over its whole parameter space for an LKJ density to be about. Every other correlated structure getslkj(1), brms's own default (seefrmtmb::set_prior()for what the density is on frmtmb's parameters).Sum-to-zero ordinal thresholds stay flat: brms puts its density on thresholds it declares before centering them, a vector frmtmb does not have.
A Student-t group term's degrees of freedom have no prior: frmtmb fixes them (
gr(dist_nu = )), where brms estimates them undergamma(2, 0.1).NONLINEAR parameters stay flat, which is what brms does: their coefficients are class
b, and the response's median and mad say nothing about a rate or a shape sitting inside a nonlinear body. Write them withset_prior(nlpar = ), as the brms nonlinear vignette does.
Overriding and opting out. A set_prior() specification takes over
the parameters it reaches. The defaults stay on every other
parameter, as in brms. set_prior(class = "sd", group = "g")
replaces the default on the standard deviations of g and keeps it
on the other groups. Adding coef = "x" narrows the prior to one
standard deviation. A default written per coefficient, such as the
lscale of each gp() term, is more specific than a class-wide
specification and stays, as in brms, and so does a more specific
specification of your own. The default-prior message names a
specification that applies to no parameter for either reason.
Write the coef to replace such a default. The named list of prior
objects keyed by internal parameter names, which earlier releases
took, is refused: it gave a parameter a second prior slot beside its
class, and a name inside a correlation block dropped the LKJ density
on the rest of the block.
A missing value under mi(). A draw of it is named as in brms,
Ymi_<resp>[k] (Ymi[k] with one response, Yl_ under mi(sd)), and
is on the response's own scale, whatever chart the sampler moves in. A
frmtmb::frm() fit leaves out of its objective a value no density reads;
sampling from that fit holds it again. laplace = TRUE, from a fit or from a
formula, integrates the values out and leaves such a value out, as frm()
does.
prior = "flat". brms has no opt-out of its priors; the nearest
brms model is "" in every slot, which is (flat): no density on the
parameter brms declares, inside the bounds it declares. That is what
prior = "flat" samples. frmtmb holds several of those parameters
through a map, so flat on brms's parameter is not flat on frmtmb's:
a standard deviation (sd, sds, sdgp, lscale) and a positive
distributional parameter (sigma, shape, nu, ...) are flat on
their own scale, not on their log; a correlation matrix (cor,
rescor, cortime, lncor) is flat on the Cholesky factor brms
declares, which beyond two dimensions is not lkj(1); the ordered
thresholds and delta are flat on themselves; and ar, ma,
cosy, lagsar and the simplexes keep their defaults, which are
already flat on brms's parameter. A coefficient, and an intercept of
a distributional parameter with a formula, is the parameter brms
declares and gets no density. Through frmtmb.sample 0.23.0 the
opt-out added nothing, which is flat on the internal parameters: on
a log standard deviation it is improper at sd = 0, and on a
first-order correlation it is (1 - rho^2)^(-3/2). Measured against
brms with "" in every slot (4 chains of 3000 draws): on
y ~ 1 + (1 | g) with 6 groups the posterior mean of sd_g__Intercept
is 1.458 in brms, 1.455 here, and was 1.197; on (1 + x | g) with 8
groups the correlation's mean is -0.28 in brms, -0.26 here, and was
-1.00 with a posterior sd below 0.0001. A parameter brms has no class
for keeps no density. One parameter brms has a class for is not yet
brms's (flat): frmtmb.spatial's car() dependence parameter, which
brms declares on (0, 1) and leaves flat, carries no density here, so
it is flat on its logit, 1 / (rho (1 - rho)) on rho, which is not
integrable at rho = 0 where a CAR likelihood stays finite: under
prior = "flat" such a posterior is improper. Give it a prior with
class "theta" on its coef, as frmtmb::default_prior() lists it.
The call warns when the model has standard
deviations: nothing but the likelihood bounds them above, and with
few groups the posterior need not be proper, which neither the
chains nor Rhat can see. On a MAP fit "flat" reaches further than
the defaults: that fit's own prior is taped into its objective, so
the penalty is rebuilt out and the call says which prior it dropped.
A bounds-only specification of your own that replaces a default,
set_prior("", class = "sd", lb = 0.1), is brms's (flat) on that
slot in the same way: flat on the standard deviation above 0.1.
Three sources, most explicit first. On a fitted model the stack is
this call's prior =, then the fit's own prior (frm(prior = )
makes a MAP fit, and that prior is part of the model the fit is),
then the defaults. Each addressed slot is settled by the most
explicit source that names it, one slot at a time: a call-level
set_prior(class = "sd") replaces the fit's sd prior and the sd
default while the Intercept default stays. prior_summary() on
the returned draws prints what the stack came to.
Prior draws
sample_prior = "yes" draws from the priors as brms does. Each
prior statement of brms's Stan program for the model is one column
of draws, independent of the likelihood: a prior on a whole class is
one column (prior_b, prior_sd_g, prior_sigma); where a
coefficient has a prior of its own, each coefficient of that vector
is a column (prior_b_x). A flat prior, the default on class "b",
has no column. A bounded parameter is drawn from its prior truncated
to the bounds: sigma, sd and sds at 0, nu at 1, a "b" prior
at its lb and ub. A correlation is one off-diagonal element of an
LKJ draw (prior_cor_g, prior_rescor), and a mo() simplex or a
mixture's weights a Dirichlet draw. prior_Intercept is drawn from
the prior of the intercept at the centered predictors, which is the
intercept brms puts that prior on, so the population-level
intercept b_Intercept has no prior draws, as in brms. An ordinal
model with a threshold prior of its own (coef = "1") has one
column per threshold, prior_Intercept__1, ..., as brms does. A
gp() term's sdgp and lscale are not drawn, nor a GEV's xi
(brms draws no prior of it either). A me() term's hyperparameters
are drawn under brms's names, which are numbered by the group of
terms that share a gr =: prior_meanme_1 and prior_sdme_1 for
a class-wide prior, prior_Lme_1 for the correlation, and, where a
term has a prior of its own, one column per term spelled as brms's
renaming spells it, prior_meanme_1_meanme_mex. With brms's
defaults only prior_Lme_1 is drawn, as in brms.
The draws are made directly from each density with R's random number generators, as brms makes them in Stan's generated quantities: they are independent, and a run of the sampler on the prior would give correlated draws of the same distribution at the cost of a second run. The saved warmup iterations get prior draws too, as in brms.
sample_prior = "only" samples every parameter from its prior with
the likelihood left out (core's prior_only objective, documented
in frmtmb::frmtmb-sampling-api): the group-level coefficients
from their group's distribution, and the latent values of a
mi(sdy) response from their measurement density, as brms does.
The draws are draws from the prior, with the parameter names of the
model, and the predictive methods give the prior predictive
distribution. Every parameter must have a proper prior: a flat one
is refused with brms's message, which names the parameter as brms's
Stan program does ('b' with brms's default priors). brms refuses
only a parameter with no prior at all and samples a coefficient
left flat beside others that have a prior; this function refuses
that too. A model with missing values of an mi() term is refused,
because without the likelihood those values have no density (brms
samples them flat). A model with me() terms needs proper priors on
"meanme" and "sdme", as in brms, which names 'meanme_1'
otherwise; the measurement and latent densities of the me()
terms stay, as brms keeps them, so the latent values and their
hyperparameters are informed by the noisy values. Draws
made with "only" have no prior draws of their own, as in brms.
Warmup
rstan saves each chain's warmup iterations unless save_warmup = FALSE is passed, and so does brms, whose fit keeps them in its
stanfit. This function keeps them there too, and the
inc_warmup = TRUE of as_draws() and its relatives, as.mcmc(),
rhat() and neff_ratio() moves them to the draws' names and
scales and puts them before each chain's draws, as brms does.
Everything else reads the post-warmup draws only: ndraws(),
summary(), variables(), the predictions and the diagnostics.
nuts_params() and log_posterior() read the stanfit, and take
bayesplot's inc_warmup.
The pre-flight against the compatibility registry
Some models cannot be sampled at all, and the reason usually
belongs to the package that supplies the feature rather than to
this one. Such a package declares the pair through
frmtmb::frmtmb_register_compat() with the status "refused",
and this function reads those rows before it does anything else:
before the rstan and tmbstan namespaces are loaded, before the
formula route calls frmtmb::frm(), and before either route
retapes the objective for its priors. as_tmbstan() runs the same
check. The refusal repeats the registry's own note and names the
row, so frm_compat("frm_sample", "<feature>") gives the same
answer without running anything.
What it can decide, and what it will not guess at. Whether a model uses a registry feature is settled by WHERE a call sits in the formula. The check decides three ways.
A response FAMILY is matched by name, which is exact, because a
family is named and never called. It does not look inside
frmtmb::mixture(), whose family name is the whole string
"mixture(gaussian, gaussian)", so a refused gaussian row does
not reach a mixture of gaussians.
An ADDITION TERM is matched by POSITION. It is written inside the
bar on the left of the response formula and nowhere else, and the
bf() object carries that, so y | se(v) ~ x is separated from
y ~ a * se(x) with a user's own se(), which frmtmb fits.
Everything else the vocabulary writes as a call is matched BY
NAME, on the display name, outside addition-term position, and
only when a formula actually writes that name and it is not also a
base R function. This is the part that is justifiable from names
alone rather than sound: a name test cannot see that a call sits
where the feature cannot occur, so a refusal registered on mo(),
ma(), cosy() or mm() would also fire on a user's own
function of that name in a nonlinear body. Nothing is registered
on those today. mi_pred(), gp_pred() and cs_pred() are
spelled in the vocabulary differently from the way a formula
writes them, which the registry itself reports, so they warn
rather than matching a call no formula makes.
A refused row naming anything else, a covariance structure, a fitting mode, another post-fit method, WARNS and samples, and so does one whose name is a base R function. Refusing would fire on correct models; silence would make a registered refusal look like a passed check.
frmtmb.ode registers the one refused row there is: an
frm_ode() fit aborts at the first warmup step inside RTMBode,
and the abort leaves rstan's autodiff arena unusable for the rest
of the session, so the refusal is worth more than the attempt.
Multiply imputed data
Given a frmtmb::frm_multiple() result, frm_sample() samples the
fit of every imputation with the same arguments and combines the
draws, as brms's brm_multiple() combines the chains of its fits:
the draws of imputation 2 follow those of imputation 1, and so on,
and the stanfit holds every chain, so nchains() is the number
of imputations times chains. A seed gives imputation j the
seed seed + j - 1, so the result is reproducible and no two
imputations share a random stream. The draws object keeps the
first imputation's fit, so a method that reads the data in sample
(posterior_predict() without newdata, pp_check(), loo())
reads the first imputed data set, which is what brms does and warns
about. Chains of different imputations sample different
posteriors, so an R-hat across them measures the imputations as
well as the sampler, as the brms vignette on missing values says;
check each imputation's chains on their own.
Multimodal posteriors
For frmtmb::mixture() fits the posterior is multimodal by construction
(label switching at minimum). Mode-centered inits, jittered or
not, leave every chain in one symmetry branch, so Rhat cannot
flag the others; use init = "random" there and inspect chains
individually.
Examples
# \donttest{
if (requireNamespace("tmbstan", quietly = TRUE) &&
requireNamespace("rstan", quietly = TRUE) &&
!frmtmb.sample:::tmbstan_build_broken()) {
set.seed(9)
dd <- data.frame(x = rnorm(80), g = factor(rep(1:8, 10)))
dd$y <- rnorm(80, 1 + 0.5 * dd$x + rnorm(8, 0, 0.5)[dd$g], 1)
fit <- frm(bf(y ~ x + (1 | g)) + gaussian(), data = dd)
ds <- frm_sample(fit, chains = 1, iter = 500, refresh = 0)
summary(ds)
fixef(ds)
hypothesis(ds, "sd_g__Intercept^2 / (sd_g__Intercept^2 + sigma^2) = 0",
class = NULL)
# the same model sampled straight from the formula, with no ML fit.
# It reports the brms default priors it chose, and prior_summary()
# gives them back.
ds2 <- frm_sample(bf(y ~ x + (1 | g)), data = dd, family = gaussian(),
chains = 1, iter = 500, refresh = 0)
prior_summary(ds2)
fixef(ds2)
# a set_prior() specification takes over the parameters it reaches
# and leaves the defaults on the rest
ds3 <- frm_sample(bf(y ~ x + (1 | g)), data = dd, family = gaussian(),
chains = 1, iter = 500, refresh = 0,
prior = set_prior("exponential(1)", class = "sd"))
prior_summary(ds3)
}
#> frm_sample(): default priors (brms 2.23 defaults; prior = "flat" puts brms's (flat) in every slot)
#> Intercept student_t(3, 1, 2.5)
#> sigma student_t(3, 0, 2.5) [natural scale]
#> sd student_t(3, 0, 2.5) [natural sd scale]
#> b (flat), as brms leaves slopes
#> Warning: The largest R-hat is 1.1, indicating chains have not mixed.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#r-hat
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#tail-ess
#> frm_sample(): default priors (brms 2.23 defaults; prior = "flat" puts brms's (flat) in every slot)
#> Intercept student_t(3, 1, 2.5)
#> sigma student_t(3, 0, 2.5) [natural scale]
#> sd student_t(3, 0, 2.5) [natural sd scale]
#> b (flat), as brms leaves slopes
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#tail-ess
#> frm_sample(): default priors (brms 2.23 defaults; prior = "flat" puts brms's (flat) in every slot)
#> Intercept student_t(3, 1, 2.5)
#> sigma student_t(3, 0, 2.5) [natural scale]
#> b (flat), as brms leaves slopes
#> Warning: There were 1 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: Examine the pairs() plot to diagnose sampling problems
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#tail-ess
#> prior class coef group resp dpar nlpar lb ub source
#> student_t(3, 1, 2.5) Intercept <NA> <NA> user
#> student_t(3, 0, 2.5) sigma <NA> <NA> user
#> exponential(1) sd <NA> <NA> user
# }