`statspai.mendelian`¶

mendelian ¶

Mendelian Randomization methods.

Uses genetic variants as instrumental variables to estimate causal effects in epidemiological/health economics studies.

Core estimators (mr): - IVW (inverse-variance weighted) - MR-Egger regression - Weighted / penalized median

Diagnostics (diagnostics): - Cochran's Q / Rucker's Q' (heterogeneity) - MR-Egger intercept test (directional pleiotropy) - Leave-one-out sensitivity - Steiger directionality test - MR-PRESSO global test + outlier correction - Radial MR

Frontier (frontier, v1.6): - MR-Lap (Burgess-Davies-Thompson 2016 sample-overlap correction) - MR-Clust (Foley et al. 2021 clustered pleiotropy) - GRAPPLE (Wang et al. 2021 profile-likelihood MR) - MR-cML-BIC (Xue-Shen-Pan 2021 constrained maximum likelihood)

MRResult ¶

Bases: ResultProtocolMixin

Results from Mendelian Randomization analysis.

Returned by :func:mendelian_randomization. Holds the per-method estimates table plus heterogeneity (Cochran's Q) and pleiotropy (MR-Egger intercept) diagnostics; render with .summary() or .plot().

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> n_snps = 15
>>> bx = rng.uniform(0.1, 0.5, n_snps)
>>> by = 0.3 * bx + rng.normal(0, 0.02, n_snps)
>>> snp_stats = pd.DataFrame({
...     "beta_x": bx, "beta_y": by,
...     "se_x": rng.uniform(0.01, 0.05, n_snps),
...     "se_y": rng.uniform(0.01, 0.05, n_snps),
... })
>>> res = sp.mendelian_randomization(
...     data=snp_stats, beta_exposure="beta_x", beta_outcome="beta_y",
...     se_exposure="se_x", se_outcome="se_y",
...     exposure_name="BMI", outcome_name="T2D", seed=0)
>>> isinstance(res, sp.MRResult)
True
>>> res.n_snps
15
>>> list(res.estimates["method"])
['IVW', 'MR-Egger', 'Weighted Median']

plot ¶

plot(ax: Any = None, **kwargs: Any) -> Any

Scatter plot of SNP effects with MR lines.

HeterogeneityResult `dataclass` ¶

Bases: ResultProtocolMixin

Container for Cochran's Q / Rücker's Q' heterogeneity diagnostics.

Returned by :func:mr_heterogeneity. Holds the Q statistic, its degrees of freedom and p-value, the I^2 percentage, and the fitting method ("ivw" or "egger").

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 12
>>> bx = rng.uniform(0.2, 0.5, n)
>>> by = 0.3 * bx + rng.normal(0, 0.02, n)
>>> sy = rng.uniform(0.02, 0.05, n)
>>> het = sp.mr_heterogeneity(bx, by, sy)
>>> isinstance(het, sp.HeterogeneityResult)
True
>>> het.method
'ivw'
>>> bool(het.Q >= 0.0)
True

PleiotropyResult `dataclass` ¶

Bases: ResultProtocolMixin

Container for the MR-Egger intercept (directional-pleiotropy) test.

Returned by :func:mr_pleiotropy_egger. Holds the fitted Egger intercept, its standard error se, and the two-sided p_value for H0: intercept == 0.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 12
>>> bx = rng.uniform(0.2, 0.5, n)
>>> by = 0.3 * bx + rng.normal(0, 0.02, n)
>>> sy = rng.uniform(0.02, 0.05, n)
>>> pleio = sp.mr_pleiotropy_egger(bx, by, sy)
>>> isinstance(pleio, sp.PleiotropyResult)
True
>>> bool(0.0 <= pleio.p_value <= 1.0)
True

LeaveOneOutResult `dataclass` ¶

Bases: ResultProtocolMixin

Container for leave-one-out IVW estimates.

Returned by :func:mr_leave_one_out. Wraps a :class:pandas.DataFrame table with one row per dropped SNP and columns dropped_snp, estimate, se, ci_lower, ci_upper and p_value.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 12
>>> bx = rng.uniform(0.2, 0.5, n)
>>> by = 0.3 * bx + rng.normal(0, 0.02, n)
>>> sy = rng.uniform(0.02, 0.05, n)
>>> loo = sp.mr_leave_one_out(bx, by, sy)
>>> isinstance(loo, sp.LeaveOneOutResult)
True
>>> len(loo.table)
12

SteigerResult `dataclass` ¶

Bases: ResultProtocolMixin

Container for the Steiger directionality test.

Returned by :func:mr_steiger. Reports whether the assumed causal direction is supported (correct_direction), the one-sided steiger_pvalue, the variance explained on each trait (r2_exposure, r2_outcome), and the effective sample sizes.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 12
>>> bx = rng.uniform(0.2, 0.5, n)
>>> by = 0.3 * bx + rng.normal(0, 0.02, n)
>>> sx = rng.uniform(0.02, 0.05, n)
>>> sy = rng.uniform(0.02, 0.05, n)
>>> st = sp.mr_steiger(bx, sx, 50000, by, sy, 50000)
>>> isinstance(st, sp.SteigerResult)
True
>>> isinstance(st.correct_direction, bool)
True

RadialResult `dataclass` ¶

Bases: ResultProtocolMixin

Container for radial-MR diagnostics.

Returned by :func:mr_radial. Wraps a per-SNP table (columns snp, W, beta_hat, q_contribution), the total_Q statistic with its Q_pvalue, and the list of Bonferroni-flagged outliers.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 12
>>> bx = rng.uniform(0.2, 0.5, n)
>>> by = 0.3 * bx + rng.normal(0, 0.02, n)
>>> sy = rng.uniform(0.02, 0.05, n)
>>> rad = sp.mr_radial(bx, by, sy)
>>> isinstance(rad, sp.RadialResult)
True
>>> len(rad.table)
12

ModeBasedResult `dataclass` ¶

Bases: ResultProtocolMixin

Result container returned by :func:mr_mode.

Holds the mode-based MR point estimate, bootstrap standard error, confidence interval, p-value, the number of SNPs, and the Silverman bandwidth used for the kernel density. Call :meth:summary for a formatted text report.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> bx = rng.uniform(0.2, 0.5, 25)
>>> by = 0.3 * bx + rng.normal(0, 0.01, 25)
>>> sx = rng.uniform(0.02, 0.04, 25)
>>> sy = rng.uniform(0.02, 0.04, 25)
>>> res = sp.mr_mode(bx, by, sx, sy, n_boot=200, seed=0)
>>> isinstance(res, sp.ModeBasedResult)
True
>>> res.method
'weighted-mode'

FStatisticResult `dataclass` ¶

Bases: ResultProtocolMixin

Result container returned by :func:mr_f_statistic.

Summarises instrument strength across SNPs: the mean / min / max per-SNP F-statistic, mean explained variance r2_mean, the full per_snp_F array, and a weak_instrument_risk flag that is True when any SNP's F falls below 10 (Staiger-Stock 1997).

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> bx = rng.uniform(0.2, 0.5, 25)
>>> sx = rng.uniform(0.02, 0.04, 25)
>>> f = sp.mr_f_statistic(bx, sx)
>>> isinstance(f, sp.FStatisticResult)
True
>>> f.per_snp_F.shape
(25,)

MVMRResult `dataclass` ¶

Bases: ResultProtocolMixin

Multivariable MR output.

Examples:

>>> import numpy as np, pandas as pd, statspai as sp
>>> rng = np.random.default_rng(0)
>>> bx1 = rng.uniform(0.1, 0.5, 30); bx2 = rng.uniform(0.1, 0.5, 30)
>>> snp = pd.DataFrame({"beta_x1": bx1, "beta_x2": bx2,
...                     "beta_y": 0.4 * bx1 + rng.normal(0, 0.02, 30),
...                     "se_y": rng.uniform(0.01, 0.05, 30)})
>>> res = sp.mr_multivariable(snp)
>>> isinstance(res, sp.MVMRResult)
True

MediationMRResult `dataclass` ¶

Bases: ResultProtocolMixin

Two-step MR output.

Examples:

>>> import numpy as np, pandas as pd, statspai as sp
>>> rng = np.random.default_rng(1)
>>> bx = rng.uniform(0.1, 0.5, 30); bm = 0.5 * bx + rng.normal(0, 0.02, 30)
>>> se = lambda: rng.uniform(0.01, 0.05, 30)
>>> snp = pd.DataFrame({"beta_x": bx, "se_x": se(), "beta_m": bm, "se_m": se(),
...     "beta_y": 0.2 * bx + 0.6 * bm + rng.normal(0, 0.02, 30), "se_y": se()})
>>> res = sp.mr_mediation(snp)
>>> isinstance(res, sp.MediationMRResult)
True

MRBMAResult `dataclass` ¶

Bases: ResultProtocolMixin

MR-BMA (Bayesian Model Averaging) output.

Examples:

>>> import numpy as np, pandas as pd, statspai as sp
>>> rng = np.random.default_rng(0)
>>> bx1 = rng.uniform(0.1, 0.5, 30); bx2 = rng.uniform(0.1, 0.5, 30)
>>> snp = pd.DataFrame({"beta_x1": bx1, "beta_x2": bx2,
...                     "beta_y": 0.4 * bx1 + rng.normal(0, 0.02, 30),
...                     "se_y": rng.uniform(0.01, 0.05, 30)})
>>> res = sp.mr_bma(snp)
>>> isinstance(res, sp.MRBMAResult)
True

MRLapResult `dataclass` ¶

Bases: ResultProtocolMixin

Output of :func:mr_lap.

Attributes:

Name	Type	Description
`estimate`	`float`	Overlap-corrected causal estimate.
`se`	`float`	First-order SE (identical to IVW SE; overlap correction does not change the variance under the Burgess-Davies-Thompson 2016 approximation).
`ci_lower, ci_upper`	`float`
`p_value`	`float`
`estimate_ivw`	`float`	The uncorrected IVW point estimate (for comparison).
`bias_correction`	`float`	`estimate_ivw - estimate`, i.e. the amount subtracted from IVW.
`overlap_fraction`	`float`
`overlap_rho`	`float`
`f_mean`	`float`	Mean first-stage F-statistic across SNPs.
`n_snps`	`int`

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> rng = np.random.default_rng(8)
>>> bx = rng.uniform(0.1, 0.5, 25)
>>> sx = rng.uniform(0.02, 0.05, 25)
>>> sy = rng.uniform(0.02, 0.05, 25)
>>> by = 0.3 * bx + rng.normal(0, sy)
>>> res = sp.mr_lap(bx, by, sx, sy, overlap_fraction=0.5, overlap_rho=0.1)
>>> type(res).__name__
'MRLapResult'
>>> res.n_snps
25
>>> bool(np.isfinite(res.estimate))
True

MRClustResult `dataclass` ¶

Bases: ResultProtocolMixin

Output of :func:mr_clust.

Attributes:

Name	Type	Description
`cluster_estimates`	`DataFrame`	One row per cluster, columns `[cluster, estimate, se, ci_lower, ci_upper, weight, n_snps]`. Cluster 0 is a 'null' component at θ=0 in the Foley et al. 2021 specification.
`assignments`	`ndarray`	Integer cluster id (0..K-1) for each SNP.
`responsibilities`	`ndarray`	(n_snps, K) matrix of posterior cluster-membership probabilities.
`K`	`int`	Selected number of clusters (including the null cluster).
`bic`	`dict`	BIC per K tried; keys are ints.
`loglik`	`float`	Final log-likelihood at the selected K.
`wald_ratios`	`ndarray`	Per-SNP Wald ratio (β_y / β_x).
`wald_se`	`ndarray`	Per-SNP Wald-ratio first-order SE.
`n_snps`	`int`

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> rng = np.random.default_rng(8)
>>> bx = rng.uniform(0.1, 0.5, 25)
>>> sx = rng.uniform(0.02, 0.05, 25)
>>> sy = rng.uniform(0.02, 0.05, 25)
>>> by = 0.3 * bx + rng.normal(0, sy)
>>> res = sp.mr_clust(bx, by, sx, sy)
>>> type(res).__name__
'MRClustResult'
>>> res.n_snps
25
>>> bool(len(res.assignments) == 25)
True

GrappleResult `dataclass` ¶

Bases: ResultProtocolMixin

Output of :func:grapple.

Attributes:

Name	Type	Description
`estimate`	`float`	Robust profile-score estimate of the causal effect β.
`se`	`float`	Sandwich standard error of the estimate (S2 Text).
`ci_lower, ci_upper`	`float`
`p_value`	`float`
`tau2`	`float`	Estimated pleiotropy variance τ² (balanced pleiotropy SD² = τ²).
`loglik`	`float`	Robust objective `-sum rho(t)` at the estimate.
`tau2_se`	`float`
`loss`	`str`
`converged`	`bool`
`n_snps`	`int`

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> rng = np.random.default_rng(8)
>>> bx = rng.uniform(0.1, 0.5, 25)
>>> sx = rng.uniform(0.02, 0.05, 25)
>>> sy = rng.uniform(0.02, 0.05, 25)
>>> by = 0.3 * bx + rng.normal(0, sy)
>>> res = sp.grapple(bx, by, sx, sy)
>>> type(res).__name__
'GrappleResult'
>>> res.n_snps
25
>>> bool(np.isfinite(res.estimate))
True

MRcMLResult `dataclass` ¶

Bases: ResultProtocolMixin

Output of :func:mr_cml.

Attributes:

Name	Type	Description
`estimate`	`float`	cML-BIC estimate (or the BIC model average when `model_average=True`).
`se`	`float`	Standard error from the Fisher information at the selected K (or the model-averaged SE).
`ci_lower, ci_upper, p_value`	`float`
`K_selected`	`int`	Number of variants flagged as invalid (pleiotropic) by BIC.
`invalid_snps`	`ndarray`	Boolean mask of the K_selected variants flagged as invalid.
`path`	`DataFrame`	One row per K in [0, K_max] with columns `[K, estimate, se, loglik, bic, converged]`.
`loglik`	`float`	Maximised constrained Gaussian log-likelihood at K_selected (Equation 2 including its normalising constant).
`n_snps`	`int`
`n_samples`	`int or None`	GWAS sample size used in the BIC penalty (`None` if not given, in which case the number of variants stood in and a warning was issued).
`model_average`	`bool`	Whether `estimate` / `se` are BIC model-averaged over K.

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> rng = np.random.default_rng(8)
>>> bx = rng.uniform(0.1, 0.5, 25)
>>> sx = rng.uniform(0.02, 0.05, 25)
>>> sy = rng.uniform(0.02, 0.05, 25)
>>> by = 0.3 * bx + rng.normal(0, sy)
>>> res = sp.mr_cml(bx, by, sx, sy, n=50_000)
>>> type(res).__name__
'MRcMLResult'
>>> res.n_snps
25
>>> bool(np.isfinite(res.estimate))
True

MRRapsResult `dataclass` ¶

Bases: ResultProtocolMixin

Output of :func:mr_raps.

Attributes:

Name	Type	Description
`estimate`	`float`	Causal effect estimate.
`se`	`float`	Asymptotic standard error: Theorem 3.3 of Zhao et al. (2020) for the model without over-dispersion, the `beta` entry of the joint `(beta, tau2)` sandwich (Theorems 4.2 / 5.2) otherwise.
`ci_lower, ci_upper`	`float`	Normal confidence interval at level `1 - alpha`.
`p_value`	`float`	Two-sided normal p-value for `beta = 0`.
`tau2`	`float`	Over-dispersion variance (0 for the model without it).
`tau2_se`	`float`	Standard error of `tau2` (`nan` without over-dispersion).
`loglik_robust`	`float`	Objective at the estimate, `-sum_j rho(t_j)` (eq. 5.1, which is the profile log-likelihood (3.2) when `rho(r) = r^2 / 2`).
`converged`	`bool`	Whether the estimating equations were solved within `niter` iterations to `tol`.
`tuning_c`	`float`	Loss tuning constant (`nan` for `loss="l2"`).
`n_snps`	`int`	Variants used, after pruning.
`loss`	`str`
`over_dispersion`	`bool`
`n_pruned`	`int`	Variants dropped because `se_exposure > 10 * median(se_exposure)`.

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> rng = np.random.default_rng(8)
>>> bx = rng.uniform(0.1, 0.5, 25)
>>> sx = rng.uniform(0.02, 0.05, 25)
>>> sy = rng.uniform(0.02, 0.05, 25)
>>> by = 0.3 * bx + rng.normal(0, sy)
>>> res = sp.mr_raps(bx, by, sx, sy)
>>> type(res).__name__
'MRRapsResult'
>>> res.n_snps
25
>>> bool(np.isfinite(res.estimate))
True

mendelian_randomization ¶

mendelian_randomization(data: Optional[DataFrame] = None, beta_exposure: Optional[str] = None, beta_outcome: Optional[str] = None, se_exposure: Optional[str] = None, se_outcome: Optional[str] = None, exposure_name: str = 'Exposure', outcome_name: str = 'Outcome', methods: Optional[List[str]] = None, alpha: float = 0.05, seed: Optional[int] = None) -> MRResult

Mendelian Randomization analysis using summary statistics.

Equivalent to R's MendelianRandomization::mr_allmethods() and TwoSampleMR::mr().

Parameters:

Name	Type	Description	Default
`data`	`DataFrame`	Summary statistics with one row per SNP/instrument.	`None`
`beta_exposure`	`str`	Column name for SNP-exposure association beta.	`None`
`beta_outcome`	`str`	Column name for SNP-outcome association beta.	`None`
`se_exposure`	`str`	Column name for SNP-exposure SE.	`None`
`se_outcome`	`str`	Column name for SNP-outcome SE.	`None`
`exposure_name`	`str`		`'Exposure'`
`outcome_name`	`str`		`'Outcome'`
`methods`	`list of str`	MR methods to use. Default: ['ivw', 'egger', 'weighted_median'].	`None`
`alpha`	`float`		`0.05`
`seed`	`int`		`None`

Returns:

Type	Description
`MRResult`	Results with .summary(), .plot(), estimates table.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n_snps = 25
>>> bx = rng.uniform(0.05, 0.25, n_snps)          # SNP-exposure effects
>>> by = 0.4 * bx + rng.normal(0, 0.01, n_snps)   # causal effect ~ 0.4
>>> snp_stats = pd.DataFrame({
...     'beta_x': bx, 'beta_y': by,
...     'se_x': rng.uniform(0.01, 0.03, n_snps),
...     'se_y': rng.uniform(0.01, 0.03, n_snps),
... })
>>> result = sp.mendelian_randomization(
...     data=snp_stats,
...     beta_exposure='beta_x', beta_outcome='beta_y',
...     se_exposure='se_x', se_outcome='se_y',
...     exposure_name='BMI', outcome_name='T2D',
... )
>>> bool(result.n_snps == 25)
True
>>> summary_text = result.summary()
>>> ax = result.plot()

mr_egger ¶

mr_egger(beta_exposure: ndarray, beta_outcome: ndarray, se_exposure: ndarray, se_outcome: ndarray, alpha: float = 0.05) -> Dict[str, float]

MR-Egger regression.

Allows for directional pleiotropy via a non-zero intercept. Variants are oriented to a positive exposure association before fitting, and the standard errors use multiplicative random effects with the residual standard error floored at 1 -- both as in MendelianRandomization and TwoSampleMR, which this matches to machine precision.

.. versionchanged:: 1.28.0 Variants were not oriented, so the result depended on which allele each variant was coded against: on the LDL-C / CHD example shipped with MendelianRandomization the slope was 14% off and the intercept -- the pleiotropy test -- 38% off. The residual standard error was also not floored at 1, which shrinks the standard errors below the fixed-effect ones under under-dispersion.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> n_snps = 15
>>> beta_exposure = rng.uniform(0.1, 0.5, n_snps)
>>> beta_outcome = 0.3 * beta_exposure + rng.normal(0, 0.02, n_snps)
>>> se_exposure = rng.uniform(0.01, 0.05, n_snps)
>>> se_outcome = rng.uniform(0.01, 0.05, n_snps)
>>> res = sp.mr_egger(beta_exposure, beta_outcome, se_exposure, se_outcome)
>>> sorted(res.keys())
['ci_lower', 'ci_upper', 'estimate', 'intercept', 'intercept_p', 'intercept_se', 'p_value', 'se']
>>> bool("intercept" in res)
True

References

bowden2015mendelian

mr_ivw ¶

mr_ivw(beta_exposure: ndarray, beta_outcome: ndarray, se_exposure: ndarray, se_outcome: ndarray, alpha: float = 0.05, model: str = 'default') -> Dict[str, float]

Inverse-Variance Weighted (IVW) MR estimator.

Weighted regression of the variant-outcome associations on the variant-exposure associations through the origin (equivalently, an inverse-variance-weighted mean of the Wald ratios).

model follows MendelianRandomization::mr_ivw:

"fixed" -- standard error sqrt(1 / sum(w bx^2));
"random" -- multiplicative random effects: the fixed-effect standard error times max(1, RSE), RSE = sqrt(Q / (n - 1));
"default" -- random with more than three variants, fixed otherwise (the package's own rule; TwoSampleMR's mr_ivw is the same multiplicative random-effects form).

.. versionchanged:: 1.28.0 The standard error used to be the fixed-effect one regardless of heterogeneity. On MendelianRandomization's own LDL-cholesterol / CHD example (28 variants, Cochran's Q p = 3e-10) that is 0.276 against the reference's 0.530 -- about half. The point estimate is unchanged.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> n_snps = 15
>>> beta_exposure = rng.uniform(0.1, 0.5, n_snps)
>>> beta_outcome = 0.3 * beta_exposure + rng.normal(0, 0.02, n_snps)
>>> se_exposure = rng.uniform(0.01, 0.05, n_snps)
>>> se_outcome = rng.uniform(0.01, 0.05, n_snps)
>>> res = sp.mr_ivw(beta_exposure, beta_outcome, se_exposure, se_outcome)
>>> sorted(res.keys())
['I2', 'Q', 'Q_df', 'Q_p', 'ci_lower', 'ci_upper', 'estimate', 'model',
 'p_value', 'rse', 'se']
>>> bool(0.0 < res["estimate"] < 0.6)
True

References

burgess2013mendelian

mr_median ¶

mr_median(beta_exposure: ndarray, beta_outcome: ndarray, se_exposure: ndarray, se_outcome: ndarray, penalized: bool = False, n_boot: int = 1000, alpha: float = 0.05, seed: Optional[int] = None, weighting: Optional[str] = None) -> Dict[str, float]

Weighted median MR estimator.

Consistent when at least 50% of the weight comes from valid instruments. Matches MendelianRandomization::mr_median in all three weightings.

.. versionchanged:: 1.28.0 Three corrections, found against MendelianRandomization:

the weighted median was a step function (first ratio reaching half the weight) instead of Bowden et al.'s interpolated definition -- 0.8% on the package's LDL-C / CHD example;
the penalty used the LOWER-tail chi-square probability, so it up-weighted exactly the heterogeneous variants it exists to suppress, had no min(1, 20 q) cap, and centred on the IVW estimate rather than the weighted median -- 38% off;
the bootstrap recomputed the weights from every draw instead of holding them fixed.

Parameters:

Name	Type	Description	Default
`weighting`	`(weighted, simple, penalized)`	`"weighted"` uses `(bx / se_y)^2`; `"simple"` equal weights; `"penalized"` multiplies the weighted weights by `min(1, 20 q_j)`, where `q_j` is the upper-tail chi-square(1) p-value of variant j's heterogeneity about the weighted median. Defaults to `"penalized"` if `penalized=True` else `"weighted"`.	`"weighted"`
`penalized`	`bool`	Kept for backward compatibility; equivalent to `weighting="penalized"`.	`False`
`n_boot`	`int`	Parametric bootstrap for the standard error: draw `bx` and `by` from their normal sampling distributions and take the weighted median of the resampled ratios with the ORIGINAL weights, as the reference does. The standard error is stochastic; the estimate is not.	`1000`
`seed`	`int`	Parametric bootstrap for the standard error: draw `bx` and `by` from their normal sampling distributions and take the weighted median of the resampled ratios with the ORIGINAL weights, as the reference does. The standard error is stochastic; the estimate is not.	`1000`

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> n_snps = 15
>>> beta_exposure = rng.uniform(0.1, 0.5, n_snps)
>>> beta_outcome = 0.3 * beta_exposure + rng.normal(0, 0.02, n_snps)
>>> se_exposure = rng.uniform(0.01, 0.05, n_snps)
>>> se_outcome = rng.uniform(0.01, 0.05, n_snps)
>>> res = sp.mr_median(beta_exposure, beta_outcome, se_exposure,
...                    se_outcome, n_boot=200, seed=0)
>>> sorted(res.keys())
['ci_lower', 'ci_upper', 'estimate', 'p_value', 'se', 'weighting']
>>> bool(0.0 < res["estimate"] < 0.6)
True

References

bowden2016consistent

mr_plot ¶

mr_plot(result: Optional[MRResult] = None, ax: Any = None, **kwargs: Any) -> Any

Scatter plot of SNP-exposure vs SNP-outcome effects with MR lines.

mr_heterogeneity ¶

mr_heterogeneity(beta_exposure: ndarray, beta_outcome: ndarray, se_outcome: ndarray, *, method: str = 'ivw', se_exposure: Optional[ndarray] = None) -> HeterogeneityResult

Cochran's Q (IVW) or Rücker's Q' (Egger) for pleiotropic heterogeneity.

Parameters:

Name	Type	Description	Default
`beta_exposure`	`array - like`		required
`beta_outcome`	`array - like`		required
`se_outcome`	`array - like`		required
`method`	`('ivw', 'egger')`		`"ivw"`
`se_exposure`	`array - like`	Required for Rücker's Q' (Egger); if omitted we default to first-order weights and raise for `method='egger'`.	`None`

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 12
>>> beta_exposure = rng.uniform(0.2, 0.5, n)
>>> beta_outcome = 0.3 * beta_exposure + rng.normal(0, 0.02, n)
>>> se_outcome = rng.uniform(0.02, 0.05, n)
>>> het = sp.mr_heterogeneity(beta_exposure, beta_outcome, se_outcome)
>>> het.method
'ivw'
>>> het.Q_df  # n - 1 for IVW
11
>>> bool(0.0 <= het.I2 <= 100.0)
True

References

[@bowden2016assessing]

mr_pleiotropy_egger ¶

mr_pleiotropy_egger(beta_exposure: ndarray, beta_outcome: ndarray, se_outcome: ndarray) -> PleiotropyResult

Test the MR-Egger intercept for directional (unbalanced) pleiotropy.

H0: intercept == 0 (no directional pleiotropy).

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 12
>>> beta_exposure = rng.uniform(0.2, 0.5, n)
>>> beta_outcome = 0.3 * beta_exposure + rng.normal(0, 0.02, n)
>>> se_outcome = rng.uniform(0.02, 0.05, n)
>>> pleio = sp.mr_pleiotropy_egger(beta_exposure, beta_outcome, se_outcome)
>>> bool(pleio.se >= 0.0)
True
>>> bool(0.0 <= pleio.p_value <= 1.0)
True

mr_leave_one_out ¶

mr_leave_one_out(beta_exposure: ndarray, beta_outcome: ndarray, se_outcome: ndarray, *, snp_ids: Optional[List[str]] = None, alpha: float = 0.05, model: str = 'default') -> LeaveOneOutResult

IVW estimate with each SNP dropped in turn.

Identifies SNPs that drive the overall estimate disproportionately.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 12
>>> beta_exposure = rng.uniform(0.2, 0.5, n)
>>> beta_outcome = 0.3 * beta_exposure + rng.normal(0, 0.02, n)
>>> se_outcome = rng.uniform(0.02, 0.05, n)
>>> loo = sp.mr_leave_one_out(beta_exposure, beta_outcome, se_outcome)
>>> len(loo.table)  # one row per dropped SNP
12
>>> sorted(loo.table.columns)
['ci_lower', 'ci_upper', 'dropped_snp', 'estimate', 'p_value', 'se']

mr_steiger ¶

mr_steiger(beta_exposure: ndarray, se_exposure: ndarray, n_exposure: SampleSizeInput, beta_outcome: ndarray, se_outcome: ndarray, n_outcome: SampleSizeInput, *, eaf: Optional[ndarray] = None, alternative: str = 'two-sided') -> SteigerResult

Steiger test for the direction of the causal effect (Hemani 2017).

Compares the R^2 of the SNPs on the exposure with their R^2 on the outcome. If R^2_exposure > R^2_outcome the assumed causal direction (exposure -> outcome) is supported. Matches TwoSampleMR::mr_steiger with r from get_r_from_bsen.

.. versionchanged:: 1.28.0 The p-value was one-sided and computed as 1 - Phi(z), which is exactly 0.0 for z beyond ~8.3. It is now two-sided by default, as the method authors' TwoSampleMR reports it, computed from the upper tail directly; pass alternative="greater" for the one-sided test.

Parameters:

Name	Type	Description	Default
`alternative`	`('two-sided', 'greater')`	`"two-sided"` is the reference's p-value; `"greater"` tests the one-sided hypothesis that the exposure R^2 is the larger.	`"two-sided"`

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 12
>>> beta_exposure = rng.uniform(0.2, 0.5, n)
>>> beta_outcome = 0.3 * beta_exposure + rng.normal(0, 0.02, n)
>>> se_exposure = rng.uniform(0.02, 0.05, n)
>>> se_outcome = rng.uniform(0.02, 0.05, n)
>>> st = sp.mr_steiger(
...     beta_exposure, se_exposure, 50000,
...     beta_outcome, se_outcome, 50000,
... )
>>> isinstance(st.correct_direction, bool)
True
>>> bool(st.r2_exposure >= 0.0 and st.r2_outcome >= 0.0)
True

References

[@hemani2017orienting]

mr_radial ¶

mr_radial(beta_exposure: ndarray, beta_outcome: ndarray, se_outcome: ndarray, *, snp_ids: Optional[List[str]] = None, alpha: float = 0.05, bonferroni: bool = True) -> RadialResult

Radial IVW MR (Bowden et al. 2018).

Reparameterizes each SNP's Wald ratio as a coordinate in a "radial" space. SNPs whose individual chi-square contribution to the Cochran Q exceeds the threshold are flagged as outliers: the Bonferroni level alpha / n by default, or alpha itself with bonferroni=False, which is RadialMR::ivw_radial's rule. The per-variant contributions and total Q match RadialMR (first-order weights) to machine precision either way.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 12
>>> beta_exposure = rng.uniform(0.2, 0.5, n)
>>> beta_outcome = 0.3 * beta_exposure + rng.normal(0, 0.02, n)
>>> se_outcome = rng.uniform(0.02, 0.05, n)
>>> rad = sp.mr_radial(beta_exposure, beta_outcome, se_outcome)
>>> len(rad.table)  # one row per SNP
12
>>> bool(0.0 <= rad.Q_pvalue <= 1.0)
True
>>> isinstance(rad.outliers, list)
True

References

[@bowden2018improving]

mr_mode ¶

mr_mode(beta_exposure: ndarray, beta_outcome: ndarray, se_exposure: ndarray, se_outcome: ndarray, *, method: str = 'weighted', n_boot: int = 1000, alpha: float = 0.05, seed: Optional[int] = None, phi: float = 1.0, refine: bool = False) -> ModeBasedResult

Mode-based MR estimator (Hartwig, Davey Smith & Bowden 2017).

Matches MendelianRandomization::mr_mbe(stderror = "simple") to machine precision in both weightings: same bandwidth rule, same 512-point density grid, same bootstrap construction.

.. versionchanged:: 1.28.0 The bandwidth was a weighted-SD / IQR rule rather than Hartwig et al.'s 0.9 min(sd, mad) n^(-1/5) on the unweighted ratios, the mode was searched on a different grid, and the bootstrap redrew the exposure associations and recomputed the weights instead of resampling the ratios with the weights held fixed. On the reference package's LDL-C / CHD example the estimate moved 7%. The standard error is now the MAD of the bootstrap draws, as in the reference, rather than their SD.

Parameters:

Name	Type	Description	Default
`method`	`('weighted', 'simple')`	`weighted` uses IVW weights when finding the mode; `simple` uses unit weights (robust but less efficient).	`"weighted"`
`phi`	`float`	Bandwidth multiplier, the reference's `phi`.	`1.0`
`refine`	`bool`	The reference estimate is a point of its 512-point density grid, whose step is `(range + 6h) / 511` -- coarse when a few ratios are extreme. `True` climbs to the continuous mode of the same density. Leave `False` to reproduce `mr_mbe`.	`False`
`n_boot`	`int`	Bootstrap replicates for SE.	`1000`

References

hartwig2017robust

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> bx = rng.uniform(0.2, 0.5, 25)
>>> by = 0.3 * bx + rng.normal(0, 0.01, 25)
>>> sx = rng.uniform(0.02, 0.04, 25)
>>> sy = rng.uniform(0.02, 0.04, 25)
>>> res = sp.mr_mode(bx, by, sx, sy, n_boot=200, seed=0)
>>> res.n_snps
25
>>> bool(0.0 < res.estimate < 1.0)  # near the true 0.3
True

mr_f_statistic ¶

mr_f_statistic(beta_exposure: ndarray, se_exposure: ndarray, *, n_samples: Optional[int] = None) -> FStatisticResult

Per-SNP F-statistic as an instrument-strength summary.

Uses the standard summary-stat approximation F_i = (beta_i / se_i)^2, which is valid for large GWAS where the first-stage regression is asymptotically Normal. Flags weak-IV risk when any SNP's F falls below 10 (Staiger-Stock 1997).

Parameters:

Name	Type	Description	Default
`n_samples`	`int`	Exposure GWAS sample size. If provided, each variant's R^2 is `F / (F + n - 2)` and `r2_mean` their mean; otherwise `r2_mean` is NaN (R^2 needs the sample size; the F-statistics do not). .. versionchanged:: 1.41 Without `n_samples` the number of variants used to stand in for the sample size.	`None`

References

staiger1997instrumental

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> bx = rng.uniform(0.2, 0.5, 25)
>>> sx = rng.uniform(0.02, 0.04, 25)  # strong instruments
>>> f = sp.mr_f_statistic(bx, sx)
>>> bool(f.f_mean > 10)  # Staiger-Stock rule of thumb
True
>>> f.weak_instrument_risk
False

mr_funnel_plot ¶

mr_funnel_plot(beta_exposure: ndarray, beta_outcome: ndarray, se_outcome: ndarray, *, snp_ids: Optional[List[str]] = None, ax: Any = None) -> Any

Funnel plot of SNP-specific Wald ratios vs. precision.

An asymmetric funnel around the IVW estimate suggests directional horizontal pleiotropy. Returns the matplotlib axis.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> bx = rng.uniform(0.2, 0.5, 25)
>>> by = 0.3 * bx + rng.normal(0, 0.01, 25)
>>> sy = rng.uniform(0.02, 0.04, 25)
>>> ax = sp.mr_funnel_plot(bx, by, sy)
>>> ax is not None
True

mr_scatter_plot ¶

mr_scatter_plot(beta_exposure: ndarray, beta_outcome: ndarray, se_exposure: ndarray, se_outcome: ndarray, *, ax: Any = None) -> Any

Classic MR scatter plot with IVW and MR-Egger lines.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> bx = rng.uniform(0.2, 0.5, 25)
>>> by = 0.3 * bx + rng.normal(0, 0.01, 25)
>>> sx = rng.uniform(0.02, 0.04, 25)
>>> sy = rng.uniform(0.02, 0.04, 25)
>>> ax = sp.mr_scatter_plot(bx, by, sx, sy)
>>> ax is not None
True

mr_multivariable ¶

mr_multivariable(snp_associations: DataFrame, *, outcome: str = 'beta_y', outcome_se: str = 'se_y', exposures: Optional[Sequence[str]] = None, alpha: float = 0.05, exposure_se: Optional[Sequence[str]] = None) -> MVMRResult

Multivariable Mendelian Randomization (Sanderson et al. 2019).

Parameters:

Name	Type	Description	Default
`snp_associations`	`DataFrame`	One row per SNP. Must contain: `outcome` — outcome-SNP association (β_Y). `outcome_se` — its standard error. exposure columns `beta_X1, beta_X2, ...` (configurable via the `exposures` argument).	required
`outcome`	`str`		`'beta_y'`
`outcome_se`	`str`		`'beta_y'`
`exposures`	`sequence of str`	Exposure-beta column names. If omitted, all `beta_*` columns other than the outcome are used.	`None`
`alpha`	`float`		`0.05`
`exposure_se`	`sequence of str`	Columns of the variant-exposure standard errors, in the order of `exposures`. Needed for the conditional F-statistics; without them those are NaN.	`None`

Returns:

Type Description

MVMRResult

Contains direct-effect estimates for each exposure and the two-sample conditional F-statistics of Sanderson, Spiller & Bowden (2021): for exposure k, Q_k = sum_j (b_kj - d' b_-k,j)^2 / s2_kj with d the least-squares regression (no intercept) of b_k on the other exposures' associations, s2_kj = se_kj^2 + sum_m d_m^2 se_mj^2 (the exposure associations taken as uncorrelated, as from separate samples), and F_k = Q_k / (L - K + 1) for L variants and K exposures: R MVMR::strength_mvmr(gencov = 0). Values below 10 flag weak instruments conditional on the other exposures.

.. versionchanged:: 1.41 conditional_f_stats used to be the F of regressing b_k on the other exposures (with an intercept and outcome weights), which is large, or infinite, exactly when the exposures are collinear and the conditional F should be small. It is now the statistic above, and NaN without exposure_se.

Notes

Model:

β_Y = Σ_j α_j β_{X_j} + ε,  ε ~ Normal(0, σ_Y²)

fitted by weighted least squares with weights 1 / se_y². The α_j are the direct causal effects, holding other exposures fixed.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n_snps = 30
>>> bx1 = rng.uniform(0.1, 0.5, n_snps)
>>> bx2 = rng.uniform(0.1, 0.5, n_snps)
>>> # True direct effects: X1 -> Y = 0.4, X2 -> Y = 0.0.
>>> by = 0.4 * bx1 + 0.0 * bx2 + rng.normal(0, 0.02, n_snps)
>>> snp = pd.DataFrame({
...     "beta_x1": bx1, "beta_x2": bx2, "beta_y": by,
...     "se_y": rng.uniform(0.01, 0.05, n_snps),
... })
>>> res = sp.mr_multivariable(snp, outcome="beta_y", outcome_se="se_y")
>>> isinstance(res, sp.MVMRResult)
True
>>> res.exposures
['beta_x1', 'beta_x2']

References

Sanderson, E., Davey Smith, G., Windmeijer, F. & Bowden, J. (2019). "An examination of multivariable Mendelian randomization in the single-sample and two-sample summary data settings." IJE 48(3). [@sanderson2019examination]

mr_mediation ¶

mr_mediation(snp_associations: DataFrame, *, beta_exposure: str = 'beta_x', se_exposure: str = 'se_x', beta_mediator: str = 'beta_m', se_mediator: str = 'se_m', beta_outcome: str = 'beta_y', se_outcome: str = 'se_y', exposure_name: str = 'exposure', mediator_name: str = 'mediator', outcome_name: str = 'outcome') -> MediationMRResult

Two-step MR mediation: decompose total effect into direct + indirect.

Parameters:

Name	Type	Description	Default
`snp_associations`	`DataFrame`	One row per SNP. Must contain SNP-associations with the exposure, mediator and outcome (both beta and SE).	required
`beta_exposure`	`str`		`'beta_x'`
`se_exposure`	`str`		`'beta_x'`
`beta_mediator`	`str`		`'beta_x'`
`se_mediator`	`str`		`'beta_x'`
`beta_outcome`	`str`		`'beta_y'`
`se_outcome`	`str`		`'beta_y'`

Returns:

Type	Description
`MediationMRResult`

Notes

Two-step approach (Burgess et al. 2015):

Step 1 — total effect (IVW): α_total = IVW(β_Y ~ β_X) (sp.mr_ivw's default model: multiplicative random effects).
Step 2a — exposure → mediator: α_XM = IVW(β_M ~ β_X).
Step 2b — direct effect via MVMR: α_direct = MVMR(β_Y ~ β_X, β_M).
Indirect: α_indirect = α_total − α_direct.

The direct effect is sp.mr_multivariable (mr_mvivw). The SE of the indirect effect combines the SEs of steps 1 and 3 as if they were independent (they are not; treat it as approximate).

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> n_snps = 30
>>> bx = rng.uniform(0.1, 0.5, n_snps)
>>> # exposure -> mediator (0.5), mediator -> outcome (0.6),
>>> # direct (0.2).
>>> bm = 0.5 * bx + rng.normal(0, 0.02, n_snps)
>>> by = 0.2 * bx + 0.6 * bm + rng.normal(0, 0.02, n_snps)
>>> snp = pd.DataFrame({
...     "beta_x": bx, "se_x": rng.uniform(0.01, 0.05, n_snps),
...     "beta_m": bm, "se_m": rng.uniform(0.01, 0.05, n_snps),
...     "beta_y": by, "se_y": rng.uniform(0.01, 0.05, n_snps),
... })
>>> res = sp.mr_mediation(snp)
>>> isinstance(res, sp.MediationMRResult)
True
>>> bool(res.indirect_effect > 0)  # positive mediated path X -> M -> Y
True

References

Burgess, S., Daniel, R. M., Butterworth, A. S. & Thompson, S. G. (2015). "Network Mendelian randomization: using genetic variants as instrumental variables to investigate mediation in causal pathways." IJE 44(2). [@burgess2015network]

mr_bma ¶

mr_bma(snp_associations: DataFrame, *, outcome: str = 'beta_y', outcome_se: str = 'se_y', exposures: Optional[Sequence[str]] = None, max_model_size: Optional[int] = None, prior_inclusion: float = 0.5, prior_sd: float = 0.5, method: str = 'bf') -> MRBMAResult

Mendelian Randomization with Bayesian Model Averaging.

Parameters:

Name	Type	Description	Default
`snp_associations`	`DataFrame`		required
`outcome`	see :func:`mr_multivariable`		`'beta_y'`
`outcome_se`	see :func:`mr_multivariable`		`'beta_y'`
`exposures`	see :func:`mr_multivariable`		`'beta_y'`
`max_model_size`	`int`	Maximum number of exposures in any one model (Zuber et al.'s `kmax`), an integer from 1 to `len(exposures)`. Defaults to `len(exposures)`. The search is exhaustive over the models of size up to this cap and is refused beyond 16,383 models (all subsets of 14 exposures); with more exposures, cap the model size.	`None`
`prior_inclusion`	`float`	Prior probability that each exposure is in the true causal model (`prior_prob`); a model of size `s` has prior `p^s (1 - p)^(k - s)`.	`0.5`
`prior_sd`	`float`	Standard deviation `sigma` of the independent normal prior on each causal effect (Zuber et al.'s `sigma`).	`0.5`
`method`	`('bf', 'bic')`	`'bf'` is MR-BMA as Zuber, Colijn, Klaver & Burgess (2020) define and implement it (their `summary_mvMR_BF`): the IVW-scaled regression `beta_Y / se_Y ~ beta_X / se_Y` without intercept, a closed-form Bayes factor under the normal prior, posterior-mean effects averaged over models. `'bic'` is the pre-1.30 behaviour: posterior weights `exp(-BIC/2)` from weighted least squares.	`'bf'`

Returns:

Type	Description
`MRBMAResult`	`marginal_inclusion` (posterior inclusion probabilities), `model_averaged_estimate` (BMA causal effects), `best_models` (top 20; `log10_bf` or `bic`), `model_priors` (posterior probability of every model, in enumeration order).

Notes

With X = beta_X / se_Y, Y = beta_Y / se_Y and a model g, Omega^-1 = sigma^-2 I + X_g'X_g, B = Omega X_g'Y and

log10 BF_g = -0.5 log10 det(Omega^-1) - |g| log10 sigma
             - (n/2) [log10(Y'Y - B' Omega^-1 B) - log10(Y'Y)].

Before 1.30 this function returned BIC-weighted BMA under the MR-BMA name (method='bic'); its inclusion probabilities differ from Zuber et al.'s and it reported no averaged effect.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(2)
>>> n_snps = 40
>>> bx1 = rng.uniform(0.1, 0.5, n_snps)
>>> bx2 = rng.uniform(0.1, 0.5, n_snps)
>>> bx3 = rng.uniform(0.1, 0.5, n_snps)
>>> by = 0.4 * bx1 + rng.normal(0, 0.02, n_snps)  # only X1 is causal
>>> snp = pd.DataFrame({
...     "beta_x1": bx1, "beta_x2": bx2, "beta_x3": bx3, "beta_y": by,
...     "se_y": rng.uniform(0.01, 0.05, n_snps),
... })
>>> res = sp.mr_bma(snp, outcome="beta_y", outcome_se="se_y")
>>> isinstance(res, sp.MRBMAResult)
True
>>> # X1 carries the most posterior inclusion mass.
>>> bool(
...     res.marginal_inclusion["beta_x1"]
...     > res.marginal_inclusion["beta_x2"]
... )
True

References

Zuber, V., Colijn, J. M., Klaver, C. & Burgess, S. (2020). "Selecting likely causal risk factors from high-throughput experiments using multivariable Mendelian randomization." Nature Communications 11, 29. [@zuber2020selecting]

mr_lap ¶

mr_lap(beta_exposure: ndarray, beta_outcome: ndarray, se_exposure: ndarray, se_outcome: ndarray, *, overlap_fraction: float = 1.0, overlap_rho: float = 0.0, alpha: float = 0.05) -> MRLapResult

Sample-overlap-corrected inverse-variance-weighted MR.

Corrects the first-order bias introduced when the exposure and outcome GWAS share participants. Under Burgess, Davies & Thompson (2016) Eq. 8, the bias in two-sample IVW satisfies

.. math::

E[\hat\beta_{IVW}] - \beta \approx \frac{p \, \rho_{obs}}{F_{mean}}

where :math:p is the overlap fraction (share of samples appearing in both GWAS), :math:\rho_{obs} is the phenotypic correlation of exposure and outcome in the overlapping sub-sample, and :math:F_{mean} is the mean first-stage F-statistic of the SNPs.

When overlap_fraction=0 (clean two-sample design) the correction vanishes and the estimate equals standard IVW.

Parameters:

Name	Type	Description	Default
`beta_exposure`	`ndarray`	Per-SNP GWAS effect sizes.	required
`beta_outcome`	`ndarray`	Per-SNP GWAS effect sizes.	required
`se_exposure`	`ndarray`	Per-SNP GWAS standard errors.	required
`se_outcome`	`ndarray`	Per-SNP GWAS standard errors.	required
`overlap_fraction`	`float`	Share of the outcome sample that also appears in the exposure sample. 0 = clean two-sample; 1 = full overlap; in between = partial overlap. Must be in [0, 1].	`1.0`
`overlap_rho`	`float`	Phenotypic correlation between exposure and outcome measurements in the overlapping sub-sample. In practice estimated by LD-score regression (`ldsc --rg`). Must be in [-1, 1].	`0.0`
`alpha`	`float`	Two-sided significance level for the CI and p-value.	`0.05`

Returns:

Type	Description
class:`MRLapResult`

Notes

This implements the closed-form Burgess-Davies-Thompson (2016) correction. The more elaborate MR-Lap of Mounier & Kutalik (2023) additionally regresses bias on a non-linear function of (F_mean, overlap_fraction) to absorb weak-instrument curvature; that refinement is deferred to a future release. When the user- supplied overlap_rho is obtained from LD-score regression this simplified correction captures most of the overlap bias on real-world GWAS scales (F > 30).

The SE is kept equal to the naive-IVW SE: the correction shifts the point estimate only. Monte-Carlo studies (Burgess 2016 Figs. 2-3) show coverage remains nominal under this conservative choice because the correction itself has second-order variance.

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> rng = np.random.default_rng(8)
>>> bx = rng.uniform(0.1, 0.5, 25)
>>> sx = rng.uniform(0.02, 0.05, 25)
>>> sy = rng.uniform(0.02, 0.05, 25)
>>> by = 0.3 * bx + rng.normal(0, sy)
>>> res = sp.mr_lap(bx, by, sx, sy,
...                 overlap_fraction=0.4, overlap_rho=0.18)
>>> type(res).__name__
'MRLapResult'
>>> bool(np.isfinite(res.estimate))
True

mr_clust ¶

mr_clust(beta_exposure: ndarray, beta_outcome: ndarray, se_exposure: ndarray, se_outcome: ndarray, *, K_range: Tuple[int, int] = (1, 5), include_null: bool = True, alpha: float = 0.05, seed: int = 0, max_iter: int = 200, tol: float = 1e-06) -> MRClustResult

Clustered Mendelian randomization via Gaussian mixture on Wald ratios.

Implements the MR-Clust estimator of Foley, Mason, Kirk & Burgess (2021) Bioinformatics 37(4). Assigns each SNP to one of K clusters based on its Wald ratio, where cluster 0 is a 'null' component (SNPs with zero effect on the outcome). Cluster means represent distinct causal-effect estimates — multiple non-null clusters are evidence of heterogeneous / pathway-specific pleiotropy.

Model

For SNP i in cluster k:

.. math::

r_i = \frac{\beta_{y,i}}{\beta_{x,i}}, \qquad r_i \mid z_i = k \ \sim\ \mathcal{N}(\theta_k, \tau_i^2)

with :math:\tau_i = s_{y,i} / |\beta_{x,i}| the first-order Wald SE. When include_null=True cluster 0 is fixed at :math:\theta_0 = 0. Parameters fitted by EM; K selected by BIC.

Parameters:

Name	Type	Description	Default
`beta_exposure`	`ndarray`		required
`beta_outcome`	`ndarray`		required
`se_exposure`	`ndarray`		required
`se_outcome`	`ndarray`		required
`K_range`	`(int, int)`	Inclusive range of cluster counts to try. K=1 with `include_null=True` is "everything null" and K=1 without is a single cluster at median ratio.	`(1, 5)`
`include_null`	`bool`	If True the first cluster is fixed at 0 (null).	`True`
`alpha`	`float`	For per-cluster CIs.	`0.05`
`seed`	`int`	EM initialisation rng.	`0`
`max_iter`	`EM controls.`		`200`
`tol`	`EM controls.`		`200`

Returns:

Type	Description
class:`MRClustResult`

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> rng = np.random.default_rng(8)
>>> bx = rng.uniform(0.1, 0.5, 25)
>>> sx = rng.uniform(0.02, 0.05, 25)
>>> sy = rng.uniform(0.02, 0.05, 25)
>>> by = 0.3 * bx + rng.normal(0, sy)
>>> res = sp.mr_clust(bx, by, sx, sy, K_range=(1, 4))
>>> type(res).__name__
'MRClustResult'
>>> res.n_snps
25

mr_cml ¶

mr_cml(beta_exposure: ndarray, beta_outcome: ndarray, se_exposure: ndarray, se_outcome: ndarray, *, K_max: Optional[int] = None, alpha: float = 0.05, max_iter: int = 100, tol: float = 1e-07, n: Optional[int] = None, model_average: bool = False) -> MRcMLResult

Constrained maximum-likelihood MR (cML-BIC / cML-MA-BIC).

Fits the constrained MLE of Xue, Shen & Pan (2021) for every number of invalid instruments K = 0, ..., K_max, selects K by BIC(K) = -2 L(K) + K log n and reports either the selected fit (cML-BIC) or the BIC-weighted average over K (cML-MA-BIC, model_average=True). Robust to both correlated and uncorrelated pleiotropy under the plurality-valid assumption.

Parameters:

Name	Type	Description	Default
`beta_exposure`	`array - like`	Per-variant associations with the exposure and the outcome.	required
`beta_outcome`	`array - like`	Per-variant associations with the exposure and the outcome.	required
`se_exposure`	`array - like`	Their standard errors (strictly positive).	required
`se_outcome`	`array - like`	Their standard errors (strictly positive).	required
`K_max`	`int`	Largest number of invalid instruments considered. Default `n_snps - 2` (the paper's candidate set `{0, ..., m - 2}`, which keeps at least two valid instruments); must lie in `[0, n_snps - 1]`.	`None`
`alpha`	`float`	Two-sided level for the normal confidence interval.	`0.05`
`max_iter`	`int`	Maximum coordinate-descent rounds per K.	`100`
`tol`	`float`	Stop when the change in theta between rounds is below `tol`.	`1e-7`
`n`	`int`	GWAS sample size for the BIC penalty; the paper recommends `min(n_exposure, n_outcome)`. When omitted the number of variants is used and a warning is issued.	`None`
`model_average`	`bool`	Report the BIC model average (Equation 5) instead of the BIC-selected fit.	`False`

Returns:

Type	Description
`MRcMLResult`

Notes

Variants with negative exposure associations are sign-flipped first (the likelihood is invariant to flipping both associations of a variant). The data-perturbation variants (cML-*-DP) and random restarts are not implemented: the fit uses the paper's default start theta = 0, b = 0.

Clean-room implementation; the reference (MendelianRandomization 0.10.0, GPL-2 | GPL-3) was used only as a black box.

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> bx = rng.uniform(0.1, 0.5, 30)
>>> sx = np.full(30, 0.02)
>>> sy = np.full(30, 0.02)
>>> by = 0.4 * bx + rng.normal(0, sy)
>>> r = sp.mr_cml(bx, by, sx, sy, n=50_000)
>>> r.K_selected <= 2
True

References

xue2021constrained

mr_raps ¶

mr_raps(beta_exposure: ndarray, beta_outcome: ndarray, se_exposure: ndarray, se_outcome: ndarray, *, loss: Optional[str] = None, over_dispersion: bool = True, tuning_c: Optional[float] = None, alpha: float = 0.05, pruning: bool = True, niter: int = 20, tol: float = _EPS_HALF, beta_init: Optional[float] = None, tau2_init: Optional[float] = None) -> MRRapsResult

Robust Adjusted Profile Score MR (Zhao et al. 2020).

Profile-likelihood MR that accounts for measurement error in the exposure associations (weak instruments), optionally for systematic pleiotropy (a random-effect variance tau2) and, with a robust loss, for a small fraction of gross pleiotropic outliers.

Parameters:

Name	Type	Description	Default
`beta_exposure`	`array - like`	Variant-exposure and variant-outcome associations.	required
`beta_outcome`	`array - like`	Variant-exposure and variant-outcome associations.	required
`se_exposure`	`array - like`	Their standard errors (strictly positive).	required
`se_outcome`	`array - like`	Their standard errors (strictly positive).	required
`loss`	`('huber', 'tukey', 'l2')`	Loss in the profile score. `None` means `'huber'`.	`'huber'`
`over_dispersion`	`bool`	Estimate the over-dispersion variance `tau2` (Sections 4-5). `False` is the profile-score estimator of Section 3 and needs `loss='l2'`.	`True`
`tuning_c`	`float`	Loss tuning constant `k`; default 1.345 (Huber) or 4.685 (Tukey). Ignored for `loss='l2'`.	`None`
`alpha`	`float`	`1 - alpha` confidence level.	`0.05`
`pruning`	`bool`	Drop variants with `se_exposure > 10 * median(se_exposure)`.	`True`
`niter`	`int`	Maximum number of alternating `tau2` / `beta` iterations.	`20`
`tol`	`float`	Convergence threshold on the change in `(beta, tau2)`.	``sqrt(machine eps)``
`beta_init`	`float`	Deprecated and ignored: the over-dispersed fit starts from the no-over-dispersion estimate and the robust fit from the L2 fit.	`None`
`tau2_init`	`float`	Deprecated and ignored: the over-dispersed fit starts from the no-over-dispersion estimate and the robust fit from the L2 fit.	`None`

Returns:

Type	Description
`MRRapsResult`

Raises:

Type	Description
`ValueError`	On mismatched lengths, fewer than 3 variants, non-positive standard errors, an unknown `loss`, a non-positive `tuning_c`, or a robust loss without over-dispersion.

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> rng = np.random.default_rng(1)
>>> bx = rng.uniform(0.1, 0.4, 30)
>>> sx = np.full(30, 0.02)
>>> sy = np.full(30, 0.02)
>>> by = 0.5 * bx + rng.normal(0, 0.02, 30)
>>> r = sp.mr_raps(bx, by, sx, sy, loss="tukey")
>>> bool(abs(r.estimate - 0.5) < 0.2)
True

Notes

Clean-room implementation of the paper's estimating equations; the reference (R mr.raps 0.4.3, GPL-3) was used only as a black box. The (R)APS root is solved to machine precision, so it can differ from the reference in the digits that the reference's own stopping rule leaves undetermined (see docs/dev/cleanroom/mr_raps.md).

References

zhao2020statistical

mr_available_methods ¶

mr_available_methods() -> list[str]

Return the full list of registered MR method names (incl. aliases).

Examples:

>>> import statspai as sp
>>> methods = sp.mr_available_methods()
>>> isinstance(methods, list)
True
>>> "ivw" in methods and "mvmr" in methods
True
>>> methods == sorted(methods)  # returned sorted
True

statspai.mendelian¶

mendelian ¶

MRResult ¶

plot ¶

HeterogeneityResult dataclass ¶

PleiotropyResult dataclass ¶

LeaveOneOutResult dataclass ¶

SteigerResult dataclass ¶

RadialResult dataclass ¶

ModeBasedResult dataclass ¶

FStatisticResult dataclass ¶

MVMRResult dataclass ¶

MediationMRResult dataclass ¶

MRBMAResult dataclass ¶

MRLapResult dataclass ¶

MRClustResult dataclass ¶

GrappleResult dataclass ¶

MRcMLResult dataclass ¶

MRRapsResult dataclass ¶

mendelian_randomization ¶

mr_egger ¶

mr_ivw ¶

mr_median ¶

mr_plot ¶

mr_heterogeneity ¶

mr_pleiotropy_egger ¶

mr_leave_one_out ¶

mr_steiger ¶

mr_radial ¶

mr_mode ¶

mr_f_statistic ¶

mr_funnel_plot ¶

mr_scatter_plot ¶

mr_multivariable ¶

mr_mediation ¶

mr_bma ¶

mr_lap ¶

mr_clust ¶

mr_cml ¶

mr_raps ¶

mr_available_methods ¶

`statspai.mendelian`¶

HeterogeneityResult `dataclass` ¶

PleiotropyResult `dataclass` ¶

LeaveOneOutResult `dataclass` ¶

SteigerResult `dataclass` ¶

RadialResult `dataclass` ¶

ModeBasedResult `dataclass` ¶

FStatisticResult `dataclass` ¶

MVMRResult `dataclass` ¶

MediationMRResult `dataclass` ¶

MRBMAResult `dataclass` ¶

MRLapResult `dataclass` ¶

MRClustResult `dataclass` ¶

GrappleResult `dataclass` ¶

MRcMLResult `dataclass` ¶

MRRapsResult `dataclass` ¶