Handling Missing Genotype Calls in Association Regression

Why Missing Genotype Calls Need Special Handling

A missing variant call is not a genotype of zero, and it is not necessarily an uncertain call. It means that the genotype was not recorded or did not pass the calling procedure for that sample at that variant. The distinction matters in association regression: dropping samples, substituting a predicted dosage, and integrating over possible genotype states use different information and rely on different assumptions.

Consider a genetic association model relating an outcome \(Y_i\) to the number of effect alleles \(G_i\), possibly with covariates \(X_i\). If some \(G_i\) values are unavailable, the analysis must address both the unobserved genotype and the process that made it unavailable. A high overall call rate is reassuring about data quantity, but it does not by itself show that the observed calls represent the intended analysis population.

Definition: Let \(G_i\in\{0,1,2\}\) denote the number of effect alleles for individual \(i\), and let \(R_i=1\) if the genotype call is observed and \(R_i=0\) if it is missing. The call rate is the proportion of individuals with \(R_i=1\). A genotype dosage is an expected allele count, commonly \(D_i=E(G_i\mid Z_i)\), where \(Z_i\) denotes genetic information available for prediction, such as nearby typed variants.

A genotype can be missing because of low-quality reads, assay failure, or filtering. An observed genotype can instead be uncertain: sequencing or imputation information may assign probabilities to several possible values of \(G_i\). Dosage and likelihood methods can use that uncertainty, whereas complete-call analysis typically treats an accepted call as known and excludes a missing call.

Complete-Call Analysis

Complete-call analysis fits the chosen association regression only to individuals with observed genotypes, usually also requiring observed outcome and covariates. It is simple and transparent, and it avoids inventing genotype values. Its limitation is that it changes the analysis sample. Whether this changes the association estimate depends on how call availability relates to the outcome, genotype, and covariates, as well as on the regression model.

For example, if the conditional outcome model is correctly specified and call availability is independent of the outcome given the genotype and modeled covariates, selection into the called sample need not change the conditional outcome relationship. The analysis also needs adequate representation of relevant genotype and covariate combinations among called individuals. If call availability depends on the outcome even after conditioning on genotype and covariates, the outcome distribution among called individuals may differ from the target distribution. A regression fitted only to those individuals can then be biased.

Worked Example: A Complete-Call Linear Regression

Suppose a hypothetical quantitative-trait study has six participants. The trait is measured in units, and the genotype is coded as the number of effect alleles. Calls are available for four participants; the other two have missing calls.

Participant	Genotype \(G_i\)	Trait \(Y_i\)	Call observed?
1	0	10	Yes
2	1	12	Yes
3	2	14	Yes
4	1	12	Yes
5	Missing	16	No
6	Missing	15	No

For the four called participants, \(\bar G=1\) and \(\bar Y=12\). The simple-regression slope is

\hat\beta_G= \frac{\sum_{i:R_i=1}(G_i-\bar G)(Y_i-\bar Y)} {\sum_{i:R_i=1}(G_i-\bar G)^2} = \frac{(-1)(-2)+(0)(0)+(1)(2)+(0)(0)} {(-1)^2+0^2+1^2+0^2} =\frac{4}{2}=2.

The fitted slope is 2 trait units per additional effect allele among the called participants. The observed called genotypes include all three genotype levels, so the denominator is nonzero and the slope is estimable in this small example. The missing participants have higher observed trait values than the called participants, which is a reason to examine call availability by phenotype. Their unobserved genotypes prevent us from determining from this table alone whether complete-call analysis recovers the association in the full six-person group.

Dosage-Based Analysis

A dosage replaces a categorical genotype with its expected allele count using available genotype-prediction information. For example, if the predicted probabilities of \(G_i=0,1,2\) are \(q_{i0},q_{i1},q_{i2}\), then the dosage is \(D_i=q_{i1}+2q_{i2}\). In practice, these probabilities may come from imputation using nearby markers or from other genetic information. When a call is wholly absent and no informative auxiliary data are available, the individual-specific dosage may be little more than an expected allele count based on a population model.

One can use \(D_i\) in place of \(G_i\) in an association regression, but this is a plug-in procedure: it treats the expected allele count as if it were the observed predictor. It is not identical to averaging the outcome model over possible genotypes. This distinction is especially important for a nonlinear model such as logistic regression. Dosage analysis is often useful and computationally convenient, but its interpretation and accuracy depend on the quality and calibration of the genotype predictions.

Worked Example: A Dosage Is Not a Weighted Outcome Probability

Suppose a hypothetical participant has genotype probabilities \(P(G=0\mid Z)=0.15\), \(P(G=1\mid Z)=0.60\), and \(P(G=2\mid Z)=0.25\), based on genetic information \(Z\). The dosage is

$$ D=0(0.15)+1(0.60)+2(0.25)=1.10. $$

For illustration, suppose the binary-outcome model is \(\operatorname{logit}P(Y=1\mid G)= -2.197+0.6931G\). Plugging the dosage into the model gives

\operatorname{logit}^{-1}(-2.197+0.6931(1.10)) =\operatorname{logit}^{-1}(-1.43459) =\frac{1}{1+\exp(1.43459)} \approx 0.1924.

If instead we average the outcome probability over the three genotype states, the calculation is

\begin{aligned} P(Y=1\mid Z) &=\sum_{g=0}^{2}P(Y=1\mid G=g)P(G=g\mid Z)\\ &\approx (0.1000)(0.15)+(0.1818)(0.60)+(0.3077)(0.25)\\ &\approx 0.2010. \end{aligned}

The two results differ: \(0.1924\) is the fitted probability at the mean dosage, while \(0.2010\) averages the fitted probabilities across genotype states. The difference arises because the inverse-logit function is nonlinear. Thus replacing \(G\) with \(D\) is not generally the same operation as accounting for genotype uncertainty in the outcome likelihood.

A common but weak substitute is to set every missing genotype to the same sample mean or allele-frequency-based mean. This removes individual-level information about genotype variation and can reduce the apparent genotype contrast. A single imputed value also hides uncertainty: two individuals assigned the same expected dosage may have very different genotype-probability distributions. If dosages are used, the source of the predictions and their uncertainty should be considered, rather than treating a convenient numeric replacement as a known genotype.

Likelihood-Based Treatment of Genotype Uncertainty

A likelihood-based approach can account for uncertainty by summing the outcome likelihood over possible genotype values, weighted by their probabilities given the available genetic information. Let \(f(Y_i\mid G_i=g,X_i;\theta)\) be the outcome model, with parameters \(\theta\). If genotype-prediction information \(Z_i\) gives probabilities \(q_{ig}=P(G_i=g\mid Z_i,X_i)\), then the contribution for an individual with observed outcome and uncertain genotype is

L_i(\theta)=\sum_{g=0}^{2} f(Y_i\mid G_i=g,X_i;\theta)\,q_{ig}.

For binary \(Y_i\), the outcome likelihood in this expression is the fitted probability of the observed outcome: \(P(Y_i=1\mid G_i=g,X_i)\) for a case and \(1-P(Y_i=1\mid G_i=g,X_i)\) for a non-case. When the genotype is known with certainty, its probability distribution puts all weight on the called genotype, and the sum reduces to the usual regression contribution. When the genotype is missing and there is no individual-level genetic information, the weights must come from an appropriate model for genotype given available covariates or population information; they cannot be learned from nothing.

More generally, sequencing data may provide read likelihoods \(P(\text{reads}_i\mid G_i=g)\), rather than genotype probabilities. A joint likelihood can combine those read likelihoods, a model for genotype frequencies, and the outcome regression. Whichever formulation is used, the genotype model and outcome model must be suitable for the data. If a missingness or calling process is informative even after conditioning on the variables included in the likelihood, it may also need to be modeled or addressed through a justified analysis strategy.

Worked Example: A Likelihood Contribution for an Uncertain Call

Use the hypothetical genotype probabilities and logistic model from the dosage example. For a participant whose outcome is \(Y=1\), the contribution that averages over genotype uncertainty is

\begin{aligned} L(\theta) &=\sum_{g=0}^{2}P(Y=1\mid G=g)P(G=g\mid Z)\\ &=(0.1000)(0.15)+(0.1818)(0.60)+(0.3077)(0.25)\\ &\approx 0.2010. \end{aligned}

For a participant with \(Y=0\) and the same genotype probabilities, the contribution is the weighted probability of a non-case:

\begin{aligned} L(\theta) &=(1-0.1000)(0.15)+(1-0.1818)(0.60)+(1-0.3077)(0.25)\\ &=(0.9000)(0.15)+(0.8182)(0.60)+(0.6923)(0.25)\\ &\approx 0.7990. \end{aligned}

If these two participants are independent conditional on the model inputs, their joint likelihood contribution is approximately \(0.2010\times0.7990=0.1606\). This calculation retains the uncertainty over genotype states instead of assigning both participants the same fixed genotype. In a real analysis, the genotype probabilities and regression parameters would be estimated or supplied by an appropriate model; these values are hypothetical and illustrate the calculation only.

How Missingness Can Threaten Validity

The call indicator \(R_i\) should be examined in relation to the observed phenotype and covariates. For example, call failure that is more common in cases than controls may signal differential data quality, batch effects, or selection related to phenotype. Differences in call rate across ancestry groups, recruitment sites, or other covariates can also change which individuals contribute to a complete-call analysis. If such variables are associated with genotype or outcome, a change in the analysis sample can matter even when the overall call rate is high.

Missing-data terminology clarifies what can and cannot be concluded. Missingness is missing at random when the probability that a genotype is missing can depend on observed information, but not on the missing genotype after conditioning on that observed information. It is missing not at random when call availability still depends on the unobserved genotype after conditioning on observed information. In practice, the missing genotype cannot be directly compared between called and uncalled individuals, so a missing-at-random assumption is not proved simply by finding no observed differences.

Call-rate summaries by phenotype and covariate strata are useful checks, not automatic cures. A regression of \(R_i\) on observed \(Y_i\) and \(X_i\), or a table of call rates across those groups, can reveal patterns to investigate. These checks cannot establish whether missingness depends on the unobserved \(G_i\). Nor does a high call rate guarantee valid inference: a small amount of missingness concentrated in a particular outcome or covariate group may be more consequential than a larger amount missing uniformly.

Practical rule: Choose the approach to match the available information. Use complete calls when exclusion is defensible for the target association and the called sample has adequate support; use well-calibrated dosages when genotype predictions are available and a plug-in analysis is appropriate; consider likelihood-based integration when genotype uncertainty can be represented by credible probabilities or read likelihoods. State what assumptions link the observed data to the intended association.

A Short Justification for Likelihood Integration

The likelihood sum follows from the law of total probability. The genotype state is latent, so the probability of the observed outcome and genetic information is obtained by considering each possible genotype and adding the corresponding joint probabilities. Under the factorization \(P(Y_i,G_i\mid Z_i,X_i)=P(Y_i\mid G_i,X_i)P(G_i\mid Z_i,X_i)\), this gives

\begin{aligned} P(Y_i\mid Z_i,X_i) &=\sum_{g=0}^{2}P(Y_i,G_i=g\mid Z_i,X_i)\\ &=\sum_{g=0}^{2}P(Y_i\mid G_i=g,X_i)P(G_i=g\mid Z_i,X_i). \end{aligned}

The factorization states a modeling assumption: conditional on genotype and covariates, the outcome model does not require additional dependence on \(Z_i\). If that assumption is unsuitable, the joint model must be expanded. The derivation also shows why the likelihood contribution is a weighted sum of probabilities, rather than the outcome probability evaluated at the weighted-average genotype.

Common Mistakes

Treating a missing call as genotype 0. Zero denotes a known genotype with no effect alleles; missing means the genotype is unknown.
Assuming a high call rate settles the issue. The pattern of missingness matters, including its relation to phenotype and covariates.
Calling a dosage a known genotype. A dosage is an expected allele count and can represent substantial uncertainty about the underlying genotype.
Equating plug-in dosage with likelihood integration. In nonlinear models, evaluating the outcome function at the mean genotype generally differs from averaging outcome probabilities over genotypes.
Assuming complete-case analysis is always unbiased. It changes the analysis sample; its validity depends on the selection pattern, model, and support in the called sample.

Check Your Understanding

Use these questions to check whether you can distinguish the methods and their assumptions.

What is the difference between a missing genotype call and an observed but uncertain genotype?
Under what kind of selection pattern can complete-call regression fail to estimate the intended conditional association?
If genotype probabilities are 0.50 for \(G=0\), 0.30 for \(G=1\), and 0.20 for \(G=2\), calculate the dosage.
Why does evaluating a logistic model at a dosage generally differ from averaging its predicted probabilities over genotype states?
Name two observed variables against which call rates should be examined, and explain why those checks cannot establish that missingness is independent of the unobserved genotype.