軟體對照：同一個模型的各種程式寫法，以及疑難排解

各章以 R 配適模型，並以一個軟體對照方框作結，提供同一批模型在其他程式中的寫法。本附錄針對兩個模型更進一步：High School and Beyond 資料中含跨層級交互作用的隨機截距與斜率模型（第 6 章的模型 M3），以及 Riesby 資料的線性成長模型（第 19 章的模型 G1）。本附錄呈現每個模型在每一種程式中的寫法，比較各程式所得的估計值，說明為什麼看起來不同的輸出往往是相同的結果，最後提供一份指南，說明多層次軟體所印出的警告與錯誤訊息。R 的配適使用 lme4 (Bates et al., 2015)、lmerTest (Kuznetsova et al., 2017)、pbkrtest (Halekoh & Højsgaard, 2014)、nlme (Pinheiro & Bates, 2000)、glmmTMB (Brooks et al., 2017) 與 brms (Bürkner, 2017)，另有一個以 statsmodels (Seabold & Perktold, 2010) 進行的 Python 配適；這些都是為本附錄實際執行的。其他程式的語法則沒有執行，因為製作本書的環境中沒有這些程式；這些語法是根據近期版本的說明文件撰寫的（SAS 9.4 的 SAS/STAT、IBM SPSS Statistics 29、Stata 18、Mplus 8 與 HLM 8），使用時應該對照手邊版本的輸出加以檢查。以專書篇幅介紹特定程式中之多層次模型的著作，包括 SAS 方面的 Singer (1998) 與 Littell et al. (2006)、涵蓋多種程式的 West et al. (2014)，以及 Mplus User’s Guide (Muthén & Muthén, 1998--2017)。

E.1各程式及其預設值#

程式之間表面上的不一致，大多來自預設值，而不是來自不同的數學。表 E.1 針對每一種程式，列出預設的估計法、固定效果檢定的預設方法、隨機效果共變數矩陣的預設結構，以及 BIC 所用的樣本數。在兩種程式中以相同的估計法與相同的共變數結構配適同一個模型，應該能重現相同的估計值，直到小數點後好幾位；若不然，請檢查資料、平減方式與預測變項的編碼，接著檢查每一次配適是否都到達了概似的最大值（第 E.6 節）。

有三個預設值造成大部分的麻煩。第一，除非要求無結構矩陣，SAS 與 SPSS 會給隨機效果一個對角（變異成分）共變數矩陣；相對地，R 套件預設會估計隨機截距與隨機斜率的共變數，Stata 假定兩者獨立，Mplus 則把它們的共變數固定為零，除非以 WITH 釋放；因此，除非明確指定結構，含隨機斜率的模型在各程式之間並不相同。第二，Stata、glmmTMB、Mplus、Julia 與 MLwiN 預設以 ML 估計，而 lme4、nlme、SAS、SPSS、HLM、JASP、jamovi 與 statsmodels 則使用 REML；兩者的變異成分有系統性的差異（第 9 章）。第三，各程式以不同的方式檢定固定效果，從大樣本 $z$ 檢定，到採 Satterthwaite 或 Kenward–Roger 自由度的 $t$ 檢定都有（第 10 章），因此同一個估計值可能伴隨不同的 $p$ 值。

E.2同一個模型的各種程式寫法：HSB 跨層級模型#

第 6 章的模型 M3 以組平減的 SES、SES 的學校平均數、公私立別（天主教學校），以及兩個跨層級交互作用來預測數學成就，並含可以共變的隨機截距與 SES 隨機斜率：

資料是附錄 C 的 HSB 合併檔，其中包含組平減的 SES（ses_cwc）與樣本學校平均數（meanses），因此每一種程式都使用相同的變項（160 所學校的 7,185 名學生）。下列語法在程式允許時，一律要求 REML 與無結構共變數矩陣。

R。 四個套件使用相同的模型公式，語法只有些微差異：

library(lme4); library(lmerTest)          # lmerTest adds Satterthwaite df
m1 <- lmer(mathach ~ (meanses + catholic) * ses_cwc + (1 + ses_cwc | school),
           data = hsb, REML = TRUE)
summary(m1)                                # summary(m1, ddf = "Kenward-Roger")
library(nlme)
m2 <- lme(mathach ~ (meanses + catholic) * ses_cwc, random = ~ ses_cwc | school,
          data = hsb, method = "REML")
library(glmmTMB)
m3 <- glmmTMB(mathach ~ (meanses + catholic) * ses_cwc + (1 + ses_cwc | school),
              data = hsb, REML = TRUE)
library(brms)                              # Bayesian, default priors (Chapter 13)
m4 <- brm(mathach ~ (meanses + catholic) * ses_cwc + (1 + ses_cwc | school),
          data = hsb, seed = 20270001)

SAS。 使用 TYPE=UN 與 Satterthwaite 自由度的 PROC MIXED：

proc mixed data=hsb method=reml covtest;
  class school;
  model mathach = meanses catholic ses_cwc meanses*ses_cwc catholic*ses_cwc
        / solution ddfm=satterthwaite;
  random intercept ses_cwc / subject=school type=un;
run;

SPSS。 使用 COVTYPE(UN) 的 MIXED；WITH 之後的變項是共變項：

MIXED mathach WITH meanses catholic ses_cwc
  /FIXED=meanses catholic ses_cwc meanses*ses_cwc catholic*ses_cwc
  /METHOD=REML
  /PRINT=SOLUTION TESTCOV
  /RANDOM=INTERCEPT ses_cwc | SUBJECT(school) COVTYPE(UN).

Stata。 使用因子變項記法（factor-variable notation）、REML 與 Satterthwaite 自由度的 mixed：

mixed mathach c.meanses c.ses_cwc i.catholic c.meanses#c.ses_cwc ///
      i.catholic#c.ses_cwc || school: ses_cwc, ///
      covariance(unstructured) reml dfmethod(satterthwaite)

Mplus。 Mplus 沒有 REML；使用 ESTIMATOR = ML 時，其估計值對應於 R 套件的 ML 配適（REML = FALSE），而預設的 MLR 只改變標準誤，此時的標準誤會與其他程式以模型為基礎的標準誤不同。變項名稱是附錄 C 匯出壓縮檔中 Mplus 檔的名稱：

DATA:     FILE = hsb.dat;
VARIABLE: NAMES = school minority female ses mathach size catholic pracad
                  disclim himinty meanses_ n_j meanses ses_cwc;
          USEVARIABLES = mathach ses_cwc meanses catholic;
          WITHIN = ses_cwc;  BETWEEN = meanses catholic;
          CLUSTER = school;  MISSING = ALL (-999);
ANALYSIS: TYPE = TWOLEVEL RANDOM;  ESTIMATOR = ML;
MODEL:    %WITHIN%
          s | mathach ON ses_cwc;
          %BETWEEN%
          mathach s ON meanses catholic;
          mathach WITH s;

HLM。 HLM 8 以選單操作。多變量資料矩陣（MDM）由一個學生檔與一個學校檔建立，兩者都依學校排序，例如匯出壓縮檔中的 SPSS 檔 hsb_pupils.sav 與 hsb_schools.sav。模型逐一方程式設定：

其中 SES 以組平減的方式進入，並啟用 $u_{1j}$。壓縮檔中的學校檔包含隨資料發布的學校平均 SES（meanses_rb），它與此處使用的樣本平均數略有不同（附錄 C）；要完全重現表 E.2，必須把樣本平均數加入學校檔。

JASP 與 jamovi。 在 JASP 中，Mixed Models 模組的 Linear Mixed Models 分析以 mathach 為依變項，以 meanses、catholic 與 ses_cwc 為固定效果變項，並在模型項中加入兩個交互作用；以 school 為隨機效果的分組因子，含 ses_cwc 的隨機斜率，並勾選 Correlations 選項。在 jamovi 中，GAMLj 模組的 Linear Mixed Model 分析使用相同的變項，隨機效果輸入為 Intercept | school 與 ses_cwc | school，並設為相關。兩者都以 lme4 配適模型（JASP 透過 afex），但它們的預設值會改變參數化：JASP 以總和為零的對比為名義（nominal）預測變項編碼，並為每一個在學校內變動的固定效果納入隨機斜率；GAMLj 則把共變項平減於其總平均數。要重現表 E.2 的 lme4 欄，請把 catholic 以 0/1 的尺度（scale）變項輸入（或選擇處理編碼〔treatment coding〕），只保留隨機截距與 ses_cwc 斜率，並在 GAMLj 中把共變項的尺度調整設為 none；當只有編碼或平減方式不同時，隨機部分與其他係數會一致，但截距、SES 斜率，以及（在總和為零的編碼下）與天主教學校有關的項則不會一致。

Python 與 Julia。 在 Python 中，statsmodels 以一個公式與一個分組變項配適模型；在 Julia 中，MixedModels 使用 lme4 風格的公式，除非要求 REML，否則以 ML 估計：

import statsmodels.formula.api as smf
m = smf.mixedlm("mathach ~ (meanses + catholic) * ses_cwc", hsb,
                groups=hsb["school"], re_formula="~ses_cwc")
fit = m.fit(reml=True)

using MixedModels, CSV, DataFrames
hsb = CSV.read("hsb.csv", DataFrame)
m = fit(MixedModel, @formula(mathach ~ (meanses + catholic) * ses_cwc +
        (1 + ses_cwc | school)), hsb; REML = true)

估計值的比較。 表 E.2 呈現 lme4、nlme、glmmTMB 與 statsmodels 的 REML 估計值，以及 brms 的後驗平均數。四個以概似為基礎的配適彼此一致：它們在任何一個估計值上的最大差異約為 0.0004，REML 對數概似在四者中也都相同（$-2\log L_R$ = 46,503.66）。lme4、nlme 與 statsmodels 的 ML 配適也一致到相同的程度（最大差異 0.0005）。glmmTMB 與 statsmodels 的標準誤對某些係數略大，涉及 SES 的項差異最大（平均 SES 與 SES 的交互作用為 0.2996，相對於 0.2989），因為這兩個程式由所有參數聯合的觀察訊息矩陣求得標準誤，其中包含共變數參數的不確定性；而 lme4 與 nlme 則由 $({\mathbf{X}}^{\prime }{\hat{\mathbf{V}}}^{-1}\mathbf{X}{)}^{-1}$ 計算標準誤，把共變數參數當作已知；忽略這項不確定性時，glmmTMB 也得到 0.2989。

貝氏配適在固定效果上非常接近，其後驗標準差也接近 REML 標準誤，但在共變數參數上則不然：斜率變異數的後驗分配很寬，而且強烈偏斜（中位數 0.064，95% 可信區間 [0.000, 0.475]），因此其平均數 0.113 高於 REML 估計值 0.101，也遠高於它自己的中位數；共變數則被拉向零（0.123，相對於 0.192），因為資料對這麼小的變異數所提供的訊息很少（第 13 章）。每一個 $\hat{R}$ 都低於 1.01，所有參數中最小的 bulk 或 tail 有效樣本數為 775。

這張表也可以用來檢查其他程式的語法：SAS、SPSS 與 Stata 在使用 REML 與無結構共變數矩陣時，以及 HLM 在使用樣本學校平均數時，應該在其印出的精確度內重現表中的第一欄；而 Mplus 在使用 ESTIMATOR = ML 時，應該重現 ML 估計值（例如斜率變異數為 0.065，而不是 0.101）。

E.3Riesby 成長模型#

第 19 章的模型 G1 是 Riesby 資料中漢氏憂鬱量表評定的線性成長模型，含可以共變的隨機截距與週數隨機斜率，以 ML 從 66 位病人觀察到的 375 次評定估計。在長格式中，每一種程式設定這個模型的方式都和 HSB 模型一樣，以測量時間點為第一層、病人為第二層：

lmer(hamd ~ week + (1 + week | id), data = riesby, REML = FALSE)

proc mixed data=riesby method=ml covtest;
  class id;
  model hamd = week / solution ddfm=satterthwaite;
  random intercept week / subject=id type=un;
run;

mixed hamd week || id: week, covariance(unstructured)

在 SPSS 中，使用 /METHOD=ML 與 /RANDOM=INTERCEPT week | SUBJECT(id) COVTYPE(UN)；在 Mplus 中，使用與 HSB 相同的兩層次設定，以 week 為組內層次變項，或使用寬格式的潛在成長模型（第 22 章）；在 HLM 中，使用兩層次模型，WEEK 在第一層不平減。表 E.3 顯示，四個套件的結果彼此相差在 0.002 以內。

E.4核對不同程式的輸出#

當兩個程式的估計值不同時，通常問五個問題就能找到原因。

估計法相同嗎？ REML 與 ML 的變異成分估計值不同，最小的變異成分差異最大：HSB 模型的 SES 斜率變異數在 REML 下為 0.101，在 ML 下為 0.065。固定效果的差異則小得多。Mplus 與 Stata 預設以 ML 估計，因此其預設輸出必須與 ML 配適比較；Mplus 預設的 MLR 也會改變標準誤。

參數相同嗎？ 各程式以不同的形式報告隨機部分：lme4、glmmTMB 與 Julia 印出變異數、標準差與相關；nlme 印出標準差與相關；SAS、SPSS、Mplus、statsmodels 與 Stata（除非指定 stddeviations 選項）印出變異數與共變數；HLM 則兩種形式都印出。把標準差平方，並把相關乘以兩個標準差，就能把一種形式轉換成另一種形式。各程式對共變數矩陣元素的標示也不同：共變數 ${\tau }_{01}$ 在 SAS 中標示為 UN(2,1)，在 SPSS 中為 UN (2,1)，在 Stata 中為 cov(ses_cwc,_cons)，在 Mplus 中是組間層次的 MATHACH WITH S，在 HLM 中則是 $\tau$ 矩陣的非對角元素。

檢定相同嗎？ 固定效果相同，但其檢定不同。以 HSB 模型中天主教學校與 SES 的交互作用為例，Satterthwaite 近似法得出 143.2 個自由度，Kenward–Roger 近似法得出 142.0；而 nlme 因為這個交互作用含有第一層變項，就給它第一層的自由度（7,022），忽略了 SES 斜率在學校之間變動的事實；HLM 給隨機斜率的係數 $J-q-1$ 個自由度，因此給它 157。至於公私立別本身的效果，nlme 使用 157，Satterthwaite 則為 149.6。有 160 所學校時，結論相同，但叢集很少時，選擇就很重要（第 10 章）。

配適統計量相同嗎？ REML 對數概似不能與 ML 對數概似比較，固定效果不同之模型的 REML 概似則完全不能相互比較（第 9 至 11 章）。對 REML 配適，lme4 報告 REML 準則（HSB 模型的 $-2\log L_R$ = 46,503.66；以 REMLcrit() 取出，而不是 deviance(fit, REML = TRUE)）；對 ML 配適，則報告離異數（$-2\log L$ = 46,496.43）。訊息準則的懲罰項也不同：就 HSB 模型的 ML 配適而言，AIC 為 46,516.43，以學生數計算的 BIC 為 46,585.23，而像 SAS 那樣以學校數計算的 BIC 則為 46,547.18；Stata 以 estat ic, n(160) 得出後者，在其他程式中則可以由對數概似計算。只有在同一個程式中、以相同定義配適之模型間的差異才有意義。

資料相同嗎？ 變項有遺漏值時，各程式刪除不完整觀察值的方式不同，而平減也必須以完全相同的方式進行；就 Riesby 資料而言，所有程式都必須使用相同的 375 次觀察到的評定。上述 HLM 的例子提醒我們，名稱相同的變項，在這裡是 SES 的學校平均數，可能有不只一種計算方式。

E.5常見模型功能的語法對照#

表 E.4 列出本書模型最常見之功能在四種程式中的語法。各章的軟體對照方框提供完整的範例；這張表只供快速查閱。

E.6疑難排解#

多層次軟體的警告與錯誤訊息可以歸為少數幾類，每一類各有其常見原因與對策（表 E.5）。表 E.6 列出本附錄的範例在 R 與 Python 中重現的 15 則訊息，並照錄本書所用版本印出的確切文字。其中四個範例使用一筆小型模擬資料：10 個叢集、每個叢集五個觀察值，完全沒有叢集變異數；這正是無法估計隨機斜率的典型情況：lme4 報告邊界（奇異）配適，glmmTMB 報告海森矩陣非正定與奇異收斂，nlme 則以錯誤中止；第四個範例使用兩個互為倍數的預測變項，在奇異配適之外又多了一則關於秩不足（rank deficiency）的訊息。

HSB 模型的另外兩個情況顯示，警告不是出問題的唯一跡象，而沒有警告也不能保證沒有問題。以 ML 估計時，glmmTMB 對這個模型得出非正定的海森矩陣，斜率變異數為 0.018（logLik() 對這次配適傳回 NA；由其目標函數得到的對數概似為 $-$23,248.24，比最大值低 0.03），而 lme4、nlme 與 statsmodels 都得到 0.065；概似在這個小變異數上是平坦的，而這則警告正確地指出，最佳化演算法在到達最大值之前就停下來了。以 nlme 的替代最佳化演算法（lmeControl(opt = "optim")）用 REML 估計時，模型停在 0.0003 的斜率變異數，其 REML 對數概似比最大值低 0.56，卻沒有任何訊息。比較不同最佳化演算法或不同程式的對數概似，只需要一行程式碼，就能發現這兩個問題，前提是這些配適最大化的是同一個概似（相同的估計法，以及在廣義模型中相同的近似方法）。只使用 L-BFGS 最佳化演算法的 statsmodels 也過早停止，但它有發出訊息（表 E.6）；它預設的一連串最佳化演算法則到達了最大值。