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(disease_state ~ target_genotypes, data=my_data, family = 'binomial') summary(regular_logistic)
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")
**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)
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?