Diagnostics And Reporting For Epistasis Studies

Why Diagnostics Matter After an Interaction Test

An interaction test answers a specific model-based question: whether the fitted interaction term or terms are compatible with a null value under the chosen model. A small p-value does not show that the model describes the data adequately, that all genotype combinations are represented well, or that the result is insensitive to analytic choices. Diagnostics and transparent reporting help readers judge those issues without confusing them with the test itself.

For a two-locus analysis, the essential checks concern the outcome model, the distribution of people across joint genotype categories, the relationship between the loci, and the way genotypes were coded. These checks do not replace the interaction test or its multiple-testing adjustment. They reveal limitations that can affect the estimate, its uncertainty, or its interpretation.

Definition: A diagnostic is a check of whether features of the observed data and fitted model are compatible with the model’s assumptions and intended use. A reporting check asks whether the analysis is described sufficiently for a reader to identify what was tested, how it was tested, and what evidence supports the conclusion.

Check the Outcome Model

For a quantitative trait, a common model includes the two locus terms, their interaction, and any prespecified covariates. In a simple linear-model notation, \(Y_i\) is the trait for person \(i\), \(G_{1i}\) and \(G_{2i}\) are the chosen genotype codes, and \(C_i\) denotes covariates:

Y_i=\beta_0+\beta_1G_{1i}+\beta_2G_{2i}+\beta_{12}G_{1i}G_{2i}+C_i^\mathsf{T}\gamma+\varepsilon_i.

The usual linear-model conditions for the stated inference include a correctly specified conditional mean, independent errors (or an appropriate model for dependence), and, for conventional small-sample t-based inference, an adequate error-distribution approximation. Constant error variance is also relevant to conventional standard errors; robust or otherwise appropriate standard errors may address some forms of heteroscedasticity, but they do not repair a misspecified mean model. Residual-versus-fitted plots can reveal curvature or changing spread, while residual distributions and influence checks can flag unusual observations. No single plot establishes that every assumption holds.

For a binary trait, a logistic model commonly specifies the log odds as a function of the two genotype terms, their interaction, and covariates. The model assumes the stated log-odds structure is adequate and that the observations are independent conditional on the model, unless clustering or other dependence is explicitly handled. Logistic residuals are not expected to look like normally distributed linear-model residuals. Instead, inspect appropriate residual or influence summaries and assess whether predicted risks agree reasonably with observed outcomes across relevant ranges.

Worked Example: Reading a Quantitative-Trait Residual Check

Suppose a hypothetical study fits a linear interaction model for a quantitative trait measured in trait units. For one participant, the fitted value is 10.2 units and the observed value is 11.0 units. The raw residual is:

e_i=Y_i-\hat{Y}_i=11.0-10.2=0.8\text{ units}.

A positive residual means this participant’s observed trait is 0.8 units above the model prediction. It does not by itself indicate a problem. In this hypothetical fit, a residual-versus-fitted plot shows residual spread increasing as fitted values rise, and a residual-versus-\(G_1G_2\) plot shows a systematic pattern across interaction values. Those patterns suggest checking the mean and variance specification rather than declaring the interaction real or spurious from the plots alone. The analysts should also inspect influential observations and use standard errors appropriate to the variance structure if needed. A diagnostic pattern is a reason for investigation, not automatic grounds for deleting observations.

Calibration for Binary Outcomes

Calibration concerns agreement between predicted risks and observed outcome frequencies. If a model assigns risk \(\hat{q}_i\) to participant \(i\), then among participants with similar predictions the average observed outcome proportion should be reasonably close to the average predicted risk. Calibration is distinct from discrimination: a model can rank higher-risk people above lower-risk people while systematically overpredicting or underpredicting their risks.

A practical summary groups predictions into bins and compares the mean predicted risk with the observed proportion in each bin. A more formal calibration assessment can fit a logistic regression of the outcome on the model’s predicted log odds; ideal calibration corresponds to intercept 0 and slope 1. Such estimates are uncertain, especially in small samples, and in-sample checks can look overly favorable because the same data were used to fit the model. Validation data or resampling can help assess out-of-sample calibration.

Worked Example: Comparing Predicted and Observed Risks

In a hypothetical binary-trait analysis, participants are divided into three prediction groups. The first group contains 40 people with mean predicted risk 0.10, of whom 3 have the outcome. The second contains 40 people with mean predicted risk 0.30, of whom 14 have the outcome. The third contains 20 people with mean predicted risk 0.60, of whom 11 have the outcome.

For the first group, the observed proportion is \(3/40=0.075\), compared with predicted risk 0.10. For the second, it is \(14/40=0.350\), compared with 0.30. For the third, it is \(11/20=0.550\), compared with 0.60. The observed number of outcomes is \(3+14+11=28\). The total expected number, using the group mean risks, is \(40(0.10)+40(0.30)+20(0.60)=4+12+12=28\).

The total expected and observed counts agree, but the group comparisons show departures in both directions: the second group has two more observed outcomes than its expected count of 12, while the other two groups each have one fewer than expected. Equality of the totals therefore does not demonstrate perfect calibration. These summaries are descriptive; they do not by themselves establish that the differences are statistically meaningful.

Joint Genotype Cells and Linkage Disequilibrium

For two biallelic loci coded by allele dosages 0, 1, and 2, a joint genotype table has up to nine cells. Report the number of participants in each cell, not only the separate genotype counts at each locus. Marginal counts can conceal a joint cell with very few observations or none. Sparse or empty cells can make some interaction contrasts unstable or non-estimable, particularly in logistic models, where separation can also occur. There is no universal cell-count threshold that guarantees reliable estimation; the implications depend on the model, outcome frequency, and contrast being estimated.

Linkage disequilibrium (LD) describes non-random association of alleles at different loci. Let \(p_A\) and \(p_B\) be the frequencies of specified alleles, and let \(p_{AB}\) be the frequency of chromosomes carrying both. Define \(D=p_{AB}-p_Ap_B\). A common standardized measure is \(r^2\), which compares squared allelic covariance with the product of allele variances.

Definition: For two polymorphic loci, so that \(0<p_A<1\) and \(0<p_B<1\), linkage disequilibrium is \(D=p_{AB}-p_Ap_B\), and \(r^2=D^2/\{p_A(1-p_A)p_B(1-p_B)\}\). Under these conditions \(0\leq r^2\leq1\). If either allele frequency is 0 or 1, the denominator is zero and \(r^2\) is undefined.

The bound follows by considering chromosome-level indicators \(X_A\) and \(X_B\), each equal to 1 when the specified allele is present and 0 otherwise. Their covariance is \(D\), and their variances are \(p_A(1-p_A)\) and \(p_B(1-p_B)\). The squared correlation is therefore \(r^2\). By the Cauchy–Schwarz inequality, the absolute correlation is at most 1, so its square is between 0 and 1. This also explains why both allele variances must be positive for the measure to be defined.

Worked Example: Joint Cell Counts and LD

Consider a hypothetical sample of 200 people. The two loci are coded as allele dosages 0, 1, and 2. The joint counts are:

Locus 1 dosage	Locus 2 dosage 0	Locus 2 dosage 1	Locus 2 dosage 2	Row total
0	40	10	0	50
1	10	80	10	100
2	0	10	40	50
Column total	50	100	50	200

The margins alone appear balanced, but the table reveals two empty cells and several cells with only 10 people. Those counts warrant scrutiny of the interaction contrasts that rely on these combinations; the table does not establish that every possible model is unestimable.

Suppose phased chromosome data from these same hypothetical participants give haplotype counts \(00=140\), \(01=60\), \(10=60\), and \(11=140\), out of 400 chromosomes. Then \(p_A=(60+140)/400=0.50\), \(p_B=(60+140)/400=0.50\), and \(p_{AB}=140/400=0.35\). Thus \(D=0.35-(0.50)(0.50)=0.10\), and:

r^2=\frac{0.10^2}{(0.50)(0.50)(0.50)(0.50)} =\frac{0.01}{0.0625}=0.16.

This is a defined \(r^2\), because both loci are polymorphic, and it indicates nonzero LD in this hypothetical sample. LD diagnostics help readers assess dependence between the loci and potential aliasing in fitted terms; they do not, by themselves, determine whether an interaction is biologically meaningful.

Check Sensitivity to Genotype Coding

An interaction coefficient depends on how the genotypes are represented. An additive dosage code assigns 0, 1, and 2 copies of an allele. Dominant and recessive codes instead combine genotype categories, and a categorical coding can represent each genotype separately. These choices test different patterns of effect. Consequently, estimates and p-values from different coding schemes are not interchangeable evidence about one identical parameter.

Worked Example: Comparing Prespecified Coding Choices

Suppose a hypothetical study reports an additive-by-additive interaction estimate of 0.50 trait units per product of allele dosages, with \(SE=0.20\). Its Wald statistic is \(z=0.50/0.20=2.50\), corresponding to a two-sided normal-approximation p-value of about 0.0124. A 95% Wald interval is \(0.50\pm1.96(0.20)\), or approximately 0.108 to 0.892 trait units per dosage-product unit.

Suppose prespecified dominant coding gives an interaction estimate of 0.18 with \(SE=0.15\), so \(z=1.20\) and the two-sided p-value is about 0.230. A recessive coding gives 0.62 with \(SE=0.31\), so \(z=2.00\) and the two-sided p-value is about 0.0455. These are hypothetical model fits, not recalculations from the first coefficient: the coding changes the contrast being estimated. Report the coding and its corresponding interpretation for each analysis. If alternative codings were explored after seeing the results, disclose that exploration rather than presenting the most favorable fit as though it were the sole prespecified analysis.

What a Transparent Report Includes

A useful report lets a reader reconstruct the analysis and distinguish a prespecified test from exploratory follow-up. State the variants and alleles, the genotype coding, the outcome scale, the interaction parameter or joint test, and all covariates. Give the analysis sample size and explain exclusions or missing genotype handling. Summarize joint genotype counts and relevant LD estimates, including the method used to estimate haplotype frequencies if phase was not directly observed.

Report the interaction estimate with its standard error or confidence interval, the test statistic and p-value, and the multiplicity procedure used for the set of pairs searched. State how many pairs were tested and how pairs entered the analysis, including any screening or biological filtering. Identify the diagnostic checks performed and describe important findings, such as sparse cells, calibration departures, influential observations, or coding sensitivity. Avoid implying that a diagnostic plot proves model validity, or that an adjusted p-value removes limitations caused by sparse data or model misspecification.

Key takeaway: An interpretable epistasis result is more than an interaction p-value. Pair the test with model-appropriate diagnostics, joint genotype counts, LD information, explicit coding choices, and a clear account of the complete tested-pair set.

Common Mistakes

Equating a significant test with a well-fitting model. The interaction test is conditional on the fitted model. Check outcome-model diagnostics and report notable limitations.
Reporting only marginal genotype frequencies. Margins do not show whether particular two-locus combinations are rare or absent. Provide the joint-cell counts.
Computing \(r^2\) when a locus is monomorphic. If either allele frequency is 0 or 1, the denominator is zero. The standardized measure is undefined, not a value obtained by dividing by zero.
Calling coding sensitivity a contradiction without explanation. Additive, dominant, and recessive codings define different interaction contrasts. Name the contrast before comparing estimates.
Reporting only the best-looking pair or coding. Describe the search process, number of tested pairs, multiplicity handling, and any exploratory analyses.

Check Your Understanding

Answer each question by referring to the diagnostic or reporting issue it tests.

In a linear interaction model, what might a residual-versus-fitted pattern of increasing spread suggest, and why is it not by itself proof that the interaction is false?
Why can equal total predicted and observed outcome counts coexist with poor calibration in individual risk groups?
A joint genotype table contains an empty cell. What should be reported, and why is it not correct to claim automatically that every interaction model is unestimable?
What conditions make \(r^2\) defined, and what is its range under those conditions?
Why should an additive interaction estimate not be interpreted as the same contrast as a dominant-coded interaction estimate?
Name at least three elements that allow a reader to assess whether an epistasis result was reported transparently.