R语言与机器学习6logistic回归.docx
- 文档编号:23951546
- 上传时间:2023-05-22
- 格式:DOCX
- 页数:19
- 大小:491.44KB
R语言与机器学习6logistic回归.docx
《R语言与机器学习6logistic回归.docx》由会员分享,可在线阅读,更多相关《R语言与机器学习6logistic回归.docx(19页珍藏版)》请在冰豆网上搜索。
R语言与机器学习6logistic回归
R语言与机器学习(6)logistic回归
写在前面的废话
2014,又到了新的一年,首先祝大家新年快乐,也感谢那些关注我的博客的人。
现在想想数据挖掘课程都是去年的事了,一直预告着,盘算着年内完工的分类算法也拖了一年了。
本来打算去年就完成分类算法,如果有人看的话也顺带提提关联分析,聚类神马的,
可是,
。
借着新年新气象的借口来补完这一系列的文章,
但是sigmoid函数去拟合蓝线确实十分合适的。
于是我们可以考虑logistic回归模型:
假定有N个观测样本Y1,Y2,…,YN,设P(Yi=1|Xi)=π(Xi)为给定条件Xi下得到结果Yi=1的条件概率;而在同
样条件下得到结果Yi=0的条件概率为P(Yi=0|Xi)=1-π(Xi),于是得到一个观测值的概率P(Yi)=π(Xi)Yi[1-π(Xi)]
1-Yi假设各观测独立,对logistic回归模型来说,其对数似然函数为:
于是便可求解出logistic模型的MLE。
二、logit还是probit?
虽说sigmoid函数对边际递减的模型拟合良好,但是我们也要知道S型函数并非仅sigmoid函数一个,绝
大多数的累积分布函数都是S型的。
于是考虑F-1(P)(F为标准正态分布的累积分布函数)也不失为一个
很好的选择。
像这样的,对概率P做一点变换,让变换后的取值范围变得合理,且变换后我们能够有办法进
行参数估计的,就涉及到广义线性模型理论中的连接函数。
在广义线性模型中我们把log(P/(1-P))称为logit,
F-1(P)(F为标准正态分布的累积分布函数)称为probit。
那么这里就涉及到一个选择的问题:
连接函数
选logit还是probit?
logistic回归认为二分类变量服从伯努利分布,应当选择logit,而且从解释的角度说,p/
(1-p)就是我们常说的oddsratio,也就是软件报告中出现的OR值。
但是probit也有它合理的一面,首先,中心极限定理告诉我们,伯努利分布在样本够多的时候就是近似正态分布的;其次,从不确定性的角度考虑,probit认为我们的线性概率模型服从正态分布,这也是更为合理的。
我们来看一下经过变换后,自变量和P的关系是什么样子的:
如果你确实想知道到底你的数据用哪一个方法好,也不是没有办法,你可以看一下你的残差到底是符合
logit函数呢还是符合probit函数,当然,凭肉眼肯定是看不出来的,因为这两个函数本来就很接近,你可以通
过函数的假定,用拟合优度检验一下。
但通常,估计不会有人非要这么较真,因为没有必要。
但是有一点是
要注意的,logit模型较probit模型而言具有厚尾的特征,这也是为什么经济学论文爱用logit的原因。
我们以鸢尾花数据中的virginica,versicolor两类数据分类为例,看看两种办法分类有无差别。
probit.predictions
versicolorvirginica
versicolor473
virginica347
logit.predictions
versicolorvirginica
versicolor473
virginica347
从上图与比较表格均可以看出,两者差别不大。
三、多项式logit与orderlogit
对于二项分类模型的一个自然而然的推广便是多项分类模型。
我们借鉴神经网络里提到的异或模型,有:
按照上面这种方法,给定一个输入向量x,获得最大hθ(x)的类就是x所分到的类。
选择最大的hθ(x)十分好理解:
在类别选择问题中,不论要选的类别是什么,每一个类别对做选择的经济个体来说都有或多或少的效用(没有效用的类别当然不会被考虑),一个类别的脱颖而出必然是因为该类别能产生出最高的效用。
我们将多项logit模型的数学表述叙述如下:
假定对于第i个观测,因变量Yi有M个取值:
1,2,…,M,自变量为Xi,则多项logit模型为:
与logistic回归的似然估计类似,我们可以很容易写出多项logit的对数似然函数:
多项Logit模型虽然好用,但从上面的叙述可以看出,多项Logit模型最大的限制在于各个类别必须是对
等的,因此在可供选择的类别中,不可有主要类别和次要类别混杂在一起的情形。
例如在研究旅游交通工具的
选择时,可将交通工具的类别粗分为航空、火车、公用汽车、自用汽车四大类,但若将航空类别再依三家航空
公司细分出三类而得到总共六个类别,则多项Logit模型就不适用,因为航空、火车、公用汽车、自用汽车均
属同一等级的主要类别,而航空公司的区别则很明显的是较次要的类别,不应该混杂在一起。
在这个例子中,
主要类别和次要类别很容易分辨,但在其他的研究中可能就不是那么容易,若不慎将不同层级的类别混在一起
,则由多项Logit模型所得到的实证结果就会有误差。
对于分类模型,我们还会遇到被解释变量中有分类变量的情形。
对于连续变量解释离散变量,且被解释的离散变量是有顺序的(这个是和多项logit最大的区别)的情形,我们就需要考虑到orderlogit模型。
其数学模型叙述如下:
其中,F(.)表示累积分布函数,当F表示正态分布的分布函数时,对应的是orderprobit;F表示
logistic分布时,变对应orderlogit。
与logistic分布类似,我们可以很容易写出其对数似然函数:
四、dummyvariable
在logistic回归中,经常会遇到解释变量为分类变量的情形,比如收入:
高、中、低;地域:
北京、上海
、广州等。
这里对分类变量而言就涉及一个问题:
要不要将分类变量设置dummyvariable?
这个问题的答案在线性模型中很显然,必须要这么做!
!
!
如果我们不设置哑变量,而是单纯地赋值:
北
京=1,上海=2,广州=3,即我们将自变量视作连续性的数值变量,但这仅仅是一个代码而己,并不意味着地域间
存在大小次序的关系,即并非代表被解释变量(响应变量)会按此顺序线性增加或减少。
即使是有序多分类变量,
如家庭收入分为高、中、低三档,各类别间的差距也是无法准确衡量的,按编码数值来分析实际上就是强行规定
为等距,这显然可能引起更大的误差。
但是在logistic回归中,由于logit(p)变化的特殊性,在解释定序变量时,为了减少自由度(即解释变量个数),我们常常将定序变量(如家庭收入分为高、中、低)视为连续的数值变量,而且经济解释可以是XX每提高一个档次,相应的概率会提高expression(delta(XX))(expression的表达式还是很复杂的,不打了)。
当然减少变量个数是以牺牲预测精度为代价的。
毕竟数据处理是一门艺术而非一门技术,如何取舍还得具体问题具体分析。
当然,非定序的分类变量是万万不可将其视为数值变量的。
五、广义线性模型的R实现
R语言提供了广义线性模型的拟合函数glm(),其调用格式如下:
glm(formula,family=gaussian,data,weights,subset,
na.action,start=NULL,etastart,mustart,offset,
control=list(...),model=TRUE,method="glm.fit",
x=FALSE,y=TRUE,contrasts=NULL,...)
参数说明:
Formula:
回归形式,与lm()函数的formula参数用法一致
Family:
设置广义线性模型连接函数的典则分布族,glm()提供正态、指数、gamma、逆高斯、Poisson、二项分
布。
我们的logistic回归使用的是二项分布族binomial。
Binomial族默认连接函数为logit,可设置为probit。
Data:
数据集
鸢尾花例子使用的R代码:
logit.fit<-glm(Species~Petal.Width+Petal.Length,
family=binomial(link='logit'),
data=iris[51:
150,])
logit.predictions<-ifelse(predict(logit.fit)>0,'virginica','versicolor')
table(iris[51:
150,5],logit.predictions)
probit.fit<-glm(Species~Petal.Width+Petal.Length,
family=quasibinomial(link='probit'),
data=iris[51:
150,])
probit.predictions<-ifelse(predict(probit.fit)>0,'virginica','versicolor')
table(iris[51:
150,5],probit.predictions)
程序包mlogit提供了多项logit的模型拟合函数:
mlogit(formula,data,subset,weights,na.action,start=NULL,
alt.subset=NULL,reflevel=NULL,
nests=NULL,un.nest.el=FALSE,unscaled=FALSE,
heterosc=FALSE,rpar=NULL,probit=FALSE,
R=40,correlation=FALSE,halton=NULL,
random.nb=NULL,panel=FALSE,estimate=TRUE,
seed=10,...)
mlogit.data(data,choice,shape=c("wide","long"),varying=NULL,
sep=".",alt.var=NULL,chid.var=NULL,alt.levels=NULL,id.var=NULL,opposite=NULL,drop.index=FALSE,
ranked=FALSE,...)
参数说明:
formula:
mlogit提供了条件logit,多项logit,混合logit多种模型,对于多项logit的估计模型应写为:
因变量~0|自变量,如:
mode~0|income
data:
使用mlogit.data函数使得数据结构符合mlogit函数要求。
Choice:
确定分类变量是什么
Shape:
如果每一行是一个观测,我们选择wide,如果每一行是表示一个选择,那么就应该选择long。
alt.var:
对于shape为long的数据,需要标明所有的选择名称
选择wide的数据示例:
选择long的数据示例:
以fishing数据为例,来说明如何使用mlogit。
library(mlogit)
data("Fishing",package="mlogit")
Fish<-mlogit.data(Fishing,shape="wide",choice="mode")
summary(mlogit(mode~0|income,data=Fish))
这个输出的结果与nnet包中的multinom()函数一致。
由于mlogit包可以做的logit模型更多,所以这里就不在对nnet
包的multinom作介绍了,可以参见《根据EconometricsinR一书,将回归方法总结一下》一文。
程序包MASS提供polr()函数可以进行orderedlogit或probit回归。
用法如下:
polr(formula,data,weights,start,...,subset,na.action,
contrasts=NULL,Hess=FALSE,model=TRUE,
method=c("logistic","probit","cloglog","cauchit"))
参数说明:
Formula:
回归形式,与lm()函数的formula参数用法一致
Data:
数据集
Method:
默认为orderlogit,选择probit时变为orderprobit模型。
以housing数据为例说明函数用法:
house.plr<-polr(Sat~Infl+Type+Cont,weights=Freq,data=housing)
house.plr
summary(house.plr,digits=3)
这些结果十分直观,易于解读,所以我们在这里省略所有的输出结果。
再看手写数字案例:
最后,我们回到最开始的那个手写数字的案例,我们试着利用多项logit重做这个案例。
(这个案例的描述与数据参见《kNN算法》一章)
特征的选择可参见《神经网络》一章。
由于手写数字的特征选取很容易导致回归系数矩阵是降秩的,所以我们使用nnet包的multinom()函数代替mlogit()。
运行下列代码:
setwd("D:
/R/data/digits/trainingDigits")
names<-list.files("D:
/R/data/digits/trainingDigits")
data<-paste("train",1:
1934,sep="")
for(iin1:
length(names))
assign(data[i],as.matrix(read.fwf(names[i],widths=rep(1,32))))
label<-factor(rep(0:
9,c(189,198,195,199,186,187,195,201,180,204)))
feature<-matrix(rep(0,length(names)*25),length(names),25)
for(iin1:
length(names)){
feature[i,1]<-sum(get(data[i])[,16])
feature[i,2]<-sum(get(data[i])[,8])
feature[i,3]<-sum(get(data[i])[,24])
feature[i,4]<-sum(get(data[i])[16,])
feature[i,5]<-sum(get(data[i])[11,])
feature[i,6]<-sum(get(data[i])[21,])
feature[i,7]<-sum(diag(get(data[i])))
feature[i,8]<-sum(diag(get(data[i])[,32:
1]))
feature[i,9]<-sum((get(data[i])[17:
32,17:
32]))
feature[i,10]<-sum((get(data[i])[1:
8,1:
8]))
feature[i,11]<-sum((get(data[i])[9:
16,1:
8]))
feature[i,12]<-sum((get(data[i])[17:
24,1:
8]))
feature[i,13]<-sum((get(data[i])[25:
32,1:
8]))
feature[i,14]<-sum((get(data[i])[1:
8,9:
16]))
feature[i,15]<-sum((get(data[i])[9:
16,9:
16]))
feature[i,16]<-sum((get(data[i])[17:
24,9:
16]))
feature[i,17]<-sum((get(data[i])[25:
32,9:
16]))
feature[i,18]<-sum((get(data[i])[1:
8,17:
24]))
feature[i,19]<-sum((get(data[i])[9:
16,17:
24]))
feature[i,20]<-sum((get(data[i])[17:
24,17:
24]))
feature[i,21]<-sum((get(data[i])[25:
32,17:
24]))
feature[i,22]<-sum((get(data[i])[1:
8,25:
32]))
feature[i,23]<-sum((get(data[i])[9:
16,25:
32]))
feature[i,24]<-sum((get(data[i])[17:
24,25:
32]))
feature[i,25]<-sum((get(data[i])[25:
32,25:
32]))
}
data1<-data.frame(feature,label)
#降秩时mlogit不可用
#data10<-mlogit.data(data1,shape="wide",choice="label")
#m1<-mlogit(label~0|X1+X2+X3+X4+X5+X6+X7+X8+X9+X10+X11+X12+X13+X14+X15+X16+X17+X18+X19+X20+X21+X22+X23+X24+X25,data=data10)
library(nnet)
m1<-multinom(label~.,data=data1)
pred<-predict(m1,data1)
table(pred,label)
sum(diag(table(pred,label)))/length(names)
setwd("D:
/R/data/digits/testDigits")
name<-list.files("D:
/R/data/digits/testDigits")
data1<-paste("train",1:
1934,sep="")
for(iin1:
length(name))
assign(data1[i],as.matrix(read.fwf(name[i],widths=rep(1,32))))
feature<-matrix(rep(0,length(name)*25),length(name),25)
for(iin1:
length(name)){
feature[i,1]<-sum(get(data1[i])[,16])
feature[i,2]<-sum(get(data1[i])[,8])
feature[i,3]<-sum(get(data1[i])[,24])
feature[i,4]<-sum(get(data1[i])[16,])
feature[i,5]<-sum(get(data1[i])[11,])
feature[i,6]<-sum(get(data1[i])[21,])
feature[i,7]<-sum(diag(get(data1[i])))
feature[i,8]<-sum(diag(get(data1[i])[,32:
1]))
feature[i,9]<-sum((get(data1[i])[17:
32,17:
32]))
feature[i,10]<-sum((get(data1[i])[1:
8,1:
8]))
feature[i,11]<-sum((get(data1[i])[9:
16,1:
8]))
feature[i,12]<-sum((get(data1[i])[17:
24,1:
8]))
feature[i,13]<-sum((get(data1[i])[25:
32,1:
8]))
feature[i,14]<-sum((get(data1[i])[1:
8,9:
16]))
feature[i,15]<-sum((get(data1[i])[9:
16,9:
16]))
feature[i,16]<-sum((get(data1[i])[17:
24,9:
16]))
feature[i,17]<-sum((get(data1[i])[25:
32,9:
16]))
feature[i,18]<-sum((get(data1[i])[1:
8,17:
24]))
feature[i,19]<-sum((get(data1[i])[9:
16,17:
24]))
feature[i,20]<-sum((get(data1[i])[17:
24,17:
24]))
feature[i,21]<-sum((get(data1[i])[25:
32,17:
24]))
feature[i,22]<-sum((get(data1[i])[1:
8,25:
32]))
feature[i,23]<-sum((get(data1[i])[9:
16,25:
32]))
feature[i,24]<-sum((get(data1[i])[17:
24,25:
32]))
feature[i,25]<-sum((get(data1[i])[25:
32,25:
32]))
}
labeltest<-factor(rep(0:
9,c(87,97,92,85,114,108,87,96,91,89)))
data2<-data.frame(feature,labeltest)
pred1<-predict(m1,data2)
table(pred1,labeltest)
sum(diag(table(pred1,labeltest)))/length(name)
经整理,输出结果如下:
(左边为训练集,右边为测试集)
Tips:
oddsratio=p/1-p相对风险指数贝努力模型中P是发生A事件的概率,1-p是不发生A事件的概率所以p/1-p是发生与不发生的相对风险。
OR值等于1,表示该因素对A事件发生不起作用;OR值大于1,表示A事件发生的可能性大于不发生的可能性;OR值小于1,表示A事件不发生的可能性大于发生的可能性。
Furtherreading:
YvesCroissant:
EstimationofmultinomiallogitmodelsinR:
ThemlogitPackages
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 语言 机器 学习 logistic 回归