15 October, 2015

## Matching Packages

MatchIt and optmatch implement matching methods, RItools implements covariate balance checking methods.

library(MatchIt)
library(optmatch)
library(RItools)
library(haven)
library(dplyr)
library(effsize)
library(kSamples)
library(weights)

## Propensity score matching:

What is it?

A method for causal inference in observational studies whereby observations are matched on known covariates to create 'treatment' and 'control' groups that simulate the effects of random assignment (with respect to measured covariates).

What types of problems is it useful for?

Causal inference problems where random assignment to treatment groups is infeasible/impossible:

Examples:

• Honors college vs non-honors college students (Honors college admissions based on covariates)
• Smokers vs Non-smokers (Can't assign some people to smoke for decades)
• The effectiveness of alcoholics anonymous (Self-selection bias)

## What are propensity scores?

Propensity scores are the conditional probability of exposure to treatment given observed covariates.

In random assignment: All individuals have a value of .5 for their propensity scores by definition Their covariates are not useful for predicting treatment assignment.
This is the benefit of random assignment, that covariate distributions are balanced between treatment and control groups.

In Observational studies: It is infeasible to assign individuals to treatments Propensity scores for different covariates and combinations thereof are not equal between groups. Causal inference is inhibited by unequal covariate balance; covariates may be related to outcomes.

Matching on propensity scores balances observed covariates in the manner that would be expected from random treatment assignment. Ability to draw causal inferences on the basis of observational data is improved.

## Propensity score estimation

How are propensity scores estimated?

Logistic regression of treatment assignment on the chosen covariates is used to determine the likelihood of selection to treatement based on the observed covariates.

Matching methods also use different distance metrics:

• Mahalanobis distance
• Euclidean distance
• Any other specifiable distance metric

Propensity score calipers often used to control quality of matches when using other distance metrics

## Types of propensity score matching techniques:

• Exact Matching (control unit is exactly the same as treatment; often not plausible)
• Nearest neighbor / Pairmatch
• K nearest neighbor / Fullmatch
• Genetic Matching
• Propensity Score Weighting

## Options for PSA in R:

Matchit (easiest to use, most user-friendly package, but nevertheless offers more advanced options, uses greedy matching algorithm by default)

From Matchit website:

Matchit enables parametric models for causal inference to work better by selecting well-matched subsets of the original treated and control groups. MatchIt implements the suggestions of Ho, Imai, King, and Stuart (2004). MatchIt implements a wide range of sophisticated matching methods, making it possible to greatly reduce the dependence of causal inferences on statistical modeling assumptions. After preprocessing with MatchIt, researchers can use whatever parametric model they would have used, but produce inferences with substantially more robustness and less sensitivity to modeling assumptions.

Optmatch (less user-friendly, offers more accessibility to functions and customizability, uses optimal matching algorithm by default)

From optmatch documentation:

Provides routines for distance based bipartite matching to reducecovariate imbalance between treatment and control groups in observational studies. Routines are provided to generate distances from GLM models (propensity score matching) and formulas (Euclidean and Mahalanobis matching), stratified matching (exact matching), and calipers. Results of the fullmatch routine are guaranteed to provide minimum average within matched set distance.

## 1-to-1 nearest neighbor matching with Matching; An Example:

Dataset: ~22000 students from large university, of whom 266 part of ASP program. Data on everything from entering characteristics to end of college outcomes. Covariates determine ASP program participation.

modeldata <- data.frame(hcdata$ASPSCLR, hcdata$Best_ACT_Comp,
hcdata$HSGPA, hcdata$incomingcreditstotal,
hcdata$gndr_flag, hcdata$ID,
hcdata$CumGPAendmostrecentUNterm, hcdata$@12yrretention)

colnames(modeldata)[1] <- "ASPSCLR"
colnames(modeldata)[5] <- "gender"
modeldataOmitNA <- na.omit(modeldata)

ppty <- glm(as.integer(ASPSCLR) ~ hcdata.Best_ACT_Comp + hcdata.HSGPA
+ hcdata.incomingcreditstotal + gender,
data = modeldataOmitNA, family = "binomial")
modeldataOmitNA$pptyScore <- ppty$fitted.values

## Nearest Neighbor Matching.

matchit function is the main function for matching.

Arguments include:

• method: matching method ("exact", "genetic", "nearest", "optimal")
• distance: distance metric ("logit", "mahalanobis", "euclidean")
matches <- matchit(ASPSCLR ~ hcdata.Best_ACT_Comp + hcdata.HSGPA +
hcdata.incomingcreditstotal + gender,
data = modeldataOmitNA)

matchedPairs <- match.data(matches)

match.data function creates new data frame with only the treatment and matched control observations, measured propensity score, distance metric (if different), and weights.

## Pair Matching Results

summary(matches)
##
## Call:
## matchit(formula = ASPSCLR ~ hcdata.Best_ACT_Comp + hcdata.HSGPA +
##     hcdata.incomingcreditstotal + gender, data = modeldataOmitNA)
##
## Summary of balance for all data:
##                             Means Treated Means Control SD Control
## distance                           0.0422        0.0133     0.0265
## hcdata.Best_ACT_Comp              28.4774       24.7969     3.6854
## hcdata.HSGPA                       3.9797        3.6045     0.3274
## hcdata.incomingcreditstotal        8.5000        3.9776     7.7014
## gender                             0.5940        0.5713     0.4949
##                             Mean Diff eQQ Med eQQ Mean eQQ Max
## distance                       0.0289  0.0263   0.0275  0.2053
## hcdata.Best_ACT_Comp           3.6805  3.0000   3.7895 14.0000
## hcdata.HSGPA                   0.3752  0.3310   0.3769  1.9670
## hcdata.incomingcreditstotal    4.5224  5.0000   4.6353 34.0000
## gender                         0.0227  0.0000   0.0226  1.0000
##
##
## Summary of balance for matched data:
##                             Means Treated Means Control SD Control
## distance                           0.0422        0.0415     0.0471
## hcdata.Best_ACT_Comp              28.4774       28.4962     2.5167
## hcdata.HSGPA                       3.9797        3.9644     0.1742
## hcdata.incomingcreditstotal        8.5000        7.1955    10.2121
## gender                             0.5940        0.6278     0.4843
##                             Mean Diff eQQ Med eQQ Mean eQQ Max
## distance                       0.0008   0.000   0.0008  0.2053
## hcdata.Best_ACT_Comp          -0.0188   1.000   1.1241  4.0000
## hcdata.HSGPA                   0.0153   0.011   0.0251  0.5450
## hcdata.incomingcreditstotal    1.3045   3.000   2.3271 15.0000
## gender                        -0.0338   0.000   0.0338  1.0000
##
## Percent Balance Improvement:
##                             Mean Diff. eQQ Med eQQ Mean eQQ Max
## distance                       97.3080 99.9864  97.0784  0.0000
## hcdata.Best_ACT_Comp           99.4893 66.6667  70.3373 71.4286
## hcdata.HSGPA                   95.9090 96.6767  93.3366 72.2928
## hcdata.incomingcreditstotal    71.1544 40.0000  49.7972 55.8824
## gender                        -49.1617  0.0000 -50.0000  0.0000
##
## Sample sizes:
##           Control Treated
## All         19158     266
## Matched       266     266
## Unmatched   18892       0
## Discarded       0       0

## Covariate Balance: before/after

Before

RItools::xBalance(as.numeric(ASPSCLR) ~ hcdata.Best_ACT_Comp +
hcdata.HSGPA + hcdata.incomingcreditstotal +
gender, data = modeldataOmitNA)
##                             strata  unstrat
##                             stat   std.diff     z
## vars
## hcdata.Best_ACT_Comp               1.00e+00 1.62e+01 ***
## hcdata.HSGPA                       1.15e+00 1.85e+01 ***
## hcdata.incomingcreditstotal        5.86e-01 9.48e+00 ***
## gender                             4.58e-02 7.42e-01

After

RItools::xBalance(as.numeric(ASPSCLR) ~ hcdata.Best_ACT_Comp +
hcdata.HSGPA + hcdata.incomingcreditstotal +
gender,
data = matchedPairs)
##                             strata  unstrat
##                             stat   std.diff     z
## vars
## hcdata.Best_ACT_Comp               -0.00947 -0.10935
## hcdata.HSGPA                       0.08974  1.03489
## hcdata.incomingcreditstotal        0.13889  1.59945
## gender                             -0.06931 -0.79958

## Checking Covariate Balance: QQ-Plots

plot(matches)

## Checking Covariate Balance: Jitter Plots

plot(matches, type = "jitter")

## [1] "To identify the units, use first mouse button; to stop, use second."
## integer(0)

## Checking Covariate Balance: Histograms

plot(matches, type = "hist")

## Comparing propensity score distributions

Anderson-Darling test for distributions of propensity scores

asp <- dplyr::filter(matchedPairs, ASPSCLR == 1)
nonasp <- dplyr::filter(matchedPairs, ASPSCLR == 0)
allcont <- dplyr::filter(modeldataOmitNA, ASPSCLR == 0)

ad.test(asp$pptyScore, allcont$pptyScore)
##
##
##  Anderson-Darling k-sample test.
##
## Number of samples:  2
## Sample sizes:  266, 19158
## Number of ties: 4912
##
## Mean of  Anderson-Darling  Criterion: 1
## Standard deviation of  Anderson-Darling  Criterion: 0.76163
##
## T.AD = ( Anderson-Darling  Criterion - mean)/sigma
##
## Null Hypothesis: All samples come from a common population.
##
## version 1: 293.68 384.28      1.8001e-161
## version 2: 294.00 384.30      6.2751e-161

## Anderson-Darling test for distributions of propensity scores

ad.test(asp$pptyScore, nonasp$pptyScore)
##
##
##  Anderson-Darling k-sample test.
##
## Number of samples:  2
## Sample sizes:  266, 266
## Number of ties: 64
##
## Mean of  Anderson-Darling  Criterion: 1
## Standard deviation of  Anderson-Darling  Criterion: 0.75865
##
## T.AD = ( Anderson-Darling  Criterion - mean)/sigma
##
## Null Hypothesis: All samples come from a common population.
##
## version 1: 0.01454 -1.2990                1
## version 2: 0.01030 -1.3045                1

## Follow-up analysis: Outcomes

T-test: Cumulative GPA, ASP vs non-ASP, unmatched data

t.test(hcdata.CumGPAendmostrecentUNterm ~ ASPSCLR, data = modeldataOmitNA)
##
##  Welch Two Sample t-test
##
## data:  hcdata.CumGPAendmostrecentUNterm by ASPSCLR
## t = -17.513, df = 286.19, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.4675145 -0.3730428
## sample estimates:
## mean in group 0 mean in group 1
##        3.097036        3.517315

T-test: Cumulative GPA, ASP vs non-ASP, matched data

t.test(hcdata.CumGPAendmostrecentUNterm ~ ASPSCLR, data = matchedPairs)
##
##  Welch Two Sample t-test
##
## data:  hcdata.CumGPAendmostrecentUNterm by ASPSCLR
## t = -0.47014, df = 513.16, p-value = 0.6385
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.08957818  0.05498344
## sample estimates:
## mean in group 0 mean in group 1
##        3.500017        3.517315

## Effect Sizes for mean differences

cohen.d(hcdata.CumGPAendmostrecentUNterm ~ ASPSCLR, modeldataOmitNA)
##
## Cohen's d
##
## d estimate: -0.6542157 (medium)
## 95 percent confidence interval:
##        inf        sup
## -0.7754087 -0.5330226
cohen.d(hcdata.CumGPAendmostrecentUNterm ~ ASPSCLR, matchedPairs)
##
## Cohen's d
##
## d estimate: -0.04076663 (negligible)
## 95 percent confidence interval:
##        inf        sup
## -0.2114449  0.1299117

## Follow-up analysis: Outcomes

Chi-squared test of 2nd year Retention Outcome

Unmatched

chisq.test(modeldataOmitNA$ASPSCLR, modeldataOmitNA$hcdata...12yrretention.)
##
##  Pearson's Chi-squared test with Yates' continuity correction
##
## data:  modeldataOmitNA$ASPSCLR and modeldataOmitNA$hcdata...12yrretention.
## X-squared = 17.024, df = 1, p-value = 3.692e-05

Matched

matchedRetention <- chisq.test(matchedPairs$ASPSCLR, matchedPairs$hcdata...12yrretention.)
## Warning in chisq.test(matchedPairs$ASPSCLR, matchedPairs$hcdata...
## 12yrretention.): Chi-squared approximation may be incorrect
matchedRetention
##
##  Pearson's Chi-squared test with Yates' continuity correction
##
## data:  matchedPairs$ASPSCLR and matchedPairs$hcdata...12yrretention.
## X-squared = 0.80759, df = 1, p-value = 0.3688
matchedRetention$expected ## ## matchedPairs$ASPSCLR   0     1
##                    0 2.5 263.5
##                    1 2.5 263.5
matchedRetention$observed ## ## matchedPairs$ASPSCLR   0   1
##                    0   4 262
##                    1   1 265

## Optmatch package

Provides a similar set of routines to MatchIt, offering nearest neighbor, 1 to k nearest neigbor, full matching (1 to n), and optimal matching (m to n)

pairmatch: nearest neighbor, 1 to k nearest neighbor, default distance metrics is Mahalanobis distance, but logit propensity scores are available, as are euclidean and user defined distance metrics.

pairs <- pairmatch(as.numeric(ASPSCLR) ~ hcdata.Best_ACT_Comp
+ hcdata.HSGPA + hcdata.incomingcreditstotal
+ gender, data = modeldataOmitNA)

## Pairmatch Output

Output is an optmatch object, with a number of accessible attributes, including the matched.distances and default vector of the groupings.

modeldataOmitNA$pairs <- pairs ordered <- dplyr::filter(modeldataOmitNA[order(pairs),], !is.na(pairs)) head(dplyr::select(ordered, ASPSCLR, pptyScore, pairs), n = 12) ## ASPSCLR pptyScore pairs ## 1 1 0.03654729 1.1 ## 2 0 0.03654729 1.1 ## 3 1 0.03190636 1.10 ## 4 0 0.03190636 1.10 ## 5 0 0.03334549 1.100 ## 6 1 0.03345875 1.100 ## 7 0 0.08494272 1.101 ## 8 1 0.11482639 1.101 ## 9 0 0.04165946 1.102 ## 10 1 0.04165946 1.102 ## 11 1 0.03592469 1.103 ## 12 0 0.03728537 1.103 attr(pairs, "matched.distances")[1:20] ## 1.1 1.10 1.100 1.101 1.102 1.103 ## 0.000000000 0.000000000 0.003575495 0.389313579 0.000000000 0.039330443 ## 1.104 1.105 1.106 1.107 1.108 1.109 ## 0.189835144 0.003575495 0.035754948 0.000000000 0.007150990 0.064358906 ## 1.11 1.110 1.111 1.112 1.113 1.114 ## 0.148375200 0.028603958 0.053632422 0.014301979 0.025028463 0.057207916 ## 1.115 1.116 ## 0.064358906 0.110840338 ## Covariate Balance Checking modeldataOmitNA$pairings <- pairs
propScores <- glm(as.numeric(ASPSCLR) ~ hcdata.Best_ACT_Comp
+ hcdata.HSGPA + hcdata.incomingcreditstotal
+ gender, data = modeldataOmitNA,
family = "binomial")
modeldataOmitNA$pptyscores <- propScores$fitted.values
pairedData <- filter(modeldataOmitNA, !is.na(pairings))

boxplot(pptyscores ~ ASPSCLR , data = modeldataOmitNA,
main = "Prop Scores all data")
boxplot(pptyscores ~ ASPSCLR , data = pairedData,
main = "Prop Scores matched data")

## Covariate Balance

xBalance(as.numeric(ASPSCLR) ~ hcdata.Best_ACT_Comp + hcdata.HSGPA
+ hcdata.incomingcreditstotal + gender, data = modeldataOmitNA)
##                             strata  unstrat
##                             stat   std.diff     z
## vars
## hcdata.Best_ACT_Comp               1.00e+00 1.62e+01 ***
## hcdata.HSGPA                       1.15e+00 1.85e+01 ***
## hcdata.incomingcreditstotal        5.86e-01 9.48e+00 ***
## gender                             4.58e-02 7.42e-01
xBalance(as.numeric(ASPSCLR) ~ hcdata.Best_ACT_Comp + hcdata.HSGPA
+ hcdata.incomingcreditstotal + gender, data = pairedData)
##                             strata  unstrat
##                             stat   std.diff    z
## vars
## hcdata.Best_ACT_Comp                0.00587 0.06780
## hcdata.HSGPA                        0.02355 0.27178
## hcdata.incomingcreditstotal         0.00354 0.04085
## gender                              0.00000 0.00000

## Pair matching with mahalanobis distance and propensity score caliper

modeldataOmitNA$ASPSCLR <- as.integer(modeldataOmitNA$ASPSCLR)

ppty <- glm(ASPSCLR ~ hcdata.Best_ACT_Comp + hcdata.HSGPA +
hcdata.incomingcreditstotal + gender,
data = modeldataOmitNA, family = binomial)

ppty.dist <- match_on(ppty)

mhd.pptyc <- caliper(ppty.dist, width = 1) +
match_on(ASPSCLR ~ hcdata.Best_ACT_Comp + hcdata.HSGPA
+ hcdata.incomingcreditstotal + gender,
data = modeldataOmitNA)

caliperMatch <- pairmatch(mhd.pptyc, data = modeldataOmitNA)
modeldataOmitNA$calipers <- caliperMatch matchedCalipers <- filter(modeldataOmitNA, !is.na(calipers)) ## Balance Checking xBalance(ASPSCLR ~ hcdata.Best_ACT_Comp + hcdata.HSGPA + hcdata.incomingcreditstotal + gender, data = modeldataOmitNA) ## strata unstrat ## stat std.diff z ## vars ## hcdata.Best_ACT_Comp 1.00e+00 1.62e+01 *** ## hcdata.HSGPA 1.15e+00 1.85e+01 *** ## hcdata.incomingcreditstotal 5.86e-01 9.48e+00 *** ## gender 4.58e-02 7.42e-01 xBalance(ASPSCLR ~ hcdata.Best_ACT_Comp + hcdata.HSGPA + hcdata.incomingcreditstotal + gender, data = matchedCalipers) ## strata unstrat ## stat std.diff z ## vars ## hcdata.Best_ACT_Comp -0.01774 -0.20477 ## hcdata.HSGPA 0.02211 0.25518 ## hcdata.incomingcreditstotal 0.00000 0.00000 ## gender -0.00765 -0.08827 ## Full Matching with Optmatch full <- fullmatch(as.numeric(ASPSCLR) ~ hcdata.Best_ACT_Comp + hcdata.HSGPA + hcdata.incomingcreditstotal + gender, data = modeldataOmitNA) modeldataOmitNA$fullMatches <- full
groupings <- data.frame(modeldataOmitNA$fullMatches, modeldataOmitNA$pptyScore, modeldataOmitNA$ASPSCLR) strata <- groupings[order(groupings$modeldataOmitNA.fullMatches),]
colnames(strata) <- c("fullMatches", "pptyScore", "ASPSCLR")
strata[1:15,]
##       fullMatches  pptyScore ASPSCLR
## 9646          1.1 0.03654729       1
## 11986         1.1 0.02602222       0
## 19219         1.1 0.03654729       0
## 4563         1.10 0.03063138       0
## 11032        1.10 0.03190636       1
## 11033        1.10 0.03190636       0
## 11348        1.10 0.03001252       0
## 11469        1.10 0.03115651       0
## 11522        1.10 0.03147577       0
## 11650        1.10 0.03256296       0
## 18790        1.10 0.03063138       0
## 8567        1.100 0.03334549       0
## 11672       1.100 0.03278471       0
## 11673       1.100 0.03278471       0
## 14933       1.100 0.03345875       1
stratumStructure(full)
##    1:1    1:2    1:3    1:4    1:5    1:6    1:7    1:8    1:9   1:10
##     17     16     13     14     13     12     13     17      9     11
##   1:11   1:12   1:13   1:14   1:15   1:16   1:17   1:18   1:19   1:20
##     11      9      6      1      3      4      5      1      3      4
##   1:21   1:22   1:23   1:24   1:26   1:27   1:28   1:29   1:30   1:31
##      4      1      1      3      2      2      1      2      2      2
##   1:32   1:33   1:34   1:35   1:37   1:38   1:39   1:40   1:44   1:45
##      4      1      1      1      1      1      2      1      2      2
##   1:47   1:48   1:50   1:53   1:57   1:60   1:61   1:62   1:63   1:64
##      1      1      1      1      1      1      1      1      1      1
##   1:69   1:70   1:77   1:92   1:93   1:94   1:99  1:102  1:105  1:107
##      1      1      1      1      1      1      1      2      1      1
##  1:109  1:110  1:129  1:142  1:161  1:192  1:214  1:220  1:227  1:266
##      2      1      1      1      1      1      1      1      1      1
##  1:325  1:328  1:359  1:431  1:446  1:494  1:498  1:509  1:640  1:751
##      1      1      1      1      1      1      1      1      1      1
##  1:803  1:889 1:1027 1:1029 1:1967 1:2723
##      1      1      1      1      1      1

## Full Matching: Balance Check

listDists <- lapply(attr(full, "matched.distances"), length)
strata <- attr(listDists, "names")
weights <- 1 / as.integer(listDists)
attr(weights, "names") <- strata
weightFrame <- data.frame(weights, attr(weights, "names"))
merged <- merge(modeldataOmitNA, weightFrame,
by.x = "fullMatches",
by.y = "attr.weights...names..")
merged[merged$ASPSCLR == 1,]$weights <- rep.int(1,nrow(merged[merged$ASPSCLR == 1,])) xBalance(as.integer(ASPSCLR) ~ hcdata.Best_ACT_Comp + hcdata.HSGPA + hcdata.incomingcreditstotal + gender, data = modeldataOmitNA) ## strata unstrat ## stat std.diff z ## vars ## hcdata.Best_ACT_Comp 1.00e+00 1.62e+01 *** ## hcdata.HSGPA 1.15e+00 1.85e+01 *** ## hcdata.incomingcreditstotal 5.86e-01 9.48e+00 *** ## gender 4.58e-02 7.42e-01 xBalance(as.integer(ASPSCLR) ~ hcdata.Best_ACT_Comp + hcdata.HSGPA + hcdata.incomingcreditstotal + gender, data = merged, strata = merged$fullMatches,
stratum.weights = weights)
##                             strata    strat
##                             stat   std.diff    z
## vars
## hcdata.Best_ACT_Comp                0.03156 1.41084
## hcdata.HSGPA                        0.01574 0.81006
## hcdata.incomingcreditstotal         0.01098 0.79063
## gender                              0.00013 0.15361

Weighted t-test:

ASP <- filter(merged, ASPSCLR == 1)
NASP <- filter(merged, ASPSCLR == 0)

t.test(ASP$hcdata.CumGPAendmostrecentUNterm, NASP$hcdata.CumGPAendmostrecentUNterm)
##
##  Welch Two Sample t-test
##
## data:  ASP$hcdata.CumGPAendmostrecentUNterm and NASP$hcdata.CumGPAendmostrecentUNterm
## t = 17.513, df = 286.19, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.3730428 0.4675145
## sample estimates:
## mean of x mean of y
##  3.517315  3.097036

## Weighted t-test:

wtd.t.test(ASP$hcdata.CumGPAendmostrecentUNterm, NASP$hcdata.CumGPAendmostrecentUNterm,
weight = ASP$weights, weighty = NASP$weights)
## Warning in wtd.t.test(ASP$hcdata.CumGPAendmostrecentUNterm, NASP ##$hcdata.CumGPAendmostrecentUNterm, : Treating data for x and y separately
## because they are of different lengths
## $test ## [1] "Two Sample Weighted T-Test (Welch)" ## ##$coefficients
##      t.value           df      p.value
##   2.10289149 465.24811943   0.03601182
##
## 0.08842003 3.51731466 3.42889464 0.04204688