frmtmb fits by maximum likelihood. This package hands the same objective to Stan’s NUTS sampler through tmbstan and gives back the posterior, under brms’s default priors and with frmtmb’s own parameter names.
Nothing about the fit changes. The model is the one
frm() already assembled; what is added is a posterior, and
a starting point.
Most chunks on this page are shown and not run, and speed is not the
reason: four chains on the model below take seconds, not minutes. The
reason is that tmbstan and rstan are Suggests, so a build
machine need not have them, and that a tmbstan built against the wrong
StanHeaders samples the wrong density in silence (the last section says
how). Numbers that depend on which Stan the page was built against are
worse than no numbers. The two chunks that need no sampler do run: the
setup below, and the get_prior() calls under The
default priors.
library(frmtmb)
#>
#> Attaching package: 'frmtmb'
#> The following object is masked from 'package:stats':
#>
#> ar
library(frmtmb.sample)
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 = 4)
summary(ds)Two routes, one posterior
frm_sample() takes either a fit or a formula.
# from a fit: the ML mode anchors the chains and shortens warmup
ds1 <- frm_sample(fit, chains = 4)
# from a formula: assembled and taped, never optimized, so Stan's own
# overdispersed initialization is the default
ds2 <- frm_sample(bf(y ~ x + (1 | g)), data = dd, family = gaussian(),
chains = 4)Both sample the same posterior under the same default priors. What the fit adds is a starting point, not a different density.
The default priors
Both routes apply brms 2.23’s weakly-informative defaults, read off
brms::default_prior() on matched models. The call reports
what it chose, one line per class, and one line per slot it deliberately
left flat, so a model that gets few defaults says so rather than looking
flat by accident.
frmtmb::get_prior() reports them too, when you ask for
this route. It has two routes, because frmtmb has two:
route = "fit" is the default and reports what
frm() applies, which is flat in every slot.
route = "sample" reports what frm_sample()
applies, and this package is what lets it answer.
# what frm() applies. Loading this package does not change it
get_prior(bf(y ~ x + (1 | g)) + gaussian(), data = dd)
#> route = "fit": the prior defaults frm() applies
#> prior class coef group resp dpar nlpar lb ub source
#> 1 (flat) b NA NA default
#> 2 (flat) b x NA NA default
#> 3 (flat) Intercept NA NA default
#> 4 (flat) sd NA NA default
#> 5 (flat) sd g NA NA default
#> 6 (flat) sd Intercept g NA NA default
#> 7 (flat) sigma NA NA default
# what frm_sample() applies, which is the brms reading of get_prior().
# Without this package loaded, the call is refused rather than
# answered "(flat)"
get_prior(bf(y ~ x + (1 | g)) + gaussian(), data = dd,
route = "sample")
#> route = "sample": the prior defaults frm_sample() applies
#> prior class coef group resp dpar nlpar lb ub source
#> 1 (flat) b NA NA default
#> 2 (flat) b x NA NA default
#> 3 student_t(3, 1, 2.5) Intercept NA NA default
#> 4 student_t(3, 0, 2.5) sd NA NA default
#> 5 (flat) sd g NA NA default
#> 6 (flat) sd Intercept g NA NA default
#> 7 student_t(3, 0, 2.5) sigma NA NA defaultThe printed table names its route on the first line, so a table copied out of a session still says which set of defaults it describes.
prior = "flat" opts out of the defaults and puts brms’s
(flat) in every slot: no density on the parameter brms
declares, inside the bounds brms declares. That is the model brms fits
with "" in every row of its prior table. A standard
deviation is then flat on its own scale and a correlation matrix flat on
its Cholesky factor, not on the log or the unbounded parameter that
frmtmb holds. It is fragile inference all the same: nothing but the
likelihood bounds a standard deviation above, so with few groups the
posterior need not be proper, and neither the chains nor Rhat can see
it. The opt-out warns when the model has standard deviations.
A set_prior() specification takes over the parameters it
reaches and leaves the defaults on the other parameters, as brms does. A
prior on one group’s standard deviations, or on one coefficient’s
(set_prior(class = "sd", group = "g", coef = "x")), keeps
the default on the rest.
ds <- frm_sample(fit, chains = 4,
prior = set_prior("exponential(1)", class = "sd"))
prior_summary(ds)Non-centering
A random-effect block has a funnel in its centered joint posterior:
the width of the prior on b is a standard deviation 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.
reparameterize = TRUE, the default, samples
z ~ N(0, I) and computes b = L(theta) z on the
tape. Each draw is mapped back through its own theta, so
every method downstream reads the same b[i] columns and
cannot tell the two routes apart.
Not every block qualifies, and the call says which did not and why. A
block is non-centered only when every parameter it has is a standard
deviation or a correlation with a registered Cholesky factor, and every
one of those parameters carries a prior. Both conditions are about not
handing the sampler a direction it can run away in: removing the funnel
without closing off the flat tail at sd = 0 first trades a
slow chain for a wrong one.
See ?frm_sample for which structures qualify, what it is
worth, and the measurements behind that.
Once the draws exist
Whether the chains explored the posterior, whether the model
describes the data, and whether the Laplace approximation of the
maximum-likelihood fit held in the first place are three separate
questions, and vignette("posterior-diagnostics") works
through them with the tools this package provides for each:
rhat(), neff_ratio() and
nuts_params() for the sampler, pp_check() and
loo() for the model, and check_laplace() for
the approximation.
The method surface
The draws object supports the post-processing surface a brmsfit does, on the same generics:
fixef(ds); ranef(ds); VarCorr(ds)
posterior_epred(ds); posterior_predict(ds); posterior_linpred(ds)
hypothesis(ds, "sd_g__Intercept^2 / (sd_g__Intercept^2 + sigma^2) = 0",
class = NULL)
conditional_effects(ds)
pp_check(ds)
log_lik(ds); loo(ds); waic(ds); bayes_R2(ds)
as_draws(ds); as.array(ds); as.mcmc(ds)
mcmc_plot(ds, type = "trace")Parameter names are brms’s: b_Intercept,
b_x, r_g[1,Intercept],
sd_g__Intercept. The log standard deviation the sampler
samples, which brms does not, keeps the name confint()
gives it on the fit, theta_1, and
sd_g__Intercept is computed from it.
variables(ds) lists them all.
Conditional effects draw the frame the fit draws
conditional_effects(ds) returns what
conditional_effects(fit) returns: the same columns in the
same order, the same grid, the same keys, and the same attributes brms’s
plot() reads. What differs is where the band comes from.
The curves are posterior expected-response draws, so
estimate__ is their median (their mean with
robust = FALSE, brms’s switch) and lower__ and
upper__ are their percentiles, and
attr(df, "band") reads "posterior".
The default display is the EXPECTED RESPONSE, not the mu
predictor. On a zero-inflated or hurdle family, a mixture, or a
trunc() response those are different curves, and the one
you want is almost always the mean. Name a dpar = to get a
predictor instead.
Three arguments mean here what they mean on a fit.
int_conditions = labels a moderator,
allow_new_levels = is accepted in both spellings, and
re_formula = NULL conditions on a NEW group: the grouping
column of the returned frame reports NA, and each posterior
draw carries that group’s effects drawn from its own covariance
parameters. The band is wider than the population band by what that
group’s variance adds, which is the whole point of asking.
seed = makes the draw reproducible.
conditional_effects(ds, effects = "x") # population
conditional_effects(ds, effects = "x", re_formula = NULL, seed = 1)
conditional_effects(ds, effects = "x:z",
int_conditions = list(z = c(low = -1, high = 1)))
conditional_effects(ds, categorical = FALSE) # ordinal: E[category]There is no method = and no band =. The
curves ARE the posterior, so there is no wald/profile/boot choice to
make. For a predictive band, quantile posterior_predict()
over a grid you build yourself.
Evidence ratios and Bayes factors
hypothesis() returns brms’s object, and its
hypothesis frame reports brms’s Evid.Ratio and
Post.Prob beside the interval. They answer two different
questions depending on how the hypothesis is written.
A DIRECTIONAL hypothesis is the simple one. Its evidence ratio is the posterior odds of the claim, counted from the draws. No prior enters it, so it is always available.
h <- hypothesis(ds, "x > 0")
h$hypothesis$Evid.Ratio # P(x > 0) / P(x < 0)
h$hypothesis$Post.Prob # P(x > 0)A POINT hypothesis is the Savage-Dickey density ratio: the posterior
density of the tested quantity at the point, divided by its prior
density there. For a point null nested inside the model, that ratio IS
the Bayes factor. As in brms, the prior density is read from draws of
the prior, so sample with sample_prior = "yes". The
parameter needs a proper prior, and frm_sample() leaves
class "b" flat by default because brms does, so write
one:
ds <- frm_sample(fit, prior = set_prior("normal(0, 1)", class = "b"),
sample_prior = "yes")
hypothesis(ds, "x = 0")$hypothesis$Evid.Ratio
# draws made without prior draws get them with update()
ds <- update(ds, sample_prior = "yes")Without prior draws for every parameter the hypothesis names, the
ratio is NA, as in brms. That is the case for draws made
without sample_prior = "yes", for a parameter with a flat
prior, and for the overall intercept b_Intercept, whose
prior is written about the intercept at the centered predictors (brms
warns about that one). The other rows of the same call still report.
Both densities are kernel densities, with
stats::density()‘s default nrd0 bandwidth over
4096 grid points, read at the point with a spline: brms’s
evidence_ratio(), so the two packages’ numbers are the same
quantity. prior_draws() returns the prior draws, under
brms’s names, and the hypothesis object carries the prior draws of each
hypothesis as prior_samples.
A density ratio deserves an error bar. The between-chain Monte Carlo error of each ratio’s numerator rides on the returned object:
h <- hypothesis(ds, "x = 0")
attr(h, "evid_ratio_mcse")The prior predictive distribution
sample_prior = "only" samples the prior of the whole
model with the likelihood left out, as brms does. Every parameter needs
a proper prior; a flat one is refused with brms’s message. The draws
then give the prior predictive distribution:
dp <- frm_sample(fit, prior = set_prior("normal(0, 1)", class = "b"),
sample_prior = "only")
yrep <- posterior_predict(dp)
pp_check(dp)Known failure: a tmbstan built against StanHeaders 2.39
stanc 2.39.0 emits two log_prob_impl overloads where
2.32 emitted one, and tmbstan’s install-time code generator patches only
the first. HMC then reads its value and gradient from the unpatched
overload, so every chain samples a standard normal instead of the model,
with the objective, priors and data all silently absent.
frm_sample() and as_tmbstan() check the
installed model source once per session and refuse such a build before
sampling, rather than returning plausible garbage after it. If you hit
that refusal, install a binary tmbstan build or reinstall tmbstan
against StanHeaders 2.32.10, and distrust any draws the affected
installation already produced.