Worked Linear-Trait Association Scan

From a Manhattan-Plot Signal to a Single-SNP Model

A Manhattan plot can identify a genomic region with a small association p-value, but the plot does not show the estimated trait change, its uncertainty, or the covariates included in the test. A linear-trait association scan supplies those details one variant at a time. At each SNP, the central question is whether the quantitative phenotype is associated with allele dosage after adjustment for specified covariates.

The model below describes an association, not necessarily a causal effect. Its coefficient is conditional on the other included predictors and on the stated linear model. The same analysis is repeated across variants in a scan; the examples here focus on the calculation and interpretation for one SNP.

Definition: For an individual $i$, let $Y_i$ be a quantitative trait, $G_i$ the dosage of a defined effect allele, and $C_{ij}$ the value of covariate $j$. A covariate-adjusted single-SNP linear model is $$ Y_i=\beta_0+\beta_G G_i+\sum_{j=1}^{q}\gamma_j C_{ij}+\varepsilon_i. $$ Here $\beta_G$ is the adjusted mean difference in the trait per one-unit increase in effect-allele dosage, holding the listed covariates fixed.

For a genotyped diploid SNP, $G_i$ is commonly coded as 0, 1, or 2 copies of the effect allele. With imputed genotype data, dosage may be fractional: it is the expected allele count given the genotype probabilities. For example, probabilities 0.1, 0.6, and 0.3 for 0, 1, and 2 copies give dosage $0(0.1)+1(0.6)+2(0.3)=1.2$. Dosage coding uses an additive model: each additional copy is assumed to contribute the same change in the conditional mean. That is a modeling choice, not a fact guaranteed by the data.

Model Assumptions and the Adjusted Coefficient

In matrix notation, write $\mathbf y=\mathbf X\boldsymbol\theta+\boldsymbol\varepsilon$, where $\mathbf X$ contains an intercept, genotype dosage, and the covariate columns. Ordinary least squares (OLS) requires a full-rank design matrix to estimate all coefficients uniquely. For the usual exact small-sample t inference, assume errors have conditional mean zero, are independent, and have common variance $\sigma^2$; normal errors give the exact t distribution. These assumptions concern the errors after conditioning on the model predictors. In large samples, approximate inference may be useful without normal errors, but dependence or incorrectly modeled variation can still invalidate ordinary standard errors.

The coefficient $\hat\beta_G$ can be understood by removing the linear contributions of the intercept and covariates from both the trait and the genotype. Let $\mathbf Z$ contain the intercept and covariates, and define $\mathbf M_Z=\mathbf I-\mathbf Z(\mathbf Z^\mathsf{T}\mathbf Z)^{-1}\mathbf Z^\mathsf{T}$, assuming $\mathbf Z$ has full column rank. Then the Frisch–Waugh–Lovell result gives

\hat\beta_G= \frac{\mathbf g^\mathsf{T}\mathbf M_Z\mathbf y} {\mathbf g^\mathsf{T}\mathbf M_Z\mathbf g}.

The numerator is the cross-product of genotype and trait after each has been linearly adjusted for $\mathbf Z$; the denominator is the remaining genotype sum of squares. To justify the expression, the OLS residual from regressing $\mathbf y$ on $\mathbf Z$ is $\mathbf M_Z\mathbf y$, and the residualized genotype is $\mathbf M_Z\mathbf g$. Regressing the former on the latter gives slope $(\mathbf g^\mathsf{T}\mathbf M_Z\mathbf g)^{-1}\mathbf g^\mathsf{T}\mathbf M_Z\mathbf y$, since $\mathbf M_Z$ is symmetric and idempotent. This also shows why adjusting only the trait, or treating genotype as orthogonal to the covariates without checking, need not give the adjusted SNP effect.

With $p$ estimated coefficients, including the intercept and SNP, the residual degrees of freedom are $df=n-p$. Under the classical OLS assumptions, estimate the residual variance by $\hat\sigma^2=\mathrm{SSE}/df$, where $\mathrm{SSE}$ is the residual sum of squares. The standard error and test statistic for $H_0:\beta_G=0$ are

SE(\hat\beta_G)= \sqrt{\frac{\hat\sigma^2}{\mathbf g^\mathsf{T}\mathbf M_Z\mathbf g}}, \qquad t=\frac{\hat\beta_G}{SE(\hat\beta_G)}.

Under the null and the model assumptions, $t$ follows a t distribution with $df$ degrees of freedom. A two-sided p-value is $2P(T_{df}\geq |t|)$, where $T_{df}$ has that t distribution. A $100(1-\alpha)\%$ confidence interval is $\hat\beta_G\pm t_{1-\alpha/2,df}SE(\hat\beta_G)$. It describes uncertainty about the adjusted mean change in trait units per additional effect-allele copy.

Worked Example: Fitting an Adjusted SNP Effect

Consider a hypothetical dataset of eight individuals. Let $G=(1,1,1,1,0,1,1,3)$ be effect-allele dosage and let $Z=(1,1,1,1,-1,-1,-1,-1)$ be one centered covariate. To make the hypothetical trait values fully specified, define $q=(0,0,0,0,-2,3,0,-1)$, $c=\sqrt{345/1064}$, and $Y_i=1731/152+(23/19)G_i+(251/152)Z_i+cq_i$. Thus the displayed outcome is generated by stated coefficients plus a specified residual vector; it is not presented as observed study data.

The residual vector is orthogonal to the three design columns: $\sum q_i=0$, $\sum G_iq_i=0$, and $\sum Z_iq_i=0$. Consequently the OLS coefficients are the coefficients in the definition of $Y_i$. To verify the adjusted genotype slope directly, the centered cross-products are $S_{GG}=39/8$, $S_{GZ}=-1$, $S_{ZZ}=8$, $S_{GY}=17/4$, and $S_{ZY}=12$. Residualizing both genotype and trait against the intercept and $Z$ gives

\mathbf g^\mathsf{T}\mathbf M_Z\mathbf g =\frac{39}{8}-\frac{(-1)^2}{8} =\frac{19}{4}, \qquad \mathbf g^\mathsf{T}\mathbf M_Z\mathbf y =\frac{17}{4}-\frac{(-1)(12)}{8} =\frac{23}{4}.

Therefore $\hat\beta_G=(23/4)/(19/4)=23/19=1.211$ trait units per effect-allele copy. The adjusted covariate coefficient is $251/152=1.651$, and the intercept is $1731/152=11.39$. Genotype is not orthogonal to either the intercept or $Z$: its sum is 9 and its centered cross-product with $Z$ is $-1$. Omitting those adjustments would therefore not justify using the unadjusted genotype sum of squares.

There are three fitted coefficients, so $df=8-3=5$. The residual sum of squares is $\mathrm{SSE}=c^2\sum q_i^2=(345/1064)(14)=345/76=4.539$, and $\hat\sigma^2=(345/76)/5=69/76=0.9079$. Thus

SE(\hat\beta_G)= \sqrt{\frac{69/76}{19/4}} =\frac{\sqrt{69}}{19} =0.4372, \qquad t=\frac{23/19}{\sqrt{69}/19} =\frac{23}{\sqrt{69}} =2.769.

Using $df=5$, the two-sided p-value is approximately 0.039, and $t_{0.975,5}=2.5706$. The 95% confidence interval is $23/19\pm2.5706(0.4372)=[0.0867,2.3344]$ trait units per effect-allele copy. Under the model, the data provide evidence against a zero adjusted linear dosage effect at the 5% level. The interval is a range of parameter values compatible with this procedure and data under its assumptions; it is not a 95% probability statement about a fixed parameter.

From Model Output to a Scan Result

A scan report should make the allele direction, trait scale, adjustment set, and uncertainty recoverable. A coefficient without allele labels is ambiguous: reversing which allele is called the effect allele changes the sign of the coefficient, though not the two-sided p-value. Covariates should be specified based on the analysis plan and scientific context, rather than chosen after inspecting which set gives the smallest p-value.

Worked Example: Reporting a Scan Result

Suppose a hypothetical scan analyzes a standardized quantitative trait in 1,203 unrelated participants, with one SNP and three covariates plus an intercept. The effect allele is A, the reference allele is G, and the fitted dosage coefficient is $0.018$ standard-deviation units per A-allele copy, with $SE=0.005$. The model has $df=1203-5=1198$, so $t=0.018/0.005=3.6$.

The two-sided t-test p-value with 1,198 degrees of freedom is approximately 0.0003312. A 95% interval, using a t critical value of about 1.962, is $0.018\pm1.962(0.005)=[0.00819,0.02781]$ standard-deviation units per A-allele copy. A concise report is: “In this hypothetical scan, each additional A-allele copy was associated with a 0.018-SD higher trait value after adjustment for the three stated covariates (SE 0.005; 95% CI 0.00819 to 0.02781; two-sided $p=0.0003312$).” This reports the estimated association and its uncertainty without implying causation.

The p-value is calculated for this single coefficient under the stated model. A genome-wide scan tests many variants, so a single-variant p-value should not be described as if it automatically accounts for the full scan. Scan-wide interpretation and thresholding are separate from the coefficient calculation shown here.

Worked Example: Dosage Scale and Effect Direction

Suppose a second hypothetical analysis reports $\hat\beta_G=-0.42$ kg with $SE=0.12$ kg per copy of the effect allele T. The genotype dosage ranges from 0 to 2, so an increase of one dosage unit corresponds to one additional expected T-allele copy. The Wald statistic is $-0.42/0.12=-3.50$. With a large residual degrees of freedom, a 95% interval using a critical value of about 1.96 is $-0.42\pm1.96(0.12)=[-0.655,-0.185]$ kg per T-allele copy.

The negative estimate means that, conditional on the included covariates and under the additive model, one additional T-allele copy is associated with an estimated 0.42 kg lower mean trait value. It does not mean that every carrier has a trait value exactly 0.42 kg lower, nor does it establish that the allele causes the difference. If the other allele were instead labeled the effect allele, the coefficient would reverse sign. A clear report must therefore name the allele and preserve the scale.

Diagnostics, Alternatives, and Limitations

Before interpreting a linear scan result, check the aspects of the analysis that can make its estimate or uncertainty misleading. Earlier material on residual diagnostics for quantitative-trait GWAS and heteroskedasticity and standard errors in linear GWAS gives fuller treatments; the essential checks here are:

Residual behavior: inspect residual-versus-fitted plots for curvature or changing spread, and check for extreme residuals or influential observations. A linear conditional mean that is substantially wrong can distort the dosage coefficient.
Genotype and design: verify dosage ranges, allele labels, missingness handling, and that genotype is not perfectly collinear with covariates. If $\mathbf g^\mathsf{T}\mathbf M_Z\mathbf g=0$, the SNP effect cannot be estimated separately from those predictors.
Dependence: ordinary OLS standard errors presume independent errors. Relatedness or other sample dependence can make them too small or otherwise miscalibrated; a standard error calculation cannot repair an unsuitable independence assumption.
Variance and distribution: strong heteroskedasticity undermines the classical OLS standard error. Heteroskedasticity-robust standard errors can address some variance misspecification, but do not fix confounding, dependence, or a wrongly specified mean model.
Model form: additive dosage imposes a linear per-copy effect. If that form is inappropriate, the coefficient is a model-dependent summary. A different genotype model should be specified for a scientific reason, not selected solely because it gives the strongest result.

The main alternative standard-error calculations include classical OLS, which uses the common-variance model, and heteroskedasticity-robust OLS, which allows unequal error variances under appropriate independence and regularity conditions. These approaches may give different uncertainty estimates; neither automatically handles related individuals. For a quantitative phenotype with substantial relatedness or other structured dependence, a model that represents that covariance may be needed. That changes the inferential model, not merely the presentation of the same OLS result.

Common pitfall: “The genotype sum of squares is 4.875, so use $SE=\sqrt{\hat\sigma^2/4.875}$.” This is wrong when genotype is associated with the intercept or covariates. In the eight-person example, accounting for the intercept and $Z$ gives $\mathbf g^\mathsf{T}\mathbf M_Z\mathbf g=4.75$, not 4.875. The correction is to residualize genotype against every non-SNP design column, or equivalently use the corresponding diagonal element of the full OLS coefficient covariance matrix.

A second frequent error is to report a coefficient “per allele” without saying which allele or which trait units are used. A third is to read a small p-value as the probability that the null hypothesis is true. The p-value is instead the probability, assuming the null and test assumptions, of a test statistic at least as extreme as the one observed. Statistical evidence, estimated magnitude, uncertainty, model adequacy, and scan-wide interpretation should be considered together.

Check Your Understanding

Try each question, then compare your reasoning with the solutions.

A SNP is coded as 0, 1, or 2 copies of allele C. What does a fitted dosage coefficient of 0.30 mean if the trait is measured in centimeters?
In a model with an intercept, genotype, and four covariates, what are the residual degrees of freedom if $n=500$?
Why is $\mathbf g^\mathsf{T}\mathbf M_Z\mathbf g$, rather than simply $\sum_i(G_i-\bar G)^2$, used in the adjusted coefficient standard error?
A result is reported as $\hat\beta_G=-0.5$, $SE=0.2$, but the effect allele and trait units are omitted. Name two pieces of information needed to interpret it.

Solutions. (1) Under the additive model and conditional on the included covariates, each additional C-allele copy is associated with an estimated 0.30 cm increase in the conditional mean trait. (2) There are six fitted coefficients, so $df=500-6=494$. (3) The adjusted uncertainty depends on genotype variation remaining after its linear association with the intercept and covariates is removed; the raw centered sum of squares does not generally equal this residualized quantity. (4) State which allele is counted as the effect allele and the trait scale and units; also give the covariate adjustment set and uncertainty measure for a complete report.

Worked Linear-Trait Association Scan

What You'll Learn