針對肝病分類中存在的特征交互的問題, 我們研究了一種分層交互lasso分類方法。首先對logistic模型添加lasso罰函數和分層凸約束, 其次采用卡羅需-庫恩-塔克條件與廣義梯度下降法相結合的凸優化方法給出模型求解方法, 最后得到主效應特征系數與交互特征系數的稀疏解, 實現模型分類。本文在兩個肝病數據集上進行實驗, 證明了特征交互對肝病分類有貢獻。實驗結果證明了分層交互lasso方法可解釋性強, 效果、效率均優于lasso方法、全特征對lasso方法以及支持向量機、最近鄰和決策樹等傳統分類方法。
引用本文： 王金甲, 盧陽. 特征交互lasso用于肝病分類. 生物醫學工程學雜志, 2015, 32(6): 1227-1232. doi: 10.7507/1001-5515.20150218 復制
0 引言
肝病在我國是一種高發病，已經成為一種嚴重影響我國人民生命健康的疾病。我國政府對肝病的防治尤為重視，從七五計劃起一直將肝病的防治作為連續攻關課題。故此本文作者著力于肝病分析研究，以尋找有效的肝病檢測治療方案，實驗數據來自UCI學習數據集。
在過去的研究工作中，眾多學者選用1990年Forsyth收集的肝臟失調樣本(BUPA Liver Disorders, BLD)數據集作為研究對象，以探索肝病的病理特征。研究者探索方向大致分為兩類。第一類為機器學習算法[1-2]，例如實例推理法、決策樹、貝葉斯分類、人工神經網絡以及基于模糊邏輯的實例學習。Neshat等[3-4]將模糊方法和霍普菲爾德神經網絡(Hopfield neural network)方法應用于肝臟異常診斷。這些傳統方法在一定程度上提高了肝病診斷的準確性，但可解釋性都較差，且忽略了肝病特征交互。另一類方法采用優化理論，如Lin等[5]采用粒子群智能優化算法得到78.18%的正確率，Polat等[6]采用監督人工免疫系統提高診斷正確率。此兩種方法雖然提高了肝病診斷的準確率同時也具有一定的可解釋性，但復雜度較高，不適用于現代醫學應用。Chien等[7]提出了平滑套索方法，得到了68.82%的正確率，建立的模型解釋性較強，但所得分類結果不夠理想，并且忽略了特征交互對病癥的影響。
2012年Ramana等收集了印度肝病數據集(Indian Liver Patient Dataset，ILPD)，它更貼近這些年來肝病的發展現狀。Karthik等[8]采用人工神經網絡分類ILPD肝病樣本，方法容易過擬合，也忽略了疾病間的交互特征。Neshat等[9]則將粒子群優化算法與案例推理結合進行特征提取和肝病分類，此方法屬于優化算法，也沒有利用疾病間的特征交互信息。Ekong等[10]通過模糊聚類的方法來分析血清蛋白，并以臨床特征作為輔助以診斷肝病。Lin[11]綜合了多種分類器、回歸樹和案例推理技術，構建了一種智能診斷模型。Ekong與Lin將智能模型與肝病病理知識結合，提高了肝病的分類正確率，但模型復雜度較高且忽略了交互信息在醫學中的重要意義。而套索(least absolute shrinkage and selection operator，lasso)方法用于ILPD的處理則未見報道。
基于上述方法模型復雜度高以及對特征交互信息利用率低的問題，本文以簡單實用的logistic回歸模型為基礎，結合文獻[12]提出的交互lasso模型和凸優化算法，研究了數據驅動的罰分類方法，即改進的lasso方法用于肝病診斷。lasso方法中響應變量是預測變量線性加權和的表示，因此模型具有很強的可解釋性。而特征間的復雜交互對于疾病診斷具有重要意義，這推動我們將特征交互思想用于疾病診斷。作者首先采用分層lasso罰分類模型，并添加凸約束條件，對主效應特征和交互特征進行收縮和選擇處理，而后通過卡羅需-庫恩-塔克最優化條件(Karush-Kuhn-Tucker，KKT)與廣義梯度下降方法結合解決最優化問題，并求解模型參數。本文的創新在于使用了分層特征交互思想，將可解釋性強的lasso分類方法推廣用于肝病診斷。最后本文對兩組肝病數據進行分類實驗，并與lasso方法、全特征對lasso方法、支持向量機(support vector machine，SVM)、線性判別分析(linear discriminant analysis，LDA)等其他傳統分類方法進行比較。新方法用于兩組數據的實驗結果表明了特征交互lasso分類方法在肝病診斷中的有效性。
1 交互lasso模型
假設模型樣本輸出為Y，輸入X為P維特征x1, …，xj, …, xP，樣本在每維輸入特征間都有交互xjxk，則考慮交互的logistic回歸模型具有以下形式
| $\begin{gathered} \log {\text{it}}\left({P\left({Y=1\left| X \right|} \right)} \right)=\sum\limits_{j=0}^P {{\beta _j}{x_j}} + \hfill \\ \frac{1}{2}\sum\limits_{j \ne k} {{\Theta _{jk}}{x_j}{x_k} + \varepsilon }, \hfill \\ \end{gathered} $ |
其中x0為1，xj表示1階主效應特征，xjxk表示2階交互特征，1階主效應特征系數為β∈RP+1，2階交互特征系數為Θ∈RP×P，Θjj=0，ε∈N(0, σ2)。可以看出式(1)是主效應特征和交互特征的線性加權和，模型可解釋性強。
假設訓練樣本為(x1, y1), …, (xi, yi), …, (xN, yN)，1階主效應特征和2階交互特征融合后特征數目為P+P(P-1)/2，我們目的是從融合特征中通過罰分類方法進行特征子集選擇，選擇出的特征子集作為相應的預測特征。這就需要估計出模型系數中非零參數值，因為它們對應被選擇的子集。定義，則兩類的概率分別如下：
| ${p_i}=p\left({{{\text{x}}_{\text{i}}}} \right)=p\left({Y=1|{{\text{x}}_{\text{i}}}} \right)=\frac{1}{{1 + \exp \left({-\sum\limits_j {{\beta _j}{x_{ij}}}-\frac{1}{2}\sum\limits_{j \ne k} {{\Theta _{jk}}{x_{ij}}{x_{ik}}} } \right)}}$ |
| $-{p_i}=1-p\left({{{\text{x}}_{\text{i}}}} \right)=p\left({Y=2|{{\text{x}}_{\text{i}}}} \right)=\frac{1}{{1 + \exp \left({-\sum\limits_j {{\beta _j}{x_{ij}}}-\frac{1}{2}\sum\limits_{j \ne k} {{\Theta _{jk}}{x_{ij}}{x_{ik}}} } \right)}}$ |
通過極大似然估計擬合模型，使N次獨立觀測的似然函數或者對數似然函數最大。只考察大的主效應特征間的交互更為高效，因此模型添加約束‖Θj‖1≤｜βj｜，即
| $\ell\left({\beta, \Theta } \right)=\sum\limits_{i=1}^N {\left[{{y_i}\log {p_i} + \left({1-{y_i}} \right)} \right]\log \left({1-{p_i}} \right)} $ |
其中Θj是Θ的第j行。如果jk≠0，那么有‖j‖1＞0，因此有j≠0。新增約束條件導致式(4)不能凸優化求解，因此使用β+, β-∈RP+1兩個向量代替β的策略使式(4)成為凸函數，并添加lasso罰函數[13]，即
| $\begin{gathered} \ell\left({{\beta ^ \pm }, \Theta } \right)=-\sum\limits_{i=1}^n {\left[{{y_i}\log {p_i} + \left({1-{y_i}} \right)\log \left({1-{p_i}} \right)} \right]} \hfill \\ + \lambda {1^T}\left({{\beta ^ + } + {\beta ^-}} \right) + \frac{\lambda }{2}{\left\| \Theta \right\|_1} \hfill \\ s.t.\; \;\left. \begin{gathered} {\left\| {{\Theta _j}} \right\|_1} \leqslant \beta _j^ + + \beta _j^-\hfill \\ \beta _j^ + \geqslant 0, \beta _j^-\geqslant 0 \hfill \\ \end{gathered} \right\}\; \;{\text{for}}\; \;j{\text{=1, }} \cdots {\text{, }}P, \hfill \\ {\beta _j}=\beta _j^ +-\beta _j^-, {\left\| {{\beta _j}} \right\|_1}=\beta _j^ + + \beta _j^-\hfill \\ \end{gathered} $ |
其中懲罰參數λ≥0，控制著訓練數據的相對重要性和稀疏性(l1罰項)。我們的目的是求解使式(5)最小化時的主效應特征系數與交互特征系數(β±, Θ)。
2 KKT條件與廣義梯度下降求解模型
我們用廣義梯度下降法[14]求解式(5)問題。
引理：設無約束優化問題有以下形式Minimize g(φ)+h(φ)，其中g和h為凸函數，則利用廣義梯度下降法可轉化為如下形式，這里k是迭代次數，t是步長：
| $\begin{gathered} {{\hat \varphi }_k}=\arg \min \frac{1}{{2t}}{\left\| {\varphi-{{\tilde \varphi }_{k-1}}} \right\|^2} + h\left(\varphi \right)=\hfill \\ \arg \min \frac{1}{{2t}}{\left\| {\varphi-\left[{{{\tilde \varphi }_{k-1}}-t\nabla g\left({{{\tilde \varphi }_{k-1}}} \right)} \right]} \right\|^2} + h\left(\varphi \right) \hfill \\ \end{gathered} $ |
證明顯而易見。
定理1：式(5)中令，兩式相加的無約束條件的全局最小值解為
| $\left\{ \begin{gathered} \tilde \beta _j^{ \pm \left({k-1} \right)}=\tilde \beta _j^{ \pm \left({k-1} \right)} \pm t\left[{\sum\limits_{i=1}^n {{x_{ij}}\left({{y_i}-p_i^{k-1}} \right) + \lambda } } \right] \hfill \\ \tilde \Theta _{jk}^{\left({k-1} \right)}=\hat \Theta _{jk}^{\left({k-1} \right)} + \frac{t}{2}\sum\limits_{i=1}^n {{x_{ij}}\left({{y_i}-p_i^{k-1}} \right){x_{ik}}} \hfill \\ \end{gathered} \right.$ |
證明：
令，對e±求βj+的偏導數，得，同理，。。
對g(β±, Θ)求βj+的偏導數，得。同理，對g(β±, Θ)求βj-的偏導數，得。
對g(β±, Θ)求Θjk的偏導數，得。
綜上：。
根據引理中可得：
| $ \begin{align} & \tilde{\beta }{{_{j}^{+}}^{\left(k-1 \right)}}\text{=}\tilde{\beta }{{_{j}^{+}}^{\left(k-1 \right)}}\pm t\left[\sum\limits_{i=1}^{n}{{{x}_{ij}}}\left({{y}_{i}}-p_{i}^{k-1} \right)+\lambda \right] \\ & \tilde{\Theta }_{jk}^{\left(k-1 \right)}\text{=}\tilde{\Theta }_{jk}^{\left(k-1 \right)}+\frac{t}{2}\sum\limits_{i=1}^{n}{{{x}_{ij}}}\left({{y}_{i}}-p_{i}^{k-1} \right)+{{x}_{ik}} \\ \end{align} $ |
根據引理中式(6)，根據特征維數可分解為P個獨立可并行解決的子問題，式(5)可轉換為如下凸優化形式：
| $\begin{align} & \underset{{{\beta }^{\pm }}\in {{R}^{p}}, \Theta \in {{R}^{p-1}}}{\mathop{Minimze}}\, \frac{1}{2t}{{\left\| \beta _{j}^{+}-\tilde{\beta }_{j}^{+} \right\|}^{2}}+\frac{1}{2t}{{\left\| \beta _{j}^{-}-\tilde{\beta }_{j}^{-} \right\|}^{2}}+ \\ & \frac{1}{2t}\left\| {{\Theta }_{j}}-{{{\tilde{\Theta }}}_{j}} \right\|_{F}^{2}+\frac{\lambda }{2}{{\left\| {{\Theta }_{j}} \right\|}_{1}} \\ & s.t.\ \ \ \ {{\left\| {{\Theta }_{j}} \right\|}_{1}}\le \beta _{j}^{+}+\beta _{j}^{-}, \beta _{j}^{+}\ge 0, \beta _{j}^{-}\ge 0, j=1, 2, \cdots, P \\ \end{align}$ |
定理2：設和為式(8)的解，則
| $\tilde{\beta }_{j}^{\pm }={{\left[\tilde{\beta }_{j}^{\pm }+t{{{\hat{\alpha }}}_{j}} \right]}_{+}}, {{{\hat{\Theta }}}_{j}}=S\left({{{\hat{\Theta }}}_{j}}, t\left(\lambda /2+{{{\hat{\alpha }}}_{j}} \right) \right)$ |
證明：式(8)拉格朗日運算定義如下：
| $ \begin{align} & L\left(\hat{\beta }_{j}^{\pm }, {{{\hat{\Theta }}}_{j}}; {{{\hat{\alpha }}}_{j}}, {{{\hat{\gamma }}}_{j\pm }} \right)=\frac{1}{2t}{{\left(\hat{\beta }_{j}^{+}-\hat{\beta }_{j}^{+} \right)}^{2}}+\frac{1}{2t}{{\left(\hat{\beta }_{j}^{-}-\hat{\beta }_{j}^{-} \right)}^{2}}+ \\ & \frac{1}{2t}{{\left({{{\hat{\Theta }}}_{j}}-{{{\tilde{\Theta }}}_{j}} \right)}^{2}}+\left(\frac{\lambda }{2}+{{{\hat{\alpha }}}_{j}} \right){{\left\| {{{\hat{\Theta }}}_{j}} \right\|}_{1}}-\left(\hat{\gamma }_{j}^{+}+{{{\hat{\alpha }}}_{j}} \right)\hat{\beta }_{j}^{+}-\left(\hat{\gamma }_{j}^{-}+{{{\hat{\alpha }}}_{j}} \right)\hat{\beta }_{j}^{-} \\ \end{align} $ |
其中是分層約束對偶特征，是非負約束對偶特征。
有KKT條件如下[9]：平穩性條件為，平穩性條件為，其中uj是絕對值函數關于的次梯度，互補松弛條件為。
現定義，其中S(·)是lasso領域常用的函數[16]。最后的KKT條件只包含。觀察發現f關于是非增的，基于其為分段線性函數，可求解出。根據以及有，將代入β±的平穩性條件可知，由平穩性條件可得。
綜上，對λ值的一組下降序列，就可以得到β，Θ的路徑解。我們采用網格法，對于固定的λ，求解系數的過程是一個循環嵌套過程。梯度下降法求解交互lasso模型算法歸結如下：已知X, Y，根據定理1中式(7)求解，根據定理2中式(9)求解，直至收斂。
3 實驗數據和結果分析
本實驗采用的數據是加州大學歐文分校(University of California Irvine, UCI)機器學習數據庫的BLD數據集以及ILPD數據集。實驗數據來自 http://archive.ics.uci.edu/ml/。實驗通過R語言編程建立模型并實現文中算法。
3.1 數據說明
BLD數據集包含樣本345個，其中健康樣本200個，肝病樣本145個。每個樣本包括6維實數特征，前5維皆來自受試者血液檢測分析，分別為平均紅細胞容積、堿性磷酸酶、丙氨酸轉氨酶、天冬氨酸轉氨酶以及谷氨酰基轉肽酶含量，最后1維為受試者每日平均飲酒量。
ILPD數據集，由三個印度教授收集自印度安得拉邦的東北部，包括樣本583個，其中肝病患者416個以及非肝病受試者167個。583個受試者中包含441名男性和142名女性。每個樣本包括10維實數特征，分別為年齡、性別、總膽紅素、直接膽紅素、堿性磷酸酶、丙氨酸轉氨酶、天冬氨酸轉氨酶、血清總蛋白、血清白蛋白、白蛋白與球蛋白之比(簡稱白球比例)。
3.2 特征子集選擇的實驗結果和分析
本文對BLD數據集以及ILPD數據集進行100次10倍交叉驗證，實驗記錄了交互lasso分類方法的非0特征系數個數、平均錯誤率、標準差、CPU時間、所選取的λ值。
表 1、2給出了交互模型在兩組肝病數據上的系數解。行表示是非零主效應特征系數，第1列表示主效應系數，其它列表示交互系數。表 1、2中特征系數非0時，即表示此特征被選擇。表 1、2中交互系數絕對值越大說明特征交互程度越高。例如，表 1中黑體數字表示堿性磷酸酶與丙氨酸轉氨酶的交互特征被選擇，且交互系數為0.649 5；同樣地，其它非0的交互特征也被選擇。
表1
分層交互模型用于BLD數據所得系數解
Table1.
Solution of BLD dataset obtained using hierarchical interactive model
在BLD數據集實驗中，表 1是λ值為2.37時得到的模型系數解。從表中可以看出模型選擇了全部主效應特征以及部分交互項，主要包括以下4個交互項：①堿性磷酸酶與丙氨酸轉氨酶交互項，此項交互不但排除了堿性磷酸酶升高原因為骨病的可能性，而且提高了肝細胞病變診斷的可信度。②堿性磷酸酶與天冬氨酸轉氨酶交互項，此項作用與前項相似。③丙氨酸轉氨酶與天冬氨酸轉氨酶交互項，因天冬氨酸在心肌含量最高，依次為肝、骨骼肌、腎、胰，只有與其他特征結合才能斷定天冬氨酸轉氨酶升高是肝病引發的。④受試者飲酒量與幾乎其他各項均存在交互，可見飲酒對肝功能影響甚大。
在ILPD數據集實驗中，表 2是最優λ值為1.52時得到模型的主效應系數解和部分交互系數解。下面討論表 2中5項主要交互：①直接膽紅素與丙氨酸轉氨酶交互項，組成交互項的兩個特征均為肝病診斷的有效識別依據，而兩種特征交互，使此交互特征更突出。②總膽紅素與直接膽紅素交互項，總膽紅素升高但直接膽紅素不升高的原因并非肝細胞病變，故排除了許多其他病因。③丙氨酸與天冬氨酸交互項，與表 1中BLD實驗相同，因天冬氨酸存在于心肌、肝、骨骼肌、腎、胰等，在其與丙氨酸交互后，強化了此交互特征在模型中的重要性。④血清白蛋白與白球比例交互項，根據數據說明中的肝病檢測先驗知識可知，白蛋白降低且白球比例降低是典型正交互項，可提高醫生的診斷把握。⑤堿性磷酸酶與丙氨酸轉氨酶交互項，此項交互不但排除了堿性磷酸酶升高原因為骨病的可能性，并且提高了肝細胞病變診斷的可信度。
表2
分層交互模型用于ILPD數據所得系數解
Table2.
Solution of ILPD dataset obtained using hierarchical interactive model
3.3 本文方法與基本lasso方法的結果比較
表 3給出了幾種lasso方法的對比結果。從表中可以看出交互lasso方法正確率最高，全特征對lasso方法正確率次之，但較lasso方法有很大提高，這說明肝病數據存在特征交互。這是因為全特征對lasso方法考慮了所有交互項，而交互lasso方法既考慮了特征交互又考慮了特征分層，從而大大降低了時間損耗。
表3
lasso方法在BLD數據集上10倍交叉驗證的實驗結果比較
Table3.
Results of four lasso models for BLD
3.4 本文方法與傳統方法的結果比較
圖 1給出了交互lasso與傳統分類方法10倍交叉驗證的對比，可知交互lasso在正確率上優于傳統分類方法。表 4給出了SVM、最近鄰、LDA、決策樹、普通lasso方法以及交互lasso方法的特異性、敏感性以及F測度，可以看出本文方法在三種評測規則中得分較高，證明了本文方法的良好性能。
圖1
本文所建立模型與其他經典方法正確率比較
Figure1.
Accuracy rate comparison between our model and other traditional methods
表4
基本方法以及lasso方法在BLD數據集上10倍交叉驗證的實驗結果比較(%)
Table4.
Comparison between results of our model and those of the traditional methods for BLD (%)
3.5 本文方法與其他文獻結果比較
作者將本文方法與其它文獻方法用于BLD數據集時的準確率進行比較，如圖 2所示，本文方法優于其他文獻結果。文獻[3]所用方法為神經網絡方法和模糊邏輯實例學習方法，文獻[15]采用基于單純函數的稀疏網格方法，文獻[16]采用了近鄰算法與SVM相結合的方法，文獻[17]采用了進化學習理論。
圖2
本文方法與其他文獻方法正確率比較
Figure2.
Accuracy rate comparison between our model and other methods shown in the papers
3.6 討論
由上述結果可以看出，本文方法有以下三方面的貢獻：第一，交互信息在肝病診斷中具有突出貢獻，交互特征的應用能夠提高診斷準確率。第二，分層約束條件的應用在很大程度上降低了時間損耗，模型學習及識別過程相對于不添加約束的方法可以省去約三分之一的時間。第三，lasso罰函數的應用規避了高維數據的特征選擇困難的問題。總之，上述實驗展現了交互特征在肝病診斷中的廣泛存在，證明交互lasso適合處理生物醫學數據。
4 結論
lasso具有特征線性加權和的特點，故可解釋性強，應該在生物醫學數據分類中廣泛應用。但是事實恰恰相反，原因之一是生物醫學數據普遍存在特征復雜交互問題，線性加權和并不能得到優異的分類結果，專家們紛紛采用復雜的非線性分類器。而非線性分類器的黑箱性質和復雜性，使它實踐應用困難，相反決策樹方法反而在實踐中應用廣泛。因此基于lasso和特征交互思想，作者提出將交互lasso分類器用于肝病分類，給出了模型定義、凸約束條件、凸優化方法和系數解。實驗結果證實了肝病特征存在交互，以及交互對于分類意義明顯。與基本機器學習方法和lasso方法相比，交互lasso的正確率、F測度等分類性能最好。新方法可以推廣到其它生物醫學數據分類問題，可以預測這種可解釋性強的交互線性模型在醫療診斷中會得到越來越多的認可和實踐驗證。接下來可進一步深入研究高階交互模型。
0 引言
肝病在我國是一種高發病，已經成為一種嚴重影響我國人民生命健康的疾病。我國政府對肝病的防治尤為重視，從七五計劃起一直將肝病的防治作為連續攻關課題。故此本文作者著力于肝病分析研究，以尋找有效的肝病檢測治療方案，實驗數據來自UCI學習數據集。
在過去的研究工作中，眾多學者選用1990年Forsyth收集的肝臟失調樣本(BUPA Liver Disorders, BLD)數據集作為研究對象，以探索肝病的病理特征。研究者探索方向大致分為兩類。第一類為機器學習算法[1-2]，例如實例推理法、決策樹、貝葉斯分類、人工神經網絡以及基于模糊邏輯的實例學習。Neshat等[3-4]將模糊方法和霍普菲爾德神經網絡(Hopfield neural network)方法應用于肝臟異常診斷。這些傳統方法在一定程度上提高了肝病診斷的準確性，但可解釋性都較差，且忽略了肝病特征交互。另一類方法采用優化理論，如Lin等[5]采用粒子群智能優化算法得到78.18%的正確率，Polat等[6]采用監督人工免疫系統提高診斷正確率。此兩種方法雖然提高了肝病診斷的準確率同時也具有一定的可解釋性，但復雜度較高，不適用于現代醫學應用。Chien等[7]提出了平滑套索方法，得到了68.82%的正確率，建立的模型解釋性較強，但所得分類結果不夠理想，并且忽略了特征交互對病癥的影響。
2012年Ramana等收集了印度肝病數據集(Indian Liver Patient Dataset，ILPD)，它更貼近這些年來肝病的發展現狀。Karthik等[8]采用人工神經網絡分類ILPD肝病樣本，方法容易過擬合，也忽略了疾病間的交互特征。Neshat等[9]則將粒子群優化算法與案例推理結合進行特征提取和肝病分類，此方法屬于優化算法，也沒有利用疾病間的特征交互信息。Ekong等[10]通過模糊聚類的方法來分析血清蛋白，并以臨床特征作為輔助以診斷肝病。Lin[11]綜合了多種分類器、回歸樹和案例推理技術，構建了一種智能診斷模型。Ekong與Lin將智能模型與肝病病理知識結合，提高了肝病的分類正確率，但模型復雜度較高且忽略了交互信息在醫學中的重要意義。而套索(least absolute shrinkage and selection operator，lasso)方法用于ILPD的處理則未見報道。
基于上述方法模型復雜度高以及對特征交互信息利用率低的問題，本文以簡單實用的logistic回歸模型為基礎，結合文獻[12]提出的交互lasso模型和凸優化算法，研究了數據驅動的罰分類方法，即改進的lasso方法用于肝病診斷。lasso方法中響應變量是預測變量線性加權和的表示，因此模型具有很強的可解釋性。而特征間的復雜交互對于疾病診斷具有重要意義，這推動我們將特征交互思想用于疾病診斷。作者首先采用分層lasso罰分類模型，并添加凸約束條件，對主效應特征和交互特征進行收縮和選擇處理，而后通過卡羅需-庫恩-塔克最優化條件(Karush-Kuhn-Tucker，KKT)與廣義梯度下降方法結合解決最優化問題，并求解模型參數。本文的創新在于使用了分層特征交互思想，將可解釋性強的lasso分類方法推廣用于肝病診斷。最后本文對兩組肝病數據進行分類實驗，并與lasso方法、全特征對lasso方法、支持向量機(support vector machine，SVM)、線性判別分析(linear discriminant analysis，LDA)等其他傳統分類方法進行比較。新方法用于兩組數據的實驗結果表明了特征交互lasso分類方法在肝病診斷中的有效性。
1 交互lasso模型
假設模型樣本輸出為Y，輸入X為P維特征x1, …，xj, …, xP，樣本在每維輸入特征間都有交互xjxk，則考慮交互的logistic回歸模型具有以下形式
| $\begin{gathered} \log {\text{it}}\left({P\left({Y=1\left| X \right|} \right)} \right)=\sum\limits_{j=0}^P {{\beta _j}{x_j}} + \hfill \\ \frac{1}{2}\sum\limits_{j \ne k} {{\Theta _{jk}}{x_j}{x_k} + \varepsilon }, \hfill \\ \end{gathered} $ |
其中x0為1，xj表示1階主效應特征，xjxk表示2階交互特征，1階主效應特征系數為β∈RP+1，2階交互特征系數為Θ∈RP×P，Θjj=0，ε∈N(0, σ2)。可以看出式(1)是主效應特征和交互特征的線性加權和，模型可解釋性強。
假設訓練樣本為(x1, y1), …, (xi, yi), …, (xN, yN)，1階主效應特征和2階交互特征融合后特征數目為P+P(P-1)/2，我們目的是從融合特征中通過罰分類方法進行特征子集選擇，選擇出的特征子集作為相應的預測特征。這就需要估計出模型系數中非零參數值，因為它們對應被選擇的子集。定義，則兩類的概率分別如下：
| ${p_i}=p\left({{{\text{x}}_{\text{i}}}} \right)=p\left({Y=1|{{\text{x}}_{\text{i}}}} \right)=\frac{1}{{1 + \exp \left({-\sum\limits_j {{\beta _j}{x_{ij}}}-\frac{1}{2}\sum\limits_{j \ne k} {{\Theta _{jk}}{x_{ij}}{x_{ik}}} } \right)}}$ |
| $-{p_i}=1-p\left({{{\text{x}}_{\text{i}}}} \right)=p\left({Y=2|{{\text{x}}_{\text{i}}}} \right)=\frac{1}{{1 + \exp \left({-\sum\limits_j {{\beta _j}{x_{ij}}}-\frac{1}{2}\sum\limits_{j \ne k} {{\Theta _{jk}}{x_{ij}}{x_{ik}}} } \right)}}$ |
通過極大似然估計擬合模型，使N次獨立觀測的似然函數或者對數似然函數最大。只考察大的主效應特征間的交互更為高效，因此模型添加約束‖Θj‖1≤｜βj｜，即
| $\ell\left({\beta, \Theta } \right)=\sum\limits_{i=1}^N {\left[{{y_i}\log {p_i} + \left({1-{y_i}} \right)} \right]\log \left({1-{p_i}} \right)} $ |
其中Θj是Θ的第j行。如果jk≠0，那么有‖j‖1＞0，因此有j≠0。新增約束條件導致式(4)不能凸優化求解，因此使用β+, β-∈RP+1兩個向量代替β的策略使式(4)成為凸函數，并添加lasso罰函數[13]，即
| $\begin{gathered} \ell\left({{\beta ^ \pm }, \Theta } \right)=-\sum\limits_{i=1}^n {\left[{{y_i}\log {p_i} + \left({1-{y_i}} \right)\log \left({1-{p_i}} \right)} \right]} \hfill \\ + \lambda {1^T}\left({{\beta ^ + } + {\beta ^-}} \right) + \frac{\lambda }{2}{\left\| \Theta \right\|_1} \hfill \\ s.t.\; \;\left. \begin{gathered} {\left\| {{\Theta _j}} \right\|_1} \leqslant \beta _j^ + + \beta _j^-\hfill \\ \beta _j^ + \geqslant 0, \beta _j^-\geqslant 0 \hfill \\ \end{gathered} \right\}\; \;{\text{for}}\; \;j{\text{=1, }} \cdots {\text{, }}P, \hfill \\ {\beta _j}=\beta _j^ +-\beta _j^-, {\left\| {{\beta _j}} \right\|_1}=\beta _j^ + + \beta _j^-\hfill \\ \end{gathered} $ |
其中懲罰參數λ≥0，控制著訓練數據的相對重要性和稀疏性(l1罰項)。我們的目的是求解使式(5)最小化時的主效應特征系數與交互特征系數(β±, Θ)。
2 KKT條件與廣義梯度下降求解模型
我們用廣義梯度下降法[14]求解式(5)問題。
引理：設無約束優化問題有以下形式Minimize g(φ)+h(φ)，其中g和h為凸函數，則利用廣義梯度下降法可轉化為如下形式，這里k是迭代次數，t是步長：
| $\begin{gathered} {{\hat \varphi }_k}=\arg \min \frac{1}{{2t}}{\left\| {\varphi-{{\tilde \varphi }_{k-1}}} \right\|^2} + h\left(\varphi \right)=\hfill \\ \arg \min \frac{1}{{2t}}{\left\| {\varphi-\left[{{{\tilde \varphi }_{k-1}}-t\nabla g\left({{{\tilde \varphi }_{k-1}}} \right)} \right]} \right\|^2} + h\left(\varphi \right) \hfill \\ \end{gathered} $ |
證明顯而易見。
定理1：式(5)中令，兩式相加的無約束條件的全局最小值解為
| $\left\{ \begin{gathered} \tilde \beta _j^{ \pm \left({k-1} \right)}=\tilde \beta _j^{ \pm \left({k-1} \right)} \pm t\left[{\sum\limits_{i=1}^n {{x_{ij}}\left({{y_i}-p_i^{k-1}} \right) + \lambda } } \right] \hfill \\ \tilde \Theta _{jk}^{\left({k-1} \right)}=\hat \Theta _{jk}^{\left({k-1} \right)} + \frac{t}{2}\sum\limits_{i=1}^n {{x_{ij}}\left({{y_i}-p_i^{k-1}} \right){x_{ik}}} \hfill \\ \end{gathered} \right.$ |
證明：
令，對e±求βj+的偏導數，得，同理，。。
對g(β±, Θ)求βj+的偏導數，得。同理，對g(β±, Θ)求βj-的偏導數，得。
對g(β±, Θ)求Θjk的偏導數，得。
綜上：。
根據引理中可得：
| $ \begin{align} & \tilde{\beta }{{_{j}^{+}}^{\left(k-1 \right)}}\text{=}\tilde{\beta }{{_{j}^{+}}^{\left(k-1 \right)}}\pm t\left[\sum\limits_{i=1}^{n}{{{x}_{ij}}}\left({{y}_{i}}-p_{i}^{k-1} \right)+\lambda \right] \\ & \tilde{\Theta }_{jk}^{\left(k-1 \right)}\text{=}\tilde{\Theta }_{jk}^{\left(k-1 \right)}+\frac{t}{2}\sum\limits_{i=1}^{n}{{{x}_{ij}}}\left({{y}_{i}}-p_{i}^{k-1} \right)+{{x}_{ik}} \\ \end{align} $ |
根據引理中式(6)，根據特征維數可分解為P個獨立可并行解決的子問題，式(5)可轉換為如下凸優化形式：
| $\begin{align} & \underset{{{\beta }^{\pm }}\in {{R}^{p}}, \Theta \in {{R}^{p-1}}}{\mathop{Minimze}}\, \frac{1}{2t}{{\left\| \beta _{j}^{+}-\tilde{\beta }_{j}^{+} \right\|}^{2}}+\frac{1}{2t}{{\left\| \beta _{j}^{-}-\tilde{\beta }_{j}^{-} \right\|}^{2}}+ \\ & \frac{1}{2t}\left\| {{\Theta }_{j}}-{{{\tilde{\Theta }}}_{j}} \right\|_{F}^{2}+\frac{\lambda }{2}{{\left\| {{\Theta }_{j}} \right\|}_{1}} \\ & s.t.\ \ \ \ {{\left\| {{\Theta }_{j}} \right\|}_{1}}\le \beta _{j}^{+}+\beta _{j}^{-}, \beta _{j}^{+}\ge 0, \beta _{j}^{-}\ge 0, j=1, 2, \cdots, P \\ \end{align}$ |
定理2：設和為式(8)的解，則
| $\tilde{\beta }_{j}^{\pm }={{\left[\tilde{\beta }_{j}^{\pm }+t{{{\hat{\alpha }}}_{j}} \right]}_{+}}, {{{\hat{\Theta }}}_{j}}=S\left({{{\hat{\Theta }}}_{j}}, t\left(\lambda /2+{{{\hat{\alpha }}}_{j}} \right) \right)$ |
證明：式(8)拉格朗日運算定義如下：
| $ \begin{align} & L\left(\hat{\beta }_{j}^{\pm }, {{{\hat{\Theta }}}_{j}}; {{{\hat{\alpha }}}_{j}}, {{{\hat{\gamma }}}_{j\pm }} \right)=\frac{1}{2t}{{\left(\hat{\beta }_{j}^{+}-\hat{\beta }_{j}^{+} \right)}^{2}}+\frac{1}{2t}{{\left(\hat{\beta }_{j}^{-}-\hat{\beta }_{j}^{-} \right)}^{2}}+ \\ & \frac{1}{2t}{{\left({{{\hat{\Theta }}}_{j}}-{{{\tilde{\Theta }}}_{j}} \right)}^{2}}+\left(\frac{\lambda }{2}+{{{\hat{\alpha }}}_{j}} \right){{\left\| {{{\hat{\Theta }}}_{j}} \right\|}_{1}}-\left(\hat{\gamma }_{j}^{+}+{{{\hat{\alpha }}}_{j}} \right)\hat{\beta }_{j}^{+}-\left(\hat{\gamma }_{j}^{-}+{{{\hat{\alpha }}}_{j}} \right)\hat{\beta }_{j}^{-} \\ \end{align} $ |
其中是分層約束對偶特征，是非負約束對偶特征。
有KKT條件如下[9]：平穩性條件為，平穩性條件為，其中uj是絕對值函數關于的次梯度，互補松弛條件為。
現定義，其中S(·)是lasso領域常用的函數[16]。最后的KKT條件只包含。觀察發現f關于是非增的，基于其為分段線性函數，可求解出。根據以及有，將代入β±的平穩性條件可知，由平穩性條件可得。
綜上，對λ值的一組下降序列，就可以得到β，Θ的路徑解。我們采用網格法，對于固定的λ，求解系數的過程是一個循環嵌套過程。梯度下降法求解交互lasso模型算法歸結如下：已知X, Y，根據定理1中式(7)求解，根據定理2中式(9)求解，直至收斂。
3 實驗數據和結果分析
本實驗采用的數據是加州大學歐文分校(University of California Irvine, UCI)機器學習數據庫的BLD數據集以及ILPD數據集。實驗數據來自 http://archive.ics.uci.edu/ml/。實驗通過R語言編程建立模型并實現文中算法。
3.1 數據說明
BLD數據集包含樣本345個，其中健康樣本200個，肝病樣本145個。每個樣本包括6維實數特征，前5維皆來自受試者血液檢測分析，分別為平均紅細胞容積、堿性磷酸酶、丙氨酸轉氨酶、天冬氨酸轉氨酶以及谷氨酰基轉肽酶含量，最后1維為受試者每日平均飲酒量。
ILPD數據集，由三個印度教授收集自印度安得拉邦的東北部，包括樣本583個，其中肝病患者416個以及非肝病受試者167個。583個受試者中包含441名男性和142名女性。每個樣本包括10維實數特征，分別為年齡、性別、總膽紅素、直接膽紅素、堿性磷酸酶、丙氨酸轉氨酶、天冬氨酸轉氨酶、血清總蛋白、血清白蛋白、白蛋白與球蛋白之比(簡稱白球比例)。
3.2 特征子集選擇的實驗結果和分析
本文對BLD數據集以及ILPD數據集進行100次10倍交叉驗證，實驗記錄了交互lasso分類方法的非0特征系數個數、平均錯誤率、標準差、CPU時間、所選取的λ值。
表 1、2給出了交互模型在兩組肝病數據上的系數解。行表示是非零主效應特征系數，第1列表示主效應系數，其它列表示交互系數。表 1、2中特征系數非0時，即表示此特征被選擇。表 1、2中交互系數絕對值越大說明特征交互程度越高。例如，表 1中黑體數字表示堿性磷酸酶與丙氨酸轉氨酶的交互特征被選擇，且交互系數為0.649 5；同樣地，其它非0的交互特征也被選擇。
表1
分層交互模型用于BLD數據所得系數解
Table1.
Solution of BLD dataset obtained using hierarchical interactive model
在BLD數據集實驗中，表 1是λ值為2.37時得到的模型系數解。從表中可以看出模型選擇了全部主效應特征以及部分交互項，主要包括以下4個交互項：①堿性磷酸酶與丙氨酸轉氨酶交互項，此項交互不但排除了堿性磷酸酶升高原因為骨病的可能性，而且提高了肝細胞病變診斷的可信度。②堿性磷酸酶與天冬氨酸轉氨酶交互項，此項作用與前項相似。③丙氨酸轉氨酶與天冬氨酸轉氨酶交互項，因天冬氨酸在心肌含量最高，依次為肝、骨骼肌、腎、胰，只有與其他特征結合才能斷定天冬氨酸轉氨酶升高是肝病引發的。④受試者飲酒量與幾乎其他各項均存在交互，可見飲酒對肝功能影響甚大。
在ILPD數據集實驗中，表 2是最優λ值為1.52時得到模型的主效應系數解和部分交互系數解。下面討論表 2中5項主要交互：①直接膽紅素與丙氨酸轉氨酶交互項，組成交互項的兩個特征均為肝病診斷的有效識別依據，而兩種特征交互，使此交互特征更突出。②總膽紅素與直接膽紅素交互項，總膽紅素升高但直接膽紅素不升高的原因并非肝細胞病變，故排除了許多其他病因。③丙氨酸與天冬氨酸交互項，與表 1中BLD實驗相同，因天冬氨酸存在于心肌、肝、骨骼肌、腎、胰等，在其與丙氨酸交互后，強化了此交互特征在模型中的重要性。④血清白蛋白與白球比例交互項，根據數據說明中的肝病檢測先驗知識可知，白蛋白降低且白球比例降低是典型正交互項，可提高醫生的診斷把握。⑤堿性磷酸酶與丙氨酸轉氨酶交互項，此項交互不但排除了堿性磷酸酶升高原因為骨病的可能性，并且提高了肝細胞病變診斷的可信度。
表2
分層交互模型用于ILPD數據所得系數解
Table2.
Solution of ILPD dataset obtained using hierarchical interactive model
3.3 本文方法與基本lasso方法的結果比較
表 3給出了幾種lasso方法的對比結果。從表中可以看出交互lasso方法正確率最高，全特征對lasso方法正確率次之，但較lasso方法有很大提高，這說明肝病數據存在特征交互。這是因為全特征對lasso方法考慮了所有交互項，而交互lasso方法既考慮了特征交互又考慮了特征分層，從而大大降低了時間損耗。
表3
lasso方法在BLD數據集上10倍交叉驗證的實驗結果比較
Table3.
Results of four lasso models for BLD
3.4 本文方法與傳統方法的結果比較
圖 1給出了交互lasso與傳統分類方法10倍交叉驗證的對比，可知交互lasso在正確率上優于傳統分類方法。表 4給出了SVM、最近鄰、LDA、決策樹、普通lasso方法以及交互lasso方法的特異性、敏感性以及F測度，可以看出本文方法在三種評測規則中得分較高，證明了本文方法的良好性能。
圖1
本文所建立模型與其他經典方法正確率比較
Figure1.
Accuracy rate comparison between our model and other traditional methods
表4
基本方法以及lasso方法在BLD數據集上10倍交叉驗證的實驗結果比較(%)
Table4.
Comparison between results of our model and those of the traditional methods for BLD (%)
3.5 本文方法與其他文獻結果比較
作者將本文方法與其它文獻方法用于BLD數據集時的準確率進行比較，如圖 2所示，本文方法優于其他文獻結果。文獻[3]所用方法為神經網絡方法和模糊邏輯實例學習方法，文獻[15]采用基于單純函數的稀疏網格方法，文獻[16]采用了近鄰算法與SVM相結合的方法，文獻[17]采用了進化學習理論。
圖2
本文方法與其他文獻方法正確率比較
Figure2.
Accuracy rate comparison between our model and other methods shown in the papers
3.6 討論
由上述結果可以看出，本文方法有以下三方面的貢獻：第一，交互信息在肝病診斷中具有突出貢獻，交互特征的應用能夠提高診斷準確率。第二，分層約束條件的應用在很大程度上降低了時間損耗，模型學習及識別過程相對于不添加約束的方法可以省去約三分之一的時間。第三，lasso罰函數的應用規避了高維數據的特征選擇困難的問題。總之，上述實驗展現了交互特征在肝病診斷中的廣泛存在，證明交互lasso適合處理生物醫學數據。
4 結論
lasso具有特征線性加權和的特點，故可解釋性強，應該在生物醫學數據分類中廣泛應用。但是事實恰恰相反，原因之一是生物醫學數據普遍存在特征復雜交互問題，線性加權和并不能得到優異的分類結果，專家們紛紛采用復雜的非線性分類器。而非線性分類器的黑箱性質和復雜性，使它實踐應用困難，相反決策樹方法反而在實踐中應用廣泛。因此基于lasso和特征交互思想，作者提出將交互lasso分類器用于肝病分類，給出了模型定義、凸約束條件、凸優化方法和系數解。實驗結果證實了肝病特征存在交互，以及交互對于分類意義明顯。與基本機器學習方法和lasso方法相比，交互lasso的正確率、F測度等分類性能最好。新方法可以推廣到其它生物醫學數據分類問題，可以預測這種可解釋性強的交互線性模型在醫療診斷中會得到越來越多的認可和實踐驗證。接下來可進一步深入研究高階交互模型。
