Logistic Regression on very rare events - R

Question

I am trying to run a logistic regression analysis in R using the speedglm package. CNVs are a type of genetic variant. We measure whether CNVs occur in cases or controls and whether genes in a pathway are "hit" by the CNV or not (Pathway.hit), and how many genes were hit by the CNV that were not in the pathway (Pathway.out). I run two models with and without the Pathway.hit covariate and compare to see if a pathway is preferentially hit by cases vs controls.

the models and comparison of said are as follows:

fit1 = speedglm(status~size+Pathway.out, data=cnvs, family=binomial('logit'))
fit2 = speedglm(status~size+Pathway.out+Pathway.hit, data=cnvs,family=binomial('logit'))
P.anova = 1-pchisq(abs(fit1$deviance - fit2$deviance), abs(fit1$df - fit2$df))

It seems to work okay for most data I throw at it, but in a few cases I get the error:

Error in solve.default(XTX, XTz, tol = tol.solve) : 
  system is computationally singular: reciprocal condition number = 1.87978e-16

After some googling around I think I found what's causing the problem:

by(cnvs$Pathway.hit, cnvs$status, summary)
cnvs$status: 1 (controls)
    0     1 
13333     0 
------------------------------------ 
cnvs$status: 2 (cases)
    0     1 
10258     2 

So here there no observations in controls and only 2 in cases.

If I use with normal glm method however, then it does not throw an error (but that of course doesn't necessarily mean the results will be meaningful). The reason I am using the speedglm package is that I have approximately 16,000 of these analyses to run, and using the base glm function for all 16,000 takes about 20 hours, where as I think speedglm can reduce it down to 8 or so.

So my question is, should I ignore those analyses which throw an error and list the results as NA as there were too few observations, or when speed glm fails should I retry with normal glm? In the above example there are 2 observations of the covariate in cases, but 0 in controls. Might this not be interesting? Would the analysis also fail if there were 0 in controls and 20 in cases - that would certainly be interesting would it not?

Answer 1

The problem is that you only have two instances of Pathway.hit and they are all cases. If you fit the univariate model, the odds ratio would be infinite. If you inspect the covariates for the fit2 model, you will probably see that they achieve unbounded values, with log odds ratios taking values of 8 or more. This is a typical issue of separation, as @Scortchi alluded to.

This is also an issue of small cell counts. On this topic there is quite a literature. One simple and easy to implement method gives you a biased estimator that usually converges. A preprocessing step is converting your data to a tabular format.

Inference on GLMs fit to tabular data should not be done with the deviance because the degrees of freedom are calculated based on the format of the data. You can input these dfs by hand, or you can perform the likelihood ratio test instead using the base function logLik. e.g. 2*pchisq(logLik(fit1)-logLik(fit2), lower.tail=F, df=1). Or just use the lrtest from the package lmtest for the nested models.

Obtain tabular data using the following:

cnvs <- as.data.frame(table(cnvs))

Or some equivalent argument from data.table.

And fit the same model fit2 <- glm(status~size+Pathway.out+Pathway.hit, data=cnvs,family=binomial('logit'),weights=Freq)

to obtain equivalent inference to the parent model. Then apply the small cell count correct by just adding 1 to the Freq variable.

fit2.scc <- glm(status~size+Pathway.out+Pathway.hit, data=cnvs,family=binomial('logit'),weights=Freq+1)

For some intuition, the former OR of 2/0 = infinity, would in this case be 3/1 = 3 which again is biased, but a reasonable measure of association, and 3 is a big odds ratio.

Inference and estimates from this model will be far more well behaved.

Now, I don't know who wrote speedglm or why. GLMs usually converge pretty fast. 24,000 cases is not a lot of data. In some rare cases I need to speed the glm up, like fitting a jackknife variance using the DF-betas with 24,000 glms of size 23,999, so I call the glm.fit workhorse. I don't know if GLMs can get any faster than that without implementing some parallel computing, and even then I think methods are lacking since it is an iterative subroutine. glm.fit almost directly calls the Fisher Scoring algorithm from LAPACK.