Batched cell-type-specific stratified LD score regression.
Reproduces ldsc's cell_type_specific for many models at once. ldsc fits each model with its Hsq class (Hsq = "h-squared", ldsc's heritability estimator in gwaslab/extension/ldsc/ldsc_regressions.py, which fits LD score regression with a block-jackknife standard error); cell_type_specific calls it with old_weights=True, meaning a single weighted least squares with no iterative reweighting and no two-step intercept estimate. Each model regresses chi-square on [target, controls, baseline] LD scores scaled by N / mean(N), plus an intercept, and reports the target coefficient with a 200-block jackknife standard error.
Models sharing a control set share the design apart from the target column, so they are processed together as a control group. In exact mode each model keeps its own regression weights (the weights depend on the model's total LD and total M_5_50, which include the target); the per-block normal equations of a chunk of models are formed with a few large matrix products rather than one regression at a time. ldsc normalizes the square-root weights to sum to one before weighting the design and response; that scalar cancels in every solve, so the effective weight here is the raw ldsc weight.
In shared_controls mode every model in a control group uses the weight of the model without a target column, so the shared Gram matrix is formed once and only the target's border changes; the target coefficient follows from the Frisch-Waugh-Lovell identity.
Dimension names used throughout this module: S: number of regression SNPs (rows of every per-SNP array). B: number of baseline LD-score columns. U: number of LD-score columns in the annotation set (targets and controls). K: number of models (one reported target column each). k: number of models in the chunk currently being processed (k <= chunk_size). q: number of shared design columns of a control group: its controls, the baseline columns, and the intercept (q = number of controls + B + 1). p: number of design columns of one model: its target plus the shared columns (p = q + 1). n_blocks: number of jackknife blocks, min(S, 200).
M_5_50: for each annotation, the sum of its values over reference-panel SNPs with minor allele frequency above 5% (for a binary annotation, the number of common SNPs it contains). In the LDSC model E[chi2_j] = N_j * sum_c tau_c * l(j, c) + N_j * a + 1, where l(j, c) is the LD score of SNP j with annotation c and tau_c is the per-SNP heritability contributed by annotation c, M_5_50 converts per-SNP effects into the heritability attributed to an annotation (h2_c = tau_c * M_5_50_c). Here the total M_5_50 over a model's annotations scales the aggregate heritability estimate that shapes the regression weights.
Classes:
-
BatchedSLDSCResult–Per-model target coefficient, jackknife SE, and one-sided p-value, each (K,), in input model
-
SLDSCModel–target_index and control_indices are columns of SLDSCRegressionInputs.annotation_ld.
-
SLDSCRegressionInputs–Regression SNPs in ldsc row order (the jackknife blocks depend on it). Dimensions S, B, U
Functions:
-
batched_sldsc_cts– -
ldsc_block_separators–ldsc's Jackknife.get_separators: floor(linspace(0, n_snps, n_blocks + 1)).
-
ldsc_hsq_weights–Regression weights of ldsc's Hsq estimator (old_weights=True path), for k models at once.
Attributes:
Per-model target coefficient, jackknife SE, and one-sided p-value, each (K,), in input model order (K = number of models).
Attributes:
target_index and control_indices are columns of SLDSCRegressionInputs.annotation_ld.
Attributes:
-
control_indices(tuple[int, ...]) – -
label(str) – -
target_index(int) –
Regression SNPs in ldsc row order (the jackknife blocks depend on it). Dimensions S, B, U are defined in the module docstring.
chi2, n, w_ld: (S,) chi-square, integer-valued sample size, weight LD score. baseline_ld: (S, B) baseline LD scores; baseline_m: (B,) baseline M_5_50. annotation_ld: (S, U) target and control LD scores (float32 allowed); annotation_m: (U,) their M_5_50.
Attributes:
-
annotation_ld(ndarray) – -
annotation_m(ndarray) – -
baseline_ld(ndarray) – -
baseline_m(ndarray) – -
chi2(ndarray) – -
n(ndarray) – -
w_ld(ndarray) –
batched_sldsc_cts(
inputs: SLDSCRegressionInputs,
models: Sequence[SLDSCModel],
*,
weighting: SLDSCWeighting,
n_blocks: int = LDSC_DEFAULT_N_BLOCKS,
chunk_size: int = DEFAULT_CHUNK_SIZE,
) -> BatchedSLDSCResult
Source code in mecfs_bio/build_system/task/sldsc/batched_sldsc_cts.py
ldsc's Jackknife.get_separators: floor(linspace(0, n_snps, n_blocks + 1)).
ldsc_hsq_weights(
ld_tot: ndarray,
w_ld: ndarray,
n: ndarray,
m_tot: ndarray,
chi2: ndarray,
) -> np.ndarray
Regression weights of ldsc's Hsq estimator (old_weights=True path), for k models at once.
ld_tot: (S, k) unscaled total LD score of each model (sum over its annotations). w_ld: (S,) weight LD score; n: (S,) sample size; chi2: (S,) chi-square. m_tot: (k,) total M_5_50 of each model.
Returns (S, k) weights het * oc, where het = 1 / (2 (1 + h2_agg N max(ld_tot, 1) / m_tot)^2) corrects for heteroskedasticity, oc = 1 / max(w_ld, 1) for over-counting of correlated SNPs, and h2_agg = clip(m_tot (mean chi2 - 1) / mean(ld_tot N), 0, 1). ldsc weights the design and response by the square root of this, so it is the effective least-squares weight.