Multilevel mediation API reference

This page exposes the complete generated reference for the post-parity multilevel mediation module. Use the method guide for workflow decisions and interpretation; use this page for signatures, parameters, return contracts, and docstrings.

multilevel_mediation

Bayesian multilevel mediation for trial-level gaze and behavioural data.

This module consumes the canonical preparation contract produced by eyeprocesspy.prepare_multilevel_mediation_data. It does not reimplement within/between decomposition.

The primary estimands are products of path coefficients on the model's linear predictor scale. For nonlinear outcome families these are not probability- scale natural indirect effects; the scale is therefore stored and reported explicitly.

create_mediation_prior_specification

create_mediation_prior_specification(*, intercept_sd: float = 2.5, coefficient_sd: float = 1.0, group_sd_scale: float = 1.0, residual_sd_scale: float = 1.0, dispersion_rate: float = 1.0, beta_precision_rate: float = 0.1) -> MediationPriorSpecification

Create explicit weakly-informative default priors for mediation paths.

Source code in src/gp3bayespy/multilevel_mediation.py

def create_mediation_prior_specification(
    *,
    intercept_sd: float = 2.5,
    coefficient_sd: float = 1.0,
    group_sd_scale: float = 1.0,
    residual_sd_scale: float = 1.0,
    dispersion_rate: float = 1.0,
    beta_precision_rate: float = 0.1,
) -> MediationPriorSpecification:
    """Create explicit weakly-informative default priors for mediation paths."""
    return MediationPriorSpecification(
        intercept_sd=_finite_positive(intercept_sd, "intercept_sd"),
        coefficient_sd=_finite_positive(coefficient_sd, "coefficient_sd"),
        group_sd_scale=_finite_positive(group_sd_scale, "group_sd_scale"),
        residual_sd_scale=_finite_positive(residual_sd_scale, "residual_sd_scale"),
        dispersion_rate=_finite_positive(dispersion_rate, "dispersion_rate"),
        beta_precision_rate=_finite_positive(beta_precision_rate, "beta_precision_rate"),
    )

specify_multilevel_gaze_mediation

specify_multilevel_gaze_mediation(prepared: Any, *, mediator_family: str = 'gaussian', outcome_family: str = 'bernoulli', priors: MediationPriorSpecification | None = None, random_slopes: Sequence[str] = ('mediator_x',), missingness_policy: str = 'error') -> MultilevelMediationSpecification

Create a backend-independent multilevel mediation specification.

missingness_policy='error' is intentionally conservative. Explicit complete_case or quality_eligible policies are available, and every excluded input row is recorded in the returned specification.

Source code in src/gp3bayespy/multilevel_mediation.py

def specify_multilevel_gaze_mediation(
    prepared: Any,
    *,
    mediator_family: str = "gaussian",
    outcome_family: str = "bernoulli",
    priors: MediationPriorSpecification | None = None,
    random_slopes: Sequence[str] = ("mediator_x",),
    missingness_policy: str = "error",
) -> MultilevelMediationSpecification:
    """Create a backend-independent multilevel mediation specification.

    ``missingness_policy='error'`` is intentionally conservative. Explicit
    ``complete_case`` or ``quality_eligible`` policies are available, and every
    excluded input row is recorded in the returned specification.
    """
    data, columns, provenance = _prepared_contract(prepared)
    mediator_family = _normalize_family(mediator_family, _MEDIATOR_FAMILIES, "mediator_family")
    outcome_family = _normalize_family(outcome_family, _OUTCOME_FAMILIES, "outcome_family")
    if missingness_policy not in _MISSING_POLICIES:
        raise GP3BayesError(
            "`missingness_policy` must be one of: " + ", ".join(sorted(_MISSING_POLICIES)) + "."
        )
    slopes = tuple(random_slopes)
    invalid_slopes = sorted(set(slopes) - _RANDOM_SLOPES)
    if invalid_slopes:
        raise GP3BayesError("Unknown random slope(s): " + ", ".join(invalid_slopes) + ".")
    if len(set(slopes)) != len(slopes):
        raise GP3BayesError("`random_slopes` must not contain duplicates.")
    priors = priors or create_mediation_prior_specification()
    if not isinstance(priors, MediationPriorSpecification):
        raise GP3BayesError("`priors` must be a MediationPriorSpecification.")

    mediator_col = str(columns["mediator"])
    outcome_col = str(columns["outcome"])
    participant_col = str(columns["participant"])
    trial_col = str(columns["trial"])
    required = [
        participant_col,
        trial_col,
        mediator_col,
        outcome_col,
        "X_within",
        "X_between",
        "M_within",
        "M_between",
    ]
    missing_cols = [c for c in required if c not in data]
    if missing_cols:
        raise GP3BayesError(
            "Prepared data is missing model columns: " + ", ".join(missing_cols) + "."
        )

    complete = data[required].notna().all(axis=1)
    eligible = data["mediation_analysis_eligible"].fillna(False).astype(bool)
    if missingness_policy == "error":
        if not bool((complete & eligible).all()):
            n_problem = int((~(complete & eligible)).sum())
            raise GP3BayesError(
                f"{n_problem} trial row(s) are incomplete or quality-ineligible. "
                "Choose an explicit `missingness_policy` after reviewing the preparation audits."
            )
        keep = pd.Series(True, index=data.index)
    elif missingness_policy == "complete_case":
        keep = complete
    else:
        keep = complete & eligible

    excluded_positions = tuple(np.flatnonzero(~keep.to_numpy()).astype(int).tolist())
    analysis = data.loc[keep].copy()
    if analysis.empty:
        raise GP3BayesError("No rows remain under the requested missingness policy.")
    if analysis[participant_col].nunique() < 2:
        raise GP3BayesError("At least two participants are required for multilevel mediation.")
    _validate_family_data(analysis, mediator_col, outcome_col, mediator_family, outcome_family)
    estimable_paths = _estimable_simple_paths(analysis)
    required_indirect = {"a_within", "b_within"}
    if not required_indirect.issubset(estimable_paths):
        missing_paths = sorted(required_indirect - set(estimable_paths))
        raise GP3BayesError(
            "The within-participant indirect effect is not estimable from the observed design; "
            "missing variation for path(s): " + ", ".join(missing_paths) + "."
        )
    slope_path = {"mediator_x": "a_within", "outcome_x": "cprime_within", "outcome_m": "b_within"}
    unsupported_slopes = [name for name in slopes if slope_path[name] not in estimable_paths]
    if unsupported_slopes:
        raise GP3BayesError(
            "Requested random slope(s) are not estimable from observed within-participant variation: "
            + ", ".join(unsupported_slopes)
            + "."
        )

    # Trial-level random slopes need repeated observations within at least some groups.
    counts = analysis.groupby(participant_col)[trial_col].nunique()
    if slopes and int((counts >= 2).sum()) < 2:
        raise GP3BayesError("Random slopes require repeated trials for at least two participants.")

    model_provenance = {
        "preparation": dict(provenance),
        "model_specification": {
            "mediator_family": mediator_family,
            "outcome_family": outcome_family,
            "random_slopes": list(slopes),
            "missingness_policy": missingness_policy,
            "priors": priors.as_dict(),
            "estimand_scale": "linear_predictor_product",
            "estimable_paths": list(estimable_paths),
        },
    }
    return MultilevelMediationSpecification(
        data=analysis,
        participant_col=participant_col,
        trial_col=trial_col,
        outcome_col=outcome_col,
        mediator_col=mediator_col,
        mediator_family=mediator_family,
        outcome_family=outcome_family,
        priors=priors,
        random_slopes=slopes,
        missingness_policy=missingness_policy,
        input_rows=len(data),
        analysis_rows=len(analysis),
        excluded_row_positions=excluded_positions,
        provenance=model_provenance,
        estimable_paths=estimable_paths,
    )

fit_multilevel_gaze_mediation

fit_multilevel_gaze_mediation(prepared: Any, *, mediator_family: str = 'gaussian', outcome_family: str = 'bernoulli', priors: MediationPriorSpecification | None = None, random_slopes: Sequence[str] = ('mediator_x',), missingness_policy: str = 'error', chains: int = 4, draws: int = 1000, tune: int = 1000, cores: int | None = None, seed: int = 1, target_accept: float = 0.95, max_treedepth: int = 12) -> MultilevelMediationFit

Fit a Bayesian trial-level multilevel mediation model with PyMC NUTS.

Source code in src/gp3bayespy/multilevel_mediation.py

def fit_multilevel_gaze_mediation(
    prepared: Any,
    *,
    mediator_family: str = "gaussian",
    outcome_family: str = "bernoulli",
    priors: MediationPriorSpecification | None = None,
    random_slopes: Sequence[str] = ("mediator_x",),
    missingness_policy: str = "error",
    chains: int = 4,
    draws: int = 1000,
    tune: int = 1000,
    cores: int | None = None,
    seed: int = 1,
    target_accept: float = 0.95,
    max_treedepth: int = 12,
) -> MultilevelMediationFit:
    """Fit a Bayesian trial-level multilevel mediation model with PyMC NUTS."""
    spec = specify_multilevel_gaze_mediation(
        prepared,
        mediator_family=mediator_family,
        outcome_family=outcome_family,
        priors=priors,
        random_slopes=random_slopes,
        missingness_policy=missingness_policy,
    )
    controls = _validate_sampling_controls(
        chains=chains,
        draws=draws,
        tune=tune,
        cores=cores,
        seed=seed,
        target_accept=target_accept,
        max_treedepth=max_treedepth,
    )

    pm = _load_pymc()
    model = _build_pymc_model(spec)
    with model:
        idata = pm.sample(
            draws=controls["draws"],
            tune=controls["tune"],
            chains=controls["chains"],
            cores=controls["cores"],
            random_seed=controls["seed"],
            target_accept=controls["target_accept"],
            nuts={"max_treedepth": controls["max_treedepth"]},
            progressbar=False,
            compute_convergence_checks=False,
            return_inferencedata=True,
            idata_kwargs={"log_likelihood": True},
        )
    names = [
        "a_within",
        "a_between",
        "b_within",
        "b_between",
        "cprime_within",
        "cprime_between",
        "indirect_within",
        "indirect_between",
        "total_within",
        "total_between",
        "participant_a_slope",
        "participant_b_slope",
        "participant_indirect_within",
    ]
    posterior = _flatten_posterior(idata, names)
    diagnostics = _sampler_diagnostics(idata, max_treedepth=int(controls["max_treedepth"]))
    fit = MultilevelMediationFit(
        specification=spec,
        posterior=posterior,
        sampler_diagnostics=diagnostics,
        backend_fit=idata,
        backend_model=model,
        package_versions={"pymc": _package_version("pymc"), "arviz": _package_version("arviz")},
    )
    convergence = check_mediation_convergence(fit)
    if not convergence.passed:
        warnings.warn(
            "Mediation fit failed one or more critical convergence checks. "
            "Indirect-effect extraction will fail by default until diagnostics are resolved.",
            RuntimeWarning,
            stacklevel=2,
        )
    return fit

check_mediation_convergence

check_mediation_convergence(fit: MultilevelMediationFit, *, rhat_max: float = 1.01, ess_bulk_min: float = 400.0, divergences_max: int = 0, treedepth_hits_max: int = 0) -> MediationConvergence

Apply explicit convergence thresholds to a mediation fit.

Source code in src/gp3bayespy/multilevel_mediation.py

def check_mediation_convergence(
    fit: MultilevelMediationFit,
    *,
    rhat_max: float = 1.01,
    ess_bulk_min: float = 400.0,
    divergences_max: int = 0,
    treedepth_hits_max: int = 0,
) -> MediationConvergence:
    """Apply explicit convergence thresholds to a mediation fit."""
    if not isinstance(fit, MultilevelMediationFit):
        raise GP3BayesError("`fit` must be a MultilevelMediationFit.")
    d = fit.sampler_diagnostics
    max_rhat = float(d.get("max_rhat", np.nan))
    min_ess = float(d.get("min_ess_bulk", np.nan))
    divergences = int(d.get("divergences", 0))
    treedepth_hits = int(d.get("treedepth_hits", 0))
    issues: list[str] = []
    if not np.isfinite(max_rhat):
        issues.append("R-hat was not available")
    elif max_rhat > rhat_max:
        issues.append(f"max R-hat {max_rhat:.3f} exceeds {rhat_max:.3f}")
    if not np.isfinite(min_ess):
        issues.append("bulk ESS was not available")
    elif min_ess < ess_bulk_min:
        issues.append(f"minimum bulk ESS {min_ess:.0f} is below {ess_bulk_min:.0f}")
    if divergences > divergences_max:
        issues.append(f"{divergences} divergence(s) detected")
    if treedepth_hits > treedepth_hits_max:
        issues.append(f"{treedepth_hits} maximum-tree-depth hit(s) detected")
    passed = len(issues) == 0
    return MediationConvergence(
        status="pass" if passed else "fail",
        passed=passed,
        max_rhat=max_rhat,
        min_ess_bulk=min_ess,
        divergences=divergences,
        treedepth_hits=treedepth_hits,
        thresholds={
            "rhat_max": float(rhat_max),
            "ess_bulk_min": float(ess_bulk_min),
            "divergences_max": float(divergences_max),
            "treedepth_hits_max": float(treedepth_hits_max),
        },
        issues=tuple(issues),
    )

posterior_indirect_effect

posterior_indirect_effect(fit: MultilevelMediationFit, *, level: str = 'within', probability: float = 0.95, require_convergence: bool = True) -> MediationEffect

Return posterior indirect-effect draws on the linear-predictor product scale.

Source code in src/gp3bayespy/multilevel_mediation.py

def posterior_indirect_effect(
    fit: MultilevelMediationFit,
    *,
    level: str = "within",
    probability: float = 0.95,
    require_convergence: bool = True,
) -> MediationEffect:
    """Return posterior indirect-effect draws on the linear-predictor product scale."""
    _require_acceptable_fit(fit, require_convergence)
    if level not in {"within", "between"}:
        raise GP3BayesError("`level` must be `within` or `between`.")
    key = f"indirect_{level}"
    draws = fit.posterior.get(key)
    if draws is None:
        a = fit.posterior.get(f"a_{level}")
        b = fit.posterior.get(f"b_{level}")
        if a is None or b is None:
            raise GP3BayesError(f"Posterior path draws required for `{key}` are unavailable.")
        draws = np.asarray(a) * np.asarray(b)
    return _effect_from_draws(key, draws, probability, "linear_predictor_product")

summarise_multilevel_mediation

summarise_multilevel_mediation(fit: MultilevelMediationFit, *, probability: float = 0.95, require_convergence: bool = True) -> pd.DataFrame

Summarise within/between a, b, direct, indirect, and total effects.

Source code in src/gp3bayespy/multilevel_mediation.py

def summarise_multilevel_mediation(
    fit: MultilevelMediationFit,
    *,
    probability: float = 0.95,
    require_convergence: bool = True,
) -> pd.DataFrame:
    """Summarise within/between a, b, direct, indirect, and total effects."""
    _require_acceptable_fit(fit, require_convergence)
    rows: list[dict[str, Any]] = []
    for level in ("within", "between"):
        for name, scale in [
            (f"a_{level}", "linear_predictor"),
            (f"b_{level}", "linear_predictor"),
            (f"cprime_{level}", "linear_predictor"),
        ]:
            if name in fit.posterior:
                effect = _effect_from_draws(name, fit.posterior[name], probability, scale)
                rows.append(_effect_row(effect))
        with suppress(GP3BayesError):
            rows.append(
                _effect_row(
                    posterior_indirect_effect(
                        fit,
                        level=level,
                        probability=probability,
                        require_convergence=False,
                    )
                )
            )
        with suppress(GP3BayesError):
            rows.append(
                _effect_row(
                    posterior_total_effect(
                        fit,
                        level=level,
                        probability=probability,
                        require_convergence=False,
                    )
                )
            )
    return pd.DataFrame(rows)

posterior_predictive_check_mediation

posterior_predictive_check_mediation(fit: MultilevelMediationFit, *, draws: int = 200, seed: int = 1) -> pd.DataFrame

Run simple posterior-predictive checks for mediator and outcome.

Source code in src/gp3bayespy/multilevel_mediation.py

def posterior_predictive_check_mediation(
    fit: MultilevelMediationFit,
    *,
    draws: int = 200,
    seed: int = 1,
) -> pd.DataFrame:
    """Run simple posterior-predictive checks for mediator and outcome."""
    draws, seed = _validate_draw_seed(draws, seed)
    if fit.backend_model is None or fit.backend_fit is None:
        raise GP3BayesError("Posterior predictive checks require a fitted backend model.")
    pm = _load_pymc()
    with fit.backend_model:
        ppc = pm.sample_posterior_predictive(
            fit.backend_fit,
            var_names=["M_obs", "Y_obs"],
            random_seed=seed,
            progressbar=False,
            predictions=False,
        )
    predictive = getattr(ppc, "posterior_predictive", None)
    if predictive is None:
        raise GP3BayesError("Backend did not return posterior predictive draws.")
    rows: list[dict[str, Any]] = []
    for var, observed_col in [
        ("M_obs", fit.specification.mediator_col),
        ("Y_obs", fit.specification.outcome_col),
    ]:
        sims = np.asarray(predictive[var]).reshape(-1, fit.specification.analysis_rows)
        if len(sims) > draws:
            sims = sims[:draws]
        observed = pd.to_numeric(fit.specification.data[observed_col], errors="raise").to_numpy(
            dtype=float
        )
        sim_means = sims.mean(axis=1)
        obs_mean = float(np.mean(observed))
        rows.append(
            {
                "variable": var,
                "observed_mean": obs_mean,
                "predictive_mean": float(np.mean(sim_means)),
                "predictive_mean_lower": float(np.quantile(sim_means, 0.025)),
                "predictive_mean_upper": float(np.quantile(sim_means, 0.975)),
                "tail_probability": float(np.mean(sim_means >= obs_mean)),
            }
        )
    return pd.DataFrame(rows)

prior_predictive_check_mediation

prior_predictive_check_mediation(specification: MultilevelMediationSpecification, *, draws: int = 200, seed: int = 1) -> pd.DataFrame

Sample prior predictive mediator/outcome distributions from a specification.

Source code in src/gp3bayespy/multilevel_mediation.py

def prior_predictive_check_mediation(
    specification: MultilevelMediationSpecification,
    *,
    draws: int = 200,
    seed: int = 1,
) -> pd.DataFrame:
    """Sample prior predictive mediator/outcome distributions from a specification."""
    if not isinstance(specification, MultilevelMediationSpecification):
        raise GP3BayesError("`specification` must be a MultilevelMediationSpecification.")
    draws, seed = _validate_draw_seed(draws, seed)
    pm = _load_pymc()
    model = _build_pymc_model(specification)
    with model:
        prior = pm.sample_prior_predictive(draws=draws, random_seed=seed)
    group = getattr(prior, "prior_predictive", None)
    if group is None:
        raise GP3BayesError("Backend did not return prior predictive draws.")
    rows = []
    for var in ("M_obs", "Y_obs"):
        sims = np.asarray(group[var]).reshape(-1, specification.analysis_rows)
        rows.append(
            {
                "variable": var,
                "draws": len(sims),
                "mean_of_means": float(np.mean(sims.mean(axis=1))),
                "minimum": float(np.min(sims)),
                "maximum": float(np.max(sims)),
            }
        )
    return pd.DataFrame(rows)

compare_multilevel_mediation_models

compare_multilevel_mediation_models(fits: Mapping[str, MultilevelMediationFit]) -> pd.DataFrame

Compare mediation fits using joint PSIS-LOO on identical observations only.

Source code in src/gp3bayespy/multilevel_mediation.py

def compare_multilevel_mediation_models(fits: Mapping[str, MultilevelMediationFit]) -> pd.DataFrame:
    """Compare mediation fits using joint PSIS-LOO on identical observations only."""
    if not isinstance(fits, Mapping) or len(fits) < 2:
        raise GP3BayesError("`fits` must contain at least two named mediation fits.")
    fit_items = list(fits.items())
    for label, fit in fit_items:
        if not isinstance(fit, MultilevelMediationFit) or fit.backend_fit is None:
            raise GP3BayesError(f"Model `{label}` is not a fitted mediation object.")
    reference_label, reference_fit = fit_items[0]
    for label, fit in fit_items[1:]:
        if not _same_comparison_observations(reference_fit.specification, fit.specification):
            raise GP3BayesError(
                "PSIS-LOO comparison requires the same mediator/outcome observations in the same "
                f"participant-trial order; `{label}` does not match `{reference_label}`. "
                "Do not compare fits produced from different missingness/exclusion sets."
            )
    try:
        az = import_module("arviz")
    except Exception as exc:
        raise BackendUnavailableError("ArviZ is required for PSIS-LOO model comparison.") from exc
    rows: list[dict[str, Any]] = []
    for label, fit in fit_items:
        ll = getattr(fit.backend_fit, "log_likelihood", None)
        if ll is None or "M_obs" not in ll or "Y_obs" not in ll:
            raise GP3BayesError(f"Model `{label}` lacks mediator/outcome log-likelihood draws.")
        joint = _joint_pointwise_log_likelihood(
            ll["M_obs"], ll["Y_obs"], expected_rows=fit.specification.analysis_rows
        )
        temp = fit.backend_fit.copy()
        try:
            temp.log_likelihood["joint"] = joint
            loo = az.loo(temp, var_name="joint")
        except Exception as exc:
            raise GP3BayesError(f"Could not compute joint PSIS-LOO for `{label}`: {exc}") from exc
        rows.append(
            {
                "model": label,
                "elpd_loo": float(loo.elpd_loo),
                "se": float(loo.se),
                "p_loo": float(loo.p_loo),
                "warning": bool(getattr(loo, "warning", False)),
            }
        )
    return pd.DataFrame(rows).sort_values("elpd_loo", ascending=False).reset_index(drop=True)

plot_indirect_effect_distribution

plot_indirect_effect_distribution(fit: MultilevelMediationFit, *, level: str = 'within', probability: float = 0.95, require_convergence: bool = True, ax: Any = None) -> Any

Plot the posterior distribution of an indirect effect.

Source code in src/gp3bayespy/multilevel_mediation.py

def plot_indirect_effect_distribution(
    fit: MultilevelMediationFit,
    *,
    level: str = "within",
    probability: float = 0.95,
    require_convergence: bool = True,
    ax: Any = None,
) -> Any:
    """Plot the posterior distribution of an indirect effect."""
    effect = posterior_indirect_effect(
        fit, level=level, probability=probability, require_convergence=require_convergence
    )
    try:
        plt = import_module("matplotlib.pyplot")
    except Exception as exc:
        raise BackendUnavailableError("Matplotlib is required for mediation plots.") from exc
    if ax is None:
        _, ax = plt.subplots()
    ax.hist(effect.draws, bins=30, density=True)
    ax.axvline(0.0, linewidth=1)
    ax.axvline(effect.mean, linewidth=1.5)
    ax.set_xlabel(f"{level.title()} indirect effect ({effect.scale})")
    ax.set_ylabel("Posterior density")
    ax.set_title(f"Posterior {level} indirect effect")
    return ax

plot_mediation_posteriors

plot_mediation_posteriors(fit: MultilevelMediationFit, *, probability: float = 0.95, require_convergence: bool = True, ax: Any = None) -> Any

Plot posterior interval summaries for mediation paths and effects.

Source code in src/gp3bayespy/multilevel_mediation.py

def plot_mediation_posteriors(
    fit: MultilevelMediationFit,
    *,
    probability: float = 0.95,
    require_convergence: bool = True,
    ax: Any = None,
) -> Any:
    """Plot posterior interval summaries for mediation paths and effects."""
    summary = summarise_multilevel_mediation(
        fit, probability=probability, require_convergence=require_convergence
    )
    try:
        plt = import_module("matplotlib.pyplot")
    except Exception as exc:
        raise BackendUnavailableError("Matplotlib is required for mediation plots.") from exc
    if ax is None:
        _, ax = plt.subplots()
    y = np.arange(len(summary))
    ax.errorbar(
        summary["mean"],
        y,
        xerr=[summary["mean"] - summary["lower"], summary["upper"] - summary["mean"]],
        fmt="o",
    )
    ax.axvline(0.0, linewidth=1)
    ax.set_yticks(y, summary["effect"])
    ax.set_xlabel("Posterior effect")
    ax.set_title("Multilevel mediation posterior intervals")
    return ax

plot_participant_mediation_effects

plot_participant_mediation_effects(fit: MultilevelMediationFit, *, ax: Any = None) -> Any

Plot participant random-effect scales when available.

Participant-specific indirect effects are not manufactured when the fitted model did not include both participant-specific a and b paths.

Source code in src/gp3bayespy/multilevel_mediation.py

def plot_participant_mediation_effects(fit: MultilevelMediationFit, *, ax: Any = None) -> Any:
    """Plot participant random-effect scales when available.

    Participant-specific indirect effects are not manufactured when the fitted
    model did not include both participant-specific a and b paths.
    """
    keys = [k for k in fit.posterior if k.startswith("participant_indirect_")]
    if not keys:
        raise GP3BayesError(
            "Participant-specific indirect effects are unavailable because the model did not estimate "
            "participant-specific a and b paths jointly."
        )
    try:
        plt = import_module("matplotlib.pyplot")
    except Exception as exc:
        raise BackendUnavailableError("Matplotlib is required for mediation plots.") from exc
    if ax is None:
        _, ax = plt.subplots()
    values = np.asarray(fit.posterior[keys[0]])
    means = values.mean(axis=0)
    ax.plot(np.arange(len(means)), means, "o")
    ax.axhline(0.0, linewidth=1)
    ax.set_xlabel("Participant")
    ax.set_ylabel("Posterior mean indirect effect")
    return ax

report_multilevel_gaze_mediation

report_multilevel_gaze_mediation(fit: MultilevelMediationFit, *, probability: float = 0.95, require_convergence: bool = True) -> str

Create a compact, conservative reporting paragraph.

Source code in src/gp3bayespy/multilevel_mediation.py

def report_multilevel_gaze_mediation(
    fit: MultilevelMediationFit,
    *,
    probability: float = 0.95,
    require_convergence: bool = True,
) -> str:
    """Create a compact, conservative reporting paragraph."""
    convergence = check_mediation_convergence(fit)
    if require_convergence and not convergence.passed:
        raise GP3BayesError("A substantive report is blocked because convergence checks failed.")
    within = posterior_indirect_effect(
        fit, level="within", probability=probability, require_convergence=False
    )
    try:
        between = posterior_indirect_effect(
            fit, level="between", probability=probability, require_convergence=False
        )
        between_text = f"{between.mean:.3f} [{between.lower:.3f}, {between.upper:.3f}]"
    except GP3BayesError:
        between_text = "not estimable from the observed design"
    spec = fit.specification
    excluded = len(spec.excluded_row_positions)
    return (
        f"A Bayesian trial-level multilevel mediation model was fit to {spec.analysis_rows} observations "
        f"from {spec.data[spec.participant_col].nunique()} participants, using a {spec.mediator_family} "
        f"mediator model and {spec.outcome_family} outcome model. The explicit missingness policy was "
        f"'{spec.missingness_policy}' ({excluded} input rows excluded). The within-participant indirect "
        f"effect on the linear-predictor product scale was {within.mean:.3f} "
        f"[{within.lower:.3f}, {within.upper:.3f}], and the between-participant indirect effect was "
        f"{between_text}. Sampler diagnostics status: "
        f"{convergence.status}. For nonlinear families these coefficient products should not be interpreted "
        f"as probability-scale natural indirect effects."
    )

simulate_multilevel_gaze_mediation

simulate_multilevel_gaze_mediation(*, n_participants: int = 100, trials_per_participant: int = 20, a_within: float = 0.7, b_within: float = 0.8, cprime_within: float = 0.2, a_between: float = 0.2, b_between: float = 0.3, seed: int = 2026) -> pd.DataFrame

Simulate the repeated AI-advice example used by tests and documentation.

Source code in src/gp3bayespy/multilevel_mediation.py

def simulate_multilevel_gaze_mediation(
    *,
    n_participants: int = 100,
    trials_per_participant: int = 20,
    a_within: float = 0.7,
    b_within: float = 0.8,
    cprime_within: float = 0.2,
    a_between: float = 0.2,
    b_between: float = 0.3,
    seed: int = 2026,
) -> pd.DataFrame:
    """Simulate the repeated AI-advice example used by tests and documentation."""
    if n_participants < 2 or trials_per_participant < 2:
        raise GP3BayesError(
            "Simulation requires at least two participants and two trials per participant."
        )
    rng = np.random.default_rng(seed)
    rows: list[dict[str, Any]] = []
    participant_intercept_m = rng.normal(0, 0.5, n_participants)
    participant_intercept_y = rng.normal(0, 0.5, n_participants)
    participant_x_propensity = rng.beta(4, 4, n_participants)
    for i in range(n_participants):
        x = rng.binomial(1, participant_x_propensity[i], trials_per_participant)
        x_between = float(np.mean(x))
        x_within = x - x_between
        mediator = (
            1.5
            + a_within * x_within
            + a_between * x_between
            + participant_intercept_m[i]
            + rng.normal(0, 0.6, trials_per_participant)
        )
        mediator = np.maximum(mediator, 0.02)
        m_between = float(np.mean(mediator))
        m_within = mediator - m_between
        eta = (
            -0.5
            + cprime_within * x_within
            + b_within * m_within
            + b_between * m_between
            + participant_intercept_y[i]
        )
        p = 1.0 / (1.0 + np.exp(-eta))
        outcome = rng.binomial(1, p)
        for j in range(trials_per_participant):
            rows.append(
                {
                    "participant_id": f"p{i + 1:03d}",
                    "trial_id": j + 1,
                    "ai_correct": int(x[j]),
                    "source_dwell": float(mediator[j]),
                    "correct_override": int(outcome[j]),
                    "valid_fraction": 1.0,
                }
            )
    return pd.DataFrame(rows)

specify_multilevel_serial_gaze_mediation

specify_multilevel_serial_gaze_mediation(prepared: Any, *, mediator2_family: str = 'gaussian', mediator_family: str = 'gaussian', outcome_family: str = 'bernoulli', priors: MediationPriorSpecification | None = None, random_slopes: Sequence[str] = (), missingness_policy: str = 'error') -> MultilevelMediationSpecification

Specify X -> M1 -> M2 -> Y using eyeprocess-prepared components.

Serial models currently use participant random intercepts. Non-empty random_slopes are rejected explicitly rather than ignored.

Source code in src/gp3bayespy/multilevel_mediation.py

def specify_multilevel_serial_gaze_mediation(
    prepared: Any,
    *,
    mediator2_family: str = "gaussian",
    mediator_family: str = "gaussian",
    outcome_family: str = "bernoulli",
    priors: MediationPriorSpecification | None = None,
    random_slopes: Sequence[str] = (),
    missingness_policy: str = "error",
) -> MultilevelMediationSpecification:
    """Specify X -> M1 -> M2 -> Y using eyeprocess-prepared components.

    Serial models currently use participant random intercepts. Non-empty
    ``random_slopes`` are rejected explicitly rather than ignored.
    """
    data, columns, provenance = _prepared_contract(prepared)
    for key in ("mediator2", "mediator2_within", "mediator2_between"):
        if key not in columns:
            raise GP3BayesError(
                "Serial mediation requires eyeprocess-prepared `mediator2`, "
                "`mediator2_within`, and `mediator2_between` semantic columns."
            )
    slopes = tuple(random_slopes)
    if slopes:
        raise GP3BayesError(
            "Serial mediation random slopes are not yet implemented; pass `random_slopes=()` "
            "instead of requesting slopes that would be silently ignored."
        )
    m1_family = _normalize_family(mediator_family, _MEDIATOR_FAMILIES, "mediator_family")
    m2_family = _normalize_family(mediator2_family, _MEDIATOR_FAMILIES, "mediator2_family")
    y_family = _normalize_family(outcome_family, _OUTCOME_FAMILIES, "outcome_family")
    priors = priors or create_mediation_prior_specification()
    if not isinstance(priors, MediationPriorSpecification):
        raise GP3BayesError("`priors` must be a MediationPriorSpecification.")

    participant_col = str(columns["participant"])
    trial_col = str(columns["trial"])
    mediator_col = str(columns["mediator"])
    outcome_col = str(columns["outcome"])
    mediator2_col = str(columns["mediator2"])
    m2w = str(columns["mediator2_within"])
    m2b = str(columns["mediator2_between"])
    required = [
        participant_col,
        trial_col,
        mediator_col,
        mediator2_col,
        outcome_col,
        "X_within",
        "X_between",
        "M_within",
        "M_between",
        m2w,
        m2b,
    ]
    missing_cols = [c for c in required if c not in data]
    if missing_cols:
        raise GP3BayesError(
            "Prepared serial data is missing model columns: " + ", ".join(missing_cols) + "."
        )
    analysis, excluded = _select_analysis_rows(
        data, required=required, missingness_policy=missingness_policy
    )
    if analysis[participant_col].nunique() < 2:
        raise GP3BayesError(
            "At least two participants are required for serial multilevel mediation."
        )
    _validate_family_data(analysis, mediator_col, outcome_col, m1_family, y_family)
    temp = analysis.rename(columns={mediator2_col: "__m2__", outcome_col: "__y__"})
    _validate_family_data(temp, "__m2__", "__y__", m2_family, y_family)

    estimable: list[str] = []
    if _has_variation(analysis["X_within"]):
        estimable += ["a1_within", "a2_within", "cprime_within"]
    if _has_variation(analysis["X_between"]):
        estimable += ["a1_between", "a2_between", "cprime_between"]
    if _has_variation(analysis["M_within"]):
        estimable += ["d_within", "b1_within"]
    if _has_variation(analysis["M_between"]):
        estimable += ["d_between", "b1_between"]
    if _has_variation(analysis[m2w]):
        estimable.append("b2_within")
    if _has_variation(analysis[m2b]):
        estimable.append("b2_between")
    estimable_paths = tuple(estimable)
    required_serial = {"a1_within", "d_within", "b2_within"}
    if not required_serial.issubset(estimable_paths):
        raise GP3BayesError(
            "The within-participant serial indirect effect is not estimable from the observed design; "
            "missing path(s): " + ", ".join(sorted(required_serial - set(estimable_paths))) + "."
        )

    model_provenance = {
        "preparation": dict(provenance),
        "model_specification": {
            "model_kind": "serial",
            "mediator_family": m1_family,
            "mediator2_family": m2_family,
            "outcome_family": y_family,
            "random_slopes": [],
            "missingness_policy": missingness_policy,
            "priors": priors.as_dict(),
            "estimand_scale": "linear_predictor_product",
            "estimable_paths": list(estimable_paths),
            "random_effects_scope": "participant_random_intercepts",
        },
    }
    return MultilevelMediationSpecification(
        data=analysis,
        participant_col=participant_col,
        trial_col=trial_col,
        outcome_col=outcome_col,
        mediator_col=mediator_col,
        mediator_family=m1_family,
        outcome_family=y_family,
        priors=priors,
        random_slopes=(),
        missingness_policy=missingness_policy,
        input_rows=len(data),
        analysis_rows=len(analysis),
        excluded_row_positions=excluded,
        provenance=model_provenance,
        model_kind="serial",
        extra_columns={
            "mediator2": mediator2_col,
            "mediator2_within": m2w,
            "mediator2_between": m2b,
            "mediator2_family": m2_family,
        },
        serial=True,
        estimable_paths=estimable_paths,
    )

fit_multilevel_serial_gaze_mediation

fit_multilevel_serial_gaze_mediation(prepared: Any, *, mediator2_family: str = 'gaussian', mediator_family: str = 'gaussian', outcome_family: str = 'bernoulli', priors: MediationPriorSpecification | None = None, random_slopes: Sequence[str] = (), missingness_policy: str = 'error', chains: int = 4, draws: int = 1000, tune: int = 1000, cores: int | None = None, seed: int = 1, target_accept: float = 0.95, max_treedepth: int = 12) -> MultilevelMediationFit

Fit X -> M1 -> M2 -> Y without re-decomposing variables in gp3bayespy.

Source code in src/gp3bayespy/multilevel_mediation.py

def fit_multilevel_serial_gaze_mediation(
    prepared: Any,
    *,
    mediator2_family: str = "gaussian",
    mediator_family: str = "gaussian",
    outcome_family: str = "bernoulli",
    priors: MediationPriorSpecification | None = None,
    random_slopes: Sequence[str] = (),
    missingness_policy: str = "error",
    chains: int = 4,
    draws: int = 1000,
    tune: int = 1000,
    cores: int | None = None,
    seed: int = 1,
    target_accept: float = 0.95,
    max_treedepth: int = 12,
) -> MultilevelMediationFit:
    """Fit X -> M1 -> M2 -> Y without re-decomposing variables in gp3bayespy."""
    spec = specify_multilevel_serial_gaze_mediation(
        prepared,
        mediator2_family=mediator2_family,
        mediator_family=mediator_family,
        outcome_family=outcome_family,
        priors=priors,
        random_slopes=random_slopes,
        missingness_policy=missingness_policy,
    )
    controls = _validate_sampling_controls(
        chains=chains,
        draws=draws,
        tune=tune,
        cores=cores,
        seed=seed,
        target_accept=target_accept,
        max_treedepth=max_treedepth,
    )
    pm = _load_pymc()
    model = _build_serial_pymc_model(spec)
    with model:
        idata = pm.sample(
            draws=controls["draws"],
            tune=controls["tune"],
            chains=controls["chains"],
            cores=controls["cores"],
            random_seed=controls["seed"],
            target_accept=controls["target_accept"],
            nuts={"max_treedepth": controls["max_treedepth"]},
            progressbar=False,
            compute_convergence_checks=False,
            return_inferencedata=True,
            idata_kwargs={"log_likelihood": True},
        )
    names = [
        "a1_within",
        "a1_between",
        "a2_within",
        "a2_between",
        "d_within",
        "d_between",
        "b1_within",
        "b1_between",
        "b2_within",
        "b2_between",
        "cprime_within",
        "cprime_between",
        "serial_indirect_within",
        "serial_indirect_between",
        "m1_indirect_within",
        "m2_indirect_within",
        "total_indirect_within",
        "total_within",
    ]
    fit = MultilevelMediationFit(
        specification=spec,
        posterior=_flatten_posterior(idata, names),
        sampler_diagnostics=_sampler_diagnostics(
            idata, max_treedepth=int(controls["max_treedepth"])
        ),
        backend_fit=idata,
        backend_model=model,
        package_versions={"pymc": _package_version("pymc"), "arviz": _package_version("arviz")},
    )
    if not check_mediation_convergence(fit).passed:
        warnings.warn(
            "Serial mediation fit failed critical convergence checks.", RuntimeWarning, stacklevel=2
        )
    return fit

specify_multilevel_moderated_gaze_mediation

specify_multilevel_moderated_gaze_mediation(prepared: Any, *, moderation_path: str = 'a', moderator_component: str = 'within', mediator_family: str = 'gaussian', outcome_family: str = 'bernoulli', priors: MediationPriorSpecification | None = None, random_slopes: Sequence[str] = (), missingness_policy: str = 'error') -> MultilevelMediationSpecification

Specify moderation of the within a or b path with a prepared moderator.

Moderated extensions currently use participant random intercepts. Requested random slopes are rejected explicitly until the extension implements them.

Source code in src/gp3bayespy/multilevel_mediation.py

def specify_multilevel_moderated_gaze_mediation(
    prepared: Any,
    *,
    moderation_path: str = "a",
    moderator_component: str = "within",
    mediator_family: str = "gaussian",
    outcome_family: str = "bernoulli",
    priors: MediationPriorSpecification | None = None,
    random_slopes: Sequence[str] = (),
    missingness_policy: str = "error",
) -> MultilevelMediationSpecification:
    """Specify moderation of the within a or b path with a prepared moderator.

    Moderated extensions currently use participant random intercepts. Requested
    random slopes are rejected explicitly until the extension implements them.
    """
    data, columns, _ = _prepared_contract(prepared)
    if moderation_path not in {"a", "b"}:
        raise GP3BayesError("`moderation_path` must be `a` or `b`.")
    if moderator_component not in {"within", "between"}:
        raise GP3BayesError("`moderator_component` must be `within` or `between`.")
    if tuple(random_slopes):
        raise GP3BayesError(
            "Moderated mediation random slopes are not yet implemented; pass `random_slopes=()` "
            "instead of requesting slopes that would be silently ignored."
        )
    key = f"moderator_{moderator_component}"
    if key not in columns:
        raise GP3BayesError(f"Moderated mediation requires eyeprocess-prepared `{key}` semantics.")
    base = specify_multilevel_gaze_mediation(
        prepared,
        mediator_family=mediator_family,
        outcome_family=outcome_family,
        priors=priors,
        random_slopes=(),
        missingness_policy=missingness_policy,
    )
    zcol = str(columns[key])
    if zcol not in data:
        raise GP3BayesError(f"Prepared moderator column `{zcol}` is missing.")
    z_missing = data[zcol].isna()
    if bool(z_missing.any()) and missingness_policy == "error":
        raise GP3BayesError(
            "Moderator contains missing values; choose an explicit missingness policy."
        )
    analysis = base.data.loc[base.data[zcol].notna()].copy()
    excluded_extra = tuple(np.flatnonzero(z_missing.to_numpy()).astype(int).tolist())
    excluded = tuple(sorted(set(base.excluded_row_positions) | set(excluded_extra)))
    if analysis.empty:
        raise GP3BayesError("No rows remain after applying moderator missingness rules.")
    if not _has_variation(analysis[zcol]):
        raise GP3BayesError("The selected moderator component has no observed variation.")
    interaction = (
        analysis["X_within"] * pd.to_numeric(analysis[zcol], errors="raise")
        if moderation_path == "a"
        else analysis["M_within"] * pd.to_numeric(analysis[zcol], errors="raise")
    )
    if not _has_variation(interaction):
        raise GP3BayesError("The requested moderated path interaction has no observed variation.")
    estimable_paths = _estimable_simple_paths(analysis)
    if not {"a_within", "b_within"}.issubset(estimable_paths):
        raise GP3BayesError(
            "The within-participant indirect effect is not estimable after moderator filtering."
        )
    provenance2 = dict(base.provenance)
    provenance2["model_specification"] = dict(provenance2["model_specification"])
    provenance2["model_specification"].update(
        {
            "model_kind": "moderated",
            "moderation_path": moderation_path,
            "moderator_component": moderator_component,
            "random_effects_scope": "participant_random_intercepts",
            "estimable_paths": list(estimable_paths),
        }
    )
    return MultilevelMediationSpecification(
        data=analysis,
        participant_col=base.participant_col,
        trial_col=base.trial_col,
        outcome_col=base.outcome_col,
        mediator_col=base.mediator_col,
        mediator_family=base.mediator_family,
        outcome_family=base.outcome_family,
        priors=base.priors,
        random_slopes=(),
        missingness_policy=base.missingness_policy,
        input_rows=base.input_rows,
        analysis_rows=len(analysis),
        excluded_row_positions=excluded,
        provenance=provenance2,
        model_kind="moderated",
        extra_columns={"moderator": zcol},
        moderation={"path": moderation_path, "component": moderator_component},
        moderated=True,
        estimable_paths=estimable_paths,
    )

fit_multilevel_moderated_gaze_mediation

fit_multilevel_moderated_gaze_mediation(prepared: Any, *, moderation_path: str = 'a', moderator_component: str = 'within', mediator_family: str = 'gaussian', outcome_family: str = 'bernoulli', priors: MediationPriorSpecification | None = None, random_slopes: Sequence[str] = (), missingness_policy: str = 'error', chains: int = 4, draws: int = 1000, tune: int = 1000, cores: int | None = None, seed: int = 1, target_accept: float = 0.95, max_treedepth: int = 12) -> MultilevelMediationFit

Fit moderation of the within a or b path using an explicitly prepared moderator.

Source code in src/gp3bayespy/multilevel_mediation.py

def fit_multilevel_moderated_gaze_mediation(
    prepared: Any,
    *,
    moderation_path: str = "a",
    moderator_component: str = "within",
    mediator_family: str = "gaussian",
    outcome_family: str = "bernoulli",
    priors: MediationPriorSpecification | None = None,
    random_slopes: Sequence[str] = (),
    missingness_policy: str = "error",
    chains: int = 4,
    draws: int = 1000,
    tune: int = 1000,
    cores: int | None = None,
    seed: int = 1,
    target_accept: float = 0.95,
    max_treedepth: int = 12,
) -> MultilevelMediationFit:
    """Fit moderation of the within a or b path using an explicitly prepared moderator."""
    spec = specify_multilevel_moderated_gaze_mediation(
        prepared,
        moderation_path=moderation_path,
        moderator_component=moderator_component,
        mediator_family=mediator_family,
        outcome_family=outcome_family,
        priors=priors,
        random_slopes=random_slopes,
        missingness_policy=missingness_policy,
    )
    controls = _validate_sampling_controls(
        chains=chains,
        draws=draws,
        tune=tune,
        cores=cores,
        seed=seed,
        target_accept=target_accept,
        max_treedepth=max_treedepth,
    )
    pm = _load_pymc()
    model = _build_moderated_pymc_model(spec)
    with model:
        idata = pm.sample(
            draws=controls["draws"],
            tune=controls["tune"],
            chains=controls["chains"],
            cores=controls["cores"],
            random_seed=controls["seed"],
            target_accept=controls["target_accept"],
            nuts={"max_treedepth": controls["max_treedepth"]},
            progressbar=False,
            compute_convergence_checks=False,
            return_inferencedata=True,
            idata_kwargs={"log_likelihood": True},
        )
    names = [
        "a_within",
        "a_between",
        "b_within",
        "b_between",
        "cprime_within",
        "cprime_between",
        "moderation",
        "indirect_within",
        "indirect_between",
    ]
    fit = MultilevelMediationFit(
        specification=spec,
        posterior=_flatten_posterior(idata, names),
        sampler_diagnostics=_sampler_diagnostics(
            idata, max_treedepth=int(controls["max_treedepth"])
        ),
        backend_fit=idata,
        backend_model=model,
        package_versions={"pymc": _package_version("pymc"), "arviz": _package_version("arviz")},
    )
    if not check_mediation_convergence(fit).passed:
        warnings.warn(
            "Moderated mediation fit failed critical convergence checks.",
            RuntimeWarning,
            stacklevel=2,
        )
    return fit

posterior_conditional_indirect_effect

posterior_conditional_indirect_effect(fit: MultilevelMediationFit, *, moderator_value: float, probability: float = 0.95, require_convergence: bool = True) -> MediationEffect

Compute the within-person indirect effect at a specified moderator value.

Source code in src/gp3bayespy/multilevel_mediation.py

def posterior_conditional_indirect_effect(
    fit: MultilevelMediationFit,
    *,
    moderator_value: float,
    probability: float = 0.95,
    require_convergence: bool = True,
) -> MediationEffect:
    """Compute the within-person indirect effect at a specified moderator value."""
    if fit.specification.model_kind != "moderated":
        raise GP3BayesError("`fit` is not a moderated mediation model.")
    _require_acceptable_fit(fit, require_convergence)
    if (
        isinstance(moderator_value, bool)
        or not isinstance(moderator_value, (int, float))
        or not math.isfinite(float(moderator_value))
    ):
        raise GP3BayesError("`moderator_value` must be one finite numeric value.")
    z = float(moderator_value)
    raw_a = fit.posterior.get("a_within")
    raw_b = fit.posterior.get("b_within")
    raw_mod = fit.posterior.get("moderation")
    if raw_a is None or raw_b is None or raw_mod is None:
        raise GP3BayesError("Required moderated path posterior draws are unavailable.")
    a = np.asarray(raw_a)
    b = np.asarray(raw_b)
    mod = np.asarray(raw_mod)
    draws = (a + mod * z) * b if fit.specification.moderation["path"] == "a" else a * (b + mod * z)
    return _effect_from_draws(
        f"conditional_indirect_within_z_{z:g}", draws, probability, "linear_predictor_product"
    )