4 Categorical exposure
4.1 Alcohol and Other Drug treatment data
require(twang)
data(AOD)
head(AOD)
table(AOD$treat)
#>
#> community metcbt5 scy
#> 200 200 200
Here,
- Outcome variable: suf12: substance use frequency at 12 month follow-up
- exposure variable: treat: possible values
- community: traditional programs (community)
- metcbt5: evidence-based MET/CBT-5 treatment protocol
- scy: Strengthening Communities for Youth
- Covariates (pre-treatment):
- illact: illicit activities scale
- crimjust: criminal justice involvement
- subprob: substance use problem scale
- subdep: substance use dependence scale
- white: non-Hispanic white youth
4.2 Step 1: Estimation of Propensity scores
We can fit vector generalized linear models (VGLMs) to fit propensity score model with our 3 category exposure variable:
require(VGAM)
<- "treat ~ illact + crimjust + subprob + subdep + white"
psFormula <- vglm(psFormula,family=multinomial, data=AOD)
ps.model $ps <- data.frame(fitted(ps.model)) AOD
require(ggplot2)
ggplot(AOD, aes(x=ps[,1], fill=factor(treat))) +
geom_density(alpha=0.5) +
scale_fill_discrete(name="Treatment Group") +
labs(title="Density Plot of Propensity Scores by Treatment Group",
x="Propensity Score",
y="Density") +
theme_minimal()
4.3 Step 1: Estimation of Propensity scores via GBM
set.seed(1235)
<- mnps(treat ~ illact + crimjust + subprob + subdep + white,
mnps.AOD.ATT data = AOD,
interaction.depth = 3,
estimand = "ATT",
treatATT = "community", # the treated
verbose = FALSE,
stop.method = "es.mean",
n.trees = 1000)
4.4 Step 2: IPW calculation for ATT
$w.ATT <- twang::get.weights(mnps.AOD.ATT, stop.method = "es.mean")
AODsummary(AOD$w.ATT)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.1494 0.7265 1.0000 0.8896 1.0000 1.8823
by(AOD$w.ATT, AOD$treat, summary)
#> AOD$treat: community
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 1 1 1 1 1 1
#> ------------------------------------------------------------
#> AOD$treat: metcbt5
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.1494 0.4939 0.7204 0.7529 0.9460 1.8386
#> ------------------------------------------------------------
#> AOD$treat: scy
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.3548 0.7619 0.9630 0.9159 1.0263 1.8823
4.5 Step 3: Balance assessment in weighted data
::bal.table(mnps.AOD.ATT,
twangdigits = 2,
collapse.to = "covariate")[,c("max.std.eff.sz",
"stop.method")]
plot(mnps.AOD.ATT, plots = 3)
4.6 Step 4: Effect estimates from weighted outcome model
require(survey)
<- svydesign(ids=~1, weights=~w.ATT, data=AOD)
design.mnps.ATT <- svyglm(suf12 ~ treat, design = design.mnps.ATT)
fit require(Publish)
publish(fit, intercept = FALSE)
#> Variable Units Coefficient CI.95 p-value
#> treat community Ref
#> metcbt5 0.20 [-0.00;0.41] 0.05129
#> scy 0.08 [-0.11;0.27] 0.41503