動態因素分析與狀態空間模型

26.1 一個形式，容納所有模型

狀態空間形式就是兩條方程式。觀察方程式 \(\mathbf{y}_t = \mathbf{Z}\,\boldsymbol{\alpha}_t + \boldsymbol{\varepsilon}_t\)，其中 \(\boldsymbol{\varepsilon}_t \sim N(\mathbf{0}, \mathbf{H})\)，說的是在時點 (occasion) \(t\) 測到的東西，是一個負荷量矩陣 (loading matrix) \(\mathbf{Z}\) 乘上一個潛在狀態 \(\boldsymbol{\alpha}_t\)，再加上測量誤差 \(\boldsymbol{\varepsilon}_t\)。狀態方程式 \(\boldsymbol{\alpha}_t = \mathbf{T}\,\boldsymbol{\alpha}_{t-1} + \boldsymbol{\eta}_t\)，其中 \(\boldsymbol{\eta}_t \sim N(\mathbf{0}, \mathbf{Q})\)，說的是潛在狀態依一個轉移矩陣 (transition matrix) \(\mathbf{T}\) 演化，並由歷程雜訊 \(\boldsymbol{\eta}_t\) 推動。這兩個雜訊項不是同一回事，而把它們分開正是整件事的重點。測量誤差 \(\boldsymbol{\varepsilon}_t\) 是施測的雜訊，是工具的不可靠，是一個情緒評分裡「反映的是題目本身含糊、而不是這個人真實狀態」的那一部分。歷程雜訊 \(\boldsymbol{\eta}_t\) 則是真正的衝擊，是這個人的狀態從一個時點到下一個時點的真實變化。前面每一章只要把一個觀察分數迴歸到它自己的延宕 (lag) 上，就默默地把這兩者熔在一起，而本章要把它們分開。圖 26.1 把這兩條方程式連同它們的心理學讀法一起畫了出來。

註：狀態方程式讓一個潛在狀態 \(\alpha_t\) 透過轉移矩陣 \(\mathbf{T}\) 演化，並受到歷程雜訊（真實的變化）擾動。觀察方程式透過負荷量矩陣 \(\mathbf{Z}\) 去測量這個狀態，而測量會被測量誤差（施測的雜訊）汙染。把這兩種雜訊分開，是這個形式在心理計量上最核心的貢獻。這張圖在第 27 與第 32 章會再次出現。

這項分離立刻帶來一個把第 24 章重新框起來的事實。一個被帶著測量誤差觀察到的一階自我迴歸 (first-order autoregressive, AR(1))，根本就不是一個自我迴歸；代數上它是一個 ARMA(1,1)，而它天真的延宕一期自我相關 (autocorrelation) 會依這個測量不可靠的程度低估真實的延續效果 (carryover)。機制正是測量誤差加在任何一個相關上的那份衰減 (attenuation)，只是這裡被用在「一條序列與它自己的過去」的相關上。圖 26.2 在 kalman_n1 上把這件事講得很具體，那是一條由潛在的 AR(1) 生成、真實延續效果為 \(0.60\)、在信度 \(0.65\) 下被帶誤差觀察到的單一受試者 (single-subject) 序列。左圖顯示天真的自我迴歸隨著信度變差而穩定下滑，掉到遠低於真值的位置；在這條序列上，天真的延宕一期自我相關是 \(0.46\)，離生成用的 \(0.60\) 差得很遠。右圖則把修補與它的代價一起呈現：一個把測量誤差另外建模的卡爾曼最大概似 (maximum likelihood) 配適還原出 \(0.61\) 的延續效果，而跨一百條模擬序列來看它平均而言是不偏的，但它在每一條序列上的估計值都很有雜訊，因為要從一條短序列裡把兩個變異數分開，是一項苛刻的要求。這項分離在期望上是成功的，而且必須連同它的不確定性 (uncertainty) 一起報告。

註：左：一條帶誤差的序列，它天真的自我迴歸隨著信度變差而掉到真實延續效果（\(0.60\)）之下。右：跨 100 條單一序列的抽樣分配 (sampling distribution)；把測量誤差另外建模的卡爾曼估計值落在真值附近，天真的估計值則被往零的方向衰減。這項分離平均而言是不偏的，但由一條序列來看很不精確。

同樣的兩條方程式，換上不同的矩陣，就成了那些熟悉的模型，而把那些矩陣填好正是本節要建立的本事。圖 26.3 列出五張卡片。一個帶測量誤差的 AR(1) 有一個純量 (scalar) 狀態，\(\mathbf{T} = \phi\)、\(\mathbf{Z} = 1\)、歷程變異數 \(\mathbf{Q} = \sigma^2_\eta\)、測量變異數 \(\mathbf{H} = \sigma^2_\varepsilon\)。一個局部線性趨勢，也就是「一個水準 (level) 隨著一個緩慢改變的斜率 (slope) 漂移」的結構時間序列 (structural time series) 模型，有一個由水準與斜率組成的二維狀態，而它的轉移矩陣把斜率帶進水準裡。一個 VAR(1) 有一個向量狀態與一個完整的轉移矩陣，也就是第 25 章的模型，在無誤差測量時 \(\mathbf{H} = \mathbf{0}\)，否則 \(\mathbf{H} \neq \mathbf{0}\)。一個動態因素模型有一個會演化的潛在因素狀態，以及一個透過數個指標 (indicator) 去測量它的負荷量矩陣。而第 19 章的潛在成長曲線，就是一個狀態為截距 (intercept) 與斜率、轉移矩陣讓它們保持不變、負荷量矩陣隨時間變動而把斜率讀到接續的各個時點上的狀態空間模型，這一點把第五與第六部分統一在同一個形式之下。表 26.1 為每個矩陣命名，並給出它的維度、心理學意義，以及它在結構方程或多層次 (multilevel) 模型中的對應物；表 26.2 則說明五個模型各自在每個矩陣裡放了什麼。

註：每一個熟悉的模型都是「填了特定矩陣的狀態空間形式」。自我迴歸、趨勢模型、向量自我迴歸、動態因素模型與潛在成長曲線之間的差別，只在於狀態裡裝什麼、轉移矩陣怎麼移動它，以及負荷量矩陣怎麼讀它。潛在成長的負荷量隨時間變動，而那正是一個靜態的成長模型成為狀態空間模型的關鍵。

26.2 卡爾曼這具引擎

卡爾曼濾波器是一個沿著序列往前走的遞迴，而它在每一個時點做兩件事。它預測 (predict)：拿前一個時點對狀態最好的估計，用轉移矩陣把它帶到對此刻狀態的一個預測上，而不確定性依歷程雜訊增長。它更新 (update)：新的觀察值抵達，它與預測所隱含的值之間的差距就是創新，而狀態估計值會依一個叫作卡爾曼增益的量往觀察值修正。圖 26.4 把這個迴路畫了出來。增益是整件事的核心，而它就是第 13 章收縮意義下的信度權重。當測量相對於狀態的不確定性而言很精確時，增益很大，估計值幾乎整個移向觀察值；當測量有雜訊時，增益很小，估計值就留在模型的預測附近。濾波器是在以最佳的方式權衡一個事前的預報與一個帶雜訊的資料點，以精確度對精確度，而那正是「把一個人的估計值往群體平均數收縮」的同一套加權平均邏輯，只是現在在時間上運作。

註：在每一個時點，濾波器先由轉移矩陣預測狀態（不確定性依歷程雜訊增長），再依卡爾曼增益往新的觀察值更新；增益是一個信度權重，測量精確時它大，測量有雜訊時它小。概似、遺漏資料的處理與標準化殘差都由同一趟掃描掉出來。

從那一趟前向掃描裡掉出三樣東西，而它們合起來就是這套做法值得付出的理由。概似透過預測誤差分解 (prediction-error decomposition) 免費取得：每一個創新以它自己的變異數標準化之後貢獻一個高斯項，而它們的和就是精確的對數概似 (log-likelihood)，於是轉移、負荷量與各變異數參數的最大概似估計，就是一連串濾波器執行結果上的最佳化。遺漏資料被精確處理而且不需要插補 (imputation)：當一個時點沒有被觀察到，更新這一步就直接跳過，預測往前帶，而狀態的不確定性持續增長直到下一個觀察值抵達，那正是第 24 章指出「天真的內插 (interpolation) 會灌大自我相關」的那個問題有原則的修補。而標準化的一步前 (one-step-ahead) 創新就是診斷用的殘差，模型對的時候它們是白的、常態的。基礎概念方塊以形式化的方式寫出概似分解與「AR 加誤差」的代數。

kalman_n1 上的卡爾曼弧線把這項分離示範了出來。以最大化濾波器概似的方式配適「AR 加誤差」模型，還原出的延續效果是 \(0.61\)，對照生成用的 \(0.60\)；平滑器接著在每一個時點產生對潛在狀態最好的回溯 (retrospective) 估計。圖 26.5 在有雜訊的觀察值底下畫出濾波與平滑的潛在路徑，連同它們的不確定帶。濾波估計值只用過去，因而可以即時 (real time) 取得；平滑估計值用了整條序列，因而更緊、也更貼近真實的潛在狀態，是回溯的最佳猜測。在兩者之間的選擇是一項估計標的 (estimand) 的選擇，表 26.3 把它列了出來：濾波是給「資料一邊進來一邊監看」用的，也就是第 33 章在及時介入 (just-in-time intervention) 上會發展的臨床用途；平滑則是給一份已完成的研究所要的歷史重建用的。

註：kalman_n1 上潛在狀態的濾波（只用過去）與平滑（整條序列）估計值，以及平滑的不確定帶，對照有雜訊的觀察值與真實的潛在路徑。兩者都把施測雜訊底下的訊號撈了出來；平滑更緊，也更貼著真值走。

遺漏資料的修補正是第 24 章承諾過的回報。把 kalman_n1 的一段遮起來再跑濾波器，缺口中的狀態估計值不是一條直線，而是一條往歷程平均數回復的曲線，它的不確定帶隨著資料的缺席而變寬，並在觀察值恢復時再度收攏。圖 26.6 把這件事拿來與線性內插對照，後者在缺口上畫出一條自信的直線，因而宣稱了一份它並不擁有的確定性。後果是量得出來的：卡爾曼平滑器在缺口中重建的潛在軌跡對照真實潛在狀態的均方根誤差 (root-mean-square error) 是 \(1.43\)，內插則差得多，來到 \(2.04\)；而在有缺口的序列上以卡爾曼法估出來的延續效果是 \(0.57\)，仍貼近完整資料的值，反觀先內插再天真估計則會灌大被觀察到的自我相關。濾波器把一個遺漏的時點當成它本來的樣子，也就是訊息的缺席，而不是一個要去發明出來的數值。下面的程式面板把整條弧線收在一起。

註：在 kalman_n1 一段被遮起來的區間上，卡爾曼平滑器的不確定帶變寬、估計值往平均數回復，並跟著真實潛在狀態走（均方根誤差 \(1.43\)）；線性內插畫出一條自信的直線，卻沒有跟上（均方根誤差 \(2.04\)）。濾波器把遺漏表達成訊息的缺席，而不是一個被發明出來的數值。

# kalman_n1 上的卡爾曼弧線：把測量誤差與歷程雜訊分開
d <- readRDS("Examples/data/kalman_n1.rds")            # y（觀察值）、eta_true（潛在）
yc <- d$y - mean(d$y)
# 在 (phi, 歷程變異數 q, 測量變異數 r) 上最大化濾波器概似
negll <- function(par) {                               # kalman_filter() 手寫於 ch26_analysis
  phi <- tanh(par[1]); q <- exp(par[2]); r <- exp(par[3])
  -kalman_filter(matrix(yc), matrix(phi), matrix(1), matrix(q), matrix(r),
                 a1 = 0, P1 = matrix(q/(1-phi^2)))$loglik
}
fit <- optim(c(atanh(.3), 0, 0), negll, method = "BFGS")
tanh(fit$par[1])                                       # phi 約 0.61（天真的延宕一期自我相關 = 0.46）
# 平滑後的潛在狀態 + 跨越被遮起來的缺口時變寬的不確定性：
kf <- kalman_filter(matrix(d$y_gap - mean(d$y_gap, na.rm=TRUE)), ...)  # NA 的時點會被跳過
ks <- kalman_smoother(kf, matrix(tanh(fit$par[1])))

26.3 動態因素模型

動態因素模型就是「以一個潛在因素當狀態」的狀態空間形式，而它正是第 24 章留下缺陷的那個 P 技術的修補。P 技術把一個人的題目跨時點作因素分析 (factor analysis)，卻把那些時點當成彼此獨立的，於是丟掉了序列依賴，而序列依賴正是一條時間序列的重點。動態因素模型保留了因素結構，也就是少數幾個負荷到觀察題目上的潛在因素，再加上動態，讓那些因素依一個向量自我迴歸演化，於是這個潛在歷程既透過負荷量矩陣被測量，也透過轉移矩陣把自己往前帶。結果就是一個帶著逐人動態的逐人測量模型，也就是個殊 (idiographic) 測量的落實。存在兩種標準形式，而這個區分值得說清楚，圖 26.7 把它畫了出來。在歷程因素形式裡，也就是 Molenaar 的傳統，潛在因素遵循那個自我迴歸，而測量是同時 (contemporaneous) 且乾淨的，於是動態完全住在因素裡。在直接自我迴歸形式裡，自我迴歸被寫在因素分數 (factor score) 上，配上延宕的負荷量，於是記憶被分散到測量之中。這兩種形式在某些條件下可以做到等價，而表 26.4 說明何時偏好哪一種、以及各自的識別要求。

註：歷程因素形式（左）：潛在因素遵循自我迴歸，而負荷量乾淨地測量當下的因素，於是全部的動態都住在因素裡。直接自我迴歸形式（右）：記憶由「從較早的因素分數延宕負荷到當下指標」承載。兩種形式在某些條件下重合；歷程因素形式在心理學中比較常見。

這個模型還原得出一個已知的結構，是在 dfm_sim 上示範的，那是一條帶著已知負荷量矩陣與已知因素轉移的雙因素序列，這裡以一個手寫的期望最大化 (expectation-maximization, EM) 演算法配適，它的期望步驟就是卡爾曼平滑器，而它的最大化步驟以封閉解更新負荷量、轉移與各變異數。圖 26.8 呈現這份還原：負荷量回到接近它們的生成值，因素轉移還原出自我迴歸的對角線與由第二個因素進入第一個因素的那條交叉延宕，而平滑後的因素分數與真實因素的相關是 \(0.94\) 與 \(0.91\)。這個估計式從有雜訊的多變量指標中重建出一個潛在的動態結構，而那正是「在真相未知的資料上使用它」所需要的驗證。

註：以卡爾曼 EM 配適到 dfm_sim 的雙因素動態因素模型，其估計出來的負荷量對照生成值（左），以及因素轉移矩陣（右，估計值與括號內的真值）。負荷量與因素 VAR 都被還原了出來，而平滑後的因素分數與真實因素的相關是 \(0.94\) 與 \(0.91\)。

個殊的承諾在真實題目上兌現。把一個單因素的動態因素模型配適到兩位 affect_ema 受試者的四個負向情緒題目上，讓因素自由地遵循它自己的自我迴歸，得到的是兩幅真正不同的圖像，如圖 26.9 所示。對第一個人而言，負向情緒這個因素既有一貫性也有持續性 (persistence)：四個題目都與它有實質的關聯，其中「覺得緊張」最強，而這個因素帶著 \(0.41\) 的中等慣性 (inertia)，於是他的負向情緒會賴著不走。對第二個人而言，這個因素幾乎沒有延續效果，慣性只有 \(0.06\)，而且指標的剖面也不一樣：「覺得緊繃」是它最強的表徵，「覺得緊張」卻幾乎與它無關。同樣的四個題目在不同人身上組裝成不同的個人結構、帶著不同的動態，那正是第 24 章的非遍歷性 (nonergodicity) 在測量這一層被具體化，也正是「一個統一的個體間因素分析，硬把同一套負荷量型態與零動態加在每個人身上」所表達不出來的東西。多人的策略沿著兩條路走：複製式單一受試者取向為每個人各自配適，再摘述他們結構與動態的分配；而目前最先進的多層次動態因素模型則把他們彙整起來，那就是第 25 章 DSEM 的測量那一層，同一個想法從因素分析這一側抵達。

註：兩位 affect_ema 受試者四個負向情緒題目的單因素動態因素模型，呈現每個題目與這個人的因素之間的相關。第 3 人的因素有一貫性也有持續性（慣性 \(0.41\)）；第 17 人的慣性接近零（\(0.06\)），指標的剖面也不同。逐人的測量與逐人的動態，正是一個統一的個體間模型捕捉不到的東西。

26.4 診斷與選擇

一個配適好的狀態空間模型，跟任何時間序列模型一樣，是靠「它是不是把白雜訊 (white noise) 留在身後」來評判的，而受檢視的對象是標準化的一步前創新，也就是濾波器的殘差除以它自己預測出來的標準差。如果模型設定正確，這些標準化創新就是獨立的標準常態，於是有三項檢查適用。它們的自我相關應該可以忽略，這由眼睛與一個 Ljung-Box 檢定共同確認，因為殘留的序列結構意味著動態被設定錯了。它們的分配應該是常態的，這由一張分位數圖 (quantile plot) 與一個常態性檢定來查，因為厚尾 (heavy tail) 是離群值或非高斯歷程的訊號。而個別的大值標示出離群值，那在一條心理學序列中往往不是要被修剪掉的麻煩，而是要被理解而不是丟棄的實質事件，一個壓力源、一次危機、一次體制 (regime) 的轉換，是通往第 32 章的指路標。圖 26.10 呈現 kalman_n1 上「AR 加誤差」配適的診斷子圖組：標準化創新在它們的帶附近散布，它們的自我相關可以忽略，Ljung-Box 的 \(p\) 為 \(0.93\)，而它們的分位數圖是直的，常態性的 \(p\) 為 \(0.23\)，於是這個模型把序列漂白了，予以保留。

註：kalman_n1 上「AR 加誤差」配適的標準化一步前創新：它們的序列與 \(\pm 2\) 的帶（左）、它們的自我相關與 Ljung-Box 結果（中），以及它們的常態分位數圖（右）。白的、常態的、沒有離群值的殘差，表示模型是適當的。

在各種設定之間作選擇，遵循的是本書一路主張的精簡原則 (parsimony) 教條，並被心理學序列長度的現實磨得更銳利。巢套 (nested) 的模型，例如「AR 加誤差」對上一個單純的 AR，或一個單因素對上一個雙因素的動態模型，以預測誤差概似作概似比 (likelihood-ratio) 檢定來比較；非巢套的則用一個訊息準則 (information criterion)。至於「把狀態弄豐富」的誘惑，也就是加上第二個因素、更高的延宕階數、一個隨時間變動的參數，在心理學所蒐集到的序列長度上必須被抵抗，因為一個參數多過資料所能識別的狀態空間模型會收斂到一個邊界、報出一個被釘在零的變異數成分 (variance component)，或者配適了這個樣本卻複製不出來，那正是第 13 章為變異數成分提出的 Heywood 案例與過度參數化 (over-parameterization) 的同一組危害。陷阱方塊把這些模型誤導人的具體方式收在一起。

26.5 軟體，以及這套做法何時划算

狀態空間形式由數個各有所長的 R 套件支援，圖 26.11 與表 26.5 把這片地景攤了出來。KFAS 套件對線性高斯模型與它們的指數族 (exponential family) 擴充既快又精確，是本章濾波與平滑的主力。MARSS 套件專攻多變量自我迴歸狀態空間模型，介面對生態學與心理學的時間序列都很好上手。dynr 套件，它「調配模型」的慣例會在第 32 章再度出現，加上了體制轉換 (regime switching) 與非線性 (nonlinear) 動態，為接下來的變化點 (change point) 與振盪子 (oscillator) 模型作了準備。而 OpenMx 在一個結構方程的環境裡以定義變項 (definition variable) 界定狀態空間模型，是熟悉 SEM 的讀者最自然的家。Mplus 則透過它的動態結構方程模型成為鄰近的生產工具，而它獨到的加法是跨人的多層次彙整；值得記住的區分是：DSEM 加的是多層次那一層，狀態空間套件加的是結構上的自由度，也就是「設定一個固定的 DSEM 模板不會開放給你的自訂轉移與測量結構」的能力。本章的分析一個套件也沒有用，全部由手寫的濾波器、平滑器與 EM 演算法建成，好讓每一個物件都能對照一個已知的真相被稽核；有生產需求的讀者應該去拿那些套件，而想弄懂那些套件在做什麼的讀者，會在隨附的腳本裡找到機制。

註：主要幾個狀態空間工具的能力矩陣。強項各不相同：KFAS 精確而快，dynr 擅長體制轉換與非線性，OpenMx 擅長 SEM 式的界定，而 Mplus DSEM 擅長多層次彙整。選擇跟著分析所需要的東西走。

狀態空間這套做法什麼時候賺得回它的額外成本，是那個實務問題，而答案很具體。當測量誤差必須與歷程雜訊分開時它划算，因為沒有任何觀察分數模型做得到這件事；當資料有遺漏、而遺漏必須被精確處理而不是被插補時它划算；當目標是一個潛在的多變量歷程、一個沒有任何單指標模型構得到的動態因素結構時它划算；當結構是自訂的、需要一個套裝模型沒有提供的轉移或測量矩陣時它也划算。它不划算的時候，是第 25 章的多層次 VAR 用更少的組裝就回答得了這個問題、而測量誤差可以忽略的時候，那時比較簡單的工具才是誠實的選擇。下一章的連續時間 (continuous-time) 模型就是同一套狀態空間架構，只是轉移改由一個微分方程 (differential equation) 而非一個離散的步伐支配，那是這個形式容納第 25 章標出的不等間隔的方式；而第 32 章的體制轉換模型讓轉移矩陣本身會改變，那是這個形式容納本章離群值所暗示的變化點的方式。

寫出一份狀態空間或動態因素分析，是把前幾章的報告紀律延伸到這個形式特有的決定上。報告應該說明配適的是哪一個狀態空間模型、狀態裡裝了什麼、轉移與測量矩陣是怎麼設定的，尤其要說明測量誤差是被建模了還是被固定為零，因為單單這一項選擇就決定了動態會不會被衰減。它應該給出估出來的動態、它們的心理學讀法與它們的不確定性，報告遺漏資料是怎麼處理的，並呈現收斂與殘差的診斷，也就是那些認證這份配適的標準化創新檢查。對一個動態因素模型而言，它應該說明識別，也就是錨定或尺度的限制與延宕結構；而對一份個殊分析而言，它應該明說結構與動態是逐人的，並不假設可以推論到別人身上。一段示範的段落是這樣的：「本研究以卡爾曼 EM 對四個負向情緒題目配適了一個單因素的動態因素模型，因素為潛在的 AR(1)，測量誤差變異數自由估計；因素的慣性在第 3 人為 \(0.41\)、在第 17 人為 \(0.06\)，兩人的負荷量剖面也不同，顯示測量與動態都是逐人的；標準化創新是白的、近似常態的。」