Skip to content

mecfs_bio.build_system.task.gwaslab.ldsc_diagnostic

The arithmetic behind the LD-score-regression diagnostic plot.

LD score regression models the per-variant chi-square as a straight line in the LD score: E[chi^2_j] = intercept + (N * h2 / M) * ld_j, where N is the sample size, M the number of reference variants, and h2 the observed-scale heritability. The diagnostic plot bins variants by LD score and shows each bin's mean chi-square against that line, so a bin that sits well off the line, or an intercept that drifts away from one, is visible by eye.

This module holds only the pure math: binning variants into equal-count LD-score bins, the fitted line, and the rescaling that turns the chi-square axis into heritability-slope units. It takes plain numpy arrays and returns plain values, so it can be tested without reading any asset or drawing anything.

Classes:

  • LdscDiagnosticBins –

    Per-bin summaries of variants grouped into equal-count LD-score bins, ordered by

Functions:

  • bin_by_ld_score –

    Group variants into n_bins equal-count bins by ascending LD score and summarize each.

  • chi2_to_heritability_units –

    Rescale chi-square to chi^2 * M / N, the units in which the fitted line's slope equals

  • fit_line_chi2 –

    The LDSC fitted line in chi-square units: intercept + (N * h2 / M) * ld.

LdscDiagnosticBins

Per-bin summaries of variants grouped into equal-count LD-score bins, ordered by ascending LD score. Every array is 1-D of length n_bins.

Attributes:

count instance-attribute 

count: ndarray

mean_chi2 instance-attribute 

mean_chi2: ndarray

mean_ld instance-attribute 

mean_ld: ndarray

se_chi2 instance-attribute 

se_chi2: ndarray

bin_by_ld_score 

bin_by_ld_score(
    chi2: ndarray, ld_score: ndarray, n_bins: int
) -> LdscDiagnosticBins

Group variants into n_bins equal-count bins by ascending LD score and summarize each.

Equal-count (rather than equal-width) bins put a comparable number of variants in every point, so each mean carries similar noise -- matching the binned plots in the original LDSC papers. Bin sizes differ by at most one when n_bins does not divide the variant count.

Source code in mecfs_bio/build_system/task/gwaslab/ldsc_diagnostic.py 
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
def bin_by_ld_score(
    chi2: np.ndarray, ld_score: np.ndarray, n_bins: int
) -> LdscDiagnosticBins:
    """Group variants into n_bins equal-count bins by ascending LD score and summarize each.

    Equal-count (rather than equal-width) bins put a comparable number of variants in every
    point, so each mean carries similar noise -- matching the binned plots in the original LDSC
    papers. Bin sizes differ by at most one when n_bins does not divide the variant count."""
    assert chi2.shape == ld_score.shape, (
        f"chi2 and ld_score must have the same shape, got {chi2.shape} and {ld_score.shape}"
    )
    assert n_bins >= 1, f"n_bins must be at least 1, got {n_bins}"

    order = np.argsort(ld_score, kind="stable")
    ld_sorted = ld_score[order]
    chi2_sorted = chi2[order]

    mean_ld = []
    mean_chi2 = []
    se_chi2 = []
    count = []
    for ld_bin, chi2_bin in zip(
        np.array_split(ld_sorted, n_bins), np.array_split(chi2_sorted, n_bins)
    ):
        c = chi2_bin.shape[0]
        mean_ld.append(float(ld_bin.mean()))
        mean_chi2.append(float(chi2_bin.mean()))
        # Sample std (ddof=1) needs at least two values; a singleton bin has no spread.
        se_chi2.append(float(chi2_bin.std(ddof=1) / np.sqrt(c)) if c > 1 else 0.0)
        count.append(c)

    return LdscDiagnosticBins(
        mean_ld=np.array(mean_ld, dtype=float),
        mean_chi2=np.array(mean_chi2, dtype=float),
        se_chi2=np.array(se_chi2, dtype=float),
        count=np.array(count, dtype=np.int64),
    )

chi2_to_heritability_units 

chi2_to_heritability_units(
    chi2: ndarray, n: float, m: float
) -> np.ndarray

Rescale chi-square to chi^2 * M / N, the units in which the fitted line's slope equals the observed-scale heritability. This is the secondary y-axis, and because M and N are divided out it is comparable across traits of different sample size.

Source code in mecfs_bio/build_system/task/gwaslab/ldsc_diagnostic.py 
86
87
88
89
90
def chi2_to_heritability_units(chi2: np.ndarray, n: float, m: float) -> np.ndarray:
    """Rescale chi-square to chi^2 * M / N, the units in which the fitted line's slope equals
    the observed-scale heritability. This is the secondary y-axis, and because M and N are
    divided out it is comparable across traits of different sample size."""
    return chi2 * m / n

fit_line_chi2 

fit_line_chi2(
    intercept: float,
    h2_obs: float,
    n: float,
    m: float,
    ld: ndarray,
) -> np.ndarray

The LDSC fitted line in chi-square units: intercept + (N * h2 / M) * ld.

Source code in mecfs_bio/build_system/task/gwaslab/ldsc_diagnostic.py 
93
94
95
96
97
98
def fit_line_chi2(
    intercept: float, h2_obs: float, n: float, m: float, ld: np.ndarray
) -> np.ndarray:
    """The LDSC fitted line in chi-square units: intercept + (N * h2 / M) * ld."""
    slope = n * h2_obs / m
    return intercept + slope * ld