Deriving the Summary-Statistic Likelihood With Linkage Disequilibrium

Why Marginal Association Statistics Must Be Read Together

When two variants are in linkage disequilibrium (LD), their genotypes tend to vary together. Consequently, their marginal association statistics are correlated, even when only one variant has a nonzero effect in the model. A large statistic at a variant can therefore reflect its correlation with another variant’s genotype, rather than a separate effect at that variant.

Fine-mapping needs a likelihood that describes this dependence. Treating each marginal statistic as an independent observation discards the LD structure and can misrepresent how well a proposed pattern of variant effects fits the data. This tutorial derives an approximate joint likelihood for marginal association statistics conditional on the locus LD matrix. The likelihood is a model for the observed statistics given effect sizes and LD; it is not, by itself, a posterior probability that any variant is causal.

Definition: Let $R$ be the $m \times m$ LD matrix for $m$ variants. Its entry $R_{jk}$ is the correlation between the standardized genotypes at variants $j$ and $k$; in particular, $R_{jj}=1$. Let $\mathbf{z}$ be the vector of marginal association $z$-statistics and let $\boldsymbol{\theta}$ be the vector of standardized effects in the joint genetic model. Under the approximation developed here, the summary-statistic likelihood is $$ \mathbf{z}\mid\boldsymbol{\theta},R \ \approx\ N\!\left(\sqrt{N}\,R\boldsymbol{\theta},\,R\right). $$ Here $N$ is the common sample size, and the covariance matrix in the normal distribution is $R$.

Model and Assumptions

The derivation is most direct for a quantitative trait measured in the same $N$ individuals at every variant. Let $X$ be the $N \times m$ genotype matrix after centering and scaling each genotype column to have sample variance one, so that $X'X/N=R$. Let $y$ be the centered trait vector, and suppose the joint linear model is

y=X\boldsymbol{\beta}+\boldsymbol{\varepsilon}, \qquad \boldsymbol{\varepsilon}\sim N(\mathbf{0},\sigma^2 I).

The model treats residuals as independent, mean-zero normal errors with common variance $\sigma^2$, conditional on the genotypes. The vector $\boldsymbol{\beta}$ contains the joint effects of the variants: $\beta_j$ is the change in trait units per one standard deviation of genotype at variant $j$, conditional on the other variants in this model. Define $\boldsymbol{\theta}=\boldsymbol{\beta}/\sigma$, so each effect is expressed in residual-standard-deviation units. Effects equal to zero contribute no term to the joint mean.

For clarity, the derivation omits covariates. In an analysis that adjusts for covariates, the relevant genotypes and trait must be residualized with respect to those covariates, and the LD matrix and sample size must match the association statistics. The approximation also presumes consistent allele orientation: if one variant’s effect or genotype is coded on the opposite allele, its signs must be harmonized before interpreting the vector.

In practice, marginal statistics often use estimated standard errors and may be $t$-statistics rather than exact normal score statistics. The multivariate normal form is therefore an approximation, typically most useful when the sample is large and the statistics and LD reference describe the same population and analysis. It can be less reliable when sample sizes differ across variants, LD is poorly matched, or unaccounted-for structure affects the statistics.

Deriving the Joint Distribution

For each variant $j$, the standardized genotype column has squared length $N$. The vector of marginal regression coefficient estimates can therefore be written, under the stated scaling, as $\widehat{\boldsymbol{\beta}}_{\mathrm{marg}}=X'y/N$. These are marginal estimates: each uses one genotype at a time, rather than conditioning on all variants.

\begin{aligned} E(\widehat{\boldsymbol{\beta}}_{\mathrm{marg}}\mid X) &= \frac{X'X}{N}\boldsymbol{\beta} =R\boldsymbol{\beta},\\ \operatorname{Var}(\widehat{\boldsymbol{\beta}}_{\mathrm{marg}}\mid X) &= \frac{1}{N^2}X'\operatorname{Var}(\boldsymbol{\varepsilon})X = \frac{\sigma^2}{N}R. \end{aligned}

The first line follows by substituting $y=X\boldsymbol{\beta}+\boldsymbol{\varepsilon}$ and using $E(\boldsymbol{\varepsilon}\mid X)=\mathbf{0}$. The second follows from $\operatorname{Var}(\boldsymbol{\varepsilon}\mid X)=\sigma^2 I$ and $X'X/N=R$. Thus, a marginal estimate at one variant can have a nonzero mean because its genotype is correlated with genotypes at variants that have nonzero joint effects.

Scale each marginal estimate by its approximate standard error $\sigma/\sqrt{N}$, defining $\mathbf{z}=\sqrt{N}\widehat{\boldsymbol{\beta}}_{\mathrm{marg}}/\sigma$. A linear transformation of a multivariate normal vector is multivariate normal, so the preceding mean and covariance imply

\mathbf{z}\mid\boldsymbol{\theta},R \ \approx\ N\!\left(\sqrt{N}\,R\boldsymbol{\theta},\,R\right).

This result has two distinct parts. The mean, $\sqrt{N}R\boldsymbol{\theta}$, describes the association pattern expected from joint effects after each variant is tested marginally. The covariance, $R$, describes how sampling noise in the marginal statistics is correlated. LD therefore affects both the expected pattern and the uncertainty around it.

The Likelihood in Explicit Form

If $R$ is positive definite, the multivariate normal density gives the likelihood of $\mathbf{z}$ for a proposed $\boldsymbol{\theta}$:

L(\boldsymbol{\theta};\mathbf{z},R) = \frac{1}{(2\pi)^{m/2}|R|^{1/2}} \exp\!\left[ -\frac{1}{2} \left(\mathbf{z}-\sqrt{N}R\boldsymbol{\theta}\right)' R^{-1} \left(\mathbf{z}-\sqrt{N}R\boldsymbol{\theta}\right) \right].

Equivalently, expanding the quadratic gives a form that makes the effect-dependent terms visible:

\log L = -\frac{m}{2}\log(2\pi)-\frac{1}{2}\log|R| -\frac{1}{2}\mathbf{z}'R^{-1}\mathbf{z} +\sqrt{N}\boldsymbol{\theta}'\mathbf{z} -\frac{N}{2}\boldsymbol{\theta}'R\boldsymbol{\theta}.

To verify the expansion, write $\mathbf{a}=\sqrt{N}R\boldsymbol{\theta}$ and expand $(\mathbf{z}-\mathbf{a})'R^{-1}(\mathbf{z}-\mathbf{a})$. The cross terms are equal because the expression is a scalar and $R$ is symmetric. Since $R^{-1}\mathbf{a}=\sqrt{N}\boldsymbol{\theta}$, the result is $\mathbf{z}'R^{-1}\mathbf{z}-2\sqrt{N}\boldsymbol{\theta}'\mathbf{z}+N\boldsymbol{\theta}'R\boldsymbol{\theta}$, as shown.

When comparing proposed effects for the same observed statistics and the same $R$, terms that do not depend on $\boldsymbol{\theta}$ can be omitted from a likelihood ratio. They must not be omitted if a complete density is required. If $R$ is singular, the displayed density involving $R^{-1}$ does not apply as written. Perfectly correlated variants, for example, can make it singular; the data then contain redundant directions, and a nonsingular approximation or a formulation on the supported subspace is needed.

Worked Example: Evaluating a Joint Likelihood

Suppose two variants have $R_{12}=0.8$, $N=10{,}000$, and the proposed standardized effects are $\boldsymbol{\theta}=(0.02,0)'$. The observed marginal statistics are $\mathbf{z}=(2.1,1.5)'$. Calculate the likelihood density under this proposed effect vector.

The LD matrix, its determinant, and its inverse are

R= \begin{pmatrix}1&0.8\\0.8&1\end{pmatrix}, \qquad |R|=1-0.8^2=0.36, \qquad R^{-1}=\frac{1}{0.36} \begin{pmatrix}1&-0.8\\-0.8&1\end{pmatrix}.

Because $\sqrt{N}=100$, the expected statistic vector is $100R(0.02,0)'=(2,1.6)'$. The residual vector is therefore $\mathbf{r}=(0.1,-0.1)'$. Its quadratic form is

\mathbf{r}'R^{-1}\mathbf{r} = \frac{0.1^2+(-0.1)^2-1.6(0.1)(-0.1)}{0.36} = \frac{0.036}{0.36} =0.1.

Substituting into the density gives

L = \frac{\exp(-0.1/2)}{2\pi\sqrt{0.36}} = \frac{\exp(-0.05)}{2\pi(0.6)} \approx 0.2523.

This is a density value for the two-dimensional statistic vector, not the probability of observing exactly $(2.1,1.5)'$. The small quadratic discrepancy indicates that the observed vector is close to its expected value under this specified effect vector, relative to the covariance $R$.

How LD Shapes the Evidence Pattern

If only variant $j$ has a nonzero standardized effect $\theta_j$, the expected vector simplifies to

E(\mathbf{z}\mid\theta_j,R) = \sqrt{N}\theta_j R_{\cdot j},

where $R_{\cdot j}$ is column $j$ of the LD matrix. Thus a variant’s marginal statistic is expected to be proportional to its correlation with the variant carrying the effect. Strongly correlated variants can show similarly large statistics; variants correlated in the opposite direction can show statistics with opposite signs. This pattern is why the magnitude of one marginal statistic alone does not identify the variant responsible for the association.

Worked Example: Comparing Two LD-Predicted Patterns

Consider two variants with correlation $0.8$, sample size $10{,}000$, and observed statistics $\mathbf{z}=(3,2.4)'$. Compare two specified models: model A has $\boldsymbol{\theta}=(0.03,0)'$, while model B has $\boldsymbol{\theta}=(0,0.03)'$.

Under model A, the expected vector is $100R(0.03,0)'=(3,2.4)'$, exactly the observed vector, so the residual quadratic form is zero. Under model B, the expected vector is $100R(0,0.03)'=(2.4,3)'$, leaving residual $\mathbf{r}=(0.6,-0.6)'$. Using the same inverse LD matrix as above,

\mathbf{r}'R^{-1}\mathbf{r} = \frac{0.6^2+(-0.6)^2-1.6(0.6)(-0.6)}{0.36} = \frac{1.296}{0.36} =3.6.

The log-likelihood under model A exceeds that under model B by $3.6/2=1.8$, since the normalizing terms are identical. The likelihood ratio for these specified effects is $\exp(1.8)\approx 6.050$ in favor of model A. This comparison is conditional on the two fixed effect vectors: it is not a posterior odds ratio, and it does not account for prior probabilities or uncertainty in the effect sizes.

Worked Example: Two Effects Contribute to Every Marginal Statistic

Suppose three variants have LD matrix

R= \begin{pmatrix} 1&0.5&0.2\\ 0.5&1&0.4\\ 0.2&0.4&1 \end{pmatrix},

and take $N=10{,}000$ with standardized effects $\boldsymbol{\theta}=(0.02,0,-0.01)'$. The expected statistic vector is obtained by multiplying each row of $R$ by the effect vector and then by $100$:

\begin{aligned} E(\mathbf{z}) &=100R \begin{pmatrix}0.02\\0\\-0.01\end{pmatrix}\\ &=100 \begin{pmatrix} 0.02-0.002\\ 0.01-0.004\\ 0.004-0.01 \end{pmatrix} = \begin{pmatrix}1.8\\0.6\\-0.6\end{pmatrix}. \end{aligned}

Although the second variant has zero joint effect in this model, its expected marginal statistic is $0.6$, because it is correlated with both variants that have nonzero effects. The third expected statistic is negative: its positive correlation with the first effect and its own negative effect contribute in opposite directions, with the negative contribution larger here. The covariance of the three statistics is also $R$, so their random deviations from these expectations are correlated.

Interpretation, Scope, and Common Pitfalls

The joint likelihood gives a coherent measure of how compatible the observed vector of marginal statistics is with a proposed vector of effects, given the assumed LD structure. In Bayesian fine-mapping, this likelihood can be combined with prior information to obtain posterior quantities. The likelihood alone does not assign posterior probabilities, choose a causal-count prior, or establish that an associated variant has a biological causal role.

The approximation depends on correctly matched data. If the LD matrix comes from a population with different allele correlations, its predicted patterns and covariance may not represent the association statistics. Allele mismatches can also reverse signs and produce an apparent lack of fit. In addition, variants tested in different samples may have different effective sample sizes and correlations among their statistics that are not captured by a single $N$ and $R$.

Mistake: treating correlated statistics as independent evidence. Multiplying separate univariate likelihoods as though their errors were independent ignores the off-diagonal entries of $R$. The correction is to use the joint covariance, or an equivalent likelihood that retains the LD structure.
Mistake: interpreting a high marginal statistic as proof that the variant itself has a joint effect. The expected marginal statistic depends on the entire vector $R\boldsymbol{\theta}$. A correlated variant can have a substantial marginal statistic even when its own joint effect is zero.
Mistake: calling the likelihood ratio a posterior probability. A likelihood compares how well specified models explain the data. Posterior probabilities additionally depend on prior probabilities and, when effects are uncertain, on integrating over effect sizes.

Key takeaway: Under the standardized linear-model approximation, marginal statistics have mean $\sqrt{N}R\boldsymbol{\theta}$ and covariance $R$. LD therefore determines both how joint effects appear in marginal tests and how much independent information those correlated tests provide.

Check Your Understanding

Use the assumptions and results in this tutorial to answer the following questions.

In the approximate likelihood, what does $R_{jk}$ represent, and where does it enter the distribution of $\mathbf{z}$?
If only variant $j$ has a nonzero joint effect, what is the expected pattern of marginal statistics across the locus?
Why can a variant with zero joint effect have a nonzero expected marginal statistic?
What does a likelihood ratio compare, and why is it not automatically a posterior odds ratio?
Name two data or modeling problems that could make the stated joint likelihood a poor approximation.