概似、REML 與經驗貝氏的推導

第 9、11、12 與 23 章陳述了一些結果，若在正文中推導，會打斷敘述。本附錄以附錄 A 的符號，針對線性混合效果模型 ${\mathbf{y}}_j={\mathbf{X}}_j\mathbf{\gamma }+{\mathbf{Z}}_j{\mathbf{u}}_j+{\mathbf{e}}_j$ 得出這些推導，其中 ${\mathbf{u}}_j\sim \mathcal{N}(0,\mathbf{T})$、${\mathbf{e}}_j\sim \mathcal{N}(0,{\sigma }^2\mathbf{I})$，邊際共變數為 ${\mathbf{V}}_j={\mathbf{Z}}_j\mathbf{T}{\mathbf{Z}}_j^{\prime }+{\sigma }^2\mathbf{I}$；接著再把推導延伸到廣義線性混合效果模型。不帶下標 $j$ 的符號指的是全部叢集合在一起：$\mathbf{y}$ 與 $\mathbf{X}$ 把各 ${\mathbf{y}}_j$ 與 ${\mathbf{X}}_j$ 堆疊起來，$\mathbf{V}$ 則是以 ${\mathbf{V}}_j$ 為區塊的區塊對角矩陣（附錄 A 第 A.5 節）；$N={\sum }_j{n}_j$ 是觀察值的總數，$p$ 是固定效果的個數。大部分結果都在文中以數值驗證，所用的是 23 所學校之 NELS:88 子樣本中以作業時間預測數學成就的隨機截距模型 (Kreft & de Leeuw, 1998)；二元結果變項的近似方法，則以第 23 章的孟加拉避孕措施資料驗證 (Huq & Cleland, 1990)；第 B.5 節的演算法與第 B.8 節的邊界檢定，則依序在第 9 與 11 章以數值說明。

B.1邊際概似及其分數方程式#

令 $\mathbf{\theta }$ 集合各變異成分，也就是 $\mathbf{T}$ 的相異元素與 ${\sigma }^2$。由於各叢集彼此獨立，且 ${\mathbf{y}}_j\sim \mathcal{N}({\mathbf{X}}_j\mathbf{\gamma },{\mathbf{V}}_j)$，對數概似就是 $J$ 個多變量常態對數密度之和（附錄 A 第 A.8 節）：

它之所以稱為邊際概似，是因為隨機效果已經透過積分消去：$f({\mathbf{y}}_j)=\int f({\mathbf{y}}_j\mid {\mathbf{u}}_j)f({\mathbf{u}}_j)d{\mathbf{u}}_j$；對常態分配而言，這個積分就是共變數為 ${\mathbf{V}}_j$ 的常態密度。以附錄 A（第 A.9 節）的規則微分，就得到分數方程式。對固定效果而言，

而對每一個變異成分 ${\theta }_k$，令 ${\dot{\mathbf{V}}}_{jk}=\partial {\mathbf{V}}_j/\partial {\theta }_k$，殘差 ${\mathbf{r}}_j={\mathbf{y}}_j-{\mathbf{X}}_j\mathbf{\gamma }$，則

對隨機截距而言，${\dot{\mathbf{V}}}_{jk}$ 對 ${\tau }_0^2$ 是 ${\mathbf{J}}_{n_j}$，對 ${\sigma }^2$ 是 ${\mathbf{I}}_{n_j}$。在 $\mathbf{\gamma }$ 已知時，式 (B.3) 使觀察到的二次式等於其期望值：ML 估計值就是使殘差在這種加權的意義下，恰好如模型所說的那樣分散的參數值。就 NELS:88 模型而言，解析的對數概似 (B.1) 與 lme4 的值相差在 $10^{-12}$ 以內，而在 lme4 的 ML 估計值（${\hat{\tau }}_0^2=20.23$、${\hat{\sigma }}^2=71.14$）處，四個分數都在 $10^{-7}$ 以內等於零。

B.2GLS 與剖面概似#

對固定的 $\mathbf{\theta }$，式 (B.2) 是 $\mathbf{\gamma }$ 的線性方程式，其解為 GLS 估計式（附錄 A 第 A.7 節）$\hat{\mathbf{\gamma }}(\mathbf{\theta })=({\sum }_j{\mathbf{X}}_j^{\prime }{\mathbf{V}}_j^{-1}{\mathbf{X}}_j{)}^{-1}{\sum }_j{\mathbf{X}}_j^{\prime }{\mathbf{V}}_j^{-1}{\mathbf{y}}_j$。把它代入式 (B.1)，就得到剖面對數概似 ${\ell }_p(\mathbf{\theta })=\ell \{\hat{\mathbf{\gamma }}(\mathbf{\theta }),\mathbf{\theta }\}$，它只依賴變異成分；對 $\mathbf{\theta }$ 最大化 ${\ell }_p$，再計算 $\hat{\mathbf{\gamma }}(\hat{\mathbf{\theta }})$，就得到所有參數的 ML 估計值。軟體會把問題進一步化簡。lme4 以相對共變數因子 ${\mathbf{\Lambda }}_{\mathbf{\theta }}$（附錄 A 第 A.4 節）把共變數矩陣參數化為 $\mathbf{T}={\sigma }^2{\mathbf{\Lambda }}_{\mathbf{\theta }}{\mathbf{\Lambda }}_{\mathbf{\theta }}^{\prime }$，並把 ${\sigma }^2$ 也剖面消去，因此最佳化演算法只需在 ${\mathbf{\Lambda }}_{\mathbf{\theta }}$ 的元素上搜尋 (Bates et al., 2015)。

剖面化也能產生以概似為基礎的區間。對單一參數而言，該參數的剖面概似是對其他所有參數取最大值，而與最大值相差在 ${\chi }_{1;.95}^2/2=1.92$ 以內的參數值，構成 95% 概似比區間。圖 B.1 呈現學校標準差 ${\tau }_0$ 的剖面對數概似與剖面限制對數概似（第 B.3 節）。ML 曲線的峰值位於 ${\hat{\tau }}_0=4.50$，得出的區間為 [3.23, 6.46]，與 lme4 計算的剖面區間（[3.23, 6.46]）一致；REML 曲線的峰值位於 4.62，得出的區間為 [3.31, 6.68]。兩個區間在估計值周圍都不對稱，這正是概似函數的偏斜形狀所隱含的，也是為什麼對變異成分而言，剖面區間優於對稱的 Wald 區間（剖面區間：第 10 章；變異數之 Wald 推論的失效：第 11 章）。

註：NELS:88 子樣本，以作業時間預測數學成就的隨機截距模型。曲線：作為 ${\tau }_0$ 之函數的最大化對數概似（虛線，ML）與限制對數概似（實線，REML），各自以其最大值為基準，而 ${\sigma }^2$（以及 ML 的固定效果）已透過最大化消去。垂直線：ML（虛線）與 REML（實線）估計值。水平點線：$-1.92$，即 95% 概似比區間的切截值。

B.3限制最大概似估計#

變異數的 ML 估計式有向下的偏誤，因為它忽略了估計固定效果所用掉的自由度；在一般迴歸中，ML 把殘差平方和除以 $N$，REML 則除以 $N-p$（第 9 章）。REML 有兩種等價的推導 (Harville, 1974; Patterson & Thompson, 1971)。

誤差對比。 選取一個滿足 ${\mathbf{K}}^{\prime }\mathbf{X}=0$、行滿秩的 $N\times (N-p)$ 矩陣 $\mathbf{K}$。轉換後的資料 ${\mathbf{K}}^{\prime }\mathbf{y}\sim \mathcal{N}(0,{\mathbf{K}}^{\prime }\mathbf{V}\mathbf{K})$ 不依賴 $\mathbf{\gamma }$，REML 就是最大化這些資料的概似。兩個矩陣恆等式說明結果與 $\mathbf{K}$ 的選擇無關：對任何這樣的 $\mathbf{K}$，$\mathbf{K}({\mathbf{K}}^{\prime }\mathbf{V}\mathbf{K}{)}^{-1}{\mathbf{K}}^{\prime }=\mathbf{P}$，即附錄 A（第 A.7 節）的那個矩陣；而且 $\log |{\mathbf{K}}^{\prime }\mathbf{V}\mathbf{K}|=\log |{\mathbf{K}}^{\prime }\mathbf{K}|+\log |\mathbf{V}|+\log |{\mathbf{X}}^{\prime }{\mathbf{V}}^{-1}\mathbf{X}|-\log |{\mathbf{X}}^{\prime }\mathbf{X}|$，因此不同的 $\mathbf{K}$ 只會使對數概似改變一個常數 $-\frac{1}{2}\log |{\mathbf{K}}^{\prime }\mathbf{K}|$。若 $\mathbf{K}$ 的各行為正規正交，則 $|{\mathbf{K}}^{\prime }\mathbf{K}|=1$，因此

其中最後一項不涉及 $\mathbf{\theta }$，軟體通常將它省略；省略之後，預測變項的尺度改變時，REML 準則就會改變（第 10 章）。就 NELS:88 資料而言，由一個 519 $\times$ 517 的正規正交 $\mathbf{K}$ 直接計算 ${\mathbf{K}}^{\prime }\mathbf{y}$ 的對數概似，與不含該常數的公式相差 6.65，這正等於 $\frac{1}{2}\log |{\mathbf{X}}^{\prime }\mathbf{X}|$；在嘗試過的每一個 $\mathbf{\theta }$ 值上都是如此（誤差在 $10^{-11}$ 以內）。

積分概似。 以平坦的先驗分配把概似對固定效果積分，也會得到同一個函數，$L_R(\mathbf{\theta })\propto \int L(\mathbf{\gamma },\mathbf{\theta })d\mathbf{\gamma }$ (Harville, 1974)。由於 $\ell (\mathbf{\gamma },\mathbf{\theta })$ 是 $\mathbf{\gamma }$ 的二次函數，$\ell (\mathbf{\gamma },\mathbf{\theta })={\ell }_p(\mathbf{\theta })-\frac{1}{2}(\mathbf{\gamma }-\hat{\mathbf{\gamma }}{)}^{\prime }{\mathbf{X}}^{\prime }{\mathbf{V}}^{-1}\mathbf{X}(\mathbf{\gamma }-\hat{\mathbf{\gamma }})$，而第二項之指數函數的積分為 $(2\pi {)}^{p/2}|{\mathbf{X}}^{\prime }{\mathbf{V}}^{-1}\mathbf{X}{|}^{-1/2}$。因子 $|{\mathbf{X}}^{\prime }{\mathbf{V}}^{-1}\mathbf{X}{|}^{-1/2}$ 就是 REML 校正：它懲罰那些使固定效果估計得很精確的 $\mathbf{\theta }$ 值，藉此抵銷 $\hat{\mathbf{\gamma }}$ 以同一筆資料配適時所產生的過小殘差。在 $\hat{\mathbf{\gamma }}$ 處，式 (B.3) 的二次式等於 ${\mathbf{y}}^{\prime }\mathbf{P}{\dot{\mathbf{V}}}_k\mathbf{P}\mathbf{y}$，其中 ${\dot{\mathbf{V}}}_k$ 是以 ${\dot{\mathbf{V}}}_{jk}$ 為區塊的區塊對角矩陣；這個二次式的期望值是 $\mathrm{tr}(\mathbf{P}{\dot{\mathbf{V}}}_k)$，而不是 $\mathrm{tr}({\mathbf{V}}^{-1}{\dot{\mathbf{V}}}_k)$。REML 的分數方程式使用的是這個跡，因此它們是在估計 $\mathbf{\gamma }$ 之後，才使二次式等於其期望值。固定效果的 REML 估計值，就是在 ${\hat{\mathbf{\theta }}}_{\mathrm{REML}}$ 處的 GLS 估計值。由於 ${\ell }_R$ 是誤差對比 ${\mathbf{K}}^{\prime }\mathbf{y}$ 的概似，而 $\mathbf{K}$ 依賴 $\mathbf{X}$，固定效果不同的模型，其 REML 概似是不同資料的概似，因此不能互相比較（第 9 與 10 章）。

B.4變異成分的訊息與標準誤#

線性混合效果模型的期望訊息矩陣是區塊對角的：固定效果與變異成分的 ML 估計式在漸近上彼此獨立。$\mathbf{\gamma }$ 的區塊為 ${\sum }_j{\mathbf{X}}_j^{\prime }{\mathbf{V}}_j^{-1}{\mathbf{X}}_j$，其反矩陣就是 GLS 共變數矩陣；變異成分的區塊則由下列元素構成：

就 NELS:88 模型而言，在 ML 估計值處，式 (B.5) 的反矩陣得出 ${\hat{\tau }}_0^2$ 的標準誤 7.06 與 ${\hat{\sigma }}^2$ 的標準誤 4.52。Delta 法（附錄 A 第 A.9 節）把它們轉換為 ${\hat{\tau }}_0$ 的標準誤 0.79，以及殘差組內相關係數 $\hat{\rho }=.22$ 的標準誤 0.062。在 REML 下，式 (B.5) 的元素變成 $\frac{1}{2}\mathrm{tr}(\mathbf{P}{\dot{\mathbf{V}}}_k\mathbf{P}{\dot{\mathbf{V}}}_l)$；許多程式以觀察訊息，也就是常以數值方法計算的負海森矩陣，取代期望訊息；兩者在大樣本中一致。這裡的「大」指的是叢集數多：叢集層次參數的漸近結果，以叢集數 $J$ 不斷增加為前提，每個叢集有更多觀察值並不能取代更多的叢集。以這些標準誤建立的 Wald 區間，假定抽樣分配為常態且對稱；在 23 所學校的情形下，第 B.2 節的剖面區間比較可靠。

B.5演算法簡述#

牛頓-拉弗森法以對數概似的分數 $\mathbf{s}$ 與海森矩陣 $\mathbf{H}$ 進行迭代 ${\mathbf{\theta }}^{(t+1)}={\mathbf{\theta }}^{(t)}-{\mathbf{H}}^{-1}\mathbf{s}$；費雪計分法以期望訊息 (B.5) 取代 $-\mathbf{H}$，期望訊息較容易計算，而且只要變異成分可以識別就是正定的，$-\mathbf{H}$ 在遠離最大值處則不一定正定；因此，費雪計分法在最初幾次迭代中往往比較穩定。Goldstein 的迭代廣義最小平方法交替進行兩個步驟：以 GLS 估計 $\mathbf{\gamma }$，以及把殘差的平方與交乘積對 $\mathbf{V}$ 的設計做 GLS 迴歸 (Goldstein, 1986)。EM 演算法把隨機效果視為遺漏資料 (Dempster et al., 1977; Laird & Ware, 1982)：E 步驟在目前的參數值下計算 $E({\mathbf{u}}_j\mid {\mathbf{y}}_j)$ 與 $\mathrm{Var}({\mathbf{u}}_j\mid {\mathbf{y}}_j)$（第 B.6 節）；M 步驟則把 $\mathbf{T}$ 更新為 $E({\mathbf{u}}_j{\mathbf{u}}_j^{\prime }\mid {\mathbf{y}}_j)$ 的平均，把 ${\sigma }^2$ 更新為第一層殘差平方之期望值的平均，並以最小平方法對 ${\mathbf{y}}_j-{\mathbf{Z}}_{j}E({\mathbf{u}}_j\mid {\mathbf{y}}_j)$ 更新 $\mathbf{\gamma }$（第 9 章），或者在一種常見的變形中改用 GLS。EM 從不使概似下降，但在最大值附近收斂得很慢；Lindstrom 與 Bates (1988) 的結論是，在大多數情況下，實作良好的牛頓-拉弗森演算法較為可取。lme4 使用稀疏 Cholesky 分解，以不需導數、有邊界限制的最佳化演算法，在 ${\mathbf{\Lambda }}_{\mathbf{\theta }}$ 的元素上最小化剖面離異數 (Bates et al., 2015)；第 9 章在 NELS:88 的虛無模型上執行牛頓-拉弗森法、費雪計分法與 EM（表 9.3），並討論估計問題及其對策。

B.6經驗貝氏預測、BLUP 與 Henderson 方程式#

隨機效果的後驗分配。 把 ${\mathbf{u}}_j\sim \mathcal{N}(0,\mathbf{T})$ 視為先驗分配、${\mathbf{y}}_j\mid {\mathbf{u}}_j\sim \mathcal{N}({\mathbf{X}}_j\mathbf{\gamma }+{\mathbf{Z}}_j{\mathbf{u}}_j,{\sigma }^2\mathbf{I})$ 視為概似，貝氏定理得出常態的後驗分配，因為兩個密度之乘積的指數是 ${\mathbf{u}}_j$ 的二次式。整理二次項與一次項，可得

由 Woodbury 恆等式（附錄 A 第 A.5 節），後驗平均數等於 $\mathbf{T}{\mathbf{Z}}_j^{\prime }{\mathbf{V}}_j^{-1}({\mathbf{y}}_j-{\mathbf{X}}_j\mathbf{\gamma })$，後驗變異數等於 $\mathbf{T}-\mathbf{T}{\mathbf{Z}}_j^{\prime }{\mathbf{V}}_j^{-1}{\mathbf{Z}}_j\mathbf{T}$，也就是附錄 A（第 A.8 節）的條件常態結果。以估計值取代 $\mathbf{\gamma }$ 與 $\mathbf{\theta }$，就得到第 12 章的經驗貝氏預測值。對隨機截距而言，後驗平均數就是收縮估計式 ${\lambda }_j({\bar{y}}_j-{\bar{\mathbf{x}}}_j^{\prime }\mathbf{\gamma })$，其中 ${\lambda }_j=n_j{\tau }_0^2/({\sigma }^2+n_j{\tau }_0^2)$。

BLUP 與混合效果模型方程式。 不假定常態性時，把後驗平均數中的 $\mathbf{\gamma }$ 換成其 GLS 估計值所得的 $\mathbf{T}{\mathbf{Z}}_j^{\prime }{\mathbf{V}}_j^{-1}({\mathbf{y}}_j-{\mathbf{X}}_j\hat{\mathbf{\gamma }})$，在 $\mathbf{\theta }$ 已知時就是最佳線性不偏預測值（best linear unbiased predictor, BLUP）：在預測誤差的平均數對每一個 $\mathbf{\gamma }$ 都為零的所有 $\mathbf{y}$ 之線性函數中，它的均方預測誤差最小 (Henderson, 1975; Robinson, 1991)。$\mathbf{\theta }$ 為估計值時，它是經驗 BLUP（第 12 章）。Henderson 藉由最大化 $\mathbf{y}$ 與 $\mathbf{u}$ 的聯合密度，同時求得 $\mathbf{\gamma }$ 的估計值與 $\mathbf{u}=({\mathbf{u}}_1^{\prime },\dots ,{\mathbf{u}}_J^{\prime }{)}^{\prime }$ 的預測值；這等價於最小化懲罰平方和 ${\sigma }^{-2}\Vert \mathbf{y}-\mathbf{X}\mathbf{\gamma }-\mathbf{Z}\mathbf{u}{\Vert }^2+{\mathbf{u}}^{\prime }{\mathbf{G}}^{-1}\mathbf{u}$，其中 $\mathbf{G}={\mathbf{I}}_J\otimes \mathbf{T}$，$\mathbf{Z}=\mathrm{diag}({\mathbf{Z}}_1,\dots ,{\mathbf{Z}}_J)$。令導數為零，就得到混合效果模型方程式

解第二個區塊得 $\hat{\mathbf{u}}=({\mathbf{Z}}^{\prime }\mathbf{Z}+{\sigma }^2{\mathbf{G}}^{-1}{)}^{-1}{\mathbf{Z}}^{\prime }(\mathbf{y}-\mathbf{X}\hat{\mathbf{\gamma }})$，代入第一個區塊，再利用 Woodbury 恆等式（附錄 A 第 A.5 節），得到 ${\mathbf{X}}^{\prime }{\mathbf{V}}^{-1}\mathbf{X}\hat{\mathbf{\gamma }}={\mathbf{X}}^{\prime }{\mathbf{V}}^{-1}\mathbf{y}$；因此 $\hat{\mathbf{\gamma }}$ 就是 GLS 估計式，而 $\hat{\mathbf{u}}$ 就是在 $\hat{\mathbf{\gamma }}$ 處計算的後驗平均數 (B.6)。乘上 ${\sigma }^2$ 之後，對隨機截距而言，懲罰項為 $({\sigma }^2/{\tau }_0^2){\sum }_j{u}_j^2$，也就是第 32 章懲罰參數為 ${\sigma }^2/{\tau }_0^2$ 的脊懲罰。式 (B.7) 中係數矩陣的反矩陣，其左上區塊包含 $\hat{\mathbf{\gamma }}$ 的共變數矩陣，右下區塊則包含預測誤差變異數 $\mathrm{Var}(\hat{\mathbf{u}}-\mathbf{u})$；與式 (B.6) 的後驗變異數不同，這些預測誤差變異數納入了 $\hat{\mathbf{\gamma }}$ 的不確定性。

就 NELS:88 模型而言，在 REML 變異成分下，式 (B.7) 的解與 lme4 的固定效果、其共變數矩陣以及學校效果預測值，差異都在 $10^{-12}$ 以內。對最小的學校（$n_j=5$），lme4 報告的後驗標準差為 2.92，預測誤差標準差為 2.99；對最大的學校（$n_j=67$），兩者為 1.01 與 1.47。對隨機截距而言，$\hat{\mathbf{\gamma }}$ 的不確定性使後驗變異數增加 ${\lambda }_j^2{\bar{\mathbf{x}}}_j^{\prime }\mathrm{Var}(\hat{\mathbf{\gamma }}){\bar{\mathbf{x}}}_j$；這個增量往往在那些以自身資料最精確地決定其效果的學校中最大，因為它們的預測值收縮得最少。

B.7廣義線性混合效果模型：拉普拉斯近似、數值積分與 PQL#

在 GLMM 中，給定 ${\mathbf{u}}_j$ 時 ${\mathbf{y}}_j$ 的條件分配可以是伯努利、二項或卜瓦松分配等；叢集 $j$ 的邊際概似

沒有封閉形式。常用的近似方法有三種（第 23 章）。

拉普拉斯近似。 把 $h_j$ 在其最大值 ${\hat{\mathbf{u}}}_j$（即條件眾數）附近展開到二階：$h_j(\mathbf{u})\approx h_j({\hat{\mathbf{u}}}_j)-\frac{1}{2}(\mathbf{u}-{\hat{\mathbf{u}}}_j{)}^{\prime }{\mathbf{H}}_j(\mathbf{u}-{\hat{\mathbf{u}}}_j)$，其中 ${\mathbf{H}}_j=-{\partial }^2{h}_j/\partial \mathbf{u}\partial {\mathbf{u}}^{\prime }$ 在 ${\hat{\mathbf{u}}}_j$ 處計算。剩下的積分是常態積分，因此若以 $q$ 表示 ${\mathbf{u}}_j$ 的維度，則

當被積函數與常態密度成比例時，這個近似是精確的；當被積函數偏斜時，近似就會變差，例如每個叢集只有少數幾個二元觀察值的情形。它的漸近偏誤在二元與次序結果變項中比在計數結果變項中大，並隨叢集大小增加而減少 (Joe, 2008)。

適應性高斯積分法。 數值積分法以 $Q$ 個節點上的加權和取代積分。非適應性積分法對每個叢集使用相同的節點，例如 mirt 的固定網格（第 30 章），而且需要很多節點。適應性積分法把節點的中心放在 ${\hat{\mathbf{u}}}_j$，並以 $\sqrt{2}{\mathbf{H}}_j^{-1/2}$ 調整其尺度，使節點落在被積函數有質量的地方；只用一個節點時，它就化約為拉普拉斯近似，而誤差隨 $Q$ 增加而迅速減少 (Pinheiro & Bates, 1995; Rabe-Hesketh et al., 2002)。函數求值的次數隨 ${\mathbf{u}}_j$ 的維度 $q$ 以 $Q^q$ 增加，而 glmer 只對含單一純量隨機效果項的模型提供數值積分法。

懲罰準概似。 PQL 在目前的預測值附近把模型線性化，對一個工作變項配適線性混合效果模型，並反覆迭代直到收斂 (Breslow & Clayton, 1993)。它的速度快，也能推廣到複雜的隨機效果結構；但在二元結果變項、叢集小且變異數大時，它會低估變異成分，並使勝算對數尺度上的固定效果往零衰減。Rodríguez 與 Goldman (1995) 已就密切相關的 MQL 指出這一點，Rodríguez 與 Goldman (2001) 則就 PQL 指出同樣的結果，第 23 章的模擬也證實了這一點。

範例。 圖 B.2 呈現避孕措施資料中第 41 號地區（26 位婦女，其中 13 位使用避孕措施）的 $\exp \{h_j(u)-h_j({\hat{u}}_j)\}$；在所有地區中，拉普拉斯近似對這個地區之概似貢獻的重現最差。即使在這裡，兩條曲線也無法分辨，因此這張圖呈現的是近似如何運作，而不是近似如何失敗：這個地區的拉普拉斯誤差為 0.3%（各地區的中位數為 0.2%），而五個適應性積分節點就能把相對誤差降到 $10^{-6}$ 以下。表 B.1 比較整個模型在各種近似下的結果。拉普拉斯近似的 $-2\log L$ 比 25 個節點之積分法的值大 0.27，其中幾乎全部是近似誤差，而不是適配上的差異（因此，不應比較以不同近似方法求得的概似；第 23 與 24 章）；地區標準差則從 0.474 變為 0.477。三個節點與收斂後的 $-2\log L$ 相差不到 0.02，而從七個節點起，在表中所示的精確度下就不再有任何變化。PQL 得出較小的標準差（0.464）與較小的都市係數（0.681 對 0.692），方向與其已知的偏誤一致。像這裡一樣，地區大小的中位數為 26、變異數又小時，拉普拉斯近似就已足夠。在每人五個二元試題的情形下，它把一個模擬之 Rasch 模型的個人標準差低估了約 5%（第 30 章）；在叢集大小為二時，對於 0.5 到 5 的真實變異數，它把變異數低估了 26% 到 33%（第 23 章）。

註：粗淺色線：$\exp \{h_j(u)-h_j({\hat{u}}_j)\}$，即一個地區之式 (B.8) 被積函數除以其最大值，在固定效果與地區變異數的拉普拉斯估計值處計算；細虛線：拉普拉斯近似的常態曲線，眾數與曲率都與前者相同。圓形：適應性高斯積分法的五個節點。

B.8邊界上的概似比檢定#

概似比統計量通常所用的卡方分配，建立在對數概似於參數空間內部某個真值附近的二次近似之上。當虛無假設把某個變異數設為零時，真值就位於邊界上，而估計式受到限制，不得為負。在大樣本中，不受限制的估計值近似以零為中心的常態分配；它有一半的機率為負，此時受限制的估計值為零，概似比統計量也為零；另有一半的機率為正，此時統計量是一個近似標準常態變項的平方。這個啟發式論證假定只有一個參數位於邊界上、其他參數都在內部，由此得到 ${\chi }_0^2$ 與 ${\chi }_1^2$ 以 50:50 組成的混合分配 (Self & Liang, 1987; Stram & Lee, 1994)。檢定一個與截距相關的隨機斜率時，它的變異數以及它與截距的共變數在虛無假設下都為零，這時得到的則是 ${\chi }_1^2$ 與 ${\chi }_2^2$ 以 50:50 組成的混合分配；臨界值與其他情形見附錄 A（第 A.12 節），應用見第 11 章。這個漸近論證同樣以叢集數為依據；在有限樣本中，混合分配往往偏保守，因此較安全的做法是對概似比統計量做參數拔靴法，或者，對只含單一變異成分的模型，使用其精確的虛無分配 (Crainiceanu & Ruppert, 2004)（第 11 章）。

延伸閱讀#

Harville (1977); Searle et al. (1992)

變異成分的 ML 與 REML 估計。

Robinson (1991)

BLUP、Henderson 方程式及其多種詮釋。

Pinheiro 與 Bates (1995); Rabe-Hesketh et al. (2002)

拉普拉斯近似與適應性積分法。

Verbeke 與 Molenberghs (2000)

用於縱貫性資料的線性混合效果模型，並附推導。