Simulation-Based Power For Epistasis Screens

Why Simulate Power for an Interaction Screen?

For a prespecified locus pair, a power calculation can use a particular model and its assumptions to predict how often a test will detect an interaction. A simulation makes those assumptions explicit by generating genotypes, phenotypes, and missing calls, then applying the planned analysis to each generated data set. This is useful when linkage disequilibrium, missingness, or the phenotype model makes a simple approximation inadequate.

The result is conditional on the scenario being simulated. It is not a guarantee that a screen will discover a biological interaction: it estimates the probability that the specified analysis rejects its null under the specified data-generating model. The test and rejection threshold must match the planned analysis. In a screen, that threshold should reflect the intended testing procedure; the simulation here does not determine how to control error across variant pairs.

Specify the Genetic and Phenotype Models

First specify allele frequencies and linkage disequilibrium (LD) for two biallelic loci, A and B. Let \(p_A\) and \(p_B\) be the frequencies of the alleles coded 1, and let \(D\) describe the difference between the haplotype frequency and the frequency expected under independence. The four haplotype probabilities are:

These probabilities must each be between 0 and 1 and sum to 1. Equivalently, \(D\) must satisfy \(\max\{-p_Ap_B,-(1-p_A)(1-p_B)\}\leq D\leq\min\{p_A(1-p_B),(1-p_A)p_B\}\). To generate one diploid individual under random mating, draw two haplotypes independently from these probabilities and count the allele coded 1 at each locus. The resulting dosages \(G_A\) and \(G_B\) are each 0, 1, or 2. This construction gives the specified LD and, under the random-union-of-haplotypes assumption, Hardy–Weinberg genotype proportions at each locus.

Next specify how the phenotype depends on those dosages. For an additive-by-additive interaction, one common quantitative-trait model is a linear mean with product term. For a binary trait, a corresponding model uses the logit of disease probability. These models define the interaction on different scales, so their coefficients and power results are not interchangeable.

For the quantitative model, this specifies normally distributed residual variation with constant variance and mean zero, independent of genotype. For the binary model, it specifies a Bernoulli outcome with the probability given by the logistic equation. In either case, the interaction alternative has \(\beta_{AB}\ne0\), and the interaction null is \(H_0:\beta_{AB}=0\). The main effects remain in both models: this is a hierarchical model. Covariates, if part of the planned analysis, should also be generated and included in both the data-generating and fitted models as appropriate.

Represent Missing Genotypes and the Planned Test

Missingness must be part of the design, not an afterthought. For example, suppose each locus has an independent probability \(q_A\) or \(q_B\) of a missing genotype call, and those events are independent of genotype and phenotype. If the planned analysis requires both calls, an individual is retained only when both are observed. The probability of retaining a generated individual is then \((1-q_A)(1-q_B)\), and the expected number retained from \(n\) generated individuals is \(n(1-q_A)(1-q_B)\). This expectation is not a fixed realized sample size: the number retained varies across simulations.

That calculation depends on the stated missingness assumptions. If missingness depends on genotype, phenotype, or other variables, simulate that dependence explicitly and use an analysis that is appropriate for the resulting missing-data process. Applying complete-case analysis under an unexamined missingness mechanism can change the analyzed population or bias estimates; simply reducing the sample size by an average percentage does not represent all consequences of informative missingness.

For an additive dosage interaction, fit the same hierarchical model in every replicate and test \(H_0:\beta_{AB}=0\). For example, a likelihood-ratio test compares the model containing the product term with the otherwise identical model omitting it. The test's degrees of freedom are 1 for this single product coefficient, provided the model is identifiable and fits successfully. The significance level \(\alpha\) must be chosen to match the planned screen analysis. A single-pair threshold and a screen-wide threshold can produce different power, even with identical simulated data.

Estimate Power and Its Monte Carlo Error

Generate \(B\) independent data sets under the alternative. In each one, generate genotypes, phenotypes, and missingness; apply the planned inclusion and analysis rules; fit the model; and record whether the interaction test rejects at the chosen threshold. If \(R_b=1\) when replicate \(b\) rejects and \(R_b=0\) otherwise, the estimated power is the mean of these indicators.

The Monte Carlo standard error (MCSE) describes simulation noise in the estimated rejection probability. It does not capture uncertainty about whether the assumed allele frequencies, effect sizes, or missingness rates are realistic. Conditional on a fixed scenario and a valid analysis, each rejection indicator is Bernoulli with probability equal to power. Independent replicates therefore have variance \(\text{Power}(1-\text{Power})\), and the sample mean has variance \(\text{Power}(1-\text{Power})/B\). Replacing the unknown power by its estimate and taking the square root gives the displayed MCSE.

A Binary-Trait Simulation and Scenario Comparisons

Binary outcomes require generating event probabilities from the fitted scale's model, then drawing each outcome from a Bernoulli distribution. The baseline coefficient should be chosen to represent a plausible baseline risk at the reference genotype, and all effect coefficients must be interpreted on the log-odds scale in a logistic model. A simulation should also track model failures: silently counting nonconverged fits as non-rejections changes the estimated performance of the planned procedure.

Validation, Interpretation, and Common Pitfalls

A simulation is only as trustworthy as its data generation and analysis. Before estimating power, check that haplotype probabilities are nonnegative and sum to one; simulated allele frequencies and LD are near their targets; and the generated phenotype distribution and missing-call rates match their specifications. Confirm that the fitted model includes the intended main effects and interaction, and that the test and threshold match the planned analysis. Report the number of replicates, estimated power, MCSE, and the assumptions that define the scenario.

It is also useful to simulate under the null, setting \(\beta_{AB}=0\), and measure the rejection rate. Under a correctly calibrated test and a sufficiently large number of replicates, this rate should be compatible with the intended \(\alpha\), allowing for Monte Carlo variation. A substantial discrepancy can reveal model mismatch, a coding error, convergence problems, or inappropriate use of the test. For small or sparse genotype cells, some replicates may be unidentifiable; report their frequency and decide in advance how the planned analysis handles them.

Power is scenario-specific. Rather than relying on one set of assumed values, a useful design study can vary plausible allele frequencies, LD, interaction effects, residual variation or baseline disease risk, and missingness rates. This reveals which assumptions most influence detection and makes clear what the estimated power does—and does not—support. Because simulation error decreases only at the rate \(1/\sqrt{B}\), more replicates are especially useful when comparing scenarios whose estimated powers are close.