Continuous-Shrinkage Scores with PRS-CS

Why Continuous Shrinkage?

A polygenic score combines many variant effects, but GWAS estimates can be noisy, especially when variants are correlated through linkage disequilibrium (LD). A strategy that discards every variant below a significance threshold can remove real but modest effects; leaving every marginal estimate unaltered can retain substantial noise. Bayesian shrinkage instead estimates effects using both the observed summary statistics and a prior that favors small effects while allowing some effects to be larger.

PRS-CS is a Bayesian method for this task. Its name refers to a continuous-shrinkage prior: unlike a prior that assigns a variant exactly zero effect with positive probability, the prior allows effects across a continuous range. A global scale controls the overall amount of shrinkage, while local scales let individual variants have different prior variances. An LD reference panel supplies the correlations needed to distinguish evidence about correlated variants.

Model and Notation

Let $m$ be the number of variants in a genomic region. Let $X$ be the $n$-by-$m$ genotype matrix, with each column centered and scaled to have variance one, and let $y$ be a centered, scaled quantitative phenotype. Write $\beta_j$ for the joint effect of variant $j$ in this model. The basic regression model is $y=X\beta+\epsilon$, with residuals modeled as independent, mean-zero normal errors of variance one. This scaling makes the effects and summary-statistic approximations below interpretable on standardized scales.

PRS-CS works with GWAS summary statistics rather than requiring individual-level genotypes from the discovery study. Denote the marginal GWAS effect estimate for variant $j$ by $\hat b_j$, and let $R$ be the LD correlation matrix for the variants, estimated from a reference panel. Under the standardized model and the usual large-sample summary-statistic approximation, the vector of marginal estimates has likelihood

\hat b \mid \beta,R \ \approx\ N\!\left(R\beta,\frac{1}{n}R\right).

Thus, a marginal estimate is not generally the joint effect: its expected value depends on the effects at correlated variants. The likelihood also represents dependence among the summary statistics. The reference panel is used to estimate $R$; it does not supply the phenotype effects. In practice, sample-size and phenotype scaling must be handled consistently with the GWAS summary statistics. The approximation presumes appropriately aligned alleles, an LD matrix on the same variants and allele coding, and a sufficiently suitable reference population. It is not a guarantee that those conditions hold in a particular analysis.

Definition: In the PRS-CS prior, each joint effect has a zero-centered normal distribution conditional on its local variance: $$ \beta_j\mid \phi,\psi_j \sim N(0,\phi\psi_j). $$ Here $\phi>0$ is the global shrinkage parameter and $\psi_j>0$ is the local variance for variant $j$. A commonly specified continuous-shrinkage construction uses the gamma-gamma hierarchy $\psi_j\mid\delta_j\sim\operatorname{Gamma}(a,\delta_j)$, $\delta_j\sim\operatorname{Gamma}(b,1)$, with gamma distributions parameterized by shape and rate and commonly fixed $a=b=1/2$. This mixture gives a continuous prior with substantial concentration near zero and heavier tails than a single normal prior.

The global scale $\phi$ adjusts the prior variance across the genome: holding local scales fixed, smaller $\phi$ makes effects a priori more concentrated near zero, and larger $\phi$ permits larger effects overall. The local scale $\psi_j$ changes the prior variance for variant $j$ relative to other variants. A small local variance gives that effect stronger shrinkage; a large local variance gives its likelihood evidence more room to support a non-negligible effect. The local scales are estimated through the hierarchical model rather than assigned from a hard significance cutoff.

Combining the Likelihood, LD, and Prior

For fixed $\phi$ and local variances $\psi=(\psi_1,\ldots,\psi_m)$, the posterior density for the effects is proportional to the summary-statistic likelihood multiplied by the normal prior densities. In compact form,

p(\beta\mid \hat b,R,\phi,\psi) \ \propto\ N\!\left(\hat b\mid R\beta,\frac{1}{n}R\right) \prod_{j=1}^{m}N(\beta_j\mid 0,\phi\psi_j).

The full PRS-CS model also accounts for uncertainty in the local scales under their gamma-gamma hierarchy; estimation therefore involves more than substituting a single fixed shrinkage factor into each GWAS coefficient. The posterior distribution combines evidence across correlated variants, with the LD matrix determining how the marginal signals relate to candidate joint effects. Posterior means of the effects are commonly used as score weights.

The formula is a model statement, not a validation check. It relies on the specified standardized regression and approximate summary-statistic likelihood, on consistent variant and effect-allele alignment, and on an LD reference matrix that represents the relevant genotype correlations. Violations can affect the resulting estimates even if the calculations under the model are performed correctly.

Worked Example: Global and Local Prior Scales

Suppose, purely hypothetically, that a PRS-CS model uses $\phi=0.01$ and considers three variants with local variances $\psi_1=0.5$, $\psi_2=4$, and $\psi_3=25$. These values are stipulated to illustrate the prior component; they are not estimates inferred from data in this calculation.

For each variant, the conditional prior variance is $\phi\psi_j$. Therefore,

\begin{aligned} \operatorname{Var}(\beta_1\mid\phi,\psi_1)&=0.01(0.5)=0.005,\\ \operatorname{Var}(\beta_2\mid\phi,\psi_2)&=0.01(4)=0.04,\\ \operatorname{Var}(\beta_3\mid\phi,\psi_3)&=0.01(25)=0.25. \end{aligned}

The corresponding conditional prior standard deviations are $\sqrt{0.005}=0.07071$, $\sqrt{0.04}=0.2000$, and $\sqrt{0.25}=0.5000$. With the same global scale, variant 1 is most tightly concentrated around zero, while variant 3 has the widest prior distribution. This comparison explains the roles of the scales; it does not establish that variant 3 has a larger true effect or that its local variance would be estimated as 25 in an actual analysis.

A Conditional Update and Its Derivation

A useful way to see how LD enters the calculation is to examine the conditional distribution for one effect, with the other effects and the local scales held fixed. Assume $R$ is positive definite for this derivation, $R_{jj}=1$, and use the likelihood and independent conditional normal priors stated above. The likelihood exponent, up to terms that do not involve $\beta$, is

-\frac{n}{2}\left(\hat b'R^{-1}\hat b-2\beta'\hat b+\beta'R\beta\right).

The prior contribution involving $\beta_j$ is $-\beta_j^2/(2\phi\psi_j)$. Collecting the likelihood terms that involve $\beta_j$, and writing $\beta_{-j}$ for all effects other than $\beta_j$, gives a quadratic with precision $nR_{jj}+1/(\phi\psi_j)$ and linear coefficient $n\{\hat b_j-\sum_{k\ne j}R_{jk}\beta_k\}$. Completing the square yields

\beta_j\mid\beta_{-j},\hat b,R,\phi,\psi \sim N\!\left( \frac{n\left(\hat b_j-\sum_{k\ne j}R_{jk}\beta_k\right)} {nR_{jj}+1/(\phi\psi_j)}, \frac{1}{nR_{jj}+1/(\phi\psi_j)} \right).

Because $R_{jj}=1$, the term $\sum_{k\ne j}R_{jk}\beta_k$ adjusts the marginal signal for the effects of correlated variants. The prior adds $1/(\phi\psi_j)$ to the likelihood precision, pulling the conditional mean toward zero. For $R$ with unit diagonal, a smaller $\phi\psi_j$ means a larger added precision and more shrinkage. This conditional distribution is one building block of posterior computation; it is not, by itself, the final posterior mean, since other effects and local scales are uncertain too.

Worked Example: One Conditional Effect Calculation

Consider a hypothetical two-variant region with $n=1{,}000$, marginal estimates $\hat b_1=0.040$ and $\hat b_2=0.030$, and LD matrix

R=\begin{pmatrix}1&0.60\\0.60&1\end{pmatrix}.

For illustration, assume the summary-statistic likelihood above is appropriate, the alleles and variant ordering are consistently aligned, and this matrix is the standardized-variant LD matrix to use. Those assumptions are stipulated for the calculation; the supplied numbers cannot check or validate them. Hold $\beta_2=0.020$ and $\phi\psi_1=0.040$ fixed. The conditional mean for $\beta_1$ is

\frac{1{,}000\left(0.040-0.60(0.020)\right)} {1{,}000+1/0.040} = \frac{1{,}000(0.028)}{1{,}025} = \frac{28}{1{,}025} \approx 0.02732.

The conditional variance is $1/1{,}025=0.0009756$, so the conditional standard deviation is approximately $\sqrt{0.0009756}=0.03123$. The LD adjustment changes the signal used for variant 1 from $0.040$ to $0.040-0.60(0.020)=0.028$, and the prior precision further pulls the conditional mean toward zero. This result is conditional on the stipulated value of $\beta_2$ and local variance; it should not be reported as the completed PRS-CS posterior estimate.

From Posterior Effects to a Score

After fitting the model, a common choice is the posterior mean $\tilde\beta_j=E(\beta_j\mid\text{summary data and model})$ as the weight for variant $j$. The score for an individual is then a weighted sum of aligned genotype dosages, using the effect-allele convention and genotype scaling chosen for the score. A PRS-CS weight is therefore a model-based estimate informed jointly by the GWAS summary statistics, LD, and the prior; it is not simply the unmodified marginal GWAS estimate.

Worked Example: Applying Hypothetical PRS-CS Weights

Suppose a completed fit for a two-variant hypothetical example reports posterior mean weights of $0.0270$ and $0.0180$, in standardized phenotype units per standardized genotype unit. These are stipulated outputs for demonstrating score arithmetic, not values derived from the previous conditional update. Assume alleles are harmonized and a person has standardized genotype values $x_1=1$ and $x_2=2$, measured on the same genotype scale used to fit the weights.

The score is the sum of genotype values multiplied by their corresponding weights:

S=x_1\tilde\beta_1+x_2\tilde\beta_2 =1(0.0270)+2(0.0180) =0.0270+0.0360 =0.0630.

For this hypothetical individual, the score is $0.0630$ standardized phenotype units on the stated scale. It is a model-based genetic prediction component, not a causal effect estimate for that person and not, on its own, a calibrated clinical risk.

Practical Interpretation and Common Pitfalls

The central advantage of continuous shrinkage is flexibility: a shared global scale regularizes the overall collection of effects, while local scales allow effects to depart to different degrees from that common level. Incorporating LD is essential because correlated variants can carry overlapping marginal signals. Consequently, the result depends not only on each variant's GWAS statistic but also on the LD structure used in the joint model.

Confusing a marginal estimate with a joint effect. The likelihood has expectation $R\beta$, not simply $\beta$. With correlated variants, the marginal estimate for one variant can reflect effects at other variants.
Calling continuous shrinkage a hard filter. The gamma-gamma prior does not create a point mass at zero or apply a p-value threshold. Shrinkage can be strong, but the model is not a rule that deletes every small signal.
Reversing the scale interpretation. Increasing $\phi\psi_j$ increases the conditional prior variance and generally relaxes prior shrinkage; decreasing it concentrates the prior more tightly around zero.
Treating the LD panel as a harmless input. The calculation assumes appropriate variant matching, allele coding, and LD representation. Correct arithmetic does not establish that these assumptions are satisfied.
Overinterpreting an intermediate update. A conditional mean with other effects fixed is not the posterior mean after integrating over uncertainty in all effects and local scales.

PRS-CS does not remove the need to inspect the inputs and modeling scale. Its estimates are meaningful only relative to the summary-statistic likelihood and the LD matrix supplied to the procedure. Keeping the distinction between an assumption, a computed conditional quantity, and a fitted posterior estimate explicit is especially important when interpreting numerical examples or reporting score weights.

Check Your Understanding

Answer each question using the model and assumptions stated above.

In the PRS-CS prior, what changes when the global scale $\phi$ decreases while all local variances remain fixed?
Why is the expected marginal estimate vector $R\beta$, rather than $\beta$, when variants are in LD?
For fixed $\phi$, which variant has the wider conditional prior in an example where $\psi_1=2$ and $\psi_2=8$? Explain using the prior variances.
In the two-variant conditional calculation, why is the term $0.60(0.020)$ subtracted from $\hat b_1$?
Why does a worked calculation using an LD matrix not demonstrate that the chosen reference panel is appropriate?
What is the distinction between the conditional mean in the update example and a posterior mean weight used in a finished score?