Predictions from a frmtmb fit

brms's predict(): a summary of the PREDICTIVE distribution of the response, with observation noise included, in the columns Estimate, Est.Error, Q2.5 and Q97.5. It is a much wider interval than fitted()'s, which summarizes the expected response.

A maximum-likelihood fit has no posterior, so the draws are simulated. Each replicate draws the outer parameter vector from its asymptotic normal law, N(theta_hat, vcov(object, full = TRUE)), evaluates every distributional parameter there, and draws a response from the family's own simulator, the one simulate() uses. The interval is the empirical quantile of those draws, so it carries the family's skew and its discreteness: a poisson predictive interval is on the counts.

Usage

# S3 method for class 'frmtmb_fit'
predict(
  object,
  newdata = NULL,
  re_formula = NULL,
  transform = NULL,
  resp = NULL,
  negative_rt = FALSE,
  ndraws = NULL,
  draw_ids = NULL,
  sort = FALSE,
  ntrys = 100L,
  cores = NULL,
  summary = TRUE,
  robust = FALSE,
  probs = c(0.025, 0.975),
  ...,
  allow_new_levels = FALSE,
  sample_new_levels = NULL,
  propagate_error = TRUE
)

Arguments

object: A frmtmb_fit.
newdata: Optional data frame to predict on. Defaults to the training data.
re_formula: Which group-level terms enter the prediction. NULL (default) keeps all of them, NA keeps none. A one-sided formula keeps the terms it names and drops the others, as in brms: ~ (1 | g) on a fit with (1 + x | g) + (1 | h) keeps the intercept of g alone. A formula that names no group-level term, such as ~0 or ~1, is NA. A term the fit does not have is an error that names it; brms drops such a term silently. See frm_linpred() for the full rule.
transform: A function applied to the draws before they are summarized, as in brms.
resp: For multivariate fits: which response, or NULL (default) for all of them, which is brms's nrow x 4 x nresp array with the responses named on the third dimension.
negative_rt: brms's sign convention for its wiener family. Refused, with the reason: this package's evidence-accumulation families return the response and the time as they declare them.
ndraws: Number of simulated replicates. Defaults to 1000. brms thins its posterior with this; here it sets how many draws are taken, so a larger number narrows the Monte Carlo error of the quantile columns and costs proportionally more.
draw_ids, sort, cores: Refused by name. There are no stored draws to index, rows always come back in the order of the data, and the simulation is not parallelized.
ntrys: Rejection-sampling attempts per row for a truncated response, brms's spelling of the family simulator's max_iter.
summary: If FALSE, the ndraws x nrow matrix of simulated draws instead of the summary.
robust: If TRUE, the median and the median absolute deviation of the draws instead of the mean and the standard deviation. brms's argument, and it is answerable here because these draws exist.
probs: Probabilities of the quantile columns.
...: Refused. An argument this method does not have is an error naming it, and type, se.fit, dpar and scale are refused with the name of the function that took over each one.
allow_new_levels: Predict rows whose grouping-factor level the fit never saw, instead of erroring. Each replicate draws that level's effect from the block's estimated covariance.
sample_new_levels: brms's argument. "gaussian" (or NULL) draws an unseen level's effect from its block's estimated covariance. "old_levels" gives each unseen level the effects of one seen level of its grouping factor, chosen at random once per call as brms chooses it, and every term of that factor reads the same seen level (within the row's by-level for gr(by = )), from each replicate's draw of the group effects. One choice serves every response of a multivariate call. Two different unseen values of one mm() term read two seen levels; brms 2.23.0 gives them one, because it numbers unseen values per member rather than by value. "uncertainty" mixes the two per posterior draw, and a maximum-likelihood fit has none, so it is refused by name.
propagate_error: Whether the error in the estimates is propagated into the interval. TRUE, the default, draws the parameters and the group effects of a level the fit saw. FALSE holds both at their estimates, so the interval carries the observation noise alone: the plug-in predictive distribution. It is the only way to include a group's own effect and still treat it as known; re_formula cannot say that, because it chooses which terms are in the prediction and not whether their error is carried.

Value

With summary = TRUE (the default) an nrow x 4 matrix with the columns Estimate, Est.Error and one per entry of probs; for an ordinal or categorical response an nrow x K matrix of simulated category proportions; for a multivariate fit an nrow x 4 x nresp array with the responses named. With summary = FALSE the ndraws x nrow matrix of draws, or ndraws x nrow x nresp for a multivariate fit.

Which uncertainty is in the interval

All three sources.

Observation noise, from the family's simulator. This is almost always the largest term.
Uncertainty in the estimates, from their asymptotic covariance.
Uncertainty in the group effects of a level the fit SAW. Each replicate draws the whole vector of group effects from its conditional law given the data and that replicate's parameters, which the joint precision of the fit describes. For a linear mixed model this is the prediction error variance of the best linear unbiased predictor, jointly with the fixed effects. brms's draws carry the posterior of each seen level's effect, and this is the frequentist analogue. Two rows of one group share that group's draw, so summary = FALSE is right across rows as well as per row.

propagate_error = FALSE drops the second and the third together and simulates at the estimates alone, which leaves the observation noise. Everything estimated is then treated as known.

What the interval covers: a new observation at a known level, with the nominal coverage AVERAGED over the groups. It is not a coverage statement for one group's realized effect. The predictor shrinks a group toward the population, so for a group far from the population the interval covers less than nominal and for a group near it more. dev/reunc-findings.md reports the measured coverage.

A new row of a me() term also has no latent value, so each replicate draws one from normal(x, sdx), per row, or per level under me(gr = ), as brms's posterior_predict() does on new data. The latent value of a new row is measurement noise in the grid, not an estimate, so propagate_error = FALSE keeps it. corme is not read, as in brms. Two responses or two dpars that read one me() term read ONE latent value, a stated divergence from brms, which draws one per predictor against its own Stan program: at the training rows of a two-response fit the residuals correlate at 0.734, one latent predicts 0.739 and brms's replicates give -0.001, and the 95 percent width of y1 - y2 at sdx = 2 is 3.18 with one latent and 15.08 with two. newdata must hold sdx. On the data the model was fitted on the latent values are the fit's own.

A level the fit did NOT see is different: there is no conditional law to draw from, so allow_new_levels = TRUE draws its effect from the block's own estimated covariance, once per replicate, at that replicate's parameters. This is brms's sample_new_levels = "gaussian" and it is the whole between-group variance, so the interval at an unseen level is wider than at a known one.

re_formula and propagate_error are different questions

The two are easy to confuse and neither can express the other.

re_formula chooses WHICH terms enter the prediction. Only re_formula = NA can say "predict for an average group", by leaving the group effects out of the linear predictor.
propagate_error chooses whether the error in the ESTIMATES is propagated into the interval. Only propagate_error = FALSE can say "include this group's own effect but treat it as known".

So a prediction for the group in front of you, with its effect held at the fitted value, is propagate_error = FALSE, and a prediction for a group you have not met is re_formula = NA. They combine: re_formula = NA with propagate_error = FALSE is the population curve with no estimation error at all.

brms has no such argument, because its draws always carry both: a posterior draw has its own parameters and its own group effects, so there is nothing to switch off. A brms user looking for the missing argument is looking for this one.

Only the terms re_formula keeps are drawn. A population smooth and a gp() curve stay at their modes under every setting.

A fit whose covariance could not be recovered from the Hessian has no second term to draw: the simulation falls back to the estimates and warns, rather than drawing from a matrix of NaN. A fit made with quadrature = TRUE has no group effects in its joint precision: its group effects stay at their modes, and a warning says so.

Under REML = TRUE (or frmtmb_control(profile = TRUE)) the fixed effects are integrated out of the outer problem, so they are not in vcov(object, full = TRUE). They are drawn from their joint covariance with the outer parameters, which the joint precision gives, and the group effects are then drawn given both.

Ordinal and categorical responses

There is no mean to summarize, so predict() returns the simulated proportion of each category, one column per category named P(Y = k), which is brms's shape. fitted() returns the MODELLED category probabilities with their standard errors instead, which is the smoother quantity.

Examples

set.seed(1)
dd <- data.frame(x = rnorm(100))
dd$y <- rpois(100, exp(0.3 + 0.4 * dd$x))
fit <- frm(bf(y ~ x) + poisson(), data = dd)

# the predictive summary: observation noise is in it
head(predict(fit, ndraws = 200))
#>      Estimate Est.Error Q2.5 Q97.5
#> [1,]    0.765 0.8738381    0 3.000
#> [2,]    1.270 1.0876650    0 4.000
#> [3,]    0.765 0.8445956    0 3.000
#> [4,]    2.430 1.6087387    0 6.000
#> [5,]    1.235 1.1072057    0 4.000
#> [6,]    0.665 0.8099724    0 2.025
# against the expected response, which is much tighter
head(fitted(fit))
#>       Estimate Est.Error      Q2.5     Q97.5
#> [1,] 0.7751024 0.1130250 0.5824215 1.0315275
#> [2,] 1.1849807 0.1126043 0.9836142 1.4275712
#> [3,] 0.6946365 0.1132732 0.5046074 0.9562283
#> [4,] 2.4828747 0.3477984 1.8867719 3.2673090
#> [5,] 1.2791025 0.1161520 1.0705573 1.5282725
#> [6,] 0.7001766 0.1132798 0.5099103 0.9614382

# the draws themselves
dim(predict(fit, summary = FALSE, ndraws = 50))
#> [1]  50 100

# the linear predictor is frm_linpred() now
head(frm_linpred(fit, type = "link"))
#>          1          2          3          4          5          6 
#> -0.2547601  0.1697265 -0.3643666  0.9094170  0.2461587 -0.3564227