Why Missingness Can Look Like Population Structure
Principal component analysis (PCA) summarizes patterns of similarity among individuals. In genomic data, those patterns can reflect ancestry, but they can also reflect technical differences in which genotypes were called. If samples from one batch have systematically different call rates, or if missing genotypes are filled in using a value that does not match the observed genotypes, the processed data can acquire sample-level patterns that PCA represents as components. Those components are real features of the processed matrix; they need not be real ancestry differences.
The key issue is not simply that some entries are absent. It is how the missing entries are distributed across samples and variants, what values replace them, and whether missingness is related to the genotype values that PCA uses. The earlier material on adjusting for ancestry with principal components motivates interpreting PCs as summaries of genetic variation. Here the focus is a technical alternative explanation for apparent structure.
A Model for Missing and Processed Genotypes
Let \(G_{ij}\) denote the allele dosage for individual \(i\) at variant \(j\), coded as 0, 1, or 2, and let \(R_{ij}=1\) when that genotype is observed and \(R_{ij}=0\) when it is missing. Let \(\mu_j\) be the value subtracted from observed genotypes at variant \(j\), and let \(c_j\) be the value used to fill a missing entry on that same centered scale. The processed matrix entry is
For a batch \(b\), define \(r_{bj}=P(R_{ij}=1\mid b)\), the call rate at variant \(j\) in that batch, and define \(m_{bj}=E(G_{ij}-\mu_j\mid R_{ij}=1,b)\), the mean centered dosage among called genotypes. Then
This identity separates two contributors to a batch mean: the mean of called, centered genotypes and the fill value for missing entries. If the called-genotype mean is common across batches, write it as \(m_j\). For batches A and B using the same fill value \(c_j\), the difference in expected processed means is then
Thus a call-rate difference produces a mean difference when the observed centered genotypes have a nonzero mean and the fill value differs from that mean. If genotypes are centered so that \(m_j=0\), and missing entries are filled with \(c_j=0\), this particular expected batch shift is zero. That does not guarantee that missingness is harmless: it can still alter variances, covariances, and which individuals are compared using which variants.
How a Batch Shift Enters PCA
PCA is based on variation in the processed matrix across individuals and variants. Suppose a variant has two equally sized batches, with expected processed values \(a\) and \(b\). After centering across all individuals, the batch means are \((a-b)/2\) and \(-(a-b)/2\). The between-batch contribution to the variance at that variant is therefore
Across many variants, if the same batch division creates shifts in a coordinated direction, those shifts can contribute a strong sample-level contrast. PCA may then assign a leading component to that contrast. The effect is especially plausible when batch or processing date aligns with ancestry groups, recruitment sites, or another sample grouping. A PC can reflect a mixture of genetic and technical variation rather than one source alone.
This argument concerns the expected mean pattern. Missingness can also change the variance of processed values. For example, if centered observed genotypes have mean zero, missing entries are filled with zero, and the probability of observing a genotype is \(r\), then under independent missingness,
The variance changes with call rate even though the expected mean does not. Such heteroskedasticity does not automatically create a PC separating batches: under a simple null with independent, mean-zero genotypes and independent missingness, it need not create a systematic between-batch mean contrast. But real genotype covariance, correlated missingness across variants, finite-sample variation, and sample-specific patterns can make the resulting PCA sensitive to missingness. It is important to distinguish a possible mechanism from a guarantee that every missingness difference produces false structure.
Worked Example: Quantifying a Batch Mean Shift
Consider a hypothetical variant with a common mean centered dosage \(m_j=0.40\) among called genotypes in two equally sized batches. Genotypes are filled with \(c_j=0\). Batch A has call rate \(r_{Aj}=0.90\), and batch B has call rate \(r_{Bj}=0.60\). Assume the called-genotype mean is indeed the same in both batches, and the fill value is identical. The expected processed means are
and
The expected difference is \(0.36-0.24=0.12\) dosage units, also given directly by \((0.90-0.60)(0.40-0)=0.12\). With equal batch sizes, the centered batch means are \(+0.06\) and \(-0.06\). Their between-batch variance contribution is \(0.5(0.06^2)+0.5((-0.06)^2)=0.0036\) dosage-squared units. A single variant may contribute little, but a coordinated shift across many variants can become visible in PCA.
The assumptions matter. The calculation would not apply as written if the called-genotype means differed by batch, or if the two batches used different fill values. Those cases can be evaluated using the more general expectation \(r_{bj}m_{bj}+(1-r_{bj})c_j\) for each batch separately.
Worked Diagnostic: Do PC Scores Track Call Rate?
A practical first check is to compare each sample's call rate with its PC scores. This is a diagnostic association, not a causal test: ancestry and batch can be correlated, and both can be associated with call rate. The calculation below shows how such a check can be made and interpreted.
Worked Example: Correlating Call Rate With a PC
Suppose six hypothetical samples have call rates \(0.60, 0.65, 0.70, 0.90, 0.95,\) and \(1.00\). Their mean call rate is \(0.80\), so the centered call rates are \(-0.20,-0.15,-0.10,0.10,0.15,\) and \(0.20\). Suppose a centered PC score vector, in the same sample order, is approximately \(-1.038,-0.931,-1.028,1.028,0.931,\) and \(1.038\). Using the unrounded values that generated these displayed scores, the relevant sums are
where \(r_i\) is sample \(i\)'s call rate, \(\bar r=0.80\), and \(t_i\) is its PC score. The sample correlation is
The strong positive correlation means that, in this hypothetical set, samples with higher call rates tend to have higher scores on this PC. It does not establish that missingness caused the component: the association could arise because call rate and genuine ancestry are aligned, or because a third factor affects both. It does warrant checking batch labels, variant-level missingness, and whether the pattern persists under defensible changes to the data processing.
Checks That Help Identify the Failure Mode
A useful diagnostic set examines missingness at both the sample and variant levels, then asks whether the structure is robust to choices that affect missing entries. These checks should be interpreted together, because no single correlation or plot distinguishes technical structure from ancestry in every setting.
Calculate call rates by sample and missingness rates by variant. Look for unusually low sample call rates, variants with high missingness, and concentration of either pattern in a batch, plate, site, or processing period.
Plot or quantify associations between sample call rates and PC scores. Examine whether clusters or gradients in PC space coincide with missingness summaries or technical labels. An association is a warning signal, not a stand-alone verdict.
Check whether variants with strong loadings on a suspicious PC have distinctive missingness rates or batch-specific call rates. A component supported disproportionately by technically problematic variants deserves scrutiny.
Compare the broad PC pattern after reasonable sample or variant quality filters and under a justified alternative handling of missing entries. Large changes in sample positions or leading directions indicate sensitivity; they do not alone identify which version is correct.
When checking sensitivity, keep the sample set, variant set, and scaling as comparable as possible. Otherwise a changed PC may reflect a changed analysis target rather than a missingness artifact. Also avoid interpreting a component solely by its sign: multiplying all scores and loadings by \(-1\) represents the same PC. The relevant evidence is the pattern of associations and its stability, not whether scores are positive or negative.
Worked Example: Interpreting a Sensitivity Check
Suppose a hypothetical analysis finds that PC1 is strongly associated with sample call rate and with sequencing batch. After removing a small set of variants flagged for unusually high missingness, the correlation between PC1 scores from the original and filtered analyses is \(0.30\), after matching the direction of the component. Assume the samples and scaling are otherwise comparable. A correlation of \(0.30\) indicates that the sample ordering along this component has changed substantially; it is not evidence that the filtered PC is automatically correct. The analyst should check whether the removed variants were genuinely unreliable, whether filtering changed the genetic variation represented, and whether batch still tracks the remaining PCs.
If, instead, PC scores remain highly aligned after the filtering and the call-rate association weakens, that supports the view that the original component was sensitive to those variants. It still does not prove that all remaining structure is ancestry. The appropriate conclusion is limited to the diagnostic evidence: the observed component was, or was not, robust to that particular missingness-related change.
Alternatives, Limitations, and Common Pitfalls
There is no universal missing-data treatment that removes every risk. Filtering samples with very low call rates or variants with high missingness can reduce technical problems, but aggressive filtering can discard useful data or change which genetic variation is represented. Mean imputation can keep a rectangular matrix for PCA, but its behavior depends on the centering and fill values, as the batch-shift derivation shows. Pairwise-complete calculations use different sets of variants for different sample pairs; the resulting similarities can be difficult to compare and need not form a coherent positive semidefinite matrix. More elaborate imputation can retain information, but it can also preserve or introduce technical patterns if the missingness mechanism or imputation model is misspecified.
A missingness indicator, \(R_{ij}\), can be analyzed to reveal technical patterns directly, but it is not a substitute for genotype PCA and may reflect both technical and biological factors. Likewise, simply adding a missingness-associated PC as a covariate is not a guaranteed correction: the component may combine technical and ancestry variation, and removing it can remove genuine genetic signal. Diagnostic checks identify concerns; handling the concern requires an explicit analysis decision and awareness of what that decision changes.
- Pitfall: “Mean imputation makes missingness irrelevant.” Mean imputation replaces values, but if call rates differ and the relevant centered mean is not the fill value, processed means can differ by batch. Even when expected means match, variances and pairwise patterns can still be affected.
- Pitfall: “A PC correlated with call rate must be technical.” Correlation is not proof of cause. Ancestry, recruitment, and technical batch may be associated with one another; inspect several diagnostics and the study context.
- Pitfall: “Zero mean shift means no missingness problem.” Zero expected shift only rules out the particular mean-shift mechanism under its assumptions. It does not establish that covariance, variance, or finite-sample effects are negligible.
- Pitfall: “Filtering until the PC disappears proves the problem is fixed.” The component may disappear because filtering removed technical noise, genuine variation, or both. Check the reason for filtering and the stability of the remaining structure.
Check Your Understanding
Try each question, then compare your reasoning with the solutions.
- For one variant, \(m_j=0.30\), \(c_j=0\), and the call rates are \(0.80\) and \(0.50\) in two batches. What is the expected processed mean difference?
- If \(m_j=c_j=0\) in both batches, does a difference in call rate necessarily create a batch mean shift? What does this conclusion not rule out?
- Why is a strong association between a PC and sample call rate a diagnostic warning rather than proof that the PC is technical?
- Two PC analyses use different variant filters and produce substantially different sample scores. Give one interpretation and one reason not to conclude immediately that the second analysis is superior.
Solutions. (1) The difference is \((0.80-0.50)(0.30-0)=0.09\) dosage units. (2) No: the expected difference from this mechanism is \((r_{Aj}-r_{Bj})(0-0)=0\). This does not rule out changes in variance or covariance, correlated missingness, or finite-sample effects. (3) Call rate may be associated with batch, ancestry, or both; correlation alone cannot identify the cause. Batch labels, variant missingness, and sensitivity checks provide additional evidence. (4) The large change indicates sensitivity to the variant set. The second analysis is not automatically superior because filtering could also remove genuine genetic variation or otherwise change the analysis target.