【Note】 因果推断(五)——基于模型的检验

在前面的三节中，Fisher 和 Neyman 的方法都将潜在结果 \(Y_i​\) 视为定值，但只是观测到了一部分即 \(Y_i^\mathrm{obs}​\)。随机性主要来自于分配机制 \(W_i​\) ；如果考虑样本是从总体抽样得到的话，这一步也会引入一定的随机性。而在回归模型以及这章中，我们将会把 \(Y_i​\) 同样视为随机变量——即使在考虑有限样本下依然如此。后面我们会看到，基于模型的方法1. 带来了相当大的灵活性——我们可以和之前一样给出平均组间效应，我们还可以给出分位数、给出 \(Y_i​\) 的方差等等。事实上，我们可以对任何估计量 \(\tau=\tau(\mathbf{Y}(0),\mathbf{Y}(1))​\) ，甚至 \(\tau=\tau(\mathbf{Y}(0),\mathbf{Y}(1),\mathbf{W},\mathbf{X})​\) 给出推断(当然，我们会要求这个 \(\tau​\) 是行可交换的)。 2. 可以很容易的扩展到总体样本的情况。3. 可以在观察实验(也就是我们无法对分配机制 \(W_i​\) 做干预) 下发挥作用。

基于模型的算法的核心只有两步

1. 给出 \(f(\mathbf{Y}^\mathrm{mis}\mid \mathbf{Y}^\mathrm{obs},\mathbf{W})\)
2. 推导 \(\tau=\tau(\mathbf{Y}(0),\mathbf{Y}(1),\mathbf{W})\) 的分布并进行推断

数据

序号	\(Y_i(0)\)	\(Y_i(1)\)	\(W_i\)	\(Y_i^\mathrm{obs}\)
1	0	?	0	0
2	?	9.9	1	9.9
3	12.4	?	0	12.4
4	?	3.6	1	3.6
5	0	?	0	0
6	?	24.9	1	24.9

我们关心的是

\[\begin{align}\tau_\mathrm{fs}&=\tau(\mathbf{Y}(0),\mathbf{Y}(1))\\&=\frac{1}{6}\sum_{i=1}^{6}\big(Y_i(1)-Y_i(0)\big)\\&=\frac{1}{6}\sum_{i=1}^{6}\big((2\cdot W_i-1)(Y_i^\mathrm{obs}-Y_i^\mathrm{mis}\big)\\&=\tilde\tau(\mathbf{Y}^\mathrm{obs},\mathbf{Y}^\mathrm{mis},\mathbf{W})\end{align}\]

由于 \(\mathbf{Y}^\mathrm{mis}\) 的缺失，我们只能给出估计值，此时为

\(\widehat{\tau}=\tilde{\tau}(\mathbf{Y}^\mathrm{obs},\widehat{\mathbf{Y}}^\mathrm{mis},\mathbf{W})\)

因此核心就在于给出 \(\widehat{\mathbf{Y}}^\mathrm{mis}\)

朴素模型一

使用各组观测到的平均值来替换缺失值，即

\[\mathbb{P}[Y_i^\mathrm{mis}=y\mid\mathbf{Y}^\mathrm{obs},\mathbf{W}]=\left\{\begin{array}\\1&W_i=0,\ y=12.8\\1&W_i=1,\ y=4.1\\0&\mathrm{o.w.}\end{array}\right.\]

此时给出的估计 \(\widehat{\tau}=12.8-4.1=8.7=\bar{Y}_t^\mathrm{obs}-\bar{Y}_c^\mathrm{obs}=\widehat{\tau}^\mathrm{dif}\)

这个方法的劣势在于，由于完全没有随机性，因此我们只能给出点估计而无法给出精确度，正如 Fisher 的方法一样。而且这样的赋值方法也没有保存下来 \(Y_i\) 的方差。比如 \(Y_i(1)\) 中我们已经观测到了 3.9,9.9,24.9，但我们的赋值湮灭了这一点。

朴素模型二

对于每一个缺失值，从已观测到的值中抽样，以本数据为例子就是

\[\mathbb{P}[Y_i^\mathrm{mis}=y\mid\mathbf{Y}^\mathrm{obs},\mathbf{W}]=\left\{\begin{array}\\1/3&W_i=0,\ y\in\{3.6,9.9,24.9\}\\1/3&W_i=1,\ y=12.4\\2/3&W_i=1,\ y=0\\0&\mathrm{o.w.}\end{array}\right.\]

显然此时有 \(3^6=729\) 种赋值方式，使用代码跑一下得到均值 8.7，标准差 3.1

这一模型比上一模型复杂，但可以给我们一个区间估计。但他并没有囊括所有的随机性——理论上我们应该从 \(Y_i\) 的精确分布中抽样然后填入缺失值，但我们这里只用了观测到的数据来替代这一行为。为了将总体的信息引入，我们使用下面的贝叶斯模型。

Y0 = c(0,0,12.4)
Y1 = c(9.9,3.6,24.9)
exacts = rep(NA, 729)
for (i in 1:729) {
  id1 = i %% 3 + 1
  id2 = round(i/3) %% 3 + 1
  id3 = round(i/9) %% 3 + 1
  id4 = round(i/27) %% 3 + 1
  id5 = round(i/81) %% 3 + 1
  id6 = round(i/243) %% 3 + 1
  y0 = c(Y0, Y0[id1], Y0[id2], Y0[id3])
  y1 = c(Y1, Y1[id4], Y1[id5], Y1[id6])
  exacts[i] = mean(y1-y0)
}
mean(exacts) # 8.666667
sd(exacts) # 3.084173
# or you can just sample to get an approximate vallue
taus = rep(NA, 10000)
set.seed(42)
for (i in 1:10000) {
  y0 = c(Y0, sample(Y0, 3, replace=T))
  y1 = c(Y1, sample(Y1, 3, replace=T))
  taus[i] = mean(y1-y0)
}
mean(taus) # 8.671788
sd(taus) # 3.069502

贝叶斯模型的建立

三个输入

1. 潜在结果的联合分布 \(f(\mathbf{Y}(0),\mathbf{Y}(1))\)

   由于行可交换，其可以写成各样本的乘积

   \(f(\mathbf{Y}(0),\mathbf{Y}(1))=\int\prod_{i=1}^{N}f(Y_i(0),Y_i(1)\mid\theta)\cdot p(\theta)\,\mathrm{d}\theta\)

   因此我们需要给出第二个输入

2. \(p(\theta)\) 即参数的先验分布

3. 在观察实验中，我们可能还需要引入 \(f(\mathbf{W}\mid\mathbf{Y}(0),\mathbf{Y}(1))\)。但在完全随机试验中，我们知道其在支撑集上其为常数 \(\binom{N}{N_t}^{-1}\)，因此无需引入。

第一步：推导 \(f(\mathbf{Y}^\mathrm{mis}\mid\mathbf{Y}^\mathrm{obs},\mathbf{W},\theta)\)

首先我们有 \(f(\mathbf{Y}(0),\mathbf{Y}(1),\mathbf{W}\mid\theta)=\mathbb{P}[\mathbf{W}\mid\mathbf{Y}(0),\mathbf{Y}(1),\theta]\cdot f(\mathbf{Y}(0),\mathbf{Y}(1)\mid\theta)\)

接着我们有 \(f(\mathbf{Y}(0),\mathbf{Y}(1)\mid\mathbf{W},\theta)=\frac{f(\mathbf{Y}(0),\mathbf{Y}(1),\mathbf{W}\mid\theta)}{\mathbb{P}[\mathbf{W}\mid\theta]}=\frac{f(\mathbf{Y}(0),\mathbf{Y}(1),\mathbf{W}\mid\theta)}{\iint f(\mathbf{Y}(0),\mathbf{Y}(1),\mathbf{W}\mid\theta)\,\mathrm{d}\mathbf{Y}(0)\mathrm{d}\mathbf{Y}(1)}\)

在完全随机试验下 \(\mathbf{W}\perp\!\!\!\perp (\mathbf{Y}(0),\mathbf{Y}(1))\) 即 \(f(\mathbf{Y}(0),\mathbf{Y}(1)\mid\mathbf{W},\theta)=f(\mathbf{Y}(0),\mathbf{Y}(1)\mid\theta)\)

注意到 \((\mathbf{Y}^\mathrm{obs},\mathbf{Y}^\mathrm{mis})\) 可以写为 \((\mathbf{Y}(0),\mathbf{Y}(1),\mathbf{W})\) 的变换，因此我们可以推导

\(f(\mathbf{Y}^\mathrm{mis}\mid\mathbf{Y}^\mathrm{obs},\mathbf{W},\theta)=\frac{f(\mathbf{Y}^\mathrm{mis},\mathbf{Y}^\mathrm{obs}\mid\mathbf{W},\theta)}{f(\mathbf{Y}^\mathrm{obs}\mid\mathbf{W},\theta)}=\frac{f(\mathbf{Y}^\mathrm{mis},\mathbf{Y}^\mathrm{obs}\mid\mathbf{W},\theta)}{\int f(\mathbf{Y}^\mathrm{mis},\mathbf{Y}^\mathrm{obs}\mid\mathbf{W},\theta)\,\mathrm{d}\mathbf{y}^\mathrm{mis}}\)

这也称为 \(\mathbf{Y}^\mathrm{mis}\) 的后验预测分布

第二步：推导 \(p(\theta\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})\)

首先我们有似然函数 \(\mathcal{L}(\theta\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})\equiv f(\mathbf{Y}^\mathrm{obs},\mathbf{W}\mid\theta)=\int f(\mathbf{y}^\mathrm{mis},\mathbf{Y}^\mathrm{obs},\mathbf{W}\mid\theta)\,\mathrm{d}\mathbf{y}^\mathrm{mis}\)

于是 \(p(\theta\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})=\frac{p(\theta)\mathcal{L}(\theta\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})}{f(\mathbf{Y}^\mathrm{obs},\mathbf{W})}=\frac{p(\theta)\mathcal{L}(\theta\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})}{\int p(\theta)\mathcal{L}(\theta\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})\,\mathrm{d}\theta}\)

第三步：推导 \(f(\mathbf{Y}^\mathrm{mis}\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})\)

上面两式合并得 \(f(\mathbf{Y}^\mathrm{mis},\theta\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})=f(\mathbf{Y}^\mathrm{mis}\mid\mathbf{Y}^\mathrm{obs},\mathbf{W},\theta)\cdot p(\theta\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})\)

于是 \(f(\mathbf{Y}^\mathrm{mis}\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})=\int f(\mathbf{Y}^\mathrm{mis},\theta\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})\,\mathrm{d}\theta\)

第四步：推导 \(f(\tau\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})\)

我们知道 \(\tau=\tau(\mathbf{Y}(0),\mathbf{Y}(1),\mathbf{W})=\tilde\tau(\mathbf{Y}^\mathrm{mis},\mathbf{Y}^\mathrm{obs},\mathbf{W})\)

而给定 \(\mathbf{Y}^\mathrm{obs},\mathbf{W}\) 后，我们可以得到 \(\mathbf{Y}^\mathrm{mis}\) 的条件分布，也就可以得到 \(\tau=\tilde{\tau}\) 的条件分布

贝叶斯模型(无协变量)的例子

三个输入

首先，我们假设给定参数 \(\theta=(\mu_c,\mu_t)\) 后我们的服从二元正态分布如下

\[\begin{pmatrix}Y_i(0)\\Y_i(1)\end{pmatrix}\Bigg|\,\theta\sim\mathcal{N}\Bigg(\begin{pmatrix}\mu_c\\\mu_t\end{pmatrix},\begin{pmatrix}100&0\\0&64\end{pmatrix}\Bigg)\]

而我们的参数的先验为

\[\theta=\begin{pmatrix}\mu_c\\\mu_t\end{pmatrix}\sim\mathcal{N}\Bigg(\begin{pmatrix}0\\0\end{pmatrix},\begin{pmatrix}10000&0\\0&10000\end{pmatrix}\Bigg)\]

分配机制 \[\mathbb{P}[\mathbf{W}=\mathbf{w}\mid\mathbf{Y}(0),\mathbf{Y}(1),\mu_c,\mu_t]=\binom{N}{N_t}^{-1},\quad\sum_{i=1}^{N}\mathbf{w}_i=N_t\]

第一步：推导 \(f(\mathbf{Y}^\mathrm{mis}\mid\mathbf{Y}^\mathrm{obs},\mathbf{W},\theta)\)

我们有 \(f(\mathbf{Y}(0),\mathbf{Y}(1)\mid\mathbf{W},\theta)=f(\mathbf{Y}(0),\mathbf{Y}(1)\mid\theta)=\prod_{i=1}^{N}f(Y_i(0),Y_i(1)\mid\theta)\)

\(\begin{pmatrix}Y_i^\mathrm{mis}\\Y_i^\mathrm{obs}\end{pmatrix}\Bigg|\mathbf{W},\theta\sim\mathcal{N}\Bigg(\begin{pmatrix}W_i\cdot\mu_c+(1-W_i)\cdot\mu_t\\(1-W_i)\cdot\mu_c+W_i\cdot\mu_t\end{pmatrix},\begin{pmatrix}W_i\cdot100+(1-W_i)\cdot64&0\\0&(1-W_i)\cdot100+W_i\cdot\64\end{pmatrix} \Bigg)\)

故 \[f(\mathbf{Y}^\mathrm{mis},\mathbf{Y}^\mathrm{obs}\mid\mathbf{W},\theta)=\prod_{i=1}^{N}f(Y_i^\mathrm{mis},Y_i^\mathrm{obs}\mid\mathbf{W},\theta)\]

则 \(Y_i^\mathrm{mis}\mid\mathbf{Y}^\mathrm{obs},\mathbf{W},\theta\sim\mathcal{N}(W_i\cdot\mu_c+(1-W_i)\cdot\mu_t,W_i\cdot100+(1-W_i)\cdot64)\)

第二步：推导 \(p(\theta\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})\)

\[\begin{align}\mathcal{L}(\theta\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})&\propto\prod_{i:W_i=0}\frac{1}{\sqrt{2\pi\cdot 100}}\exp\left\{-\frac{1}{2}\left(\frac{1}{100}(Y_i^\mathrm{obs}-\mu_c)^2\right)\right\}\\&\quad\times\prod_{i:W_i=1}\frac{1}{\sqrt{2\pi\cdot 64}}\exp\left\{-\frac{1}{2}\left(\frac{1}{64}(Y_i^\mathrm{obs}-\mu_t)^2\right)\right\}\end{align}\]

于是 \[p(\theta\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})\propto\mathcal{L}(\theta\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})p(\theta)\]，合并指数项后可以推出

\[\theta\mid\mathbf{Y}^\mathrm{obs},\mathbf{W}\sim\mathcal{N}\Bigg(\begin{pmatrix}\bar{Y}_c^\mathrm{obs}\cdot\frac{N_c\cdot 10,000}{N_c\cdot 10,000+100}\\\bar{Y}_t^\mathrm{obs}\cdot\frac{N_t\cdot 10,000}{N_t\cdot 10,000+64}\end{pmatrix},\begin{pmatrix}(N_c/100+1/10000)^{-1}&0\\0&(N_t/64+1/10000)^{-1}\end{pmatrix}\Bigg)\]

代入数据得

\[\theta\mid\mathbf{Y}^\mathrm{obs},\mathbf{W}\sim\mathcal{N}\Bigg(\begin{pmatrix}4.1\\12.8\end{pmatrix},\begin{pmatrix}5.8^2&0\\0&4.6^2\end{pmatrix}\Bigg)\]

第三步：推导 \(f(\mathbf{Y}^\mathrm{mis}\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})\)

\(f(\mathbf{Y}^\mathrm{mis},\theta\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})=f(\mathbf{Y}^\mathrm{mis}\mid\mathbf{Y}^\mathrm{obs},\mathbf{W},\theta)\cdot p(\theta\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})\)

注意到后面这两项都是正态分布，可以推出 \(f(\mathbf{Y}^\mathrm{mis}\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})\) 也服从正态分布

\[\mathbb{E}[Y_i^\mathrm{mis}\mid\mathbf{Y}^\mathrm{obs},\mathbf{W},\mu_c,\mu_t]=W_i\cdot\mu_c+(1-W_i)\cdot\mu_t\]

\[\mathbb{E}[Y_i^\mathrm{mis}\mid\mathbf{Y}^\mathrm{obs},\mathbf{W}]=W_i\cdot\left(\bar{Y}_c^\mathrm{obs}\cdot\frac{N_c\cdot 10000}{N_c\cdot 10000+100}\right)+(1-W_i)\cdot\left(\bar{Y}_t^\mathrm{obs}\cdot\frac{N_t\cdot 10000}{N_t\cdot 10000+100}\right)\]

\[\begin{align}\mathbb{V}[Y_i^\mathrm{mis}\mid \mathbf{Y}^\mathrm{obs},\mathbf{W}]&=\mathbb{E}[\mathbb{V}[Y_i^\mathrm{mis}\mid\mathbf{Y}^\mathrm{obs},\mathbf{W},\theta]\mid\mathbf{Y}^\mathrm{obs},\mathbf{W}]+\mathbb{V}[\mathbb{E}[Y_i^\mathrm{mis}\mid\mathbf{Y}^\mathrm{obs},\mathbf{W},\mu_c,\mu_t]\mid\mathbf{Y}^\mathrm{obs},\mathbf{W}]\\&=W_i\cdot 100+(1-W_i)\cdot 64+W_i\cdot\frac{1}{N_c/100+1/10,000}+(1-W_i)\cdot\frac{1}{M_t/64+1/10,000}\end{align}\]

\[\begin{align}\mathbb{C}[Y_i^\mathrm{mis},Y_{i'}^\mathrm{mis}\mid\mathbf{Y}^\mathrm{obs},\mathbf{W}]&=\mathbb{E}[\mathbb{C}[Y_i^\mathrm{mis},Y_{i'}^\mathrm{mis}\mid\mathbf{Y}^\mathrm{obs},\mathbf{W},\mu_c,\mu_t]\mid\mathbf{Y}^\mathrm{obs},\mathbf{W}]\\&\quad+\mathbb{C}[\mathbb{E}[Y_i^\mathrm{mis}\mid\mathbf{Y}^\mathrm{obs},\mathbf{W},\mu_c,\mu_t],\mathbb{E}[Y_{i'}^\mathrm{mis}\mid\mathbf{Y}^\mathrm{obs},\mathbf{W},\mu_c,\mu_t]\mid\mathbf{Y}^\mathrm{obs},\mathbf{W}]\\&=0+\mathbb{C}[W_i\cdot\mu_c+(1-W_i)\cdot\mu_t,W_{i'}\cdot\mu_c+(1-W_{i'})\cdot\mu_t\mid\mathbf{Y}^\mathrm{obs},\mathbf{W}]\\&=W_i\cdot W_{i'}\cdot\frac{1}{N_c/100+1/10000}+(1-W_i)\cdot (1-W_{i'})\cdot\frac{1}{N_t/64+1/10000}\end{align}\]

利用这些信息即得分布

\[\begin{align}\left.\begin{pmatrix}Y_1^\mathrm{mis}\\Y_2^\mathrm{mis}\\Y_3^\mathrm{mis}\\Y_4^\mathrm{mis}\\Y_5^\mathrm{mis}\\Y_6^\mathrm{mis}\end{pmatrix}\right|\mathbf{Y}^\mathrm{obs},\mathbf{W}\sim\mathcal{N}\left(\begin{pmatrix}12.8\\4.1\\12.8\\4.1\\12.8\\4.1\\\end{pmatrix},\begin{pmatrix}85.3&0&21.3&0&21.3&0\\0&133.2&0&33.2&0&33.2\\21.3&0&85.3&0&21.3&0\\0&0&0&133.2&0&33.2\\21.3&0&21.3&0&85.3&0\\0&33.2&0&33.2&0&133.2\\\end{pmatrix}\right)\end{align}\]

第四步：推导 \(f(\tau\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})\)

\(\tau_\mathrm{fs}=\tau(\mathbf{Y}(0),\mathbf{Y}(1),\mathbf{W})=\frac{1}{N}\sum_{i=1}^{N}(Y_i(1)-Y_i(0))\)

\(\tau_\mathrm{fs}=\tilde{\tau}(\mathbf{Y}^\mathrm{mis},\mathbf{Y}^\mathrm{obs},\mathbf{W})=\frac{1}{N}\sum_{i=1}^{N}(1-2\cdot W_i)\cdot Y_i^\mathrm{mis}+\frac{1}{N}\sum_{i=1}^{N}(2\cdot W_i-1)\cdot Y_i^\mathrm{obs}\)

\[\begin{align}\mathbb{E}[\tau_\mathrm{fs}\mid\mathbf{Y}^\mathrm{obs},\mathbf{W}]&=\frac{N_t\bar{Y}_t^\mathrm{obs}-N_c\cdot\bar{Y}_c^\mathrm{obs}}{N}+\frac{1}{N}\sum_{i=1}^{N}(1-2\cdot W_i)\cdot\mathbb{E}[Y_i^\mathrm{obs}\mid\mathbf{Y}^\mathrm{obs},\mathbf{Y}]\\&=\bar{Y}_t^\mathrm{obs}\cdot\frac{10000N_t+64N_t/N}{10000N_t+64}-\bar{Y}_c^\mathrm{obs}\cdot\frac{10000N_c+100N_c/N}{10000N_c+100}\end{align}\]

\[\begin{align}\mathbb{V}(\tau_\mathrm{fs}\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})&=\frac{1}{N^2}\sum_{i=1}^{N}\mathbb{V}[(1-2\cdot W_i)\cdot Y_i^\mathrm{mis}\mid\mathbf{Y}^\mathrm{obs},\mathbf{W}]+\frac{1}{N^2}\sum_{i=1}^{N}\sum_{i'\neq i}\mathbb{C}[(1-2\cdot W_i)\cdot Y_i^\mathrm{mis},(1-2\cdot W_{i'})\cdot Y_{i'}^\mathrm{mis}]\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})]\end{align}\]

代入数据得 \(\tau_\mathrm{fs}\mid\mathbf{Y}^\mathrm{obs},\mathbf{W}\sim\mathcal{N}(8.7,5.2^2)\)

可以看到均值和之前保持一致，而方差比之前的模型增大了。这是十分合理的——因为我们现在的模型在估计 \(Y_i^\mathrm{mis}\) 的时候引入了额外的随机性。

数学模拟

学过贝叶斯推断的话应该知道，很多时候我们很难像上面的推导过程一样给出精确的解析解：其难点都集中在后验分布 (如 \(p(\theta\mid y), p(\tilde{y}\mid y)\)) 的导出——往往我们能给出其非归一化的概率密度，即 \(p(\theta\mid y)\propto p(\theta)p(y\mid\theta)\)，但归一化因子 \(p(y)\) 一般不好得出。事实上贝叶斯推断针对这一难点给出了相当多的解决方案——比如经典的 MCMC 方法。

仍以该题为例子说明，由于我们已经显式得得到了 \(p(theta\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})\) 和 \(f(\mathbf{Y}^\mathrm{mis}\mid \mathbf{Y}^\mathrm{obs},\mathbf{W},\mu_c,\mu_t)\)这两个后验分布，因此我们的数值模拟是十分简单和直接的：先从后验分布采样出 \(\theta_{(i)}\)，然后采样出 \(\mathbf{Y}^\mathrm{mis}_{(i)}\)，接着计算得出 \(\widehat{tau}_{(i)}\)，并根据模拟结果推得均值和方差。一份简单的模拟如下

set.seed(42)
taus = rep(NA, 10000)
y_obs_0_mean = mean(c(0, 12.4, 0))
y_obs_1_mean = mean(c(9.9, 3.6, 24.9))
tau_obs = y_obs_1_mean - y_obs_0_mean
for (i in 1:10000) {
  mu_c = rnorm(1, 4.1, 5.8)
  mu_t = rnorm(1, 12.8, 4.6)
  y_mis_0_mean = mean(rnorm(3, mu_c, 10))
  y_mis_1_mean = mean(rnorm(3, mu_t, 8))
  tau_mis = y_mis_1_mean - y_mis_0_mean
  taus[i] = (tau_obs + tau_mis) / 2
}
mean(taus) # 8.748153
sd(taus) # 5.195695

如果没有简单的形式，则采样的时候需要使用其他的方法从非归一化的概率密度中进行采样。详细的介绍在之前写过的 【Series】 贝叶斯数据分析 中的采样章节有着详细的介绍。

变形——去除潜在结果间的独立性假设

在上面的推到中，我们假设 \(Y_i(0),Y_i(1)\) 是独立(至少不相关)的，即

\(Y_i(0),Y_i(1)\mid\theta\sim\mathcal{N}\left(\begin{pmatrix}\mu_c\\\mu_t\end{pmatrix},\begin{pmatrix}\sigma_c^2&0\\0&\sigma_t^2\end{pmatrix}\right)\)

现在我们将其放宽至

\(Y_i(0),Y_i(1)\mid\theta\sim\mathcal{N}\left(\begin{pmatrix}\mu_c\\\mu_t\end{pmatrix},\begin{pmatrix}\sigma_c^2&\rho\sigma_c\sigma_t\\\rho\sigma_c\sigma_t&\sigma_t^2\end{pmatrix}\right)\)

为此，\(\theta=(\mu_c,\mu_t,\sigma_c^2,\sigma_t^2,\rho)\) ,此时 \(p(\theta)=p(\rho)\cdot p(\mu_c,\mu_t,\sigma_c^2,\sigma_t^2)\)

但如果推导一下 \(f(Y_i^\mathrm{obs}\mid\mathbf{W},\theta)\)，就会发现其和 \(\rho\) 无关，因此我们的似然函数 \(\mathcal{L}\) 也和 \(\rho\) 无关。（这是十分自然的：我们每次只拿到了 \(Y_i(0),Y_i(1)\) 的一个，凭此数据不可能得到 \(\rho\) 的有效推断——因为数据中完全不包含这一信息，事实上，贝叶斯推断里也专门针对此问题有过讨论），因此参数 \(\rho\) 的后验分布也会等于先验分布，因为数据无法提供信息进行更新。

现在考虑一个具体而极端的例子——\(\rho=1\)

\(Y_i(0),Y_i(1)\mid\theta\sim\mathcal{N}\left(\begin{pmatrix}\mu_c\\\mu_t\end{pmatrix},\begin{pmatrix}100&80\\80&64\end{pmatrix}\right)\)

可以导出

\(\begin{pmatrix}Y_i^\mathrm{mis}\\Y_i^\mathrm{obs}\end{pmatrix}\Bigg|\mathbf{W},\theta\sim\mathcal{N}\Bigg(\begin{pmatrix}W_i\cdot\mu_c+(1-W_i)\cdot\mu_t\\(1-W_i)\cdot\mu_c+W_i\cdot\mu_t\end{pmatrix},\begin{pmatrix}W_i\cdot100+(1-W_i)\cdot64&80\\80&(1-W_i)\cdot100+W_i\cdot64\end{pmatrix} \Bigg)\)

由于此时协方差不为零，导致 \(Y_i^\mathrm{mis}\) 的边缘分布发生了变化，为

\(Y_i^\mathrm{mis}\mid Y_i^\mathrm{obs},\mathbf{W},\theta\sim\mathcal{N}(W_i\cdot(\mu_c+\frac{80}{64}(Y_i^\mathrm{obs}-\mu_t)+(1-W_i)\cdot(\mu_t+\frac{80}{100}(Y_i^\mathrm{obs}-\mu_t),0)\)

这里方差为 \(0\) 可以直接推导，也可以由 \(\rho=1\) 直接给出。

注意到参数后验 \(p(\theta\mid\mathbf{Y}^\mathrm{obs},\mathbf{W})\) 不因为 \(\rho\) 的引入而改变，因此和之前的推导一致。

类似的，最后可以得到 \(Y_i^\mathrm{mis}\) 的后验期望、方差和协方差，并得到 \(\tau_\mathrm{fs}\)的方差和协方差。代入数据最后得到结果

\(\tau\mid\mathbf{Y}^\mathrm{obs},\mathbf{W}\sim\mathcal(8.7,7.7^2)\)

均值没变——这自然符合预期。方差进一步增大——这可能有点困惑。但如果读者之前看到 Neyman 里面的方差和 \(\rho_{tc}\) 的关系后，应该能够理解这一现象。

变形——引入协变量

此时模型变形为

\[f(Y_i(0),Y_i(1),X\,\mid\,\theta_{Y\mid X},\theta_X)=f(Y_i(0),Y_i(1)\mid X,\theta_{Y\mid X})\cdot f(X\mid\theta_X)\]

引入的先验为 \(p(\theta_{Y\mid X},\theta_X)=p(\theta_{Y\mid X})\cdot p(\theta_X)\)

这一独立性假设在实践中通常会被采用，但其不总是无害的——比如协变量包含之前观测的时序数据，特别是和结果 \(Y\) 高度相关的，那么我们的参数 \(\theta_X\) 可能会含有很强的控制组的分布的信息，但这一假设通常会大大简化我们的模型——我们只需 对 \(f(Y_i(0),Y_i(1)\mid X_i,\theta_{Y\mid X})\) 进行建模，仍以二元正态为例，此时为

\(Y_i(0),Y_i(1)\mid X_i,\theta\sim\mathcal{N}\left(\begin{pmatrix}X_i\beta_c\\X_i\beta_t\end{pmatrix},\begin{pmatrix}\sigma_c^2&0\\0&\sigma_t^2\end{pmatrix}\right)\)

其中 \(\theta=(\beta_c,\beta_t,\sigma_c^2,\sigma_t^2)\)

而后续的步骤完全一样。

变形——考察全体分布下的估计

之前讨论的都是有限样本下的估计，现在我们将样本视为从一个全体分布抽样得到的结果，再次考察 \(\tau_\mathrm{sp}=\mathbb{E}_\mathrm{sp}[Y_i(1)-Y_i(0)]\)

当我们的模型已经建立好时，通常 \(\tau_\mathrm{sp}\) 可以视为参数 \(\theta\) 的一个函数，即

\[\tau_\mathrm{sp}=\tau(\theta)=\mathbb{E}_\mathrm{sp}[Y_i(1)-Y_i(0)\mid\theta]=\iint (Y(1)-Y(0))f(Y(1),Y(0)\mid\theta)\,\mathrm{d}Y(1)\mathrm{d}Y(0)\]

在我们之前的模型种，则就是 \(\tau_t-\tau_c\)

如果有协变量，则 \(\tau_\mathrm{sp}=\mathbb{E}_\mathrm{sp}[\tau(\theta,\mathbf{X})]\) 其中 \(\tau(\theta,\mathbf{X})=\mathbb{E}_\mathrm{sp}[Y_i(1)-Y_i(0)\mid\mathbf{X},\theta]\)

实际计算时，我们只需从 \(\theta\) 的后验分布中采取大量样本，然后计算 \(\tau(\theta)\)，即可得到均值和方差，而无需再次从 \(f(\mathbf{Y}^\mathrm{mis}\mid\mathbf{Y}^\mathrm{obs},\mathbf{W},\theta)\) 中抽样再计算。

回到我们之前的例子，由于

\[\theta\mid\mathbf{Y}^\mathrm{obs},\mathbf{W}\sim\mathcal{N}\Bigg(\begin{pmatrix}4.1\\12.8\end{pmatrix},\begin{pmatrix}5.8^2&0\\0&4.6^2\end{pmatrix}\Bigg)\]

我们可以得到 \(\tau_\mathrm{sp}=\mu_t-\mu_c\mid\mathbf{Y}^\mathrm{obs},\mathbf{W}\sim\mathcal{N}(12.8-4.1,5.8^2+4.6^2)\sim\mathcal{N}(8.7,7.4^2)\)

同样均值一样，但是方差比独立下的方差来得大 (\(5.2^2\))，因为我们引入了样本的随机性，即使我们获得了我们拿到的样本的所有潜在结果并精确的算出了 \(\tau_\mathrm{fs}\)，我们仍然对 \(\tau_\mathrm{sp}\) 抱有一定的随机性；但比 \(\rho=1\) 的极端情况来的小 \(7.7^2\)。而 \(\rho=1\) 则相当于给出了一个最坏的、最保守的一个估计，但同时也给出了 \(\tau_\mathrm{sp}\) 的一个无偏估计。

另外值得注意的是，我们关注的量 \(\tau_\mathrm{sp}\) 不依赖于 \(\rho\)，因为似然函数不包含 \(\rho\)，则\(\tau\) 的后验分布也不会包含 \(\rho\)，如果 \(\rho\) 和 \(\mu_c,\mu_t\) 独立的话。