diff --git a/NAMESPACE b/NAMESPACE index 42468a082..0f263ce95 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -99,6 +99,8 @@ S3method(expected_response_defined,vwm) S3method(fit_info,brmsfit) S3method(fit_info,brmsfit_list) S3method(fitted,bmmfit) +S3method(generated_parameters,default) +S3method(generated_parameters,mpt) S3method(get_resp_vars,bmmodel) S3method(has_intercept,formula) S3method(has_intercept,terms) @@ -144,9 +146,11 @@ S3method(pp_simulate,ddm) S3method(pp_simulate,default) S3method(pp_simulate,ezdm_3par) S3method(pp_simulate,ezdm_4par) +S3method(ppp,bmmfit) S3method(print,bmm_data_check) S3method(print,bmm_models) S3method(print,bmm_parameters) +S3method(print,bmm_ppp) S3method(print,bmm_prior_info) S3method(print,bmm_sdt_auc) S3method(print,bmm_sdt_latent) @@ -280,6 +284,7 @@ export(mixture2p) export(mixture3p) export(model_citation) export(mpt) +export(mpt_data_reshape) export(mpt_from_eqn) export(mpt_from_string) export(mpt_tree) @@ -297,6 +302,8 @@ export(pp_check) export(pp_check_vars) export(pp_observables) export(pp_simulate) +export(ppp) +export(print_model_details) export(print_pretty_models_md) export(prior_info) export(psdm) diff --git a/NEWS.md b/NEWS.md index ca8b8f91d..3881c4413 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,18 +7,20 @@ * Add **Confidence Rating Signal Detection Theory** (`sdt_rating`) for old/new recognition and detection tasks with ordered confidence ratings, aggregated to counts per rating category. It estimates sensitivity (`d`), the criterion (as `criterion` or `bias`, as in `sdt_yn()`) and the confidence thresholds, and optionally the unequal-variance ratio (`sdratio`); `threshold_type` chooses how the thresholds are parameterized. Also adds `dsdt_rating()` and `rsdt_rating()`. See `?sdt_rating` for the parameters, the threshold types, the default priors and how `d` relates to the sensitivity `sdt_yn()` reports. Thanks to @GidonFrischkorn * Add two versions of `sdt_rating()`, selected with `version`: `"dpsdt"`, the dual-process model (Yonelinas, 1994) with recollection of old items (`Ro`) and recall-to-reject of new items (`Rn`; as in Yonelinas, 2024), and `"metad"`, the meta-d' model, which estimates metacognitive efficiency as the log M-ratio `logmratio` (Maniscalco & Lau, 2012). Both reduce to the standard model while their extra parameters are off. Also adds `mratio()`, `dsdt_dpsdt()`, `rsdt_dpsdt()`, `dsdt_metad()` and `rsdt_metad()`. See `?sdt_rating`. Thanks to @GidonFrischkorn * Add **Continuous Dual-Process Signal Detection Theory** (`sdt_cdp`) for old/new recognition with Remember/Know judgments (Wixted & Mickes, 2010). It estimates familiarity (`dfam`) and recollection (`drec`) sensitivity, the old/new `criterion`, the confidence thresholds and the Remember criterion (`rcrit`) from aggregated response counts; the recollection SD (`sigmar`), the familiarity-recollection correlation (`rho`) and a Know/Guess criterion (`kcrit`) can be freed through the formula. Also adds `aggregate_sdt_cdp_data()`, `dsdt_cdp()` and `rsdt_cdp()`. See `?sdt_cdp` for the parameters, the response format and how `drec` relates to the sensitivity `sdt_yn()` reports. Thanks to @GidonFrischkorn -* Add **Multinomial Processing Tree models** (`mpt`) for tasks with categorical responses that come from discrete processing states, such as the two-high-threshold model of recognition. 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. Parameters take a logit or probit link and can be fixed in the formula (`bmf(g = 0.5)`) or given non-linear formulas. Also adds `dmpt()` and `rmpt()`. See `?mpt` (#255). Thanks to @GidonFrischkorn -* `mpt()` gains `restrictions`, which equates parameters or fixes them to constants in the syntax of MPTinR and TreeBUGS (`c("Dn = Do", "g = 0.5")`), and `simplex`, which declares groups of parameters that sum to 1, such as guessing rates across response options. Blank and `#` comment lines in the restrictions are skipped, so `readLines()` of an MPTinR restrictions file of equalities and constants works. Printing a model lists its restrictions and simplex groups. See `?mpt` (#255). Thanks to @GidonFrischkorn -* By default, a `simplex` group gets a uniform prior over all ways of splitting the probability among its members, so each of K members has prior mean 1/K whatever the order they are listed in. Printing the model states this prior. `?mpt` shows how to get the per-branch uniform priors of TreeBUGS instead, under which the first-listed member has prior mean 0.5 (#255). -* `mpt()` gains `covariates`, for data columns used in branch expressions, such as guessing rates fixed by the design; `bmm()` stops when the branches of a row do not sum to 1 with that row's values. `mpt_tree()` gains `impossible`, for response categories a tree cannot produce, such as distractor responses on trials without distractors. `dmpt()` and `rmpt()` take covariate values too. See `?mpt` and `?mpt_tree` (#255). Thanks to @GidonFrischkorn -* New functions `mpt_from_string()` and `mpt_from_eqn()` build an `mpt()` model from MPTinR-style model text (a string or `readLines()` output) or from an EQN file. Both take `impossible` for response categories a tree cannot produce. `mpt_from_eqn()` renames names `mpt()` cannot use, such as `D_o`, reports the renaming, takes restrictions in the file's names, and maps numeric tree labels through `tree_names`. MPTinR and TreeBUGS skip the first line of an EQN file; bmm reads it unless it is the equation count, so delete a title line or start it with `#`. See `?mpt_from_eqn` (#255). Thanks to @GidonFrischkorn -* `plot()` draws the processing-tree diagram of an `mpt_tree()`, or of every tree of an `mpt()` model, from the branch expressions. See `?plot.mpt_tree` (#255). Thanks to @GidonFrischkorn +* Add **Multinomial Processing Tree models** (`mpt`) for categorical responses from discrete processing states, such as the two-high-threshold model of recognition. Write each tree's branches with `mpt_tree()`, as formulas or as a named list of strings, name a shared parameter alike in every tree, and fit the model to response counts, one row per participant and tree, with `tree_id` naming each row's tree. `links` sets a logit or probit link for all or single parameters; parameters can be fixed in the formula (`bmf(g = 0.5)`) or given non-linear formulas. Also adds `dmpt()` and `rmpt()`. See `?mpt` (#255). Thanks to @GidonFrischkorn +* `mpt()` takes `restrictions` that equate parameters or fix them to constants, in MPTinR and TreeBUGS syntax (`c("Dn = Do", "g = 0.5")`, or `readLines()` of such a restrictions file), and `simplex` for parameters that sum to 1, such as guessing rates. A simplex group's default uniform prior gives each of K members prior mean 1/K whatever their order; `?mpt` shows how to get the per-branch uniform priors of TreeBUGS, which give the first member prior mean 0.5. Printing the model lists the restrictions, the simplex groups with their prior and an identifiability check (#255). Thanks to @GidonFrischkorn +* Branch expressions can use data columns declared in `covariates`, for probabilities the design fixes, such as guessing rates that depend on the number of response options; `bmm()` stops when a row's values make the branches of a tree not sum to 1. Response categories a tree cannot produce, such as distractor responses on trials without distractors, are declared with `impossible` in `mpt_tree()` and in the importers. `dmpt()` and `rmpt()` take covariate values too. See `?mpt_tree` (#255). Thanks to @GidonFrischkorn +* New functions `mpt_from_string()` and `mpt_from_eqn()` build an `mpt()` model from MPTinR-style model text (a string or `readLines()` output) or from an EQN file (multiTree, MPTinR, TreeBUGS); both take `restrictions` as their second argument. `mpt_from_eqn()` renames and reports names `mpt()` cannot use, such as `D_o`; it takes restrictions in the file's names, maps per-tree category labels onto shared categories with `categories`, and numeric tree labels with `tree_names`. Unlike MPTinR and TreeBUGS, bmm reads an EQN file's first line unless it is the equation count, so start a title line with `#` or delete it. See `?mpt_from_eqn` (#255). Thanks to @GidonFrischkorn +* New function `mpt_data_reshape()` reshapes the wide frequency tables of MPTinR and TreeBUGS, one row per participant and one column per category, into the one-row-per-participant-and-tree layout `bmm()` expects, with 0 for the categories a tree cannot produce. Columns that are no counts, such as a group factor, are carried to every row; `columns = "mptinr"` reads the unnamed, sorted columns of MPTinR data. See `?mpt_data_reshape`. Thanks to @GidonFrischkorn +* `plot()` draws the processing-tree diagram of an `mpt_tree()`, or of every tree of an `mpt()` model, from the branch expressions, and marks restricted parameters on their edges. See `?plot.mpt_tree` (#255). Thanks to @GidonFrischkorn +* New article *Multinomial Processing Tree (MPT) models* walks through specifying, fitting and interpreting `mpt()` models, constraining parameters, design-fixed covariates, a fit to real data, non-linear parameters, importing model files and migrating from MPTinR or TreeBUGS (#255). Thanks to @GidonFrischkorn ### New datasets * Add **`broeder_schuetz_2009_e3`**, binary old/new recognition data from Broeder & Schuetz (2009, Exp. 3), with five base-rate conditions from 40 subjects. See `?broeder_schuetz_2009_e3`. * Add **`meyer_grant_jakob_2025`**, ranking data from Meyer-Grant & Jakob (2025): 60 subjects ranking 3, 4 or 5 faces, aggregated to the target's rank counts per subject and set size. See `?meyer_grant_jakob_2025`. ### New features +* New function `ppp()` computes the goodness-of-fit statistics T1 (category frequencies) and T2 (their covariances across participants) of Klauer (2010) with their posterior predictive p-values, the checks TreeBUGS' `PPP()` reports, for `mpt()`, `m3()` and the SDT count models. `level = "participant"` gives one T1 per participant. See `?ppp` (#497). * `m3()` accepts `choice_rule = "gaussian"`, a Thurstonian rule: every candidate gets independent normal noise and the strongest one is chosen (`"softmax"` is the same with Gumbel noise). Activation estimates are on a smaller scale than under `"softmax"`, and conclusions about `c` against `a`, or about conditions that change the number of candidates such as set size, can differ between the rules. Fits take much longer than with `"softmax"`. The Gaussian rule requires whole numbers of options. See `choice_rule` in `?m3`. Thanks to @GidonFrischkorn * New function `bmm_setup()` checks whether your machine can fit models: the C++ toolchain, `cmdstanr` with CmdStan, `rstan`, and the backend `bmm()` will use. It then compiles and samples a small test model and gives one fix for every check that failed. When `bmm()` fails because the compiler or the Stan backend is missing or broken, the error now points to `bmm_setup()` (#431). * New post-processing functions for the signal detection models: `roc_sdt()`, `roc_observed()` and `auc_sdt()` return the model-implied and observed ROC and the area under it for `sdt_yn` and `sdt_rating`; `latent_sdt()` the latent evidence distributions of all four models; `sdt_thresholds()` the confidence thresholds of `sdt_rating`; and `sdt_sensitivity()` the sensitivity on the noise, signal or root-mean-square SD scale. All but `roc_observed()` print a summary; `plot()` draws the ROC with the observed points, the AUC and the latent distributions; `summary()` says when `d` is *d_a*. The new article *Signal Detection Theory models* walks through all four models and these functions (#370). Thanks to @GidonFrischkorn diff --git a/R/bmmformula.R b/R/bmmformula.R index aecc15ee2..501287efa 100644 --- a/R/bmmformula.R +++ b/R/bmmformula.R @@ -251,11 +251,24 @@ bmf2bf.default <- function(model, formula) { NULL } +# The parameters whose formulas a model's check_formula method generates, so +# that no user formula is expected for them. Their placeholders are still added +# below, because the generating method overwrites them by name. +generated_parameters <- function(model) { + UseMethod("generated_parameters") +} + +#' @exportS3Method +generated_parameters.default <- function(model) { + character(0) +} + add_missing_parameters <- function(model, formula = bmmformula()) { formula_pars <- names(formula) model_pars <- setdiff(names(model$parameters), unused_alternative_parameters(model, formula)) fixed_pars <- names(model$fixed_parameters) + generated_pars <- generated_parameters(model) estimable_formula_pars <- setdiff(formula_pars, fixed_pars) missing_pars <- setdiff(model_pars, estimable_formula_pars) @@ -264,7 +277,7 @@ add_missing_parameters <- function(model, formula = bmmformula()) { add <- stats::as.formula(paste(mpar, "~ 1")) if (mpar %in% fixed_pars) { attr(add, "constant") <- TRUE - } else { + } else if (!mpar %in% generated_pars) { message2("No formula for parameter {mpar} provided. Only a fixed intercept will be estimated.") } formula[mpar] <- add diff --git a/R/distributions.R b/R/distributions.R index fc70851a7..2b205fae1 100644 --- a/R/distributions.R +++ b/R/distributions.R @@ -846,7 +846,11 @@ rm3 <- function(n, size, pars, m3_model, act_funs = NULL, unpack = FALSE, #' categories, giving the number of observed responses per category. If #' `x` is a named vector, its names are matched to the response categories #' and the vector is reordered accordingly; an unnamed vector is read in -#' the order of the branch names of the selected tree. +#' the order of the model's response categories, +#' `mpt_model$resp_vars$resp_cats`, which is the order of the first tree +#' given to [mpt()] ([mpt()] stores the branches of every tree in that +#' order, whatever order they were written in). `rmpt()` returns its +#' columns in the same order. #' @param n Integer. Number of observations to generate data for. #' @param size The total number of observations across all response categories. #' @param pars A named vector or list with the values of at least the latent @@ -884,23 +888,23 @@ rm3 <- function(n, size, pars, m3_model, act_funs = NULL, unpack = FALSE, #' and `rmpt` generates random response counts for the response categories. #' #' @examples -#' tree_old <- mpt_tree("old", list( +#' tree_target <- mpt_tree("target", list( #' old = "D + (1 - D) * g", #' new = "(1 - D) * (1 - g)" #' )) -#' tree_new <- mpt_tree("new", list( +#' tree_lure <- mpt_tree("lure", list( #' old = "(1 - D) * g", #' new = "D + (1 - D) * (1 - g)" #' )) -#' model <- mpt(list(tree_old, tree_new), tree_id = "item_type") +#' model <- mpt(list(tree_target, tree_lure), tree_id = "item_type") #' #' dmpt( #' x = c(35, 15), pars = c(D = 0.7, g = 0.5), -#' mpt_model = model, tree = "old" +#' mpt_model = model, tree = "target" #' ) #' rmpt( #' n = 10, size = 50, pars = c(D = 0.7, g = 0.5), -#' mpt_model = model, tree = "new" +#' mpt_model = model, tree = "lure" #' ) #' @export dmpt <- function(x, pars, mpt_model, tree = NULL, log = TRUE, diff --git a/R/helpers-model.R b/R/helpers-model.R index 19363ae96..964dff841 100644 --- a/R/helpers-model.R +++ b/R/helpers-model.R @@ -486,6 +486,7 @@ print.bmmodel <- function(x, ...) { #' #' @return Called for its printed output; returns `NULL` invisibly #' @keywords internal developer +#' @export print_model_details <- function(model, ...) { UseMethod("print_model_details") } diff --git a/R/model_mpt.R b/R/model_mpt.R index 32ef20fcf..6a2086e1c 100644 --- a/R/model_mpt.R +++ b/R/model_mpt.R @@ -12,18 +12,34 @@ #' #' @param name Character. Label for the tree. For models with multiple trees, #' the label identifies the tree in the data: the column named by the -#' `tree_id` argument of [mpt()] holds one such label per observation. -#' @param branches A named list of character strings. Each element gives the -#' branch probability expression for one response category, and the element -#' names are the response categories. The expressions can use latent -#' parameters (e.g., `"D + (1 - D) * g"`), numeric constants, and declared -#' covariates (see [mpt()]). +#' `tree_id` argument of [mpt()] holds one such label per observation. A +#' tree stands for one kind of trial, such as targets or lures in a +#' recognition test, not for a response; the responses are the branch +#' names. +#' @param branches Either a named list of character strings or the first of +#' several formulas. Each list element (or formula) gives the branch +#' probability expression for one response category; the element names (or +#' the left-hand sides) are the response categories. The expressions can use +#' latent parameters (e.g., `D + (1 - D) * g`), numeric constants, and +#' declared covariates (see [mpt()]). +#' @param ... Further formulas of the form `category ~ expression`, one per +#' response category, when `branches` is a formula; they take no argument +#' names. A call uses either +#' formulas or a named list of strings, not a mix; both spellings give the +#' same tree (see the examples). #' @param impossible Character vector. Response categories that cannot occur in #' this tree, for example distractor responses when no distractors were #' presented. Such categories get no branch expression, and the remaining -#' branches must sum to 1 on their own. +#' branches must sum to 1 on their own. It follows `...`, so it has to be +#' named in full. #' -#' @details Constant arithmetic such as `1/4` or `1/(2*2)` is folded into a +#' @details Nothing in the formulas is evaluated: the left-hand side must be a +#' single response category name, the right-hand side is kept as an +#' expression (with constant arithmetic folded, see below), and the formula +#' environment is dropped. A one-sided formula, a left-hand side that is not +#' a name, and a mix of formulas and strings are errors. +#' +#' Constant arithmetic such as `1/4` or `1/(2*2)` is folded into a #' decimal literal (`0.25`) when the tree is created. Stan compiles a bare #' integer fraction as integer division (`1/4 == 0`), which would silently #' corrupt the likelihood. @@ -37,37 +53,43 @@ #' #' @return An object of class `mpt_tree` #' -#' @keywords transform +#' @keywords bmmodel #' #' @examples -#' tree_old <- mpt_tree( -#' name = "old", +#' # the target tree of the two-high-threshold model: targets (studied items) +#' # end in the responses "old" and "new" +#' tree_target <- mpt_tree( +#' "target", +#' old ~ D + (1 - D) * g, +#' new ~ (1 - D) * (1 - g) +#' ) +#' tree_target +#' +#' # the same tree with the branches as a named list of strings +#' tree_target_strings <- mpt_tree( +#' name = "target", #' branches = list( #' old = "D + (1 - D) * g", #' new = "(1 - D) * (1 - g)" #' ) #' ) -#' tree_old +#' identical(tree_target, tree_target_strings) #' #' # a three-alternative recognition task in which lures are shown on some #' # trials only: `lure` responses cannot occur on the others #' tree_lures <- mpt_tree( -#' name = "lures", -#' branches = list( -#' target = "D * b + (1 - D) * (1/3)", -#' lure = "D * (1 - b) + (1 - D) * (1/3)", -#' new = "(1 - D) * (1/3)" -#' ) +#' "lures", +#' target ~ D * b + (1 - D) * (1/3), +#' lure ~ D * (1 - b) + (1 - D) * (1/3), +#' new ~ (1 - D) * (1/3) #' ) #' #' # the branches of the lure-free tree sum to 1 without the `lure` category; #' # no probability is reserved for it #' tree_nolures <- mpt_tree( -#' name = "nolures", -#' branches = list( -#' target = "D + (1 - D) * (1/2)", -#' new = "(1 - D) * (1/2)" -#' ), +#' "nolures", +#' target ~ D + (1 - D) * (1/2), +#' new ~ (1 - D) * (1/2), #' impossible = "lure" #' ) #' tree_nolures @@ -77,16 +99,50 @@ #' # e.g. bmf(D ~ 0 + cond) #' model <- mpt(list(tree_lures, tree_nolures), tree_id = "tree") #' @export -mpt_tree <- function(name, branches, impossible = NULL) { +mpt_tree <- function(name, branches, ..., impossible = NULL) { stop_missing_args() stopif( !is.character(name) || length(name) != 1L || !nzchar(name), "The tree name must be a single non-empty character string." ) + stopif( + any(!nzchar(as.character(as.list(substitute(list(...)))[-1]))), + "An argument of mpt_tree() is empty, probably from a trailing comma after \\ + the last branch." + ) + extra <- list(...) + # name, branches and impossible match partially, so b = b ~ y would take the + # slot of branches and reorder the categories, and imp = binds to ... + inexact_names <- setdiff(names(sys.call()), c("", "name", "branches", "impossible")) + stopif( + length(extra) > 0 && length(inexact_names) > 0, + "Write the arguments name, branches and impossible out in full when \\ + several branches are given, and give the formulas after branches without \\ + a name: the left-hand side of a formula is its response category. \\ + Offending: {collapse_comma(inexact_names)}" + ) + from_formulas <- inherits(branches, "formula") + mixed <- length(extra) > 0 && + (!from_formulas || !all(vapply(extra, inherits, logical(1), "formula"))) + stopif( + mixed, + "Cannot mix formulas and other branch expressions in tree '{name}'. Write \\ + every branch as a formula, category ~ expression, or pass a named list of \\ + character strings. The impossible argument has to be named, \\ + impossible = \"c\"." + ) + stopif( + is.list(branches) && any(vapply(branches, inherits, logical(1), "formula")), + "Pass the branch formulas as separate arguments, \\ + mpt_tree(name, a ~ expr_a, b ~ expr_b), not as a list." + ) + if (from_formulas) { + branches <- .mpt_formula_branches(c(list(branches), extra), name) + } stopif( !is.list(branches) || length(branches) == 0L, "The branches must be a named list with one branch probability expression \\ - per response category." + per response category, or formulas of the form category ~ expression." ) resp_cats <- names(branches) stopif( @@ -116,15 +172,25 @@ mpt_tree <- function(name, branches, impossible = NULL) { ) branches[] <- lapply(resp_cats, function(resp_cat) { expr <- branches[[resp_cat]] + if (from_formulas) { + parsed <- expr + expr <- deparse1(parsed) + } else { + stopif( + !is.character(expr) || length(expr) != 1L || !nzchar(expr), + "Branch probability expressions must be single character strings \\ + (tree '{name}', category '{resp_cat}')." + ) + parsed <- try(str2lang(expr), silent = TRUE) + stopif( + is_try_error(parsed), + "Cannot parse the branch expression '{expr}' in tree '{name}'." + ) + } stopif( - !is.character(expr) || length(expr) != 1L || !nzchar(expr), - "Branch probability expressions must be single character strings \\ - (tree '{name}', category '{resp_cat}')." - ) - parsed <- try(str2lang(expr), silent = TRUE) - stopif( - is_try_error(parsed), - "Cannot parse the branch expression '{expr}' in tree '{name}'." + !is.language(parsed) && !(is.numeric(parsed) && length(parsed) == 1L), + "The branch of category '{resp_cat}' in tree '{name}' must be an \\ + expression in parameters and numbers, such as D + (1 - D) * g, not {expr}." ) parsed <- .mpt_fold_numeric_division(parsed) scientific <- .mpt_scientific_constants(parsed) @@ -166,6 +232,9 @@ print.mpt_tree <- function(x, ...) { simplex = NULL, restrictions = NULL, unrestricted_trees = NULL, links = "logit", default_priors = NULL, call = NULL, ...) { + # a named list sets single parameters on top of the model-wide default link + parameter_links <- if (is.list(links)) links + links <- if (is.list(links)) "logit" else links trees <- .mpt_as_tree_list(trees) if (length(trees)) names(trees) <- vapply(trees, `[[`, character(1), "name") covariates <- covariates %||% character(0) @@ -275,7 +344,8 @@ print.mpt_tree <- function(x, ...) { call = call ) out$default_priors[names(default_priors)] <- default_priors - set_links(out, NULL) + out <- set_links(out, parameter_links) + .mpt_match_link_priors(out) } # matched priors: effects get the intercept's own prior, because a narrower @@ -321,9 +391,21 @@ settable_links.mpt <- function(model) { names(defaults)[defaults != "identity"] } +.mpt_link_functions <- c("logit", "probit") + #' @exportS3Method settable_link_functions.mpt <- function(model) { - c("logit", "probit") + .mpt_link_functions +} + +# the derived member of each simplex group and the stick-breaking components +# get their formulas from check_formula.mpt(), so none is expected from the user +#' @exportS3Method +generated_parameters.mpt <- function(model) { + c( + vapply(model$other_vars$simplex, function(grp) grp[length(grp)], character(1)), + unname(model$other_vars$simplex_raw) + ) } # user facing alias @@ -336,8 +418,9 @@ settable_link_functions.mpt <- function(model) { #' @description #' Multinomial Processing Tree (MPT) models measure the probabilities of #' discrete latent cognitive states from categorical responses. A model is -#' specified as a set of trees (one per distinct branch structure) created -#' with [mpt_tree()]. Branch expressions that terminate in the same response +#' specified as a set of trees created with [mpt_tree()], one per distinct +#' branch structure, such as one tree for targets and one for lures in a +#' recognition test. Branch expressions that terminate in the same response #' category are summed within each tree, and parameters shared across trees #' are equated by giving them the same name. All latent probability #' parameters are estimated on an unconstrained latent scale (logit or @@ -350,7 +433,7 @@ settable_link_functions.mpt <- function(model) { #' @param tree_id Character. Name of the data column whose values name the tree #' each observation belongs to; the values must match the tree names. Can be #' omitted for single-tree models. Effects of experimental conditions belong -#' in the parameter formulas, not here — see Details. +#' in the parameter formulas, not here; see Details. #' @param covariates Character vector. Names of data columns that appear in #' branch expressions but are not latent parameters, for example #' design-fixed guessing rates. Covariates pass into the model formulas @@ -383,12 +466,17 @@ settable_link_functions.mpt <- function(model) { #' restrictions file, in `bmm` it holds the restrictions. To use such a file, #' pass `readLines(path)`; blank lines and lines starting with `#` are #' skipped, and a `#` after a restriction starts a comment, as in MPTinR. -#' Restrictions are substituted into the branch expressions before the -#' parameters are identified, so a restricted parameter is not part of the -#' model and cannot be a simplex member. Order constraints (`"Do > Dn"`) are -#' not supported here; see Details. -#' @param links Character. The link function for all latent probability -#' parameters: `"logit"` (default) or `"probit"`. +#' A restricted parameter is removed from the model; the Details say how +#' this differs from fixing a parameter in the formula. Order constraints +#' (`"Do > Dn"`) are not supported here; see Details. +#' @param links The link function of the latent probability parameters: +#' `"logit"` (default) or `"probit"`. A single string sets every parameter; +#' a named list, `links = list(D = "probit")`, sets the named parameters +#' and leaves the others at `"logit"`. Members of a `simplex` group and +#' their `praw` components take no link of their own. A parameter switched +#' to the other link, in the list or afterwards with +#' `model$links$D <- "probit"`, takes that link's default priors (see +#' Details). #' @param ... used internally for testing, ignore it #' #' @details `r model_docs(.model_mpt(), components = c('domain', 'task', 'name', 'citation'))` @@ -399,11 +487,13 @@ settable_link_functions.mpt <- function(model) { #' like the branch names given to [mpt_tree()] (`old` and `new` in the #' example below); `bmm()` finds them by these names, so no response #' variable is passed. The `tree_id` column names the tree of each row, and -#' the counts of a row must come from trials of that tree. +#' the counts of a row must come from trials of that tree. Wide tables with +#' one column per category of an EQN file, as MPTinR and TreeBUGS take them, +#' are reshaped with [mpt_data_reshape()]. #' -#' A separate tree is needed only when the *branch expressions* differ — that +#' A separate tree is needed only when the *branch expressions* differ, that #' is, when the trial determines which latent processes apply or which -#' response categories are reachable (old versus new probes, a response set +#' response categories are reachable (targets versus lures, a response set #' with or without distractors, inclusion versus exclusion instructions). #' Experimental factors whose *effect on the parameters* you want to estimate #' belong in the parameter formulas, exactly as in every other `bmm` model: @@ -412,23 +502,36 @@ settable_link_functions.mpt <- function(model) { #' well. When several levels of a factor share a branch structure, add a #' column naming the tree and keep the factor for the formulas. #' -#' Parameters that receive a non-linear predictor formula (a formula whose -#' right-hand side references other formula parameters, e.g. -#' `D ~ Dmax * (1 - exp(-rate * ptime))`) are not transformed by the link -#' function: the user-supplied expression must keep the parameter within -#' (0, 1). Sub-parameters of such formulas (e.g., `Dmax` and `rate`) are -#' estimated on the identity scale with `normal(0, 1)` default priors. Their -#' random-effect SDs keep brms's default prior, `student_t(3, 0, 2.5)`, -#' because no single rate fits every scale such a formula can give them. +#' **Equating and fixing parameters.** Parameters are equated across trees +#' by giving them the same name. Beyond that, there are two mechanisms, and +#' they differ in when the constraint is decided. A *restriction* +#' (`restrictions = c("Dn = Do", "g = 0.5")`) is applied when the model is +#' built: it is substituted into the branch expressions before the +#' parameters are identified, so the restricted parameter is not part of +#' the model, gets no formula and no prior, cannot be a simplex member, and +#' is listed by `print(model)` under `Restrictions:` but not by +#' [parameter_info()]. A *constant in the formula* (`bmf(D ~ 1, g = 0.5)`) +#' is applied at fit time: the parameter stays part of the model, is listed +#' by `print()` under `Fixed:` and by [parameter_info()] with its value, and +#' is freed again by giving it a formula. The value is mapped to the latent +#' scale when the constant prior is built, so [default_prior()] shows +#' `constant(0)` for a guessing rate of 0.5 under the logit link, and +#' `summary()` lists it as `g_Intercept 0.00` under "Constant Parameters"; +#' [prior_info()] and [parameter_info()] show the probability, 0.5. Members +#' of a `simplex` group and their `praw` components cannot be fixed this +#' way; fix parameters outside the group instead. As a rule of thumb, a +#' constant that defines the model goes into `restrictions`, a constant you +#' compare against a free version goes into the formula; the section +#' "Constraining parameters" of the MPT article walks through both. #' -#' Under the logit link, intercepts and regression coefficients get a -#' `logistic(0, 1)` prior and random-effect SDs get `exponential(1)`. Under -#' the probit link, intercepts and coefficients get `normal(0, 1)` and SDs -#' get `exponential(2)`. Group-level correlation matrices get `lkj(2)`. -#' Coefficients of the sub-parameters above get `normal(0, 0.5)`. A parameter -#' switched to the other link after construction (`model$links$D <- -#' "probit"`) takes that link's priors unless its default prior was -#' replaced. [default_prior()] lists the rows a given formula produces. +#' **Default priors.** Under the logit link, intercepts and regression +#' coefficients get a `logistic(0, 1)` prior and random-effect SDs get +#' `exponential(1)`. Under the probit link, intercepts and coefficients get +#' `normal(0, 1)` and SDs get `exponential(2)`. Group-level correlation +#' matrices get `lkj(2)`. A parameter switched to the other link after +#' construction (`model$links$D <- "probit"`) takes that link's priors +#' unless its default prior was replaced. [default_prior()] lists the rows +#' a given formula produces. #' #' The default prior of a simplex group of K members is the uniform #' (symmetric Dirichlet(1)) distribution on the simplex: each member has @@ -451,32 +554,22 @@ settable_link_functions.mpt <- function(model) { #' factor, leaving their ratios unchanged. Read such effects on the probability scale with #' [native_parameters()]. #' -#' A parameter can also be fixed to a probability at fit time, in the -#' formula: `bmf(D ~ 1, g = 0.5)`. Unlike a restriction, the parameter stays -#' part of the model and can be freed again by giving it a formula. The -#' value is mapped to the latent scale when the constant prior is built, so -#' [default_prior()] shows `constant(0)` for a guessing rate of 0.5 under -#' the logit link. For the same reason, `summary()` lists a fixed `g = 0.5` -#' as `g_Intercept 0.00` under "Constant Parameters"; [prior_info()] and -#' [parameter_info()] show the probability, 0.5. Members of a `simplex` -#' group and their `praw` components cannot be fixed this way; fix parameters -#' outside the group instead. -#' -#' Printing the model ends with an identifiability check for intercept-only -#' formulas. It compares the number of free parameters with the degrees of -#' freedom (response categories minus 1, summed over trees) and computes the -#' rank of the Jacobian of the category probabilities with respect to the -#' free parameters, from exact derivatives at five interior test values (and -#' five test values of each covariate); parameters fixed in the formula enter -#' at their values. A rank below the number of free parameters means that a -#' combination of the named parameters is not identified by the branch -#' expressions, even when the count passes. The check is local: it holds at -#' the test values, not at the boundaries of the parameter space. A test -#' value with a non-finite derivative is left out, and the printout says how -#' many; if [stats::D()] cannot differentiate a branch expression, the rank -#' is not computed and the count decides. The rank sees the tree parameters -#' only: after a non-linear formula, the rank line names that formula's -#' parameter as not analysed. +#' **Identifiability.** Printing the model ends with an identifiability +#' check for intercept-only formulas. It compares the number of free +#' parameters with the degrees of freedom (response categories minus 1, +#' summed over trees) and computes the rank of the Jacobian of the category +#' probabilities with respect to the free parameters, from exact derivatives +#' at five interior test values (and five test values of each covariate); +#' parameters fixed in the formula enter at their values. A rank below the +#' number of free parameters means that a combination of the named +#' parameters is not identified by the branch expressions, even when the +#' count passes. The check is local: it holds at the test values, not at the +#' boundaries of the parameter space. A test value with a non-finite +#' derivative is left out, and the printout says how many; if [stats::D()] +#' cannot differentiate a branch expression, the rank is not computed and +#' the count decides. The rank sees the tree parameters only: after a +#' non-linear formula, the rank line names that formula's parameter as not +#' analysed. #' #' `bmm()` warns when the rank is deficient and every formula is #' intercept-only. When a formula has a predictor or is non-linear, it only @@ -487,10 +580,10 @@ settable_link_functions.mpt <- function(model) { #' `Do`, `Dn` and `g` cannot be estimated; a second condition that varies #' the response bias identifies them. #' -#' `summary()` reports intercepts and regression coefficients on the latent -#' (logit or probit) scale. [native_parameters()] returns the posterior -#' draws of every parameter on the probability scale for each observed -#' combination of the predictors, e.g. +#' **Reading the estimates.** `summary()` reports intercepts and regression +#' coefficients on the latent (logit or probit) scale. [native_parameters()] +#' returns the posterior draws of every parameter on the probability scale +#' for each observed combination of the predictors, e.g. #' `native_parameters(fit, re_formula = NA, summary = TRUE)`; a difference #' between conditions is the difference of these draws. With #' `re_formula = NA` these are the probabilities of the median participant @@ -505,20 +598,31 @@ settable_link_functions.mpt <- function(model) { #' #' For posterior predictive checks, [brms::posterior_predict()] returns #' simulated counts (draws x rows x categories) and [brms::posterior_epred()] -#' the expected counts; the MPT article computes the T1 statistic and its -#' posterior predictive p-value from them. With `newdata`, +#' the expected counts; [ppp()] computes the T1 and T2 statistics of Klauer +#' (2010) and their posterior predictive p-values from them. With `newdata`, #' [brms::posterior_predict()] needs the columns `nTrials`, `Idx_` and #' `Poss_` that `check_data()` builds; `fit$data` contains them. #' +#' **Non-linear formulas.** Parameters that receive a non-linear predictor +#' formula (a formula whose right-hand side references other formula +#' parameters, e.g. `D ~ Dmax * (1 - exp(-rate * ptime))`) are not +#' transformed by the link function: the user-supplied expression must keep +#' the parameter within (0, 1). Sub-parameters of such formulas (e.g., +#' `Dmax` and `rate`) are estimated on the identity scale with +#' `normal(0, 1)` default priors on their intercepts and `normal(0, 0.5)` +#' on their regression coefficients. Their random-effect SDs keep brms's +#' default prior, `student_t(3, 0, 2.5)`, because no single rate fits every +#' scale such a formula can give them. +#' #' Order constraints between parameters (`Do > Dn`) are expressed by #' reparameterizing the larger parameter in the model formula, e.g. #' `bmf(Do ~ Dn + (1 - Dn) * inv_logit(phi), Dn ~ 1, phi ~ 1)`; the section #' "Order constraints" of the MPT article walks through the #' recipe. #' -#' Parameter and response category names must start with a letter and may -#' contain only letters and digits, because brms does not allow underscores -#' or dots in non-linear parameter names. +#' **Names.** Parameter and response category names must start with a +#' letter and may contain only letters and digits, because brms does not +#' allow underscores or dots in non-linear parameter names. #' #' All latent parameters are non-linear parameters of the underlying brms #' model, and brms has no separate `Intercept` prior class for those: their @@ -535,17 +639,19 @@ settable_link_functions.mpt <- function(model) { #' @examples #' \dontrun{ #' # two-high-threshold (2HTM) model of recognition memory -#' 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)" -#' )) +#' # one tree per probe type: targets (studied items) and lures (new items), +#' # each ending in the responses "old" and "new" +#' tree_target <- mpt_tree("target", +#' old ~ D + (1 - D) * g, +#' new ~ (1 - D) * (1 - g) +#' ) +#' tree_lure <- mpt_tree("lure", +#' old ~ (1 - D) * g, +#' new ~ D + (1 - D) * (1 - g) +#' ) #' #' model <- mpt( -#' trees = list(tree_old, tree_new), +#' trees = list(tree_target, tree_lure), #' tree_id = "item_type" #' ) #' @@ -553,7 +659,7 @@ settable_link_functions.mpt <- function(model) { #' # one row per participant x tree; count columns named after the branches #' data <- data.frame( #' id = rep(1:20, each = 2), -#' item_type = rep(c("old", "new"), 20) +#' item_type = rep(c("target", "lure"), 20) #' ) #' counts <- t(sapply(data$item_type, function(cond) { #' rmpt( @@ -587,11 +693,37 @@ mpt <- function(trees, tree_id = NULL, covariates = NULL, simplex = NULL, restrictions = NULL, links = "logit", ...) { call <- match.call() stop_missing_args() - links <- match.arg(links, c("logit", "probit")) + # NULL means the default, as in the other model constructors + links <- links %||% "logit" + # a named list is validated against the parameters once the model exists + stopif( + !is.list(links) && !is.character(links), + "The links argument must be a single link function for all latent \\ + parameters, links = \"probit\", or a named list with one link per \\ + parameter, links = list(D = \"probit\")." + ) + stopif( + is.character(links) && length(links) != 1L, + "The links argument takes one link function for the whole model. Links \\ + for single parameters go in a named list, links = list(D = \"probit\")." + ) + stopif( + is.character(links) && !links %in% .mpt_link_functions, + "Unknown link function(s): {collapse_comma(links)}. mpt() takes \\ + {collapse_comma(.mpt_link_functions)}" + ) trees <- .mpt_as_tree_list(trees) covariates <- covariates %||% character(0) simplex <- .mpt_as_simplex_list(simplex) restrictions <- .mpt_parse_restrictions(restrictions) + # the link matcher would repair a removed name to a parameter one edit away + restricted_links <- intersect(names(links), names(restrictions)) + stopif( + length(restricted_links) > 0, + "The link target(s) {collapse_comma(restricted_links)} are removed by the \\ + restrictions and take no link. Set the link of the parameter a \\ + restriction equates them to instead." + ) stopif( length(trees) == 0L || !all(vapply(trees, inherits, logical(1), "mpt_tree")), @@ -821,13 +953,14 @@ print_model_details.mpt <- function(model, ...) { "\n" ) } - # the links field lists the sticks as identity because their generated - # formulas apply the model link; the display says which one + # the links field lists the members and sticks as identity because their + # generated formulas apply the model link; the display says which one for (grp in model$other_vars$simplex) { cat(glue( "Simplex: {collapse_comma(grp)} via stick-breaking; the components \\ {collapse_comma(model$other_vars$simplex_raw[grp[-length(grp)]])} pass \\ - through {inv_link('x', model$other_vars$link)[[1]]}()" + through {inv_link('x', model$other_vars$link)[[1]]}() in the generated \\ + formulas (hence the identity links above)" ), "\n") cat(strrep(" ", 12), glue( "default prior: symmetric Dirichlet(1), prior mean 1/{length(grp)} = \\ @@ -1273,7 +1406,7 @@ check_data.mpt <- function(model, data, formula) { unused_trees <- setdiff(tree_names, unique(tree_values)) warnif( length(unused_trees) > 0, - "The data contain no observations for tree(s): {collapse_comma(unused_trees)}" + "The data contain no rows for tree(s): {collapse_comma(unused_trees)}" ) # brms needs the indicators as data columns, so they are handed on through # the data rather than through the attribute bridge @@ -1362,7 +1495,7 @@ check_data.mpt <- function(model, data, formula) { stopif( length(observed) > 0, "The response category '{resp_cat}' is declared impossible for \\ - {length(observed)} observation(s), but the data record responses in it \\ + {length(observed)} row(s), but the data record responses in it \\ (first: row {observed[1]}, count {data$Y[observed[1], resp_cat]}). \\ Please check the impossible argument of mpt_tree() and the response counts." ) diff --git a/R/model_mpt_data.R b/R/model_mpt_data.R new file mode 100644 index 000000000..2e4986751 --- /dev/null +++ b/R/model_mpt_data.R @@ -0,0 +1,333 @@ +#' Reshape wide MPT frequency tables to the long format of `bmm()` +#' +#' MPTinR and TreeBUGS take one row per participant and one column per +#' response category of the EQN file. [bmm()] takes one row per participant +#' and tree, a column naming the tree, and one count column per *shared* +#' response category (see [mpt()]). `mpt_data_reshape()` reshapes the first +#' layout into the second for a given model. +#' +#' @param data A data frame or matrix with one row per participant. The count +#' columns are read as described under `columns`; with `columns = "names"`, +#' all other columns except `id` are carried to the long data. +#' @param model An `mpt` model from [mpt()], [mpt_from_eqn()] or +#' [mpt_from_string()]. For a model from [mpt_from_eqn()] the `mpt_renaming` +#' attribute maps the labels of the EQN file onto the response categories +#' of the model. For any other model, the columns of `data` are named +#' `.`, e.g. `target.old`, with the tree names and response +#' categories of the model. +#' @param id Character. The name of the participant column in `data`, or +#' `NULL` (the default) to number the rows of `data`, as `id = seq_len(nrow(data))` +#' in a column named `id`. The participants must be unique. +#' @param columns Character. How the count columns of `data` are identified: +#' `"names"` (the default) by their names, as in TreeBUGS data; `"mptinr"` +#' by their position, as in MPTinR data. See Details. +#' +#' @details +#' **Named columns (`columns = "names"`).** For a model from +#' [mpt_from_eqn()], a column is matched to a category by the EQN label +#' (`hit`), or, when trees share labels, by the tree name and label +#' (`target.hit`, with the bmm name of the tree); the qualified name wins +#' where both exist. A label used in several trees must be qualified. +#' Labels that the `categories` argument of [mpt_from_eqn()] merged onto one +#' response category are summed. All columns of `data` that are neither +#' count columns nor `id` are carried unchanged to every long row, e.g. +#' between-participant covariates. +#' +#' **Positional columns (`columns = "mptinr"`).** MPTinR data for an EQN +#' model have no column names. Every column of `data` except `id` is read as +#' a count, in the order of `check.mpt()` in MPTinR (version 1.14.1), +#' elements `eqn.order.trees` and `eqn.order.categories`: trees sorted by +#' their label in the EQN file and, within each tree, the category labels of +#' the file sorted. No other column can be carried; add covariates +#' afterwards, e.g. with [merge()] on the participant column. The help page +#' of `check.mpt()` states only that these elements give the order; the +#' sorting rule is that of its internal `.get.EQN.model.order()`, which +#' sorts the columns that `read.table()` returns for the file: numerically +#' if all tree labels (all category labels) are numbers, else as character +#' strings in the collation of the R session. `mpt_data_reshape()` applies the +#' same rule, so run it in the session where MPTinR's order was determined +#' if the labels mix upper and lower case. This mode needs a model from +#' [mpt_from_eqn()]: MPTinR's own model format orders the columns by the +#' lines of the model text, and [mpt()] reorders the branches of every tree +#' to the categories of the first, so the line order cannot be recovered +#' from the model. For other models, name the columns `.` +#' and use `columns = "names"`. +#' +#' **Impossible categories.** A response category that a tree cannot +#' produce (`impossible` in [mpt_tree()], [mpt_from_string()] and +#' [mpt_from_eqn()]) has no column in `data` and gets 0 in the long rows of +#' that tree, which [bmm()] accepts. +#' +#' **Counts.** The count columns must be numeric, without missing values, +#' and hold non-negative whole numbers; the long data hold them as integers. +#' +#' **Checks.** The function cannot tell a correct reshape from one that +#' swaps columns of equal meaning. Run [bmm_data_check()] on the result, and +#' compare one participant's long rows with the wide row (see the Examples +#' of [mpt_from_eqn()]). +#' +#' @return A data frame with one row per participant and tree, participants +#' first. Its columns are the participant column (`id`), the tree column +#' named by the `tree_id` of the model (absent in a model without it), the +#' carried columns, and one integer count column per response category of +#' the model. +#' +#' @seealso [mpt_from_eqn()] and [mpt_from_string()] for importing the model +#' itself; [mpt()] for the data layout [bmm()] expects; [bmm_data_check()] +#' to check the result before fitting. +#' +#' @keywords transform +#' +#' @examples +#' eqn_file <- tempfile(fileext = ".eqn") +#' writeLines(c( +#' "6", +#' "target hit D_o", +#' "target hit (1-D_o)*g", +#' "target miss (1-D_o)*(1-g)", +#' "lure fa (1-D_n)*g", +#' "lure cr D_n", +#' "lure cr (1-D_n)*(1-g)" +#' ), eqn_file) +#' model <- mpt_from_eqn( +#' eqn_file, +#' restrictions = "D_n = D_o", +#' categories = c(hit = "yes", fa = "yes", miss = "no", cr = "no"), +#' tree_id = "item_type" +#' ) +#' +#' # TreeBUGS data: one named column per category label of the EQN file. +#' # Columns that are no counts, here a between-participant factor, are +#' # carried to the long data +#' wide <- data.frame( +#' subject = c("s1", "s2"), group = c("a", "b"), +#' hit = c(30, 25), miss = c(10, 15), fa = c(8, 12), cr = c(32, 28) +#' ) +#' long <- mpt_data_reshape(wide, model, id = "subject") +#' long +#' bmm_data_check(bmf(Do ~ 1, g ~ 1), long, model) +#' +#' # MPTinR data: unnamed columns, sorted by tree, then by label (cr, fa for +#' # tree lure, then hit, miss for tree target) +#' wide_mptinr <- rbind(c(32, 8, 30, 10), c(28, 12, 25, 15)) +#' identical( +#' mpt_data_reshape(wide_mptinr, model, columns = "mptinr"), +#' mpt_data_reshape(wide[c("cr", "fa", "hit", "miss")], model) +#' ) +#' +#' # trees with response categories that only one tree can produce +#' model_pc <- mpt( +#' list( +#' mpt_tree("pairs", list(E1 = "c * r", E2 = "1 - c * r"), +#' impossible = "F1"), +#' mpt_tree("singles", list(F1 = "u", E2 = "1 - u"), +#' impossible = "E1") +#' ), +#' tree_id = "item_type" +#' ) +#' wide_pc <- data.frame( +#' pairs.E1 = c(20, 22), pairs.E2 = c(20, 18), +#' singles.F1 = c(25, 24), singles.E2 = c(15, 16) +#' ) +#' mpt_data_reshape(wide_pc, model_pc) +#' @export +mpt_data_reshape <- function(data, model, id = NULL, + columns = c("names", "mptinr")) { + stop_missing_args() + columns <- match.arg(columns) + stopif( + !inherits(model, "mpt"), + "The model must be an mpt model, created with mpt(), mpt_from_eqn() or \\ + mpt_from_string()." + ) + stopif( + !is.data.frame(data) && !is.matrix(data), + "The wide data must be a data frame or a matrix." + ) + stopif( + !is.null(id) && !(is.character(id) && length(id) == 1 && !is.na(id)), + "The id argument must be NULL or the name of a single column of the wide \\ + data." + ) + data <- as.data.frame(data, stringsAsFactors = FALSE) + stopif( + anyDuplicated(names(data)) > 0, + "The wide data have duplicated column names: \\ + {collapse_comma(unique(names(data)[duplicated(names(data))]))}. Each \\ + column is read once, so a second column of the same name would be dropped." + ) + stopif(nrow(data) == 0, "The wide data contain no rows.") + stopif( + !is.null(id) && !id %in% names(data), + "The id column '{id}' is not in the wide data. Columns: \\ + {collapse_comma(names(data))}" + ) + + participants <- if (is.null(id)) seq_len(nrow(data)) else data[[id]] + stopif( + anyNA(participants) || anyDuplicated(participants) > 0, + "The participant column '{id %||% 'id'}' must identify each row of the \\ + wide data: it has missing or repeated values." + ) + id_name <- id %||% "id" + + layout <- .mpt_long_layout(model, columns) + non_id <- setdiff(names(data), id) + count_cols <- if (columns == "names") { + .mpt_named_count_columns(layout, non_id) + } else { + .mpt_positional_count_columns(layout, non_id) + } + counts <- data[count_cols] + .mpt_check_wide_counts(counts) + + carried <- data[setdiff(non_id, count_cols)] + tree_id <- model$other_vars$tree_id + resp_cats <- model$resp_vars$resp_cats + taken <- c(id_name, tree_id, resp_cats) + collisions <- intersect(names(carried), taken) + stopif( + length(collisions) > 0, + "The wide data hold column(s) {collapse_comma(collisions)} that the long \\ + data need for the participant, the tree or a response category. Please \\ + rename them. A participant column is passed as id = 'name'." + ) + + tree_names <- names(model$other_vars$trees) + blocks <- lapply(tree_names, function(tree) { + in_tree <- layout$tree == tree + tree_counts <- matrix(0, nrow(data), length(resp_cats), + dimnames = list(NULL, resp_cats)) + summed <- t(rowsum( + t(as.matrix(counts[in_tree])), layout$category[in_tree], reorder = FALSE + )) + tree_counts[, colnames(summed)] <- summed + storage.mode(tree_counts) <- "integer" + data.frame( + c( + setNames(list(participants), id_name), + if (!is.null(tree_id)) setNames(list(rep(tree, nrow(data))), tree_id), + as.list(carried), + as.list(as.data.frame(tree_counts)) + ), + check.names = FALSE + ) + }) + long <- do.call(rbind, blocks) + long <- long[order( + rep(seq_len(nrow(data)), times = length(tree_names)), + rep(seq_along(tree_names), each = nrow(data)) + ), ] + rownames(long) <- NULL + long +} + +# one row per count column of the wide data: the tree, the response category +# of the model, and for EQN models the label of the file. For columns = +# "mptinr" the rows are in the order of MPTinR's columns. +.mpt_long_layout <- function(model, columns) { + renaming <- attr(model, "mpt_renaming") + if (is.null(renaming)) { + stopif( + columns == "mptinr", + "columns = 'mptinr' needs a model from mpt_from_eqn(), because the \\ + column order of MPTinR follows the labels of the EQN file. For other \\ + models, name the columns '.' and use columns = 'names'." + ) + trees <- model$other_vars$trees + return(data.frame( + tree = rep(names(trees), lengths(lapply(trees, `[[`, "branches"))), + category = unlist(lapply(trees, function(tree) names(tree$branches)), + use.names = FALSE), + label = NA_character_ + )) + } + rows <- renaming[renaming$kind == "category", ] + layout <- data.frame( + tree = rows$tree, category = rows$bmm_name, label = rows$file_name + ) + if (columns == "mptinr") { + layout <- layout[.mpt_eqn_column_order(layout, renaming), ] + } + layout +} + +# MPTinR sorts the tree labels and, within a tree, the category labels, each +# after read.table() has typed the file's column: numbers sort numerically, +# anything else as character strings +.mpt_eqn_column_order <- function(layout, renaming) { + trees <- renaming[renaming$kind == "tree", ] + tree_key <- utils::type.convert(trees$file_name, as.is = TRUE) + order( + tree_key[match(layout$tree, trees$bmm_name)], + utils::type.convert(layout$label, as.is = TRUE) + ) +} + +.mpt_named_count_columns <- function(layout, available) { + qualified <- paste(layout$tree, ifelse( + is.na(layout$label), layout$category, layout$label + ), sep = ".") + shared_label <- !is.na(layout$label) & !duplicated(layout$label) & + !duplicated(layout$label, fromLast = TRUE) + count_cols <- ifelse( + qualified %in% available, qualified, + ifelse(shared_label & layout$label %in% available, layout$label, NA) + ) + missing_cols <- is.na(count_cols) + stopif( + any(missing_cols), + "The wide data lack count columns for the following response categories \\ + (expected column names in parentheses): \\ + {paste0(layout$category[missing_cols], ' in tree ', layout$tree[missing_cols], + ' (', qualified[missing_cols], + ifelse(shared_label[missing_cols], + paste0(' or ', layout$label[missing_cols]), ''), ')', + collapse = ', ')}" + ) + # a bare label of one tree can equal the qualified name of a label of another + shared_cols <- unique(count_cols[duplicated(count_cols)]) + stopif( + length(shared_cols) > 0, + "The wide column(s) {collapse_comma(shared_cols)} would be read as the \\ + counts of more than one response category. Name the count columns of \\ + these labels by tree and label, {collapse_comma(qualified[count_cols %in% \\ + shared_cols])}." + ) + count_cols +} + +.mpt_positional_count_columns <- function(layout, available) { + stopif( + length(available) != nrow(layout), + "With columns = 'mptinr', every column of the wide data except the \\ + participant column is read as a count, in the order of MPTinR. The model \\ + expects {nrow(layout)} columns ({collapse_comma(layout$category)} in the \\ + trees {collapse_comma(unique(layout$tree))}), but the wide data have \\ + {length(available)}." + ) + available +} + +.mpt_check_wide_counts <- function(counts) { + non_numeric <- names(counts)[!vapply(counts, is.numeric, logical(1))] + stopif( + length(non_numeric) > 0, + "The count columns must be numeric. Not numeric: \\ + {collapse_comma(non_numeric)}" + ) + values <- as.matrix(counts) + affected <- function(problem) collapse_comma(colnames(values)[colSums(problem) > 0]) + stopif( + anyNA(values), + "The count columns must not contain missing values. Missing values in: \\ + {affected(is.na(values))}" + ) + not_counts <- values < 0 | values != round(values) | is.infinite(values) + stopif( + any(not_counts), + "The count columns must hold non-negative whole numbers. Other values in: \\ + {affected(not_counts)}" + ) + invisible(NULL) +} diff --git a/R/model_mpt_expressions.R b/R/model_mpt_expressions.R index 82b41299e..0fee638bd 100644 --- a/R/model_mpt_expressions.R +++ b/R/model_mpt_expressions.R @@ -22,6 +22,36 @@ unlist(lapply(as.list(expr)[-1], .mpt_scientific_constants)) %||% numeric(0) } +# mpt_tree() formulas, category ~ expression, become the same named list of +# language objects that the string form produces; the right-hand sides are +# never evaluated and the formula environments are dropped +.mpt_formula_branches <- function(formulas, tree_name) { + is_formula <- vapply(formulas, inherits, logical(1), "formula") + stopif( + !all(is_formula), + "Cannot mix formulas and other branch expressions in tree '{tree_name}'. \\ + Write every branch as a formula, category ~ expression, or pass a named \\ + list of character strings." + ) + one_sided <- lengths(formulas) != 3L + stopif( + any(one_sided), + "Branch formulas must be two-sided, category ~ expression (tree \\ + '{tree_name}'). One-sided: {collapse_comma(vapply(formulas[one_sided], deparse1, ''))}" + ) + lhs_is_symbol <- vapply(formulas, function(f) is.symbol(f[[2]]), logical(1)) + stopif( + !all(lhs_is_symbol), + "The left-hand side of a branch formula must be a single response category \\ + name (tree '{tree_name}'). Offending: \\ + {collapse_comma(vapply(formulas[!lhs_is_symbol], deparse1, ''))}" + ) + setNames( + lapply(formulas, function(f) f[[3]]), + vapply(formulas, function(f) as.character(f[[2]]), character(1)) + ) +} + # the constants as the user wrote them in the branch string, each followed by # the form brms writes when the two differ; a value folded from constant # arithmetic (1/10000) has no written form and is shown as brms writes it @@ -670,8 +700,14 @@ too_low <- if (point <= length(interior)) probs <= 0 else probs < -tolerance outside <- names(probs)[too_low | probs > 1 + tolerance][1] if (!is.na(outside)) { + # interpolated into the glue() below, which leaves the line + # continuations of an interpolated value in place zero_hint <- if (probs[[outside]] == 0) { - " A zero probability makes the likelihood undefined." + paste( + " A zero probability makes the likelihood undefined. Declare a", + "response category that a tree cannot produce with", + "mpt_tree(impossible = ) instead of giving it a zero branch." + ) } else { "" } diff --git a/R/model_mpt_import.R b/R/model_mpt_import.R index 8789bc16c..f52c06632 100644 --- a/R/model_mpt_import.R +++ b/R/model_mpt_import.R @@ -14,6 +14,10 @@ #' #' @param text Character. The model definition, either as a single string or #' as a vector of lines (e.g., from [readLines()]). +#' @param restrictions Parameter restrictions as in [mpt()], written with the +#' parameter names of the model text, e.g. `"g = 0.5"` or `list(g = 0.5)`. +#' They mirror the `restrictions` argument of MPTinR and TreeBUGS. +#' Restriction files are not read: pass their lines (see Examples). #' @param tree_names Character vector with one name per tree block. The names #' are matched to the blocks by position: the first name goes to the first #' block, and so on. bmm cannot check this order, so names given in a @@ -36,8 +40,8 @@ #' each observation belongs to (see [mpt()]). Add this column to the data #' yourself: for each row it holds the name, as given in `tree_names`, of #' the tree the counts come from. -#' @param covariates,simplex,restrictions,links passed to [mpt()]. -#' `print(model)` lists the symbols that became parameters. +#' @param covariates,simplex,links passed to [mpt()]. `print(model)` lists +#' the symbols that became parameters. #' @param impossible A named list with the response categories that cannot #' occur in a tree, one element per such tree: #' `list( = c("", ...))` (see the `impossible` argument of @@ -77,9 +81,8 @@ #' vector for aggregated counts) and one unnamed column per line of the #' model file, in file order, because each line is one response category. #' (For EQN files MPTinR sorts the columns instead; see [mpt_from_eqn()].) -#' Name the columns by tree and category first, then stack one block per -#' tree (see Examples); when a tree lacks some categories, start each tree -#' from zeros, as in the Examples of [mpt_from_eqn()]. +#' Name the columns `.`; [mpt_data_reshape()] then stacks the +#' trees and gives 0 to the categories a tree cannot produce (see Examples). #' #' **Provenance.** The returned model carries an `mpt_source` attribute, #' `list(text = )`, with the model definition as one string. @@ -91,7 +94,7 @@ #' [mpt()] and [mpt_tree()] for building the trees in R, including the data #' layout `bmm()` expects. #' -#' @keywords transform +#' @keywords bmmodel #' #' @examples #' model_2htm <- " @@ -104,7 +107,7 @@ #' " #' model <- mpt_from_string( #' model_2htm, -#' tree_names = c("old", "new"), +#' tree_names = c("target", "lure"), #' tree_id = "item_type" #' ) #' model @@ -112,7 +115,7 @@ #' # fix the guessing rate as in MPTinR / TreeBUGS #' mpt_from_string( #' model_2htm, -#' tree_names = c("old", "new"), +#' tree_names = c("target", "lure"), #' tree_id = "item_type", #' restrictions = "g = 0.5" #' ) @@ -149,21 +152,17 @@ #' restr <- trimws(sub("#.*", "", readLines(restr_file))) #' model <- mpt_from_string( #' readLines(model_file), -#' tree_names = c("old", "new"), +#' tree_names = c("target", "lure"), #' tree_id = "item_type", #' restrictions = restr[nzchar(restr)] #' ) #' #' # MPTinR data: one unnamed column per line of the model file, in file -#' # order. Name the columns by tree and category, then stack one block per -#' # tree to get bmm's one row per participant and tree. +#' # order. Name the columns by tree and category, then let mpt_data_reshape() +#' # stack the trees to get bmm's one row per participant and tree. #' wide <- rbind(c(30, 10, 8, 32), c(25, 15, 12, 28)) -#' colnames(wide) <- c("old.old", "old.new", "new.old", "new.new") -#' long <- do.call(rbind, lapply(c("old", "new"), function(tree) { -#' counts <- wide[, paste(tree, c("old", "new"), sep = "."), drop = FALSE] -#' colnames(counts) <- c("old", "new") -#' data.frame(id = seq_len(nrow(wide)), item_type = tree, counts) -#' })) +#' colnames(wide) <- c("target.old", "target.new", "lure.old", "lure.new") +#' long <- mpt_data_reshape(wide, model) #' long #' # the row sums are the trials per row (40 here), but they cannot show #' # swapped columns or tree labels: compare one participant with the wide @@ -174,9 +173,9 @@ #' wide[1, ] #' bmm_data_check(bmf(D ~ 1), long, model) #' @export -mpt_from_string <- function(text, tree_names, categories = NULL, - tree_id = NULL, covariates = NULL, - simplex = NULL, restrictions = NULL, +mpt_from_string <- function(text, restrictions = NULL, tree_names, + categories = NULL, tree_id = NULL, + covariates = NULL, simplex = NULL, links = "logit", impossible = NULL) { stop_missing_args() # a vector of lines keeps its "" elements, which separate the trees @@ -376,7 +375,7 @@ mpt_from_string <- function(text, tree_names, categories = NULL, #' `impossible` does not cover the difference, the error names the #' categories each tree lacks. #' @param tree_names A named character vector mapping tree labels of the EQN -#' file onto tree names in bmm, e.g. `c("1" = "old", "2" = "new")`. bmm tree +#' file onto tree names in bmm, e.g. `c("1" = "target", "2" = "lure")`. bmm tree #' names must start with a letter and contain only letters, digits, or #' underscores, so files with numeric tree labels need this mapping. Labels #' not listed keep their name. @@ -412,9 +411,11 @@ mpt_from_string <- function(text, tree_names, categories = NULL, #' parameter names, parameter and response category names are sanitized: a #' leading `_` prefix is removed from category names that are not #' mapped with `categories`, and all remaining underscores and dots are -#' stripped. Names that would clash after this are an error; names are -#' case-sensitive, so `P_G1` and `P_g1` stay distinct (`PG1`, `Pg1`). Tree -#' labels are kept, unless `tree_names` maps them. The formulas passed to +#' stripped. Two labels of the file, or two different values within one +#' tree, that become the same name are an error; mapped values that become +#' the same name in different trees are that one shared category. Names +#' are case-sensitive, so `P_G1` and `P_g1` stay distinct (`PG1`, `Pg1`). +#' Tree labels are kept, unless `tree_names` maps them. The formulas passed to #' [bmm()] use the bmm names (`bmf(Do ~ 1)` for a file parameter `D_o`; #' `bmf()` is short for [bmmformula()]); the `parameter` rows of the #' renaming map translate them, and `print(model)` lists them. Of the @@ -434,33 +435,24 @@ mpt_from_string <- function(text, tree_names, categories = NULL, #' one count column per shared response category, named by `bmm_name` in #' the renaming map (see [mpt()]). A category that is `impossible` in a #' tree gets 0 (or `NA`) in the rows of that tree. TreeBUGS data are wide: -#' one row per participant, one column per category label of the EQN file, -#' and no participant column, so add an `id` column first. The `category` -#' rows of the renaming map hold every tree with its labels and bmm names, -#' so the wide data can be reshaped by name: for each tree, start from a -#' block of zeros over all shared categories and fill in the columns of -#' that tree (see Examples). MPTinR data for an EQN model are unnamed, and -#' their columns do not follow the file: the trees and, within each tree, -#' the category labels are sorted (alphabetically, or numerically when they -#' are numbers); `MPTinR::check.mpt("model.eqn")` returns the order in its -#' elements `eqn.order.trees` and `eqn.order.categories`. (For MPTinR's own -#' model format the columns follow the lines instead; see -#' [mpt_from_string()].) Name the columns by tree and label, -#' `paste(tree, label, sep = ".")` with the bmm tree names, because trees -#' may share labels, and select them with the same key (see Examples). -#' -#' If `categories` merged several labels of a tree onto one name, sum their -#' columns in the wide data into the first label's column and drop the map -#' rows of the other labels (`map[!duplicated(map[c("tree", "bmm_name")]), ]`) -#' before reshaping; otherwise the recipe fails with "duplicate subscripts". +#' one row per participant, one column per category label of the EQN file. +#' MPTinR data for an EQN model are unnamed, and their columns do not follow +#' the file: the trees and, within each tree, the category labels are sorted +#' (alphabetically, or numerically when all labels in the file are numbers); +#' `MPTinR::check.mpt("model.eqn")` returns the order in its elements +#' `eqn.order.trees` and `eqn.order.categories`. (For MPTinR's own model +#' format the columns follow the lines instead; see [mpt_from_string()].) +#' [mpt_data_reshape()] reshapes both layouts with the renaming map, sums the +#' labels that `categories` merged onto one response category, and gives 0 to +#' the categories a tree cannot produce (see Examples). #' #' bmm cannot tell a reshape that swaps tree labels or category columns -#' from a correct one: such data fit without an error. The row sums, -#' `rowSums(long[cats])`, are the trials per row; they cannot show swapped -#' columns, and they show swapped tree labels only when the trees differ -#' in their number of trials. Comparing the response proportions per tree -#' with the design (old items must get more "yes" responses than new -#' items) catches swapped tree labels and response columns swapped in +#' from a correct one: such data fit without an error. The row sums of the +#' count columns are the trials per row; they cannot show swapped columns, +#' and they show swapped tree labels only when the trees differ in their +#' number of trials. Comparing the response proportions per tree +#' with the design (targets must get more "yes" responses than lures) +#' catches swapped tree labels and response columns swapped in #' every tree, but not always a column slipped or a map entry wrong in one #' tree, nor a swap between two trees of the same kind; so also compare #' one participant's long rows with the wide row. [bmm_data_check()] runs @@ -479,18 +471,18 @@ mpt_from_string <- function(text, tree_names, categories = NULL, #' [mpt()] and [mpt_tree()] for building the trees in R, including the data #' layout `bmm()` expects. #' -#' @keywords transform +#' @keywords bmmodel #' #' @examples #' eqn_file <- tempfile(fileext = ".eqn") #' writeLines(c( #' "6", -#' "old hit D_o", -#' "old hit (1-D_o)*g", -#' "old miss (1-D_o)*(1-g)", -#' "new fa (1-D_n)*g", -#' "new cr D_n", -#' "new cr (1-D_n)*(1-g)" +#' "target hit D_o", +#' "target hit (1-D_o)*g", +#' "target miss (1-D_o)*(1-g)", +#' "lure fa (1-D_n)*g", +#' "lure cr D_n", +#' "lure cr (1-D_n)*(1-g)" #' ), eqn_file) #' # restrictions use the file's names; the formulas use the bmm names #' model <- mpt_from_eqn( @@ -502,50 +494,29 @@ mpt_from_string <- function(text, tree_names, categories = NULL, #' attr(model, "mpt_renaming") #' #' # TreeBUGS data: one row per participant, columns named by the category -#' # labels of the EQN file, no participant column. bmm needs one row per -#' # participant and tree, columns named by bmm_name: reshape by name through -#' # the map, starting each tree from zeros so that categories the tree -#' # cannot produce get 0 +#' # labels of the EQN file, no participant column. mpt_data_reshape() reshapes +#' # them to bmm's one row per participant and tree, through the map; the +#' # categories a tree cannot produce get 0 #' wide <- data.frame( #' hit = c(30, 25), miss = c(10, 15), fa = c(8, 12), cr = c(32, 28) #' ) -#' wide$id <- seq_len(nrow(wide)) -#' map <- attr(model, "mpt_renaming") -#' map <- map[map$kind == "category", ] -#' cats <- unique(map$bmm_name) -#' long <- do.call(rbind, lapply(split(map, map$tree), function(tree_map) { -#' counts <- as.data.frame(matrix(0, nrow(wide), length(cats))) -#' names(counts) <- cats -#' counts[tree_map$bmm_name] <- wide[tree_map$file_name] -#' data.frame(id = wide$id, item_type = tree_map$tree[1], counts) -#' })) +#' long <- mpt_data_reshape(wide, model) #' long #' # the row sums are the trials per row (40 here). The proportions catch -#' # swapped tree labels and swapped response columns (old items must get -#' # more "yes" than new items), but not always a column slipped in one +#' # swapped tree labels and swapped response columns (targets must get +#' # more "yes" than lures), but not always a column slipped in one #' # tree: also compare one participant's long rows with the wide row -#' rowSums(long[cats]) -#' tapply(long$yes / rowSums(long[cats]), long$item_type, mean) +#' rowSums(long[c("yes", "no")]) +#' tapply(long$yes / rowSums(long[c("yes", "no")]), long$item_type, mean) #' long[long$id == 1, ] #' wide[1, ] #' bmm_data_check(bmf(Do ~ 1, g ~ 1), long, model) #' #' # MPTinR data for an EQN model are unnamed and sorted by tree, then by #' # label: MPTinR::check.mpt(eqn_file)$eqn.order.categories lists the -#' # labels (here cr, fa for tree new, then hit, miss for tree old). Name the -#' # columns by tree and label, so that a label used in several trees stays -#' # apart, and select with the same key +#' # labels (here cr, fa for tree new, then hit, miss for tree old) #' wide_mptinr <- rbind(c(32, 8, 30, 10), c(28, 12, 25, 15)) -#' colnames(wide_mptinr) <- c("new.cr", "new.fa", "old.hit", "old.miss") -#' long2 <- do.call(rbind, lapply(split(map, map$tree), function(tree_map) { -#' counts <- as.data.frame(matrix(0, nrow(wide_mptinr), length(cats))) -#' names(counts) <- cats -#' key <- paste(tree_map$tree, tree_map$file_name, sep = ".") -#' counts[tree_map$bmm_name] <- wide_mptinr[, key] -#' data.frame(id = seq_len(nrow(wide_mptinr)), item_type = tree_map$tree[1], -#' counts) -#' })) -#' identical(long2, long) +#' identical(mpt_data_reshape(wide_mptinr, model, columns = "mptinr"), long) #' @export mpt_from_eqn <- function(file, restrictions = NULL, categories = NULL, tree_id = NULL, covariates = NULL, simplex = NULL, @@ -843,10 +814,23 @@ mpt_from_eqn <- function(file, restrictions = NULL, categories = NULL, removed. Rename them in the EQN file or map them with the categories \\ argument." ) - cat_renaming <- .mpt_sanitized_names(unique(category_clean)) + # a mapped value is the user's choice of a shared name, so two of them that + # sanitize alike in different trees are one category; two file labels that + # do so are a clash in any tree, as are two different values in one tree + labels <- unique(data.frame( + tree = eqn$tree, original = eqn$category, clean = category_clean, + mapped = user_mapped, bmm = gsub("[._]", "", category_clean) + )) + in_tree <- unique(labels[c("tree", "clean", "bmm")]) + file_labels <- unique(labels[!labels$mapped, c("clean", "bmm")]) + clashing <- unique(c( + in_tree$clean[.mpt_repeated(in_tree[c("tree", "bmm")])], + file_labels$clean[.mpt_repeated(file_labels["bmm"])] + )) # the message names the labels as the EQN file writes them, before the # prefix strip - names(cat_renaming) <- eqn$category[match(names(cat_renaming), category_clean)] + cat_renaming <- .mpt_sanitized_names(clashing) + names(cat_renaming) <- labels$original[match(clashing, labels$clean)] .mpt_check_sanitized_names(cat_renaming, for_categories = TRUE) eqn$bmm_category <- gsub("[._]", "", category_clean) eqn$prefix <- ifelse(has_tree_prefix, paste0(eqn$tree, "_"), "") @@ -972,7 +956,7 @@ mpt_from_eqn <- function(file, restrictions = NULL, categories = NULL, .mpt_model_from_branches <- function(branches, impossible, ...) { trees <- Map(function(tree_name, tree_branches) { - mpt_tree(tree_name, tree_branches, impossible[[tree_name]]) + mpt_tree(tree_name, tree_branches, impossible = impossible[[tree_name]]) }, names(branches), branches) mpt(trees = trees, ...) } @@ -1009,8 +993,11 @@ mpt_from_eqn <- function(file, restrictions = NULL, categories = NULL, # the rows whose bmm category is shared by another row of the same tree .mpt_merged_rows <- function(rows) { - merge_key <- rows[c("tree", "bmm_category")] - rows[duplicated(merge_key) | duplicated(merge_key, fromLast = TRUE), ] + rows[.mpt_repeated(rows[c("tree", "bmm_category")]), ] +} + +.mpt_repeated <- function(key) { + duplicated(key) | duplicated(key, fromLast = TRUE) } .mpt_format_merges <- function(merged) { diff --git a/R/model_mpt_plot.R b/R/model_mpt_plot.R index e8217f977..f9800e326 100644 --- a/R/model_mpt_plot.R +++ b/R/model_mpt_plot.R @@ -69,19 +69,19 @@ #' empirical review of multinomial process tree modeling. Psychonomic #' Bulletin & Review, 6(1), 57-86. https://doi.org/10.3758/BF03210812 #' -#' @keywords transform +#' @keywords bmmodel #' #' @examples -#' tree_old <- mpt_tree("old", list( +#' tree_target <- mpt_tree("target", list( #' old = "D + (1 - D) * g", #' new = "(1 - D) * (1 - g)" #' )) -#' plot(tree_old) +#' plot(tree_target) #' #' model <- mpt( #' list( -#' tree_old, -#' mpt_tree("new", list( +#' tree_target, +#' mpt_tree("lure", list( #' old = "(1 - D) * g", #' new = "D + (1 - D) * (1 - g)" #' )) diff --git a/R/pp_check.R b/R/pp_check.R index 402b6b734..4220338f6 100644 --- a/R/pp_check.R +++ b/R/pp_check.R @@ -57,7 +57,9 @@ #' `group = NA` to pool them into a single profile. The [sdt_ranking()] model #' shows the rank positions as the categories of the proportion profile; for #' mixed set sizes, pass the set-size column as `group` so each facet has a -#' homogeneous denominator. The custom-family +#' homogeneous denominator. For [mpt()] fits, the `tree_id` column can be +#' named even if no parameter formula uses it, giving one facet per tree; +#' the default pools the trees. The custom-family #' models [sdt_yn()] and [sdt_mafc()] delegate to [brms::pp_check()]; for #' their aggregated counts an overlaid empirical CDF (`type = "ecdf_overlay"`, #' adding `group = "stimulus"` for [sdt_yn()]) is usually clearer than the @@ -144,6 +146,7 @@ pp_check.bmmfit <- function(object, type = NULL, ndraws = NULL, if (identical(family(object)$family, "multinomial")) { group <- .pp_check_resolve_group(object, group) object <- .pp_check_restore_set_size(object, group) + object <- .pp_check_restore_tree(object, group) return(.pp_check_multinomial(object, type = type, ndraws = ndraws %||% 100L, group = group, ...)) } @@ -203,6 +206,22 @@ pp_check.bmmfit <- function(object, type = NULL, ndraws = NULL, } +# brms keeps the tree column of an mpt fit only when a formula names it; +# otherwise the one-hot tree indicators are the record of each row's tree +.pp_check_restore_tree <- function(object, group) { + model <- object$bmm$model + indicators <- model$other_vars$indicators$tree + tree_id <- model$other_vars$tree_id + if (!inherits(model, "mpt") || is.null(indicators) || + !identical(group, tree_id) || tree_id %in% names(object$data)) { + return(object) + } + tree <- max.col(as.matrix(object$data[indicators]), ties.method = "first") + object$data[[tree_id]] <- factor(names(indicators)[tree], levels = names(indicators)) + object +} + + .auto_grouped_type <- function(type) { grouped <- paste0(type, "_grouped") if (endsWith(type, "_grouped") || is.null(.ppc_fun(grouped))) type else grouped diff --git a/R/ppp.R b/R/ppp.R new file mode 100644 index 000000000..76af6a816 --- /dev/null +++ b/R/ppp.R @@ -0,0 +1,471 @@ +############################################################################# ! +# PPP.R #### +# Posterior predictive p-values (T1, T2) for the count models. #### +############################################################################# ! + +#' Posterior predictive p-values for count models +#' +#' Computes the goodness-of-fit statistics \eqn{T_1} and \eqn{T_2} of Klauer +#' (2010) and their posterior predictive p-values (PPP) for models that fit +#' aggregated response counts: [mpt()], [m3()], [sdt_yn()], [sdt_rating()], +#' [sdt_mafc()], [sdt_ranking()] and [sdt_cdp()]. TreeBUGS' `PPP()` reports the +#' same two statistics; Details says where their values differ. +#' +#' @details +#' Both statistics are computed once per posterior draw, for the observed counts +#' and for counts simulated from that draw with [brms::posterior_predict()], +#' each against the expected counts of the same draw from +#' [brms::posterior_epred()]. The simulated data are new data from the same +#' participants. The PPP is the proportion of draws in which the simulated +#' counts deviate more from the expectation than the observed counts do. A PPP +#' near 0 flags misfit; a PPP away from 0 means that the model reproduces the +#' aspect of the data the statistic measures. +#' +#' **\eqn{T_1}: mean structure.** The rows of the data are summed into cells, +#' and \eqn{T_1} is Pearson's \eqn{X^2} over all cells and response categories, +#' \deqn{T_1 = \sum_{c,k} \frac{(n_{ck} - \hat{n}_{ck})^2}{\hat{n}_{ck}},} +#' where \eqn{n_{ck}} is the observed and \eqn{\hat{n}_{ck}} the expected count +#' of category \eqn{k} in cell \eqn{c}. Klauer (2010, Eq. 17) and TreeBUGS +#' average the counts over participants instead of summing them. With the same +#' number \eqn{T} of participants in every cell, the two give the same PPP and +#' \eqn{T_1} here is \eqn{T} times theirs; otherwise the mean gives the cells +#' with fewer participants more weight. +#' +#' **\eqn{T_2}: covariance structure.** For every participant, the counts are +#' summed per cell, giving one vector per participant over all cells and +#' categories. \eqn{T_2} compares the covariance matrix \eqn{s} of these vectors +#' across participants with the covariance matrix \eqn{\sigma} the model +#' expects (Klauer, 2010, Eq. 18): +#' \deqn{T_2 = \sum_{i,j} \frac{(s_{ij} - \sigma_{ij})^2}{\sqrt{\sigma_{ii}\sigma_{jj}}}.} +#' \eqn{\sigma} is the covariance of the expected counts across participants +#' plus the average multinomial covariance within a participant, which is the +#' expected value of \eqn{s} given the draw. Both covariances divide by +#' \eqn{T - 1}, where \eqn{T} is the number of participants. Klauer (2010) +#' multiplies the multinomial part by \eqn{(T - 1) / T}, the expected value +#' for covariances that divide by \eqn{T}; that version of \eqn{T_2} is +#' \eqn{(T - 1) / T} times this one and has the same PPP. TreeBUGS combines +#' the factor with covariances that divide by \eqn{T - 1}, so its \eqn{\sigma} +#' falls short of the expected \eqn{s} by \eqn{1 / T} of the multinomial part. +#' The sum runs over +#' the full matrix, so the covariances between trees or conditions, which +#' carry the correlations between the person parameters, are part of the +#' check. When participants differ in the cells they have, as in a +#' between-participants design, they are grouped by their set of cells, +#' \eqn{T_2} is computed per group, and the results are summed. A participant +#' who shares the set of cells with nobody else has no covariance and is left +#' out, with a message. +#' +#' **Cells.** By default, a cell is a combination of the tree ([mpt()]), the +#' stimulus (the SDT models) or the set size ([sdt_ranking()]), and the +#' categorical population-level predictors of the model formula. Numeric +#' predictors are not part of the default, because every distinct value would +#' form a cell of its own. Rows of a cell are summed, so the trees or stimuli +#' must be part of `by` for the statistic to compare each tree's own +#' categories. +#' +#' **Impossible categories.** Categories with an expected count of about 0 in +#' every draw, such as the impossible categories of an [mpt()] or [m3()] model +#' or the ranks beyond the set size of an [sdt_ranking()] model, are left out +#' of both statistics. A cell with a single possible category, such as an +#' [m3()] set size without other candidates, is left out of \eqn{T_2}: the +#' number of trials fixes its count. +#' +#' @param fit A `bmmfit` object returned by [bmm()]. +#' @param statistic Character. `"T1"`, `"T2"`, or both (the default). If the +#' participants cannot be identified, the default computes `"T1"` only. +#' @param by Character vector naming the columns of the model data that define +#' the cells. `NULL` (default) uses the cells described in Details; +#' `character(0)` pools all rows into one cell. For [mpt()] fits, the +#' tree column can be named even if it is not a predictor. +#' @param ndraws Integer. Number of posterior draws to use. `NULL` (default) +#' uses all draws. +#' @param level Character. `"group"` (default) for the statistics of the whole +#' data set; `"participant"` for one \eqn{T_1} and PPP per participant, as +#' TreeBUGS reports them, to find the participants the model fits poorly. +#' @param participant Character. The column that identifies participants, +#' needed for \eqn{T_2} and for `level = "participant"`. It must be a +#' variable of the model formula, because the model data contain only +#' those; a model without random effects has one only if a predictor such +#' as `0 + id` names the participants. `NULL` (default) +#' uses the grouping variable of the random effects; when the model has more +#' than one (e.g. participants and items), name it here. +#' @param ... Unused. +#' +#' @return A data frame of class `bmm_ppp` with one row per statistic (and per +#' participant for `level = "participant"`) and the columns `statistic`, +#' `observed` and `replicated` (the means of the statistic over draws, for the +#' observed and the simulated counts, as reported by Klauer, 2010, and +#' TreeBUGS), `ppp` and `ndraws`. The draws of the statistics are stored in +#' the attribute `"draws"`. +#' +#' @references Klauer, K. C. (2010). Hierarchical multinomial processing tree +#' models: A latent-trait approach. *Psychometrika, 75*(1), 70–98. +#' https://doi.org/10.1007/s11336-009-9141-0 +#' +#' @seealso [pp_check.bmmfit()] for the graphical check. +#' @export +#' @examplesIf isTRUE(Sys.getenv("BMM_EXAMPLES")) +#' trees <- list( +#' mpt_tree("target", list(old = "D + (1 - D) * g", new = "(1 - D) * (1 - g)")), +#' mpt_tree("lure", list(old = "(1 - D) * g", new = "D + (1 - D) * (1 - g)")) +#' ) +#' dat <- expand.grid(id = factor(1:20), item_type = c("target", "lure"), +#' stringsAsFactors = FALSE) +#' dat$old <- rbinom(nrow(dat), 50, ifelse(dat$item_type == "target", 0.8, 0.2)) +#' dat$new <- 50 - dat$old +#' +#' fit <- bmm( +#' bmf(D ~ 1 + (1 | id), g ~ 1 + (1 | id)), dat, mpt(trees, "item_type"), +#' cores = 4, backend = "cmdstanr" +#' ) +#' ppp(fit) +#' ppp(fit, level = "participant") +ppp <- function(fit, ...) { + UseMethod("ppp") +} + +#' @rdname ppp +#' @export +ppp.bmmfit <- function(fit, statistic = c("T1", "T2"), by = NULL, ndraws = NULL, + level = c("group", "participant"), participant = NULL, + ...) { + model <- fit$bmm$model + stopif(!inherits(model, .ppp_models), + "ppp() is available for the count models ({collapse_comma(.ppp_models)}), \\ + not for the {utils::tail(class(model), 1)} model.") + statistic_given <- !missing(statistic) + statistic <- match.arg(statistic, several.ok = TRUE) + level <- match.arg(level) + stopif(level == "participant" && statistic_given && "T2" %in% statistic, + "T2 is a group-level statistic; level = 'participant' computes T1 only.") + stopif(!is.null(participant) && + (!is.character(participant) || length(participant) != 1L || + !participant %in% names(fit$data)), + "'participant' must name a column of the model data, which contain \\ + only the variables of the formula.") + stopif(!is.null(ndraws) && + (!is.numeric(ndraws) || length(ndraws) != 1L || ndraws %% 1 != 0 || + ndraws < 1 || ndraws > brms::ndraws(fit)), + "'ndraws' must be a whole number between 1 and {brms::ndraws(fit)}.") + + # one set of draws for posterior_predict() and posterior_epred(), so that + # each replicated data set is compared with its own expectation + draw_ids <- if (!is.null(ndraws)) sample.int(brms::ndraws(fit), ndraws) + counts <- .ppp_counts(model, fit, draw_ids) + by <- by %||% .ppp_default_by(model, fit, counts$data) + stopif(!is.character(by) || !all(by %in% names(counts$data)), + "'by' must name columns of the model data. Available: \\ + {collapse_comma(setdiff(names(counts$data), 'Y'))}.") + + needs_participant <- level == "participant" || "T2" %in% statistic + participant <- participant %||% if (needs_participant) .ppp_participant(fit) + if (needs_participant && length(participant) == 0L) { + requested <- if (level == "participant") "level = 'participant'" else "T2" + stopif(level == "participant" || statistic_given, + "{requested} needs to know which rows belong to which participant, \\ + and the model has no random effects to tell. Name the column with \\ + 'participant'; the model data contain only the variables of the \\ + formula, so it must be one of them.") + message2("The model has no random effects to identify participants, so \\ + ppp() computes T1 only. T2 needs the participant column named \\ + with 'participant', and the model data contain only the \\ + variables of the formula.") + statistic <- "T1" + } + person <- if (length(participant) == 1L) as.factor(counts$data[[participant]]) + + cell <- .ppp_index(counts$data[by]) + out <- if (level == "participant") { + .ppp_by_participant(counts, cell, person) + } else { + .ppp_by_group(counts, cell, statistic, person) + } + structure(out$table, class = c("bmm_ppp", "data.frame"), draws = out$draws, + level = level, by = by, n_cells = max(cell)) +} + +#' @export +print.bmm_ppp <- function(x, digits = 3, ...) { + by <- attr(x, "by") + cells <- if (length(by) == 0L) { + "all rows pooled into one cell" + } else { + glue("{collapse_comma(by)} ({attr(x, 'n_cells')} cells)") + } + cat("Posterior predictive p-values (Klauer, 2010)\n") + cat("Cells: ", cells, "\n\n", sep = "") + print(as.data.frame(unclass(x)), digits = digits, row.names = FALSE) + invisible(x) +} + +.ppp_models <- c("mpt", "m3", "sdt_yn", "sdt_rating", "sdt_mafc", "sdt_ranking", + "sdt_cdp") + + +############################################################################# ! +# COUNTS PER MODEL #### +############################################################################# ! + +# observed: rows x categories; replicated and expected: draws x rows x +# categories, on the count scale; data: the model data the cells come from +.ppp_counts <- function(model, fit, draw_ids) { + UseMethod(".ppp_counts") +} + +.ppp_counts.default <- function(model, fit, draw_ids) { + list( + observed = unclass(fit$data$Y), + replicated = brms::posterior_predict(fit, draw_ids = draw_ids), + expected = brms::posterior_epred(fit, draw_ids = draw_ids), + data = fit$data + ) +} + +# the response column holds one category; the other is trials - Y +.ppp_counts.sdt_yn <- function(model, fit, draw_ids) { + response <- model$resp_vars$response + trials <- fit$data[[model$other_vars$n_trials]] + y <- fit$data[[response]] + list( + observed = matrix(c(y, trials - y), ncol = 2L, + dimnames = list(NULL, c(response, "other"))), + replicated = .ppp_complement(brms::posterior_predict(fit, draw_ids = draw_ids), + trials), + expected = .ppp_complement(brms::posterior_epred(fit, draw_ids = draw_ids), + trials), + data = fit$data + ) +} + +# draws x rows counts of one category -> draws x rows x 2 with trials - count +.ppp_complement <- function(draws, trials) { + array(c(draws, rep(trials, each = nrow(draws)) - draws), c(dim(draws), 2L)) +} + +.ppp_counts.sdt_mafc <- .ppp_counts.sdt_yn + +# the tree column enters no formula, so brms leaves it out of the model data; +# the tree indicators that do enter it identify the tree of each row +.ppp_counts.mpt <- function(model, fit, draw_ids) { + out <- NextMethod() + indicators <- model$other_vars$indicators$tree + tree_id <- model$other_vars$tree_id + if (is.null(indicators) || tree_id %in% names(out$data)) { + return(out) + } + tree <- max.col(as.matrix(out$data[indicators]), ties.method = "first") + out$data[[tree_id]] <- factor(names(indicators)[tree], levels = names(indicators)) + out +} + + +############################################################################# ! +# CELLS #### +############################################################################# ! + +.ppp_default_by <- function(model, fit, data) { + UseMethod(".ppp_default_by") +} + +.ppp_default_by.default <- function(model, fit, data) { + predictors <- intersect(unique(unlist(.population_vars(fit$bmm$user_formula))), + names(data)) + is_categorical <- vapply(data[predictors], function(x) { + is.factor(x) || is.character(x) || is.logical(x) + }, logical(1)) + predictors[is_categorical] +} + +.ppp_default_by.mpt <- function(model, fit, data) { + tree <- if (!is.null(model$other_vars$indicators$tree)) model$other_vars$tree_id + unique(c(tree, NextMethod())) +} + +.ppp_default_by.sdt <- function(model, fit, data) { + unique(c(model$other_vars$stimulus, NextMethod())) +} + +.ppp_default_by.sdt_ranking <- function(model, fit, data) { + unique(c("max_rank", NextMethod())) +} + +# Only plain grouping variables can identify a participant: an interaction +# such as id:session groups the rows of a participant further +.ppp_participant <- function(fit) { + candidates <- grep(":", re_group_vars(fit$bmm$user_formula), fixed = TRUE, + invert = TRUE, value = TRUE) + candidates <- intersect(candidates, names(fit$data)) + stopif(length(candidates) > 1L, + "The model has several grouping variables ({collapse_comma(candidates)}). \\ + Name the one that identifies participants with 'participant'.") + candidates +} + +.ppp_index <- function(df) { + if (ncol(df) == 0L) { + return(rep(1L, nrow(df))) + } + as.integer(interaction(df, drop = TRUE, lex.order = TRUE)) +} + +# sums the rows of a draws x rows x categories array within each group +.ppp_rowsum <- function(x, group) { + d <- dim(x) + summed <- rowsum(matrix(aperm(x, c(2L, 1L, 3L)), d[2]), group) + aperm(array(summed, c(nrow(summed), d[1], d[3])), c(2L, 1L, 3L)) +} + +# A category is possible in a cell if its expected share is not negligible in +# at least one draw. Impossible categories sit at exp(-100) relative to the +# others, far below any probability a fitted category reaches. +.ppp_possible <- function(observed, expected) { + share <- expected / as.vector(apply(expected, c(1L, 2L), sum)) + possible <- apply(!is.na(share) & share > 1e-10, c(2L, 3L), any) + impossible_seen <- observed > 0 & !possible + stopif(any(impossible_seen), + "The model gives probability 0 to observed responses in categories \\ + {collapse_comma(unique(colnames(observed)[col(observed)[impossible_seen]]))}.") + possible +} + + +############################################################################# ! +# STATISTICS #### +############################################################################# ! + +.ppp_by_group <- function(counts, cell, statistic, person) { + draws <- list() + if ("T1" %in% statistic) { + draws$T1 <- .ppp_t1(counts, cell) + } + if ("T2" %in% statistic) { + draws$T2 <- .ppp_t2(counts, cell, person) + } + table <- do.call(rbind, lapply(names(draws), function(stat) { + .ppp_summary(draws[[stat]][, "observed", drop = FALSE], + draws[[stat]][, "replicated", drop = FALSE], stat) + })) + list(table = table, draws = draws) +} + +# one T1 per participant over all of the participant's cells, as TreeBUGS' +# individual T1; the draws are draws x participants matrices +.ppp_by_participant <- function(counts, cell, person) { + person_cell <- .ppp_index(data.frame(person, cell)) + terms <- .ppp_t1_terms(counts, person_cell) + pc_person <- person[match(seq_len(max(person_cell)), person_cell)] + per_person <- function(x) t(rowsum(t(apply(x, c(1L, 2L), sum)), pc_person)) + draws <- list(T1 = list(observed = per_person(terms$observed), + replicated = per_person(terms$replicated))) + table <- data.frame( + participant = factor(colnames(draws$T1$observed), levels(person)), + .ppp_summary(draws$T1$observed, draws$T1$replicated, "T1") + ) + list(table = table, draws = draws) +} + +# observed and replicated are draws x statistics matrices +.ppp_summary <- function(observed, replicated, statistic) { + data.frame( + statistic = statistic, + observed = colMeans(observed), + replicated = colMeans(replicated), + ppp = colMeans(replicated > observed), + ndraws = nrow(observed), + row.names = NULL + ) +} + +.ppp_t1 <- function(counts, cell) { + terms <- .ppp_t1_terms(counts, cell) + ndraws <- dim(terms$observed)[1] + cbind(observed = rowSums(matrix(terms$observed, ndraws)), + replicated = rowSums(matrix(terms$replicated, ndraws))) +} + +# the X^2 terms per draw x group x category, 0 for impossible categories +.ppp_t1_terms <- function(counts, group) { + observed <- rowsum(counts$observed, group) + expected <- .ppp_rowsum(counts$expected, group) + replicated <- .ppp_rowsum(counts$replicated, group) + ndraws <- dim(expected)[1] + impossible <- !rep(.ppp_possible(observed, expected), each = ndraws) + terms <- function(n) { + x <- (n - expected)^2 / expected + x[impossible] <- 0 + x + } + list(observed = terms(rep(observed, each = ndraws)), + replicated = terms(replicated)) +} + +.ppp_t2 <- function(counts, cell, person) { + possible <- .ppp_possible(rowsum(counts$observed, cell), + .ppp_rowsum(counts$expected, cell)) + cell_sets <- tapply(cell, person, function(x) paste(sort(unique(x)), collapse = "-")) + patterns <- split(names(cell_sets), cell_sets) + usable <- lengths(patterns) >= 2L + if (any(!usable)) { + message2("T2 leaves out {sum(lengths(patterns[!usable]))} participant(s) \\ + whose set of cells no other participant shares.") + } + if (!any(usable)) { + warning2("No two participants share the same set of cells, so T2 cannot \\ + be computed.") + missing_draws <- rep(NA_real_, dim(counts$expected)[1]) + return(cbind(observed = missing_draws, replicated = missing_draws)) + } + Reduce(`+`, lapply(patterns[usable], function(members) { + rows <- which(person %in% members) + .ppp_t2_pattern(counts, rows, cell[rows], droplevels(person[rows]), possible) + })) +} + +# T2 for participants who have the same cells: the participants are the rows +# and the possible cell x category pairs the columns of the count matrices +.ppp_t2_pattern <- function(counts, rows, cell, person, possible) { + n_people <- nlevels(person) + cells <- sort(unique(cell)) + # the trials fix the count of a cell's only possible category, so its + # variance is 0 and would sit in the denominator of T2 + possible <- possible[cells, , drop = FALSE] + columns <- which(possible & rowSums(possible) > 1L, arr.ind = TRUE) + col_cell <- columns[, "row"] + col_cat <- columns[, "col"] + # every participant has every cell, so participant p's row sum for the j-th + # cell is group (p - 1) * length(cells) + j + person_cell <- (as.integer(person) - 1L) * length(cells) + match(cell, cells) + lookup <- cbind(c(outer(seq_len(n_people) - 1L, col_cell, + function(p, j) p * length(cells) + j)), + rep(col_cat, each = n_people)) + wide <- function(per_group) matrix(per_group[lookup], n_people) + + observed <- counts$observed[rows, , drop = FALSE] + observed_cov <- stats::cov(wide(rowsum(observed, person_cell))) + expected <- .ppp_rowsum(counts$expected[, rows, , drop = FALSE], person_cell) + replicated <- .ppp_rowsum(counts$replicated[, rows, , drop = FALSE], person_cell) + weight <- 1 / sqrt(rowSums(observed)) + weight[!is.finite(weight)] <- 0 + row_col <- match(cell, cells) + + t2 <- function(s, sigma) { + sum((s - sigma)^2 / sqrt(outer(diag(sigma), diag(sigma)))) + } + t(vapply(seq_len(dim(expected)[1]), function(d) { + expected_wide <- wide(expected[d, , ]) + # the average multinomial covariance of a participant: diag(n) - n n' / N + # summed over the participant's rows, which are independent + within <- diag(colMeans(expected_wide), length(col_cell)) + for (j in seq_along(cells)) { + in_cell <- col_cell == j + scaled <- matrix(counts$expected[d, rows[row_col == j], col_cat[in_cell]], + ncol = sum(in_cell)) * weight[row_col == j] + within[in_cell, in_cell] <- within[in_cell, in_cell] - crossprod(scaled) / n_people + } + # the expected value of stats::cov() of the counts (divisor T - 1); the + # (T - 1) / T of Klauer (2010) belongs to the divisor T + sigma <- stats::cov(expected_wide) + within + c(observed = t2(observed_cov, sigma), + replicated = t2(stats::cov(wide(replicated[d, , ])), sigma)) + }, numeric(2))) +} diff --git a/R/summary.R b/R/summary.R index f6bfa6151..dc3406ff3 100644 --- a/R/summary.R +++ b/R/summary.R @@ -177,13 +177,16 @@ summarise_links <- function(links) { } # the formula as text, one element per parameter; with a model, also the -# parameters the user left out and the constants the model fixes them to +# parameters the user left out and the constants the model fixes them to, but +# not the formulas the model generates for parameters the user never wrote summarise_formula <- function(formula, model = NULL) { if (!is.null(model)) { + generated <- setdiff(generated_parameters(model), names(formula)) formula <- suppressMessages(add_missing_parameters(model, formula)) fixpars <- model$fixed_parameters fixpars <- fixpars[names(fixpars) %in% names(model$parameters)] formula[names(fixpars)] <- fixpars + formula <- formula[setdiff(names(formula), generated)] } formula_lines(formula) } diff --git a/README.Rmd b/README.Rmd index dc284c850..2034bdd61 100644 --- a/README.Rmd +++ b/README.Rmd @@ -94,10 +94,22 @@ articles. bmm::print_pretty_models_md(group = "Choices and response times") ``` +**Processing-tree models.** Participants give categorical responses, and +you describe the discrete processes that lead to each response as a +multinomial processing tree: the two-high-threshold model of recognition, +the pair-clustering model of free recall, or a tree of your own. The model is +fit to response counts per participant and condition. See the +[MPT article](https://popov-lab.github.io/bmm/articles/bmm_mpt.html). + +```{r, results = "asis", echo = FALSE} +bmm::print_pretty_models_md(group = "Processing-tree models") +``` + ```{r other-groups, results = "asis", echo = FALSE} known <- c( "Continuous reproduction", "Categorical recall and n-AFC decisions", - "Detection, recognition and confidence judgments", "Choices and response times" + "Detection, recognition and confidence judgments", "Choices and response times", + "Processing-tree models" ) other <- setdiff(unique(bmm:::model_registry()$group), known) if (length(other) > 0) { diff --git a/README.md b/README.md index f8dc319ec..a65151c55 100644 --- a/README.md +++ b/README.md @@ -13,7 +13,7 @@ badge](https://popov-lab.r-universe.dev/badges/bmm)](https://popov-lab.r-univers [![R-CMD-check](https://github.com/popov-lab/bmm/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/popov-lab/bmm/actions/workflows/R-CMD-check.yaml) [![test-coverage](https://github.com/popov-lab/bmm/actions/workflows/test-coverage.yaml/badge.svg)](https://github.com/popov-lab/bmm/actions/workflows/test-coverage.yaml) [![downloads](https://cranlogs.r-pkg.org/badges/bmm)](https://cran.r-project.org/package=bmm) -[![Dependencies](https://img.shields.io/badge/dependencies-14/25-orange?style=flat)](#) +[![Dependencies](https://img.shields.io/badge/dependencies-15/25-orange?style=flat)](#) *Bayesian Measurement Models for behavioral research in R* @@ -103,6 +103,17 @@ articles. - [`ezdm()`](https://popov-lab.github.io/bmm/reference/ezdm.html): EZ-Diffusion Model +**Processing-tree models.** Participants give categorical responses, and +you describe the discrete processes that lead to each response as a +multinomial processing tree: the two-high-threshold model of +recognition, the pair-clustering model of free recall, or a tree of your +own. The model is fit to response counts per participant and condition. +See the [MPT +article](https://popov-lab.github.io/bmm/articles/bmm_mpt.html). + +- [`mpt()`](https://popov-lab.github.io/bmm/reference/mpt.html): + Multinomial Processing Tree (MPT) models + The tasks bmm covers and the models for each `bmm_models()` prints the same list in R, and `?modelname` (for example diff --git a/_pkgdown.yml b/_pkgdown.yml index 57ede022c..f75407410 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -46,6 +46,7 @@ reference: - "conditional_effects.bmmfit" - "emmeans-bmmfit" - "pp_check.bmmfit" + - "ppp" - "posterior_epred.bmmfit" - title: "Model evaluation & post-processing" desc: "ROC/AUC, quantile-transformed ROC, and latent-distribution analysis for signal detection models" @@ -109,6 +110,7 @@ articles: desc: Fitting models to categorical decision tasks (e.g. cued or free recall, n-AFC, etc) contents: - articles/bmm_m3 + - articles/bmm_mpt - title: Signal detection & recognition desc: Fitting signal detection models and analysing ROC/AUC diff --git a/man/figures/README-task-map.png b/man/figures/README-task-map.png index 3f2d12c87..64fdf12f6 100644 Binary files a/man/figures/README-task-map.png and b/man/figures/README-task-map.png differ diff --git a/man/mpt.Rd b/man/mpt.Rd index 5d74227d9..0cd6541f3 100644 --- a/man/mpt.Rd +++ b/man/mpt.Rd @@ -22,7 +22,7 @@ declared impossible.} \item{tree_id}{Character. Name of the data column whose values name the tree each observation belongs to; the values must match the tree names. Can be omitted for single-tree models. Effects of experimental conditions belong -in the parameter formulas, not here — see Details.} +in the parameter formulas, not here; see Details.} \item{covariates}{Character vector. Names of data columns that appear in branch expressions but are not latent parameters, for example @@ -58,13 +58,18 @@ files are not read: in MPTinR and TreeBUGS a character vector names a restrictions file, in \code{bmm} it holds the restrictions. To use such a file, pass \code{readLines(path)}; blank lines and lines starting with \verb{#} are skipped, and a \verb{#} after a restriction starts a comment, as in MPTinR. -Restrictions are substituted into the branch expressions before the -parameters are identified, so a restricted parameter is not part of the -model and cannot be a simplex member. Order constraints (\code{"Do > Dn"}) are -not supported here; see Details.} +A restricted parameter is removed from the model; the Details say how +this differs from fixing a parameter in the formula. Order constraints +(\code{"Do > Dn"}) are not supported here; see Details.} -\item{links}{Character. The link function for all latent probability -parameters: \code{"logit"} (default) or \code{"probit"}.} +\item{links}{The link function of the latent probability parameters: +\code{"logit"} (default) or \code{"probit"}. A single string sets every parameter; +a named list, \code{links = list(D = "probit")}, sets the named parameters +and leaves the others at \code{"logit"}. Members of a \code{simplex} group and +their \code{praw} components take no link of their own. A parameter switched +to the other link, in the list or afterwards with +\code{model$links$D <- "probit"}, takes that link's default priors (see +Details).} \item{...}{used internally for testing, ignore it} } @@ -74,8 +79,9 @@ An object of class \code{bmmodel} \description{ Multinomial Processing Tree (MPT) models measure the probabilities of discrete latent cognitive states from categorical responses. A model is -specified as a set of trees (one per distinct branch structure) created -with \code{\link[=mpt_tree]{mpt_tree()}}. Branch expressions that terminate in the same response +specified as a set of trees created with \code{\link[=mpt_tree]{mpt_tree()}}, one per distinct +branch structure, such as one tree for targets and one for lures in a +recognition test. Branch expressions that terminate in the same response category are summed within each tree, and parameters shared across trees are equated by giving them the same name. All latent probability parameters are estimated on an unconstrained latent scale (logit or @@ -98,11 +104,13 @@ column per response category. The count columns must be named exactly like the branch names given to \code{\link[=mpt_tree]{mpt_tree()}} (\code{old} and \code{new} in the example below); \code{bmm()} finds them by these names, so no response variable is passed. The \code{tree_id} column names the tree of each row, and -the counts of a row must come from trials of that tree. +the counts of a row must come from trials of that tree. Wide tables with +one column per category of an EQN file, as MPTinR and TreeBUGS take them, +are reshaped with \code{\link[=mpt_data_reshape]{mpt_data_reshape()}}. -A separate tree is needed only when the \emph{branch expressions} differ — that +A separate tree is needed only when the \emph{branch expressions} differ, that is, when the trial determines which latent processes apply or which -response categories are reachable (old versus new probes, a response set +response categories are reachable (targets versus lures, a response set with or without distractors, inclusion versus exclusion instructions). Experimental factors whose \emph{effect on the parameters} you want to estimate belong in the parameter formulas, exactly as in every other \code{bmm} model: @@ -111,22 +119,36 @@ column can serve both when a factor happens to change the structure as well. When several levels of a factor share a branch structure, add a column naming the tree and keep the factor for the formulas. -Parameters that receive a non-linear predictor formula (a formula whose -right-hand side references other formula parameters, e.g. -\code{D ~ Dmax * (1 - exp(-rate * ptime))}) are not transformed by the link -function: the user-supplied expression must keep the parameter within -(0, 1). Sub-parameters of such formulas (e.g., \code{Dmax} and \code{rate}) are -estimated on the identity scale with \code{normal(0, 1)} default priors. Their -random-effect SDs keep brms's default prior, \code{student_t(3, 0, 2.5)}, -because no single rate fits every scale such a formula can give them. +\strong{Equating and fixing parameters.} Parameters are equated across trees +by giving them the same name. Beyond that, there are two mechanisms, and +they differ in when the constraint is decided. A \emph{restriction} +(\code{restrictions = c("Dn = Do", "g = 0.5")}) is applied when the model is +built: it is substituted into the branch expressions before the +parameters are identified, so the restricted parameter is not part of +the model, gets no formula and no prior, cannot be a simplex member, and +is listed by \code{print(model)} under \verb{Restrictions:} but not by +\code{\link[=parameter_info]{parameter_info()}}. A \emph{constant in the formula} (\code{bmf(D ~ 1, g = 0.5)}) +is applied at fit time: the parameter stays part of the model, is listed +by \code{print()} under \verb{Fixed:} and by \code{\link[=parameter_info]{parameter_info()}} with its value, and +is freed again by giving it a formula. The value is mapped to the latent +scale when the constant prior is built, so \code{\link[=default_prior]{default_prior()}} shows +\code{constant(0)} for a guessing rate of 0.5 under the logit link, and +\code{summary()} lists it as \verb{g_Intercept 0.00} under "Constant Parameters"; +\code{\link[=prior_info]{prior_info()}} and \code{\link[=parameter_info]{parameter_info()}} show the probability, 0.5. Members +of a \code{simplex} group and their \code{praw} components cannot be fixed this +way; fix parameters outside the group instead. As a rule of thumb, a +constant that defines the model goes into \code{restrictions}, a constant you +compare against a free version goes into the formula; the section +"Constraining parameters" of the MPT article walks through both. -Under the logit link, intercepts and regression coefficients get a -\code{logistic(0, 1)} prior and random-effect SDs get \code{exponential(1)}. Under -the probit link, intercepts and coefficients get \code{normal(0, 1)} and SDs -get \code{exponential(2)}. Group-level correlation matrices get \code{lkj(2)}. -Coefficients of the sub-parameters above get \code{normal(0, 0.5)}. A parameter -switched to the other link after construction (\code{model$links$D <- "probit"}) takes that link's priors unless its default prior was -replaced. \code{\link[=default_prior]{default_prior()}} lists the rows a given formula produces. +\strong{Default priors.} Under the logit link, intercepts and regression +coefficients get a \code{logistic(0, 1)} prior and random-effect SDs get +\code{exponential(1)}. Under the probit link, intercepts and coefficients get +\code{normal(0, 1)} and SDs get \code{exponential(2)}. Group-level correlation +matrices get \code{lkj(2)}. A parameter switched to the other link after +construction (\code{model$links$D <- "probit"}) takes that link's priors +unless its default prior was replaced. \code{\link[=default_prior]{default_prior()}} lists the rows +a given formula produces. The default prior of a simplex group of K members is the uniform (symmetric Dirichlet(1)) distribution on the simplex: each member has @@ -149,32 +171,22 @@ on the first member therefore scales all later members by the same factor, leaving their ratios unchanged. Read such effects on the probability scale with \code{\link[=native_parameters]{native_parameters()}}. -A parameter can also be fixed to a probability at fit time, in the -formula: \code{bmf(D ~ 1, g = 0.5)}. Unlike a restriction, the parameter stays -part of the model and can be freed again by giving it a formula. The -value is mapped to the latent scale when the constant prior is built, so -\code{\link[=default_prior]{default_prior()}} shows \code{constant(0)} for a guessing rate of 0.5 under -the logit link. For the same reason, \code{summary()} lists a fixed \code{g = 0.5} -as \verb{g_Intercept 0.00} under "Constant Parameters"; \code{\link[=prior_info]{prior_info()}} and -\code{\link[=parameter_info]{parameter_info()}} show the probability, 0.5. Members of a \code{simplex} -group and their \code{praw} components cannot be fixed this way; fix parameters -outside the group instead. - -Printing the model ends with an identifiability check for intercept-only -formulas. It compares the number of free parameters with the degrees of -freedom (response categories minus 1, summed over trees) and computes the -rank of the Jacobian of the category probabilities with respect to the -free parameters, from exact derivatives at five interior test values (and -five test values of each covariate); parameters fixed in the formula enter -at their values. A rank below the number of free parameters means that a -combination of the named parameters is not identified by the branch -expressions, even when the count passes. The check is local: it holds at -the test values, not at the boundaries of the parameter space. A test -value with a non-finite derivative is left out, and the printout says how -many; if \code{\link[stats:D]{stats::D()}} cannot differentiate a branch expression, the rank -is not computed and the count decides. The rank sees the tree parameters -only: after a non-linear formula, the rank line names that formula's -parameter as not analysed. +\strong{Identifiability.} Printing the model ends with an identifiability +check for intercept-only formulas. It compares the number of free +parameters with the degrees of freedom (response categories minus 1, +summed over trees) and computes the rank of the Jacobian of the category +probabilities with respect to the free parameters, from exact derivatives +at five interior test values (and five test values of each covariate); +parameters fixed in the formula enter at their values. A rank below the +number of free parameters means that a combination of the named +parameters is not identified by the branch expressions, even when the +count passes. The check is local: it holds at the test values, not at the +boundaries of the parameter space. A test value with a non-finite +derivative is left out, and the printout says how many; if \code{\link[stats:D]{stats::D()}} +cannot differentiate a branch expression, the rank is not computed and +the count decides. The rank sees the tree parameters only: after a +non-linear formula, the rank line names that formula's parameter as not +analysed. \code{bmm()} warns when the rank is deficient and every formula is intercept-only. When a formula has a predictor or is non-linear, it only @@ -185,10 +197,10 @@ test give three free parameters and rank 2 of 3, so one combination of \code{Do}, \code{Dn} and \code{g} cannot be estimated; a second condition that varies the response bias identifies them. -\code{summary()} reports intercepts and regression coefficients on the latent -(logit or probit) scale. \code{\link[=native_parameters]{native_parameters()}} returns the posterior -draws of every parameter on the probability scale for each observed -combination of the predictors, e.g. +\strong{Reading the estimates.} \code{summary()} reports intercepts and regression +coefficients on the latent (logit or probit) scale. \code{\link[=native_parameters]{native_parameters()}} +returns the posterior draws of every parameter on the probability scale +for each observed combination of the predictors, e.g. \code{native_parameters(fit, re_formula = NA, summary = TRUE)}; a difference between conditions is the difference of these draws. With \code{re_formula = NA} these are the probabilities of the median participant @@ -203,20 +215,31 @@ corresponds to the latent-trait MPT model. For posterior predictive checks, \code{\link[brms:posterior_predict]{brms::posterior_predict()}} returns simulated counts (draws x rows x categories) and \code{\link[brms:posterior_epred]{brms::posterior_epred()}} -the expected counts; the MPT article computes the T1 statistic and its -posterior predictive p-value from them. With \code{newdata}, +the expected counts; \code{\link[=ppp]{ppp()}} computes the T1 and T2 statistics of Klauer +(2010) and their posterior predictive p-values from them. With \code{newdata}, \code{\link[brms:posterior_predict]{brms::posterior_predict()}} needs the columns \code{nTrials}, \verb{Idx_} and \verb{Poss_} that \code{check_data()} builds; \code{fit$data} contains them. +\strong{Non-linear formulas.} Parameters that receive a non-linear predictor +formula (a formula whose right-hand side references other formula +parameters, e.g. \code{D ~ Dmax * (1 - exp(-rate * ptime))}) are not +transformed by the link function: the user-supplied expression must keep +the parameter within (0, 1). Sub-parameters of such formulas (e.g., +\code{Dmax} and \code{rate}) are estimated on the identity scale with +\code{normal(0, 1)} default priors on their intercepts and \code{normal(0, 0.5)} +on their regression coefficients. Their random-effect SDs keep brms's +default prior, \code{student_t(3, 0, 2.5)}, because no single rate fits every +scale such a formula can give them. + Order constraints between parameters (\code{Do > Dn}) are expressed by reparameterizing the larger parameter in the model formula, e.g. \code{bmf(Do ~ Dn + (1 - Dn) * inv_logit(phi), Dn ~ 1, phi ~ 1)}; the section "Order constraints" of the MPT article walks through the recipe. -Parameter and response category names must start with a letter and may -contain only letters and digits, because brms does not allow underscores -or dots in non-linear parameter names. +\strong{Names.} Parameter and response category names must start with a +letter and may contain only letters and digits, because brms does not +allow underscores or dots in non-linear parameter names. All latent parameters are non-linear parameters of the underlying brms model, and brms has no separate \code{Intercept} prior class for those: their @@ -230,17 +253,19 @@ formula produces. \examples{ \dontrun{ # two-high-threshold (2HTM) model of recognition memory -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)" -)) +# one tree per probe type: targets (studied items) and lures (new items), +# each ending in the responses "old" and "new" +tree_target <- mpt_tree("target", + old ~ D + (1 - D) * g, + new ~ (1 - D) * (1 - g) +) +tree_lure <- mpt_tree("lure", + old ~ (1 - D) * g, + new ~ D + (1 - D) * (1 - g) +) model <- mpt( - trees = list(tree_old, tree_new), + trees = list(tree_target, tree_lure), tree_id = "item_type" ) @@ -248,7 +273,7 @@ model <- mpt( # one row per participant x tree; count columns named after the branches data <- data.frame( id = rep(1:20, each = 2), - item_type = rep(c("old", "new"), 20) + item_type = rep(c("target", "lure"), 20) ) counts <- t(sapply(data$item_type, function(cond) { rmpt( diff --git a/man/mpt_data_reshape.Rd b/man/mpt_data_reshape.Rd new file mode 100644 index 000000000..3ecedc877 --- /dev/null +++ b/man/mpt_data_reshape.Rd @@ -0,0 +1,145 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/model_mpt_data.R +\name{mpt_data_reshape} +\alias{mpt_data_reshape} +\title{Reshape wide MPT frequency tables to the long format of \code{bmm()}} +\usage{ +mpt_data_reshape(data, model, id = NULL, columns = c("names", "mptinr")) +} +\arguments{ +\item{data}{A data frame or matrix with one row per participant. The count +columns are read as described under \code{columns}; with \code{columns = "names"}, +all other columns except \code{id} are carried to the long data.} + +\item{model}{An \code{mpt} model from \code{\link[=mpt]{mpt()}}, \code{\link[=mpt_from_eqn]{mpt_from_eqn()}} or +\code{\link[=mpt_from_string]{mpt_from_string()}}. For a model from \code{\link[=mpt_from_eqn]{mpt_from_eqn()}} the \code{mpt_renaming} +attribute maps the labels of the EQN file onto the response categories +of the model. For any other model, the columns of \code{data} are named +\verb{.}, e.g. \code{target.old}, with the tree names and response +categories of the model.} + +\item{id}{Character. The name of the participant column in \code{data}, or +\code{NULL} (the default) to number the rows of \code{data}, as \code{id = seq_len(nrow(data))} +in a column named \code{id}. The participants must be unique.} + +\item{columns}{Character. How the count columns of \code{data} are identified: +\code{"names"} (the default) by their names, as in TreeBUGS data; \code{"mptinr"} +by their position, as in MPTinR data. See Details.} +} +\value{ +A data frame with one row per participant and tree, participants +first. Its columns are the participant column (\code{id}), the tree column +named by the \code{tree_id} of the model (absent in a model without it), the +carried columns, and one integer count column per response category of +the model. +} +\description{ +MPTinR and TreeBUGS take one row per participant and one column per +response category of the EQN file. \code{\link[=bmm]{bmm()}} takes one row per participant +and tree, a column naming the tree, and one count column per \emph{shared} +response category (see \code{\link[=mpt]{mpt()}}). \code{mpt_data_reshape()} reshapes the first +layout into the second for a given model. +} +\details{ +\strong{Named columns (\code{columns = "names"}).} For a model from +\code{\link[=mpt_from_eqn]{mpt_from_eqn()}}, a column is matched to a category by the EQN label +(\code{hit}), or, when trees share labels, by the tree name and label +(\code{target.hit}, with the bmm name of the tree); the qualified name wins +where both exist. A label used in several trees must be qualified. +Labels that the \code{categories} argument of \code{\link[=mpt_from_eqn]{mpt_from_eqn()}} merged onto one +response category are summed. All columns of \code{data} that are neither +count columns nor \code{id} are carried unchanged to every long row, e.g. +between-participant covariates. + +\strong{Positional columns (\code{columns = "mptinr"}).} MPTinR data for an EQN +model have no column names. Every column of \code{data} except \code{id} is read as +a count, in the order of \code{check.mpt()} in MPTinR (version 1.14.1), +elements \code{eqn.order.trees} and \code{eqn.order.categories}: trees sorted by +their label in the EQN file and, within each tree, the category labels of +the file sorted. No other column can be carried; add covariates +afterwards, e.g. with \code{\link[=merge]{merge()}} on the participant column. The help page +of \code{check.mpt()} states only that these elements give the order; the +sorting rule is that of its internal \code{.get.EQN.model.order()}, which +sorts the columns that \code{read.table()} returns for the file: numerically +if all tree labels (all category labels) are numbers, else as character +strings in the collation of the R session. \code{mpt_data_reshape()} applies the +same rule, so run it in the session where MPTinR's order was determined +if the labels mix upper and lower case. This mode needs a model from +\code{\link[=mpt_from_eqn]{mpt_from_eqn()}}: MPTinR's own model format orders the columns by the +lines of the model text, and \code{\link[=mpt]{mpt()}} reorders the branches of every tree +to the categories of the first, so the line order cannot be recovered +from the model. For other models, name the columns \verb{.} +and use \code{columns = "names"}. + +\strong{Impossible categories.} A response category that a tree cannot +produce (\code{impossible} in \code{\link[=mpt_tree]{mpt_tree()}}, \code{\link[=mpt_from_string]{mpt_from_string()}} and +\code{\link[=mpt_from_eqn]{mpt_from_eqn()}}) has no column in \code{data} and gets 0 in the long rows of +that tree, which \code{\link[=bmm]{bmm()}} accepts. + +\strong{Counts.} The count columns must be numeric, without missing values, +and hold non-negative whole numbers; the long data hold them as integers. + +\strong{Checks.} The function cannot tell a correct reshape from one that +swaps columns of equal meaning. Run \code{\link[=bmm_data_check]{bmm_data_check()}} on the result, and +compare one participant's long rows with the wide row (see the Examples +of \code{\link[=mpt_from_eqn]{mpt_from_eqn()}}). +} +\examples{ +eqn_file <- tempfile(fileext = ".eqn") +writeLines(c( + "6", + "target hit D_o", + "target hit (1-D_o)*g", + "target miss (1-D_o)*(1-g)", + "lure fa (1-D_n)*g", + "lure cr D_n", + "lure cr (1-D_n)*(1-g)" +), eqn_file) +model <- mpt_from_eqn( + eqn_file, + restrictions = "D_n = D_o", + categories = c(hit = "yes", fa = "yes", miss = "no", cr = "no"), + tree_id = "item_type" +) + +# TreeBUGS data: one named column per category label of the EQN file. +# Columns that are no counts, here a between-participant factor, are +# carried to the long data +wide <- data.frame( + subject = c("s1", "s2"), group = c("a", "b"), + hit = c(30, 25), miss = c(10, 15), fa = c(8, 12), cr = c(32, 28) +) +long <- mpt_data_reshape(wide, model, id = "subject") +long +bmm_data_check(bmf(Do ~ 1, g ~ 1), long, model) + +# MPTinR data: unnamed columns, sorted by tree, then by label (cr, fa for +# tree lure, then hit, miss for tree target) +wide_mptinr <- rbind(c(32, 8, 30, 10), c(28, 12, 25, 15)) +identical( + mpt_data_reshape(wide_mptinr, model, columns = "mptinr"), + mpt_data_reshape(wide[c("cr", "fa", "hit", "miss")], model) +) + +# trees with response categories that only one tree can produce +model_pc <- mpt( + list( + mpt_tree("pairs", list(E1 = "c * r", E2 = "1 - c * r"), + impossible = "F1"), + mpt_tree("singles", list(F1 = "u", E2 = "1 - u"), + impossible = "E1") + ), + tree_id = "item_type" +) +wide_pc <- data.frame( + pairs.E1 = c(20, 22), pairs.E2 = c(20, 18), + singles.F1 = c(25, 24), singles.E2 = c(15, 16) +) +mpt_data_reshape(wide_pc, model_pc) +} +\seealso{ +\code{\link[=mpt_from_eqn]{mpt_from_eqn()}} and \code{\link[=mpt_from_string]{mpt_from_string()}} for importing the model +itself; \code{\link[=mpt]{mpt()}} for the data layout \code{\link[=bmm]{bmm()}} expects; \code{\link[=bmm_data_check]{bmm_data_check()}} +to check the result before fitting. +} +\keyword{transform} diff --git a/man/mpt_from_eqn.Rd b/man/mpt_from_eqn.Rd index 0be10fde5..57a409ef6 100644 --- a/man/mpt_from_eqn.Rd +++ b/man/mpt_from_eqn.Rd @@ -59,7 +59,7 @@ categories. When the trees do not share their categories and categories each tree lacks.} \item{tree_names}{A named character vector mapping tree labels of the EQN -file onto tree names in bmm, e.g. \code{c("1" = "old", "2" = "new")}. bmm tree +file onto tree names in bmm, e.g. \code{c("1" = "target", "2" = "lure")}. bmm tree names must start with a letter and contain only letters, digits, or underscores, so files with numeric tree labels need this mapping. Labels not listed keep their name.} @@ -109,9 +109,11 @@ declare it with the \code{impossible} argument instead. parameter names, parameter and response category names are sanitized: a leading \verb{_} prefix is removed from category names that are not mapped with \code{categories}, and all remaining underscores and dots are -stripped. Names that would clash after this are an error; names are -case-sensitive, so \code{P_G1} and \code{P_g1} stay distinct (\code{PG1}, \code{Pg1}). Tree -labels are kept, unless \code{tree_names} maps them. The formulas passed to +stripped. Two labels of the file, or two different values within one +tree, that become the same name are an error; mapped values that become +the same name in different trees are that one shared category. Names +are case-sensitive, so \code{P_G1} and \code{P_g1} stay distinct (\code{PG1}, \code{Pg1}). +Tree labels are kept, unless \code{tree_names} maps them. The formulas passed to \code{\link[=bmm]{bmm()}} use the bmm names (\code{bmf(Do ~ 1)} for a file parameter \code{D_o}; \code{bmf()} is short for \code{\link[=bmmformula]{bmmformula()}}); the \code{parameter} rows of the renaming map translate them, and \code{print(model)} lists them. Of the @@ -131,33 +133,24 @@ and tree: a column named by \code{tree_id} holding the bmm tree names, and one count column per shared response category, named by \code{bmm_name} in the renaming map (see \code{\link[=mpt]{mpt()}}). A category that is \code{impossible} in a tree gets 0 (or \code{NA}) in the rows of that tree. TreeBUGS data are wide: -one row per participant, one column per category label of the EQN file, -and no participant column, so add an \code{id} column first. The \code{category} -rows of the renaming map hold every tree with its labels and bmm names, -so the wide data can be reshaped by name: for each tree, start from a -block of zeros over all shared categories and fill in the columns of -that tree (see Examples). MPTinR data for an EQN model are unnamed, and -their columns do not follow the file: the trees and, within each tree, -the category labels are sorted (alphabetically, or numerically when they -are numbers); \code{MPTinR::check.mpt("model.eqn")} returns the order in its -elements \code{eqn.order.trees} and \code{eqn.order.categories}. (For MPTinR's own -model format the columns follow the lines instead; see -\code{\link[=mpt_from_string]{mpt_from_string()}}.) Name the columns by tree and label, -\code{paste(tree, label, sep = ".")} with the bmm tree names, because trees -may share labels, and select them with the same key (see Examples). - -If \code{categories} merged several labels of a tree onto one name, sum their -columns in the wide data into the first label's column and drop the map -rows of the other labels (\code{map[!duplicated(map[c("tree", "bmm_name")]), ]}) -before reshaping; otherwise the recipe fails with "duplicate subscripts". +one row per participant, one column per category label of the EQN file. +MPTinR data for an EQN model are unnamed, and their columns do not follow +the file: the trees and, within each tree, the category labels are sorted +(alphabetically, or numerically when all labels in the file are numbers); +\code{MPTinR::check.mpt("model.eqn")} returns the order in its elements +\code{eqn.order.trees} and \code{eqn.order.categories}. (For MPTinR's own model +format the columns follow the lines instead; see \code{\link[=mpt_from_string]{mpt_from_string()}}.) +\code{\link[=mpt_data_reshape]{mpt_data_reshape()}} reshapes both layouts with the renaming map, sums the +labels that \code{categories} merged onto one response category, and gives 0 to +the categories a tree cannot produce (see Examples). bmm cannot tell a reshape that swaps tree labels or category columns -from a correct one: such data fit without an error. The row sums, -\code{rowSums(long[cats])}, are the trials per row; they cannot show swapped -columns, and they show swapped tree labels only when the trees differ -in their number of trials. Comparing the response proportions per tree -with the design (old items must get more "yes" responses than new -items) catches swapped tree labels and response columns swapped in +from a correct one: such data fit without an error. The row sums of the +count columns are the trials per row; they cannot show swapped columns, +and they show swapped tree labels only when the trees differ in their +number of trials. Comparing the response proportions per tree +with the design (targets must get more "yes" responses than lures) +catches swapped tree labels and response columns swapped in every tree, but not always a column slipped or a map entry wrong in one tree, nor a swap between two trees of the same kind; so also compare one participant's long rows with the wide row. \code{\link[=bmm_data_check]{bmm_data_check()}} runs @@ -173,12 +166,12 @@ it was built from. eqn_file <- tempfile(fileext = ".eqn") writeLines(c( "6", - "old hit D_o", - "old hit (1-D_o)*g", - "old miss (1-D_o)*(1-g)", - "new fa (1-D_n)*g", - "new cr D_n", - "new cr (1-D_n)*(1-g)" + "target hit D_o", + "target hit (1-D_o)*g", + "target miss (1-D_o)*(1-g)", + "lure fa (1-D_n)*g", + "lure cr D_n", + "lure cr (1-D_n)*(1-g)" ), eqn_file) # restrictions use the file's names; the formulas use the bmm names model <- mpt_from_eqn( @@ -190,54 +183,33 @@ model <- mpt_from_eqn( attr(model, "mpt_renaming") # TreeBUGS data: one row per participant, columns named by the category -# labels of the EQN file, no participant column. bmm needs one row per -# participant and tree, columns named by bmm_name: reshape by name through -# the map, starting each tree from zeros so that categories the tree -# cannot produce get 0 +# labels of the EQN file, no participant column. mpt_data_reshape() reshapes +# them to bmm's one row per participant and tree, through the map; the +# categories a tree cannot produce get 0 wide <- data.frame( hit = c(30, 25), miss = c(10, 15), fa = c(8, 12), cr = c(32, 28) ) -wide$id <- seq_len(nrow(wide)) -map <- attr(model, "mpt_renaming") -map <- map[map$kind == "category", ] -cats <- unique(map$bmm_name) -long <- do.call(rbind, lapply(split(map, map$tree), function(tree_map) { - counts <- as.data.frame(matrix(0, nrow(wide), length(cats))) - names(counts) <- cats - counts[tree_map$bmm_name] <- wide[tree_map$file_name] - data.frame(id = wide$id, item_type = tree_map$tree[1], counts) -})) +long <- mpt_data_reshape(wide, model) long # the row sums are the trials per row (40 here). The proportions catch -# swapped tree labels and swapped response columns (old items must get -# more "yes" than new items), but not always a column slipped in one +# swapped tree labels and swapped response columns (targets must get +# more "yes" than lures), but not always a column slipped in one # tree: also compare one participant's long rows with the wide row -rowSums(long[cats]) -tapply(long$yes / rowSums(long[cats]), long$item_type, mean) +rowSums(long[c("yes", "no")]) +tapply(long$yes / rowSums(long[c("yes", "no")]), long$item_type, mean) long[long$id == 1, ] wide[1, ] bmm_data_check(bmf(Do ~ 1, g ~ 1), long, model) # MPTinR data for an EQN model are unnamed and sorted by tree, then by # label: MPTinR::check.mpt(eqn_file)$eqn.order.categories lists the -# labels (here cr, fa for tree new, then hit, miss for tree old). Name the -# columns by tree and label, so that a label used in several trees stays -# apart, and select with the same key +# labels (here cr, fa for tree new, then hit, miss for tree old) wide_mptinr <- rbind(c(32, 8, 30, 10), c(28, 12, 25, 15)) -colnames(wide_mptinr) <- c("new.cr", "new.fa", "old.hit", "old.miss") -long2 <- do.call(rbind, lapply(split(map, map$tree), function(tree_map) { - counts <- as.data.frame(matrix(0, nrow(wide_mptinr), length(cats))) - names(counts) <- cats - key <- paste(tree_map$tree, tree_map$file_name, sep = ".") - counts[tree_map$bmm_name] <- wide_mptinr[, key] - data.frame(id = seq_len(nrow(wide_mptinr)), item_type = tree_map$tree[1], - counts) -})) -identical(long2, long) +identical(mpt_data_reshape(wide_mptinr, model, columns = "mptinr"), long) } \seealso{ \code{\link[=mpt_from_string]{mpt_from_string()}} for model text in MPTinR's own format; \code{\link[=mpt]{mpt()}} and \code{\link[=mpt_tree]{mpt_tree()}} for building the trees in R, including the data layout \code{bmm()} expects. } -\keyword{transform} +\keyword{bmmodel} diff --git a/man/mpt_from_string.Rd b/man/mpt_from_string.Rd index 2fe17ce48..6aac672b1 100644 --- a/man/mpt_from_string.Rd +++ b/man/mpt_from_string.Rd @@ -6,12 +6,12 @@ \usage{ mpt_from_string( text, + restrictions = NULL, tree_names, categories = NULL, tree_id = NULL, covariates = NULL, simplex = NULL, - restrictions = NULL, links = "logit", impossible = NULL ) @@ -20,6 +20,11 @@ mpt_from_string( \item{text}{Character. The model definition, either as a single string or as a vector of lines (e.g., from \code{\link[=readLines]{readLines()}}).} +\item{restrictions}{Parameter restrictions as in \code{\link[=mpt]{mpt()}}, written with the +parameter names of the model text, e.g. \code{"g = 0.5"} or \code{list(g = 0.5)}. +They mirror the \code{restrictions} argument of MPTinR and TreeBUGS. +Restriction files are not read: pass their lines (see Examples).} + \item{tree_names}{Character vector with one name per tree block. The names are matched to the blocks by position: the first name goes to the first block, and so on. bmm cannot check this order, so names given in a @@ -45,8 +50,8 @@ each observation belongs to (see \code{\link[=mpt]{mpt()}}). Add this column to yourself: for each row it holds the name, as given in \code{tree_names}, of the tree the counts come from.} -\item{covariates, simplex, restrictions, links}{passed to \code{\link[=mpt]{mpt()}}. -\code{print(model)} lists the symbols that became parameters.} +\item{covariates, simplex, links}{passed to \code{\link[=mpt]{mpt()}}. \code{print(model)} lists +the symbols that became parameters.} \item{impossible}{A named list with the response categories that cannot occur in a tree, one element per such tree: @@ -101,9 +106,8 @@ this format are wide and positional: one row per participant (or one vector for aggregated counts) and one unnamed column per line of the model file, in file order, because each line is one response category. (For EQN files MPTinR sorts the columns instead; see \code{\link[=mpt_from_eqn]{mpt_from_eqn()}}.) -Name the columns by tree and category first, then stack one block per -tree (see Examples); when a tree lacks some categories, start each tree -from zeros, as in the Examples of \code{\link[=mpt_from_eqn]{mpt_from_eqn()}}. +Name the columns \verb{.}; \code{\link[=mpt_data_reshape]{mpt_data_reshape()}} then stacks the +trees and gives 0 to the categories a tree cannot produce (see Examples). \strong{Provenance.} The returned model carries an \code{mpt_source} attribute, \code{list(text = )}, with the model definition as one string. @@ -119,7 +123,7 @@ D + (1 - D) * (1 - g) # new " model <- mpt_from_string( model_2htm, - tree_names = c("old", "new"), + tree_names = c("target", "lure"), tree_id = "item_type" ) model @@ -127,7 +131,7 @@ model # fix the guessing rate as in MPTinR / TreeBUGS mpt_from_string( model_2htm, - tree_names = c("old", "new"), + tree_names = c("target", "lure"), tree_id = "item_type", restrictions = "g = 0.5" ) @@ -164,21 +168,17 @@ writeLines(c("# fix the guessing rate", "g = 0.5"), restr_file) restr <- trimws(sub("#.*", "", readLines(restr_file))) model <- mpt_from_string( readLines(model_file), - tree_names = c("old", "new"), + tree_names = c("target", "lure"), tree_id = "item_type", restrictions = restr[nzchar(restr)] ) # MPTinR data: one unnamed column per line of the model file, in file -# order. Name the columns by tree and category, then stack one block per -# tree to get bmm's one row per participant and tree. +# order. Name the columns by tree and category, then let mpt_data_reshape() +# stack the trees to get bmm's one row per participant and tree. wide <- rbind(c(30, 10, 8, 32), c(25, 15, 12, 28)) -colnames(wide) <- c("old.old", "old.new", "new.old", "new.new") -long <- do.call(rbind, lapply(c("old", "new"), function(tree) { - counts <- wide[, paste(tree, c("old", "new"), sep = "."), drop = FALSE] - colnames(counts) <- c("old", "new") - data.frame(id = seq_len(nrow(wide)), item_type = tree, counts) -})) +colnames(wide) <- c("target.old", "target.new", "lure.old", "lure.new") +long <- mpt_data_reshape(wide, model) long # the row sums are the trials per row (40 here), but they cannot show # swapped columns or tree labels: compare one participant with the wide @@ -194,4 +194,4 @@ bmm_data_check(bmf(D ~ 1), long, model) \code{\link[=mpt]{mpt()}} and \code{\link[=mpt_tree]{mpt_tree()}} for building the trees in R, including the data layout \code{bmm()} expects. } -\keyword{transform} +\keyword{bmmodel} diff --git a/man/mpt_tree.Rd b/man/mpt_tree.Rd index 8bcca40ee..3e5ce4d89 100644 --- a/man/mpt_tree.Rd +++ b/man/mpt_tree.Rd @@ -4,23 +4,34 @@ \alias{mpt_tree} \title{Create a tree for a Multinomial Processing Tree (MPT) model} \usage{ -mpt_tree(name, branches, impossible = NULL) +mpt_tree(name, branches, ..., impossible = NULL) } \arguments{ \item{name}{Character. Label for the tree. For models with multiple trees, the label identifies the tree in the data: the column named by the -\code{tree_id} argument of \code{\link[=mpt]{mpt()}} holds one such label per observation.} +\code{tree_id} argument of \code{\link[=mpt]{mpt()}} holds one such label per observation. A +tree stands for one kind of trial, such as targets or lures in a +recognition test, not for a response; the responses are the branch +names.} -\item{branches}{A named list of character strings. Each element gives the -branch probability expression for one response category, and the element -names are the response categories. The expressions can use latent -parameters (e.g., \code{"D + (1 - D) * g"}), numeric constants, and declared -covariates (see \code{\link[=mpt]{mpt()}}).} +\item{branches}{Either a named list of character strings or the first of +several formulas. Each list element (or formula) gives the branch +probability expression for one response category; the element names (or +the left-hand sides) are the response categories. The expressions can use +latent parameters (e.g., \code{D + (1 - D) * g}), numeric constants, and +declared covariates (see \code{\link[=mpt]{mpt()}}).} + +\item{...}{Further formulas of the form \code{category ~ expression}, one per +response category, when \code{branches} is a formula; they take no argument +names. A call uses either +formulas or a named list of strings, not a mix; both spellings give the +same tree (see the examples).} \item{impossible}{Character vector. Response categories that cannot occur in this tree, for example distractor responses when no distractors were presented. Such categories get no branch expression, and the remaining -branches must sum to 1 on their own.} +branches must sum to 1 on their own. It follows \code{...}, so it has to be +named in full.} } \value{ An object of class \code{mpt_tree} @@ -32,6 +43,12 @@ response category. Trees created with \code{mpt_tree()} are combined into a full model specification with \code{\link[=mpt]{mpt()}}. } \details{ +Nothing in the formulas is evaluated: the left-hand side must be a +single response category name, the right-hand side is kept as an +expression (with constant arithmetic folded, see below), and the formula +environment is dropped. A one-sided formula, a left-hand side that is not +a name, and a mix of formulas and strings are errors. + Constant arithmetic such as \code{1/4} or \code{1/(2*2)} is folded into a decimal literal (\code{0.25}) when the tree is created. Stan compiles a bare integer fraction as integer division (\code{1/4 == 0}), which would silently @@ -45,34 +62,40 @@ rather than a literal zero; no probability has to be taken away from the remaining branches to compensate. } \examples{ -tree_old <- mpt_tree( - name = "old", +# the target tree of the two-high-threshold model: targets (studied items) +# end in the responses "old" and "new" +tree_target <- mpt_tree( + "target", + old ~ D + (1 - D) * g, + new ~ (1 - D) * (1 - g) +) +tree_target + +# the same tree with the branches as a named list of strings +tree_target_strings <- mpt_tree( + name = "target", branches = list( old = "D + (1 - D) * g", new = "(1 - D) * (1 - g)" ) ) -tree_old +identical(tree_target, tree_target_strings) # a three-alternative recognition task in which lures are shown on some # trials only: `lure` responses cannot occur on the others tree_lures <- mpt_tree( - name = "lures", - branches = list( - target = "D * b + (1 - D) * (1/3)", - lure = "D * (1 - b) + (1 - D) * (1/3)", - new = "(1 - D) * (1/3)" - ) + "lures", + target ~ D * b + (1 - D) * (1/3), + lure ~ D * (1 - b) + (1 - D) * (1/3), + new ~ (1 - D) * (1/3) ) # the branches of the lure-free tree sum to 1 without the `lure` category; # no probability is reserved for it tree_nolures <- mpt_tree( - name = "nolures", - branches = list( - target = "D + (1 - D) * (1/2)", - new = "(1 - D) * (1/2)" - ), + "nolures", + target ~ D + (1 - D) * (1/2), + new ~ (1 - D) * (1/2), impossible = "lure" ) tree_nolures @@ -82,4 +105,4 @@ tree_nolures # e.g. bmf(D ~ 0 + cond) model <- mpt(list(tree_lures, tree_nolures), tree_id = "tree") } -\keyword{transform} +\keyword{bmmodel} diff --git a/man/mptdist.Rd b/man/mptdist.Rd index 8bf5b3fb8..1195612bb 100644 --- a/man/mptdist.Rd +++ b/man/mptdist.Rd @@ -24,7 +24,11 @@ rmpt( categories, giving the number of observed responses per category. If \code{x} is a named vector, its names are matched to the response categories and the vector is reordered accordingly; an unnamed vector is read in -the order of the branch names of the selected tree.} +the order of the model's response categories, +\code{mpt_model$resp_vars$resp_cats}, which is the order of the first tree +given to \code{\link[=mpt]{mpt()}} (\code{\link[=mpt]{mpt()}} stores the branches of every tree in that +order, whatever order they were written in). \code{rmpt()} returns its +columns in the same order.} \item{pars}{A named vector or list with the values of at least the latent parameters appearing in the branch expressions of the selected tree, on @@ -72,23 +76,23 @@ Unlike the densities of the circular models in this package (\code{dsdm()}, \code{dmpt()} returns the log density by default, matching \code{\link[=dm3]{dm3()}}. } \examples{ -tree_old <- mpt_tree("old", list( +tree_target <- mpt_tree("target", list( old = "D + (1 - D) * g", new = "(1 - D) * (1 - g)" )) -tree_new <- mpt_tree("new", list( +tree_lure <- mpt_tree("lure", list( old = "(1 - D) * g", new = "D + (1 - D) * (1 - g)" )) -model <- mpt(list(tree_old, tree_new), tree_id = "item_type") +model <- mpt(list(tree_target, tree_lure), tree_id = "item_type") dmpt( x = c(35, 15), pars = c(D = 0.7, g = 0.5), - mpt_model = model, tree = "old" + mpt_model = model, tree = "target" ) rmpt( n = 10, size = 50, pars = c(D = 0.7, g = 0.5), - mpt_model = model, tree = "new" + mpt_model = model, tree = "lure" ) } \references{ diff --git a/man/plot.mpt_tree.Rd b/man/plot.mpt_tree.Rd index c65993101..2f93965a0 100644 --- a/man/plot.mpt_tree.Rd +++ b/man/plot.mpt_tree.Rd @@ -77,16 +77,16 @@ labels near the root can be clipped earlier, when the leaf labels leave little room. } \examples{ -tree_old <- mpt_tree("old", list( +tree_target <- mpt_tree("target", list( old = "D + (1 - D) * g", new = "(1 - D) * (1 - g)" )) -plot(tree_old) +plot(tree_target) model <- mpt( list( - tree_old, - mpt_tree("new", list( + tree_target, + mpt_tree("lure", list( old = "(1 - D) * g", new = "D + (1 - D) * (1 - g)" )) @@ -100,4 +100,4 @@ Batchelder, W. H., & Riefer, D. M. (1999). Theoretical and empirical review of multinomial process tree modeling. Psychonomic Bulletin & Review, 6(1), 57-86. https://doi.org/10.3758/BF03210812 } -\keyword{transform} +\keyword{bmmodel} diff --git a/man/pp_check.bmmfit.Rd b/man/pp_check.bmmfit.Rd index c219fd094..1ced21e2f 100644 --- a/man/pp_check.bmmfit.Rd +++ b/man/pp_check.bmmfit.Rd @@ -41,7 +41,9 @@ signal and noise rating distributions are checked separately; pass \code{group = NA} to pool them into a single profile. The \code{\link[=sdt_ranking]{sdt_ranking()}} model shows the rank positions as the categories of the proportion profile; for mixed set sizes, pass the set-size column as \code{group} so each facet has a -homogeneous denominator. The custom-family +homogeneous denominator. For \code{\link[=mpt]{mpt()}} fits, the \code{tree_id} column can be +named even if no parameter formula uses it, giving one facet per tree; +the default pools the trees. The custom-family models \code{\link[=sdt_yn]{sdt_yn()}} and \code{\link[=sdt_mafc]{sdt_mafc()}} delegate to \code{\link[brms:pp_check]{brms::pp_check()}}; for their aggregated counts an overlaid empirical CDF (\code{type = "ecdf_overlay"}, adding \code{group = "stimulus"} for \code{\link[=sdt_yn]{sdt_yn()}}) is usually clearer than the diff --git a/man/ppp.Rd b/man/ppp.Rd new file mode 100644 index 000000000..193cebd29 --- /dev/null +++ b/man/ppp.Rd @@ -0,0 +1,148 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ppp.R +\name{ppp} +\alias{ppp} +\alias{ppp.bmmfit} +\title{Posterior predictive p-values for count models} +\usage{ +ppp(fit, ...) + +\method{ppp}{bmmfit}( + fit, + statistic = c("T1", "T2"), + by = NULL, + ndraws = NULL, + level = c("group", "participant"), + participant = NULL, + ... +) +} +\arguments{ +\item{fit}{A \code{bmmfit} object returned by \code{\link[=bmm]{bmm()}}.} + +\item{...}{Unused.} + +\item{statistic}{Character. \code{"T1"}, \code{"T2"}, or both (the default). If the +participants cannot be identified, the default computes \code{"T1"} only.} + +\item{by}{Character vector naming the columns of the model data that define +the cells. \code{NULL} (default) uses the cells described in Details; +\code{character(0)} pools all rows into one cell. For \code{\link[=mpt]{mpt()}} fits, the +tree column can be named even if it is not a predictor.} + +\item{ndraws}{Integer. Number of posterior draws to use. \code{NULL} (default) +uses all draws.} + +\item{level}{Character. \code{"group"} (default) for the statistics of the whole +data set; \code{"participant"} for one \eqn{T_1} and PPP per participant, as +TreeBUGS reports them, to find the participants the model fits poorly.} + +\item{participant}{Character. The column that identifies participants, +needed for \eqn{T_2} and for \code{level = "participant"}. It must be a +variable of the model formula, because the model data contain only +those; a model without random effects has one only if a predictor such +as \code{0 + id} names the participants. \code{NULL} (default) +uses the grouping variable of the random effects; when the model has more +than one (e.g. participants and items), name it here.} +} +\value{ +A data frame of class \code{bmm_ppp} with one row per statistic (and per +participant for \code{level = "participant"}) and the columns \code{statistic}, +\code{observed} and \code{replicated} (the means of the statistic over draws, for the +observed and the simulated counts, as reported by Klauer, 2010, and +TreeBUGS), \code{ppp} and \code{ndraws}. The draws of the statistics are stored in +the attribute \code{"draws"}. +} +\description{ +Computes the goodness-of-fit statistics \eqn{T_1} and \eqn{T_2} of Klauer +(2010) and their posterior predictive p-values (PPP) for models that fit +aggregated response counts: \code{\link[=mpt]{mpt()}}, \code{\link[=m3]{m3()}}, \code{\link[=sdt_yn]{sdt_yn()}}, \code{\link[=sdt_rating]{sdt_rating()}}, +\code{\link[=sdt_mafc]{sdt_mafc()}}, \code{\link[=sdt_ranking]{sdt_ranking()}} and \code{\link[=sdt_cdp]{sdt_cdp()}}. TreeBUGS' \code{PPP()} reports the +same two statistics; Details says where their values differ. +} +\details{ +Both statistics are computed once per posterior draw, for the observed counts +and for counts simulated from that draw with \code{\link[brms:posterior_predict]{brms::posterior_predict()}}, +each against the expected counts of the same draw from +\code{\link[brms:posterior_epred]{brms::posterior_epred()}}. The simulated data are new data from the same +participants. The PPP is the proportion of draws in which the simulated +counts deviate more from the expectation than the observed counts do. A PPP +near 0 flags misfit; a PPP away from 0 means that the model reproduces the +aspect of the data the statistic measures. + +\strong{\eqn{T_1}: mean structure.} The rows of the data are summed into cells, +and \eqn{T_1} is Pearson's \eqn{X^2} over all cells and response categories, +\deqn{T_1 = \sum_{c,k} \frac{(n_{ck} - \hat{n}_{ck})^2}{\hat{n}_{ck}},} +where \eqn{n_{ck}} is the observed and \eqn{\hat{n}_{ck}} the expected count +of category \eqn{k} in cell \eqn{c}. Klauer (2010, Eq. 17) and TreeBUGS +average the counts over participants instead of summing them. With the same +number \eqn{T} of participants in every cell, the two give the same PPP and +\eqn{T_1} here is \eqn{T} times theirs; otherwise the mean gives the cells +with fewer participants more weight. + +\strong{\eqn{T_2}: covariance structure.} For every participant, the counts are +summed per cell, giving one vector per participant over all cells and +categories. \eqn{T_2} compares the covariance matrix \eqn{s} of these vectors +across participants with the covariance matrix \eqn{\sigma} the model +expects (Klauer, 2010, Eq. 18): +\deqn{T_2 = \sum_{i,j} \frac{(s_{ij} - \sigma_{ij})^2}{\sqrt{\sigma_{ii}\sigma_{jj}}}.} +\eqn{\sigma} is the covariance of the expected counts across participants +plus the average multinomial covariance within a participant, which is the +expected value of \eqn{s} given the draw. Both covariances divide by +\eqn{T - 1}, where \eqn{T} is the number of participants. Klauer (2010) +multiplies the multinomial part by \eqn{(T - 1) / T}, the expected value +for covariances that divide by \eqn{T}; that version of \eqn{T_2} is +\eqn{(T - 1) / T} times this one and has the same PPP. TreeBUGS combines +the factor with covariances that divide by \eqn{T - 1}, so its \eqn{\sigma} +falls short of the expected \eqn{s} by \eqn{1 / T} of the multinomial part. +The sum runs over +the full matrix, so the covariances between trees or conditions, which +carry the correlations between the person parameters, are part of the +check. When participants differ in the cells they have, as in a +between-participants design, they are grouped by their set of cells, +\eqn{T_2} is computed per group, and the results are summed. A participant +who shares the set of cells with nobody else has no covariance and is left +out, with a message. + +\strong{Cells.} By default, a cell is a combination of the tree (\code{\link[=mpt]{mpt()}}), the +stimulus (the SDT models) or the set size (\code{\link[=sdt_ranking]{sdt_ranking()}}), and the +categorical population-level predictors of the model formula. Numeric +predictors are not part of the default, because every distinct value would +form a cell of its own. Rows of a cell are summed, so the trees or stimuli +must be part of \code{by} for the statistic to compare each tree's own +categories. + +\strong{Impossible categories.} Categories with an expected count of about 0 in +every draw, such as the impossible categories of an \code{\link[=mpt]{mpt()}} or \code{\link[=m3]{m3()}} model +or the ranks beyond the set size of an \code{\link[=sdt_ranking]{sdt_ranking()}} model, are left out +of both statistics. A cell with a single possible category, such as an +\code{\link[=m3]{m3()}} set size without other candidates, is left out of \eqn{T_2}: the +number of trials fixes its count. +} +\examples{ +\dontshow{if (isTRUE(Sys.getenv("BMM_EXAMPLES"))) withAutoprint(\{ # examplesIf} +trees <- list( + mpt_tree("target", list(old = "D + (1 - D) * g", new = "(1 - D) * (1 - g)")), + mpt_tree("lure", list(old = "(1 - D) * g", new = "D + (1 - D) * (1 - g)")) +) +dat <- expand.grid(id = factor(1:20), item_type = c("target", "lure"), + stringsAsFactors = FALSE) +dat$old <- rbinom(nrow(dat), 50, ifelse(dat$item_type == "target", 0.8, 0.2)) +dat$new <- 50 - dat$old + +fit <- bmm( + bmf(D ~ 1 + (1 | id), g ~ 1 + (1 | id)), dat, mpt(trees, "item_type"), + cores = 4, backend = "cmdstanr" +) +ppp(fit) +ppp(fit, level = "participant") +\dontshow{\}) # examplesIf} +} +\references{ +Klauer, K. C. (2010). Hierarchical multinomial processing tree +models: A latent-trait approach. \emph{Psychometrika, 75}(1), 70–98. +https://doi.org/10.1007/s11336-009-9141-0 +} +\seealso{ +\code{\link[=pp_check.bmmfit]{pp_check.bmmfit()}} for the graphical check. +} diff --git a/tests/local/README.md b/tests/local/README.md index 51a87090a..6290f7afb 100644 --- a/tests/local/README.md +++ b/tests/local/README.md @@ -31,9 +31,11 @@ Files: - `helper-mpt.R`: sources `tests/testthat/helper-mpt.R`, so the MPT fits use the same tree and data builders as the unit tests. - `test-*.R`: the tests, among them `test-kfold-update.R` (refits through - `update()` on a slice of the fit's data, grouped `kfold()`) and + `update()` on a slice of the fit's data, grouped `kfold()`), `test-mpt-fits.R` (`mpt()` fits read through `conditional_effects()`, and - the prior-only fit behind the default simplex prior). Keep each fit + the prior-only fit behind the default simplex prior) and + `test-ppp-fits.R` (`ppp()` on well-specified, misspecified and + impossible-category `mpt()` fits). Keep each fit short (one or two chains, a few hundred iterations) and assert with `expect_range()` or a tolerance, never with exact equality to a stored fit. - `make-fixtures.R`: regenerates four fixtures in `tests/testthat/assets/`: diff --git a/tests/local/test-ppp-fits.R b/tests/local/test-ppp-fits.R new file mode 100644 index 000000000..6e74bb8fb --- /dev/null +++ b/tests/local/test-ppp-fits.R @@ -0,0 +1,43 @@ +test_that("ppp() does not flag a well-specified mpt model and flags a misspecified one", { + withr::local_seed(2) + withr::local_options("bmm.silent" = 2) + n_id <- 15 + dat <- expand.grid(id = factor(seq_len(n_id)), item_type = c("old", "new"), + stringsAsFactors = FALSE) + D <- stats::plogis(stats::qlogis(0.6) + stats::rnorm(n_id, 0, 0.5))[dat$id] + p_old <- ifelse(dat$item_type == "old", D + (1 - D) * 0.5, (1 - D) * 0.5) + dat$old <- stats::rbinom(nrow(dat), 60, p_old) + dat$new <- 60L - dat$old + + fit <- bmm( + bmf(D ~ 1 + (1 | id), g ~ 1 + (1 | id)), dat, mpt(mpt_2htm_trees(), "item_type"), + backend = "cmdstanr", chains = 2, iter = 1000, refresh = 0, seed = 2 + ) + res <- ppp(fit) + expect_equal(attr(res, "by"), "item_type") + expect_true(all(res$ppp > 0.05)) + expect_equal(nrow(ppp(fit, level = "participant")), n_id) + + misfit <- bmm( + bmf(D ~ 1, g = 0.9), dat, mpt(mpt_2htm_trees(), "item_type"), + backend = "cmdstanr", chains = 2, iter = 1000, refresh = 0, seed = 2 + ) + withr::local_options("bmm.silent" = 1) + expect_message(res_misfit <- ppp(misfit), "computes T1 only") + expect_lt(res_misfit$ppp, 0.01) +}) + +test_that("ppp() handles the impossible categories of a fitted mpt model", { + withr::local_seed(3) + withr::local_options("bmm.silent" = 2) + dat <- mpt_impossible_data() + dat$corr <- dat$corr + stats::rbinom(nrow(dat), 5, 0.5) + + fit <- bmm( + bmf(Pm ~ 1 + (1 | id), Pb ~ 1), dat, mpt(mpt_impossible_trees(), "tree"), + backend = "cmdstanr", chains = 2, iter = 1000, refresh = 0, seed = 3 + ) + res <- ppp(fit) + expect_equal(attr(res, "by"), "tree") + expect_true(all(is.finite(c(res$observed, res$replicated, res$ppp)))) +}) diff --git a/tests/testthat/test-links.R b/tests/testthat/test-links.R index dd9ce1944..9979d0c07 100644 --- a/tests/testthat/test-links.R +++ b/tests/testthat/test-links.R @@ -8,6 +8,10 @@ test_that("an unrecognized link target is refused, not appended", { links = list(zzz = "log")), "Unrecognized link target" ) + expect_error( + mpt(mpt_2htm_trees(), tree_id = "item_type", links = list(zzz = "probit")), + "Unrecognized link target" + ) }) test_that("a link target one edit from a parameter is read as that parameter", { @@ -23,6 +27,13 @@ test_that("a link target one edit from a parameter is read as that parameter", { "'bond' read as 'bound'" ) expect_equal(model$links$bound, "softplus") + + expect_warning( + model <- mpt(mpt_2htm_trees(), tree_id = "item_type", links = list(gg = "probit")), + "'gg' read as 'g'" + ) + expect_equal(model$links$g, "probit") + expect_false("gg" %in% names(model$links)) }) test_that("a link target that abbreviates a parameter is read as that parameter", { @@ -90,6 +101,25 @@ test_that("a link the model cannot honour is refused at the constructor", { settable_link_functions(m3_ss(NULL)), c("log", "softplus", "logit", "probit", "identity") ) + + # mpt applies its links in the generated category formulas, so it honours + # the two latent scales it has default priors for, as a list or as a string + expect_error( + mpt(mpt_2htm_trees(), tree_id = "item_type", links = list(D = "cauchit")), + "Unknown link function\\(s\\): 'cauchit'\\. mpt\\(\\) takes 'logit', 'probit'" + ) + expect_error( + mpt(mpt_2htm_trees(), tree_id = "item_type", links = "cauchit"), + "Unknown link function\\(s\\): 'cauchit'\\. mpt\\(\\) takes 'logit', 'probit'" + ) + expect_error( + mpt(mpt_2htm_trees(), tree_id = "item_type", links = c("logit", "probit")), + "one link function for the whole model" + ) + expect_equal( + settable_link_functions(mpt(mpt_2htm_trees(), tree_id = "item_type")), + c("logit", "probit") + ) }) test_that("a parameter fixed for scaling is refused with its own message", { @@ -173,8 +203,20 @@ test_that("a link the model does not pass on to the fit is refused", { links = list(kappa = "softplus")), "cannot be changed in mixture3p\\(\\)" ) + # a simplex member and its stick-breaking component take the link of the + # generated stick-breaking formulas, not one of their own + simplex_tree <- mpt_tree("t", list(A = "gA", B = "gB", C = "gC")) + expect_error( + mpt(simplex_tree, simplex = c("gA", "gB", "gC"), links = list(gA = "probit")), + "cannot be changed in mpt\\(\\)" + ) + expect_error( + mpt(simplex_tree, simplex = c("gA", "gB", "gC"), links = list(gAraw = "probit")), + "cannot be changed in mpt\\(\\)" + ) # naming the link the model already uses asks for no change expect_silent(sdm(resp_error = "y", links = list(kappa = "log"))) + expect_silent(mpt(mpt_2htm_trees(), tree_id = "item_type", links = list(D = "logit"))) }) test_that("the links a model applies are exactly the settable ones", { diff --git a/tests/testthat/test-model_mpt.R b/tests/testthat/test-model_mpt.R index 93609582c..06d7dc333 100644 --- a/tests/testthat/test-model_mpt.R +++ b/tests/testthat/test-model_mpt.R @@ -19,6 +19,146 @@ test_that("mpt_tree validates its inputs", { expect_error(mpt_tree("t", list(a = 0.5, b = "x")), "character strings") }) +test_that("mpt_tree takes formulas and strings interchangeably", { + mpt_formula_tree_pairs <- list( + twohtm_old = list( + formulas = mpt_tree("old", + old ~ D + (1 - D) * g, + new ~ (1 - D) * (1 - g) + ), + strings = mpt_tree("old", list( + old = "D + (1 - D) * g", + new = "(1 - D) * (1 - g)" + )) + ), + pcm = list( + formulas = mpt_tree("study", + C ~ cp + (1 - cp) * rp * rp, + E ~ 2 * (1 - cp) * rp * (1 - rp), + U ~ (1 - cp) * (1 - rp) * (1 - rp) + ), + strings = mpt_tree("study", list( + C = "cp + (1 - cp) * rp * rp", + E = "2 * (1 - cp) * rp * (1 - rp)", + U = "(1 - cp) * (1 - rp) * (1 - rp)" + )) + ), + source_monitoring = list( + formulas = mpt_tree("sourceA", + A ~ dA + (1 - dA) * gA, + B ~ (1 - dA) * gB, + New ~ (1 - dA) * gNew + ), + strings = mpt_tree("sourceA", list( + A = "dA + (1 - dA) * gA", + B = "(1 - dA) * gB", + New = "(1 - dA) * gNew" + )) + ), + folded_fractions = list( + formulas = mpt_tree("nodist", + corr ~ Pm * Pb + (1 - Pm) * (1/5), + npl ~ Pm * (1 - Pb) + (1 - Pm) * (4/5), + impossible = "dist" + ), + strings = mpt_tree("nodist", list( + corr = "Pm * Pb + (1 - Pm) * (1/5)", + npl = "Pm * (1 - Pb) + (1 - Pm) * (4/5)" + ), impossible = "dist") + ) + ) + for (pair in mpt_formula_tree_pairs) { + expect_identical(pair$formulas, pair$strings) + } + expect_equal(mpt_formula_tree_pairs$folded_fractions$formulas$branches$corr[[3]][[3]], 0.2) + + formula_model <- mpt( + list( + mpt_tree("old", old ~ D + (1 - D) * g, new ~ (1 - D) * (1 - g)), + mpt_tree("new", old ~ (1 - D) * g, new ~ D + (1 - D) * (1 - g)) + ), + tree_id = "item_type" + ) + string_model <- mpt(mpt_2htm_trees(), tree_id = "item_type") + expect_equal(formula_model, string_model, ignore_attr = "call") + dat <- mpt_2htm_data() + expect_identical( + stancode(bmf(D ~ 1, g ~ 1), data = dat, model = formula_model), + stancode(bmf(D ~ 1, g ~ 1), data = dat, model = string_model) + ) +}) + +test_that("mpt_tree validates formula input", { + expect_error(mpt_tree("t", ~ D), "two-sided") + expect_error(mpt_tree("t", a ~ D, ~ 1 - D), "two-sided") + expect_error(mpt_tree("t", "a" ~ D, b ~ 1 - D), "single response category name") + expect_error(mpt_tree("t", a ~ D, c(a, b) ~ 1 - D), "single response category name") + expect_error(mpt_tree("t", a ~ D, a ~ 1 - D), "unique within a tree") + expect_error(mpt_tree("t", a ~ D, "1 - D"), "mix") + expect_error(mpt_tree("t", a ~ D, 0.5), "mix") + expect_error(mpt_tree("t", list(a = "D", b = "1 - D"), c ~ 1), "Cannot mix") + expect_error(mpt_tree("t", list(a = "D", b = "1 - D"), "a"), "impossible argument has to be named") + expect_error( + mpt_tree("t", a ~ D, b ~ 1 - D, "c"), + "Cannot mix formulas.*impossible argument has to be named, impossible = \"c\"" + ) + expect_error(mpt_tree("t", list(a ~ D, b ~ 1 - D)), "separate arguments") + expect_error(mpt_tree("t", a ~ D, b ~ 1 - D, impossible = "a"), "both impossible") + expect_error(mpt_tree("t", a ~ 1e-5 * D, b ~ 1 - D), "scientific notation") +}) + +test_that("mpt_tree gives the same tree for a formula as for the same string", { + expressions <- c( + "D", "-D", "-(1 - D)", "(D)", "((D))", "D^2", "D * 2^-1", "2 * (1 - D) * g", + "pnorm(-x) + exp(-(y + 1))", "`my par` * (1 - `my par`)", "`a.b` + 0.5", + "1/3 + D * (1/3)", "(1 - D) * (4/5)", "-1 + 2 * D", "D / (1 + D)" + ) + for (expression in expressions) { + from_string <- mpt_tree("t", list(a = expression, b = "1")) + from_formula <- mpt_tree("t", as.formula(paste("a ~", expression)), b ~ 1) + expect_identical(from_formula, from_string, label = expression) + } +}) + +test_that("mpt_tree stores nothing from the environment of a formula", { + formula_with_large_environment <- function() { + large_object <- numeric(1e5) + a ~ D + } + tree <- mpt_tree("t", formula_with_large_environment(), b ~ 1 - D) + string_tree <- mpt_tree("t", list(a = "D", b = "1 - D")) + expect_identical(length(serialize(tree, NULL)), length(serialize(string_tree, NULL))) +}) + +test_that("mpt_tree does not let a partially matched name reorder or drop formulas", { + expect_error(mpt_tree("t", a ~ D, b = b ~ 1 - D), "in full") + expect_error(mpt_tree("t", a ~ D, foo = b ~ 1 - D), "in full") + expect_error(mpt_tree("t", list(a = "D"), imp = "b"), "in full") + expect_identical( + mpt_tree("t", br = list(a = "D", b = "1 - D")), + mpt_tree("t", list(a = "D", b = "1 - D")) + ) + expect_identical( + names(mpt_tree(name = "t", a ~ D, b ~ 1 - D)$branches), + c("a", "b") + ) +}) + +test_that("mpt_tree rejects branch values that are not expressions", { + expect_error(mpt_tree("t", a ~ D, b ~ NULL), "must be an expression") + expect_error(mpt_tree("t", a ~ D, b ~ "1 - D"), "must be an expression") + expect_error(mpt_tree("t", a ~ D, eval(bquote(b ~ .(c(0.2, 0.3))))), "must be an expression") + expect_error(mpt_tree("t", list(a = "D", b = "NULL")), "must be an expression") + expect_error(mpt_tree("t", list(a = "D", b = "'1 - D'")), "must be an expression") + expect_s3_class(mpt_tree("t", a ~ 0.5, b ~ 0.5), "mpt_tree") +}) + +test_that("mpt_tree names a trailing comma instead of reporting a missing argument", { + expect_error(mpt_tree("t", a ~ D, b ~ 1 - D, ), "trailing comma") + expect_error(mpt_tree("t", a ~ D, , b ~ 1 - D), "trailing comma") + expect_error(mpt_tree("t", list(a = "D", b = "1 - D"), ), "trailing comma") +}) + test_that("mpt stores its derived state once", { model <- mpt(mpt_impossible_trees(), tree_id = "tree") expect_equal(model$other_vars$link, "logit") @@ -121,6 +261,16 @@ test_that("mpt errors on inconsistent trees and a missing tree_id", { }) test_that("mpt refuses a zero branch and points to impossible categories", { + expect_error( + mpt(mpt_tree("a", list(y = "D", z = "0 * D", w = "1 - D"))), + "is 0 at the test values.*mpt_tree\\(impossible = \\)" + ) + # the hint reads as one sentence, not with the source's line breaks + expect_error( + mpt(mpt_tree("a", list(y = "D", z = "0 * D", w = "1 - D"))), + "undefined. Declare a response category that a tree cannot produce with mpt_tree(impossible = ) instead", + fixed = TRUE + ) zero_padded <- list( mpt_tree("a", list(x = "D", y = "1 - D", z = "0")), mpt_tree("b", list(x = "1 - 1", y = "g", z = "1 - g")) @@ -207,6 +357,86 @@ test_that("check_formula uses Phi for the probit link", { expect_false(grepl("pnorm", old_rhs, fixed = TRUE)) }) +test_that("links = list() sets one parameter's link and matched prior at construction", { + model <- mpt(mpt_2htm_trees(), tree_id = "item_type", links = list(D = "probit")) + expect_equal(model$links, list(D = "probit", g = "logit")) + expect_equal(model$default_priors$D, bmm:::.mpt_latent_prior("probit")) + expect_equal(model$default_priors$g, bmm:::.mpt_latent_prior("logit")) + expect_output(print(model), "Links:\\s+D = probit; g = logit") + + pr <- suppressMessages(default_prior(bmf(D ~ 1, g ~ 1), mpt_2htm_data(), model)) + intercepts <- pr$class == "b" & pr$coef == "Intercept" + expect_equal(pr$prior[intercepts & pr$nlpar == "D"], "normal(0, 1)") + expect_equal(pr$prior[intercepts & pr$nlpar == "g"], "logistic(0, 1)") + + # the single string remains the shorthand for every latent parameter + all_probit <- mpt(mpt_2htm_trees(), tree_id = "item_type", links = "probit") + expect_equal(all_probit$links, list(D = "probit", g = "probit")) + expect_equal(all_probit$default_priors$g, bmm:::.mpt_latent_prior("probit")) + expect_equal(all_probit$other_vars$link, "probit") +}) + +test_that("the links argument of mpt() refuses what set_links() refuses", { + trees <- mpt_2htm_trees() + expect_warning( + model <- mpt(trees, tree_id = "item_type", links = list(gg = "probit")), + "'gg' read as 'g'" + ) + expect_equal(model$links$g, "probit") + expect_error( + mpt(trees, tree_id = "item_type", links = list(zzz = "probit")), + "Unrecognized link target\\(s\\): 'zzz'" + ) + expect_error( + mpt(trees, tree_id = "item_type", links = list(D = "cauchit")), + "Unknown link function\\(s\\): 'cauchit'\\. mpt\\(\\) takes 'logit', 'probit'" + ) + expect_error( + mpt(trees, tree_id = "item_type", links = "cauchit"), + "Unknown link function\\(s\\): 'cauchit'\\. mpt\\(\\) takes 'logit', 'probit'" + ) + expect_error( + mpt(trees, tree_id = "item_type", links = c("logit", "probit")), + "one link function for the whole model" + ) + expect_error(mpt(trees, tree_id = "item_type", links = 1), "named list") + # NULL takes the default link, as in the other model constructors + expect_equal( + mpt(trees, tree_id = "item_type", links = NULL)$links, + list(D = "logit", g = "logit") + ) + + # simplex members and their stick-breaking components carry identity links + tree <- mpt_tree("t", list(A = "gA", B = "gB", C = "gC")) + expect_error( + mpt(tree, simplex = c("gA", "gB", "gC"), links = list(gA = "probit")), + "The link of 'gA' cannot be changed in mpt\\(\\)" + ) + expect_error( + mpt(tree, simplex = c("gA", "gB", "gC"), links = list(gAraw = "probit")), + "The link of 'gAraw' cannot be changed in mpt\\(\\)" + ) +}) + +test_that("mpt() refuses a link for a parameter the restrictions removed", { + trees <- list( + mpt_tree("old", list(old = "Do + (1 - Do) * g", new = "(1 - Do) * (1 - g)")), + mpt_tree("new", list(old = "(1 - Dn) * g", new = "Dn + (1 - Dn) * (1 - g)")) + ) + expect_error( + mpt(trees, tree_id = "item_type", restrictions = "Dn = Do", links = list(Dn = "probit")), + "'Dn' .*removed by the restrictions" + ) + # 'g' is one edit from 'D', so the link matcher used to read it as D + expect_error( + mpt(mpt_2htm_trees(), tree_id = "item_type", restrictions = "g = 0.5", + links = list(g = "probit")), + "'g' .*removed by the restrictions" + ) + model <- mpt(trees, tree_id = "item_type", restrictions = "Dn = Do", links = list(Do = "probit")) + expect_equal(model$links, list(Do = "probit", g = "logit")) +}) + test_that("formulas for response categories are rejected", { model <- mpt(mpt_2htm_trees(), tree_id = "item_type") dat <- mpt_2htm_data() @@ -546,7 +776,10 @@ test_that("printing an mpt model lists trees, restrictions and the identifiabili simplex = c("gA", "gB", "gC") ) expect_output(print(simplex), "2 free parameter\\(s\\), 2 degrees of freedom") - expect_output(print(simplex), "'gAraw', 'gBraw' pass through inv_logit\\(\\)") + expect_output( + print(simplex), + "pass through inv_logit\\(\\) in the generated formulas \\(hence the identity links above\\)" + ) expect_output(print(simplex), "symmetric Dirichlet\\(1\\), prior mean 1/3 = 0.333") expect_match( parameter_info(simplex)$description[4], @@ -1245,6 +1478,37 @@ test_that("stick-breaking members form a simplex for any raw values", { ) }) +test_that("generated simplex formulas are neither announced as missing nor listed in the summary", { + tree <- mpt_tree("t", list(A = "gA", B = "gB", New = "gNew")) + model <- mpt(tree, simplex = c("gA", "gB", "gNew")) + dat <- data.frame(A = 10L, B = 10L, New = 10L) + expect_no_message( + fit <- bmm( + bmf(gA ~ 1, gB ~ 1), dat, model, backend = "mock", mock_fit = 1, rename = FALSE + ), + message = "No formula" + ) + expect_equal(summarise_formula(fit$bmm$user_formula, fit$bmm$model), c("gA ~ 1", "gB ~ 1")) + # the header block print(summary(fit)) uses; a mock fit has no draws to summarise + withr::local_options(bmm.color_summary = FALSE) + expect_output( + bmm:::print_model_header(fit$bmm$model, fit$bmm$user_formula, "dat", 1L), + "Formula: gA ~ 1\\s+gB ~ 1\\s+Data:" + ) + + # a parameter the user leaves out of a plain model is still announced, and + # the summary still lists the constants the user fixed + plain <- mpt(mpt_2htm_trees(), tree_id = "item_type") + expect_message( + bmm(bmf(D ~ 1), mpt_2htm_data(), plain, backend = "mock", mock_fit = 1, rename = FALSE), + "No formula for parameter g" + ) + fixed <- suppressMessages(bmm( + bmf(D ~ 1, g = 0.5), mpt_2htm_data(), plain, backend = "mock", mock_fit = 1, rename = FALSE + )) + expect_equal(summarise_formula(fixed$bmm$user_formula, fixed$bmm$model), c("D ~ 1", "g = 0.5")) +}) + test_that("predictors on the derived simplex parameter are rejected", { trees <- list(mpt_tree("t", list(A = "gA", B = "gB", New = "gNew"))) model <- mpt(trees, simplex = c("gA", "gB", "gNew")) @@ -1505,7 +1769,16 @@ test_that("check_data builds possibility indicators for impossible categories", dat_observed$dist[dat_observed$cond == "same"][1] <- 3L expect_error( check_data(model, dat_observed, bmf(Pm ~ 1, Pb ~ 1)), - "declared impossible for 1 observation" + "declared impossible for 1 row\\(s\\)" + ) +}) + +test_that("check_data warns about a tree without rows in units of rows", { + model <- mpt(mpt_2htm_trees(), tree_id = "item_type") + dat <- mpt_2htm_data() + expect_warning( + check_data(model, dat[dat$item_type == "old", ], bmf(D ~ 1, g ~ 1)), + "no rows for tree\\(s\\): 'new'" ) }) diff --git a/tests/testthat/test-model_mpt_data.R b/tests/testthat/test-model_mpt_data.R new file mode 100644 index 000000000..af833ce97 --- /dev/null +++ b/tests/testthat/test-model_mpt_data.R @@ -0,0 +1,325 @@ +eqn_2htm_model <- function(...) { + eqn_file <- withr::local_tempfile(fileext = ".eqn", .local_envir = parent.frame()) + writeLines(c( + "6", + "old hit D_o", + "old hit (1-D_o)*g", + "old miss (1-D_o)*(1-g)", + "new fa (1-D_n)*g", + "new cr D_n", + "new cr (1-D_n)*(1-g)" + ), eqn_file) + suppressMessages(mpt_from_eqn( + eqn_file, categories = c(hit = "yes", fa = "yes", miss = "no", cr = "no"), + tree_id = "item_type", restrictions = "D_n = D_o", ... + )) +} + +wide_2htm <- function() { + data.frame(hit = c(30, 25), miss = c(10, 15), fa = c(8, 12), cr = c(32, 28)) +} + +test_that("mpt_data_reshape gives the same long data for TreeBUGS and MPTinR columns", { + model <- eqn_2htm_model() + by_name <- mpt_data_reshape(wide_2htm(), model) + expect_equal( + by_name, + data.frame( + id = c(1L, 1L, 2L, 2L), item_type = c("old", "new", "old", "new"), + yes = c(30L, 8L, 25L, 12L), no = c(10L, 32L, 15L, 28L) + ) + ) + # MPTinR sorts the trees (new, old) and then the labels (cr, fa, hit, miss) + mptinr <- unname(as.matrix(wide_2htm()[c("cr", "fa", "hit", "miss")])) + expect_identical(mpt_data_reshape(mptinr, model, columns = "mptinr"), by_name) + expect_no_error(suppressMessages( + bmm_data_check(bmf(Do ~ 1, g ~ 1), by_name, model) + )) +}) + +test_that("mpt_data_reshape sorts MPTinR columns like check.mpt()", { + eqn_text <- function(first_label) { + c( + "8", + "10 12 a", "10 2 (1-a)*b", "10 1 (1-a)*(1-b)", + paste("9", first_label, "a"), "9 4 1-a", + "2 3 a", "2 11 (1-a)*b", "2 20 (1-a)*(1-b)" + ) + } + build_model <- function(first_label) { + eqn_file <- withr::local_tempfile(fileext = ".eqn", .local_envir = parent.frame()) + writeLines(eqn_text(first_label), eqn_file) + suppressMessages(mpt_from_eqn( + eqn_file, + categories = setNames( + c("x", "y", "z", "x", "y", "x", "y", "z"), + c("12", "2", "1", first_label, "4", "3", "11", "20") + ), + tree_names = c("10" = "t10", "9" = "t9", "2" = "t2"), + impossible = list(t9 = "z"), tree_id = "tree" + )) + } + # one distinct count per label, so that a wrong order changes the long data + counts <- function(labels) { + as.data.frame( + as.list(setNames(seq_along(labels) * 10, labels)), check.names = FALSE + ) + } + + # MPTinR 1.14.1, check.mpt(): numbers sort numerically if all labels of the + # file are numbers (3 11 20 | 4 5 | 1 2 12) ... + numeric_model <- build_model("5") + numeric_order <- c("3", "11", "20", "4", "5", "1", "2", "12") + expect_identical( + mpt_data_reshape( + unname(as.matrix(counts(numeric_order))), numeric_model, columns = "mptinr" + ), + mpt_data_reshape(counts(numeric_order), numeric_model) + ) + # ... and as character strings if one label is not (11 20 3 | 4 a | 1 12 2) + mixed_model <- build_model("a") + mixed_order <- c("11", "20", "3", "4", "a", "1", "12", "2") + expect_identical( + mpt_data_reshape( + unname(as.matrix(counts(mixed_order))), mixed_model, columns = "mptinr" + ), + mpt_data_reshape(counts(mixed_order), mixed_model) + ) +}) + +test_that("mpt_data_reshape sets the categories a tree cannot produce to 0", { + model <- mpt( + list( + mpt_tree("pairs", list(E1 = "c * r", E2 = "1 - c * r"), impossible = "F1"), + mpt_tree("singles", list(F1 = "u", E2 = "1 - u"), impossible = "E1") + ), + tree_id = "item_type" + ) + wide <- data.frame( + pairs.E1 = c(20, 22), pairs.E2 = c(20, 18), + singles.F1 = c(25, 24), singles.E2 = c(15, 16) + ) + long <- mpt_data_reshape(wide, model) + expect_equal( + long, + data.frame( + id = c(1L, 1L, 2L, 2L), item_type = rep(c("pairs", "singles"), 2), + E1 = c(20L, 0L, 22L, 0L), E2 = c(20L, 15L, 18L, 16L), + F1 = c(0L, 25L, 0L, 24L) + ) + ) + expect_error( + mpt_data_reshape(unname(as.matrix(wide)), model, columns = "mptinr"), + "needs a model from mpt_from_eqn" + ) + expect_no_error(suppressMessages( + bmm_data_check(bmf(c ~ 1, r ~ 1, u ~ 1), long, model) + )) +}) + +test_that("mpt_data_reshape fills the impossible categories of an EQN model with 0", { + eqn_file <- withr::local_tempfile(fileext = ".eqn") + writeLines(c( + "P P_E1 c*r", + "P P_E2 (1-c)*u*u", + "P P_E3 2*(1-c)*u*(1-u)", + "P P_E4 c*(1-r)+(1-c)*(1-u)*(1-u)", + "S F1 u", + "S F2 1-u" + ), eqn_file) + model <- suppressMessages(mpt_from_eqn( + eqn_file, tree_id = "tree", + impossible = list(P = c("F1", "F2"), S = c("E1", "E2", "E3", "E4")) + )) + # the labels keep the file's prefix; the long data use the bmm names + wide <- data.frame( + P_E1 = c(5, 1), P_E2 = c(3, 2), P_E3 = c(4, 3), P_E4 = c(8, 4), + F1 = c(6, 5), F2 = c(4, 6) + ) + long <- mpt_data_reshape(wide, model) + expect_equal( + long, + data.frame( + id = c(1L, 1L, 2L, 2L), tree = c("P", "S", "P", "S"), + E1 = c(5L, 0L, 1L, 0L), E2 = c(3L, 0L, 2L, 0L), + E3 = c(4L, 0L, 3L, 0L), E4 = c(8L, 0L, 4L, 0L), + F1 = c(0L, 6L, 0L, 5L), F2 = c(0L, 4L, 0L, 6L) + ) + ) + # MPTinR sorts the trees (P, S) and the labels within a tree + expect_identical( + mpt_data_reshape(unname(as.matrix(wide)), model, columns = "mptinr"), long + ) + expect_no_error(suppressMessages( + bmm_data_check(bmf(c ~ 1, r ~ 1, u ~ 1), long, model) + )) +}) + +test_that("mpt_data_reshape carries non-count columns and uses the id column", { + model <- eqn_2htm_model() + wide <- cbind( + subject = c("s1", "s2"), group = c("a", "b"), age = c(21, 34), wide_2htm() + ) + long <- mpt_data_reshape(wide, model, id = "subject") + expect_named(long, c("subject", "item_type", "group", "age", "yes", "no")) + expect_equal(long$group, c("a", "a", "b", "b")) + expect_equal(long$age, c(21, 21, 34, 34)) + expect_equal(long$subject, c("s1", "s1", "s2", "s2")) +}) + +test_that("mpt_data_reshape reads a label shared by trees from the qualified name", { + eqn_file <- withr::local_tempfile(fileext = ".eqn") + writeLines(c("4", "old 1 D", "old 2 1-D", "new 1 g", "new 2 1-g"), eqn_file) + model <- suppressMessages(mpt_from_eqn( + eqn_file, categories = c("1" = "yes", "2" = "no"), tree_id = "tree" + )) + wide <- data.frame(old.1 = 30, old.2 = 10, new.1 = 8, new.2 = 32) + expect_equal(mpt_data_reshape(wide, model)$yes, c(30L, 8L)) + expect_error( + mpt_data_reshape(data.frame(`1` = 1, `2` = 2, check.names = FALSE), model), + "old.1" + ) +}) + +test_that("mpt_data_reshape stops when one wide column would feed two categories", { + eqn_file <- withr::local_tempfile(fileext = ".eqn") + writeLines(c("4", "old hit p", "old miss 1-p", "new old.hit q", "new old.miss 1-q"), eqn_file) + model <- suppressMessages(mpt_from_eqn( + eqn_file, + categories = c(hit = "yes", miss = "no", old.hit = "yes", old.miss = "no"), + tree_id = "tree" + )) + # 'old.hit' is the label of tree new and the qualified name of hit in tree old + wide <- data.frame(hit = 1, miss = 2, old.hit = 30, old.miss = 40) + expect_error(mpt_data_reshape(wide, model), "more than one response category.*old.hit") + qualified <- data.frame(old.hit = 1, old.miss = 2, new.old.hit = 30, new.old.miss = 40) + expect_equal(mpt_data_reshape(qualified, model)$yes, c(1L, 30L)) +}) + +test_that("mpt_data_reshape sums the labels that categories merged", { + eqn_file <- withr::local_tempfile(fileext = ".eqn") + writeLines(c( + "5", "old hit D", "old miss (1-D)*g", "old miss2 (1-D)*(1-g)", + "new fa g", "new cr 1-g" + ), eqn_file) + model <- suppressMessages(mpt_from_eqn( + eqn_file, + categories = c(hit = "yes", fa = "yes", miss = "no", miss2 = "no", cr = "no"), + tree_id = "tree" + )) + wide <- data.frame(hit = 30, miss = 4, miss2 = 6, fa = 8, cr = 32) + expect_equal(mpt_data_reshape(wide, model)$no, c(10L, 32L)) +}) + +test_that("mpt_data_reshape returns integer counts and keeps the types of carried columns", { + model <- eqn_2htm_model() + wide <- data.frame( + subject = factor(c("s1", "s2")), when = as.Date("2026-01-01") + 0:1, + flag = c(TRUE, NA), group = factor(c("a", "b"), levels = c("b", "a", "unused")), + wide_2htm() + ) + long <- mpt_data_reshape(wide, model, id = "subject") + expect_type(long$yes, "integer") + expect_type(long$no, "integer") + expect_identical(long$subject, factor(c("s1", "s1", "s2", "s2"))) + expect_identical(long$when, rep(wide$when, each = 2)) + expect_identical(long$flag, c(TRUE, TRUE, NA, NA)) + expect_identical(long$group, factor(c("a", "a", "b", "b"), levels = levels(wide$group))) +}) + +test_that("mpt_data_reshape takes the qualified name where a bare label also exists", { + model <- eqn_2htm_model() + wide <- cbind(wide_2htm(), old.hit = c(5, 6)) + long <- mpt_data_reshape(wide, model) + expect_equal(long$yes[long$item_type == "old"], c(5L, 6L)) +}) + +test_that("mpt_data_reshape without a tree_id column in a single-tree model", { + model <- mpt(mpt_tree("t", list(a = "p", b = "1 - p"))) + long <- mpt_data_reshape(data.frame(t.a = c(3, 4), t.b = c(7, 6)), model) + expect_equal(long, data.frame(id = 1:2, a = c(3L, 4L), b = c(7L, 6L))) +}) + +test_that("mpt_data_reshape reports unusable input", { + model <- eqn_2htm_model() + expect_error(mpt_data_reshape(wide_2htm()), "model") + expect_error(mpt_data_reshape(model = model), "data") + expect_error(mpt_data_reshape(wide_2htm(), "model"), "mpt model") + expect_error(mpt_data_reshape(list(1), model), "data frame or a matrix") + expect_error(mpt_data_reshape(wide_2htm()[0, ], model), "no rows") + expect_error(mpt_data_reshape(wide_2htm(), model, columns = "other"), "should be one of") + + expect_error(mpt_data_reshape(wide_2htm(), model, id = "subject"), "id column 'subject'") + expect_error(mpt_data_reshape(wide_2htm(), model, id = c("a", "b")), "id argument") + repeated <- cbind(subject = c("s1", "s1"), wide_2htm()) + expect_error(mpt_data_reshape(repeated, model, id = "subject"), "repeated values") +}) + +test_that("mpt_data_reshape refuses duplicated column names", { + model <- eqn_2htm_model() + repeated <- cbind(wide_2htm(), hit = c(1, 1)) + expect_error(mpt_data_reshape(repeated, model), "duplicated column names.*hit") + positional <- unname(as.matrix(wide_2htm())) + colnames(positional) <- c("a", "a", "b", "b") + expect_error( + mpt_data_reshape(positional, model, columns = "mptinr"), + "duplicated column names.*a.*b" + ) +}) + +test_that("mpt_data_reshape names every missing count column", { + model <- eqn_2htm_model() + err <- expect_error( + mpt_data_reshape(wide_2htm()[c("hit", "fa")], model), "lack count columns" + ) + expect_match(conditionMessage(err), "old.miss or miss", fixed = TRUE) + expect_match(conditionMessage(err), "new.cr or cr", fixed = TRUE) +}) + +test_that("mpt_data_reshape rejects counts that are not counts", { + model <- eqn_2htm_model() + with_hit <- function(value) { + wide <- wide_2htm() + wide$hit[1] <- value + wide + } + expect_error(mpt_data_reshape(with_hit(-1), model), "non-negative whole numbers.*hit") + expect_error(mpt_data_reshape(with_hit(2.5), model), "non-negative whole numbers.*hit") + expect_error(mpt_data_reshape(with_hit(Inf), model), "non-negative whole numbers.*hit") + expect_error(mpt_data_reshape(with_hit(NA), model), "missing values.*hit") + wide <- wide_2htm() + wide$hit <- as.character(wide$hit) + expect_error(mpt_data_reshape(wide, model), "must be numeric.*hit") +}) + +test_that("mpt_data_reshape checks the column count in MPTinR mode", { + model <- eqn_2htm_model() + expect_error( + mpt_data_reshape(matrix(1:6, 2), model, columns = "mptinr"), + "expects 4 columns.*have 3" + ) + expect_error( + mpt_data_reshape( + cbind(subject = 1:2, group = 1:2, unname(as.matrix(wide_2htm()))), + model, id = "subject", columns = "mptinr" + ), + "expects 4 columns.*have 5" + ) +}) + +test_that("mpt_data_reshape refuses carried columns that the long data need", { + model <- eqn_2htm_model() + expect_error( + mpt_data_reshape(cbind(id = 1:2, wide_2htm()), model), + "id = 'name'" + ) + expect_error( + mpt_data_reshape(cbind(item_type = "x", wide_2htm()), model), + "item_type" + ) + expect_error( + mpt_data_reshape(cbind(yes = 1, wide_2htm()), model), + "yes" + ) + expect_no_error(mpt_data_reshape(cbind(id = 1:2, wide_2htm()), model, id = "id")) +}) diff --git a/tests/testthat/test-model_mpt_import.R b/tests/testthat/test-model_mpt_import.R index 580d16e1e..128ff3977 100644 --- a/tests/testthat/test-model_mpt_import.R +++ b/tests/testthat/test-model_mpt_import.R @@ -35,6 +35,13 @@ test_that("mpt_from_string parses MPTinR-style model definitions", { ) }) +test_that("the importers take restrictions right after the model source", { + expect_equal(names(formals(mpt_from_eqn))[1:2], c("file", "restrictions")) + expect_equal(names(formals(mpt_from_string))[1:2], c("text", "restrictions")) + # tree_names stays required behind the optional restrictions + expect_error(mpt_from_string("D # a\n1 - D # b"), "tree_names") +}) + test_that("mpt_from_string says when it sums lines that share a category", { expect_message( mpt_from_string( @@ -252,6 +259,46 @@ test_that("mpt_from_eqn errors when sanitizing merges categories across trees", ) }) +test_that("mapped values that sanitize alike across trees are one category", { + model <- suppressMessages(mpt_from_eqn( + write_eqn(eqn_2htm_lines), + categories = c(hit = "y_es", fa = "yes", miss = "no", cr = "no"), + tree_id = "item_type" + )) + expect_setequal(model$resp_vars$resp_cats, c("yes", "no")) + renaming <- attr(model, "mpt_renaming") + expect_equal( + renaming$bmm_name[renaming$file_name %in% c("hit", "fa")], c("yes", "yes") + ) + expect_equal( + deparse1(model$other_vars$trees$old$branches$yes), "(Do) + ((1 - Do) * g)" + ) + + # a mapped value also meets the file label of another tree + with_label <- suppressMessages(mpt_from_eqn( + write_eqn(c("old hit D", "old no 1-D", "new yes 1-D", "new no D")), + categories = c(hit = "y_es"), tree_id = "item_type" + )) + expect_setequal(with_label$resp_vars$resp_cats, c("yes", "no")) +}) + +test_that("mapped values that sanitize alike within one tree are a clash", { + eqn_file <- write_eqn(c("t a D", "t b (1-D)*g", "t c (1-D)*(1-g)")) + clash <- tryCatch( + mpt_from_eqn(eqn_file, categories = c(a = "y_es", b = "yes", c = "no")), + error = conditionMessage + ) + expect_match(clash, "produces duplicated names: 'a', 'b' become 'yes'") + expect_match(clash, "map them with the categories argument") + # written identically, the two values merge instead + expect_match( + paste(testthat::capture_messages( + mpt_from_eqn(eqn_file, categories = c(a = "yes", b = "yes", c = "no")) + ), collapse = ""), + "tree t: 'a', 'b' -> yes" + ) +}) + test_that("mpt_from_eqn errors when a sanitized parameter name equals a covariate", { clash_file <- write_eqn(c("t yes a_b", "t no 1-a_b", "u yes ab*g", "u no 1-ab*g")) expect_error( diff --git a/tests/testthat/test-pp_check.R b/tests/testthat/test-pp_check.R index 0b162fed2..46adfff47 100644 --- a/tests/testthat/test-pp_check.R +++ b/tests/testthat/test-pp_check.R @@ -512,3 +512,66 @@ test_that(".resolve_pp_conditions() keeps a user column named Poss_*", { expect_true("Poss_block" %in% .resolve_pp_conditions(mock_poss_fit("Poss_block"))) expect_true("Poss_load" %in% .resolve_pp_conditions(mock_sdm_poss_fit())) }) + + +# The tree column of an mpt fit is a valid group even when no formula uses it + +# no parameter formula names item_type, so brms keeps only the Idx_ columns +mock_2htm_fit <- function() { + suppressMessages(bmm( + bmf(D ~ 1, g ~ 1), mpt_2htm_data(), + mpt(mpt_2htm_trees(), tree_id = "item_type"), + backend = "mock", mock_fit = 1, rename = FALSE + )) +} + +# a mock fit has no draws; the check only needs replicates of the right shape +local_mock_posterior_predict <- function(env = parent.frame()) { + local_mocked_bindings( + posterior_predict = function(object, ndraws = 100L, ...) { + n_obs <- nrow(object$data) + cats <- colnames(object$data$Y) + array(stats::rbinom(ndraws * n_obs * length(cats), 50L, 0.5), + c(ndraws, n_obs, length(cats)), dimnames = list(NULL, NULL, cats)) + }, + .package = "brms", .env = env + ) +} + +test_that(".pp_check_restore_tree() rebuilds the tree column from the Idx_ columns", { + fit <- mock_2htm_fit() + expect_false("item_type" %in% names(fit$data)) + restored <- .pp_check_restore_tree(fit, "item_type") + expect_identical(levels(restored$data$item_type), c("old", "new")) + expect_identical(restored$data$item_type == "old", restored$data$Idx_old == 1) + expect_identical(.pp_check_restore_tree(fit, "id"), fit) +}) + +test_that("pp_check() facets an mpt fit by its tree column", { + fit <- mock_2htm_fit() + local_mock_posterior_predict() + expect_no_warning(p <- pp_check(fit, group = "item_type", ndraws = 5)) + expect_s3_class(p$facet, "FacetWrap") + expect_setequal(as.character(p$data$item_type), c("old", "new")) + expect_equal(nrow(p$data), 4L) + old_tree <- fit$data$Idx_old == 1 + expect_equal( + p$data$proportion[p$data$item_type == "old" & p$data$category == "old"], + sum(fit$data$Y[old_tree, "old"]) / sum(fit$data$nTrials[old_tree]) + ) +}) + +test_that("pp_check() pools the trees of an mpt fit by default", { + fit <- mock_2htm_fit() + local_mock_posterior_predict() + p <- pp_check(fit, ndraws = 5) + expect_s3_class(p$facet, "FacetNull") + expect_equal(nrow(p$data), 2L) +}) + +test_that("pp_check() still warns for a column that is neither predictor nor tree", { + fit <- mock_2htm_fit() + local_mock_posterior_predict() + expect_warning(p <- pp_check(fit, group = "id", ndraws = 5), "not a predictor") + expect_s3_class(p$facet, "FacetNull") +}) diff --git a/tests/testthat/test-ppp.R b/tests/testthat/test-ppp.R new file mode 100644 index 000000000..7855d32f2 --- /dev/null +++ b/tests/testthat/test-ppp.R @@ -0,0 +1,279 @@ +# A bmmfit stand-in for ppp() with fixed draws: .ppp_counts() is mocked, so +# only the model class, the formula and the data are read from the fit +ppp_toy_fit <- function(data, formula = bmf(D ~ 1, g ~ 1)) { + structure( + list(bmm = list(model = mpt(mpt_2htm_trees(), "item_type"), user_formula = formula), + data = data), + class = c("bmmfit", "brmsfit") + ) +} + +# draws given as a list of rows x categories matrices, one per draw +ppp_toy_draws <- function(...) { + aperm(simplify2array(list(...)), c(3L, 1L, 2L)) +} + +local_ppp_counts <- function(observed, expected, replicated, data, + env = parent.frame()) { + counts <- list(observed = observed, expected = expected, + replicated = replicated, data = data) + local_mocked_bindings(.ppp_counts = function(model, fit, draw_ids) counts, + .env = env) +} + +# Klauer (2010), Eq. 18, written out: cells is a list of participants x +# categories matrices of expected counts, one per cell, with trials n per +# participant and cell +klauer_t2 <- function(observed, cells, n) { + n_people <- nrow(observed) + blocks <- lapply(seq_along(cells), function(j) { + Reduce(`+`, lapply(seq_len(n_people), function(t) { + e <- cells[[j]][t, ] + diag(e, length(e)) - tcrossprod(e) / n[t, j] + })) / n_people + }) + within <- matrix(0, ncol(observed), ncol(observed)) + offset <- 0 + for (b in blocks) { + idx <- offset + seq_len(ncol(b)) + within[idx, idx] <- b + offset <- offset + ncol(b) + } + # under cov()'s divisor T - 1 the expected sample covariance is + # cov(expected) + within; Klauer's factor (T - 1) / T goes with the divisor T + sigma <- stats::cov(do.call(cbind, cells)) + within + sum((stats::cov(observed) - sigma)^2 / sqrt(outer(diag(sigma), diag(sigma)))) +} + +test_that("ppp() computes T1 by hand on a two-row, two-category example", { + data <- data.frame(row = c("a", "b")) + local_ppp_counts( + observed = matrix(c(3, 6, 7, 4), 2, dimnames = list(NULL, c("old", "new"))), + expected = ppp_toy_draws(matrix(c(4, 5, 6, 5), 2), matrix(c(6, 5, 4, 5), 2)), + replicated = ppp_toy_draws(matrix(c(5, 5, 5, 5), 2), matrix(c(3, 7, 7, 3), 2)), + data = data + ) + fit <- ppp_toy_fit(data) + + # pooled: observed (9, 11); expected (9, 11) and (11, 9); replicated (10, 10) + pooled <- ppp(fit, statistic = "T1", by = character(0)) + t1_obs <- c(0, 2^2 / 11 + 2^2 / 9) + t1_rep <- c(1 / 9 + 1 / 11, 1 / 11 + 1 / 9) + expect_equal(attr(pooled, "draws")$T1[, "observed"], t1_obs) + expect_equal(attr(pooled, "draws")$T1[, "replicated"], t1_rep) + expect_equal(pooled$observed, mean(t1_obs)) + expect_equal(pooled$replicated, mean(t1_rep)) + expect_equal(pooled$ppp, 0.5) + expect_equal(pooled$ndraws, 2L) + + # one cell per row + per_row <- ppp(fit, statistic = "T1", by = "row") + expect_equal( + attr(per_row, "draws")$T1[, "observed"], + c(1 / 4 + 1 / 6 + 1 / 5 + 1 / 5, 3^2 / 6 + 3^2 / 4 + 1 / 5 + 1 / 5) + ) +}) + +test_that("ppp() computes T2 as in Klauer (2010), across cells", { + # three participants with two cells each; trials 10 per cell + data <- data.frame(id = factor(rep(1:3, 2)), item_type = rep(c("old", "new"), each = 3)) + observed <- matrix(c(7, 5, 9, 2, 4, 3, 3, 5, 1, 8, 6, 7), 6) + expected <- matrix(c(6, 6, 8, 3, 3, 2, 4, 4, 2, 7, 7, 8), 6) + replicated <- matrix(c(5, 7, 8, 3, 2, 3, 5, 3, 2, 7, 8, 7), 6) + local_ppp_counts(observed, ppp_toy_draws(expected), ppp_toy_draws(replicated), data) + fit <- ppp_toy_fit(data, bmf(D ~ 1 + (1 | id), g ~ 1 + (1 | id))) + + res <- ppp(fit, statistic = "T2", by = "item_type") + # cells are sorted: "new" (rows 4:6) before "old" (rows 1:3) + wide <- function(m) cbind(m[4:6, ], m[1:3, ]) + cells <- list(expected[4:6, ], expected[1:3, ]) + n <- matrix(10, 3, 2) + expect_equal(unname(attr(res, "draws")$T2[, "observed"]), klauer_t2(wide(observed), cells, n)) + expect_equal(unname(attr(res, "draws")$T2[, "replicated"]), klauer_t2(wide(replicated), cells, n)) +}) + +test_that("T2 is 0 when the observed covariance is the one the model implies", { + # E[cov(counts)] by enumeration: 3 participants, 3 binomial trials each + n_trials <- 3 + p <- c(0.2, 0.5, 0.7) + outcomes <- as.matrix(expand.grid(0:n_trials, 0:n_trials, 0:n_trials)) + prob <- apply(outcomes, 1, function(k) prod(stats::dbinom(k, n_trials, p))) + expected_cov <- Reduce(`+`, lapply(seq_len(nrow(outcomes)), function(o) { + prob[o] * stats::cov(cbind(outcomes[o, ], n_trials - outcomes[o, ])) + })) + # observed counts whose sample covariance is exactly expected_cov, which has + # rank 1 because the two counts of a participant sum to n_trials + eig <- eigen(expected_cov, symmetric = TRUE) + centred <- outer(c(1, -1, 0), sqrt(eig$values[1]) * eig$vectors[, 1]) + observed <- sweep(centred, 2, n_trials * c(0.5, 0.5), "+") + colnames(observed) <- c("old", "new") + + data <- data.frame(id = factor(1:3), item_type = "old") + expected <- cbind(n_trials * p, n_trials * (1 - p)) + local_ppp_counts(observed, ppp_toy_draws(expected), ppp_toy_draws(expected), data) + res <- ppp(ppp_toy_fit(data, bmf(D ~ 1 + (1 | id), g ~ 1)), statistic = "T2") + expect_equal(unname(attr(res, "draws")$T2[, "observed"]), 0) +}) + +test_that("ppp() leaves participants without a matching set of cells out of T2", { + data <- data.frame(id = factor(c(1, 2, 3, 1, 2)), item_type = c(rep("old", 3), "new", "new")) + observed <- matrix(c(7, 5, 9, 2, 4, 3, 5, 1, 8, 6), 5) + local_ppp_counts(observed, ppp_toy_draws(observed + 0.5), ppp_toy_draws(observed), data) + fit <- ppp_toy_fit(data, bmf(D ~ 1 + (1 | id), g ~ 1)) + + expect_message(res <- ppp(fit, statistic = "T2"), "leaves out 1 participant") + expect_true(is.finite(res$observed)) +}) + +test_that("ppp() skips impossible categories without Inf or NaN", { + data <- data.frame(id = factor(1:3), item_type = "old") + observed <- cbind(a = c(4, 5, 6), b = c(6, 5, 4), c = 0) + expected <- cbind(c(5, 5, 5), c(5, 5, 5), 1e-43) + replicated <- cbind(c(5, 4, 6), c(5, 6, 4), 0) + local_ppp_counts(observed, ppp_toy_draws(expected), ppp_toy_draws(replicated), data) + fit <- ppp_toy_fit(data, bmf(D ~ 1 + (1 | id), g ~ 1)) + + res <- ppp(fit) + expect_true(all(is.finite(res$observed)) && all(is.finite(res$replicated))) + expect_equal(unname(attr(res, "draws")$T1[, "observed"]), 0) + # the same counts without the impossible category give the same statistics + local_ppp_counts(observed[, 1:2], ppp_toy_draws(expected[, 1:2]), + ppp_toy_draws(replicated[, 1:2]), data) + expect_equal(unclass(ppp(fit))$observed, res$observed) +}) + +test_that("ppp() leaves out categories whose expected count is exactly 0", { + data <- data.frame(id = factor(1:3), item_type = "old") + observed <- cbind(a = c(4, 5, 6), b = c(6, 5, 4), c = 0) + expected <- cbind(c(5, 4, 6), c(5, 6, 4), 0) + replicated <- cbind(c(6, 3, 6), c(4, 7, 4), 0) + fit <- ppp_toy_fit(data, bmf(D ~ 1 + (1 | id), g ~ 1)) + local_ppp_counts(observed, ppp_toy_draws(expected), ppp_toy_draws(replicated), data) + group <- attr(ppp(fit), "draws") + person <- attr(ppp(fit, level = "participant"), "draws") + + local_ppp_counts(observed[, 1:2], ppp_toy_draws(expected[, 1:2]), + ppp_toy_draws(replicated[, 1:2]), data) + expect_equal(group, attr(ppp(fit), "draws")) + expect_equal(person, attr(ppp(fit, level = "participant"), "draws")) +}) + +test_that("T2 leaves out a cell in which a single category is possible", { + # e.g. an m3 set size without "other" candidates: every response is correct + data <- data.frame(id = factor(rep(1:3, 2)), item_type = rep(c("old", "new"), each = 3)) + observed <- cbind(old = c(10, 10, 10, 4, 6, 5), new = c(0, 0, 0, 6, 4, 5)) + expected <- cbind(old = c(10, 10, 10, 5, 4, 6), new = c(1e-43, 1e-43, 1e-43, 5, 6, 4)) + replicated <- cbind(old = c(10, 10, 10, 5, 3, 6), new = c(0, 0, 0, 5, 7, 4)) + fit <- ppp_toy_fit(data, bmf(D ~ 1 + (1 | id), g ~ 1)) + local_ppp_counts(observed, ppp_toy_draws(expected), ppp_toy_draws(replicated), data) + with_forced_cell <- attr(ppp(fit, statistic = "T2"), "draws")$T2 + + rows <- 4:6 + local_ppp_counts(observed[rows, ], ppp_toy_draws(expected[rows, ]), + ppp_toy_draws(replicated[rows, ]), data[rows, ]) + expect_equal(with_forced_cell, attr(ppp(fit, statistic = "T2"), "draws")$T2) +}) + +test_that("ppp() refuses observed counts in an impossible category", { + data <- data.frame(item_type = "old") + observed <- cbind(a = 4, b = 5, c = 1) + expected <- cbind(5, 5, 1e-43) + local_ppp_counts(observed, ppp_toy_draws(expected), ppp_toy_draws(expected), data) + expect_error(ppp(ppp_toy_fit(data), statistic = "T1"), "probability 0 to observed responses in categories 'c'") +}) + +test_that("ppp() asks for the participant column when it cannot find one", { + data <- data.frame(id = factor(1:2), item = factor(1:2), item_type = "old") + observed <- matrix(c(4, 6, 6, 4), 2) + local_ppp_counts(observed, ppp_toy_draws(observed), ppp_toy_draws(observed), data) + + no_re <- ppp_toy_fit(data) + expect_message(res <- ppp(no_re), "computes T1 only") + expect_equal(res$statistic, "T1") + expect_error(ppp(no_re, statistic = "T2"), "Name the column with 'participant'") + expect_error(ppp(no_re, level = "participant"), "level = 'participant' needs") + expect_no_message(ppp(no_re, statistic = "T2", participant = "id")) + + crossed <- ppp_toy_fit(data, bmf(D ~ 1 + (1 | id) + (1 | item), g ~ 1)) + expect_error(ppp(crossed), "several grouping variables \\('id', 'item'\\)") + expect_no_error(ppp(crossed, statistic = "T1")) +}) + +test_that("ppp() says where the participant column has to be", { + # brms keeps only formula variables, so a model without person effects has + # no participant column in its data + data <- data.frame(item_type = c("old", "old")) + observed <- matrix(c(4, 6, 6, 4), 2) + local_ppp_counts(observed, ppp_toy_draws(observed), ppp_toy_draws(observed), data) + no_re <- ppp_toy_fit(data) + expect_match(tryCatch(ppp(no_re), message = conditionMessage), + "contain only the variables of the formula") + expect_match(tryCatch(ppp(no_re, participant = "id"), error = conditionMessage), + "contain only the variables of the formula") +}) + +test_that("ppp() checks its arguments", { + data <- data.frame(id = factor(1:2), item_type = "old") + observed <- matrix(c(4, 6, 6, 4), 2) + local_ppp_counts(observed, ppp_toy_draws(observed), ppp_toy_draws(observed), data) + fit <- ppp_toy_fit(data, bmf(D ~ 1 + (1 | id), g ~ 1)) + + expect_error(ppp(fit, level = "participant", statistic = "T2"), "group-level statistic") + expect_error(ppp(fit, participant = "subject"), "'participant' must name a column") + expect_error(ppp(fit, by = "cond"), "'by' must name columns") + expect_error(ppp(fit, statistic = "T3"), "should be one of") + + not_counts <- fit + not_counts$bmm$model <- mixture2p(resp_error = "y") + expect_error(ppp(not_counts), "available for the count models") +}) + +test_that("ppp() returns one T1 per participant at the participant level", { + data <- data.frame(id = factor(c(1, 2, 1, 2)), item_type = rep(c("old", "new"), each = 2)) + observed <- matrix(c(7, 5, 2, 4, 3, 5, 8, 6), 4) + expected <- matrix(c(6, 6, 3, 3, 4, 4, 7, 7), 4) + local_ppp_counts(observed, ppp_toy_draws(expected), ppp_toy_draws(expected), data) + fit <- ppp_toy_fit(data, bmf(D ~ 1 + (1 | id), g ~ 1)) + + res <- ppp(fit, level = "participant") + expect_equal(as.character(res$participant), c("1", "2")) + # participant 1: rows 1 and 3, each its own cell + expect_equal(res$observed[1], 1 / 6 + 1 / 4 + 1 / 3 + 1 / 7) + expect_equal(res$ppp, c(0, 0)) + expect_output(print(res), "'item_type' \\(2 cells\\)") +}) + +test_that(".ppp_complement() adds the second category of a binomial count", { + draws <- matrix(c(3, 4, 7, 8), 2) + out <- .ppp_complement(draws, trials = c(10, 20)) + expect_equal(dim(out), c(2L, 2L, 2L)) + expect_equal(out[, , 1], draws) + expect_equal(out[, , 2], matrix(c(7, 6, 13, 12), 2)) +}) + +test_that("ppp() compares each replicated data set with the expectation of its own draw", { + path <- test_path("assets/bmmfit_m3_ppcheck.rds") + skip_if_not(file.exists(path), "M3 fixture not available (excluded by .Rbuildignore)") + fit <- readRDS(path) + withr::local_seed(1) + res <- ppp(fit, statistic = "T1", by = character(0), ndraws = 5) + expect_s3_class(res, "bmm_ppp") + + # the same random numbers as ppp(): the draw ids first, then the simulation + withr::local_seed(1) + draw_ids <- sample.int(brms::ndraws(fit), 5) + yrep <- brms::posterior_predict(fit, draw_ids = draw_ids) + mu <- brms::posterior_epred(fit, draw_ids = draw_ids) + pooled <- function(x) apply(x, c(1L, 3L), sum) + x2 <- function(n) rowSums((n - pooled(mu))^2 / pooled(mu)) + observed <- matrix(colSums(unclass(fit$data$Y)), 5, dim(mu)[3], byrow = TRUE) + expect_equal(unname(attr(res, "draws")$T1[, "observed"]), x2(observed)) + expect_equal(unname(attr(res, "draws")$T1[, "replicated"]), x2(pooled(yrep))) +}) + +test_that("ppp() takes ndraws as a whole number of draws", { + path <- test_path("assets/bmmfit_m3_ppcheck.rds") + skip_if_not(file.exists(path), "M3 fixture not available (excluded by .Rbuildignore)") + fit <- readRDS(path) + expect_error(ppp(fit, ndraws = 2.5), "'ndraws' must be a whole number") +}) diff --git a/vignettes/articles/REFERENCES.bib b/vignettes/articles/REFERENCES.bib index e6f44b564..a88613924 100644 --- a/vignettes/articles/REFERENCES.bib +++ b/vignettes/articles/REFERENCES.bib @@ -261,6 +261,94 @@ @article{ratcliff_EstimatingParametersDiffusion_2002 } +@article{rieferBatchelder1988, + title = {Multinomial modeling and the measurement of cognitive processes}, + author = {Riefer, David M. and Batchelder, William H.}, + year = {1988}, + journal = {Psychological Review}, + volume = {95}, + number = {3}, + pages = {318--339}, + doi = {10.1037/0033-295X.95.3.318} +} + +@article{batchelderRiefer1999, + title = {Theoretical and empirical review of multinomial process tree modeling}, + author = {Batchelder, William H. and Riefer, David M.}, + year = {1999}, + journal = {Psychonomic Bulletin \& Review}, + volume = {6}, + number = {1}, + pages = {57--86}, + doi = {10.3758/BF03210812} +} + +@article{erdfelderEtAl2009, + title = {Multinomial processing tree models: {A} review of the literature}, + author = {Erdfelder, Edgar and Auer, Tina-Sarah and Hilbig, Benjamin E. and A{\ss}falg, Andr{\'e} and Moshagen, Morten and Nadarevic, Lena}, + year = {2009}, + journal = {Zeitschrift f{\"u}r Psychologie / Journal of Psychology}, + volume = {217}, + number = {3}, + pages = {108--124}, + doi = {10.1027/0044-3409.217.3.108} +} + +@article{oberauerBindings2019, + title = {Working memory capacity limits memory for bindings}, + author = {Oberauer, Klaus}, + year = {2019}, + journal = {Journal of Cognition}, + volume = {2}, + number = {1}, + pages = {40}, + doi = {10.5334/joc.86} +} + +@article{alhadhramiEtAl2025, + title = {A multinomial model-based analysis of bindings in working memory}, + author = {Al Hadhrami, Suaad S. and Bartsch, Lea M. and Oberauer, Klaus}, + year = {2025}, + journal = {Psychological Review}, + volume = {132}, + number = {4}, + pages = {828--856}, + doi = {10.1037/rev0000540} +} + +@article{singmannKellen2013, + title = {{MPTinR}: {Analysis} of multinomial processing tree models in {R}}, + author = {Singmann, Henrik and Kellen, David}, + year = {2013}, + journal = {Behavior Research Methods}, + volume = {45}, + number = {2}, + pages = {560--575}, + doi = {10.3758/s13428-012-0259-0} +} + +@article{heckEtAl2018, + title = {{TreeBUGS}: {An} {R} package for hierarchical multinomial-processing-tree modeling}, + author = {Heck, Daniel W. and Arnold, Nina R. and Arnold, Denis}, + year = {2018}, + journal = {Behavior Research Methods}, + volume = {50}, + number = {1}, + pages = {264--284}, + doi = {10.3758/s13428-017-0869-7} +} + +@article{klauer2010, + title = {Hierarchical multinomial processing tree models: {A} latent-trait approach}, + author = {Klauer, Karl Christoph}, + year = {2010}, + journal = {Psychometrika}, + volume = {75}, + number = {1}, + pages = {70--98}, + doi = {10.1007/s11336-009-9141-0} +} + @book{Green_Swets_1966, title = {Signal Detection Theory and Psychophysics}, author = {Green, David M. and Swets, John A.}, diff --git a/vignettes/articles/bmm.Rmd b/vignettes/articles/bmm.Rmd index 99ec79117..5d3172002 100644 --- a/vignettes/articles/bmm.Rmd +++ b/vignettes/articles/bmm.Rmd @@ -213,10 +213,26 @@ responses are the topic of the knitr::kable(bmm:::model_overview(group = "Choices and response times")) ``` +## Processing-tree models + +Participants give categorical responses, and you describe the discrete +processes that lead to each response as a multinomial processing tree, for +example the two-high-threshold model of recognition or the pair-clustering +model of free recall. You write each tree's branch probabilities as +expressions, and `mpt()` fits them to response counts. The +[MPT article](https://popov-lab.github.io/bmm/articles/bmm_mpt.html) covers +specifying trees, importing model files from MPTinR and TreeBUGS, and +constraining parameters. + +```{r, echo = FALSE} +knitr::kable(bmm:::model_overview(group = "Processing-tree models")) +``` + ```{r other-groups, results = "asis", echo = FALSE} known <- c( "Continuous reproduction", "Categorical recall and n-AFC decisions", - "Detection, recognition and confidence judgments", "Choices and response times" + "Detection, recognition and confidence judgments", "Choices and response times", + "Processing-tree models" ) registry <- bmm:::model_registry() other <- setdiff(unique(registry$group), known) diff --git a/vignettes/articles/bmm_mpt.Rmd b/vignettes/articles/bmm_mpt.Rmd new file mode 100644 index 000000000..470a3fa57 --- /dev/null +++ b/vignettes/articles/bmm_mpt.Rmd @@ -0,0 +1,716 @@ +--- +title: "Multinomial Processing Tree (MPT) models" +output: bookdown::html_document2 +author: + - Gidon Frischkorn +bibliography: REFERENCES.bib +header-includes: + - \usepackage{amsmath} +vignette: > + %\VignetteIndexEntry{Multinomial Processing Tree (MPT) models} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +pkgdown: + as_is: true +--- + +```{=html} + +``` +```{r, include = FALSE} +options(crayon.enabled = TRUE) +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + dev = "jpeg", + dpi = 100, + fig.asp = 0.8, + fig.width = 5, + out.width = "80%", + fig.align = "center" +) +fansi::set_knit_hooks(knitr::knit_hooks, which = c("output", "message", "error")) +``` + +# Introduction to the model + +Multinomial Processing Tree (MPT) models are a family of measurement models for tasks with categorical responses [@rieferBatchelder1988; @batchelderRiefer1999; @erdfelderEtAl2009]. They assume that each observed response results from a sequence of discrete latent processing states, for example detecting an item in memory and, if detection fails, guessing. Each latent process succeeds with some probability, and the branches of the resulting processing tree map every combination of latent outcomes onto one observable response category. The probability of a response category is the sum of the probabilities of all branches that terminate in it, and each branch probability is the product of the process probabilities along its path. + +The classic example is the two-high-threshold model (2HTM) of recognition memory. In an old/new recognition test, participants see two kinds of probe, targets (items studied before) and lures (new items), and respond "old" or "new" to each. The 2HTM assumes different processes for the two kinds of probe. A target is detected as old with probability $D$; if detection fails (probability $1 - D$), the participant guesses "old" with probability $g$. A lure is detected as new with probability $D$ (the two detection probabilities can also be separated), and otherwise guessed on in the same way. The model therefore has one tree per probe type, not one per response: the target tree and the lure tree each end in both responses, "old" and "new". + +```{r fig-tree-2htm-intro, echo = FALSE, fig.width = 8, fig.asp = 0.45, out.width = "100%"} +bmm::mpt( + trees = list( + bmm::mpt_tree("target", old ~ D + (1 - D) * g, new ~ (1 - D) * (1 - g)), + bmm::mpt_tree("lure", old ~ (1 - D) * g, new ~ D + (1 - D) * (1 - g)) + ), + tree_id = "item_type" +) |> + plot() +``` + +Each panel is one tree, titled with the probe type it describes. Its leaves are the response categories, "old" and "new", which both trees share. A response category collects the probability of every leaf that ends in it, which gives: + +$$ +\begin{aligned} +P(\text{"old"} \mid \text{target}) &= D + (1 - D) \cdot g \\ +P(\text{"new"} \mid \text{target}) &= (1 - D) \cdot (1 - g) \\ +P(\text{"old"} \mid \text{lure}) &= (1 - D) \cdot g \\ +P(\text{"new"} \mid \text{lure}) &= D + (1 - D) \cdot (1 - g) +\end{aligned} +$$ + +In `bmm`, MPT models are fit as hierarchical Bayesian models: all latent probability parameters are estimated on an unconstrained latent scale (logit by default, or probit), where they can be predicted by fixed and random effects using the standard `bmmformula` syntax. With a probit link and correlated random effects across all parameters, this is the latent-trait approach of @klauer2010 as implemented, for example, in TreeBUGS [@heckEtAl2018]. + +The article follows the order of a typical analysis: specifying the trees, preparing and simulating data, fitting, interpreting the estimates, constraining parameters, and two extensions (covariates inside branch expressions, non-linear parameter formulas), followed by a real-data example. If you already have a model file for multiTree, MPTinR or TreeBUGS, jump to [Importing existing model definitions]; if you are coming from MPTinR or TreeBUGS, read [Specifying an MPT model] first and then [Migrating from MPTinR or TreeBUGS], which maps the functions of those packages onto `bmm`. + +# Specifying an MPT model + +An MPT model is specified as a set of trees. Each tree is created with `mpt_tree()`: it has a name, and one branch-probability expression per response category, written as a formula `category ~ expression`: + +```{r} +library(bmm) + +tree_target <- mpt_tree("target", + old ~ D + (1 - D) * g, + new ~ (1 - D) * (1 - g) +) +tree_lure <- mpt_tree("lure", + old ~ (1 - D) * g, + new ~ D + (1 - D) * (1 - g) +) +tree_target +``` + +The first argument names the probe type the tree describes; the left-hand sides are the response categories, and every tree lists the same ones. The formulas only carry the expressions and the response category on the left-hand side; nothing in them is evaluated. The expressions can equally be given as a named list of character strings, which is the form that `mpt_from_string()` and `mpt_from_eqn()` produce. The two spellings give identical trees: + +```{r} +tree_target_strings <- mpt_tree("target", list( + old = "D + (1 - D) * g", + new = "(1 - D) * (1 - g)" +)) +identical(tree_target, tree_target_strings) +``` + +The trees are combined into a model object with `mpt()`. The `tree_id` argument names the data column whose values say which tree each observation belongs to; here `item_type`, holding `"target"` and `"lure"`, which must match the tree names. The response counts go in columns named after the response categories, `old` and `new` (see [Data format and simulation]). Printing the model lists the trees and then checks identifiability in two ways: it compares the number of free parameters with the degrees of freedom the trees provide (response categories minus 1, summed over trees), the classical bound that MPTinR's `check.mpt()` reports, and it computes the rank of the Jacobian of the category probabilities at interior parameter values: + +```{r} +model_2htm <- mpt( + trees = list(tree_target, tree_lure), + tree_id = "item_type" +) +model_2htm +``` + +The two checks answer different questions. The classical bound counts parameters against the information the response categories carry; the rank asks whether the branch expressions actually separate the parameters they contain, which can fail even when the counting works out. *Locally identified* means the rank is full at the test values, so no direction in parameter space leaves every category probability unchanged there. Both checks assume intercept-only formulas, which is why the printed line says so. `bmm()` follows two rules: it warns when the rank is deficient and every formula is intercept-only, and it only says that predictors or non-linear formulas may identify the parameters, without checking, when any formula has a predictor or is non-linear. It does not look at the data. `parameter_info()` lists the parameters with their links and any fixed values. + +`tree_id` is about *structure*, not about experimental manipulations. A separate tree is needed only when the branch expressions themselves differ, that is, when the trial determines which latent processes apply or which response categories are reachable at all. Effects of experimental factors on the parameters go in the parameter formulas, exactly as in every other `bmm` model. The two roles are independent, and the same column can serve both when a factor happens to change the structure as well. + +The `plot()` method draws the processing trees implied by the branch expressions: each product in an expression is one root-to-leaf path, and paths sharing a common prefix share the corresponding nodes. This is a quick way to verify that the expressions encode the intended tree; the figure in the [Introduction to the model] was drawn with this call: + +```{r, eval = FALSE} +plot(model_2htm) +``` + +The latent parameters are identified automatically from the branch expressions: every symbol that is not a declared covariate (see [Design-fixed covariates in branch expressions]) is treated as a parameter. Equality constraints follow directly from the naming. Because both trees use the same symbol `D`, detection is equated for targets and lures. To estimate separate detection probabilities, write `Do` and `Dn` instead: + +```{r} +tree_target_sep <- mpt_tree("target", list( + old = "Do + (1 - Do) * g", + new = "(1 - Do) * (1 - g)" +)) +tree_lure_sep <- mpt_tree("lure", list( + old = "(1 - Dn) * g", + new = "Dn + (1 - Dn) * (1 - g)" +)) +model_2htm_sep <- mpt( + trees = list(tree_target_sep, tree_lure_sep), + tree_id = "item_type" +) +model_2htm_sep +``` + +The printed note is the well-known result that a 2HTM with separate detection parameters is not identified from a single old/new test: three free parameters, two degrees of freedom, and a Jacobian of rank 2, so one combination of the three cannot be estimated. The model becomes identifiable with a second condition that varies the response bias (different base rates or payoffs), which adds two degrees of freedom and one parameter, or with a restriction (see [Constraining parameters]). The order constraint $D_o > D_n$ shown there is a different kind of remedy: it keeps the two detection parameters ordered in every posterior draw, but it does not add information, so the rank deficit remains and $D_n$ stays weakly determined. + +Two things to keep in mind when writing branch expressions: + +- Parameter and category names must consist of letters and digits only (no underscores or dots), because they become non-linear parameters in **brms**. +- Numeric fractions such as `1/4` are converted to decimal literals (`0.25`) when the tree is created. In Stan, a bare integer fraction would be evaluated as integer division. + +The default link is the logit. `links = "probit"` selects the probit link that latent-trait MPT models traditionally use for every parameter, and switches the default priors accordingly (see [Default priors]); a named list, `links = list(D = "probit")`, switches single parameters and their priors: + +```{r} +model_2htm_probit <- mpt( + trees = list(tree_target, tree_lure), + tree_id = "item_type", + links = "probit" +) +``` + +# Data format and simulation + +The data contain one column of aggregated response counts per response category, named after the branch names of the trees, plus the `tree_id` column. Aggregated counts are sufficient whenever parameters and covariates are constant within participant and condition. Trial-level data (one row per trial with a 1 in the chosen category's column) also work and are only needed when predictors vary from trial to trial. + +The distribution functions `dmpt()` and `rmpt()` provide the multinomial density and random generation for any `mpt()` model. `rmpt()` draws the counts for one tree and one parameter vector; a short function around it simulates a hierarchical data set in the shape of TreeBUGS' `genTraitMPT()`: population means on the probability scale, participant heterogeneity as standard deviations on the latent scale (here the logit scale; the `sigma` of `genTraitMPT()` is on the probit scale, see [Default priors] for the conversion). + +```{r} +simulate_2htm <- function(n_subj, n_items, mean, sd) { + latent <- data.frame( + D = rnorm(n_subj, qlogis(mean[["D"]]), sd[["D"]]), + g = rnorm(n_subj, qlogis(mean[["g"]]), sd[["g"]]) + ) + dat <- expand.grid( + id = factor(seq_len(n_subj)), + item_type = c("target", "lure"), + stringsAsFactors = FALSE + ) + counts <- t(sapply(seq_len(nrow(dat)), function(i) { + subj <- as.integer(dat$id[i]) + rmpt( + n = 1, size = n_items, + pars = plogis(c(D = latent$D[subj], g = latent$g[subj])), + mpt_model = model_2htm, tree = dat$item_type[i] + )[1, ] + })) + cbind(dat, counts) +} + +set.seed(124) +dat <- simulate_2htm( + n_subj = 40, n_items = 100, + mean = c(D = 0.7, g = 0.5), sd = c(D = 0.4, g = 0.3) +) +head(dat) +``` + +Each participant contributes two rows, one per probe type, and the counts in the `old` and `new` columns sum to 100 within a row. + +# Fitting the model + +Model parameters are predicted with the usual `bmmformula` syntax. Here both parameters get a population intercept and a random intercept per participant. Written as `(1 | id)` in each formula, the two person effects are independent of each other, which matches the simulation above: it draws `D` and `g` independently. A shared label, `(1 |p| id)` in both formulas, would estimate their correlation as well, which is the latent-trait model; [Migrating from MPTinR or TreeBUGS] discusses these options. + +```{r fit-2htm, results = "hide", message = FALSE, warning = FALSE} +fit_2htm <- bmm( + formula = bmf( + D ~ 1 + (1 | id), + g ~ 1 + (1 | id) + ), + data = dat, + model = model_2htm, + backend = "cmdstanr", + cores = 4, + refresh = 0, + file = "fits/mpt_2htm_v2" +) +``` + +```{r} +summary(fit_2htm) +``` + +All estimates are on the logit scale; the next section converts them to probabilities. + +## Default priors + +The default priors are matched to the link function: with the logit link, intercepts and regression effects both get a `logistic(0, 1)` prior, which for an intercept is a uniform distribution on the probability scale. Effects get the same width as intercepts because a narrower prior pulls condition differences toward zero when the formula is written as `1 + cond`. With `links = "probit"`, both get `normal(0, 1)`, matching the priors typically used in latent-trait MPT models. Random-effects standard deviations get `exponential(1)` with the logit link and `exponential(2)` with the probit link, because a probit SD is about 1/1.7 of the logit SD it corresponds to; the rate of 1.7 this implies is rounded to 2. + +```{r} +default_prior(bmf(D ~ 1, g ~ 1), data = dat, model = model_2htm) +``` + +A prior given as `prior("normal(0, 3.1623)", class = "b", nlpar = "D")` replaces the defaults of every coefficient of `D`, its intercept included; add `coef = "Intercept"` to set the intercept alone. + +# Interpreting the estimates + +## Population-level probabilities + +`native_parameters()` turns the latent-scale estimates into probabilities. With `re_formula = NA` it sets all person effects to zero, and `summary = TRUE` summarizes the transformed draws: + +```{r} +native_parameters(fit_2htm, re_formula = NA, summary = TRUE) +``` + +The function applies the inverse link to the linear predictor of each combination of predictors, not to individual coefficients. That is why it also works for models with predictors, where the inverse link applied to an effect coefficient is not a probability, and for the probit link. Summarizing after the transformation matters for the point estimate: transforming the summary of the latent estimates instead returns correct interval bounds, because quantiles survive a monotone transformation, but a wrong `Estimate`, because the mean of a transformed parameter is not the transform of its mean. + +These are the probabilities of a participant at the latent mean, that is, of the median participant. They are not the population mean of the participants' probabilities: because the inverse link is non-linear, that mean lies closer to 0.5, and the two differ noticeably when the latent standard deviation is near 1. The median participant is the quantity that TreeBUGS reports as `"mean"` in `getParam()`. `?native_parameters` (section "Group-level effects") explains the distinction. + +## Participant-level probabilities + +`native_parameters()` returns a parameter on the probability scale, population effects and random effects combined, for every combination of its predictors in the data. `D` depends on the participant only, so this is one row per participant; `summary = TRUE` summarizes the draws after the transformation, as above: + +```{r} +D_subj <- native_parameters(fit_2htm, pars = "D", summary = TRUE) +head(D_subj) +``` + +With `newdata`, the same call gives the probabilities for new participants or conditions, and `re_formula = NA` leaves out the random effects. For work on the latent scale, `brms::posterior_linpred()` with the `nlpar` argument returns the latent value of one parameter for every row of the data and every posterior draw. + +## Posterior predictive checks + +Because the model uses a multinomial family, `pp_check()` produces a response-proportion plot: the observed proportion of each response category against the posterior predictive distribution. + +```{r fig-ppc-2htm, fig.width = 7, fig.asp = 0.5, out.width = "100%"} +pp_check(fit_2htm) +``` + +The `group` argument splits the check by a predictor variable of the model. The tree column of this fit is not a predictor (no parameter formula uses it), so the check pools both probe types here; the real-data example in [Fitting real data] groups by its condition factor, which also selects the tree. Condition effects, hypothesis tests, a goodness-of-fit statistic and model comparison are shown there as well. + +# Constraining parameters + +MPT analyses routinely fix parameters or set them equal. `bmm` offers two places to do so, and the choice is a matter of when the constraint is decided: at construction with `restrictions`, or at fit time in the formula. + +## Equality restrictions and constants at construction + +The `restrictions` argument of `mpt()` accepts the string syntax of MPTinR and TreeBUGS. A restriction either equates a parameter with another one or fixes it to a constant, and chains such as `"G1 = G2 = G3"` map all earlier names onto the last. Restrictions are substituted into the branch expressions before the parameters are identified, so a restricted parameter is not part of the model at all. The separate-detection 2HTM with `Dn = Do` is the shared-detection model again: + +```{r} +model_2htm_restricted <- mpt( + trees = list(tree_target_sep, tree_lure_sep), + tree_id = "item_type", + restrictions = "Dn = Do" +) +model_2htm_restricted +``` + +A constant works the same way. `dmpt()` and the branch-sum check see the restricted expressions, and `plot()` labels the edges of `g` as `g = 0.5`: + +```{r} +model_2htm_g5 <- mpt( + trees = list(tree_target, tree_lure), + tree_id = "item_type", + restrictions = "g = 0.5" +) +names(model_2htm_g5$parameters) +dmpt( + x = c(old = 80, new = 20), pars = c(D = 0.6), + mpt_model = model_2htm_g5, tree = "target", log = FALSE +) +plot(model_2htm_g5) +``` + +Order constraints (`"Do > Dn"`) are not accepted by `restrictions`; they are expressed through the formula, as shown below. + +## Constants at fit time + +A parameter can also be fixed in the model formula, `bmf(D ~ 1 + (1 | id), g = 0.5)`. The parameter then stays part of the model, and freeing it again is a matter of giving it a formula. The rule of thumb: put a constant that defines the model into `restrictions`, and a constant that you compare against a free version into the formula. Fixed values are written on the probability scale in both cases; the formula mechanism maps them to the latent scale when the constant prior is built, so `default_prior()` shows `constant(0)` for `g = 0.5` under the logit link: + +```{r} +default_prior(bmf(D ~ 1 + (1 | id), g = 0.5), data = dat, model = model_2htm) +``` + +## Order constraints + +Requiring detection of targets to exceed detection of lures ($D_o > D_n$) needs no dedicated syntax: derive the larger parameter from the smaller one with a non-linear formula such that the constraint holds by construction. The recommended formulation interpolates on the probability scale: + +```{r, eval = FALSE} +fit_ordered <- bmm( + bmf( + Do ~ Dn + (1 - Dn) * inv_logit(phi), + Dn ~ 1 + (1 | id), + g ~ 1 + (1 | id), + phi ~ 1 + (1 | id) + ), + data = dat, + model = model_2htm_sep +) +``` + +Here `Dn` remains an ordinary logit-scale parameter with its usual priors, and `phi` is the (logit of the) proportion of the remaining probability space above $D_n$ that $D_o$ covers: $D_o = D_n + (1 - D_n) \cdot \text{logit}^{-1}(\phi)$, which lies in $(D_n, 1)$ for every posterior draw. Because `Dn` appears in a non-linear formula, its link transformation is applied automatically: you write plain `Dn` and get $\text{logit}^{-1}(D_n)$ in the model. Sub-parameters such as `phi` default to a `normal(0, 1)` prior; for a uniform prior on the covered proportion, replace that default with `prior("logistic(0, 1)", class = "b", coef = "Intercept", nlpar = "phi")`. + +An equivalent alternative shifts on the latent scale, as latent-trait MPT implementations usually do: `Dn ~ inv_logit(lDn)` and `Do ~ inv_logit(lDn + exp(delta))`, with `lDn` and `delta` receiving the predictor formulas. Both formulations enforce the order in every posterior draw. In a validation on simulated 2HTM data (40 participants, 100 trials per probe type, 10 data sets per gap size), the interpolation sampled without divergences at the default `adapt_delta`, whereas the latent-offset version produced divergences in every fit. A closer look at two of these fits tied the divergences to small values of the random-effect standard deviation of the log-gap parameter `delta`, not to a closing gap. Raising `adapt_delta` to 0.99 removed them; 0.95 did not. With a true gap of 0.20 between the two detection parameters, the 95% intervals of both detection parameters contained the generating values in all 10 data sets, under both formulations. With a gap of 0.04 they missed them in two or three of 10 data sets, and both formulations overestimated the gap. + +What an order constraint does not do is restore identifiability. Replacing `Do` by `phi` leaves three free parameters against the two degrees of freedom of a single old/new test, so `bmm()` still reports the rank deficit of the branch expressions, as a message here because the formula for `Do` is non-linear. The constraint makes the posterior proper and the ordering certain, and the hierarchical structure contributes what the aggregate counts cannot, but the lower parameter is identified mainly by the prior and the order: in the validation above, the 95% credible interval of $D_n$ averaged 0.3 to 0.4 in width, two to four times that of $D_o$ (0.1 to 0.2), and exceeded half the unit interval in 6 of the 20 interpolation fits. Read the interval before interpreting the smaller parameter, and prefer a design with a second response-bias condition when both detection parameters matter. + +One pitfall to avoid: do **not** write the latent offset against the parameter name directly, as in `Do ~ inv_logit(Dn + exp(delta))`. The automatic link substitution turns the inner `Dn` into $\text{logit}^{-1}(D_n)$, a probability, so the expression adds a probability and a latent-scale offset, and the resulting $D_o$ is *not* guaranteed to exceed $D_n$ (for example, $D_n = 0.9$ with a small offset yields $D_o \approx 0.71$). Constraint expressions must keep all terms on one scale: either work with probabilities (the interpolation above) or introduce explicit latent-scale sub-parameters like `lDn`. + +## Simplex-constrained parameter groups + +When several guessing probabilities must sum to 1, the parameters form a probability simplex. A source-monitoring model illustrates this: items come from source A or source B, or are new, and a participant who neither detects the item nor its source guesses one of the three responses with probabilities `gA`, `gB` and `gNew` that sum to 1. + +```{r} +tree_sourceA <- mpt_tree("sourceA", list( + A = "dA + (1 - dA) * gA", + B = "(1 - dA) * gB", + New = "(1 - dA) * gNew" +)) +tree_sourceB <- mpt_tree("sourceB", list( + A = "(1 - dB) * gA", + B = "dB + (1 - dB) * gB", + New = "(1 - dB) * gNew" +)) +tree_newitem <- mpt_tree("newitem", list( + A = "gA", + B = "gB", + New = "gNew" +)) + +model_smm <- mpt( + trees = list(tree_sourceA, tree_sourceB, tree_newitem), + tree_id = "source", + simplex = c("gA", "gB", "gNew") +) +model_smm +``` + +Declaring the group via the `simplex` argument applies a stick-breaking reparameterization: the last parameter of each group is derived as 1 minus the sum of the others, and each remaining parameter `p` is estimated through an unconstrained component `praw`, which is what the `gAraw` and `gBraw` entries in the parameter list are. Predictor formulas are written for the simplex parameters themselves and are applied to the corresponding components; the derived parameter (here `gNew`) cannot receive its own predictors. Several independent simplex groups can be declared by passing a list of character vectors. + +By default, the group gets a uniform prior over the simplex, a symmetric Dirichlet(1): every way of splitting the guessing probability among `gA`, `gB` and `gNew` is equally likely a priori, so each has prior mean 1/3, whatever the order in which the group lists them. Under stick breaking this is a Beta(1, 2) prior on the first stick (`gA`'s share) and a Beta(1, 1) prior on the second (`gB`'s share of what `gA` leaves), which `default_prior()` shows on the latent scale as `logitbeta(1, 2)` for `gAraw` and `logistic(0, 1)` for `gBraw`. The 1/3 holds in every cell when the members get cell-means formulas (`gA ~ 0 + cond`); under `gA ~ 1 + cond` it holds in the reference cell only. TreeBUGS has no simplex construct and puts a uniform prior on each branch parameter instead, which here would give prior means of 0.50, 0.25 and 0.25. To match it, replace the first stick's prior: + +```{r, eval = FALSE} +prior_treebugs <- brms::set_prior("logistic(0, 1)", nlpar = "gAraw") +``` + +```{r, eval = FALSE} +fit_smm <- bmm( + bmf(dA ~ 1 + (1 | id), dB ~ 1 + (1 | id), gA ~ 1, gB ~ 1), + data = source_data, + model = model_smm +) +``` + +# Design-fixed covariates in branch expressions + +Some MPT models contain probabilities that are fixed by the experimental design rather than estimated, for example guessing rates determined by the number of response options. Such quantities can be included in branch expressions as *covariates*: data columns that are declared in the `covariates` argument and passed through to the model unchanged. + +The model @oberauerBindings2019 used for n-AFC recognition data illustrates this. The probability of a correct response depends on binding memory (`Pb`) and item memory (`Pi`), while three guessing rates are computed from the size of the response set in each condition: + +```{r} +tree_ss <- mpt_tree("ss", list( + correct = "Pb + (1 - Pb) * Pi * GcorrPi + (1 - Pb) * (1 - Pi) * GcorrNoPi", + other = "(1 - Pb) * Pi * (1 - GcorrPi) + + (1 - Pb) * (1 - Pi) * (1 - GcorrNoPi) * GotherNoPi", + npl = "(1 - Pb) * (1 - Pi) * (1 - GcorrNoPi) * (1 - GotherNoPi)" +)) + +model_ss <- mpt( + trees = tree_ss, + covariates = c("GcorrPi", "GcorrNoPi", "GotherNoPi") +) +model_ss +``` + +Only `Pb` and `Pi` are latent parameters; the `Gcorr*` symbols must exist as numeric columns in the data (for example, `GcorrPi = 1 / rsizeList`). This mechanism goes beyond what classical MPT software supports, because the guessing rates can vary continuously across trials or conditions. When the data are checked, the branch sums are re-evaluated with the observed covariate values row by row, and `bmm()` stops when covariate columns violate their joint constraints. + +# Fitting real data + +The `oberauer_lewandowsky_2019_e1` data set that the [m3 article](https://popov-lab.github.io/bmm/articles/bmm_m3.html) analyzes with activation-based measurement models can equally be described by a processing-tree model, which makes for an instructive real-data example. Participants recalled a list item from a set of candidates in three conditions: in the `new distractors` condition the response set additionally contained five distractors, while in the `old reordered` and `old same` conditions no distractors were present. Each row records how often a participant chose the correct item (`corr`, 1 candidate), another list item (`other`, 4 candidates), a distractor (`dist`, 5 or 0 candidates), or a not-presented lure (`npl`, 5 or 10 candidates). The participant column of this data set is called `ID`: + +```{r} +data_ol <- oberauer_lewandowsky_2019_e1 +data_ol$cond <- factor( + data_ol$cond, + levels = c("new distractors", "old reordered", "old same"), + labels = c("newdist", "oldreord", "oldsame") +) +head(data_ol[, 1:6]) +``` + +We describe the responses by three processes, ordered by how much they constrain the response. First, with probability `Pi`, the list items themselves are remembered. Given item memory, with probability `Pb` the item bound to the probed position is retrieved, yielding a correct response. If the binding is not available, distractor activation may still persist (probability `Pd`) and capture the response; this stage exists only in the distractor condition, otherwise the response is a guess among the five known list candidates. Without item memory, the response is a guess across the full response set of 15 candidates. This yields two tree structures, and illustrates a parameter that exists only in some conditions: `Pd` appears only in the branch expressions of the `newdist` tree, so it is estimated exclusively from that condition, without any indicator-variable bookkeeping. + +One response category is impossible in the no-distractor conditions: `dist` responses cannot occur there. Listing it in the `impossible` argument of `mpt_tree()` removes it from that tree's branch list (the remaining branches then have to sum to 1 among themselves), and `bmm` gives the category a large negative linear predictor for the rows concerned. Because the model is estimated on the log scale and the multinomial family renormalizes across categories, this is an exact zero for practical purposes, and no probability has to be taken away from the other branches to compensate. Written with candidate-count fractions, which `mpt_tree()` folds into decimal literals, the trees are: + +```{r} +tree_newdist <- mpt_tree("newdist", list( + corr = "Pi * Pb + Pi * (1 - Pb) * (1 - Pd) * (1/5) + (1 - Pi) * (1/15)", + other = "Pi * (1 - Pb) * (1 - Pd) * (4/5) + (1 - Pi) * (4/15)", + dist = "Pi * (1 - Pb) * Pd + (1 - Pi) * (5/15)", + npl = "(1 - Pi) * (5/15)" +)) +tree_nodist <- mpt_tree("nodist", list( + corr = "Pi * Pb + Pi * (1 - Pb) * (1/5) + (1 - Pi) * (1/15)", + other = "Pi * (1 - Pb) * (4/5) + (1 - Pi) * (4/15)", + npl = "(1 - Pi) * (10/15)" +), impossible = "dist") + +data_ol$tree <- ifelse(data_ol$cond == "newdist", "newdist", "nodist") + +model_ol <- mpt( + trees = list(tree_newdist, tree_nodist), + tree_id = "tree" +) +``` + +Note that `tree` is a coarser column than `cond`: three experimental conditions, two branch structures. Both no-distractor conditions get the same tree, and the difference between them is estimated where it belongs, in the parameter formulas (`Pi ~ 0 + cond` below). Writing the same branch list twice under two names would fit the same model but cost real work, because the category probabilities are built as a sum over all trees and every duplicate is evaluated for every row. If all conditions shared one structure, `tree_id` could be dropped altogether and a single tree given. + +The two tree structures: + +```{r fig-tree-ol, fig.width = 9, fig.asp = 0.44, out.width = "100%"} +plot(model_ol) +``` + +## Condition effects + +`Pi` and `Pb` get one population parameter per condition (`0 + cond`) and, per participant, one uncorrelated random effect per condition (`||`); `Pd` exists in one condition only and takes a plain random intercept. Without a shared `|p|` label the three parameters' person effects are also uncorrelated with each other. This keeps the model to standard deviations: correlating all seven person effects would add 21 correlations. Both fits in this section draw 3000 post-warmup samples per chain instead of the default 1000, which gives every parameter, the random effects included, a bulk effective sample size above 400: + +```{r fit-ol, results = "hide", message = FALSE, warning = FALSE} +fit_ol <- bmm( + formula = bmf( + Pi ~ 0 + cond + (0 + cond || ID), + Pb ~ 0 + cond + (0 + cond || ID), + Pd ~ 1 + (1 | ID) + ), + data = data_ol, + model = model_ol, + backend = "cmdstanr", + cores = 4, + iter = 4000, + warmup = 1000, + refresh = 0, + file = "fits/mpt_ol" +) +``` + +```{r} +native_parameters(fit_ol, re_formula = NA, summary = TRUE) +``` + +`Pd` has no predictor, so its single value is repeated for each condition; it is estimated from the `newdist` condition only. Item memory is high in all three conditions (not-presented lures are rarely chosen), while binding memory is clearly higher when the list is re-presented in its original order than in the other two conditions. When item memory succeeds but the binding is unavailable in the distractor condition, distractor activation captures the response about half of the time. + +## Testing differences between conditions + +`brms::hypothesis()` evaluates a contrast in every posterior draw. Parameters are addressed by their **brms** names without the `b_` prefix, that is, `_`; `variables(fit_ol)` lists them. The contrast is evaluated on the latent scale, where the difference of two logits is a log odds ratio: + +```{r} +brms::hypothesis(fit_ol, "Pb_condoldsame > Pb_condoldreord") +``` + +`Post.Prob` is the posterior probability that the difference is positive. The difference on the probability scale, if that is the quantity to report, is the difference of the probabilities in the two conditions, draw by draw: + +```{r} +Pb_draws <- native_parameters(fit_ol, pars = "Pb", re_formula = NA) +Pb_draws <- Pb_draws[order(Pb_draws$cond, Pb_draws$.draw), ] +Pb_by_cond <- split(Pb_draws$value, Pb_draws$cond) +brms::posterior_summary(Pb_by_cond$oldsame - Pb_by_cond$oldreord) +``` + +## Posterior predictive check + +The posterior predictive distribution reproduces the response-category proportions within each condition, including the near-zero `dist` proportions where no distractors were shown: + +```{r fig-ppc-ol, fig.width = 7, fig.asp = 0.7, out.width = "100%"} +pp_check(fit_ol, group = "cond") +``` + +## Goodness of fit: T1 and T2 + +The plot compares proportions by eye. A posterior predictive p-value (PPP) puts a number on it. @klauer2010 proposed two statistics for hierarchical MPT models. $T_1$ compares the observed frequency of each response category in each cell of the design with the frequency the model expects: + +$$ +T_1 = \sum_{c,k} \frac{(n_{ck} - \hat{n}_{ck})^2}{\hat{n}_{ck}}, +$$ + +summed over all categories $k$ of all cells $c$, where $n_{ck}$ is the observed and $\hat{n}_{ck}$ the expected frequency. $T_2$ compares the covariances of the category frequencies across participants with the covariances the model expects, so it checks whether the model reproduces how participants differ, not only the average. Each statistic is computed once for the observed counts and once for counts simulated from the model, both against the expected counts of the same posterior draw. The PPP is the proportion of draws in which the simulated data deviate more than the observed data; a PPP near 0 flags misfit. `ppp()` computes both: + +```{r ppp-ol} +ppp_ol <- ppp(fit_ol) +ppp_ol +``` + +By default, a cell combines the tree column with the categorical predictors of the formula, here `tree` and `cond`, which gives the three conditions; the `by` argument changes this. The participants for $T_2$ come from the grouping variable of the random effects, `ID`. `ppp()` sums the counts of a cell over participants, where @klauer2010 and TreeBUGS average them; with the same number of participants in every cell, as here, both give the same PPP. + +The $T_1$ PPP of `r round(ppp_ol$ppp[1], 2)` shows that the model reproduces the category frequencies of each condition. The $T_2$ PPP of `r round(ppp_ol$ppp[2], 3)` is small: data simulated from the model rarely deviate from the expected covariances as much as the observed data do. The model keeps all person effects uncorrelated, so a participant with a high `Pi` in one condition is not expected to have a high `Pi` in the other conditions as well. Correlated person effects, `(0 + cond | ID)` or the `|p|` syntax described under [Migrating from MPTinR or TreeBUGS], are the first extension to try. + +With `level = "participant"`, `ppp()` computes $T_1$ and its PPP for each participant separately, as TreeBUGS does, to find the participants the model fits poorly: + +```{r ppp-ol-participant} +ppp_ol_id <- ppp(fit_ol, level = "participant") +head(ppp_ol_id) +``` + +The smallest of the `r nrow(ppp_ol_id)` participant PPPs is `r round(min(ppp_ol_id$ppp), 2)`. + +## Model comparison + +Because a `bmm` fit is a `brmsfit`, the model-comparison tools of **brms** apply to it. Which of them can be trusted depends on the structure of the data: here each row aggregates all trials of one participant in one condition. As an example, take a model in which binding memory does not vary across conditions: + +```{r fit-ol-restricted, results = "hide", message = FALSE, warning = FALSE} +fit_ol_pb1 <- bmm( + formula = bmf( + Pi ~ 0 + cond + (0 + cond || ID), + Pb ~ 1 + (1 | ID), + Pd ~ 1 + (1 | ID) + ), + data = data_ol, + model = model_ol, + backend = "cmdstanr", + cores = 4, + iter = 4000, + warmup = 1000, + refresh = 0, + file = "fits/mpt_ol_pb1" +) +``` + +```{r} +loo_ol <- brms::loo(fit_ol) +loo_ol_pb1 <- brms::loo(fit_ol_pb1) +brms::loo_compare(loo_ol, loo_ol_pb1) +``` + +`loo()` is a quick screen here, not a result. It flags `r sum(loo::pareto_k_values(loo_ol) > 0.7)` of the `r length(loo::pareto_k_values(loo_ol))` rows of `fit_ol` and `r sum(loo::pareto_k_values(loo_ol_pb1) > 0.7)` of `fit_ol_pb1` with a Pareto-k above 0.7. That is expected with aggregated data: leaving out a row removes all trials of one participant in one condition, which shifts that participant's random effects considerably, and the importance-sampling approximation breaks down. With this many flagged rows, neither the elpd differences nor their standard errors are trustworthy, and the output above should not be reported as evidence for or against the restriction. The hypothesis test above answers the same question directly. + +The prediction target that matches the data is a new participant, so cross-validation should leave out all rows of a participant at once. In **brms** this is `kfold(fit, folds = "grouped", group = "ID")`, and it runs on a `bmmfit` of an MPT model. It refits the model once per fold, so the call costs K full fits where `loo()` costs none. That is why `loo()` stays the cheap screen and grouped `kfold()` is the check to run when the screen flags many rows. The article renders from cached fits, so the call is shown but not evaluated: + +```{r, eval = FALSE} +kfold_ol <- brms::kfold(fit_ol, folds = "grouped", group = "ID", K = 10) +``` + +Bayes factors between nested models, for example with and without person effects or with and without a condition effect, are available through `brms::bayes_factor()`, which runs bridge sampling [@gronauBridgesamplingPackageEstimating2020] on both fits. This needs two changes to the fitting call. First, `save_pars = save_pars(all = TRUE)`, because the sampler needs the draws of all parameters, including the person effects. Second, many more post-warmup draws than a posterior summary needs, so that the estimate is stable when repeated. The fits in this article were fitted without these settings, so the sampler cannot be run on them here. The [hypothesis testing article](https://popov-lab.github.io/bmm/articles/bmm_hypothesis_testing.html) shows the full workflow. + +## Relation to the m3 model + +This MPT is not merely inspired by the [m3 analysis of the same data](https://popov-lab.github.io/bmm/articles/bmm_m3.html): under the `simple` choice rule the two models are *the same likelihood* in different parameterizations. With the m3 activations `corr ~ b + a + c`, `other ~ b + a`, `dist ~ b + d`, `npl ~ b` (candidate counts 1, 4, 5, 5 and $S = 15b + 5a + c + 5d$), the mapping is + +$$ +P_i = \frac{5a + c + 5d}{S}, \qquad +P_b = \frac{c}{5a + c + 5d}, \qquad +P_d = \frac{d}{a + d}, +$$ + +and both parameter sets span exactly the same category probabilities (the identity is verified to machine precision in the package's test suite). What remains different is where the priors and random effects live: on the log scale of activations in the m3, on the logit scale of the process probabilities here. + +# Time-varying and other non-linear parameters + +Because MPT parameters are ordinary `bmm` parameters, they can be given non-linear predictor formulas. For example, a detection probability that grows exponentially with presentation time, given a numeric column `ptime` in the data: + +```{r, eval = FALSE} +fit_tv <- bmm( + bmf( + D ~ inv_logit(Dmax) * (1 - exp(-exp(rate) * ptime)), + Dmax ~ 1 + (1 | id), + rate ~ 1 + (1 | id) + ), + data = dat_ptime, + model = model_2htm +) +``` + +When a parameter receives a non-linear formula, the automatic link transformation is switched off for it: the user-supplied expression is responsible for keeping the parameter within (0, 1). The sub-parameters of the expression (`Dmax` and `rate` above) are estimated on the identity scale with `normal(0, 1)` default priors. Tip: when the random-effect standard deviations of such sub-parameters are close to zero, the posterior can develop a funnel geometry, and increasing `adapt_delta` (for example to 0.95) helps. + +# Importing existing model definitions + +Models written for other MPT software can be imported directly. `mpt_from_eqn()` reads EQN files (used by multiTree, MPTinR and TreeBUGS). Each line holds the tree, the response category and one branch expression; lines that end in the same category are summed. Parameter names containing underscores, which are common in published EQN files but not allowed in **brms**, are sanitized automatically, with the complete renaming map reported and stored in the `mpt_renaming` attribute: + +```{r} +eqn_file <- tempfile(fileext = ".eqn") +writeLines(c( + "6", + "target target_hit D_o", + "target target_hit (1-D_o)*g", + "target target_miss (1-D_o)*(1-g)", + "lure lure_fa (1-D_n)*g", + "lure lure_cr D_n", + "lure lure_cr (1-D_n)*(1-g)" +), eqn_file) + +model_eqn <- mpt_from_eqn( + eqn_file, + categories = c(target_hit = "old", lure_fa = "old", + target_miss = "new", lure_cr = "new"), + tree_id = "item_type" +) +model_eqn +``` + +Note the `categories` argument: classical EQN files often label response categories per tree (`hit`/`miss` for targets, `fa`/`cr` for lures), while in `bmm` every tree shares the same response categories, the actual response options. The mapping states which tree-specific labels refer to the same response. Lines whose category carries a `_` prefix have the prefix removed automatically. Two more arguments cover the remaining differences between files: a tree that has no branch for some of the shared categories declares them with `impossible`, as in `mpt_tree()`, and numeric tree labels are renamed with `tree_names`, e.g. `tree_names = c("1" = "target", "2" = "lure")`, because `bmm` tree names start with a letter. + +The file above starts with the equation count, the first line of classical EQN files, which `mpt_from_eqn()` skips. TreeBUGS and MPTinR skip the first line whatever it holds, while `mpt_from_eqn()` reads any other first line as an equation. A file without the count line therefore gives `bmm` one equation more than those programs read; the details of `?mpt_from_eqn` list the cases. + +Restrictions are written with the names used in the EQN file and renamed together with the parameters: + +```{r} +mpt_from_eqn( + eqn_file, + categories = c(target_hit = "old", lure_fa = "old", + target_miss = "new", lure_cr = "new"), + tree_id = "item_type", + restrictions = c("D_n = D_o", "g = 0.5") +) +``` + +`mpt_from_string()` reads MPTinR-style model definitions [@singmannKellen2013], where each line is a branch expression, blank lines separate trees, and inline comments name the response categories: + +```{r} +model_string <- " +D + (1 - D) * g # old +(1 - D) * (1 - g) # new + +(1 - D) * g # old +D + (1 - D) * (1 - g) # new +" +mpt_from_string( + model_string, + tree_names = c("target", "lure"), + tree_id = "item_type" +) +``` + +As a real-world example, the 12- and 18-parameter hierarchical binding models of @alhadhramiEtAl2025 (their unitization and hybrid models; 3 trees, 10 response categories, TreeBUGS latent-trait structure) can be imported from the model text in their Experiment 1 TreeBUGS scripts: saved to a file without its title line, it is in EQN layout for `mpt_from_eqn()`. Fitting them took two steps beyond the import: the nine guessing probabilities that the design fixes went into `restrictions`, and the wide frequency table was reshaped as shown under [Migrating from MPTinR or TreeBUGS]. With `links = "probit"` and correlated random effects across all parameters, the posterior means of the 18-parameter model lay inside the paper's 95% credible intervals for all 18 group-level estimates of Experiment 1 (Table 5). The paper reports none for the 12-parameter model, whose posterior means lay inside the 95% intervals of a TreeBUGS fit of the same data. + +# Migrating from MPTinR or TreeBUGS + +Importing the model structure is the easy part; what usually needs translating is the hierarchical structure, the priors, and the reporting conventions. The table maps the functions and arguments of MPTinR [@singmannKellen2013] and TreeBUGS [@heckEtAl2018] onto their `bmm` counterparts. + +| MPTinR / TreeBUGS | `bmm` | +|---|---| +| model file (EQN, MPTinR text) | `mpt_from_eqn()`, `mpt_from_string()`, or `mpt_tree()` + `mpt()` | +| wide data, one row per participant and one column per EQN category (`genTraitMPT()$data`) | long data, one row per participant and tree, with a `tree_id` column and one count column per shared response category ([Data format and simulation]) | +| reshaping that wide data: named columns (TreeBUGS), or unnamed columns sorted by tree and label (MPTinR, `check.mpt()$eqn.order.categories`) | `mpt_data_reshape(wide, model)`, or `mpt_data_reshape(wide, model, columns = "mptinr")` | +| `restrictions = list("Dn = Do", "g = 0.5")` (TreeBUGS; MPTinR takes the same list, or a file name, as `restrictions.filename`) | `mpt(..., restrictions = c("Dn = Do", "g = 0.5"))`; a fit-time constant is `bmf(g = 0.5)`. TreeBUGS' `FE` has no restriction counterpart: leave the random effect out of that parameter's formula | +| `check.mpt()` | `print(model)`: parameters, trees, restrictions, the identifiability bound and a Jacobian rank check ([Specifying an MPT model]) | +| `fit.mpt()`, `betaMPT()`, `traitMPT()` | `bmm()` with the random-effect structure written in the formula (see below) | +| `traitMPT(..., predStructure = )`, `covData` | predictors in the parameter formulas, e.g. `D ~ 0 + cond + (1 | id)` | +| `transformedParameters` | `brms::hypothesis()`, or arithmetic on the posterior draws | +| `getParam(fit, "theta")` | `native_parameters(fit, pars = "D", summary = TRUE)` | +| `getParam(fit, "mean")` | `native_parameters(fit, re_formula = NA, summary = TRUE)`: the probability of the median participant, as in TreeBUGS ([Population-level probabilities]) | +| `probitInverse()` | `?native_parameters`, section "Group-level effects" | +| `genTraitMPT()`, `genBetaMPT()` | `rmpt()` inside a small simulation function ([Data format and simulation]) | +| `PPP()` | `ppp()`: $T_1$ and $T_2$ with their posterior predictive p-values, and with `level = "participant"` the $T_1$ of each participant ([Goodness of fit: T1 and T2]); `pp_check()` for the plot | +| `WAIC()`, `select.mpt()` | `loo()`, `waic()` and `bayes_factor()` from **brms**; for aggregated counts see [Model comparison] before using any of them | +| `testHetChi()`, `testHetPerm()` | these test raw frequencies before any model is fitted: run them in TreeBUGS on the wide data. In `bmm`, compare models with and without person effects ([Model comparison]) | +| `BayesFactorMPT()`, `BayesFactorSlope()` | `brms::bayes_factor()` on fits stored with `save_pars = save_pars(all = TRUE)` ([Model comparison]) | + +**The data change from wide to long.** A frame in the layout of `genTraitMPT()`'s `data`, one row per participant and one column per EQN category, becomes one row per participant and tree, with the counts under the shared response categories. `mpt_data_reshape()` does the reshape for the imported model, matching the columns by the category labels of the EQN file; columns that are no counts, such as a between-participant factor, are carried to every row, and `id = ` names the participant column (without it, the rows are numbered): + +```{r} +wide <- data.frame(target_hit = c(40, 35), target_miss = c(10, 15), + lure_fa = c(8, 12), lure_cr = c(42, 38)) +long <- mpt_data_reshape(wide, model_eqn) +long +``` + +MPTinR data for an EQN model are unnamed. Their columns are sorted by tree and then by category label, the order of `check.mpt()$eqn.order.categories`; `columns = "mptinr"` reads them in that order. The row sums are the trials per row, but they cannot show swapped columns, so compare one participant's long rows with the wide row and run `bmm_data_check()` before fitting (the identifiability warning it prints is the single-test 2HTM result discussed above): + +```{r} +wide_mptinr <- unname(as.matrix(wide[c("lure_cr", "lure_fa", "target_hit", "target_miss")])) +identical(mpt_data_reshape(wide_mptinr, model_eqn, columns = "mptinr"), long) +bmm_data_check(bmf(Do ~ 1, Dn ~ 1, g ~ 1), long, model_eqn) +``` + +For a model from `mpt_from_string()` or `mpt()`, name the columns `.` instead, e.g. `target.old`. + +**The random-effect structure is written as a formula, not chosen from a menu.** The two families of hierarchical MPT models map onto the `|` and `||` operators, and the choice is per parameter: + +- `D ~ 0 + condition + (1 | participant:condition)`: independent person effects per condition with one shared standard deviation. +- `D ~ 0 + condition + (0 + condition || participant)`: one standard deviation per condition, uncorrelated across conditions. The beta-MPT of TreeBUGS also keeps person effects independent, but across parameters and with a beta distribution on the probability scale; `bmm` has no exact counterpart to it. +- `D ~ 0 + condition + (0 + condition | participant)`: one standard deviation per condition plus the full correlation matrix, that is, the latent-trait approach of @klauer2010. + +The last of these, combined with `links = "probit"`, is the model TreeBUGS fits with `traitMPT()`. Correlations *across* parameters, the defining feature of the latent-trait model, are requested with the `|p|` syntax of **brms**: `D ~ 1 + (1 |p| id)` and `g ~ 1 + (1 |p| id)` share one correlation matrix. Nothing forces `D` and `g` to share a random-effect structure, and parameters that only exist in one condition can simply take `(1 | participant)`. + +**Priors need converting from precision to standard deviation.** JAGS parameterizes the normal distribution by precision, Stan by standard deviation, so a JAGS `dnorm(0, 0.1)` on a latent-scale parameter becomes `normal(0, 3.1623)` in `bmm` ($1/\sqrt{0.1}$). Gamma priors carry over unchanged: JAGS `dgamma(1, 0.1)` on a random-effect standard deviation is `gamma(1, 0.1)` here, since both use the shape/rate parameterization. + +**Group-level estimates are reported on the latent scale.** TreeBUGS reports group means as probabilities, namely the probability of the median participant. To get there from a `bmm` fit, use `native_parameters(fit, re_formula = NA, summary = TRUE)`, as shown under [Population-level probabilities]; it applies the inverse link to the linear predictor, which also covers the probit link and models with predictors. + +# References