Skip to content

Use an unambiguous missing-value sentinel in marker regressions - #111

Open
thierrygosselin wants to merge 1 commit into
bcm-uga:masterfrom
thierrygosselin:fix-regression-missing-sentinel
Open

thierrygosselin wants to merge 1 commit into
bcm-uga:masterfrom
thierrygosselin:fix-regression-missing-sentinel

Conversation

@thierrygosselin

Copy link
Copy Markdown
Contributor

Summary

This PR prevents legitimate standardized genotype values equal to 3 from being mistaken for missing genotypes during marker-level regression.

Problem

multLinReg() currently constructs a scaled genotype accessor using 3 as its missing-value sentinel:

matAccScaled macc(obj, ind_col, af, ploidy, 3);

It then detects missing observations with:

not_missing[i] = (x[i] != 3);

However, x[i] is no longer a raw genotype. It is the centred and scaled genotype value:

z = (g - ploidy × p) / sqrt(ploidy × p × (1 - p))

A legitimate observed genotype can therefore equal 3.

For diploid data, an alternate homozygote at allele frequency p = 2/11 gives exactly:

(2 - 2 × 2/11) / sqrt(2 × 2/11 × 9/11) = 3

This frequency is realizable, for example, with four alternate alleles among 22 observed chromosomes.

The affected alternate homozygotes are consequently removed from the marker regression as though their genotypes were missing.

Demonstrated effect

Using 11 diploid individuals with two observed alternate homozygotes and nine reference homozygotes gives an allele frequency of 2/11.

In a controlled two-component example, the current sentinel collision changed the marker statistics from:

Using all observed genotypes:
-0.636317, -2.119022

to:

After legitimate values equal to 3 were treated as missing:
0.103821, 1.424499

Both component statistics changed substantially and reversed sign.

Changes

  • Represent missing scaled genotypes with R’s numerical missing value rather than the ordinary number 3.
  • Detect missing values using NumericVector::is_na().
  • Apply the correction to both matrix-backed and BED-backed genotype accessors.
  • Preserve the existing marker-regression calculation otherwise.
  • Add a regression test using the exactly realizable p = 2/11 collision.
  • Add a separate test confirming that genuinely missing genotypes remain excluded.

Scope

This PR only corrects the ambiguous sentinel.

It does not change the statistical treatment of missing genotypes or resolve the separate question of whether marker-specific regression should use the full cross-product matrix of the observed PC scores.

Validation

  • 42 package tests passed.
  • The package compiled successfully.
  • R CMD check --no-manual: 0 errors, 0 warnings and 0 notes when unavailable suggested packages were not forced.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant