基于惩罚函数的Beta回归变量选择方法.doc

基于惩罚函数的Beta回归变量选择方法本文受到国家自然科学基金青年项目（71201139，71303200）、国家社科基金（13&ZD148，13CTJ001）和教育部人文社科项目(12YJC790263)的资助。

内容提要：当因变量取值在(0,1)区间的比例数据时，Beta回归是常用的模型，该模型基于因变量服从Beta分布的假设，通过连接函数将样本服从Beta分布的均值与自变量联系起来。本文主要研究基于惩罚函数的Beta回归模型变量选择方法，研究了在LASSO SCAD和MCP三种惩罚方法下的极大似然估计，同时实现了变量选择和回归系数的估计。本文还研究了其渐进性质与Oracle性质，并且给出了计算机模拟。最后，将该方法应用到股息率的影响因素研究中。

关键词：Beta回归；变量选择；惩罚函数；Oracle性质
中图分类号：0212 文献标识码：A
Variable Selection for Beta Regression based on Penalized Function
Abstract: Beta regression is the popular model when the dependent variable is proportional data in interval (0,1), which assumes the dependent variable is beta distribution, connecting the mean of beta distribution and independent variables via link function. This article main focus on the variable selection for beta regression based on penalized function, giving maximum likelihood estimation with LASSO, SCAD and MCP penalty methods, the variable selection and estimation of coefficients can be conducted at the same time. We also give the proof of its asymptotical and oracle properties and the results of simulations. Finally, we apply this method to select significant factors for dividend rate.

Key Words: Beta Regression; Variable selection; Penalized function; Oracle Properties

一、引言
线性回归模型是最常用的回归模型，其假设因变量的取值区间为，但对因变量取值限制在区间的比例数据时，直接线性回归建模就不合适了，因为其拟合值很可能会落在区间外。而比例数据在实际中是很常见的，比如股息率、升学率、电视收视率、环境效率等。对于这种比例数据的因变量，传统的方法是将因变量转换到实数区间上，用自变量与转换后的因变量的均值建模。然而，这种方法的回归系数难以解释，而且比例数据通常具有非对称性，基于正态分布假设下做出的推断可能会导致预测失误。为解决传统回归所存在的问题，Ferrari 和 Cribari-Neto（2004）提出了Beta回归。由于Beta分布的密度函数的形状变化十分灵活，能够描述各种各样的分布，且回归系数具有较好的解释意义，因此广受欢迎。Beta回归模型的参数估计和统计诊断可以通过R包betareg实现（Cribari-Neto 和Zeileis，2010）。
随着计算机技术的飞速发展，极大地便利了数据的获取和存储，很多部门积累了大量的数据，比如股票市场的逐笔交易记录、商业银行交易记录等，同时这些数据的维度也越来越高，高维模型中噪声变量多，一个最重要的问题是模型选择，在许多情形下，这等价于如何选择自变量的问题。若不对它们加以筛选，引入过多变量，从理论上来说，一方面会导致模型不稳健，极大地降低了估计和预测精度；另一方面会加大模型的复杂度，无法突出最重要的自变量。
在回归分析中，往往先引入大量的自变量，然后筛选对因变量影响显著的自变量，这样可以减少模型的设定偏差。关于变量的选择，计量经济学中常用的方法有最优子集法（subset selection）、逐步（向前、向后）回归法(stepwise)、逐段法（stagewise）等（方匡南等，2014）。但是这些方法存在许多缺点，Fan和Li（2001）指出最优子集法遍历所有子集搜索最优解，因此计算成本高，特别是很大时，基本不可能实现求解; Breiman（1996）指出逐步回归法缺乏稳定性，对数据的微小变动非常敏感。惩罚函数法是目前高维数据研究领域广受欢迎的一类方法,它通过对模型的回归系数进行约束，使部分系数压缩为零，实现变量选择，克服了传统方法计算量大和稳定性差等缺点。常用的基于最小二乘方法的惩罚回归有LASSO(least absolute shrinkage and selection operator)（Tibshirani,1996），SCAD( the smoothly clipped absolute deviation)（Fan and Li,2001）和MCP (minimax concave penalty) （Zhang，2010）等。这种方法随后也被推广到极大似然估计中。与传统的变量选择方法不同，基于惩罚函数的极大似然估计具有很好的理论性质，当调和参数（tuning parameter）选择适当时，基于非凹惩罚函数的极大似然估计在变量选择方面具有oracle性质，即当真实的回归系数为0时，参数估计值依概率收敛到0，而非0系数的估计值与只用非0系数的变量进行回归时具有相同的表现。简而言之，基于非凹惩罚函数的极大似然估计估计好像提前知道了那些系数为0，哪些系数不为0，该种性质被称为Oracle性质。
本文的主要目的是研究基于惩罚函数下的Beta回归变量选择方法及其理论性质，并探讨该方法在实际中的应用。本文的结构安排如下：第二部分介绍了Beta回归模型；第三部分介绍基于惩罚函数的极大似然估计及其渐进性质和oracle性质；第四部分是计算机模拟；第五部分我们将该方法运用到利息率的影响因素分析中；最后是论文的小结。

二、Beta回归模型
Beta回归假设因变量服从Beta分布，其概率密度函数为：
，
其中，为gamma函数。的均值和方差为
， 。
当时，Beta分布为均匀分布。为了便于研究因变量的均值和精度的具体形式，我们对Beta分布进行参数化，令，，也就是，，因此，，其中。所以为因变量的均值，可以被解释为因变量的精度，对于固定的，越大，因变量的方差越小。的概率密度函数可以被改写为
，
其中。
在本文中，我们只考虑因变量落在区间的情况下，而Beta回归同样也适用于因变量落在区间的情况，此时我们只需要将因变量转换为即可。
设相互独立，其中服从均值为，方差为的Beta分布。模型假设的均值满足

其中是未知的回归系数向量，为个自变量的观测值，在模型中被假设为固定且已知。为严格单调且二阶可导的连接函数，可以将区间映射到实数域上。注意到的方差是和的函数，因此，因变量的方差变换也被考虑到了模型中。
连接函数有几种形式可供选择，这里我们采用logit连接函数，即，其它的连接函数参见Atkinson (1985)，McCullagh and Nelder (1989)等。在采用logit连接函数时，均值可以写成，其中。在这里回归系数有重要的解释意义，假设当其他自变量不变时，第个自变量增加了个单位，容易看出，也就是说，等于比例值之比，其中表示自变量取新的值时的均值，而表示自变量取原来的值时的均值。
个观测值的对数似然函数为，其中
，
极大似然估计通过令似然函数的一阶导为0获得，即和。此处极大似然函数没有显式解，我们需要通过非线性最优化算法来获取数值，如牛顿算法，高斯牛顿算法，详见Nocedal and Wright (1999)。最优化算法需要提供迭代的初始值，根据Ferrari and Cribari-Neto(2004)的建议，采用转换形式的因变量 对进行普通线性回归来所得的系数作为的初始值，即，其中。我们还需要的初始值，如前所述，由此。因为
。
因此，采用作为的初始值，其中通过将因变量 对进行回归的拟合值带入得到，也就是，而，其中为因变量对进行回归的残差。

三、基于惩罚函数的极大似然估计及其性质
基于惩罚函数的极大似然估计需在似然函数后添加惩罚项，即为
，
其中为惩罚函数，LASSO的惩罚函数为，SCAD的惩罚函数为，
其中且，MCP的惩罚函数为
，且
惩罚函数中的参数称为调整参数，在进行极大似然估计之前，调整参数必须事先选择确定。对带惩罚项的似然函数求极大值，所得的系数估计值就是基于惩罚函数的极大似然估计。通过选择适当的惩罚函数，对似然函数进行惩罚，可以使某些回归系数缩减为0，从而将该系数所对应的自变量从模型中剔除。
接下来，我们研究基于惩罚函数的极大似然估计量的渐进性质。令。不失一般性，假设。设为知道时的费雪信息阵。则基于惩罚函数的极大似然估计有以下的渐进性质：
(1).基于惩罚的极大似然估计以的速度收敛，其中，这表示如果，则基于惩罚函数的极大似然估计收敛（容易看出，则）。
(2).如果，收敛的估计值一定满足，且服从方差为的渐进正态分布，这表示基于惩罚函数的极大似然估计的性质表现与提前知道时相同。
接下来我们将证明上述结论。为了证明上述结论，我们先假设概率密度函数满足三个正则条件。
条件1：相互独立，的密度函数为。有共同的支撑。且满足和

条件2：费雪信息阵在处正定且有限。
条件3：对于几乎所有的，存在一个的开子集包含了真实参数，使得对于所有的，所有的三阶导数存在。并且存在函数使得对于所有的，。
此外，我们还需引入定理1、定理2以及引理1。设，。
定理1：相互独立，的密度函数 满足三个正则条件，如果，那么在的范围内，存在一个的局部最大值，其中。
证明：令，我们将证明对于任意的，总存在一个足够大的常数，使得，这表示在球内，至少有的概率存在一个局部最优解。因此存在一个最优解，使得最优解落在。
由于，我们有
，其中为的元素个数。
令为导数向量，对在处进行泰勒展开，有
，由于，注意为个样本的似然函数，为单个样本的似然函数，因此
。
对在处进行泰勒展开，有
。
而且
，
因为，所以，
所以，
，
，
由于，
所以。
因此当足够大时，在时是四者中最大的,因此，。
因此，根据定理1，如果选择合适的，如果存在一个收敛的带惩罚项的极大似然估计。
接下来我们证明带惩罚项的极大似然估计具有稀疏性质，即，证明如下。
引理1：相互独立，的密度函数 满足附录A中的条件1到条件3，假设，如果且，对任意满足的，和任意的常数，。
证明：对于任意满足的，和较小的，

为了证明上述引理，我们对进行泰勒展开，
其中位于和之间，由于，且，因此
对于满足的，我们有
，当且时，导数的符号完全由决定，因此

定理2：记，
，
为的元素个数，相互独立，的密度函数满足三个正则条件，假设，如果且，定理1中的收敛局部最优解一定满足：
(1)稀疏性：
(2)渐进正态性：
，
为知道 时的费雪信息阵。
证明：根据定理1，我们知道存在为的最优解，满足，因此

收敛于，
，，收敛于
根据Slutsky定理

收敛于。
因此

收敛于，其中，
。
如果，则。根据定理2，当时，带SCAD和MCP惩罚项的极大似然估计具有Oracle性质。对于LASSO来说，，由于收敛要求，所以，而定理2要求的条件，由于Lasso不能同时满足两个条件，因此Lasso不满足Oracle性质。

四、计算机模拟
接下来，我们通过计算机模拟进一步研究基于惩罚函数的极大似然估计量的性质。首先，我们需要确定惩罚函数中的调整参数，即和的值。我们通过五折交叉验证（fivefold cross-validation）来确定调整参数（Breiman，1995；Tibshirani 1996；Fu,1998)。交叉验证的具体过程为：记为全部的样本数据，将样本分为五部分，对于每个，分别用作为训练样本，作为测试样本。对于每个和，我们有参数估计值。交叉验证准则为
，
使CV最小时的即为我们要选择的调整参数。由于在不同的下，不同方法的相对好坏并不会发生变化，且模拟结果相差不大。为了表达简洁，我们在下面只给出取常用的3.7(SCAD)和 1.5(MCP)时的模拟结果。
1. 模拟1
在该例中，我们模拟100组样本，每组样本中有个观测值，三种情形，回归系数设为，精度设为三种情形。即，即模型中有个自变量的真实系数为0。自变量是服从均值为0的多元正态分布生成的，变量之间的相关系数设为，误差设为0.01和0.001两种，即当系数估计值小于时，认为系数为0。
表1左侧给出了当精度=120，样本量时100次模拟中系数被正确估计为0的平均个数。可以看出随着样本量的增大，三种惩罚方法的变量选择的精度都在提高。当样本量为50时，三种方法的变量选择效果差不多，被正确估计为0的变量平均有3.3个，比真实的0系数变量数少；当样本量为100时，LASSO和SCAD方法被正确估计为0的变量平均有3.54个，而MCP比前两者都要好，平均个数为3.92；在样本量为200时，SCAD的效果比LASSO好，而MCP的效果又比SCAD的效果好。说明MCP的方法会优于SCAD和LASSO，随着样本量的增大，SCAD的方法也将优于LASSO。
表1 模拟结果
误差
0.01
0.001
0.01
0.001

正确
错误
正确
错误
正确
错误
正确
错误
=50, =120

=100, =60

MLE
1.18
0
0.1
0
1.2
0.01
0.1
0
LASSO
3.3
0.01
2.45
0.01
3.74
0.02
2.98
0
SCAD
3.3
0.01
2.45
0.01
4.06
0.04
3.52
0
MCP
3.34
0.01
2.37
0.01
4.1
0.06
3.51
0
=100, =120

=100, =120

MLE
1.67
0
0.18
0
1.67
0
0.18
0
LASSO
3.54
0
2.51
0
3.54
0
2.51
0
SCAD
3.54
0
2.56
0
3.54
0
2.56
0
MCP
3.92
0
3.02
0
3.92
0
3.02
0
=200, =120

=100, =240

MLE
2.26
0
0.21
0
2.12
0
0.22
0
LASSO
4.28
0
2.98
0
3.19
0
1.6
0
SCAD
4.49
0
3.46
0
3.14
0
1.57
0
MCP
4.63
0
3.86
0
3.18
0
1.47
0

表1-2右边给出了样本量，精度时100次模拟中被正确估计为0的平均个数。可以看出当精度时，三种方法的效果差不多；在时，LASSO和SCAD的效果差不多，而MCP比前两者的效果要好；在精度时，SCAD的效果比LASSO好，而MCP的效果又比SCAD的效果好。由此可以说明，在精度较大时，三种方法的选择效果差不多的；随着精度的降低，MCP的方法会率先优于前两者。
2. 模拟2
在该例中，假设除了截距项外，所有变量的系数均为0的特殊情形，即。设样本量，精度，重复模拟100次。自变量的生成过程与模拟1相同。表2给出了100次模拟中系数被正确估计为0的平均个数。从表2可以看出，MLE的变量选择效果很差，三种惩罚方法的变量选择效果都很好，系数被正确估计为0的平均个数都很接近真实模型中0系数个数。说明这三种惩罚方法都可以很好地识别出不显著的变量。
表2 模拟2的结果
误差
0.01
0.001

正确
错误
正确
错误
MLE
1.02
0.2
0.07
0.2
LASSO
2.94
0.17
2.75
0.17
SCAD
2.94
0.17
2.75
0.17
MCP
2.96
0.18
2.69
0.18

3. 模拟3
在该例中，模拟了模型中仅截距项为0，其他回归系数均不为0的特殊情形，即。设样本量，精度，重复模拟100次。自变量的生成过程与模拟1相同。表3给出了100次模拟中系数被正确估计为0的平均个数。从表3可以看出，三种惩罚方法的变量选择效果都很好，系数被正确估计为0的平均个数都很接近真实模型中0系数个数，相对来讲，MCP的方法略好于LASSO和SCAD方法。说明这三种惩罚方法都可以很好地识别出不显著的变量。
表3 模拟3的结果
误差
0.01
0.001

正确
错误
正确
错误
MLE
0.46
0.2
0.06
0.2
LASSO
0.94
0.18
0.86
0.18
SCAD
0.94
0.18
0.86
0.18
MCP
0.98
0.17
0.94
0.17

五、股息率影响因素分析
接下来，我们利用本文提出的变量选择方法对我国上市公司的股息率影响因素进行分析。选取了沪深300指数的300只成分股作为我们的研究对象。之所以这样选择，是因为沪深300指数覆盖了沪深市场六成左右的市值，具有良好的市场代表性。数据来源于wind金融数据库。由于2013年上市公司财务报表还没有披露完整，因此我们选取2012年的沪深300指数的300只成分股的相关数据。经过计算，发现股息率的范围在（0,0.07），分布极其不均匀，而且有30家公司在2012年没有分红。对于股息率恰好等于0时，通过进行转换，使其取值范围在（0,1）间。对于这种分布不对称的比例数据，传统的方法就不适用了，需要使用beta回归。
影响上市公司股息率的因素很多，根据财务金融相关文献（阎大颖，2004；王化成，2010），我们选取了19个变量（详见表4），这些变量从上市公司的基本信息、市场表现、成长能力、营运能力、风险系数、偿债能力等角度反映公司的情况。由于变量众多，不一定所有的变量都和股息率有关，包含过多变量的beta模型往往会因为多重共线性导致部分变量的检验统计量不显著，这往往会降低模型的解释性并且影响模型的预测准确性，此外，模型一旦确定并选入一些无关的自变量，不仅干扰了对变量间关系的理解，而且浪费人力物力搜集这些变量信息，因此我们有必要进行变量筛选。利用本文提出的基于三种惩罚方法下的Beta回归变量选择方法筛选对股息率有影响的变量。
表4 变量解释及说明
编码
变量
计算说明
经济含义
Y
股息率
每股现金红利/股票年末价格
被解释变量
X1
总股本
股份总数
资本规模
X2
年成交量
2012年股票成交总量
股票流动性
X3
年成交额
2012年股票成交总额
市场表现
X4
年涨跌幅(%)
年末收盘价/年初开盘价
价格波动
X5
市净率
每股股价/每股净资产
投资潜力
X6
股票风险系数
利用股票的beta系数衡量
风险系数
X7
机构评级
2012年机构评级
市场评价
X8
每股收益
税后利润/股本总数
盈利能力
X9
每股净资产
股东权益/股本总数
X10
每股营业收入
营业收入/股本总数
X11
净资产收益率(%)
净利润/平均股东权益
X12
资产负债率(%)
总负债/总资产
偿债能力
X13
总资产周转率(%)
业务收入净额/平均资产总额
营运能力
X14
每股收益同比增长率
（2012年每股收益-2011年每股收益）/2011年每股收益
成长能力
X15
净资产同比增长率
（2012年净资产-2011年净资产）/2011年净资产
X16
税后每股股利
税后股利总额/普通股股数
分红政策
X17
年度分红总额
2012年分红总额
X18
3年累计分红占比(%)
截止指定年度连续三年的派息总额/截止指定年度连续三年年均归属母公司股东净利润合计×100%
再融资条件
X19
每股现金净额
现金流净额/股本总数
派现能力
表5给出了三种方法的变量选择和参数估计结果。从表5中可以看出，三种方法的选择结果比较一致。市净率、风险系数、综合评级、每股收益、税后每股股利、年度分红总额等变量对股息率的影响较大；而年成交量、年成交额、每股营业收入、每股收益同比增长率、净资产同比增长率、三年累计分红占比、每股现金流量净额、总股本等变量对股息率的影响不显著。
表5 回归结果

LASSO
SCAD
MCP

截距项
-0.1522
-0.1104
-0.1353
总股本
x1
0.0003
0.0004
0.0004
年成交量
x2
0
0
0
年成交额
x3
0
0.0002
0
年涨跌幅
x4
-0.0035
-0.0023
-0.0025
市净率
x5
-0.4343
-0.4670
-0.4529
风险系数
x6
-0.1662
-0.1281
-0.1281
综合评级
x7
-0.1434
-0.1469
-0.1496
每股收益
x8
-1.0554
-0.7849
-1.0332
每股净资产
x9
0.0533
0
0.0521
每股营业收入
x10
0.0013
0
0
净资产收益率
x11
0.1009
0.0907
0.1003
资产负债率
x12
0.0024
0.0003
0
总资产周转率
x13
-0.0707
0
-0.0549
每股收益同比增长率
x14
0
0
0
净资产同比增长率
x15
-0.0004
0
0
税后每股股利
x16
-1.8595
-2.3138
-2.2156
年度分红总额
x17
3.5354
3.7687
3.8671
三年累计分红占比
x18
0.0058
0.0056
0.0055
每股现金流量净额
x19
0.006309
0
0

精度
1.9627
1.9010
1.9452

六、小结
本文主要构建了基于惩罚函数的Beta回归变量选择方法。在极大似然函数基础上分别添加LASSO、SCAD和MCP三种惩罚项，同时实现变量选择和参数估计。用五折交叉验证（CV）方法选择调整参数，并从理论上证明了基于LASSO方法的极大似然估计不具有Oracle性质，而基于SCAD和MCP方法的极大似然估计具有Oracle性质。计算机模拟表明，
随着样本量的增大，三种惩罚方法的变量选择的精度都在提高，MCP的方法会优于SCAD和LASSO，并且随着样本量的增大，SCAD的方法也将优于LASSO；在精度较大时，三种方法的选择效果差不多的，随着精度的降低，MCP的方法会率先优于前两者。对于仅截距项为0和仅截距项非0的特殊情形，三种方法的变量选择效果都很好。最后，我们将本文提出的方法应用到我国上市公司的股息率影响因素分析中。

参考文献
[1]. Ferrari S. and Cribari-Neto F., Beta regression for modeling rates and proportions[J], Journal of Applied Statistics, 2004，7(31)：799–815.
[2]. Cribari-Neto F., Zeileis A., Beta regression in R[J], Journal of Statistical Software. 2010, 2(34)：1–24.
[3]. J. Fan and R. Li, Variable selection via nonconcave penalized likelihood and its oracle properties[J], Journal of American Statistical Association. 2001,456(96)：1348–1360.
[4]. Breiman, L., Heuristics of Instability and Stabilization in Model Selection[J], The Annals of Statistics, 1996, 6(24)：2350–2383.
[5]. Tibshirani R., Regression shrinkage and selection via the LASSO[J], Journal Royal Statistical Society. Series B, 1996,1(58)：267–288.
[6]. Zhang C., Nearly unbiased variable selection under minimax concave penalty[J], Annals of Statistics, 2010,2(38)：894–942.
[7]. Atkinson, A.C., Plots, Transformations and Regression: An Introduction to Graphical Methods of Diagnostic Regression Analysis[M]. New York: Oxford University Press, 1985.
[8]. McCullagh, P., Nelder, J.A., Generalized Linear Models[M], 2nd ed. London: Chapman and Hall, 1989.
[9]. Nocedal, J., Wright, S.J., Numerical Optimization[M]. New York: Springer–Verlag, 1999.
[10]. Breiman,L., Better Subset Regression Using the Nonnegative Garrote[J], Technometrics, 1995, 4(37)： 373–384.
[11]. Fu, W. J., Penalized Regression: The Bridge Versus the LASSO[J], Journal of Computational and Graphical Statistics, 1998,3(7), 397–416.
[12]. 方匡南,章贵军,张慧颖.基于lasso-logistic的个人信用风险预警方法[J]，数量经济技术经济研究，2014(2)：125-136.
[13]. 王化成,李春玲,卢闯.控股股东对上市公司现金股利政策影响的实证研究[J]，管理世界，2007(1):122-136
[14]. 阎大颖.中国上市公司控股股东价值取向对股利政策影响的证研究 [J]，南开经济研究，2004(6) : 94-105.