Default Step 2 Settings (is_noadjCov=TRUE) produce incorrect p-values for variants with AF > 0.5

[dannydelr0550]

Hello SAIGE team,

This is closely related to #197 and shares the same origin; thanks to Fred for bringing this up and helping me fix a problem in my GWAS. The consequences are pretty important -- if you select parameters that route you through the scoreTestFast_noadjCov() branch, this will lead to inaccurate effect size and p-value estimates for variants with AF > 0.05. A notable application I have found is that it can deflate the overall GC lambda of your GWAS to far below one, while the GC lambda for variants with very high MAF (> 0.4) is unaffected.

In short, there appears to be a bug in the branch of code that does not residualize the genotype against the covariates before calculating the score statistic. This is triggered by is_noadjCov=TRUE, which was implemented in July 2025 and is the default for SAIGE v1.5.0.

The bug is triggered by either of the following scenarios, assuming is_noadjCov is left as true:

• You use a full GRM
• You use a sparse GRM with fast_test = TRUE

Currently, imputeGenoAndFlip() flips alternate alleles with AF > 0.5 so that counts are given in terms of the minor allele. It flips the allele frequency along with it, giving the correct allele frequency for the minor allele:

if(t_altFreq > 0.5){
    flip = true;
    t_GVec = 2 - t_GVec;
    t_altFreq = 1 - t_altFreq; 
  }

but then at the end of imputeGenoAndFlip(), the allele frequency is flipped back:

if(flip){
	t_altFreq = 1 - t_altFreq;
        t_altCount = 2*dosagesSize - t_altCount;	
	//t_altCount = 2 * t_altFreq * dosagesSize;
     }	     

The allele frequency now corresponds to the original allele, but the genotype is flipped.

Later, the working score variance is calculated with this estimate of the allele frequency:

 double var2_a = dot(m_mu21,pow(g1,2));
    double var2_b = dot(m_mu21, 2*2*t_altFreq*g1);
    double var2_c = arma::accu(m_mu2)*pow(2*t_altFreq, 2);
    var2 = var2_a - var2_b + var2_c;

Which is just the following, expanded out:

$$ Var(S) = \sum_{i=1}^{m} \hat{\mu_i} ( 1 - \hat{\mu_i})(G_i - 2AF)^2 $$

If AF corresponds to the original allele and G_i is a count of the flipped allele, G_i will often be way off its estimated mean of 2*AF. As AF gets larger (e.g. 0.60 --> 0.80), this discrepancy will widen.

The simple fix is --is_noadjCov=FALSE. In my GWAS, this flag moved the GC lambda from 0.82 to 0.98, which was closely aligned with a simple logistic regression model I fit. This seems to mostly affect variants without large effect sizes to begin with (not sure why), so this deflation in the lower half of the p-value distribution is the most prominent result. Fred (issue 197) showed that on a variant with AF = 0.9515:
is_fastTest=TRUE: var=717.7, BETA=-0.018, p=0.64
is_fastTest=FALSE: var=20.1, BETA=-0.626, p=0.005

As a >30x increase in variance. In his example, he shows this with is_fastTest=TRUE (the Sparse GRM version of the path), but same idea.

Here is a reproducible example on simulated data, using SAIGE directly. You can first run simulate.R to make a plink binary fileset for a null GWAS for steps 1 and 2, run saige through run_saige.sh, and examine the results using compare.R . When 20% of the variants have AF > 0.5, the GC lambda goes from 0.77 ( --is_noadjCov=TRUE) to 0.99 (--is_noadjCov=FALSE).

saige_AF_05_sim.zip

I could be off base here, but I am pretty sure that this is accurate, especially given that it is a replication of a previous issue identified independently. Could someone please take a look? Thank you.

Best,
Danny