連續時間模型

27.1 間隔的醜聞

推動整章的那個事實令人不安：離散時間的動態係數，離開它的間隔就沒有意義。設想一個連續的歷程(process)，壓力隨時間把負向情緒往上推。三支研究團隊研究它，一支每半個時間單位取樣一次，一支每兩個單位，一支每四個單位；每一支都配適一個交叉延宕(cross-lag)模型，並報告壓力對後來負向情緒的效果。圖 27.1 是他們的發現：交叉延宕效果分別是 \(0.10\)、\(0.19\) 與 \(0.14\)。三支團隊看起來彼此不合，比對這三篇論文的讀者可能會斷定這個效果很脆弱，或者這些研究互相牴觸。但三者都對，都是同一個底層系統在不同取樣率下的正確估計值。離散的交叉延宕不是這個歷程單獨的性質，而是歷程乘上間隔的性質；報告它卻不報告間隔，就像報告一個速度卻不給時間單位。

註：三份研究看的是同一個連續時間系統，取樣間隔分別為 \(0.5\)、\(2\) 與 \(4\) 個時間單位，報出來的交叉延宕效果卻差了將近兩倍。沒有一份是錯的；每一份都是在不同的間隔上讀同一個底層歷程。離散的交叉延宕離開它的間隔就沒有意義。

產生這種行為的關係就是本章的主方程式，先在數值上見一面，再談形式。底層的連續動態若由漂移矩陣(drift matrix) \(\mathbf{A}\) 支配，跨越間隔 \(\Delta t\) 的離散轉移矩陣(transition matrix)就是矩陣指數(matrix exponential) \(\boldsymbol{\Phi}(\Delta t) = e^{\mathbf{A}\Delta t}\)；轉移矩陣裝的正是離散模型在那個間距上會估出來的自我迴歸與交叉延宕。把單一個漂移矩陣餵進這個關係，跑過一整排間隔，就得到圖 27.2 的效果對間隔曲線。那是本章的核心圖，也是要內化下來的那一張。自效果(auto-effect)是對角線上的元素，隨著間隔變長而由 \(1\) 單調衰減(decay)到零，因為一個變項對它自己遙遠過去的依賴會淡去。交叉效果(cross-effect)是非對角線上的元素，行為更有意思：極短的間隔上它從接近零開始，影響還沒有時間累積(accumulate)；中等的間隔上升到高峰，這裡靠近 \(\Delta t = 2\)；再長就隨著整個系統被沖淡而下降。在極短間隔上取樣的研究幾乎找不到交叉效果，可能斷定它不存在；落在高峰上的研究會發現它很強；長間隔的研究則會發現它又變弱了。存在一個讓效果最容易被偵測到的最適間隔(optimal interval)，這一點由 Dormann and Griffin (2015) 發展出來；而一份離散研究報告的那個單一數字，只是這條曲線上的一個抽樣(draw)。帶有回饋(feedback)的系統裡，兩條交叉效果還會隨間隔改變大小順序，於是哪一個變項看起來領先竟取決於間距，這是 Kuiper and Ryan (2018) 分析過的一項危害。第 20 章為潛在變化分數提出的間隔告誡、第 21 章的跨研究可比較性問題、第 25 章離散化(discretization)帶來的粗化(coarsening)，全都是同一個底層真相的不同側面。

註：一個固定的漂移矩陣經由 \(\boldsymbol{\Phi}(\Delta t) = e^{\mathbf{A}\Delta t}\) 所隱含的離散自效果與交叉效果，畫成間隔 \(\Delta t\) 的函數。自效果單調衰減；交叉效果升到靠近 \(\Delta t = 2\) 的高峰之後下降。任何一份離散研究都只報出這些曲線上的一個點。這是本章的反射圖。

27.2 連續時間模型

這些曲線底下的模型是一個隨機微分方程(stochastic differential equation)，它的符號帶著撐得過形式化的心理學意義。以 \(\boldsymbol{\eta}(t)\) 表示連續時間 \(t\) 上的潛在狀態(latent state)向量，方程式 \(d\boldsymbol{\eta}(t) = \big(\mathbf{A}\,\boldsymbol{\eta}(t) + \mathbf{b}\big)\,dt + \mathbf{G}\,d\mathbf{W}(t)\) 有三個部分。漂移項 \(\mathbf{A}\,\boldsymbol{\eta}(t)\,dt\) 是決定性(deterministic)的那份拉力：任何一個瞬間，系統都被拉往漂移矩陣作用在當下狀態上所設定的方向；歷程若會回復到平衡點，這份拉力就是往平衡點回去。截距 \(\mathbf{b}\) 決定那個平衡點的位置。擴散項 \(\mathbf{G}\,d\mathbf{W}(t)\) 是持續不斷的隨機輸入，是創新(innovation)在連續時間中的對應物。其中 \(\mathbf{W}(t)\) 是一個 Wiener 歷程(Wiener process)，也就是連續時間中的隨機漫步(random walk)，它的增量獨立且高斯，是那些讓系統一直活著、永不停歇的小衝擊(shock)的來源。要用這個模型並不需要測度論(measure theory)：漂移是系統被拉往哪裡，擴散是它一路上被搖晃得多厲害。

漂移矩陣是要拿來解讀的那個物件，而它的對角線正是本章交得出來最可解讀的通貨。一個為負的對角元素 \(a_{ii}\) 就是一個回復速率：它越負，這個變項受擾動之後被拉回平衡點(equilibrium)的速度就越快。那個速率翻譯成一個半衰期(half-life)，也就是一次擾動衰減到一半大小所需要的時間，由 \(\ln 2 / |a_{ii}|\) 給出。這正是第 24 章對自我迴歸做過的那個詮釋動作，只是改在連續時間裡進行，而且不依賴任何取樣間隔。一個 \(-0.40\) 的漂移自效果對應到 \(1.73\) 個時間單位的半衰期，一個 \(-0.60\) 的則對應到 \(1.16\)。非對角的元素是連續時間的交叉效果，是一個變項對另一個變項的瞬時耦合，它們帶著與第 25 章離散交叉延宕相同的詮釋紀律：它們是條件的(conditional)、嵌在系統裡的量，不是孤立的因果(causal)箭號，而且沒有標準化(standardization)就不可以跨「量尺不同的變項」比較。表 27.1 把連續與離散兩套詞彙並排，表 27.2 則連同半衰期一起解讀那些漂移參數。

連續時間能夠涵蓋任何間距，理由在那個精確的離散解。跨越一個間隔 \(\Delta t\)，這個歷程精確地滿足 \(\boldsymbol{\eta}_{t} = e^{\mathbf{A}\Delta t}\,\boldsymbol{\eta}_{t-1} + \mathbf{w}_t\)；\(\mathbf{w}_t\) 是高斯的，它的共變數是擴散在這個間隔上的一個積分，用基礎概念方塊記下的 Van Loan 方法算出來。同一個漂移矩陣因此在每一個間隔上都生成正確的離散動態。於是一份間距混雜的資料，日間空隙短、跨夜空隙長，處理起來不必假裝那些空隙相等。測量那一層就是第 26 章的狀態空間觀察方程式：連續的狀態透過一個負荷量矩陣(loading matrix)被觀察到，並帶著測量誤差(measurement error)。配適這個模型的引擎，就是那一章的卡爾曼濾波器(Kalman filter)，只是轉移矩陣在每一個時點依它自己的間隔重算一次。Ornstein-Uhlenbeck 歷程是那個純量(scalar)的情形，也就是連續時間的一階自我迴歸(first-order autoregressive, AR(1))，是建立直覺最好的起點。圖 27.3 呈現它的樣本路徑：每一條都以自效果所設定的速率被拉回平衡點，永遠被擴散推開，也永遠被拉回來，而半衰期標示出一次偏離衰減得多快。它的定態(stationary)分配有一個變異數，基礎概念方塊會給出。整幅圖像把第 9 與第 20 章預告過的向量場(vector field)與吸引子(attractor)直覺形式化了。

註：一個 Ornstein-Uhlenbeck 歷程的四條樣本路徑，那就是連續時間的 AR(1)。每一條都以自效果所設定的速率被拉向平衡點（藍線），並被擴散搖晃；陰影帶是一個定態標準差。半衰期，也就是那支紅色箭號，是回復速率可解讀的摘要。

整個系統的性格由漂移矩陣的特徵值(eigenvalue)決定，那正是連續時間與第 32 章振盪子(oscillator)接上頭的地方。圖 27.4 呈現兩個系統。特徵值若是實數且為負，系統在受擾動之後平滑地衰減回平衡點，那就是工作範例裡回復到平衡點的行為。特徵值若是複數、實部為負而虛部不為零，系統就一邊回去一邊振盪(oscillate)，越過平衡點再以一種阻尼(damped)的節奏盪回來；那是第 24 章遇過的二階振盪在連續時間裡的面貌，也是通往第 32 章微分方程振盪子的橋。特徵值的實部(real part)支配系統安頓下來的速度，虛部(imaginary part)支配它振盪的速度。圖 27.5 把同樣的動態畫成一個向量場：每一支箭都是漂移矩陣在狀態空間中該一點上施加的瞬時拉力，全部朝平衡點流動，並被交叉效果扳彎。這就把本書自第 9 章起以預告方式呈現的那些吸引子圖像形式化了。

註：左：一個衰減系統（實數、為負）與一個振盪系統（複數）的特徵值在平面上的位置；陰影區標示穩定（實部為負）。右：兩者各自產生的軌跡，平滑的回復對照阻尼的振盪。特徵值的虛部是振盪的來源，也是通往第 32 章的橋。

註：每一支箭都是漂移在該狀態上施加的瞬時拉力 \(\mathbf{A}\boldsymbol{\eta}\)；全部朝平衡點（紅點）流動。那份不對稱，也就是流被扳彎的地方，正是壓力到情緒的交叉效果。這是第 9 與第 20 章那些吸引子預告的正式版本。

27.3 配適連續時間模型

這些模型的生產工具是 ctsem 套件。它設定一個連續時間模型，用頻率論(frequentist)的卡爾曼最大概似(maximum likelihood)配適，或用一個以 Stan 為底的階層貝氏(hierarchical Bayesian)模式配適，再回傳帶著不確定性的漂移矩陣。本章的分析自己寫了一具連續時間引擎：轉移用矩陣指數，歷程共變數用 Van Loan 積分，外加一個隨時間變動版本的第 26 章卡爾曼濾波器，在每一個時點依它自己的間隔把兩者重算一次。這樣做是為了讓整套做法能對照一個已知的真相逐項稽核；應用工作上該去拿的工具仍然是 ctsem，表 27.3 記下它的各種模式。工作階梯的第一階是還原。模擬一個帶著已知漂移矩陣的雙變項連續時間歷程，在不規則的間隔上取樣，再以最大概似配適連續時間模型，漂移就還原了出來：自效果回來成 \(-0.45\) 與 \(-0.56\)，生成用的是 \(-0.40\) 與 \(-0.60\)；壓力到情緒的交叉效果回來成 \(0.29\)，生成用的是 \(0.25\)；半衰期回來成 \(1.55\) 與 \(1.25\) 個單位，真值則是 \(1.73\) 與 \(1.16\)，以單一條序列而言算是接近了。這個估計式從不等間隔的資料中還原出與間隔無關的動態，那正是整件事的重點。

# 延宕翻譯器與一次連續時間配適（手寫引擎）
library(Matrix)
A <- matrix(c(-0.40, 0.00,          # 漂移：為負的對角線就是回復速率
               0.25, -0.60), 2, 2, byrow = TRUE)
Phi <- function(dt) as.matrix(expm(A * dt))       # 間隔 dt 上的離散效果
Phi(0.5); Phi(2); Phi(4)            # 效果對間隔曲線，一次一個 dt
log(2) / abs(diag(A))               # 半衰期：1.73 與 1.16 個時間單位
# 以轉移的最大概似對不等間隔資料配適一個 CT-VAR（分析腳本中的 ct_fit_bivar）：
# 先為每一個相異的間隔預先算好 Phi(dt) 與 Van Loan 的 Q*(dt)，
# 再在 A 與擴散上最大化高斯的轉移概似。
fit <- ct_fit_bivar(Y = as.matrix(ct_sim[, c("stress","na")]), dtv = ct_sim$dt)
fit$A; fit$half_life                # 還原出來的漂移與半衰期

第二階是人。對每一個個體的序列配適一個連續時間 AR(1)，就還原出一個逐人的回復速率，因而也還原出一個逐人的半衰期；而這些半衰期的分配(distribution)，是一項以動態為通貨的個別差異(individual differences)發現。圖 27.6 同時呈現真實資料的分配與它的驗證。在 affect_ema 的 \(120\) 人中選了 \(80\) 人（執行時取的子集），用每個人實際的提示(beep)時戳(timestamp)，也就是日間兩個半小時、跨夜十一個半小時的間隔，負向情緒的半衰期散得很開：中位數(median)靠近三小時，第十到第九十百分位數(percentile)的範圍由不到一小時到將近八小時。人與人在「情緒回到基線(baseline)有多快」上差別很大，這是一項逐人(person-specific)的參數。這樣的異質性(heterogeneity)還原得出來，由模擬確認：生成一批回復速率已知而各不相同的人，再逐一重新估計，還原出來的半衰期與真值相關 \(0.63\)，在一份日記研究的序列長度下並不完美，但顯然帶著訊息。在 ctsem 的階層貝氏模式裡，這些逐人參數會連同它們的母體(population)分配一起估計，那正是第 25 章隨機效果那套跨人借力；而隨機自我迴歸在這裡再度以「來自一個漂移矩陣分配的抽樣」的身分登場。

註：左：對 affect_ema 作連續時間 AR(1) 配適得到的逐人負向情緒半衰期分配，用的是每個人真實而不規則的時戳；中位數靠近三小時。右：一項模擬驗證，估出來的半衰期對照跨人的真實半衰期，以 \(0.63\) 的相關還原出這份異質性。

第三階是刻意的示範：連續時間不是只給密集資料用的，它也修補追蹤研究。設想一份只有少數幾波(wave)的追蹤(panel)研究，那幾波在不同的人身上落在真正不同的實際時間上，有些人隔一個月測、有些人隔三個月測。這正是離散模型會處理不好的情形，因為它把同一個延宕係數套在並不相等的間隔上。圖 27.7 呈現時點的柵狀圖(raster)，每個人的施測各自落在自己的時間上，後果就放在旁邊。假裝間距相等的離散模型只報出單一個 \(0.14\) 的交叉延宕；真實的交叉效果卻隨著資料裡實際存在的那些間隔而變動，由最短間距上的大約 \(0.10\) 到最長間距上的 \(0.18\)。那份假裝把這一整段範圍壓成單一個數字，而這一次模擬報出的那個數字對應到大約 \(0.84\) 個時間單位的間距。\(0.84\) 只是一次抽樣的落點：同一份 \(200\) 人、\(5\) 波的模擬重複 \(300\) 次，合併估計值的平均是 \(0.1601\)，圖上這一次的 \(0.14\) 落在第 \(22\) 百分位數。系統性的答案要看估計標的：定態之下，這個假裝等距的估計式的機率極限是轉移矩陣在間隔分配上的期望值，不是平均間隔上的轉移矩陣，而交叉效果在這份設計的間隔範圍上是凹的，由 Jensen 不等式得 \(0.1576\) 對照 \(0.1753\)；把它放回效果對間隔曲線解出來的隱含間距是 \(1.08\)，短於這份設計裡實際的平均間隔 \(1.42\)。連續時間模型還原出把上述全部正確生成出來的那個漂移矩陣。對照第 21 章的隨機截距(random intercept)交叉延宕追蹤模型，讀法是：連續時間加上去的，恰恰是跨人不等間距所要求的東西。所以任何一次連續時間配適，應用上的產出從來都不是那個原始的漂移本身，而是它的翻譯：讀者關心的那些間隔上隱含的(implied)離散矩陣、帶著不確定區間(interval)的半衰期，以及效果對間隔曲線，好讓這份結果與其他間距的研究比較得起來。

註：左：一份五波追蹤研究中逐人施測時點的柵狀圖；每個人的施測都落在不同的實際時間上。右：一個假裝間距相等的離散模型報出單一個交叉延宕（\(0.14\)），而真實的交叉效果在資料實際存在的那些間隔上橫跨 \(0.10\) 到 \(0.18\)。連續時間如實對待那份不相等的間距。

27.4 限制、誠實與選擇

連續時間很有力，但不是萬能，而它最重要的限制正是第 1 章以疊頻(aliasing)之名預告過的那一項，這裡改用矩陣指數這具透鏡重新看一次。歷程若運作得比取樣速度快得多，就還原不出來：兩次觀察之間，那些快的動態早已跑完，觀察值看起來幾乎彼此獨立，不帶任何關於回復速率的訊息。圖 27.8 把這件事量化。在一整排取樣間隔上估計一個快速歷程的回復速率，只要取樣間隔短於歷程的半衰期，半衰期就還原得不錯；間隔一旦超過半衰期，不確定性就炸開；等到間隔是半衰期的好幾倍，點估計本身也垮掉，漂移實際上就識別不出來了。這兩個位置量化得出來：不確定性在比值 \(3.9\) 附近炸開，點估計的中心則要到比值 \(7.8\) 才失效，中間那一段最危險，因為點估計看起來還好好的。圖上比值 \(5.19\) 那一格印的是一次抽樣；在同一個比值上重複估計，中位數只偏離真值 \(8\%\)，垮掉的是離散程度，不是中心。實務上的讀法是：取樣率必須配上你希望研究的那個歷程的時間尺度，而那是一項事後任何分析都補救不了的設計決定。另外兩項限制界定了適用範圍。這裡呈現的模型是線性的；真正帶著非線性(nonlinear)動態、閾值(threshold)或依狀態而定之回復的歷程，要用第 32 章的非線性擴充。這裡也假設歷程在一個固定的平衡點附近是定態的；讓平衡點本身移動的連續時間成長模型是存在的，本章也引用了，但這套主力做法假設的是一個穩定的吸引子。

註：一個快速歷程估出來的半衰期，對照「取樣間隔比上真實半衰期」的比值。比值小於一（間隔在半衰期之下）時還原良好；越過之後不確定性，也就是那條陰影帶，先炸開，點估計要到比值數倍時才垮掉。圖上每一格印的都是一次抽樣；重複估計時，點估計的中心要到比值 \(7.8\) 附近才失效。取樣率必須配上歷程的時間尺度。

要不要付出連續時間那份概念成本，這個決定跟著資料的間隔結構走，圖 27.9 把它攤了出來。人之內或人之間的間隔變動很大時，連續時間值得，因為替代方案會誤述動態。間隔相等或幾乎相等時，一個離散模型就夠用，也比較簡單；而第 25 章動態結構方程模型的 TINTERVAL 離散化是一條中間路線，它把中度不等的間隔對齊到一格細網格上，不必動用完整的連續時間那一套做法。等距的假裝只有在間隔變異數小的時候才勉強可以接受，而本章的模擬把「什麼時候不可以」量化了出來：以追蹤範例那種很寬的間距變動來說，這份假裝顯著地錯過了那些依間隔而異的效果。更大的願景是可比較性，那也是這份額外成本值得學的理由。漂移矩陣與間隔無關，所以它是比較設計不同的研究、以及對動態作後設分析(meta-analysis)的天然共同通貨；後者是一個活躍的研究方向，目標是像一般後設分析彙整效果量(effect size)那樣，跨研究彙整漂移矩陣。表 27.4 以表格形式給出這個決定，表 27.5 則給出報告檢核表，而它不可讓步的那一項就是隱含間隔的翻譯。

註：間隔相等或幾乎相等時，一個離散模型是可以的；中度不等的間隔可以用第 25 章的 TINTERVAL 離散化處理；顯著不等的間隔，或者跨研究比較的需要，就要用連續時間的漂移矩陣。代價是概念上的額外成本，回報是可比較性。

寫一份連續時間分析，有一項要求不可讓步，那也正是它有別於任何離散報告的地方：結果必須先以與間隔無關的形式給出，再翻譯成讀者用得上的那些間隔。報告要呈現帶著不確定性的漂移矩陣，把自效果轉成帶區間的半衰期，並在給出交叉效果時附上告誡，說明它們是連續的耦合，不是離散的交叉延宕。接著要翻譯：一到兩個具實質意義的間隔上隱含的離散轉移矩陣，最好還有效果對間隔曲線本身，好讓這項發現與其他間距的研究比較得起來。報告也要說明被建模的觀察時點(observation time)結構，註明用的是頻率論還是階層模式以及它的先驗，並誠實面對取樣率的限制，承認這份設計解析得了與解析不了哪些時間尺度(timescale)。示範的段落是這樣的：「本研究對不規則間隔的資料配適了一個雙變項連續時間模型，漂移矩陣以最大概似估計。負向情緒以 \(1.25\) 個時間單位的半衰期回復到平衡點（在這份日記資料中約為三小時），壓力的半衰期則為 \(1.55\)；壓力對負向情緒有一個為正的連續時間效果，它隱含的離散交叉延宕在接近兩個單位的間隔上達到高峰。因為這個效果依賴間隔，我們一併報告在一個單位與四個單位間距上隱含的離散矩陣，好讓這些估計值能與使用那些延宕的研究比較。」

習題

1. 27.1 延宕翻譯器。用表 27.2 的漂移矩陣算出 \(\Delta t\) 為 \(0.5\)、\(2\) 與 \(4\) 上的轉移矩陣，重現圖 27.1 那三個看似互相牴觸的交叉延宕；再掃過一整排間隔畫出圖 27.2 的效果對間隔曲線，指出交叉效果的高峰落在哪個間距。
2. 27.2 還原與翻譯。用 ctsem 對 ct_sim 的不等間隔資料配適雙變項連續時間模型，把還原出來的自效果與半衰期對照生成的真值，再依表 27.5 寫出結果段落，一併交出一個單位與四個單位間距上隱含的離散矩陣。
3. 27.3 逐人半衰期。用每個人真實的提示時戳，對 affect_ema 逐人配適連續時間 AR(1)，畫出負向情緒半衰期的分配並與圖 27.6 對照；說明只取 \(120\) 人中的 \(80\) 人這件事，會怎麼限制你對這份異質性所能作的陳述。
4. 27.4 等距的假裝。對圖 27.7 那份波次間距不等的追蹤資料，分別配適一個假裝間隔相等的離散模型與一個連續時間模型，把兩者的估計標的對照清楚，並解釋那個單一的 \(0.14\) 為什麼只對應到 \(0.84\) 個時間單位，遠短於設計裡 \(1.42\) 的平均間隔。
5. 27.5 手算與說明。以 \(a = -0.40\) 手算間隔為 \(2\) 時隱含的離散自我迴歸與 \(\ln 2 / |a|\) 的半衰期，再寫一份給非技術審查人的備忘錄，說明一個 \(0.25\) 的漂移交叉效果為什麼在任何間距上都不是「一個 \(0.25\) 的交叉延宕」。
6. 27.6 取樣率的誠實。重現圖 27.8 的疊頻模擬，指出不確定性在哪個比值上先炸開、點估計又要到哪裡才垮掉；接著依表 27.4 判斷自己手上資料的間隔結構，並寫出第 27.4 節所要求的那段取樣率告白。