Abstract

Background: Propensity score methods have become a cornerstone of modern causal inference, enabling researchers to approximate the conditions of randomized experiments in observational studies. Despite their widespread adoption, most established propensity score approaches were originally developed for two-group comparisons, leaving a notable methodological gap for multi-group data commonly encountered in clinical trials, public health interventions, and comparative effectiveness research. Methods: We conducted a comparative evaluation of several propensity score methods in balancing confounders and estimating treatment effects using Monte Carlo simulation. Datasets of varying sample sizes were generated under two distinct hybrid data-generating structures. Propensity scores were estimated using both generalized linear models (GLM) and generalized boosting models (GBM), and were subsequently applied via inverse probability of treatment weighting (IPTW), overlap weighting (OW), and matching. Five specific method combinations were evaluated: GLM-IPTW, GLM-OW, GLM-matching, GBM-IPTW, and GBM-OW. Covariate balance was assessed using standardized mean differences (SMD), while treatment effect estimation performance was evaluated based on point estimate accuracy and root mean square error (RMSE). Results: Across simulation scenarios with both linear and non-linear underlying relationships, the GLM-matching approach generally outperformed other methods. GLM-OW and GBM-OW demonstrated superior performance in achieving covariate balance, while GLM-IPTW and GBM-IPTW yielded more accurate point estimates of the treatment effect. Conclusion: When the relationship between covariates and outcome is relatively simple and treatment assignment follows a linear model, the GLM-matching method proved particularly advantageous. It produced estimates closer to the true value and exhibited a stronger ability to balance covariates compared to the other methods considered.

Keywords

Multi-Group Data, Generalized Boosting Model, Generalized Linear Model, Overlap Weighting

Share and Cite:

1. Background

Both observational studies and randomized controlled trials (RCTs) are essential methodologies in medical research. In RCTs, randomization is typically employed to assign participants to experimental or control groups, thereby minimizing confounding factors and ensuring balanced distribution of known and unknown confounders across groups in order to enable unbiased estimation of treatment effects. In contrast, observational studies are susceptible to selection bias due to the non-random allocation of subjects. To address the above issue, propensity score (PS) methods have emerged as a widely adopted approach for confounding adjustment in real-world studies [1]. By balancing the distributions of observed covariates across treatment groups, PS techniques reduce selection bias and strengthen the validity of effect estimates.

Propensity score methods have become a cornerstone in modern causal inference, enabling researchers to approximate randomized experiment conditions in observational settings. Several methods have been developed for propensity score estimation, such as generalized boosted models (GBM) and Generalized Linear Models (GLM) [2]. Once estimated, propensity scores can be incorporated into analyses through matching, weighting, or stratification to balance covariates across groups. Common implementations involve estimating propensity scores using either GBM or GLM, with the subsequent application of balancing techniques, such as Inverse Probability of Treatment Weighting (IPTW), Overlap Weighting (OW), or matching. Consequently, five distinct methodologies can be utilized: GBM with IPTW (GBM-IPTW), GBM with OW (GBM-OW), GLM with IPTW (GLM-IPTW), GLM with OW (GLM-OW), and GLM with matching (GLM-matching). Despite their widespread adoption, most established PS methodologies were tailored to two-group comparisons, leaving a significant gap for multi-group data commonly encountered in clinical trials, public health interventions, and comparative effectiveness research.

This study aims to address these methodological gaps by conducting a comprehensive comparison of five PS approaches in multi-group contexts. We will evaluate their relative performance in achieving covariate balance and producing accurate treatment effect estimates under varying conditions, including different sample sizes, functional forms of covariate-outcome relationships, and model specifications. These findings will contribute to the growing methodological literature on causal inference and provide evidence-based recommendations for applied researchers.

2. Methodology

2.1. Generalized Augmented Model (GBM)

The GBM estimates the propensity score p(x) by iteratively forming a collection of simple regression tree models and then combining them into a robust model [3].

$g\left(x\right)=\log \left(p\left(x\right)/\left(1-p\left(x\right)\right)\right)$ (1)

In the GBM model (Formula (1)), initial propensity scores are estimated using a logistic regression model to estimate the initial propensity scores g(x), and these initial propensity scores are used as the starting point for the GBM. Then, the adjustment function h(x) is found to add to g(x). In each iteration, the GBM model trains a new h(x) in this regression tree based on the current propensity score, observes the data, and adds it to the model. This process continues several times until a predetermined number of iterations is reached or a convergence condition is met [4]. In each iteration, the GBM weights the propensity score to update it based on the performance of the previous model and other factors. This process can help the model better adapt to the characteristics of the data and improve the accuracy of the estimation. The GBM is the algorithm for estimating the greatest likelihood of the logarithmic dominance of the propensity score values g(x) [5].

2.2. Generalized Linear Model (GLM)

The Generalized Linear Model (GLM) can be modeled in various ways to estimate propensity scores, and we used a polynomial logistic regression model in this study. When modeling the GLM, appropriate covariates need to be selected to predict the probability that an individual will be assigned to a treatment group, conditional on the specified covariates:

$e_j\left(x\right)=Pr\left(j|x\right)$ (2)

where x is the specified covariate and j = {1, 2, …, J} represents the different treatment groups (J > 2) [6]. A GLM is used to fit the data, with the fitting process determining the model parameters by maximizing the likelihood function or minimizing the loss function to obtain optimal parameter estimates [7]. Ordinary least squares (OLS) is generally used, or a machine algorithm is employed to find the optimal form of a weighted combination of several prediction algorithms through data-driven cross-validation [8], where the prediction algorithms include both parametric and nonparametric forms.

2.3. Weighted Methods

Weighting is the use of propensity scores to calculate weighting coefficients to weight the study population through different weighting methods in the study. The objects are assigned weights to make the confounders of the two groups converge. The weights are calculated using the formula [9].

$W=\frac{h\left(x\right)}{A\times e\left(x\right)+\left(1-A\right)\times \left(1-e\left(x\right)\right)}$ (3)

The inverse probability weighted target population includes all study subjects across different treatment groups, with the inverse probability weight being the inverse of the propensity score of the respective treatment groups, using h(x) = 1 in the weighting formula. The overlap weighting method, proposed by Li [10], defines the overlap weighting target population as a group with similar covariate characteristics. After identifying the overlap region between treatment groups, the density or probability of the propensity scores within this overlap region is used as the weight to adjust for individual subjects. Overlap weighting increases the relative weight (i.e., the ratio of individual sample weights to the overall weight) of the overlapping portion of the propensity score distribution [11].

2.4. Matching Method

Matching methods include nearest neighbor matching, caliper matching, and Mahalanobis distance matching. We used caliper matching in our study. In caliper matching, the difference in propensity scores of research subjects within each treatment group is matched within a predefined range, with this predefined range referred to as the caliper value [12]. Starting with the first subject in the first group, the subject from the second group with the closest propensity score is selected. If there is more than one subject in the second group with a propensity score identical to that of the first group’s subject, one is randomly chosen for matching. After matching the first two groups, the differences between the propensity score values of the third group’s research subjects and those of the first two groups are calculated separately, and the subject with the smallest sum of these differences is selected for matching [13]. The caliper value is set during the matching process to control the maximum allowable matching distance, thereby ensuring the accuracy of the matching.

2.5. Balance Analysis

Propensity score balance analysis is designed to assess whether the distributions of baseline characteristics are similar between the matched treatment and control groups, thereby ensuring the reliability and validity of subsequent comparisons. Statistically, this involves testing whether the differences in baseline characteristics between groups are significant. A common metric for this assessment is the Standardized Mean Difference (SMD) [14], which quantifies the magnitude of the difference between two group means, standardized by the pooled standard deviation. This standardization allows for the comparison of balance across variables measured on different scales. The SMD is calculated using the following formula:

$\text{SMD}=\frac{{\overline{X}}_1-\overline{X}}{{\text{SD}}_{\text{pooled}}}$(4)

where ${\overline{X}}_1$ and $\overline{X}$ correspond to the mean of group 1 and the overall mean, respectively, and ${\text{SD}}_{\text{pooled}}$ is the pooled standard deviation. A reduction in the SMD value reflects improved covariate balance. When extending this evaluation to multi-group data, the maximum absolute value of the SMD from all possible group pairs for a given covariate is used as the summary metric. This strategy is grounded in the logic that demonstrating balance between the most dissimilar groups provides sufficient evidence to infer balance among all groups, thereby justifying the pooling of data for subsequent analysis [15].

3. Simulation Study Design

A comprehensive simulation study was designed to evaluate and compare the performance of the five propensity score methodologies—GBM-IPTW, GBM-OW, GLM-IPTW, GLM-OW, and GLM-matching—under a range of conditions typically encountered in clinical research. The study was structured around the following key components:

3.1. Data Generation Mechanism

Multi-group clinical trial data were simulated with varying sample sizes (e.g., from 500 to 5000) to assess performance in both small-sample and large-sample settings. Baseline covariates were generated to include both continuous and binary variables. Crucially, the relationship between these covariates and the hypothetical outcome variable was manipulated to span two primary scenarios: a simple linear relationship and a more complex nonlinear relationship.

Firstly, six total sample sizes were established for the simulation: 500, 1000, 2000, 3000, 4000, and 5000. For each sample, an unordered trichotomous grouping variable T_j (j = 1, 2, 3) was randomly generated. The sample allocation ratios for Groups 1, 2, and 3 were set as pertr1, pertr2, and 1-pertr1-pertr2, respectively. In this study, these ratios were fixed at 1:2:7. Corresponding dummy variables for the three groups were then included in the subsequent models.

$j=\{\begin{array}{ll}1\hfill & x_1=1,x_2=0\hfill \\ 2\hfill & x_1=0,x_2=1\hfill \\ 3\hfill & x_1=0,x_2=0\hfill \end{array}$

Second, to emulate a clinical setting where baseline characteristics influence treatment assignment [16], we simulated eight covariates with imbalanced distributions across the treatment groups. This set comprised four continuous variables generated from normal distributions and four binary variables from binomial distributions. Furthermore, the continuous outcome variable was generated to depend on both these covariates and the treatment assignments. To comprehensively evaluate methodological performance, we specifically designed two distinct simulation scenarios.

Scene 1:

$\begin{array}{l}y=1.6+c_1\times 1.5+c_2\times 2+c_3\times 3+c_4\times 4\\ +b_1\times 4+b_2\times 3+b_3\times 2+b_4\times 1.5+x_1+x_2+b-\end{array}$

Scene 2:

$\begin{array}{l}y=0.5\times c_1\times c_2+c_2\times c_4+1.8\times c_3^2+c_4\times c_3\\ +b_1\times 4+b_2\times 3+b_3\times 2+b_4\times 1.5+x_1+x_2+ℰ\end{array}$

The error term $ℰ$ follows a normal distribution with a mean of 0 and a standard deviation of 2. Two distinct simulation scenarios were defined based on the relationship between the outcome and covariates: Scenario 1 assumed a simple linear relationship, whereas Scenario 2 incorporated a complex nonlinear relationship. Using the simulated dataset, the grouping variable was treated as the response, with all eight covariates included as independent variables in the propensity score model. Confounder adjustment was performed using each of the five propensity score methods. For the GLM-matching approach, 1:1:1 caliper matching without replacement was implemented, with a caliper width set to 0.1 [17]. The performance of alternative caliper widths was not explored and could be considered in future work. All simulations and analyses were conducted in the R software environment. To ensure robust performance evaluation, the entire data generation and analysis process was repeated 3000 times for each scenario, enabling a precise comparison of the five methods in handling multi-group propensity scores.

3.2. Performance Metrics

Following successful covariate balance via propensity score adjustment, the analysis proceeds analogously to a randomized controlled trial, wherein treatment effects are estimated by comparing the weighted means of the outcome across groups. To evaluate the performance of the five methods, we employed two primary classes of metrics. The first assessed covariate balance, quantified by the maximum absolute Standardized Mean Difference (SMD) across all pairwise group comparisons for each covariate. An SMD below the threshold of 0.1 was considered indicative of adequate balance [18]. The second assessed estimation accuracy, primarily measured using the Root Mean Square Error (RMSE) of the treatment effect estimate [19], as provided by the R software output. The RMSE, in conjunction with the point estimate, was used to determine each method’s ability to recover the true, pre-specified parameter value and to evaluate their relative strengths and weaknesses.

4. Results

4.1. Comparison of Effect Estimates of PS Methods under Scenario 1

Under a simple linear confounding structure (Scenario 1), all five methods produced similar Root Mean Square Error (RMSE) values. In terms of RMSE, GLM-matching exhibited the smallest variability in effect estimates across sample sizes, followed by GLM-OW and GBM-OW, while GLM-IPTW and GBM-IPTW showed notably higher variability (Table 1). Similarly, with respect to point estimation accuracy, GLM-matching provided the most precise estimates, with GLM-OW and GBM-OW also performing reasonably well, and GLM-IPTW and GBM-IPTW demonstrating larger deviations from the true value (Table 2). As sample size increased, the RMSE consistently decreased across all methods, and their point estimates progressively converged toward the true parameter value.

Table 1. RMSE of the five propensity scores under Scenario 1 (linear relationship only).

Table 2. Point estimation of the five propensity scores under Scenario 1 (linear relationship only).

Both evaluation metrics lead to a consistent conclusion: under a simple linear confounding structure, GLM-matching performed best, closely followed by GLM-OW and GBM-OW, whereas GLM-IPTW and GBM-IPTW were less accurate. All methods benefited from increased sample size, showing improved estimation accuracy as sample size grew (see Table 1 and Table 2).

4.2. Comparison of Effect Estimates of PS Methods under Scenario 2

Under a nonlinear relationship between covariates and the outcome variable (Scenario 2), GLM-matching, GLM-IPTW, and GBM-IPTW exhibited similar variability in effect estimates across all sample sizes, which was consistently lower than that of GLM-OW and GBM-OW (Table 3). As the sample size increased, the RMSE values of all five methods progressively declined. In terms of point estimation accuracy, GLM-matching outperformed the other methods under small sample sizes, followed by GLM-OW and GBM-OW, while GLM-IPTW and GBM-IPTW showed considerable deviation from the true treatment effect (Table 4). However, as the sample size grew, the accuracy of GLM-IPTW and GBM-IPTW improved substantially. At a sample size of 5000, GLM-IPTW achieved the highest precision, followed by GLM-OW and GLM-matching, whereas GBM-IPTW and GBM-OW displayed the largest deviations. Overall, increasing the sample size enhanced the estimation accuracy across all five methods (Table 3 and Table 4).

Table 3. RMSE of the five propensity scores under Scenario 2 (nonlinear relationships).

Table 4. Point estimation of the five propensity scores under Scenario 2 (nonlinear relationships).

4.3. Comparison of Covariate Equilibrium of PS Methods under Scenario 1

In large-sample studies, covariates are more likely to be balanced across groups. In contrast, with relatively small samples, achieving covariate balance—even after randomization—becomes challenging, making balance testing necessary to ensure comparability between groups. Standardized mean differences (SMD) were used to assess covariate balance in small-sample settings. This study specifically examined the performance of five covariate balancing methods under small-sample conditions.

When sample sizes were small (n ≤ 700), regardless of whether the outcome variable had a simple linear relationship or a complex nonlinear relationship with the covariate, both GLM-OW and GLM-matching outperformed the other three methods in balancing covariate effects, with all SMD values below 0.1—indicating that each covariate achieved adequate balance (Table 5 and Table 6).

Table 5. Equilibrium analysis of the five propensity scores under sample Scenario 1 (linear relationship only).

Table 6. Equilibrium analysis of the five propensity scores under sample Scenario 2 (nonlinear relationships).

In cases where only a linear relationship existed between the covariates and the outcome, GBM-OW demonstrated better covariate balance compared to GBM-IPTW and GLM-IPTW (Table 5). However, under more complex scenarios involving nonlinear relationships—particularly with smaller sample sizes—the performance of GBM-OW in balancing covariates was inferior to that of GBM-IPTW and GLM-IPTW (Table 6).

As sample size increased, the ability of GBM-OW, GBM-IPTW, and GLM-IPTW to balance covariates improved. Notably, the performance of GBM-OW gradually surpassed that of GBM-IPTW and GLM-IPTW. When the sample size exceeded 1000, all five methods achieved SMD values below 0.1 for each covariate, indicating satisfactory balance. For example, taking C₁ as an illustration, the SMD decreased with increasing sample size, and GBM-OW, GBM-IPTW, and GLM-IPTW performed similarly to GLM-OW and GLM-matching in balancing this covariate (Figure 1 and Figure 2).

Figure 1. Comparison of C1’s MD of five propensity score methods under Scenario 1.

Figure 2. Comparison of C1’SMD of five propensity score methods under Scenario 2.

5. Discussions

In recent years, propensity score methods have gained increasing popularity in clinical trial design for their ability to emulate randomized experiments and improve comparability between study groups. By addressing imbalances in the distribution of baseline covariates between treatment and control groups, these methods help reduce confounding and enhance causal inference. However, while the majority of existing applications focus on two-group comparisons, many real-world clinical trials involve multi-arm interventions, raising important methodological questions regarding the appropriate use of propensity score approaches in multi-group settings.

Existing PS estimation methods each present distinct advantages and limitations. Parametric approaches such as GLM rely on distributional assumptions that, when violated, could lead to model misspecification and biased estimates [20]. Machine learning alternatives like GBM offer flexibility in capturing complex relationships and handling data complexities, yet might introduce estimation instability [21] [22]. Recent innovations such as overlap weighting demonstrate superior theoretical properties in achieving exact covariate balance, though their practical performance across diverse multi-group settings warrants further investigation.

To address this limitation, data-driven approaches have been proposed, including the use of cross-validation to optimize weighted combinations of multiple prediction algorithms. Another alternative, the generalized boosted model (GBM), employed an iterative ensemble of regression trees. This method automatically selected covariates, captured complex nonlinear relationships, handled missing data directly, and resisted overfitting. However, in practice, GBM might still exhibit instability in weight estimation and yield imprecise effect estimates. In simpler clinical scenarios, both GLM and GBM demonstrated comparable performance in covariate balancing and estimation accuracy. With the methodological advance introduced in the overlap weighting (OW) approach, it could ensure exact balance of all covariate moments across groups, driving standardized mean differences toward zero. Moreover, OW was robust to extreme propensity scores and consistently achieved superior covariate balance compared to other weighting schemes, such as inverse probability of treatment weighting (IPTW) [23]. Our findings aligned with this theoretical advantage: both GBM-OW and GLM-OW yielded better covariate balance than their IPTW counterparts, reinforcing the value of overlap weighting in clinical trial settings where balance is prioritized. The inferior performance of IPTW in small samples or under nonlinear data-generating mechanisms could be attributed to two primary factors. First, model misspecification in the treatment model (e.g., using a GLM for a nonlinear relationship) lead to biased propensity scores and, consequently, biased effect estimates. Second, IPTW was particularly susceptible to extreme weights in these scenarios, which increased variance and can destabilize estimates [24].

When the relationship between covariates and the outcome was linear or mildly nonlinear, GLM-based matching emerged as a favorable option, offering optimal balance and estimation accuracy combined with analytical simplicity. In such cases, GLM-matching not only produced estimates closer to the true treatment effect but also maintained strong balancing properties. Meanwhile, GLM-OW and GBM-OW provided competitive balance performance, whereas GLM-IPTW and GBM-IPTW tended to yield more accurate estimates under correct model specification. Thus, in conventional clinical trial contexts with straightforward covariate-outcome relationships and well-specified treatment models, GLM-matching represented a straightforward and effective strategy for robust inference.

6. Limitations

This study had several limitations. First, the simulations considered only two confounders (related to both treatment assignment and the outcome) and a linear treatment-covariate relationship—idealized conditions under which GLM and GBM performed similarly. Real-world clinical data often involve more complex structures, where their performance might differ. Second, the degree of propensity score overlap was not varied, limiting insight into its influence, particularly in small samples. Third, only continuous outcomes were examined; performance with binary, time-to-event, or other outcome types remains to be evaluated. Third, the 1:2:7 allocation ratio was chosen to mimic a clinical trial with one smaller treatment arm. Different allocation ratios, especially those with more extreme imbalance, may affect the relative performance of the methods, particularly those reliant on weighting, and represent an important area for future research.

Funding

Research reported in this publication was supported by College Students’ Innovation and Entrepreneurship Training Program (2023, Southern Medical University, S202312121136).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References