Pointwise log-likelihood of posterior draws

The ndraws x nobs matrix of per-observation log-densities, each row evaluated at one draw's own parameter vector. Every leave-one-out and information-criterion quantity in this package is a function of this matrix, and frmtmb::loo(), frmtmb::waic() and loo::psis() take it directly.

Usage

# S3 method for class 'frmtmb_draws'
log_lik(
  object,
  newdata = NULL,
  re_formula = arg_unset(),
  resp = NULL,
  ndraws = NULL,
  draw_ids = NULL,
  pointwise = FALSE,
  combine = TRUE,
  add_point_estimate = FALSE,
  cores = NULL,
  point_estimate = NULL,
  ndraws_point_estimate = 1,
  ...
)

# S3 method for class 'frmtmb_draws'
logLik(
  object,
  newdata = NULL,
  re_formula = NULL,
  resp = NULL,
  ndraws = NULL,
  draw_ids = NULL,
  pointwise = FALSE,
  combine = TRUE,
  cores = NULL,
  ...
)

Arguments

object: A frmtmb_draws from frm_sample().
newdata, re_formula: Accepted in brms's own second and third positions and refused: see What is conditioned on. The refusal names the reason and the replacement.
resp: For a multivariate model, the responses whose density to report: the default is every response, summed without rescor and joint with it (see Multivariate models).
ndraws: Number of draws to use, evenly spaced through the matrix (default: all of them).
draw_ids: The draws to use, by row index, instead of the evenly spaced subsample ndraws takes. Give one or the other.
pointwise, combine, add_point_estimate, cores: brms's slots, in brms's positions. pointwise = TRUE (brms's function of one observation) and add_point_estimate = TRUE are refused by name; combine = FALSE is accepted on a univariate model, on one response, and on a set_rescor(TRUE) model, where there is nothing to combine; cores is accepted and unused, because the matrix is built in this process.
point_estimate, ndraws_point_estimate: brms's arguments: collapse the draws to their "mean" or "median" first and use that one parameter vector, repeated ndraws_point_estimate times.
...: Refused: an argument the method does not have is an error naming it, rather than silently changing nothing.

Value

A numeric matrix with one row per draw and one column per observation (the rows the model was fitted on).

Details

logLik() of draws is log_lik(), as brms's logLik() of a fit is: the draws by observations matrix, not stats' logLik object, which describes one parameter vector.

What is conditioned on

The density is CONDITIONAL on the draw's own group-level values: frm_sample() samples b alongside everything else, so each row of the draws matrix is a complete parameter vector and no integration is left to do. This is exactly brms's convention, where the Stan model also samples the group-level parameters. The consequence is worth stating plainly: loo() on such a matrix is leave-one-OBSERVATION-out with the groups held fixed, not leave-one-group-out, and for a model with few observations per group the two differ.

Addition terms enter exactly as they enter the fitted objective, because they enter through the same code: cens() replaces a row's density by the matching CDF difference, trunc() divides by the window mass, weights() multiplies the row's contribution, and trials() is the family's own argument. brms composes them in the same order (log_lik_censor(), log_lik_truncate(), log_lik_weight()).

Autocorrelation, and what a row conditions on

brms's default cov = FALSE form of ar(), ma() and arma() (frmtmb::frmtmb-autocor) has a column per observation, and it is the column brms's own log_lik() returns. A row's mu is the one-step conditional mean: the linear predictor plus a regression on the OBSERVED earlier residuals of that row's group. A group starts with no residuals behind it, so lag i first reaches the row at within-group position i + 1: the FIRST row of each group takes the unshifted mean, a row at position k carries only the lags up to k - 1, and from position max(p, q) + 1 on a row carries all of them. Only the first row is unshifted, at every order. So the density of row t is conditional on y at the rows before it, and the columns multiply to the whole likelihood exactly as they do for an independent model.

Say what that makes the leave-one-out quantity. PSIS-LOO over these columns leaves out one CONDITIONAL density, the density of y_t given the observed past, and not the series: the retained columns still read y_t through their own lagged residuals, so it is not a forecast of an unseen point and it is not leave-one-group-out either. It is brms's number, computed the way brms computes it, and it answers "how well does the one-step-ahead density fit this row" rather than "how well would the model predict a held-out series". For the latter, refit without the held-out group.

Residual covariance terms

A residual covariance term (cov = TRUE, cosy(), unstr(), fcor(), and frmtmb.spatial's sar()) makes the residuals of a group one joint draw, so a row has no density of its own. A column is then brms's leave-one-out conditional density: the density of the row given the other rows of its group at the draw's parameters, normal with precision \(Q_{ii}\) around \(y_i - (Qe)_i / Q_{ii}\), or for student() the conditional t, after the loo package's vignette on non-factorized models (frmtmb::frmtmb-autocor, "The pointwise log-density"). The values are brms's log_lik(), and they do NOT sum to the log-likelihood: each conditions on the rest of its group. The cost per draw is one small inverse per pattern of time points for the ar() family, and the nonzeros of the weights matrix for sar(); brms inverts a dense covariance over each group per draw.

A group that skips a time point inside its range differs from brms under ar(), ma() and arma(): frmtmb reads the lag from the time levels and brms from the position within the group, so the two fit different models there and their columns differ. In a measured case with one such group (22 rows, 3 draws), the largest difference over the matrix was 17 percent of brms's largest value, and on the rows of the groups without a gap the two agreed to 8e-16. cosy() has one correlation for every pair and brms's unstr() reads the time levels too, so those agree with gaps.

Multivariate models

With set_rescor(TRUE) a column is the joint density of the row's response VECTOR, so the matrix keeps one column per observation. resp naming several responses gives the joint density of those responses, and naming one gives that response's own marginal density (the normal, or for student() the t with the shared nu), as brms's log_lik(resp =) computes them. combine = FALSE has nothing to combine there and returns the same matrix, as in brms. Without rescor, the responses are independent given the predictors and the default sums their log-densities per row; pass resp to get one response's contribution alone.

Likelihoods with no per-observation column

A model whose smallest independent unit is a group has no per-observation column to leave out, and this refuses rather than inventing one: a frmtmb.latent::hmm() sequence, and a group-level mixture (mixture(groups = )). An frmtmb.latent::lca() subject is one row, so its column is well defined and is not refused. In-model imputation (mi(), me()) is refused for the same kind of reason: a latent value is a parameter, not an observation. Use AIC() on the maximum-likelihood fits or frmtmb::frm_bootstrap() for those.

A family that declares how its likelihood factorizes (frmtmb::frmtmb_structure(loglik_group = ) or (loglik_row = )) is not in that position and is not refused. Its columns are the pieces it declares, at the COARSEST granularity it gives: a family that groups its likelihood is saying that its rows are not independently droppable, which is exactly the question a leave-one-out column asks. The matrix then carries attr(x, "unit") naming what a column is. loo::loo.matrix() never sees that attribute and its printout says only "Computed from N by K log-likelihood matrix", which reads exactly like a per-observation one, so loo(), waic() and psis() emit a message naming the unit and the column count when they are handed such a matrix. The elpd is then leave-one-UNIT-out, and the number itself carries no mark of that.

Examples

# \donttest{
if (requireNamespace("tmbstan", quietly = TRUE) &&
    requireNamespace("rstan", quietly = TRUE) &&
    !frmtmb.sample:::tmbstan_build_broken()) {
  set.seed(9)
  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)), family = gaussian(), data = dd)
  ds <- frm_sample(fit, chains = 1, iter = 500, refresh = 0)

  ll <- log_lik(ds)
  dim(ll)
  # the column means are the per-observation expected log-densities
  head(colMeans(ll))
}
#> frm_sample(): default priors (brms 2.23 defaults; prior = "flat" puts brms's (flat) in every slot)
#>   Intercept          student_t(3, 0.8, 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: There were 4 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
#> [1] -0.8735093 -2.1721647 -1.2481160 -2.5636748 -0.8763703 -3.1446238
# }