Matching Weights Tutorial
Kazuki Yoshida
2016-07-26
References
Aim
This document provides a step-by-step guide for implementation of matching weight method in practice. The example is in the three-group setting. However, the essentially the same code can be used in the two-group setting or settings where there are more than three groups. The example is written in R, but it can be implemented in any statistical environment that has (multinomial) logistic regression and weighted data analysis capabilities.
Dataset
The tutoring dataset included in the TriMatch R package is used. The exposure is the treat variable, which takes one of Treat1, Treat2, and Control. These represent the tutoring method each student received. The outcome is the Grade ordinal variable, which takes one of 0, 1, 2, 3, or 4. Pre-treatment potential confounders include gender, ethnicity, military service status of the student, non-native English speaker status, education level of the subject’s mother (ordinal), education level of the subject’s father (ordinal), age of the student, employment status (no, part-time, full-time), household income (ordinal), number of transfer credits, grade point average. The dataset does not contain any missing values. See ?tutoring for details. The employment categorical variable is coded numerically. Thus, it is converted to a factor.
## Load data
library(TriMatch)
data(tutoring)
summary(tutoring)## treat Course Grade Gender Ethnicity Military ESL
## Control:918 Length:1142 Min. :0.000 FEMALE:627 Black:211 Mode :logical Mode :logical
## Treat1 :134 Class :character 1st Qu.:2.000 MALE :515 Other:193 FALSE:783 FALSE:1049
## Treat2 : 90 Mode :character Median :4.000 White:738 TRUE :359 TRUE :93
## Mean :2.891 NA's :0 NA's :0
## 3rd Qu.:4.000
## Max. :4.000
## EdMother EdFather Age Employment Income Transfer GPA
## Min. :1.000 Min. :1.000 Min. :20.00 Min. :1.000 Min. :1.000 Min. : 3.00 Min. :0.000
## 1st Qu.:3.000 1st Qu.:3.000 1st Qu.:30.00 1st Qu.:3.000 1st Qu.:3.000 1st Qu.: 36.66 1st Qu.:2.890
## Median :3.000 Median :3.000 Median :37.00 Median :3.000 Median :5.000 Median : 48.31 Median :3.215
## Mean :3.785 Mean :3.684 Mean :36.92 Mean :2.667 Mean :5.059 Mean : 52.12 Mean :3.166
## 3rd Qu.:5.000 3rd Qu.:5.000 3rd Qu.:43.00 3rd Qu.:3.000 3rd Qu.:7.000 3rd Qu.: 65.00 3rd Qu.:3.518
## Max. :8.000 Max. :9.000 Max. :65.00 Max. :3.000 Max. :9.000 Max. :126.00 Max. :4.000
## GradeCode Level ID
## Length:1142 Lower:988 Min. : 1.0
## Class :character Upper:154 1st Qu.: 286.2
## Mode :character Median : 571.5
## Mean : 571.5
## 3rd Qu.: 856.8
## Max. :1142.0## Make employment categorical
tutoring$Employment <- factor(tutoring$Employment, levels = 1:3,
labels = c("no","part-time","full-time"))Pre-weighting balance assessment
The tableone package can be utilized for covariate balance assessment using standardized mean differences (SMD). SMD greater than 0.1 is often regarded as a substantial imbalance. The SMD shown in the table is the average of all possible pairwise SMDs.
## Examine covariate balance
library(tableone)
covariates <- c("Gender", "Ethnicity", "Military", "ESL",
"EdMother", "EdFather", "Age", "Employment",
"Income", "Transfer", "GPA")
tab1Unadj <- CreateTableOne(vars = covariates, strata = "treat", data = tutoring)
print(tab1Unadj, test = FALSE, smd = TRUE)## Stratified by treat
## Control Treat1 Treat2 SMD
## n 918 134 90
## Gender = MALE (%) 449 (48.9) 38 (28.4) 28 (31.1) 0.287
## Ethnicity (%) 0.095
## Black 166 (18.1) 24 (17.9) 21 (23.3)
## Other 157 (17.1) 23 (17.2) 13 (14.4)
## White 595 (64.8) 87 (64.9) 56 (62.2)
## Military = TRUE (%) 309 (33.7) 32 (23.9) 18 (20.0) 0.208
## ESL = TRUE (%) 76 ( 8.3) 8 ( 6.0) 9 (10.0) 0.100
## EdMother (mean (sd)) 3.80 (1.49) 3.78 (1.51) 3.67 (1.54) 0.057
## EdFather (mean (sd)) 3.68 (1.65) 3.66 (1.73) 3.78 (1.73) 0.044
## Age (mean (sd)) 36.75 (8.95) 37.10 (9.41) 38.41 (9.49) 0.119
## Employment (%) 0.248
## no 95 (10.3) 24 (17.9) 18 (20.0)
## part-time 75 ( 8.2) 20 (14.9) 11 (12.2)
## full-time 748 (81.5) 90 (67.2) 61 (67.8)
## Income (mean (sd)) 5.10 (2.24) 5.04 (2.60) 4.69 (2.51) 0.111
## Transfer (mean (sd)) 51.40 (24.38) 57.37 (25.10) 51.61 (26.39) 0.158
## GPA (mean (sd)) 3.16 (0.58) 3.16 (0.46) 3.24 (0.58) 0.097## Examine all pairwise SMDs
ExtractSmd(tab1Unadj)## average 1 vs 2 1 vs 3 2 vs 3
## Gender 0.28718081 0.431825669 0.369462797 0.06025398
## Ethnicity 0.09475231 0.004540496 0.137619463 0.14209699
## Military 0.20773590 0.217301900 0.312032587 0.09387322
## ESL 0.09955894 0.089842245 0.059753148 0.14908142
## EdMother 0.05735067 0.014182489 0.086066827 0.07180268
## EdFather 0.04433253 0.007919139 0.059274560 0.06580389
## Age 0.11889003 0.038429226 0.179969129 0.13827175
## Employment 0.24838203 0.332479394 0.324590337 0.08807636
## Income 0.11113230 0.025003114 0.171951403 0.13644238
## Transfer 0.15777889 0.241327245 0.008454888 0.22355453
## GPA 0.09651297 0.009213587 0.128230886 0.15209444Propensity score modeling
As the exposure is a three-category variable, the propensity score model can be modeled using multinomial logistic regression. In R, the VGAM (vector generalized linear and additive models) package provides a flexible framework for this. Because the sample size of the treatment 2 group is small, making flexible modeling difficult, the ordinal variables are used only as linear terms. Predicting the “response” gives predicted probabilities of each treatment as a (sample size) \(\times\) 3 matrix, which then can be added to the dataset. The following AddGPS function can be used to ease this process. Three propensity scores (one for each treatment category) are added to the dataset.
## Function to add generalized PS to dataset
AddGPS <- function(data, formula, family = multinomial(), psPrefix = "PS_") {
library(VGAM)
## Fit multinomial logistic regression
resVglm <- vglm(formula = formula, data = data, family = family)
## Calculate PS
psData <- as.data.frame(predict(resVglm, type = "response"))
names(psData) <- paste0(psPrefix, names(psData))
cbind(data, psData)
}
tutoring <- AddGPS(data = tutoring, # dataset
## Propensity score model for multinomial regression
formula = treat ~ Gender + Ethnicity + Military +
ESL + EdMother + EdFather + Age +
Employment + Income + Transfer + GPA)Weight creation
As mentioned in the manuscript, the matching weight is defined as follows.
\[\begin{align*} MW_{i} &= \frac{{Smallest~ PS}}{{PS~ of~ assigned~ treatment}}\\ &= \frac{{min}(e_{1i},...,e_{Ki})}{\sum^K_{k=1} I(Z_{i} = k) e_{ki}}\\ \end{align*}\]where \(e_{ki}\) is the \(i\)-th individual’s probability of being assigned to the \(k\)-th treatment category given the covariate pattern, \(Z_i \in \{1,...,K\}\) is the categorical variable indicating the \(i\)-th individual’s treatment assignment.
The following function can be used to add matching weight to the dataset. Individuals’ matching weights have a range of [0,1], where as the inverse probability treatment weights have a range of [1,\(\infty\)].
## Function to add matching weight as mw to dataset
AddMwToData <- function(data, txVar, txLevels, psPrefix = "PS_") {
## Treatment indicator data frame (any number of groups allowed)
dfAssign <- as.data.frame(lapply(txLevels, function(tx_k) {
as.numeric(data[txVar] == tx_k)
}))
## Name of PS variables
psVars <- paste0(psPrefix, txLevels)
## Pick denominator (PS for assigned treatment)
data$PS_assign <- rowSums(data[psVars] * dfAssign)
## Pick numerator
data$PS_min <- do.call(pmin, data[psVars])
## Calculate the IPTW
data$iptw <- 1 / data$PS_assign
## Calculate the matching weight
data$mw <- exp(log(data$PS_min) - log(data$PS_assign))
## Return the whole data
data
}
## Add IPTW and MW
tutoring <- AddMwToData(data = tutoring, # dataset
txVar = "treat", # treatment variable name
tx = c("Control", "Treat1", "Treat2")) # treatment levels
## Check how weights are defined
head(tutoring[c("treat","PS_Control","PS_Treat1","PS_Treat2","PS_assign","PS_min","iptw","mw")], 20)## treat PS_Control PS_Treat1 PS_Treat2 PS_assign PS_min iptw mw
## 3 Control 0.8192816 0.11440448 0.06631388 0.81928164 0.06631388 1.220581 0.08094149
## 4 Control 0.8313205 0.10516348 0.06351606 0.83132046 0.06351606 1.202906 0.07640383
## 11 Control 0.6346235 0.22597339 0.13940309 0.63462352 0.13940309 1.575737 0.21966266
## 12 Control 0.7203265 0.11853269 0.16114082 0.72032649 0.11853269 1.388259 0.16455412
## 14 Control 0.6759314 0.15931947 0.16474916 0.67593137 0.15931947 1.479440 0.23570361
## 16 Treat1 0.7278386 0.18054526 0.09161616 0.18054526 0.09161616 5.538777 0.50744155
## 17 Control 0.7963014 0.09228518 0.11141339 0.79630143 0.09228518 1.255806 0.11589227
## 18 Control 0.7963014 0.09228518 0.11141339 0.79630143 0.09228518 1.255806 0.11589227
## 19 Control 0.4011609 0.29293705 0.30590201 0.40116094 0.29293705 2.492765 0.73022327
## 23 Control 0.7980564 0.14170696 0.06023666 0.79805638 0.06023666 1.253044 0.07547920
## 28 Treat2 0.7696177 0.11208565 0.11829667 0.11829667 0.11208565 8.453323 0.94749620
## 31 Treat1 0.7876534 0.11912070 0.09322587 0.11912070 0.09322587 8.394847 0.78261688
## 32 Control 0.7602112 0.13218394 0.10760486 0.76021120 0.10760486 1.315424 0.14154600
## 34 Treat2 0.6994628 0.12694918 0.17358797 0.17358797 0.12694918 5.760768 0.73132478
## 38 Treat1 0.6359332 0.24401948 0.12004734 0.24401948 0.12004734 4.098034 0.49195804
## 39 Control 0.7523881 0.15006473 0.09754713 0.75238814 0.09754713 1.329101 0.12965001
## 40 Control 0.8281320 0.11921012 0.05265789 0.82813199 0.05265789 1.207537 0.06358635
## 49 Treat1 0.7963180 0.09950924 0.10417277 0.09950924 0.09950924 10.049318 1.00000000
## 50 Control 0.8929612 0.06199434 0.04504442 0.89296124 0.04504442 1.119869 0.05044387
## 51 Control 0.6910650 0.16455995 0.14437500 0.69106505 0.14437500 1.447042 0.20891666## Check weight distribution
summary(tutoring[c("mw","iptw")])## mw iptw
## Min. :0.01025 Min. : 1.052
## 1st Qu.:0.05546 1st Qu.: 1.154
## Median :0.09410 Median : 1.258
## Mean :0.21706 Mean : 3.066
## 3rd Qu.:0.17721 3rd Qu.: 1.465
## Max. :1.00000 Max. :46.446Post-weighting balance assessment
All analyses afterward should be proceeded as weighted analyses. In R, this is most easily achieved by using the survey package. Firstly, a survey design object must be created with svydesign function. The resulting object is then used as the dataset. The weighted covariate table can be constructed with the tableone package. All SMDs are less than 0.1 after weighting, indicating better covariate balance.
## Created weighted data object
library(survey)
tutoringSvy <- svydesign(ids = ~ 1, data = tutoring, weights = ~ mw)
## Weighted table with tableone
tab1Mw <- svyCreateTableOne(vars = covariates, strata = "treat", data = tutoringSvy)
print(tab1Mw, test = FALSE, smd = TRUE)## Stratified by treat
## Control Treat1 Treat2 SMD
## n 82.8 82.6 82.5
## Gender = MALE (%) 24.9 (30.1) 25.0 (30.3) 24.4 (29.6) 0.010
## Ethnicity (%) 0.010
## Black 18.9 (22.9) 19.2 (23.3) 18.8 (22.8)
## Other 11.7 (14.1) 11.3 (13.7) 11.6 (14.1)
## White 52.2 (63.0) 52.1 (63.1) 52.0 (63.0)
## Military = TRUE (%) 17.2 (20.8) 19.7 (23.8) 17.4 (21.1) 0.048
## ESL = TRUE (%) 6.1 ( 7.4) 6.4 ( 7.7) 8.1 ( 9.8) 0.056
## EdMother (mean (sd)) 3.66 (1.49) 3.65 (1.47) 3.65 (1.55) 0.006
## EdFather (mean (sd)) 3.71 (1.70) 3.66 (1.75) 3.73 (1.70) 0.024
## Age (mean (sd)) 38.13 (9.68) 38.21 (9.63) 38.01 (9.38) 0.014
## Employment (%) 0.041
## no 16.3 (19.7) 15.6 (18.9) 15.2 (18.4)
## part-time 10.2 (12.3) 9.2 (11.2) 10.5 (12.7)
## full-time 56.3 (68.0) 57.7 (69.9) 56.8 (68.9)
## Income (mean (sd)) 4.76 (2.35) 4.72 (2.47) 4.80 (2.47) 0.023
## Transfer (mean (sd)) 52.46 (24.04) 51.39 (25.02) 53.48 (26.19) 0.055
## GPA (mean (sd)) 3.21 (0.49) 3.21 (0.45) 3.21 (0.59) 0.004## All pairwise SMDs
ExtractSmd(tab1Mw)## average 1 vs 2 1 vs 3 2 vs 3
## Gender 0.010336859 0.004393687 0.0111115330 0.0155053556
## Ethnicity 0.009595945 0.013881066 0.0006174629 0.0142893048
## Military 0.047738733 0.071609306 0.0067821033 0.0648247896
## ESL 0.055666107 0.010019487 0.0834804231 0.0734984115
## EdMother 0.005765913 0.008755059 0.0082762793 0.0002663992
## EdFather 0.023721214 0.024874520 0.0107632204 0.0355259006
## Age 0.013982735 0.008033386 0.0128645704 0.0210502478
## Employment 0.040896810 0.043102022 0.0330741322 0.0465142771
## Income 0.023351441 0.019691181 0.0157469189 0.0346162234
## Transfer 0.055073782 0.043293809 0.0406028456 0.0813246930
## GPA 0.003834104 0.006104611 0.0018132523 0.0035844491Visualizing the covariate balance before and after weighting can sometimes be helpful. Extracted SMD data can be fed to a plotting function (here ggplot2).
## Create SMD data frame
dataPlot <- data.frame(variable = rownames(ExtractSmd(tab1Unadj)),
Unadjusted = ExtractSmd(tab1Unadj)[,"average"],
Weighted = ExtractSmd(tab1Mw)[,"average"])
## Reshape to long format
library(reshape2)
dataPlotMelt <- melt(data = dataPlot,
id.vars = "variable",
variable.name = "method",
value.name = "SMD")
## Variables names ordered by unadjusted SMD values
varsOrderedBySmd <- rownames(dataPlot)[order(dataPlot[,"Unadjusted"])]
## Reorder factor levels
dataPlotMelt$variable <- factor(dataPlotMelt$variable,
levels = varsOrderedBySmd)
dataPlotMelt$method <- factor(dataPlotMelt$method,
levels = c("Weighted","Unadjusted"))
## Plot
library(ggplot2)
ggplot(data = dataPlotMelt, mapping = aes(x = variable, y = SMD, group = method, linetype = method)) +
geom_line() +
geom_point() +
geom_hline(yintercept = 0, size = 0.3) +
geom_hline(yintercept = 0.1, size = 0.1) +
coord_flip() +
theme_bw() + theme(legend.key = element_blank())
Outcome analysis
The outcome analyses should also be proceeded as weighted analyses. Most functions in the survey package is named svy* with * being the name of the unweighted counterpart.
The outcome was handled as a continuous outcome for simplicity. In weighted linear regression, both treatments appear superior to the control without tutoring regarding the course grade assuming the propensity score model was correctly specified. The mean difference was 0.45 [0.23, 0.67] for treatment 1 vs control and 0.67 [0.45, 0.89] for treatment 2 vs control.
## Weighted group means of Grade
svyby(formula = ~ Grade, by = ~ treat, design = tutoringSvy, FUN = svymean)## treat Grade se
## Control Control 2.792759 0.06648740
## Treat1 Treat1 3.244832 0.09179853
## Treat2 Treat2 3.463329 0.09070431## Group difference tested in weighted regression
modelOutcome1 <- svyglm(formula = Grade ~ treat, design = tutoringSvy)
summary(modelOutcome1)##
## Call:
## svyglm(formula = Grade ~ treat, design = tutoringSvy)
##
## Survey design:
## svydesign(ids = ~1, data = tutoring, weights = ~mw)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.79276 0.06649 42.004 < 2e-16 ***
## treatTreat1 0.45207 0.11335 3.988 0.00007076303 ***
## treatTreat2 0.67057 0.11246 5.963 0.00000000331 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 1.394533)
##
## Number of Fisher Scoring iterations: 2## ShowRegTable in tableone may come in handy
ShowRegTable(modelOutcome1, exp = FALSE)## coef [confint] p
## (Intercept) 2.79 [2.66, 2.92] <0.001
## treatTreat1 0.45 [0.23, 0.67] <0.001
## treatTreat2 0.67 [0.45, 0.89] <0.001Bootstrapping
As discussed in the manuscript, bootstrapping may provide better variance estimates than model-based inference. The boot package is a general purpose bootstrapping package. The following context-specific wrapper functions can be used to simplify the process. In this specific example, the bootstrap confidence intervals for the treatment effects were somewhat narrower.
## Define a function for each bootstrap step
BootModelsConstructor <- function(formulaPs, formulaOutcome, OutcomeRegFun, ...) {
## Obtain treatment variable name
txVar <- as.character(formulaPs[[2]])
## Return a function
function(data, i) {
## Obtain treatment levels
txLevels <- names(table(data[,txVar]))
## Add generalized propensity scores
dataB <- AddGPS(data = data[i,], formula = formulaPs)
## Add matching weight
dataB <- AddMwToData(data = dataB, txVar = txVar, txLevels = txLevels)
## Weighted analysis (lm() ok as only the estimates are used)
lmWeighted <- OutcomeRegFun(formula = formulaOutcome, data = dataB,
weights = mw, ...)
## Extract coefs
coef(lmWeighted)
}
}
## Define a function to summarize bootstrapping
BootSummarize <- function(data, R, BootModels, level = 0.95, ...) {
## Use boot library
library(boot)
## Run bootstrapping
outBoot <- boot(data = data, statistic = BootModels, R = R, ...)
out <- outBoot$t
colnames(out) <- names(outBoot$t0)
## Confidence intervals
lower <- apply(out, MARGIN = 2, quantile, probs = (1 - level) / 2)
upper <- apply(out, MARGIN = 2, quantile, probs = (1 - level) / 2 + level)
outCi <- cbind(lower = lower, upper = upper)
## Variance of estimator
outVar <- apply(out, MARGIN = 2, var)
outSe <- sqrt(outVar)
## Return as a readable table
cbind(est = outBoot$t0, outCi, var = outVar, se = outSe)
}
## Construct a custom bootstrap function with specific formulae
## formulaPs is propensity score model
BootModels <- BootModelsConstructor(formulaPs = treat ~ Gender + Ethnicity + Military +
ESL + EdMother + EdFather + Age +
Employment + Income + Transfer + GPA,
## Outcome model
formulaOutcome = Grade ~ treat,
## Regression function for outcome model
OutcomeRegFun = lm)
## Use a clean dataset without PS and weight variables
data(tutoring)
## Make employment categorical
tutoring$Employment <- factor(tutoring$Employment, levels = 1:3,
labels = c("no","part-time","full-time"))
## Run bootstrap
set.seed(201508131)
system.time(bootOut1 <- BootSummarize(data = tutoring, R = 2000, BootModels = BootModels))## user system elapsed
## 228.096 6.023 251.345bootOut1## est lower upper var se
## (Intercept) 2.7927593 2.6130814 2.9872607 0.008972568 0.09472364
## treatTreat1 0.4520730 0.2325361 0.6577786 0.011831058 0.10877067
## treatTreat2 0.6705692 0.4626595 0.8484488 0.009776627 0.09887683## Show naive confidence interval again
ShowRegTable(modelOutcome1, exp = FALSE, digits = 7)## coef [confint] p
## (Intercept) 2.7927593 [2.6624464, 2.9230722] <0.001
## treatTreat1 0.4520730 [0.2299169, 0.6742290] <0.001
## treatTreat2 0.6705692 [0.4501465, 0.8909920] <0.001- github: https://github.com/kaz-yos