【python因果推断库7】使用 pymc 模型的工具变量建模 (IV)2

news2024/11/15 19:30:36

与普通最小二乘法 (OLS) 的比较

应用理论：政治制度与GDP

拟合模型：贝叶斯方法

多变量结果和相关性度量

结论

与普通最小二乘法 (OLS) 的比较

simple_ols_reg = sk_lin_reg().fit(X.reshape(-1, 1), y)

print("Intercept:", simple_ols_reg.intercept_, "Beta:", simple_ols_reg.coef_[0])

Intercept: 0.5677845035965572 Beta: 4.427701928515228

我们可以看到，由于处理变量的内生性，OLS 过高地估计了焦点参数的值，而 IV 回归则更接近真实值。正是这种偏误，工具变量设计的目的就是要缓解。

应用理论：政治制度与GDP

提醒一下，我们想要建模以下关系：

iv_df = cp.load_data("risk")
iv_df[["longname", "loggdp", "risk", "logmort0"]].head()

当我们观察到政治制度在这个增长系统中是内生的时候，问题就出现了。这意味着我们需要以某种方式控制测量误差和偏误，如果我们简单地拟合一个OLS模型的话。他们继续论证说，我们可以使用一个工具变量，这个变量仅通过政治制度的程度与GDP相关联，通过使用工具变量回归。他们最终建议使用欧洲定居者在那个时期的死亡率作为工具变量，因为较高的死亡率会导致较少的移民和对该地区的投资，这应该会减少在殖民地建立的政治制度。他们可以使用军事记录来收集这些数据。

我们可以手动估计两阶段最小二乘法 (2SLS) 的处理效应如下：

X = iv_df.risk.values.reshape(-1, 1)
Z = iv_df.logmort0.values.reshape(-1, 1)
t = iv_df.risk.values
y = iv_df.loggdp.values

simple_ols_reg = sk_lin_reg().fit(X, y)
first_stage_reg = sk_lin_reg().fit(Z, t)
fitted_risk_values = first_stage_reg.predict(Z)


second_stage_reg = sk_lin_reg().fit(X=fitted_risk_values.reshape(-1, 1), y=y)

print(
    "Simple OLS Parameters: Intercept and Beta Coeff",
    simple_ols_reg.intercept_,
    simple_ols_reg.coef_,
)
print(
    "First Stage Parameters: Intercept and Beta Coeff",
    first_stage_reg.intercept_,
    first_stage_reg.coef_,
)
print(
    "Second Stage Parameters Intercept and Beta Coeff",
    second_stage_reg.intercept_,
    second_stage_reg.coef_,
)

Simple OLS Parameters: Intercept and Beta Coeff 4.687414702305412 [0.51618698]
First Stage Parameters: Intercept and Beta Coeff 9.365894904697788 [-0.61328925]
Second Stage Parameters Intercept and Beta Coeff 1.9942956864448975 [0.92948966]

请注意，朴素的OLS估计值0.515与2SLS估计值0.92在处理效应上的显著差异。这与论文中报告的结果相符。

在这个笔记中，我们不会进一步讨论弱工具变量和强工具变量的问题，也不会讨论如何找到并测试工具变量的强度，但我们将会展示如何在贝叶斯设置下拟合这类模型。我们还将讨论贝叶斯方法如何在幕后将焦点（第二阶段）回归和工具（第一阶段）回归建模为具有明确相关性的多元随机变量。想法是将这两个结果一起建模，并带有明确的相关性。这种方法的好处是我们可以获得关于“工具”和结果之间关系的额外见解。

$\begin{aligned}\begin{pmatrix}y\\t\end{pmatrix}&\sim\text{MultiNormal}(\mu,\Sigma)\\\mu&=\begin{pmatrix}\mu_y\\\mu_t\end{pmatrix}=\begin{pmatrix}\beta_{00}+\beta_{01}t\\\beta_{10}+\beta_{11}z\end{pmatrix}\end{aligned}$

在这个实现的选择上，我们遵循了Juan Orduz博客的例子，该例子又借鉴了Jim Savage的工作。这样做有一个好处，那就是能够明确地表达出我们对处理变量和工具变量联合分布的兴趣。

拟合模型：贝叶斯方法

我们使用CausalPy来处理我们的数据，具体如下：

sample_kwargs = {"tune": 1000, "draws": 2000, "chains": 4, "cores": 4}
instruments_formula = "risk  ~ 1 + logmort0"
formula = "loggdp ~  1 + risk"
instruments_data = iv_df[["risk", "logmort0"]]
data = iv_df[["loggdp", "risk"]]
iv = InstrumentalVariable(
    instruments_data=instruments_data,
    data=data,
    instruments_formula=instruments_formula,
    formula=formula,
    model=InstrumentalVariableRegression(sample_kwargs=sample_kwargs),
)


az.plot_trace(iv.model.idata, var_names=["beta_z", "beta_t"]);

az.summary(iv.model.idata, var_names=["beta_t", "beta_z"])[
    ["mean", "sd", "hdi_3%", "hdi_97%", "r_hat"]
]

hdi_prob = 0.94
ax = az.plot_posterior(
    data=iv.model.idata,
    var_names=["beta_z"],
    hdi_prob=hdi_prob,
)

ax[0].axvline(
    iv.ols_beta_params["Intercept"],
    label="Naive OLS Intercept \n Estimate",
    color="red",
)
ax[1].axvline(
    iv.ols_beta_params[iv.instrument_variable_name],
    label="Naive OLS Treatment \n Estimate",
    color="red",
)
ax[0].axvline(
    iv.ols_beta_second_params[0], label="MLE 2SLS Intercept \n Estimate", color="purple"
)
ax[1].axvline(
    iv.ols_beta_second_params[1], label="MLE 2SLS Treatment \n Estimate", color="purple"
)
ax[0].legend()
ax[1].legend();

多变量结果和相关性度量

正如我们上面所述，贝叶斯方法的一个好处是我们可以直接测量工具变量和处理变量之间的双变量关系。我们可以看到（在二维空间中）估计的处理系数差异如何扭曲预期结果的表示。

az.summary(iv.model.idata, var_names=["chol_cov_corr"])[
    ["mean", "sd", "hdi_3%", "hdi_97%", "r_hat"]
]

fig, axs = plt.subplots(1, 3, figsize=(20, 8))

diffs = (
    iv.model.idata["posterior"]["beta_z"].sel(covariates=[iv.instrument_variable_name])
    - iv.ols_beta_params[iv.instrument_variable_name]
)
axs[0].hist(diffs.values.flatten(), bins=30, ec="black", color="blue", alpha=0.4)
axs[0].axvline(
    np.mean(diffs.values.flatten()),
    label="Expected Diff \n In Treatment Effect \n Estimate",
    color="magenta",
)
axs[0].set_xlabel("Difference")
axs[0].legend()

intercepts = iv.model.idata["posterior"]["beta_z"].sel(covariates=["Intercept"])
betas = iv.model.idata["posterior"]["beta_z"].sel(
    covariates=[iv.instrument_variable_name]
)


raw_df = pd.DataFrame(iv.X, columns=iv.labels)
x = np.linspace(0, 10, 10)
uncertainty = [intercepts.values.flatten() + betas.values.flatten() * i for i in x]
uncertainty = pd.DataFrame(uncertainty).T

ols = [
    iv.ols_beta_params["Intercept"]
    + iv.ols_beta_params[iv.instrument_variable_name] * i
    for i in x
]

custom_lines = [
    Line2D([0], [0], color="orange", lw=4),
    Line2D([0], [0], color="black", lw=4),
]

uncertainty.sample(500).T.plot(legend=False, color="orange", alpha=0.4, ax=axs[1])
axs[1].plot(x, ols, color="black", label="OLS fit")
axs[1].set_title("OLS versus Instrumental Regression Fits", fontsize=20)
axs[1].legend(custom_lines, ["IV fits", "OlS fit"])
axs[1].set_xlabel("Treatment Scale/ Risk")
axs[1].set_ylabel("Outcome Scale/ Log GDP")

axs[0].set_title(
    "Posterior Differences between \n OLS and IV beta coefficients", fontsize=20
)

corr = az.extract(data=iv.model.idata, var_names=["chol_cov_corr"])[0, 1, :]
axs[2].hist(corr, bins=30, ec="black", color="C2", label="correlation")
axs[2].set_xlabel("Correlation Measure")
axs[2].set_title("Correlation between \n Outcome and Treatment", fontsize=20);