最大似然估計(jì)（ML）和REML

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

當(dāng)利用混合效應(yīng)建模時闰蚕，你會遇到諸如REML和ML這樣的詞匯。不像線性回歸模型雷酪，就算你不知道背后的數(shù)學(xué)理論也可以照樣使用阿宅。但是在混合效應(yīng)建模中嚼酝，你必須懂得一些相關(guān)的數(shù)學(xué)知識莺琳。所以REML是什么意思嚷缭？它能干什么蝶桶？
第一個問題比較簡單，REML限制性最大似然估計(jì)英文首字母的縮寫爱榔。但是對于第二個問題被环，大多數(shù)教材在這點(diǎn)上就變得相當(dāng)專業(yè)，或者解釋得比較粗略详幽，只是提到它是矯正自由度的一個神秘方法【ML沒有考慮到估計(jì)固定效應(yīng)帶來自由度的損失筛欢，造成參數(shù)的低估】^[1]。而我們則選擇嘗試更為詳細(xì)地解釋它唇聘，所以需要利用矩陣代數(shù)的知識版姑。但是為了理解REML，首先需要理解最大似然估計(jì)的原理迟郎，我們從這兒開始剥险。如果你不熟悉矩陣代數(shù)，或者如果這一部分對數(shù)學(xué)水平的要求太高宪肖，我們?nèi)越ㄗh你跳過這一部分表制。
我們首先回顧用于線性回歸的最大似然，然后給出REML是如何用于矯正方差估計(jì)量的控乾。
假設(shè)有一個線性回歸模型
$Y_{i}=\alpha+\beta \times X_{i}+\varepsilon_{i}$ 么介，其中 $\varepsilon_{i}\sim N(0,\sigma^2)$ 。模型中有3個未知參數(shù) $\alpha$ 蜕衡、 $\beta$ 和 $\sigma$ 壤短。為了簡便，令 $\boldsymbol{\theta}=(\alpha, \beta, \sigma)$ 慨仿。普通最小二乘是估計(jì) $\boldsymbol{\theta}$ 的一個方法久脯，它給出 $\boldsymbol{\theta}$ 中每一個元素的表達(dá)式。利用線性回歸獲取方差估計(jì)的表達(dá)式是：
$\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)$
我們給參數(shù)加了一個帽子^镰吆，表示它是估計(jì)值桶现， $n$ 是觀測值的個數(shù)《︽ⅲ可以證明 $\hat{\sigma}$ 是 $\sigma$ 的無偏估計(jì)，意味著 $E[\hat{\sigma}]=\sigma$ ∠嗔蓿現(xiàn)在讓我們看一看最大似然估計(jì)方法相寇。假設(shè) $Y_i$ 服從正態(tài)分布，其密度函數(shù)為：
$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àn)槲覀円布僭O(shè)了 $Y_i$ 是獨(dú)立的钮科，可以將 $Y_{1}, Y_{2}, \dots, Y_{n}$ 的聯(lián)合密度函數(shù)寫成單個密度曲線 $f_{1}, f_{2}, \dots, f_{n}$ 乘積的形式唤衫。這個乘積就叫做似然函數(shù)L。它是數(shù)據(jù)和 $\boldsymbol{\theta}$ 的一個函數(shù)绵脯。問題是如何選擇 $\boldsymbol{\theta}$ 使L最大佳励。為了簡化休里，對L取自然對數(shù)，得到如下的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)$
我們需要最大化這個式子赃承，問題就變成了對每一個參數(shù)偏導(dǎo)妙黍，令偏導(dǎo)數(shù)=0并求解。因?yàn)槲覀兒苋菀子?jì)算瞧剖，這些偏導(dǎo)數(shù)=0的式子稱之為封閉解拭嫁。對于廣義線性混合模型我們將會看到開放解，意思是參數(shù)沒有直接的解抓于。
這里沒有給出 $\alpha$ 和 $\beta$ 估計(jì)量的式子做粤，但對于方差我們得到：
$\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）非常相似。實(shí)際上捉撮，受到因子 $(n-2) / n$ 的影響怕品，利用最大似然得到的方差估計(jì)量是有偏的（回歸分析為什么誤差方差中自由度是n-2？）巾遭。如果線性回歸模型含有p個解釋變量肉康，那么偏度是 $(n-p) / n$ 。最大似然是有偏的的原因是它忽略了截距和斜率也被估計(jì)的事實(shí)恢总。所以我們需要更好的ML估計(jì)量迎罗，而這正是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$ 厢岂。正態(tài)性假設(shè)意味著
$Y_{i} \sim N\left(\mathbf{X}_{i} \times \beta, \sigma^{2}\right)\quad\quad(5.18)$
用ML估計(jì)量的問題是我們不得不估計(jì)式子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估計(jì)咖熟。正交指的是如果 $\mathbf{A}$ 與 $\mathbf{X}$ 相乘圃酵，結(jié)果是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}$ 」停現(xiàn)在 $\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}$ 進(jìn)行似然估計(jì)就會得到 $\sigma^{2}$ 的無偏估計(jì)量（5.14）∪贩校現(xiàn)在我們討論REML如何應(yīng)用到混合線性模型捌锭。我們的起點(diǎn)是邊際模型
$\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似然式子观谦。未知參數(shù)是 $\boldsymbol{\beta}$ 和 $\mathbf{D}$ 及 $\mathbf{\Sigma}_{i}$ 中的元素拉盾，依然用 $\boldsymbol{\theta}$ 表示。似然函數(shù)：
$\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}$ 求偏導(dǎo)并=0解方程盾剩。如之前討論的例子，得到的參數(shù)是有偏的替蔬，因此我們需要REML告私。
總之，REML承桥，就是用一個特殊的矩陣乘以Y驻粟，這樣X×β就消去。然后用ML估計(jì)得到的參數(shù)估計(jì)子就是無偏的凶异，并且與特定的矩陣相乘無關(guān)蜀撑。因此， $\boldsymbol{\beta}$ 的REML估計(jì)子與ML的估計(jì)子不同剩彬。如果相對于觀測值的個數(shù)酷麦，固定協(xié)變量的個數(shù)很少，就沒有太大的不同喉恋，相反有許多的固定協(xié)變量沃饶，情況就大不相同。

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

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者

人面猴
序言：七十年代末轻黑，一起剝皮案震驚了整個濱河市糊肤，隨后出現(xiàn)的幾起案子，更是在濱河造成了極大的恐慌氓鄙，老刑警劉巖馆揉，帶你破解...
沈念sama閱讀 221,820評論 6贊 515
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件，死亡現(xiàn)場離奇詭異抖拦，居然都是意外死亡升酣，警方通過查閱死者的電腦和手機(jī)，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 94,648評論 3贊 399
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進(jìn)店門态罪，熙熙樓的掌柜王于貴愁眉苦臉地迎上來噩茄，“玉大人，你說我怎么就攤上這事向臀。” “怎么了诸狭？”我有些...
開封第一講書人閱讀 168,324評論 0贊 360
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵券膀，是天一觀的道長君纫。經(jīng)常有香客問我，道長芹彬，這世上最難降的妖魔是什么蓄髓？我笑而不...
開封第一講書人閱讀 59,714評論 1贊 297
?港島之戀（遺憾婚禮）
正文為了忘掉前任，我火速辦了婚禮舒帮，結(jié)果婚禮上会喝，老公的妹妹穿的比我還像新娘。我一直安慰自己玩郊，他們只是感情好肢执，可當(dāng)我...
茶點(diǎn)故事閱讀 68,724評論 6贊 397
惡毒庶女頂嫁案：這布局不是一般人想出來的
文/花漫我一把揭開白布。她就那樣靜靜地躺著译红，像睡著了一般预茄。火紅的嫁衣襯著肌膚如雪。梳的紋絲不亂的頭發(fā)上侦厚，一...
開封第一講書人閱讀 52,328評論 1贊 310
城市分裂傳說
那天耻陕，我揣著相機(jī)與錄音，去河邊找鬼刨沦。笑死诗宣，一個胖子當(dāng)著我的面吹牛，可吹牛的內(nèi)容都是我干的想诅。我是一名探鬼主播召庞，決...
沈念sama閱讀 40,897評論 3贊 421
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開眼，長吁一口氣：“原來是場噩夢啊……” “哼侧蘸！你這毒婦竟也來了裁眯？” 一聲冷哼從身側(cè)響起，我...
開封第一講書人閱讀 39,804評論 0贊 276
萬榮殺人案實(shí)錄
序言：老撾萬榮一對情侶失蹤讳癌，失蹤者是張志新（化名）和其女友劉穎穿稳，沒想到半個月后，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體晌坤，經(jīng)...
沈念sama閱讀 46,345評論 1贊 318
?護(hù)林員之死
正文獨(dú)居荒郊野嶺守林人離奇死亡逢艘，尸身上長有42處帶血的膿包…… 初始之章·張勛以下內(nèi)容為張勛視角年9月15日...
茶點(diǎn)故事閱讀 38,431評論 3贊 340
?白月光啟示錄
正文我和宋清朗相戀三年，在試婚紗的時候發(fā)現(xiàn)自己被綠了骤菠。大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片它改。...
茶點(diǎn)故事閱讀 40,561評論 1贊 352
活死人
序言：一個原本活蹦亂跳的男人離奇死亡，死狀恐怖商乎，靈堂內(nèi)的尸體忽然破棺而出央拖，到底是詐尸還是另有隱情，我是刑警寧澤，帶...
沈念sama閱讀 36,238評論 5贊 350
?日本核電站爆炸內(nèi)幕
正文年R本政府宣布鲜戒，位于F島的核電站专控，受9級特大地震影響，放射性物質(zhì)發(fā)生泄漏遏餐。R本人自食惡果不足惜伦腐，卻給世界環(huán)境...
茶點(diǎn)故事閱讀 41,928評論 3贊 334
男人毒藥：我在死后第九天來索命
文/蒙蒙一、第九天我趴在偏房一處隱蔽的房頂上張望失都。院中可真熱鬧柏蘑，春花似錦、人聲如沸粹庞。這莊子的主人今日做“春日...
開封第一講書人閱讀 32,417評論 0贊 24
一樁弒父案，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽信粮。三九已至黔攒，卻和暖如春，著一層夾襖步出監(jiān)牢的瞬間强缘，已是汗流浹背督惰。一陣腳步聲響...
開封第一講書人閱讀 33,528評論 1贊 272
情欲美人皮
我被黑心中介騙來泰國打工，沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留旅掂，地道東北人赏胚。一個月前我還...
沈念sama閱讀 48,983評論 3贊 376
代替公主和親
正文我出身青樓，卻偏偏與公主長得像商虐，于是被迫代替她去往敵國和親觉阅。傳聞我的和親對象是個殘疾皇子，可洞房花燭夜當(dāng)晚...
茶點(diǎn)故事閱讀 45,573評論 2贊 359

最大似然估計(jì)（ML）和REML

推薦閱讀更多精彩內(nèi)容