1 Эксперименты

1.1 Короткое напоминание теории

1.1.1 Потенциальные исходы

Обозначения:

\(X_i\) – независимые переменные (covariates)
\(T_i\) – бинарная переменная воздействия (treatment variable): \[\begin{equation*} T_i = \begin{cases} 1, &\text{воздействие на объект i оказано}\\ 0, &\text{воздействие на объект i не оказано} \end{cases} \end{equation*}\]
\(Y_{i1}\), \(Y_{i0}\) – потенциальные исходы (potential outcomes)

Наблюдаемые исходы \(Y_i\) отличаются от потенциальных исходов. Потенциальные исходы являются гипотетическими случайными величинами, когда наблюдаемые исходы являются фактическими случайными величинами. Наблюдаемые исходы являются функцией от потенциальных исходов:

\(Y_i = T_i \cdot Y_{i1} + (1-T_i) \cdot Y_{i0}\)

Действительно при \(T=1\) мы получим \(Y_i = Y_{i1}\), а при \(T=0\) получим \(Y_i = Y_{i0}\).

Тогда наш эффект воздействия для конкретного наблюдения равен разнице между двумя состояниями мира для этого наблюдения (потенциальными исходами):

\(\tau_i = Y_{i1} - Y_{i0}\)

Фундаментальная проблема причинного вывода (Fundamental problem of causal inference):

Чтобы оценить эффект воздействия для конкретного индивида, мы должны знать потенциальные исходы сразу для двух его состояний мира, но реально мы наблюдаем только одно из них – либо \(Y_{i1}\), если индивид подвергся воздействию, либо \(Y_{i0}\), если он ему не подвергался. Оценка индивидуального эффекта требует доступа к данным, которых у нас физически не может быть.

1.1.2 Средние эффекты

Если с распределением индивидуального эффекта воздействия (treatment effect) работать не получается, будем довольствоваться средними величинами. Например, попробуем рассчитать средний эффект воздействия (average treatment effect):

\(ATE = \mathbb{E}[\tau_i] = \mathbb{E}[Y_{i1} - Y_{i0}] = \mathbb{E}[Y_{i1}] - \mathbb{E}[Y_{i0}] = \frac{1}{N_1}\sum \limits_{i=1}^{N_1}Y_{i1} - \frac{1}{N_0}\sum \limits_{i=1}^{N_0}Y_{i0}\)

Следующий шаг, который мы предпримем, попробуем рассчиать величину эффекта только для тритмент группы – средний эффект воздействия на задействованных (average treatment effect for the treatment group):

\(ATT = \mathbb{E}[\tau_i|T_i=1] = \mathbb{E}[Y_{i1} - Y_{i0}|T_i=1] = \mathbb{E}[Y_{i1}|T_i=1] - \mathbb{E}[Y_{i0}|T_i=1]\)

А теперь то же самое, но для контрольной группы – средний эффект воздействия на незадействованных (average treatment on the non-treated):

\(ATnT = \mathbb{E}[\tau_i|T_i=0] = \mathbb{E}[Y_{i1} - Y_{i0}|T_i=0] = \mathbb{E}[Y_{i1}|T_i=0] - \mathbb{E}[Y_{i0}|T_i=0]\)

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

Самая простая идея для оценки ATE, которая всем придет в голову, взять простую разницу в средних: \(\mathbb{E}[Y_1|T=1] - \mathbb{E}[Y_0|T=0]\).

Но тут всё не так просто, после небольших преобразований, с которыми можно ознакомиться в части 4.1.3 учебника мы получим следующее:

\(\mathbb{E}[Y_1|T=1] - \mathbb{E}[Y_0|T=0] = \underbrace{\mathbb{E}[Y_1] - \mathbb{E}[Y_0]}_{\text{ATE}} + \underbrace{\mathbb{E}[Y_0|T=1] - \mathbb{E}[Y_0|T=0]}_{\text{Selection Bias}} + \underbrace{(1-\pi)(ATT - ATnT)}_{\text{Heterogeneous treatment effect bias}}\)

ATE – интересующий нас эффект
Selection Bias – смещение, возникающее из-за того, что контрольная группа и группа воздействия различались, даже если бы на них не было оказано воздействие, то есть имеет место некоторый дисбаланс
Heterogeneous treatment effect bias – различие в интенсивности эффекта для тритмент и контрольной группы, взвешенное на долю выборки \((1-\pi)\), которая попала в контрольную группу

Далее мы рассмотрим предпосылки, которые позволяют нивелировать влияние этих двух смещений и получить оценку ATE.

1.1.3 Предпосылки

Экзогенность воздействия (Independence assumption)

Экзогенность воздействия (Independence assumption) означает, что распределение объекта в тритмент или контрольную группы осуществляется случайно и независимо от его изначальных характеристик. Данная предпосылка обычно обозначается следующим образом \((Y_1, Y_0, X)_i \perp T_i\)

Технически для нас это значит следующее:

\(\mathbb{E}[Y_0|T=1] - \mathbb{E}[Y_0|T=0] = 0 \Rightarrow \text{Selection Bias}=0\)
\(\mathbb{E}[Y_1|T=1] - \mathbb{E}[Y_1|T=0] = 0\)
- \((1-\pi)(ATT - ATnT) = (1-\pi)\left[(\mathbb{E}[Y_1|T=1]-\mathbb{E}[Y_0|T=1])-(\mathbb{E}[Y_1|T=0]-\mathbb{E}[Y_0|T=0])\right] = 0 \Rightarrow \text{Heterogeneous treatment effect bias}=0\)

То есть хорошая рандомизация, а следовательно, и выполнение предпосылок, позволяет нам очистить эффект воздействия от двух типов смещения, в этом случае:

\(ATE = \mathbb{E}[Y_1] - \mathbb{E}[Y_0] = \mathbb{E}[Y_1|T=1] - \mathbb{E}[Y_0|T=0] \xrightarrow{p} \frac{1}{N_1}\sum \limits_{i=1}^{N_1}Y_{i1} - \frac{1}{N_0}\sum \limits_{i=1}^{N_0}Y_{i0}\)

Отсутствие “внешних эффектов” воздействия (SUTVA – Stable unit treatment value assumption)

Эта предпосылка подразумевает выполнение двух вещей. Во-первых, воздействие оказывается только на один объект и внешние эффекты у него отсутствуют. Во-вторых, воздействие гомогенно – существует только один тип тритмента.

1.2 Работа со случайными числами в R

В рамках курса нам неоднократно придется прибегнуть к генерации случайных чисел из какого-нибудь закона распределения, как правило, из нормального. Для работы со случайными числами из нормального распределения в R используется несколько функций: dnorm() для плотности вероятности, pnorm() для функции распределения, qnorm() для квантилей распределения и rnorm() для генерации случайных чисел. Почитать подробнее про них и посмотреть примеры можно тут или тут.

Если вы не уверены в синтаксисе какой-то функции, вы можете воспользоваться встроенной справкой в RStudio (справочная информация появится во вкладке Help):

?rnorm

Если вам когда-нибудь понадобятся функции для других распределений, вы можете, конечно, их просто загуглить или ввести в окне Help запрос “Distributions”:

?Distributions

Теперь, когда мы точно знаем, как устроены аргументы у функции rnorm(n;mean;sd), давайте попробуем достать 5 чисел из стандартного нормального распределения \(\mathcal{N} \sim \left(0;1\right)\):

rnorm(5,0,1)

[1]  0.8028321  0.1841247 -0.1678564  0.6213981 -1.4385412

Если вы повторяли за мной, то у вас, вероятно, получились другие значения. А теперь я тоже попробую ещё раз:

rnorm(5,0,1)

[1] -1.91234816 -0.05973835 -1.86892906 -0.40398381  1.19708204

Как и ожидалось, у меня тоже не совпало с предыдущим результатом. Существует по меньшей мере две причины, почему нам хотелось бы получать одинаковые результаты:

Работа на семинарах, когда нам с вами было бы удобно сверяться
Воспроизводимые научные исследования, когда любой читатель может реплицировать и проверить на корректность ваш результат

В R существует несколько алгоритмов, позволящих генерировать случайные числа (Random Number Generator (RNG), подробнее – ?RNGkind). По умолчанию используется Mersenne twister (Вихрь Мерсенна), поскольку он работает наиболее быстро. Однако все эти алгоритмы на самом деле детерминированы, поэтому генерируемые ими числа корректнее называть “псевдослучайными”.

Генератор псевдослучайных чисел начинает свою работу с определенной точки в пространстве возможных чисел. Эту точку приянто называть начальное число или seed. Начальное число – это число (или вектор), используемое для инициализации генератора псевдослучайных чисел. Если вы хотите получить одинаковые результаты с помощью генератора случайных чисел, важно установить начальное число. Для этого в R используется функция set.seed(). По умолчанию начальное число не установлено, если ничего не указать, то создается новое из текущего времени и идентификатора процесса на вашем устройстве. Следовательно, по умолчанию разные сеансы дают разные результаты моделирования.

Давайте снова попробуем сгенерировать 5 чисел из случайного нормального распределения, но в этот раз зафиксируем seed:

set.seed(123)
rnorm(5,0,1)

[1] -0.56047565 -0.23017749  1.55870831  0.07050839  0.12928774

И еще разок:

set.seed(123)
rnorm(5,0,1)

[1] -0.56047565 -0.23017749  1.55870831  0.07050839  0.12928774

Получилось! В качестве seed можно использовать любое число, главное, чтобы при репликации оно каждый раз было одинаковым.

Почитать про set.seed можно тут на русском и тут на английском.

1.3 Мини-симуляция

Симуляции – это “игрушечные” примеры, которые позволят нам на протяжении курса разбираться с разными методами. Игрушечные они потому, что за ними не стоят реальные данные. Данные для симуляций мы будем специальным образом заранее моделировать. Это бывает удобно, когда идеально подходящих данных нет, или они не лежат в открытом доступе. К тому же, живые данные часто могут быть зашумлены из-за других факторов, на которые нам не всегда будет удобно отвлекаться на занятиях. Но с живыми данными мы тоже обязательно будем работать! В том числе, они будут доступны вам в домашних заданиях :)

Теперь смоделируем небольшую симуляцию по мотивам изученного материала.

1.3.1 Подготовка данных

Давайте смоделируем гипотетическую ситуацию. Мы хотим оценить величину эффекта от использование сайта с расписанием cacs.ws на свободное время студента.

Предположим, что наша экспериментальная выборка состоит из 1000 человек:

N <- 1000

Для простоты у них будет всего две характеристики (\(X\) – возраст и \(Z\) – время в пути от дома до ЭФ), имеющих влияние на потенциальный исход (\(Y_0\) и \(Y_1\) – свободное время). При этом предположим, что реальный эффект воздействия равен \(\tau=15\) минутам. Будем считать, что реальная зависимость потенциального исхода от ковариатов и тритмента выглядит следующим образом:

\(Y_0 = 240 - 3 \cdot X - Z + 15 \cdot \underbrace{T}_{= 0} + \varepsilon = 120 - 3 \cdot X - Z + \varepsilon\) – потенциальный исход, если не было оказано воздействие
\(Y_1 = 240 - 3 \cdot X - Z + 15 \cdot \underbrace{T}_{= 1} + \varepsilon = 120 - 3 \cdot X - Z + 15 + \varepsilon\) – потенциальный исход, если было оказано воздействие

Сгенерируем \(X\) и \(Z\) случайным образом:

set.seed(123)
X <- runif(N, 18, 25)
Z <- rnorm(N, 60, 20)
# Посмотрим на наши данные, нарисуем сразу 2 графика -- диаграмму разброса и гистограмму
par(mfrow=c(2,2)) # читаем буквально как матрицы -- 2 строки и 2 столбца
plot(X)
hist(X)
plot(Z)
hist(Z)

Также создадим случайные ошибки для каждого из типов индивидов:

set.seed(123)
e <- rnorm(N, 0, 1)
plot(e, col="red")

Вспомогательно пронумеруем наши наблюдения:

number <- 1:N

И соберем всё в общий датасет:

data <- data.frame(number = number, X=X, Z=Z) # сшиваем столбики в data frame
head(data, 10) # смотрим только первые несколько строк датасета

   number        X        Z
1       1 20.01304 47.96214
2       2 23.51814 40.12603
3       3 20.86284 80.53570
4       4 24.18112 75.02123
5       5 24.58327 29.81667
6       6 18.31890 58.09705
7       7 21.69674 42.08104
8       8 24.24693 18.58498
9       9 21.86005 63.00240
10     10 21.19630 58.41577

1.3.2 Рандомизация

Что касается рандомизации, то наша глобальная задача, используя разные функции придумать способ, как случайным образом разбить выборку на 2 группы. Существует огромное множество вариантов, как это можно сделать. Мы обсудим лишь несколько из них. На практике вы можете пользоваться как этими вариантами, так и придумать что-то своё.

1.3.2.1 Использование порядкового номера

Первый вариант – просто взять половину наблюдений по их порядку расположения в датасете. Пусть воздействие будет оказано только на первую половину нашей выборки.

T1 <- c(rep(1,N/2), rep(0,N/2))
mean(T1) # проверяем долю тритмента, должна получиться ровно половина

[1] 0.5

data$T1 <- T1 # добавляем переменную в датасет

Тут мы воспользовались функцией c(...), которая позволяет объединить величины внутри неё в вектор, и функцией rep(x;times) которая повторяет величину \(x\) столько раз, сколько указано в величине \(times\).

Второй вариант, который практически дублирует первый, но немного отличается – мы можем взять каждое второе наблюдение:

T2 <- rep(0:1, N/2)
mean(T2) # проверяем долю тритмента, должна получиться ровно половина

[1] 0.5

data$T2 <- T2 # добавляем переменную в датасет

1.3.2.2 Использование свойств распределений

Логично, что если мы хотим случайным образом присвоить наблюдениям разные группы, то мы можем начать с того, что присвоим им случайные числа, которые мы потом разными способами переведем в бинарный формат, который уже будет соотвествовать конкретной группе. Протестируем этот способ на примере трех распределений – нормального, равномерного и биномиального.

1.3.2.2.1 Равномерное распределение

set.seed(123)
V1 <- runif(N) # генерим вспомогательную переменную из равномерного распределения
T3 <- as.numeric(V1 < median(V1)) #генерим тритмент, отсекая половину выборки по медиане; as.numeric используется, чтобы перейти от логического типа (true, false) к численному (1 и 0)
mean(T3) # проверяем долю тритмента, должна получиться ровно половина

[1] 0.5

data$T3 <- T3 # добавляем переменную в датасет

1.3.2.2.2 Нормальное распределение

set.seed(123)
V2 <- rnorm(N) # генерим вспомогательную переменную из нормального распределения
T4 <- as.numeric(V2>0) #генерим тритмент, отсекая половину выборки по знаку вспомогательной переменной; as.numeric используется, чтобы перейти от логического типа (true, false) к численному (1 и 0)
summary(T4) # проверяем долю тритмента, должна получиться примерно половина

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   1.000   0.505   1.000   1.000

data$T4 <- T4 # добавляем переменную в датасет

1.3.2.2.3 Биномиальное распределение

set.seed(123)
T5 <- rbinom(N, 1, 0.5) #генерим тритмент, отсекая половину выборки по знаку вспомогательной переменной; as.numeric используется, чтобы перейти от логического типа (true, false) к численному (1 и 0)
summary(T5) # проверяем долю тритмента, должна получиться примерно половина

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   0.000   0.493   1.000   1.000

data$T5 <- T5 # добавляем переменную в датасет

1.3.2.3 Использование готовых пакетов

1.3.2.3.1 Пакет `experiment`

library('experiment') 
set.seed(123)
rand <- randomize(data, group= c("Treat", "Control"))
T6 <- as.numeric(rand$treatment == "Treat") # преобразовываем переменную в численный формат 
data$T6 <- T6 # добавляем переменную в датасет
head(data, 10)

   number        X        Z T1 T2 T3 T4 T5 T6
1       1 20.01304 47.96214  1  0  1  0  0  1
2       2 23.51814 40.12603  1  1  0  0  1  1
3       3 20.86284 80.53570  1  0  1  1  0  1
4       4 24.18112 75.02123  1  1  0  1  1  0
5       5 24.58327 29.81667  1  0  0  1  1  1
6       6 18.31890 58.09705  1  1  1  1  0  0
7       7 21.69674 42.08104  1  0  0  1  1  0
8       8 24.24693 18.58498  1  1  0  0  1  1
9       9 21.86005 63.00240  1  0  0  0  1  1
10     10 21.19630 58.41577  1  1  1  0  0  1

1.3.2.3.2 Пакет `radomizr`

library('randomizr') 
set.seed(123)
T7 <- complete_ra(N=N, m=N/2)
data$T7 <- T7 # добавляем переменную в датасет
head(data, 10)

   number        X        Z T1 T2 T3 T4 T5 T6 T7
1       1 20.01304 47.96214  1  0  1  0  0  1  0
2       2 23.51814 40.12603  1  1  0  0  1  1  0
3       3 20.86284 80.53570  1  0  1  1  0  1  0
4       4 24.18112 75.02123  1  1  0  1  1  0  1
5       5 24.58327 29.81667  1  0  0  1  1  1  0
6       6 18.31890 58.09705  1  1  1  1  0  0  1
7       7 21.69674 42.08104  1  0  0  1  1  0  1
8       8 24.24693 18.58498  1  1  0  0  1  1  0
9       9 21.86005 63.00240  1  0  0  0  1  1  0
10     10 21.19630 58.41577  1  1  1  0  0  1  0

1.3.2.4 Хэш-функции

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

Хэш-функция (hash function) – функция, осуществляющая преобразование массива входных данных произвольной длины в выходную битовую строку установленной длины, выполняемое определённым алгоритмом. Преобразование, производимое хэш-функцией, называется хэшированием.

Хэш-функции применяются в следующих случаях:

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

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

Пример хэширования с помощью пакета digest:

library("digest")
hash <- digest('econometrics', algo="murmur32")
hash

[1] "24b6de5a"

Пакет digest применяет хеш-функцию к произвольным объектам R. В пакете реализовано много разных алгоритмов преобразований, однако мы будем использовать алгоритм “murmur32” (подробнее почитать можно тут). Этот алгоритм совершенно не подходит для криптографических целей, но зато идеально подходит нам, поскольку он 32-битный, что позволяет нам перевести полученный хэш-код в целое число.

Прежде чем начать разбирать пример, немного отвлечемся на техническую полезную вещь – функцию sapply().

Про семейство apply функций рекомендую почитать подробнее на русском тут или тут, а на английском тут.

Прелесть функции sapply() состоит в том, что она позволяет нам избежать громоздких циклов при написании кода и ускорить вычисления благодаря “векторной ориентированности” языка R. Например, функция позволяет нам найти минимальное и максимальное значение для каждой из ковариат:

sapply(list(X,Z), min)

[1] 18.003257  3.804506

sapply(list(X,Z), max)

[1]  24.99583 127.80742

То есть функция проводит однотипную операцию (по сути вложенную функцию) над каждым элементом списка. Если кто-то хорошо владеет python, то аналогичной функцией там является map().

Вернемся к нашей симуляции. Изначально мы это все затеяли, чтобы сделать хорошую рандомизацию. Сначала хэшируем номера наших наблюдений. Так делать не очень хорошо, обычно, для хэширования данных используют ФИО и/или СНИЛС, но мы не будем ради этого громоздить генерацию еще одной переменной. Для иллюстрации просто воспользуемся номерами, это тоже сработает.

hashes <- sapply(data$number, function(x) {digest(x, algo="murmur32")}) # хэшируем строчки
hashes[1:20]

 [1] "5e6216f3" "8b8b3789" "ad1bf356" "b2df92b2" "1d1b8e3a" "98d09cab"
 [7] "2e366c5a" "4d110361" "91d6bfe0" "f568e543" "08bd74cd" "e979af4b"
[13] "a82e5827" "485cbd08" "bc4e5bd8" "95f68430" "7180d744" "da414596"
[19] "70aa2885" "afd0ffd2"

Далее нужно перевести получившиеся строчки в цифровой числовой. Функция strtoi() конвертирует строковое представление числа, которое хранится в строке, в длинное целое.

result <- strtoi(substring(hashes, 2), base=16) 
result[1:20]

 [1] 241309427 193673097 219935574  48206514 219909690 147889323 238447706
 [8] 219218785  30851040  90760515 146633933 158969675 137254951 140295432
[15] 206461912 100041776  25220932 172049814  11151493 265355218

Далее, используя эти числа, нужно как-то разбить выборку на тритмент и контроль. Будем смотреть на последнюю цифру, если она меньше 5, то наблюдение попадет в тритмент, если больше – в контроль.

T8 <- result %% 10 < 5 # берем остаток от деления на 10 (получится последняя цифра) и сравниваем его с 5
data$T8 <- as.numeric(T8) # переводим из логического типа в числовой и добавляем в датасет
data$T8[1:20]

 [1] 0 0 1 1 1 1 0 0 1 0 1 0 1 1 1 0 1 1 1 0

Проверим нашу рандомизацию на сбалансированность:

summary(data)

     number             X               Z                 T1            T2     
 Min.   :   1.0   Min.   :18.00   Min.   :  3.805   Min.   :0.0   Min.   :0.0  
 1st Qu.: 250.8   1st Qu.:19.78   1st Qu.: 46.232   1st Qu.:0.0   1st Qu.:0.0  
 Median : 500.5   Median :21.43   Median : 60.576   Median :0.5   Median :0.5  
 Mean   : 500.5   Mean   :21.48   Mean   : 60.239   Mean   :0.5   Mean   :0.5  
 3rd Qu.: 750.2   3rd Qu.:23.23   3rd Qu.: 73.101   3rd Qu.:1.0   3rd Qu.:1.0  
 Max.   :1000.0   Max.   :25.00   Max.   :127.807   Max.   :1.0   Max.   :1.0  
       T3            T4              T5              T6            T7     
 Min.   :0.0   Min.   :0.000   Min.   :0.000   Min.   :0.0   Min.   :0.0  
 1st Qu.:0.0   1st Qu.:0.000   1st Qu.:0.000   1st Qu.:0.0   1st Qu.:0.0  
 Median :0.5   Median :1.000   Median :0.000   Median :0.5   Median :0.5  
 Mean   :0.5   Mean   :0.505   Mean   :0.493   Mean   :0.5   Mean   :0.5  
 3rd Qu.:1.0   3rd Qu.:1.000   3rd Qu.:1.000   3rd Qu.:1.0   3rd Qu.:1.0  
 Max.   :1.0   Max.   :1.000   Max.   :1.000   Max.   :1.0   Max.   :1.0  
       T8       
 Min.   :0.000  
 1st Qu.:0.000  
 Median :1.000  
 Mean   :0.508  
 3rd Qu.:1.000  
 Max.   :1.000

1.3.3 Реализация исходов

Посчитаем потенциальные исходы:

Y0 <- 240 - 3*X - Z + 15*0 + e # T=0
Y1 <- 240 - 3*X - Z + 15*1 + e # T=1
data$Y0 <- Y0 # добавляем переменные в датасет
data$Y1 <- Y1

А также наблюдаемые исходы:

data$Y <- T1*Y1 + (1-T1)*Y0

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

Функция для генерации данных

generate_data <- function(N){
  N=N
  X <- runif(N, 18, 25)
  Z <- rnorm(N, 60, 20)
  e0 <- rnorm(N/2, 0, 1)
  e1 <- rnorm(N/2, 0, 1)
  number <- 1:N
  Y0 <- 240 - 3*X - Z + 15*0 + e0 
  Y1 <- 240 - 3*X - Z + 15*1 + e1 
  data <- data.frame(number = number, X=X, Z=Z, Y0=Y0, Y1=Y1)
}

Функция для рандомизации

randomization <- function(data, type){
  N <- nrow(data)
  if (type=='in order'){
    T <- c(rep(1,N/2), rep(0,N/2))
    data$T <- T
  } else if (type=='by turns'){
    T <- rep(0:1, N/2)
    data$T <- T
  } else if (type=='unif'){
    V <- runif(N)
    T <- as.numeric(V1 < median(V1))
    data$T <- T
  } else if (type=='norm'){
    V <- rnorm(N) 
    T <- as.numeric(V2>0)
    data$T <- T
  } else if (type=='binom'){
    T <- rbinom(N, 1, 0.5)
  }
  return(data)
}

Функция для расчета наблюдаемого исхода

outcome <- function(data){
  observed_data <- data
  observed_data$Y <- data$T*observed_data$Y1 + (1 - data$T)*observed_data$Y0
  return(observed_data)
}

А теперь посмотрим, сколько строчек займет наша симуляция с использованием функций:

set.seed(123)
data <- generate_data(1000)
set.seed(123)
data <- randomization(data=data, type='in order')
data <- outcome(data=data)
tail(data,10)

     number        X         Z        Y0        Y1 T         Y
991     991 18.61054  75.36015 107.71110 125.93429 0 107.71110
992     992 19.56646  37.83344 144.39269 158.58181 0 144.39269
993     993 22.00504  44.27528 129.95640 144.11565 0 129.95640
994     994 20.80118 105.68233  71.17735  87.99218 0  71.17735
995     995 21.95826  38.13398 134.71123 149.89166 0 134.71123
996     996 23.80737  64.28959 104.36495 120.01488 0 104.36495
997     997 22.49480  77.85142  94.91935 111.10506 0  94.91935
998     998 20.74049  80.37516  97.68081 112.19320 0  97.68081
999     999 22.96706  81.78224  89.85344 105.76786 0  89.85344
1000   1000 18.76177  56.73742 126.51679 142.61883 0 126.51679

1.3.4 Считаем эффект

Когда данные подготовлены, мы забываем, что что-то о них знали :)

С данного момента мы действуем согласно предпосылке, что мы знаем только \(Y\), \(X\) и \(T\).

Поскольку мы знаем, что тритмент распределялся среди наблюдений случайно, мы можем сделать вывод, что предпосылка о независимости воздействия выполнена. Также мы знаем, что данные устроены так, что SUTVA тоже выполнена. Следовательно, можем использовать обычную разницу в средних, чтобы оценить эффнект воздействия:

\(\widehat{ATE} = \overline{Y_1} - \overline{Y_0}\)

ATE_hat <- mean(data[which(data$T == 1),]$Y) - mean(data[which(data$T == 0),]$Y)
ATE_hat

[1] 15.53169

Из курса ЭКМ-2 помним, что то же самое можно было бы получить, с помощью обычного парного МНК:

model1 <- lm(Y ~ T, data=data)
summary(model1)


Call:
lm(formula = Y ~ T, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-78.525 -14.782  -0.268  14.040  64.919 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 115.0500     0.9369  122.80   <2e-16 ***
T            15.5317     1.3249   11.72   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 20.95 on 998 degrees of freedom
Multiple R-squared:  0.121, Adjusted R-squared:  0.1202 
F-statistic: 137.4 on 1 and 998 DF,  p-value: < 2.2e-16

Если обе предпосылки выполнены, то оценка должна быть несмещенной даже без контрольных переменных (ковариат). Попробуем их добавить и сравним оценки эффекта:

model2 <- lm(Y ~ T + X + Z, data=data)
summary(model2)


Call:
lm(formula = Y ~ T + X + Z, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.97985 -0.63057 -0.01498  0.69005  3.04683 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 240.123033   0.347927   690.2   <2e-16 ***
T            14.874547   0.061883   240.4   <2e-16 ***
X            -2.988124   0.015384  -194.2   <2e-16 ***
Z            -1.005280   0.001546  -650.4   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9783 on 996 degrees of freedom
Multiple R-squared:  0.9981,    Adjusted R-squared:  0.9981 
F-statistic: 1.732e+05 on 3 and 996 DF,  p-value: < 2.2e-16

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

Далее в качестве упраженения посчитаем средние эффекты на разных группах:

ATT <- mean(data[which(data$T == 1),]$Y1) - mean(data[which(data$T == 1),]$Y0)
ATT

[1] 14.87723

ATnT <- mean(data[which(data$T == 0),]$Y1) - mean(data[which(data$T == 0),]$Y0)
ATnT

[1] 14.87723

Получили, то что и должны были получить \(ATT = ATnT\). Если мы вернемся в раздел, где мы подготовливали данные, то увидим, что в двух потенциальных исходах мы “зашифровали” одинаковый эффект:

\(Y_0 = 120 - 3 \cdot X - Z + 15 \cdot \underbrace{T}_{= 0} + \varepsilon = 120 - 3 \cdot X - Z + \varepsilon\)
\(Y_1 = 120 - 3 \cdot X - Z + 15 \cdot \underbrace{T}_{= 1} + \varepsilon = 120 - 3 \cdot X - Z + 15 + \varepsilon\)

Содержательно это иллюстрирует то, что эффект воздействия гомогенен, то есть \(\text{Heterogeneous treatment effect bias} = (1-\pi)(ATT - ATnT) = 0\).

1.3.5 Баланс ковариатов

Теперь проверим насколько наша рандомизация оказалась хорошей, то есть проверим контрольную и тритмент группу на “похожесть”. Сравним средние значения ковариат в двух группах.

1.3.5.1 Сравнение средних

mean(data[which(data$T == 1),]$X) - mean(data[which(data$T == 0),]$X)

[1] -0.02791759

mean(data[which(data$T == 1),]$Z) - mean(data[which(data$T == 0),]$Z)

[1] -0.5707093

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

1.3.5.2 Тест Стьюдента

А теперь проведем тест на разницу в средних с помощью формального t-теста Стьюдента:

\(H_0:\) среднее значение параметра в группах одинаковое, группы не различаются
\(H_1:\) среднее значение параметра в группах разное, группы различаются

t.test(data[which(data$T == 1),]$X, data[which(data$T == 0),]$X)


    Welch Two Sample t-test

data:  data[which(data$T == 1), ]$X and data[which(data$T == 0), ]$X
t = -0.21924, df = 997.51, p-value = 0.8265
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.2777933  0.2219581
sample estimates:
mean of x mean of y 
 21.46699  21.49490

t.test(data[which(data$T == 1),]$Z, data[which(data$T == 0),]$Z)


    Welch Two Sample t-test

data:  data[which(data$T == 1), ]$Z and data[which(data$T == 0), ]$Z
t = -0.45032, df = 997.68, p-value = 0.6526
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -3.057656  1.916238
sample estimates:
mean of x mean of y 
 59.95331  60.52402

Согласно значению p-value мы принимаем \(H_0\) и делаем вывод, что группы одинаковые.

1.3.5.3 Использование готовых пакетов

library("tableone")
table <- CreateTableOne(vars=c("X", "Z"), strata="T", data=data, test=TRUE)
print(table)

               Stratified by T
                0             1             p      test
  n               500           500                    
  X (mean (SD)) 21.49 (2.04)  21.47 (1.99)   0.827     
  Z (mean (SD)) 60.52 (19.86) 59.95 (20.22)  0.653

1.4 Про важность предпосылок Гаусса-Маркова

1.4.1 Предпосылки Гаусса-Маркова

Линейная модель: \(Y_i=\tau T_i+\beta X_i+\varepsilon_i\)
\(T_i\) не выражается линейно через \(X_i\)
\(\mathbb{E}\left(\varepsilon_i\right)=0\)
Гомоскедастичность, некоррелированность ошибок \(\mathbb{V}\left(\varepsilon_i\right)=\sigma I\)
Экзогенность \(\operatorname{Cov}\left(X_i, \varepsilon_i\right)=0\)
Оцениваем методом наименьших квадратов: \[ \left(Y_i-\hat{\tau} T_i-\hat{\beta} X_i\right)^2 \rightarrow \min _{\hat{\tau}, \hat{\beta}} \]
Свойства оценки (теорема Гаусса-Маркова): несмещенность, эффективность, состоятельность

1.4.2 Пример

Мы располагаем данными из 1000 наблюдений, про которые известно:

\(X_i \sim U[0,1]\) - контрольная переменная: успехи команды в прошлом году
\(T_i \mid X_i \sim \operatorname{Bernoulli}\left(0.1+0.8 X_i\right)\) - бинарная переменная переменная интереса: решение тренера взять команду
\(Y_i=f\left(X_i, T_i\right)+\varepsilon_i\) - переменная исхода: успехи команды в этом году
В соответствии с предпосылками оценим \(Y_i \sim T_i+X_i\)
То же ли самое получится, если \(T_i \sim \operatorname{Bernoulli}(0.5)\) и оценим \(Y_i \sim T_i\)?

1.4.2.1 Генерация данных

Сгенерим данные согласно условию:

library(stargazer)

set.seed(123)

X <- runif(1000, 0, 1)

Y0 <- rnorm(1000)
Y1 <- 10 * X + rnorm(1000)

1.4.2.2 Генерация тритмента

Случай 1: эндогенный тритмент

T_one = rbinom(1000, 1, 0.1 + 0.8*X)
Y_one = Y1 * T_one + Y0 * (1 - T_one)

data_one = data.frame(X=X, Y=Y_one, T=T_one)

Случай 2: экзогенный тритмент

T_two = rbinom(1000, 1, 0.5)
Y_two = Y1 * T_two + Y0 * (1 - T_two)

data_two = data.frame(X=X, Y=Y_two, T=T_two)

1.4.2.3 Оценка эффекта

Сравним оценки:

model_one_1 <- lm(Y ~ T, data=data_one)
model_one_2 <- lm(Y ~ T + X, data=data_one)
model_two_1 <- lm(Y ~ T, data=data_two)
model_two_2 <- lm(Y ~ T + X, data=data_two)

#stargazer(model_one_1, model_one_2, model_two_1, model_two_2, type = 'text') 

# В обычной жизни вы просто можете пользоваться старгейзером, как в комментраии выше
# Чтобы конкретно на этой странице влезло 4 регрессии, я воспользуюсь другой функцией msummary
m_list <- list(endogenous_without_control = model_one_1, 
               endogenous_with_control = model_one_2,
               exogenous_without_control = model_two_1, 
               exogenous_with_control = model_two_2)
library('modelsummary')
msummary(m_list, output = 'markdown', stars = TRUE, title = 'Сравнение моделей')

Сравнение моделей
	endogenous_without_control	endogenous_with_control	exogenous_without_control	exogenous_with_control
(Intercept)	0.022	-1.847***	0.030	-2.234***
	(0.094)	(0.102)	(0.096)	(0.120)
T	6.260***	4.859***	5.069***	5.007***
	(0.131)	(0.114)	(0.137)	(0.109)
X		5.206***		4.612***
		(0.199)		(0.189)
Num.Obs.	1000	1000	1000	1000
R2	0.695	0.820	0.577	0.735
R2 Adj.	0.695	0.819	0.577	0.735
AIC	4299.4	3778.1	4390.7	3924.6
BIC	4314.2	3797.7	4405.4	3944.3
Log.Lik.	-2146.714	-1885.056	-2192.343	-1958.311
F	2279.232	2264.086	1363.774	1385.366
RMSE	2.07	1.59	2.17	1.71

Note: ^^ + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

Выводы:

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

1.4.2.4 Баланс ковариатов

\(H_0:\) среднее значение параметра в группах одинаковое, группы не различаются
\(H_1:\) среднее значение параметра в группах разное, группы различаются

table1 <- CreateTableOne(vars=c("X", "Z"), strata="T", data=data_one, test=TRUE)

Warning in ModuleReturnVarsExist(vars, data): The data frame does not have: Z
Dropped

print(table1)

               Stratified by T
                0           1           p      test
  n              486         514                   
  X (mean (SD)) 0.36 (0.25) 0.63 (0.26) <0.001

table2 <- CreateTableOne(vars=c("X", "Z"), strata="T", data=data_two, test=TRUE)

Warning in ModuleReturnVarsExist(vars, data): The data frame does not have: Z
Dropped

print(table2)

               Stratified by T
                0           1           p      test
  n              514         486                   
  X (mean (SD)) 0.49 (0.29) 0.50 (0.29)  0.463

Согласно значению p-value мы принимаем \(H_0\) для второго случая, где тритмент экзогенный, и делаем вывод, что группы одинаковые; для первого случая с эндогенным тритментом баланс ковариатов не выполняется.