Вы можете вручную выполнить сопоставление, используя 13 строк R-кодов.

Задний план

Проблемой многих обсервационных исследований для получения причинного эффекта является самоотбор — во многих случаях люди выбирают лечение по каким-то причинам, и, следовательно, леченных людей нельзя напрямую сравнивать с нелечеными людьми. Например, учащиеся из богатых семей с большей вероятностью предпочтут учебу в частном колледже. Если мы хотим знать причинно-следственный эффект посещения частного колледжа на заработок, мы должны исключить (или контролировать) влияние семейного происхождения.

Оценка сопоставления широко используется в таких дисциплинах, как статистика, экономика, социология, политология и т. д., для оценки причинно-следственных связей. Это квазиэкспериментальный метод, целью которого является поиск контрфактической единицы, сравнимой с обработанной единицей, среди множества необработанных единиц.

Сопоставление простое, и технически оно во многом связано с обработкой данных. Тем не менее, я видел много материалов, которые затрудняют понимание простого метода. Это расстраивает! Поэтому я решил пройтись по методу и показать вам, что вы можете реализовать метод менее чем с 15 строками R-кода. Обещаю!

R-коды для сопоставления оценок предрасположенности к запуску

Здесь я сосредоточусь на сопоставлении показателей склонности, потому что это популярный метод сопоставления. Другие методы сопоставления аналогичны этому.

Во-первых, давайте смоделируем некоторые данные, используя следующие R-коды. Это несущественные. Просто скопируйте и вставьте. В этих кодах я указываю две ковариаты: W1 и W2. A — бинарное лечение, Y — наблюдаемый результат. Y.0 и Y.1 — это потенциальные исходы, поэтому истинный эффект лечения — это просто разница между Y.0 и Y.1.

# Make the libraries ready
library(dplyr)
library(Matching)
generateData<- function(n) {
 # Generate baseline covariates:
 W1 <- runif(n, min = 0.02, max = 0.7) 
 W2 <- rnorm(n, mean = (0.2 + 0.125*W1), sd = 1) 
 # Generate binary treatment indicator
 A <- rbinom(n, 1, plogis(-.7 + 1.8*W1–1.4*W2))
 # Generate the potential outcomes
 Y.0 <- rnorm(n, mean = (-.5 + 2*poly(W1,2) — W2), sd=1)
 Y.1 <- rnorm(n, mean = (-.5 + 2*poly(W1,2) — W2 + 1.5 + 2*poly(W1,2) — W2), sd=1)
 # Generate the observed outcome
 Y <- ifelse(A == 1, Y.1, Y.0)
 return(data.frame(W1,W2,A,Y,Y.1,Y.0))
}
set.seed(101)
data <- generateData(n = 1000) # Generate 1000 data.

Коды R для сопоставления:

Шаг 1. Оцените логит-модель, чтобы получить оценки склонности:

logitmodel <- glm(A ~ W1 + W2, data = data, family=’binomial’) # L1

Шаг 2. Поиск среди необработанных единиц, которые имеют наименьшие различия в оценочных показателях склонности с обработанными единицами. Затем найденные единицы используются в качестве контрфактических единиц для условного исчисления потенциальных результатов обработанных единиц.

# Add the estimated ps to the data. 
data <- mutate(data, pscore = logitmodel$fitted.values) #L2
data1 <- filter(data, A == 1) # Treated units data #L3
data0 <- filter(data, A == 0) # Untreated units data #L4
match.data <- data1 #Make a copy of treated units data #L5
for(i in 1:nrow(match.data)){#Now find counterfactual units #L6
 temp <- data0 %>% # Search among untreated units #L7
 mutate(pscorei = data1$pscore[i], # PS of the treated unit i #L8
      dif.score = abs(pscorei — pscore)) %>% # Score difference
      arrange(dif.score) %>% # Arrange the data  #L10
      slice(1) %>% # Choose the top one with lowest score dif
      dplyr::select(!!colnames(match.data)) # Keep needed cols# L11
 match.data[i,] <- temp[1,] # Replace with the found unit #L12
}

Шаг 3. Рассчитайте среднюю разницу между наблюдаемым Y и согласованным Y. Эта средняя разница представляет собой средний эффект лечения на леченных (ATT).

mean(data1$Y — match.data$Y) #L13

Конец

Я просто покажу вам, как выполнить сопоставление, используя 13 строк кода R. Существуют альтернативы для ускорения кодов, расчета стандартных ошибок, оценки эффекта лечения на необработанных и т. д. Но это не цели данной статьи. Следите за моим следующим постом!

На всякий случай, если вы мне не верите, вы можете запустить следующие коды из пакета Matching, чтобы получить такое же значение ATT .

m.out <- Match(Y = data$Y, Tr = data$A, X = data$pscore, distance.tolerance = 0)
summary(m.out)
mean(data1$Y.1 - data1$Y.0) # True ATT