Subject-wise maximum likelihood estimates of a bmm model

Fits a model with no pooling — one independent set of parameters per subject — and returns the estimates in the same tibble shape a hierarchical fit gives, so that the two can be scored against one truth and compared. It exists to measure what hierarchical estimation buys, not to recommend maximum likelihood for inference.

Usage

fit_ml(
  model,
  data,
  formula = NULL,
  by = "id",
  method = c("stan", "optim"),
  prior = c("flat", "default"),
  ci_level = 0.95,
  draws = 1000,
  max_abs_link = 20,
  start = NULL,
  nll = NULL,
  file = NULL,
  refit = c("on_change", "never", "always"),
  seed = NULL,
  ...,
  .fitter = NULL
)

Arguments

model: A bmmodel, as bmm::mixture2p() returns.
data: The data to fit.
formula: An already-reduced bmmformula with no group-level terms. NULL, the default, writes p ~ 0 + <by> for every free parameter. Supply one for a design fit_ml() does not write itself. method = "optim" builds no formula and does not accept one. For a model with alternative parameterizations (criterion or bias of sdt_yn), the member the formula names is estimated; the default formula writes bmm's default, criterion.
by: The column defining one independent set of parameters. A single column in this version.
method: "stan" optimises bmm's own generated likelihood through algorithm = "laplace", re-implementing nothing, and needs CmdStan. "optim" maximises bmm's R density for the model directly with stats::optim() and needs no compiler. See the two-routes section.
prior: "flat", the default, puts an empty prior on every free parameter so the mode is the maximum-likelihood estimate. "default" keeps bmm's priors, which makes the mode a penalised estimate rather than an ML one. method = "optim" has no priors at all and accepts "flat" only.
ci_level: The interval level. The interval is the Laplace approximation's on the Stan route (ci_method = "laplace") and a Wald interval from the optimiser's Hessian on the other (ci_method = "wald"). Both routes also return the central 50 % interval of the same kind in ci_low_50 and ci_high_50, whatever ci_level is, and in post_mean_link and post_sd_link the point estimate and its standard error (on the Stan route, the mean and SD of the Laplace draws), so that a posterior z-score can be computed for ML rows too.
draws: How many draws to take from the Laplace approximation. method = "stan" only.
max_abs_link: The largest abs(estimate) on the link scale that still counts as converged.
start: Named link-scale start values for method = "optim", one per free parameter. NULL, the default, starts at zero. On 40 subjects of mixture2p at 100 trials, a start at zero reached the same optimum as a start at the mean of the Laplace draws of the "stan" route, to 4.3e-07 on the link scale. Under bmm::ezdm(version = "4par") it does not: see the section "Four-parameter ezdm, one subject at a time". Its names also choose the member of an alternative parameterization that is estimated.
nll: A function of (pars, data, model) returning the negative log-likelihood of one subject's data as a single number, where pars is a named list on the natural scale with the model's fixed parameters included. Supply it for a model bmmtools carries no density for. method = "optim" only.
file, refit: Passed to fit_cached(). NULL uses a temporary file, so nothing is cached between sessions. method = "optim" does not cache and ignores both: it produces no fit object, and at well under a second per call a cache would add an invalidation surface without saving anything.
seed: Passed to the fitter. method = "optim" is deterministic and ignores it.
...: Passed to fit_cached() and on to the fitter.
.fitter: As in fit_cached(). method = "stan" only.

Value

A bmmtools_ml: a tibble satisfying the estimates contract, one row per subject and term, with level = "subject" and estimator = "ml". attr(x, "ml_cells") holds one row per subject with its term and convergence counts.

Two routes

method = "stan" optimises bmm's own generated likelihood with algorithm = "laplace". It re-implements nothing, so it works for any bmm model, and it needs a working CmdStan.

method = "optim" maximises the likelihood directly with stats::optim() on bmm's R density for the model. It needs no compiler and is roughly twenty times faster — 0.35 s against 8 s for 40 subjects of mixture2p — but it is capped at the models bmmtools carries a density for, listed in the error it raises otherwise, and nll is the escape hatch for anything else.

The two routes return the same columns, the same class and the same estimator, and differ only in ci_method. Measured on 40 subjects of mixture2p at 100 trials, they agree to Monte Carlo error on the estimate and their Wald and Laplace standard errors agree to a mean ratio of 1.000 (range 0.95 to 1.05), giving the same coverage and calibration_slope to three decimals. The route is not a moderator of a scored number.

One fit, not one per subject (`method = "stan"`)

Under p ~ 0 + <by> every free parameter belongs to one subject and bmm's remaining parameters are constants, so the log posterior is a sum of per-subject terms and the joint mode is the vector of per-subject modes. The Stan route therefore runs a single fit for the whole data set. This was checked rather than assumed: over 40 subjects, the joint Stan mode and independent optimisation of bmm's own density agreed to 3.0e-06 on the link scale.

The cost is that convergence is no longer separable by subject — one subject at a boundary can stall the whole optimisation — so the per-cell verdict in ml_cells is the range check described below. The optim route does not have this property: it is one independent optimisation per subject.

What counts as converged

converged is TRUE where the estimate is finite and abs(estimate) <= max_abs_link on the link scale, and, on the optim route, where optim() also returned its success code. brms keeps no optimiser return code on a Laplace fit, so on the Stan route the range check is the only verdict available — and not a weak one: fitting a single subject, Stan reported convergence at a point 0.064 log-likelihood short of the mode.

A row that fails is kept with estimate = NA rather than dropped, because metric_bias() and metric_rmse() drop incomplete pairs silently: dropping would compare the Bayesian estimator on every subject against maximum likelihood on the easiest ones and report the difference as the estimator's. Report n alongside any comparison.

On the optim route a Hessian that cannot be inverted gives a row with no interval: ci_low and ci_high are NA while the estimate and converged stand, so the subject counts in bias, rmse and r and drops out of coverage only.

rhat, ess_bulk and ess_tail are always NA. Neither route produces MCMC draws, so the values posterior would return describe the approximation and not the estimate.

Four-parameter `ezdm`, one subject at a time

Under bmm::ezdm(version = "4par") one subject's row holds five statistics — the mean and variance of the response times at each boundary, and n_upper — for four parameters. A boundary reached fewer than twice leaves its two RT summaries NA, three statistics remain for four parameters, and no-pooling cannot identify that subject. On the optim route such a subject nearly always fails to converge: of 30 subjects simulated at drift 1.5, bound 1.5, ndt 0.3 and zr 0.7 with no spread between them (seed = 1), 23 failed at 25 trials and 6 at 100, exactly the subjects with fewer than two lower-boundary responses. Two or more responses at each boundary are necessary for identification but not sufficient for a good fit. At 400 trials every subject reached both boundaries and one still failed: from the default start at zero, optim() used up its 100 BFGS iterations without a success code, although it stood within 5e-05 on the link scale of the optimum that a start at the true values reaches. Nor does converged vouch for every subject that passes it. In the same run, subject 17 (15 lower-boundary responses) passed with drift -0.002 against a true 1.5 and log bound 6.73 against a true 0.41, at a negative log-likelihood of 117.26 against -6.28 at the true values. That point is not an optimum. The gradient there is -50 on drift and 15 on log bound, and lowering log bound by 1 lowers the negative log-likelihood to 107.56: BFGS stalled short of the optimum in an ill-conditioned valley and returned its success code all the same. Nelder-Mead started from that point with maxit = 1000 reaches -8.10, the optimum that a start at the true values reaches; with the default maxit of 500 it stops at -6.24 without a success code. At 25 trials with SDs between subjects of 0.3 on drift, 0.2 on bound, 0.2 on ndt and 0.3 on zr (seed = 1), one subject with no lower-boundary response passed at a drift of -1.36 against a true 1.31, with Wald intervals 3308 wide on bound and on zr. SDs, estimates and widths are on the link scale. Give each subject enough trials to reach both boundaries and pass start near plausible values: started at the true values, all 30 subjects of the 400-trial run converged and subject 17 returned a drift of 1.35. Compare the estimates with a hierarchical fit by bmm::bmm(), which pools across subjects. Both routes of fit_ml() are no-pooling fits.

Examples

if (FALSE) { # \dontrun{
model <- bmm::mixture2p(resp_error = "y")
sim <- simulate_recovery(
  model,
  pars = c(kappa = 2, thetat = 1), sds = c(kappa = 0.3, thetat = 0.5),
  n_subjects = 40, n_trials = 100, seed = 1
)
ml <- fit_ml(model, sim$data)
summary(recover_subjects(ml, sim$truth$subjects))

# the same thing without a compiler, and without Stan
ml2 <- fit_ml(model, sim$data, method = "optim")
} # }