假设有三个烤饼：一个烤饼两面都进行了烤制；一个只烤制了一面；一个两面都没有进行烤制。现在你买了一个烤饼，服务员给你放到了餐桌上。此时，你发现你看到的烤饼的一面是烤制过的，那么此时你猜测它的另一面也是烤制的概率是多少？
可能有不少人会回答：1/2。原因很简单，因为3个烤饼的6个面，有一半是烤制过的，另一半是没有烤制的，所以概率是1/2。
但是实际上，这儿犯了一个错误，忽略了现有的条件，因为你已经知道了1个面是烤制过的，所以该事件的全集不是所有可能性，而是在新条件成立的情况下的所有事件发生的可能性。这也是著名的Bertrand's box悖论。其实也是条件概率的一个例子：
所以我们上面的例子可以写成：
分子是两面同时烤制的概率，即1/3； 分母是上述三种情况出现B面烤制的概率，即(1/3)*1+(1/3)*0.5+(1/3)*0 = 1/2
所以，最终我们得到的概率是：
如果你不用条件概率来解释，可以采用枚举法，把所有事件发生的可能路径都罗列出来，然后数一数各条件是发生路径的多少即可。概率理论并不是很晦涩的东西，本质上就是计数，数一数各个事件发生方式/路径的多少而已。
贝叶斯方法中大量使用了条件概率，不仅仅包括观测数据，还包括在模型参数中。我们建立一个模型，提供一些数据，剩下的则交给条件概率去做就可以了。也就是”假设-演绎推断“的逻辑思路。
下面进入正题，本文探讨”假设-演绎推断“策略在常用的两种情况：一种是在模型中纳入测量误差；另一种是通过贝叶斯归因进行缺失数据的估计。
一、测量误差
之前文章中介绍过美国各个州的离婚率和结婚率的例子。该数据有50条观测（即50个州），每一个州有一个调查统计数据（行数据），包括人口数、中位结婚年龄、平均结婚率、结婚率的标准误、平均离婚率、离婚率的标准误等。
在之前的例子中，我们用来说明各个因变量之间的虚假相关性和如何通过多元回归来识别校正。但是我们当时使用的是结婚率和离婚率的平均值（因为数据中没有具体各个观测值，只有平均值），但是我们的数据中还提供了两个变量的标准误差。显然，如果一个模型仅仅使用点平均值来代替整个分布，而没有考虑标准误的话，数据所提供的信息并没有得到充分的挖掘。
所以，这儿在我们的新模型中，不仅考虑了数据的均值，还要考虑数据的标准误差。假设离婚观测数据来自一个正态分布，该正态分布有两个参数：均值
和标准差
也就是每一个州，离婚率的观测数据（这儿的观测数据实际上数据集中给出的各个州离婚率的均值）是来自不同的正态分布的。因为真实的分布标准差无法得知，我们可以使用数据中提供的测量误差来代替（在没有其他信息证据的前提下，这种代替能够使得熵值极大化）。
通过上述方法将测量误差纳入模型之后，我们再来探究结婚中位年龄(A)、结婚率(R)这两个自变量和离婚率(D)的关系。
其中蓝色部分是我们加入了标准误差之后导致的模型变化。最大的区别是模型的因变量变成了一个分布的参数！测量的不确定性会影响模型线性部分的参数，而这些线性部分的参数也会影响测量的不确定性。
实现上述模型如下：
library(rethinking)
data(WaffleDivorce)
d <- WaffleDivorce
dlist <- list(
div_obs = d$Divorce,
div_sd = d$Divorce.SE,
R = d$Marriage,
A = d$MedianAgeMarriage
)
m14.1 <- map2stan(
alist(
div_est ~ dnorm(mu, sigma),
mu <- a + bA*A + bR*R,
div_obs ~ dnorm(div_est, div_sd),
c(a,bA,bR) ~ dnorm(0,10),
sigma ~ dcauchy(0,0.25)
),
data = dlist,
start = list(div_est=dlist$div_obs),
WAIC=FALSE,iter=5000,warmup = 1000, chains=2, cores=2,
control=list(adapt_delta=0.95)
)
precis(m14.1, depth=2)
输出结果：
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## div_est[1] 11.72 0.66 10.68 12.78 6178 1
## div_est[2] 11.16 1.03 9.54 12.78 9604 1
## div_est[3] 10.44 0.61 9.43 11.36 9045 1
## ...
## div_est[49] 8.47 0.52 7.69 9.35 8511 1
## div_est[50] 11.56 1.06 9.87 13.22 6578 1
## a 21.65 6.47 10.87 31.54 2902 1
## bA -0.56 0.21 -0.89 -0.22 2924 1
## bR 0.13 0.07 0.01 0.24 3623 1
## sigma 1.07 0.20 0.74 1.37 1975 1
上述模型我们一共得到54个参数，其中50个是属于各个州的离婚率均值，另外4个是：a, bA, bR 和sigma。
在之前的模型中，我们没有在模型中包含测量误差，估计的中位结婚年龄的斜率bA大小约为-1，这儿的斜率明显减小了。也就是说包含了测量误差之后，模型预测的离婚率和结婚年龄之间的相关关系减弱了。这不是一个个例，普遍来说，如果模型忽略测量误差的存在，会夸大自变量和因变量之间的相关性。
在本例中，结婚年龄很小和很大的这些极端州的测量误差也相对较大，所以对这些州拟合的模型有较大的灵活度，使得模型这些州的预测值可以与观察值存在较大差异。也就我们之前说的模型收缩shrinkage。
如上图左，每个点表示一个州的离婚率的测量误差（横坐标）和预测值与观测值之间的差异（纵坐标），可以看出，测量误差越大的州，预测和观测之间的差异越大。
上图右，表示自变量结婚年龄和离婚率之间的关系，灰色是之前没有考虑测量误差的模型作出的预测，蓝色是本部分在模型中加入测量误差的模型作出的模型预测。可以看到，加入测量误差之后，模型斜率变小了。对于测量误差比较大的州，模型的斜率很大程度上影响了这些州的预测值；对于测量误差比较小的州，观测数据很大程度上影响了模型的斜率。
换句话说，相当于给各个州的观测值增加了一个权重，测量误差小的州，权重系数比较大；测量误差大的州，权重系数比较小。
自变量和因变量都存在测量误差
上面的例子我们考虑了因变量离婚率存在测量误差，而实际上自变量也可以存在测量误差。即每一个自变量的观测数据都是来自一个分布总体的抽样。比如上述数据中，因变量结婚率同样是一个含有测量误差的变量。为了把结婚率的测量误差也纳入模型中，我们建立了一下模型：
其中模型蓝色部分就是加入的自变量结婚率的测量误差。实现上述模型如下：
dlist <- list(
div_obs = d$Divorce,
div_sd = d$Divorce.SE,
mar_obs = d$Marriage,
mar_sd = d$Marriage.SE,
A = d$MedianAgeMarriage
)
m14.2 <- map2stan(
alist(
div_est ~ dnorm(mu, sigma),
mu <- a + bA*A + bR*mar_est[i],
div_obs ~ dnorm(div_est, div_sd),
mar_obs ~ dnorm(mar_est, mar_sd),
c(a,bA,bR) ~ dnorm(0,10),
sigma ~ dcauchy(0,2.5)
),
data=dlist,
start = list(div_est = dlist$div_obs,mar_est = dlist$mar_obs),
WAIC=FALSE, iter=5000, warmup=1000, chains=3, cores=3,
control=list(adapt_delta=0.95)
)
precis(m14.2, depth=2)
输出结果：
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## div_est[1] 11.78 0.68 10.66 12.84 14368 1
## div_est[2] 11.28 1.10 9.55 13.05 13099 1
## div_est[3] 10.47 0.62 9.51 11.49 19656 1
## ...
## div_est[49] 8.46 0.50 7.64 9.26 18487 1
## div_est[50] 11.51 1.16 9.68 13.36 11946 1
## mar_est[1] 20.54 1.28 18.48 22.55 17740 1
## mar_est[2] 26.29 2.90 21.60 30.87 20871 1
## ...
## mar_est[49] 17.15 0.77 15.96 18.42 22794 1
## mar_est[50] 29.82 3.77 23.62 35.63 16516 1
## a 20.83 6.69 10.16 31.31 5307 1
## bA -0.54 0.21 -0.87 -0.19 5909 1
## bR 0.14 0.08 0.01 0.28 4848 1
## sigma 1.10 0.21 0.78 1.44 4273 1
不过和上一个模型比较，你会发现，加入自变量结婚率的测量误差之后，模型的参数好像变化不大，对模型的预测和推断没有太大影响。 当然，上述模型中还有一个自变量是结婚年龄，实际上它也存在测量误差，只是我们的数据集中并没有给出误差大小。如果数据集给出了每个州结婚年龄的测量误差，或者其分布的参数，我们同样可以把结婚年龄的测量误差加入到模型中。 和前面章节不同，这部分的数据是含有“观测均数”和标准误差的数据，如果我们仅仅使用“观测均数”，会导致数据中信息没有被充分挖掘。这也提示我们，平时在建模、分析数据的时候，一定要注意当你使用均值的时候，你其实把一个数据压缩成一个点数据。而理想的做法，是用分布来描述一组数据，而不仅仅是一个点！ 所以：不要均值！
二、缺失数据
上面的一些例子我们考虑的是存在测量误差的情况，此时，用分布来代替观测，给模型一个足够的灵活度，使得各个不同的观测之间可以形成“信息流”，通过回归模型，不同观测之间相互影响，相互“借鉴学习”。 但是有些时候，我们的数据中存在很多缺失值，使得观测不完整。这时应该怎么做？直接把这些含有缺失值的观测都删掉吗？但是，删掉这些观测会浪费掉很多信息。其实和前面的测量误差一样，各个观测之间的“信息流”可以相互影响，可以通过现有完整观测的数据来推测缺失观测的数据。
大脑新皮质层归因模型
之前的文章中还介绍过一个灵长类动物乳汁中含有能量多少的例子，模型中以不同灵长类物种乳汁能量密度为因变量，大脑新皮质层的百分比和体型大小为两个自变量。一共29个观测，但是其中12个观测是缺失数据（大脑新皮质层百分比数据缺失），但是我们的做法是使用complete.cases()函数，剔除了含有缺失数据的观测，最后只有17个观测数据用在了模型中。
现在，我们不剔除缺失数据，而是使用贝叶斯归因对这些数据进行预测。当然，这些预测得到的其实是一个后验概率分布。
需要说明的是，这儿的缺失数据我们认为是随机缺失，而不是系统性存在偏倚的缺失。比如，如果有些体型小的灵长类动物，其大脑皮层数据可能不容易测得，而造成体型小的动物的缺失数据较多。这种缺失数据不是本研究的范畴。
含有缺失数据的模型实际上要同时做两件事，对缺失数据建模，同时还要对结果变量建模。已存在的数据能够为缺失数据提供一个先验分布，然后根据自变量和因变量的关系对缺失数据进行预测。所以每一个缺失数据预测得到的都是一个后验概率分布。
在灵长类乳汁能量数据中，大脑新皮质变量存在12个缺失值。下面我们是建立了一个包含12个缺失数据的模型，来研究乳汁能量密度和大脑皮质百分比以及体型大小的关系。
实际上，对于大脑新皮质变量，其变成了一个既含有数据又含有参数的变量。如下：
其中
表示缺失值，每一个缺失值都会有一个后验分布。
上面模型中，蓝色部分是对含有缺失值的大脑新皮质层百分比N的建模。当
不是缺失值，而是观测数据的时候，模型第三行就是该观测数据的似然函数v，认为该数据来自均值为，标准差为
的正态分布；而如果
是缺失值的时候，第三行就被解释为一个先验分布。
实现上述模型如下：
library(rethinking)
data(milk)
d <- milk
d$neocortex.prop <- d$neocortex.perc/100
d$logmass <- log(d$mass)
data_list <- list(
kcal = d$kcal.per.g,
neocortex = d$neocortex.prop,
logmass = d$logmass
)
m14.3 <- map2stan(
alist(
kcal ~ dnorm(mu, sigma),
mu <- a + bN*neocortex + bM*logmass,
neocortex ~ dnorm(nu, sigma_N),
a ~ dnorm(0,100),
c(bN,bM) ~ dnorm(0,10),
nu ~ dnorm(0.5,1),
c(sigma, sigma_N) ~ dcauchy(0,1)
),
data=data_list, iter=1e4, chains=2
)
precis(m14.3, depth=2)
输出结果：
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## neocortex_impute[1] 0.63 0.05 0.55 0.71 8789 1
## neocortex_impute[2] 0.62 0.05 0.54 0.70 7923 1
## neocortex_impute[3] 0.62 0.05 0.54 0.70 7562 1
## neocortex_impute[4] 0.65 0.05 0.58 0.73 9213 1
## neocortex_impute[5] 0.70 0.05 0.62 0.78 10238 1
## neocortex_impute[6] 0.66 0.05 0.58 0.73 8556 1
## neocortex_impute[7] 0.69 0.05 0.61 0.76 10760 1
## neocortex_impute[8] 0.70 0.05 0.62 0.77 9855 1
## neocortex_impute[9] 0.71 0.05 0.64 0.79 8440 1
## neocortex_impute[10] 0.64 0.05 0.57 0.73 9230 1
## neocortex_impute[11] 0.66 0.05 0.58 0.73 9130 1
## neocortex_impute[12] 0.69 0.05 0.62 0.77 8721 1
## a -0.54 0.47 -1.31 0.19 2456 1
## bN 1.92 0.73 0.80 3.13 2421 1
## bM -0.07 0.02 -0.11 -0.04 3319 1
## nu 0.67 0.01 0.65 0.69 6945 1
## sigma 0.13 0.02 0.10 0.17 3995 1
## sigma_N 0.06 0.01 0.04 0.08 4837 1
通过上述的结果可以看到，这12个缺失值，每一个缺失值都相当于一个参数，有各自的估计值。 下面是我们之前把缺失值剔除掉之后拟合模型的做法，我们可以比较一下剔除缺失值和不剔除缺失值的差别。
dcc <- d[complete.cases(d$neocortex.prop),]
data_list_cc <- list(
kcal = dcc$kcal.per.g,
neocortex = dcc$neocortex.prop,
logmass = dcc$logmass
)
m14.3cc <- map2stan(
alist(
kcal ~ dnorm(mu, sigma),
mu <- a + bN*neocortex + bM*logmass,
a ~ dnorm(0,100),
c(bN,bM) ~ dnorm(0,10),
sigma ~ dcauchy(0,1)
),
data = data_list_cc, iter=1e4, chains=2
)
precis(m14.3cc)
输出结果:
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a -1.10 0.59 -2.03 -0.15 2205 1
## bN 2.82 0.92 1.36 4.28 2168 1
## bM -0.10 0.03 -0.14 -0.05 2583 1
## sigma 0.14 0.03 0.09 0.18 3345 1
比较两个模型中bN和bM的参数估计值发现，前面保留缺失值后的模型，模型估计的斜率都变小了。这个结果似乎让人很失望，自变量和因变量之间的关系变得不那么明显了。为什么会出现这种情况？
我们通过散点图来看一下缺失值的预测情况：
上图中空心点是缺失值的预测，黑线表示其89%后验分布区间。蓝色实心点表示非缺失的观测值。左图是自变量新皮质百分比和因变量乳汁能量密度的关系，我们可以发现虽然缺失值的预测分布区间很宽，但是整体上还是有向回归线靠拢的趋势，毕竟我们是通过非缺失值来引导估计缺失值的。
右图是两个自变量之间的关系。我们可以发现非缺失值（实心点）中两个自变量存在正相关关系，而缺失值中，两个变量的却看不出相关关系。这也就是问题所在，在上述模型中，我们虽然考虑到了自变量自身和因变量的关系，但是却忽略了两个自变量之间存在的相关关系！贝叶斯模型是逻辑，不是魔法。模型中没有加入两个自变量之间的关系，最后的后验概率分布中，自然也不会体现出他们之间的相关关系。
所以，我们需要对模型进行改进。
改进的归因模型
这儿我们要在模型中加入两个变量之间的相关关系。做法是把上述模型中的第三行：
改写成：
这儿
，
和用来描述两个两个自变量之间的线性相关关系。
m14.4 <- map2stan(
alist(
kcal ~ dnorm(mu, sigma),
mu <- a + bN*neocortex + bM*logmass,
neocortex ~ dnorm(nu, sigma_N),
nu <- a_N + gM*logmass,
a ~ dnorm(0,100),
c(bN,bM,gM) ~ dnorm(0,10),
a_N ~ dnorm(0.5,1),
sigma_N ~ dcauchy(0,1),
sigma ~ dcauchy(0,1)
),
data=data_list, iter=1e4, chains=2
)
precis(m14.4, depth=2)
输出结果
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## neocortex_impute[1] 0.63 0.04 0.58 0.69 7983 1
## neocortex_impute[2] 0.63 0.04 0.57 0.69 9208 1
## neocortex_impute[3] 0.62 0.04 0.56 0.68 8807 1
## neocortex_impute[4] 0.65 0.03 0.59 0.70 9442 1
## neocortex_impute[5] 0.66 0.04 0.61 0.72 7972 1
## neocortex_impute[6] 0.63 0.04 0.57 0.68 8621 1
## neocortex_impute[7] 0.68 0.03 0.62 0.73 9925 1
## neocortex_impute[8] 0.70 0.03 0.64 0.75 9691 1
## neocortex_impute[9] 0.71 0.03 0.65 0.76 9840 1
## neocortex_impute[10] 0.66 0.04 0.61 0.72 8718 1
## neocortex_impute[11] 0.68 0.03 0.62 0.73 8996 1
## neocortex_impute[12] 0.74 0.04 0.68 0.80 9325 1
## a -0.86 0.48 -1.63 -0.11 2749 1
## bN 2.42 0.75 1.21 3.60 2719 1
## bM -0.09 0.02 -0.12 -0.05 3608 1
## gM 0.02 0.01 0.01 0.03 6208 1
## a_N 0.64 0.01 0.62 0.66 5024 1
## sigma_N 0.04 0.01 0.03 0.05 5040 1
## sigma 0.13 0.02 0.09 0.16 5067 1
从拟合结果来看，模型中的两个自变量成正相关gM。而此时两个bN和bM两个变量和因变量的关系也变得更为显著了。我们再看一下缺失值的预测情况。
左图看出，各个缺失值的预测分布区间变窄了（横线变短了）。这也很好理解，因为我们考虑了缺失值和另一个自变量的相关关系，所以能够获得更多与该缺失变量的信息，估计的分布区间自然会变窄。
右图更明显的表现出了缺失变量预测值和另一个自变量的相关关系，这和我们前一个拟合的模型有明显的区别。此时模型对缺失值的预测要远比之前的模型更加准确。
三、总结
本章讨论了两种特殊数据形式的贝叶斯建模，第一种数据并非原始数据，而是经过汇总的，包括均值和标准误差。既往我们仅仅对均值感兴趣，只将均值纳入模型，而忽略了标准误差。本章则将标准误差也纳入到模型中。 第二种是数据中存在缺失值，既往我们可能会直接将含有缺失值的观测删除，但是本章保留了这些缺失值。含有缺失值的变量在模型中既是参数、又是数据。通过模型我们可以估计出各个缺失值的后验分布。对这两种数据形式的建模本质上都是尽最大可能压榨数据中的信息。 最后最后，还是要提醒一下，对一个变量或一组数据的正确描述应该是用分布来描述，不应该仅仅用一个均值，这会丢失掉太多信息！应该摒弃这一传统做法。正如全国人均收入xxx元类似的说法，不应该出现在学术资料中。
【谢谢阅读，欢迎转发分享】
资料来源：《statistical rethinking》
