Caveats on Using Firth's Penalization in the Model-Based Regression Standardization for Rare Diseases.

Introduction

Model‐based regression standardization [1], also known as parametric g‐formula or g‐computation for fixed exposures [2], is a popular method for the estimation of marginal effect measures in both randomized and observational studies [3,4]. For example, logistic regression models for binary outcomes provide not only the covariate‐conditional odds ratios by exponentiated regression coefficients but also the marginal “counterfactual” risks by predicted outcomes under exposed and unexposed conditions, from which risk differences and ratios can be derived as marginal effect measures.

In observational studies of rare diseases, we may encounter a small number of events relative to the number of confounding variables which should be adjusted for in the models. In such cases, a problem known as (quasi‐)complete separation may occur, and the nonexistence of the maximum likelihood estimates of the logistic regression models can be problematic [5,6,7]. Even if statistical software or packages provide finite estimates in datasets with “separated” outcomes, the obtained values are contingent on the stopping criterion of the iterative computation in each program [8]. Furthermore, the estimated standard errors of the regression parameters increase to near infinity [7,9]. If such separation matters, then model‐based regression standardization would also provide invalid effect estimates.

A convenient solution to this convergence issue in generalized linear models, including logistic regression models, is to use Firth's penalized likelihood estimates [6,10]. This approach may yield finite estimates, even when the maximum likelihood estimates would become infinite owing to separation [6,9]. Originally designed to correct for first‐order bias in maximum likelihood estimates [10], Firth's penalization is applicable even in situations with sufficient event counts, regardless of the presence of separation [6,8,9,11,12,13]. Hence, it may be attractive to apply Firth's method to the estimation of regression parameters in model‐based standardization estimators in rare disease settings without practical concerns. However, the use of Firth's penalization in model‐based standardization has a problem similar to those previously reported by Puhr et al. [14] regarding biases in the estimation of predicted probabilities. Specifically, Firth's penalized likelihood estimates tend to pull the predicted probabilities toward 0.5, leading to a non‐negligible discrepancy between the predicted and observed event rates when events are extremely rare or common [14]. As model‐based standardization estimators are based on the predicted outcomes under different exposures, Firth's method may inflate rather than reduce bias in rare‐event situations, even if the regression coefficients in the model are accurately estimated. Our study examines the potential biases introduced by using Firth's penalized likelihood for estimating marginal effects, within the framework of a counterfactual model in causal inference, thereby underscoring the novelty of this research.

In this study, we illustrate the possible bias using Firth's method for model‐based regression standardization estimators, which has received little attention in the literature. In Section2, we introduce model‐based regression standardization estimators and Firth's method to estimate the underlying logistic regression models. Section3compares regression standardization estimators with or without following Firth's penalization, with illustration using stratified data from a study of surgical site infection (SSI) in orthopedic surgeries. This section also presents twoad hocmethods of Firth's penalized estimates that correct the predicted probabilities such that their average equals the observed event rate in each exposure group. These regression standardization methods are evaluated along with propensity score‐based methods using simulation studies in Section4in terms of bias and convergence rates. Finally, in Section5, the association of SSI (which is known as a rare disease) and smoking (a known risk factor for the infection) was analyzed in the database of orthopedic surgeries.

Regression Standardization Estimators Following Firth's Penalization

Notation and Assumptions

Let$Y_{i}$be the outcome (1: event occurred, 0: event did not occur),$Z_{i}$be the exposure (1: exposed, 0: unexposed), and$\mathbf{\mathit{L}}_{i}$be the set of confounding variables for patient$i$($i = 1 , \ldots , n$). Using a potential outcome$Y_{i} \left(z\right)$under exposure$Z_{i} = z$with thecausal consistencyassumption, that is,$Y_{i} = Y_{i} \left(z\right)$if$Z_{i} = z$[15], the marginal causal effect can be defined as the contrast of the potential outcome means

$$E \left[Y \left(1\right)\right] - E \left[Y \left(0\right)\right] .$$

Hereafter, we omit the subscript$i$if it is unnecessary. This marginal “mean” effect (1) is identifiable from the probability distribution of observed data ($Y_{i} , Z_{i} , \mathbf{\mathit{L}}_{i}$) if the following two assumptions are additionally met [15]: themean conditional exchangeability, i.e.,$E \left[\right. Y \left(z\right) \left|\right. Z = z , \mathbf{\mathit{L}} \left]\right. = E \left[\right. Y \left(z\right) \left|\right. \mathbf{\mathit{L}} \left]\right.$, and theconditional positivity, i.e.,$0 < P \left(\right. Z = z \left|\right. \mathbf{\mathit{L}} \left.\right)$, across possible support of$\mathbf{\mathit{L}}$and for$z = 0 , 1$. Under these assumptions, we obtain the following identification formula:

$$E \left[Y \left(1\right)\right] - E \left[Y \left(0\right)\right] = \int R \left(1 , \mathbf{\mathit{l}}\right) f_{\mathbf{\mathit{L}}} \left(\mathbf{\mathit{l}}\right) d \mathbf{\mathit{l}} - \int R \left(0 , \mathbf{\mathit{l}}\right) f_{\mathbf{\mathit{L}}} \left(\mathbf{\mathit{l}}\right) d \mathbf{\mathit{l}} ,$$

where$R \left(z , \mathbf{\mathit{l}}\right) = E \left[\right. Y \left|\right. Z = z , \mathbf{\mathit{L}} = \mathbf{\mathit{l}} \left]\right.$is regression of$Y$on$\left(Z , \mathbf{\mathit{L}}\right)$evaluated at$\left(z , \mathbf{\mathit{l}}\right)$, and$f_{\mathbf{\mathit{L}}} \left(\mathbf{\mathit{l}}\right)$is the marginal density of$\mathbf{\mathit{L}}$at$\mathbf{\mathit{l}}$.

Model‐Based Regression Standardization Estimators With and Without Penalized Likelihood

To estimate the identification formula (2) from the observed data, we may substitute$f_{\mathbf{\mathit{L}}} \left(\mathbf{\mathit{l}}\right)$with the empirical distribution (i.e., the probability distribution that places probability$1 / n$at each observed value$\mathbf{\mathit{L}}_{i}$) and$R \left(z , \mathbf{\mathit{l}}\right)$with the estimates of the parametric regression models, including the unknown parameters$\mathbf{\mathit{\beta}}$. The model‐based regression standardization estimator [1], or the parametric g‐formula for a fixed exposure effect [2], is given by

$$\hat{E} \left[Y \left(1\right)\right] - \hat{E} \left[Y \left(0\right)\right] = \frac{1}{n} \sum_{i = 1}^{n} R \left(1 , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}\right) - \frac{1}{n} \sum_{i = 1}^{n} R \left(0 , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}\right) ,$$

where$R \left(z , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}\right)$is the outcome predicted from the estimated parametric model for patient$i$by replacing the exposure status$Z_{i}$with$z = 0 , 1$.

The popular parametric model for the binary outcome$Y_{i}$is a logistic regression model. For example,

$$logit R \left(z , \mathbf{\mathit{l}} ; \mathbf{\mathit{\beta}}\right) = \beta_{0} + \beta_{1} z + \mathbf{\mathit{\beta}}_{2}^{\top} \mathbf{\mathit{l}} + \mathbf{\mathit{\beta}}_{3}^{\top} z \mathbf{\mathit{l}} ,$$

where$logit r = log \left\{r / \left(1 - r\right)\right\}$for$r \in \left(0 , 1\right)$. We typically estimate the parameters$\mathbf{\mathit{\beta}} = \left(\beta_{0} , \beta_{1} , \mathbf{\mathit{\beta}}_{2}^{\top} , \mathbf{\mathit{\beta}}_{3}^{\top}\right)^{\top}$by the maximum likelihood, which is equivalent to solving the following score equations:

$$\sum_{i = 1}^{n} \left\{Y_{i} - R \left(Z_{i} , \mathbf{\mathit{L}}_{i} ; \mathbf{\mathit{\beta}}\right)\right\} \mathbf{\mathit{X}}_{i} = 0 ,$$

where$\mathbf{\mathit{X}}_{i} = \left(1 , Z_{i} , \mathbf{\mathit{L}}_{i}^{\top} , Z_{i} \mathbf{\mathit{L}}_{i}^{\top}\right)^{\top}$is the shorthand notation for regressors in the logistic model. Using the resulting maximum likelihood estimator$\hat{\mathbf{\mathit{\beta}}}_{\text{ML}}$, the model‐based standardization estimator can be expressed as follows:

$$\begin{aligned} & \frac{1}{n} \sum_{i = 1}^{n} R \left(1 , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{ML}}\right) - \frac{1}{n} \sum_{i = 1}^{n} R \left(0 , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{ML}}\right) \\ & = \frac{1}{n} \sum_{i = 1}^{n} expit \left(\hat{\beta}_{0 , \text{ML}} + \hat{\beta}_{1 , \text{ML}} + \hat{\mathbf{\mathit{\beta}}}_{2 , \text{ML}}^{\top} \mathbf{\mathit{L}}_{i} + \hat{\mathbf{\mathit{\beta}}}_{3 , \text{ML}}^{\top} \mathbf{\mathit{L}}_{i}\right) \\ & - \frac{1}{n} \sum_{i = 1}^{n} expit \left(\hat{\beta}_{0 , \text{ML}} + \hat{\mathbf{\mathit{\beta}}}_{2 , \text{ML}}^{\top} \mathbf{\mathit{L}}_{i}\right) ,\end{aligned}$$

where$expit \left(u\right) = e^{u} / \left(1 + e^{u}\right)$is the inverse function of the$logit$function.

When the number of events is extremely small or large, the maximum likelihood estimates$\hat{\mathbf{\mathit{\beta}}}_{\text{ML}}$would suffer from sparse data bias or even cannot be obtained because of the (quasi‐)complete separation [6,8,9,11,12,13]. In such cases, the standardization estimators would also become invalid. Firth's method circumvents such a difficulty by penalizing the likelihood function of logistic models, corresponding to solving the following “penalized” score equations [14]:

$$\sum_{i = 1}^{n} \left\{Y_{i} - R \left(Z_{i} , \mathbf{\mathit{L}}_{i} ; \mathbf{\mathit{\beta}}\right) + h_{i} \left(\frac{1}{2} - R \left(Z_{i} , \mathbf{\mathit{L}}_{i} ; \mathbf{\mathit{\beta}}\right)\right)\right\} \mathbf{\mathit{X}}_{i} = 0 ,$$

where$h_{i} = h \left(\mathbf{\mathit{X}}_{i} , \mathbf{\mathit{\beta}}\right)$represents the$i$‐th diagonal element of the “hat” matrix$\mathbf{\mathit{H}}_{\mathbf{\mathit{\beta}}}$such that$\mathbf{\mathit{H}}_{\hat{\mathbf{\mathit{\beta}}}_{\text{FML}}} \left(Y_{1} , \ldots , Y_{n}\right)^{\top} = \left(R \left(Z_{1} , \mathbf{\mathit{L}}_{1} ; \hat{\mathbf{\mathit{\beta}}}_{\text{FML}}\right) , \ldots , R \left(Z_{n} , \mathbf{\mathit{L}}_{n} ; \hat{\mathbf{\mathit{\beta}}}_{\text{FML}}\right)\right)^{\top}$, producing predicted outcomes at the Firth's maximum penalized likelihood estimates satisfying the equations (6).

The regression standardization estimators are obtained simply by replacing$\hat{\mathbf{\mathit{\beta}}}_{\text{ML}}$with$\hat{\mathbf{\mathit{\beta}}}_{\text{FML}}$in (5). Although there is minimal bias in$\hat{\mathbf{\mathit{\beta}}}_{\text{FML}}$for$\mathbf{\mathit{\beta}}$in a correctly specified model [16], the standardization estimators based on Firth's penalized likelihood estimates may be biased owing to the known property of predicted outcome probabilities$R \left(Z_{i} , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{FML}}\right)$[14]. The predicted probabilities tend to be pulled toward 0.5, as described in the next section.

Discrepancy Between Standardization Estimators Following the Standard and Penalized Maximum Likelihood Estimation

Systematic Change in Predicted Probabilities Due to Firth's Penalization

Unlike the standard maximum likelihood method, the average predicted probability using Firth's method is not equal to the observed event rate [14]. Namely, (1) theaverageof the predicted probabilities$R \left(Z_{i} , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{FML}}\right)$is in between the average of predicted probabilities$R \left(Z_{i} , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{ML}}\right)$(which is equal to the event rate in a sample) and 0.5 and (2)eachpredicted probability$R \left(Z_{i} , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{FML}}\right)$tends to be pulled to 0.5 from$R \left(Z_{i} , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{ML}}\right)$[14,17]. To explain these results, we rewrite the penalized score equations (6) as follows:

$$\begin{aligned} & \sum_{i = 1}^{n} \left\{\left(Y_{i} - R \left(Z_{i} , \mathbf{\mathit{L}}_{i} ; \mathbf{\mathit{\beta}}\right)\right)\right) \\ & \left(+ \frac{h_{i}}{2} \left(1 - R \left(Z_{i} , \mathbf{\mathit{L}}_{i} ; \mathbf{\mathit{\beta}}\right)\right) + \frac{h_{i}}{2} \left(0 - R \left(Z_{i} , \mathbf{\mathit{L}}_{i} ; \mathbf{\mathit{\beta}}\right)\right)\right\} \mathbf{\mathit{X}}_{i} = 0 .\end{aligned}$$

This indicates that equations (6) augment the standard score equations (4) by introducing two artificial “copies”$\left(y , Z_{i} , \mathbf{\mathit{L}}_{i}\right)$of the data$\left(Y_{i} , Z_{i} , \mathbf{\mathit{L}}_{i}\right)$, where the outcomes are replaced with$y = 1$and$y = 0$, and weighted by$h_{i} / 2$. Note that$0 < h_{i} \leq 1$such that$\sum_{i = 1}^{n} h_{i} = dim \left(\mathbf{\mathit{\beta}}\right)$if the design matrix of a logistic model is full‐rank [14,17]. Hence, it is clear that the event rate in this “augmented” dataset, which is equal to average$R \left(Z_{i} , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{FML}}\right)$, is closer to 0.5 than the event rate in the original dataset, which is equal to average$R \left(Z_{i} , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{ML}}\right)$. Moreover, the above expression also provides insight that each predicted probability$R \left(Z_{i} , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{FML}}\right)$tends to be pulled toward 0.5 because one event ($y = 1$) and one non‐event ($y = 0$) with the same$\left(Z_{i} , \mathbf{\mathit{L}}_{i}\right)$are augmented with the same weight$h_{i} / 2$in the penalized equations (7). In fact, as illustrated in the next subsection, a simple condition exists that ensures that each$R \left(Z_{i} , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{FML}}\right)$is closer to 0.5 than$R \left(Z_{i} , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{ML}}\right)$if saturated regression models are used.

Hereafter, we assume the event rates to be much smaller than 0.5 in both exposure groups. Furthermore, because$\mathbf{\mathit{X}}_{i}$includes an exposure indicator$Z_{i}$, the relationship above for the average probabilities remains true for each exposure group, leading to the following inequalities:

$$\begin{aligned}\sum_{i = 1}^{n} Z_{i} R \left(Z_{i} , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{ML}}\right) & < \sum_{i = 1}^{n} Z_{i} R \left(Z_{i} , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{FML}}\right) < 0 . 5 , \\ \sum_{i = 1}^{n} \left(1 - Z_{i}\right) R \left(Z_{i} , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{ML}}\right) & < \sum_{i = 1}^{n} \left(1 - Z_{i}\right) R \left(Z_{i} , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{FML}}\right) < 0 . 5 .\end{aligned}$$

The standardization estimators also rely on the predicted outcomes under “opposite” exposure status, that is,$R \left(1 - Z_{i} , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}\right)$. The standardization estimator$\hat{E} \left[Y \left(z\right)\right]$based on Firth's penalized maximum likelihood estimator is

$$\begin{aligned} & \hat{E}_{\text{FML}} \left[Y \left(1\right)\right] = \frac{1}{n} \sum_{i = 1}^{n} R \left(1 , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{FML}}\right) \\ & = \frac{1}{n} \sum_{i = 1}^{n} \left\{Z_{i} R \left(Z_{i} , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{FML}}\right) + \left(1 - Z_{i}\right) R \left(1 - Z_{i} , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{FML}}\right)\right\} \\ & \hat{E}_{\text{FML}} \left[Y \left(0\right)\right] = \frac{1}{n} \sum_{i = 1}^{n} R \left(0 , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{FML}}\right) \\ & = \frac{1}{n} \sum_{i = 1}^{n} \left\{\left(1 - Z_{i}\right) R \left(Z_{i} , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{FML}}\right) + Z_{i} R \left(1 - Z_{i} , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{FML}}\right)\right\} \\ & \end{aligned}$$

The first terms in the summation of RHSs of these equations are larger than their counterparts based on the maximum likelihood estimates$\hat{\mathbf{\mathit{\beta}}}_{\text{ML}}$. Unfortunately, however, we do not have further general results regarding the relationship between standardization estimators with and without Firth's penalization because the average “counterfactual” predictions (i.e., the second terms) may or may not become greater by Firth's penalization than that from maximum likelihood estimates in each exposure group.

However, as the next subsection demonstrates, we can derive a simple sufficient condition for$R \left(z , \mathbf{\mathit{l}} ; \hat{\mathbf{\mathit{\beta}}}_{\text{FML}}\right)$to be larger than$R \left(z , \mathbf{\mathit{l}} ; \hat{\mathbf{\mathit{\beta}}}_{\text{ML}}\right)$for the given combinations of$\left(z , \mathbf{\mathit{l}}\right)$if we employ saturated logistic regression models. We defer the quantitative evaluation in more practical settings with parametric models, as well as the examination of whether the discrepancy between these standardization estimators leads to a bias, to the next section.

Illustration Using Stratified Data

Table1represents the data from the Society for Orthopedic Surgical Site Infection (OSSI), which consists of seven hospitals in the metropolitan area of Japan [18]. Briefly, the Society for OSSI prospectively collected the preoperative, intraoperative, and postoperative variables of consecutive patients undergoing orthopedic surgery in a clean surgical environment between 2013 and 2016. Overall, 10 943 eligible patients were enrolled, and 148 variables were measured.

Table: Smoking–surgical site infection (SSI) association stratified by American Society of Anesthesiologists (ASA) classification among patients undergoing arthropathy surgeries in the data obtained from the Society of OSSI (n=2717).

Prognostic factors for SSI events include age, sex, American Society of Anesthesiologists (ASA) classification of severity, and diabetes, with smoking being particularly important [19]. In the present study, we examined the association between smoking habits ($Z = 1$: smoker,$0$: non‐smoker) as an exposure and incidence within 30 days after surgery ($Y = 1$: occurrence,$0$: no occurrence) as an outcome, adjusting for ASA classification ($L = 1$:$\geq 3$,$0$:$\leq 2$) as a confounder. Postoperative complications, including SSI, are generally more common in smokers than in nonsmokers, and a previous report revealed that perioperative smoking cessation reduces SSIs [20]. In this database, most non‐arthritic surgeries required perioperative hospitalization, which led patients to stop smoking. Hence, Table1presents only$n = 2717$arthropathic surgeries as an illustrative example as these patients were not necessarily hospitalized before surgery.

We modeled the conditional probabilities for$Y = 1$in the stratified data in Table1using the saturated logistic model:

$$logit P \left(\right. Y = 1 \left|\right. Z = z , L = l \left.\right) = \beta_{0} + \beta_{1} z + \beta_{2} l + \beta_{3} z l .$$

Using the estimated conditional probabilities from the estimates of the saturated model, the maximum likelihood–based standardization estimates (3) yielded:

$$\begin{aligned}\hat{E}_{\text{ML}} \left[Y \left(1\right)\right] & = \frac{1}{6} \left(\frac{114}{2717}\right) + \frac{2}{189} \left(\frac{2603}{2717}\right) = 1 . 713 \% , \\ \hat{E}_{\text{ML}} \left[Y \left(0\right)\right] & = \frac{1}{108} \left(\frac{114}{2717}\right) + \frac{11}{2414} \left(\frac{2603}{2717}\right) = 0 . 475 \% .\end{aligned}$$

Although Firth's penalization can be applied to the saturated logistic model using statistical software, it is mathematically equivalent to obtaining the maximum likelihood estimates from a contingency table with pseudo‐counts of$+ 0 . 5$added to each cell of Table1[6]. As Firth's penalization affects only the estimation of conditional probabilities$R \left(z , l\right) = P \left(\right. Y = 1 \left|\right. Z = z , L = l \left.\right)$, the resulting standardization estimates are as follows:

$$\begin{aligned}\hat{E}_{\text{FML}} \left[Y \left(1\right)\right] & = \frac{1 . 5}{7} \left(\frac{114}{2717}\right) + \frac{2 . 5}{190} \left(\frac{2603}{2717}\right) = 2 . 160 \% , \\ \hat{E}_{\text{FML}} \left[Y \left(0\right)\right] & = \frac{1 . 5}{109} \left(\frac{114}{2717}\right) + \frac{11 . 5}{2415} \left(\frac{2603}{2717}\right) = 0 . 514 \% .\end{aligned}$$

Hence,$\hat{E}_{\text{FML}} \left[Y \left(z\right)\right] > \hat{E}_{\text{ML}} \left[Y \left(z\right)\right]$for both$z = 0 , 1$in this dataset.

There is a sufficient condition for Firth's penalization to inflate the maximum likelihood‐based standardization estimates using saturated regression models. As demonstrated above, the predicted outcome probability$R \left(z , l ; \hat{\mathbf{\mathit{\beta}}}_{\text{ML}}\right)$from the saturated regression models is equal to the exposure$z$‐specific observed event rates in confounder stratum$L = l$. We also demonstrated that Firth's penalization in the likelihood of the saturated regression model was equivalent to the maximum likelihood estimation applied to the augmented data, where all stratified cell counts are increased by 0.5. Let$N_{z , k}$denote the number of patients exposed to$Z = z$in the confounder stratum$L = l_{k}$($k = 1 , \ldots , K$; for example, in Table1,$K = 2$such that$l_{1} = 1$and$l_{2} = 0$), of which$A_{z , k}$patients develop events. The predicted outcome under$Z = z$following Firth's penalization$R \left(z , l_{k} ; \hat{\mathbf{\mathit{\beta}}}_{\text{FML}}\right)$is larger than the maximum likelihood prediction,$R \left(z , l_{k} ; \hat{\mathbf{\mathit{\beta}}}_{\text{ML}}\right)$if and only if$\left(A_{z , k} + 0 . 5\right) / \left(N_{z , k} + 1\right) - A_{z , k} / N_{z , k} > 0$, or$A_{z , k} < N_{z , k} - A_{z , k}$. In other words, this inequality suggests that the event rate in the exposure$Z = z$is lower than 0.5, in stratum$l_{k}$. Hence,$\hat{E}_{\text{FML}} \left[Y \left(z\right)\right]$is necessarily greater than$\hat{E}_{\text{ML}} \left[Y \left(z\right)\right]$if the observed event rates in the exposure$Z = z$are smaller than 0.5 across all strata$l_{1} , \ldots , l_{K}$.

Ad Hoc Modiications of Standardization Estimators

To address the overestimation (relative to the standard maximum likelihood) in regression standardization estimators following Firth's maximum penalized likelihood, we apply two modified algorithms proposed in the context of prediction modeling: Firth's logistic regression with intercept correction (FLIC) and added covariate (FLAC) methods [14]. Both methods are implemented using the R packagelogistf(https://cran.r‐project.org/web/packages/logistf/index.html) [21]. The following simple modifications match the average predicted probabilities with the observed event proportions:

Both the maximum likelihood fit for FLIC and the (weighted) maximum likelihood fit with an added dummy variable$G_{i , j}$with weight$W_{i , j}$in the augmented data for FLAC ensure the equivalence of the average predicted outcome and observed event rates in the original dataset [14]. To apply them to the model‐based regression standardization estimators, we replace$R \left(z , L_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{ML}}\right)$in Equation (5) with the corresponding outcome predictions:

$$\begin{aligned} & \hat{E}_{\text{FLIC}} \left[Y \left(z\right)\right] = \frac{1}{n} \sum_{i = 1}^{n} expit \left\{\hat{\beta}_{\text{FLIC}} + logit R \left(z , \mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\beta}}}_{\text{FML}}\right)\right\} , \\ & \hat{E}_{\text{FLAC}} \left[Y \left(z\right)\right] \\ & = \frac{1}{n} \sum_{i = 1}^{n} expit \left(\hat{\beta}_{0 , \text{FLAC}} + \hat{\beta}_{1 , \text{FLAC}} z + \hat{\mathbf{\mathit{\beta}}}_{2 , \text{FLAC}}^{\top} \mathbf{\mathit{L}}_{i} + \hat{\mathbf{\mathit{\beta}}}_{3 , \text{FLAC}}^{\top} z \mathbf{\mathit{L}}_{i}\right) .\end{aligned}$$

Table2depicts the results of the regression standardization based on the saturated logistic model presented in Table1. All Firth's estimates were obtained from logistic model fits using the R packagelogistf. While applying Firth's penalization can result in higher standardization estimates than those obtained using maximum likelihood estimation, FLIC and FLAC can mitigate this overestimation even when Firth's penalization is applied. However, in this analysis, we did not evaluate whether the deviation in Firth's method (2nd row) from the others should be considered a bias; this issue is addressed in the next section using simple simulations.

Table: Estimates of model‐based standardization in the stratified data from Table1using saturated logistic models with different fitting methods.

A Simulation Study

Data Generation

We evaluated the bias of model‐based regression standardization using distinct estimates$\hat{\mathbf{\mathit{\beta}}}$(i.e., maximum likelihood, Firth's penalized maximum likelihood, FLIC, and FLAC) for estimating the marginal causal risks$E \left[Y \left(z\right)\right]$($z = 0 , 1$), marginal risk difference$E \left[Y \left(1\right)\right] - E \left[Y \left(0\right)\right]$, and the log marginal risk ratio$log \left\{E \left[Y \left(1\right)\right] / E \left[Y \left(0\right)\right]\right\}$through a simple simulation study, assuming situations with rare events and the possibility of separation occurring. Each simulation was repeated 1000 times. The simulations were performed using R version 4.2.2.

For data$i = 1 , \ldots , n = \left(500 , 2500\right)$in each repetition, we generated 10‐dimensional latent variables from a multivariate normal distribution with mean vector fixed at zero and covariance matrix$\sum$, where diagonal elements were set to 1 and off‐diagonal elements were set to$\rho \in \left\{0 . 5 , 0 . 25 , 0\right\}$. These latent variables were then converted into binary variables$L_{i} = \left(L_{i 1} , \ldots , L_{i 10}\right)$using an indicator function for values greater than zero. Additional details of the data‐generating parameters for the exposure and outcome models are provided in theSupporting Information, with the full set of parameters listed inSupporting InformationTables1and2.

$$logit P \left(\right. Z = 1 \left|\right. \mathbf{\mathit{L}} \left.\right) = \alpha_{0} + \sum_{k = 1}^{10} \alpha_{k} L_{k}$$

and outcome follows the logistic model

$$logit P \left(\right. Y = 1 \left|\right. Z , \mathbf{\mathit{L}} \left.\right) = \beta_{0} + \sum_{k = 1}^{10} \beta_{k} L_{k} + \beta_{Z} Z ,$$

where$\beta_{Z}$was set as$log \left(0 . 8\right)$.

Exposure probability was set to$10 \%$, and the event rates were 20%, 10%, 5%, 3%, 1%, and 0.5%. For each event rate, the frequency of separation was controlled by the correlation between the latent normal variables behind the confounders (0, 0.25, and 0.5). We manipulated the intercepts$\alpha_{0}$and$\beta_{0}$to achieve the desired event rates.

As the analytical representation of true marginal risks and their contrasts (difference and log ratio) is not feasible, we numerically approximated the targeted values by directly generating$Y_{i} \left(1\right)$and$Y_{i} \left(0\right)$with a sample size of$n = 10 , 0 , 0$.

While our simulations were conducted under conditions that avoided separation in the propensity score model, we emphasize that the absence of separation does not a sufficient condition for satisfying the positivity assumption and/or ensuring adequate covariate overlap between the groups. Positivity violations can be either deterministic (precluding nonparametricidentifiability) or stochastic (affecting finite‐sampleestimability); separation in a propensity score model can be viewed as a manifestation of a lack of estimability, which may or may not result from (stochastic or deterministic) positivity violation [22]. Hence, in practical applications, the (deterministic) positivity assumption required for nonparametric identifiability should be assessed based on subject‐matter knowledge, and potential near‐violations of positivity should be diagnosed beyond checking for separation to ensure valid inference. To this end, the overlap of the covariates and estimated propensity score distributions between exposure groups should be examined visually and statistically [23]. Researchers should also be alert to sparsity in covariate strata, where random positivity violations can occur because of small sample sizes or high‐dimensional data. In practice, such sparsity may be identified by examining contingency tables of key categorical covariates, where empty or near‐empty cells may indicate limited support. In such cases, researchers may consider simple remedies such as collapsing categories or treating them as continuous covariates in regression and propensity score models. When inverse probability weights are used, the distribution of the resulting weights should also be examined. For example, estimated weights with very extreme values may indicate nonpositivity or model misspecification [24]. When extreme weights affect the stability of the analysis,ad hocmethods such as weight truncation may offer a better tradeoff between bias and precision. In our simulation, these issues were mitigated by design, but such diagnostics are essential when applying confounding‐adjustment methods to real‐world data.

Evaluated Methods and Performance Measures

Six estimation methods were compared. In addition to the maximum likelihood with or without the Firth method, including the FLIC and FLAC modifications, two propensity score‐based estimators were employed to reduce the number of adjustment variables in the outcome regression models. The first is the model‐based regression standardization estimator based on the logistic regression model, which adjusts for the estimated propensity scores instead of confounders$\mathbf{\mathit{L}}_{i}$. The second method is an inverse probability weighted (IPW) estimator. The maximum‐likelihood estimates of propensity scores were obtained through a logistic model for$e \left(\mathbf{\mathit{L}}\right) = P \left(\right. Z = 1 \left|\right. \mathbf{\mathit{L}} \left.\right)$, namely$e \left(\mathbf{\mathit{L}}_{i} ; \hat{\mathbf{\mathit{\alpha}}}_{\text{ML}}\right) = expit \left(\hat{\alpha}_{0 , \text{ML}} + \hat{\mathbf{\mathit{\alpha}}}_{1 , \text{ML}}^{\top} \mathbf{\mathit{L}}_{i}\right)$. Both propensity score‐based methods have been demonstrated to circumvent instability owing to sparse data problems in rare disease settings [25,26,27,28,29,30].

Separation was detected by monitoring the standard error estimates of the parameters as established in [31]. We used the R packagedetectseparationto determine whether separation occurs by specifyingmethod = “detect_separation”in theglmfunction [32]. In addition, if no event occurs in either exposure group, maximum likelihood estimates without Firth's penalization and propensity score‐based estimates cannot be obtained from such datasets. We recorded them as “Not Applicable (NA)” and did not include them in the summary measures.

We estimated the marginal risks$E \left[Y \left(1\right)\right]$and$E \left[Y \left(0\right)\right]$, along with their differences and log‐transformed ratios. The asymptotic standard errors were estimated using multivariate delta methods (for regression standardization estimators) or the (conservative) sandwich estimator (for IPW estimators). The estimates were evaluated using bias (i.e., estimated − true effect), relative bias (i.e., bias divided by true values inSupporting InformationTable3), Monte Carlo standard error (MCSE), mean estimated standard error (MESE), and coverage proportion of the 95% confidence interval. To assess the reliability of the simulation results, we evaluated Monte Carlo standard errors (MCSE) and mean estimated standard errors (MESE). Across all scenarios, MCSE values were small relative to the estimated effects, and MESE closely matched the empirical standard deviations, indicating that the stable and reliable inference. As the true effect is substantially small, the relative bias is reported in the main text to clarify the extent of bias against the true effect, whereas the bias itself is reported in theSupporting Information. The R code for the simulation experiments is provided in theSupporting Information.

Results

Among 1000 replications, separation occurred only in the scenario with an event rate of 0.005 under a sample size of 2500; the separation rates were$24 \%$(between‐covariate correlation: 0),$83 \%$(correlation: 0.25), and$90 \%$(correlation: 0.5), resulting in “NA” values for maximum‐likelihood estimation without Firth's penalization (seeSupporting InformationTable2). For propensity score‐based estimators, we excluded replications in which at least one cell in the$\left(Y , Z\right)$contingency table had zero counts (seeSupporting InformationTable2). Exclusions were negligible in most scenarios but increased sharply when events were rare. At an event rate of$0 . 5 \%$, exclusions ranged from$0 . 0 \%$to$49 . 3 \%$of replications (0–493 out of 1000), depending on the sample size and the between‐covariate correlation. Accordingly, PS‐based results were summarized using a scenario‐specific denominator of$1000 - m$, where$m$denotes the number of excluded replications.

Figures1,2,3,4shows the relative bias of each method, and the relative bias of the unadjusted risk estimates is shown inSupporting InformationFigure1. As expected, the relative bias was larger at smaller event rates for the regression standardization methods, especially for the maximum likelihood and Firth's estimates. This trend was not observed in the propensity score‐based methods. Increasing the correlation between covariates did not necessarily increase the bias for all methods.

Relative bias forE[Y(1)]andEY(0). The compared methods are regression standardization with maximum likelihood estimates (ML), ML not adjusted for confounding variables (Unadj), regression standardization following Firth's method (FML), regression standardization following Firth's method with corresponding modification (FLIC, FLAC), regression standardization with propensity score‐adjusted model (PS‐adj), and inverse probability weighted estimates (IPW).

Relative bias forE[Y(1)]−E[Y(0)]andlog(E[Y(1)])−log(E[Y(0)])(n=500). The compared methods are ML, Unadj, FML, FLIC, FLAC, PS‐adj, and IPW.

Relative bias forE[Y(1)]andEY(0). The compared methods are ML, Unadj, FML, FLIC, FLAC, PS‐adj, and IPW.

Relative bias forE[Y(1)]−E[Y(0)]andlog(E[Y(1)])−log(E[Y(0)])(n=2500). The compared methods are ML, Unadj, FML, FLIC, FLAC, PS‐adj, and IPW.

When we focused on rare‐event scenarios prone to separation (event rate 0.005 with nonzero between‐covariate correlation), the relative bias of$E \left[Y \left(z\right)\right]$at$n = 2500$ranged from$10 \%$to$45 \%$for maximum likelihood estimation and from$50 \%$to$95 \%$for Firth's method across$z = 0$and$z = 1$; FLIC and FLAC exhibited smaller ranges. A similar pattern was observed for the marginal risk differences. At$n = 500$, the qualitative patterns were consistent with those at$n = 2500$; the relative bias increased as the event rate decreased for regression‐standardization methods, particularly for maximum likelihood estimation and Firth's method, whereas FLIC and FLAC showed comparatively smaller biases. The ranges were wider owing to more frequent separation and, for propensity score‐based methods, exclusions based on zero cells in the$\left(Y , Z\right)$table; however, the overall conclusions regarding relative performance and coverage remained unchanged.

Relative biases in the propensity score‐based methods were generally larger than those in the four multivariate‐adjusted regression model‐based standardization methods. Specifically, propensity score‐adjusted regression standardization exhibited a non‐negligible bias, which may be explained by model misspecification due to including propensity scores rather than confounders in the outcome regression models. The IPW estimator provides relatively large biases, especially under large correlations between confounders. This is possible because correlation may have caused instability in the propensity score estimates, which is known to result in finite‐sample bias in weighted estimators. In contrast, propensity score‐based methods demonstrated no consistent monotonic trends across event rates. For propensity score‐based estimators, we excluded replications in which at least one$\left(Y , Z\right)$cell had zero counts, according to our pre‐specified rule (seeSupporting InformationTable2). As a result, summaries are reported using a scenario‐specific denominator of$1000 - m$, where$m$denotes the number of excluded replications. In rare‐event settings (event rate 0.005), such exclusions ranged from$0 . 0 \%$to$49 . 3 \%$of replications, depending on the sample size and the correlation structure.

The MCSE and MESE (seeSupporting InformationFigures2and3) were comparable between the methods across all scenarios. Across all scenarios and estimators, the ratios of MCSE to MESE were generally close to 1 (ranging approximately from 0.9 to 1.5), indicating that the estimated standard errors were reasonably well calibrated. No systematic pattern of severe miscalibration was observed. The ratios of MCSE to MESE are provided inSupporting InformationTables4to7for completeness. We defined “approximately nominal”$95 \%$coverage as$95 \% \pm 2 \times \text{MCSE}$(with 1000 replications, MCSE is approximately 0.7 percentage points, corresponding to a range of$93 . 6 \%$–$96 . 4 \%$). Under this criterion, the$95 \%$coverage rates of all methods, except for the maximum likelihood estimator, were nominal in most scenarios (Supporting InformationFigure3). The undercoverage of the maximum likelihood estimator was concentrated in settings with frequent separation, which reduced the yield of finite estimates (separation rates of up to$90 \%$in the rare‐event scenarios).

Analysis of the Data From the Society of OSSI

We analyzed the data from the Society for OSSI (Section3.2) in a more realistic manner by adjusting for multiple confounders using parametric models. Note that maximum likelihood estimates under the absence of events in at least one level of each categorical covariate may lead to non‐convergent estimates or produce “NA” results dependent on the “stopping” criteria of a fitting procedure used in statistical packages. We used theglmfunction in the R package to fit the models using maximum likelihood.

The outcome of interest was the occurrence of SSI within 30 days of surgery, with smoking as the exposure variable. We estimated the marginal causal risks$E \left[Y \left(1\right)\right]$and$E \left[Y \left(0\right)\right]$and their difference and ratio using the methods compared in the simulation experiments. Based on a prior publication from the OSSI study [18], we adjusted for the two distinct sets of covariates as linear terms. Adjustment Set 1 included age category (20–29, 30–39,…, 90–99,$\geq$100), gender, ASA classification ($\leq$2,$\geq$3), and presence or absence of diabetes. Adjustment Set 2 additionally included BMI category (< 25,25–< 30,$30 \leq$), total surgical time (< 60, 60–150,$150 \leq$min), postoperative drainage duration (none, < 48, 48 h or more), presence of rheumatoid arthritis, timing of prophylactic antibiotics (none, < 24, 24 h or more), highest postoperative blood glucose level (< 200 mg/dL,$\geq$200 mg/dL, not measured), surgical duration (hours), and blood loss (ml). These variables were selected based on their clinical relevance.Supporting InformationTable8presents the distribution of variables among smokers and nonsmokers.

Table3shows a series of estimates. The results after adjusting for Adjustment Set 1 are similar to those in Table2, which stratifies the data only by diabetes. Note that parametric modeling by logistic regression with linear terms circumvents the quasi‐complete separation of outcomes by model parameters, which is observed in the stratified data (i.e., the saturated model) in Table1. With Adjustment Set 2, the maximum likelihood estimates could not be obtained owing to quasi‐complete separation. Firth's penalized likelihood estimates provide larger estimates for both marginal risks,$E \left[Y \left(1\right)\right]$and$E \left[Y \left(0\right)\right]$. Correcting these estimates using FLIC and FLAC decreases risk by approximately$1 / 3 - 1 / 2$. The propensity score‐based IPW estimates provide a slightly lower risk under the exposed$E \left[Y \left(1\right)\right]$, resulting in a smaller risk ratio than that obtained by the other methods.

Table: Estimates of the marginal risk difference and ratio by adjusting for confounding variables by logistic regression model or propensity score model in arthritic patients from the Society of OSSI data (n=2717).

Discussion

Using real and simulated datasets, this study points out caveats in the application of Firth's method to regression standardization estimators, which have received little attention in the context of causal inference, and evaluates a simple remedy for them. In particular, Firth's method has a non‐negligible bias when estimating marginal counterfactual risks (and their difference) in situations where outcome separation occurs. These biases can be reduced byad hocmodifications such as FLIC and FLAC.

Our simulations revealed that in situations where outcome separation occurred, the regression standardization estimators following Firth's penalized likelihood estimation of logistic models yielded approximately$50 \% - - 66 \%$relative bias. These biases do not necessarily cancel out when estimating the risk differences or log risk ratios, as they may be introduced at different magnitudes in the exposed and unexposed groups. The biases were mitigated using FLIC and FLAC without inflating the standard error. The confidence intervals were also estimated sufficiently correctly using the typical multivariate delta method by replacing the parameter estimates with those from FLIC or FLAC. We also confirmed that the bias in regression standardization using Firth's method is negligible in situations where separation is absent, and the event proportion is not extremely low. However, we have not yet precisely determined the event‐rate threshold at which the FML begins to exhibit a notable bias. Moreover, various factors, such as exposure prevalence, outcome‐model misspecification, magnitude of the exposure effect, absolute number of events, and events per variable may influence this threshold. Therefore, future investigations should focus on the event‐rate threshold and these influencing factors.

In targeting the marginal log–marginal risk ratio as an effect measure, all regression standardization estimators, including Firth's method, suffer from a relative bias of up to 3%. In rare diseases, marginal risk ratios are approximated by marginal odds ratios, which are further approximated by the conditional odds ratios modeled in logistic regression. As Firth's method was originally proposed to reduce bias in typical maximum likelihood estimates of regression coefficients in logistic models [14,33], the marginal risk ratio can be estimated in an approximately unbiased manner by estimating an exposure coefficient, which represents the conditional odds ratio, with rare events.

In the primary simulations, we assumed that the regression models were correctly specified. However, additional robustness checks under two model misspecification settings yielded qualitatively similar conclusions. Because correctly specifying all working models is challenging in practical data analyses, the sensitivity of the methods to model misspecification remains an important topic for applied use. In practice, however, the problem of quasi‐complete separation of regression models may be mitigated by deliberate model misspecification by merging covariate categories, dropping higher‐order terms, omitting covariates that strongly predict exposure, and so froth. These tradeoffs between different sources of bias will be a topic for future studies.

The present study aimed to highlight the bias introduced into regression standardization by the easy‐to‐use Firth's method, we focused only on the limited approach for addressing (quasi‐)complete separation. For rare event analysis, a wide range of methods exists to handle “sparse‐data bias,” including Bayesian methods such as log‐F(1,1) prior [34] or the Cauchy prior [35], as well as penalized likelihood methods, including lasso or ridge. Moreover, “iterative” adjustment procedures for model fitting can enhance convergence stability [36]. In principle, these methods can be directly incorporated into regression standardization procedures. Evaluating potential bias and its relative performance in rare‐event settings is beyond the scope of this study.

Finally, a wide variety of models can be applied using Firth's method. For parametric models other than logistic models, an iterative algorithm can be used by adjusting the score functions [37]. In practice, we often choose alternative parametric models such as log‐linear Poisson or linear normal models, accompanied by robust sandwich variance estimators. With iterative fitting with Firth's correction of these alternative models, it is unknown whether the problems highlighted in this paper occur in rare event data and modifications such as FLIC/FLAC. Multinomial logit models, which utilize the framework of Poisson log‐linear models, have the potential to expand Firth's method beyond logistic regression [38].

Author Contributions

Conceptualization and methodology: Tomohiro Shinozaki; Code implementation: Hashibe S; Validation: Wataru Hongo; Writing – original draft: Sotaro Hashibe; Writing – review and editing: Tomohiro Shinozaki; Supervision:Tomohiro Shinozaki. All authors have read and approved the final manuscript.

Funding

This work was supported by the Japan Society for the Promotion of Science (Grant No. 24K14864) and the Japan Agency for Medical Research and Development (Grant No. JP25mk0121297).

Disclosure

The authors have nothing to report.

Conflicts of Interest

The authors declare no conflicts of interest.

Acknowledgments

We thank all researchers of the Society of OSSI, especially Dr. Koji Yamada (OrthoSupport/Nakanoshima Orthopedics), for their permission to use the data. This work was supported by JSPS KAKENHI Grant Number 24K14864 and AMED under Grant Number JP25mk0121297.

Data Availability Statement

Data sharing not applicable to this article as no datasets were generated or analysed during the current study.

Associated Data

Data Availability Statement

Data sharing not applicable to this article as no datasets were generated or analysed during the current study.

References

1. Greenland S., “Model‐Based Estimation of Relative Risks and Other Epidemiologic Measures in Studies of Common Outcomes and in Case‐Control Studies,” American Journal of Epidemiology 160, no. 4 (2004): 301–305. doi.org/10.1093/aje/kwh221
2. Hernán M. A. and Robins J. M., “Estimating Causal Effects From Epidemiological Data,” Journal of Epidemiology and Community Health 60, no. 7 (2006): 578–586. doi.org/10.1136/jech.2004.029496
3. Bartlett J. W., “Covariate Adjustment and Estimation of Mean Response in Randomised Trials,” Pharmaceutical Statistics 17, no. 5 (2018): 648–666. doi.org/10.1002/pst.1880
4. Rothman K. J., Greenland S., Lash T. L., et al., Modern Epidemiology, 3rd ed. (Wolters Kluwer Health/Lippincott Williams & Wilkins, 2008).
5. Albert A. and Anderson J. A., “On the Existence of Maximum Likelihood Estimates in Logistic Regression Models,” Biometrika 71, no. 1 (1984): 1–10.
6. Heinze G. and Schemper M., “A Solution to the Problem of Separation in Logistic Regression,” Statistics in Medicine 21, no. 16 (2002): 2409–2419. doi.org/10.1002/sim.1047
7. Van Smeden M., Groot d J A., Moons K. G., et al., “No Rationale for 1 Variable Per 10 Events Criterion for Binary Logistic Regression Analysis,” BMC Medical Research Methodology 16 (2016): 1–12. doi.org/10.1186/s12874-016-0267-3
8. Mansournia M. A., Geroldinger A., Greenland S., and Heinze G., “Separation in Logistic Regression: Causes, Consequences, and Control,” American Journal of Epidemiology 187, no. 4 (2018): 864–870. doi.org/10.1093/aje/kwx299
9. Heinze G., “A Comparative Investigation of Methods for Logistic Regression With Separated or Nearly Separated Data,” Statistics in Medicine 25, no. 24 (2006): 4216–4226. doi.org/10.1002/sim.2687
10. Firth D., “Bias Reduction of Maximum Likelihood Estimates,” Biometrika 80, no. 1 (1993): 27–38.
11. Heinze G., “The Application of Firth's Procedure to Cox and Logistic Regression,” Technical Report 10/1999, update in January 2001, Section of Clinical Biometrics, Department of Medical Computer Sciences, University of Vienna (2001).
12. Greenland S., Mansournia M. A., and Altman D. G., “Sparse Data Bias: A Problem Hiding in Plain Sight,” BMJ 352 (2016): i1981. doi.org/10.1136/bmj.i1981
13. Heinze G. and Puhr R., “Bias‐Reduced and Separation‐Proof Conditional Logistic Regression With Small or Sparse Data Sets,” Statistics in Medicine 29, no. 7–8 (2010): 770–777. doi.org/10.1002/sim.3794
14. Puhr R., Heinze G., Nold M., Lusa L., and Geroldinger A., “Firth's Logistic Regression With Rare Events: Accurate Effect Estimates and Predictions?,” Statistics in Medicine 36, no. 14 (2017): 2302–2317. doi.org/10.1002/sim.7273
15. Hernán M. A. and Robins J. M., CAUSAL INFERENCE What if (Chapman & Hall/CRC, 2020).
16. Ross R. K., Cole S. R., and Richardson D. B., “Decreased Susceptibility of Marginal Odds Ratios to Finite‐Sample Bias,” Epidemiology 32, no. 5 (2021): 648–652. doi.org/10.1097/EDE.0000000000001370
17. Kosmidis I. and Firth D., “Jeffreys‐Prior Penalty, Finiteness and Shrinkage in Binomial‐Response Generalized Linear Models,” Biometrika 108, no. 1 (2021): 71–82.
18. Yamada K., Nakajima K., Nakamoto H., et al., “Association Between Normothermia at the End of Surgery and Postoperative Complications Following Orthopedic Surgery,” Clinical Infectious Diseases 70, no. 3 (2020): 474–482. doi.org/10.1093/cid/ciz213
19. Mangram A. J., Horan T. C., Pearson M. L., Silver L. C., Jarvis W. R., and The Hospital Infection Control Practices Advisory Committee , “Guideline for Prevention of Surgical Site Infection, 1999,” Infection Control and Hospital Epidemiology 20, no. 4 (1999): 247–280. doi.org/10.1086/501620
20. Sørensen L. T., “Wound Healing and Infection in Surgery: The Clinical Impact of Smoking and Smoking Cessation: A Systematic Review and Meta‐Analysis,” Archives of Surgery 147, no. 4 (2012): 373–383. doi.org/10.1001/archsurg.2012.5
21. Heinze G., Ploner M., and Jiricka L., “logistf: Firth's Bias‐Reduced Logistic Regression,” R package version 1.24.1 (2022).
22. Zivich P. N., Cole S. R., and Westreich D., “Positivity: Identifiability and Estimability,” arXiv Preprint arXiv:2207.05010 (2022).
23. Imbens G. W. and Rubin D. B., Causal Inference for Statistics, Social, and Biomedical Sciences: An Introduction (Cambridge University Press, 2015).
24. Cole S. R. and Hernán M. A., “Constructing Inverse Probability Weights for Marginal Structural Models,” American Journal of Epidemiology 168, no. 6 (2008): 656–664. doi.org/10.1093/aje/kwn164
25. Braitman L. E. and Rosenbaum P. R., “Rare Outcomes, Common Treatments: Analytic Strategies Using Propensity Scores,” Annals of Internal Medicine 137, no. 8 (2002): 693–695. doi.org/10.7326/0003-4819-137-8-200210150-00015
26. Pirracchio R., Resche‐Rigon M., and Chevret S., “Evaluation of the Propensity Score Methods for Estimating Marginal Odds Ratios in Case of Small Sample Size,” BMC Medical Research Methodology 12 (2012): 1–10. doi.org/10.1186/1471-2288-12-70
27. Gagne J. J., Thompson L., O'Keefe K., and Kesselheim A. S., “Innovative Research Methods for Studying Treatments for Rare Diseases: Methodological Review,” BMJ 349 (2014): 349. doi.org/10.1136/bmj.g6802
28. Elze M. C., Gregson J., Baber U., et al., “Comparison of Propensity Score Methods and Covariate Adjustment: Evaluation in 4 Cardiovascular Studies,” Journal of the American College of Cardiology 69, no. 3 (2017): 345–357. doi.org/10.1016/j.jacc.2016.10.060
29. Cenzer I., Boscardin W. J., and Berger K., “Performance of Matching Methods in Studies of Rare Diseases: A Simulation Study,” Intractable & Rare Diseases Research 9, no. 2 (2020): 79–88. doi.org/10.5582/irdr.2020.01016
30. Almaghlouth I., Pullenayegum E., Gladman D. D., Urowitz M. B., and Johnson S. R., “Propensity Score Methods in Rare Disease: A Demonstration Using Observational Data in Systemic Lupus Erythematosus,” Journal of Rheumatology 48, no. 3 (2021): 321–325. doi.org/10.3899/jrheum.200254
31. Lesaffre E. and Albert A., “Partial Separation in Logistic Discrimination,” Journal of the Royal Statistical Society, Series B: Statistical Methodology 51, no. 1 (1989): 109–116.
32. Kosmidis I., Schumacher D., and Schwendinger F., “detectseparation: Detect and Check for Separation and Infinite Maximum Likelihood Estimates,” R Package Version 0.2 (2021).
33. Zhang J. and Kai F. Y., “What's the Relative Risk? A Method of Correcting the Odds Ratio in Cohort Studies of Common Outcomes,” JAMA 280, no. 19 (1998): 1690–1691. doi.org/10.1001/jama.280.19.1690
34. Greenland S. and Mansournia M. A., “Penalization, Bias Reduction, and Default Priors in Logistic and Related Categorical and Survival Regressions,” Statistics in Medicine 34, no. 23 (2015): 3133–3143. doi.org/10.1002/sim.6537
35. Gelman A., Jakulin A., Pittau M. G., and Su Y. S., “A Weakly Informative Default Prior Distribution for Logistic and Other Regression Models,” Annals of Applied Statistics 2, no. 4 (2008): 1360–1383.
36. Kosmidis I., “On Iterative Adjustment of Responses for the Reduction of Bias in Binary Regression Models,” CRiSM Working Paper Series 9, no. 36 (2009): 1–9.
37. Kosmidis I. and Firth D., “A Generic Algorithm for Reducing Bias in Parametric Estimation,” Electronic Journal of Statistics 4 (2010): 1097–1112.
38. Kosmidis I. and Firth D., “Multinomial Logit Bias Reduction via Poisson Log‐Linear Model,” Biometrika 98, no. 3 (2011): 755–759.

---

Republished from the open web under CC-BY. Authors: Hashibe S, Hongo W, Shinozaki T. Read the original.