為個體內變異性建模：混合效果的位置尺度模型

16.1 變異性就是訊號

本章的前提是一組在平均數上完全相同、在其他一切要緊的事情上都不同的人。圖 16.1 顯示四個人，他們的負向情緒 (negative affect) 平均數相同，個體內變異性卻天差地遠，從幾乎持平到劇烈擺盪。在臨床情境中，這是四個不同的病人：那個穩定的與那個高度易變的可能共享同一個診斷與同一個平均症狀水準，卻面對不同的風險、對不同的治療有反應，也需要不同的監測。研究這類差異的傳統，貫穿了精神病理學中關於情緒不穩定的工作 (Ebner-Priemer et al., 2009; Jahng et al., 2008)、把表現不一致當成認知衰退 (cognitive decline) 早期標記 (early marker) 的研究 (Hultsch et al., 2002)，以及把短期個體內變異性當成一個穩定的個別差異構念 (individual-differences construct) 的更大方案 (Wang et al., 2012)。在每一種情形下，主張都是同一個：變異數就是現象，而一個只為平均數建模的方法把它丟掉了。

註：四個模擬出來的人，負向情緒的平均數相同，個體內變異性由穩定遞增到高度易變。平均數分辨不出他們，個體內標準差可以，而後者經常正是臨床上感興趣的量。本圖延伸第 7 章的個體內變異性圖示，並加上明確的實質解讀。

第 7 章以個體內標準差 (intraindividual standard deviation, iSD) 測量這種變異性，一次算一個人，並警告這個指標很脆弱。這項警告現在可以寫成一份四點的缺陷清查。第一，兩階段的標準差在估計時帶著誤差，而誤差大小取決於時點數，因此一個被測了二十次提示 (beep) 的人，其估計值遠比一個被測了五十次的人來得不穩，圖 16.2 顯示了由此產生的散布。第二，沒有收縮 (shrinkage)：每個人的估計值各自為政，於是那些極端值，而它們又不成比例地落在時點數最少的人身上，被照單全收，而不是被拉向團體。第三，只要水準與變異性有關聯，這個指標就把兩者混淆 (confounding) 在一起，而圖 16.2 的右圖顯示的正是一個有界 (bounded) 或容易碰到地板 (floor) 的量尺所引發的、平均數與標準差之間的正向耦合 (coupling)，於是變異性上的差異可能只是水準差異的假象。第四，當這個指標接著被帶進第二個分析當預測變項時，它的估計誤差就成了變項含誤差 (error-in-variables)，會使下游的係數產生偏誤，這是「把一個帶雜訊的估計值當成已知量」時一再出現的代入謬誤 (plug-in fallacy) 的一個實例。一個以模型為基礎的做法必須把四項都修好：它必須跨人借力 (borrow strength)、必須容許兩個層次都能放共變項 (covariate)、必須為水準與變異性聯合建模，並且必須把屬於個人的變異性中的不確定性，傳遞到任何用到它的地方去。表 16.1 把兩種做法並排放在一起。

註：左：原始的個體內標準差對上提示次數，顯示時點越少，估計越吵、越極端，而且沒有任何收縮。右：原始的個體內標準差對上個人的平均數，顯示結果變項有地板時，水準與變異性混淆在一起。這兩項缺陷都由一個「帶收縮地估計個人波動性、並把它與平均數分開」的模型修補。

16.2 分三步組裝這個模型

混合效果的位置尺度模型最好分三步組裝，每一步都是前一步一眼認得出來的延伸。第一步是異質的 (heterogeneous) 第一層變異數。一般的混合模型假定每個人共用單一個殘差變異數 (residual variance) \(\sigma^2\)；第一次鬆綁讓這個變異數隨組別或隨共變項而不同，也就是第 14 章殘差變異數函數 (variance function) 那塊熟悉的地盤。第二步把變異數變成一個迴歸。與其為每一組各給一個變異數，一個對數線性的尺度次模型 (log-linear scale submodel) 把個體內變異數的對數寫成預測變項的線性函數，\(\log \sigma^2_{it} = \alpha_0 + \alpha_1 x_{it}\)，於是一個共變項可以抬高或壓低當下的變異性。對數連結 (log link) 不是偶然的：它讓變異數在係數取任何值時都保持為正，也讓係數有一個乘性 (multiplicative) 的讀法，因為預測變項變動一個單位，變異數就乘以 \(e^{\alpha_1}\)。第三步在尺度次模型上加一個隨機效果。人們在殘差波動性上的差異，超出共變項所能解釋的範圍，而一個隨機尺度效果 (random scale effect) \(\omega_i\) 就捕捉這份差異，\(\log \sigma^2_{it} = \alpha_0 + \alpha_1 x_{it} + \omega_i\)，於是每個人都有一個屬於自己的基線波動性，正如位置次模型給每個人一個屬於自己的基線水準。

完整的模型把兩個次模型耦合起來，而且關鍵在於容許它們的隨機效果相關。位置次模型是 \(y_{it} = \gamma_{0} + \gamma_{1} x_{it} + u_{i} + \varepsilon_{it}\)，其中 \(\varepsilon_{it} \sim N(0, \sigma^2_{it})\)；尺度次模型是 \(\log \sigma^2_{it} = \alpha_0 + \alpha_1 x_{it} + \omega_i\)；而兩個屬於個人的效果是一起抽出來的，\((u_i, \omega_i)' \sim N(\mathbf{0}, \mathbf{\Sigma})\)，兩者之間有一個相關 \(\rho\)。圖 16.3 畫出這個架構。相關 \(\rho\) 是最有實質分量的參數：它問的是平均起來比較高的人是不是也比較容易起伏，例如平均負向情緒較高的病人是不是在情緒上也比較波動，而這是任何兩階段分析都問不出來的問題，因為它從來沒有在同一個模型裡估計過這兩個量。表 16.2 是這些參數以及各自讀法的一份詞彙表。

註：位置次模型在平均數上放一個屬於個人的隨機截距 \(u_i\)；尺度次模型在個體內變異數的對數上放一個屬於個人的隨機截距 \(\omega_i\)；一個個體內的預測變項可以進任一個次模型。兩個隨機效果由一個帶相關 \(\rho\) 的聯合分配抽出，而這個位置尺度相關把一個人的水準與這個人的波動性連在一起。

這個模型比一般的混合模型難估計，因為概似 (likelihood) 要對兩個相關的隨機效果積分，而其中一個還坐在一個變異數裡面，且這個積分沒有封閉解。有兩條路可走。以適應性求積 (adaptive quadrature) 做最大概似 (maximum likelihood)，實作在 Hedeker 的專用程式裡 (Hedeker et al., 2008; Hedeker & Nordgren, 2013)，這些程式為生態瞬時評估 (ecological momentary assessment, EMA) 資料開創了這個模型。以馬可夫鏈蒙地卡羅 (Markov chain Monte Carlo) 做貝氏估計，也就是本書採用的路線，則能更優雅地推廣 (generalize) 到隨機尺度斜率 (random scale slope)、兩層以上的結構，以及非高斯 (non-Gaussian) 的位置次模型 (Rast et al., 2012; Williams et al., 2019)。下面的實作分析以二維 Gauss-Hermite 求積 (quadrature) 的最大概似配適，這是一個能還原同一個估計標的 (estimand)、且透明的實作 (implementation)；而實務上建議使用的工具是貝氏的那一個，其做法由第 17 章展開。這是刻意安排的交棒：位置尺度模型促成了轉向貝氏的一步，而下一章把它交出來。

16.3 配適這個模型

實作範例問的是：在控制水準之後，當下的壓力能不能預測負向情緒的波動性？用的是模擬的日誌資料集 ema_lability，其中 \(120\) 個人被測到的提示 (beep) 次數不等，而資料生成歷程 (data-generating process) 中含有一個真正的隨機尺度效果與一個正向的位置尺度相關。個體內的壓力同時進入兩個次模型：進位置次模型是為了讓壓力能抬高負向情緒的水準，進尺度次模型是為了讓它能抬高波動性。壓力的位置效果 (location effect) 被還原為 \(0.53\)，接近它的真值，但有意思的發現在尺度次模型。壓力的尺度係數是 \(0.41\)，標準誤 (standard error) \(0.03\)，在乘性的尺度上意思是：一個人當下的壓力每上升一個單位，他負向情緒的個體內變異數就乘以 \(e^{0.41} = 1.50\)。圖 16.4 把這件事畫成一條隨壓力上升的個體內標準差預測值，並附上不確定性帶 (uncertainty band)。壓力不只抬高負向情緒的平均數，它還讓負向情緒變得不穩定，而尺度次模型正是讓第二個、與第一個截然不同的主張變得可估計的東西。

註：由尺度次模型算出的負向情緒個體內標準差預測值，隨個體內壓力變化，附 90% 的不確定性帶。當下的壓力每上升一個單位，殘差變異數就乘以約 \(1.5\)。尺度次模型把這個「使之不穩定」的效果，與壓力對平均水準的效果分開估計。

模型對兩階段指標的修補，是把屬於個人的波動性估計值收縮起來。圖 16.5 把每個人原始的個體內對數標準差，與模型的經驗貝氏 (empirical Bayes) 估計值對照，橫軸是提示次數。原始估計值很極端，時點數最少的人尤其如此，而模型依它們的不可靠 (unreliable) 程度按比例把它們拉向團體，正是第 13 章的收縮，如今施用在一個變異數上。收縮後的估計值不只是比較平滑而是真的比較好，這件事在這裡是查得出來的，因為真值已知：模型對個人波動性的估計值與真正的隨機尺度效果相關 \(0.90\)，原始指標則是 \(0.81\)。位置尺度相關這個實質目標估計為 \(0.40\)，標準誤 \(0.10\)，而圖 16.6 顯示它背後那些屬於個人的效果的聯合分配：平均負向情緒較高的人，在情緒上確實傾向於比較波動。那個很寬的標準誤本身就是一課。位置尺度相關的估計精確度遠低於一個平均數效果，因為它是一個二階 (second-order) 的量，取決於各人波動性的分散程度，所以以它為估計目標的研究，需要比只瞄準平均數的研究更多的人與更多的時點，這是一項與第 4 章相接的檢定力 (power) 事實，也是誠實的報告必須承認的。

註：每個人原始的個體內對數標準差（灰）與模型的經驗貝氏估計值（藍），以一條線段相連，橫軸為提示次數。模型把時點少、雜訊大的估計值向團體平均收縮，而因為資料生成的真值已知，可以顯示收縮後的估計值比原始指標更能還原真正的個人波動性（相關 \(0.90\) 對 \(0.81\)）。

註：經驗貝氏的位置隨機效果（個人平均數）與尺度隨機效果（個人的對數波動性）的聯合分配。那條正斜率就是估計出來的位置尺度相關，這裡約為 \(0.40\)，標準誤 \(0.10\)。這個相關是任何兩階段分析都取用不到的，而它通常正是一項位置尺度研究的實質目標。

16.4 延伸與它的親戚

尺度次模型接受任何迴歸都接受的那些精緻化，其中兩項特別有用。變異性的預測變項可以是時間本身，這使模型變成一項關於「波動性在研究期間如何改變」的陳述。圖 16.7 顯示一條變異數軌跡，其中個體內波動性隨治療的週次遞減，由一個時間變動的尺度模型還原出來：變異性在這裡本身就是一項治療結果，是症狀的穩定化 (stabilization) 而不是症狀水準的降低，而一項治療可能在第二件事上成功，即使它在第一件事上進展緩慢。個體間與個體內變異數之分是第二項精緻化。隨機尺度效果為波動性上的個體間差異建模，尺度共變項則為它在個體內、時點與時點之間 (occasion-to-occasion) 的改變建模，而 Hedeker 的架構明確地把兩者分開 (Hedeker et al., 2012)。為個體內變異的個別差異而設的多層次做法早於這個完全耦合的模型，至今仍是一個有用的特例 (Hoffman, 2007)。

註：一個時間變動的尺度模型，其中個體內殘差標準差隨治療週次遞減。配適出來的衰退（藍）貼著真正的生成衰退（虛線）。當變異性本身就是結果變項時，它在一項介入期間的軌跡就是研究發現，與平均水準上的任何改變截然不同。

位置尺度模型坐落在一族分配模型 (distributional model) 之中，這族模型讓平均數以外的參數也隨共變項變化。表 16.3 畫出這一帶的地圖。雙重廣義線性模型 (double generalized linear model) 在廣義線性模型的離散度上放一個次模型；為位置、尺度與形狀而設的廣義加法模型 (generalized additive models for location, scale and shape) 讓一個有彈性的分配的每一個參數，包括偏態 (skewness) 與峰度 (kurtosis)，都能是共變項的平滑函數 (Rigby & Stasinopoulos, 2005)；變異數函數混合模型與現代軟體的離散度 (dispersion) 公式能快速配適固定的 (fixed) 尺度次模型；而多層次教育研究中的變異數建模傳統，則從另一個領域處理同一個結構 (Leckie et al., 2014)。它們共享同一個前提：值得建模的不只是中心，還有分散程度 (spread)。

最後一項延伸關乎把估計出來的變異性拿到下游使用。當一個人的尺度效果 \(\omega_i\) 要在後續分析中當預測變項或結果變項時，它不可以被縮成一個點估計值 (point estimate) 再代進去，因為那會丟掉它的不確定性、使第二階段的結果產生偏誤，正是任何「把收縮後的估計值當成資料」時都會犯的那個代入謬誤。正確的做法是把 \(\omega_i\) 的完整後驗分配帶到下游，途徑是把兩個階段聯合配適，或是傳遞後驗抽樣 (posterior draws)，好讓一個波動性被估得很不準的人，按比例地貢獻得少一些。這也是貝氏路線在這裡格外自然的又一個理由，因為它產生的後驗抽樣，讓誠實的傳遞變得直截了當。

16.5 驗證、陷阱與報告

一個位置尺度模型必須在它下主張的地方受檢查，也就是在變異數結構上，而工具是一個瞄準變異數特徵 (variance feature)、而不是瞄準平均數的後驗預測檢查 (posterior predictive check)。圖 16.8 把觀察到的、每個人的個體內標準差分配，與由配適好的模型複製出來的分配對照；觀察到的分配穩穩地落在複製出來的分散範圍之內，這是模型重現了它本來就是為了解釋而建的那個波動性型態的證據。這一類檢查，瞄準的是 iSD 跨人的分散程度或屬於個人的變異數是否校準 (calibration) 得宜，是一般後驗預測檢查在位置尺度模型上的對應物，應該放進任何一份報告裡。表 16.4 是報告檢核表，它在第 13 章的混合模型檢核表之上，再加上貝氏配適所需要、而由第 17 章完整說明的估計品質診斷 (diagnostics)。

註：觀察到的、每個人的個體內標準差分配（紅），對上由配適好的位置尺度模型複製出來的四十條分配（灰）。觀察到的分配落在複製出來的分散範圍之內，表示模型重現了波動性的個體間型態。位置尺度模型的後驗預測檢查應該瞄準變異數特徵，而不是平均數。

有三個陷阱一再出現。把尺度係數當成原始標準差上的效果來解讀，忽略了對數連結，會低估這個效果；尺度係數是變異數上的乘數 (multiplier)。從一條有界的量尺讀出「低變異性」可能是假象，因為一個靠近地板或天花板的人本來就沒什麼空間可以變動，所以表面上波動性的差異可能是喬裝過的水準差異，這與第 7 章討論相對變異性指數時提出的告誡相同。而把一個人波動性的點估計值當成結果變項使用，等於把它的不確定性洗掉 (launder away)，並使下游分析產生偏誤。三者都能以這個模型本身提供的紀律避開。

16.6 在 R 中配適位置尺度模型

brms 的分配模型語法把兩個次模型寫在一條式子裡，第一個給平均數，第二個冠上 sigma 前綴、給對數標準差，兩者各有自己的固定效果 (fixed effects) 與隨機效果。

library(brms)
# 位置：na ~ stress + (1 | person)；尺度：log sigma ~ stress + (1 | person)
f <- bf(na    ~ stress_wc + (1 | p | person),
        sigma ~ stress_wc + (1 | p | person))   # 共用的 "p"：讓兩個 RE 相關
m <- brm(f, data = ema, chains = 4, cores = 4)   # 先驗 + MCMC；見第 17 章
summary(m)                                        # 尺度效果在 log-SD 上

共用的分組標籤 | p | 就是讓位置與尺度兩個隨機效果相關起來的東西，由此得到位置尺度相關。一個只有固定效果、沒有隨機尺度效果的尺度模型，可以更快地透過 nlme 的殘差變異數函數或 glmmTMB 的離散度公式取得，當成踏腳石 (stepping stone)、或在不需要隨機尺度時很有用。

library(nlme)   # 固定的對數線性尺度：殘差 SD ~ exp(theta * stress)
m_fixed <- lme(na ~ stress_wc, random = ~1 | person, data = ema,
               weights = varExp(form = ~ stress_wc))
# glmmTMB 的替代做法：glmmTMB(na ~ stress_wc + (1|person),
#                              dispformula = ~ stress_wc, data = ema)

隨附的腳本 ch16_analysis_V01.R 收錄了完整的分析，包括兩階段缺陷的示範、以二維 Gauss-Hermite 最大概似配適並取經驗貝氏個人效果的完整位置尺度模型、尺度效果與位置尺度的圖示、變異數軌跡，以及瞄準變異數的後驗預測檢查；圖形則由中文版的 ch16_figures_zh_V01.R 繪出，資料集由 gen_ema_lability_V01.R 產生。該腳本為了透明與可攜性而以最大概似配適；brms 以馬可夫鏈蒙地卡羅在單一次呼叫中配適同一個模型，是應用工作建議的工具，尤其當尺度次模型自己還帶隨機斜率時。

16.7 常見的迷思

關於個體內變異性，有幾個信念會誤導讀者。第一個是一個估得夠好的個體內標準差就跟一個模型一樣好；並不是，因為即使時點數很多，它仍然沒有收縮、無法承載共變項，而且只要水準與變異性有耦合，它就把兩者混淆。第二個是尺度次模型是一種異質變異數的修正；在這裡它就是假設本身，而個體內變異數是感興趣的結果變項，不是要掃掉的干擾項。第三個是關於變異數的研究發現本質上比較難複製；真相是它們對設計的要求比較高，要更多的人與更多的時點，而它們表面上的脆弱 (fragility) 通常是一個檢定力問題，由第 4 章那套誠實的規劃來處理。第四個是五波的追蹤資料就能支撐個體內變異數的建模；不能，因為每人寥寥幾個時點根本估不出那個人的波動性，而這套方法屬於密集設計 (intensive design)。

接下來要去哪裡

本章刻意停在貝氏估計的門檻上。位置尺度模型在這裡為了透明而以最大概似配適，但它自然的引擎是馬可夫鏈蒙地卡羅，它能處理隨機尺度斜率、多層次的延伸，以及最大概似很吃力才給得出來的後驗傳遞。第 17 章完整提供這套做法，而想先拿到食譜再回頭補理論的讀者，現在就可以用 brms 配適本章的模型，之後再回到第 17 章求理解。個體內變異數會在本書密集縱貫的部分一再出現：第 25 章為一個動態模型的創新變異數 (innovation variance) 中的異質性建模，那是隨機尺度效果的動態版表親；第 33 章與第 35 章則把情緒不穩定部署成臨床與人格上的構念。「值得建模的不只是中心、還有分散程度」這一課，會出現在每一個被測量得夠密、密到足以顯示它擺盪多大的歷程上。

習題

1. 16.1 缺陷示範。模擬一批時點數不等、但真實波動性相同的人，並顯示依原始個體內標準差為人排序是不穩定的，而位置尺度模型收縮後的估計值則不然。
2. 16.2 配適尺度次模型。在一個日誌結果變項上，配適一個兩個次模型都放進同一個個體內預測變項的位置尺度模型，並寫出尺度次模型的結果段落，以乘性的變異數尺度報告效果。
3. 16.3 一個尺度隨機斜率。把模型延伸成允許預測變項對波動性的效果隨人而異，誠實地評估收斂，並報告資料是否支持這份額外的複雜度。
4. 16.4 瞄準變異數的檢查。提出並實作兩項瞄準變異數結構的後驗預測檢查，並解讀各自能揭露尺度次模型的哪一種設定錯誤。
5. 16.5 位置、尺度，還是兩者。針對三個研究問題，判定各自關切的是水準、變異性，還是兩者的關係，並設定能回答它的次模型。