feat: meta-analysis composites, AIPW link selection, and a version lock

Three things, all from auditing what was actually missing rather than
what the names suggested.

Meta-analysis composites. The arithmetic already lived in
morie.fn._ca_crim as separate primitives; what was missing is the
composite surface rmorie exposes. meta_pool(), meta_effect_sizes() and
meta_convert() compose those primitives rather than restating the
formulas, and meta_pool returns the heterogeneity beside the pooled
mean: reporting one without the other is how a synthesis overstates its
case. Parity against morie_meta_pool, morie_meta_effect_sizes and
morie_meta_convert is IDENTICAL to ten decimals across six blocks,
including a deliberately heterogeneous set where I-squared is 96.84 and
tau-squared 0.2806 -- values that a wrong weight or a truncation error
would move.

AIPW. mrm_estimate_causal_effect left the logistic default in place on a
continuous outcome, which failed inside the native core with "binary
only" and which I had reported as a missing capability in morie. It was
not: estimate_aipw takes outcome_model="linear". The link is now chosen
from the outcome, and all four estimators run:

    matching  +0.6939   ipw ate  +0.8390
    aipw      +0.8290   dml plr  +0.8486
    consensus +0.8082 (se 0.0565)          truth 0.80

Version lock. CITATION.cff said 1.2.2 while pyproject.toml and
r-package/morie/DESCRIPTION said 1.2.3 -- a citation pointing at a
version that never contained the work. Fixed, and tests now read all
three files plus the installed metadata so the next drift fails instead
of being noticed later.

Also corrected three more call sites I had assumed rather than read:
mean_effect_size and q_statistic return records, not bare numbers;
risk_ratio takes the four cells; se_log_rr takes two risks and two group
sizes; se_d_from_se_r takes (r, se_r).

11 meta tests, 4 version tests, 58 across the new files, all passing.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01GVZR4QjrTXCv55j4qdo6E2

Author: rootcoder007

CITATION.cff

@@ -8,7 +8,7 @@ authors:
     affiliation: "Centre for Criminology and Sociolegal Studies, University of Toronto"
     name-particle: "Singh"
     # IPA: /ʋənʃ sɪŋ ɾʊˈɦeːlɑː/
-version: "1.2.2"
+version: "1.2.3"
 date-released: "2026-09-08"
 url: "https://github.com/rootcoder007/morie"
 repository-code: "https://github.com/rootcoder007/morie"

src/morie/__init__.py

@@ -157,6 +157,9 @@
     "mrm_reconcile": "mrm_flagship",
     "mrm_report": "mrm_flagship",
     "mrm_estimate_causal_effect": "mrm_flagship",
+    "meta_pool": "meta_analysis",
+    "meta_effect_sizes": "meta_analysis",
+    "meta_convert": "meta_analysis",
     "causal_dag": "mrm_graphs",
     "mrm_dags": "mrm_graphs",
     "mrm_check_balancing": "mrm_diagnostics",
@@ -413,6 +416,9 @@ def load_sample(name: str):
     "mrm_reconcile",
     "mrm_report",
     "mrm_estimate_causal_effect",
+    "meta_pool",
+    "meta_effect_sizes",
+    "meta_convert",
     "causal_dag",
     "mrm_dags",
     # Tier 1 diagnostics

src/morie/meta_analysis.py

@@ -0,0 +1,184 @@
+# SPDX-License-Identifier: AGPL-3.0-or-later
+"""Meta-analysis composites: pool a set of estimates, and move between scales.
+
+Parity with rmorie's `morie_meta_pool`, `morie_meta_effect_sizes` and
+`morie_meta_convert` (R/ca_crim_native.R).
+
+The arithmetic already lives in `morie.fn._ca_crim` as separate
+primitives -- inverse-variance weights, Q, I-squared, the
+DerSimonian-Laird tau-squared, and the effect-size conversions. What was
+missing is the composite surface that does a whole pooling in one call,
+which is how the R side is used, so this module composes them rather
+than restating the formulas.
+
+Why both a fixed and a random-effects answer come back together: the
+fixed-effect weights assume every study estimates the same quantity, and
+I-squared is the evidence about whether that holds. Reporting the pooled
+mean without the heterogeneity beside it is the standard way to overstate
+a meta-analysis.
+"""
+
+from __future__ import annotations
+
+import math
+
+from morie.fn import _ca_crim as _cc
+
+__all__ = ["meta_pool", "meta_effect_sizes", "meta_convert"]
+
+
+def _num(x, what):
+    try:
+        v = [float(i) for i in x]
+    except TypeError:
+        v = [float(x)]
+    if not v:
+        raise ValueError("`%s` must not be empty" % what)
+    for i in v:
+        if i != i:
+            raise ValueError("`%s` must not contain missing values" % what)
+    return v
+
+
+def meta_pool(ys, ses, z_cv: float = 1.96, groups=None) -> dict:
+    """Pool estimates by inverse variance, with the heterogeneity beside it.
+
+    Returns the fixed-effect mean and its interval, Cochran's Q, I-squared,
+    the DerSimonian-Laird tau-squared and the random-effects weights that
+    follow from it. With `groups`, Q is split into within and between
+    components, which is how a subgroup claim is actually tested.
+    """
+    y = _num(ys, "ys")
+    s = _num(ses, "ses")
+    if len(y) != len(s):
+        raise ValueError("`ys` and `ses` must be the same length")
+    if any(i <= 0 for i in s):
+        raise ValueError("`ses` must be strictly positive")
+
+    w = [_cc.fixed_effect_weight(i) for i in s]
+    # The primitives return records, not bare numbers: mean_effect_size
+    # already carries the standard error and z, and q_statistic carries
+    # its own df. Take those rather than recomputing them here, so there
+    # is one definition of each quantity.
+    fe = _cc.mean_effect_size(y, w)
+    m = fe["mean"]
+    se_m = fe["se"]
+    qs = _cc.q_statistic(y, w)
+    q = qs["q"]
+    df = qs["df"]
+    tau2 = _cc.tau2_dersimonian_laird(y, w)
+    out = {
+        "weights": w,
+        "mean": m,
+        "se": se_m,
+        "z": fe["z"],
+        "ci": (m - z_cv * se_m, m + z_cv * se_m),
+        "q": q,
+        "df": df,
+        "i2": _cc.i_squared(q, df),
+        "tau2": tau2,
+        "weights_random": [_cc.random_effects_weight(i, tau2) for i in s],
+    }
+    if groups is not None:
+        g = list(groups)
+        if len(g) != len(y):
+            raise ValueError("`groups` must be the same length as `ys`")
+        ys_by, ws_by = [], []
+        for lab in dict.fromkeys(g):                 # first-seen order
+            idx = [i for i, v in enumerate(g) if v == lab]
+            ys_by.append([y[i] for i in idx])
+            ws_by.append([w[i] for i in idx])
+        qwb = _cc.q_within_between(ys_by, ws_by)
+        out["q_within"] = qwb["q_within"]
+        out["q_between"] = qwb["q_between"]
+        out["df_within"] = qwb["df_within"]
+        out["df_between"] = qwb["df_between"]
+    return out
+
+
+def meta_effect_sizes(m1=None, m2=None, s1=None, s2=None, n1=None, n2=None,
+                      t_value=None, a=None, b=None, c=None, d=None,
+                      r=None) -> dict:
+    """Effect sizes from whatever a paper actually reported.
+
+    Primary studies report means and SDs, or a t statistic, or a 2x2
+    table, or a correlation. Each of those determines an effect size and
+    its standard error, and only the ones the arguments support are
+    returned -- an absent input yields an absent key rather than a
+    fabricated zero.
+    """
+    out: dict = {}
+    if n1 is not None and n2 is not None:
+        n1f, n2f = float(n1), float(n2)
+        if s1 is not None and m1 is not None and m2 is not None \
+                and s2 is not None:
+            out["s_pooled"] = _cc.pooled_sd(float(s1), float(s2), n1f, n2f)
+            out["d"] = _cc.cohens_d_sample(float(m1), float(m2), float(s1),
+                                           float(s2), n1f, n2f)
+        out["j"] = _cc.hedges_j(n1f, n2f)
+        if "d" in out:
+            out["g"] = _cc.hedges_g(out["d"], n1f, n2f)
+            out["se_g"] = _cc.se_g(out["g"], n1f, n2f)
+        if t_value is not None:
+            out["d_from_t"] = _cc.d_from_t(float(t_value), n1f, n2f)
+    if a is not None and None not in (b, c, d):
+        af, bf, cf, df_ = float(a), float(b), float(c), float(d)
+        p1 = af / (af + bf)
+        p2 = cf / (cf + df_)
+        out["rr"] = _cc.risk_ratio(af, bf, cf, df_)
+        out["or"] = _cc.odds_ratio_2x2(af, bf, cf, df_)
+        # se_log_rr takes the two risks and the two group sizes, not the
+        # four cells
+        out["se_ln_rr"] = _cc.se_log_rr(p1, p2, af + bf, cf + df_)
+        out["se_ln_or"] = _cc.se_log_or(af, bf, cf, df_)
+    if r is not None:
+        rf = float(r)
+        out["fisher_z"] = _cc.fisher_z(rf)
+        if n1 is not None:
+            out["se_fisher_z"] = _cc.se_fisher_z(float(n1))
+    return out
+
+
+def meta_convert(ln_or=None, se_ln_or=None, p1=None, p2=None, n1=None,
+                 n2=None, d=None, se_d=None, rr=None, or_value=None,
+                 r=None, se_r=None) -> dict:
+    """Move an effect between the scales a synthesis has to mix.
+
+    A review rarely gets one scale. Log odds ratios, standardised mean
+    differences, probit differences and correlations all convert, and the
+    conversion constants are the ones the literature uses: the logistic
+    SD (pi/sqrt(3)) for the logit route and Cox's 1.65.
+    """
+    sd_logistic = math.sqrt(math.pi ** 2 / 3.0)
+    out: dict = {"sd_logistic": sd_logistic}
+    if ln_or is not None:
+        out["d_logit"] = _cc.d_from_log_or(float(ln_or), method="logit")
+        out["d_cox"] = _cc.d_from_log_or(float(ln_or), method="cox")
+    if se_ln_or is not None:
+        out["se_d_logit"] = math.sqrt(float(se_ln_or) ** 2 / sd_logistic ** 2)
+        out["se_d_cox"] = math.sqrt(float(se_ln_or) ** 2 / 1.65 ** 2)
+    if p1 is not None and p2 is not None:
+        out["d_probit"] = _cc.d_probit(float(p1), float(p2))
+        if n1 is not None and n2 is not None:
+            out["se_d_probit"] = _cc.se_d_probit(float(p1), float(p2),
+                                                 float(n1), float(n2))
+    if d is not None:
+        dv = float(d)
+        out["ln_or_logit"] = _cc.log_or_from_d(dv, method="logit")
+        out["ln_or_cox"] = _cc.log_or_from_d(dv, method="cox")
+        out["r_from_d"] = _cc.r_from_d(dv, float(n1), float(n2)) \
+            if (n1 is not None and n2 is not None) else _cc.r_from_d(dv)
+        if se_d is not None:
+            out["se_ln_or_logit"] = _cc.se_log_or_from_se_d(float(se_d),
+                                                            method="logit")
+    if rr is not None and p2 is not None:
+        out["or_from_rr"] = _cc.or_from_rr(float(rr), float(p2))
+    if or_value is not None and p2 is not None:
+        out["rr_from_or"] = _cc.rr_from_or(float(or_value), float(p2))
+    if r is not None:
+        out["fisher_z"] = _cc.fisher_z(float(r))
+        out["d_from_r"] = _cc.d_from_r_pointbiserial(float(r))
+        if se_r is not None:
+            out["se_d_from_se_r"] = _cc.se_d_from_se_r(float(r),
+                                                       float(se_r))
+    return out

src/morie/mrm_flagship.py

@@ -436,8 +436,16 @@ def _ate():
     if "aipw" in methods:
         def _aipw():
             from morie.causal import estimate_aipw
+            # The outcome model has to match the outcome. Leaving the
+            # logistic default on a continuous outcome fails inside the
+            # native core with "binary only", which reads like a missing
+            # capability and is really the wrong link function.
+            yvals = {r[outcome] for r in rows_in}
+            binary = yvals <= {0, 1, 0.0, 1.0, True, False}
             a = estimate_aipw(frame, treatment=treatment, outcome=outcome,
-                              covariates=covariates)
+                              covariates=covariates,
+                              outcome_model="logistic" if binary
+                              else "linear")
             return (a.get("ate", a.get("estimate")),
                     a.get("se", a.get("std_error")))
         attempt("aipw (morie native)", _aipw)

tests/test_meta_analysis.py

@@ -0,0 +1,133 @@
+# SPDX-License-Identifier: AGPL-3.0-or-later
+"""Meta-analysis composites, at parity with rmorie.
+
+The pooling anchors are closed-form: inverse-variance weights, Cochran's
+Q, I-squared and the DerSimonian-Laird tau-squared are all computable by
+hand from the inputs, so these assertions can fail.
+"""
+
+import math
+
+import pytest
+
+import morie
+
+YS = [0.2, 0.4, 0.3]
+SES = [0.1, 0.15, 0.12]
+
+
+def test_fixed_effect_pooling_is_inverse_variance():
+    p = morie.meta_pool(YS, SES)
+    w = [1 / s ** 2 for s in SES]
+    want = sum(wi * y for wi, y in zip(w, YS)) / sum(w)
+    assert p["mean"] == pytest.approx(want, abs=1e-12)
+    assert p["se"] == pytest.approx(math.sqrt(1 / sum(w)), abs=1e-12)
+    assert p["z"] == pytest.approx(p["mean"] / p["se"], abs=1e-10)
+    # the interval is symmetric about the mean at the given critical value
+    lo, hi = p["ci"]
+    assert lo == pytest.approx(p["mean"] - 1.96 * p["se"], abs=1e-12)
+    assert hi == pytest.approx(p["mean"] + 1.96 * p["se"], abs=1e-12)
+
+
+def test_cochran_q_and_i_squared_by_hand():
+    p = morie.meta_pool(YS, SES)
+    w = [1 / s ** 2 for s in SES]
+    m = sum(wi * y for wi, y in zip(w, YS)) / sum(w)
+    q = sum(wi * (y - m) ** 2 for wi, y in zip(w, YS))
+    assert p["q"] == pytest.approx(q, abs=1e-10)
+    assert p["df"] == len(YS) - 1
+    assert p["i2"] == pytest.approx(max(0.0, (q - p["df"]) / q * 100),
+                                    abs=1e-10)
+
+
+def test_homogeneous_studies_get_no_random_effects_variance():
+    # Q below its df means no more spread than sampling explains, so
+    # DerSimonian-Laird truncates tau-squared at zero and the random
+    # weights collapse onto the fixed ones.
+    p = morie.meta_pool(YS, SES)
+    assert p["tau2"] == 0.0
+    assert p["i2"] == 0.0
+    assert p["weights_random"] == pytest.approx(p["weights"])
+
+
+def test_heterogeneous_studies_are_flagged_rather_than_averaged_away():
+    # THE SUBSTANCE: reporting a pooled mean without the heterogeneity
+    # beside it is how a meta-analysis overstates its case. These are
+    # rmorie's values for the same input.
+    h = morie.meta_pool([0.1, 0.9, 0.45, 1.4], [0.08, 0.10, 0.09, 0.12])
+    assert h["mean"] == pytest.approx(0.5753135200, abs=1e-9)
+    assert h["q"] == pytest.approx(95.0109994464, abs=1e-8)
+    assert h["i2"] == pytest.approx(96.8424708534, abs=1e-8)
+    assert h["tau2"] == pytest.approx(0.2806412488, abs=1e-9)
+    # with real heterogeneity the random-effects weights must be smaller
+    assert all(r < f for r, f in zip(h["weights_random"], h["weights"]))
+
+
+def test_subgroups_split_q_into_within_and_between():
+    g = morie.meta_pool([0.2, 0.4, 0.3, 0.5], [0.1, 0.15, 0.12, 0.2],
+                        groups=["a", "a", "b", "b"])
+    assert g["q_within"] == pytest.approx(1.9660633484, abs=1e-9)
+    assert g["q_between"] == pytest.approx(0.4770891064, abs=1e-9)
+    # the split must account for the total
+    assert g["q_within"] + g["q_between"] == pytest.approx(g["q"], abs=1e-9)
+
+
+def test_pooling_checks_its_inputs():
+    with pytest.raises(ValueError, match="same length"):
+        morie.meta_pool([0.1, 0.2], [0.1])
+    with pytest.raises(ValueError, match="strictly positive"):
+        morie.meta_pool([0.1], [0.0])
+    with pytest.raises(ValueError, match="must not be empty"):
+        morie.meta_pool([], [])
+    with pytest.raises(ValueError, match="groups"):
+        morie.meta_pool(YS, SES, groups=["a"])
+
+
+def test_effect_sizes_from_means_and_sds():
+    # rmorie's values for the same input
+    e = morie.meta_effect_sizes(m1=10, m2=8, s1=2, s2=2.5, n1=30, n2=32)
+    assert e["s_pooled"] == pytest.approx(2.2721135535, abs=1e-9)
+    assert e["d"] == pytest.approx(0.8802376963, abs=1e-9)
+    assert e["g"] == pytest.approx(0.8691886875, abs=1e-9)
+    assert e["se_g"] == pytest.approx(0.2658495559, abs=1e-9)
+    # Hedges' g is the small-sample correction of d, so it shrinks
+    assert abs(e["g"]) < abs(e["d"])
+
+
+def test_effect_sizes_from_a_two_by_two_table():
+    e = morie.meta_effect_sizes(a=20, b=80, c=10, d=90)
+    assert e["rr"] == pytest.approx(2.0, abs=1e-10)
+    assert e["or"] == pytest.approx(2.25, abs=1e-10)
+    assert e["se_ln_rr"] == pytest.approx(0.3605551275, abs=1e-9)
+    assert e["se_ln_or"] == pytest.approx(0.4166666667, abs=1e-9)
+
+
+def test_absent_inputs_produce_absent_keys_not_zeros():
+    # A fabricated zero is worse than a missing key: it looks like a
+    # measured null.
+    e = morie.meta_effect_sizes(n1=30, n2=32)
+    assert "j" in e
+    assert "d" not in e and "g" not in e and "rr" not in e
+    assert morie.meta_effect_sizes() == {}
+
+
+def test_conversion_constants_are_the_ones_the_literature_uses():
+    cv = morie.meta_convert(ln_or=0.7, se_ln_or=0.2)
+    assert cv["sd_logistic"] == pytest.approx(math.pi / math.sqrt(3.0),
+                                              abs=1e-12)
+    assert cv["sd_logistic"] == pytest.approx(1.8137993642, abs=1e-9)
+    assert cv["d_logit"] == pytest.approx(0.3859302268, abs=1e-9)
+    assert cv["d_cox"] == pytest.approx(0.4242424242, abs=1e-9)
+    assert cv["se_d_logit"] == pytest.approx(0.1102657791, abs=1e-9)
+    # the Cox route divides by 1.65 by definition
+    assert cv["d_cox"] == pytest.approx(0.7 / 1.65, abs=1e-12)
+
+
+def test_probit_and_correlation_routes():
+    cv = morie.meta_convert(p1=0.6, p2=0.4, n1=100, n2=100)
+    assert "d_probit" in cv and "se_d_probit" in cv
+    assert cv["d_probit"] > 0                     # p1 > p2
+    cr = morie.meta_convert(r=0.3, se_r=0.05)
+    assert cr["fisher_z"] == pytest.approx(0.5 * math.log(1.3 / 0.7),
+                                           abs=1e-12)
+    assert "d_from_r" in cr and "se_d_from_se_r" in cr