Использованные данные

В 1988 году в Калифорнии был принят знаменитый «Закон о налогообложении табака и здравоохранении», который получил известность под названием «Предложение 99». «Основной эффект от этого закона заключался во взимании акцизного налога в размере 25 центов с каждой проданной пачки сигарет в пределах штата Калифорния, с примерно эквивалентными акцизными налогами на розничную продажу другой табачной продукции, в частности сигар и жевательного табака. Среди дополнительных ограничений, налагаемых на продажу табака — запрет на установку вендинговых автоматов с табачной продукцией в общественных местах, открытых для детей, а также запрет на поштучную продажу сигарет».

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

import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
import statsmodels.formula.api as smf
from sklearn import linear_model

%matplotlib inline

cig = (pd.read_csv("/content/drive/MyDrive/smoking.csv")
         .drop(columns=["lnincome","beer", "age15to24"]))

style.use('fivethirtyeight') # You can change this to your preference

cig.query("california").head()

Фрагмент датасета

Здесь state — индекс/номер штата, у Калифорнии номер 3. В качестве ковариат берём розничную стоимость сигарет и подушевые продажи сигарет в пересчёте на пачки. Переменная, интересующая нас на выходе, — cigsale. Наконец, у нас есть вспомогательные булевы переменные, по которым можно судить о штате Калифорния и периоде после вмешательства. Вот что получится, если построить график продаж сигарет в Калифорнии и график средних продаж сигарет в других штатах в зависимости от времени.


В течение того периода, данными за который мы располагаем, калифорнийцы явно покупали меньше сигарет, чем в целом по стране. Кроме того, просматривается тренд на снижение потребления табака после 80-х. Представляется, что после принятия «Предложения 99» этот тренд на снижение в Калифорнии усилился (в сравнении с другими штатами), но с уверенностью этого утверждать нельзя. Это просто версия, которую нам требуется проверить по графику.

Чтобы ответить на вопрос о том, повлияло ли «Предложение 99» на потребление сигарет, изучим период до вмешательства и выстроим синтетический контроль. Мы скомбинируем данные по другим штатам, чтобы собрать на их основе фиктивный штат, тренды в котором очень напоминают калифорнийские. Затем посмотрим, как картина в этом синтетическом контрольном штате изменится после вмешательства.

Математическая нотация

Допустим, у нас J+1 единиц. Далее, без ущерба нашим обобщениям, предположим, что единицу 1 затрагивает рассматриваемое вмешательство, то есть «Предложение 99». Единицы j=2,…, J+1 — это коллекция незатронутых единиц (штатов), которые мы в совокупности назовём «пулом доноров». Также предположим, что имеющиеся у нас данные охватывают T периодов с T0 периодов до вмешательства. Для каждой единицы j и каждого интервала t наблюдается результат Yjt. Для каждой единицы j и каждого периода t определим YNjt как потенциальный результат, наступающий при отсутствии вмешательства, и YIjt — как потенциальный результат, наступающий после вмешательства.


Тогда эффект затронутой единицы j=1 во время t для t>T0 определяется как

Поскольку j=1 затронута вмешательством, YIjt наблюдается на практике, а YNjt — нет. Теперь возникает вопрос: как нам оценить YNjt. Обратите внимание: эффект вмешательства определяется для каждого периода, то есть со временем может меняться. Он не должен проявляться мгновенно, а может накапливаться или рассеиваться. Чтобы показать это нагляднее, укажем, что задача оценки эффекта вмешательства сводится к следующему: оценить, что произошло бы, если бы выход от единицы j=1 был получен без вмешательства.

Для оценки YNjt не забудем, что комбинация единиц в пуле доноров может аппроксимировать характеристики затрагиваемой единицы гораздо лучше, чем любая незатронутая вмешательством единица как таковая. Следовательно, синтетический контроль определяется как взвешенное среднее от всех единиц в контрольном пуле. Имея веса W=(w2,…,wJ+1), получаем через синтетический контроль оценку YNjt равной

Наглядное объяснение

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

При применении синтетического контроля у нас не так много единиц, зато много временных отрезков. В таком случае мы просто переворачиваем входную матрицу. Тогда единицы становятся «переменными», и мы представляем итог как средневзвешенное значение от единиц — см. следующее перемножение матриц.

Если у нас имеется более чем по одному признаку на период времени, то можно группировать признаки следующим образом. Здесь самое важное добиться, чтобы регрессия «пыталась предсказать» затронутую единицу 1, используя другие единицы. Поэтому можно подобрать веса оптимальным образом, чтобы достичь нужной нам степени приближения. Мы можем даже по-разному отмасштабировать признаки, присвоив им разную степень важности.

Реализация

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

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

Для этого нам первым делом требуется преобразовать единицы (в нашем случае — штаты) в столбцы, а время — в строки. Поскольку у нас 2 признака, cigsale и retprice, расположим их друг над другом, как было сделано на предыдущей картинке. Соберём синтетический контрольный штат, который до вмешательства выглядел примерно как Калифорния, и посмотрим, что с ним произойдёт в период после вмешательства. Поэтому важно выбрать лишь период до вмешательства. Здесь признаки откладываются в похожем масштабе, поэтому мы их не нормализуем. Если бы масштабы у признаков отличались, например, один исчислялся тысячами, а другой десятками, то при минимизации разницы нам был бы наиболее важен больший признак. Чтобы этого избежать, важно сначала отмасштабировать признаки.

features = ["cigsale", "retprice"]

inverted = (cig.query("~after_treatment") # отфильтруем период до вмешательства
            .pivot(index='state', columns="year")[features] # предусмотрим по столбцу на год и по строке на штат 
            .T) # перевернём таблицу так, чтобы у нас было по столбцу на штат 

inverted.head()

Теперь можно определить переменную Y как штат Калифорния, а переменную X как другие штаты.

y = inverted[3].values # штат Калифорния
X = inverted.drop(columns=3).values  # другие штаты

Далее выполним лассо-регрессию. Лассо-регрессия, также именуемая L1, используется для того, чтобы избежать переобучения наших данных по штатам. Для этого также можно воспользоваться ридж-регрессией. Регрессия вернёт набор весов, минимизирующих квадратичную разницу между рассматриваемой единицей и единицами из пула доноров.

from sklearn.linear_model import Lasso
weights_lr = Lasso(fit_intercept=False).fit(X, y).coef_
weights_lr.round(3)

По этим весам видно, как построить синтетический контроль. Умножим итог штата 1 на 0,566, штата 3 на 0,317, штата 4 на 0,158 и т. д. Этого можно добиться скалярным произведением между матрицей штатов в пуле и весами.

calif_synth_lr = (cig.query("~california")
                  .pivot(index='year', columns="state")["cigsale"]
                  .values.dot(weights_lr))

Теперь, когда синтетический контроль готов, можно нанести его на график вместе с итоговой переменной для штата Калифорния.

plt.figure(figsize=(10,6))
plt.plot(cig.query("california")["year"], cig.query("california")["cigsale"], label="California")
plt.plot(cig.query("california")["year"], calif_synth_lr, label="Synthetic Control")
plt.vlines(x=1988, ymin=40, ymax=140, linestyle=":", lw=2, label="Proposition 99")
plt.ylabel("Gap in per-capita cigarette sales (in packs)")
plt.legend();

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

plt.figure(figsize=(10,6))
plt.plot(cig.query("california")["year"], cig.query("california")["cigsale"] - calif_synth_lr,
         label="California Effect")
plt.vlines(x=1988, ymin=-30, ymax=7, linestyle=":", lw=2, label="Proposition 99")
plt.hlines(y=0, xmin=1970, xmax=2000, lw=2)
plt.title("State - Synthetic Across Time")
plt.ylabel("Gap in per-capita cigarette sales (in packs)")
plt.legend();

Представляется, что к 2000 году благодаря «Предложению 99» подушевые продажи сигарет снизились на 25 пачек. Теперь нужно выяснить, является ли этот показатель статистически значительным.

Логический вывод

Здесь вооружимся точным тестом Фишера. Интуитивно он совершенно понятен. Мы основательно переставляем единицы в затронутой и контрольной части. Поскольку в затронутой части у нас всего одна единица, это будет означать, что каждую единицу мы будем считать обрабатываемой, а все остальные в то же время — контрольными.

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

Эта функция возвращает фрейм данных, где один столбец соответствует штату, один – году, один — итогу cigsale и один — синтетическому итогу для данного штата.

Вот что получается при применении функции к первому штату:

def synthetic_control(state, pool, data) -> np.array:
    
    features = ["cigsale", "retprice"]
    
    inverted = (data.query("~after_treatment")
                .pivot(index='state', columns="year")[features]
                .T)
    
    y = inverted[state].values # treated
    X = inverted.drop(columns=state).values # donor pool

    weights = Lasso(fit_intercept=False).fit(X, y).coef_.round(3)
    synthetic = (data.query(f"~(state=={state})")
                 .pivot(index='year', columns="state")["cigsale"]
                 .values.dot(weights))

    return (data
            .query(f"state=={state}")[["state", "year", "cigsale", "after_treatment"]]
            .assign(synthetic=synthetic))

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

sinthetic_states = [synthetic_control(state, control_pool, cig) for state in control_pool]
plt.figure(figsize=(12,7))
for state in sinthetic_states:
    plt.plot(state["year"], state["cigsale"] - state["synthetic"], color="C5",alpha=0.4)

plt.plot(cigar.query("california")["year"], cigar.query("california")["cigsale"] - calif_synth_lr,
        label="California");

plt.vlines(x=1988, ymin=-50, ymax=120, linestyle=":", lw=2, label="Proposition 99")
plt.hlines(y=0, xmin=1970, xmax=2000, lw=3)
plt.ylabel("Gap in per-capita cigarette sales (in packs)")
plt.title("State - Synthetic Across Time")
plt.legend();

------

control_pool = cig["state"].unique()

sinthetic_states = [synthetic_control(state, control_pool, cig) for state in control_pool]

plt.figure(figsize=(12,7))
for state in sinthetic_states:
    plt.plot(state["year"], state["cigsale"] - state["synthetic"], color="C5",alpha=0.4)

plt.plot(cig.query("california")["year"], cig.query("california")["cigsale"] - calif_synth_lr,
        label="California");

plt.vlines(x=1988, ymin=-50, ymax=120, linestyle=":", lw=2, label="Proposition 99")
plt.hlines(y=0, xmin=1970, xmax=2000, lw=3)
plt.ylabel("Gap in per-capita cigarette sales (in packs)")
plt.title("State - Synthetic Across Time")
plt.legend();

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

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

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

def pre_treatment_error(state):
    pre_treat_error = (state.query("~after_treatment")["cigsale"] 
                       - state.query("~after_treatment")["synthetic"]) ** 2
    return pre_treat_error.mean()

plt.figure(figsize=(12,7))
for state in sinthetic_states:
    
    # удалить единицы со средней ошибкой выше 5.
    if pre_treatment_error(state) < 5:
        plt.plot(state["year"], state["cigsale"] - state["synthetic"], color="C5",alpha=0.4)

plt.plot(cig.query("california")["year"], cig.query("california")["cigsale"] - calif_synth_lr,
        label="California");

plt.vlines(x=1988, ymin=-50, ymax=120, linestyle=":", lw=2, label="Proposition 99")
plt.hlines(y=0, xmin=1970, xmax=2000, lw=3)
plt.ylabel("Gap in per-capita cigarette sales (in packs)")
plt.title("Distribution of Effects")
plt.title("State - Synthetic Across Time (Large Pre-Treatment Errors Removed)")
plt.legend();

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


Эффект от вмешательства в штате Калифорния в 2000 году составляет -31.419, то есть в результате вмешательства подушевое потребление сигарет уменьшилось почти на 31 пачку.

Заключение

Пользуясь данными по другим штатам, отражающими период до вмешательства, мы построили модель на основе лассо-регрессии, в которой присвоили фиксированные веса каждому контрольному штату и получили средневзвешенное значение, которое очень напоминает характеристику курения для штата Калифорния до того, как было принято «Предложение 99».

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

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

Книга

Causal Inference for The Brave and True by Matheus Facure