0
$\begingroup$

I am performing lots of association tests between genotypes and binary diseases. The genotypes can be very rare, and my tests often have extreme case control imbalance (e.g. 500 controls for every case among 50K samples). When running logistic regression using glm() in R, my test often produce the error: glm.fit: fitted probabilities numerically 0 or 1 occurred.

Here is some dummy data that replicates the issue. There are more predictors than a single "target genotype" in my real data (multiple genotypes at a multi-allelic site), but I leave them out here for simplicity:

set.seed(001)
reference_genotypes_control <- sample(c(rep(0,5),rep(1,495), rep(2,39500)))
target_genotypes_controls <- sample(c(rep(0,39750),rep(1,249), rep(2,1)))
disease_state_controls <- rep(0,40000)

controls <- data.frame(reference_genotypes=reference_genotypes_control,
                       target_genotypes=target_genotypes_controls,
                       disease_state=disease_state_controls)

reference_genotypes_cases <- sample(c(rep(1,5),rep(2,495)))
target_genotypes_cases <- rep(0,500)
disease_state_cases <- rep(1,500)

cases <- data.frame(reference_genotypes=reference_genotypes_cases,
                       target_genotypes=target_genotypes_cases,
                       disease_state=disease_state_cases)

my_data <- rbind(controls,cases)

Genotypes can be 0,1, or 2. I am interested in the impact of the target genotype on the binary disease state. I leave reference genotype out of the model specification as reference.

Regular old logistic regression gives the glm.fit error:

glm(disease_state ~ target_genotypes, data=my_data, family = 'binomial')
summary(regular_logistic)

which gives:

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.1581  -0.1581  -0.1581  -0.1581   2.9625  

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)        -4.376      0.045 -97.235   <2e-16 ***
target_genotypes  -13.180    249.390  -0.053    0.958    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 5388.3  on 40499  degrees of freedom
Residual deviance: 5382.0  on 40498  degrees of freedom
AIC: 5386

Number of Fisher Scoring iterations: 16

Am I having issues because the target genotype is always 0 in the cases and this is considered complete separation issue? My coefficients and standard error are very large and p value is close to 1 for target_genotypes, which I have read is a side effect of complete separation. I also checked complete separation with brglm2.

library("brglm2")
glm(disease_state ~ target_genotypes,data=my_data,family=binomial("logit"),method="detect_separation")

which gives:

**Separation: TRUE** 
Existence of maximum likelihood estimates
     (Intercept) target_genotypes 
               0             -Inf 
0: finite value, Inf: infinity, -Inf: -infinity

I read that logistf can work with data that produces the separation problem, so I tried to run:

logistf(disease_state ~ target_genotypes, data = my_data, firth=TRUE,pl=TRUE)

which gives:

Model fitted by Penalized ML
Coefficients:
                      coef   se(coef) lower 0.95  upper 0.95   Chisq          p
(Intercept)      -4.374772 0.04498013  -4.464213 -4.28785813     Inf 0.00000000
target_genotypes -1.836535 1.41644734  -6.670910  0.07811345 3.41453 0.06462481
                 method
(Intercept)           2
target_genotypes      2

Method: 1-Wald, 2-Profile penalized log-likelihood, 3-None

Likelihood ratio test=3.41453 on 1 df, p=0.06462481, n=40500
Wald test = 1.681115 on 1 df, p = 0.1947764>

Is the logistf approach I show here an appropriate way to deal with my apparent separation problem?

$\endgroup$
5
  • 1
    $\begingroup$ Since yu seem to be fitting predictors one-at-a-time and your interest seems to be limited to p-values (why?) you can ignore the need for penalization and just compute the score test or likelihood ratio $\chi^2$ for association and be done with it. If the predictor is binary the score test is the ordinary Pearson $\chi^2$ test. $\endgroup$ Nov 27 '20 at 17:02
  • $\begingroup$ There are more predictors in the model, I just put the one for simplicity to demonstrate the concept. I updated the post to indicate this. $\endgroup$
    – curious
    Nov 27 '20 at 17:05
  • $\begingroup$ Thanks. Some of that still applies. As long as you stay away from Wald tests (i.e., you don't use the standard error in any way) you will get valid p-values without penalization (and it is hard to get p-values with penalization; it's even hard to define them). $\endgroup$ Nov 27 '20 at 17:18
  • $\begingroup$ >it is hard to get p-values with penalization I am getting p values in my penalized logistf example above, are you just saying that they are not readily interpretable? Also lets say I wanted to compare two nested models and get valid p values even if fit1 had a complete separation problem, would this be the score-based approach? input <- mtcars[,c("am","cyl","hp","wt")] fit1 = glm(formula = am ~ cyl + hp + wt, data = input, family = binomial) fit2 = glm(formula = am ~ cyl + hp, data = input, family = binomial) anova(fit1, fit2, test="Rao") Thank you. $\endgroup$
    – curious
    Nov 29 '20 at 16:34
  • $\begingroup$ p-value are intended for situations where you do not impose external knowledge and don't necessarily work when the model is intentionally biased (here, towards the null). If you set your sights as low as wanting p-values, it's best to compute them without penalization. $\endgroup$ Nov 30 '20 at 12:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.