最大似然估计（ML）和REML

A.F. Zuur et al., Mixed Effects Models and Extensions in Ecology with R, p116-119, chr 5.6

当利用混合效应建模时，你会遇到诸如REML和ML这样的词汇。不像线性回归模型，就算你不知道背后的数学理论也可以照样使用。但是在混合效应建模中，你必须懂得一些相关的数学知识。所以REML是什么意思？它能干什么？
第一个问题比较简单，REML限制性最大似然估计英文首字母的缩写。但是对于第二个问题，大多数教材在这点上就变得相当专业，或者解释得比较粗略，只是提到它是矫正自由度的一个神秘方法【ML没有考虑到估计固定效应带来自由度的损失，造成参数的低估】^[1]。而我们则选择尝试更为详细地解释它，所以需要利用矩阵代数的知识。但是为了理解REML，首先需要理解最大似然估计的原理，我们从这儿开始。如果你不熟悉矩阵代数，或者如果这一部分对数学水平的要求太高，我们仍建议你跳过这一部分。
我们首先回顾用于线性回归的最大似然，然后给出REML是如何用于矫正方差估计量的。
假设有一个线性回归模型
$Y_{i}=\alpha+\beta \times X_{i}+\varepsilon_{i}$ ，其中 $\varepsilon_{i}\sim N(0,\sigma^2)$ 。模型中有3个未知参数 $\alpha$ 、 $\beta$ 和 $\sigma$ 。为了简便，令 $\boldsymbol{\theta}=(\alpha, \beta, \sigma)$ 。普通最小二乘是估计 $\boldsymbol{\theta}$ 的一个方法，它给出 $\boldsymbol{\theta}$ 中每一个元素的表达式。利用线性回归获取方差估计的表达式是：
$\hat{\sigma}^{2}=\frac{1}{n-2} \sum_{i=1}^{n}\left(Y_{i}-\hat{\alpha}-\hat{\beta} \times X_{i}\right)^{2} \quad\quad(5.14)$
我们给参数加了一个帽子^，表示它是估计值， $n$ 是观测值的个数。可以证明 $\hat{\sigma}$ 是 $\sigma$ 的无偏估计，意味着 $E[\hat{\sigma}]=\sigma$ 。现在让我们看一看最大似然估计方法。假设 $Y_i$ 服从正态分布，其密度函数为：
$f_{i}\left(Y_{i}, X_{i}, \alpha, \beta, \sigma\right)=\frac{1}{\sigma \sqrt{2 \pi}} \mathrm{e}^{\frac{\left(V_{i}-\alpha-\beta \times X_{i}\right)^{2}}{2 \sigma^{2}}}\quad\quad(5.15)$

因为我们也假设了 $Y_i$ 是独立的，可以将 $Y_{1}, Y_{2}, \dots, Y_{n}$ 的联合密度函数写成单个密度曲线 $f_{1}, f_{2}, \dots, f_{n}$ 乘积的形式。这个乘积就叫做似然函数L。它是数据和 $\boldsymbol{\theta}$ 的一个函数。问题是如何选择 $\boldsymbol{\theta}$ 使L最大。为了简化，对L取自然对数，得到如下的log似然方程式：
$\ln L\left(Y_{i}, X_{i}, \alpha, \beta, \sigma\right)=-\frac{n}{2} \ln 2 \pi-\frac{n}{2} \ln \sigma^{2}-\frac{1}{2 \sigma^{2}} \sum_{i=1}^{n}\left(Y_{i}-\alpha-\beta \times X_{i}\right)^{2}\quad\quad(5.16)$
我们需要最大化这个式子，问题就变成了对每一个参数偏导，令偏导数=0并求解。因为我们很容易计算，这些偏导数=0的式子称之为封闭解。对于广义线性混合模型我们将会看到开放解，意思是参数没有直接的解。
这里没有给出 $\alpha$ 和 $\beta$ 估计量的式子，但对于方差我们得到：
$\hat{\sigma}^{2}=\frac{1}{n} \sum_{i=1}^{n}\left(Y_{i}-\hat{\alpha}-\hat{\beta} \times X_{i}\right)^{2}\quad\quad(5.17)$
注意这个式子与我们利用普通最小二乘得到的式子（5.14）非常相似。实际上，受到因子 $(n-2) / n$ 的影响，利用最大似然得到的方差估计量是有偏的（回归分析为什么误差方差中自由度是n-2？）。如果线性回归模型含有p个解释变量，那么偏度是 $(n-p) / n$ 。最大似然是有偏的的原因是它忽略了截距和斜率也被估计的事实。所以我们需要更好的ML估计量，而这正是REML所做的事情。

REML的工作如下：有线性回归模型 $Y_{i}=\alpha+\beta \times X_{i}+\varepsilon_{i}$ 可以写成 $Y_{i}=\mathbf{X}_{i} \times \boldsymbol{\beta}+\varepsilon_{i}$ 。这是简单的矩阵形式， $\mathbf{X}_{i}=\left(1 X_{i}\right)$ ， $\boldsymbol{\beta}$ 的第一个元素是截距，第二个元素是原始的 $\beta$ 。正态性假设意味着
$Y_{i} \sim N\left(\mathbf{X}_{i} \times \beta, \sigma^{2}\right)\quad\quad(5.18)$
用ML估计量的问题是我们不得不估计式子5.18中 $\boldsymbol{\beta}$ 中的截距和斜率。显然，如果没有 $\boldsymbol{\beta}$ ，就能解决问题。为了消除 $\boldsymbol{\beta}$ ，可以找到一个维度 $n \times (n – 2)$ 的特殊矩阵 $\mathbf{A}$ ，特殊指的是“与 $\mathbf{A}^\prime$ 正交”，然后用这个矩阵乘以 $\mathbf{Y}=\left(Y_{1}, \dots,Y_{n} \right)^{\prime}$ 之后再用ML估计。正交指的是如果 $\mathbf{A}$ 与 $\mathbf{X}$ 相乘，结果是0。因此，我们得到 $\mathbf{A}^{\prime} \times \mathbf{Y}=\mathbf{A}^{\prime} \times \mathbf{X} \times \boldsymbol{\beta}+\mathbf{A}^{\prime} \times \boldsymbol{\varepsilon}=\mathbf{0}+\mathbf{A}^{\prime} \times \boldsymbol{\varepsilon}=\mathbf{A}^{\prime} \times \boldsymbol{\varepsilon}$ 。现在 $\mathbf{A}^{\prime} \times \mathbf{Y}$ 的分布是
$\mathbf{A}^{\prime} \times \mathbf{Y} \sim N\left(\mathbf{0}, \sigma^{2} \times \mathbf{A}^{\prime} \times \mathbf{A}\right)$
而不再依赖 $\boldsymbol{\beta}$ 。那么对 $\mathbf{A}^{\prime} \times \mathbf{Y}$ 进行似然估计就会得到 $\sigma^{2}$ 的无偏估计量（5.14）。现在我们讨论REML如何应用到混合线性模型。我们的起点是边际模型
$\mathbf{Y}_{i} \sim N\left(\mathbf{X}_{i} \times \boldsymbol{\beta}, \mathbf{V}_{i}\right) \quad \mathbf{V}_{i}=\mathbf{Z}_{i} \times \mathbf{D} \times \mathbf{Z}_{i}^{\prime}+\mathbf{\Sigma}_{i}$
故事又重新开始，如之前，我们可以写一个略微不同的log似然式子。未知参数是 $\boldsymbol{\beta}$ 和 $\mathbf{D}$ 及 $\mathbf{\Sigma}_{i}$ 中的元素，依然用 $\boldsymbol{\theta}$ 表示。似然函数：
$\ln L\left(\mathbf{Y}_{i}, \mathbf{X}_{i}, \boldsymbol{\theta}\right)=-\frac{n}{2} \ln 2 \pi-\frac{1}{2} \sum_{i=1}^{n} \ln \left|\mathbf{V}_{i}\right|-\frac{1}{2} \sum_{i=1}^{n}\left(\mathbf{Y}_{i}-\mathbf{X}_{i} \times \boldsymbol{\beta}\right)^{\prime} \times \mathbf{V}_{i}^{-1} \times\left(\mathbf{Y}_{i}-\mathbf{X}_{i} \times \boldsymbol{\beta}\right)$
$|\mathbf{V}_{i}|$ 是 $\mathbf{V}_{i}$ 的行列式。对 $\boldsymbol{\beta}$ 求偏导并=0解方程。如之前讨论的例子，得到的参数是有偏的，因此我们需要REML。
总之，REML，就是用一个特殊的矩阵乘以Y，这样X×β就消去。然后用ML估计得到的参数估计子就是无偏的，并且与特定的矩阵相乘无关。因此， $\boldsymbol{\beta}$ 的REML估计子与ML的估计子不同。如果相对于观测值的个数，固定协变量的个数很少，就没有太大的不同，相反有许多的固定协变量，情况就大不相同。

Lindstrom MJ, Bates DM (1988) Newton—Raphson and EM Algorithms for Linear Mixed-Effects Models for Repeated-Measures Data. J Am Stat Assoc 83:1014–1022. ↩

最后编辑于：2019.04.10 12:06:32

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 203,456评论 5赞 477
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 85,370评论 2赞 381
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 150,337评论 0赞 337
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 54,583评论 1赞 273
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 63,596评论 5赞 365
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 48,572评论 1赞 281
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 37,936评论 3赞 395
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,595评论 0赞 258
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 40,850评论 1赞 297
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,601评论 2赞 321
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,685评论 1赞 329
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,371评论 4赞 318
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 38,951评论 3赞 307
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 29,934评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 31,167评论 1赞 259
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 43,636评论 2赞 349
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 42,411评论 2赞 342

最大似然估计（ML）和REML

推荐阅读更多精彩内容