合成控制法（Synthetic Control Method）笔记及R/Python实现

kkkkkpaboo

已于 2024-04-30 21:50:25 修改

阅读量2.9k

点赞数 33

文章标签：笔记 python r语言

于 2024-04-30 21:44:43 首次发布

本文链接：https://blog.csdn.net/kkkkkangel/article/details/137747827

版权

最近在写本科毕业论文，想要探究中国区域政策对某个商品销量的具体影响，于是有幸了解到合成控制法（DID做不到但是它做到了！），记录一下学习笔记，以免之后忘记。

本文主要翻译至Abadie 等（2010）提出的合成控制法https://doi.org/10.1198/jasa.2009.ap08746

R代码引用GitHub - edunford/tidysynth: A tidy implementation of the synthetic control method in R

Python代码引用tom-beer/Synthetic-Control-Examples: Step by step code for applying synthetic control to interesting use cases (github.com)

1.引入

合成控制法（Synthetic Control Method, SCM）是一种用于因果推断的统计学方法，适用于评估政策干预或事件对特定单位（如国家、企业等）影响的情况。该方法通过构造一个“合成控制组”，模拟政策干预或事件未发生时的情况，从而估计干预的因果效应。

那它与双重差分法（DID）的区别在哪里？或者说优势在哪里？双重差分法通常对处理组（受到干预的单元）和控制组（未受到干预的单元）有一定的要求：平行假设等。但是在一些情况下，我们无法寻找到在干预前具有相同趋势的处理组和控制组（比如：我的本科毕业论文就不适用555），而合成控制法利用控制组中各个单元的权重加和去构建一个“合成控制组”去拟合出一个与控制组相似的单元。（思路打开芜湖）

2.原理

（开始硬核翻译论文）

假设：在T时期有J+1个可观测的样本，其中不妨假设第1个样本在 $T_0(1\leq T_0\leq T )$ 时刻受到政策干预，即为处理组，剩余J个样本为控制组。

令 $Y_{it}^{N}$ 为在没有干预下样本i在t时点的观测值， $Y_{it}^{I}$ 为在干预下样本i在t时点的观测值；

我们假设在受到干预前，干预对样本的观测值没有影响，即 $Y_{it}^{I}=Y_{it}^{N},i=1,2,...,J+1,t=1,2,..,T_0$ 。

于是，样本i在t时点受到干预的影响可以表示为 $\alpha _{it}=Y_{it}^{I}-Y_{it}^{N},t=T_0+1,...,T$ 。其中 $Y_{it}^{I}$ 是可观测的， $Y_{it}^{N}$ 未知，于是作者通过控制组中各个单元的权重组合得到。

下面讲讲为什么这样可以。

假设：

$Y_{it}^N=\delta_t+\theta_tZ_i+\lambda_t\mu_i+\! \varepsilon_{it} (1)$

其中 $\delta_t$ 为未知的固定个体效应， $Z_i$ 为一组可观测且不受干预影响的协变量构成的 $(r\times1)$ 维向量， $\theta_t$ 为未知系数构成的 $(1\times r)$ 维向量， $\lambda _t$ 为未知的公共因子构成的 $(1\times F)$ 维向量， $\mu _i$ 为因子载荷构成的 $(F\times1)$ 维向量， $\varepsilon_{it}$ 为不可观测的均值为0的随机误差项。

假设存在 $(J\times1)$ 维权重向量 $W=(w_2,..,w_{J+1})'$ 满足 $w_j\geq 0,j=2,...,J+1$ 且 $w_2+...+w_{J+1}=1$ 。此时，加权后的公式(1)可以表示为：

$\sum_{j=2}^{J+1}w_jY_{jt}=\delta_t+\theta_t\sum_{j=2}^{J+1}w_jZ_j+\lambda_t\sum_{j=2}^{J+1}w_j\mu_j+\sum_{j=2}^{J+1}w_j\varepsilon_{jt}$

假设存在 $W^*=(w_2^*,...,w_{J+1}^*)'$ 满足：

$\sum_{j=2}^{J+1}w_j^*Y_{j1}=Y_{11},\sum_{j=2}^{J+1}w_j^*Y_{j2}=Y_{12},...,\sum_{j=2}^{J+1}w_j^*Y_{jT_0}=Y_{1T_0},$ $\sum_{j=2}^{J+1}w_j^*Z_{j}=Z_{1}$

当 $\sum_{t=1}^{T_0}\lambda_t'\lambda_t$ 是非奇异，干预前的时间段大于干预后的时间段或者干预前只有1期数据时，下图公式右边趋近于0。

因此，我们可以利用 $Y_{1t}^I-\sum_{j=2}^{J+1}w_j^*Y_{jt},t=T_0+1,...,T$ 去估计 $\alpha_{1t}$ ！

OK那我们现在要解决的就是怎么得到 $W^*=(w_2^*,...,w_{J+1}^*)'$ ，即怎么求解最优权重使得控制组在干预前的观测值最接近处理组。

设 $X_1$ 为处理组受到干预前协变量的值构成的 $(k\times 1)$ 维向量， $X_0$ 为控制组受到干预前协变量的值构成的 $(k\times J)$ 维矩阵（即公式(1)中的Z）， $V$ 为 $(k\times k)$ 的对角线代表不同协变量的相对重要性的半正定对称矩阵， $Z_1$ 为处理组受到干预前的结果变量构成的 $(T_0\times 1)$ 维向量， $Z_0$ 为控制组受到干预前的结果变量构成的 $(T_0\times J)$ 维矩阵（即公式(1)中的Y，应用中的人均香烟销售量），则

$W^*=argmin\left \| X_1-X_0W_V \right \|=argmin\sqrt{(X_1-X_0W)'V(X_1-X_0W)}$

$V^*=argmin(Z_1-Z_0W^*)'(Z_1-Z_0W^*)$

以下R和python代码求解最优权重W*和V*的具体方法如下：

给定权重矩阵V初始值，通过最小化 $(X_1-X_0W)'V(X_1-X_0W)$ 求解控制组权重W*（V）
将上述求出的W*（V）代入 $(Z_1-Z_0W)'V(Z_1-Z_0W)$
通过最小化 $(Z_1-Z_0W^*(v)))'V(Z_1-Z_0W^*(V)))$ 求解最优变量权重矩阵V*
将最优变量权重矩阵V*代入 $(X_1-X_0W)'V(X_1-X_0W)$ ，通过最小化 $(X_1-X_0W)'V^*(X_1-X_0W)$ 得到最优控制组权重矩阵W*

3.应用

1988年11月，加州通过了著名的99号提案：烟草税与健康保护法。其主要内容是对加州的卷烟销售征收每包25美分的州消费税，对其他商业类烟草产品的销售也征收约同等的消费税。烟草销售的其他限制包括禁止在青少年可进入的公共区域安装自动售烟机，禁止个人销售香烟。同时，该法案产生的收入专门用于各种环境和保健计划以及反烟草广告。

由于数据的可获取性和2000年开始大部分州实施反烟草措施，作者收集1970-2000年39个州面板数据探究99号提案对人均卷烟销售量的具体影响。其中，作者保证控制组（38个州）在1989年-2000年没有实施类似法案（我们在使用合成控制法的时候也要注意这一点，去选择合适的控制单元）。

由上图可以看出99号提案实施之后，加州相比于其他州人均卷烟销量下降幅度更大。为了评估99号提案对加州人均卷烟销量的影响，核心问题在于如果没有实施99号提案，1988年后加州的卷烟消费将如何变化。合成控制方法则提供了一种系统的方法来估计这种反事实。

数据链接：

'smoking.rda' data taken from mixtape/data at master · johnson-shuffle/mixtape · GitHub

'prop99.csv' taken from tslib/tests/testdata at master · jehangiramjad/tslib · GitHub

#也可以从R中tidysynth包直接导入
require(tidysynth)
data("smoking")

3.1 R语言实现（推荐）

# 安装包
install.packages('tidysynth')
# 导入数据并查看
require(tidysynth)
data("smoking")
smoking %>% dplyr::glimpse()
#输出：
## Rows: 1,209
## Columns: 7
## $ state     <chr> "Rhode Island", "Tennessee", "Indiana", "Nevada", "Louisiana…
## $ year      <dbl> 1970, 1970, 1970, 1970, 1970, 1970, 1970, 1970, 1970, 1970, …
## $ cigsale   <dbl> 123.9, 99.8, 134.6, 189.5, 115.9, 108.4, 265.7, 93.8, 100.3,…
## $ lnincome  <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ beer      <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ age15to24 <dbl> 0.1831579, 0.1780438, 0.1765159, 0.1615542, 0.1851852, 0.175…
## $ retprice  <dbl> 39.3, 39.9, 30.6, 38.9, 34.3, 38.4, 31.4, 37.3, 36.7, 28.8, …

变量解释：state为州名，year为年份，cigsale为卷烟销售，lnincome为对数人均可支配收入，beer为人均啤酒消费量，age15to24为年龄在15-24岁之间的人口比例，retprice为卷烟平均零售价格。

#合成控制
smoking_out <-
  
  smoking %>%
  
  # initial the synthetic control object
  synthetic_control(outcome = cigsale, # 结果变量Y
                    unit = state, # 面板数据中的单元
                    time = year, # 时间
                    i_unit = "California", # 处理组
                    i_time = 1988, # 干预实施时间
                    generate_placebos=T # 后续是否做安慰剂检验
                    ) %>%
  
  # Generate the aggregate predictors used to fit the weights

  # 设置协变量为：1980-1988年人均收入的对数变换、卷烟平均销售价格、年龄为15-24岁人口比例、
  # 1984-1988年人均啤酒消费量、1975年卷烟销量、1980年卷烟销量、1988年卷烟销量
  # average log income, retail price of cigarettes, and proportion of the
  # population between 15 and 24 years of age from 1980 - 1988
  generate_predictor(time_window = 1980:1988,
                     ln_income = mean(lnincome, na.rm = T),
                     ret_price = mean(retprice, na.rm = T),
                     youth = mean(age15to24, na.rm = T)) %>%
  
  # average beer consumption in the donor pool from 1984 - 1988
  generate_predictor(time_window = 1984:1988,
                     beer_sales = mean(beer, na.rm = T)) %>%
  
  # Lagged cigarette sales 
  generate_predictor(time_window = 1975,
                     cigsale_1975 = cigsale) %>%
  generate_predictor(time_window = 1980,
                     cigsale_1980 = cigsale) %>%
  generate_predictor(time_window = 1988,
                     cigsale_1988 = cigsale) %>%
  
  
  # Generate the fitted weights for the synthetic control
  generate_weights(optimization_window = 1970:1988, # 干预前时间段
                   margin_ipop = .02,sigf_ipop = 7,bound_ipop = 6 # optimizer options
  ) %>%
  
  # 开始利用控制组的协变量数据合成加州
  generate_control()

# 输出干预前真实加州、合成加州、控制组的协变量表格
smoking_out %>% grab_balance_table()

输出表格第一列、第二列、第三列分别是真实加州、合成加州、38个控制单元平均的协变量值。由表中数据可以看出相比于第三列，第二列的值与第一列的值更接近。这可以说明控制单元的加权组合更好的再现了干预前的加州，因此这个加权组合在干预后的表现可以看作没有实施干预的加州的表现。

# 输出协变量的相对重要性
smoking_out %>% grab_predictor_weights()

以上结果为协变量的相对重要性，即使干预前真是加州和合成加州人均卷烟销售量最小的矩阵V的对角线值。可以看出，人均收入的对数变换对人均卷烟销量的预测能力较低，重要性较低。

# 输出每个控制单元的权重
smoking_out %>%grab_unit_weights()

上图结果为最优控制单元权重矩阵W*。在99号提案通过之前，加利福尼亚州的吸烟趋势主要由科罗拉多州、康涅狄格州、蒙大拿州、内华达州和犹他州的组合来重现，其他控制单元权重都为0。

# 绘制真实加州和合成加州的人均卷烟销量趋势图
smoking_out %>% plot_trends()

上图即为2中的 $Y_{it}^{I}$ 和 $Y_{it}^{N}$ ，可以看出在99号提案通过之前合成加州的人均卷烟销量与真实加州几乎重合，这可以说明控制单元的加权组合是1970-2000年真实加州的一个合理近似。并且在99号提案通过之后加州的人均卷烟销量迅速下降，可以说明99号提案对人均卷烟销量具有负面影响。此外，我们可以通过以下函数输出 $\alpha _{it}=Y_{it}^{I}-Y_{it}^{N},t=T_0+1,...,T$ 具体值

#真实加州和合成加州差值趋势图
smoking_out %>% plot_differences()

为检验上述结果（99号提案对人均卷烟销量具有负面影响）并非偶然，作者使用了安慰剂检验，即迭代地将控制组中的单元作为处理组，原本的处理组加州作为控制单元进行合成控制，然后绘制差值趋势图。代码和结果如下：

# 安慰剂检验（默认删除干预前拟合较差的安慰剂检验）
smoking_out %>% plot_placebos()

图中紫线为加州作为处理组的合成控制结果，灰线为安慰剂检验结果。可以看出相较于灰线，紫线下降趋势更为明显。

此外，作者还提出通过计算加州和控制组干预后/干预前MSPE检验99号提案对人均卷烟销量具有负面影响并非偶然。

# 绘制干预后/前MSPE
smoking_out %>% plot_mspe_ratio()
# 干预后/前MSPE具体数值及P值
# smoking_out %>% grab_significance()

加州干预后/前MSPE远超控制单元，且出现这么大比值的概率为1/39=0.026，小于显著性水平0.05。因此，上述所有结果可以证明99号提案对人均卷烟销量具有负面影响并非偶然。

3.2 Python语言实现

（此部分不含安慰剂检验，如果有请评论补充）

准备数据：

from matplotlib import pyplot as plt
import pyreadr
import numpy as np
import pandas as pd
import copy
from scipy.optimize import fmin_slsqp
from sklearn.metrics import mean_squared_error

START_TIME = 1970 #开始时间
INTERVENTION_TIME = 1989 #干预时间
STOP_TIME = 2001 #结束时间

# 导入数据
df_outcome_raw = pd.read_csv('../Data/prop99.csv')
df_outcome_raw = df_outcome_raw[df_outcome_raw['SubMeasureDesc'] == 'Cigarette Consumption (Pack Sales Per Capita)']
df_outcome = pd.DataFrame(df_outcome_raw.pivot_table(values='Data_Value', index='LocationDesc', columns=['Year']).to_records()) #结果变量

rda_predictors = pyreadr.read_r('../Data/smoking.rda')
df_predictors = pd.DataFrame(list(rda_predictors.values())[0]) #协变量

# 按照论文3.2部分筛选出符合条件的州作为控制组和处理组
print(f'In the original dataset there are {df_outcome.LocationDesc.unique().shape[0]} states')
# Section 3.2 in the paper
bad_states = ['Massachusetts', 'Arizona', 'Oregon', 'Florida', 'Alaska', 'Hawaii', 'Maryland', 
              'Michigan', 'New Jersey', 'New York', 'Washington', 'District of Columbia']

df_outcome.drop(df_outcome[df_outcome['LocationDesc'].isin(bad_states)].index, inplace=True)
ca_id = df_outcome[df_outcome['LocationDesc'] == 'California'].index.item()
df_outcome = df_outcome.reset_index()
df_outcome = df_outcome.rename(columns={'index': 'org_index'})
print(f'After filtering out some states, we are left with {df_outcome.LocationDesc.unique().shape[0]} states (including California):')
df_outcome.head()

df_outcome_ca = df_outcome.loc[df_outcome['LocationDesc'] == 'California', :]
df_outcome_control = df_outcome.loc[df_outcome['LocationDesc'] != 'California', :]

ca_outcomes_pre = df_outcome_ca.loc[:,[str(i) for i in list(range(START_TIME, INTERVENTION_TIME))]].values.reshape(-1,1)
control_outcomes_pre = df_outcome_control.loc[:,[str(i) for i in list(range(START_TIME, INTERVENTION_TIME))]].values.transpose()

ca_outcomes_post = df_outcome_ca.loc[:,[str(i) for i in list(range(INTERVENTION_TIME, STOP_TIME))]].values.reshape(-1,1)
control_outcomes_post = df_outcome_control.loc[:,[str(i) for i in list(range(INTERVENTION_TIME, STOP_TIME))]].values.transpose()

Z0 = control_outcomes_pre  #干预之前的控制组的结果变量
Z1 = ca_outcomes_pre #干预之前的处理组的结果变量
Y0 = control_outcomes_post #干预之后的控制组的结果变量
Y1 = ca_outcomes_post #干预之后的处理组的结果变量

# 绘制控制组和试验组的结果变量
mean_outcomes = np.vstack([Z0, Y0]).mean(axis=1)
CA_outcomes = np.vstack([Z1, Y1]).flatten()
fig = plt.figure(figsize=(7.5,4.5))
plt.plot(range(START_TIME,STOP_TIME),mean_outcomes, 'r--', label="rest of the U.S.");
plt.plot(range(START_TIME,STOP_TIME),CA_outcomes, 'b-', label="California");

plt.ylabel('per-capita cigarette sales (in packs)')
plt.xlabel('year')
plt.legend(loc='upper right')
plt.title("Figure 1: Trends in per-capita cigarette sales: California vs. the rest of the United States")
plt.axvline(INTERVENTION_TIME)
plt.text(x=INTERVENTION_TIME+0.2, y=30, s='Passage of Proposition 99')
plt.xlim([START_TIME, STOP_TIME-1])
plt.ylim([0, 150])
plt.grid()
plt.show()
fig.savefig("name", dpi=300)

# 生成协变量
def extract_predictor_vec(state):
    df_outcome_state = df_outcome[df_outcome['LocationDesc'] == state]
    cigsale88_predictor = df_outcome_state['1988'].item()
    cigsale80_predictor = df_outcome_state['1980'].item()
    cigsale75_predictor = df_outcome_state['1975'].item()
    
    state_id_predictors_df = df_outcome[df_outcome['LocationDesc'] == 'California'].index.item() + 1
    df_predictors_state = df_predictors[df_predictors['state'] == state_id_predictors_df]
    beer_predictor = df_predictors_state.loc[(df_predictors_state['year'] >= 1984) & (df_predictors_state['year'] < INTERVENTION_TIME), 'beer'].mean()
    age15to24_predictor = df_predictors_state.loc[(df_predictors_state['year'] >= 1980) & (df_predictors_state['year'] < INTERVENTION_TIME), 'age15to24'].mean()*100  # Should I turn multiply by 100? In table 1 it looks like it
    retprice_predictor = df_predictors_state.loc[(df_predictors_state['year'] >= 1980) & (df_predictors_state['year'] < INTERVENTION_TIME), 'retprice'].mean()
    lnincome_predictor = df_predictors_state.loc[(df_predictors_state['year'] >= 1980) & (df_predictors_state['year'] < INTERVENTION_TIME), 'lnincome'].mean()
    
    return np.array([lnincome_predictor, age15to24_predictor, retprice_predictor, beer_predictor,  
                     cigsale88_predictor, cigsale80_predictor, cigsale75_predictor]).reshape(-1,1)


control_predictors = []
for state in df_outcome['LocationDesc'].unique():
    print(state)
    state_predictor_vec = extract_predictor_vec(state)
    print(state_predictor_vec)
    if state == 'California':
        ca_predictors = state_predictor_vec
    else:
        control_predictors += [state_predictor_vec]

control_predictors = np.hstack(control_predictors)

X0 = control_predictors #控制组预测变量
X1 = ca_predictors #处理组预测变量

定义优化函数，求解最优权重W*,V*

def w_mse(w, v, x0, x1): return mean_squared_error(x1, x0.dot(w), sample_weight=v)

# 约束控制单元权重之和为1
def w_constraint(w, v, x0, x1): return np.sum(w) - 1 

# 约束协变量权重之和为1
def v_constraint(V, W, X0, X1, Z0, Z1): return np.sum(V) - 1

def fun_w(w, v, x0, x1): return fmin_slsqp(w_mse, w, bounds=[(0.0, 1.0)]*len(w), f_eqcons=w_constraint, args=(v, x0, x1), disp=False, full_output=True)[0]

def fun_v(v, w, x0, x1, z0, z1): return mean_squared_error(z1, z0.dot(fun_w(w, v, x0, x1)))

def solve_synthetic_control(X0, X1, Z0, Z1, Y0):
    k,j = X0.shape
    V0 = 1/k*np.ones(k)
    W0 = 1/j*np.zeros(j).transpose()
    V = fmin_slsqp(fun_v, V0, args=(W0, X0, X1, Z0, Z1), bounds=[(0.0, 1.0)]*len(V0), disp=True, f_eqcons=v_constraint, acc=1e-6)
    W = fun_w(W0, V, X0, X1)
    return V, W

V, W = solve_synthetic_control(X0, X1, Z0, Z1, Y0)

# 绘制真实加州和合成加州的趋势图
SC_outcomes = np.vstack([Z0, Y0]).dot(W)
CA_outcomes = np.vstack([Z1, Y1]).flatten()
fig = plt.figure(figsize=(7.5,4.5)) 

plt.plot(range(START_TIME,STOP_TIME),SC_outcomes, 'r--', label="Synthetic California");
plt.plot(range(START_TIME,STOP_TIME),CA_outcomes, 'b-', label="California");

plt.ylabel('per-capita cigarette sales (in packs)')
plt.xlabel('year')
plt.legend(loc='upper right')
plt.title("Figure 2: Trends in per-capita cigarette sales: California vs. synthetic California")
plt.axvline(INTERVENTION_TIME)
plt.text(x=INTERVENTION_TIME+0.2, y=30, s='Passage of Proposition 99')
plt.xlim([START_TIME, STOP_TIME-1])
plt.ylim([0, 140])
plt.grid()
plt.show()
fig.savefig("prop99_figure2", dpi=300)