Skip to content

Add multinomial processing tree models with mpt() (MPT 1/6, #255) - #484

Merged
GidonFrischkorn merged 62 commits into
developfrom
255-mpt-core
Oct 14, 2026
Merged

GidonFrischkorn merged 62 commits into
developfrom
255-mpt-core

Conversation

@GidonFrischkorn

@GidonFrischkorn GidonFrischkorn commented Oct 3, 2026 •

Copy link
Copy Markdown
Collaborator

Part of #255.

Base branch: develop (rebased on 6f04c69, 2026-10-04; synced with develop 186250e by merge 036de49, 2026-10-08; develop 222a7ef merged in 2026-10-08; head e7c8504, 2026-10-09). This is PR 1 of 6 in the MPT stack. Each PR targets the branch below it, so each diff shows one feature group only:

  1. 255-mpt-core → develop (this PR): the model itself
  2. 255-mpt-constraints → 255-mpt-core: restrictions = and simplex =
  3. 255-mpt-design → 255-mpt-constraints: covariates = and mpt_tree(impossible =)
  4. 255-mpt-import → 255-mpt-design: mpt_from_string() and mpt_from_eqn()
  5. 255-mpt-plot → 255-mpt-import: plot() for trees
  6. 255-mpt-article → 255-mpt-plot: the bmm_mpt article

This PR: 22 files, +3063/-8.

What this adds

mpt() fits multinomial processing tree models. You write each tree's branch probabilities with mpt_tree(), give a parameter the same name in every tree that shares it, and fit the model to response counts. tree_id names the data column that says which tree a row belongs to.

tree_old <- mpt_tree("old", list(
  old = "D + (1 - D) * g",
  new = "(1 - D) * (1 - g)"
))
tree_new <- mpt_tree("new", list(
  old = "(1 - D) * g",
  new = "D + (1 - D) * (1 - g)"
))
model <- mpt(list(tree_old, tree_new), tree_id = "item_type")

fit <- bmm(bmf(D ~ 1 + (1 | id), g ~ 1 + (1 | id)), data, model)
  • The latent parameters are read from the branch expressions. Shared names are equality constraints.
  • links = "logit" (default) or "probit". Parameters are estimated on the latent scale and predicted with ordinary bmf() formulas. Probit with (1 |p| id) is the TreeBUGS latent-trait structure.
  • A parameter can be fixed in the formula: bmf(D ~ 1, g = 0.5).
  • Non-linear parameter formulas, such as D ~ inv_logit(Dmax) * (1 - exp(-exp(rate) * ptime)), bypass the link. Their sub-parameters are collected as in m3 custom models.
  • print(model) lists the trees, the count of free parameters against response categories minus 1 (summed over trees), and the rank of the Jacobian of all category probabilities, which names the parameters of a non-identified model (see "Local identifiability").
  • dmpt() and rmpt() mirror dm3() and rm3(), including the log = TRUE default and unpack =. They accept parameter values strictly between 0 and 1, as bmm() does for fixed values.

Framework touches

  • print.bmmodel() now calls a new internal generic print_model_details(), with a no-op default. Classes are ordered general to specific, so a print.mpt() method would never dispatch.
  • response_annotations(): mpt shares the m3 annotation for count columns.
  • uses_aggregate_data(), model_groups, data_column_roles and revert_check_data() gain mpt entries (next section).
  • update.bmmfit() re-checks the stored model when formula. is given, for every model (4224afc). It never called check_model(), so a model checked with the old formula kept what that formula derived. For mpt that broke update() outright (see "Decided", D3). Other models: check_model() on an already checked model is a no-op for the same formula (regex replacement is guarded by regex_applied; check_model.m3_custom() only adds parameters the new formula introduces). The full suite passes at this level (FAIL 0, SKIP 1, PASS 6092).

Rebased onto develop 6f04c69 (2026-10-04)

The stack was cut on develop 2b25ab5. Develop's 47 commits since then renamed functions this PR calls and added contracts; each is ported in its own commit:

Synced with develop 186250e (2026-10-08)

Develop merged in once (036de49), including #508, #509 and #512. One conflict, in .ce_compute_and_transform(): this PR's inherits(x$bmm$model, c("m3", "mpt")) bypass gives way to develop's route for multinomial families (#512), which mpt reaches through its brms::multinomial() family. The bypass function no longer exists on develop. On a real 2HTM fit the route gives the same estimates and intervals as the bypass at robust = TRUE; the default estimate is now the posterior median, as in brms. conditions, re_formula and spaghetti now apply, and categorical and method = "posterior_predict" stop with a message naming mpt. Calls without effects print In log(old) : NaNs produced: brms builds a grid for the interaction of the tree indicator columns, which develop's route computes and then discards. The estimates are unaffected; the warning comes from develop's route, not from this PR.

Two of develop's new rules are guarded by tests that list their models by hand, so mpt gets its own (4dfa85d; both fail on 8dc5e4c): the fit keeps the name of the user's data (#508), and (1 |p| id), which ?mpt recommends, is not read as a data column (#509).

How mpt() meets the contracts develop added since the original branch

The original branch (new_model/mpt, e2afd98) predates several framework contracts on develop. I merged develop (2b25ab5) and adapted mpt() to each contract in its own commit.

  • Links (links argument: an unrecognized parameter name is appended instead of rejected, in every model constructor and in the model template #420). .model_mpt() assigns links through set_links(). settable_links.mpt() returns the latent probability parameters, and settable_link_functions.mpt() allows "logit" and "probit". A parameter switched to the other link after construction gets that link's default prior, unless the user replaced its prior (.mpt_match_link_priors()). check_model.mpt() sets identity links for parameters with a non-linear formula. It records these in the links_checked attribute, so a later check_links() does not read them as user changes. Mutation check on the full tree: 5 of 5 mutants killed.
  • other_vars and field order. All model state (trees, tree_id, link, generated indicator names) lives in other_vars. The top-level fields follow the canonical order. MPT is unreleased, so no saved fit depends on the old layout.
  • update() via revert_check_data(). revert_check_data.mpt() unpacks Y into the count columns and drops nTrials and the Idx_ tree indicators. brms keeps the tree column only when a formula names it, so the method rebuilds that column from the indicators and appends it to the "rebuilt" attribute. It then chains NextMethod(). Three stored-frame cases cover two trees, the tree column as a predictor, and a single tree. Mutation check on the full tree: 5 of 6 mutants killed. The survivor (assign instead of append "rebuilt") is equivalent for mpt.
  • sd default priors. exponential(1) with the logit link, the rate of the mixture weights. exponential(2) with the probit link (decided, D1 below). Mutation check: 2 of 2 killed.
  • effects default priors. See the next section.
  • bmm_data_check(). MPT rows are aggregated counts. uses_aggregate_data() now includes mpt, so the "fewer than 10 observations per cell" finding no longer fires on them.
  • Overview table (Data columns, key parameters and minimum trials in the task-to-model overview #467). mpt joins model_groups as "Processing-tree models". data_column_roles$mpt labels trees (the count columns) and tree_id. The README and overview-article sections for the new group come with PR 6.
  • No change needed: lkj(2) on correlations, the default create_initfun(), pp_check() (multinomial path, as for m3) and native_transform(), which reads model$links.

The effects prior

The original branch gave effects half the intercept's scale: logistic(0, 0.5) with logit, normal(0, 0.5) with probit. Under 1 + cond coding that pulls condition contrasts toward zero. A posterior-mode sweep on expected 2HTM counts measured how much. For realistic designs (contrast up to 3 logits, at least 50 old and 50 new items per condition), the largest shrinkage of the condition-B detection probability was as below. z is the shrinkage divided by the posterior SD under a flat prior.

link effects prior max shrinkage max z
logit logistic(0, 0.5) (old) 0.043 1.91
logit logistic(0, 1) (new) 0.018 0.83
probit normal(0, 0.5) (old) 0.048 1.11
probit normal(0, 1) (new) 0.015 0.34

On oberauer_lewandowsky_2019_e1, pooled over participants, every coding gives modes within 0.0016 of the flat-prior mode. Per participant, the median "old same" minus "new distractors" contrast in Pb is 0.1245 under 1 + cond with the old default, 0.1416 under 1 + cond with logistic(0, 1), and 0.1454 under 0 + cond with logistic(0, 1).

So effects now get the intercept's prior: logistic(0, 1) with logit, normal(0, 1) with probit. With that, 1 + cond and 0 + cond put the same prior on every coefficient. Sub-parameters of non-linear formulas keep normal(0, 1) for the intercept and normal(0, 0.5) for effects. The test "condition effects under reference coding get the intercept's prior scale" failed before the change. The change enters no validation fit: all 8 cached validation fits have identical Stan code and prior tables before and after.

Design decisions

  • Native link machinery, no reparameterization layer. The exploration prototype emitted nlf(D ~ inv_logit(lD)) + lf(lD ~ ...) per parameter. Here check_formula.mpt() prepends generated category formulas, like construct_m3_act_funs(), and apply_links() substitutes inv_logit() or Phi(). Parameters keep their names and live on the latent scale, like the log-link parameters of m3. A unit test shows the category probabilities equal the production m3 likelihood to < 1e-10 through the analytic simple-span bijection.
  • One multinomial path. Two-category models use the same multinomial(refcat = NA) family as models with more categories. brms supports two-category multicolumn families. This keeps one code path and gives the multinomial pp_check() for free.
  • Safety guards. Every variable-free arithmetic constant in a branch expression is folded to its value, because Stan compiles 1/4 as integer division (0). Parameter and category names are checked against the brms nlpar rules at construction. Branch sums are checked at five numeric test points at construction, where every parameter takes its own value; a tree whose branches do not sum to 1, a branch outside (0, 1] at a test point (2 * a, or 0 * D), or a branch that is the constant 0, is an error. User formulas for response categories are rejected, and so is a response category on the right-hand side of any formula (Dn ~ old), which brms otherwise refuses as a "circular dependency". Numeric constants that brms would write in scientific notation are rejected: brms's operator spacing turns 6.7e-05 into the subtraction 6.7e - 05, and the Stan error is very hard to trace. Every spelling of such a value is affected (0.00010 is 1e-04), so the message quotes the constant as written, with the form brms would emit, and suggests a rescaling: (1 + 9999 * p) / 10000 for a small constant, p / 1000 / 100 for a large one.
  • Fixed values stay probabilities. bmf(g = 0.5) stays on the probability scale in the model object, which is what print() and parameter_info() show. configure_prior.mpt() maps the value through the link when it builds the constant() prior, so check_model() is idempotent.
  • tree_id, not condition. A separate tree is needed only when the branch expressions differ, for example old versus new probes. Effects of experimental factors belong in the parameter formulas, as in every other bmm model. The two roles are independent, and one column can serve both. Naming the argument condition invited the reading that each factor level needs its own tree. That costs real work: the category probabilities sum over all trees, so every duplicate tree is evaluated for every row. tree_id also matches the EQN vocabulary (multiTree, MPTinR and TreeBUGS call column 1 the tree).
  • Order constraints by documentation. Do > Dn is written with the existing non-linear formula mechanism: Do ~ Dn + (1 - Dn) * inv_logit(phi). A latent-scale offset (TreeBUGS style) works too but needs adapt_delta = 0.99. ?mpt points to the recipe.
  • No init_ranges. mpt() declares none and uses develop's default init function. Every validation fit below ran with it.

Decided (GF, 2026-10-03)

  • D1. Probit sd rate stays exponential(2). Logit uses exponential(1) (median 0.69); a probit SD is about the logit SD / 1.7; the round rate matches the SDT rate scheme. No code change.
  • D2. Sub-parameters of non-linear formulas keep brms's default sd prior, documented (4cc4087). Measured with default_prior() (brms 2.23.0): student_t(3, 0, 2.5) under both links. The message check_model.mpt() prints for them and ?mpt now name it; a test ties the message to the prior default_prior() reports (red before the change).
  • D3. Review finding 7 was real, fixed here (4224afc). A checked model re-used with a formula in which a formerly non-linear parameter is linear did not just keep the identity link without a prior: it failed in brms (The parameter 'a' is not a valid distributional or non-linear parameter), on a real update() and on default_prior()/stancode() with fit$bmm$model. check_model.mpt() now starts from the constructor's parameters, the links the parameters had before the bypass (kept as an attribute, so a link switched by the user survives) and the model's latent prior, which it then matches to a switched link; update.bmmfit() re-checks the stored model when formula. is given (see "Framework touches"). Tests: "a parameter made linear again in a re-check gets back its link and prior", "update() to a linear formula restores the link and prior of a non-linear parameter" (red before the change). Mutation check: 4 of 4 killed; the two that remove the reset or the update() call kill the update() test by the brms error itself. Real fits after the fix: the update samples (max Rhat 1.0033, 0 divergences), D is back on logit with logistic(0, 1), and its Stan code is identical to a fresh model's. The 2HTM and PCM Stan code did not change.

Changes since the first review (2026-10-04 to 2026-10-06)

An internal review (2026-10-04): merge after fixes, nothing blocking. Each fix is its own commit with a test that failed before it.

  • Branch expressions. Every variable-free constant is folded, so 1/(2*2), -1/4 and 1/2^2 no longer reach Stan as integer division (0d49ed9). A tree whose branches do not sum to 1 and a branch that is the constant 0 are errors at mpt(), not warnings or Stan failures (6dbe773, 16bbdc5, 45f72d8). The test values follow a golden-ratio sequence, because the four-value cycle let a typo between parameters four places apart pass (1a7a242).
  • Formulas and data. A symbol in a non-linear formula that is neither a data column nor a formula entry is an error (50e82e7). check_data() warns on missing counts and refuses user columns Y and nTrials (4ed2d83).
  • dmpt() and rmpt(). A named count vector is matched to the categories by name (7bbb6d0). Parameters the selected tree does not use are ignored (08111f3).
  • Count line. It counts only latent parameters that are not fixed, and says it holds for intercept-only formulas (d9be1c6).
  • ?mpt. The default priors, the data layout, native_parameters() for probabilities, correlated person effects with |p|, posterior_predict()/posterior_epred() for goodness of fit, and why summary() shows a fixed g = 0.5 as 0 on the latent scale (993bca7, d29e4a5, 7e97f0d, 75f693e, 3306830, df0967e).

Local identifiability (0de1416)

The count of free parameters against degrees of freedom misses models that pass it but are not identified. The unrestricted unitization and hybrid EQN files of Al Hadhrami et al. (2025) pass the count, yet their Jacobians have rank 19 of 21 and 21 of 27; with the design constants fixed the ranks are 12 of 12 and 18 of 18. Fits of such models sample the prior along the redundant directions without any message.

print() now reports the rank of the Jacobian of all category probabilities with respect to the free parameters. The derivatives are exact (stats::D() on the stored branch calls), at five deterministic interior test points; parameters fixed in the formula enter at their values. A deficit names the parameters involved and, separately, parameters whose derivatives vanish (they cancel from the likelihood). check_model() warns about a deficit unless a formula could resolve it. When any formula has a data predictor, the predictor can identify the model across design cells, also from a parameter outside the deficit, so the deficit is only a message. A non-linear formula that defines or reads a parameter of the deficit (Dn ~ Do, Do ~ Dn + (1 - Dn) * inv_logit(phi)) also makes it a message: the rank does not read the formula, and print() of the checked model says the formula was not analysed. A non-linear formula elsewhere (h ~ g next to a deficit in D and r) leaves the warning. Sub-parameters (a ~ inv_logit(u + v)) are not checked either; the message about them says so. A branch that D() cannot differentiate leaves the rank uncomputed, and the printout says so. The count line stays.

Verified in four rounds by an independent agent: the rank agrees with exact derivatives on 488 random models and a 19-model zoo, chains are exact to depth 32, and a print takes about 9 ms. PR 2 extends the check to simplex groups and PR 3 to covariates.

Decided 2026-10-06, implemented

  • conditional_effects() works for mpt fits. It stopped with "Please set 'categorical' to TRUE" on every fit. Since the develop sync, mpt fits go through develop's route for multinomial families (conditional_effects() fails on SDT rating, CDP and ranking fits and ignores conditions on m3 fits #512); this PR's class bypass is gone (see "Synced with develop 186250e"). With unequal trial counts per row the prediction grid held the mean trial count, which brms refused, so check_data.mpt() stores nTrials as an integer. A real cmdstanr fit with unequal counts matches native_parameters() (item 7, be49641).
  • The unused internal .mpt_constructor_args() is deleted (item 8, 342ba51).

Changes since the merge-readiness read (2026-10-07)

Each fix has a test that failed before it.

  • Formula ties are not reported as a rank deficit (246ef97). bmf(Do ~ 1, Dn ~ Do, g ~ 1) on the 2HTM warned "The model is not identified", and a ~ inv_logit(u + v) printed "locally identified". A deficit under a non-linear formula is now a message, like a data predictor. This also turns the warning for Do ~ inv_logit(phi) into a message, although that deficit is real. The branch-sum error no longer says restrictions do not belong in the formula. Narrowed in the next section.
  • Branches outside (0, 1] are refused at mpt() (bdaa62f). 2 * a and 1 - 2 * a sum to 1 but reached Stan as log(2 * a). A branch that is 0 for every value, such as 0 * D or (1 - a) * 0, is refused the same way; decision 1 had left it uncaught.
  • The constant message says what is wrong and what works (4b47a54); see "Safety guards".
  • dmpt()/rmpt() refuse 0 and 1 (ed4c138), matching check_model() for fixed values.
  • Refactors, output unchanged (5c21a73, 05013d0, 9c3d029): one helper for the tree parameters, the count line takes its numbers from the rank check, four single-use variables inlined. Stancode of 2HTM and PCM and print() of seven models are byte-identical before and after.

Changes since the re-review of the identifiability policy (2026-10-07)

Fixes were implemented by one agent and read by a second, independent one. Each fix has a test that failed before it.

  • A non-linear formula no longer hides an unrelated deficit (0f9f9e9, 00fb817). After 246ef97, h ~ inv_logit(k) or h ~ g turned the warning for a deficit in D and r (D * r) into a message, although nothing in the formula can separate them. The message now needs a data predictor anywhere, or a non-linear formula that defines or reads a parameter of the deficit. The first version also warned when a predictor sat outside the deficit. The independent read showed that such a predictor can still identify the model (in h * r + (1 - h) * D, h ~ cond separates D and r), so any data predictor keeps the message (00fb817).
  • The message names design cells only when predictors make them (73166a4). Under a non-linear formula alone it says "not identified by the branch expressions alone", as print() does.
  • A response category on the right-hand side of a formula is an error at check_model() (8a1b54b). Dn ~ old previously printed the identifiability message and then died in brms with "Cannot handle circular dependency structures".
  • The message about non-linear sub-parameters says that their identification is not checked (8339e91). With a ~ inv_logit(u + v), bmm() printed nothing about u and v.
  • The scientific-constant message fits the constant (1f2e416). It quotes the constant as written and shows the scientific form only when it differs; large constants get their own rewrite example.
  • A comment on refused zero branches now covers underflow at a test value (3a8e445; (1 - D)^n is refused from n = 392), and one lint fix (8dc5e4c).

Validation

Validation ran on a git archive snapshot of the full tree at an earlier head of this stack, with fresh fits. It was not rerun per level. What ties it to this level is the stancode gate: the 2HTM and the single-tree pair-clustering model (PCM) compile to Stan code identical to the final tree's. After that head, the R code changed only in a parameter description, comments, overview labels and inlined variables in revert_check_data.mpt(), and in the checks and messages of the 2026-10-07 round, none of which changes the Stan code of a valid model.

  • Parameter recovery. Hierarchical 2HTM, simulated by hand, independent of the package's distribution functions: mu_lD 0.874 [0.746, 0.997] (true 0.837), mu_lg 0.109 [0.000, 0.217] (true 0.026), sd_lD 0.420 [0.325, 0.534] (true 0.461), sd_lg 0.252 [0.092, 0.392] (true 0.274); max Rhat 1.0059. PCM, simulated with rmpt() and fitted with brms-generated Stan code: cp 0.601 [0.527, 0.662] (true 0.600), rp 0.752 [0.705, 0.795] (true 0.750); max Rhat 1.0129. All true values covered. On the original branch (2026-09-04) both max Rhats were 1.0071.
  • Ordered constraints. Recommended pattern B: Do 0.772 [0.699, 0.833] (true 0.750), Dn 0.483 [0.080, 0.672] (true 0.550), g 0.358 [0.185, 0.518] (true 0.400). The order holds in all 4000 draws (min difference 0.0083); max Rhat 1.0110. Latent-offset pattern A: 85 of 4000 (2.0%) divergences at default settings, 0 with adapt_delta = 0.99 (max Rhat 1.0060). This moved since 2026-09-04, when pattern A had 263 of 4000 (7.0%) divergences and max Rhat 1.1092 before the refit. Candidate causes are develop's starting step size 0.01 and the exponential sd prior; I did not separate them.
  • update() without newdata. On a real cmdstanr fit of the two-tree 2HTM, update() reproduces the standata exactly. The validation script reports FAIL only because of an extra stancode check. The only differing line is the version header, which update() writes twice. mixture2p shows the same, so this is develop behaviour.

Since 2026-10-08 (develop sync 2, open items)

Gates on 255-mpt-core (e7c8504)

  • Full suite (NOT_CRAN=true): 0 failed, 0 errors, 1 skipped, 7024 passed.
  • check() on an rsync copy: Status OK, 0 errors, 0 warnings, 0 notes.
  • Stancode of 2HTM and PCM identical, and print() of the seven reference models identical except the parameter_info() line of e7c8504, against 9f140ee.
  • document(): no drift. git merge-tree against develop 222a7ef: clean.
  • CI at e7c8504: runs after the push. The last run, at 9f140ee, passed file-size and R-CMD-check on ubuntu, macOS and Windows (release).

Review focus

  • The pipeline methods: check_model.mpt(), check_data.mpt(), check_formula.mpt(), bmf2bf.mpt(), configure_model.mpt().
  • The generated category formulas (.mpt_category_formulas()) and how they meet apply_links().
  • The matched priors and .mpt_match_link_priors().
  • The latent-scale constant() prior in configure_prior.mpt(). It is combined last, so it overrides the probability-scale constant from fixed_pars_priors().
  • revert_check_data.mpt(), which is new since the original branch.
  • .mpt_undo_formula_state() and the update.bmmfit() change (D3): the attribute mpt_bypassed_links, and whether re-checking a stored model on update(formula. = ) changes anything for the other models.

Deviations from the split plan at this level:

  • No test here for the "tree identifier column not present" error. The tree-column tests use the impossible-category fixture and arrive in PR 3.
  • dmpt() and rmpt() keep ....
  • ?mpt refers to the section "Ordered parameter constraints" of the MPT article, which arrives in PR 6.
  • Docs and messages that mention features of later PRs are reworded here. Each later PR restores its own wording.

Known noise

Fits with a non-linear parameter formula print the generic "Non-linear transformations of model parameters detected ... flat priors" warning from validate_default_priors(). The MPT parameters do have priors. PR 2 restricts the warning to formulas that replace a default prior (cf1f259).

Not fixed here

  • nlpar intercept priors. Every mpt parameter is a brms non-linear parameter, and bmm stores its default intercept prior as class = "b", coef = "Intercept". combine_prior() matches priors by exact key, so a user prior given as class = "b", nlpar = "D" does not replace it. This affects mixture2p, mixture3p, imm, m3 and mpt alike, so it belongs in its own issue. ?mpt documents the coef = "Intercept" workaround.

Deferred

Update 2026-10-10: identifiability messaging simplified (48378e8)

This update supersedes the per-case wording described under "Local identifiability" and "Changes since the re-review of the identifiability policy". The rank computation is unchanged.

  • print() ends with one rank line: "locally identified", "the combination of {params} is not identified by the branch expressions", or "Jacobian rank not computed: {reason}". A non-linear formula on a tree parameter is named as not analysed, so a model is never called locally identified without that qualifier.
  • bmm() applies two rules. It warns when the rank is deficient and every formula is intercept-only. Otherwise it gives one fixed message: the branch expressions alone do not identify {params}, predictors or non-linear formulas may identify them, and this was not checked. When no rank can be computed, free parameters beyond the degrees of freedom decide.
  • Behaviour that changes: a non-linear formula on a parameter outside the deficit now gives the message, not the warning. The "zero up to rounding" sentence is gone.
  • Removed code: the tied-rank and rank-deficit text helpers and the logic that decided which parameters a formula reaches.
  • The identifiability paragraph of ?mpt is rewritten (12 lines plus the separate-detection 2HTM example, rank 2 of 3).
  • Tests: every rank-value test is kept. One test per rule (warning, message, rank not computed) replaces the wording tests.
  • Verified at this head: NOT_CRAN=true devtools::test() FAIL 0, WARN 1, SKIP 1, PASS 7043; devtools::check() 0 errors, 0 warnings, 0 notes.

Update 2026-10-13: develop through #546, MPT fits moved to tests/local (head 3bdab5b)

This update replaces the "Gates" section above.

Measured at this head: test-conventions.R and test-model_mpt.R 0 failed (the known warning at test-model_mpt.R:279 remains); the moved fit passes from tests/local; document() no drift. The full suite and R CMD check were run at the tip of the stack (#489).

CI on this head: R-CMD-check passes on macOS, Ubuntu and Windows (Windows check step 9 of its 30 minutes), tests-source passes within the time budget (suite 594 s, cap 800 s), and file-size passes.

@GidonFrischkorn
GidonFrischkorn added this pull request to stack #491 October 3, 2026 11:32
mpt_tree() and mpt(trees, tree_id, links) specify the model; dmpt()/rmpt() give its distribution. Parameters fixed in the formula and non-linear parameter formulas are supported, print() lists the trees and the identifiability bound, and the model meets develop's link, update(), field-order, overview, sd-prior and data-check contracts.
…info() (#255)

develop renamed the internal model_info() to model_docs() (9ea0f9d). ?mpt
called the old name in @details, which document() could not evaluate; the
generated mpt.Rd is unchanged.
…dels (#255)

supported_models() became bmm_models() on develop (#481), and the group list
that the MPT branch had extended in ?supported_models now lives in
?bmm_models.
…un, as the model template does (#255)

use_model_template() (#490) writes version = "NA" for an unversioned model,
as every unversioned constructor on develop stores it, and puts the alias
example under \dontrun{}. MPT stored "" and used the @examplesIf condition
that is never true. model_docs() and model_citation() read "" and "NA"
alike, so nothing printed or returned changes.
…ula get (#255)

Sub-parameters of a non-linear formula live on whatever scale the formula
gives them, so mpt() sets no sd default for them and brms's
student_t(3, 0, 2.5) applies (measured with default_prior(), brms 2.23.0).
The message check_model() prints for them and ?mpt now say so; a test ties
the message to the prior default_prior() reports.
)

check_model.mpt() switches off the link of a parameter with a non-linear
formula and adds that formula's sub-parameters to the model. A checked model
re-used with another formula kept both: update(fit, formula. = ) or
default_prior()/stancode() with fit$bmm$model and a formula in which that
parameter is linear again failed in brms ('The parameter 'a' is not a valid
distributional or non-linear parameter'), and the parameter would have kept
the identity link without a prior.

check_model.mpt() now starts from the constructor's parameters, the links
the parameters had before the bypass (kept as an attribute, so a link the
user switched survives) and the model's latent prior, which it then matches
to a switched link. update.bmmfit() did not call check_model() at all; it
now re-checks the stored model whenever formula. is given, for every model.
Review finding 7 of the phase-1 review.
)

develop defines the expected response per model (#475). On a real 2HTM fit,
posterior_epred() returns draws x rows x categories of expected counts that
sum to the trials of each row and match the counts computed from the latent
parameters to 4e-16; fitted() agrees. The default (defined) is right for mpt.
…essions (#255)

Only num / num was folded, so 1/(2*2), -1/4 and 1/2^2 reached the Stan code as integer division and evaluated to 0 without an error. Every subexpression built from numeric literals with + - * / ^ and parentheses is now folded to its value before the scientific-notation screen runs.
dmpt() read x positionally, so a named vector in another order than the first tree's branches gave a wrong density silently. Names are now matched to the categories; a mismatch errors. An unnamed x is read positionally as before.
The page named no prior for the latent parameters, the random-effect SDs or the correlations.
…y line (#255)

print() of a checked model counted the sub-parameters of non-linear formulas and parameters fixed in the formula as free, so an identified model printed as unidentified. The count now uses the constructor's parameters minus the fixed ones, and the line says the bound holds for intercept-only formulas.
A tree whose branches summed to something other than 1 only warned at mpt()
and the fit then sampled renormalised probabilities. mpt() now errors, and
the message says that restrictions between parameters belong in the branch
expressions.
A branch written as 0 (or folding to 0) passed mpt() and the fit died inside
Stan with a non-finite gradient. mpt() now errors and says that impossible
categories arrive with mpt_tree(impossible = ) in a later release; that
release removes this guard.
#255)

A column missing from the data (or mistyped) in a non-linear formula became a
free parameter with an implicit intercept and only a message. check_model()
now errors and asks for a formula for the new parameter or for the column.
#255)

A missing count silently became 0, and user columns named Y or nTrials were
silently replaced by the response matrix and trial counts bmm builds.
check_data() now warns on missing counts and errors on those two columns.
…se (#255)

Both refused a parameter vector carrying names the selected tree does not
use, which made per-row simulation over trees with different parameters
awkward. Extra names are now ignored; missing ones still error and are
listed first in the message.
The help page never said that the data hold one row per participant and
tree with one count column per response category, named after the branches;
only the example showed it. A paragraph in Details and a comment in the
example now say so.
…es (#255)

summary() reports latent intercepts only and the help page never said how to
read the parameters as probabilities. Details now name native_parameters()
and the example calls it after summary().
The example used (1 | id) per parameter, which gives independent person
effects, and the help page never mentioned (1 |p| id). Details now explain
the shared label and the example uses it.
)

Nothing on the help page said where simulated and expected counts come
from or where a T1 posterior predictive p-value is computed. Details now
name posterior_predict() and posterior_epred() and refer to the MPT article.
…ale (#255)

A fixed g = 0.5 appears as g_Intercept 0.00 under "Constant Parameters" in
summary(), which reads like a different value. The fixed-parameter paragraph
now explains this and names prior_info() and parameter_info() for the
probability.
With several zero branches the message read "The branch probability of 'z'
in tree 'a', 'x' in tree 'b' is the constant 0". The list now comes after a
colon.
…heck (#255)

The four test points cycled four values, so parameters four places apart
in first-appearance order always took the same value, and mpt() accepted
a typo between them (e.g. (1 - A) for (1 - F) in a five-parameter
chain). A golden-ratio sequence now gives each symbol its own value.
The parameters-versus-degrees-of-freedom count missed models that pass the
count but are not identified (the unrestricted Al Hadhrami unitization and
hybrid files: rank 19 of 21 and 21 of 27), and such fits sample the prior
along the redundant directions without any message.

print() now reports the rank of the Jacobian of all category probabilities
(exact derivatives at five deterministic test points) and names the
parameters involved, including parameters that cancel from the likelihood.
check_model() warns for intercept-only formulas and gives a message when
predictors may identify the model across design cells. The count line stays.
The helper had no caller at any level of the stack and was exercised only by a round-trip test. Removing it changes nothing a user can see.
… fits (#255)

The develop sync brought two fixes whose tests name their models by hand:
f9e24ac/383b802e stamp the user's data name after check_data(), which
rebuilds the frame for mpt as for m3, and a0d7ee9 stops reading the
correlation ID of (1 |p| id), the formula ?mpt recommends, as a data
column. Both tests fail on the pre-sync core 8dc5e4c (data_name NULL,
`p` listed as a predictor) and pass on the merged tree.
@GidonFrischkorn GidonFrischkorn changed the title Add multinomial processing tree models (MPT 1/6, #255) Add multinomial processing tree (MPT) models (MPT 1/6, #255) Oct 8, 2026
@GidonFrischkorn GidonFrischkorn changed the title Add multinomial processing tree (MPT) models (MPT 1/6, #255) Add multinomial processing tree models with mpt() (MPT 1/6, #255) Oct 8, 2026
print() of every model ended with 'Use parameter_info() for more details.'
Under the MPT identifiability block this read as a pointer to rank detail,
which parameter_info() does not give. The line now names what it lists:
descriptions, links and fixed values. Raised as C1 in the #489 review.
…526)

Replaces the coef = "Intercept" workaround in ?mpt with the behaviour
after the combine_prior() fix, and tests it on the 2HTM.
The generic carried @Keywords internal without a title, so roxygen wrote
no help page for it.
…nt (#255)

update() re-runs check_model() on the stored model whenever the formula
changes; covers m3 (custom), sdm, imm (regex), ddm, cswald and sdt_rating.
@GidonFrischkorn

Copy link
Copy Markdown
Collaborator Author

Shared-code audit of the MPT stack

I compared develop with the stack head for every model in bmm_models(), 22 cases in all: every version of imm, m3 (ss, cs, custom), sdt_rating (sdt, dpsdt, metad), plus mixture2p, mixture3p (with and without regex), sdm, ddm, cswald, ezdm, sdt_yn, sdt_mafc, sdt_ranking and sdt_cdp. For each case it records print(model), default_prior(), stancode(), the warnings and messages of bmm(backend = "mock"), and the priors, Stan code, parameters and conditions after update() with a changed formula, the same formula and a non-linear formula. Neither branch errors on any case.

The diff has two kinds of lines:

Difference Models Verdict
Last line of print(model): "Use parameter_info() for more details." became "Use parameter_info() for parameter descriptions, links and fixed values." all intended
No "Non-linear transformations ..." warning for b ~ exp(lb) in bmm() and update() m3 (any version) intended, see below

Nothing else differs: priors, Stan code, model parameters and all other warnings and messages are identical, including after update().

update() re-runs check_model(). Re-checking the stored model returns it unchanged and silent for m3 (custom), sdm, imm (regex), ddm, cswald, sdt_rating and mpt, in 1 to 6 ms. The m3 custom default-prior warning appears once at bmm() and not again on update(). New test: check_model() leaves a checked model unchanged and silent.

validate_default_priors(). The warning now fires only when a default prior is dropped. Every free parameter of every built-in model has a default prior, so the two conditions disagree only for parameters that are fixed by default and have none. A non-linear formula on one of those frees it: m3's b, and also mu2/kappa2 (mixtures and imm) and mu (ddm, cswald, ezdm), which are fixed in the same way. Only m3's b was probed. With b ~ exp(lb) the old condition warned, the new one does not. lb gets a flat prior either way, exactly as b does with b ~ 1, which never warned. So the behaviour is now consistent between linear and non-linear formulas. No existing test expectation changed. New test: freeing m3's b through b ~ 1 or b ~ exp(lb) is silent.

print_model_details(). Output of every non-mpt model is unchanged apart from the last line. The generic had @keywords internal but no title, so roxygen wrote no help page; it now has one (@keywords internal developer, like check_model()). It does not collide with #442. #442 should list it in the method contract as an optional method, and export it, since a model in another package cannot register a method for an unexported generic without registerS3method().

Filters (uses_aggregate_data(), response_annotations(), data_column_roles, model_groups, pp_check, conditional_effects). The Poss_<cat> names are exact, anchored matches taken from model$other_vars$indicators$possible, which is NULL for every other model, and check_data() refuses data columns that already carry those names. No user column of another model can be matched.

#526. Fixed in its own PR against develop and merged into 255-mpt-core; ?mpt and the article no longer describe the workaround.

print() gives one rank line: locally identified, the parameters that the
branch expressions leave unidentified, or that the rank was not computed.
Non-linear formulas on tree parameters are named as not analysed.

bmm() warns about a rank deficit when every formula is intercept-only and
otherwise says once, in fixed text, that predictors or non-linear formulas
may identify the parameters and that this was not checked. When no rank can
be computed, free parameters beyond the degrees of freedom decide.
…prior

update() fed every row of the fit's stored prior back in as a user prior,
including the flat default brms fills in for the class of an intercept-only
non-linear parameter. A formula. that added predictors to such a parameter
gave the new effects that flat prior instead of the bmm default. Rows brms
filled in itself are now left out, as brms::update.brmsfit() does.
print() names a parameter with a non-linear formula as not analysed in its
rank line; the identifiability paragraph of ?mpt did not mention it.
@GidonFrischkorn

Copy link
Copy Markdown
Collaborator Author

Update 2026-10-12: review fixes landed, develop merged (head b4a1212)

An independent review of the changes made at this level since the last review (the #526 prior fix, the update() re-check of the model, the simplified identifiability messaging, the print_model_details() generic) found the four changes correct and added two commits:

  • update(fit, formula. = …) gave the coefficients a changed formula introduces a flat prior, and with 0 + x the whole parameter had no prior, because the stored fit's brms default rows were carried as if the user had set them. Fixed in shared code; the Stan code of an updated fit now equals that of a fresh bmm() call for all 25 model cases probed (14 differed before). The same fix is lifted into its own PR against develop, because it reaches mixture2p, mixture3p, imm and m3 independently of MPT.
  • ?mpt says that the rank check does not analyse non-linear formulas, matching what print() reports.

develop (#530, the #526 fix) is merged in. Measured at this head: parse of every R file, document() without drift, byte-identical Stan code for the eight reference models before and after, and the focused test files (mpt, distributions, conditional_effects, pp_check, update, data_check, prior and model helpers) with no failure; the one warning in that run is the pre-existing "flat priors will be used" message of a non-linear-formula test, recorded as a 1.5.0 follow-up.

Deferred to 1.5.0 by the review: storing the user's own prior on the fit so that update() can separate it from defaults in every case; an mpt entry in the stored-frame update tests once the kfold fix (#536) is on develop; log_lik()/posterior_predict() with newdata in the bmm() layout for count models; a silent constant covariate; the misleading flat-prior warning.

The macOS and Windows pull-request jobs ship no CmdStan since #546, so
the conditional_effects() fit opens with skip_unless_stan_tests().
test-conventions.R registered the counts of the stack's tip; a present
file must match, so core registers its one cmdstanr fit and no
set.seed() in test-model_mpt_expressions.R.
ADR 0004 R1: nothing in tests/testthat samples. The fit ran on the
Ubuntu pull-request jobs and would run on every push, compiling on
Windows. tests/local/helper-mpt.R sources the testthat builders so the
trees stay defined once. test-conventions.R drops the MPT entries
registered ahead (test-model_mpt.R, test-ppp.R): the stack adds no
real-backend site.
@GidonFrischkorn

Copy link
Copy Markdown
Collaborator Author

I am merging this PR and the five above it in the stack (#485 to #489) into develop without a second approving review, using the admin override on the branch rule.

Why. I wrote the stack and am the only maintainer who reviewed it, so a required review cannot be met from my side. The stack did get an independent evaluation; its four blocking issues are closed (details in the #489 description). I would rather record the override here than wait on the branch rule.

What was green before the merge:

  • All six heads are pushed and match the PR heads. develop (the test: tier the suite, four workers and a time budget (ADR 0004) #546 merge) is contained in each of them.
  • CI passes on every PR: R-CMD-check on macOS, Ubuntu and Windows, the tests-source job and the file-size check. The Windows check step took 9 to 13 minutes against a 30-minute limit.
  • Full test suite at the article tip with Stan on: 1907 tests, 0 failed, 0 warnings. R CMD check on a clean export: 0 errors, 0 warnings, 0 notes. No document() drift on any level.
  • Parameter recovery for the new models ran at the final head. Scenario 2 has Rhat 1.0129. I accept it as borderline; the script is illustration only.

How the merge runs. Bottom-up, one merge commit per PR (no squash), each level merged only after its own CI is green on the pushed head. #485 to #489 merge under this same decision.

@GidonFrischkorn
GidonFrischkorn merged commit 6f7cde2 into develop Oct 14, 2026
6 checks passed
@GidonFrischkorn
GidonFrischkorn deleted the 255-mpt-core branch October 14, 2026 04:42
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant