Skip to content

Bayesian Analysis

metasignal.sdtbayes fits hierarchical Bayesian metacognition models and analyses posteriors with ArviZ. Two backends are in use, depending on the approach:

  • cmdstanpy (recommended, no R required) — the primary backend, used by fit_meta_formula and the approaches redirected to it below.
  • brms via brmspy (requires R) — used by the plain-regression approaches below.

See Tutorial 7 — Bayesian Hierarchical Meta-d' for a hands-on walkthrough using the canonical HMeta-d dataset and a side-by-side comparison with metadpy.

Installation

The Bayesian submodule is an optional extra — install it alongside the base package:

pip install "metasignal[sdtbayes] @ git+https://github.com/saurabhr/metasignal.git"  # from GitHub
pip install ".[sdtbayes]"       # from a local clone

For the cmdstanpy-backed approaches (3, 4, and Approach 10's full path), install CmdStan once:

import cmdstanpy
cmdstanpy.install_cmdstan()

For the brms-backed approaches (1, 2, 7, 10's two-stage path, 12), run the one-time R/brms runtime setup instead (downloads Stan and installs brms, ~5 min):

from metasignal.sdtbayes import setup_runtime
setup_runtime()

Modelling approaches

# Approach Backend Status Primary function(s) Description
1 Ordered logistic brms ✅ working fit_hierarchical_metad Trial-level; confidence ratings as cumulative outcome. Fast and extensible.
2 Two-stage Bayesian brms ✅ working fit_two_stage_group MLE per participant → Bayesian on log M-ratio. Fast Stage 2; no Stage-1 uncertainty propagation.
3 Full HMeta-d cmdstanpy ✅ working fit_full_metad Single-stage port of Fleming (2017), delegating to fit_meta_formula. Numerically stable log-space likelihood; most accurate.
4 Subject-level Bayesian cmdstanpy ✅ working fit_subject_level Single-participant model. Matches metadpy.hmetad(nR_S1, nR_S2) API; reproduces metadpy's reference values to within MCMC noise.
7 Beta AUC brms ✅ working fit_beta_auc_group Non-parametric alternative to meta-d'. Models Type-2 AUC directly with a Beta likelihood; no Gaussian SDT assumption.
10 Meta-regression brms + cmdstanpy ✅ working fit_two_stage_regression, fit_full_metad_regression Regresses log M-ratio on participant-level covariates via either the two-stage (brms) or full hierarchical (cmdstanpy) path.
12 Within-subject comparison brms ✅ working fit_within_subject_comparison Paired model for within-subject designs: same participants in two conditions; participant random intercepts absorb stable individual differences.

Approaches 5 (robust HMeta-d), 6 (variational inference), 8 (Gaussian mixture), 9 (bivariate hierarchical), and 11 (state-space) were removed — each was permanently blocked by an upstream brmspy/rpy2 limitation with no working implementation. See the roadmap for status.

Approaches 1–3 and 7 have matching two-group comparison counterparts: fit_group_comparison, fit_two_stage_comparison, fit_full_metad_comparison, fit_beta_auc_comparison.

API Reference

Approach 1 — Ordered Logistic

fit_hierarchical_metad(participants, n_ratings, chains=4, n_iter=2000, warmup=1000, seed=42, items=None, **kwargs)

Fit a hierarchical ordered-logistic model of confidence on accuracy.

Uses an ordered logistic (cumulative) regression where confidence ratings are the outcome and correct is the key predictor. Random slopes per participant allow metacognitive discrimination to vary across the group.

The population-level coefficient b_correct is a monotonically-related index of group-mean metacognitive discrimination on the model's own log-odds scale; the random slopes in r_participant capture individual deviations from the group mean. b_correct is not meta-d' and is not on the same scale as Fleming's (2017) or Maniscalco & Lau's (2012) meta-d' — see the module docstring for why. Use :func:fit_full_metad or :func:fit_subject_level for an estimate on that scale.

Parameters:

Name Type Description Default
participants list[tuple[ndarray, ndarray, ndarray]]

List of (stim, resp, conf) tuples, one per participant. stim and resp are binary (0/1); conf is an integer from 1 to n_ratings.

 required
n_ratings int

Number of confidence rating categories.

 required
chains int

Number of MCMC chains (default 4).

 4
n_iter int

Total iterations per chain including warmup (default 2000).

 2000
warmup int

Warmup (burn-in) iterations per chain (default 1000).

 1000
seed int

Random seed for reproducibility (default 42).

 42
items list[ndarray] | None

Optional list of item-ID arrays, one per participant (same length as participants). Each array gives the stimulus item identity for every trial of that participant. When supplied, an item-level random intercept (1 | item) is added, producing a cross-classified model that separates participant-level metacognition from item-level difficulty effects.

 None
**kwargs Any

Additional arguments forwarded to brmspy.brms.brm.

 {}

Returns:

Type Description
Any

FitResult with:
Any
  • .idata — ArviZ InferenceData object containing posterior samples, log-likelihood, and prior predictive (if requested).
Any
  • .r — Lightweight R object handle for downstream brms calls.

Raises:

Type Description
ImportError

If brmspy is not installed.
ValueError

If items is provided but has a different length from participants.

Example::

import numpy as np
from metasignal.sdtbayes import fit_hierarchical_metad, posterior_summary

rng = np.random.default_rng(0)
N, n_items = 20, 50
participants = [
    (rng.integers(0, 2, 200), rng.integers(0, 2, 200), rng.integers(1, 5, 200))
    for _ in range(N)
]

# Without item effects
fit = fit_hierarchical_metad(participants, n_ratings=4)

# With crossed item random effects
item_ids = [rng.integers(0, n_items, 200) for _ in range(N)]
fit_crossed = fit_hierarchical_metad(participants, n_ratings=4, items=item_ids)
print(posterior_summary(fit_crossed, var_names=["b_correct", "sd_item__Intercept"]))
Source code in src/metasignal/sdtbayes/hierarchical.py 
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
def fit_hierarchical_metad(
    participants: list[tuple[np.ndarray, np.ndarray, np.ndarray]],
    n_ratings: int,
    chains: int = 4,
    n_iter: int = 2000,
    warmup: int = 1000,
    seed: int = 42,
    items: list[np.ndarray] | None = None,
    **kwargs: Any,
) -> Any:
    """Fit a hierarchical ordered-logistic model of confidence on accuracy.

    Uses an ordered logistic (cumulative) regression where confidence
    ratings are the outcome and ``correct`` is the key predictor.
    Random slopes per participant allow metacognitive discrimination to
    vary across the group.

    The population-level coefficient ``b_correct`` is a monotonically-related
    index of group-mean metacognitive discrimination on the model's own
    log-odds scale; the random slopes in ``r_participant`` capture individual
    deviations from the group mean. ``b_correct`` is not meta-d' and is not
    on the same scale as Fleming's (2017) or Maniscalco & Lau's (2012)
    meta-d' — see the module docstring for why. Use :func:`fit_full_metad`
    or :func:`fit_subject_level` for an estimate on that scale.

    Args:
        participants: List of ``(stim, resp, conf)`` tuples, one per participant.
            ``stim`` and ``resp`` are binary (0/1); ``conf`` is an integer
            from 1 to ``n_ratings``.
        n_ratings: Number of confidence rating categories.
        chains: Number of MCMC chains (default 4).
        n_iter: Total iterations per chain including warmup (default 2000).
        warmup: Warmup (burn-in) iterations per chain (default 1000).
        seed: Random seed for reproducibility (default 42).
        items: Optional list of item-ID arrays, one per participant (same
            length as ``participants``).  Each array gives the stimulus item
            identity for every trial of that participant.  When supplied, an
            item-level random intercept ``(1 | item)`` is added, producing a
            cross-classified model that separates participant-level
            metacognition from item-level difficulty effects.
        **kwargs: Additional arguments forwarded to ``brmspy.brms.brm``.

    Returns:
        ``FitResult`` with:

        - ``.idata`` — ArviZ ``InferenceData`` object containing posterior
            samples, log-likelihood, and prior predictive (if requested).
        - ``.r`` — Lightweight R object handle for downstream brms calls.

    Raises:
        ImportError: If ``brmspy`` is not installed.
        ValueError: If ``items`` is provided but has a different length from
            ``participants``.

    Example::

        import numpy as np
        from metasignal.sdtbayes import fit_hierarchical_metad, posterior_summary

        rng = np.random.default_rng(0)
        N, n_items = 20, 50
        participants = [
            (rng.integers(0, 2, 200), rng.integers(0, 2, 200), rng.integers(1, 5, 200))
            for _ in range(N)
        ]

        # Without item effects
        fit = fit_hierarchical_metad(participants, n_ratings=4)

        # With crossed item random effects
        item_ids = [rng.integers(0, n_items, 200) for _ in range(N)]
        fit_crossed = fit_hierarchical_metad(participants, n_ratings=4, items=item_ids)
        print(posterior_summary(fit_crossed, var_names=["b_correct", "sd_item__Intercept"]))
    """
    brms = require_brms()

    if items is not None and len(items) != len(participants):
        msg = (
            f"'items' has {len(items)} entries but 'participants' has "
            f"{len(participants)}.  They must have the same length."
        )
        raise ValueError(msg)

    df = _trials_to_dataframe(participants, n_ratings, items=items)

    formula_str = "conf ~ correct + (correct | participant)"
    if items is not None:
        formula_str += " + (1 | item)"
    formula = brms.bf(formula_str)

    priors = [
        brms.prior("normal(0, 2)", class_="b"),
        brms.prior("normal(0, 2)", class_="Intercept"),
        brms.prior("exponential(1)", class_="sd"),
        brms.prior("lkj(2)", class_="cor"),
    ]

    _result = brms.brm(
        formula=formula,
        data=df,
        family="cumulative",
        priors=priors,
        chains=chains,
        iter=n_iter,
        warmup=warmup,
        seed=seed,
        **kwargs,
    )
    return FitResult(idata=_result.idata, r=_result.r)

fit_group_comparison(group_a, group_b, n_ratings, chains=4, n_iter=2000, warmup=1000, seed=42, items_a=None, items_b=None, **kwargs)

Bayesian comparison of metacognition between two groups.

Extends the hierarchical ordered-logistic model with a group predictor and a correct × group interaction. The interaction coefficient b_correct:group1 represents the difference in metacognitive discrimination between groups on the model's own log-odds scale (not meta-d' units — see the module docstring), with a full posterior distribution.

Parameters:

Name Type Description Default
group_a list[tuple[ndarray, ndarray, ndarray]]

Participants in group A — list of (stim, resp, conf) tuples.

 required
group_b list[tuple[ndarray, ndarray, ndarray]]

Participants in group B — list of (stim, resp, conf) tuples.

 required
n_ratings int

Number of confidence rating categories.

 required
chains int

Number of MCMC chains (default 4).

 4
n_iter int

Total iterations per chain including warmup (default 2000).

 2000
warmup int

Warmup iterations per chain (default 1000).

 1000
seed int

Random seed (default 42).

 42
items_a list[ndarray] | None

Optional item-ID arrays for group A (parallel to group_a). When provided alongside items_b, adds (1 | item) to the formula for crossed stimulus random effects.

 None
items_b list[ndarray] | None

Optional item-ID arrays for group B (parallel to group_b).

 None
**kwargs Any

Additional arguments forwarded to brmspy.brms.brm.

 {}

Returns:

Type Description
Any

FitResult — same structure as :func:fit_hierarchical_metad.
Any

The parameter of interest is b_correct:group1 (group B minus
Any

group A difference in metacognitive discrimination, log-odds scale).

Raises:

Type Description
ImportError

If brmspy is not installed.
ValueError

If items_a / items_b lengths do not match groups, or if either group has fewer than 3 participants.

Example::

fit = fit_group_comparison(healthy, patient, n_ratings=4)
# Posterior probability that group B has lower metacognitive discrimination than group A:
import arviz as az
post = az.extract(fit.idata)["b_correct:group1"].values
print(f"P(group B < group A): {(post < 0).mean():.3f}")
Source code in src/metasignal/sdtbayes/hierarchical.py 
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
def fit_group_comparison(
    group_a: list[tuple[np.ndarray, np.ndarray, np.ndarray]],
    group_b: list[tuple[np.ndarray, np.ndarray, np.ndarray]],
    n_ratings: int,
    chains: int = 4,
    n_iter: int = 2000,
    warmup: int = 1000,
    seed: int = 42,
    items_a: list[np.ndarray] | None = None,
    items_b: list[np.ndarray] | None = None,
    **kwargs: Any,
) -> Any:
    """Bayesian comparison of metacognition between two groups.

    Extends the hierarchical ordered-logistic model with a ``group``
    predictor and a ``correct × group`` interaction.  The interaction
    coefficient ``b_correct:group1`` represents the difference in
    metacognitive discrimination between groups on the model's own
    log-odds scale (not meta-d' units — see the module docstring), with a
    full posterior distribution.

    Args:
        group_a: Participants in group A — list of ``(stim, resp, conf)`` tuples.
        group_b: Participants in group B — list of ``(stim, resp, conf)`` tuples.
        n_ratings: Number of confidence rating categories.
        chains: Number of MCMC chains (default 4).
        n_iter: Total iterations per chain including warmup (default 2000).
        warmup: Warmup iterations per chain (default 1000).
        seed: Random seed (default 42).
        items_a: Optional item-ID arrays for group A (parallel to ``group_a``).
            When provided alongside ``items_b``, adds ``(1 | item)`` to the
            formula for crossed stimulus random effects.
        items_b: Optional item-ID arrays for group B (parallel to ``group_b``).
        **kwargs: Additional arguments forwarded to ``brmspy.brms.brm``.

    Returns:
        ``FitResult`` — same structure as :func:`fit_hierarchical_metad`.
        The parameter of interest is ``b_correct:group1`` (group B minus
        group A difference in metacognitive discrimination, log-odds scale).

    Raises:
        ImportError: If ``brmspy`` is not installed.
        ValueError: If ``items_a`` / ``items_b`` lengths do not match groups,
            or if either group has fewer than 3 participants.

    Example::

        fit = fit_group_comparison(healthy, patient, n_ratings=4)
        # Posterior probability that group B has lower metacognitive discrimination than group A:
        import arviz as az
        post = az.extract(fit.idata)["b_correct:group1"].values
        print(f"P(group B < group A): {(post < 0).mean():.3f}")
    """
    brms = require_brms()

    if items_a is not None and len(items_a) != len(group_a):
        msg = f"'items_a' has {len(items_a)} entries but 'group_a' has {len(group_a)}."
        raise ValueError(msg)
    if items_b is not None and len(items_b) != len(group_b):
        msg = f"'items_b' has {len(items_b)} entries but 'group_b' has {len(group_b)}."
        raise ValueError(msg)
    for participants, label in ((group_a, "A"), (group_b, "B")):
        if len(participants) < 3:
            msg = (
                f"Group {label}: only {len(participants)} participants — need at "
                "least 3 to fit a hierarchical group comparison."
            )
            raise ValueError(msg)

    import pandas as pd
    df_a = _trials_to_dataframe(group_a, n_ratings, group_label=0, pid_prefix="A", items=items_a)
    df_b = _trials_to_dataframe(group_b, n_ratings, group_label=1, pid_prefix="B", items=items_b)
    df = pd.concat([df_a, df_b], ignore_index=True)
    df["group"] = df["group"].astype("category")

    use_items = (items_a is not None) and (items_b is not None)
    formula_str = "conf ~ correct * group + (correct | participant)"
    if use_items:
        formula_str += " + (1 | item)"
    formula = brms.bf(formula_str)

    priors = [
        brms.prior("normal(0, 2)", class_="b"),
        brms.prior("normal(0, 2)", class_="Intercept"),
        brms.prior("exponential(1)", class_="sd"),
        brms.prior("lkj(2)", class_="cor"),
    ]

    _result = brms.brm(
        formula=formula,
        data=df,
        family="cumulative",
        priors=priors,
        chains=chains,
        iter=n_iter,
        warmup=warmup,
        seed=seed,
        **kwargs,
    )
    return FitResult(idata=_result.idata, r=_result.r)

Approach 2 — Two-Stage Bayesian

fit_two_stage_group(participants, n_ratings, chains=4, n_iter=2000, warmup=1000, seed=42, **kwargs)

Two-stage hierarchical Bayesian model for group-level M-ratio.

Runs MLE estimation per participant (Stage 1) then fits a hierarchical Bayesian model over log M-ratio across participants (Stage 2). The key group-level parameters are:

  • b_Intercept — posterior mean log M-ratio (exp gives group mean M-ratio)
  • sigma — between-subject SD on the log scale

Parameters:

Name Type Description Default
participants list[tuple[ndarray, ndarray, ndarray]]

List of (stim, resp, conf) tuples, one per participant.

 required
n_ratings int

Number of confidence rating categories.

 required
chains int

MCMC chains (default 4).

 4
n_iter int

Total iterations per chain including warmup (default 2000).

 2000
warmup int

Warmup iterations (default 1000).

 1000
seed int

Random seed (default 42).

 42
**kwargs Any

Forwarded to brmspy.brms.brm.

 {}

Returns:

Type Description
Any

Tuple of FitResult and pd.DataFrame.
'pd.DataFrame'

The FitResult has .idata (ArviZ InferenceData) and .r
tuple[Any, 'pd.DataFrame']

(R handle). The key population-level parameter is b_Intercept
tuple[Any, 'pd.DataFrame']

(group mean log M-ratio).
tuple[Any, 'pd.DataFrame']

The pd.DataFrame contains per-participant Stage-1 MLE estimates
tuple[Any, 'pd.DataFrame']

(columns: participant, dprime, c, meta_da, da,
tuple[Any, 'pd.DataFrame']

m_ratio, log_m_ratio).

Raises:

Type Description
ImportError

If brmspy is not installed.

Example::

import numpy as np
from metasignal.sdtbayes import fit_two_stage_group, posterior_summary

rng = np.random.default_rng(0)
participants = [
    (rng.integers(0, 2, 200), rng.integers(0, 2, 200), rng.integers(1, 5, 200))
    for _ in range(20)
]

fit, mle_df = fit_two_stage_group(participants, n_ratings=4)
print(mle_df[["participant", "m_ratio", "log_m_ratio"]])
print(posterior_summary(fit, var_names=["b_Intercept", "sigma"]))
Source code in src/metasignal/sdtbayes/two_stage.py 
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
def fit_two_stage_group(
    participants: list[tuple[np.ndarray, np.ndarray, np.ndarray]],
    n_ratings: int,
    chains: int = 4,
    n_iter: int = 2000,
    warmup: int = 1000,
    seed: int = 42,
    **kwargs: Any,
) -> tuple[Any, "pd.DataFrame"]:
    """Two-stage hierarchical Bayesian model for group-level M-ratio.

    Runs MLE estimation per participant (Stage 1) then fits a hierarchical
    Bayesian model over log M-ratio across participants (Stage 2).  The key
    group-level parameters are:

    - ``b_Intercept`` — posterior mean log M-ratio (exp gives group mean M-ratio)
    - ``sigma`` — between-subject SD on the log scale

    Args:
        participants: List of ``(stim, resp, conf)`` tuples, one per participant.
        n_ratings: Number of confidence rating categories.
        chains: MCMC chains (default 4).
        n_iter: Total iterations per chain including warmup (default 2000).
        warmup: Warmup iterations (default 1000).
        seed: Random seed (default 42).
        **kwargs: Forwarded to ``brmspy.brms.brm``.

    Returns:
        Tuple of ``FitResult`` and ``pd.DataFrame``.

        The ``FitResult`` has ``.idata`` (ArviZ InferenceData) and ``.r``
        (R handle). The key population-level parameter is ``b_Intercept``
        (group mean log M-ratio).

        The ``pd.DataFrame`` contains per-participant Stage-1 MLE estimates
        (columns: ``participant``, ``dprime``, ``c``, ``meta_da``, ``da``,
        ``m_ratio``, ``log_m_ratio``).

    Raises:
        ImportError: If ``brmspy`` is not installed.

    Example::

        import numpy as np
        from metasignal.sdtbayes import fit_two_stage_group, posterior_summary

        rng = np.random.default_rng(0)
        participants = [
            (rng.integers(0, 2, 200), rng.integers(0, 2, 200), rng.integers(1, 5, 200))
            for _ in range(20)
        ]

        fit, mle_df = fit_two_stage_group(participants, n_ratings=4)
        print(mle_df[["participant", "m_ratio", "log_m_ratio"]])
        print(posterior_summary(fit, var_names=["b_Intercept", "sigma"]))
    """
    brms = require_brms()

    mle_df = _compute_participant_estimates(participants, n_ratings)

    valid = mle_df.dropna(subset=["log_m_ratio"])
    if len(valid) < 3:
        msg = f"Only {len(valid)} participants have valid MLE estimates — need at least 3."
        raise ValueError(msg)

    formula = brms.bf("log_m_ratio ~ 1")
    priors = [
        brms.prior("normal(0, 1)", class_="Intercept"),
        brms.prior("exponential(1)", class_="sigma"),
    ]

    fit = brms.brm(
        formula=formula,
        data=valid,
        family="student",
        priors=priors,
        chains=chains,
        iter=n_iter,
        warmup=warmup,
        seed=seed,
        **kwargs,
    )
    return FitResult(idata=fit.idata, r=fit.r), mle_df

fit_two_stage_comparison(group_a, group_b, n_ratings, chains=4, n_iter=2000, warmup=1000, seed=42, **kwargs)

Two-stage Bayesian comparison of M-ratio between two groups.

MLE is run per participant (Stage 1); a Bayesian regression with a group predictor is fit on log M-ratio (Stage 2). The coefficient b_group1 is the posterior difference in log M-ratio (group B − group A).

Parameters:

Name Type Description Default
group_a list[tuple[ndarray, ndarray, ndarray]]

Participants in group A.

 required
group_b list[tuple[ndarray, ndarray, ndarray]]

Participants in group B.

 required
n_ratings int

Number of confidence rating categories.

 required
chains int

MCMC chains (default 4).

 4
n_iter int

Total iterations per chain including warmup (default 2000).

 2000
warmup int

Warmup iterations (default 1000).

 1000
seed int

Random seed (default 42).

 42
**kwargs Any

Forwarded to brmspy.brms.brm.

 {}

Returns:

Type Description
Any

Tuple of FitResult and per-participant MLE DataFrame (with a
'pd.DataFrame'

group column, 0 = A, 1 = B).

Raises:

Type Description
ImportError

If brmspy is not installed.

Example::

fit, mle_df = fit_two_stage_comparison(healthy, patient, n_ratings=4)

# Posterior probability that group B has lower M-ratio than group A
import arviz as az
post = az.extract(fit.idata)["b_group1"].values
print(f"P(group B < group A): {(post < 0).mean():.3f}")
Source code in src/metasignal/sdtbayes/two_stage.py 
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
def fit_two_stage_comparison(
    group_a: list[tuple[np.ndarray, np.ndarray, np.ndarray]],
    group_b: list[tuple[np.ndarray, np.ndarray, np.ndarray]],
    n_ratings: int,
    chains: int = 4,
    n_iter: int = 2000,
    warmup: int = 1000,
    seed: int = 42,
    **kwargs: Any,
) -> tuple[Any, "pd.DataFrame"]:
    """Two-stage Bayesian comparison of M-ratio between two groups.

    MLE is run per participant (Stage 1); a Bayesian regression with a
    ``group`` predictor is fit on log M-ratio (Stage 2).  The coefficient
    ``b_group1`` is the posterior difference in log M-ratio (group B − group A).

    Args:
        group_a: Participants in group A.
        group_b: Participants in group B.
        n_ratings: Number of confidence rating categories.
        chains: MCMC chains (default 4).
        n_iter: Total iterations per chain including warmup (default 2000).
        warmup: Warmup iterations (default 1000).
        seed: Random seed (default 42).
        **kwargs: Forwarded to ``brmspy.brms.brm``.

    Returns:
        Tuple of ``FitResult`` and per-participant MLE DataFrame (with a
        ``group`` column, 0 = A, 1 = B).

    Raises:
        ImportError: If ``brmspy`` is not installed.

    Example::

        fit, mle_df = fit_two_stage_comparison(healthy, patient, n_ratings=4)

        # Posterior probability that group B has lower M-ratio than group A
        import arviz as az
        post = az.extract(fit.idata)["b_group1"].values
        print(f"P(group B < group A): {(post < 0).mean():.3f}")
    """
    brms = require_brms()

    import pandas as pd

    df_a = _compute_participant_estimates(group_a, n_ratings)
    df_a["group"] = 0
    df_b = _compute_participant_estimates(group_b, n_ratings)
    df_b["group"] = 1
    # Re-index participant IDs so they are unique across groups
    df_b["participant"] = df_b["participant"] + len(group_a)
    mle_df = pd.concat([df_a, df_b], ignore_index=True)
    mle_df["group"] = mle_df["group"].astype("category")

    valid = mle_df.dropna(subset=["log_m_ratio"])

    for g, label in ((0, "A"), (1, "B")):
        n_valid = int((valid["group"] == g).sum())
        if n_valid < 3:
            total = int((mle_df["group"] == g).sum())
            msg = (
                f"Group {label}: only {n_valid} of {total} participants have valid "
                "MLE estimates — need at least 3."
            )
            raise ValueError(msg)

    formula = brms.bf("log_m_ratio ~ group")
    priors = [
        brms.prior("normal(0, 1)", class_="Intercept"),
        brms.prior("normal(0, 1)", class_="b"),
        brms.prior("exponential(1)", class_="sigma"),
    ]

    fit = brms.brm(
        formula=formula,
        data=valid,
        family="student",
        priors=priors,
        chains=chains,
        iter=n_iter,
        warmup=warmup,
        seed=seed,
        **kwargs,
    )
    return FitResult(idata=fit.idata, r=fit.r), mle_df

Approach 3 — Full HMeta-d

fit_full_metad(participants, n_ratings, chains=4, n_iter=2000, warmup=1000, seed=42, tol=1e-07, **kwargs)

Full hierarchical meta-d' (HMeta-d) — Fleming (2017) model.

This delegates to :func:metasignal.sdtbayes.fit_meta_formula with parameterization="mratio" and backend="stan" (cmdstanpy), which implements the same Fleming (2017) HMeta-d model with a numerically stable log-space likelihood, hard-by-construction ordered Type-2 criteria, and a Type-1 binomial likelihood tying d'/c' directly to the observed hit/false-alarm totals.

This function previously injected custom Stan code into brms via brmspy; that path is no longer used because the brmspy/rpy2 bridge cannot reliably convert a list of multiple stanvar() objects (a block != "data" stanvar round-trips to Python as a differently-shaped object than a data-carrying one, so a combined list fails rpy2's homogeneous-type dispatch — this is an upstream brmspy limitation, not fixable from calling code).

The group-level parameter of primary interest is alpha_logMratio (aliased as mu_logMratio), the posterior mean log M-ratio. exp(alpha_logMratio) gives the group mean M-ratio.

Key parameters in the posterior

  • alpha_logMratio / mu_logMratio — group mean log M-ratio
  • sigma_logMratio — between-subject SD on the log M-ratio scale
  • Mratio[s] — per-subject M-ratio (transformed parameter)
  • meta_d[s] — per-subject meta-d' (transformed parameter)
  • d1[s], c1[s] — per-subject Type-1 d' and criterion
  • mu_d1, mu_c1 — group mean Type-1 parameters

Parameters:

Name Type Description Default
participants list[tuple[ndarray, ndarray, ndarray]]

List of (stim, resp, conf) tuples, one per participant.

 required
n_ratings int

Number of confidence rating categories.

 required
chains int

MCMC chains (default 4).

 4
n_iter int

Total iterations per chain including warmup (default 2000).

 2000
warmup int

Warmup iterations (default 1000).

 1000
seed int

Random seed (default 42).

 42
tol float

Minimum probability floor for multinomial cells (default 1e-7).

 1e-07
**kwargs Any

Forwarded to :func:fit_meta_formula (e.g. cmdstan_path).

 {}

Returns:

Type Description
Any

FitResult with .idata (ArviZ InferenceData).

Raises:

Type Description
ImportError

If cmdstanpy is not installed.

Example::

import numpy as np
from metasignal.sdtbayes import fit_full_metad, posterior_summary

rng = np.random.default_rng(0)
participants = [
    (rng.integers(0, 2, 200), rng.integers(0, 2, 200), rng.integers(1, 5, 200))
    for _ in range(20)
]

fit = fit_full_metad(participants, n_ratings=4)
print(posterior_summary(fit, var_names=["mu_logMratio", "sigma_logMratio",
                                         "mu_d1", "Mratio"]))
Source code in src/metasignal/sdtbayes/full_metad.py 
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
def fit_full_metad(
    participants: list[tuple[np.ndarray, np.ndarray, np.ndarray]],
    n_ratings: int,
    chains: int = 4,
    n_iter: int = 2000,
    warmup: int = 1000,
    seed: int = 42,
    tol: float = 1e-7,
    **kwargs: Any,
) -> Any:
    """Full hierarchical meta-d' (HMeta-d) — Fleming (2017) model.

    This delegates to :func:`metasignal.sdtbayes.fit_meta_formula` with
    ``parameterization="mratio"`` and ``backend="stan"`` (cmdstanpy), which
    implements the same Fleming (2017) HMeta-d model with a numerically
    stable log-space likelihood, hard-by-construction ordered Type-2 criteria,
    and a Type-1 binomial likelihood tying d'/c' directly to the observed
    hit/false-alarm totals.

    This function previously injected custom Stan code into brms via
    ``brmspy``; that path is no longer used because the brmspy/rpy2 bridge
    cannot reliably convert a list of multiple ``stanvar()`` objects (a
    ``block != "data"`` stanvar round-trips to Python as a differently-shaped
    object than a data-carrying one, so a combined list fails rpy2's
    homogeneous-type dispatch — this is an upstream brmspy limitation, not
    fixable from calling code).

    The group-level parameter of primary interest is **``alpha_logMratio``**
    (aliased as ``mu_logMratio``), the posterior mean log M-ratio.
    ``exp(alpha_logMratio)`` gives the group mean M-ratio.

    Key parameters in the posterior
    --------------------------------
    - ``alpha_logMratio`` / ``mu_logMratio`` — group mean log M-ratio
    - ``sigma_logMratio`` — between-subject SD on the log M-ratio scale
    - ``Mratio[s]`` — per-subject M-ratio (transformed parameter)
    - ``meta_d[s]`` — per-subject meta-d' (transformed parameter)
    - ``d1[s]``, ``c1[s]`` — per-subject Type-1 d' and criterion
    - ``mu_d1``, ``mu_c1`` — group mean Type-1 parameters

    Args:
        participants: List of ``(stim, resp, conf)`` tuples, one per participant.
        n_ratings: Number of confidence rating categories.
        chains: MCMC chains (default 4).
        n_iter: Total iterations per chain including warmup (default 2000).
        warmup: Warmup iterations (default 1000).
        seed: Random seed (default 42).
        tol: Minimum probability floor for multinomial cells (default 1e-7).
        **kwargs: Forwarded to :func:`fit_meta_formula` (e.g. ``cmdstan_path``).

    Returns:
        ``FitResult`` with ``.idata`` (ArviZ InferenceData).

    Raises:
        ImportError: If ``cmdstanpy`` is not installed.

    Example::

        import numpy as np
        from metasignal.sdtbayes import fit_full_metad, posterior_summary

        rng = np.random.default_rng(0)
        participants = [
            (rng.integers(0, 2, 200), rng.integers(0, 2, 200), rng.integers(1, 5, 200))
            for _ in range(20)
        ]

        fit = fit_full_metad(participants, n_ratings=4)
        print(posterior_summary(fit, var_names=["mu_logMratio", "sigma_logMratio",
                                                 "mu_d1", "Mratio"]))
    """
    import pandas as pd

    from metasignal.sdtbayes.formula import fit_meta_formula

    pred_df = pd.DataFrame({"participant": range(len(participants))})
    return fit_meta_formula(
        participants=participants,
        n_ratings=n_ratings,
        formula="~ 1",
        data=pred_df,
        parameterization="mratio",
        backend="stan",
        chains=chains,
        n_iter=n_iter,
        warmup=warmup,
        seed=seed,
        tol=tol,
        **kwargs,
    )

fit_full_metad_comparison(group_a, group_b, n_ratings, chains=4, n_iter=2000, warmup=1000, seed=42, tol=1e-07, **kwargs)

Full HMeta-d comparison between two groups (Fleming 2017 model).

Delegates to :func:fit_meta_formula with a binary group covariate (0 = group A, 1 = group B) and parameterization="mratio". See :func:fit_full_metad for why this no longer uses the brmspy stanvar injection path.

Because the covariate is mean-centred internally (0/1 → -0.5/+0.5, spanning exactly one unit), the fitted slope beta_logMratio[1] is the full group B − group A difference in log M-ratio — the same quantity the old delta_logMratio parameter represented.

Parameters:

Name Type Description Default
group_a list[tuple[ndarray, ndarray, ndarray]]

Participants in group A.

 required
group_b list[tuple[ndarray, ndarray, ndarray]]

Participants in group B.

 required
n_ratings int

Number of confidence rating categories.

 required
chains int

MCMC chains (default 4).

 4
n_iter int

Total iterations per chain including warmup (default 2000).

 2000
warmup int

Warmup iterations (default 1000).

 1000
seed int

Random seed (default 42).

 42
tol float

Minimum probability floor (default 1e-7).

 1e-07
**kwargs Any

Forwarded to :func:fit_meta_formula.

 {}

Returns:

Type Description
Any

FitResult. Key parameters: alpha_logMratio (group A mean
Any

log M-ratio, at the covariate mean), beta_logMratio (length-1
Any

vector — group B − A difference).

Raises:

Type Description
ImportError

If cmdstanpy is not installed.

Example::

fit = fit_full_metad_comparison(healthy, patient, n_ratings=4)

import arviz as az
delta = az.extract(fit.idata)["beta_logMratio"].values[:, 0]
print(f"P(group B > group A): {(delta > 0).mean():.3f}")
Source code in src/metasignal/sdtbayes/full_metad.py 
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
def fit_full_metad_comparison(
    group_a: list[tuple[np.ndarray, np.ndarray, np.ndarray]],
    group_b: list[tuple[np.ndarray, np.ndarray, np.ndarray]],
    n_ratings: int,
    chains: int = 4,
    n_iter: int = 2000,
    warmup: int = 1000,
    seed: int = 42,
    tol: float = 1e-7,
    **kwargs: Any,
) -> Any:
    """Full HMeta-d comparison between two groups (Fleming 2017 model).

    Delegates to :func:`fit_meta_formula` with a binary ``group`` covariate
    (0 = group A, 1 = group B) and ``parameterization="mratio"``.  See
    :func:`fit_full_metad` for why this no longer uses the brmspy stanvar
    injection path.

    Because the covariate is mean-centred internally (0/1 → -0.5/+0.5,
    spanning exactly one unit), the fitted slope ``beta_logMratio[1]`` is
    the full group B − group A difference in log M-ratio — the same
    quantity the old ``delta_logMratio`` parameter represented.

    Args:
        group_a: Participants in group A.
        group_b: Participants in group B.
        n_ratings: Number of confidence rating categories.
        chains: MCMC chains (default 4).
        n_iter: Total iterations per chain including warmup (default 2000).
        warmup: Warmup iterations (default 1000).
        seed: Random seed (default 42).
        tol: Minimum probability floor (default 1e-7).
        **kwargs: Forwarded to :func:`fit_meta_formula`.

    Returns:
        ``FitResult``.  Key parameters: ``alpha_logMratio`` (group A mean
        log M-ratio, at the covariate mean), ``beta_logMratio`` (length-1
        vector — group B − A difference).

    Raises:
        ImportError: If ``cmdstanpy`` is not installed.

    Example::

        fit = fit_full_metad_comparison(healthy, patient, n_ratings=4)

        import arviz as az
        delta = az.extract(fit.idata)["beta_logMratio"].values[:, 0]
        print(f"P(group B > group A): {(delta > 0).mean():.3f}")
    """
    import pandas as pd

    from metasignal.sdtbayes.formula import fit_meta_formula

    participants = list(group_a) + list(group_b)
    pred_df = pd.DataFrame({
        "participant": range(len(participants)),
        "group": [0] * len(group_a) + [1] * len(group_b),
    })
    return fit_meta_formula(
        participants=participants,
        n_ratings=n_ratings,
        formula="~ group",
        data=pred_df,
        parameterization="mratio",
        backend="stan",
        chains=chains,
        n_iter=n_iter,
        warmup=warmup,
        seed=seed,
        tol=tol,
        **kwargs,
    )

Approach 4 — Subject-Level Bayesian (metadpy API)

fit_subject_level(nR_S1, nR_S2, chains=4, n_iter=2000, warmup=1000, seed=42, tol=1e-05, **kwargs)

Subject-level Bayesian meta-d' — matches metadpy's hmetad().

Estimates d1, c1, and meta_d for a single participant using the same fixed-prior model as metadpy's PyMC hmetad() implementation, but running via cmdstanpy (delegates to :func:metasignal.sdtbayes.fit_meta_formula with parameterization="mratio").

This function previously injected custom Stan code into brms via brmspy; that path is no longer used because the brmspy/rpy2 bridge cannot reliably convert a list of multiple stanvar() objects into R (see :func:metasignal.sdtbayes.fit_full_metad docstring for details).

The posterior for meta_d should reproduce metadpy's meta_d parameter to within MCMC noise. For the canonical dataset from Fleming (2017) / metadpy tutorial 3::

nR_S1 = [52, 32, 35, 37, 26, 12, 4, 2]
nR_S2 = [2, 5, 15, 22, 33, 38, 40, 45]

both implementations return meta_d ≈ 1.57 ± 0.20 and d1 ≈ 1.53 ± 0.14.

Parameters:

Name Type Description Default
nR_S1 ndarray

Rating distribution for S1 (noise) trials, length 2 * nratings. Elements [0..nR-1] are correct rejections (high→low confidence); elements [nR..2*nR-1] are false alarms (low→high confidence).

 required
nR_S2 ndarray

Rating distribution for S2 (signal) trials, same format. Elements [0..nR-1] are misses; elements [nR..2*nR-1] are hits.

 required
chains int

MCMC chains (default 4).

 4
n_iter int

Total iterations per chain including warmup (default 2000).

 2000
warmup int

Warmup iterations (default 1000).

 1000
seed int

Random seed (default 42).

 42
tol float

Minimum probability floor for multinomial cells (default 1e-5, matching metadpy's Tol).

 1e-05
**kwargs Any

Forwarded to :func:fit_meta_formula.

 {}

Returns:

Type Description
Any

FitResult with .idata (ArviZ InferenceData) containing
Any

posterior samples for d1, c1, meta_d, Mratio.

Raises:

Type Description
ImportError

If cmdstanpy is not installed.
ValueError

If nR_S1 and nR_S2 have different lengths or the length is not even.

Example::

import numpy as np
from metasignal.sdtbayes import fit_subject_level, posterior_summary

# Canonical dataset from Fleming (2017) and metadpy tutorial
nR_S1 = np.array([52, 32, 35, 37, 26, 12, 4, 2])
nR_S2 = np.array([2, 5, 15, 22, 33, 38, 40, 45])

fit = fit_subject_level(nR_S1, nR_S2)
print(posterior_summary(fit, var_names=["d1", "c1", "meta_d"]))

# Expected output (comparable to metadpy):
#   d1      mean ≈  1.53,  sd ≈ 0.14
#   c1      mean ≈ -0.01,  sd ≈ 0.07
#   meta_d  mean ≈  1.57,  sd ≈ 0.20
Source code in src/metasignal/sdtbayes/subject_level.py 
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
def fit_subject_level(  # pylint: disable=invalid-name
    nR_S1: np.ndarray,
    nR_S2: np.ndarray,
    chains: int = 4,
    n_iter: int = 2000,
    warmup: int = 1000,
    seed: int = 42,
    tol: float = 1e-5,
    **kwargs: Any,
) -> Any:
    """Subject-level Bayesian meta-d' — matches metadpy's ``hmetad()``.

    Estimates ``d1``, ``c1``, and ``meta_d`` for a single participant using
    the same fixed-prior model as metadpy's PyMC ``hmetad()`` implementation,
    but running via **cmdstanpy** (delegates to
    :func:`metasignal.sdtbayes.fit_meta_formula` with
    ``parameterization="mratio"``).

    This function previously injected custom Stan code into brms via
    ``brmspy``; that path is no longer used because the brmspy/rpy2 bridge
    cannot reliably convert a list of multiple ``stanvar()`` objects into R
    (see :func:`metasignal.sdtbayes.fit_full_metad` docstring for details).

    The posterior for ``meta_d`` should reproduce metadpy's ``meta_d``
    parameter to within MCMC noise.  For the canonical dataset from Fleming
    (2017) / metadpy tutorial 3::

        nR_S1 = [52, 32, 35, 37, 26, 12, 4, 2]
        nR_S2 = [2, 5, 15, 22, 33, 38, 40, 45]

    both implementations return ``meta_d ≈ 1.57 ± 0.20`` and
    ``d1 ≈ 1.53 ± 0.14``.

    Args:
        nR_S1: Rating distribution for S1 (noise) trials, length ``2 * nratings``.
            Elements ``[0..nR-1]`` are correct rejections (high→low confidence);
            elements ``[nR..2*nR-1]`` are false alarms (low→high confidence).
        nR_S2: Rating distribution for S2 (signal) trials, same format.
            Elements ``[0..nR-1]`` are misses; elements ``[nR..2*nR-1]`` are hits.
        chains: MCMC chains (default 4).
        n_iter: Total iterations per chain including warmup (default 2000).
        warmup: Warmup iterations (default 1000).
        seed: Random seed (default 42).
        tol: Minimum probability floor for multinomial cells (default 1e-5,
            matching metadpy's ``Tol``).
        **kwargs: Forwarded to :func:`fit_meta_formula`.

    Returns:
        ``FitResult`` with ``.idata`` (ArviZ ``InferenceData``) containing
        posterior samples for ``d1``, ``c1``, ``meta_d``, ``Mratio``.

    Raises:
        ImportError: If ``cmdstanpy`` is not installed.
        ValueError: If ``nR_S1`` and ``nR_S2`` have different lengths or the
            length is not even.

    Example::

        import numpy as np
        from metasignal.sdtbayes import fit_subject_level, posterior_summary

        # Canonical dataset from Fleming (2017) and metadpy tutorial
        nR_S1 = np.array([52, 32, 35, 37, 26, 12, 4, 2])
        nR_S2 = np.array([2, 5, 15, 22, 33, 38, 40, 45])

        fit = fit_subject_level(nR_S1, nR_S2)
        print(posterior_summary(fit, var_names=["d1", "c1", "meta_d"]))

        # Expected output (comparable to metadpy):
        #   d1      mean ≈  1.53,  sd ≈ 0.14
        #   c1      mean ≈ -0.01,  sd ≈ 0.07
        #   meta_d  mean ≈  1.57,  sd ≈ 0.20
    """
    nR_S1 = np.asarray(nR_S1, dtype=int)
    nR_S2 = np.asarray(nR_S2, dtype=int)

    if len(nR_S1) != len(nR_S2):
        raise ValueError("nR_S1 and nR_S2 must have the same length.")
    if len(nR_S1) % 2 != 0:
        raise ValueError("Length of nR_S1 must be even (= 2 * nratings).")

    import pandas as pd

    from metasignal.sdtbayes.formula import fit_meta_formula

    n_ratings = len(nR_S1) // 2

    # Expand counts back to trial-level (stim, resp, conf) arrays — exact
    # inverse of metasignal.stdpy.core.trials_to_counts.
    stims: list[int] = []
    resps: list[int] = []
    confs: list[int] = []
    for idx in range(n_ratings):
        conf = n_ratings - idx  # high confidence at idx=0
        n = int(nR_S1[idx])         # CR: S1 stim, S1 resp
        stims += [0] * n; resps += [0] * n; confs += [conf] * n
        n = int(nR_S2[idx])         # Miss: S2 stim, S1 resp
        stims += [1] * n; resps += [0] * n; confs += [conf] * n
    for idx in range(n_ratings):
        conf = idx + 1               # low confidence at idx=0
        n = int(nR_S1[n_ratings + idx])  # FA: S1 stim, S2 resp
        stims += [0] * n; resps += [1] * n; confs += [conf] * n
        n = int(nR_S2[n_ratings + idx])  # Hit: S2 stim, S2 resp
        stims += [1] * n; resps += [1] * n; confs += [conf] * n

    participants = [(np.array(stims), np.array(resps), np.array(confs))]
    pred_df = pd.DataFrame({"participant": [0]})
    return fit_meta_formula(
        participants=participants,
        n_ratings=n_ratings,
        formula="~ 1",
        data=pred_df,
        parameterization="mratio",
        backend="stan",
        chains=chains,
        n_iter=n_iter,
        warmup=warmup,
        seed=seed,
        tol=tol,
        **kwargs,
    )

Approach 7 — Beta Regression on Type-2 AUC

fit_beta_auc_group(participants, n_ratings, chains=4, n_iter=2000, warmup=1000, seed=42, **kwargs)

Hierarchical Beta regression on Type-2 AUC for a single group.

Computes AUC2 per participant (Stage 1) and models the group-level distribution using a Beta likelihood (Stage 2). Unlike meta-d' approaches this makes no equal-variance Gaussian SDT assumption — AUC2 is a non-parametric measure of metacognitive performance.

The group mean AUC2 is expit(b_Intercept).

Parameters:

Name Type Description Default
participants list[tuple[ndarray, ndarray, ndarray]]

List of (stim, resp, conf) tuples, one per participant.

 required
n_ratings int

Number of confidence rating categories.

 required
chains int

MCMC chains (default 4).

 4
n_iter int

Total iterations per chain including warmup (default 2000).

 2000
warmup int

Warmup iterations (default 1000).

 1000
seed int

Random seed (default 42).

 42
**kwargs Any

Forwarded to brmspy.brms.brm.

 {}

Returns:

Type Description
FitResult

FitResult with posterior for b_Intercept (logit-scale group mean
FitResult

AUC2) and phi (precision).

Raises:

Type Description
ImportError

If brmspy is not installed.
ValueError

If fewer than 3 participants have valid AUC2 estimates.

Example::

import numpy as np
from scipy.special import expit
import arviz as az
from metasignal.sdtbayes import fit_beta_auc_group

rng = np.random.default_rng(0)
participants = [
    (rng.integers(0, 2, 200), rng.integers(0, 2, 200), rng.integers(1, 5, 200))
    for _ in range(20)
]
fit = fit_beta_auc_group(participants, n_ratings=4)

intercept = az.extract(fit.idata)["b_Intercept"].values
print(f"Group mean AUC2: {expit(intercept).mean():.3f}")
Source code in src/metasignal/sdtbayes/beta_auc.py 
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
def fit_beta_auc_group(
    participants: list[tuple[np.ndarray, np.ndarray, np.ndarray]],
    n_ratings: int,
    chains: int = 4,
    n_iter: int = 2000,
    warmup: int = 1000,
    seed: int = 42,
    **kwargs: Any,
) -> FitResult:
    """Hierarchical Beta regression on Type-2 AUC for a single group.

    Computes AUC2 per participant (Stage 1) and models the group-level
    distribution using a Beta likelihood (Stage 2).  Unlike meta-d'
    approaches this makes no equal-variance Gaussian SDT assumption —
    AUC2 is a non-parametric measure of metacognitive performance.

    The group mean AUC2 is ``expit(b_Intercept)``.

    Args:
        participants: List of ``(stim, resp, conf)`` tuples, one per participant.
        n_ratings: Number of confidence rating categories.
        chains: MCMC chains (default 4).
        n_iter: Total iterations per chain including warmup (default 2000).
        warmup: Warmup iterations (default 1000).
        seed: Random seed (default 42).
        **kwargs: Forwarded to ``brmspy.brms.brm``.

    Returns:
        ``FitResult`` with posterior for ``b_Intercept`` (logit-scale group mean
        AUC2) and ``phi`` (precision).

    Raises:
        ImportError: If ``brmspy`` is not installed.
        ValueError: If fewer than 3 participants have valid AUC2 estimates.

    Example::

        import numpy as np
        from scipy.special import expit
        import arviz as az
        from metasignal.sdtbayes import fit_beta_auc_group

        rng = np.random.default_rng(0)
        participants = [
            (rng.integers(0, 2, 200), rng.integers(0, 2, 200), rng.integers(1, 5, 200))
            for _ in range(20)
        ]
        fit = fit_beta_auc_group(participants, n_ratings=4)

        intercept = az.extract(fit.idata)["b_Intercept"].values
        print(f"Group mean AUC2: {expit(intercept).mean():.3f}")
    """
    brms = require_brms()

    df = _compute_auc2_per_participant(participants, n_ratings)
    valid = _clip_auc2(df.dropna(subset=["auc2"]))

    if len(valid) < 3:
        msg = f"Only {len(valid)} participants have valid AUC2 — need at least 3."
        raise ValueError(msg)

    priors = [
        brms.prior("normal(0, 2)", class_="Intercept"),
        brms.prior("gamma(0.1, 0.1)", class_="phi"),
    ]

    _result = brms.brm(
        formula=brms.bf("auc2 ~ 1"),
        data=valid,
        family="beta",
        priors=priors,
        chains=chains,
        iter=n_iter,
        warmup=warmup,
        seed=seed,
        **kwargs,
    )
    return FitResult(idata=_result.idata, r=_result.r)

fit_beta_auc_comparison(group_a, group_b, n_ratings, chains=4, n_iter=2000, warmup=1000, seed=42, **kwargs)

Beta regression comparison of Type-2 AUC between two groups.

The key posterior parameter b_group1 is the logit-scale difference in AUC2 between group B and group A. Posterior probability that group B has higher AUC2::

post = az.extract(fit.idata)["b_group1"].values
p_b_higher = (post > 0).mean()

Parameters:

Name Type Description Default
group_a list[tuple[ndarray, ndarray, ndarray]]

Participants in group A — list of (stim, resp, conf) tuples.

 required
group_b list[tuple[ndarray, ndarray, ndarray]]

Participants in group B — list of (stim, resp, conf) tuples.

 required
n_ratings int

Number of confidence rating categories.

 required
chains int

MCMC chains (default 4).

 4
n_iter int

Total iterations per chain including warmup (default 2000).

 2000
warmup int

Warmup iterations (default 1000).

 1000
seed int

Random seed (default 42).

 42
**kwargs Any

Forwarded to brmspy.brms.brm.

 {}

Returns:

Type Description
FitResult

FitResult with posterior for b_Intercept (group A logit-AUC2),
FitResult

b_group1 (group B − A difference), and phi.

Raises:

Type Description
ImportError

If brmspy is not installed.
ValueError

If either group has fewer than 3 participants with valid AUC2 estimates.

Example::

fit = fit_beta_auc_comparison(healthy, patient, n_ratings=4)
import arviz as az
delta = az.extract(fit.idata)["b_group1"].values
print(f"P(patient AUC2 < healthy): {(delta < 0).mean():.3f}")
Source code in src/metasignal/sdtbayes/beta_auc.py 
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
def fit_beta_auc_comparison(
    group_a: list[tuple[np.ndarray, np.ndarray, np.ndarray]],
    group_b: list[tuple[np.ndarray, np.ndarray, np.ndarray]],
    n_ratings: int,
    chains: int = 4,
    n_iter: int = 2000,
    warmup: int = 1000,
    seed: int = 42,
    **kwargs: Any,
) -> FitResult:
    """Beta regression comparison of Type-2 AUC between two groups.

    The key posterior parameter ``b_group1`` is the logit-scale difference in
    AUC2 between group B and group A.  Posterior probability that group B has
    higher AUC2::

        post = az.extract(fit.idata)["b_group1"].values
        p_b_higher = (post > 0).mean()

    Args:
        group_a: Participants in group A — list of ``(stim, resp, conf)`` tuples.
        group_b: Participants in group B — list of ``(stim, resp, conf)`` tuples.
        n_ratings: Number of confidence rating categories.
        chains: MCMC chains (default 4).
        n_iter: Total iterations per chain including warmup (default 2000).
        warmup: Warmup iterations (default 1000).
        seed: Random seed (default 42).
        **kwargs: Forwarded to ``brmspy.brms.brm``.

    Returns:
        ``FitResult`` with posterior for ``b_Intercept`` (group A logit-AUC2),
        ``b_group1`` (group B − A difference), and ``phi``.

    Raises:
        ImportError: If ``brmspy`` is not installed.
        ValueError: If either group has fewer than 3 participants with valid
            AUC2 estimates.

    Example::

        fit = fit_beta_auc_comparison(healthy, patient, n_ratings=4)
        import arviz as az
        delta = az.extract(fit.idata)["b_group1"].values
        print(f"P(patient AUC2 < healthy): {(delta < 0).mean():.3f}")
    """
    brms = require_brms()
    import pandas as pd

    df_a = _compute_auc2_per_participant(group_a, n_ratings)
    df_a["group"] = 0
    df_b = _compute_auc2_per_participant(group_b, n_ratings)
    df_b["group"] = 1
    df = pd.concat([df_a, df_b], ignore_index=True)
    df["group"] = df["group"].astype("category")
    valid = _clip_auc2(df.dropna(subset=["auc2"]))

    for g, label in ((0, "A"), (1, "B")):
        n_valid = int((valid["group"] == g).sum())
        if n_valid < 3:
            total = int((df["group"] == g).sum())
            msg = (
                f"Group {label}: only {n_valid} of {total} participants have valid "
                "AUC2 estimates — need at least 3."
            )
            raise ValueError(msg)

    priors = [
        brms.prior("normal(0, 2)", class_="Intercept"),
        brms.prior("normal(0, 1)", class_="b"),
        brms.prior("gamma(0.1, 0.1)", class_="phi"),
    ]

    _result = brms.brm(
        formula=brms.bf("auc2 ~ group"),
        data=valid,
        family="beta",
        priors=priors,
        chains=chains,
        iter=n_iter,
        warmup=warmup,
        seed=seed,
        **kwargs,
    )
    return FitResult(idata=_result.idata, r=_result.r)

Approach 10 — Bayesian Meta-Regression

fit_two_stage_regression(participants, n_ratings, covariates, formula=None, chains=4, n_iter=2000, warmup=1000, seed=42, **kwargs)

Two-stage Bayesian meta-regression of log M-ratio on covariates.

Stage 1 computes per-participant MLE estimates; Stage 2 fits a brms regression of log M-ratio on the supplied covariates.

Parameters:

Name Type Description Default
participants list[tuple[ndarray, ndarray, ndarray]]

List of (stim, resp, conf) tuples, one per participant. Participant indices must align with the "participant" column in covariates (0-based integers).

 required
n_ratings int

Number of confidence rating categories.

 required
covariates 'pd.DataFrame'

:class:pandas.DataFrame with a "participant" column (integer) and one or more numeric predictor columns. Passed directly to brms after merging with MLE estimates.

 required
formula str | None

brms formula string. Default: "log_m_ratio ~ col1 + col2" using all non-participant columns in covariates.

 None
chains int

MCMC chains (default 4).

 4
n_iter int

Total iterations per chain including warmup (default 2000).

 2000
warmup int

Warmup iterations (default 1000).

 1000
seed int

Random seed (default 42).

 42
**kwargs Any

Forwarded to brmspy.brms.brm.

 {}

Returns:

Type Description
FitResult

Tuple of:
'pd.DataFrame'
  • FitResult with Stage-2 posterior.
tuple[FitResult, 'pd.DataFrame']
  • pd.DataFrame with per-participant MLE estimates and covariates.

Raises:

Type Description
ImportError

If brmspy is not installed.
ValueError

If "participant" column is missing from covariates or fewer than 3 participants have valid log M-ratio estimates.

Example::

import numpy as np
import pandas as pd
from metasignal.sdtbayes import fit_two_stage_regression

rng = np.random.default_rng(0)
N = 30
participants = [
    (rng.integers(0, 2, 200), rng.integers(0, 2, 200), rng.integers(1, 5, 200))
    for _ in range(N)
]
covs = pd.DataFrame({
    "participant": range(N),
    "age": rng.normal(35, 10, N),
    "score": rng.normal(50, 15, N),
})
fit, mle_df = fit_two_stage_regression(participants, n_ratings=4, covariates=covs)
print(fit.posterior_summary(var_names=["b_Intercept", "b_age", "b_score"]))
Source code in src/metasignal/sdtbayes/meta_regression.py 
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
def fit_two_stage_regression(
    participants: list[tuple[np.ndarray, np.ndarray, np.ndarray]],
    n_ratings: int,
    covariates: "pd.DataFrame",
    formula: str | None = None,
    chains: int = 4,
    n_iter: int = 2000,
    warmup: int = 1000,
    seed: int = 42,
    **kwargs: Any,
) -> tuple[FitResult, "pd.DataFrame"]:
    """Two-stage Bayesian meta-regression of log M-ratio on covariates.

    Stage 1 computes per-participant MLE estimates; Stage 2 fits a brms
    regression of log M-ratio on the supplied covariates.

    Args:
        participants: List of ``(stim, resp, conf)`` tuples, one per participant.
            Participant indices must align with the ``"participant"`` column in
            ``covariates`` (0-based integers).
        n_ratings: Number of confidence rating categories.
        covariates: :class:`pandas.DataFrame` with a ``"participant"`` column
            (integer) and one or more numeric predictor columns.  Passed
            directly to brms after merging with MLE estimates.
        formula: brms formula string.  Default: ``"log_m_ratio ~ col1 + col2"``
            using all non-``participant`` columns in ``covariates``.
        chains: MCMC chains (default 4).
        n_iter: Total iterations per chain including warmup (default 2000).
        warmup: Warmup iterations (default 1000).
        seed: Random seed (default 42).
        **kwargs: Forwarded to ``brmspy.brms.brm``.

    Returns:
        Tuple of:

        - ``FitResult`` with Stage-2 posterior.
        - ``pd.DataFrame`` with per-participant MLE estimates and covariates.

    Raises:
        ImportError: If ``brmspy`` is not installed.
        ValueError: If ``"participant"`` column is missing from ``covariates``
            or fewer than 3 participants have valid log M-ratio estimates.

    Example::

        import numpy as np
        import pandas as pd
        from metasignal.sdtbayes import fit_two_stage_regression

        rng = np.random.default_rng(0)
        N = 30
        participants = [
            (rng.integers(0, 2, 200), rng.integers(0, 2, 200), rng.integers(1, 5, 200))
            for _ in range(N)
        ]
        covs = pd.DataFrame({
            "participant": range(N),
            "age": rng.normal(35, 10, N),
            "score": rng.normal(50, 15, N),
        })
        fit, mle_df = fit_two_stage_regression(participants, n_ratings=4, covariates=covs)
        print(fit.posterior_summary(var_names=["b_Intercept", "b_age", "b_score"]))
    """
    brms = require_brms()
    import pandas as pd

    if "participant" not in covariates.columns:
        msg = "'covariates' must contain a 'participant' column."
        raise ValueError(msg)

    mle_df = _compute_participant_estimates(participants, n_ratings)
    merged = mle_df.merge(covariates, on="participant", how="left")
    valid = merged.dropna(subset=["log_m_ratio"])

    if len(valid) < 3:
        msg = f"Only {len(valid)} participants have valid estimates — need at least 3."
        raise ValueError(msg)

    if formula is None:
        pred_cols = [c for c in covariates.columns if c != "participant"]
        formula = "log_m_ratio ~ " + " + ".join(pred_cols)

    priors = [
        brms.prior("normal(0, 1)", class_="Intercept"),
        brms.prior("normal(0, 1)", class_="b"),
        brms.prior("exponential(1)", class_="sigma"),
    ]

    _result = brms.brm(
        formula=brms.bf(formula),
        data=valid,
        family="student",
        priors=priors,
        chains=chains,
        iter=n_iter,
        warmup=warmup,
        seed=seed,
        **kwargs,
    )
    return FitResult(idata=_result.idata, r=_result.r), merged

fit_full_metad_regression(participants, n_ratings, covariates, chains=4, n_iter=2000, warmup=1000, seed=42, tol=1e-07, **kwargs)

Full hierarchical meta-d' with covariate regression on log M-ratio.

A single-stage fully Bayesian model that extends :func:fit_full_metad with participant-level covariates. Covariates are incorporated into the SDT likelihood directly, so estimation uncertainty propagates through the regression.

The model replaces the scalar group hyperprior with a linear predictor::

logMratio[s] ~ Normal(alpha + X[s] @ beta, sigma_logMratio)

where alpha is the expected log M-ratio at the covariate mean and beta is the vector of slopes. Covariates are mean-centred automatically before fitting.

Parameters:

Name Type Description Default
participants list[tuple[ndarray, ndarray, ndarray]]

List of (stim, resp, conf) tuples, one per participant.

 required
n_ratings int

Number of confidence rating categories.

 required
covariates 'np.ndarray'

:class:numpy.ndarray of shape (N,) for a single continuous covariate, or (N, p) for multiple covariates. Must have the same length as participants. Covariates are mean-centred automatically.

 required
chains int

MCMC chains (default 4).

 4
n_iter int

Total iterations per chain including warmup (default 2000).

 2000
warmup int

Warmup iterations (default 1000).

 1000
seed int

Random seed (default 42).

 42
tol float

Minimum probability floor for multinomial cells (default 1e-7).

 1e-07
**kwargs Any

Forwarded to brmspy.brms.brm.

 {}

Returns:

Type Description
FitResult

FitResult. Key parameters:
FitResult
  • alpha_logMratio — intercept (log M-ratio at covariate mean).
FitResult
  • beta_logMratio — slope(s), one per covariate column.
FitResult
  • sigma_logMratio — residual between-subject SD.
FitResult
  • Mratio[s], meta_d[s] — covariate-adjusted per-subject estimates.

Raises:

Type Description
ImportError

If cmdstanpy is not installed.
ValueError

If covariates length does not match participants.

Example::

import numpy as np
from metasignal.sdtbayes import fit_full_metad_regression

rng = np.random.default_rng(0)
N = 25
participants = [
    (rng.integers(0, 2, 200), rng.integers(0, 2, 200), rng.integers(1, 5, 200))
    for _ in range(N)
]
age = rng.normal(35, 10, N)   # single continuous covariate

fit = fit_full_metad_regression(participants, n_ratings=4, covariates=age)
print(fit.posterior_summary(var_names=["alpha_logMratio", "beta_logMratio",
                                        "sigma_logMratio"]))
Notes

Delegates to :func:metasignal.sdtbayes.fit_meta_formula (parameterization="mratio", backend="stan") rather than injecting custom Stan code into brms via brmspy — the brmspy/rpy2 bridge cannot reliably convert a list of multiple stanvar() objects into R (see :func:metasignal.sdtbayes.fit_full_metad docstring for details).

Source code in src/metasignal/sdtbayes/meta_regression.py 
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
def fit_full_metad_regression(
    participants: list[tuple[np.ndarray, np.ndarray, np.ndarray]],
    n_ratings: int,
    covariates: "np.ndarray",
    chains: int = 4,
    n_iter: int = 2000,
    warmup: int = 1000,
    seed: int = 42,
    tol: float = 1e-7,
    **kwargs: Any,
) -> FitResult:
    """Full hierarchical meta-d' with covariate regression on log M-ratio.

    A single-stage fully Bayesian model that extends :func:`fit_full_metad`
    with participant-level covariates.  Covariates are incorporated into the
    SDT likelihood directly, so estimation uncertainty propagates through the
    regression.

    The model replaces the scalar group hyperprior with a linear predictor::

        logMratio[s] ~ Normal(alpha + X[s] @ beta, sigma_logMratio)

    where ``alpha`` is the expected log M-ratio at the covariate mean and
    ``beta`` is the vector of slopes.  Covariates are mean-centred
    automatically before fitting.

    Args:
        participants: List of ``(stim, resp, conf)`` tuples, one per participant.
        n_ratings: Number of confidence rating categories.
        covariates: :class:`numpy.ndarray` of shape ``(N,)`` for a single
            continuous covariate, or ``(N, p)`` for multiple covariates.
            Must have the same length as ``participants``.  Covariates are
            mean-centred automatically.
        chains: MCMC chains (default 4).
        n_iter: Total iterations per chain including warmup (default 2000).
        warmup: Warmup iterations (default 1000).
        seed: Random seed (default 42).
        tol: Minimum probability floor for multinomial cells (default 1e-7).
        **kwargs: Forwarded to ``brmspy.brms.brm``.

    Returns:
        ``FitResult``.  Key parameters:

        - ``alpha_logMratio`` — intercept (log M-ratio at covariate mean).
        - ``beta_logMratio`` — slope(s), one per covariate column.
        - ``sigma_logMratio`` — residual between-subject SD.
        - ``Mratio[s]``, ``meta_d[s]`` — covariate-adjusted per-subject estimates.

    Raises:
        ImportError: If ``cmdstanpy`` is not installed.
        ValueError: If ``covariates`` length does not match ``participants``.

    Example::

        import numpy as np
        from metasignal.sdtbayes import fit_full_metad_regression

        rng = np.random.default_rng(0)
        N = 25
        participants = [
            (rng.integers(0, 2, 200), rng.integers(0, 2, 200), rng.integers(1, 5, 200))
            for _ in range(N)
        ]
        age = rng.normal(35, 10, N)   # single continuous covariate

        fit = fit_full_metad_regression(participants, n_ratings=4, covariates=age)
        print(fit.posterior_summary(var_names=["alpha_logMratio", "beta_logMratio",
                                                "sigma_logMratio"]))

    Notes:
        Delegates to :func:`metasignal.sdtbayes.fit_meta_formula`
        (``parameterization="mratio"``, ``backend="stan"``) rather than
        injecting custom Stan code into brms via ``brmspy`` — the brmspy/rpy2
        bridge cannot reliably convert a list of multiple ``stanvar()``
        objects into R (see :func:`metasignal.sdtbayes.fit_full_metad`
        docstring for details).
    """
    import pandas as pd

    from metasignal.sdtbayes.formula import fit_meta_formula

    cov_mat = np.atleast_2d(np.asarray(covariates, dtype=float))
    if cov_mat.shape[0] == 1 and cov_mat.shape[1] == len(participants):
        cov_mat = cov_mat.T
    if cov_mat.shape[0] != len(participants):
        msg = (
            f"covariates has {cov_mat.shape[0]} rows but participants has "
            f"{len(participants)} entries."
        )
        raise ValueError(msg)

    p_cov = cov_mat.shape[1]
    col_names = [f"cov{i}" for i in range(p_cov)]
    pred_df = pd.DataFrame(cov_mat, columns=col_names)
    pred_df.insert(0, "participant", range(len(participants)))
    formula = "~ " + " + ".join(col_names)

    return fit_meta_formula(
        participants=participants,
        n_ratings=n_ratings,
        formula=formula,
        data=pred_df,
        parameterization="mratio",
        backend="stan",
        chains=chains,
        n_iter=n_iter,
        warmup=warmup,
        seed=seed,
        tol=tol,
        **kwargs,
    )

Approach 12 — Within-Subject Condition Comparison

fit_within_subject_comparison(condition_a, condition_b, n_ratings, chains=4, n_iter=2000, warmup=1000, seed=42, **kwargs)

Two-stage within-subject comparison of metacognitive efficiency across conditions.

Estimates the change in log M-ratio between two conditions tested within the same participants. Because individual participants serve as their own controls, this design is statistically more efficient than between-subject group comparisons and controls for stable differences in metacognitive ability.

The key estimand is b_condition1, the posterior distribution of the within-person condition-B − condition-A difference in log M-ratio.

Parameters:

Name Type Description Default
condition_a list[tuple[ndarray, ndarray, ndarray]]

Per-participant (stim, resp, conf) tuples for condition A.

 required
condition_b list[tuple[ndarray, ndarray, ndarray]]

Per-participant (stim, resp, conf) tuples for condition B. Must have the same length as condition_a and list participants in the same order.

 required
n_ratings int

Number of confidence rating categories.

 required
chains int

MCMC chains (default 4).

 4
n_iter int

Total iterations per chain including warmup (default 2000).

 2000
warmup int

Warmup iterations (default 1000).

 1000
seed int

Random seed (default 42).

 42
**kwargs Any

Forwarded to brmspy.brms.brm.

 {}

Returns:

Type Description
FitResult

Tuple of (FitResult, mle_df).
'pd.DataFrame'
  • FitResult.idata — ArviZ InferenceData. Key parameter: b_condition1 = condition B − condition A difference.
tuple[FitResult, 'pd.DataFrame']
  • mle_df — long-format pd.DataFrame with Stage-1 MLE estimates for all participants × conditions.

Raises:

Type Description
ImportError

If brmspy is not installed.
ValueError

If condition_a and condition_b have different lengths, or if fewer than 3 participants have valid estimates in both conditions.

Example::

import numpy as np
import arviz as az
from metasignal.sdtbayes import fit_within_subject_comparison

rng = np.random.default_rng(0)
n_part = 25

# Same 25 participants, two experimental conditions
condition_a = [
    (rng.integers(0, 2, 200), rng.integers(0, 2, 200), rng.integers(1, 5, 200))
    for _ in range(n_part)
]
condition_b = [
    (rng.integers(0, 2, 200), rng.integers(0, 2, 200), rng.integers(1, 5, 200))
    for _ in range(n_part)
]

fit, mle_df = fit_within_subject_comparison(condition_a, condition_b, n_ratings=4)

# Posterior probability that condition B has higher metacognition
post = az.extract(fit.idata)["b_condition1"].values
print(f"P(M-ratio higher in B): {(post > 0).mean():.3f}")

# Condition-B minus condition-A M-ratio ratio
print(f"Median M-ratio ratio B/A: {np.exp(np.median(post)):.3f}")
Source code in src/metasignal/sdtbayes/within_subject.py 
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
def fit_within_subject_comparison(
    condition_a: list[tuple[np.ndarray, np.ndarray, np.ndarray]],
    condition_b: list[tuple[np.ndarray, np.ndarray, np.ndarray]],
    n_ratings: int,
    chains: int = 4,
    n_iter: int = 2000,
    warmup: int = 1000,
    seed: int = 42,
    **kwargs: Any,
) -> tuple[FitResult, "pd.DataFrame"]:
    """Two-stage within-subject comparison of metacognitive efficiency across conditions.

    Estimates the change in log M-ratio between two conditions tested within
    the same participants.  Because individual participants serve as their own
    controls, this design is statistically more efficient than between-subject
    group comparisons and controls for stable differences in metacognitive
    ability.

    The key estimand is ``b_condition1``, the posterior distribution of the
    within-person condition-B − condition-A difference in log M-ratio.

    Args:
        condition_a: Per-participant ``(stim, resp, conf)`` tuples for condition A.
        condition_b: Per-participant ``(stim, resp, conf)`` tuples for condition B.
            Must have the same length as ``condition_a`` and list participants in
            the same order.
        n_ratings: Number of confidence rating categories.
        chains: MCMC chains (default 4).
        n_iter: Total iterations per chain including warmup (default 2000).
        warmup: Warmup iterations (default 1000).
        seed: Random seed (default 42).
        **kwargs: Forwarded to ``brmspy.brms.brm``.

    Returns:
        Tuple of ``(FitResult, mle_df)``.

        - ``FitResult.idata`` — ArviZ InferenceData.  Key parameter:
            ``b_condition1`` = condition B − condition A difference.
        - ``mle_df`` — long-format ``pd.DataFrame`` with Stage-1 MLE estimates
            for all participants × conditions.

    Raises:
        ImportError: If ``brmspy`` is not installed.
        ValueError: If ``condition_a`` and ``condition_b`` have different lengths,
            or if fewer than 3 participants have valid estimates in both conditions.

    Example::

        import numpy as np
        import arviz as az
        from metasignal.sdtbayes import fit_within_subject_comparison

        rng = np.random.default_rng(0)
        n_part = 25

        # Same 25 participants, two experimental conditions
        condition_a = [
            (rng.integers(0, 2, 200), rng.integers(0, 2, 200), rng.integers(1, 5, 200))
            for _ in range(n_part)
        ]
        condition_b = [
            (rng.integers(0, 2, 200), rng.integers(0, 2, 200), rng.integers(1, 5, 200))
            for _ in range(n_part)
        ]

        fit, mle_df = fit_within_subject_comparison(condition_a, condition_b, n_ratings=4)

        # Posterior probability that condition B has higher metacognition
        post = az.extract(fit.idata)["b_condition1"].values
        print(f"P(M-ratio higher in B): {(post > 0).mean():.3f}")

        # Condition-B minus condition-A M-ratio ratio
        print(f"Median M-ratio ratio B/A: {np.exp(np.median(post)):.3f}")
    """
    brms = require_brms()

    if len(condition_a) != len(condition_b):
        msg = (
            f"condition_a has {len(condition_a)} participants but condition_b has "
            f"{len(condition_b)}.  Both must list the same participants in the same order."
        )
        raise ValueError(msg)

    mle_df = _compute_paired_estimates(condition_a, condition_b, n_ratings)
    mle_df["condition"] = mle_df["condition"].astype("category")

    valid = mle_df.dropna(subset=["log_m_ratio"])

    # Require each participant to have at least one valid condition
    valid_participants = valid["participant"].nunique()
    if valid_participants < 3:
        msg = (
            f"Only {valid_participants} participants have any valid MLE estimates — "
            "need at least 3."
        )
        raise ValueError(msg)

    formula = brms.bf("log_m_ratio ~ condition + (1 | participant)")
    priors = [
        brms.prior("normal(0, 1)", class_="Intercept"),
        brms.prior("normal(0, 1)", class_="b"),
        brms.prior("exponential(1)", class_="sd"),
        brms.prior("exponential(1)", class_="sigma"),
    ]

    _result = brms.brm(
        formula=formula,
        data=valid,
        family="student",
        priors=priors,
        chains=chains,
        iter=n_iter,
        warmup=warmup,
        seed=seed,
        **kwargs,
    )
    return FitResult(idata=_result.idata, r=_result.r), mle_df

Formula Interface

fit_meta_formula(participants, n_ratings, formula, data, *, parameterization='mratio', backend='stan', chains=4, n_iter=2000, warmup=1000, seed=42, tol=1e-07, n_gh=20, cmdstan_path=None, **kwargs)

Fit a hierarchical Bayesian metacognition model using a brms-style formula.

Parameters:

Name Type Description Default
participants list[tuple[ndarray, ndarray, ndarray]]

List of (stim, resp, conf) tuples, one per participant. stim and resp are 0/1 integer arrays; conf is an integer array with values 1 … n_ratings.

 required
n_ratings int

Number of confidence rating levels.

 required
formula str

Wilkinson formula string. The LHS is ignored. Examples::

"~ sensory_cond * perform_cond + vviq_c"
"~ sensory_cond + perform_cond + (1 | group)"   # brms only
required
data 'pd.DataFrame'

:class:pandas.DataFrame with one row per participant containing all predictor columns named in the formula. Row order must match participants, or include a "participant" integer column.

 required
parameterization str

Generative model to use:

  • "mratio" (default) — Fleming (2017) HMeta-d. Primary group parameter: alpha_logMratio (log M-ratio).
  • "meta_noise" — Guggenmos (2022) / Maniscalco & Lau (2014). Primary group parameter: mu_logSigmaMeta (log σ_meta).
  • "casandre" — Boundy-Singer et al. (2023) CASANDRE. Primary group parameter: mu_logPhi (log meta-uncertainty φ).
 'mratio'
backend str

"stan" (cmdstanpy, default) or "brms" (two-stage brmspy, parameterization="mratio" only).

 'stan'
chains int

MCMC chains (default 4).

 4
n_iter int

Total iterations per chain including warmup (default 2000).

 2000
warmup int

Warmup iterations (default 1000).

 1000
seed int

Random seed (default 42).

 42
tol float

Probability floor for multinomial cells (default 1e-7). Stan only.

 1e-07
n_gh int

Number of Gauss-Hermite quadrature points for the "casandre" model (default 20; increase for more accuracy at higher cost).

 20
cmdstan_path str | None

Path to CmdStan installation. If None, auto-detected from cmdstanpy config or ~/.cmdstan/cmdstan-*. Stan only.

 None
**kwargs Any

Forwarded to cmdstanpy.CmdStanModel.sample (Stan) or brmspy.brms.brm (brms).

 {}

Returns:

Type Description
FitResult | tuple[FitResult, 'pd.DataFrame']

FitResult — Stan backend.
FitResult | tuple[FitResult, 'pd.DataFrame']

(FitResult, merged_df) — brms backend.

Examples

M-ratio model (default)::

fit = fit_meta_formula(
    participants, n_ratings=4,
    formula="~ sensory_cond * perform_cond + vviq_c",
    data=df,
    parameterization="mratio",
)
# Key posteriors: alpha_logMratio, beta_logMratio, group_Mratio, Mratio[s]

Meta-noise model::

fit = fit_meta_formula(
    participants, n_ratings=4,
    formula="~ perform_cond + vviq_c",
    data=df,
    parameterization="meta_noise",
)
# Key posteriors: mu_logSigmaMeta, sigma_meta[s], group_sigma_meta

CASANDRE meta-uncertainty model::

fit = fit_meta_formula(
    participants, n_ratings=4,
    formula="~ sensory_cond + perform_cond",
    data=df,
    parameterization="casandre",
    n_gh=20,
)
# Key posteriors: mu_logPhi, phi[s], group_phi, log_theta[1..K-1]
Source code in src/metasignal/sdtbayes/formula.py 
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
def fit_meta_formula(
    participants: list[tuple[np.ndarray, np.ndarray, np.ndarray]],
    n_ratings: int,
    formula: str,
    data: "pd.DataFrame",
    *,
    parameterization: str = "mratio",
    backend: str = "stan",
    chains: int = 4,
    n_iter: int = 2000,
    warmup: int = 1000,
    seed: int = 42,
    tol: float = 1e-7,
    n_gh: int = 20,
    cmdstan_path: str | None = None,
    **kwargs: Any,
) -> FitResult | tuple[FitResult, "pd.DataFrame"]:
    """Fit a hierarchical Bayesian metacognition model using a brms-style formula.

    Args:
        participants: List of ``(stim, resp, conf)`` tuples, one per participant.
            ``stim`` and ``resp`` are 0/1 integer arrays; ``conf`` is an integer
            array with values ``1 … n_ratings``.
        n_ratings: Number of confidence rating levels.
        formula: Wilkinson formula string.  The LHS is ignored.  Examples::

                "~ sensory_cond * perform_cond + vviq_c"
                "~ sensory_cond + perform_cond + (1 | group)"   # brms only

        data: :class:`pandas.DataFrame` with one row per participant containing
            all predictor columns named in the formula.  Row order must match
            ``participants``, or include a ``"participant"`` integer column.
        parameterization: Generative model to use:

            - ``"mratio"`` *(default)* — Fleming (2017) HMeta-d.
              Primary group parameter: ``alpha_logMratio`` (log M-ratio).
            - ``"meta_noise"`` — Guggenmos (2022) / Maniscalco & Lau (2014).
              Primary group parameter: ``mu_logSigmaMeta`` (log σ_meta).
            - ``"casandre"`` — Boundy-Singer et al. (2023) CASANDRE.
              Primary group parameter: ``mu_logPhi`` (log meta-uncertainty φ).

        backend: ``"stan"`` (cmdstanpy, default) or ``"brms"`` (two-stage brmspy,
            ``parameterization="mratio"`` only).
        chains: MCMC chains (default 4).
        n_iter: Total iterations per chain including warmup (default 2000).
        warmup: Warmup iterations (default 1000).
        seed: Random seed (default 42).
        tol: Probability floor for multinomial cells (default 1e-7). Stan only.
        n_gh: Number of Gauss-Hermite quadrature points for the ``"casandre"``
            model (default 20; increase for more accuracy at higher cost).
        cmdstan_path: Path to CmdStan installation.  If ``None``, auto-detected
            from cmdstanpy config or ``~/.cmdstan/cmdstan-*``. Stan only.
        **kwargs: Forwarded to ``cmdstanpy.CmdStanModel.sample`` (Stan) or
            ``brmspy.brms.brm`` (brms).

    Returns:
        ``FitResult`` — Stan backend.
        ``(FitResult, merged_df)`` — brms backend.

    Examples
    --------
    M-ratio model (default)::

        fit = fit_meta_formula(
            participants, n_ratings=4,
            formula="~ sensory_cond * perform_cond + vviq_c",
            data=df,
            parameterization="mratio",
        )
        # Key posteriors: alpha_logMratio, beta_logMratio, group_Mratio, Mratio[s]

    Meta-noise model::

        fit = fit_meta_formula(
            participants, n_ratings=4,
            formula="~ perform_cond + vviq_c",
            data=df,
            parameterization="meta_noise",
        )
        # Key posteriors: mu_logSigmaMeta, sigma_meta[s], group_sigma_meta

    CASANDRE meta-uncertainty model::

        fit = fit_meta_formula(
            participants, n_ratings=4,
            formula="~ sensory_cond + perform_cond",
            data=df,
            parameterization="casandre",
            n_gh=20,
        )
        # Key posteriors: mu_logPhi, phi[s], group_phi, log_theta[1..K-1]
    """
    if parameterization not in _VALID_PARAMETERIZATIONS:
        raise ValueError(
            f"Unknown parameterization {parameterization!r}. "
            f"Choose from: {_VALID_PARAMETERIZATIONS}"
        )

    if backend == "brms" and parameterization != "mratio":
        raise ValueError(
            "The brms two-stage backend only supports parameterization='mratio'. "
            "For meta_noise or casandre use backend='stan'."
        )

    # Strip and clean formula RHS
    rhs = formula.split("~", 1)[1].strip() if "~" in formula else formula.strip()

    if backend == "stan":
        if "|" in rhs:
            warnings.warn(
                "Random-effects terms (|) are not supported in the Stan backend "
                "and will be ignored. The hierarchical Stan model already pools "
                "over participants. Use backend='brms' (mratio only) for additional "
                "random-effects structure.",
                UserWarning,
                stacklevel=2,
            )
            rhs = " + ".join(t for t in rhs.split("+") if "|" not in t)
        return _fit_stan(
            participants, n_ratings, rhs, data,
            parameterization=parameterization,
            chains=chains, n_iter=n_iter, warmup=warmup, seed=seed,
            tol=tol, cmdstan_path=cmdstan_path, n_gh=n_gh, **kwargs,
        )

    if backend == "brms":
        return _fit_brms(
            participants, n_ratings, formula, data,
            chains=chains, n_iter=n_iter, warmup=warmup, seed=seed, **kwargs,
        )

    raise ValueError(f"Unknown backend {backend!r}. Choose 'stan' or 'brms'.")

Runtime and Result Types

setup_runtime(use_prebuilt=True)

Install the brms/Stan toolchain (run once before fitting any model).

Downloads a prebuilt Stan runtime and installs the brms R package into an isolated environment managed by brmspy. Internet access is required on first call; subsequent calls are no-ops if the runtime is already present.

Parameters:

Name Type Description Default
use_prebuilt bool

Use a prebuilt Stan binary (faster). Set False to compile from source (requires a C++ toolchain).

 True

Example::

from metasignal.sdtbayes import setup_runtime
setup_runtime()   # ~5 min on first run
Source code in src/metasignal/sdtbayes/_runtime.py 
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
def setup_runtime(use_prebuilt: bool = True) -> None:
    """Install the brms/Stan toolchain (run once before fitting any model).

    Downloads a prebuilt Stan runtime and installs the brms R package into
    an isolated environment managed by brmspy.  Internet access is required
    on first call; subsequent calls are no-ops if the runtime is already
    present.

    Args:
        use_prebuilt: Use a prebuilt Stan binary (faster). Set ``False`` to
            compile from source (requires a C++ toolchain).

    Example::

        from metasignal.sdtbayes import setup_runtime
        setup_runtime()   # ~5 min on first run
    """
    try:
        from brmspy import brms
    except ImportError as e:
        msg = (
            "brmspy is not installed. Run:\n"
            '    pip install "metasignal[sdtbayes] @ git+https://github.com/saurabhr/metasignal.git"'
        )
        raise ImportError(msg) from e

    with brms.manage() as ctx:
        ctx.install_brms(use_prebuilt=use_prebuilt)

FitResult dataclass

Container returned by every sdtbayes fitting function.

Wraps the brmspy/Stan fit with ArviZ diagnostics as methods, so that every approach exposes the same interface::

fit = fit_hierarchical_metad(participants, n_ratings=4)
fit.posterior_summary()
fit.plot_trace(var_names=["b_correct"])
fit.convergence_diagnostics()

Attributes:

Name Type Description
idata Any

ArviZ InferenceData with posterior samples, log-likelihood, and (if requested) prior predictive.
r Any

Lightweight R object handle for downstream brms calls (e.g. brms.hypothesis(fit.r, "b_correct > 0")).
Source code in src/metasignal/sdtbayes/diagnostics.py 
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
@dataclass
class FitResult:
    """Container returned by every sdtbayes fitting function.

    Wraps the brmspy/Stan fit with ArviZ diagnostics as methods, so that
    every approach exposes the same interface::

        fit = fit_hierarchical_metad(participants, n_ratings=4)
        fit.posterior_summary()
        fit.plot_trace(var_names=["b_correct"])
        fit.convergence_diagnostics()

    Attributes:
        idata: ArviZ ``InferenceData`` with posterior samples, log-likelihood,
            and (if requested) prior predictive.
        r: Lightweight R object handle for downstream brms calls (e.g.
            ``brms.hypothesis(fit.r, "b_correct > 0")``).
    """

    idata: Any
    r: Any

    # ------------------------------------------------------------------
    # Diagnostics — one implementation, shared by all four approaches
    # ------------------------------------------------------------------

    def posterior_summary(
        self,
        var_names: Sequence[str] | None = None,
        *,
        hdi_prob: float = 0.94,
    ) -> Any:
        """Posterior summary table (mean, SD, HDI, R-hat, ESS).

        Args:
            var_names: Parameter names to include. ``None`` returns all.
            hdi_prob: Highest density interval probability mass (default 0.94).

        Returns:
            ``pandas.DataFrame`` with one row per parameter.

        Example::

            summary = fit.posterior_summary()
            print(summary[["mean", "sd", "hdi_3%", "hdi_97%", "r_hat"]])
        """
        az = _require_arviz()
        return az.summary(self.idata, var_names=var_names, hdi_prob=hdi_prob)

    def convergence_diagnostics(self) -> Any:
        """R-hat and ESS for all parameters, with a ``converged`` flag.

        R-hat > 1.01 or ESS < 400 indicates sampling problems.

        Returns:
            ``pandas.DataFrame`` with columns ``r_hat``, ``ess_bulk``,
            ``ess_tail``, ``converged``.

        Example::

            diag = fit.convergence_diagnostics()
            print(diag[~diag["converged"]])   # non-converged parameters
        """
        az = _require_arviz()
        # round_to="none": az.summary()'s default "auto" rounding formats
        # numeric columns (including r_hat) as display strings, which breaks
        # the numeric comparisons below.
        summary = az.summary(self.idata, round_to="none")
        diag = summary[["r_hat", "ess_bulk", "ess_tail"]].astype(float).copy()
        diag["converged"] = (
            (diag["r_hat"] <= 1.01)
            & (diag["ess_bulk"] >= 400)
            & (diag["ess_tail"] >= 400)
        )
        return diag

    def plot_trace(
        self,
        var_names: Sequence[str] | None = None,
        **kwargs: Any,
    ) -> Any:
        """MCMC trace and rank plots for convergence inspection.

        Args:
            var_names: Parameters to plot. ``None`` plots all (can be slow).
            **kwargs: Forwarded to ``arviz.plot_trace``.

        Returns:
            Array of ``matplotlib.axes.Axes``.

        Example::

            fit.plot_trace(var_names=["b_correct", "sd_participant__correct"])
        """
        az = _require_arviz()
        return az.plot_trace(self.idata, var_names=var_names, **kwargs)

    def plot_posterior(
        self,
        var_names: Sequence[str] | None = None,
        *,
        hdi_prob: float = 0.94,
        ref_val: float | None = None,
        **kwargs: Any,
    ) -> Any:
        """Marginal posterior density plots with HDI.

        Args:
            var_names: Parameters to plot.
            hdi_prob: HDI probability mass to shade (default 0.94).
            ref_val: Optional reference line (e.g. ``0`` for a null).
            **kwargs: Forwarded to ``arviz.plot_posterior``.

        Returns:
            Array of ``matplotlib.axes.Axes``.

        Example::

            fit.plot_posterior(var_names=["b_correct:group1"], ref_val=0)
        """
        az = _require_arviz()
        kw: dict[str, Any] = {"hdi_prob": hdi_prob}
        if ref_val is not None:
            kw["ref_val"] = ref_val
        kw.update(kwargs)
        return az.plot_posterior(self.idata, var_names=var_names, **kw)

    def plot_forest(
        self,
        var_names: Sequence[str] | None = None,
        *,
        hdi_prob: float = 0.94,
        r_hat: bool = True,
        ess: bool = True,
        **kwargs: Any,
    ) -> Any:
        """Forest plot of parameter estimates with HDI intervals.

        Args:
            var_names: Parameters to include. Defaults to all.
            hdi_prob: HDI probability mass (default 0.94).
            r_hat: Show R-hat values (default True).
            ess: Show effective sample size (default True).
            **kwargs: Forwarded to ``arviz.plot_forest``.

        Returns:
            Array of ``matplotlib.axes.Axes``.

        Example::

            fit.plot_forest(var_names=["b_correct", "r_participant"])
        """
        az = _require_arviz()
        return az.plot_forest(
            self.idata,
            var_names=var_names,
            hdi_prob=hdi_prob,
            r_hat=r_hat,
            ess=ess,
            **kwargs,
        )

convergence_diagnostics()

R-hat and ESS for all parameters, with a converged flag.

R-hat > 1.01 or ESS < 400 indicates sampling problems.

Returns:

Type Description
Any

pandas.DataFrame with columns r_hat, ess_bulk,
Any

ess_tail, converged.

Example::

diag = fit.convergence_diagnostics()
print(diag[~diag["converged"]])   # non-converged parameters
Source code in src/metasignal/sdtbayes/diagnostics.py 
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
def convergence_diagnostics(self) -> Any:
    """R-hat and ESS for all parameters, with a ``converged`` flag.

    R-hat > 1.01 or ESS < 400 indicates sampling problems.

    Returns:
        ``pandas.DataFrame`` with columns ``r_hat``, ``ess_bulk``,
        ``ess_tail``, ``converged``.

    Example::

        diag = fit.convergence_diagnostics()
        print(diag[~diag["converged"]])   # non-converged parameters
    """
    az = _require_arviz()
    # round_to="none": az.summary()'s default "auto" rounding formats
    # numeric columns (including r_hat) as display strings, which breaks
    # the numeric comparisons below.
    summary = az.summary(self.idata, round_to="none")
    diag = summary[["r_hat", "ess_bulk", "ess_tail"]].astype(float).copy()
    diag["converged"] = (
        (diag["r_hat"] <= 1.01)
        & (diag["ess_bulk"] >= 400)
        & (diag["ess_tail"] >= 400)
    )
    return diag

plot_forest(var_names=None, *, hdi_prob=0.94, r_hat=True, ess=True, **kwargs)

Forest plot of parameter estimates with HDI intervals.

Parameters:

Name Type Description Default
var_names Sequence[str] | None

Parameters to include. Defaults to all.

 None
hdi_prob float

HDI probability mass (default 0.94).

 0.94
r_hat bool

Show R-hat values (default True).

 True
ess bool

Show effective sample size (default True).

 True
**kwargs Any

Forwarded to arviz.plot_forest.

 {}

Returns:

Type Description
Any

Array of matplotlib.axes.Axes.

Example::

fit.plot_forest(var_names=["b_correct", "r_participant"])
Source code in src/metasignal/sdtbayes/diagnostics.py 
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
def plot_forest(
    self,
    var_names: Sequence[str] | None = None,
    *,
    hdi_prob: float = 0.94,
    r_hat: bool = True,
    ess: bool = True,
    **kwargs: Any,
) -> Any:
    """Forest plot of parameter estimates with HDI intervals.

    Args:
        var_names: Parameters to include. Defaults to all.
        hdi_prob: HDI probability mass (default 0.94).
        r_hat: Show R-hat values (default True).
        ess: Show effective sample size (default True).
        **kwargs: Forwarded to ``arviz.plot_forest``.

    Returns:
        Array of ``matplotlib.axes.Axes``.

    Example::

        fit.plot_forest(var_names=["b_correct", "r_participant"])
    """
    az = _require_arviz()
    return az.plot_forest(
        self.idata,
        var_names=var_names,
        hdi_prob=hdi_prob,
        r_hat=r_hat,
        ess=ess,
        **kwargs,
    )

plot_posterior(var_names=None, *, hdi_prob=0.94, ref_val=None, **kwargs)

Marginal posterior density plots with HDI.

Parameters:

Name Type Description Default
var_names Sequence[str] | None

Parameters to plot.

 None
hdi_prob float

HDI probability mass to shade (default 0.94).

 0.94
ref_val float | None

Optional reference line (e.g. 0 for a null).

 None
**kwargs Any

Forwarded to arviz.plot_posterior.

 {}

Returns:

Type Description
Any

Array of matplotlib.axes.Axes.

Example::

fit.plot_posterior(var_names=["b_correct:group1"], ref_val=0)
Source code in src/metasignal/sdtbayes/diagnostics.py 
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
def plot_posterior(
    self,
    var_names: Sequence[str] | None = None,
    *,
    hdi_prob: float = 0.94,
    ref_val: float | None = None,
    **kwargs: Any,
) -> Any:
    """Marginal posterior density plots with HDI.

    Args:
        var_names: Parameters to plot.
        hdi_prob: HDI probability mass to shade (default 0.94).
        ref_val: Optional reference line (e.g. ``0`` for a null).
        **kwargs: Forwarded to ``arviz.plot_posterior``.

    Returns:
        Array of ``matplotlib.axes.Axes``.

    Example::

        fit.plot_posterior(var_names=["b_correct:group1"], ref_val=0)
    """
    az = _require_arviz()
    kw: dict[str, Any] = {"hdi_prob": hdi_prob}
    if ref_val is not None:
        kw["ref_val"] = ref_val
    kw.update(kwargs)
    return az.plot_posterior(self.idata, var_names=var_names, **kw)

plot_trace(var_names=None, **kwargs)

MCMC trace and rank plots for convergence inspection.

Parameters:

Name Type Description Default
var_names Sequence[str] | None

Parameters to plot. None plots all (can be slow).

 None
**kwargs Any

Forwarded to arviz.plot_trace.

 {}

Returns:

Type Description
Any

Array of matplotlib.axes.Axes.

Example::

fit.plot_trace(var_names=["b_correct", "sd_participant__correct"])
Source code in src/metasignal/sdtbayes/diagnostics.py 
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
def plot_trace(
    self,
    var_names: Sequence[str] | None = None,
    **kwargs: Any,
) -> Any:
    """MCMC trace and rank plots for convergence inspection.

    Args:
        var_names: Parameters to plot. ``None`` plots all (can be slow).
        **kwargs: Forwarded to ``arviz.plot_trace``.

    Returns:
        Array of ``matplotlib.axes.Axes``.

    Example::

        fit.plot_trace(var_names=["b_correct", "sd_participant__correct"])
    """
    az = _require_arviz()
    return az.plot_trace(self.idata, var_names=var_names, **kwargs)

posterior_summary(var_names=None, *, hdi_prob=0.94)

Posterior summary table (mean, SD, HDI, R-hat, ESS).

Parameters:

Name Type Description Default
var_names Sequence[str] | None

Parameter names to include. None returns all.

 None
hdi_prob float

Highest density interval probability mass (default 0.94).

 0.94

Returns:

Type Description
Any

pandas.DataFrame with one row per parameter.

Example::

summary = fit.posterior_summary()
print(summary[["mean", "sd", "hdi_3%", "hdi_97%", "r_hat"]])
Source code in src/metasignal/sdtbayes/diagnostics.py 
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
def posterior_summary(
    self,
    var_names: Sequence[str] | None = None,
    *,
    hdi_prob: float = 0.94,
) -> Any:
    """Posterior summary table (mean, SD, HDI, R-hat, ESS).

    Args:
        var_names: Parameter names to include. ``None`` returns all.
        hdi_prob: Highest density interval probability mass (default 0.94).

    Returns:
        ``pandas.DataFrame`` with one row per parameter.

    Example::

        summary = fit.posterior_summary()
        print(summary[["mean", "sd", "hdi_3%", "hdi_97%", "r_hat"]])
    """
    az = _require_arviz()
    return az.summary(self.idata, var_names=var_names, hdi_prob=hdi_prob)

Diagnostics and Posteriors

posterior_summary(fit, var_names=None, *, hdi_prob=0.94)

Module-level alias for :meth:FitResult.posterior_summary.

Source code in src/metasignal/sdtbayes/diagnostics.py 
181
182
183
184
185
186
187
188
def posterior_summary(
    fit: FitResult,
    var_names: Sequence[str] | None = None,
    *,
    hdi_prob: float = 0.94,
) -> Any:
    """Module-level alias for :meth:`FitResult.posterior_summary`."""
    return fit.posterior_summary(var_names=var_names, hdi_prob=hdi_prob)

plot_trace(fit, var_names=None, **kwargs)

Module-level alias for :meth:FitResult.plot_trace.

Source code in src/metasignal/sdtbayes/diagnostics.py 
196
197
198
199
200
201
202
def plot_trace(
    fit: FitResult,
    var_names: Sequence[str] | None = None,
    **kwargs: Any,
) -> Any:
    """Module-level alias for :meth:`FitResult.plot_trace`."""
    return fit.plot_trace(var_names=var_names, **kwargs)

plot_posterior(fit, var_names=None, *, hdi_prob=0.94, ref_val=None, **kwargs)

Module-level alias for :meth:FitResult.plot_posterior.

Source code in src/metasignal/sdtbayes/diagnostics.py 
205
206
207
208
209
210
211
212
213
214
def plot_posterior(
    fit: FitResult,
    var_names: Sequence[str] | None = None,
    *,
    hdi_prob: float = 0.94,
    ref_val: float | None = None,
    **kwargs: Any,
) -> Any:
    """Module-level alias for :meth:`FitResult.plot_posterior`."""
    return fit.plot_posterior(var_names=var_names, hdi_prob=hdi_prob, ref_val=ref_val, **kwargs)

plot_forest(fit, var_names=None, *, hdi_prob=0.94, r_hat=True, ess=True, **kwargs)

Module-level alias for :meth:FitResult.plot_forest.

Source code in src/metasignal/sdtbayes/diagnostics.py 
217
218
219
220
221
222
223
224
225
226
227
def plot_forest(
    fit: FitResult,
    var_names: Sequence[str] | None = None,
    *,
    hdi_prob: float = 0.94,
    r_hat: bool = True,
    ess: bool = True,
    **kwargs: Any,
) -> Any:
    """Module-level alias for :meth:`FitResult.plot_forest`."""
    return fit.plot_forest(var_names=var_names, hdi_prob=hdi_prob, r_hat=r_hat, ess=ess, **kwargs)

convergence_diagnostics(fit)

Module-level alias for :meth:FitResult.convergence_diagnostics.

Source code in src/metasignal/sdtbayes/diagnostics.py 
191
192
193
def convergence_diagnostics(fit: FitResult) -> Any:
    """Module-level alias for :meth:`FitResult.convergence_diagnostics`."""
    return fit.convergence_diagnostics()

Quick-start example (Option A — two-stage)

import numpy as np
from metasignal.sdtbayes import (
    fit_two_stage_group,
    posterior_summary,
    plot_posterior,
    convergence_diagnostics,
)

rng = np.random.default_rng(0)

# 20 participants: (stim, resp, conf) per participant
participants = [
    (rng.integers(0, 2, 200), rng.integers(0, 2, 200), rng.integers(1, 5, 200))
    for _ in range(20)
]

# Fit
fit, mle_df = fit_two_stage_group(participants, n_ratings=4)

# Check convergence
diag = convergence_diagnostics(fit)
print(diag[~diag["converged"]])   # empty if all converged

# Group mean M-ratio posterior (exp converts log scale → M-ratio)
print(posterior_summary(fit, var_names=["b_Intercept", "sigma"]))

# Plot
plot_posterior(fit, var_names=["b_Intercept"], ref_val=0)

Full HMeta-d example (Option B, cmdstanpy)

from metasignal.sdtbayes import fit_full_metad, plot_trace

fit = fit_full_metad(participants, n_ratings=4)

# Key parameters (mu_logMratio is an alias for alpha_logMratio)
print(posterior_summary(fit,
      var_names=["mu_logMratio", "sigma_logMratio", "mu_d1", "Mratio"]))

plot_trace(fit, var_names=["mu_logMratio", "sigma_logMratio"])

Group comparison example

from metasignal.sdtbayes import fit_full_metad_comparison
import arviz as az

fit_cmp = fit_full_metad_comparison(group_a, group_b, n_ratings=4)

# beta_logMratio[0] = group B - group A difference in log M-ratio
# (the covariate is mean-centred 0/1 -> -0.5/+0.5, spanning exactly one unit)
delta = az.extract(fit_cmp.idata)["beta_logMratio"].values[:, 0]
print(f"P(group B > group A): {(delta > 0).mean():.3f}")