Fitting Weighted Stratified LD Score Regression

Why Fit the Regression Carefully?

Stratified LD score regression connects GWAS association statistics to annotation-specific genetic contributions through a regression on LD scores. Constructing the annotation-specific LD scores supplies the predictors; fitting the regression estimates their coefficients from the observed GWAS statistics. The details matter: a mismatch in sample-size scaling changes the coefficient scale, and an inappropriate intercept choice can push baseline inflation into the annotation slopes.

The response is the observed chi-square statistic for each SNP. The predictors summarize how much annotated genetic variation that SNP tags through linkage disequilibrium. The regression describes the expected trend across SNPs; an individual chi-square value need not lie on that trend.

Definition: For SNP \(j\) and annotation \(k\), the annotation-specific LD score is \(\ell(j,k)=\sum_i r_{ji}^2 a_{ik}\), where \(r_{ji}\) is the genotype correlation between SNPs \(j\) and \(i\), and \(a_{ik}\) records SNP \(i\)'s value for annotation \(k\). In a GWAS with sample size \(N\), the stratified LD score regression model is

E[\chi_j^2]=a+N\sum_{k=1}^{K}\tau_k\ell(j,k).

Here \(a\) is the regression intercept, \(K\) is the number of annotations, and \(\tau_k\) is the coefficient for annotation \(k\). The intercept is often near 1 when there is no excess baseline inflation, but it may be estimated rather than fixed at 1.

Model and Fitting Choices

The regression uses the GWAS summary statistics and the annotation-specific LD scores calculated from an appropriate LD reference. Its basic assumptions are that the expected chi-square statistics follow the stated linear mean model, the LD scores represent the GWAS population and variant set adequately, and the sample-size scale is appropriate for the phenotype and study design. For a binary trait, this may require the effective sample size used in earlier material on case-control sampling, rather than simply the total number of participants.

The sample size multiplies every annotation term. Thus, when the same per-allele genetic architecture is studied in a larger sample, the expected genetic contribution to chi-square statistics grows with \(N\). The intercept is not multiplied by \(N\) in this parameterization: it represents the baseline level, including any excess inflation captured by the intercept. Keep that convention consistent throughout the fit.

The intercept can be handled in either of two ways. One may estimate \(a\) along with the annotation coefficients, allowing the fit to accommodate baseline inflation beyond 1. Alternatively, one may fix \(a=1\) if that constraint is justified by the analysis goal and assumptions. Fixing it when the true baseline differs from 1 can force the slopes to absorb that discrepancy. Estimating it does not guarantee that every source of bias has been removed; it gives the regression a baseline parameter.

In weighted least squares, SNP \(j\) receives a positive weight \(w_j\). The fit chooses coefficients to minimize the weighted sum of squared residuals between observed and fitted chi-square statistics. A high weight makes a SNP's residual contribute more to that objective; it does not mean that the SNP has a larger biological effect.

Regression Weights

Chi-square statistics are heteroskedastic: their variability generally increases with their expected value. A common working approximation for constructing regression weights is \(\operatorname{Var}(\chi_j^2)\approx 2\mu_j^2\), where \(\mu_j=E[\chi_j^2]\). Inverse-variance weighting then gives \(w_j\propto 1/(2\mu_j^2)\). Since the mean depends on unknown model parameters, practical weights are based on estimated or smoothed expected values, often using LD-score information.

These are working weights, not a claim that the chi-square statistics are independent or that the variance approximation is exact for every SNP. Linkage disequilibrium can make statistics at nearby SNPs correlated. The weighted point fit is still defined by the stated objective, but the diagonal weights alone do not represent all dependence when quantifying uncertainty.

Worked Example: How Sample Size Scales the Mean

Suppose a hypothetical GWAS has \(N=30{,}000\), two annotation-specific LD scores for one SNP of \(\ell(j,1)=40\) and \(\ell(j,2)=15\), coefficients \(\tau_1=2.0\times10^{-6}\) and \(\tau_2=1.0\times10^{-6}\), and fitted intercept \(a=1.05\). Use the model to calculate the expected chi-square statistic.

First calculate the annotation contribution before sample-size scaling:

\tau_1\ell(j,1)+\tau_2\ell(j,2) =(2.0\times10^{-6})(40)+(1.0\times10^{-6})(15) =80\times10^{-6}+15\times10^{-6} =95\times10^{-6}.

Now multiply by the sample size and add the intercept:

E[\chi_j^2]=1.05+(30{,}000)(95\times10^{-6}) =1.05+2.85=3.90.

The fitted model therefore predicts a mean chi-square statistic of 3.90 at this SNP. If the sample size were 15,000 while the coefficients and LD scores remained the same, the annotation contribution would be \(15{,}000(95\times10^{-6})=1.425\), giving a mean of \(1.05+1.425=2.475\). The intercept stays fixed in this model; the genetic contribution scales with \(N\).

Weighted Least-Squares Coefficients

For \(m\) SNPs, let \(\mathbf y\) be the vector of observed chi-square statistics. Form a design matrix \(\mathbf X\) whose first column is 1 for the intercept and whose remaining columns are \(N\ell(j,k)\), one for each annotation. Let \(\boldsymbol\beta=(a,\tau_1,\ldots,\tau_K)^\mathsf{T}\), and let \(\mathbf W\) be diagonal with entries \(w_j\). The weighted least-squares estimate minimizes \((\mathbf y-\mathbf X\boldsymbol\beta)^\mathsf{T}\mathbf W(\mathbf y-\mathbf X\boldsymbol\beta)\).

When \(\mathbf X^\mathsf{T}\mathbf W\mathbf X\) is invertible, the minimizer is

\widehat{\boldsymbol\beta} =(\mathbf X^\mathsf{T}\mathbf W\mathbf X)^{-1} \mathbf X^\mathsf{T}\mathbf W\mathbf y.

Invertibility is a necessary condition for this particular inverse formula. If all weights are strictly positive and \(\mathbf X\) has full column rank, it holds: for any nonzero vector \(\mathbf v\), \(\mathbf v^\mathsf{T}\mathbf X^\mathsf{T}\mathbf W\mathbf X\mathbf v=(\mathbf X\mathbf v)^\mathsf{T}\mathbf W(\mathbf X\mathbf v)>0\). If weights are allowed to be zero, full column rank of \(\mathbf X\) alone is not enough; the direct condition is that \(\mathbf W^{1/2}\mathbf X\) has full column rank.

A Short Derivation of the Estimate

The weighted objective is \(Q(\boldsymbol\beta)=(\mathbf y-\mathbf X\boldsymbol\beta)^\mathsf{T}\mathbf W(\mathbf y-\mathbf X\boldsymbol\beta)\). Differentiating with respect to \(\boldsymbol\beta\) and setting the derivative to zero gives the weighted normal equations:

\frac{\partial Q}{\partial\boldsymbol\beta} =-2\mathbf X^\mathsf{T}\mathbf W(\mathbf y-\mathbf X\boldsymbol\beta)=\mathbf 0, \qquad \mathbf X^\mathsf{T}\mathbf W\mathbf X\boldsymbol\beta =\mathbf X^\mathsf{T}\mathbf W\mathbf y.

When \(\mathbf X^\mathsf{T}\mathbf W\mathbf X\) is invertible, multiplying both sides by its inverse yields the estimate above. The objective is a convex quadratic because \(\mathbf W\) is nonnegative; under the stated invertibility condition it is strictly convex and the solution is unique.

Worked Example: Estimating a Slope and Intercept

Consider a small hypothetical one-annotation regression with three SNPs. To keep the arithmetic readable, define the scaled predictor \(q_j=N\ell(j)/1000\), so that the fitted model is \(y_j=a+\beta q_j\), where \(y_j=\chi_j^2\) and \(\beta=1000\tau\). Suppose the observations and positive illustrative weights are:

SNP	\(q_j\)	\(y_j\)	\(w_j\)
1	1	1.8	1
2	2	2.1	2
3	3	2.7	1

For a weighted regression with an intercept, the weighted means are \(\bar q_w=\sum_jw_jq_j/\sum_jw_j\) and \(\bar y_w=\sum_jw_jy_j/\sum_jw_j\). Here \(\sum_jw_j=4\), so

\bar q_w=\frac{1(1)+2(2)+1(3)}{4}=2, \qquad \bar y_w=\frac{1(1.8)+2(2.1)+1(2.7)}{4} =\frac{8.7}{4}=2.175.

The weighted slope is the weighted covariance numerator divided by the weighted predictor sum of squares:

\widehat\beta= \frac{\sum_j w_j(q_j-\bar q_w)(y_j-\bar y_w)} {\sum_j w_j(q_j-\bar q_w)^2} = \frac{1(-1)(-0.375)+2(0)(-0.075)+1(1)(0.525)} {1(1)+2(0)+1(1)} =\frac{0.900}{2}=0.450.

Thus \(\widehat a=\bar y_w-\widehat\beta\bar q_w=2.175-(0.450)(2)=1.275\). The fitted values are 1.725, 2.175, and 2.625; their residuals are 0.075, \(-0.075\), and 0.075. Because \(q=N\ell/1000\), the coefficient in the unscaled model is \(\widehat\tau=\widehat\beta/1000=0.000450\). These arbitrary positive weights illustrate the calculation; they are not presented as calibrated LDSC weights.

Worked Example: Constructing Relative Weights

Suppose, for illustration, that a hypothetical single-annotation model has \(N=10{,}000\), \(M=100{,}000\) variants, \(h_g^2=0.20\), and intercept 1. Use the rough mean approximation \(\mu_j=1+Nh_g^2\ell_j/M\) to compare the working inverse-variance weights for SNPs with LD scores 50 and 200.

The scaling factor is \(Nh_g^2/M=(10{,}000)(0.20)/100{,}000=0.020\). The estimated means are therefore \(\mu_1=1+(0.020)(50)=2\) and \(\mu_2=1+(0.020)(200)=5\). Using \(w_j=1/(2\mu_j^2)\), the weights are

w_1=\frac{1}{2(2^2)}=0.125, \qquad w_2=\frac{1}{2(5^2)}=0.020.

Their relative weights are \(w_1/w_2=0.125/0.020=6.25\). Under this working approximation, the higher-LD-score SNP receives less weight because its expected chi-square statistic, and hence its approximate variance, is larger. In an applied fit, weights are constructed from the specified regression-weight procedure rather than treating these illustrative values as universal constants.

Practical Checks and Common Pitfalls

Before interpreting a fit, verify that the chi-square statistics and LD scores refer to compatible variants and that the sample-size definition matches the analysis. Check whether the intercept was estimated or fixed and report that choice. Also confirm that the design matrix and weights permit a unique fit. Since the regression is a joint model, a coefficient is conditional on the other predictors included; its detailed interpretation can depend on how annotations overlap.

Confusing total sample size with the regression scale. The model multiplies annotation terms by \(N\); changing or omitting this factor changes the coefficient scale.
Fixing the intercept automatically at 1. If the baseline is elevated, fixing it can shift that discrepancy into the annotation slopes. State and justify the intercept choice.
Assuming higher weights mean stronger biological effects. Weights control contributions to the fitting objective; they are not effect sizes.
Using full column rank as the only invertibility check. The weighted design must have full column rank. With zero weights, the unweighted design can have full rank while the weighted matrix is singular.
Treating fitted values as observations or causal effects. The model concerns expected chi-square statistics and estimates an annotation-based association pattern; it does not establish that an annotation causes greater heritability.

Key takeaway: Fit the observed chi-square statistics against sample-size-scaled annotation LD scores using positive, justified weights. Make the intercept choice explicit, and ensure the weighted design identifies a unique coefficient vector before interpreting the regression.

Check Your Understanding

Use the model and fitting ideas above to answer the following questions.

In \(E[\chi_j^2]=a+N\sum_k\tau_k\ell(j,k)\), which part is multiplied by sample size, and what does the intercept represent?
If the same coefficient and LD score are used but \(N\) doubles, how does the modeled annotation contribution change?
Why can a high-LD-score SNP receive a smaller working regression weight?
In the weighted example, what would be the unscaled coefficient \(\tau\) if the predictor were \(q=N\ell/1000\) and the fitted slope were 0.45?
Why is full column rank of \(\mathbf X\) not sufficient for the inverse formula if some diagonal weights are zero?