目录

工具变量回归与使用 pymc 验证工具变量

回归机制与局部平均处理效应

旁白：从多元正态分布中采样

---

import arviz as az
import daft
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import scipy
from matplotlib.lines import Line2D
from matplotlib.patches import Ellipse
from sklearn.preprocessing import StandardScaler

import causalpy as cp
from causalpy.pymc_experiments import InstrumentalVariable
from causalpy.pymc_models import InstrumentalVariableRegression
%config InlineBackend.figure_format = 'retina'
seed = 42
np.random.seed(seed)

对于这个笔记，我们还安装了 jax 和 numpyro。这使得我们在 `InstrumentalVariable` 类中对后验预测分布的采样速度大大加快。您可以在 Jupyter 笔记本中运行以下命令来安装这两个库：

!conda install jax
!conda install numpyro

工具变量回归与使用 pymc 验证工具变量

决定学业成功和后续就业的因素是多样的。有许多社会人口统计特征和个人能力方面的因素在起作用。因此，如何干净地评估教育的影响存在合理的担忧。如果不设法考虑能力和教育之间的混淆关系，我们可能会使结果出现偏差。在这个例子中，我们将探讨著名的经济学问题之一：投资于教育所能获得的经济回报。我的教育程度对我一生财富的预期增加有何影响？经济学文献的建议是利用工具变量回归来分解教育对未来工资获取的影响。

具体来说，建议我们可以使用一个人距离四年制大学的（假设如同随机的）距离作为一个工具变量。回想一下，这个工具变量 `nearcollege` 必须满足：(i) 相关性，即它会影响我们感兴趣的因果路径，例如靠近大学被认为会对个人的教育产生因果影响；(ii) 独立性，即分配给工具变量应尽可能随机，从而减轻原始的选择效应偏差；(iii) 需要满足排除限制，即工具变量应该仅仅通过我们感兴趣的处理变量（教育）来影响结果（工资）。如果这些条件得到满足，那么我们可以认为工具变量回归技术能够提炼出处理变量的因果影响。

GRID_UNIT = 2.0
DPI = 200
NODE_EC = "none"

pgm = daft.PGM(dpi=DPI, grid_unit=GRID_UNIT, node_ec=NODE_EC)

pgm.add_node("z", "$Z$", 0, 0)
pgm.add_node("t", "$T$", 1, 0)
pgm.add_node("y", "$Y$", 2, 0)
pgm.add_node("x", "$\mathbf{X}$", 1.5, 0.75)
pgm.add_edge("z", "t")
pgm.add_edge("x", "t")
pgm.add_edge("x", "y")
pgm.add_edge("t", "y")


pgm.add_node("nearcollege", "$Near$", 3, 0)
pgm.add_node("education", "$Edu$", 4, 0)
pgm.add_node("wage", "$Wage$", 5, 0)
pgm.add_node("ability", "$\mathbf{ability}$", 4.5, 0.75)
pgm.add_edge("nearcollege", "education")
pgm.add_edge("ability", "education")
pgm.add_edge("ability", "wage")
pgm.add_edge("education", "wage")

pgm.render();

在许多方面，工具变量回归是经济学和社会科学中可信度革命的一个典范方法。对社会科学中可信度的需求源自复制危机。下面我们将探讨工具变量设计的假设，并强调在哪些条件下它可以合理地应用。首先，我们将概述工具变量回归中因果估计量的本质。然后，我们将探讨这些估计程序是如何根据上述假设工作的。我们将展示不同的方法来检验这些模型，并在类似的模型之间进行比较，以论证结果的可信度。

回归机制与局部平均处理效应

工具变量问题的结构帮助我们探究因果影响的问题，但它并不针对典型的平均处理效应<ATE>因果估计量。相反，我们必须满足于估计局部平均处理效应<LATE>量，在我们的设计中的“遵从者”群体中评估受限的处理效应。这是有道理的，因为我们隐含地评估了一个没有遵守保证的群体中的处理效应。更具体地说，我们的工具变量，由于其独立性，旨在在一个响应我们工具变量的子人群中诱导出如同随机分配的处理。所以我们不能做出更广泛的声称已经评估了平均变化。相反，它是对那些被诱导遵守我们设计的子群体中变化的评估。这些个体的处理（教育）受到了我们的工具变量（nearcollege）的影响。

LATE估计量的理论是非常丰富的。有多种方法可以隔离LATE量，也有多种LATE估计量。这里我们展示了一种经典的表示方法，即在我们工具变量的实现下的结果和处理的条件期望。

最后的表述突出了我们对工具变量 `Z` 引起的两个变化成分感兴趣的事实。这是一种看待我们试图在一个结构方程模型中估计两个随机变量的联合分布的方式。但更重要的是，我们通过相关性假设认为工具变量 `Z` 的实现与处理变量 `T` 之间存在一些非平凡的相关性。因此，选择与感兴趣的子群体相关的工具变量或预期遵守率高的工具变量尤为重要。

旁白：从多元正态分布中采样

我们如何衡量工具变量与处理变量之间的相关性？我们应该期望多少相关性？在 CausalPy 的工具变量回归实现中，我们使用 LKJ Cholesky 先验对多元正态分布建模，以显式地模拟这种相关性。在这里值得做一个简短的旁白，展示在不同先验下从这种分布中采样如何塑造联合分布的相关性。下面我们将展示这为我们提供了一种机制，以对我们关于工具变量之间关系的信念施加约束。

sd_dist = pm.Exponential.dist(1.0, shape=2)
lkj = pm.LKJCholeskyCov.dist(
    eta=2,
    n=2,
    sd_dist=sd_dist,
)
lkj1 = pm.LKJCholeskyCov.dist(
    eta=10,
    n=2,
    sd_dist=sd_dist,
)
chol, corr, sigmas = pm.draw(lkj, draws=100)
chol1, corr1, sigmas = pm.draw(lkj1, draws=100)

fig, axs = plt.subplots(4, 1, figsize=(9, 12))
axs = axs.flatten()
corrs = []
corrs1 = []
for i in range(100):
    xy = pm.draw(pm.MvNormal.dist(0, chol=chol[i]), 1000)
    xy1 = pm.draw(pm.MvNormal.dist(0, chol=chol1[i]), 1000)
    corrs.append(np.corrcoef(xy[:, 0], xy[:, 1])[0][1])
    corrs1.append(np.corrcoef(xy1[:, 0], xy1[:, 1])[0][1])

    axs[0].scatter(xy[:, 0], xy[:, 1], alpha=0.1, color="blue", rasterized=True)
    axs[1].scatter(xy1[:, 0], xy1[:, 1], alpha=0.1, color="red", rasterized=True)
    axs[2].scatter(xy[:, 0], xy[:, 1], alpha=0.1, color="blue", rasterized=True)
    axs[2].scatter(xy1[:, 0], xy1[:, 1], alpha=0.1, color="red", rasterized=True)
axs[3].hist(
    np.array(corrs),
    bins=30,
    ec="black",
    label="Correlation of Observations LKJ(eta=2)",
    color="blue",
    alpha=0.5,
    rasterized=True,
)
axs[3].hist(
    np.array(corrs1),
    bins=30,
    ec="black",
    label="Correlation of Observations LKJ(eta=10)",
    color="red",
    alpha=0.5,
    rasterized=True,
)
axs[3].legend()
axs[3].set_title("Correlation Coefficient Distribution")
axs[0].set_title("Realisations from LKJ(eta=2)")
axs[2].set_title("Overlapping Realisations")
axs[1].set_title("Realisations from LKJ(eta=10)");

 在上面的一系列图中，我们从 LKJ 先验概率分布的不同参数化中进行了采样。这个分布是对多元正态概率分布的期望值和协方差结构的先验。我们展示了先验的不同参数化会导致所采样的双变量正态分布的相关性更多或更少。在这里我们可以看到，增加 LKJ 先验中的 eta 参数会缩小可接受的相关性参数的范围。默认情况下，CausalPy 中 `InstrumentalVariable` 类的实现是从处理变量和结果变量的双变量正态分布中采样，eta 参数设置为 2。如果你的模型使得这样的潜在相关性不太可能发生，这是值得注意的。下面我们将展示如何在工具变量的背景下对这些参数应用先验。