机器算法验证 - 从识别到估计 - 吾爱随笔录

从识别到估计

机器算法验证估计参考因果关系

2022-03-11 02:47:11

我目前正在阅读 Pearl 关于因果关系的文章（Pearl，2009 年，第 2 版），并努力建立模型的非参数识别与实际估计之间的联系。不幸的是，珀尔本人对这个话题非常沉默。

举个例子，我有一个带有因果路径的简单模型， $x \rightarrow z \rightarrow y$ ，以及影响所有变量的混杂因素 $w \rightarrow x$ , $w \rightarrow z$ 和 $w \rightarrow y$ . 此外， $x$ 和 $y$ 与未观察到的影响有关， $x \leftarrow \rightarrow y$ . 根据 do-calculus 的规则，我现在知道干预后（离散）概率分布由下式给出：

P (y ∣ d o (x)) = \sum_{w, z} [P (z ∣ w, x) P (w) \sum_{x} [P (y ∣ w, x, z) P (x ∣ w)]] .

$P(y \mid do(x)) = \sum_{w,z}\bigl[P(z\mid w,x)P(w)\sum_{x}\bigl[P(y\mid w,x,z)P(x\mid w)\bigr]\bigr].$

我想知道如何估计这个数量（非参数或通过引入参数假设）？尤其是当 $w$ 是一组几个混杂变量，感兴趣的数量是连续的。在这种情况下，估计数据的联合干预前分布似乎非常不切实际。有人知道处理这些问题的珍珠方法的应用吗？我会很高兴得到一个指针。

1个回答

这个问题问得好。首先让我们验证您的公式是否正确。您提供的信息对应于以下因果模型：

正如你所说，我们可以推导出 $P(Y|do(X))$ 使用do-calculus的规则。在 R 中，我们可以使用 package 轻松做到这一点causaleffect。我们首先加载igraph以使用您提出的因果图创建一个对象：

library(igraph)
g <- graph.formula(X-+Y, Y-+X, X-+Z-+Y, W-+X, W-+Z, W-+Y, simplify = FALSE)
g <- set.edge.attribute(graph = g, name = "description", index = 1:2, value = "U")

其中前两项X-+Y, Y-+X代表未观察到的混杂因素 $X$ 和 $Y$ 其余术语代表您提到的有向边。

然后我们要求我们的估计：

library(causaleffect)
cat(causal.effect("Y", "X", G = g, primes = TRUE, simp = T, expr = TRUE))

\sum_{W, Z} (\sum_{X^{'}} P (Y | W, X^{'}, Z) P (X^{'} | W)) P (Z | W, X) P (W)

$\sum_{W,Z}\left(\sum_{X^{\prime}}P(Y|W,X^{\prime},Z)P(X^{\prime}|W)\right)P(Z|W,X)P(W)$

这确实与您的公式不谋而合——带有观察到的混杂因素的前门案例。

现在让我们进入估计部分。如果您假设线性（和正态性），事情就会大大简化。基本上你想要做的是估计路径的系数 $X \rightarrow Z \rightarrow Y$ .

让我们模拟一些数据：

set.seed(1)
n <- 1e3
u <- rnorm(n) # y -> x unobserved confounder
w <- rnorm(n)
x <- w + u + rnorm(n)
z <- 3*x + 5*w + rnorm(n)
y <- 7*z + 11*w + 13*u + rnorm(n)

注意在我们的模拟中，变化的真实因果效应 $X$ 在 $Y$ 是 21。您可以通过运行两个回归来估计这一点。第一的 $Y \sim Z+W+X$ 获得效果 $Z$ 在 $Y$ 进而 $Z \sim X + W$ 获得效果 $X$ 在 $Z$ . 您的估计将是两个系数的乘积：

yz_model <- lm(y ~ z + w + x)
zx_model <- lm(z ~ x + w)

yz <- coef(yz_model)[2]
zx <- coef(zx_model)[2]
effect <- zx*yz
effect
       x 
21.37626

并且为了推断，您可以计算产品的（渐近）标准误差：

se_yz <- coef(summary(yz_model))[2, 2]
se_zx <- coef(summary(zx_model))[2, 2]
se <- sqrt(yz^2*se_zx^2 + zx^2*se_yz^2)

您可以将其用于测试或置信区间：

c(effect - 1.96*se, effect + 1.96*se) # 95% CI
       x        x 
19.66441 23.08811

您还可以执行（非/半）参数估计，稍后我将尝试更新此答案，包括其他程序。

其它你可能感兴趣的问题

上一篇比例数据 - beta 分布与具有二项分布和 logit 链接的 GLM 下一篇直观地总结一堆有向线段