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

目录

使用 pymc 模型的工具变量建模 (IV)

使用 pymc 模型的工具变量建模 (IV)

这份笔记展示了一个使用工具变量模型（Instrumental Variable, IV）的例子。我们将会遵循 Acemoglu, Johnson 和 Robinson (2001) 的一个案例研究，该研究尝试解开强大的政治机构对于以国内生产总值（GDP）衡量的经济生产力的影响。本示例借鉴了Hansen的《计量经济学》以及Acemoglu等人关于“殖民主义起源、繁荣及持久影响”的讨论。

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from matplotlib.lines import Line2D
from sklearn.linear_model import LinearRegression as sk_lin_reg

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

这项工作的目的是考察各个被殖民国家在以对数GDP形式的不同结果，作为政治制度稳健性的一种度量函数，这里的稳健性度量为：风险。这个变量是对抗财产剥夺的法律保护程度的评分。评分越高，表示可用的法律保护越多。在实验情境中，你可能会期待处理变量是二元的，但在这里它是连续的，就像剂量一样。想法是法律和政治制度结构具有长期的好处，因为各个殖民地在不同的殖民力量下受到不同的对待（具有不同程度的制度发展），所以论点是我们有可能利用这类数据来估计政治制度对GDP的影响。

用一张图来表示：

在这里，协变量包含一个处理变量（risk），它通过未测量的历史（unmeasured history）这一中介值与结果（log-gdp）相关联，违反了OLS假设中的独立性。想法是通过引入我们的工具变量logmort0来恢复无偏的处理效应，该变量假设只通过其对（risk）的影响与（unmeasured history）相关。

这种相对简单的图在政策采纳的情境中相当常见，其中我们关心的是某个处理/政策（）对人口的影响，其中（）是一些影响采纳的度量。在医学情境中，第一阶段（）被称为“意向治疗”回归。由于这种与政策实施的关系，工具变量回归在行业中往往非常有用。

我们现在将具体展示如何在 CausalPy 中估计这类回归，以便在这种由DAG特征化的情境中恢复准确的参数。

N = 100
e1 = np.random.normal(0, 3, N)
e2 = np.random.normal(0, 1, N)
Z = np.random.uniform(0, 1, N)
## Ensure the endogeneity of the the treatment variable
X = -1 + 4 * Z + e2 + 2 * e1
y = 2 + 3 * X + 3 * e1

test_data = pd.DataFrame({"y": y, "X": X, "Z": Z})

sample_kwargs = {
    "tune": 1000,
    "draws": 2000,
    "chains": 4,
    "cores": 4,
    "target_accept": 0.99,
}
instruments_formula = "X  ~ 1 + Z"
formula = "y ~  1 + X"
instruments_data = test_data[["X", "Z"]]
data = test_data[["y", "X"]]
iv = InstrumentalVariable(
    instruments_data=instruments_data,
    data=data,
    instruments_formula=instruments_formula,
    formula=formula,
    model=InstrumentalVariableRegression(sample_kwargs=sample_kwargs),
)

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

我们可以看到，beta_z 参数针对 X 变量接近真实值 3。