Fitting the Null Model Before Genome-Wide Testing

Why fit a null model first?

In a quantitative-trait genome-wide association study (GWAS), the null model describes the trait without including the variant currently being tested. It accounts for measured covariates and for correlation among individuals due to genetic relatedness. Fitting this model before genome-wide testing provides a common estimate of the background mean and covariance structure. The null model is not a claim that genes have no effect: its random genetic effect represents aggregate relatedness-associated variation, while the variant-specific effect is left out.

The earlier material on mixed-model covariance decomposition and restricted maximum likelihood (REML) supplies the main ingredients. Here the emphasis is on putting those ingredients together into a fit-ready model and checking whether the fit is credible. A successful numerical fit does not, by itself, establish that the covariates, sample, or relationship matrix were specified correctly.

Model and assumptions

Definition: For $n$ individuals, the quantitative-trait null linear mixed model is $$ y=X\beta+u+\varepsilon,\qquad u\sim N(0,\sigma_g^2K),\qquad \varepsilon\sim N(0,\sigma_e^2I_n), $$ where $y$ is the $n$-vector of trait values, $X$ is the $n\times p$ fixed-effect design matrix, $\beta$ is its coefficient vector, $K$ is the $n\times n$ genomic relationship matrix, and $I_n$ is the identity matrix. The random genetic effect $u$ and residual error $\varepsilon$ are assumed independent. Thus $$ y\sim N(X\beta,V),\qquad V=\sigma_g^2K+\sigma_e^2I_n. $$ The null model contains no tested-variant term.

The intercept is ordinarily a column of $X$. Other columns may encode prespecified covariates such as age, sex, technical batch, or ancestry principal components. Their choice and coding should reflect the study design and scientific rationale. A column that is constant or exactly redundant with other columns makes $X$ rank deficient; the fixed-effect coefficients then cannot all be uniquely estimated.

The model assumes that the mean is represented adequately by $X\beta$, that $K$ is symmetric and positive semidefinite, and that $\sigma_g^2\geq 0$ and $\sigma_e^2>0$. The latter condition, together with positive semidefinite $K$, makes $V$ positive definite. The model also assumes that the relationship structure represented by $K$ is appropriate for these individuals and that the random effects have the stated covariance. These are modeling assumptions, not facts guaranteed by an optimizer.

Every object must refer to the same individuals in the same order. If the trait or a covariate is missing for an individual, the analysis set must be defined consistently, and $K$ must be subset and reordered to match it. A correctly sized but misordered matrix is still wrong: it assigns covariances to the wrong pairs.

What the fit estimates

For given variance components, the generalized least-squares estimate of the fixed effects is $\hat{\beta}=(X^\mathsf{T}V^{-1}X)^{-1}X^\mathsf{T}V^{-1}y$, provided $X$ has full column rank. The variance components are commonly estimated by REML. REML accounts for the degrees of freedom used to estimate the fixed effects; it is not the same objective as ordinary maximum likelihood.

For a full-rank $X$, define $$ P=V^{-1}-V^{-1}X(X^\mathsf{T}V^{-1}X)^{-1}X^\mathsf{T}V^{-1}. $$ Up to a constant that does not depend on the variance components, the restricted log-likelihood is $$ \ell_R=-\frac{1}{2}\left\{\log|V|+\log|X^\mathsf{T}V^{-1}X|+y^\mathsf{T}Py\right\}. $$ REML chooses nonnegative variance-component estimates that maximize this objective. In practical fitting, software evaluates the objective and its derivatives numerically; convergence should be assessed rather than inferred from the fact that an estimate was returned.

Build the fit before running it

Define the analysis individuals.
Choose the individuals with usable trait and covariate data. Record their identifiers and order; use that same order for every model input.

Construct the fixed-effect design.
Include an intercept and the prespecified covariate columns. Check their coding, ranges, missingness, and linear independence.

Align and inspect $K$.
Subset and reorder the relationship matrix to the analysis individuals. Check its dimensions, symmetry, diagonal, and eigenvalues within numerical tolerance.

Fit and verify.
Estimate $\beta$ and the variance components under the null model. Review convergence, estimates, covariance validity, and sensitivity to reasonable alternative starting values.

The null model is fitted without one column per tested variant. Adding a tested genotype would make the fit a variant-specific alternative model rather than the shared null fit described here. This distinction does not mean that genotype information is irrelevant to constructing $K$: the relationship matrix summarizes genome-wide similarity according to its specified construction.

Worked Example: Specifying the Null Model and Checking Its Inputs

Suppose a hypothetical analysis includes $n=5$ individuals with trait values, in identifier order A through E, of $y=(2.1,\,3.0,\,2.7,\,4.2,\,3.5)^\mathsf{T}$ grams. The planned fixed effects are an intercept, age in years, and one ancestry principal component. The age values are $10,\,12,\,11,\,14,\,13$, and the principal-component values are $-1,\,0,\,0.5,\,1,\,-0.5$. The genotype data for the variants to be tested are not included as fixed effects in this null fit.

The design matrix is $$ X= \begin{pmatrix} 1&10&-1\\ 1&12&0\\ 1&11&0.5\\ 1&14&1\\ 1&13&-0.5 \end{pmatrix}. $$ It has 5 rows and 3 columns. The columns are not constant multiples of one another, and these five rows provide more observations than columns; in an actual fit, full column rank should be checked directly. The coefficient units are trait units for the intercept, trait units per year for the age coefficient, and trait units per unit of the principal component for its coefficient.

Suppose $K$ is supplied as a $5\times5$ symmetric positive-semidefinite matrix in the same A-to-E order. Then the null covariance is $V=\sigma_g^2K+\sigma_e^2I_5$. Before fitting, one should verify that row and column labels both read A through E, rather than relying only on matrix dimensions. The numerical example specifies the model and input checks; it does not claim that these five observations support stable genome-wide variance-component estimation. With so few individuals, estimates would be particularly uncertain.

Interpreting and checking a fitted model

Variance-component estimates describe the fitted covariance under the chosen model. $\hat{\sigma}_g^2$ is the estimated scale of variation associated with $K$, and $\hat{\sigma}_e^2$ is the residual variance scale. Neither estimate alone demonstrates that the model is well specified. A nonzero estimate is not evidence that every pattern of relatedness or confounding has been captured.

Convergence checks should include more than a software flag. Review the optimizer’s termination message and, where available, the size of the gradient at the reported solution. Check that the objective is finite, $V$ is positive definite, and estimates satisfy their constraints. Repeat the fit from more than one reasonable starting value: materially different solutions or log-likelihoods suggest that the reported result may not be a stable optimum. Boundary estimates require particular attention, as discussed in the earlier material on boundary estimates and singular fits.

Also inspect the inputs and fitted model in context. Confirm that the trait and covariates have plausible ranges and units, that the chosen covariates are not redundant, and that $K$ has the expected scale and sample correspondence. Examine residuals with the covariance structure in mind: raw residuals need not be independent under this model. Persistent mean patterns against a covariate can signal a misspecified mean, while extreme observations or unexplained structure can warrant investigation. Such checks diagnose possible problems; they do not automatically prescribe a model change.

Worked Example: Verifying a Reported Fit Numerically

Consider a hypothetical three-individual null-model fit with only an intercept, $K$ having diagonal entries 1 and off-diagonal entries 0.25, and reported variance estimates $\hat{\sigma}_g^2=0.4$ trait-units squared and $\hat{\sigma}_e^2=0.6$ trait-units squared. The trait values are $y=(2,3,5)^\mathsf{T}$ units. These values are stipulated as a fit output so that the covariance and fixed-effect calculations can be checked; they are not claimed to be a REML optimum for this toy data.

The fitted covariance has diagonal entries $0.4(1)+0.6=1.0$ and off-diagonal entries $0.4(0.25)=0.1$: $$ \hat V= \begin{pmatrix} 1&0.1&0.1\\ 0.1&1&0.1\\ 0.1&0.1&1 \end{pmatrix}. $$ Its eigenvalue in the direction of the all-ones vector is $1+2(0.1)=1.2$, and each of the two orthogonal eigenvalues is $1-0.1=0.9$. All are positive, so $\hat V$ is positive definite.

For an intercept-only design, the generalized least-squares estimate is $\hat\beta=(\mathbf{1}^\mathsf{T}\hat V^{-1}y)/(\mathbf{1}^\mathsf{T}\hat V^{-1}\mathbf{1})$. Because each row of $\hat V$ sums to 1.2, $\hat V^{-1}\mathbf{1}=\mathbf{1}/1.2$. Thus the numerator is $(2+3+5)/1.2=8.333$, and the denominator is $3/1.2=2.500$. Therefore $\hat\beta=8.333/2.500=3.333$ trait units. This checks the fixed-effect calculation and positive definiteness, but it does not establish convergence of the variance-component optimizer.

Worked Example: Distinguishing Convergence from a Boundary Warning

Suppose a hypothetical fitting program reports the following results for the same null model from two starting values. Run A terminates with $\hat{\sigma}_g^2=0.32$, $\hat{\sigma}_e^2=0.68$, log-restricted-likelihood $-120.40$, and a small gradient. Run B terminates with $\hat{\sigma}_g^2=0.321$, $\hat{\sigma}_e^2=0.679$, log-restricted-likelihood $-120.40$, and a similarly small gradient. The values are in trait-units squared and the likelihood values are hypothetical.

The variance components differ by $0.001$ each, and the reported likelihoods agree to the displayed precision. Agreement from separate starts, alongside small gradients and a valid covariance matrix, supports the conclusion that the optimizer found a stable interior solution. It is evidence about numerical stability, not proof that the model assumptions are true.

Now suppose a third run reports $\hat{\sigma}_g^2=0$, $\hat{\sigma}_e^2=1.00$, a warning that the estimate is on the boundary, and a higher log-restricted-likelihood of $-119.90$. Since $-119.90>-120.40$, this candidate has a larger restricted likelihood and cannot be dismissed just because it is a boundary solution. The fit should be investigated as a boundary estimate, including the optimizer’s constraints and local behavior. Reporting only the first interior solution would not be justified.

Common failures and alternative approaches

A common mistake is to interpret an optimizer’s “converged” message as a complete validation. The incorrect reasoning is: “The software returned estimates, so the null model is ready.” The correction is to verify alignment, design-matrix rank, the relationship matrix, parameter constraints, gradients, and sensitivity to starting values. Numerical convergence cannot repair a sample-order mismatch or an unsuitable mean model.

Another mistake is to fit ordinary least squares and treat its residuals as independent despite relatedness. That analysis uses a different covariance assumption, $V=\sigma_e^2I_n$, and can misrepresent uncertainty when the actual data contain covariance associated with $K$. A mixed model is appropriate when that covariance structure is part of the intended analysis and is adequately specified. Conversely, adding a relationship matrix does not guarantee protection from every source of confounding or misspecification.

REML is commonly used for variance components when comparing covariance specifications with the same fixed-effect design. Ordinary maximum likelihood is another fitting approach, but it does not make the same adjustment for fixed effects. Estimates and likelihoods from the two objectives should not be treated as interchangeable. The choice of fitting objective should be made deliberately and reported clearly.

The null fit has limits. It depends on the chosen $K$, covariates, and analysis sample; it can be sensitive to poor scaling, redundant covariates, or an incorrect relationship matrix. A boundary estimate may reflect weak information for separating variance components, rather than a computational defect. Finally, fitting a null model is a preparation step, not a substitute for checking the calibration and interpretation of the subsequent association analysis.

Key takeaway: A trustworthy null-model fit requires a clearly specified mean and covariance, correctly aligned inputs, and evidence of numerical stability. A convergence flag is one diagnostic, not the conclusion.

Check Your Understanding

Answer each question, then compare with the solutions below.

Why does the null model omit the currently tested SNP while still allowing a genomic relationship matrix?
A relationship matrix has the correct number of rows and columns, but its individual order differs from the phenotype order. Is it ready to use? Explain.
In the model $V=\sigma_g^2K+\sigma_e^2I_n$, what conditions on $K$ and the variance components ensure positive definiteness of $V$?
Two optimizer runs return nearly identical variance components, likelihoods, and small gradients. What does this support, and what does it not establish?

Solutions. (1) The null model represents the trait without the variant-specific fixed effect, while $K$ represents genome-wide relatedness covariance; these are distinct model components. (2) No. Reorder or subset the matrix so every row and column corresponds to the same individual in the same order as the phenotype and covariates; dimensions alone do not establish correct pairing. (3) $K$ must be symmetric positive semidefinite, $\sigma_g^2\geq0$, and $\sigma_e^2>0$. Then $\sigma_g^2K$ is positive semidefinite and adding $\sigma_e^2I_n$ makes $V$ positive definite. (4) The agreement supports numerical stability across those starts. It does not prove that the model assumptions, covariate choices, or relationship structure are scientifically correct.