第11章 不依从与工具变量
企业向现有客户群体提供产品或服务的情况并不少见。例如，零售商可以推出一项基于订阅的计划，客户能享受免费配送服务。流媒体公司可以额外收费提供无广告版的服务。或者，银行可以为消费超过一定门槛的客户提供一张附带诸多优惠的高端信用卡。
在所有这些例子中，客户必须选择加入额外服务，这使得推断该服务的影响颇具挑战性。由于是否参与的选择权在客户手中，这种选择往往会干扰对服务影响的评估；毕竟，选择加入的客户和不选择加入的客户很可能有着不同的 \(Y_0\)（注：此处 \(Y_0\) 为专业符号，可结合具体研究场景理解，一般指未接受处理时的结果变量）。即便企业将服务或产品的可获得性进行随机分配，也无法强迫客户接受。这被称为 “不依从”，即并非所有被分配到处理组的人都会接受处理。本章将引导你思考这个问题，以及当你想要设计一个存在不依从情况的实验时需要考虑哪些方面。
不依从
不依从（Noncompliance）这一概念源于制药科学（不过，一些用于应对它的方法工具来自经济学）。想象一下，你正在开展一项实验，以测试一种新药对某种疾病的疗效。每个实验对象都会被分配到一个处理组：用药组或者安慰剂组。但这些实验对象都是普通人，有时会忘记服药。结果就是，并非所有被分配到用药组的人都会真正服药。此外，病情严重的人可能会发现自己被分配到了安慰剂组，然后想方设法获取药物治疗。也就是说，如果你把 “处理分配” 和 “处理接受” 区分开来，最终会得到四个组：
依从者：接受分配给他们的处理的人。
总是接受者：无论分配如何，总是接受处理的人。
从不接受者：无论分配如何，从不接受处理的人。
违背者：接受与分配相反的处理的人。
这里的难题在于，你不知道谁属于每一个组。
你也可以在有向无环图（DAG）中表示不依从情况（见图 11-1），其中Z是处理分配（在这种情况下是随机的），T是处理，Y是结果，U是混淆处理选择和结果的隐藏因素。Z就是所谓的工具变量：一个满足以下两个条件的变量：（1）以无混淆的方式影响处理；（2）不会影响结果，除非通过处理。
由于依从组和处理分配确定性地导致处理的接受，你也可以把 U 看作是导致依从组出现的未知因素。从图 11-1 的有向无环图（DAG）中可以看到，在没有进一步假设的情况下，由于存在通过 U 的开放后门路径，你无法识别处理对 Y 的效应。正如你很快会看到的，对该效应的识别将涉及对 Z 的巧妙运用。
为了让事情更具体，我们来看一个典型的行业场景：一家银行想了解向其客户提供高端信用卡的影响。由于这项高端服务成本很高，银行会收取少量费用，但这不足以覆盖其所有成本。不过，如果那些高端客户的购买额 —— 即使用该卡的总消费额 —— 至少增加 500 美元，那么提供这项服务就是值得的。因此，银行想知道这张高端信用卡能在多大程度上提高客户的购买额。
银行设法开展了一项实验，将高端信用卡的可获得性（prime_eligible）随机分配给
10,000 名客户，每位客户有 50% 的概率获得资格，50%
的概率处于对照组。当然，银行无法强迫客户选择这张卡，这就使得这项实验存在不依从的情况。
如果你把这些变量对应到图 11-1 中，购买额就是 Y，高端信用卡的可获得性就是 Z，拥有高端信用卡就是 T。所有这些信息都存储在下面的数据框中。银行还掌握客户的年龄、收入和信用评分等信息，但现在我们先不考虑这些变量。此外，我还添加了关于高端信用卡对购买额（PV）的真实效应 \(\tau\)（tau），以及客户所属组别的信息。要记住，在现实中，依从类别和 \(\tau\) 是你无法获取的。我在这里使用它们只是为了让一些例子更容易理解：
from toolz.curried import *
import pandas as pd
import numpy as np
from scipy.special import expit
import statsmodels.formula.api as smf
import seaborn as sns
from matplotlib import pyplot as plt
import matplotlib
from cycler import cycler
color=['0.0', '0.4', '0.8']
default_cycler = (cycler(color=color))
linestyle=['-', '--', ':', '-.']
marker=['o', 'v', 'd', 'p']
plt.rc('axes', prop_cycle=default_cycler)
from graphviz import Digraph
gr = Digraph(format="png", graph_attr={"rankdir":"LR"})
gr.edge("U", "T")
gr.edge("U", "Y")
gr.edge("Z", "T")
gr.edge("T", "Y")import pandas as pd
import numpy as np
df = pd.read_csv("./data/prime_card.csv")
df.head()## age income credit_score ... pv tau categ
## 0 37.7 9687.0 822.0 ... 4913.79 700.0 complier
## 1 46.0 13731.0 190.0 ... 5637.66 200.0 never-taker
## 2 43.1 2839.0 214.0 ... 2410.45 700.0 complier
## 3 36.0 1206.0 318.0 ... 1363.06 700.0 complier
## 4 39.7 4095.0 430.0 ... 2189.80 700.0 complier
##
## [5 rows x 8 columns]拓展潜在结果
为了更精准地研究不依从问题并推进效应识别，你必须拓展潜在结果的表示方法。由于Z会影响T，现在你可以定义一个潜在处理\(T_z\)。此外，潜在结果有了关于工具变量Z的新反事实情况\(Y_{z,t}\)。
在高端信用卡的例子中，Z是随机分配的，这意味着Z对Y的效应 —— 也被称为意向治疗效应（ITTE）—— 是相当容易识别的：
\(ITTE = E[Y \vert Z = 1] - E[Y \vert Z = 0] = E[Y_{1,t} - Y_{0,t}]\)
你可以用一个简单的线性回归来估计它：
m = smf.ols("pv~prime_elegible", data=df).fit()
print(m.summary().tables[1])## ==================================================================================
## coef std err t P>|t| [0.025 0.975]
## ----------------------------------------------------------------------------------
## Intercept 2498.3618 24.327 102.701 0.000 2450.677 2546.047
## prime_elegible 321.3880 34.321 9.364 0.000 254.113 388.663
## ==================================================================================意向治疗效应（ITTE）本身就是一个很有价值的指标，因为它衡量了分配一项处理的影响，比如在这个场景中提供高端信用卡。对银行来说，这个数值表明，通过将高端信用卡纳入其产品组合，每位客户预期能增加的购买额（PV）。然而，必须注意的是，意向治疗效应和处理效应并不相同。银行的主要目标是确定高端信用卡的收益是否超过其成本。因此，银行需要识别选择该信用卡的处理效应，而不是仅仅依赖意向治疗效应。
在这个特定的例子中，银行完全控制着哪些客户可以获得高端信用卡。因此，存在单向不依从，因为没有资格获得高端信用卡的客户无法得到它，但有资格获得该卡的客户仍然可以选择不使用它。这就使得总是接受者（always takers）处于依从状态，而违背者（defiers）成为从不接受者（never takers），从而将依从组的数量从四个减少到两个。
既然你已经了解了这个场景，那我们来思考如何识别高端信用卡的效应。一个很直接的想法是用意向治疗效应（ITTE）作为信用卡效应的替代指标。或许它们终究并没有那么大的差别。那么，意向治疗效应到底是什么呢？
由于处理分配是随机的，意向治疗效应可以通过比较被分配到处理组的人和被分配到对照组的人来获得。但你很快会发现，这种比较会让你得到一个向零偏倚的处理效应估计值（见图 11-2）。这是因为，那些被分配到处理组的人中，有一部分实际上接受的是对照处理，这就减小了两组之间可感知到的差异。
为了证明这一点，你可以利用我添加到数据集中的那些\(\tau\)值。平均处理效应要比意向治疗效应（ITTE）大得多：
df["tau"].mean()## 413.45好的，看来那条路走不通了。不过，直接对处理组和对照组进行简单的均值比较，即\(E[Y \vert T = 1] - E[Y \vert T = 0]\)，会怎么样呢？也许随机分配会让这个结果成为你所关注的效应估计的一个良好替代指标。嗯…… 那我们来估计一下看看：
m = smf.ols("pv~prime_card", data=df).fit()
m.summary().tables[1]| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 2534.4947 | 19.239 | 131.740 | 0.000 | 2496.783 | 2572.206 |
| prime_card | 588.1388 | 41.676 | 14.112 | 0.000 | 506.446 | 669.831 |
现在，测量到的效应比真实效应大得多（见图 11-3）。原因在于，在这个例子中，偏差是向上的，\(E[Y_0 \vert T = 1] > E[Y_0 \vert T = 0]\)，这意味着选择高端信用卡的客户，无论是否使用该信用卡，消费都更多。换句话说，从不接受者（never takers）的\(Y_0\)比依从者（compliers）的\(Y_0\)更低，这就人为地降低了未处理组的平均结果。
工具变量识别假设
让我们把我之前给你们展示的有向无环图（DAG）再在这里呈现一遍，以便更好地阅读。你们会看到，识别所需的最初几个假设已经在那个有向无环图中明确说明了：
from graphviz import Digraph
gr = Digraph(format="png", graph_attr={"rankdir":"LR"})
gr.edge("U", "T")
gr.edge("U", "Y")
gr.edge("Z", "T")
gr.edge("T", "Y")
gr## <graphviz.graphs.Digraph object at 0x000001B5FF366340>
它们如下：
你需要的第一个假设是独立性；即Z和T之间不存在未测量的混杂因素，\(T_z \perp Z \vert X\)，且Z和Y之间也不存在未测量的混杂因素，\(Y(Z, T_z) \perp Z \vert X\)。这个假设表明工具变量就如同是随机分配的。这个假设无法检验，但可以通过实验设计使其更具合理性。在你的例子中，你大概可以说这个假设是成立的，因为银行对高端信用卡的可获得性进行了随机分配。
第二个假设是排除限制，即\(Y_{z,t} = Y_t\)，这意味着不存在从Z到Y且不经过处理T的路径。换句话说，它表明工具变量只能通过处理对结果产生影响。这个假设更 tricky（棘手）。即使Z是随机分配的，它也可能通过其他渠道影响结果。例如，假设客户弄清楚了自己被分配到了哪一组，对照组的客户对银行非常生气，于是决定关闭账户。在这种情况下，随机分配通过非处理的渠道影响了结果。
第三个假设是相关性，即\(E[T_1 - T_0] \neq 0\)，也就是存在从Z到T的箭头。这个假设表明工具变量必须对处理有影响。幸运的是，这个假设是可检验的，因为你可以估计工具变量对处理的效应。
第四个也是最后一个假设在有向无环图（DAG）中没有说明。它主要是对因果模型施加的一种函数形式假设：单调性，即\(T_{i1} \geq T_{i0}\)（或者相反）。这可能看起来令人困惑，但它只是表明工具变量仅在一个方向上改变处理情况。它要么增加了所有获得工具变量的人接受处理的概率，这相当于假设不存在违背者（defiers）；要么降低了该概率，这相当于假设不存在依从者（compliers）。在你的例子中，这也是一个合理的假设，因为对照组的客户无法强行获得高端信用卡。因此，违背者被归到了从不接受者（never takers）中。
现在，让我们看看如何利用这些假设来进行识别。这里的目标是从意向治疗效应（ITTE）入手，看看我们能否得到类似平均处理效应的结果。首先，让我们把结果展开为潜在结果。回想一下，你可以用处理作为一个开关来做这件事 ——\(Y = Y_1 T + Y_0 (1 - T)\)：
\(\begin{align*} E[Y \vert Z = 1] - E[Y \vert Z = 0] &= E[Y_{1,1} T_1 + Y_{1,0} (1 - T_1) \vert Z = 1] \\ &\quad - E[Y_{0,1} T_0 + Y_{0,0} (1 - T_0) \vert Z = 0] \end{align*}\)
现在，由于排除限制，你可以去掉\(Y_{z,t}\)的工具变量下标：
\(E[Y_1 T_1 + Y_0 (1 - T_1) \vert Z = 1] - E[Y_1 T_0 + Y_0 (1 - T_0) \vert Z = 0]\)
接下来，利用独立性假设，你可以合并这两个期望：
\(E[Y_1 T_1 + Y_0 (1 - T_1) - Y_1 T_0 - Y_0 (1 - T_0)]\)
你可以将其简化为：
\(E[(Y_1 - Y_0)(T_1 - T_0)].\)
接下来，让我们利用单调性假设，将这个期望展开为可能的情况，即\(T_1 > T_0\)和\(T_1 = T_0\)：
\(\begin{align*} & E[(Y_1 - Y_0)(T_1 - T_0) \vert T_1 > T_0] * P(T_1 > T_0) \\ + & E[(Y_1 - Y_0)(T_1 - T_0) \vert T_1 = T_0] * P(T_1 = T_0) \end{align*}\)
而且，由于如果\(T_1 = T_0\)，\(T_1 - T_0\)就为 0，所以只剩下第一项。因为Z是二元的，\(T_1 - T_0 = 1\)，并且：
\(E[Y \vert Z = 1] - E[Y \vert Z = 0] = E[Y_1 - Y_0 \vert T_1 > T_0] * P(T_1 > T_0).\)
现在是停下来看看你取得了什么成果的好时机。首先，要注意到\(T_1 > T_0\)的是依从者，也就是工具变量将处理从 0 转变为 1 的那部分人群。这意味着最后这个等式告诉你，工具变量对结果的效应等于依从者的处理效应乘以依从率。这就解释了为什么意向治疗效应（ITTE）是对该效应的一个向零偏倚的估计值：你是用一个介于 0 和 1 之间的比率去乘以它。如果你们能估计出\(P(T_1 > T_0)\)，那么你们就能修正之前的估计量……
不过，等一下！由于工具变量是随机分配的，你可以估计它对处理的影响，即\(E[T_1 - T_0]\)。而且，由于对于依从者来说\(T_1 - T_0 = 1\)，其他情况下为 0（由于单调性假设），这个影响就是依从率：
\(E[T_1 - T_0] = P(T_1 > T_0)\)
把这些都综合起来，这意味着你可以通过用依从率（也就是工具变量对处理的影响）去放大工具变量对结果的影响，从而识别出依从者的平均处理效应：
\(E[Y_1 - Y_0 \vert T_1 > T_0] = \frac{E[Y \vert Z = 1] - E[Y \vert Z = 0]}{E[T \vert Z = 1] - E[T \vert Z = 0]}\)
这就是你在存在不依从情况时，如何利用工具变量来识别效应的方法。好消息是，你能够识别出那种效应。坏消息是，它不是平均处理效应（ATE），而只是依从者的效应，这通常被称为局部平均处理效应（LATE）。遗憾的是，你无法识别平均处理效应。但这可能不是什么问题。在你的信用卡例子中，局部平均处理效应会是对那些在高端信用卡可供选择时选择了它的客户的效应。现在，银行想知道额外购买额（PV）方面的效应是否能弥补高端信用卡的成本，而这两者都只发生在选择了高端信用卡的客户身上。在这种情况下，了解局部平均处理效应就足够了。银行并不关心那些永远不会选择高端信用卡的客户的效应。
在学习了关于工具变量识别的理论之后，现在是时候看看如何在实践中应用它了。
第一阶段
工具变量分析的第一步是进行所谓的 “第一阶段回归”，也就是将处理变量对工具变量进行回归。在这个步骤中，你可以检验相关性假设 —— 如果与工具变量相关的参数估计值较大且具有统计显著性，你就有充分的理由相信该假设成立：
first_stage = smf.ols("prime_card ~ prime_elegible", data=df).fit()
print(first_stage.summary().tables[1])## ==================================================================================
## coef std err t P>|t| [0.025 0.975]
## ----------------------------------------------------------------------------------
## Intercept -5.041e-17 0.005 -1.02e-14 1.000 -0.010 0.010
## prime_elegible 0.4242 0.007 60.536 0.000 0.410 0.438
## ==================================================================================在这个例子中，依从率估计约为 42%，且具有统计显著性（95% 置信区间为 [0.410, 0.438]）。由于我已经包含了每位客户所属的真实组别，你甚至可以检验这是否确实是实际的依从率：
df.groupby("categ").size()/len(df)## categ
## complier 0.4269
## never-taker 0.5731
## dtype: float64简化式
第二步被称为简化式。在这个阶段，你将结果对工具变量进行回归，以获得意向治疗效应：
red_form = smf.ols("pv ~ prime_elegible", data=df).fit()
print(red_form.summary().tables[1])## ==================================================================================
## coef std err t P>|t| [0.025 0.975]
## ----------------------------------------------------------------------------------
## Intercept 2498.3618 24.327 102.701 0.000 2450.677 2546.047
## prime_elegible 321.3880 34.321 9.364 0.000 254.113 388.663
## ==================================================================================一旦你完成了第一阶段和简化式的回归，你可以用第一阶段的参数估计值除以简化式的参数估计值，从而得到局部平均处理效应的估计值：
late = (red_form.params["prime_elegible"] / first_stage.params["prime_elegible"])
late## 757.6973795343862如你所见，这个效应是意向治疗效应（ITTE）的两倍多。这是符合预期的，因为依从率低于 50%。它也比平均处理效应（ATE）大，但这是因为依从者的效应高于平均水平。事实上，如果你计算依从者的效应，你会发现你的局部平均处理效应（LATE）估计值与它相当接近：
df.groupby("categ")["tau"].mean()## categ
## complier 700.0
## never-taker 200.0
## Name: tau, dtype: float64不过，仍然存在一些差异。如果你不把那个点估计值放在置信区间内，就很难判断这是否正确。你可以用自助法（bootstrap）来做这件事，但我认为有必要研究一下工具变量（IV）估计值的标准误的实际公式。要做到这一点，你必须学习另一种估计局部平均处理效应（LATE）的方法。
两阶段最小二乘法
如果你仔细看有向无环图（DAG）中处理部分，会发现它由两个部分导致：首先，有一个随机部分，也就是随机化的工具变量。其次，有一个U部分，混杂偏倚就来自这个部分：
gr = Digraph(format="png", graph_attr={"rankdir":"LR"})
gr.edge("U", "T")
gr.edge("U", "Y")
gr.edge("Z", "T")
gr## <graphviz.graphs.Digraph object at 0x000001B5922CFF10>
回想一下，第一阶段是将处理变量对工具变量进行回归，这本质上意味着你在估计路径\(Z \to T\)。但还有更多内容。由于这就是第一阶段所估计的，你可以把它的预测值\(\widehat{T}\)看作是处理变量的一个无偏版本。这进而意味着，如果你将结果对这些预测值进行回归，你会得到与之前相同的工具变量（IV）估计值：
gr = Digraph(format="png", graph_attr={"rankdir":"LR"})
gr.edge("U", "T")
gr.edge("U", "Y")
gr.edge("Z", "T")
gr## <graphviz.graphs.Digraph object at 0x000001B5923613D0>iv_regr = smf.ols(
"pv ~ prime_card",
data=df.assign(prime_card=first_stage.fittedvalues)).fit()
iv_regr.summary().tables[1]| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 2498.3618 | 24.327 | 102.701 | 0.000 | 2450.677 | 2546.047 |
| prime_card | 757.6974 | 80.914 | 9.364 | 0.000 | 599.091 | 916.304 |
这种方法被称为两阶段最小二乘法（2SLS）。好的，但它为什么有用呢？首先，因为它能让你正确计算标准误；其次，因为它使得添加更多工具变量和协变量变得像在回归模型中添加变量一样容易。现在，让我们依次来谈谈这些方面。
标准误
第二阶段预测的残差可以定义如下：
\(\hat{e}_{IV} = Y - \hat{\beta}_{IV}T\)
注意，这与你用第二阶段的.resid方法得到的残差不同，因为后者会是\(Y -
\hat{\beta}_{IV}\hat{T}\)。你需要的残差使用的是处理变量的原始版本，而不是预测版本。
利用该残差，你可以计算工具变量（IV）估计值的标准误：
\(SE\left( \widehat{\beta}_{IV} \right) = \frac{\sigma\left( \hat{e}_{IV} \right)}{\widehat{\beta}_{z,1st} \sigma\left( Z \right) \sqrt{n}}\)
其中，\(\sigma(.)\)代表标准差函数，\(\widehat{\beta}_{z,1st}\)是估计的依从率，它是你从第一阶段得到的：
Z = df["prime_elegible"]
T = df["prime_card"]
n = len(df)
# not the same as iv_regr.resid!
e_iv = df["pv"] - iv_regr.predict(df)
compliance = np.cov(T, Z)[0, 1]/Z.var()
se = np.std(e_iv)/(compliance*np.std(Z)*np.sqrt(n))
print("SE IV:", se)## SE IV: 80.52861026141971print("95% CI:", [late - 2*se, late + 2*se])## 95% CI: [596.6401590115468, 918.7546000572256]为了再次检查你的结果，你可以使用 Python
包linearmodels中的2SLS模块。借助它，你可以将第一阶段（如[T~Z]）包含在模型的公式中，并拟合一个工具变量（IV）模型。如你所见，使用这个包不仅能得到与你之前得到的相同的局部平均处理效应（LATE）估计值，还能得到与你刚刚计算的相同的标准误：
from linearmodels import IV2SLS
formula = 'pv ~ 1 + [prime_card ~ prime_elegible]'
iv_model = IV2SLS.from_formula(formula, df).fit(cov_type="unadjusted")
iv_model.summary.tables[1]| Parameter | Std. Err. | T-stat | P-value | Lower CI | Upper CI | |
|---|---|---|---|---|---|---|
| Intercept | 2498.4 | 24.211 | 103.19 | 0.0000 | 2450.9 | 2545.8 |
| prime_card | 757.70 | 80.529 | 9.4090 | 0.0000 | 599.86 | 915.53 |
无论你使用哪种方法，都可以看到这是一个相当大的置信区间，即便它确实包含了依从者的真实平均处理效应（ATE），也就是 700。更重要的是，我认为标准误公式可以让我们对不依从实验的挑战有所了解。首先，注意分母中的\(\sigma(Z)\)。由于Z是一个二元变量，\(\sigma(Z)\)的最大值是 0.5。这与具有二元处理的普通最小二乘法（OLS）并没有太大不同。（回想一下，在那种情况下，标准误是\(\sigma(\hat{e}) / (\sigma(T)\sqrt{n})\)）。它只是表明，你可以通过以 50%-50% 的方式对处理进行随机化来最大化检验的效力。
但现在分母中多了一个项：依从率\(\widehat{\beta}_{z,1st}\)。不出所料，如果依从率是 100%，那么\(Z = T\)，\(\widehat{\beta}_{z,1st} = 1\)，你就回到了普通最小二乘法（OLS）的标准误。但在存在不依从的情况下，标准误会增大，因为\(\widehat{\beta}_{z,1st} < 1\)。例如，当依从率为 50% 时，工具变量（IV）的标准误会是 OLS 标准误的两倍。因此，对于依从率为 50% 的实验，所需的样本量是依从率为 100% 时所需样本量的 4 倍。
工具变量的偏倚
事实证明，工具变量（IV）估计量是一致的，但不是无偏的。也就是说，\(E[\beta_{IV}] \neq \beta\)。这主要是由于第一阶段存在抽样误差。因为你没有无限的数据，处理变量的拟合值\(\widehat{T}\)会是Z和U的函数，这意味着并非所有的偏倚都会消失。当你收集更多数据时，\(\widehat{T}\)对U的依赖性会越来越小。这就是该估计量具有一致性的原因，即\(\text{plim}_{n \to +\infty} \beta_{IV} = \beta\)。
下面的图表比较了在假设不同估计依从率的情况下，你的局部平均处理效应（LATE）参数估计值的置信区间大小（第一张图）。它还展示了，考虑到多种依从率，对于存在不依从情况的检验，你还需要多少额外的样本：
se_formula_iv = lambda compliance: np.std(e_iv)/(compliance*np.std(Z)*np.sqrt(n))
x = np.linspace(0.01, 1, 50)
effect = iv_regr.params["prime_card"]
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
ax1.plot(x, effect-se_formula_iv(x)*2, label="SE($\\beta_{IV}$) x2", ls=":", color="0")
ax1.plot(x, effect+se_formula_iv(x)*2, ls=":", color="0")
ax1.hlines(effect, 0, 1, ls="-.", label="LATE")
ax1.hlines(0, 0, 1)
ax1.set_xlabel("Compliance")
ax1.set_ylim(-(effect+100), (effect+100)*2)## (-857.6973795343919, 1715.3947590687837)ax1.legend(loc="lower right")
ax1.set_title("95% CI around LATE");
x = np.linspace(0.2, 1, 50)
ax2.plot(x, 1/(x**2))
ax2.hlines(10, 0.2, 1, ls=":", label="10")
ax2.hlines(4, 0.2, 1, ls="-.", label="4")
ax2.set_xlabel("Compliance")
ax2.set_ylabel("$N_{IV}$/N")
ax2.set_title("Required Sample Size Ratio")
ax2.legend()
50% 的依从率仍然算高的。在大多数应用场景中，只有一小部分客户会选择优质服务或产品，这使得估计局部平均处理效应（LATE）更加困难。例如，如果依从率低至 30%，你需要的样本量将是不存在依从问题时所需样本量的 10 倍。收集如此大的样本量，即便不是完全不可能，往往也是不切实际的。但如果你遇到这样的问题，也并非毫无办法。仍然有一些技巧可以用来降低工具变量（IV）的标准误。要做到这一点，你需要在分析中纳入额外的协变量。
额外的控制变量和工具变量
还记得高端信用卡数据除了处理变量、工具变量和结果变量之外，还有三个额外的协变量吗？它们是客户的收入、年龄和信用评分。现在，假设描述它们与T、Z和Y之间关系的因果图如下：
gr = Digraph(format="png", graph_attr={"rankdir":"LR"})
gr.edge("U", "T")
gr.edge("U", "Y")
gr.edge("Z", "T")
gr.edge("T", "Y")
gr.edge("Income", "Y")
gr.edge("Age", "T")
gr.edge("Age", "Y")
gr.edge("Score", "T")
gr## <graphviz.graphs.Digraph object at 0x000001B5973F7D90>
也就是说，收入能高度预测结果，但不能预测依从性；信用评分能预测依从性，但不能预测结果；而年龄能预测这两者（是一个混杂因素）。如果你能巧妙运用这些变量，通过使用所有这些变量，你可以降低标准误。
首先，来看信用评分。信用评分会导致依从性，但不会导致结果。这意味着它可以被当作一个额外的工具变量。从有向无环图（DAG）中，你可以看到它和Z一样，满足前三个工具变量（IV）假设。你只需要假设正性（positivity）。在你的两阶段最小二乘法（2SLS）模型中加入这个作为额外工具变量的变量，会显著降低局部平均处理效应（LATE）参数的标准误：
formula = 'pv ~ 1 + [prime_card ~ prime_elegible + credit_score]'
iv_model = IV2SLS.from_formula(formula, df).fit()
iv_model.summary.tables[1]| Parameter | Std. Err. | T-stat | P-value | Lower CI | Upper CI | |
|---|---|---|---|---|---|---|
| Intercept | 2519.4 | 21.168 | 119.02 | 0.0000 | 2477.9 | 2560.9 |
| prime_card | 659.04 | 58.089 | 11.345 | 0.0000 | 545.19 | 772.90 |
现在，你在这里必须要小心。如果不是把它当作工具变量，而是对它进行控制，也把它加入到第二阶段中，误差会增大。但你已经知道这一点了，从第 4 章中，你了解到对那些导致处理但不导致结果的变量进行控制，会增加方差。
这里更大的问题是，除非你了解工具变量的分配机制（就像对Z那样），否则很难知道排除限制是否成立。例如，你假设信用评分不会导致结果，主要是因为我这么告诉你，而且你信任我，因为我是生成数据的人。但在现实生活中，很难找到这样的工具变量。大多数情况下，一个协变量会同时影响依从性和结果，就像这里的年龄情况一样。因此，一种更有希望降低工具变量（IV）估计方差的方法是纳入对结果有高度预测性的控制变量。例如，在这个例子中，客户收入对购买量有很强的预测性，所以将其作为额外的控制变量纳入会大幅降低标准误：
formula = '''pv ~ 1
+ [prime_card ~ prime_elegible + credit_score]
+ income + age'''
iv_model = IV2SLS.from_formula(formula, df).fit(cov_type="unadjusted")
iv_model.summary.tables[1]| Parameter | Std. Err. | T-stat | P-value | Lower CI | Upper CI | |
|---|---|---|---|---|---|---|
| Intercept | 210.62 | 37.605 | 5.6008 | 0.0000 | 136.91 | 284.32 |
| income | 0.3998 | 0.0008 | 471.04 | 0.0000 | 0.3981 | 0.4014 |
| age | 9.7444 | 0.8873 | 10.982 | 0.0000 | 8.0053 | 11.483 |
| prime_card | 693.12 | 12.165 | 56.978 | 0.0000 | 669.28 | 716.96 |
对于像年龄这样既影响结果又影响依从性的变量，其对标准误的影响会更加微妙。就像在回归的情况中，如果它对处理的解释程度远超过对结果的解释程度，最终可能会增加方差。
手动实现两阶段最小二乘法
由于你并非总能用到工具变量相关的软件，我认为有必要学习如何手动实现两阶段最小二乘法，尤其是因为它一点也不复杂。如果你有不止一个工具变量以及额外的协变量，你可以通过以下方式将它们纳入你的模型：
进行第一阶段，将处理变量对工具变量和额外的协变量进行回归，即\(T \sim Z + X\)。
进行第二阶段，将结果变量对第一阶段得到的处理变量拟合值以及额外的协变量进行回归，即\(Y \sim T_{hat} + X\)：
formula_1st = "prime_card ~ prime_elegible + credit_score + income+age"
first_stage = smf.ols(formula_1st, data=df).fit()
iv_model = smf.ols(
"pv ~ prime_card + income + age",
data=df.assign(prime_card=first_stage.fittedvalues)).fit()
iv_model.summary().tables[1]| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 210.6177 | 40.832 | 5.158 | 0.000 | 130.578 | 290.657 |
| prime_card | 693.1207 | 13.209 | 52.474 | 0.000 | 667.229 | 719.013 |
| income | 0.3998 | 0.001 | 433.806 | 0.000 | 0.398 | 0.402 |
| age | 9.7444 | 0.963 | 10.114 | 0.000 | 7.856 | 11.633 |
矩阵实现
这会给你与用linearmodels得到的完全相同的工具变量（IV）估计值，但标准误会有偏差。如果你想要准确的标准误，你可能最好使用两阶段最小二乘法（2SLS）的矩阵实现。要做到这一点，你必须将额外的协变量添加到处理矩阵和工具变量矩阵中。然后，你可以按如下方式求出
IV 参数：
\(\widehat{X} = Z(Z^{\prime}Z)^{-1}Z^{\prime}X\)
\(\beta_{IV} = \left( \widehat{X}^{\prime}\widehat{X} \right)^{-1}\widehat{X}^{\prime}Y\)
在编写代码时，你只需要注意大样本量N的情况。\(Z(Z^{\prime}Z)^{-1}Z^{\prime}\)会是一个巨大的\(N \times N\)矩阵，不过，如果你先计算\((Z^{\prime}Z)^{-1}Z^{\prime}X\)，然后再将其与Z左乘，就可以避免这种情况：
Z = df[["prime_elegible", "credit_score", "income", "age"]].values
X = df[["prime_card", "income", "age"]].values
Y = df[["pv"]].values
def add_intercept(x):
return np.concatenate([np.ones((x.shape[0], 1)), x], axis=1)
Z_ = add_intercept(Z)
X_ = add_intercept(X)
# pre-multiplying Z_.dot(...) last is important to avoid
# creating a huge NxN matrix
X_hat = Z_.dot(np.linalg.inv(Z_.T.dot(Z_)).dot(Z_.T).dot(X_))
b_iv = np.linalg.inv(X_hat.T.dot(X_hat)).dot(X_hat.T).dot(Y)
b_iv[1]## array([693.12072518])再一次，你得到了和之前完全相同的系数。一旦得到了这个系数，你就可以计算工具变量（IV）的残差和方差：
\(\widehat{Var}\left( \widehat{\beta}_{IV} \right) = \sigma^2 \left( \hat{e}_{iv} \right) diag\left( \left( \widehat{X}^{\prime}\widehat{X} \right)^{-1} \right)\)
e_hat_iv = (Y - X_.dot(b_iv))
var = e_hat_iv.var()*np.diag(np.linalg.inv(X_hat.T.dot(X_hat)))
np.sqrt(var[1])## 12.164694395031463由于矩阵符号的存在，这个方差公式有点难以解释，但在没有额外协变量的情况下，你可以近似得到更符合之前结果的式子：
\(SE \left( \widehat{\beta}_{IV} \right) \approx \frac{\sigma \left( \hat{e}_{IV} \right)}{\sigma \left( \widetilde{T} \right) \sqrt{n R_{1st}^2}}\)
这里，\(\widetilde{T}\)是处理变量对额外协变量（而非工具变量）回归的残差，\(R_{1st}^2\)是第一阶段的\(R^2\)：
t_tilde = smf.ols("prime_card ~ income + age", data=df).fit().resid
e_hat_iv.std()/(t_tilde.std()*np.sqrt(n*first_stage.rsquared))## 12.156252763191322通过这个公式，你可以看到，除了增加样本量之外，你有三个途径来降低标准误：
提高第一阶段的\(R^2\)。这可以通过寻找强工具变量来实现，强工具变量是善于预测依从性且满足排除限制（不会导致结果）的变量。
移除对T高度预测的变量，以增加\(\sigma \left( \widetilde{T} \right)\)。
减小第二阶段残差的大小，这可以通过寻找对结果高度预测的变量来实现。
在这三个途径中，我承认我只喜欢最后一个。正如我之前所说，在实际中很难找到工具变量（IV）。而且，为了减小\(\sigma \left( \widetilde{T} \right)\)，你能移除的变量也是有限的。这就使得找到善于预测结果的变量，成为了你降低方差的唯一可靠方法。
不连续性设计
除了传统的工具变量和不依从设计之外，回归不连续性设计（RDD）是另一种值得提及的设计。尽管 RDD 在学术界被广泛使用，但其在工业界的应用可能更有限。RDD 利用处理分配中的人为不连续性来识别处理效应。例如，假设政府实施了一项资金转移计划，向贫困家庭提供每月相当于 200 美元当地货币的支票，但只有收入低于 50 美元的家庭才有资格。这在该计划的分配上于 50 美元处产生了一个不连续性，使得研究人员能够比较刚好高于和刚好低于阈值的家庭，以衡量该计划的有效性，前提是这两组家庭是相似的。
除了资金转移项目的例子之外，回归不连续性设计（RDD）还可应用于许多其他情况。不连续性是普遍存在的，这使得 RDD 对研究人员极具吸引力。例如，为了了解大学的影响，研究人员可以比较在入学考试中分数刚好高于和刚好低于及格阈值的人。为了评估女性对政治的影响，研究人员可以比较女性候选人以微弱劣势落败的城市和女性候选人以微弱优势获胜的城市。其应用是无穷无尽的。
回归不连续性设计（RDD）在工业界也可能有用，但应用程度较低。例如，假设一家银行为其所有客户提供信用卡，但会向账户余额低于 5000 美元的客户收取费用。这就在信用卡的提供方式上造成了一种不连续性，账户余额高于阈值的客户更有可能选择高端信用卡，而余额低于阈值的客户则不太可能。因此，只要阈值上下的客户在其他方面是相似的，RDD 就可以用于比较拥有高端信用卡和普通信用卡的影响。
关于不连续性设计在工业界的相关性，我认为它的适用性较低，因为企业可以很容易地进行实验来随机化资格，正如我们之前所讨论的。不过，为了这个例子，让我们假设进行这样的实验会很耗时。也许是因为由于依从性低，所需的样本量太大。
相比之下，相关银行已经拥有遵循之前描述的不连续性设计的数据。因此，银行可以利用这些数据来确定高端信用卡的效果。银行如何为此目的利用这种不连续性呢？基本思路是认识到这个阈值可以被理解为一个工具变量，因为越过它会增加接受处理的可能性。
在下面的图像中，你可以看到不连续性设计与工具变量是如何关联的。下方部分展示了按账户余额划分的反事实处理情况。由于工具变量是跨越 5000 美元的阈值，当余额小于 5000 时，你可以观察到\(T_0\)，否则观察到\(T_1\)。此外，由于该工具变量增加了获得处理（高端信用卡）的可能性，一旦越过这个阈值，\(P(T = 1)\)就会出现一个跃升。图表的上方部分反映了处理概率的这些变化是如何影响结果的。
即使在阈值以上，处理的概率仍小于 1，这使得你观察到的结果小于真实的潜在结果\(Y_1\)。同理，你在阈值以下观察到的结果高于真实的潜在结果\(Y_0\)。这就使得阈值处的处理效应看起来比实际的要小，因此你必须使用工具变量（IV）来对此进行修正。
不连续性设计的假设
除了工具变量（IV）的假设之外，不连续性设计还需要关于潜在结果和潜在处理函数平滑性的进一步假设。我们定义一个运行变量R，使得在阈值\(R = c\)处，处理概率是该变量的不连续函数。在你的银行业例子中，R是账户余额，\(c = 5000\)。
现在，你需要假设：
\(\begin{align*} \lim_{r \to c^-} E[Y_t | R = r] &= \lim_{r \to c^+} E[Y_t | R = r] \\ \lim_{r \to c^-} E[T_z | R = r] &= \lim_{r \to c^+} E[T_z | R = r] \end{align*}\)
换句话说，在不连续点\(R = c\)处，潜在结果\(Y_t\)和潜在处理\(T_z\)，当从左侧或右侧趋近该点时是相同的。
有了这些假设，你可以推导出不连续性设计的局部平均处理效应估计量：
\(\begin{align*} LATE &= \frac{\lim_{r \to c^+} E[Y|R = r] - \lim_{r \to c^-} E[Y|R = r]}{\lim_{r \to c^+} E[T|R = r] - \lim_{r \to c^-} E[T|R = r]} \\ &= E[Y_1 - Y_0 | T_1 > T_0, R = c] \end{align*}\)
重要的是，这个估计量在两个意义上是局部的。首先，它是局部的，因为它只给出阈值\(R = c\)处的处理效应，这是不连续性设计的局部性。其次，它是局部的，因为它只估计依从者的处理效应，这是工具变量（IV）的局部性。
意向治疗效应
在图 11-4 的上方部分，阈值处观测结果的跃变就是意向治疗效应，因为它衡量了当你改变工具变量时，结果是如何变化的。现在，让我们看看如何估计它，因为这将是你最终工具变量（IV）估计值的分子。为了做到这一点，我们首先读取包含客户账户余额、是否选择高端信用卡以及他们的购买量等信息的数据：
df_dd = pd.read_csv("./data/prime_card_discontinuity.csv")
df_dd.head()## balance prime_card pv tau categ
## 0 12100.0 1 356.472 300.0 always-takers
## 1 4400.0 1 268.172 300.0 always-takers
## 2 4600.0 1 668.896 300.0 always-takers
## 3 3500.0 1 428.094 300.0 always-takers
## 4 12700.0 1 1619.793 700.0 complier接下来，你需要将结果变量对运行变量R（余额）与一个表示是否高于阈值（\(R > c\)）的虚拟变量的交互项进行回归：
\(y_i = \beta_0 + \beta_1 r_i + \beta_2 \mathbb{1}(r_i > c) + \beta_3 \mathbb{1}(r_i > c) r_i\)
与跨越阈值相关的参数估计\(\widehat{\beta_2}\)，可以被解释为意向治疗效应。
m = smf.ols(f"pv~balance*I(balance>0)", df_dd.assign(balance = lambda d: d["balance"] - 5000)).fit()
m.summary().tables[1]| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 559.2980 | 8.395 | 66.621 | 0.000 | 542.843 | 575.753 |
| I(balance > 0)[T.True] | 261.0699 | 10.128 | 25.777 | 0.000 | 241.218 | 280.922 |
| balance | 0.0616 | 0.005 | 11.892 | 0.000 | 0.051 | 0.072 |
| balance:I(balance > 0)[T.True] | -0.0187 | 0.005 | -3.488 | 0.000 | -0.029 | -0.008 |
请注意，这本质上是拟合两条回归直线：一条在阈值之上，一条在阈值之下。如果依从性不是问题，也就是说，阈值之上的每个人都会接受处理，阈值之下的每个人都会处于控制组，你仍然可以使用这种方法。如果是这种情况，依从性将是 100%，意向治疗效应（ITTE）就已经是平均处理效应（ATE）了。
plt_df = df_dd.round({"balance": -2}).assign(size=1).groupby("balance").agg({"pv":"mean", "size": "sum"}).reset_index()
plt.figure(figsize=(8,4))
sns.scatterplot(data=plt_df, y="pv", x="balance", size="size", color="C5")
plt.plot(plt_df.query("balance<5000")["balance"], m.predict(plt_df.query("balance<5000").assign(balance = lambda d: d["balance"] - 5000)), color="C0", lw=2, label="regression")
plt.plot(plt_df.query("balance>5000")["balance"], m.predict(plt_df.query("balance>5000").assign(balance = lambda d: d["balance"] - 5000)), color="C0", lw=2)
plt.legend(fontsize=14)
工具变量估计
由于依从性并非 100%，你需要将意向治疗效应除以依从率。在不连续性设计的背景下，这指的是当你跨越阈值时，治疗概率的变化程度。为了估计这个数值，你可以简单地重复之前的步骤，将结果变量 “pv” 替换为治疗变量 “prime_card”。这里有一个用于计算不连续性设计中工具变量估计值的简单函数。它会估计意向治疗效应和依从率，然后将前者除以后者：
def rdd_iv(data, y, t, r, cutoff):
centered_df = data.assign(**{r: data[r]-cutoff})
compliance = smf.ols(f"{t}~{r}*I({r}>0)", centered_df).fit()
itte = smf.ols(f"{y}~{r}*I({r}>0)", centered_df).fit()
param = f"I({r} > 0)[T.True]"
return itte.params[param]/compliance.params[param]
rdd_iv(df_dd, y="pv", t="prime_card", r="balance", cutoff=5000)## 732.8534752298925将这个函数应用到你的数据上，会得到一个与真实局部平均处理效应（LATE）非常接近的估计值。要记住，你可以对这一点进行检验，因为这个数据集包含了存储在 “tau” 列中的个体水平处理效应以及依从类别。
最后，尽管你可以推导出一个公式来计算该估计量的置信区间，但最简单的方法是直接将整个函数套用到自助法（bootstrap）流程中。我会把这部分代码隐藏起来，因为它相当重复，但你可以在这里看到得出的区间：
(df_dd
.round({"balance":-2}) # round to nearest hundred
.query("balance==5000 & categ=='complier'")["tau"].mean())## 700.0聚堆现象
在结束本章之前，我只想提一下不连续性设计识别中存在的一个潜在问题。如果个体（在你的例子中是客户）能够操纵运行变量，他们也能自行选择进入处理组。在高端信用卡的例子中，客户可能会决定增加存款，直到存款额刚好达到 5000，这样他们就能免费获得高端信用卡。这会违反关于潜在结果平滑性的假设，因为刚好在阈值之上的那些个体将不再与刚好在阈值之下的个体具有可比性。
一种简单且直观的检查这种情况是否发生的方法是绘制阈值附近的密度图。如果个体在自行选择进入处理组，你会预期在阈值处的密度会有一个巨大的峰值。幸运的是，这个数据似乎不存在这种情况：
核心要点
在本章中，你了解到当人们可以选择不接受某种处理时，不依从就会成为一个问题。这在行业中相当常见，因为企业往往有一系列可选的产品或服务。在这些情况下，即便企业可以对产品或服务的可获得性进行随机化处理，客户的选择也会混淆产品或服务的效果。
你还了解了依从性群体或类型：
依从者：接受分配给他们的处理的人。
总是接受者：无论分配如何，总是接受处理的人。
从不接受者：无论分配如何，从不接受处理的人。
违抗者：接受与分配相反的处理的人。
而且你还了解了如何使用工具变量来应对不依从问题。也就是说，工具变量Z是这样一种变量：（1）以无混淆的方式影响处理；（2）不会直接影响结果，除非通过处理。
除此之外，如果你假设工具变量在单一方向上改变处理（单调性假设），你可以用它来识别依从者的平均处理效应：
\(E\left[Y_1 - Y_0 \big| T_1 > T_0\right] = \frac{E[Y \big| Z = 1] - E[Y \big| Z = 0]}{E[T \big| Z = 1] - E[T \big| Z = 0]}\)
换句话说，你所要做的就是用依从率来对意向治疗效应进行标准化，而如果工具变量是随机化的，这两者都很容易识别。
然而，在方差方面仍然需要付出代价。如果依从性低，工具变量估计量的方差将大大超过普通最小二乘法（OLS）的方差。特别是，如果依从性是 50%，要达到与无依从性问题（100% 依从性）时相同的标准误，你需要 4 倍多的样本量。有一些额外的技巧可以降低方差，但最有希望的似乎是找到善于预测结果的变量，这和普通最小二乘法的情况大致相同。
此外，你还了解到数据中的不连续性也可以被当作工具变量。一般来说，你可能不需要依赖它们，因为在行业中开展实验是相当常见且容易的。不过，在实验不可行的情况下，你可以利用这些不连续性来识别局部平均处理效应。
实际案例：出生季度工具变量
正如我之前所说，在现实中很难找到有效的工具变量，但出生季度可能是其中之一。在美国，出生在最后一个季度意味着你可能会有更多的学龄，因为你会在人生中更早地入学。如果出生季度不影响收入（除非通过教育），并且几乎是随机的，经济学家就可以用它来识别教育对收入的影响。
通过这种方法，经济学家估计，平均而言，多一年的教育应该能使工资提高 8.5%：
