inst/doc/propensity.R

In propensity: A Toolkit for Calculating and Working with Propensity Scores

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(propensity)

## ----simulate-data------------------------------------------------------------
set.seed(42)
n <- 100
x1 <- rnorm(n)
x2 <- rnorm(n)
z <- rbinom(n, 1, plogis(0.5 * x1 + 0.3 * x2))
y <- rbinom(n, 1, plogis(-0.5 + 0.8 * z + 0.3 * x1 + 0.2 * x2))
dat <- data.frame(x1, z, y, x2)

## ----step1--------------------------------------------------------------------
ps_mod <- glm(z ~ x1 + x2, data = dat, family = binomial())

## ----step2--------------------------------------------------------------------
wts <- wt_ate(ps_mod)
outcome_mod <- glm(y ~ z, data = dat, family = binomial(), weights = wts)

## ----psw-inspect--------------------------------------------------------------
estimand(wts)
is_stabilized(wts)

## ----data-frame-input---------------------------------------------------------
ps <- fitted(ps_mod)
wt_ate(ps, dat$z)

## ----step3--------------------------------------------------------------------
result <- ipw(ps_mod, outcome_mod)
result

## ----switching-estimands------------------------------------------------------
wts_ate <- wt_ate(ps_mod)
wts_att <- wt_att(ps_mod)
wts_ato <- wt_ato(ps_mod)

## ----diagnose-weights---------------------------------------------------------
summary(wts_ate)

## ----overlap-estimands--------------------------------------------------------
summary(wt_ato(ps_mod))
summary(wt_atm(ps_mod))

## ----trim-ps------------------------------------------------------------------
ps_trimmed <- ps_trim(ps, method = "ps")

## ----trim-adaptive------------------------------------------------------------
ps_trimmed_adapt <- ps_trim(ps, method = "adaptive")

## ----trim-diagnostics---------------------------------------------------------
# Confirm the object has been trimmed
is_ps_trimmed(ps_trimmed)

# Which observations were removed?
sum(is_unit_trimmed(ps_trimmed))

# View trimming metadata (method, cutoffs, etc.)
ps_trim_meta(ps_trimmed)

## ----trim-subset--------------------------------------------------------------
retained <- !is_unit_trimmed(ps_trimmed)
dat_trimmed <- dat[retained, ]

## ----refit--------------------------------------------------------------------
ps_refitted <- ps_refit(ps_trimmed, ps_mod)
is_refit(ps_refitted)

## ----weights-from-refit-------------------------------------------------------
wts_trimmed <- wt_ate(ps_refitted, dat$z)
summary(wts_trimmed)

## ----truncate-----------------------------------------------------------------
ps_truncated <- ps_trunc(ps, lower = 0.05, upper = 0.95)

## ----trunc-diagnostics--------------------------------------------------------
is_ps_truncated(ps_truncated)
sum(is_unit_truncated(ps_truncated))
ps_trunc_meta(ps_truncated)

## ----weights-from-trunc-------------------------------------------------------
wts_truncated <- wt_ate(ps_truncated, dat$z)
summary(wts_truncated)

## ----interpret-binary---------------------------------------------------------
result

## ----as-data-frame------------------------------------------------------------
as.data.frame(result)

## ----exponentiate-------------------------------------------------------------
as.data.frame(result, exponentiate = TRUE)

## ----continuous-outcome-------------------------------------------------------
y_cont <- 2 + 0.8 * z + 0.3 * x1 + 0.2 * x2 + rnorm(n)
dat$y_cont <- y_cont
outcome_cont <- lm(y_cont ~ z, data = dat, weights = wts)
ipw(ps_mod, outcome_cont)

## ----continuous-exposure, eval = FALSE----------------------------------------
# # Fit a model for the continuous exposure
# ps_cont <- glm(continuous_exposure ~ x1 + x2, data = dat, family = gaussian())
# 
# # Stabilized weights (strongly recommended for continuous exposures)
# wts_cont <- wt_ate(ps_cont, stabilize = TRUE)

## ----categorical-exposure, eval = FALSE---------------------------------------
# # Multinomial propensity scores (one column per treatment level)
# ps_matrix <- predict(multinom_model, type = "probs")
# wt_ate(ps_matrix, exposure, exposure_type = "categorical")
# 
# # ATT and ATU require specifying the focal level
# wt_att(ps_matrix, exposure, .focal_level = "treated")

## ----calibrate, eval = FALSE--------------------------------------------------
# ps_calibrated <- ps_calibrate(ps, dat$z, method = "logistic", smooth = FALSE)
# is_ps_calibrated(ps_calibrated)
# 
# wts_calibrated <- wt_ate(ps_calibrated, dat$z)

## ----censoring-weights, eval = FALSE------------------------------------------
# # Model the probability of being uncensored
# cens_mod <- glm(uncensored ~ x1 + x2, data = dat, family = binomial())
# wts_cens <- wt_cens(cens_mod)
# 
# # Censoring weights use the same formula as ATE weights
# estimand(wts_cens) # "uncensored"