【python因果推断库12】工具变量回归与使用 pymc 验证工具变量5
目录
可信推断与可信度革命
比较模型推断
在结果空间中评估模型
可信推断与可信度革命
我们可以将这一想法推进到什么程度？我们的推断对模型误指定有多脆弱？我们能否通过尝试强烈的先验来压力测试参数估计？IV 方法论明显是关于论证——对提出的机制的可信论证。鉴于此，贝叶斯建模方法应用于 IV 的一个好处是，我们可以在模型设计中表达并扩展机制的可信度。我们可以通过尝试将争议性的信念作为模型的先验，并观察推断在多大程度上由我们的先验规范锚定，以及数据在多大程度上使我们远离不太可信的假设，来压力测试推断的可信度。
在这里，我们将重新拟合我们最初的 IV 模型，但是我们会将 `experience_1` 变量缩放，使其均值为 0。这样我们就可以为所有变量放置大致相同规模的先验。
instruments_formula = """education ~ experience_1 + experience_2 + ethnicity_indicator + south_indicator +
smsa_indicator + nearcollege_indicator
"""
formula = "log_wage ~ 1 + education + experience_1 + experience_2 + ethnicity_indicator + south_indicator + smsa_indicator"
instruments_data = df[
[
"education",
"nearcollege_indicator",
"experience_1",
"experience_2",
"ethnicity_indicator",
"smsa_indicator",
"south_indicator",
]
]
data = df[
[
"log_wage",
"education",
"experience_1",
"experience_2",
"ethnicity_indicator",
"smsa_indicator",
"south_indicator",
]
]
scaler = StandardScaler()
data["experience_1"] = scaler.fit_transform(data[["experience_1"]])
instruments_data["experience_1"] = scaler.fit_transform(
instruments_data[["experience_1"]]
)
data[["log_wage", "education", "experience_1", "experience_2"]].head(5)
然后我们将修改 LATE 估计的先验，将其向上调整至 0.2。出于论证的目的，尝试将推断向上拉动。
modified = iv.ols_beta_second_params
modified[1] = 0.20
modified[2] = 0
modified[3] = 0
modified[4.486489832009674, 0.2, 0, 0, -0.13080191375573674, -0.10490054155737207, 0.13132367504470194]
modified_first = iv.ols_beta_first_params
modified_first[1] = 0
modified_first[13.098074632594448, 0, -0.8633797925737376, -1.0061382678610173, -0.29146401700569985, 0.40387687667134015, 0.33732078008876376]
此外，我们还将通过提高 eta 参数并减小多元正态分布的标准差来约束我们认为可接受的相关性的可能范围。
iv2 = InstrumentalVariable(
instruments_data=instruments_data,
data=data,
instruments_formula=instruments_formula,
formula=formula,
model=InstrumentalVariableRegression(sample_kwargs=sample_kwargs),
priors={
"mus": [modified_first, modified],
"sigmas": [1, 1],
"eta": 10,
"lkj_sd": 0.5,
},
)
az.summary(iv2.idata, var_names=["beta_t", "beta_z"])[
["mean", "sd", "hdi_3%", "hdi_97%", "r_hat"]
]
我们在这里看到 LATE 估计 beta_z[education] 从之前模型记录的值有所缩小，但仍一致高于通过简单控制回归得到的估计。
比较模型推断
为了使情况更加清晰，我们提取并绘制每个模型所暗示的参数估计的置信区间。在这里我们可以看到每个模型配置带来的不同含义。请注意，所有的 IV 模型都将置信区间从更简单的 OLS 类型模型的估计中拉开。
az.plot_forest(
[idata_reg, iv.idata, iv1.idata, iv2.idata],
var_names=["beta_z"],
coords={
"covariates": [
"education",
"ethnicity_indicator",
"smsa_indicator",
"south_indicator",
]
},
combined=True,
model_names=["OLS", "IV", "IV1", "IV2"],
figsize=(8, 7),
);
我们最终模型中的不确定性增加是由相关参数的不确定性驱动的。但是尽管我们试图将 LATE 估计向上拉，贝叶斯更新的过程还是将我们的置信区间的高密度区间 (HDI) 拉入了一个与我们第一次 IV 模型拟合有很大重叠的合理范围内。
az.plot_forest(
[iv.idata, iv1.idata, iv2.idata],
var_names=["chol_cov_corr"],
coords={"chol_cov_corr_dim_1": [1], "chol_cov_corr_dim_0": [0]},
combined=True,
model_names=["IV", "IV1", "IV2"],
figsize=(8, 4),
);
最后一个模型中所隐含的相关性不确定性某种程度上削弱了这个模型的设定。如果我们关于工具变量的论证要有说服力，我们期望相关性假设成立。通过降低相关性来削弱相关性的模型设定可能过于极端。实际上，我们降低了工具变量的相关性，但仍得到了 `beta_z[education]` 的强烈正效应。这里的关键不是争论参数设置，而是展示需要考虑多个模型，并且在验证工具变量设计时总是需要进行一些敏感性测试。
iv1.model.sample_predictive_distribution(ppc_sampler="jax")在结果空间中评估模型
比较每个模型可信度的另一种方式是将模型的含义转化为结果空间，并思考参数估计的差异真正意味着什么。下面我们将可视化每个模型下因教育效应而导致的对数工资的影响。
def make_compare_plot(iv, y, model_name, covariates):
data = df.copy()
data["Intercept"] = 1
covariate_df = data[covariates]
params = az.extract(iv["posterior"]["beta_z"])["beta_z"]
fig, axs = plt.subplots(2, 1, figsize=(10, 9))
axs = axs.flatten()
for i in range(100):
if i == 99:
axs[0].hist(
np.exp(np.dot(covariate_df, params[:, i].T)),
alpha=0.3,
bins=30,
ec="black",
color="grey",
rasterized=True,
label="Predicted Histogram",
)
else:
axs[0].hist(
np.exp(np.dot(covariate_df, params[:, i].T)),
alpha=0.3,
bins=30,
ec="black",
color="grey",
rasterized=True,
)
axs[0].set_title(
f"Predicted Wage Distribution {model_name} \n Observed Values",
fontsize=12,
fontweight="bold",
)
axs[0].hist(np.exp(y), bins=30, alpha=0.2, color="blue", label="Observed Histogram")
axs[0].set_xlabel("Predicted Wage", fontsize=9)
axs[0].legend()
for c, ed in zip(["red", "purple", "blue"], [2, 10, 18]):
temp = covariate_df.copy()
temp["education"] = ed
means = []
for i in range(100):
dist = np.exp(np.dot(temp, params[:, i].T))
means.append(np.mean(dist))
if i == 99:
axs[1].hist(
dist,
alpha=0.3,
bins=30,
ec="black",
color=c,
rasterized=True,
label=f"Edu: {ed}, Mean {np.round(np.mean(means))} ",
)
else:
axs[1].hist(
dist, alpha=0.3, bins=30, ec="black", color=c, rasterized=True
)
axs[1].axvline(np.mean(means), color="k", rasterized=True)
axs[1].set_title(
f"Predicted Wage Distribution {model_name} \n Counterfactually Set Education Values",
fontsize=12,
fontweight="bold",
)
axs[1].set_xlabel("Predicted Counterfactual Wage Distribution", fontsize=9)
axs[1].legend()
covariates = list(idata_reg["posterior"]["beta_z"]["covariates"].values)
make_compare_plot(idata_reg, data["log_wage"], "OLS", covariates)
在这里，我们在反事实的情境下展示了简单回归的含义，即数据集中的每个人将其受教育年限分别设置为 2 年、10 年和 18 年。然后绘制期望值的差异。接下来，我们展示了候选 IV 模型的相同情况。
covariates = list(iv.idata["posterior"]["beta_z"]["covariates"].values)
make_compare_plot(iv.idata, data["log_wage"], "IV", covariates)
在这里，我们可以看到所暗示分布形状之间存在着明显的差异，以及在后验预测图中对观测数据有了更好的重现。
covariates = list(iv1.idata["posterior"]["beta_z"]["covariates"].values)
make_compare_plot(iv1.idata, data["log_wage"], "IV1", covariates)
我们的第二个 IV 模型与第一个模型讲述了一个类似的故事，但在覆盖对数工资分布的长尾特征方面有所改进。
covariates = list(iv2.idata["posterior"]["beta_z"]["covariates"].values)
make_compare_plot(iv2.idata, data["log_wage"], "IV2", covariates)
我们的最终模型显示了预测分布和反事实分布有些混乱的实现，表明它实际上并不是我们数据的好模型。我们也可以通过它们的预测能力来正式比较这些模型。
compare_df = az.compare({"IV": iv.idata, "IV1": iv1.idata, "IV2": iv2.idata})
compare_df[["rank", "elpd_loo", "p_loo", "elpd_diff", "weight"]]
最终分析应由原始工资尺度上的含义的合理性驱动。但了解特定模型拟合与其它合理候选模型在预测能力方面的比较是有用的。