文/中央研究院統計科學研究所博士後研究員 周芷妤
矩陣型態資料在當今生物醫學研究中相當常見。以基因與環境交互作用的偵測為例,傳統統計分析方法通常先將矩陣資料向量化,再建立統計模型 (例如廣義線性模型;Generalized Linear Model, GLM) 進行推論。然而,這類策略往往面臨高維度的問題,即參數數量龐大而樣本數相對不足,因此,難以獲得穩定且可靠的分析結果。在本研究中[3],我們利用矩陣資料的結構特性建立統計模型,並發展一套基於低秩 (low-rank) 結構的估計與檢定方法,以處理高維度問題所帶來的挑戰。此方法的優點之一在於,能夠像傳統方法一樣同時進行估計和假設檢定,進而提升估計效率和檢測能力。另一個優點是,我們提出的方法避免了現有方法中必須對模型參數施加可識別性 (identifiability) 限制的問題。
模型的建立
-
Matrix-covariate regression
令Y∈R為反應變數, M∈R^(p×q) 為矩陣型態資料,及 Z∈R^m 為其他可能的干擾變數 (confounding factors)。我們的研究重點在於探討 Y 與 M 之間的關聯,具體而言,包括:(1) Y 與 M 之間是否存在關聯;以及 (2) M 如何影響 Y 。接下來,我們透過統計模型的建立來回答這些問題。為了保留資料的矩陣結構,我們將廣義線性模型改寫為矩陣形式:g{E(Y|Z,M)}=γ+ξ^⊤ Z+vec(η)^⊤ vec(M) ,其中 g(∙) 為連結函數 (link function), γ 為截距項, ξ 為 Z 的迴歸係數,而 η 為對應於矩陣資料 M 的迴歸係數矩陣; vec(M) 表示將 M 的各行依序堆疊而成的 pq 維向量。在上述矩陣形式的GLM 架構下,回答問題 (1) 與 (2) 的核心在於對矩陣參數 η 進行估計,並檢定其是否為零矩陣,即 H_0 :η=0_(p×q) - Rank-r GLM我們所提出方法的基本理念是假設 M 的大多數列與/或行特徵在解釋Y時並不具有實質作用,這表示 η 中與這些特徵對應的大多數元素為零。基於此稀疏 (sparse) 結構假設,η 可合理視為一個低秩矩陣,因此,採用秩為 r 的 GLM作為建模工具是合理的,即 g{E(Y|Z,M)}=γ+ξ^⊤ Z+vec(η)^⊤ vec(M),η=AB^⊤ ,其中 A∈R^(p×r) , B∈R^(q×r) ,且 1≤r≤min{p,q} ,使得 rank(η)=r 。此外, (A,B) 本身並不具有可識別性。由於我們感興趣的是 η ,而非 (A,B),因此,不需要對 (A,B) 施加額外限制;只需要 AB^⊤ 具有可識別性,即可建立對 η 的推論程序。此低秩模型的一項優勢在於其參數化方式具有精簡性 (parsimony);當 r 較小時,利用此模型進行配適通常可提升估計效率。
估計及檢定流程
- 參數的估計
給定資料 {(Y_i,Z_i,M_i )}_(i=1)^n ,在 rank-r GLM 下,我們感興趣的參數為 β=β(θ)=(γ, ξ^⊤,vec(AB^⊤ )^⊤ )^⊤ ,其中 θ=(γ,ξ^⊤,vec(A)^⊤,vec(B)^⊤ )^⊤ 為模型參數。為了讓估計過程更為穩定,常見且廣泛採用的作法是以懲罰最大概似估計法 (penalized maximum likelihood estimation, penalized MLE) 來估計 θ 。在本研究中,我們使用 ‖A‖_F^2 ‖B‖_F^2 作為懲罰項,即 θ ̂=argmax┬θ{l(θ)-λ/2 ‖A‖_F^2 ‖B‖_F^2 } ,其中 l(θ) 為對數概似函數 (log-likelihood function), λ≥0 用以控制懲罰項的影響程度,而 ‖∙‖_F 表示 Frobenius 範數 (norm)。得到 θ 的估計值後,便可進一步得到的估計值 β ̂=β(θ ̂)=(γ ̂, ξ ̂^⊤,vec(η ̂ )^⊤ )^⊤ ,η ̂=A ̂B ̂^⊤ 。
我們採用交替法 (alternating method) 求解 θ ,其詳細作法請參見文獻[3]。 - 低秩檢定統計量及 p-值
針對虛無假設 (null hypothesis) H_0:η=0_(p×q) ,我們提出三種檢定統計量,以因應 η 在不同結構下可能呈現的訊號特性。
1. Wald 檢定統計量: T_wald=vec(η ̂ )^⊤ {[Σ ̂ ]_η/n}^+ vec(η ̂ ) ,其中 Σ ̂ 表示 β ̂ 的漸近共變異矩陣, [Σ ̂ ]_η 表示 Σ ̂ 中對應於 vec(η ̂ ) 的漸近共變異矩陣,而 (∙)^+ 表示 Moore-Penrose 廣義反矩陣。由於此統計量是對 η ̂-0 的加權總和,因此,當 η 較為密集 (dense) 時,使用 T_wald 較為適合;若 η 具有稀疏結構,則部分訊號可能在加總時被平均掉,進而降低檢定力。
2. Max 檢定統計量: T_max=max┬(l∈{m+2,…,1+m+pq})〖(β ̂_l^2)/([Σ ̂ ]_(β_l )/n)〗,其中最大值取自 η 各參數的估計量。此統計量著重於最顯著的單一參數,因此,當 η 為稀疏矩陣時,通常具有較佳的表現。
3. 乘積型檢定統計量: T=T_wald∙T_max∙T_gesat ,其中 T_gesat [5]利用 variance-component test [4] 的概念來處理高維度問題;當 η 的效果較弱 (weak effect) 時, T_gesat 預期表現較佳。由於實際應用中通常無法事先得知 η 的結構,而 T_wald 、 T_max 與 T_gesat 各有其適用情境,因此,我們將三種統計量結合,並提出此乘積型統計量。當 T 的值越大時,越傾向拒絕虛無假設 H_0 。
除了檢定整體假設 (overall hypothesis) 外,也可以透過 pq 個個別假設 (individual hypotheses),進一步檢定 M(j,k) 中對反應變數具有顯著影響的元素: H_0^((jk)):η_jk=0,1≤j≤p,1≤k≤q 。
由於 T 的虛無分配 (null distribution) 不易直接推導,我們根據資料的結構考慮兩種重抽樣方法 (re-sampling procedures)。在 rank-r GLM 下,我們採用參數化拔靴法 (parametric bootstrap) [2] 來計算其 p-值 (p-value);而當模型可簡化為 g{E(Y|M)}=γ+vec(η)^⊤ vec(M), η=AB^⊤ 時 (亦即不含干擾變數 Z ),則可進一步建構精確檢定 (exact test)。在此情況下,可透過隨機排列 (randomly permuting) Y_i 來破壞其與 M_i 之間的關聯,並以此產生虛無假設下的資料 (null data)。此外,這兩種重抽樣方法也可以用於計算 T_wald 、 T_max 與個別檢定的p-值。
應用於分析腦波圖 (Electroencephalography, EEG)
接下來,我們將所提出的估計及檢定方法應用於EEG資料分析。EEG資料來自 UCI Machine Learning Repository的EEG Database [1],用以探討腦波訊號與遺傳性酒精成癮之間的關聯。此資料共有122位受試者納入分析,其中 77 位為酒精成癮者 (Y=1) ,45位為非酒精成癮者 (Y=1) 。每位受試者在三種不同情況下 (single stimulus、two matched stimuli、two non-matched stimuli),共完成120次腦波檢查;在每次檢查中,皆記錄64個通道 (channels) 於256個時間點的電壓值。在本研究中,我們僅考慮single stimulus的情況,並以每位受試者於各次檢查中所測得電壓值的平均作為分析資料,因此,每位受試者對應一個256×64的矩陣型態資料。
我們的研究目的在於回答以下問題:(1) 透過整體假設檢定,評估腦波訊號是否與酒精成癮狀態具有顯著關聯;(2) 透過個別假設檢定,進一步辨識具有顯著影響的變數 (即特定時間區段下的特定電極通道)。
在分析中,我們先對資料進行前處理,將每8個時間點的電壓值取平均,經過此前處理後,每位受試者對應一個32×64的矩陣型變數 M 。此外,我們再對 M 的各元素進行標準化,使每個 M(j,k) 皆具有平均數0及變異數1。接著,以 r=1 的rank-r GLM對 (Y,M) 進行模型配適,並計算統計量 T 的 p-值。結果顯示,p-值小於 10^(-3) ,表示在給定顯著水準為0.05下,腦波訊號與酒精成癮狀態之間具有顯著的關聯。在控制整體型一誤差率 (family-wise error rate, FWER) 為0.05下,我們進一步進行個別檢定,來辨識 M 中具有顯著影響的變數,分析結果如表一所示。根據辨識結果,顯著通道主要分布於頂葉區 (Parietal) 與枕葉區 (Occipital),這兩個腦區分別與感覺功能及視覺功能有關,且在酒精成癮狀態中扮演重要角色。此外,大多數被辨識出的顯著通道集中於第 11 個時間區段,顯示大腦在刺激初期即產生明顯反應。
表一:整體型一誤差為0.05下,個別檢定所辨識出EEG訊號中的顯著變數。
|
Region |
Channel |
Time period |
|
|
Parietal |
20 |
(CP6) |
11 |
|
20 |
(CP6) |
18 |
|
|
22 |
(CP2) |
11 |
|
|
50 |
(CP4) |
11 |
|
|
24 |
(P4) |
11 |
|
|
26 |
(P8) |
11 |
|
|
27 |
(P7) |
11 |
|
|
27 |
(P7) |
12 |
|
|
51 |
(P5) |
11 |
|
|
52 |
(P6) |
11 |
|
|
Occipital |
28 |
(PO2) |
11 |
|
55 |
(PO7) |
11 |
|
|
56 |
(PO8) |
11 |
|
|
31 |
(O1) |
11 |
|
|
Central |
42 |
(C6) |
11 |
|
Temporal |
46 |
(TP8) |
11 |

圖一:人類大腦區塊及64個EEG電極通道。
總結
從統計方法的發展到實際應用,我們所關注的不僅是資料分析本身,更重要的是如何透過資料科學提升生物醫學與公共衛生研究的解釋力與決策力。以本研究的 EEG 分析為例,所提出的方法不僅能辨識腦波訊號與酒精成癮狀態之間的關聯,也有助於找出具有潛在生物意義的重要腦區與關鍵時間區段。這類方法的發展與應用,將有助於提升疾病風險辨識與早期偵測的能力,並進一步將統計研究成果轉化為具體可行的行動,實踐 Public Health in Action 的核心精神。
參考文獻
- Begleiter, H. (1995). EEG Database [Dataset]. UCI Machine Learning Repository. https://doi.org/10.24432/C5TS3D.
- Bůžková, P., Lumely, T., and Rice., K. (2011). Permutation and parametric bootstrap tests for gene-gene and gene-environment interactions. Annals of Human Genetics 75, 36-45.
- Hung, H., and Jou, Z. Y. (2019). A low rank-based estimation-testing procedure for matrix-covariate regression. Statistica Sinica, 29(2), 1025-1046.
- Lin, X. (1997). Variance component testing in generalised linear models with random effects. Biometrika 84, 309-326.
- Lin, X., Lee, S., Christiani, D. C., and Lin, X. (2013). Test for interactions between a genetic marker set and environment in generalized linear models. Biostatistics 14, 667-681.
----------------------------------------------------------------------
公衛四季:臺大公衛電子報第021期【2026春季號】文章 延伸閱讀