Estimate propensity scores via ordinal (cumulative logit) logistic regression

Fits a proportional-odds cumulative logit model via MASS::polr() for an ordered treatment variable and appends one probability column per treatment level to the dataset. The model is equivalent to the default PROC LOGISTIC behaviour for ordinal responses in SAS, standard errors included: they come from the expected information, as fit_logistic() describes under "Ordinal standard errors".

Usage

ps_ordinal(
  formula,
  data,
  treatment_col = NULL,
  id_col = "id",
  imputation_col = NULL,
  score_col_prefix = "prob",
  covariates = NULL,
  treatment_levels = NULL,
  smd_pairs = "reference"
)

Arguments

formula: A formula with the ordered factor treatment on the LHS, e.g. nyha_grp ~ age + female + ef. The treatment must be an ordered factor, whose level order is the order the model uses, unless treatment_levels gives that order. Any other treatment is an error: alphabetical order is not a safe guess for an ordinal scale.
data: A data frame (single or stacked MI).
treatment_col: Name of the ordinal treatment column. NULL (default) extracts from formula.
id_col: Patient identifier column. Required when imputation_col is set. Default "id".
imputation_col: Imputation-index column for stacked MI data. NULL (default) assumes a single complete dataset.
score_col_prefix: Prefix for the per-level output columns. Columns are named <prefix>_<level> for each level of the treatment. Default "prob".
covariates: Covariate columns for the balance table. NULL (default) uses every numeric column other than the treatment, score, identifier and strata columns.
treatment_levels: Complete treatment levels, lowest first. NULL (default) takes the levels of the treatment, which must then be an ordered factor.
smd_pairs: Which pairs of treatment levels the balance table compares. "reference" (default) compares each level with the first level; "adjacent" compares each level with the one below it; "all" compares every pair. Give more than one to combine them, for example c("reference", "adjacent").

Value

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

$data: data (the lowest-numbered imputation when imputation_col is set, without that column) with one <score_col_prefix>_<level> column per treatment level, the marginal probability of that level, and quintile and decile, which rank the probability of the highest level. The treatment column is returned as an ordered factor.
$meta: Named list: formula, treatment_col, id_col, imputation_col, score_cols, score_col_prefix, levels, smd_pairs, treatment_levels, bundle_version, model_family ("ordinal"), package_versions, cumulative_direction, method ("ordinal-logistic" or "ordinal-logistic-MI"), n_imputations, n_total.
$tables: Named list: smd, group_counts (group, n), and the model tables of fit_logistic(): estimates, covariance, by_imputation and fit_status. smd has one row per covariate and pair of levels: variable, level, versus, smd.
$models: The fitted models, one per imputation, named by imputation index: MASS::polr() fits, or stats::glm() fits for a two-level treatment.

Details

When imputation_col is supplied, models are fit on each imputed dataset separately and per-level probabilities are averaged per patient across imputations, on the probability scale. This is the "across" approach of Leyrat et al. (2019), not Rubin's rules; ps_logistic() describes it and the "within" alternative. $data and the smd table use the lowest-numbered imputation's covariates.

Package dependency

Requires the MASS package (a recommended package shipped with base R). Install with install.packages("MASS") if missing.

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

Examples

# Mirrors the SAS template workflow:
#   PROC LOGISTIC (cumulative logit, default for ordinal response)
#   followed by decomposition of cumulative to marginal probabilities
#   and rank-based quintile / decile assignment.
# \donttest{
if (requireNamespace("MASS", quietly = TRUE)) {
  dta <- sample_ps_data_ordinal(n = 300, seed = 42)
  obj <- ps_ordinal(
    nyha_grp ~ age + female + ef + diabetes,
    data = dta
  )
  print(obj)

  # Each level gets its own probability column (marginal, not cumulative).
  # The SAS template computes: p1=col1; p2=col2-col1; p3=1-col2.
  # ps_ordinal() performs this decomposition internally.
  head(obj$data[, c("id", "nyha_grp", "prob_I", "prob_II", "prob_III")])

  # Quintile and decile columns are appended, ordered by p(highest level).
  table(obj$data$quintile)
}
#> <ps_ordinal>
#>   N total     : 900
#>   Treatment   : nyha_grp (3 levels: I < II < III)
#>   Score cols  : prob_I, prob_II, prob_III
#>   Method      : ordinal-logistic
#>   Tables      : smd, group_counts, estimates, covariance, by_imputation, fit_status 
#> 
#>   1   2   3   4   5 
#> 180 180 180 180 180 
# }