Single-Variant Association with Polygenic Random Effects

Why Test a SNP in a Mixed Model?

Individuals in a genetic cohort may be related or share genome-wide ancestry. As a result, their phenotypes can be correlated even after measured covariates are included. A standard regression that treats all observations as independent may then give an inaccurate standard error for a candidate SNP effect. A linear mixed model addresses this by representing polygenic similarity through a genomic relationship matrix, while estimating the candidate SNP effect directly as a fixed effect.

The target is an association between dosage at one candidate variant and a quantitative phenotype, conditional on the covariates and the covariance model. The SNP effect is not itself the polygenic random effect: it is a specific, interpretable regression coefficient. The random effect accounts for residual genetic resemblance among individuals that is not represented by that candidate coefficient.

Definition: For $n$ individuals, let $y$ be the $n$-vector of quantitative phenotypes, $C$ an $n \times q$ matrix of covariates (including an intercept), and $x$ the $n$-vector of candidate-SNP dosages. The single-variant linear mixed model is $$ y = C\gamma + x\beta + g + e, $$ where $\gamma$ is a vector of covariate effects, $\beta$ is the fixed candidate-SNP effect, $g$ is a polygenic random effect, and $e$ is residual error. Assume $g \sim N(0,\sigma_g^2 K)$, $e \sim N(0,\sigma_e^2 I)$, and $g$ and $e$ are independent. Here $K$ is a specified genomic relationship matrix and $I$ is the $n \times n$ identity matrix. Thus $y$, conditional on the design, has covariance $V=\sigma_g^2K+\sigma_e^2I$.

Dosage usually counts copies of a designated effect allele, with values 0, 1, or 2; imputed dosages can take fractional values between 0 and 2. Under this coding, $\beta$ is the expected phenotype difference for a one-unit increase in dosage, holding covariates and the modeled covariance structure fixed. If the phenotype is measured in millimetres, for example, $\beta$ is measured in millimetres per additional effect-allele copy.

The model assumes that the conditional mean is linear in the covariates and dosage, that the random effects have the stated mean-zero Gaussian distributions, and that the covariance matrix $K$ is suitable for representing the genetic covariance. It also assumes the residual component is independent across individuals under $\sigma_e^2I$, and independent of $g$. The genotype and covariates are treated as fixed design values for this analysis. These assumptions describe the model; they do not guarantee that a candidate-SNP association is causal.

Estimating and Testing the Candidate Effect

If the variance components are known, generalized least squares (GLS) estimates the fixed effects using the inverse covariance $V^{-1}$. Write $A=[C\ \ x]$ for the design matrix containing the covariates and candidate dosage, and $\theta=(\gamma^\mathsf{T},\beta)^\mathsf{T}$. Then the GLS estimate and its covariance are

\hat{\theta}=(A^\mathsf{T}V^{-1}A)^{-1}A^\mathsf{T}V^{-1}y, \qquad \operatorname{Var}(\hat{\theta}\mid A)= (A^\mathsf{T}V^{-1}A)^{-1}.

In practice, $\sigma_g^2$ and $\sigma_e^2$ are usually estimated from the data, commonly by REML for variance-component estimation. The fitted covariance is then used to obtain the fixed-effect estimate and its standard error. The standard error for $\hat{\beta}$ is the square root of its estimated variance. Standard Wald inference treats the fitted variance components as plug-in estimates and is generally an approximation, especially when the sample is small.

The usual association hypotheses are $H_0:\beta=0$ against $H_a:\beta\ne 0$. A two-sided Wald statistic is $Z=\hat{\beta}/SE(\hat{\beta})$; under the null and the usual large-sample approximation, it is compared with a standard normal distribution. A confidence interval is $\hat{\beta}\pm z_{1-\alpha/2}SE(\hat{\beta})$. A likelihood-ratio test instead compares the maximized likelihood under the null and alternative fixed-effect models. Its one-degree-of-freedom reference distribution is asymptotically chi-square under the regularity conditions for this fixed-effect comparison.

Worked Example: GLS With a Simple Relatedness Pattern

Suppose four hypothetical individuals have phenotypes $y=(1,2,4,5)^\mathsf{T}$, candidate dosages $x=(0,0,1,1)^\mathsf{T}$, and an intercept as the only covariate. For a transparent calculation, suppose the variance components are known and $$ V= \begin{pmatrix} 2&1&0&0\\ 1&2&0&0\\ 0&0&2&1\\ 0&0&1&2 \end{pmatrix}. $$ This covariance has positive within-pair covariance for individuals 1 and 2 and for individuals 3 and 4; different pairs have zero covariance. It is positive definite because each two-by-two block has eigenvalues 1 and 3. The example treats this as the fitted covariance implied by a genomic relationship structure and positive variance components.

For either pair, the inverse block is $\frac{1}{3}\begin{pmatrix}2&-1\\-1&2\end{pmatrix}$. Thus $V^{-1}$ consists of two such blocks. Let $A$ have columns for the intercept and dosage. The resulting GLS quantities are $$ A^\mathsf{T}V^{-1}A= \begin{pmatrix} 4/3&2/3\\ 2/3&2/3 \end{pmatrix}, \qquad A^\mathsf{T}V^{-1}y= \begin{pmatrix} 4\\3 \end{pmatrix}. $$

The determinant of the first matrix is $4/9$, so its inverse is $$ (A^\mathsf{T}V^{-1}A)^{-1} = \begin{pmatrix} 3/2&-3/2\\ -3/2&3 \end{pmatrix}. $$ Multiplying by the second vector gives $\hat{\theta}=(1.5,3)^\mathsf{T}$. Therefore the fitted intercept is 1.5 phenotype units, and the candidate-SNP effect is $\hat{\beta}=3$ phenotype units per dosage unit. With $V$ treated as known, the standard error is $\sqrt{3}=1.732$, giving $Z=3/1.732=1.732$ and an approximate two-sided normal p-value of 0.0833. This small hypothetical example does not provide strong evidence against a zero SNP effect at the 0.05 level. The calculation also shows that the covariance affects both the GLS weights and the uncertainty, rather than merely changing the reported interpretation after an ordinary regression.

Wald Inference in Phenotype Units

A test statistic is most useful when paired with the effect estimate and its uncertainty. The p-value addresses how incompatible the data are with the specified null under the model; it does not measure the size or clinical importance of the effect. The estimate and confidence interval give the effect's direction and scale, in the phenotype's units per dosage unit. A confidence interval procedure has its stated long-run coverage under the model and approximation; it does not assign a probability to the fixed parameter after the observed interval has been calculated.

Worked Example: Interpreting an Estimated Effect and Interval

Suppose a hypothetical mixed-model analysis uses a phenotype measured in centimetres and codes the candidate effect allele by dosage. The fitted effect is $\hat{\beta}=-0.180$ cm per allele copy, with $SE(\hat{\beta})=0.0600$ cm per copy. Assume the large-sample normal approximation is appropriate for this illustration.

The Wald statistic is $$ Z=\frac{-0.180}{0.0600}=-3.00. $$ The two-sided p-value is $2P(Z_0\ge 3.00)\approx 0.00270$, where $Z_0$ is standard normal. A 95% confidence interval is $$ -0.180 \pm 1.96(0.0600) =-0.180 \pm 0.1176 =(-0.2976,-0.0624)\text{ cm per allele copy}. $$

Under this model, the estimated mean phenotype is 0.180 cm lower for each additional effect-allele copy, conditional on the included covariates and the modeled genetic covariance. The interval excludes zero, consistent with the two-sided test at the 0.05 level. It describes uncertainty in the estimated association; it does not show that changing a person's genotype would cause the phenotype to change by that amount.

Likelihood-Ratio Testing

For a likelihood-ratio test (LRT), the null model includes the covariates and genomic random effect but omits the candidate dosage. The alternative adds $x\beta$. Fit both models with the same individuals, phenotype, covariates, and covariance structure. For comparing models that differ in fixed effects, maximize the ordinary likelihood under each model (using ML), rather than comparing REML likelihoods from different fixed-effect design matrices. Let $\ell_0$ and $\ell_1$ denote the maximized log likelihoods under the null and alternative, respectively. The statistic is $2(\ell_1-\ell_0)$. Under $H_0$, with a single added fixed-effect parameter and standard regularity conditions, it is approximately $\chi^2_1$. The variance components are nuisance parameters estimated in both fits; this is not a test that a variance component equals zero.

Worked Example: A One-Parameter Likelihood-Ratio Test

Suppose a hypothetical analysis gives maximized ML log likelihoods $\ell_0=-120.4$ for the model without the candidate dosage and $\ell_1=-117.9$ for the model with dosage. Both fits use the same 300 individuals, covariates, and random-effect covariance specification. The alternative has one additional fixed coefficient, $\beta$.

The test statistic is $$ 2(\ell_1-\ell_0)=2[-117.9-(-120.4)]=2(2.5)=5.00. $$ Comparing this with a $\chi^2_1$ distribution gives $p\approx 0.0253$. Under the stated model and asymptotic reference distribution, this is evidence against $H_0:\beta=0$ at the 0.05 level. The result supports an association conditional on the specified covariates and covariance model; it does not establish causation. The LRT p-value alone does not give the effect's size, so the estimate and its standard error or interval should also be reported.

Why the Covariance Model Matters

The key distinction is between a mean effect and a covariance effect. The candidate dosage enters the conditional mean through the fixed coefficient $\beta$. The genomic relationship matrix enters the covariance through $V=\sigma_g^2K+\sigma_e^2I$. By modeling covariance, the analysis avoids treating genetically similar individuals as though their phenotypes necessarily supplied fully independent evidence. The SNP coefficient remains a population-level conditional association estimate, not an individual-specific prediction of the polygenic effect.

A valid analysis still depends on an appropriate model and data handling. The genotype coding and effect allele should be stated; the phenotype and dosage must be aligned to the same individuals; and covariates should be chosen for the scientific design. The relationship matrix must be scaled and constructed consistently with the variance-component interpretation. Missingness, poor genotype quality, or an inappropriate covariance structure can undermine inference. A small p-value cannot repair those problems.

Common Mistakes

Treating the SNP as the random effect. The candidate SNP is a fixed dosage predictor with coefficient $\beta$; the polygenic effect $g$ is random and has covariance $\sigma_g^2K$.
Interpreting $\beta$ without its units or coding. State which allele dosage is counted and the phenotype change per one-unit dosage increase.
Using ordinary-regression uncertainty after fitting a relatedness model. The standard error must reflect the fitted covariance; otherwise the reported test may not match the mixed model.
Reading a p-value as an effect size or probability that the null is true. A p-value is calculated under the null model and measures tail-area extremeness, not the probability that $H_0$ is true.
Calling an association causal. The coefficient describes a conditional association under the specified model; it does not by itself identify a causal variant or mechanism.

Key takeaway: Single-variant mixed-model association analysis puts the candidate dosage in the fixed-effect mean and uses the genomic relationship matrix to model polygenic covariance. Report the dosage coding, effect in phenotype units, uncertainty, and test result together.

Check Your Understanding

Try these questions without looking back.

In the model $y=C\gamma+x\beta+g+e$, which term represents the candidate-SNP association and which term represents polygenic covariance?
If dosage counts copies of an effect allele and $\hat{\beta}=0.4$ mm, what is the estimated phenotype difference associated with a one-copy increase, in the model's conditional interpretation?
Why does the covariance matrix affect the standard error of the SNP effect?
A Wald test gives $Z=2.2$. What is the two-sided p-value formula, and what does that p-value not mean?
Why should ML likelihoods, rather than REML likelihoods from different fixed-effect designs, be compared for the fixed-effect LRT described here?