面向多元時間序列的群體因果關系發現算法

蔡瑞初，伍運金，陳薇，郝志峰，2

（1.廣東工業大學 計算機學院，廣州 510006；2.汕頭大學 理學院，廣東 汕頭 515063）

0 概述

因果發現旨在從觀測數據中發現變量之間的因果關系，可以揭露數據的生成機制，幫助人們理解數據，輔助人們進行干預和決策［1］。近年來，因果關系在深度學習［2］、金融經濟［3］、神經科學［4］、生物信息學［5］、社會科學［6］等領域受到了廣泛關注。

當觀測數據是多元時間序列數據時，現有的時序因果發現算法［7-10］通常認為個體之間是獨立的，為每一個個體的多元時間序列數據單獨學習一個因果關系作為該個體背后的因果關系，而個體間的因果關系學習過程是彼此無關的。然而在實際中，個體之間可能存在相同的因果關系。舉例來說，在電商場景下，同一群體中的個體的購買行為可能具有相同的思維方式，如對于家庭群體的個體是否購買某個商品的影響變量是商品的質量，而對于普通家庭群體的個體是否購買某個商品的影響變量則是商品的質量和價格。

因此，來源于不同群體的個體的數據背后會具有不同的因果關系（產生機制），而相同群體的個體的數據背后會具有相同的因果關系。如果能知道哪些個體屬于同一群體，便能利用群體內多個獨立同分布的個體數據一起學習該群體共同的因果關系。然而，在實際中無法預先知道不同個體是否屬于相同群體，也無法判斷總共有多少個群體。

針對上述問題，本文提出一種面向多元時間序列數據的群體因果關系發現算法。首先基于因果關系的相似性，將所有個體劃分成多個群體且無須指定群體的個數。對于每一個群體，使用變分推斷方法學習群體因果關系，從而充分利用多個個體數據。當所有個體均在一個群體時，該算法利用所有個體數據共同學習一個因果關系。當一個群體內只有一個個體時，該算法與現有時序因果關系發現算法類似，僅利用單個個體數據單獨學習一個因果關系。

1 相關工作

因果關系發現算法按照觀測數據的類型可以分為基于非時序數據的因果發現方法和基于時序數據的因果發現方法［11］。

基于非時序數據的因果發現方法中包括基于約束的方法［12］、基于評分的方法［13］和基于函數的方法［14-16］。基于約束的方法利用（條件）獨立性檢驗來判斷變量之間是否存在因果關系，而基于評分的方法通過給DAG 打分并尋找得分最高的DAG 作為變量間的因果關系，但2 種方法都存在馬爾可夫等價類的問題。為了解決這個問題，學者們提出了基于函數的方法，此類方法從數據產生機制出發，假設原因變量與結果變量存在函數映射，以及存在與原因變量獨立的噪聲變量，通過原因與噪聲的獨立性來識別因果關系。基于函數的方法包括線性非高斯無環模型（Linear Non-Gaussan Acyclic Model，LiNGAM）［14］、非線性加性噪聲模型（nonlinear Additive Noise Model，ANM）［15］和后非線性因果模型（Post-NonLinear causal model，PNL）［16］。

上述基于非時序數據的因果發現方法也被拓展到了時序數據上，如同樣基于約束的PCMCI算法［8］、基于評分的DYNOTEARS算法［9］和基于函數的VAR-LiNGAM算法［10］。PCMCI 算法基于條件獨立性測試框架，使用PC 算法［17］發現變量的馬爾可夫等價類集合從而縮短條件集，進一步用瞬時條件獨立性（Momentary Conditional Independence，MCI）檢驗降低誤發現率。DYNOTEARS算法從優化問題的角度出發，最小化一個帶無環約束的損失函數，從而學習變量之間的瞬時影響和時延影響。VAR-LiNGAM 算法結合自回歸模型和線性非高斯無環模型來識別變量的瞬時和時延因果關系的影響權重，并通過非高斯性假設保證了算法的可識別性。

目前一些工作［4，18］也考慮到了不同個體的樣本背后的因果關系可能存在一定共性，并嘗試從這種混雜樣本中將樣本劃分成不同類別并學習因果關系。文獻［4］與本文工作同樣是面向多元時間序列數據，但其需要指定群體個數，并在此基礎上學習個體特定的因果關系和共性的因果關系，再基于共性的因果關系計算某個個體屬于不同群體的概率。文獻［4］與本文工作的不同之處在于，其學習的是個體個性的因果關系以及個體間共性的因果關系，認為不同個體背后的因果關系仍然是不同的，僅是存在一定的共性因果關系，并基于共性進行聚類。本文工作則是在個體間的因果關系可能相同的場景下，對多個個體聚類并學習群體因果關系。此外，在實際應用中，群體個數往往是未知的，而本文所提出的算法無須指定群體個數。文獻［18］考慮的則是二元變量之間的因果關系，在為每個樣本識別了因果關系后，基于每個樣本所對應的因果關系參數進行K-Means［19］聚類，因此也面臨著需要指定群體個數的問題。

2 問題定義

本節對所研究的問題進行符號化定義和說明。定義n個個體的多元時間序列數據集X={X1，X2，…，Xn}，其中第s個個體的多元時間序列數據Xs={，，…，}，且每個個體的變量 數均為m，時間序列長度均為T。將最長的因果關系時間間隔記為k，多元時間序列數據的因果關系表示為k+1個m×m的矩陣{B0，B1，…，Bk}，Bτ[i，j]≠0 表示第t時刻的變量xi(t)受到第t-τ時刻的變量xj(t-τ)的因果影響，且時間間隔為τ，τ∈{0，1，…，k}，其中τ=0時Bτ表示瞬時因果關系，τ＞0時Bτ表示時延因果關系。

本文所考慮的問題是：給定多元時間序列數據集X，如何基于數據背后的因果關系將n個個體劃分成c個群體，且c無須人工先驗指定，并學習每個群體的因果關系。

3 LEAD 算法

針對所研究問題，本文提出面向多元時間序列的群體因果關系發現算法LEAD。該算法分為2 個階段：第一階段是基于因果關系對全體的多元時間序列數據進行個體級別的聚類，將具有相同因果關系的個體聚集成一個群體，并且無須指定群體的個數；第二階段基于前一階段的聚類結果，利用每個群體中所有個體的多元時間序列數據進行群體因果關系發現。圖1 展示了LEAD 算法的整體框架。

圖1 LEAD 算法框架Fig.1 Framework of LEAD algorithm

下面將先討論如何基于變分推斷的方法，利用個體的多元時間序列數據學習該個體背后的因果關系，該方法在LEAD 算法的2 個階段中都會涉及，隨后將對2 個階段的具體流程進行介紹。

3.1 基于變分推斷的個體因果關系發現算法

本文主要考慮變量之間的因果關系服從線性非高斯的情況，并采用函數模型進行建模。第s個個體的多元時間序列數據的生成機制如下：

其中：N(·)表示一個高斯分布。因此，主要的推斷任務是估計真實后驗分布p(Bs|Xs)。然而直接計算該后驗分布十分困難，本小節采用變分推斷的技巧，利用一個易于計算的變分分布q(Bs)去近似真實后驗分布p(Bs|Xs)。首先給出第s個個體的觀測數據Xs的對數似然度：

從上述推導中可以看出，觀測數據Xs的對數似然度存在一個下界，當且僅當KL(q(Bs)||p(Bs|Xs))=0時，該下界等式成立。因此，只需要最大化這個下界，就可以找到一個變分分布q(Bs)足夠近似真實后驗分布p(Bs|Xs)，由此得到以下目標函數：

通過優化變分分布q(Bs)來最大化式（4）即可近似真實后驗分布，從而得到第s個個體背后因果關系。

為了計算式（4），對其進行分解如下：

將等式最后的第一項記為Lell(q(Bs)，Xs)，第二項記為Lkl(q，p)，即：

可以看出：Lell(q(Bs)，Xs)是給定變分分布q(Bs)下觀測數據Xs的對數似然度的期望；Lkl(q，p)是近似后驗分布q(Bs)與先驗分布p(Bs)的KL 散度的負數。因此，最大化式（4）意味著最大化觀測數據Xs的對數似然度且同時最小化q(Bs)與先驗分布的差距。

式（4）的計算可以通過式（6）和式（7）實現，因此，對變分分布q(Bs)做變分推斷中常見的平均場假設，即：

下面給出lnp(Xs|Bs)的具體計算方式。由數據的產生機制（式（1））可知：

進一步用混合高斯分布來對噪聲項進行建模，即：

對于式（7），先驗分布p(Bs)和變分分布q(Bs)都是高斯分布，因此，可以通過解析式直接計算2 個高斯分布的KL 散度。基于上述推導分析，可以對目標函數式（4）進行估計并計算梯度，進而采用梯度上升的方式對變分分布q(Bs)和噪聲分布pe(es(t))中的參數進行優化，從而最大化目標函數，得到第s個個體背后的因果關系Bs。

3.2 基于因果關系的聚類算法

上節給出了如何基于變分推斷的方法，利用個體的多元時間序列數據學習該個體背后的因果關系，本節將介紹如何基于學習到的個體因果關系進行聚類。

設計聚類算法的目的是為了將具有相同或者相似因果關系的個體聚成一個群體，從而利用群體的數據來學習群體背后的因果關系。因此，本文從因果關系的角度來度量個體之間的相似性并設計個體聚類過程，提出基于因果關系的時間序列聚類算法，如算法1 所示。

算法1基于因果關系的時間序列聚類算法

算法2Split(Cv，Lv)函數算法

算法1 的設計主要借鑒了分裂層次聚類的思想，初始時將全部個體作為一個群體，隨機選擇一個個體Xs作為群體代表進行因果關系學習，以學習得到的變分分布q(Bs)作為群體共同的因果關系，并計算每個個體在群體因果關系下的似然度（算法第1～7 行）。個體似然度高意味著該個體的因果關系與該個體所在群體對應的群體共同因果關系相似；個體似然度低意味著該個體的因果關系與該個體所在群體對應的群體共同因果關系不相似。隨后算法進入循環過程（算法第8～18 行），在每一次循環中，需要從當前已有的群體中選擇出群體內部個體似然度方差最大的群體。群體內部個體似然度的方差越大，意味著該群體內部個體的數據背后越可能是不同的因果關系，導致在相同的群體因果關系下，有些個體似然度高，有些個體似然度低。在此基礎上，基于選出的群體，利用算法2 對其做分裂，并根據劃分前后群體所有個體的似然度之和是否增加來決定是否接受這次分裂，如果總體似然度之和沒有增加，則拒絕分裂且跳出循環，此時算法1 返回聚類結果并結束（算法第19 和20 行）。

由于分裂后的2 個群體分別學習各自的因果關系并在對應群體因果關系下計算該群體內的個體的似然度，無法直接比較分裂前后群體方差的變化，同時本文希望個體似然度高，因此算法1 以群體內部的所有個體似然度之和是否增加作為是否分裂的條件。同時，算法1 在拒絕分裂后隨即跳出循環結束算法流程具有2 點好處：一方面緩解了經典分裂層次聚類算法需要分裂到單個個體作為一個簇才停止的情況，降低了時間復雜度；另一方面無須為算法指定群體的個數，循環停止時即得到了群體的個數，避免了經典層次聚類算法結束后仍需要選擇簇個數的問題。

3.3 群體因果關系發現算法

上節基于因果關系將全部個體劃分成了多個群體，并且無須指定群體的個數，本小節將基于前一階段的聚類結果，利用每個群體中的所有個體的多元時間序列數據進行群體因果關系發現。

本節所提出的群體因果發現算法主要基于貝葉斯思想，即利用數據更新先驗分布得到后驗分布并循環這一更新過程，算法偽代碼如算法3 所示。

算法3群體因果關系發現算法

算法3對于每個群體進行遍歷，在同一個群體內，先初始化該群體因果關系的先驗分布（算法第2 和3 行），再對該群體內的個體進行遍歷，對于一個個體，利用該個體數據對式（4）最大化得到近似后驗分布，并以此近似后驗分布作為下一次的先驗分布，重復這個過程直到所有個體都參與過群體因果關系分布的更新（算法第4 和7 行），以最終得到的近似后驗分布的期望作為該群體的群體因果關系并保存結果（算法第8～9 行），直到遍歷完所有群體。通過上述步驟就可以學習出各個群體共享的群體因果關系。得益于3.1 節中基于變分推斷的個體因果關系發現算法與貝葉斯思想，算法3 能夠充分地利用群體內所有個體的數據，并且學習出群體背后的因果關系。

4 實驗

為了驗證本文提出的LEAD 算法，本節將使用仿真數據和真實數據對算法進行實驗評估。

4.1 仿真數據集

設計4 組控制變量實驗，具體如下：

1）群體個數c={2，3，4，5，6}，個體總數n=60，個體的時間序列長度T=60，變量個數m=5。

2）群體個數c=2，個體總數n={20，30，40，50，60}，個體的時間序列長度T=60，變量個數m=5。

3）群體個數c=2，個體總數n=30，個體的時間序列長度T={40，60，80，100，120}，變量個數m=5。

4）群體個數c=2，個體總數n=30，個體的時間序列長度T=60，變量個數m={6，8，10，12，14}。

本文所提出的LEAD 算法對多個個體進行聚類并學習個體背后的因果關系（即該個體所在群體的群體因果關系），因此，本節將其與主流聚類算法對比聚類的效果，與時間序列因果發現算法對比因果發現的效果。

聚類對比算法如下：1）K-Means（Euclidean），基于歐氏距離的K-Means 算法［19］；2）K-Means（DTW），基于DTW 距離［22］的K-Means 算法；3）DBSCAN（DTW），基于DTW 距離的DBSCAN 算法［23］；4）OPTICS（DTW），基于DTW 距離的OPTICS 算法［24］。上述聚類對比算法采用Python 包scikit-learn 中的實現方法，參數設置為各算法實現的默認參數。其中，K-Means 算法需要指定簇的個數K，本小節實驗將按真實值來指定K，注意K在實際中往往是未知的。

因果關系發現對比算法如下：1）基于約束的PCMCI 算法［8］；2）基于評分的DYNOTEARS 算法［9］；3）基于函數的VAR-LiNGAM 算法［10］。上述時間序列因果發現對比算法均采用官方開源代碼實現，并按真實值設置對比算法的最大因果關系時間間隔參數，其余參數采用原算法實現的默認參數。

在評價指標方面，本實驗選擇ARI（Adjusted Rand Index）指標［25］對聚類算法進行評價，選擇AUC（Area Under Curve）指標對因果關系發現算法進行評價。ARI 計算公式為：

其中：n是給定數據集中個體總數；nij是在聚類結果R中屬于簇i同時在聚類結果R*中屬于簇j的個體數量；ni是在聚類結果R中屬于簇i的個體數量；nj是在聚類結果R*中屬于簇j的個體數量；

ARI 的取值范圍是［-1，1］，隨機劃分下ARI 為0，ARI 越大說明R和R*結果越一致。AUC 是ROC曲線下的面積，ROC 曲線是以假正率（FPR）為橫軸，真正率（TPR）為縱軸，在多個閾值下得到一系列點（FPR，TPR）所構成的曲線，其中FPR 和TPR 計算公式為：

其中：將預測因果結構與實際因果結構相比；FP是預測有邊但實際無邊的數量；TN是預測無邊實際也無邊的數量；TP是預測有邊實際也有邊的數量；FN是預測無邊但實際有邊的數量，AUC 的取值范圍是［0，1］，AUC 越接近1，預測因果結構與真實因果結構越接近。

在仿真實驗中，相同參數的實驗都在10 個不同隨機種子下進行數據生成并運行LEAD 算法和對比算法，分別計算評價指標并取10 次結果的平均值作為算法的最終結果，結果如圖2～圖5 所示。

圖2 針對群體個數的控制實驗Fig.2 Control experiments on the number of groups

圖2 顯示了個體總數相同但背后群體個數不同情況下的實驗結果。可以看出本文提出的LEAD 算法的因果發現效果要優于對比算法。但隨著群體個數的增加，LEAD 算法的因果發現效果下降，這是因為實際群體個數越多，聚類的難度會越大，特別是在群體個數未知的情況下難度會更大，因此LEAD 算法第一階段的聚類效果變差，導致第二階段的因果發現效果也下降。事實上，圖2（b）中的聚類方法都隨著群體個數的增加而效果變差，其中K-Means（Euclidean）和K-Means（DTW）由于給定了真實群體個數作為算法參數，因此效果要優于其他方法。而本文提出的LEAD 算法在沒有給定真實群體個數的情況下，也取得了相近的聚類效果。

圖3 顯示了群體個數相同但個體總數不同情況下的實驗結果。可以看出LEAD 算法的因果發現能力遠高于對比算法，這是因為LEAD 算法可以利用群體內多個個體數據。隨著個體總數增加，聚類的難度增大，LEAD 算法的聚類效果有所影響，但其因果發現效果基本沒有影響，甚至因為數據增加而略微改善。

圖3 針對個體總數的控制實驗Fig.3 Control experiments on the total number of individuals

圖4 顯示了不同時間序列長度下的實驗結果。隨著序列長度增加，圖4 中的因果發現算法和聚類算法效果都有所改善，且LEAD 算法無論是因果發現還是聚類的效果均要優于對比算法，說明LEAD算法不僅具有優異的因果發現能力，而且還具有很好的多元時間序列聚類能力。

圖4 針對時間序列長度的控制實驗Fig.4 Control experiments on the length of time series

圖5顯示了不同變量個數下的實驗結果。從圖5（b）可以看出，隨著變量個數增加，個體間的差異信息越來越豐富，聚類的難度會下降，聚類算法的效果均明顯改善。但與此同時，變量間的因果關系也更復雜，因果發現的難度增大，因此圖5（a）中的因果發現算法均隨著變量個數增加而效果變差，但LEAD 算法仍然優于對比算法，驗證了算法的有效性。

圖5 針對變量個數的控制實驗Fig.5 Control experiments on the number of variables

4.2 真實數據集

本節選擇在真實的Sachs 數據集［5］進行測試。Sachs 數據集在不同的干預措施下測量了細胞中11種磷酸化蛋白質和磷脂分子的濃度變化，不同干預措施下這11 種分子的因果關系會發生變化［4］，本節使用干預cd3cd28+U0126 和干預cd3cd28+aktinhib下的數據，2 種干預下的數據分別屬于2 個群體，并以30 作為一個個體的時間序列長度對數據進行切分，得到每個25 個個體的數據。

將LEAD 算法和聚類算法應用到該數據集中，并展示LEAD 算法所學習出的因果關系。由表1 可以看出，本文提出的LEAD 算法的聚類效果要優于對比算法，因為LEAD 算法能很好地捕捉到因果關系上的變化而其他距離度量方式不容易捕捉到這種變化。進一步分析聚類結果發現，LEAD 算法很好地將屬于干預cd3cd28+aktinhib 組的25 個個體劃分在同一群體（記為群體a1）中，但是將屬于干預cd3cd28+U0126 組的25 個個體分成了2 個群體，其中一個群體（記為群體U1）有21 個個體，另一個群體（記為群體U2）有4 個個體。

表1 在真實數據集中的聚類性能 Table 1 The clustering performance on real dataset

本文算法學習的因果關系如圖6 所示，其中有向邊a1 為群體a1 中的因果關系，有向邊U1 為群體U1 中的因果關系，有向邊U2 為群體U2 中的因果關系，僅顯示每個群體中影響權重前十的因果關系。可以看出，群體a1 的重要因果關系與群體U1、群體U2有較大區別，而群體U1、U2的重要因果關系有5條相同，說明算法能夠區分在不同干預下的個體。此外，在群體U1、U2 中都發現了因果關系Erk →Akt，該因果關系與文獻［5］中報告的一致，但在群體a1 中并未被發現，這是因為在群體a1 中Akt 被干預，從而切斷了Akt 的原因變量，驗證了算法的有效性。

圖6 真實數據集中學習到的3 個群體因果關系Fig.6 Three collective causal relations learned from real dataset

5 結束語

現有因果關系發現方法在多個個體數據背后有相同因果關系的情況下樣本利用不足。針對該問題，本文提出一種面向多元時間序列的群體因果發現算法。該算法通過基于因果關系的聚類，將多個個體按照因果關系的相同或不同劃分成多個群體，且算法無須指定群體個數。在此基礎上，為每個群體通過群體因果發現算法學習群體因果關系。實驗結果表明，該算法可以充分利用具有相同因果關系的多個個體數據，因果發現能力優于對比算法，并且同時具有和對比算法相近的多元時間序列聚類能力。本文關注于因果關系是線性的情況，下一步將考慮在因果關系滿足非線性的情況下如何進行群體因果發現，進一步提升算法的適用范圍。