Why Analysis Order Matters
At a genetic locus, several correlated variants may each show an association with a phenotype when considered one at a time. A joint regression asks a different question: what is the association estimate for each included variant when the other included variants are held fixed? In a fixed model, changing the order in which the genotype columns are entered does not change those estimates. A sequential procedure, however, repeatedly changes the conditioning set as it adds, tests, or stops at variants. Its path—and sometimes its final selected set—can depend on the order of analysis.
These are distinct meanings of “order.” Reordering columns in one simultaneous regression is a relabeling operation. Changing which variant is conditioned on first in a sequential procedure changes the intermediate model and may change the next test. The distinction matters when interpreting conditional association results: an order-dependent selection path is not evidence that the final coefficients of a fixed, full model are themselves order-dependent.
Fixed Joint Models and Column Order
Suppose a quantitative phenotype \(Y\) is modeled using an intercept, a fixed set of covariates, and \(m\) genotype variables \(G_1,\ldots,G_m\). With additive dosage coding, a genotype variable records the number of copies of a specified allele, usually 0, 1, or 2. The linear model is
Here \(\beta_j\) is the conditional association coefficient for variant \(j\), holding the other included genotypes and covariates fixed. Its units are phenotype units per additional allele copy when the dosage is unstandardized. For ordinary least-squares estimates and their usual exact \(t\)-tests, assume a correctly specified linear conditional mean, independent errors with mean zero and constant variance, and normally distributed errors for exact finite-sample \(t\)-inference. The coefficient-order result below is an algebraic property of the fitted design matrix; it does not itself require normal errors.
Let \(X\) be the design matrix, including the intercept and all covariates and genotypes, and let \(\mathbf{y}\) be the outcome vector. When \(X\) has full column rank, the least-squares estimate is \(\hat{\boldsymbol{\beta}}=(X^{\mathsf T}X)^{-1}X^{\mathsf T}\mathbf{y}\). If \(P\) is a permutation matrix that reorders columns, the reordered design is \(X^*=XP\). The coefficient vector must be reordered in the opposite way to keep the same fitted values.
The equality uses \(P^{\mathsf T}P=I\), the defining property of a permutation matrix. Thus the coefficient attached to each named variant, its fitted contribution, and the model’s fitted values are unchanged. Only the positions of the coefficients in the printed vector change. Standard errors and tests for the same named coefficients in the same model are unchanged as well, apart from numerical rounding.
Worked Example: The Same Joint Estimates After Reordering
Consider a hypothetical sample with centered genotype columns \(G_A\) and \(G_B\), and a centered quantitative phenotype \(Y\). Suppose the cross-products are \(G_A^{\mathsf T}G_A=40\), \(G_B^{\mathsf T}G_B=30\), \(G_A^{\mathsf T}G_B=12\), \(G_A^{\mathsf T}Y=16\), and \(G_B^{\mathsf T}Y=9\). The model includes both genotypes simultaneously. Since the variables are centered, the intercept is zero in this centered calculation.
The normal equations for the two genotype coefficients are
The determinant is \(40(30)-12(12)=1056\), so the inverse cross-product matrix is \(1056^{-1}\begin{pmatrix}30&-12\\-12&40\end{pmatrix}\). Therefore,
If the columns are entered as \(G_B,G_A\), the normal equations and coefficient positions are permuted, but the estimate for \(G_A\) remains \(0.3523\) and the estimate for \(G_B\) remains \(0.1591\). If \(Y\) is measured in a particular unit and the genotypes are unstandardized allele dosages, these coefficients are in that phenotype unit per allele copy. This calculation assumes the same individuals and the same included variables in both fits.
Sequential Conditioning Changes the Model Along the Way
A sequential conditioning procedure does not merely rearrange a fixed design matrix. It begins with a conditioning set, tests or scans other variants conditional on that set, and may then enlarge the set. At each step, the question is conditional on the variables already included. If the sequence changes, the intermediate conditioning sets may change too.
For a candidate genotype \(G\) tested after a set of variables \(S\), the incremental test asks whether the coefficient of \(G\) is zero conditional on \(S\). If the order changes \(S\), it changes the null hypothesis being tested and the residual variation used in that test. Consequently, even with identical data, different sequences can pass or fail different tests at a chosen threshold.
A crucial qualification is that order dependence is not inevitable. If a procedure eventually includes exactly the same predictors and fits the same ordinary least-squares model to the same observations, its final joint estimates are invariant to entry order. Order can affect which predictors are selected on the way to that fit, especially when the method stops early or applies a threshold at each step.
Worked Example: A Sequential Stop Can Depend on the Starting Variant
Suppose a hypothetical sample contains \(n=40\) unrelated individuals with complete data. Let \(Y\) be a standardized quantitative phenotype and \(G_A,G_B\) standardized genotype dosages. Assume their sample correlations are \(r_{YA}=0.50\), \(r_{YB}=0.40\), and \(r_{AB}=0.30\). These correlations are jointly feasible: their correlation matrix has determinant \(1+2(0.50)(0.40)(0.30)-0.50^2-0.40^2-0.30^2=0.62>0\). Consider the rule “start with the specified variant, test the other variant conditional on it, and include the second variant only if its two-sided test has \(p<0.05\).” For the usual exact linear-model tests, assume independent, homoskedastic, normally distributed errors conditional on the genotypes.
For two standardized predictors, the partial correlation of \(Y\) with \(G_B\) given \(G_A\) is
With one variable controlled, the partial-correlation test has \(df=n-3=37\). Its statistic is \(t=r\sqrt{df/(1-r^2)}\), giving
The two-sided \(p\)-value is approximately \(0.061\), so \(G_B\) is not added at the \(0.05\) threshold when the sequence starts with \(G_A\). Reversing the starting variant gives
This test has a two-sided \(p\)-value of approximately \(0.0057\), so \(G_A\) is added when the sequence starts with \(G_B\). Under the stated stop rule, the first sequence ends with \(G_A\) alone, whereas the reversed sequence includes both variants. This example does not say that a fixed two-variant regression has different estimates in the two orders; it shows that the sequential decision path can end in different models.
Joint Estimates and Sequential Decisions Answer Different Questions
In a simultaneous model, each coefficient answers a conditional association question for a specified set of included predictors. In a sequential analysis, a test at a given step asks whether a candidate contributes evidence beyond the current conditioning set. Those are not interchangeable questions. A marginal association, a conditional test after one selected variant, and a coefficient from a larger joint model need not have the same value or interpretation.
The distinction is especially important when communicating a locus analysis. “Variant \(B\) was not selected after conditioning on \(A\)” describes a method’s result under a specific sequence and threshold. It is not equivalent to saying that \(B\)'s coefficient in every possible joint model is zero. Likewise, changing the candidate order does not establish that one variant has a biological effect and another does not. These regressions describe statistical associations under their coding and model assumptions; association alone does not establish causation.
Worked Example: Translating a Joint Estimate Into Context
Return to the hypothetical centered cross-products in the first example. Suppose the phenotype is measured in millimeters of a laboratory measurement and each genotype column is allele dosage. The simultaneous estimates were \(\hat{\beta}_A=0.3523\) and \(\hat{\beta}_B=0.1591\) millimeters per allele copy, conditional on the other included genotype. If the analyst first enters \(G_B\) and then \(G_A\), the same two-variant fit must still report \(0.3523\) for \(G_A\) and \(0.1591\) for \(G_B\).
If instead an analysis stops after fitting \(G_B\) alone, it has not estimated \(\beta_B\) from the two-variant model. It has estimated the association in a different, one-genotype model. Adding \(G_A\) changes the conditioning set, so the coefficient being estimated is different by definition. The correct report therefore names the model: for example, “the joint dosage coefficient for \(G_B\), adjusted for \(G_A\), was \(0.1591\) mm per allele copy,” rather than describing that number as an order-specific estimate.
A Practical Way to Keep the Two Ideas Separate
When reviewing locus results, first identify whether the analysis is a fixed joint fit or a sequence of model-building steps. For a joint fit, record the complete predictor set and verify that the same individuals, phenotype, coding, and covariates are used before comparing results under a different column order. For sequential analysis, record the starting variant, each conditioning set, the test used at each step, the entry threshold, and the stopping rule.
Write down the phenotype, genotype coding, covariates, and included variants for a simultaneous regression. For a sequential procedure, state the starting set and its selection rule.
A coefficient in a joint model is conditional on all the other included predictors. An incremental test is conditional on the predictors present at that particular step.
If a fixed model gives different coefficients merely because columns were reordered, check whether the model, data, coding, or fitting procedure actually changed. If the selected set differs, inspect the sequence and stopping rule before interpreting it.
Common Mistakes
- Calling joint coefficients order-dependent. Reordering the columns of the same full-rank design matrix only reorders the coefficient vector. The estimate attached to each named variant remains the same.
- Assuming a sequential path is just a reordered joint fit. Each sequential step may use a different conditioning set and test a different conditional null hypothesis. Its selection decisions can therefore depend on the path.
- Comparing numbers without naming their models. A coefficient from a one-variant model and one from a two-variant model answer different conditional questions. State the included variants and units.
- Treating failure to enter as proof of no association. A variant that fails a sequential threshold did not meet that procedure’s rule at that step. This alone does not prove its coefficient is zero in a prespecified joint model.
- Reading a conditional association as a causal effect. The regression coefficient is an association under the stated model and assumptions; it does not by itself identify a causal genetic effect.
Check Your Understanding
Answer each question by distinguishing a fixed joint model from a sequential procedure.
- Two genotype columns are swapped in an ordinary least-squares fit, with no other changes. What happens to the fitted values and to the coefficient attached to each named variant?
- Why can a sequential procedure’s test for a candidate variant change when the preceding conditioning set changes?
- A variant is not added because its conditional \(p\)-value exceeds the procedure’s threshold. What conclusion is justified, and what stronger conclusion is not justified?
- What conditions should be held fixed before attributing a difference in joint coefficients to anything other than numerical rounding?
- In a report, what information should accompany a sequentially selected signal so readers can interpret how it was obtained?