This is a second post in a series on splitting samples. In this case, say you have a very small sub-group of a large sample. You want to look at that subgroup and controls, but you don’t want your sample to be 90% controls. Instead, you want the subgroup and a sub-sample of controls matched on some demographic variables. As a further complication, lets make one variable (age) continuous, and lets make age and sex correlated with subgroup membership. This example is heavily cribbed from a post by Norbert Köhler.

library(sn)
library(tidyverse); library(ggthemes); theme_set(theme_tufte())
library(ggExtra)
library(pander)
library(MatchIt)
library(simstudy)


# Simulation

set.seed(31453)

simdef<-defData(varname="age",
dist="uniformInt",
formula="120;144") #age in months between 10-12
simdef<-defData(simdef,varname="sex",
dist="binary",
formula=".5")
simdef<-defData(simdef,varname="parent.ed",
dist="categorical",
formula=genCatFormula(n=6))
simdef<-defData(simdef,varname="missingdata",
dist="binary",
formula=".2")
simdef<-defData(simdef,varname="inSubGroup",
dist="binary",
formula=".005/12 * (age-132) + .005*sex + .0175")

df<-genData(12000,simdef)

df$income<-rsn(nrow(df),alpha=3) numbers_of_bins = 10 df <- df %>% mutate( # bin i: i.bin = cut(income, breaks = unique(quantile( income, probs = seq.int(0, 1, by = 1 / numbers_of_bins) )), include.lowest = TRUE, labels=FALSE) ) df<-as.data.frame(df) factorialize<-c("sex","missingdata","parent.ed","inSubGroup") df[factorialize] <- lapply(df[factorialize], factor) levels(df$inSubGroup)<-c("control","treatment")


id age sex parent.ed missingdata inSubGroup income i.bin
1 128 0 4 1 control 1.392 9
2 138 1 5 0 control -0.02661 1
3 130 1 6 0 control 1.911 10
4 127 1 5 0 control 0.3791 4
5 121 0 5 0 control 0.909 7
6 132 1 1 0 control 1.695 10
pander(table(df$inSubGroup))  control treatment 11761 239 # Is there an imbalance? We only need to match on age, sex, and one other categorical variable. imbalance_model <- matchit( inSubGroup ~ age +sex + i.bin, data = df, method = NULL, distance = "glm" ) summary(imbalance_model)  ## ## Call: ## matchit(formula = inSubGroup ~ age + sex + i.bin, data = df, ## method = NULL, distance = "glm") ## ## Summary of Balance for All Data: ## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean ## distance 0.0202 0.0199 0.1269 0.9939 0.0346 ## age 132.4812 131.8617 0.0850 1.0289 0.0266 ## sex0 0.4686 0.5021 -0.0671 . 0.0335 ## sex1 0.5314 0.4979 0.0671 . 0.0335 ## i.bin 5.3096 5.5039 -0.0672 1.0122 0.0195 ## eCDF Max ## distance 0.0856 ## age 0.0600 ## sex0 0.0335 ## sex1 0.0335 ## i.bin 0.0367 ## ## ## Sample Sizes: ## Control Treated ## All 11761 239 ## Matched 11761 239 ## Unmatched 0 0 ## Discarded 0 0  Yes, age and sex are imbalanced (which we simulated). So is income! # Nearest Neighbor Matching Note that it is important to code variable type correctly, i.e. that factors are factors and not numeric. ## Sub-sampling We want 2 controls per ‘treatment’ participant. matching_model <- matchit( inSubGroup ~ age + sex + i.bin, data = df, method = "nearest", distance = "glm", ratio= 2 ) summary(matching_model,un=FALSE)  ## ## Call: ## matchit(formula = inSubGroup ~ age + sex + i.bin, data = df, ## method = "nearest", distance = "glm", ratio = 2) ## ## Summary of Balance for Matched Data: ## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean ## distance 0.0202 0.0202 0 1.0021 0 ## age 132.4812 132.4812 0 1.0021 0 ## sex0 0.4686 0.4686 0 . 0 ## sex1 0.5314 0.5314 0 . 0 ## i.bin 5.3096 5.3096 0 1.0021 0 ## eCDF Max Std. Pair Dist. ## distance 0 0 ## age 0 0 ## sex0 0 0 ## sex1 0 0 ## i.bin 0 0 ## ## Sample Sizes: ## Control Treated ## All 11761 239 ## Matched 478 239 ## Unmatched 11283 0 ## Discarded 0 0  The distribution parameters should be similar, and the control n should be twice the treated n. Then we save the new data frame: df.match<-match.data(matching_model)[,1:ncol(df)]  ## Checking df.match$inSubGroup<-factor(df.match\$inSubGroup,labels=c("Control","Treatment"))
histbygroup<-function(plotdata,xvar="i") {
ggplot(plotdata,aes_string(x=xvar)) +
geom_density() +
facet_grid(~inSubGroup)
}

histbygroup(df.match,"i.bin")


histbygroup(df.match,"age")


barbygroup<-function(plotdata,xvar="i") {
ggplot(plotdata,aes_string(x=xvar)) +
geom_bar() +
facet_grid(~inSubGroup)
}

barbygroup(df.match,"sex")


Nearest neighbor is doing pretty well! Here’s another method for comparison:

# Optimal Matching

## Sub-sampling

We want 2 controls per ‘treatment’ participant.

matching_model2 <-
matchit(
inSubGroup ~ age + sex + i.bin,
data = df,
method = "optimal",
distance = "glm",
ratio= 2
)

summary(matching_model2,un=FALSE)

##
## Call:
## matchit(formula = inSubGroup ~ age + sex + i.bin, data = df,
##     method = "optimal", distance = "glm", ratio = 2)
##
## Summary of Balance for Matched Data:
##          Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance        0.0202        0.0201          0.0546     1.0202    0.0177
## age           132.4812      132.3891          0.0126     1.0956    0.0176
## sex0            0.4686        0.5209         -0.1048          .    0.0523
## sex1            0.5314        0.4791          0.1048          .    0.0523
## i.bin           5.3096        5.2573          0.0181     1.0224    0.0073
##          eCDF Max Std. Pair Dist.
## distance   0.0607          0.2131
## age        0.0460          0.7938
## sex0       0.0523          0.7085
## sex1       0.0523          0.7085
## i.bin      0.0167          0.9767
##
## Sample Sizes:
##           Control Treated
## All         11761     239
## Matched       478     239
## Unmatched   11283       0


The distribution parameters should be similar, and the control n should be twice the treated n. Then we save the new data frame:

df.match2<-match.data(matching_model2)[,1:ncol(df)]


## Checking

histbygroup<-function(plotdata,xvar="i") {
ggplot(plotdata,aes_string(x=xvar)) +
geom_density() +
facet_grid(~inSubGroup)
}

histbygroup(df.match2,"i.bin")


histbygroup(df.match2,"age")


barbygroup<-function(plotdata,xvar="i") {
ggplot(plotdata,aes_string(x=xvar)) +
geom_bar() +
facet_grid(~inSubGroup)
}

barbygroup(df.match2,"sex")


Nearest neighbor appears to do a better job