紅斑鱗狀皮膚病是6種常見皮膚科疾病的統稱,其診斷一直是皮膚科的難題。本文提出了基于名詞性變量虛擬編碼的彈性網罰多值logistic回歸方法。該方法從皮膚病數據的變量屬性出發,首先將名詞性變量進行虛擬編碼,避免了標稱值直接計算帶來的不合理性;繼而通過帶彈性網罰的多值logistic回歸模型擬合特征與疾病分類間的關系;最后通過坐標下降法求得模型的參數估計。實驗中采用10折交叉驗證方法并達到了98.34%±0.0027%的診斷正確率,與其他方法相比,本文方法正確率相當,且步驟簡單,穩定性很強。
引用本文： 王金甲, 李慧. 基于虛擬編碼和彈性網多值Logistic回歸的紅斑鱗狀皮膚病診斷方法. 生物醫學工程學雜志, 2015, 32(4): 757-762. doi: 10.7507/1001-5515.20150138 復制
引言
紅斑鱗狀皮膚病是一組以慢性皮膚病為主的常見皮膚科疾病，包括牛皮癬、脂溢性皮炎、扁平苔蘚、玫瑰糠疹、克羅尼克皮炎和毛發紅糠疹。患者不僅面臨身體上的痛苦，更受到精神上的困擾。正確的診斷是治療的前提，然而其診斷卻存在兩方面的困境：一方面，這組疾病在臨床上都表現出紅斑、脫屑等特點，且在發病初期，某一種疾病的獨有特征并不會表現出來，各疾病癥狀之間差異很小；另一方面，對皮膚樣本的病理組織學數據分析發現，這組疾病也有許多特征是相似的，增加了診斷的難度。
國內外眾多學者對該組疾病進行了深入研究，并取得了顯著的成績。最初，該數據集的提供者Guvenir等[1-2]提出特征間隔投票(voting feature intervals,VFI)分類算法取得了96.2%的正確率，并通過遺傳算法的特征加權策略使正確率達到99.2%。之后，部分學者從分類器設計的角度提出了多種解決方案，包括Ubeyli等[3-5]的自適應模糊神經網絡系統、帶有誤差修正輸出碼的多類支持向量機(support vector machine,SVM)、聯合神經網絡模型；Luukka等[6-7]的模糊相似分類器和基于Yu’s范數的相似性分類器；Polat等[8-9]的基于模糊加權預處理和最近鄰加權預處理的決策樹分類器、基于C4.5決策樹和一對多策略的混合智能方法，都取得了較高的分類準確度。Ubeyli等[10]還提出了基于k-Means聚類的自動檢測方法。近年來，更多學者從特征選擇與分類器結合的角度提出了紅斑鱗狀皮膚病診斷的新方法，包括Nanni[11]的基于拉格朗日支持向量機(Lagrange support vector machine,LSVM)的最優子空間方法；Liu等[12]的動態互信息特征選擇結合樸素貝葉斯、最近鄰、C4.5決策樹和重復增量修枝(repeated incremental pruning to produce error reduction,RIPPER)分類器的方法；Karabatak等[13]的基于關聯規則和神經網絡的特征選擇算法；Xie等[14]的基于改進F-score的順序前向浮動特征選擇方法；Ozcift等[15-16]的基于旋轉森林系統的魯棒特征選擇策略、遺傳算法與貝葉斯網絡結合的特征選擇方法；Abdi等[17]的關聯規則特征選擇結合基于粒子群參數優化的SVM分類器方法，使該疾病的診斷正確率得到了進一步提高。
本文提出了基于名詞性變量虛擬編碼的彈性網罰多值logistic回歸方法。首先將皮膚病數據中的名詞性變量進行虛擬編碼，繼而通過帶彈性網罰的多值logistic回歸模型擬合特征與疾病分類間的關系，采用10折交叉驗證方法并計算診斷正確率，以期為紅斑鱗狀皮膚病的診斷提供一種新的途徑，也對名詞性變量的處理方法進行有益探索。
1 數據說明
本文實驗采用的數據是UCI數據庫中的dermatology數據集，該組數據共包含366個樣本(其中8個年齡數據缺失，實驗中將其刪除)，每個樣本特征數為34維，12維是臨床醫學特征，通過對患者的觀察和詢問獲得，22維是病理組織學特征，通過對患者皮膚樣本的顯微觀察獲得。其中，患者年齡是整數值；家族病史取值為1表示家族中有人患過上述任意一種疾病，否則取值為0；其他的特征(臨床的和病理組織學的)都給定了0到3的取值范圍，0表示該特征不存在，3表示可能的最大程度，1、2表示相對的中間值。
2 名詞性變量虛擬編碼
樣本的特征是通過變量描述的，變量一般根據變量屬性的定義進行取值，因此變量值所能進行的計算也不盡相同。根據變量屬性，可將其劃分為名詞性的(定性的)和數值性的(定量的)[18]，具體區分如表 1所示。
表1
變量屬性的劃分
Table1.
Division of the properties of variable
可見，在皮膚病數據集的全部34維特征中，只有年齡特征是數值變量，其他33維變量都是名詞性的標稱值，其取值僅僅表示了患者是否出現該特征或其表現出的程度，并不真正具有數值上的意義。顯然，對名詞性變量和數值性變量不加區分地直接進行計算既不科學也不符合數據代表的實際意義。
針對上述情況，我們采用虛擬編碼的方法來處理名詞性變量。所謂虛擬編碼就是采用類似編碼的方式將不同取值的名詞性變量轉換到不同的維度，避免各個原始取值之間直接進行計算，該過程相當于映射f(·)，將名詞性變量轉換為矢量變量。假設名詞性變量的取值集合為{u∈Z:1≤u≤H}，虛擬編碼是將每一個取值u映射為一個矢量x∈RH，其中第u維賦值為1，其余維度賦值為0。以取值為1、2、3的特征為例：
f(3)=(0,0,1);f(2)=(0,1,0);f(1)=(1,0,0)
對皮膚病數據進行虛擬編碼，原始特征的33維標稱值變量轉換成了129維(第11維轉換為2維；第13維轉換為3維，其余轉換為4維)。第34維是年齡特征，實驗中我們采用了不同的處理方式并對實驗結果進行了比較。
將名詞性變量的不同取值轉化到相應維度上，雖然增加了樣本的總體維數，且使樣本變得稀疏，但避免了直接對名詞性變量計算帶來的不合理性，后續實驗結果證明，該方法顯著提高了樣本分類的準確率，是處理名詞性變量的有效途徑。
3 彈性網多值分類器
經過虛擬編碼后的樣本數據是高維稀疏的，傳統的分類器在處理這樣的數據時效果并不理想，因此，我們選擇采用多值logistic回歸模型擬合特征與分類間的關系，并加入彈性網罰，同時實現了特征選擇和樣本分類。
3.1 多值logistic回歸模型
多值logistic回歸模型是二值logistic回歸模型的推廣，是常用的廣義線性模型。定義g表示響應變量，取值g∈{1,2,…,K}，預測變量xi∈Rp，有N對觀測(xi,gi),i=1,…,N。
可以建立模型：
| $\Pr \left( {g = \left. \ell \right|x} \right) = \frac{{{{\rm{e}}^{{\beta _{0\ell }} + {x^T}{\beta _\ell }}}}}{{\sum\limits_{k = 1}^K {{{\rm{e}}^{{\beta _{0k}} + {x^T}{\beta _k}}}} }}\ell = 1, \cdots ,K$ |
其中β是p維的系數向量，β0是截矩，選擇g=K作為參照組，可以得到K-1個logistic模型：
| $ln = \frac{{\Pr \left( {g = \left. \ell \right|x} \right)}}{{\Pr \left( {g = \left. K \right|x} \right)}} = {\beta _{0\ell }} + {x^T}{\beta _\ell }\ell = 1, \cdots ,K - 1$ |
這里，有聯系函數gK(x)=β0K+xTβK=0。
用極大似然估計擬合上述模型。實驗數據集可以視為N個相互獨立的觀測值，其聯合概率分布為N個邊際概率的乘積，即N個觀測的似然函數為
| $\ell \left( {{\beta _0},\beta } \right) = \prod\limits_{i = 1}^N {\prod\limits_{\ell = 1}^K {{p_\ell }{{\left( {{x_i}} \right)}^{{y_{i\ell }}}}} } ,$ |
其中條件概率p(xi)=Pr(gi=｜xi)，對一個特定觀測。
現在的問題是求出最有可能得到現有觀測樣本集合的回歸模型，也就是使上述似然函數取值最大的估計參數。為便于計算，取對數似然函數，則有：
| $L\left( {{\beta _0},\beta } \right) = \frac{1}{N}\ln \left[ {\ell \left( {{\beta _0},\beta } \right)} \right] = \frac{1}{N}\sum\limits_{i = 1}^N {\sum\limits_{\ell = 1}^K {{y_{i\ell }}\ln {p_\ell }\left( {{x_i}} \right)} } $ |
注意到：，則有：
| $L\left( {{\beta _0},\beta } \right) = \frac{1}{N}\sum\limits_{i = 1}^N {\left\{ {\sum\limits_{\ell = 1}^K {{y_{i\ell }}\left[ {\ln {p_\ell }\left( {{x_i}} \right) - \ln {p_K}\left( {{x_i}} \right) + \ln {p_K}\left( {{x_i}} \right)} \right]} } \right\} = } \frac{1}{N}\sum\limits_{i = 1}^N {\left[ {\sum\limits_{\ell = 1}^K {{y_{i\ell }}\left( {{\beta _{0\ell }} + x_i^T{\beta _\ell }} \right) - \ln \left( {\sum\limits_{k = 1}^K {{e^{{\beta _{0k}} + x_i^T{\beta _k}}}} } \right)} } \right]} $ |
由于對數似然函數是似然函數的單調函數，所以問題等價于求對數似然函數的極大值，本文采用的方法是加入彈性網罰并通過坐標下降法求解。
3.2 彈性網
以線性回歸模型為例，假設預測變量X∈Rp，響應變量Y∈R，線性回歸模型為E(Y｜X=x)=β0+xTβ，有N對觀測值(xi,yi)，且xij是標準化的，則彈性網問題可以描述如下[19]：
| $\mathop {\min }\limits_{\left( {{\beta _0},\beta } \right) \in {R^{p + 1}}} {R_\lambda }\left( {\beta 0,\beta } \right) = \mathop {\min }\limits_{\left( {{\beta _0},\beta } \right) \in {R^{p + 1}}} \left[ {\frac{1}{{2N}}\sum\limits_{i = 1}^N {{{\left( {{y_i} - {\beta _0} - x_i^T\beta } \right)}^2} + \lambda {P_\alpha }\left( \beta \right)} } \right],$ |
其中
| ${P_\alpha }\left( \beta \right) = \left( {1 - \alpha } \right)\frac{1}{2}\left\| \beta \right\|_{{\ell _2}}^2 + \alpha {\left\| \beta \right\|_{{\ell _1}}}$ |
可見，彈性網方法是在最小二乘估計的基礎上加入了彈性網罰函數，是對響應變量的有偏估計。彈性網罰是脊回歸罰和最小絕對值收縮和選擇算子(least absolute shrinkage and selection operator,lasso)罰的凸結合，當α=0，彈性網罰變為脊回歸罰；當α=1，彈性網罰變為lasso罰。當α從0增加到1，對于給定的λ，式(6)解的稀疏值(系數等于0的數)從0單調增加到lasso解的稀疏值。
彈性網罰是脊回歸罰和lasso罰的加權平衡，因此它能夠綜合兩種方法的優點，既能收縮系數，提高模型穩定性，又能進行變量選擇，建立可解釋的模型，進而提高預測精度。
對于上述式(5)所示的求對數似然函數極大值的問題，我們先做一項轉換，常用的牛頓解法等價于迭代再加權最小二乘[20]，假設參數當前估計是{,}1K，每次只允許一類{,}變化，在{,}點Taylor展開，得到對數似然函數的二次逼近：
| ${L_Q}\left( {{\beta _{0\ell }},{\beta _\ell }} \right) = - \frac{1}{{2N}}\sum\limits_{i = 1}^N {{\omega _{i\ell }}{{\left( {{z_{i\ell }} - {\beta _{0\ell }} - x_i^T{\beta _\ell }} \right)}^2} + C{{\left( {\left\{ {{{\tilde \beta }_{0k}},{{\tilde \beta }_k}} \right\}_1^K} \right)}^2},} $ |
其中，，是根據當前參數計算所得的值，針對固定的{,}1K，最后一項是常數。
那么，加入彈性網罰之后，需要求解的問題可以描述為
| $\mathop {\min }\limits_{\left( {{\beta _{0\ell }},{\beta _\ell }} \right) \in {R^{p + 1}}} \left\{ { - {L_Q}\left( {{\beta _{0\ell }},{\beta _\ell }} \right) + \lambda {P_\alpha }\left( {{\beta _\ell }} \right)} \right\}$ |
3.3 坐標下降法
與著名的最小角回歸(least angle regression,LARS)算法相比，坐標下降法[21]在求解lasso等一類凸優化問題時更簡單、有效，且容易推廣到其他相關問題，如彈性網等。其基本思想是每次只優化一維變量，這樣優化問題可以很快完成，且相關方程可以在變量循環中更新；而通常情況下，最小化變量在眾多參數中不隨循環而改變，因此整個迭代過程將很快完成。需要注意的是，該方法的前提條件是多個預測變量之間互不相關，這樣才能將多變量問題理解為多個單變量問題的組合。
對式(9)的問題，記：
| $L = - {L_Q}\left( {{\beta _{0\ell }},{\beta _\ell }} \right) + \lambda {P_\alpha }\left( {{\beta _\ell }} \right) = \frac{1}{{2N}}\sum\limits_{i = 1}^N {{\omega _{i\ell }}\left( {{{\left( {{z_{i\ell }} - {{\tilde g}_\ell }{{\left( {{x_i}} \right)}^{\left( j \right)}} - {x_{ij}}{\beta _{j\ell }}} \right)}^2}} \right. + C{{\left( {\left\{ {{{\tilde \beta }_{0k}},\tilde \beta } \right\}_1^K} \right)}^2} + \lambda {P_\alpha }\left( {{\beta _\ell }} \right)} ,$ |
其中，是除去xij一項的聯系函數。
針對當前的類別，每次只優化系數β的一維，其他維視為常數。假設當前估計是{,}1K，對系數β的第j維βj求導，當βj>0時：
| $\begin{array}{l} \frac{{\partial L}}{{\partial {\beta _{j\ell }}}} = \frac{1}{N}\sum\limits_{i = 1}^N {{\omega _{i\ell }}\left( {{z_{i\ell }} - {{\tilde g}_\ell }{{\left( {{x_i}} \right)}^{\left( j \right)}} - {x_{ij}}{\beta _{j\ell }}} \right)\left( { - {x_{ij}}} \right) + \lambda \left( {1 - \alpha } \right){\beta _{j\ell }} + \lambda \alpha } \\ = - \frac{1}{N}\sum\limits_{i = 1}^N {{\omega _{i\ell }}{x_{ij}}\left( {{z_{i\ell }} - {{\tilde g}_\ell }{{\left( {{x_i}} \right)}^{\left( j \right)}}} \right) + \lambda \alpha + \frac{1}{N}\sum\limits_{i = 1}^N {{\omega _{i\ell }}x_{ij}^2} {\beta _{j\ell }} + \lambda \left( {1 - \alpha } \right){\beta _{j\ell }}} \end{array}$ |
令上式等于0，可以求得βj的坐標更新形式：
| ${\beta _{j\ell }} = \frac{{\frac{1}{N}\sum\limits_{i = 1}^N {{\omega _{i\ell }}{x_{ij}}\left( {{z_{i\ell }} - {{\tilde g}_\ell }{{\left( {{x_i}} \right)}^{\left( j \right)}}} \right) - \lambda \alpha } }}{{\frac{1}{N}\sum\limits_{i = 1}^N {{\omega _{i\ell }}x_{ij}^2 + \lambda \left( {1 - \alpha } \right)} }}\left( {{\beta _{j\ell }} > 0} \right)$ |
當βj<0時，同樣可以得到類似的表達式。其他情況下βj=0，即有軟閾值形式：
| ${\beta _{j\ell }} = \frac{{S\left( {\frac{1}{N}\sum\limits_{i = 1}^N {{\omega _{i\ell }}{x_{ij}}\left( {{z_{i\ell }} - {{\tilde g}_\ell }{{\left( {{x_i}} \right)}^{\left( j \right)}}} \right),\lambda \alpha } } \right)}}{{\frac{1}{N}\sum\limits_{i = 1}^N {{\omega _{i\ell }}x_{ij}^2 + \lambda \left( {1 - \alpha } \right)} }}\left( {{\beta _{j\ell }} > 0} \right)$ |
總的來說，對λ值的序列，可以得到β的一組路徑解。對于每一個λ，根據當前觀測值，建立多值logistic回歸模型，即求解系數β的過程是一系列循環嵌套的迭代求解過程，每個循環嵌套于上一個循環中：
循環1：循環∈{1,2,…,K,1,2,…}，直到β收斂；
循環2：使用當前參數{,}1K更新LQ(β0,β)二次逼近；
循環3：通過坐標下降算法求解式(9)的罰加權最小二乘問題，得到βj。
4 實驗結果及分析
實驗中，我們采用帶彈性網罰的多值logistic模型進行學習和分類，為提高結果的科學性，采用10折交叉驗證方法，并重復100次得到測試錯誤率的均值和標準差。
本文方法中影響診斷結果的參數主要有兩個：彈性網罰系數λ和彈性網罰中的正則化參數α。λ起到控制模型自由度的作用，實驗中給定100個λ值的序列，將訓練數據送入彈性網多值分類器，將得到100種參數估計結果，選擇最低分類錯誤率對應的參數估計作為訓練結果，因此系數λ的優化在程序中自動完成。α控制參數估計的稀疏值，實驗中對α的優化過程是通過比較實現的，其不同取值的實驗結果如表 2所示，該實驗結果是利用了名詞性變量的虛擬編碼和歸一化的年齡特征，樣本維數為130維。可以看出，α取值為0.01時能取得更低的錯誤率，且標準差較小。表中列出的時間是一次10折交叉驗證的時間。
表2
不同α值實驗結果比較
Table2.
Experimental results with different values of α
在一次10折交叉驗證中，每一折訓練過程所用數據均有差異，因此模型最終所選變量并不完全相同，但與傳統的特征子集選擇相比，該方法選擇的特征子集和模型分類性能的穩定性更高。以α取值為0.01時的100次特征選擇結果為例，我們把與第6類對應的100個特征子集(對應β，=6)進行統計(按被選擇的次數排列特征，橫軸標度僅代表特征數，不代表特征標號)，如圖 1所示。
圖1
100次特征選擇結果統計
Figure1.
Statistic of 100 feature selection results
由圖 1可見，有70維左右的特征幾乎每次都被選擇，而其余60維左右的特征從未被選擇，說明彈性網罰logistic回歸方法在子集選擇上具有很高的穩定性。
實驗中采用了多種方式處理名詞性變量和數值性變量。對于名詞性變量，對比了原始特征和虛擬編碼特征，對于數值性變量年齡，采用了4種處理方法，分別是：歸一化的年齡變量；離散化的年齡變量；虛擬編碼的年齡變量；去除年齡特征。其對應關系如表 3所示，不同處理方法的實驗結果對比如圖 2所示。
表3
年齡特征不同處理方法間的對應關系
Table3.
Correspondence of the feature of age based on different handing methods
圖 2中：A組采用原始特征，特征維數34；B組采用歸一化的原始特征，特征維數34；C組采用名詞性變量的虛擬編碼，年齡變量做歸一化處理，特征維數130；D組僅采用33維名詞性變量的虛擬編碼，去掉了年齡特征，特征維數129；E組采用名詞性變量的虛擬編碼和離散的年齡特征，特征維數130；F組采用所有34維變量的虛擬編碼，特征維數133。
圖2
特征不同處理方法的實驗結果比較
A,B,C,D,E,F分別表示不同的特征組合
Figure2. Comparison of results based on different handing methods of the featuresA,B,C,D,E and F represent different combinations,respectively
可以看出，與原始特征的實驗結果相比較，采用名詞性變量的虛擬編碼能顯著提高結果的預測精度，說明該方法是一種有效的名詞性變量處理方法。另外，E組在α取值為0.01時能取得更高的正確率，為98.34%，且標準差較小，說明在各種變量處理方法中，名詞性變量虛擬編碼和數值性變量離散化是一種可行的途徑。
我們采用上述E組實驗的特征處理方法結合多種分類器進行了訓練和測試，同樣采用100次的10折交叉驗證，實驗結果對比如圖 3所示。
圖3
多種分類器實驗結果比較
Figure3.
Comparison of experimental results based on different classifiers
實驗中我們對比了最近鄰(KNN-1)、3-近鄰(KNN-3)、K-近鄰(KNN)、SVM、貝葉斯分類(BC)、Logistic、Fisher、感知器(Perceptron)、偏最小二乘(PLS)等多種分類模型，各種模型對數據的處理方式相同，均將名詞性特征做虛擬編碼處理，數值性特征做離散化處理。可以看出，本文提出的方法對皮膚病數據集的分析起到了很好的效果，取得了較高的正確率(98.34%±0.002 7%)，且算法簡單穩定，易于編程實現，預測結果能夠為疾病的診斷提供較為可靠的參考。
5 結論
皮膚病數據集是一個同時包含了數值性變量和名詞性變量的數據集，傳統的處理方法沒有對兩種類型的變量區別對待。本文提出的名詞性變量虛擬編碼方法從特征取值的實際意義入手，將不同取值轉換到不同的維度，實驗中對比了采用原始數據與虛擬編碼的不同結果，證明該方法能夠顯著提高樣本分類的準確率，是一種處理名詞性變量的有效途徑。另外，針對虛擬編碼后的高維稀疏數據，采用帶彈性網罰的多值logistic模型擬合特征與疾病分類間的關系，能夠較好地反映數據的內在聯系，得到預測精度較高的模型，且加入彈性網罰起到了變量選擇的作用，省去了特征選擇的步驟，使整體處理過程簡單明了，易于實現。
本文提出的方法在實驗中取得了98.34%的平均正確率，100次實驗結果的標準差僅為0.002 7%，也就是說，對全部358個樣本的分類平均出錯不會超過6個，且算法步驟簡單，穩定性很強，可以考慮推廣到相似的分類問題中。同時，對名詞性變量和數值性變量不同處理方法的探索也具有較強的實際意義，值得進一步深入研究。
引言
紅斑鱗狀皮膚病是一組以慢性皮膚病為主的常見皮膚科疾病，包括牛皮癬、脂溢性皮炎、扁平苔蘚、玫瑰糠疹、克羅尼克皮炎和毛發紅糠疹。患者不僅面臨身體上的痛苦，更受到精神上的困擾。正確的診斷是治療的前提，然而其診斷卻存在兩方面的困境：一方面，這組疾病在臨床上都表現出紅斑、脫屑等特點，且在發病初期，某一種疾病的獨有特征并不會表現出來，各疾病癥狀之間差異很小；另一方面，對皮膚樣本的病理組織學數據分析發現，這組疾病也有許多特征是相似的，增加了診斷的難度。
國內外眾多學者對該組疾病進行了深入研究，并取得了顯著的成績。最初，該數據集的提供者Guvenir等[1-2]提出特征間隔投票(voting feature intervals,VFI)分類算法取得了96.2%的正確率，并通過遺傳算法的特征加權策略使正確率達到99.2%。之后，部分學者從分類器設計的角度提出了多種解決方案，包括Ubeyli等[3-5]的自適應模糊神經網絡系統、帶有誤差修正輸出碼的多類支持向量機(support vector machine,SVM)、聯合神經網絡模型；Luukka等[6-7]的模糊相似分類器和基于Yu’s范數的相似性分類器；Polat等[8-9]的基于模糊加權預處理和最近鄰加權預處理的決策樹分類器、基于C4.5決策樹和一對多策略的混合智能方法，都取得了較高的分類準確度。Ubeyli等[10]還提出了基于k-Means聚類的自動檢測方法。近年來，更多學者從特征選擇與分類器結合的角度提出了紅斑鱗狀皮膚病診斷的新方法，包括Nanni[11]的基于拉格朗日支持向量機(Lagrange support vector machine,LSVM)的最優子空間方法；Liu等[12]的動態互信息特征選擇結合樸素貝葉斯、最近鄰、C4.5決策樹和重復增量修枝(repeated incremental pruning to produce error reduction,RIPPER)分類器的方法；Karabatak等[13]的基于關聯規則和神經網絡的特征選擇算法；Xie等[14]的基于改進F-score的順序前向浮動特征選擇方法；Ozcift等[15-16]的基于旋轉森林系統的魯棒特征選擇策略、遺傳算法與貝葉斯網絡結合的特征選擇方法；Abdi等[17]的關聯規則特征選擇結合基于粒子群參數優化的SVM分類器方法，使該疾病的診斷正確率得到了進一步提高。
本文提出了基于名詞性變量虛擬編碼的彈性網罰多值logistic回歸方法。首先將皮膚病數據中的名詞性變量進行虛擬編碼，繼而通過帶彈性網罰的多值logistic回歸模型擬合特征與疾病分類間的關系，采用10折交叉驗證方法并計算診斷正確率，以期為紅斑鱗狀皮膚病的診斷提供一種新的途徑，也對名詞性變量的處理方法進行有益探索。
1 數據說明
本文實驗采用的數據是UCI數據庫中的dermatology數據集，該組數據共包含366個樣本(其中8個年齡數據缺失，實驗中將其刪除)，每個樣本特征數為34維，12維是臨床醫學特征，通過對患者的觀察和詢問獲得，22維是病理組織學特征，通過對患者皮膚樣本的顯微觀察獲得。其中，患者年齡是整數值；家族病史取值為1表示家族中有人患過上述任意一種疾病，否則取值為0；其他的特征(臨床的和病理組織學的)都給定了0到3的取值范圍，0表示該特征不存在，3表示可能的最大程度，1、2表示相對的中間值。
2 名詞性變量虛擬編碼
樣本的特征是通過變量描述的，變量一般根據變量屬性的定義進行取值，因此變量值所能進行的計算也不盡相同。根據變量屬性，可將其劃分為名詞性的(定性的)和數值性的(定量的)[18]，具體區分如表 1所示。
表1
變量屬性的劃分
Table1.
Division of the properties of variable
可見，在皮膚病數據集的全部34維特征中，只有年齡特征是數值變量，其他33維變量都是名詞性的標稱值，其取值僅僅表示了患者是否出現該特征或其表現出的程度，并不真正具有數值上的意義。顯然，對名詞性變量和數值性變量不加區分地直接進行計算既不科學也不符合數據代表的實際意義。
針對上述情況，我們采用虛擬編碼的方法來處理名詞性變量。所謂虛擬編碼就是采用類似編碼的方式將不同取值的名詞性變量轉換到不同的維度，避免各個原始取值之間直接進行計算，該過程相當于映射f(·)，將名詞性變量轉換為矢量變量。假設名詞性變量的取值集合為{u∈Z:1≤u≤H}，虛擬編碼是將每一個取值u映射為一個矢量x∈RH，其中第u維賦值為1，其余維度賦值為0。以取值為1、2、3的特征為例：
f(3)=(0,0,1);f(2)=(0,1,0);f(1)=(1,0,0)
對皮膚病數據進行虛擬編碼，原始特征的33維標稱值變量轉換成了129維(第11維轉換為2維；第13維轉換為3維，其余轉換為4維)。第34維是年齡特征，實驗中我們采用了不同的處理方式并對實驗結果進行了比較。
將名詞性變量的不同取值轉化到相應維度上，雖然增加了樣本的總體維數，且使樣本變得稀疏，但避免了直接對名詞性變量計算帶來的不合理性，后續實驗結果證明，該方法顯著提高了樣本分類的準確率，是處理名詞性變量的有效途徑。
3 彈性網多值分類器
經過虛擬編碼后的樣本數據是高維稀疏的，傳統的分類器在處理這樣的數據時效果并不理想，因此，我們選擇采用多值logistic回歸模型擬合特征與分類間的關系，并加入彈性網罰，同時實現了特征選擇和樣本分類。
3.1 多值logistic回歸模型
多值logistic回歸模型是二值logistic回歸模型的推廣，是常用的廣義線性模型。定義g表示響應變量，取值g∈{1,2,…,K}，預測變量xi∈Rp，有N對觀測(xi,gi),i=1,…,N。
可以建立模型：
| $\Pr \left( {g = \left. \ell \right|x} \right) = \frac{{{{\rm{e}}^{{\beta _{0\ell }} + {x^T}{\beta _\ell }}}}}{{\sum\limits_{k = 1}^K {{{\rm{e}}^{{\beta _{0k}} + {x^T}{\beta _k}}}} }}\ell = 1, \cdots ,K$ |
其中β是p維的系數向量，β0是截矩，選擇g=K作為參照組，可以得到K-1個logistic模型：
| $ln = \frac{{\Pr \left( {g = \left. \ell \right|x} \right)}}{{\Pr \left( {g = \left. K \right|x} \right)}} = {\beta _{0\ell }} + {x^T}{\beta _\ell }\ell = 1, \cdots ,K - 1$ |
這里，有聯系函數gK(x)=β0K+xTβK=0。
用極大似然估計擬合上述模型。實驗數據集可以視為N個相互獨立的觀測值，其聯合概率分布為N個邊際概率的乘積，即N個觀測的似然函數為
| $\ell \left( {{\beta _0},\beta } \right) = \prod\limits_{i = 1}^N {\prod\limits_{\ell = 1}^K {{p_\ell }{{\left( {{x_i}} \right)}^{{y_{i\ell }}}}} } ,$ |
其中條件概率p(xi)=Pr(gi=｜xi)，對一個特定觀測。
現在的問題是求出最有可能得到現有觀測樣本集合的回歸模型，也就是使上述似然函數取值最大的估計參數。為便于計算，取對數似然函數，則有：
| $L\left( {{\beta _0},\beta } \right) = \frac{1}{N}\ln \left[ {\ell \left( {{\beta _0},\beta } \right)} \right] = \frac{1}{N}\sum\limits_{i = 1}^N {\sum\limits_{\ell = 1}^K {{y_{i\ell }}\ln {p_\ell }\left( {{x_i}} \right)} } $ |
注意到：，則有：
| $L\left( {{\beta _0},\beta } \right) = \frac{1}{N}\sum\limits_{i = 1}^N {\left\{ {\sum\limits_{\ell = 1}^K {{y_{i\ell }}\left[ {\ln {p_\ell }\left( {{x_i}} \right) - \ln {p_K}\left( {{x_i}} \right) + \ln {p_K}\left( {{x_i}} \right)} \right]} } \right\} = } \frac{1}{N}\sum\limits_{i = 1}^N {\left[ {\sum\limits_{\ell = 1}^K {{y_{i\ell }}\left( {{\beta _{0\ell }} + x_i^T{\beta _\ell }} \right) - \ln \left( {\sum\limits_{k = 1}^K {{e^{{\beta _{0k}} + x_i^T{\beta _k}}}} } \right)} } \right]} $ |
由于對數似然函數是似然函數的單調函數，所以問題等價于求對數似然函數的極大值，本文采用的方法是加入彈性網罰并通過坐標下降法求解。
3.2 彈性網
以線性回歸模型為例，假設預測變量X∈Rp，響應變量Y∈R，線性回歸模型為E(Y｜X=x)=β0+xTβ，有N對觀測值(xi,yi)，且xij是標準化的，則彈性網問題可以描述如下[19]：
| $\mathop {\min }\limits_{\left( {{\beta _0},\beta } \right) \in {R^{p + 1}}} {R_\lambda }\left( {\beta 0,\beta } \right) = \mathop {\min }\limits_{\left( {{\beta _0},\beta } \right) \in {R^{p + 1}}} \left[ {\frac{1}{{2N}}\sum\limits_{i = 1}^N {{{\left( {{y_i} - {\beta _0} - x_i^T\beta } \right)}^2} + \lambda {P_\alpha }\left( \beta \right)} } \right],$ |
其中
| ${P_\alpha }\left( \beta \right) = \left( {1 - \alpha } \right)\frac{1}{2}\left\| \beta \right\|_{{\ell _2}}^2 + \alpha {\left\| \beta \right\|_{{\ell _1}}}$ |
可見，彈性網方法是在最小二乘估計的基礎上加入了彈性網罰函數，是對響應變量的有偏估計。彈性網罰是脊回歸罰和最小絕對值收縮和選擇算子(least absolute shrinkage and selection operator,lasso)罰的凸結合，當α=0，彈性網罰變為脊回歸罰；當α=1，彈性網罰變為lasso罰。當α從0增加到1，對于給定的λ，式(6)解的稀疏值(系數等于0的數)從0單調增加到lasso解的稀疏值。
彈性網罰是脊回歸罰和lasso罰的加權平衡，因此它能夠綜合兩種方法的優點，既能收縮系數，提高模型穩定性，又能進行變量選擇，建立可解釋的模型，進而提高預測精度。
對于上述式(5)所示的求對數似然函數極大值的問題，我們先做一項轉換，常用的牛頓解法等價于迭代再加權最小二乘[20]，假設參數當前估計是{,}1K，每次只允許一類{,}變化，在{,}點Taylor展開，得到對數似然函數的二次逼近：
| ${L_Q}\left( {{\beta _{0\ell }},{\beta _\ell }} \right) = - \frac{1}{{2N}}\sum\limits_{i = 1}^N {{\omega _{i\ell }}{{\left( {{z_{i\ell }} - {\beta _{0\ell }} - x_i^T{\beta _\ell }} \right)}^2} + C{{\left( {\left\{ {{{\tilde \beta }_{0k}},{{\tilde \beta }_k}} \right\}_1^K} \right)}^2},} $ |
其中，，是根據當前參數計算所得的值，針對固定的{,}1K，最后一項是常數。
那么，加入彈性網罰之后，需要求解的問題可以描述為
| $\mathop {\min }\limits_{\left( {{\beta _{0\ell }},{\beta _\ell }} \right) \in {R^{p + 1}}} \left\{ { - {L_Q}\left( {{\beta _{0\ell }},{\beta _\ell }} \right) + \lambda {P_\alpha }\left( {{\beta _\ell }} \right)} \right\}$ |
3.3 坐標下降法
與著名的最小角回歸(least angle regression,LARS)算法相比，坐標下降法[21]在求解lasso等一類凸優化問題時更簡單、有效，且容易推廣到其他相關問題，如彈性網等。其基本思想是每次只優化一維變量，這樣優化問題可以很快完成，且相關方程可以在變量循環中更新；而通常情況下，最小化變量在眾多參數中不隨循環而改變，因此整個迭代過程將很快完成。需要注意的是，該方法的前提條件是多個預測變量之間互不相關，這樣才能將多變量問題理解為多個單變量問題的組合。
對式(9)的問題，記：
| $L = - {L_Q}\left( {{\beta _{0\ell }},{\beta _\ell }} \right) + \lambda {P_\alpha }\left( {{\beta _\ell }} \right) = \frac{1}{{2N}}\sum\limits_{i = 1}^N {{\omega _{i\ell }}\left( {{{\left( {{z_{i\ell }} - {{\tilde g}_\ell }{{\left( {{x_i}} \right)}^{\left( j \right)}} - {x_{ij}}{\beta _{j\ell }}} \right)}^2}} \right. + C{{\left( {\left\{ {{{\tilde \beta }_{0k}},\tilde \beta } \right\}_1^K} \right)}^2} + \lambda {P_\alpha }\left( {{\beta _\ell }} \right)} ,$ |
其中，是除去xij一項的聯系函數。
針對當前的類別，每次只優化系數β的一維，其他維視為常數。假設當前估計是{,}1K，對系數β的第j維βj求導，當βj>0時：
| $\begin{array}{l} \frac{{\partial L}}{{\partial {\beta _{j\ell }}}} = \frac{1}{N}\sum\limits_{i = 1}^N {{\omega _{i\ell }}\left( {{z_{i\ell }} - {{\tilde g}_\ell }{{\left( {{x_i}} \right)}^{\left( j \right)}} - {x_{ij}}{\beta _{j\ell }}} \right)\left( { - {x_{ij}}} \right) + \lambda \left( {1 - \alpha } \right){\beta _{j\ell }} + \lambda \alpha } \\ = - \frac{1}{N}\sum\limits_{i = 1}^N {{\omega _{i\ell }}{x_{ij}}\left( {{z_{i\ell }} - {{\tilde g}_\ell }{{\left( {{x_i}} \right)}^{\left( j \right)}}} \right) + \lambda \alpha + \frac{1}{N}\sum\limits_{i = 1}^N {{\omega _{i\ell }}x_{ij}^2} {\beta _{j\ell }} + \lambda \left( {1 - \alpha } \right){\beta _{j\ell }}} \end{array}$ |
令上式等于0，可以求得βj的坐標更新形式：
| ${\beta _{j\ell }} = \frac{{\frac{1}{N}\sum\limits_{i = 1}^N {{\omega _{i\ell }}{x_{ij}}\left( {{z_{i\ell }} - {{\tilde g}_\ell }{{\left( {{x_i}} \right)}^{\left( j \right)}}} \right) - \lambda \alpha } }}{{\frac{1}{N}\sum\limits_{i = 1}^N {{\omega _{i\ell }}x_{ij}^2 + \lambda \left( {1 - \alpha } \right)} }}\left( {{\beta _{j\ell }} > 0} \right)$ |
當βj<0時，同樣可以得到類似的表達式。其他情況下βj=0，即有軟閾值形式：
| ${\beta _{j\ell }} = \frac{{S\left( {\frac{1}{N}\sum\limits_{i = 1}^N {{\omega _{i\ell }}{x_{ij}}\left( {{z_{i\ell }} - {{\tilde g}_\ell }{{\left( {{x_i}} \right)}^{\left( j \right)}}} \right),\lambda \alpha } } \right)}}{{\frac{1}{N}\sum\limits_{i = 1}^N {{\omega _{i\ell }}x_{ij}^2 + \lambda \left( {1 - \alpha } \right)} }}\left( {{\beta _{j\ell }} > 0} \right)$ |
總的來說，對λ值的序列，可以得到β的一組路徑解。對于每一個λ，根據當前觀測值，建立多值logistic回歸模型，即求解系數β的過程是一系列循環嵌套的迭代求解過程，每個循環嵌套于上一個循環中：
循環1：循環∈{1,2,…,K,1,2,…}，直到β收斂；
循環2：使用當前參數{,}1K更新LQ(β0,β)二次逼近；
循環3：通過坐標下降算法求解式(9)的罰加權最小二乘問題，得到βj。
4 實驗結果及分析
實驗中，我們采用帶彈性網罰的多值logistic模型進行學習和分類，為提高結果的科學性，采用10折交叉驗證方法，并重復100次得到測試錯誤率的均值和標準差。
本文方法中影響診斷結果的參數主要有兩個：彈性網罰系數λ和彈性網罰中的正則化參數α。λ起到控制模型自由度的作用，實驗中給定100個λ值的序列，將訓練數據送入彈性網多值分類器，將得到100種參數估計結果，選擇最低分類錯誤率對應的參數估計作為訓練結果，因此系數λ的優化在程序中自動完成。α控制參數估計的稀疏值，實驗中對α的優化過程是通過比較實現的，其不同取值的實驗結果如表 2所示，該實驗結果是利用了名詞性變量的虛擬編碼和歸一化的年齡特征，樣本維數為130維。可以看出，α取值為0.01時能取得更低的錯誤率，且標準差較小。表中列出的時間是一次10折交叉驗證的時間。
表2
不同α值實驗結果比較
Table2.
Experimental results with different values of α
在一次10折交叉驗證中，每一折訓練過程所用數據均有差異，因此模型最終所選變量并不完全相同，但與傳統的特征子集選擇相比，該方法選擇的特征子集和模型分類性能的穩定性更高。以α取值為0.01時的100次特征選擇結果為例，我們把與第6類對應的100個特征子集(對應β，=6)進行統計(按被選擇的次數排列特征，橫軸標度僅代表特征數，不代表特征標號)，如圖 1所示。
圖1
100次特征選擇結果統計
Figure1.
Statistic of 100 feature selection results
由圖 1可見，有70維左右的特征幾乎每次都被選擇，而其余60維左右的特征從未被選擇，說明彈性網罰logistic回歸方法在子集選擇上具有很高的穩定性。
實驗中采用了多種方式處理名詞性變量和數值性變量。對于名詞性變量，對比了原始特征和虛擬編碼特征，對于數值性變量年齡，采用了4種處理方法，分別是：歸一化的年齡變量；離散化的年齡變量；虛擬編碼的年齡變量；去除年齡特征。其對應關系如表 3所示，不同處理方法的實驗結果對比如圖 2所示。
表3
年齡特征不同處理方法間的對應關系
Table3.
Correspondence of the feature of age based on different handing methods
圖 2中：A組采用原始特征，特征維數34；B組采用歸一化的原始特征，特征維數34；C組采用名詞性變量的虛擬編碼，年齡變量做歸一化處理，特征維數130；D組僅采用33維名詞性變量的虛擬編碼，去掉了年齡特征，特征維數129；E組采用名詞性變量的虛擬編碼和離散的年齡特征，特征維數130；F組采用所有34維變量的虛擬編碼，特征維數133。
圖2
特征不同處理方法的實驗結果比較
A,B,C,D,E,F分別表示不同的特征組合
Figure2. Comparison of results based on different handing methods of the featuresA,B,C,D,E and F represent different combinations,respectively
可以看出，與原始特征的實驗結果相比較，采用名詞性變量的虛擬編碼能顯著提高結果的預測精度，說明該方法是一種有效的名詞性變量處理方法。另外，E組在α取值為0.01時能取得更高的正確率，為98.34%，且標準差較小，說明在各種變量處理方法中，名詞性變量虛擬編碼和數值性變量離散化是一種可行的途徑。
我們采用上述E組實驗的特征處理方法結合多種分類器進行了訓練和測試，同樣采用100次的10折交叉驗證，實驗結果對比如圖 3所示。
圖3
多種分類器實驗結果比較
Figure3.
Comparison of experimental results based on different classifiers
實驗中我們對比了最近鄰(KNN-1)、3-近鄰(KNN-3)、K-近鄰(KNN)、SVM、貝葉斯分類(BC)、Logistic、Fisher、感知器(Perceptron)、偏最小二乘(PLS)等多種分類模型，各種模型對數據的處理方式相同，均將名詞性特征做虛擬編碼處理，數值性特征做離散化處理。可以看出，本文提出的方法對皮膚病數據集的分析起到了很好的效果，取得了較高的正確率(98.34%±0.002 7%)，且算法簡單穩定，易于編程實現，預測結果能夠為疾病的診斷提供較為可靠的參考。
5 結論
皮膚病數據集是一個同時包含了數值性變量和名詞性變量的數據集，傳統的處理方法沒有對兩種類型的變量區別對待。本文提出的名詞性變量虛擬編碼方法從特征取值的實際意義入手，將不同取值轉換到不同的維度，實驗中對比了采用原始數據與虛擬編碼的不同結果，證明該方法能夠顯著提高樣本分類的準確率，是一種處理名詞性變量的有效途徑。另外，針對虛擬編碼后的高維稀疏數據，采用帶彈性網罰的多值logistic模型擬合特征與疾病分類間的關系，能夠較好地反映數據的內在聯系，得到預測精度較高的模型，且加入彈性網罰起到了變量選擇的作用，省去了特征選擇的步驟，使整體處理過程簡單明了，易于實現。
本文提出的方法在實驗中取得了98.34%的平均正確率，100次實驗結果的標準差僅為0.002 7%，也就是說，對全部358個樣本的分類平均出錯不會超過6個，且算法步驟簡單，穩定性很強，可以考慮推廣到相似的分類問題中。同時，對名詞性變量和數值性變量不同處理方法的探索也具有較強的實際意義，值得進一步深入研究。
