R语言cem包进行粗化精确匹配coarsened exact matching (CEM)

粗化精确匹配coarsened exact matching (CEM)是一种进行数据匹配，减少基线间不平衡的办法，也有的叫广义精确匹配。

其实简单的理解就是可以先设定切点，然后把协变量进行分层，把分层的协变量进行匹配，如果我们不设置k2k这个选项，就是随机匹配，如果设置了可以做成近邻匹配或者其他匹配。
下面我们进行演示，先导入数据和R包

library(cem)
Le<-read.csv("E:/r/test/peixun.csv",sep=',',header=TRUE)

这是一个培训数据(公众号回复：培训数据，可以获得数据)，就是观察培训后是否对培训人员收入有改变，treated：是否培训过，age年龄，education：受教育年龄，也就是读了多少年书，black：是不是黑人，married;是否结婚，nodegree：是否有学位，re74：1974年收入指标，re75：1975年收入指标，re78：1978年收入指标，u74：1974年失业率指标，u75：1975年失业率指标，u78：1978年失业率指标，q1：培训前的调查表强烈同意，同意，中立、强烈反对、不同意。
CEM的算法包括三个步骤：

1. 暂时粗化X中的每个控制变量， 为了匹配的目的。例如，受教育年限可以粗略地分为小学、初中、高中、大学和研究生。学校、初中、高中、大学、研究生。
2. 将所有单位分类为阶层，每个阶层都有相同的粗化的X值。从数据集中剔除任何层中不包括至少一个治疗单元和一个控制单元的单元。
   3.从数据集中删去任何不包括至少一个治疗单元和一个控制单元的分层。

我们先删除数据中的缺失值，

Le <- data.frame(na.omit(Le))

先把培训和没有培训的数据提取出来，等会要进行平衡性的检验。

tr <- which(Le$treated==1)
ct <- which(Le$treated==0)
ntr <- length(tr)
nct <- length(ct)

我们比较一下 1978年培训组和没有培训组的平均收入是否有差别，差别还是很明显的

mean(Le$re78[tr]) - mean(Le$re78[ct])

下面我们先定义一下要比较的协变量

vars <- c("age", "education", "black", "married", "nodegree", "re74","re75", "hispanic", "u74", "u75","q1")

观察两组之间协变量的差异

imbalance(group=Le$treated, data=Le[vars])

请注意红箭头标注的L1，若L1=0，则说明两组数据完全平衡，若L1=1，则说明两组数据完全不平衡，越接近1则说明不平衡程度越大。一般来说，匹配后 L1 较匹配前的 L1 有所下降，CEM 的匹配效果较好。CEM 匹配后两组数据的样本量可能会不相等，因此 CEM 匹配过程中会产生权重变量(Weight)，以此来平衡每层中参与组和控制组的人数。
在上表中我们可以看到，变量re74和re75在原始数据中是不平衡的，变量年龄在平均值上是平衡的，但在两个分布的量级上并不平衡。
下面开始进行cem匹配，结果变量re78不能放进去, keep.all=TRUE会返回粗化后的数据集

mat <- cem(treatment = "treated", data = Le, drop = "re78",keep.all=TRUE)
mat

结果显示g0组配对成功95个，g1组配对成功84个，两组个数并不相等。

想看权重的话可以

mat$w

我们把权重加入数据中，

Le$w<-mat[["w"]]

我们把匹配后的数据提取出来

mle<-Le[mat[["matched"]],]

再做一次平衡性，我们可以看到，L1下降了很多

imbalance(group=mle$treated, data=mle, drop=c(todrop,"w"))

刚才的切点是程序自己设定的，下面我们来手动设定一下，分组变量和连续变量设定是不同的，可以看到q1是分成如下组的

levels(Le$q1)

对于有序的分组变量也要按有序排列

q1.grp <- list(c("strongly agree", "agree"), c("neutral","no opinion"), c("strongly disagree","disagree"))

假设我们要对连续变量education（受教育年龄）进行分割，我们先看看变量分布

table(Le$education)

我们可以把它做成小学、中学、高中、大学等进行切割

educut <- c(0, 6.5, 8.5, 12.5, 17)

最后进行匹配

mat1 <- cem(treatment = "treated", data = Le, drop = "re78",
            cutpoints = list(education=educut), grouping=list(q1=q1.grp))

对于教育这个变量，我们看看程序自己设定的切点和我们自己设定的切点

mat$breaks$education
mat1$breaks$education

看一下匹配数，匹配数量进一步增多了

mat1

把数据提取出来，

mle1<-Le[mat1[["matched"]],]
imbalance(group=mle1$treated, data=mle1, drop=todrop)

可以看到L1比之前增大了，不如程序自动分配的好.
接下来我们看一下匹配前和匹配后的效果

allVars <-c("age", "education", "black", "married", "nodegree", "re74","re75", "hispanic", "u74", "u75","q1")###所有变量名
fvars<-c("married", "black","nodegree","q1")#分类变量定义为fvars
tab2 <- CreateTableOne(vars = allVars, strata = "treated" , data = Le, factorVars=fvars,
                       addOverall = TRUE )###绘制基线表
print(tab2,smd = TRUE)#输出表格
tab3 <- CreateTableOne(vars = allVars, strata = "treated" , data = mle, factorVars=fvars,
                        addOverall = TRUE )###绘制基线表
print(tab3,smd = TRUE)#输出表格

Nodegree这个变量改善非常明显，匹配效果还是很好的。