[Stata] Firth Logistic Regression for Rare Outcomes and Separation: firthlogit
Logistic regression runs into trouble when an outcome is rare or a subgroup is small. If every person in some category of a predictor has the same outcome, the predictor separates the outcome perfectly, and the maximum likelihood estimate for that category does not exist; its coefficient heads toward infinity. Stata’s logit handles this by dropping the category and every observation in it, which quietly changes the sample and the comparison. Even without complete separation, maximum likelihood estimates from small samples with few events are biased away from zero.
Firth’s penalized likelihood (Firth, 1993) adds a small penalty to the likelihood that removes most of this small-sample bias and always yields finite estimates, and Heinze and Schemper (2002) showed that it solves the separation problem in logistic regression. The Stata command firthlogit by Joseph Coveney fits these models.
This comes up often in research on smaller language and immigrant groups, where a subgroup analysis can leave a few hundred respondents and a handful of events.
Installation
ssc install firthlogitThe current version is 1.2 (April 2021), and it runs on Stata 13.1 or newer.
Data
The examples use the 2024 California Health Interview Survey (CHIS) adult public use file, which has 24,810 adult respondents. The CHIS public use files are free and available through a brief application on the UCLA Center for Health Policy Research site. Every post in this series builds its variables with the same short do-file, 00_chis_prep.do:
use "chis_adult_2024.dta", clear
* outcomes
gen k6 = distress if distress >= 0 // K6 score, past month (0-24)
gen nousual = usual == 2 // no usual source of care
gen delay = ah22 == 1 // delayed or did not get medical care, past 12 months
gen erv = er == 1 // any ER visit, past 12 months
gen uninsured = ins == 2
gen needhelp = af81 == 1 if af81 > 0 // needed help for emotional/mental or substance problem
gen fairpoor = ab1 >= 4 // fair or poor self-rated health
gen spd = dstrs12 == 1 // likely serious psychological distress, past year
* groups
gen lep = spk_eng == 3 // speaks English not well or not at all
gen foreign = ah33new == 2
gen female = srsex == 2
* covariates
gen age = srage_p1
recode ahedc_p1 (1 2 = 1 "Less than high school") (3 = 2 "High school") ///
(4 5 6 = 3 "Some college") (7 8 9 = 4 "Bachelor's or more"), gen(educ)
gen pov = povll
label define povllf 4 "300%+ FPL", modify
label values pov povllf
recode ombsrr_p1 (1 = 1 "Latino") (2 = 2 "White") (3 = 3 "Black") (5 = 4 "Asian") ///
(4 6 = 5 "Other or multiracial"), gen(race)
* CHIS design: final weight plus 80 jackknife replicate weights
svyset [pw=rakedw0], jkrweight(rakedw1-rakedw80, multiplier(1)) vce(jackknife) mseLimited English proficiency (lep) means the respondent reported speaking English not well or not at all; 1,520 respondents (6.1 percent) meet that definition. Commands that accept the svy prefix use the jackknife replicate weights; commands that do not accept it use the final weight rakedw0 as a probability weight, which gives the right point estimates and robust standard errors that ignore the replicate design.
A small subgroup with a rare outcome
The example uses adults who speak only Korean at home and models being uninsured.
keep if lnghm_p1 == 5
count
tab uninsured
tab uninsured pov 191
Uninsured | Freq. Percent Cum.
------------+-----------------------------------
No | 180 94.24 94.24
Yes | 11 5.76 100.00
------------+-----------------------------------
Total | 191 100.00
| Family income (% FPL)
Uninsured | 0-99% FPL 100-199% 200-299% 300%+ FPL | Total
-----------+--------------------------------------------+----------
No | 22 25 19 114 | 180
Yes | 0 4 2 5 | 11
-----------+--------------------------------------------+----------
Total | 22 29 21 119 | 191There are 191 respondents and 11 uninsured. None of the 22 respondents below 100 percent of the poverty line is uninsured, so that income category perfectly predicts the outcome.
Standard logit
logit uninsured age female i.educ ib4.pov, ornote: 1.pov != 0 predicts failure perfectly;
1.pov omitted and 22 obs not used.
Logistic regression Number of obs = 169
-------------------------------------------------------------------------------------
uninsured | Odds ratio Std. err. z P>|z| [95% conf. interval]
--------------------+----------------------------------------------------------------
age | .9480172 .0239344 -2.11 0.034 .9022483 .9961078
female | 2.589541 1.754344 1.40 0.160 .6863676 9.76987
|
educ |
High school | .834411 1.40074 -0.11 0.914 .0310781 22.40301
Some college | .7891978 1.338941 -0.14 0.889 .0283824 21.94431
Bachelor's or more | .9151599 1.257381 -0.06 0.949 .0619429 13.5208
|
pov |
0-99% FPL | 1 (empty)
100-199% FPL | 4.742422 3.765863 1.96 0.050 1.00018 22.48653
200-299% FPL | 3.119656 2.87634 1.23 0.217 .5120285 19.00725
|
_cons | .4136718 .6572788 -0.56 0.579 .0183735 9.313673
-------------------------------------------------------------------------------------logit drops the lowest income category and the 22 people in it, so the model runs on 169 of the 191 respondents. The note at the top is the only sign of this, and it is easy to miss in a long log or a table built with esttab.
Firth logit
firthlogit does not take factor-variable notation, so the categorical predictors go in as indicator variables. The reference categories below match the logit model: less than high school and 300 percent of the poverty line or higher.
gen ed_hs = educ == 2
gen ed_sc = educ == 3
gen ed_ba = educ == 4
gen pov_lt100 = pov == 1
gen pov_100199 = pov == 2
gen pov_200299 = pov == 3
firthlogit uninsured age female ed_hs ed_sc ed_ba pov_lt100 pov_100199 pov_200299
firthlogit, or Number of obs = 191
Wald chi2(8) = 11.04
Penalized log likelihood = -30.212476 Prob > chi2 = 0.1996
------------------------------------------------------------------------------
uninsured | Odds ratio Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
age | .9560259 .0215278 -2.00 0.046 .9147497 .9991647
female | 2.344795 1.45867 1.37 0.171 .6927592 7.936472
ed_hs | .808314 1.131371 -0.15 0.879 .0520213 12.5597
ed_sc | .7150446 1.006092 -0.24 0.812 .045359 11.27204
ed_ba | .6434588 .7592273 -0.37 0.709 .0637052 6.499306
pov_lt100 | .6367047 .9614604 -0.30 0.765 .0330046 12.2829
pov_100199 | 4.228696 3.105513 1.96 0.050 1.002518 17.83695
pov_200299 | 3.166552 2.627904 1.39 0.165 .6225567 16.10624
_cons | .506679 .7081795 -0.49 0.627 .0327354 7.842375
------------------------------------------------------------------------------All 191 respondents remain in the model, and the separated category gets a finite odds ratio (0.64) with a wide interval that says, accurately, that 22 people with no events tell us little. The other estimates move toward 1: the odds ratio for age goes from 0.948 to 0.956 and for women from 2.59 to 2.34, which is the bias reduction at work.
Notes
Report the method. “Firth’s penalized likelihood logistic regression (firthlogit in Stata)” in the methods section, with the reason: small sample, rare outcome, or separation.
Confidence intervals. The default intervals are Wald intervals. Heinze and Schemper recommend profile penalized likelihood intervals, which are asymmetric and more accurate with few events; firthlogit‘s help file shows how to obtain likelihood ratio tests with constraint().
Weights and survey designs. firthlogit accepts frequency weights only, so the model above is unweighted, and it does not take the svy prefix. For a weighted or design-based analysis of a rare outcome, the options are to report the unweighted Firth model alongside a weighted model of the larger sample, or to use exact logistic regression (exlogistic) for very small samples.
Separation is information. A category with no events is a finding worth reporting in a table of counts. Firth’s method makes the model estimable; it does not create information that the data do not have.
References
Firth, D. (1993). Bias reduction of maximum likelihood estimates. Biometrika, 80(1), 27–38. https://doi.org/10.1093/biomet/80.1.27
Heinze, G., & Schemper, M. (2002). A solution to the problem of separation in logistic regression. Statistics in Medicine, 21(16), 2409–2419. https://doi.org/10.1002/sim.1047
