Skip to contents

Diagnosing a sampled fit

Source: vignettes/posterior-diagnostics.Rmd
posterior-diagnostics.Rmd

vignette("diagnostics", package = "frmtmb") asks three questions of a maximum-likelihood fit: did the optimizer converge, does the family describe the data, and are the reported uncertainties trustworthy.

A sampled fit reorders that list. There is no optimizer to interrogate, so the first question becomes whether the sampler explored the posterior it was given. The third question changes character: the posterior IS the uncertainty statement, so the interesting comparison is against the maximum-likelihood fit that the same objective produced.

Every chunk here is shown and not run, because every one needs a sampler this package only suggests. vignette("sampling") says more about why, and its own prior tables do run.

library(frmtmb)
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)

Did the sampler explore the posterior?

summary(ds) puts Rhat, Bulk_ESS and Tail_ESS beside each estimate, which is where to start. Read the pattern and not only the worst entry. One parameter with a low effective sample size and a clean Rhat is a slow chain. Several at once usually means the geometry is wrong rather than the run too short, and sampling longer will not fix it. A high Bulk_ESS with a low Tail_ESS says the centre of the posterior is well explored and the extreme quantiles are not, so an interval is less trustworthy than a mean.

summary(ds)
rhat(ds)
neff_ratio(ds)

These are brms’s diagnostics throughout, computed by the posterior package on these draws and reported under brms’s names, which variables(ds) lists. rhat() is the rank-normalized split-R-hat, neff_ratio() is min(ess_bulk, ess_tail) / ndraws, and summary(ds) reports the same Rhat with the two effective sizes it takes the minimum of. The three columns and the two accessors are one set of numbers: summary(ds)[, "Rhat"] IS rhat(ds) on those rows.

None of them is what the sampler itself recorded. Stan’s own classic split-R-hat and n_eff are a different definition of the same idea, and they are in rstan::summary(ds$stanfit)$summary if you want to compare. Expect them to disagree by about the size of the excess over 1 that either R-hat reports.

The sampler’s own record is in nuts_params(), with the divergent transitions, the tree depths and the step size per iteration, and the log posterior is log_posterior(). Both feed bayesplot’s mcmc_nuts_*() displays.

np <- nuts_params(ds)
sum(np$Value[np$Parameter == "divergent__"])   # zero is the target

mcmc_plot(ds, type = "trace")
mcmc_plot(ds, type = "rank_overlay")
pairs(ds, variable = c("theta_1", "r_g[1,Intercept]"))

Divergences concentrate in the funnel of a random-effect block, so the pairs() display to draw is a group’s log standard deviation against something it multiplies.

Three moves, in the order worth trying:

  1. Give the variance parameters a prior. reparameterize = TRUE is the default, but a block qualifies for the non-centered form only when every one of its variance parameters carries a prior, and the call reports each block it left centered and why. Under prior = "flat" a standard deviation is flat on its own scale, as brms’s (flat) is, which closes the tail at sd = 0, so a block qualifies there too. prior = set_prior("exponential(1)", class = "sd") is the usual fix, and it fixes two things at once: it closes the flat tail at sd = 0 and it unlocks the reparameterization that removes the funnel.
  2. Sample longer. iter = and chains = pass straight to tmbstan. A low Bulk_ESS with Rhat near 1 is answered here.
  3. Tighten the sampler. frm_sample(fit, control = list(adapt_delta = 0.99)) is brms’s spelling and brms’s meaning: the list goes to rstan through tmbstan. Raising adapt_delta shortens the step size the adaptation settles on, which is what removes the divergences a residual funnel produces, and the cost is wall-clock time. max_treedepth = 12 answers the companion warning. Read the step size back with nuts_params(ds) to see what the change bought. The fit-time options of frmtmb_control() are fit_control in this function.

One check comes before all of these. A tmbstan built against StanHeaders 2.39 samples a standard normal instead of the model, and every diagnostic on this page looks perfect when it does. frm_sample() refuses such a build before sampling; the sampling vignette has the detail.

Does the model describe the data?

pp_check() overlays the observed response with responses drawn from the posterior predictive distribution. This is the check brms users already run, on the same bayesplot displays, and it is a stronger statement than the maximum-likelihood version in core, because the draws carry parameter uncertainty into the simulated responses.

pp_check(ds)
pp_check(ds, type = "stat", stat = "sd")

loo() is the real leave-one-out here, not the AIC() substitute core offers. It runs Pareto-smoothed importance sampling on the log_lik() matrix and returns the loo package’s own object, so print() and loo::pareto_k_table() work on it unchanged.

loo(ds)

# a comparison is two sampled models, each with its own elpd
ds0 <- frm_sample(bf(y ~ 1 + (1 | g)) + gaussian(), data = dd, chains = 4)
loo_compare(ds, ds0)

A high Pareto k is a real finding here, because the refit-based remedies are absent: reloo(), kfold(), loo_moment_match() and loo_subsample() all refuse, and say so, since each needs to refit the model on modified data. What to do instead depends on why k is high. Under prior = "flat" the elpd is unregularized and a model with many group-level parameters will produce those warnings by construction, so the first move is to sample under the default priors and read the number again. If the warnings survive that, the observations loo names are genuinely influential: compare the models by AIC() on the maximum-likelihood fits, or measure the influence directly with frmtmb::frm_bootstrap().

Read every one of these numbers as a posterior quantity. Both of frm_sample()’s routes apply brms’s default priors, so an elpd is regularized the way brms regularizes it unless the call opted out.

Does the Laplace approximation hold?

frmtmb::frm() integrates random effects out with a Laplace approximation and reports Wald standard errors. Both assumptions degrade in the same places: variance components from few groups, binary data in small clusters. Core routes around that with confint(method = "profile"), frm_bootstrap() and frm(quadrature = TRUE). This package measures it instead.

check_laplace() samples the density the fit maximized, and no other. It hands Stan fit$obj as it stands, which is why it opts out of frm_sample()’s default priors explicitly rather than inheriting them: a prior nothing in the fit ever saw would change the very thing being measured. For an ordinary fit that objective is the likelihood, so the run is flat. For a MAP fit (frm(prior = )) the penalty is taped into the objective, so the run carries it, which is right, because the mode and the Wald standard errors it compares against come from the same penalized Hessian.

cl <- check_laplace(fit, chains = 2, iter = 1000)
cl

One row per outer parameter, with the maximum-likelihood estimate (ml), the posterior mean (post_mean), the two uncertainty statements (wald_se, post_sd), and the two comparisons worth reading:

  • z_shift is (post_mean - ml) / post_sd, the distance from the mode to the posterior mean in posterior standard deviations. It measures asymmetry the quadratic approximation cannot see.
  • sd_ratio is post_sd / wald_se. Above one, the Wald interval is too narrow.
# the parameters whose Wald interval to replace with a profile or
# bootstrap one
cl[abs(cl$z_shift) > 0.3 | cl$sd_ratio > 1.3, ]

Those thresholds are working values, not a test. A variance component estimated from eight groups will fail them and should: the log standard deviation has a skewed likelihood at that sample size, and the skew is the finding.

Because the defaults are off, this function samples the centered parameterization on an ordinary fit. That costs it nothing. The non-centered form is offered only to blocks whose variance parameters carry priors, and the flat prior that makes this comparison meaningful is exactly the one that leaves a tail at sd = 0 for a non-centered chain to walk into. Give the variance parameters a prior through prior = and the run does non-center, but it is then measuring the Laplace approximation of a different posterior, which is usually not the question.

Which page to read

question where
did the optimizer converge vignette("diagnostics", package = "frmtmb")
does the family describe the data (ML) same page: dharma_residuals(), residuals(type = "osa")
did the chains converge here
does the model describe the data (posterior) here: pp_check(), loo()
is the Laplace approximation good enough here: check_laplace()

vignette("sampling", package = "frmtmb.sample") is the page on getting the draws in the first place: the two routes, the default priors, and the non-centered parameterization.