酷应用

不想去健身房的我，最后被贝叶斯分析说服了...

百家作者：大数据文摘 2018-07-27 15:20:39

大数据文摘出品

编译：Zhifu、JonyKai、汤圆、夏雅薇

可能经常你会听到一些很主观的评价比如“你太瘦了”或者“你怎么那么高”，但这里瘦或者高都是基于评价者的主观判断和视觉记忆做出的评述，并没有严格的参照。

作者从小被人说体型瘦小，于是他用了贝叶斯分析最终得出了自己的体重确实低于本国平均水平的结论。没有比直白的数字更有说服力了，想要说服自己健身的小伙伴不妨也试试作者的统计学暴击法！

这篇文章将叙述一个在线性回归理论中应用贝叶斯分析的有趣试验（一个小秘密：我在这篇文章中使用了公制计量单位呦）。

如文章标题所述，我将会对自己的体格进行一番科学的研究。

在开始之前，我希望大家多了解我一些：我在越南出生，在新加坡上高中，现在在美国上大学。我经常因为身形瘦小被人们取笑，说我应该多去健身房锻炼增肌，拥有更强壮的体魄。这些评价我一般都一笑置之。对于一个身高169厘米（5尺6寸），体重58公斤（127磅）的人来说，我的BMI指数(20.3)几近完美。

仔细一想，大家可能没说错：我比一般的越南男性要高，但却只有平均体重（维基百科里越南男性的平均体重是58公斤,平均身高是162厘米），“看起来”可能是要稍微瘦一些。

这里“看起来”是关键：背后的逻辑很清楚，不是吗？如果你把自己抻长一些，体重不变，那确实应该看起来苗条一些。我把这个看作是严肃的科学问题，并准备深入研究。

对于一个169厘米高的越南男性来说，我到底轻了多少呢？

我们需要一种有理有据的方式来研究这个问题。有个好方法是尽可能多地找到越南男子身高和体重的数据，来判断我在这个样本中的位置。

越南人口数据

在浏览各种网页后，我找到了一份调查研究数据，包含超过10,000名越南人的人口统计信息。将抽样条件设置为年龄18-29岁的越南男性，从而得到数量为383的样本，这个样本足以用来进行接下来的分析啦。

首先，通过人口体重的直方图，看看我在年轻越南男性中的体重分布位置。

红线表示样本的中位数，而橙色线表示平均值

这张直方图表明我的体重略低于383名年轻越南男子的体重平均值和中位数。看起来是与我们要研究的相关呀！然而问题并不在于我的体重与这个样本本身的比较，是假设这383个人可以代表越南男性，在身高169厘米的情况下，我的体重与整个越南人口相比处于一个什么位置。为此，我们需要进行回归分析。

首先绘制一个身高和体重的二维散点图

好吧，看起来我处在平均水平。但是如果我们只看身高169厘米的数据（想象一条垂直x轴于169厘米这个刻度并穿过红点的直线），我的体重在他们之中处于下游。

另一个重要的发现是越南男性身高和体重呈正相关。我们将进行定量分析来进一步了解这种关系。

首先，让我们快速添加“普通最小二乘”线。我稍后会详细介绍这一点，现在先在图上展示出来。

最小二乘线可以表示为y = -86.32 + 0.889x，这表明通常情况下，我这个年龄的越南男性，每增加1厘米的身高，体重会增加0.88千克。

但是，这并没有解决我们的问题；身高169厘米，体重58公斤到底是太沉，太轻还是刚刚好呢？要以定量的方式进一步解释这个问题，如果有所有身高1米68的人的体重分布，那么我的体重排在前25％，50％或75％的几率是多少？要弄清这一点，我们需要深入学习并理解回归背后的理论。

线性回归理论

在线性回归模型中，Y变量（在我们的例子中，是人的体重）是x（身高）的线性函数。在这个线性关系中截距和斜率分别为β0和β1；也就是说，假设E（Y | X = x）=β0+β1x。我们不知道β0和β1是多少，所以将它们视为未知参数。

在大多数标准线性回归模型中，我们进一步假设给定X = x的情况下，Y的条件分布是正态分布的。

这就是基本的线性回归模型：

可以被改写成：

注意，在许多模型中，我们可以用精度参数τ替换方差参数σ，其中τ= 1 /σ。

总结：因变量Y遵循由平均数μi和精度参数τ决定的正态分布。μi是由β0和β1决定的X的线性函数。

最后，我们还需假设未知方差不依赖于x；这种假设称为同方差性。

涉及的内容很多，都涵括在下面这张图里啦。

图像来源：Joseph Chang（耶鲁大学）

在实际的数据分析问题中，我们掌握的只是黑点 - 数据。使用这些数据，我们的目标是推断不知道的事情，包括β0，β1（在图片中的蓝色虚线）和σ（它决定了在给定一个y值的时候，红色正态分布密度的宽度）。注意，每个黑点周围的正态分布看起来完全相同。这是同方差性的本质。

估算参数

现在，有几种方法可以估算β0和β1。如果你使用普通最小二乘估计这样的模型，你不必担心概率公式，因为你正在寻找β0和β1的最佳值，而这是通过最小化拟合值与预测值的平方误差得到的。

另一方面，你可以使用最大似然估计（MLE）来估计此类模型：通过最大化似然函数来寻找参数的最佳值。

注意：一个有趣的结果（我没有放上具体的数学证明）是，如果我们假设误差项也属于正态分布的话，那么最小二乘估计的结果也是最大似然估计的结果。

贝叶斯方法的线性回归

不同于最大化似然函数，贝叶斯方法会假设参数服从一个先验分布。根据贝叶斯公式计算出参数后验概率：

贝叶斯方法的似然函数同上面的类似，不同之处在于加入了对估计参数β0，β1，τ的先验分布：

等等，什么是先验分布，为什么这会让公式看上去更加复杂？

请相信我，先验分布虽然看上去很奇怪，但实际上很直观。事实上，我们对未知的参数（例子中的β0，β1，τ）分配一个看上去很随意的先验分布，这里存在着很强的哲学缘由。

先验分布能够反映出在没看见数据之前我们对数据的假设理解。在观察过一些数据之后，我们应用贝叶斯公式，就得到了同时考虑到了先验和数据的未知参数后验分布。根据后验分布，我们就能预测出未来的数据的分布。

最终的参数估计虽然取决于数据和先验分布，但是如果数据中包含的信息越多，那先验的影响就越小。

那么我该如何选择先验分布

这是个好问题，因为这里存在着无数种可能。理论上只有存在一个正确的先验能够能够反应出你的先验假设。但是在实际中，先验分布的选取是相当的主观，甚至有时可以是任意的。

我们可以选择一个有着较大标准方差（意味着精度很低）的正态分布先验。比如，我们假设参数β0和β1是服从均值为0标准方差为10000的正态分布。这种分布是毫无信息的分布，因为分布十分平坦（这意味着，参数在任意区间的取值概率几乎相同）。

如果选取了这种类型的先验分布，那么我们就不用考虑在这类分布中哪种分布更好，因为分布几乎都很平坦，在每个地方的概率都可以忽略不计。此外，后验分布不会受这种分布的影响。

同样，对于精度τ，因为其必须是非负的，所以需要选取一个取值限定在非负范围的分布。比如，在这里可以选取一个带有较小形状和尺度参数的伽马分布。

另外一种很有用的不附带信息的分布是均匀分布。如果对 σ 或 τ选择了均匀分布，那么你最终将会得到这样的模型。如下图所示，由John K. Kruschke绘制。

John K. Kruschke绘制的模型

用R语言和JAGS模拟数据

到目前为止，我们仍只停留在理论阶段。大多数情况下。后验分布并不能直接得到（想想正态分布和伽马分布有多复杂，然后还要再将他们乘起来）。Markov Chain Monte Carlo方法常常用来估计模型的参数。利用JAGS就能帮我们完成。

“等一下！！！这个工具真能够帮我们解决这些复杂的公式么？”

JAGS模型是一个基于Markov Chain Monte Carlo（MCMC）的仿真过程，它能够生成出参数空间θ=(β0;β1;τ)的许多次迭代。由参数空间中每个参数生成的样本分布会估计出参数最有可能的分布

为什们是这样呢？这个解释起来十分复杂，已经超出了本篇介绍的范围。直接来说就是：MCMC通过构建一个马尔可夫链产生了服从后验分布的样本，这个马尔可夫链有着同样的目标后验分布！？

这个过程很没意思。最快的方法就是：不去解等式（2），因为通常不可能得到解析解。我们能做的就是聪明的采样，而在数学上证明了这些样本确实服从以β0，β1，τ为参数的分布。

那么挥别了数学之后，我该怎么使用JAGS？

我们按下面的步骤在R语言中运行JAGS

首先以文本的形式写下模型

然后，我们让JAGS执行仿真模拟。这里我使用JAGs对参数空间θ进行10000次模拟

抽样之后，我们就得到了θ=(β0;β1;τ)的抽样数据，如下表所示：

“看上去好酷，那又怎样呢？”

现在我们对参数空间θ进行10000次迭代，根据等式

这就意味着，如果用x=169cm替代每个迭代值，我们将会得到10000个体重值，也就是身高169cm情况下体重的分布。

我们都对以我身高为条件下体重的百分比分布很感兴趣。为了达到这个目的，需要找到基于我身高的体重分布。

上面这张图表明我的体重（给定169的身高）最有可能在模拟越南人口中的后30%左右。

比如，我们能发现我的体重在前40%甚至更少的位置

因此大量证据表明，在身高169的条件下，体重58kg会让我处于越南人口的较低百分比处。我确实需要去健身房锻炼并增加些体重了。毕竟如果你不信详尽的贝叶斯统计分析，还能相信什么呢？

我有一份包含了美国8169名年轻男性和女性身高体重的数据（国家寿命调查1997），你能做同样的分析么，看看你会得到什么样的结论？