Estimate propensity scores via binary logistic regression

Fits stats::glm(..., family = binomial()) and appends the estimated propensity score, its logit, and a matching weight to the dataset. Supports multiply-imputed data in a stacked "long" format.

Usage

ps_logistic(
  formula,
  data,
  treatment_col = NULL,
  id_col = "id",
  imputation_col = NULL,
  score_col = "prob_t",
  logit_col = "logit_t",
  weight_col = "mt_wt",
  covariates = NULL,
  treatment_levels = NULL,
  treated_level = NULL
)

Arguments

formula: A formula with the binary treatment on the left-hand side, e.g. tavr ~ age + female + ef + diabetes.
data: A data frame. If imputation_col is set, data must be a stacked "long" MI data frame as produced by mice::complete(mids, action = "long").
treatment_col: Name of the binary treatment column: 0/1 or logical, or any two values when treated_level or treatment_levels names them. If NULL (default), extracted from the LHS of formula.
id_col: Patient identifier column. Required when imputation_col is set. Default "id".
imputation_col: Name of the imputation-index column in a stacked MI data frame. NULL (default) means a single complete dataset is supplied.
score_col: Output column name for the propensity score (probability of treatment). Default "prob_t".
logit_col: Output column name for the logit of the propensity score, log(p / (1 - p)). Default "logit_t".
weight_col: Output column name for the Li-Greene matching weight (ATM). Default "mt_wt".
covariates: Character vector of covariate column names for SMD balance diagnostics, used as given. If NULL (default), every numeric column that varies is used, except treatment_col, score_col, logit_col, weight_col, id_col and the fixed columns this package appends (quintile, decile, iptw, match, pair_id and the like). A predictor named in formula is always kept.
treatment_levels: Complete binary treatment levels. NULL takes them from the data: the two values present when treated_level is given, otherwise 0/1 or FALSE/TRUE.
treated_level: Level whose probability is the propensity score. When treatment_levels is supplied and this is NULL, its last value is used. Given alone, the treatment may hold any two values (a factor such as "TAVR"/"SAVR", for example) and this names the treated one.

Value

An object of class c("ps_logistic", "ps_data") with:

$data: data (the lowest-numbered imputation when imputation_col is set, without that column) with score_col, logit_col, weight_col, quintile, and decile columns appended. score_col is the probability of treated_level. quintile and decile rank score_col among patients with a treatment. weight_col holds the Li-Greene matching weight, min(p, 1 - p) / p for treated patients and min(p, 1 - p) / (1 - p) for controls, which targets the ATM estimand (average treatment effect among the matched; Li & Greene, 2013, doi:10.1515/ijb-2012-0030 ).
$meta: Named list: formula, treatment_col, id_col, imputation_col, score_col, logit_col, weight_col, treatment_levels, treated_level, bundle_version, model_family ("binary"), package_versions, method ("logistic" or "logistic-MI"), n_imputations, n_total, n_missing_treatment.
$tables: Named list: smd (variable, smd), group_counts, and the model tables of fit_logistic(): estimates, covariance, by_imputation and fit_status. group_counts has one row per treatment level, named as declared, with a logical treated flag and the count n.
$models: The fitted stats::glm() models, one per imputation, named by imputation index.

Details

Multiple imputation. When imputation_col is supplied, the model is fit separately on each imputed dataset, and each patient's predicted probabilities are averaged across the imputations on the probability scale, not the logit scale. logit_col is the logit of that average. Leyrat et al. (2019) call this the "across" approach. It is not Rubin's rules, which pool coefficients and their variances; $tables$estimates holds the coefficients pooled that way. The averaging mirrors the SAS template workflow of PROC LOGISTIC ... BY _IMPUTATION_ followed by PROC SUMMARY mean=.

$data is built from the lowest-numbered imputation, and its covariate values are the ones the smd balance table measures. That table reports balance in one imputed dataset, not across all of them.

In Leyrat et al.'s simulations of inverse probability of treatment weighting, only the "within" approach was unbiased in every scenario they studied, and they recommend it with the outcome included in the imputation model. To follow it, call ps_logistic() once per imputed dataset, without imputation_col; match or weight within each dataset; estimate the treatment effect in each; then pool those estimates and their variances with Rubin's rules.

The matching weight appended as weight_col is min(p, 1-p) / (p * trt + (1-p) * (1-trt)), which equals the mt_wt column produced in the SAS templates (Li & Greene, 2013). It targets the ATM estimand (average treatment effect among the matched). Neither ps_match() nor ps_weight() reads it; both work from score_col. Pass it as weight_col to ps_stddiff() or ps_mw_var(), or use it as the weight in an outcome model.

A patient with a missing treatment is left out of the model fit, as stats::glm() drops the row, and out of the strata, the balance tables and the group counts. The model needs no treatment to predict, so that patient still gets score_col and logit_col; weight_col, quintile and decile are missing. One warning states how many such patients there are, and $meta$n_missing_treatment keeps the count. A patient missing a predictor gets no score; $tables$fit_status counts those rows.

References

Leyrat C, Seaman SR, White IR, Douglas I, Smeeth L, Kim J, Resche-Rigon M, Carpenter JR, Williamson EJ (2019). Propensity score analysis with partially observed covariates: how should multiple imputation be used? Statistical Methods in Medical Research, 28(1), 3–19. doi:10.1177/0962280217713032

Li L, Greene T (2013). A weighting analogue to pair matching in propensity score analysis. The International Journal of Biostatistics, 9(2), 215–234. doi:10.1515/ijb-2012-0030

Examples

dta <- sample_ps_data(n = 200, seed = 42)
dta$prob_t <- NULL

# --- Single complete dataset (mirrors tp.lm.logistic_propensity_score.nomi.sas)
# Equivalent to: PROC LOGISTIC data=built descending; model tavr = ...;
obj <- ps_logistic(
  tavr ~ age + female + ef + diabetes + hypertension,
  data = dta
)
print(obj)
#> <ps_logistic>
#>   N total     : 400
#>   Treatment   : tavr
#>   PS column   : prob_t
#>   Weight col  : mt_wt
#>   Method      : logistic
#>   Tables      : smd, group_counts, estimates, covariance, by_imputation, fit_status 
summary(obj)
#> Summary of <ps_logistic>
#> 
#> Smd:
#>                  variable     smd
#> age                   age -0.7994
#> female             female  0.3989
#> ef                     ef  0.6139
#> diabetes         diabetes -0.0554
#> hypertension hypertension -0.0310
#> 
#> Group counts:
#>   group treated   n
#> 1     0   FALSE 200
#> 2     1    TRUE 200
#> 
#> Estimates:
#>           term    estimate  std.error  statistic  df      p.value    conf.low
#> 1  (Intercept)  4.24504327 1.30086827  3.2632384 Inf 1.101468e-03  1.69538831
#> 2          age -0.10638803 0.01575522 -6.7525593 Inf 1.452597e-11 -0.13726769
#> 3       female  0.75688483 0.23165843  3.2672449 Inf 1.085997e-03  0.30284264
#> 4           ef  0.06194245 0.01223150  5.0641747 Inf 4.101735e-07  0.03796915
#> 5     diabetes -0.06529637 0.25459818 -0.2564683 Inf 7.975892e-01 -0.56429962
#> 6 hypertension -0.01084919 0.23719636 -0.0457393 Inf 9.635180e-01 -0.47574552
#>     conf.high odds_ratio pooled
#> 1  6.79469822 69.7587783  FALSE
#> 2 -0.07550838  0.8990757  FALSE
#> 3  1.21092701  2.1316255  FALSE
#> 4  0.08591575  1.0639011  FALSE
#> 5  0.43370689  0.9367898  FALSE
#> 6  0.45404713  0.9892094  FALSE
#> 
#> Covariance:
#>               (Intercept)           age        female            ef
#> (Intercept)   1.692258249 -1.724556e-02  7.125339e-03 -6.768683e-03
#> age          -0.017245556  2.482268e-04 -3.290512e-04 -1.950182e-05
#> female        0.007125339 -3.290512e-04  5.366563e-02 -8.526194e-05
#> ef           -0.006768683 -1.950182e-05 -8.526194e-05  1.496096e-04
#> diabetes     -0.012736931 -1.199707e-04 -1.767623e-03  1.059612e-04
#> hypertension -0.039021927 -6.723910e-06 -1.699813e-03  1.024544e-04
#>                   diabetes  hypertension
#> (Intercept)  -0.0127369307 -3.902193e-02
#> age          -0.0001199707 -6.723910e-06
#> female       -0.0017676234 -1.699813e-03
#> ef            0.0001059612  1.024544e-04
#> diabetes      0.0648202310 -2.589114e-03
#> hypertension -0.0025891136  5.626211e-02
#> 
#> By imputation:
#>              imputation         term    estimate  std.error
#> (Intercept)           1  (Intercept)  4.24504327 1.30086827
#> age                   1          age -0.10638803 0.01575522
#> female                1       female  0.75688483 0.23165843
#> ef                    1           ef  0.06194245 0.01223150
#> diabetes              1     diabetes -0.06529637 0.25459818
#> hypertension          1 hypertension -0.01084919 0.23719636
#> 
#> Fit status:
#>   imputation converged n_input n_analyzed n_excluded
#> 1          1      TRUE     400        400          0
#> 

# The function appends:
#   prob_t   -- propensity score (p-hat; SAS _p_ / _propen_)
#   logit_t  -- log(p/(1-p))    (SAS _logit_)
#   mt_wt    -- matching weight  (SAS mt_wt = min(p,1-p)/(p*trt+(1-p)*(1-trt)))
#   quintile -- rank-based quintile, ceiling(rank * 5 / n)
#   decile   -- rank-based decile
head(obj$data[, c("id", "tavr", "prob_t", "logit_t", "mt_wt",
                  "quintile", "decile")])
#>   id tavr    prob_t    logit_t     mt_wt quintile decile
#> 1  1    0 0.1112484 -2.0780524 0.1251738        1      1
#> 2  2    0 0.5270896  0.1084644 1.0000000        3      6
#> 3  3    0 0.3951929 -0.4255358 0.6534196        2      4
#> 4  4    0 0.2836606 -0.9263756 0.3959863        2      3
#> 5  5    0 0.3751620 -0.5101344 0.6004149        2      4
#> 6  6    0 0.3148051 -0.7777498 0.4594387        2      3

# Pass the scored data to ps_match() for downstream matching
matched <- ps_match(obj$data, score_col = obj$meta$score_col, seed = 42)
nrow(matched$data[matched$data$match == 1L, ])
#> [1] 400

# --- Multiply-imputed data (mirrors tp.lm.logistic_propensity_score.sas)
# Equivalent to:
#   PROC LOGISTIC data=built descending; BY _IMPUTATION_; model tavr = ...;
#   PROC SUMMARY data=decile; class ccfid; var _p_; output out=... mean=_propen;
# \donttest{
# Simulate a stacked MI dataset (2 imputations, column "_Imputation_")
dta_mi <- rbind(
  cbind(dta, `_Imputation_` = 1L),
  cbind(dta, `_Imputation_` = 2L)
)
names(dta_mi)[names(dta_mi) == "_Imputation_"] <- "imp"

obj_mi <- ps_logistic(
  tavr ~ age + female + ef + diabetes + hypertension,
  data           = dta_mi,
  imputation_col = "imp",
  id_col         = "id"
)
print(obj_mi)
#> <ps_logistic>
#>   N total     : 400
#>   Treatment   : tavr
#>   PS column   : prob_t
#>   Weight col  : mt_wt
#>   Method      : logistic-MI (2 imputations)
#>   Tables      : smd, group_counts, estimates, covariance, by_imputation, fit_status 
# Per-patient PS is the average across the two imputed-dataset predictions,
# matching the PROC SUMMARY mean= step in the SAS template.
head(obj_mi$data[, c("id", "tavr", "prob_t")])
#>   id tavr    prob_t
#> 1  1    0 0.1112484
#> 2  2    0 0.5270896
#> 3  3    0 0.3951929
#> 4  4    0 0.2836606
#> 5  5    0 0.3751620
#> 6  6    0 0.3148051
# }