因果推断系列17 - 合成控制法

1. 神奇的数学技巧

研究双重差分法(DID)时，我们有来自两个不同城市（Porto Alegre和Florianopolis）的多个客户的数据。数据跨越了两个不同的时间段：在Porto Alegre进行市场干预之前和之后以提高客户存款。为了估计处理效应，我们进行了一次回归，得到了DID估计量及其标准误差。

在这种情况下，我们有很多样本，数据是个体层面的。但是如果我们只有城市层面的聚合数据呢？例如，我们只有干预前后两个城市的平均存款水平。如下表所示：

city	before	after
FL	171.64	206.16
POA	46.01	87.06

我们仍然可以计算差异法估计量。

$(E[Y(1)|D=1]-E[Y(1)|D=0])-(E[Y(0)|D=1]-E[Y(0)|D=0])=(87.06-206.16)-(46.01-171.64)=6.53$

然而，这里的样本量只有4，这也是DID模型中的参数数量。在这种情况下，标准误差没有明确定义，那么我们该怎么办呢？另一个问题是，Florianopolis可能与Porto Alegre的相似性不如我们所希望的那样。例如，Florianopolis以其美丽的海滩和随和的人民而闻名，而Porto Alegre则以其烧烤和草原而更为著名。问题在于，您无法确定是否使用了适当的对照组。

为了解决这个问题，我们将使用所谓的“近年来政策评估文献中最重要的创新”，即合成控制。它基于一个简单而强大的想法。我们不需要找到未经处理的单个单位与受处理单位非常相似。相反，我们可以将多个未经处理的单位组合在一起，创建一个有效的合成控制。合成控制非常有效而直观，它的提出是在华盛顿邮报上发表了一篇文章，而不是在科学期刊上发表。

import warnings
warnings.filterwarnings('ignore')

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

%matplotlib inline

pd.set_option("display.max_columns", 6)
style.use("fivethirtyeight")

为了看到它的实际效果，思考一下对香烟税对其消费量的影响的问题。为了给出一些背景，这是一个长期以来在经济学领域一直存在争议的问题。争论的一方认为税收将增加香烟的成本，从而降低需求。另一方认为，由于香烟会导致成瘾，它们的价格变动对需求的影响不会很大。用经济学的术语来说，香烟的需求对价格是缺乏弹性的，增加税收只是以牺牲吸烟者为代价增加政府收入的一种方式。为了解决这个问题，我们将查看一些与此问题相关的美国数据。

1988年，加利福尼亚通过了一项著名的烟草税和健康保护法案，被称为99号提案。它的主要内容是"对加利福尼亚州内销售的烟草香烟征收每包25美分的州消费税，对雪茄和咀嚼烟等其他商业烟草产品的零售销售同样征收大致相等的消费税。对烟草销售的其他限制措施包括禁止在年轻人能够进入的公共区域放置香烟自动贩卖机，以及禁止单独销售散装香烟。该法案所产生的收入被用于各种环境和医疗保健项目，以及反烟广告。"

为了评估其效果，收集了来自多个州多年的香烟销售数据。在这个案例中，收集了39个州1970年至2000年的数据。其他州也有类似的烟草控制计划，它们在分析中被排除在外。数据样本如下，smoking.csv数据下载。

cigar = pd.read_csv("data/smoking.csv").drop(columns=["lnincome","beer", "age15to24"])

cigar.query("california").head()

	state	year	cigsale	retprice	california	after_treatment
62	3	1970	123.000000	38.799999	True	False
63	3	1971	121.000000	39.700001	True	False
64	3	1972	123.500000	39.900002	True	False
65	3	1973	124.400002	39.900002	True	False
66	3	1974	126.699997	41.900002	True	False

我们将state作为州索引，其中加利福尼亚州编号为3。我们的协变量是retprice，香烟零售价格，以及cigsale，人均销售的香烟包数。我们感兴趣的结果变量是cigsale。最后，我们有布尔类型的辅助变量来表示加利福尼亚州和干预后的时期。如果我们将加利福尼亚州和其他州的香烟销售随时间的变化绘制成图表，将会得到以下结果。

ax = plt.subplot(1, 1, 1)

(cigar
 .assign(california = np.where(cigar["california"], "California", "Other States"))
 .groupby(["year", "california"])
 ["cigsale"]
 .mean()
 .reset_index()
 .pivot("year", "california", "cigsale")
 .plot(ax=ax, figsize=(10,5)))

plt.vlines(x=1988, ymin=40, ymax=140, linestyle=":", lw=2, label="Proposition 99")
plt.ylabel("Cigarette Sales Trend")
plt.title("Gap in per-capita cigarette sales (in packs)")
plt.legend()

数据展示的这段时间内，显然加州人们购买的香烟比全国平均水平要少。此外，80年代后似乎有减少香烟消费的趋势。看起来在99号提案通过后，与其他州相比，加州的减少趋势加速了，但我们不能确定这一点。这只是通过观察数据图形猜测得出的结论。

为了回答99号提案是否对香烟消费产生了影响，我们将利用干预前期建立一个合成对照组。合并其他州的数据，构建一个非常接近加州趋势的虚拟州。然后，观察这个合成对照组在干预后的行为。

2. 时间变量

为了更正式一些，假设我们有$J+1$个单位。不失一般性，假设单位1是受到干预影响的单位。单位$j=2,...,J+1$是一组未接受干预的单位，我们将其称为“捐赠池”。还假设我们的数据跨越了T个时间段，其中有$T_0$个时间段在干预之前。对于每个单位j和每个时间t，我们观察到结果$Y_{jt}$。对于每个单位j和每个时间段t，定义$Y_{jt}^N$为没有干预的潜在结果，$Y_{jt}^I$为有干预的潜在结果。那么，对于接受处理的单位$j=1$在时间t上，对于$t>T_0$，效应被定义为：

$${\tau }_{1t}=Y_{jt}^I-Y_{jt}^N$$

由于单位$j=1$是接受处理的单位，$Y_{jt}^I$是事实的，但$Y_{jt}^N$不是。那么难点就在于如何估计$Y_{jt}^N$。请注意，处理效应是针对每个时期定义的，这意味着它可以随时间变化。它不需要是瞬时的，可以累积或消散。用一张图片来描述，估计处理效应的问题归结为估计如果单位$j=1$没有接受处理，其结果会发生什么。

为了估计$Y_{jt}^N$，我们记住，与未经处理的单位相比，捐赠者池中单位的组合可能更好地近似处理单位的特征。因此，合成控制被定义为控制池中单位的加权平均值。给定权重$W=(w_2,...,w_{J+1})$，$Y_{jt}^N$的合成控制估计为：

$${\hat{Y}}_{jt}^N=\sum _{j=2}^{J+1}{w}_j{Y}_{jt}$$

如果这些数学让你头疼，别担心，我们可以使它更直观。首先，合成控制可以看作是一种倒过来做回归的方法。我们知道，线性回归也是通过变量的加权平均值来进行预测的方法。现在，想象一下那些回归，比如DID示例中每个变量都是一个时间段的哑变量。在这种情况下，回归可以表示为以下矩阵乘法:

在合成控制案例中，我们没有大量的个体，但我们有大量的时间段。因此，我们要翻转输入矩阵。然后，个体变成了"变量"，我们将结果表示为个体的加权平均值，如以下矩阵乘法所示：

如果每个时间段有多个特征，可以将这些特征堆叠在一起。重要的是要使回归试图通过使用其他个体"预测"处理个体1。这样，可以以某种最佳方式选择权重，以实现我们想要的接近度。甚至可以以不同的比例缩放特征，以赋予它们不同的重要性。

因此，如果将合成控制视为线性回归，这意味着可以使用OLS估计其权重，对吗？是的！事实上，我们现在就可以尝试这样做。

3.合成控制vs线性回归

为了估计合成控制的处理效应，我们将尝试构建一个在干预期之前类似于受处理个体的“虚拟个体”。然后，我们将观察这个“虚拟个体”在干预后的行为。合成控制与它所模拟的单位之间的差异就是处理效应。

为了用线性回归实现这一点，我们将使用最小二乘法来找到权重。我们将最小化给定的筹集池中单位的加权平均值与干预前期受治疗单位之间的平方距离。

为了实现这一点，我们首先需要将单位（在我们的例子中是州）转换为列，时间转换为行。由于我们有两个特征，cigsale和retprice，我们将它们堆叠在一起，如上面的图片所示。我们将构建一个在干预期前期看起来非常像加利福尼亚州的合成控制，并观察它在干预期后期的行为。因此，选择仅包括干预期前期是很重要的。在这里，特征似乎具有相似的尺度，所以我们不需要对它们进行任何处理。如果特征具有不同的尺度，一个在千位数级别，另一个在小数位数级别，那么在最小化差异时，较大的特征将起到更重要的作用。为了避免这种情况，需要首先对它们进行缩放。

features = ["cigsale", "retprice"]

inverted = (cigar.query("~after_treatment") # filter pre-intervention period
            .pivot(index='state', columns="year")[features] # make one column per year and one row per state
            .T) # flip the table to have one column per state

inverted.head()

	state	1	2	3	...	37	38	39
	year
cigsale	1970	89.800003	100.300003	123.000000	...	114.500000	106.400002	132.199997
	1971	95.400002	104.099998	121.000000	...	111.500000	105.400002	131.699997
	1972	101.099998	103.900002	123.500000	...	117.500000	108.800003	140.000000
	1973	102.900002	108.000000	124.400002	...	116.599998	109.500000	141.199997
	1974	108.199997	109.699997	126.699997	...	119.900002	111.800003	145.800003

5 rows × 39 columns

现在，我们可以将变量Y定义为加利福尼亚州的状态，将X定义为其他州。

y = inverted[3].values # state of california
X = inverted.drop(columns=3).values  # other states

然后，我们进行回归分析。有一个截距相当于添加另一个状态，其中每一行都是1。你可以这样做，但这会更复杂，我会将其省略。回归分析将返回一组权重，使得处理单元与捐助池中的个体之间平方差最小。

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

array([-0.436, -1.038,  0.679,  0.078,  0.339,  1.213,  0.143,  0.555,
       -0.295,  0.052, -0.529,  1.235, -0.549,  0.437, -0.023, -0.266,
       -0.25 , -0.667, -0.106, -0.145,  0.109,  0.242, -0.328,  0.594,
        0.243, -0.171, -0.02 ,  0.14 , -0.811,  0.362,  0.519, -0.304,
        0.805, -0.318, -1.246,  0.773, -0.055, -0.032])

这些权重告诉我们如何构建合成对照组。我们将州1的结果乘以-0.436，州2的结果乘以-1.038，州4的结果乘以0.679，依此类推。我们可以通过合成对照组中各个状态的矩阵与权重的点积来实现这一点。

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

现在有了合成对照组，我们可以将其与加利福尼亚州的结果变量绘制在一起。

plt.figure(figsize=(10,6))
plt.plot(cigar.query("california")["year"], cigar.query("california")["cigsale"], label="California")
plt.plot(cigar.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();

OK…有些不对劲。注意上图，首先，在干预后，合成对照组的香烟销量比加利福尼亚州要高。这表明干预措施成功地降低了香烟需求。其次，注意到干预前期的完美拟合。合成对照组能够与加利福尼亚州完全匹配。这表明我们的合成对照模型可能对数据过拟合了。另一个迹象是干预后合成对照组的结果变量存在巨大的方差。注意它不遵循平滑的模式，而是上下波动。

思考一下为什么会发生这种情况，请记住捐赠池中有38个州。因此，线性回归有38个参数可以使用，以使预处理池尽可能接近处理组。即使T很大，N也很大的情况下，这给了我们的线性回归模型太大的灵活性。如果您熟悉正则化模型，可以使用岭回归或Lasso回归来解决这个问题。在这里，我们介绍另一种更传统的避免过拟合的方法。

4. 外推？

假设您有如下表格中的数据，并被要求使用控制单元的任何线性组合来构建一个合成对照组，以复制处理组个体。

unit	sales	price
control 1	8	8
control 2	8	4
control 3	4	5
treated	2	10

由于有3个单元和2个属性要匹配，因此存在多个精确解决方案，但一个好的解决方案是将第一个控制单元乘以2.25，将第二个单元乘以-2，然后相加。注意第二次乘法创建了一个销售额为-16、价格为-8的虚拟个体。这种乘法将控制单元2外推到一个数据区域，这个区域几乎没有多少意义，因为负价格和销售额几乎是不可能的。第一次乘法也是外推，因为它将第一个单元带到了销售额和价格为18的区域。这些数字远高于我们数据中的任何值，因此是在进行外推。

这就是当我们要求回归创建一个合成对照组时所做的。外推在技术上并没有错，但在实践中是危险的。假设我们从未见过的数据的行为与我们已有的数据相似。

为了更加安全，一种方法是限制我们的合成对照组只进行内插。为此，我们将限制权重为正数并总和为1。现在，合成对照组将是捐赠池中个体的凸组合。在进行内插时，我们将处理单元投影到由未处理单元定义的凸包中，就像下面的图片中所示。

请注意两点。首先，在这种情况下，插值无法完美匹配处理单元。这是因为处理单元是销售数量最少、价格最高的单元。凸组合只能精确复制介于对照单元之间的特征。另一个要注意的是，插值是稀疏的。我们将处理单元投影到凸包的一面上，而这一面只由几个单元定义。因此，插值将对许多单元分配零权重。

这是总体思路，现在让我们稍微正式化一下。合成控制仍然被定义为：

$${\hat{Y}}_{jt}^N=\sum _{j=2}^{J+1}{w}_j{Y}_{jt}$$

但现在，我们将使用权重$W=(w_2,...,w_{J+1})$来最小化

$$||X_1-X_{0}W||=\left(\sum _{h=1}^k{v}_h\right(X_{h1}-\sum _{j=2}^{J+1}{w}_j{X}_{hj}{)}^2{)}^{\frac{1}{2}}$$

并且$w_2,...,w_{J+1}$受到正值限制且总和为1。注意当最小化处理单元和合成控制之间的差异时，$v_h$反映了每个变量的重要性。不同的$v$会给出不同的最佳权重。选择$V$的一种方法是使每个变量的均值为0，方差为1。更复杂的方法是选择$V$，使帮助更好地预测$Y$的变量具有更高的重要性。由于我们希望保持代码简单，所以简单地给予每个变量相同的重要性。

要实现这一点，首先定义上述损失函数。

from typing import List
from operator import add
from toolz import reduce, partial

def loss_w(W, X, y) -> float:
    return np.sqrt(np.mean((y - X.dot(W))**2))

由于我们对每个特征使用相同的重要性，不需要担心$v$。

现在，为了获得最佳权重，我们将使用scipy的二次规划优化方法。我们将通过以下约束使权重总和为1：

lambda x: np.sum(x) - 1

此外，我们将将优化边界设置为介于0和1之间。

from scipy.optimize import fmin_slsqp

def get_w(X, y):
    
    w_start = [1/X.shape[1]]*X.shape[1]

    weights = fmin_slsqp(partial(loss_w, X=X, y=y),
                         np.array(w_start),
                         f_eqcons=lambda x: np.sum(x) - 1,
                         bounds=[(0.0, 1.0)]*len(w_start),
                         disp=False)
    return weights

完成这一步后，便获得了定义合成对照组的权重。

calif_weights = get_w(X, y)
print("Sum:", calif_weights.sum())
np.round(calif_weights, 4)

Sum: 1.000000000000236


array([0.    , 0.    , 0.    , 0.0852, 0.    , 0.    , 0.    , 0.    ,
       0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    ,
       0.    , 0.    , 0.    , 0.113 , 0.1051, 0.4566, 0.    , 0.    ,
       0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    ,
       0.2401, 0.    , 0.    , 0.    , 0.    , 0.    ])

通过这些权重，我们将州索引1、2和3乘以零，州索引4乘以0.0852，依此类推。注意权重是稀疏的，正如我们预测的那样。而且，所有权重的总和为1，并且介于0和1之间，满足我们的凸组合约束条件。

现在，为了得到合成对照组，我们可以将这些权重与状态相乘，就像我们之前用回归权重做的那样。

calif_synth = cigar.query("~california").pivot(index='year', columns="state")["cigsale"].values.dot(calif_weights)

如果我们现在绘制合成对照组的结果，我们会得到一个更平滑的趋势。还要注意，在干预前的时期，合成对照组不再完全重现处理组的情况。这是一个好迹象，表明没有过拟合。

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

有了合成对照组的数据，我们可以估计处理效应，即处理组和合成对照组结果之间的差距。

$${\tau }_{1t}=Y_{jt}^I-Y_{jt}^N$$

在我们的香烟案例中，随着时间的推移，效应越来越大。

plt.figure(figsize=(10,6))
plt.plot(cigar.query("california")["year"], cigar.query("california")["cigsale"] - calif_synth,
         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包。这非常好，但你可能会问自己：这个结果是否具有统计学意义？

5. 推断

由于我们的样本规模非常小（39个观测值），如果要确定结果是否具有统计学意义而不仅仅是运气时，需要更加精巧。此时使用费舍尔精确检验的思想。它的思想非常简单：对处理组和对照组进行穷举的置换。由于我们只有一个处理单元，这意味着对于每个单元，假装它是处理组，而其他单元则是对照组。

iteration	1	2	…	39
1	treated	0	0	0
2	0	treated	0	0
…	0	0	0	0
39	0	0	0	treated

最终，我们将为每个州得到一个合成对照和效果估计值。所以这样做的目的是假设治疗实际上发生在另一个州而不是加利福尼亚，并观察对于这个未发生治疗的估计效果会是什么。然后，我们比较加利福尼亚州的治疗效果是否显著大于其他虚拟治疗效果。这个想法是对于实际上没有接受治疗的州，一旦我们假设它们接受了治疗，我们将无法找到任何显著的治疗效果。

为了实现这一点，我编写了这个函数，它以一个州作为输入，估计该州的合成对照。该函数返回一个数据框，其中包含一个列用于州，一个列用于年份，一个列用于cigsale结果以及该州的合成结果。

def synthetic_control(state: int, data: pd.DataFrame) -> 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 = get_w(X, y)
    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))

以下是将该函数应用于编号为1的州时的结果。

synthetic_control(1, cigar).head()

	state	year	cigsale	after_treatment	synthetic
0	1	1970	89.800003	False	95.029419
1	1	1971	95.400002	False	99.118199
2	1	1972	101.099998	False	101.881329
3	1	1973	102.900002	False	103.938655
4	1	1974	108.199997	False	107.038474

为了得到所有州的结果，我们将计算并行化到4个进程中。如果您的计算机核心更多或更少，可以使用不同的数值。此代码将返回一个类似上述数据框的数据框列表。

from joblib import Parallel, delayed

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

parallel_fn = delayed(partial(synthetic_control, data=cigar))

synthetic_states = Parallel(n_jobs=4)(parallel_fn(state) for state in control_pool)

synthetic_states[0].head()

	state	year	cigsale	after_treatment	synthetic
0	1	1970	89.800003	False	95.029419
1	1	1971	95.400002	False	99.118199
2	1	1972	101.099998	False	101.881329
3	1	1973	102.900002	False	103.938655
4	1	1974	108.199997	False	107.038474

通过为所有州获得合成对照，以便估计所有州的合成对照和真实州之间的差距。对于加利福尼亚州，这是处理效应。对于其他州，这类似于安慰剂效应，我们估计了处理实际上没有发生时的合成对照处理效应。如果我们将所有安慰剂效应以及加利福尼亚州的处理效应绘制在一起，可得到以下图像。

plt.figure(figsize=(12,7))
for state in synthetic_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,
        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()

这张图中有两个显著的方面。首先，干预后的方差高于干预前的方差。这是预期的，因为合成控制被设计为最小化干预前期的差异。另一个有趣的方面是，即使在干预前期，也有一些单位无法很好地拟合。这也是可以预期的。例如，如果某些州的香烟消费非常高，那么其他州的凸组合也无法与其匹配。

由于这些单位的拟合效果非常差，将它们从分析中移除是一个好主意。一种客观的做法是设置一个干预前期误差的阈值，

$MSE=\frac{1}{N}\sum (Y_t-{\hat{Y}}_t^{Synth}{)}^2$

并移除那些具有较高误差的单元。如果我们按照这种方式进行，并绘制相同的图形，会得到以下结果。

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 synthetic_states:
    
    # remove units with mean error above 80.
    if pre_treatment_error(state) < 80:
        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,
        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();

去除噪声后，我们可以看到加利福尼亚州的效果值是多么极端。这张图告诉我们，如果假设处理发生在任何其他州，几乎不太可能得到像加利福尼亚州那样极端的效果。

这张图片本身就是一种推断的形式，但我们还可以从这些结果中推导出p值。我们只需要看看我们得到的效果有多少次低于加利福尼亚州的效果。

calif_number = 3

effects = [state.query("year==2000").iloc[0]["cigsale"] - state.query("year==2000").iloc[0]["synthetic"]
           for state in synthetic_states
           if pre_treatment_error(state) < 80] # filter out noise

calif_effect = cigar.query("california & year==2000").iloc[0]["cigsale"] - calif_synth[-1] 

print("California Treatment Effect for the Year 2000:", calif_effect)
np.array(effects)

California Treatment Effect for the Year 2000: -24.830159744091695


array([  5.79715892,   0.8945899 , -24.83015974,  -7.16628128,
       -10.92204859,  37.11640554, -15.06971737,  -0.49805125,
       -18.45795064,  21.13366433,  12.57782758,  -1.47547827,
        10.49627385, -11.67012334,   4.29850825,   8.04811409,
        14.02322406,   8.25002726,   0.32576355,  -8.40826864,
        -2.12402711,  -7.42865049,   2.96157493,  24.10478129,
         4.25211769, -17.75844593,   7.93334018,   2.81640125,
        12.64955948, -17.47677513, -25.16040936, -12.26469128,
        24.69067348,  10.3629958 ,  -8.5988038 ])

如果我们想要测试加利福尼亚州效果低于零的单侧假设，我们可以将p值估计为加利福尼亚州效果大于所有估计效果的比例。

P V = 1 N ∑ 1 { τ ^ C a l i f > τ ^ j } PV=\frac{1}{N}\sum \mathcal{1}\{\hat{\tau}_{Calif} > \hat{\tau}_j\} PV=N1​∑1{τ^Calif​>τ^j​}

结果发现，2000年加利福尼亚州的处理效应为-24.8，意味着干预减少了近25包香烟的消费量。在我们估计的其他34个安慰剂效应中，只有一个比我们在加利福尼亚州发现的效果更高。因此，p值将为1/35。

np.mean(np.array(effects) < calif_effect)

0.02857142857142857

最后，我们可以展示效果分布，以了解加利福尼亚州效果值的极端程度。

_, bins, _ = plt.hist(effects, bins=20, color="C5", alpha=0.5);
plt.hist([calif_effect], bins=bins, color="C0", label="California")
plt.ylabel("Frquency")
plt.title("Distribution of Effects")
plt.legend();

小结

通过本文了解到，如果只有关于城市或州等实体的聚合级别数据，DID无法进行推断。此外，它还有一些其他限制，因为它必须定义一个控制单元，而一个单一的控制单元可能不是处理单元反事实很好的代表。

为了纠正这一点，我们学会了构建一个合成控制单元，将多个控制单元组合起来使它们类似于待处理单元。通过这个合成控制单元，能够看到在没有处理的情况下，我们的待处理单元会发生什么。

最后，我们看到了如何使用费舍尔精确检验在合成控制中进行推断。即假设非处理单元实际上是处理单元，并计算它们的效果。这些是安慰剂效应：即使没有处理，我们也会观察到的效果。我们使用这些效应来判断我们估计的处理效应是否具有统计学意义。