【python因果推断库10】工具变量回归与使用 pymc 验证工具变量3

本文链接：https://blog.csdn.net/qq_32146369/article/details/140756584

验证模型

我们从简单的回归情境开始。这有两个目的：(i) 我们可以探索在简单的回归中如何度量教育的影响，并且 (ii) 我们可以基准测试我们的工具变量 `nearcollege_indicator` 在试图预测教育时的有效性。这些回归实际上是我们的工具变量的诊断测试。接下来，我们将分别查看

(i) LATE 估计的第一阶段回归，
(ii) 减少形式回归，
(iii) 从 (i) 和 (ii) 得到的 LATE 量的 Wald 估计，
(iv) 将 (iii) 与简单回归调整估计的教育效果进行比较。

第一阶段

通过拟合第一阶段回归，我们可以评估相关性假设。这个假设是可以检验的，因为我们想要确保分配给我们的工具变量的系数权重不为零，从而证明工具变量的加入可以帮助预测处理变量的部分变异。通常会包括类似于 F 检验的检查，以帮助论证你的工具变量确实值得被纳入你对处理的理解之中。下面我们将展示如何执行这项练习。

让我们首先看一下当我们在回归中加入人口统计控制变量时，这如何运作。

covariate_df = df[
    [
        "experience_1",
        "experience_2",
        "ethnicity_indicator",
        "smsa_indicator",
        "south_indicator",
        "nearcollege_indicator",
    ]
].copy()
covariate_df["Intercept"] = 1
covariate_df = covariate_df[
    [
        "Intercept",
        "experience_1",
        "experience_2",
        "ethnicity_indicator",
        "smsa_indicator",
        "south_indicator",
        "nearcollege_indicator",
    ]
]


def make_reg_model(covariate_df, target="log_wage", prior_beta=[0, 1]):
    coords = {"covariates": list(covariate_df.columns)}
    with pm.Model(coords=coords) as reg_model:
        X = pm.Data("X", covariate_df)
        beta = pm.Normal("beta_z", prior_beta[0], prior_beta[1], dims="covariates")
        mu = pm.Deterministic("mu", pm.math.dot(X, beta))
        sigma = pm.TruncatedNormal("sigma", 1, lower=0.005)
        _ = pm.Normal("likelihood", mu, sigma, observed=df[target].values)
        idata_reg = pm.sample_prior_predictive()
        idata_reg.extend(
            pm.sample(
                nuts_sampler="numpyro",
                target_accept=0.99,
                idata_kwargs={"log_likelihood": True},
            )
        )  ## requires Jax and Numpyro install
    return reg_model, idata_reg


first_stage_model, idata_first_stage = make_reg_model(covariate_df, target="education")

现在我们可以提取第一阶段模型的参数估计，显示出一个正的 beta 系数 beta_z[nearcollege_indicator]，这表明我们的工具变量确实有助于预测教育指标。

az.summary(idata_first_stage, var_names=["beta_z"])[
    ["mean", "sd", "hdi_3%", "hdi_97%", "r_hat"]
]

在这里，我们应用类似于 F 检验的逻辑来评估我们的模型所解释的变异比例与仅包含截距（或零模型）的模型相比是否有显著差异。这个检验之所以成立，是因为当我们的模型所暗示的 F 统计量是基于零假设的 F 分布中的极端异常值时，我们可以拒绝零假设。

def make_F_test(idata, df, df_mod=6, df_res=3003, param="beta_z"):
    covariates = list(idata["posterior"][param]["covariates"].values)
    temp = df.copy()
    temp["Intercept"] = 1
    covariate_df = temp[covariates]
    params = az.extract(idata["posterior"][param])[param]
    preds = np.dot(covariate_df, params.mean(axis=1))
    resids = df["education"] - preds
    ss_res = np.sum(np.square(resids))
    ## Intercept Only model
    ss_tot = np.sum((df["education"] - np.mean(df["education"])) ** 2)
    ss_mod = ss_tot - ss_res
    ss_mod = ss_mod / df_mod
    ss_res = ss_res / df_res
    F_stat = ss_mod / ss_res
    p_value = np.round(scipy.stats.f.sf(F_stat, df_mod, df_res), 5)
    test = {
        True: "less than the Alpha value of 0.05 meaning we can reject the Null Hypothesis",
        False: "greater than the Alpha value of 0.05 meaning we cannot reject the Null Hypothesis",
    }
    result = f"""The F-Stat: {F_stat} implies a p-values of {p_value} \nwhich is {test[p_value < 0.05]}"""
    print(result)
    return F_stat, p_value


F_stat, p_value = make_F_test(idata_first_stage, df)

The F-Stat: 133.703914026679 implies a p-values of 0.0 
which is less than the Alpha value of 0.05 meaning we can reject the Null Hypothesis

这至少表明了一些证据，即工具变量加上控制变量与处理变量相关。您可能会担心解释的变异中有太多是由于控制变量造成的，因此我们可以单独提取工具变量，并对单个工具变量执行同样的检验。

first_stage_model_instrument, idata_first_stage_instrument = make_reg_model(
    covariate_df[["Intercept", "nearcollege_indicator"]], target="education"
)
az.summary(idata_first_stage_instrument, var_names=["beta_z"])[
    ["mean", "sd", "hdi_3%", "hdi_97%", "r_hat"]
]

F_stat, p_value = make_F_test(idata_first_stage_instrument, df)

The F-Stat: 10.451566987162975 implies a p-values of 0.0

在这里，我们仍然可以在任何水平上拒绝零假设。文献中基于模拟研究提出了各种阈值，用于判断何时我们认为工具变量“足够强”，但这些都是粗略的经验法则，不应依赖它们来辩护你的设计。我们希望通过理论来证明我们的工具变量，并确保它在第一阶段回归中得到了良好识别。在贝叶斯框架下，我们想知道先验/后验更新步骤如何改变了我们的分布。

az.plot_dist_comparison(
    idata_first_stage,
    var_names=["beta_z"],
    coords={"covariates": ["nearcollege_indicator"]},
    figsize=(8, 4),
);

贝叶斯更新通过将数据中的证据纳入考虑，有意义地改变了从先验实现到后验分布的转变。我们将继续构建包括这些人人口统计控制变量的模型，但关注的重点仍然是工具变量及其合理性。

减少形式

虽然在相关性假设的情况下，我们可以像上面那样努力证明我们的假设，但对于排除限制假设，我们无法直接测试这个主张。记住我们正在论证的是，我们的工具变量 `Z` 对结果 `Y` 的影响仅通过处理变量 `T` 实现。在我们的情况下，这意味着如果大学的接近程度对收入的影响独立于学校教育的中介效应，则排除限制假设会被违反。这并非不可能。为了使 `nearcollege` 通过与教育无关的机制影响工资，必须假定存在某种工资增长或下降的动力，这种动力受到教育接近程度的影响。你可能会认为商业投资的密度集中在教育中心附近，因此即使你没有上过大学，没有直接接受教育，你的工资也会受到当地同伴教育水平的影响。这种情况可以接受吗？这是否符合大多数城市中心内的大学分布情况？在评估我们的工具变量和结果之间的关系时，我们需要关注这些问题。

接下来，我们拟合减少形式模型，该模型基于我们的控制变量和工具变量来预测结果，但移除了处理变量。

reduced_form_model, idata_reduced_form = make_reg_model(covariate_df, target="log_wage")
az.summary(idata_reduced_form["posterior"]["beta_z"])[
    ["mean", "sd", "hdi_3%", "hdi_97%", "r_hat"]
]

az.plot_dist_comparison(
    idata_reduced_form,
    var_names=["beta_z"],
    coords={"covariates": ["nearcollege_indicator"]},
    figsize=(8, 4),
);

在这里，减少形式回归指示了一个接近于零的参数估计，这表明至少在预测对数工资的任务中，当我们将教育变量包含在我们的工具变量模型中时，我们应该看到接近度的预测影响被吸收到了教育变量中。

Wald 估计

这些探索性回归让我们理解了工具变量设置中的动态，但我们还可以（免费地）通过 Wald 估计程序计算 LATE 量。

reduced = az.extract(idata_reduced_form["posterior"])["beta_z"].sel(
    covariates="nearcollege_indicator"
)
first = az.extract(idata_first_stage["posterior"])["beta_z"].sel(
    covariates="nearcollege_indicator"
)
estimate = reduced / first
estimate.mean().data