Support for Firth Logistic Regression (for Rare Events)

294 views
Skip to first unread message

Parin Sripakdeevong

unread,
Mar 16, 2017, 8:52:38 PM3/16/17
to pystatsmodels
Hi All,

I am a research scientist at Dana-Farber Cancer Institute.

I have been using the pystatsmodels to run logistic regression models (specifically the functions sm.Logit() and sm.GLM(,,families.Binomial())

I am working with data where the # of positive outcomes in certain explanatory group has very low prevalence (e.g. ~10 positives, ~10,000 negatives).

In this situation, the maximum likelihood estimation in convention logistic model is known to suffer from small-sample bias. For discussion, see http://statisticalhorizons.com/logistic-regression-for-rare-events.

I searched online and it seem that the solution is to use 'Firth Logistic Regression'. This model appear to be supported in both SAS and R. For example, see https://www.r-bloggers.com/example-8-15-firth-logistic-regression/.

However, it appear that Firth Logistic Regression is currently not support by pystatsmodel. Is this correct?

If so, can anyone give me any tips/instructions on how I could extend/modify the existing functionality of pystatsmodel to support Firth Logistic Regression?

Thank you!

Sincerely,
Parin S.


josef...@gmail.com

unread,
Mar 16, 2017, 9:35:02 PM3/16/17
to pystatsmodels
Hi,

sorry for the delay in posting, I got the notification about pending
messages just a short time ago.

On Wed, Mar 15, 2017 at 12:14 PM, Parin Sripakdeevong
<pa...@jimmy.harvard.edu> wrote:
> Hi All,
>
> I am a research scientist at Dana-Farber Cancer Institute.
>
> I have been using the pystatsmodels to run logistic regression models
> (specifically the functions sm.Logit() and sm.GLM(,,families.Binomial())
>
> I am working with data where the # of positive outcomes in certain
> explanatory group has very low prevalence (e.g. ~10 positives, ~10,000
> negatives).
>
> In this situation, the maximum likelihood estimation in convention logistic
> model is known to suffer from small-sample bias. For discussion, see
> http://statisticalhorizons.com/logistic-regression-for-rare-events.
>
> I searched online and it seem that the solution is to use 'Firth Logistic
> Regression'. This model appear to be supported in both SAS and R. For
> example, see
> https://www.r-bloggers.com/example-8-15-firth-logistic-regression/.
>
> However, it appear that Firth Logistic Regression is currently not support
> by pystatsmodel. Is this correct?

correct, it's not supported
we have some related open issues
https://github.com/statsmodels/statsmodels/issues/2293
https://github.com/statsmodels/statsmodels/issues/2282


>
> If so, can anyone give me any tips/instructions on how I could extend/modify
> the existing functionality of pystatsmodel to support Firth Logistic
> Regression?

I only know about Firth for handling perfect separation cases, and
didn't know about the bias correction part. I also never looked at
what kind of penalization Firth Logit actually uses.
We have now elastic net L1/L2 penalization for GLM and more structured
L2 penalization waiting in a PR. However, based on a quick look, Firth
seems to use the log determinant of the information matrix (jeffreys'
prior) for the penalization term which is not covered by any of the
existing or work in progress penalization terms.

Based on my brief skimming, there might be a way to implement this by
data augmentation, otherwise we would need to use either a special
case IRLS or a generic penalization method.


So for now I don't know what the best way of implementing this is.

Josef

Parin Sripakdeevong

unread,
Mar 17, 2017, 1:32:24 PM3/17/17
to pystatsmodels
Hi Josef,

Thank you for looking into this issue!

Sincerely,
Parin S.

josef...@gmail.com

unread,
Mar 17, 2017, 2:38:43 PM3/17/17
to pystatsmodels
According to Firth 1993, it looks like it would be enough to add h_i /
2 to the count of successes and failures, where h_i is the diagonal of
the (symmetric) hat matrix, and then use standard MLE. This would work
for GLM Binomial but not for Logit which doesn't take counts.

I won't have time to compare with other packages to check whether the
standard errors and other results are still correct in that case.

All other methods to estimate a penalized log-likelihood will be more
work, at least for results beyond params.

Josef




On Fri, Mar 17, 2017 at 1:32 PM, Parin Sripakdeevong

josef...@gmail.com

unread,
Mar 18, 2017, 3:12:38 PM3/18/17
to pystatsmodels
Here's a just published article that compares different penalties for
effect/parameter estimation and prediction

Puhr, Rainer, Georg Heinze, Mariana Nold, Lara Lusa, and Angelika
Geroldinger. 2017. “Firth’s Logistic Regression with Rare Events:
Accurate Effect Estimates and Predictions?” Statistics in Medicine,
January, n/a-n/a. doi:10.1002/sim.7273.

I opened an issue to track this
https://github.com/statsmodels/statsmodels/issues/3561

I haven't looked at inference yet, Firth and similar recommend profile
likelihood, no idea which methods would be best for small sample
hypothesis tests and confidence intervals, (analogous to hypothesis
tests and confidence intervals for proportions without regression.)

PRs, examples and test cases would be welcome.

Josef
(how many new topics can we study in half a year?)
Reply all
Reply to author
Forward
0 new messages