gp3bayespy.duration

12 public functions in this module.

← API reference hub

Refit the approved duration model under declared prior-scale multipliers.

Source code in src/gp3bayespy/duration.py

def assess_duration_prior_sensitivity(
    fit: DurationFit,
    scale_multipliers: Mapping[str, float] | None = None,
    chains: int | None = None,
    iter: int | None = None,
    warmup: int | None = None,
    cores: int | None = None,
    seed: int | None = None,
    adapt_delta: float | None = None,
    max_treedepth: int | None = None,
    refresh: int = 0,
    maximum_standardized_shift: float = 0.25,
    review_standardized_shift: float = 0.50,
    retain_fits: bool = False,
):
    """Refit the approved duration model under declared prior-scale multipliers."""
    if not isinstance(fit, DurationFit):
        raise GP3BayesError("`fit` must be a duration gp3bayes fit.")
    from .sensitivity import PriorSensitivity, _classify_upper, _worst_status

    multipliers = (
        {"tighter": 0.5, "wider": 2.0} if scale_multipliers is None else dict(scale_multipliers)
    )
    if not multipliers or any(
        not str(k) or not math.isfinite(float(v)) or float(v) <= 0 for k, v in multipliers.items()
    ):
        raise GP3BayesError("`scale_multipliers` must be named unique positive finite values.")
    if maximum_standardized_shift < 0 or review_standardized_shift < maximum_standardized_shift:
        raise GP3BayesError("Sensitivity shift thresholds are invalid.")
    controls = fit.sampling
    chains = int(controls["chains"] if chains is None else chains)
    iter = int(controls["iter"] if iter is None else iter)
    warmup = int(controls["warmup"] if warmup is None else warmup)
    cores = int(controls["cores"] if cores is None else cores)
    seed = int(controls["seed"] + 2000 if seed is None else seed)
    adapt_delta = float(controls["adapt_delta"] if adapt_delta is None else adapt_delta)
    max_treedepth = int(controls["max_treedepth"] if max_treedepth is None else max_treedepth)
    reference_summary = summarise_duration_posterior(fit).table
    reference = reference_summary.loc[
        reference_summary["variable"].astype(str).str.match(r"^b_|^sigma$"),
        ["variable", "median", "sd"],
    ].copy()
    reference_diagnostics = diagnose_duration_fit(fit)
    base = fit.specification.priors
    rows = []
    alternative_fits = {}
    for index, (label, mvalue) in enumerate(multipliers.items()):
        multiplier = float(mvalue)
        cor_rows = base.table.loc[base.table["parameter_class"].eq("cor")]
        priors = create_prior_specification(
            fit.specification.contract,
            baseline=base.baseline,
            intercept_scale=_duration_prior_scale_from_table(base, "Intercept") * multiplier,
            coefficient_scale=_duration_prior_scale_from_table(base, "b") * multiplier,
            group_sd_scale=_duration_prior_scale_from_table(base, "sd") * multiplier,
            residual_scale=_duration_prior_scale_from_table(base, "sigma") * multiplier,
            correlation_eta=float(cor_rows["shape"].iloc[0]) if not cor_rows.empty else 2.0,
            student_df=float(base.table.loc[base.table["parameter_class"].eq("sd"), "df"].iloc[0]),
        )
        specification = replace(fit.specification, priors=priors)
        alternative_fit = fit_duration_model(
            specification,
            chains=chains,
            iter=iter,
            warmup=warmup,
            cores=cores,
            seed=seed + index,
            adapt_delta=adapt_delta,
            max_treedepth=max_treedepth,
            refresh=refresh,
        )
        alternative_diagnostics = diagnose_duration_fit(alternative_fit)
        alternative_summary = summarise_duration_posterior(alternative_fit).table
        alternative = alternative_summary.loc[
            alternative_summary["variable"].astype(str).str.match(r"^b_|^sigma$"),
            ["variable", "median"],
        ].rename(columns={"median": "alternative_median"})
        merged = reference.merge(alternative, on="variable", how="inner", sort=False)
        merged["scenario"] = str(label)
        merged["scale_multiplier"] = multiplier
        merged["median_shift"] = merged["alternative_median"] - merged["median"]
        merged["standardized_shift"] = merged["median_shift"].abs() / np.maximum(
            merged["sd"].abs(), np.finfo(float).eps
        )
        merged["shift_status"] = [
            _classify_upper(float(v), maximum_standardized_shift, review_standardized_shift)
            for v in merged["standardized_shift"]
        ]
        merged["diagnostic_status"] = alternative_diagnostics.status
        merged["status"] = [
            _worst_status([v, alternative_diagnostics.status]) for v in merged["shift_status"]
        ]
        rows.append(
            merged[
                [
                    "scenario",
                    "scale_multiplier",
                    "variable",
                    "median",
                    "alternative_median",
                    "median_shift",
                    "standardized_shift",
                    "shift_status",
                    "diagnostic_status",
                    "status",
                ]
            ]
        )
        if retain_fits:
            alternative_fits[str(label)] = alternative_fit
    comparison = pd.concat(rows, ignore_index=True) if rows else pd.DataFrame()
    scenario_rows = []
    for scenario, frame in comparison.groupby("scenario", sort=False):
        scenario_rows.append(
            {
                "scenario": scenario,
                "scale_multiplier": float(frame["scale_multiplier"].iloc[0]),
                "maximum_standardized_shift": float(frame["standardized_shift"].max()),
                "diagnostic_status": str(frame["diagnostic_status"].iloc[0]),
                "status": _worst_status(frame["status"].astype(str).tolist()),
            }
        )
    scenario_status = pd.DataFrame(scenario_rows)
    overall = _worst_status(
        [
            reference_diagnostics.status,
            *scenario_status.get("status", pd.Series(dtype=str)).astype(str).tolist(),
        ]
    )
    return PriorSensitivity(
        "0.1",
        "duration",
        multipliers,
        comparison,
        scenario_status,
        reference_diagnostics.status,
        overall,
        alternative_fits if retain_fits else None,
    )

Check selected duration posterior-predictive summaries conservatively.

Source code in src/gp3bayespy/duration.py

def check_duration_posterior_predictive(
    fit: DurationFit,
    draws: int = 500,
    seed: int = 1,
    pass_probability: float = 0.80,
    review_probability: float = 0.95,
):
    """Check selected duration posterior-predictive summaries conservatively."""
    if not isinstance(fit, DurationFit):
        raise GP3BayesError("`fit` must inherit from `gp3bayes_fit`.")

    from .ppc import (
        _check_table,
        _duration_result,
        _duration_summary,
        _replicated_table,
        _validate_controls,
    )
    from .predictive import extract_expected_predictions, extract_posterior_predictions

    draw_count, seed_value, pass_value, review_value = _validate_controls(
        draws, seed, pass_probability, review_probability
    )
    prepared = cast(DurationPrepared, fit.specification.prepared)
    data = prepared.data
    contract = fit.specification.contract
    outcome_col = cast(str, contract.mappings["outcome"])
    participant_col = cast(str, contract.mappings["participant"])
    condition_col = contract.mappings["condition"]
    item_col = contract.mappings["item"]

    y = pd.to_numeric(data[outcome_col], errors="raise").to_numpy(dtype=float)
    participant = data[participant_col].to_numpy(copy=True)
    condition = None if condition_col is None else data[condition_col].to_numpy(copy=True)
    item = None if item_col is None else data[item_col].to_numpy(copy=True)

    yrep = extract_posterior_predictions(
        fit,
        newdata=data,
        include_group_effects=True,
        allow_new_levels=False,
        ndraws=draw_count,
        seed=seed_value,
    )
    if not np.isfinite(yrep).all() or np.any(yrep <= 0):
        raise GP3BayesError("Posterior predictive draws must be finite and strictly positive.")
    replicated = _replicated_table(
        yrep,
        summary_function=_duration_summary,
        condition=condition,
        participant=participant,
        item=item,
    )
    observed = _duration_summary(
        y,
        condition=condition,
        participant=participant,
        item=item,
    )
    checks = _check_table(
        observed,
        replicated,
        pass_probability=pass_value,
        review_probability=review_value,
    )
    expected = extract_expected_predictions(
        fit,
        newdata=data,
        include_group_effects=True,
        allow_new_levels=False,
        ndraws=draw_count,
    )
    predicted_mean = np.mean(expected, axis=0)
    log_scale_rmse = float(np.sqrt(np.mean((np.log(y) - np.log(predicted_mean)) ** 2)))
    return _duration_result(
        outcome_unit=fit.outcome_unit,
        draws=int(yrep.shape[0]),
        seed=seed_value,
        observed=observed,
        replicated=replicated,
        checks=checks,
        log_scale_rmse=log_scale_rmse,
    )

Simulate lognormal prior predictive outcomes without fitting a model.

Source code in src/gp3bayespy/duration.py

def check_duration_prior_predictive(
    specification: DurationModelSpecification,
    draws: int = 500,
    seed: int = 1,
    plausible_median: Sequence[float] | None = None,
    maximum_q99: float | None = None,
    maximum_cv: float = 5.0,
    maximum_condition_ratio: float = 10.0,
    maximum_extreme_probability: float = 0.25,
) -> DurationPriorPredictiveCheck:
    """Simulate lognormal prior predictive outcomes without fitting a model."""
    if not isinstance(specification, DurationModelSpecification):
        raise GP3BayesError(
            "`specification` must inherit from `gp3bayes_duration_model_specification`."
        )
    prepared = specification.prepared
    if not isinstance(prepared, DurationPrepared):
        raise GP3BayesError(
            "`specification$prepared` must inherit from `gp3bayes_duration_prepared`."
        )
    validate_prior_specification(specification.priors, specification.contract)
    draws = _integer(draws, "draws", minimum=50)
    seed = _integer(seed, "seed", minimum=0)
    baseline = float(specification.priors.baseline)
    plausible = _positive_pair(plausible_median, baseline)
    q99_limit = _numeric_scalar(
        baseline * 50 if maximum_q99 is None else maximum_q99,
        "maximum_q99",
        lower=0,
        lower_open=True,
    )
    maximum_cv = _numeric_scalar(maximum_cv, "maximum_cv", lower=0, lower_open=True)
    maximum_condition_ratio = _numeric_scalar(
        maximum_condition_ratio,
        "maximum_condition_ratio",
        lower=1,
        lower_open=True,
    )
    maximum_extreme_probability = _numeric_scalar(
        maximum_extreme_probability,
        "maximum_extreme_probability",
        lower=0,
        upper=1,
    )

    data = prepared.data
    contract = prepared.contract
    model_matrix, _ = _fixed_model_matrix(data, contract)
    participant_column = cast(str, contract.mappings["participant"])
    participant = data[participant_column].astype(str).to_numpy()
    _, participant_index = np.unique(participant, return_inverse=True)
    participant_count = int(participant_index.max()) + 1

    item: np.ndarray | None = None
    item_index: np.ndarray | None = None
    item_count = 0
    item_column = contract.mappings["item"]
    if item_column is not None:
        item = data[item_column].astype(str).to_numpy()
        _, item_index = np.unique(item, return_inverse=True)
        item_count = int(item_index.max()) + 1

    condition: np.ndarray | None = None
    condition_column = contract.mappings["condition"]
    if condition_column is not None:
        condition = pd.to_numeric(data[condition_column], errors="raise").to_numpy(dtype=float)

    intercept_prior = _prior_row(specification.priors, "Intercept")
    coefficient_prior = _prior_row(specification.priors, "b")
    group_prior = _prior_row(specification.priors, "sd")
    sigma_prior = _prior_row(specification.priors, "sigma")
    correlation_prior = _prior_row(specification.priors, "cor") if contract.random_slope else None

    rng = np.random.RandomState(seed)
    rows: list[dict[str, float]] = []
    for _ in range(draws):
        intercept = rng.normal(
            loc=float(intercept_prior["location"]),
            scale=float(intercept_prior["scale"]),
        )
        coefficient_count = model_matrix.shape[1] - 1
        linear_predictor = np.full(len(data), intercept, dtype=float)
        if coefficient_count > 0:
            coefficients = rng.normal(
                loc=float(coefficient_prior["location"]),
                scale=float(coefficient_prior["scale"]),
                size=coefficient_count,
            )
            linear_predictor = linear_predictor + model_matrix[:, 1:] @ coefficients

        group_scale = float(group_prior["scale"])
        group_df = float(group_prior["df"])
        participant_intercept_sd = abs(rng.standard_t(group_df)) * group_scale
        participant_slope_sd = (
            abs(rng.standard_t(group_df)) * group_scale if contract.random_slope else 0.0
        )
        z_intercept = rng.normal(size=participant_count)
        z_slope = rng.normal(size=participant_count)
        correlation = 0.0
        if contract.random_slope:
            assert correlation_prior is not None
            shape = float(correlation_prior["shape"])
            correlation = 2 * rng.beta(shape, shape) - 1
        participant_intercept = participant_intercept_sd * z_intercept
        participant_slope = participant_slope_sd * (
            correlation * z_intercept + math.sqrt(1 - correlation**2) * z_slope
        )
        linear_predictor = linear_predictor + participant_intercept[participant_index]
        if contract.random_slope:
            assert condition is not None
            linear_predictor = linear_predictor + participant_slope[participant_index] * condition

        if item is not None:
            assert item_index is not None
            item_effect_sd = abs(rng.standard_t(group_df)) * group_scale
            item_effect = rng.normal(loc=0.0, scale=item_effect_sd, size=item_count)
            linear_predictor = linear_predictor + item_effect[item_index]

        sigma = abs(rng.standard_t(float(sigma_prior["df"]))) * float(sigma_prior["scale"])
        y = rng.lognormal(mean=linear_predictor, sigma=sigma, size=len(data))
        rows.append(_duration_summary(y, condition, participant, item))

    summaries = pd.DataFrame(rows)
    median_violation = float(
        np.mean((summaries["median"] < plausible[0]) | (summaries["median"] > plausible[1]))
    )
    q99_violation = float(np.mean(summaries["q99"] > q99_limit))
    cv_values = summaries["coefficient_of_variation"].to_numpy(dtype=float)
    cv_violation = (
        1.0
        if np.isnan(cv_values).all()
        else float(np.mean(cv_values[~np.isnan(cv_values)] > maximum_cv))
    )
    ratios = summaries["condition_median_ratio"].to_numpy(dtype=float)
    if np.isnan(ratios).all():
        condition_violation = math.nan
    else:
        finite = ratios[~np.isnan(ratios)]
        condition_violation = float(
            np.mean((finite > maximum_condition_ratio) | (finite < 1 / maximum_condition_ratio))
        )
    nonfinite_violation = float(np.mean(summaries["nonfinite_fraction"] > 0))

    probabilities = [
        median_violation,
        q99_violation,
        cv_violation,
        condition_violation,
        nonfinite_violation,
    ]
    statuses = [
        "not_applicable"
        if math.isnan(probability)
        else ("pass" if probability <= maximum_extreme_probability else "fail")
        for probability in probabilities
    ]
    checks = pd.DataFrame(
        {
            "check": [
                "overall_median",
                "upper_tail_q99",
                "coefficient_of_variation",
                "condition_median_ratio",
                "nonfinite_predictions",
            ],
            "violation_probability": probabilities,
            "maximum_probability": [maximum_extreme_probability] * 5,
            "status": statuses,
        }
    )
    adequate = bool(checks["status"].isin(["pass", "not_applicable"]).all())
    return DurationPriorPredictiveCheck(
        check_version="0.1",
        family="duration",
        draws=draws,
        seed=seed,
        outcome_unit=prepared.outcome_unit,
        summaries=summaries,
        checks=checks,
        thresholds={
            "plausible_median": plausible,
            "maximum_q99": q99_limit,
            "maximum_cv": maximum_cv,
            "maximum_condition_ratio": maximum_condition_ratio,
            "maximum_extreme_probability": maximum_extreme_probability,
        },
        adequate=adequate,
    )

Write an explicit Markdown inventory for a fitted duration model.

Source code in src/gp3bayespy/duration.py

def create_duration_model_report(
    fit: DurationFit,
    diagnostics: Any = None,
    posterior_summary: Any = None,
    posterior_predictive: Any = None,
    prior_sensitivity: Any = None,
    recovery: Any = None,
    file: str = "",
    overwrite: bool = False,
):
    """Write an explicit Markdown inventory for a fitted duration model."""
    from pathlib import Path

    from .sensitivity import ModelReport

    if not isinstance(fit, DurationFit):
        raise GP3BayesError("`fit` must be a duration gp3bayes fit.")
    path = Path(file)
    if path.suffix.lower() != ".md" or not file:
        raise GP3BayesError("`file` must end in `.md`.")
    if path.exists() and not overwrite:
        raise GP3BayesError("The report file already exists. Use `overwrite=True`.")
    if not path.parent.exists():
        raise GP3BayesError("The report parent directory does not exist.")
    diagnostics = diagnose_duration_fit(fit) if diagnostics is None else diagnostics
    posterior_summary = (
        summarise_duration_posterior(fit) if posterior_summary is None else posterior_summary
    )
    lines = [
        "# gp3bayes duration model report",
        "",
        f"- Formula: `{fit.translation.formula_text}`",
        "- Family: lognormal",
        f"- Outcome unit: {fit.outcome_unit}",
        f"- Interface: {fit.backend_interface}",
        f"- Sampling backend: {fit.sampling_backend}",
        "",
        "## Sampling diagnostics",
        "",
        f"**Threshold status: {diagnostics.status}**",
        "",
        _duration_markdown_table(diagnostics.component_table),
        "",
        "Numerical thresholds were assessed, but no automatic convergence or posterior-adequacy claim is made.",
        "",
        "## Posterior summaries",
        "",
        _duration_markdown_table(posterior_summary.table),
    ]
    registry = [
        {"section": "sampling_diagnostics", "status": diagnostics.status},
        {"section": "posterior_summary", "status": "reported"},
    ]
    for name, obj in (
        ("posterior_predictive", posterior_predictive),
        ("prior_sensitivity", prior_sensitivity),
        ("recovery", recovery),
    ):
        if obj is not None:
            status = getattr(obj, "status", "reported")
            lines += ["", f"## {name.replace('_', ' ').title()}", "", f"Status: {status}"]
            registry.append({"section": name, "status": status})
    lines += [
        "",
        "## Interpretation boundary",
        "",
        "This report does not automatically establish convergence, posterior adequacy, robustness, substantive validity, or causal identification.",
    ]
    path.write_text("\n".join(lines) + "\n", encoding="utf-8")
    return ModelReport("duration", str(path.resolve()), pd.DataFrame(registry))

Apply the frozen sampling-diagnostic thresholds to a duration fit.

Source code in src/gp3bayespy/duration.py

def diagnose_duration_fit(
    fit: DurationFit,
    rhat_pass: float = 1.01,
    rhat_fail: float = 1.05,
    ess_per_chain_pass: float = 100,
    ess_per_chain_fail: float = 50,
    maximum_treedepth_fraction: float = 0.01,
    ebfmi_pass: float = 0.30,
    ebfmi_fail: float = 0.20,
):
    """Apply the frozen sampling-diagnostic thresholds to a duration fit."""
    if not isinstance(fit, DurationFit):
        raise GP3BayesError("`fit` must inherit from `gp3bayes_fit`.")
    return _diagnose_fit(
        fit,
        family="duration",
        rhat_pass=rhat_pass,
        rhat_fail=rhat_fail,
        ess_per_chain_pass=ess_per_chain_pass,
        ess_per_chain_fail=ess_per_chain_fail,
        maximum_treedepth_fraction=maximum_treedepth_fraction,
        ebfmi_pass=ebfmi_pass,
        ebfmi_fail=ebfmi_fail,
    )

Fit the approved lognormal duration model with optional PyMC NUTS.

Source code in src/gp3bayespy/duration.py

def fit_duration_model(
    specification: DurationModelSpecification,
    chains: int = 4,
    iter: int = 2000,
    warmup: int = 1000,
    cores: int | None = None,
    seed: int = 1,
    adapt_delta: float = 0.95,
    max_treedepth: int = 12,
    refresh: int = 0,
) -> DurationFit:
    """Fit the approved lognormal duration model with optional PyMC NUTS."""
    specification = _validate_duration_model_specification(specification)
    controls = _validate_sampling_controls(
        chains=chains,
        iter=iter,
        warmup=warmup,
        cores=cores,
        seed=seed,
        adapt_delta=adapt_delta,
        max_treedepth=max_treedepth,
        refresh=refresh,
    )
    translation = translate_duration_model_to_brms(specification)
    _require_pymc("fit a duration model through the approved Python sampling backend")
    backend_model, backend_fit = _run_duration_pymc(specification, controls.as_dict())
    return DurationFit(
        fit_version="0.1",
        family="duration",
        model_family="hierarchical_lognormal_duration",
        specification=specification,
        translation=translation,
        backend_fit=backend_fit,
        backend_model=backend_model,
        outcome_unit=specification.outcome_unit,
        backend_interface="pymc",
        sampling_backend="pymc",
        algorithm="NUTS",
        sampling=controls.as_dict(),
        package_versions=_backend_versions(),
    )

Validate, explicitly convert, scale, and readiness-gate duration data.

Source code in src/gp3bayespy/duration.py

def prepare_hierarchical_duration_data(
    data: pd.DataFrame,
    contract: ModelContract,
    condition_levels: Sequence[object] | None = None,
    condition_coding: Sequence[float] = (-0.5, 0.5),
    scale_predictors: Sequence[str] | str = (),
    scale_time: bool = False,
    outcome_multiplier: float = 1.0,
    converted_unit: str | None = None,
    missing: str = "error",
) -> DurationPrepared:
    """Validate, explicitly convert, scale, and readiness-gate duration data."""
    if not isinstance(data, pd.DataFrame):
        raise GP3BayesError("`data` must be a data frame.")
    _validate_duration_contract(contract)
    if missing not in {"error", "drop"}:
        raise GP3BayesError('`missing` must be either "error" or "drop".')
    scale_predictors_tuple = _character_vector(scale_predictors, "scale_predictors")
    scale_time = _flag(scale_time, "scale_time")
    outcome_multiplier = _numeric_scalar(
        outcome_multiplier, "outcome_multiplier", lower=0, lower_open=True
    )
    if any(value not in contract.predictors for value in scale_predictors_tuple):
        raise GP3BayesError("Every scaled predictor must be declared in the model contract.")

    if not math.isclose(outcome_multiplier, 1.0):
        if not isinstance(converted_unit, str) or not converted_unit:
            raise GP3BayesError(
                "`converted_unit` must be supplied when `outcome_multiplier` is not one."
            )
    elif converted_unit is not None and (not isinstance(converted_unit, str) or not converted_unit):
        raise GP3BayesError("`converted_unit` must be NULL or one non-empty character value.")

    required = _required_columns(contract)
    absent = [value for value in required if value not in data.columns]
    if absent:
        raise GP3BayesError("Required duration columns were not found: " + ", ".join(absent) + ".")

    required_frame = cast(pd.DataFrame, data.loc[:, list(required)])
    complete_rows = ~required_frame.isna().any(axis=1)
    dropped_positions = [index + 1 for index, complete in enumerate(complete_rows) if not complete]
    if dropped_positions and missing == "error":
        raise GP3BayesError(
            "Missing values were found in required duration columns. "
            'Use `missing = "drop"` only after an explicit exclusion decision.'
        )
    working = (
        cast(pd.DataFrame, data.loc[complete_rows]).copy() if dropped_positions else data.copy()
    )

    outcome_column = cast(str, contract.mappings["outcome"])
    raw_outcome = pd.to_numeric(working[outcome_column], errors="coerce")
    values = raw_outcome.to_numpy(dtype=float)
    if not np.isfinite(values).all():
        raise GP3BayesError("The duration outcome must contain only finite numeric values.")
    values = values * outcome_multiplier
    if np.any(values <= 0):
        raise GP3BayesError(
            "The duration outcome must be strictly positive. "
            "Zero and negative values are unsupported."
        )
    working[outcome_column] = values

    analysis_contract = (
        replace(contract, outcome_unit=converted_unit) if converted_unit is not None else contract
    )
    analysis_unit = cast(str, analysis_contract.outcome_unit)
    transformations: dict[str, Any] = {
        "outcome": {
            "source_unit": contract.outcome_unit,
            "analysis_unit": analysis_unit,
            "multiplier": outcome_multiplier,
            "strictly_positive": True,
            "finite": True,
            "censored": False,
        },
        "condition": None,
        "scaled_columns": {},
        "missing": {
            "action": missing,
            "dropped_row_positions": tuple(dropped_positions),
        },
    }

    condition_column = analysis_contract.mappings["condition"]
    if condition_column is not None:
        coded, source_levels, coding = _code_condition(
            working[condition_column], condition_levels, condition_coding
        )
        working[condition_column] = coded
        transformations["condition"] = {
            "source_levels": source_levels,
            "coding": coding,
        }

    scaling_columns = list(scale_predictors_tuple)
    time_column = analysis_contract.mappings["time"]
    if scale_time and time_column is not None and time_column not in scaling_columns:
        scaling_columns.append(time_column)
    scaled_registry = cast(dict[str, dict[str, float | str]], transformations["scaled_columns"])
    for column in scaling_columns:
        series = working[column]
        if not pd.api.types.is_numeric_dtype(series.dtype):
            raise GP3BayesError(
                f"Only numeric predictors can be scaled. `{column}` is not numeric."
            )
        numeric = pd.to_numeric(series, errors="coerce").to_numpy(dtype=float)
        center = float(np.mean(numeric))
        scale_value = _sample_sd(numeric)
        if not math.isfinite(scale_value) or scale_value <= 0:
            raise GP3BayesError(
                f"The declared scaling column `{column}` has zero or undefined standard deviation."
            )
        working[column] = (numeric - center) / scale_value
        scaled_registry[column] = {
            "action": "centre_and_scale",
            "centre": center,
            "scale": scale_value,
        }

    audit = audit_model_readiness(working, analysis_contract)
    if not audit.ready:
        raise GP3BayesError("The prepared duration data did not pass the readiness gate.")

    fixed_formula = _fixed_formula_text(analysis_contract)
    model_matrix, matrix_columns = _fixed_model_matrix(working, analysis_contract)
    if np.linalg.matrix_rank(model_matrix) < model_matrix.shape[1]:
        raise GP3BayesError("The prepared fixed-effects design matrix is rank deficient.")

    condition_transform = transformations["condition"]
    if condition_transform is None:
        condition_value = "not applicable"
    else:
        condition_map = cast(Mapping[str, float], condition_transform["coding"])
        condition_value = ", ".join(f"{key}={value:g}" for key, value in condition_map.items())
    decision_log = pd.DataFrame(
        {
            "decision": [
                "outcome_validation",
                "unit_conversion",
                "missing_values",
                "condition_coding",
                "predictor_scaling",
            ],
            "value": [
                "strictly positive finite uncensored",
                f"{contract.outcome_unit} x {outcome_multiplier:g} -> {analysis_unit}",
                f"{missing}; rows removed = {len(dropped_positions)}",
                condition_value,
                ", ".join(scaling_columns) if scaling_columns else "none",
            ],
        }
    )

    return DurationPrepared(
        preparation_version="0.1",
        family="duration",
        data=working,
        source_contract=contract,
        contract=analysis_contract,
        audit=audit,
        transformations=transformations,
        decision_log=decision_log,
        fixed_formula=fixed_formula,
        fixed_formula_text=fixed_formula,
        model_matrix_columns=matrix_columns,
        n_input_rows=len(data),
        n_analysis_rows=len(working),
        rows_removed=len(dropped_positions),
        outcome_unit=analysis_unit,
    )

Run simulation-based duration parameter recovery with the restricted fitter.

Source code in src/gp3bayespy/duration.py

def run_duration_recovery(
    repetitions: int = 20,
    n_participants: int = 30,
    trials_per_participant: int = 16,
    n_items: int = 12,
    include_items: bool = True,
    random_slope: bool = True,
    baseline_median: float = 500.0,
    outcome_unit: str = "milliseconds",
    seed: int = 2001,
    chains: int = 4,
    iter: int = 1500,
    warmup: int = 750,
    cores: int | None = None,
    adapt_delta: float = 0.95,
    max_treedepth: int = 12,
    refresh: int = 0,
    interval_probability: float = 0.95,
    minimum_repetitions: int = 20,
    maximum_standardized_bias: float = 0.25,
    minimum_coverage: float = 0.80,
    minimum_diagnostic_pass_fraction: float = 0.80,
    continue_on_error: bool = True,
):
    """Run simulation-based duration parameter recovery with the restricted fitter."""
    from .contracts import create_model_contract
    from .sensitivity import RecoveryResult, _worst_status

    repetitions = _integer(repetitions, "repetitions", minimum=2)
    estimates = []
    fit_rows = []
    for rep in range(repetitions):
        rep_seed = int(seed) + rep
        try:
            sim = simulate_hierarchical_duration_data(
                n_participants=n_participants,
                trials_per_participant=trials_per_participant,
                n_items=n_items,
                include_items=include_items,
                random_slope_sd=0.15 if random_slope else 0.0,
                baseline_median=baseline_median,
                outcome_unit=outcome_unit,
                seed=rep_seed,
            )
            contract = create_model_contract(
                family="duration",
                outcome_col="duration",
                participant_col="participant_id",
                item_col="item_id" if include_items else None,
                trial_col="trial_id",
                condition_col="condition",
                predictors=("participant_covariate", "trial_covariate"),
                interaction=("condition", "participant_covariate"),
                random_slope=random_slope,
                outcome_unit=outcome_unit,
            )
            prepared = prepare_hierarchical_duration_data(
                sim.data, contract, converted_unit=outcome_unit
            )
            specification = specify_duration_model(prepared, baseline=baseline_median)
            fitted = fit_duration_model(
                specification,
                chains=chains,
                iter=iter,
                warmup=warmup,
                cores=cores,
                seed=rep_seed,
                adapt_delta=adapt_delta,
                max_treedepth=max_treedepth,
                refresh=refresh,
            )
            diagnostics = diagnose_duration_fit(fitted)
            table = summarise_duration_posterior(fitted, probability=interval_probability).table
            truth = {
                "b_Intercept": sim.truth["fixed_effects"]["(Intercept)"],
                "sigma": sim.truth["residual_sd"],
            }
            truth.update(
                {f"b_{k}": v for k, v in sim.truth["fixed_effects"].items() if k != "(Intercept)"}
            )
            for variable, value in truth.items():
                row = table.loc[table["variable"].eq(variable)]
                if row.empty:
                    continue
                r = row.iloc[0]
                estimates.append(
                    {
                        "repetition": rep + 1,
                        "variable": variable,
                        "truth": float(value),
                        "median": float(r["median"]),
                        "lower": float(r["lower"]),
                        "upper": float(r["upper"]),
                        "covered": bool(float(r["lower"]) <= float(value) <= float(r["upper"])),
                    }
                )
            fit_rows.append(
                {
                    "repetition": rep + 1,
                    "completed": True,
                    "diagnostic_status": diagnostics.status,
                    "message": "",
                }
            )
        except Exception as exc:
            if not continue_on_error:
                raise
            fit_rows.append(
                {
                    "repetition": rep + 1,
                    "completed": False,
                    "diagnostic_status": "error",
                    "message": str(exc),
                }
            )
    estimates_df = pd.DataFrame(estimates)
    fit_df = pd.DataFrame(fit_rows)
    summaries = []
    if not estimates_df.empty:
        for variable, frame in estimates_df.groupby("variable", sort=False):  # type: ignore[assignment]
            error = frame["median"].to_numpy(float) - frame["truth"].to_numpy(float)
            sd = (
                float(np.std(frame["median"].to_numpy(float), ddof=1)) if len(frame) > 1 else np.nan
            )
            bias = float(np.mean(error))
            standardized = abs(bias) / max(sd, np.finfo(float).eps) if np.isfinite(sd) else np.nan
            coverage = float(frame["covered"].mean())
            rmse = float(np.sqrt(np.mean(error**2)))
            summaries.append(
                {
                    "variable": variable,
                    "bias": bias,
                    "standardized_bias": standardized,
                    "coverage": coverage,
                    "rmse": rmse,
                    "n": len(frame),
                    "status": "pass"
                    if len(frame) >= minimum_repetitions
                    and np.isfinite(standardized)
                    and standardized <= maximum_standardized_bias
                    and coverage >= minimum_coverage
                    else "review",
                }
            )
    parameter_summary = pd.DataFrame(summaries)
    diagnostic_pass = (
        float(np.mean(fit_df["diagnostic_status"].eq("pass"))) if not fit_df.empty else 0.0
    )
    status = _worst_status(
        parameter_summary.get("status", pd.Series(["review"])).astype(str).tolist()
    )
    status = "review" if diagnostic_pass < minimum_diagnostic_pass_fraction else status
    return RecoveryResult(
        "0.1", "duration", parameter_summary, estimates_df, fit_df, repetitions, status
    )

Simulate governed positive uncensored lognormal duration data.

Source code in src/gp3bayespy/duration.py

def simulate_hierarchical_duration_data(
    n_participants: int = 40,
    trials_per_participant: int = 20,
    n_items: int = 20,
    baseline_median: float = 500.0,
    condition_effect: float = math.log(1.15),
    participant_covariate_effect: float = math.log(1.08),
    trial_covariate_effect: float = math.log(1.04),
    interaction_effect: float = math.log(1.05),
    participant_sd: float = 0.35,
    item_sd: float = 0.20,
    random_slope_sd: float = 0.15,
    random_slope_cor: float = 0.0,
    residual_sd: float = 0.40,
    condition_probability: float = 0.5,
    balanced_condition: bool = True,
    include_items: bool = True,
    outcome_unit: str = "milliseconds",
    seed: int = 1,
) -> DurationSimulation:
    """Simulate governed positive uncensored lognormal duration data."""
    n_participants = _integer(n_participants, "n_participants", minimum=2)
    trials_per_participant = _integer(trials_per_participant, "trials_per_participant", minimum=2)
    n_items = _integer(n_items, "n_items", minimum=2)
    baseline_median = _numeric_scalar(baseline_median, "baseline_median", lower=0, lower_open=True)
    condition_effect = _numeric_scalar(condition_effect, "condition_effect")
    participant_covariate_effect = _numeric_scalar(
        participant_covariate_effect, "participant_covariate_effect"
    )
    trial_covariate_effect = _numeric_scalar(trial_covariate_effect, "trial_covariate_effect")
    interaction_effect = _numeric_scalar(interaction_effect, "interaction_effect")
    participant_sd = _numeric_scalar(participant_sd, "participant_sd", lower=0)
    item_sd = _numeric_scalar(item_sd, "item_sd", lower=0)
    random_slope_sd = _numeric_scalar(random_slope_sd, "random_slope_sd", lower=0)
    random_slope_cor = _numeric_scalar(
        random_slope_cor,
        "random_slope_cor",
        lower=-1,
        upper=1,
        lower_open=True,
        upper_open=True,
    )
    residual_sd = _numeric_scalar(residual_sd, "residual_sd", lower=0, lower_open=True)
    condition_probability = _numeric_scalar(
        condition_probability,
        "condition_probability",
        lower=0,
        upper=1,
        lower_open=True,
        upper_open=True,
    )
    balanced_condition = _flag(balanced_condition, "balanced_condition")
    include_items = _flag(include_items, "include_items")
    if not isinstance(outcome_unit, str) or not outcome_unit:
        raise GP3BayesError("`outcome_unit` must be one non-empty character value.")
    seed = _integer(seed, "seed", minimum=0)

    rng = np.random.RandomState(seed)
    participant_levels = np.array(
        [f"p{index:03d}" for index in range(1, n_participants + 1)], dtype=object
    )
    participant_id = np.repeat(participant_levels, trials_per_participant)
    trial_id = np.tile(np.arange(1, trials_per_participant + 1), n_participants)
    n_rows = len(participant_id)

    if balanced_condition:
        base = np.resize(np.array([-0.5, 0.5], dtype=float), trials_per_participant)
        condition_code = np.concatenate([rng.permutation(base) for _ in range(n_participants)])
    else:
        condition_code = np.where(rng.uniform(size=n_rows) < condition_probability, 0.5, -0.5)

    participant_covariate_by_id = _standardize(rng.normal(size=n_participants))
    participant_index = pd.Categorical(participant_id, categories=participant_levels).codes
    participant_covariate = participant_covariate_by_id[participant_index]
    trial_covariate = _standardize(rng.normal(size=n_rows))

    z_intercept = rng.normal(size=n_participants)
    z_slope = rng.normal(size=n_participants)
    participant_intercept = participant_sd * z_intercept
    participant_slope = random_slope_sd * (
        random_slope_cor * z_intercept + math.sqrt(1 - random_slope_cor**2) * z_slope
    )

    item_id: np.ndarray | None = None
    item_effect = np.zeros(n_rows, dtype=float)
    item_effect_by_id = np.array([], dtype=float)
    item_levels = np.array([], dtype=object)
    if include_items:
        item_levels = np.array([f"i{index:03d}" for index in range(1, n_items + 1)], dtype=object)
        participant_offsets = np.repeat(np.arange(n_participants), trials_per_participant)
        item_index = (trial_id + participant_offsets - 1) % n_items
        item_id = item_levels[item_index]
        item_effect_by_id = rng.normal(loc=0.0, scale=item_sd, size=n_items)
        item_effect = item_effect_by_id[item_index]

    linear_predictor = (
        math.log(baseline_median)
        + condition_effect * condition_code
        + participant_covariate_effect * participant_covariate
        + trial_covariate_effect * trial_covariate
        + interaction_effect * condition_code * participant_covariate
        + participant_intercept[participant_index]
        + participant_slope[participant_index] * condition_code
        + item_effect
    )
    duration = rng.lognormal(mean=linear_predictor, sigma=residual_sd, size=n_rows)

    data_dict: dict[str, Any] = {
        "participant_id": participant_id,
        "trial_id": trial_id.astype(int),
        "condition": pd.Categorical(
            np.where(condition_code < 0, "control", "treatment"),
            categories=["control", "treatment"],
            ordered=False,
        ),
        "participant_covariate": participant_covariate,
        "trial_covariate": trial_covariate,
        "duration": duration,
        "true_median": np.exp(linear_predictor),
        "true_mean": np.exp(linear_predictor + residual_sd**2 / 2),
    }
    if include_items:
        data_dict["item_id"] = item_id
    preferred = [
        "participant_id",
        *(["item_id"] if include_items else []),
        "trial_id",
        "condition",
        "participant_covariate",
        "trial_covariate",
        "duration",
        "true_median",
        "true_mean",
    ]
    data = cast(pd.DataFrame, pd.DataFrame(data_dict).loc[:, preferred])

    truth: dict[str, Any] = {
        "fixed_effects": {
            "(Intercept)": math.log(baseline_median),
            "condition": condition_effect,
            "participant_covariate": participant_covariate_effect,
            "trial_covariate": trial_covariate_effect,
            "condition:participant_covariate": interaction_effect,
        },
        "baseline_median": baseline_median,
        "participant_sd": participant_sd,
        "item_sd": item_sd if include_items else 0.0,
        "random_slope_sd": random_slope_sd,
        "random_slope_cor": random_slope_cor,
        "residual_sd": residual_sd,
        "condition_coding": {"control": -0.5, "treatment": 0.5},
        "condition_probability": condition_probability,
        "balanced_condition": balanced_condition,
        "outcome_unit": outcome_unit,
        "seed": seed,
    }
    participant_re = pd.DataFrame(
        {
            "participant_id": participant_levels,
            "intercept": participant_intercept,
            "condition_slope": participant_slope,
            "participant_covariate": participant_covariate_by_id,
        }
    )
    item_re = (
        pd.DataFrame({"item_id": item_levels, "intercept": item_effect_by_id})
        if include_items
        else None
    )
    design = {
        "n_participants": n_participants,
        "trials_per_participant": trials_per_participant,
        "n_items": n_items if include_items else 0,
        "include_items": include_items,
        "random_slope": random_slope_sd > 0,
        "row_count": n_rows,
        "outcome_unit": outcome_unit,
        "censored": False,
    }
    return DurationSimulation(
        simulation_version="0.1",
        family="duration",
        data=data,
        truth=truth,
        random_effects={"participant": participant_re, "item": item_re},
        design=design,
    )

Combine prepared duration data with approved lognormal priors.

Source code in src/gp3bayespy/duration.py

def specify_duration_model(
    prepared: DurationPrepared,
    baseline: float,
    intercept_scale: float = 1.0,
    coefficient_scale: float = 0.5,
    group_sd_scale: float = 1.0,
    residual_scale: float = 1.0,
    correlation_eta: float = 2.0,
    student_df: float = 3.0,
) -> DurationModelSpecification:
    """Combine prepared duration data with approved lognormal priors."""
    if not isinstance(prepared, DurationPrepared):
        raise GP3BayesError("`prepared` must inherit from `gp3bayes_duration_prepared`.")
    _validate_duration_contract(prepared.contract)
    if not prepared.audit.ready:
        raise GP3BayesError("`prepared$audit` is not ready for specification.")

    priors = create_prior_specification(
        prepared.contract,
        baseline=baseline,
        intercept_scale=intercept_scale,
        coefficient_scale=coefficient_scale,
        group_sd_scale=group_sd_scale,
        residual_scale=residual_scale,
        correlation_eta=correlation_eta,
        student_df=student_df,
    )
    core = create_model_specification(prepared.contract, prepared.audit, priors)
    core_values = {field.name: getattr(core, field.name) for field in fields(ModelSpecification)}
    return DurationModelSpecification(
        **core_values,
        duration_workflow_version="0.1",
        prepared=prepared,
        fixed_formula=prepared.fixed_formula,
        fixed_formula_text=prepared.fixed_formula_text,
        model_matrix_columns=prepared.model_matrix_columns,
        outcome_unit=prepared.outcome_unit,
    )

Summarise duration posterior parameters and median-ratio transforms.

Source code in src/gp3bayespy/duration.py

def summarise_duration_posterior(
    fit: DurationFit,
    probability: float = 0.95,
    variables: Sequence[str] | str | None = None,
):
    """Summarise duration posterior parameters and median-ratio transforms."""
    if not isinstance(fit, DurationFit):
        raise GP3BayesError("`fit` must inherit from `gp3bayes_fit`.")
    return _summarise_duration(fit, probability=probability, variables=variables)

Translate an approved duration specification to the Python backend plan.

Source code in src/gp3bayespy/duration.py

def translate_duration_model_to_brms(
    specification: DurationModelSpecification,
) -> DurationBackendSpecification:
    """Translate an approved duration specification to the Python backend plan."""
    specification = _validate_duration_model_specification(specification)
    parameter_table = _translation_parameter_table(
        specification.priors,
        include_sigma=True,
        random_slope=specification.contract.random_slope,
    )
    prior_text = {
        str(row.parameter_class): str(row.prior) for row in parameter_table.itertuples(index=False)
    }
    return DurationBackendSpecification(
        translation_version="0.1",
        family="duration",
        model_family="hierarchical_lognormal_duration",
        formula=specification.formula,
        formula_text=specification.formula_text,
        family_object={"family": "lognormal", "link": "identity"},
        priors=prior_text,
        prior_text=prior_text,
        validated_priors=parameter_table.copy(),
        parameter_table=parameter_table,
        specification=specification,
        outcome_unit=specification.outcome_unit,
        backend_interface="pymc",
        sampling_backend="pymc",
        algorithm="NUTS",
        backend_available=_pymc_available(),
    )