昔日因，今日意 | 统计之都

飞帅云：“三十功名尘与土，八千里路云和月。莫等闲，白了少年头，空悲切。” 可我在耶鲁两年多了，基本一事无成。既没有像当年那样死磕 Lasso 和 Boosting，也没有能追随 Deep Learning 的浪潮。曾经真的以为人生就这样了，平静的心拒绝再有浪潮。斩了千次的情丝却断不了，百转千折她将我围绕。有人问我她究竟是哪里好？我想我是鬼迷心窍。

她就是 LMM，我给她起了一个美丽的中文名:“林妹妹”。

对我这种工科男，与林妹妹相知相识，是需要一段奇缘。从在浙大本科自动化专业入学，到港科大的电子系博士快毕业，曾经有且仅有一次机会与她相识，还是被很傻很天真的我错过了。现在不管我怎么念 “菠萝菠萝蜜”，时光还是不会倒流的。我只是想，如果上天可以给我一个机会再来一次的话，我会对她说八个字：“我们好像在哪见过？” 然而，有缘人终归是有缘人，奇妙的感觉就在点火的那一刹那。

2010 年，夏。

香港，清水湾。

海浪拍打着沙滩，伴随着风的欢笑，涌出蔚蓝的一片。

那一天，香港科大来了一位远方的客人，便是我现在耶鲁的老板。在他的演讲中，提到了一个故事。本人平生听了很多故事，但是我认为这个是最精彩的。据一些专家估算，人类身高的差异（统计学上用方差来描述），70% 左右由遗传决定。另一批专家利用最新科技做了全基因组扫描，发现了大约 100 多个遗传变异点与人类身高显著的相关。这批专家就用这些显著的变异点去解释身高的方差时，惊讶地发现它们只能解释 5%！从 5% 到 70%，我靠，这中间的差距也太大了吧！“到底哪批专家是砖家，或者都是？” 这个问题在我脑子里油然而生。然而以我的智慧，只猜中了开头，却猜不到结尾。林妹妹的出现让剧情峰回路转……

如今，专家们基本已达成共识：不要只是去关注那些显著的变异点！虽然这些不显著的变异点每一个的作用都比较小，但是他们总的作用却不能忽略！如果把那些不显著变异点的作用一起考虑进去，就能解释身高方差的 45%，如果再考虑上那些没有被直接观察到的变异点的影响，就基本上接近 70%。如何能把那么多不显著变异点的作用都优雅地考虑进去呢？这里就需要林妹妹了。

前面讲的故事是当今生命科学中最重要的课题之一。2009 年的时候，科学家们还专门给这个故事起了一个名字–“missing heritability”，用今年流行的语言翻译过来就是 “遗传物质都去哪儿了?” 身高只是其中一个例子，对很多复杂疾病，比如糖尿病，高血压，精神分裂症等，科学家们也发现类似的情况。这个 missing heritability 类似于物理学上的暗物质，感觉它存在却看不到它。林妹妹的出现让我们真实地测量到遗传学中的“暗物质”，并确认它的存在。

好了，是时候停止卖萌，进入主题了。LMM 全称是 Linear Mixed Model（混合线性模型）。她血统高贵，与现代统计学之父 Ronald Fisher 提出的随机效应一脉相承。上个世纪 50 年代，Charles Henderson 为她打造了国际一流的统计性质（BLUE and BLUP）1，他的学生 Shayle Searle 更是为她配上了 “黑客帝国（Matrix）” 的装备，从此她的名字将永远记入统计学的史册。1991 年，statistical science 上有一篇很经典的文章 “That BLUP is a Good Thing: The Estimation of Random Effects”2，里面谈到了她许多超一流的品质。事实上，我们在实践中已经用到了她很多好的性质，只不过我们以前不知道罢了。

现在从她的一副黑客帝国装备说起，因为这副装备低调奢华有内涵：

y=Xβ+Zu+e,u∼N(0,σ2uI),e∼N(0,σ2eI)

这里y∈Rn是回归问题中的因变量，X∈Rn×d和Z∈Rn×p分别是固定效应和随机效应的设计矩阵，e∈Rn是随机误差，β∈Rd和u∈Rp分别是固定效应向量与随机效应向量，n,d,p分别是样本数目，固定效应的个数以及随机效应的个数。这里y,X,Z是给定的，需要估计的是β,u,σ2u,σ2e。 看到这里，我知道大部分童鞋已经有点晕了：啥叫固定效应，啥叫随机效应？先解释什么叫随机效应。 我相信大家都理解随机数，简单点说，他们就是从某一个分布里面随机抽出来的数，这些数不是固定的，但是他们总体上服从某种规律（即某种分布），比如服从正态分布。之所以用 “随机效应” 而不是用“随机数”，是为了描述设计矩阵的每一列所对应的变量对因变量y的作用，比如在上述模型中的u是一个p维的向量，它的每个元素即uj,j=1,…,p都来自于正态分布N(0,σ2u)，uj即是Z的第j列对y的效应。现在来解释啥叫固定效应，一句话，固定效应就是非随机效应。当固定效应与随机效应在一起的时候，就是所谓的 mixed model。 注意千万不要把 mixed model 与 mixture model 混为一谈！因为前者是被动在一起的，后者则是主动在一起的，想分都分不开。

既然是被动在一起，把二者拆开就比较容易。如果只看固定效应那一部分，

y=Xβ+e,e∼N(0,σ2eI),

这就是最基本的多元线性回归，β的最优解由最小二乘法 (Least square) 给出，即ˆβLS，它的统计性质由高斯马尔科夫定理 (Gauss–Markov Theorem) 保证，比如ˆβLS是所有无偏估计中方差最小的。如果只看随机效应那一部分，

y=Zu+e,u∼N(0,σ2uI),e∼N(0,σ2eI),

相信大家对上式也不会陌生。比如，看过 Pattern recognition and machine learning (by Bishop, C.) 的第三章 Linear Models for Regression 的童鞋应该会发现：当λ=σ2eσ2u时，求解上式与下式得到的解则是完全等价。

minu‖y−Zu‖2+λ‖u‖2,

不同的是，这里的λ通常由交叉验证确定，而估计前面式子中的参数 (σ2u与σ2e) 则另有办法。在机器学习中，大家把它叫做 Evidence approximation，统计学里面把它叫经验贝叶斯 (Empirical Bayes)。

现在可以再把二者合在一起了，但这里涉及到一个重要的问题。举一个简单的例子，有n个数据点(x1,x2,…,xn)，每个数据点xi∈R都独立地来自N(μ,σ2)。对σ2的最大似然估计是˜σ2=∑i(xi−ˉx)2/n，这里ˉx=∑ixi/n是均值。但是，正如大家所知道的，对σ2的无偏估计应该是ˆσ2=∑i(xi−ˉx)2/(n−1)，因为在估计均值的时候已经消耗掉一个数据了。为了补偿在有若干个固定效应情况下对方差估计的偏差，有一个办法应运而生，它的名字叫 REML (REstricted Maximum Likelihood，约束极大似然或限制极大似然)。

前面提到的探索遗传学中的暗物质是对 LMM 一次高端大气上档次的运用。遗传的变异会引起身高的差异，那么身高有多大程度上是由遗传因素决定的？翻译成统计学语言：用遗传变异点数据究竟能解释身高方差的百分之多少？回过头来再看，y中是n个样本 (∼5000) 的身高数据，X的每一列对应一个协变量，比如年龄、性别，基因组变异点的数据都放到Z中，其中每一列对应一个变异点的数据。为了写出y的边际分布，需要对u和e积分，

y∼N(Xβ,ZZTσ2u+σ2eI).

注意Xβ并不影响y的方差。Heritability 定义为：

h2=pσ2upσ2u+σ2e,

这里pσ2u是遗传因素解释的方差（p是变异点个数，大约是 50 万到 100 万这个范围），σ2e是非遗传因素造成的方差。启动 REML 以后，就能得到 (ˆσ2u,ˆσ2e)，然后算出 heritability。

3 水榭听香，指点群豪戏

“脚步不能达到的地方，眼光可以到达。” 抬望眼，满城尽是 LMM，如下图所示。由于篇幅所限，我只能简单地介绍一部分。

LMM 与很多经典方法的关系图

第一，LMM 与 JSE-Ridge Regression 的关系最为明显。当没有固定效应，Z变为单位矩阵的时候，LMM 就变为了 JSE（这个时候需要σ2e 是已知的，不然会有可辨识性的问题。在 JSE(James-Stein Estimator) 的问题中，σ2e=1，更详细的描述请参考我的《那些年我们一起追的 EB》）。LMM 与 Ridge 的关系，前面已经讲过了。

第二，RVM(relevance vector machine) 如下：

y=∑jZjuj+e,uj∼N(0,σ2uj),e∼N(0,σ2eI),

RVM 与 LMM 的差别是：RVM 允许每个随机效应uj有自己的方差σ2uj，而 LMM 中所有的随机效应 (random effects) 具有同样的σ2u。 从 RVM 到 Lasso(Least Absolute Shrinkage and Selection Operator)，只需要假设{σ2uj}来自指数分布，详情参见 Bayesian Lasso。Lasso，LARS 和 Boosting 已经成为统计学与计算科学史上的一段佳话，最好的文献当然是 Efron 教授的 Least Angle Regression，喜欢看故事的童鞋可以看我的《统计学习那些事》。

第三，LMM 与 PCA(Principal Component Analysis) 的联系似乎不是那么直接，因为这里已经从监督学习走向了非监督学习。然而，当 PCA 被赋予概率的解释后，天堑变通途。这篇里程碑式的文章就是 Probabilistic principal component analysis by Tipping and Bishop。PCA 与 clustering 的亲密关系暴露在本世纪初，clustering 和 mixture models 的关系嘛，应该是不言而喻的。

第四，大家都知道，Mixture models 里面是有一个隐状态，比如在做 clustering 的时候，这个隐状态就用来表明数据点与 cluster 的隶属关系。当这些隐状态不再是独立等分布的时候，比如，后一个状态取决于前一个状态的时候，HMM 便应运而生。HMM 与卡曼滤波 (Kalman filter)基本上可以看做孪生兄弟，一个为离散状态而生，一个为连续状态而来。LMM 与卡曼滤波的关系在 “BLUP is a good thing” 这篇雄文中早有讨论。当年学控制的我与卡曼滤波有过初步接触，但是却与 LMM 失之交臂，还好在耶鲁与 LMM 再续前缘。

第五，如今的 Matrix factorization 已经是令人眼花缭乱了，因为这里加入了很多 sparse(包括 low-rank) 与 smoothing 的技术。但不可否认，PCA 依然是矩阵分解中最重要的一种，奇异值分解依然是这里最重要的数学基石。

面对如此波澜壮阔的模型表演，不知道大家会如何感想？这里我先引用 Terry Speed 在 “BLUP is a good thing” 的评论里的最后一段话

In closing these few remarks, I cannot resist paraphrasing I.J. Good’s memorable aphorism: ‘To a Bayesian, all things are Bayesian.’ How does ‘To a non-Bayesian, all things are BLUPs’ sound as a summary of this fine paper?

大师的话值得久久回味…… 我自己总结的话，来点通俗易懂的，还是这句 “天下武功，若说邪的，那是各有各的邪法，若说正的，则都有一种‘天下武功出少林’的感觉”。

4 杏子林中，商略平生義

“眼光不能到达的地方，精神可以飞到。”

“随机还是非随机？”是一个问题，甚至是一个哲学问题。或许，我们参一生也参不透这道难题。爱因斯坦说:“上帝不玩骰子。”然而，麦克斯韦却说：“这个世界真正的逻辑就是概率的计算。” 电影《美丽心灵》的纳什也在追问 “到底什么才是真正的逻辑”。最后他在获得诺奖时说：“ 我一直以来都坚信数字，不管是方程还是逻辑都引导我们去思考。但是在如此追求了一生后，我问自己：‘逻辑到底是什么？谁决定原由？’我的探索让我从形而下到形而上，最后到了妄想症，就这样来回走了一趟。在事业上我有了重大突破，在生命中我也找到了最重要的人：只有在这种神秘的爱情方程中，才能找到逻辑或原由来。” 这是我听到的最美的答案。

如果回到工程实践的话，或许我们应该追问：“为什么引入随机效应后会有如此神奇的疗效？”Efron 教授在他的一篇文章中称赞 James-Stein Estimator：“This is the single most striking result of post-World War II statistical theory”。 我想，我们应该可以从 JSE 中寻找到一些蛛丝马迹。JSE 的原问题是：现已观察到N个z值，即[z1,z2,…,zN]，还知道zi独立地来自以μi为均值，1 为方差的正态分布，即zi|μi∼N(μi,1), i=1,2,…,N. 问题是：如何从观察到的z=[z1,z2,…,zN]估计μ=[μ1,μ2,…,μN]？最大似然法和 James-Stein Estimator 给出解答分别是（更多细节参考《那些年，我们一起追的 EB》）：

ˆμML=z,ˆμJS=(1−N−2||z||2)z.

对 JSE 的理解有很多不同的角度，个人觉得从下面的这个角度看过去是非常精彩的。如果我们把μi,i=1,…,N看做是随机，那么我们可以认为他们来自某一个分布G(μ)，随着N的增大，我们对G(μ)的估计就会越准确。原来看似独立的zi却能通过这样一个分层结构（见下图）来共享信息。虽然N很小的时候，G(μ)是没法估计准确的，但幸运的是，这里的N并不要求太大。可以证明，只要N≥3，JSE 就比 MLE 好。这就是 Efron 教授所说的 “learning from others”。如果用更加数学的语言来刻画信息共享，其实就是 Bias-Variance trade-off。当信息共享的时候，偏差增加了少许但方差却大大降低。

5 须倾英雄泪

James Watson 在他的《双螺旋》一书的序言中写道：“科学的发现很少会像门外汉所想像的那样，按照直截了当、合乎逻辑的方式进行。事实上，科学的进步（有时是倒退）往往是人为事件。在这些事件中，人性以及文化传统都起着巨大的作用。”

“庾信平生最萧瑟，暮年诗赋动江关。”这是张益唐教授为 “孪生素数猜想” 作出巨大贡献后接受采访时引用的诗句。听罢，令人感慨万千。

“胡汉恩仇”，须倾英雄泪。

最后附上本文的 pdf 原稿供下载：LMM_Yale

杨灿于 2011 年在香港科技大学电子计算机工程系获得博士学位，2011-2012 为耶鲁大学生物统计系博士后，现为耶鲁大学副研究员。

敬告各位友媒，如需转载，请与统计之都小编联系（直接留言或发至邮箱：[email protected]），获准转载的请在显著位置注明作者和出处（转载自：统计之都），并在文章结尾处附上统计之都微信二维码。

← 在 R 中使用管道操作 失联搜救中的统计数据分析 →

发表 / 查看评论