红斑鳞状皮肤病是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个，且算法步骤简单，稳定性很强，可以考虑推广到相似的分类问题中。同时，对名词性变量和数值性变量不同处理方法的探索也具有较强的实际意义，值得进一步深入研究。
