Firth Logistic Regression for Sparse Association Data

Why Penalize Logistic Likelihood?

In genetic association analyses, a genotype group or a combination of genotype and covariates can contain very few cases. Ordinary logistic regression may then give highly biased estimates or no finite maximum-likelihood estimate at all. In particular, if a covariate pattern perfectly separates cases from controls, a coefficient can grow without bound while the likelihood approaches its supremum. Firth logistic regression modifies the likelihood to reduce first-order estimation bias and can provide finite estimates in many such sparse or separated settings.

The method does not create information that the data do not contain. Sparse estimates can remain uncertain, model assumptions still matter, and finite estimates are not by themselves evidence of a reliable association. The central ideas are the exact form of the penalty, its effect on the logistic score, and careful calibration of tests and intervals.

Model and Penalized Likelihood

Let \(Y_i\) be the binary case indicator for individual \(i\), and let \(x_i\) be the vector of model covariates, including an intercept and the genotype coding being studied. Under the logistic regression model, the observations are conditionally independent Bernoulli variables with \(\operatorname{logit}(p_i)=x_i^\mathsf{T}\beta\), where \(p_i=P(Y_i=1\mid x_i)\). The usual log likelihood is \(\ell(\beta)=\sum_i\{y_i\log(p_i)+(1-y_i)\log(1-p_i)\}\).

Write \(X\) for the design matrix, \(W\) for the diagonal matrix with entries \(p_i(1-p_i)\), and \(I(\beta)=X^\mathsf{T}WX\) for the model’s Fisher information. Assume the model is correctly specified for the target association analysis, observations are independent conditional on the included covariates, and the parameters being estimated are identifiable, so \(I(\beta)\) is nonsingular in the region of interest.

Definition: Firth logistic regression maximizes the bias-reduced (or penalized) log likelihood \(\ell^*(\beta)=\ell(\beta)+\tfrac12\log|I(\beta)|\). The added term is the log of the square root of the determinant of the Fisher information. Equivalently, it is the Jeffreys-information penalty for this model. It is a mathematical penalty, not a set of literal extra individuals or observations.

The penalty depends on the model’s information, so it must be evaluated for the model being fitted. It discourages estimates from moving into regions where the fitted probabilities make the information nearly singular. Unlike ordinary maximum likelihood, Firth estimation often remains finite under complete or quasi-complete separation, provided the design is identifiable and the data/model permit a finite bias-reduced solution. It cannot repair rank deficiency, coding errors, or a misspecified association model.

How the Adjustment Enters the Score

The ordinary logistic score is \(U(\beta)=X^\mathsf{T}(y-p)\). Differentiating the penalty gives a leverage adjustment to each observation’s residual. Define \(H=W^{1/2}X(X^\mathsf{T}WX)^{-1}X^\mathsf{T}W^{1/2}\), and let \(h_i\) be its \(i\)th diagonal element. Then the adjusted score is

U^*(\beta)=X^\mathsf{T}\{y-p+h\mathbin{\circ}(1/2-p)\},

where \(\mathbin{\circ}\) denotes elementwise multiplication. Thus, observation \(i\) contributes \(y_i-p_i+h_i(1/2-p_i)\) to the residual-like quantity. The adjustment is largest for observations with high leverage; it pulls fitted probabilities away from the boundaries 0 and 1. The estimate solves \(U^*(\beta)=0\), ordinarily by iterative numerical optimization.

Justification of the adjustment: Since \(\partial I/\partial\beta_j=X^\mathsf{T}\operatorname{diag}\{(1-2p_i)p_i(1-p_i)x_{ij}\}X\), the derivative of \(\tfrac12\log|I|\) with respect to \(\beta_j\) is one-half the trace of \(I^{-1}\partial I/\partial\beta_j\). The trace simplifies to \(\sum_i h_i(1/2-p_i)x_{ij}\). Adding this derivative to the ordinary score gives the displayed adjusted score.

This derivation also explains why Firth fitting is not equivalent to adding one-half to every genotype-by-case cell in general. That simple correction is available for a saturated two-group table, but with covariates the adjustment depends on the fitted information and the leverage of each observation.

Worked Example: A Genotype Group With No Cases

Suppose hypothetical data compare noncarriers with carriers in a binary case-control sample. There are 12 noncarriers: 0 cases and 12 controls. There are 12 carriers: 4 cases and 8 controls. The ordinary maximum-likelihood estimate of the carrier log odds ratio diverges because the noncarrier group has no cases.

Fit the saturated logistic model with an intercept and a carrier indicator. In this special two-group setting, maximizing the Firth penalized likelihood gives fitted group probabilities \((y_g+0.5)/(n_g+1)\). For noncarriers this is \(0.5/13=0.03846\); for carriers it is \(4.5/13=0.3462\). This formula is the exact Firth solution for this saturated table, not a general rule for covariate-adjusted models.

The estimated odds ratio for carriers versus noncarriers is therefore \((0.3462/0.6538)/(0.03846/0.9615)=13.24\). Equivalently, the corrected cross-product calculation is \((4.5\times12.5)/(8.5\times0.5)=56.25/4.25=13.24\). The corresponding log odds ratio is \(\log(13.24)=2.583\). The estimate is finite, but the small table still signals substantial uncertainty; finiteness alone does not establish a precise or robust effect.

The half-count expression is a way to state this particular maximizer. It does not mean that the analysis observed half a case or half a control in either group.

What Happens in a Sparse Covariate Pattern?

When covariates are included, Firth’s penalty applies to the whole design matrix. A rare genotype-by-covariate pattern can have high leverage, and its score adjustment may therefore matter even when the overall case-control sample is not tiny. The following intercept-only calculation makes the adjustment visible before considering a larger model.

Worked Example: A One-Case Sample and the Adjusted Score

Consider hypothetical data with 20 independent individuals and one case. Fit an intercept-only logistic model, so every individual has probability \(p\). At \(\beta=0\), \(p=0.5\), each leverage is \(h_i=1/20=0.05\), and the information is \(I=20(0.5)(0.5)=5\).

The ordinary score at zero is \(1-20(0.5)=-9\). Because \(p=0.5\), the adjustment \(h_i(1/2-p_i)\) is zero at this particular starting value, so the adjusted score is also \(-9\). A Fisher-scoring step is \(0+(-9)/5=-1.8\). Further iterations are needed; this first step is not the final estimate.

For the exact solution, the penalized log likelihood, up to an additive constant, is \(\ell^*(p)=\log(p)+19\log(1-p)+\tfrac12\log\{20p(1-p)\}\). Its derivative with respect to \(p\) is \(1.5/p-19.5/(1-p)\). Setting it to zero gives \(1.5(1-p)=19.5p\), hence \(\hat p=1.5/21=0.07143\). The fitted intercept is \(\log\{0.07143/(1-0.07143)\}=\log(1/13)=-2.565\), finite despite there being only one case.

This example isolates the finite-estimation mechanism; it is not a genotype-effect estimate. In an association model, the same penalty is computed from the full information matrix for the intercept, genotype, and included covariates.

Testing an Association Effect

A genotype effect can be tested with a penalized likelihood-ratio procedure, but the construction of the null fit is essential. Let the full model include the genotype parameter and any prespecified nuisance covariates. To test that the genotype parameter is zero, maximize the full model’s penalized likelihood without restriction, and compare it with the maximum under the constraint that the genotype parameter equals zero. In that constrained maximization, retain the penalty from the full model’s information matrix, evaluated at the constrained parameter values. This is a profile penalized likelihood comparison.

Under regularity conditions and an adequately informative sample, twice the difference between these unrestricted and constrained maximized values is commonly compared with a chi-square distribution whose degrees of freedom equal the number of genotype parameters constrained. In sparse data, the chi-square approximation can be inaccurate; an appropriate calibration or alternative justified test may be needed. One must not instead subtract separately maximized penalties from models of different dimensions: their information determinants scale differently with sample size, so that comparison can add an artificial dimension-dependent term.

Worked Example: Testing the Carrier Effect With the Full-Model Penalty

Use the hypothetical table from the first example: noncarriers have 0 cases and 12 controls; carriers have 4 cases and 8 controls. Use the intercept-plus-carrier full model. For group probabilities \(p_0,p_1\), the determinant of the full model’s information is \(12p_0(1-p_0)\times12p_1(1-p_1)\).

In the unrestricted fit, the probabilities are \(p_0=0.5/13=1/26\) and \(p_1=4.5/13=9/26\). The ordinary log likelihood, omitting binomial constants, is \(12\log(25/26)+4\log(9/26)+8\log(17/26)=-8.116\). The information determinant is \(\{12(25/676)\}\{12(153/676)\}=1.205\), so the penalty is \(\tfrac12\log(1.205)=0.0934\). Thus the maximized full penalized log likelihood is \(-8.023\).

Under the null constraint, the carrier coefficient is zero, so both groups have a common probability \(p\). Crucially, retain the two-parameter full-model information penalty: it is \(\tfrac12\log\{144[p(1-p)]^2\}=\log\{12p(1-p)\}\). The constrained penalized log likelihood is therefore \(4\log p+20\log(1-p)+\log\{12p(1-p)\}\). Its maximum occurs at \(p=5/26=0.1923\), as follows by differentiating: \(5/p-21/(1-p)=0\). Substitution gives a constrained maximum of \(-10.244\).

The profile penalized likelihood-ratio statistic is \(2\{-8.023-(-10.244)\}=4.441\). If a one-degree-of-freedom chi-square reference is judged appropriate, this gives an approximate two-sided \(p\)-value of 0.035. The sample is sparse, however, so the numerical chi-square calibration should not be treated as guaranteed. This computation illustrates the correctly constructed statistic; it does not establish that its asymptotic reference is adequate for every sparse study.

Intervals, Interpretation, and Limits

Firth estimates can be reported as log odds ratios or exponentiated odds ratios, with the genotype coding and covariate adjustment stated. An odds ratio compares the odds of the case outcome across the specified genotype contrast, conditional on modeled covariates; it is not a risk ratio and does not imply a causal effect. An interval from a profile penalized likelihood is obtained by examining parameter values whose constrained penalized likelihood remains within a chosen likelihood-ratio cutoff of its maximum. Such intervals can be more useful than symmetric Wald intervals when the likelihood is asymmetric, but they still require a defensible calibration.

Firth’s bias reduction is a first-order property under the model and regularity conditions; it is not a guarantee of negligible bias in every tiny sample. The fitted effect can remain sensitive to a few observations, omitted covariates, phenotype or genotype errors, and the chosen genetic model. For genome-wide testing, the penalty does not itself address multiple testing or ensure correct null calibration. Those questions require their own analysis.

Specify the model.
State the binary outcome, genotype coding, covariates, and sampling unit. Check that the design is identifiable.

Fit with the information penalty.
Maximize \(\ell(\beta)+\tfrac12\log|I(\beta)|\), using the full design matrix for adjusted analyses.

Choose inference deliberately.
For a penalized likelihood-ratio test, use the full-model penalty in both unrestricted and null-constrained fits, and assess whether the reference calibration is adequate.

Common Mistakes

Calling the penalty added data. The saturated two-group table has a half-count representation, but the general method changes the likelihood mathematically; it does not add observed cases or controls.
Assuming every finite result is trustworthy. Finiteness addresses a boundary problem, not uncertainty, misspecification, or influential observations.
Comparing separately penalized models of different dimensions. The determinants of their information matrices scale differently. For the profile penalized likelihood-ratio construction, constrain the genotype term while retaining the full model’s penalty.
Reporting a sparse-data chi-square test as exact. The reference is asymptotic. Sparse samples may need an alternative or additional calibration justified for the setting.

Key takeaway: Firth logistic regression adds one-half the log determinant of the Fisher information to the logistic log likelihood. Its leverage-adjusted score can yield finite, bias-reduced estimates in sparse settings. Valid testing requires a carefully defined null fit and a reference distribution appropriate to the data.

Check Your Understanding

Try these questions without looking back.

What is the Jeffreys-information penalty added to the logistic log likelihood, and what does \(I(\beta)\) represent?
Why is a half-count correction exact for the saturated two-group table but not a general recipe for covariate-adjusted Firth regression?
In a profile penalized likelihood-ratio test of a genotype coefficient, which penalty must be used in the constrained fit?
What does a finite Firth estimate address, and what important sources of uncertainty or bias does it not automatically remove?
Why should a chi-square reference for a sparse-data penalized likelihood-ratio statistic be treated with care?