如何在rspss中将多个变量合并变量名字改为x,y

以下试题来自:
多项选择题当相关系数r=-1时,变量x和y的相关关系为(  )。A.函数关系B.不完全相关关系C.完全正相关关系D.不相关关系E.完全负相关关系
为您推荐的考试题库
你可能感兴趣的试题
1A.加权算术平均数B.简单算术平均数C.极差D.众数E.中位数2A.经济增长率B.物价指数C.股票价格D.土地面积E.商品零售额3A.引用数据时一定要注明数据来源B.要评估第二手统计数据的可用价值C.指标的含义、口径、计算方法是否具有可比性D.对不完整的历史数据要根据需要和可能设法进行适当的补充E.不能纠正存在问题的历史数据4A.没有数字的单元格应空白B.左右两边不封口C.表中数据一般是右对齐D.列标题之间一般用竖线隔开E.行标题之间不必用横线隔开5A.国际收支顺差B.对外贸易持续逆差C.通货膨胀率比他国低D.提高利率E.总需求增长快于总供给
热门相关试卷
最新相关试卷科目:高中数学
两个具有线性相关关系的变量的一组数据为:
&&& 将以上数据,以x为自变量,y为因变量,得回归方程为=bx+a;将y为自变量,x为因变量,得回归方程为=b′y+a′. &&& 定义两个变量的相关关系数r的计算公式为r=,它可表示两个变量线性关系的强弱. &&& 试问r能否用上述两方程中的b,a与b′,a′表示?如能,怎样表示?
科目:高中数学
来源:2011年湖北省武汉市高二上学期期末考试数学理卷
题型:填空题
对两个具有线性相关关系的变量进行回归分析时,得到一个回归方程为,,则&&&&&&& .&
科目:高中数学
来源:学年湖北省武汉市武昌区高二(下)期末数学试卷(理科)(解析版)
题型:填空题
对两个具有线性相关关系的变量进行回归分析时,得到一个回归方程为,x∈{1,5,7,13,14},则=&&& .
精英家教网新版app上线啦!用app只需扫描书本条形码就能找到作业,家长给孩子检查作业更省心,同学们作业对答案更方便,扫描上方二维码立刻安装!
请输入姓名
请输入手机号图形--S语言教程--R语言教程--统计软件教程
S-PLUS有很强的图形功能,它可以用简单的函数调用迅速作出数据的各种图形,当你
熟悉了S图形的技术之后也可以指定许多图形选项按自己的要求定制图形。它的另一个特色是
同一个绘图函数对不同的数据对象可以作出不同的图形。例如,用7.1.2读入的cl数据框:
& plot(cl)
& plot(cl[,1])
第一个plot()调用绘制cl中三个列的散点图矩阵,第二个plot()调用绘制身高的散点图(
纵轴为身高值,横轴为下标)。
	最常用的绘图函数为plot(),用plot()作两个变量x与y的散点图,使用如下例的方法:
& attach(cl)
& plot(Height, Weight, main="体重对身高的回归",
+ xlab="身高", ylab="体重")
上例在R中运行成功(见图2),在S-PLUS中请将汉字字符串改为英文。上例也演示了S中
如何输入较长的语句:只要语句明显地未完成(比如,缺右括号),系统将给出一个加号作
为续行提示。如果输入“x &- 1+2”时要拆行,可以在赋值号后拆,可以在加号与2之间
拆,但是如果在x后拆则只能显示x的当前值,如果在1与加号之间拆只能把1赋给x。
	为了绘制连线图,只要在plot()函数中加type="l"选项,如:
& plot((1:50)/50, log((1:50)/50), type="l")
	可以绘制变量的茎叶图,如:
& stem(Weight)
The decimal point is 1 digit(s) to the right of the |
8 | 3445508
10 | 0332233
	绘制一个变量的盒形图,如:
& boxplot(Weight)
结果见图3。可以绘制几个变量并排的盒形图,比如先计算用上面的回归拟合结果存入p1
,然后绘制并排盒形图:
& p1 &- predict(fit1, cl)[,"predictor"]
& boxplot(Weight, p1)
	用hist()函数可以绘制直方图。例如:
& hist(Weight)
	用qqnorm()函数绘制正态概率图,如:
& qqnorm(Weight)
高级图形函数
S的图形函数分为两类:高级图形函数DD直接绘制图形并可自动生成坐标轴等附属
图形元素;低级图形函数DD可以修改已有的图形或者为绘图规定一些选择项。高级图形函
数总是开始一个新图。下面我们介绍常用的高级图形函数,以及用来修饰这些高级图形函数
的常用可选参数。
最常用的是plot()函数。比如,plot(x,y)(其中x,y是向量)对两个变量画散点图。用plot(z)
(其中z是一个定义了x变量和y变量的列表,或者一个两列的矩阵)也可以达到同样目的。如
果x是一个时间序列对象(时间序列对象用ts()函数生成),plot(x)绘制时间序列曲线图。
如果x是一个普通向量,则绘制x的值对其下标的散点图。如果x是复数向量则绘制虚部对实部
的散点图。如果f是一个因子,则plot(f)绘制f的条形图(每个因子水平的个数)。如果f是
因子,y是同长度的数值向量,则plot(f,y)对f的每一因子水平绘制y中相应数值的盒形图。
如果d是一个数据框,则plot(d)对d的每两个变量之间作图(散点图等)。
如果X是一个数值型矩阵或数据框,用pairs(X)可以绘制每两列之间的散点图矩阵。这在
变量个数不太多时可以同时看到多个变量的两两关系,变量太多时则难以绘制。
协同图(coplot)是一种多变量的探索性分析图形。其形式为coplot(y ~ x | z),其中x
和y是数值型向量,z是同长度的因子。对z的每一水平,绘制相应组的x和y的散点图。如:
& attach(cl)
& coplot(Weight ~ Height | Sex)
产生图6。对不同性别分别绘制了体重对身高的散点图。如果z是一个数值型变量,则coplot()
先对z的取值分组,然后对z的每一组取值分别绘图。甚至可以用如coplot(y~x | x1+x2)表示
对x1和x2的每一水平组合绘图。coplot()和pairs()函数缺省绘制散点图,但可以用一个panel=
参数指定其它的低级绘图函数,如lines,panel.smooth等。
	tsplot(x)绘制时间序列曲线图。多个参数时tsplot(x1, x2, ...)表示绘制多条曲线
,自动统一曲线取值范围。如果参数非时间序列对象则以下标1,2,3等为横坐标绘图。
	qqnorm(x), qqline(x), qqplot(x,y)作分位数-分位数图。qqnorm(x)对向量x作正
态概率(纵轴为次序统计量值,横轴为对应该次序统计量的标准正态分布分位数值)。qqline(x)
除作qqnorm(x)图之外还画一条拟合曲线。qqplot(x,y)把x和y的次序统计量分别画在x轴和y
轴以比较两个变量的分布。
	hist(x)作向量x的直方图。缺省时自动确定分组,也可以用nclass=参数指定分组个
数,或者用breaks=参数指定一个分组点向量。如果指定了prob=T则纵轴显示密度估计。
	S也可以作三维图或等值线图,函数为persp()和contour(),见图7和图8。
高级图形函数的常用选项
高级图形函数有一些共同的选项,作为函数的可选参数(自变量)。例如:
& plot(x, main="Graph of x")
其中的main就是一个可选参数,用来指定图形的标题。没有此选项时图形就没有标题。这
样的选项还有:
使函数向低级图形函数那样不是开始一个新图形而是在原图基础上添加。
暂不画坐标轴,随后可以用axis()函数更精确地规定坐标轴的画法。缺省值是axes=T,即
有坐标轴。
把x轴,y轴或两个坐标轴用对数刻度绘制。
规定绘图方式:
绘点并在中间用线连接
绘点并画线穿过各点
从点到横轴画垂线
阶梯函数;左连续
阶梯函数;右连续
不画任何点、线,但仍画坐标轴并建立坐标系,适用于后面用低级图形函数作图。
xlab="字符串"
ylab="字符串"
main="字符串"
sub="字符串"
定义x轴和y轴的标签。缺省时使用对象名。
图形的标题。
图形的小标题,用较小字体画在x轴下方。
低级图形函数
高级图形函数可以迅速简便地绘制常见类型的图形,但是,某些情况下你可能希望绘
制一些有特殊要求的图形。比如,你希望坐标轴按照自己的设计绘制,在已有的图上增加另
一组数据,在图中加入一行文本注释,绘出多个曲线代表的数据的标签,等等。低级图形函
数让你在已有的图的基础上进行添加。
常用的低级图形函数罗列如下:
points(x,y)
lines(x,y)
在当前图形上叠加一组点或线。可以使用plot()的type=参数来指定绘制方法,缺省时points()
画点,lines()画线。
text(x,y, labels, ...)
在由坐标x和y给出的位置标出由labels指定的字符串。labels可以是数值型或字符型的向
量,labels[i]在x[i],y[i]处标出。
abline(a, b)
abline(h=y)
abline(v=x)
在当前图形上画一条直线。两个参数a, b分布给出截距和斜率。指定h=参数时绘制水平线
,指定v=参数时绘制垂直线。以一个最小二乘拟合结果lm.obj作为参数时由lm.obj的$coefficients
成员给出直线的截距和斜率。
polygon(x, y, ...)
以由向量x给出的横坐标和向量y给出的纵坐标为顶点绘制多边形。可以用col=参数指定一
个颜色填充多边形内部。
legend(x, y,
legend, ...)
density=v)
legend( , fill=v)
legend(, col=v)
legend(, lty=v)
legend(, pch=v)
legend函数用来在当前图形的指定坐标位置绘制图例。图例的说明文字由向量legend提供
。至少下面的v值要给出以确定要对什么图例进行说明,v是长度与legend相同的向量。
angle参数指定几种阴影斜角。
density参数指定几种阴影密度。
fill参数指定几种填充颜色。
col参数指定几种颜色。
lty参数指定几种线型。
pch参数指定几种散点符号。为字符型向量。
marks参数也指定几种散点符号,但使用散点符号数值代号,为数值型向量。
title(main, sub)
绘制由main指定的标题和由sub指定的小标题。
axis(side, ...)
绘制一条坐标轴。这之前的绘图函数必须已经用axes=F选项抑制了自动的坐标轴。参数side
指定在哪一边绘制坐标轴,取值为1到4,1为下边,然后逆时针数。可以用at=参数指定刻度
位置,用labels参数指定刻度处的标签。
低级图形函数一般需要指定位置信息,其中的坐标指的是所谓用户坐标,即前面的高级图
形函数所建立的坐标系中的坐标。坐标可以用两个向量x和y给出,也可以由一个两列的矩阵
给出。如果交互作图可以用下面介绍的locator()函数来交互地从图形中直接输入坐标位置。
交互图形函数
S的低级图形函数可以在已有图形的基础上添加新内容,另外,S还提供了两个函数locator
和identify可以让用户通过在图中用鼠标点击来确定位置。
函数locator(n, type)运行时会停下来等待用户在图中点击,然后返回图形中鼠标点击的
位置的坐标。等待点击时用鼠标中键点击可以选择停止等待,立即返回。参数n指定点击多少
次后自动停止,缺省为500次;参数type如果使用则可指定绘点类型,与plot()函数中的type
参数用法相同,在鼠标点击处绘点(线、垂线,等等)。locator()的返回值是一个列表,有
两个变量(元素)x和y,分别保存点击位置的横坐标和纵坐标。
例如,为了在已经绘制的曲线图中找一个空地方标上一行文本,只要使用如下程序:
& text(locator(1), "Normal density", adj=0)
text()函数的adj参数用一个数字表示文本串相对于给定的坐标的画法,adj=0表示给定坐
标为文本串左侧的坐标,adj=1表示给定坐标为文本串右侧的坐标,adj=0.5表示给定坐标为
文本串中间的坐标。
函数identify(x, y, labels)在运行时也会停下来等待用户点击,直到按了鼠标中键,任
何返回用户在图形中用鼠标点击的点的序号,点击时对点击的点加标签。参数x和y给出要识
别的各个点的坐标。labels参数指定点击某个点时要在旁边绘制的文本标签,缺省时标出此
点的序号,如果只需要返回值而不想画任何标记则可以在调用此函数时加一个plot=F参数。
注意identify()与locator()不同,locator()返回图中任意点击位置的坐标,而identify()
只返回离点击位置最近的点的序号。
例如,我们在向量x和y中有若干个点的坐标,运行如下程序:
& plot(x, y)
& identify(x, y)
这时显示转移到图形窗口,进入等待状态,用户可以点击图中特别的点,该点的序号就会
在旁边标出。为了结束,只要单击鼠标中键或单击右键并选择停止。返回结果为你点击的各
个点的序号:
[1] 4 6 7 8
图形参数的使用
前面我们已经看到了如何用main=,xlab=等参数来规定高级图形函数的一些设置。在
实际绘图,特别是绘制用于演示或出版的图形时,S用缺省设置绘制的图形往往不能满足我们
的要求。但是,S提供了一系列所谓图形参数,通过使用图形参数可以修改图形显示的所有各
方面的设置。图形参数包括关于线型、颜色、图形排列、文本对齐方式等各种设置。每个图
形参数有一个名字,比如col代表颜色,取一个值,比如col="red"是红色。每个图形设备有
一套单独的图形参数。
设置图形参数分为两种:永久设置与临时设置。永久设置使用par()函数进行设置,设置
后在退出前一直保持有效;临时设置则是在图形函数中加入图形参数,如上面的例子:
& text(locator(1), "Normal density", adj=0)
中的adj参数。
par()函数用来访问或修改当前图形设备的图形参数。如果不带参数调用,如:
………………
结果为一个列表,列表的各元素名为图形参数的名字,元素值为相应图形参数的取值。
如果调用时指定一个图形参数名的向量作为参数,则只返回被指定的图形参数的列表:
& par(c("col", "lty"))
[1] "black"
[1] "solid"
	调用时指定名字为图形参数名的有名参数,则修改指定的图形参数,并返回原值的列
& oldpar &- par(col=4, lty=2)
[1] "black"
[1] "solid"
	因为用par()修改图形参数是保持到退出以前都有效的,而且即使是在函数内此修改
仍是全局的,所以我们可以利用如下的惯用法,在完成任务后恢复原来的图形参数:
& oldpar &- par(col=4, lty=2)
………………(需要修改图形参数的绘图任务)
& par(oldpar)
# 恢复原始的图形参数
	除了象上面那样用par()函数永久修改图形参数,我们还可以在几乎任何图形函数中
指定图形参数作为有名参数,这样的修改是临时的,只对此函数起作用。例如:
& plot(x, y, pch="+")
就用图形参数pch指定了绘散点的符号为加号。这个设定只对这一张图有效,对以后的图
形没有影响。
图形参数详解
鉴于绘制有特殊需要的图形是S的一个强项,而使用图形参数是完成此类任务的重要
手段,我们在这里较详细地介绍S的各种图形参数。这些图形参数可以大体上分为以下的几个
大类,我们将分别介绍:
图形元素控制
坐标轴与坐标刻度
一、图形元素
图形由点、线、文本、多边形等元素构成。下列的图形参数用来控制图形元素的绘制
指定用于绘制散点的符号。绘制的点往往略高于或低于指定的坐标位置,只有pch="."没
有这个问题。
如果pch的值为从0到18之间的一个数字,将使用特殊的绘点符号。下例可以显示所有特殊
绘点符号:
& plot(c(0, 100), c(0, 100), type="n", axes=F, xlab='', ylab='')
& legend(10,90, as.character(0:9), pch=0:9)
& legend(50,90, as.character(10:18), pch=10:18)
指定画线用的线型。缺省值lty=1是实线。从2开始是各种虚线。
指定线粗细,以标准线粗细为单位。这个参数影响数据曲线的线宽以及坐标轴的线宽。下
例绘制正弦曲线图:
& oldpar &- par(lwd=2)
& x &- (0:100)/100*2*pi
& plot(x, sin(x), type="l", axes=F)
& abline(h=0)
& abline(v=0)
& par(oldpar)
指定颜色,可应用于绘点、线、文本、填充区域、图象。颜色值也可以用象"red","blue"
这样的颜色名指定。
用来指定字体的整数。一般font=1是正体,2是
分别用来指定坐标刻度、坐标轴标签、标题、小标题所用的字体。
指定文本相对于给定坐标的对齐方式。取0表示左对齐,取1表示右对齐,取0.5表示居中
。此参数的值实际代表的是出现在给定坐标左边的文本的比例,所以adj=-0.1的效果是文本
出现在给定坐标位置的右边并空出相当于文本10%长度的距离。
指定字符放大倍数。
二、坐标轴与坐标刻度
	许多高级图形带有坐标轴,还可以先不画坐标轴然后用axis()单独加。函数box()
用来画坐标区域四周的框线。
	坐标轴包括三个部件:轴线(用lty可以控制线型),刻度线,刻度标签。它们可以
用如下的图形参数来控制:
lab=c(5, 7, 12)
第一个数为x轴希望画几个刻度线,第二个数为y轴希望画几个刻度线,这两个数是建议性
的;第三个数是坐标刻度标签的宽度为多少个字符,包括小数点,这个数太小会使刻度标签
四舍五入成一样的值。
坐标刻度标签的方向。0表示总是平行于坐标轴,1表示总是水平,2表示总是垂直于坐标
mgp=c(3,1,0)
坐标轴各部件的位置。第一个元素为坐标轴位置到坐标轴标签的距离,以文本行高为单位
。第二个元素为坐标轴位置到坐标刻度标签的距离。第三个元素为坐标轴位置到实际画的坐
标轴的距离,通常是0。
坐标轴刻度线长度,单位是绘图区域大小,值为占绘图区域的比例。tck小于0.5时x轴和y
轴的刻度线将统一到相同的长度。取1时即画格子线。取负值时刻度线画在绘图区域的外面。
控制x轴和y轴的画轴方法。
取值为"s"(即standard)或"e"(即extended)的时候数据范围控制在最小刻度和最大刻
度之间。取"e"时如果有数据点十分靠近边缘轴的范围会略微扩大。这种画轴方式有时会在轴
的一边留下太大的空白。
取值为"i"(即internal)或"r"(此为缺省)使得刻度线都落在数据范围内部,而"r"方
式所留的边空较小。
取值设为"d"时会锁定此坐标轴,后续的图形都使用与它完全相同的坐标轴,这在要生成
一系列可比较的图形的时候是有用的。要解除锁定需要把这个图形参数设为其它值。
三、图形边空
	S中一个单独的图由绘图区域(绘图的点、线等画在这个区域中)和包围绘图区
域的边空组成,边空中可以包含坐标轴标签、坐标轴刻度标签、标题、小标题等,绘图区域
一般被坐标轴包围。见图9。
	边空的大小由mai参数或mar参数控制,它们都是四个元素的向量,分别规定下方、左
方、上方、右方的边空大小,其中mai取值的单位是英寸,而mai的取值单位是文本行高度。
& par(mai=c(1, 0.5, 0.5, 0))
& par(mar=c(4, 2, 2, 1))
	这两个图形参数不是独立的,设定一个会影响另一个。S缺省的图形边空常常太大,
以至于有时图形窗口较小时边空占了整个图形的很大一部分。通常我们可以取消右边空,并
且在不用标题时可以大大缩小上边空。例如下例可以生成十分紧凑的图形:
& oldpar &- par(mar=c(2,2,1,0.2))
& plot(x,y)
	在一个页面上画多个图时边空自动减半,但我们往往还需要进一步减小边空才能使多
个图有意义。
四、一页多图
	R可以在同一页面开若干个按行、列排列的窗格,在每个窗格中可以作一幅图。
每个图有自己的边空,而所有图的外面可以包一个“外边空”,见图10。
	一页多图用mfrow参数或mfcol参数规定,如:
& par(mfrow=c(3,2))
表示同一页有三行两列共六个图,而且次序为按行填放。类似地,
& par(mfcol=c(3,2))
规定相同的窗格结构,但是次序为按列填放,即先填满第一列的三个再填第二列。要取消
一页多图只要再运行
& par(mfrow=c(1,1))
	缺省时无外边空。为了规定外边空大小,可以用omi参数或oma参数。omi参数使用英
寸为单位,oma参数以文本行高为单位,两个参数均为四个元素的向量,分别给出下、左、上
、右方的边空大小。如:
& par(oma=c(2,0,3,0))
函数mtext用来在外边空加文字标注。其用法为
mtext(text, side = 3, line = 0, outer = FALSE)
其中text为要加的文本内容,side表示在哪一边写(1为下,2为左,3为上,4为右),line
表示边空从里向外数的第几行,最里面的一行是第0号,outer=TRUE时使用外边空,否则会使
用当前图的边空。例如:
& par(mfrow=c(2,2), oma=c(0,0,3,0), mar=c(2,1,1,0.1))
& plot(x);plot(y);boxplot(list(x=x,y=y));plot(x,y)
& mtext("Simulation Data", outer=T, cex=1.5)
	在多图环境中还可以用mfg参数来直接跳到某一个窗格,比如
& par(mfg=c(2,2,3,2))
表示在三行两列的多图环境中直接跳到第二行第二列位置。mfg参数的后两个表示多图环
境的行、列数,前两个表示要跳到的位置。
	可以不使用多图环境而直接在页面中的任意位置产生一个窗格来绘图,参数为fig,
& par(fig=c(4,9,1,4)/10)
此参数为一个向量,分别给出窗格的左、右、下、上边缘的位置,取值为占全页面的比例
,比如上面的例子在页面的右下方开一个窗格作图。
S作图支持各种图形设备,其中常用的是显示器和PostScript打印机。在一个S运行期
间可以有多个图形设备同时存在。在R中,用
打开图形窗口绘图,在S-PLUS中,用
& win.graph()
打开图形窗口绘图。再次调用这样的函数将打开第二个图形窗口。用
& dev.list()
可显示以打开的图形设备的列表。
	要关闭一个图形设备,用
& dev.off()
这可以使得图形得以完成,例如对于postscript设备关闭设备时可完成打印或存盘。用graphics.off()
函数可以关闭所有打开的图形设备。
	MS Windows下的R可以把显示窗口中的图形复制到剪贴板或存为各种格式的图形文件
,包括WMF、PostScript、PNG、BMP、JPEG,这样我们可以用R生成所需图形然后存为需要的
格式。MS Windows下的S-PLUS也具有类似功能。
	各版本的R和S-PLUS都支持生成PostScript图形的功能,生成的图形可以直接用于LaTeX
排版。如果用MS Word排版则可把屏幕图形存为WMF等格式。生成PostScript文件的设备可以
用如下函数打开:
& postscript(file="result1.ps", horizontal=FALSE, width=5, height=3)
这时用图形命令生成一个页面的图形,然后用dev.off()关闭设备,则可生成文件result1.ps
。postscript()函数中horizotal参数指定是否将图旋转90度使得x轴平行于纸的长边,width
和height规定图的宽和高,单位是英寸。
	在打开了多个图形设备后可以用dev.set()函数来选择当前设备,dev.next()和dev.prev()
分别返回下一个和上一个图形设备。比如dev.set(dev.prev())选择上一个图形设备。在日常学习或工作中经常会使用线性回归模型对某一事物进行预测,例如预测房价、身高、GDP、学生成绩等,发现这些被预测的变量都属于连续型变量。然而有些情况下,被预测变量可能是二元变量,即成功或失败、流失或不流失、涨或跌等,对于这类问题,线性回归将束手无策。这个时候就需要另一种回归方法进行预测,即Logistic回归。
在实际应用中,Logistic模型主要有三大用途:
1)寻找危险因素,找到某些影响因变量的"坏因素",一般可以通过优势比发现危险因素;
2)用于预测,可以预测某种情况发生的概率或可能性大小;
3)用于判别,判断某个新样本所属的类别。
Logistic模型实际上是一种回归模型,但这种模型又与普通的线性回归模型又有一定的区别:
1)Logistic回归模型的因变量为二分类变量;
2)该模型的因变量和自变量之间不存在线性关系;
3)一般线性回归模型中需要假设独立同分布、方差齐性等,而Logistic回归模型不需要;
4)Logistic回归没有关于自变量分布的假设条件,可以是连续变量、离散变量和虚拟变量;
5)由于因变量和自变量之间不存在线性关系,所以参数(偏回归系数)使用最大似然估计法计算。
logistic回归模型概述
& & & & 广义线性回归是探索“响应变量的期望”与“自变量”的关系,以实现对非线性关系的某种拟合。这里面涉及到一个“连接函数”和一个“误差函数”,“响应变量的期望”经过连接函数作用后,与“自变量”存在线性关系。选取不同的“连接函数”与“误差函数”可以构造不同的广义回归模型。当误差函数取“二项分布”而连接函数取“logit函数”时,就是常见的“logistic回归模型”,在0-1响应的问题中得到了大量的应用。
&&&&&&&& Logistic回归主要通过构造一个重要的指标:发生比来判定因变量的类别。在这里我们引入概率的概念,把事件发生定义为Y=1,事件未发生定义为Y=0,那么事件发生的概率为p,事件未发生的概率为1-p,把p看成x的线性函数;
&&&&&&&& 回归中,最常用的估计是最小二乘估计,因为使得p在[0,1]之间变换,最小二乘估计不太合适,有木有一种估计法能让p在趋近与0和1的时候变换缓慢一些(不敏感),这种变换是我们想要的,于是引入Logit变换,对p/(1-p)也就是发生与不发生的比值取对数,也称对数差异比。经过变换后,p对x就不是线性关系了。
logistic回归的公式可以表示为:
其中P是响应变量取1的概率,在0-1变量的情形中,这个概率就等于响应变量的期望。
这个公式也可以写成:
可以看出,logistic回归是对0-1响应变量的期望做logit变换,然后与自变量做线性回归。参数估计采用极大似然估计,显著性检验采用似然比检验。
建立模型并根据AIC准则选择模型后,可以对未知数据集进行预测,从而实现分类。模型预测的结果是得到每一个样本的响应变量取1的概率,为了得到分类结果,需要设定一个阈值p0——当p大于p0时,认为该样本的响应变量为1,否则为0。阈值大小对模型的预测效果有较大影响,需要进一步考虑。首先必须明确模型预测效果的评价指标。
对于0-1变量的二分类问题,分类的最终结果可以用表格表示为:
其中,d是“实际为1而预测为1”的样本个数,c是“实际为1而预测为0”的样本个数,其余依此类推。
显然地,主对角线所占的比重越大,则预测效果越佳,这也是一个基本的评价指标——总体准确率(a+d)/(a+b+c+d)。
准确(分类)率=正确预测的正反例数/总数&&&&&&Accuracy=(a+d)/(a+b+c+d)
误分类率=错误预测的正反例数/总数&&&&&&&&Error rate=(b+c)/(a+b+c+d)=1-Accuracy
正例的覆盖率=正确预测到的正例数/实际正例总数
Recall(True Positive Rate,or Sensitivity)=d/(c+d)
正例的命中率=正确预测到的正例数/预测正例总数
Precision(Positive Predicted Value,PV+)=d/(b+d)
负例的命中率=正确预测到的负例个数/预测负例总数
Negative predicted value(PV-)=a/(a+c)
通常将上述矩阵称为“分类矩阵”。一般情况下,我们比较关注响应变量取1的情形,将其称为Positive(正例),而将响应变量取0的情形称为Negative(负例)。常见的例子包括生物实验的响应、营销推广的响应以及信用评分中的违约等等。针对不同的问题与目的,我们通常采用ROC曲线与lift曲线作为评价logistic回归模型的指标。
1)ROC曲线
设置了两个相应的指标:TPR与FPR。
TPR:True Positive Rate(正例覆盖率),将实际的1正确地预测为1的概率,d/(c+d)。
FPR:False Positive Rate,将实际的0错误地预测为1的概率,b/(a+b)。
TPR也称为Sensitivity(即生物统计学中的敏感度),也可以称为“正例的覆盖率”——将实际为1的样本数找出来的概率。覆盖率是重要的指标,例如若分类的目标是找出潜在的劣质客户(响应变量取值为1),则覆盖率越大表示越多的劣质客户被找出。
类似地,1-FPR其实就是“负例的覆盖率”,也就是把负例正确地识别为负例的概率。
TPR与FPR相互影响,而我们希望能够使TPR尽量地大,而FPR尽量地小。影响TPR与FPR的重要因素就是上文提到的“阈值”。当阈值为0时,所有的样本都被预测为正例,因此TPR=1,而FPR=1。此时的FPR过大,无法实现分类的效果。随着阈值逐渐增大,被预测为正例的样本数逐渐减少,TPR和FPR各自减小,当阈值增大至1时,没有样本被预测为正例,此时TPR=0,FPR=0。
由上述变化过程可以看出,TPR与FPR存在同方向变化的关系(这种关系一般是非线性的),即,为了提升TPR(通过降低阈值),意味着FPR也将得到提升,两者之间存在类似相互制约的关系。我们希望能够在牺牲较少FPR的基础上尽可能地提高TPR,由此画出了ROC曲线。
ROC曲线的全称为“接受者操作特性曲线”(receiver operating characteristic),其基本形式为:
当预测效果较好时,ROC曲线凸向左上角的顶点。平移图中对角线,与ROC曲线相切,可以得到TPR较大而FPR较小的点。模型效果越好,则ROC曲线越远离对角线,极端的情形是ROC曲线经过(0,1)点,即将正例全部预测为正例而将负例全部预测为负例。ROC曲线下的面积可以定量地评价模型的效果,记作AUC,AUC越大则模型效果越好。
当我们分类的目标是将正例识别出来时(例如识别有违约倾向的信用卡客户),我们关注TPR,此时ROC曲线是评价模型效果的准绳。
2)lift曲线
在营销推广活动中,我们的首要目标并不是尽可能多地找出那些潜在客户,而是提高客户的响应率。客户响应率是影响投入产出比的重要因素。此时,我们关注的不再是TPR(覆盖率),而是另一个指标:命中率。
回顾前面介绍的分类矩阵,正例的命中率是指预测为正例的样本中的真实正例的比例,即d/(b+d),一般记作PV。
在不使用模型的情况下,我们用先验概率估计正例的比例,即(c+d)/(a+b+c+d),可以记为k。
定义提升值lift=PV/k。
lift揭示了logistic模型的效果。例如,若经验告诉我们10000个消费者中有1000个是我们的潜在客户,则我们向这10000个消费者发放传单的效率是10%(即客户的响应率是10%),k=(c+d)/(a+b+c+d)=10%。通过对这10000个消费者进行研究,建立logistic回归模型进行分类,我们得到有可能比较积极的1000个消费者,b+d=1000。如果此时这1000个消费者中有300个是我们的潜在客户,d=300,则命中率PV为30%。此时,我们的提升值lift=30%/10%=3,客户的响应率提升至原先的三倍,提高了投入产出比。
为了画lift图,需要定义一个新的概念depth深度,这是预测为正例的比例,(b+d)/(a+b+c+d)。
与ROC曲线中的TPR和FPR相同,lift和depth也都受到阈值的影响。
当阈值为0时,所有的样本都被预测为正例,因此depth=1,而PV=d/(b+d)=(0+d)/(0+b+0+d)=k,于是lift=1,模型未起提升作用。随着阈值逐渐增大,被预测为正例的样本数逐渐减少,depth减小,而较少的预测正例样本中的真实正例比例逐渐增大。当阈值增大至1时,没有样本被预测为正例,此时depth=0,而lift=0/0。
由此可见,lift与depth存在相反方向变化的关系。在此基础上作出lift图:
&与ROC曲线不同,lift曲线凸向(0,1)点。我们希望在尽量大的depth下得到尽量大的lift(当然要大于1),也就是说这条曲线的右半部分应该尽量陡峭。
至此,我们对ROC曲线和lift曲线进行了描述。这两个指标都能够评价logistic回归模型的效果,只是分别适用于不同的问题:
如果是类似信用评分的问题,希望能够尽可能完全地识别出那些有违约风险的客户(不使一人漏网),我们需要考虑尽量增大TPR(覆盖率),同时减小FPR(减少误杀),因此选择ROC曲线及相应的AUC作为指标;
如果是做类似数据库精确营销的项目,希望能够通过对全体消费者的分类而得到具有较高响应率的客户群,从而提高投入产出比,我们需要考虑尽量提高lift(提升度),同时depth不能太小(如果只给一个消费者发放传单,虽然响应率较大,却无法得到足够多的响应),因此选择lift曲线作为指标。
相关R应用包
普通二分类 logistic 回归 用系统的 glm
因变量多分类 logistic 回归
有序分类因变量:用 MASS 包里的 polrb
无序分类因变量:用 nnet 包里的 multinom
条件logistic回归,用 survival 包里的 clogit
R运行逻辑回归
R可以让逻辑回归建模过程变得很简单。我们可以使用glm()函数进行逻辑回归建模,而且,它的拟合过程和我们之前所做的线性回归没有很大的区别。在这篇博客中,我们将介绍如何一步一步的进行逻辑回归建模。
我们用r来看一下曲线的形状:
x &- seq(from = -10, to = 10, by = 0.01)y = exp(x)/(1+exp(x))library(ggplot2)p &- ggplot(data = NULL, mapping = aes(x = x,y = y))p + geom_line(colour = 'blue')+ annotate('text', x = 1, y = 0.3, label ='y==e^x / 1-e^x', parse = TRUE)+ ggtitle('Logistic曲线')
p在0和1附近变化不敏感,ps:这条曲线帮了我一个大忙,最近做RTB竞价算法时,需要按照CTR设定一个CTR指导价,然后以CPC做曲线的修正,出于成本考虑,我们的出价不会高于某个固定值K,所以把这条曲线变换一下就派上大用场了。(扯远了)
这条曲线的主要特性是θ函数取值可以在-∞变到+∞,p的取值在[0,1],且极值端变化不敏感。这就是我们想要的。值得注意的是,经过logit变换,已经不是线性模型。
二、相关应用例子:Binary Logistic(因变量只能取两个值1和0虚拟因变量)
案例一:本文用例来自于John Maindonald所著的《Data Analysis and Graphics Using&R》一书,其中所用的数据集是anesthetic,数据集来自于一组医学数据,其中变量conc表示麻醉剂的用量,move则表示手术病人是否有所移动,而我们用nomove做为因变量,因为研究的重点在于conc的增加是否会使nomove的概率增加。
首先载入数据集并读取部分文件,为了观察两个变量之间关系,我们可以利cdplot函数来绘制条件密度图
install.packages("DAAG")
library(lattice)
library(DAAG)
head(anesthetic)
&&& move conc & logconc nomove
1 & &0 &1.0 0.0000000 & & &1
2 & &1 &1.2 0.1823216 & & &0
3 & &0 &1.4 0.3364722 & & &1
4 & &1 &1.4 0.3364722 & & &0
5 & &1 &1.2 0.1823216 & & &0
6 & &0 &2.5 0.9162907 & & &1
cdplot(factor(nomove)~conc,data=anesthetic,main='条件密度图',ylab='病人移动',xlab='麻醉剂量')
从图中可见,随着麻醉剂量加大,手术病人倾向于静止。下面利用logistic回归进行建模,得到intercept和conc的系数为-6.47和5.57,由此可见麻醉剂量超过1.16(6.47/5.57)时,病人静止概率超过50%。
anes1=glm(nomove~conc,family=binomial(link='logit'),data=anesthetic)
summary(anes1)
结果显示:
glm(formula = nomove ~ conc, family = binomial(link = "logit"),&
& & data = anesthetic)
Deviance Residuals:&
& & &Min & & & &1Q & &Median & & & &3Q & & & Max &
-1.76666 &-0.74407 & 0.03413 & 0.68666 & 2.06900 &
Coefficients:
& & & & & &&&&&&&& Estimate&Std. Error&&& z value&&&&&&& Pr(&|z|) &&
(Intercept) & -6.469 & & &2.418 &&&&& -2.675 &&&&&&& 0.00748 **
conc & & & & &&& 5.567 &&&&&2.044 &&&&&&&2.724&&&&&&&&&&0.00645 **
Signif. codes: &0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
& & Null deviance: 41.455 &on 29 &degrees of freedom
Residual deviance: 27.754 &on 28 &degrees of freedom
AIC: 31.754
Number of Fisher Scoring iterations: 5
&下面做出模型的ROC曲线
anes1=glm(nomove~conc,family=binomial(link='logit'),data=anesthetic)
对模型做出预测结果
pre=predict(anes1,type='response')
将预测概率pre和实际结果放在一个数据框中
data=data.frame(prob=pre,obs=anesthetic$nomove)
将预测概率按照从低到高排序
data=data[order(data$prob),]
n=nrow(data)
tpr=fpr=rep(0,n)
根据不同的临界值threshold来计算TPR和FPR,之后绘制成图
for (i in 1:n){
&&&&& threshold=data$prob[i]
&&&& tp=sum(data$prob&threshold&data$obs==1)
&&&& fp=sum(data$prob&threshold&data$obs==0)
&&&& tn=sum(data$prob
&&&& fn=sum(data$prob
&&&& tpr[i]=tp/(tp+fn)& #真正率
&&&& fpr[i]=fp/(tn+fp)& #假正率
plot(fpr,tpr,type='l')
abline(a=0,b=1)
R中也有专门绘制ROC曲线的包,如常见的ROCR包,它不仅可以用来画图,还能计算ROC曲线下面面积AUC,以评价分类器的综合性能,该数值取0-1之间,越大越好。
library(ROCR)
pred=prediction(pre,anesthetic$nomove)
performance(pred,'auc')@y.values
perf=performance(pred,'tpr','fpr')
plot(perf)
还可以使用更加强大的pROC包,它可以方便的比较两个分类器,并且能自动标出最优临界点,图形看起来比较漂亮:
install.packages("pROC")
library(pROC)
modelroc=roc(anesthetic$nomove,pre)
plot(modelroc,print.auc=TRUE,auc.polygon=TRUE,grid=c(0.1,0.2),grid.col=c("green","red"),max.auc.polygon=TRUE,auc.polygon.col="blue",print.thres=TRUE)
上面的方法是使用原始的0-1数据进行建模,即每一行数据均表示一个个体,另一种是使用汇总数据进行建模,先将原始数据按下面步骤进行汇总
anestot=aggregate(anesthetic[,c('move','nomove')],by=list(conc=anesthetic$conc),FUN=sum)
结果如下:
&&&&& conc move nomove
1&&&& & 0.8&&&& 6 &&& 1
2&&&& & 1.0 &&& 4&&&& 1
3 &&&&& 1.2 &&& 2&& & 4
4 &&&& 1.4 &&&& 2&&&& 4
5 &&&& 1.6 &&&& 0&&&& 4
6 &&&& 2.5 &&& 0 &&&& 2
anestot$conc=as.numeric(as.character(anestot$conc))
anestot$total=apply(anestot[,c('move','nomove')],1,sum)
anestot$total
[1] 7 5 6 6 4 2
anestot$prop=anestot$nomove/anestot$total
anestot$prop
[1] 0......0000000
对于汇总数据,有两种方法可以得到同样的结果,一种是将两种结果的向量合并做为因变量,如anes2模型。另一种是将比率做为因变量,总量做为权重进行建模,如anes3模型。这两种建模结果是一样的。
anes2=glm(cbind(nomove,move)~conc,family=binomial(link='logit'),data=anestot)
summary(anes2)
结果显示如下:
glm(formula = cbind(nomove, move) ~ conc, family = binomial(link = "logit"),
data = anestot)
Deviance Residuals:
1 2 3 4 5 6
0.267 0.500 0.26
Coefficients:
&&&&&&&&&&&&&&&&& Estimate& Std. Error& z value&&& Pr(&z)
(Intercept) -6.469 &&& & 2.419 && -2.675&& & 0.00748 **
conc&&&&&&&&& 5.567 &&&&&& 2.044&&& & 2.724&&&&& 0.00645 **
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 15.4334 on 5 degrees of freedom
Residual deviance: 1.7321 on 4 degrees of freedom
AIC: 13.811
Number of Fisher Scoring iterations: 5
anes3=glm(prop~conc,family=binomial(link='logit'),weights=total,data=anestot)
结果和上面的一样。
根据logistic模型,我们可以使用predict函数来预测结果,下面根据上述模型来绘图
x=seq(from=0,to=3,length.out=30)
&y=predict(anes1,data.frame(conc=x),type='response')
plot(prop~conc,pch=16,col='red',data=anestot,xlim=c(0.5,3),main='Logistic回归曲线图',ylab='病人静止概率',xlab='麻醉剂量')
lines(y~x,lty=2,col='blue')
案例二:利用iris数据集,进行逻辑回归二分类测试,该数据集是R语言自带得数据集,包括四个属性,和三个分类。逻辑回归我们用glm函数实现,该函数提供了各种类型的回归,如:提供正态、指数、gamma、逆高斯、Poisson、二项。我们用的logistic回归使用的是二项分布族binomial。
index &- which(iris$Species == 'setosa')
将种类为setosa的数据排除出我们需要的数据集
ir &- iris[- index,]
levels(ir$Species)[1] &- ''
生成训练集
split &- sample(100,100*(2/3))
ir_train &- ir[split,]
生成测试集ir_test &- ir[-split,]
通过训练集建立模型model &- glm(Species ~.,family=binomial(link='logit'),data=ir_train)
summary(model)
模型运行结果:
glm(formula = Species ~ ., family = binomial(link = "logit"),data = ir_train)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.339e-04 -2.100e-08 2.100e-08 2.100e-08 1.059e-04
Coefficients:
Estimate Std. Error z value Pr(&z)
(Intercept) -247.01 -0.004 0.997
Sepal.Length 12.45 .000 1.000
Sepal.Width -285.61 .003 0.998
Petal.Length 154.76
0.001 0.999
Petal.Width 869.60
0.004 0.997
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 9.0949e+01 on 65 degrees of freedom
Residual deviance: 4.0575e-08 on 61 degrees of freedom
Number of Fisher Scoring iterations: 25
通过anova()函数 对模型进行方差分析
anova(model, test="Chisq")
方差分析如下:
Analysis of Deviance Table
Model: binomial, link: logit
Response: Species
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(&Chi)
NULL 65 90.949
Sepal.Length 1 18.934 64 72.015 1.353e-05 ***
Sepal.Width 1 0.131 63 71.884 0.7176
Petal.Length 1 51.960 62 19.924 5.665e-13 ***
Petal.Width 1 19.924 61 0.000 8.058e-06 ***
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
下面通过McFadden R2指标进一步对模型进行分析
install.packages("pscl")
library(pscl)
pR2(model)
llh&&&&&&&&&&&&&&&&&&&&&&&& llhNull &&&&&&&&&&&&&& G2 McFadden&&&&& r2ML&&&&&&&&& r2CU
-2. -4. 9. 1. 7. 1.
为了得到分类结果,需要设定一个阈值p0——当p大于p0时,认为该样本的响应变量为1,否则为0。阈值大小对模型的预测效果有较大影响,需要进一步考虑。首先必须明确模型预测效果的评价指标。
求解训练模型的最佳阀值
对模型做出预测结果
model &- glm(Species ~.,family=binomial(link='logit'),data=ir_train)
pre1=predict(model,type='response')
将预测概率pre1和实际结果放在一个数据框中
data1=data.frame(prob=pre1,obs=ifelse(ir_train$Species=="virginica",1,0))
将预测概率按照从低到高排序
data1=data1[order(data1$prob),]
n=nrow(data1)
tpr=fpr=rep(0,n)
根据不同的临界值threshold来计算TPR和FPR,之后绘制成图
for (i in 1:n){
threshold=data1$prob[i]
tp=sum(data1$prob&threshold&data1$obs==1)
fp=sum(data1$prob&threshold&data1$obs==0)
tn=sum(data$prob
fn=sum(data$prob
tpr[i]=tp/(tp+fn) #真正率
fpr[i]=fp/(tn+fp) #假正率
plot(fpr,tpr,type='l')
abline(a=0,b=1)
下面通过pROC包自动标出最优临界点(0.506)
install.packages("pROC")
library(pROC)
modelroc1=roc(ifelse(ir_train$Species=="virginica",1,0),pre1)
plot(modelroc1,print.auc=TRUE,auc.polygon=TRUE,grid=c(0.1,0.2),grid.col=c("green","red"),max.auc.polygon=TRUE,auc.polygon.col="skyblue",print.thres=TRUE)
评估模型的预测效果
predict &- predict(model,type='response',newdata=ir_test)
predict.results &- ifelse( predict& 0.506,"virginica","versicolor")
misClasificError &- mean(predict.results != ir_test$Species)
print(paste('Accuracy',1-misClasificError))
[1] "Accuracy 1"
最后一步,我们将通过画ROC曲线,并计算其AUC面积,作为评估二类分类效果的一个典型测量
install.packages("ROCR")
library(gplots)
library(ROCR)
p &- predict(model,type='response',newdata=ir_test)
p.results &- ifelse( p& 0.5,1,0)
pr &- prediction(p.results, ifelse(ir_test$Species=="virginica",1,0))
prf &- performance(pr, measure = "tpr", x.measure = "fpr")
auc &- performance(pr, measure = "auc")
auc &-&auc@y.values[[1]]
&real &- ir_test$Species
data.frame(real,predict)
res &- data.frame(real,predict =ifelse(predict&0.5,'virginca','versicorlor'))
查看模型效果
基于R的案例
以下这个例子主要参考《Introduction to statistical learning with R》,强烈推荐这本书。尤其是回归部分,讲解的很透彻。
数据以Smarket为例,该数据包含年标准普尔500指数1250个观测值,9个变量。9个变量分别为:年份,前5个交易日的收益率,前一个交易日和交易额(单位:billions),今天的回报率和走势(上升或下降)。
df=read.csv("house.csv",header=TRUE)
data.frame': 1250 obs. of &9 variables:
&$ Year &&&&: num &01
&$ Lag1 &&&&: num &0.381 0.959 1.032 -0.623 0.614 ...
&$ Lag2 &&&&: num &-0.192 0.381 0.959 1.032 -0.623 ...
&$ Lag3 &&&&: num &-2.624 -0.192 0.381 0.959 1.032 ...
&$ Lag4 &&&&: num &-1.055 -2.624 -0.192 0.381 0.959 ...
&$ Lag5 &&&&: num &5.01 -1.055 -2.624 -0.192 0.381 ...
&$ Volume &&: num &1.19 1.3 1.41 1.28 1.21 ...
&$ Today &&&: num &0.959 1.032 -0.623 0.614 0.213 ...
&$ Direction: Factor w/ 2 levels "Down","Up": 2 2 1 2 2 2 1 2 2 2 ...
#以前5个交易日的收益率和前一个工作日的交易额为自变量,今天的走势做因变量做Logistic回归;
glm.fit=glm(Direction~Lag1+Lag2+Lag3+Lag4+Lag5+Volume,
&&&&&&&&&&&&data=Smarket,family="binomial")
summary(glm.fit)
glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 +&
&&&&Volume, family = "binomial", data = Smarket)
Deviance Residuals:&
&&&Min &&&&&1Q &Median &&&&&3Q &&&&Max &
-1.446 &-1.203 &&1.065 &&1.145 &&1.326 &
Coefficients:
&&&&&&&&&&&&&Estimate Std. Error z value Pr(&|z|)
(Intercept) -0.126000 &&0.240736 &-0.523 &&&0.601
Lag1 &&&&&&&-0.073074 &&0.050167 &-1.457 &&&0.145
Lag2 &&&&&&&-0.042301 &&0.050086 &-0.845 &&&0.398
Lag3 &&&&&&&&0.011085 &&0.049939 &&0.222 &&&0.824
Lag4 &&&&&&&&0.009359 &&0.049974 &&0.187 &&&0.851
Lag5 &&&&&&&&0.010313 &&0.049511 &&0.208 &&&0.835
Volume &&&&&&0.135441 &&0.158360 &&0.855 &&&0.392
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1731.2 &on 1249 &degrees of freedom
Residual deviance: 1727.6 &on 1243 &degrees of freedom
AIC: 1741.6
Number of Fisher Scoring iterations: 3
有时候,预测股市就是这么不靠谱。很多自变量没通过验证,接下来我们基于AIC准则用逐步回归做一下变量筛选;
注:AIC一般公式为 AIC=2k-ln(L),其中k为参数的个数,L为估计模型最大化的似然函数值;
step(glm.fit)
Start: &AIC=1741.58
Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume
&&&&&&&&&Df Deviance &&&AIC
- Lag4 &&&1 &&9.6
- Lag5 &&&1 &&9.6
- Lag3 &&&1 &&9.6
- Lag2 &&&1 &&0.3
- Volume &1 &&0.3
&none& &&&&&&&1.6
- Lag1 &&&1 &&1.7
#此处略去一百字;
Direction ~ 1
Call: &glm(formula = Direction ~ 1, family = "binomial", data = Smarket)
Coefficients:
(Intercept) &
&&&&0.07363 &
Degrees of Freedom: 1249 Total (i.e. Null); &1249 Residual
Null Deviance:&&&&&1731&
Residual Deviance: 1731&AIC: 1733&
这个结果足以让你崩溃,扔硬币都比这靠谱。原来不用任何变量预测是最靠谱的。
#模型的评估
glm.probs =predict(glm.fit,type ="response")
glm.probs[1:10]
#生成哑变量
contrasts(Smarket$Direction)
glm.pred=rep ("Down " ,1250)
glm.pred[glm .probs &.5]=" Up"
table(glm .pred ,Direction )
mean(glm.pred== Direction )
[1] 0.5216
通过混淆矩阵计算分类的准确性仅有52%,很不乐观;
挖掘本质上是挖掘有趣的的模式,以上数据说明我们白忙活了一场,但是我们通过实例学习了Logistic模型。
当然我们还可以通过数据变换或构造新的变量来提升拟合降低AIC,最终,模型讲究的还是泛化能力;
有时候:拟合很丰满,泛化很骨感。
因变量多分类 logistic 回归
1、概述:多元Logistic回归模型被用来建立有多个输出变量的模型,且这些预测变量通过一个线性组合变成为一个最终的预测变量。Multinomial Logistic 回归模型中因变量可以取多个值.例如一个典型的例子就是分类一部电影,它是“有趣”、“一般般”还是“无聊”。
所需的应用包:
library(foreign)
library(nnet)
library(ggplot2)
library(reshape2)
我们将以Titanic数据集作为实例进行分析。关于这个数据集,在网上流传了好几个免费的版本,而我则建议大家使用在Kaggle上的那个版本,因为它已经被用过很多次了(如果你要下载这个数据集,你需要在Kaggle注册一个账号)。这个(训练)数据集是关于一些乘客的信息(精确的数值为889),而这个比赛的目标就是基于诸如服务等级、性别、年龄等一些因素来预测一个顾客是否幸存(1表示乘客幸存,0则表示死亡)。在这里,你已经看到了,这个实例即有分类变量,也有连续变量。
数据清洗过程
当我们拿到一个真实的数据集的时候,我们需要确认这个数据集有哪些缺失值或者异常值。因此,我们需要做这方面的准备工作以便我们后面对此进行数据分析。第一步则是我们使用read.csv()来读取数据集,可以在这里下载数据集。确认参数na.strings相当于c(“”),因此,每个缺失值可以用NA表示。这可以帮助我们进行下一步的工作。
training.data.raw &- read.csv('train.csv',header=T,na.strings=c(""))
现在,我们要查阅一下缺失值,并观察一下每个变量存在多少异常值,使用sapply()函数来展示数据框的每一列的值。
sapply(training.data.raw,function(x) sum(is.na(x)))
PassengerId
sapply(training.data.raw, function(x) length(unique(x)))
PassengerId
在这里对缺失值进行可视化会有帮助:Amelia包有一个特别的作图函数missmap()来针对我们的数据集进行作图,并且对缺失值做标记。
library(Amelia)
missmap(training.data.raw, main = "Missing values vs observed")
这个变量有非常多的缺失值,我们都不会使用他们。我们也会把Passengerld去掉,因为它只有索引和票。使用subset()函数,我们可以从原始数据集中只选择与之相关的列。
data &- subset(training.data.raw,select=c(2,3,5,6,7,8,10,12))
观察一下缺失值
现在,我们需要观察一下别的缺失值。R可以通过在拟合函数内设定一个参数集来产生拟合一个广义线性模型。不过,我个人认为,我们应当尽可能的手动去除这些缺失值。这里有很多种方法来说实现,而一种常用的方法就是以均值,或中位数,或者是一个众数来代替缺失值。这里,我们会使用均值。data$Age[is.na(data$Age)] &- mean(data$Age,na.rm=T)当我们把它看成是分类变量的时候,使用read.table()函数或read.csv()函数,而默认情况下是分类变量被转成因子类型。一个因子就是R处理分类变量的方式。我们可以使用下面的代码来查阅它:
is.factor(data$Sex)
is.factor(data$Embarked)
为了能更容易的理解R是怎么处理分类变量的,我们可以使用contrasts()函数。这个函数向我们展示了R是如何优化变量的,并把它转译成模型。
contrasts(data$Sex)
contrasts(data$Embarked)
比如,你可以看见在变量sex,female可被作为参考。在乘船上的缺失值,由于这里只有两行,我们可以把这两行删掉(我们也可以用众数取代缺失值,然后保留这些数据点)。
data &- data[!is.na(data$Embarked),]
rownames(data) &- NULL
在继续拟合过程之前,让我们回顾一下数据清洗和重构的重要性。这样的预处理通常对于获得一个好而且预测能力强的模型是非常重要的。
我们把数据分成两个部分,训练和测试集。训练数据集可用于模型的拟合,而测试数据集则对此进行测试。
train &- data[1:800,]
test &- data[801:889,]
现在,让我们拟合一下这个模型。确认我们在glm()函数使用参数family=binomial。
model &- glm(Survived ~.,family=binomial(link='logit'),data=train)
现在,使用summary()函数来获得这个模型的结果。
summary(model)
glm(formula = Survived ~ ., family = binomial(link = "logit"),
data = train)
Deviance Residuals:
Coefficients:
Estimate Std. Error z value Pr(&|z|)
(Intercept)
& 2e-16 ***
-7.192 6.40e-13 ***
0.212026 -13.002
& 2e-16 ***
-4.547 5.43e-06 ***
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1065.39
degrees of freedom
Residual deviance:
degrees of freedom
AIC: 727.39
Number of Fisher Scoring iterations: 5
解读逻辑回归模型的结果
现在,我们可以分析这个拟合的模型,阅读一下它的结果。首先,我们可以看到SibSp,Fare和Embarked并不是重要的变量。对于那些比较重要的变量,sex的p值最低,这说明sex与乘客的生存几率的关系是最强的。这个负系数预测变量表明其它变量都一样,男乘客生存的机率更低。记住,逻辑模型的因变量是对数机率:ln(odds) = ln(p/(1-p)) = ax1 + bx2 + … + z*xn。male是一个优化变量,男性的生还机率下载2.75个对数机率,而年龄的增大则下降0.037个对数机率。现在,我们可以使用anova()函数来分析这个数据表的偏差了:
anova(model, test="Chisq")
Analysis of Deviance Table
Model: binomial, link: logit
Response: Survived
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev
981.79 & 2.2e-16 ***
741.77 & 2.2e-16 ***
724.28 2.881e-05 ***
0.000992 ***
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
空偏差和残差的不同展示了我们的模型是空模型(只有截距)。这个沟阅读,效果越好。分析这个图表,我们可以看到每增加一个变量,同时也会降低其中的偏差。再一次,添加Pclass、Sex和Age能显著的降低残差。其它变量对于模型的提高没什么用,尽管SibSp的p值很低。一个p值较大的值预示着一个没有变量的模型多少解释了模型的变化。最终,你想要看到的就是偏差的降低和AIC。然而,这里并不存在线性回归的残差R^2,McFaddenR^2指数可用于模型的拟合。
library(pscl)
pR2(model)
-354.6950111 -532.6961008
356.0021794
评估模型的预测能力
在上述的步骤中,我们简单的评估一下模型的拟合效果。现在,我们看一下使用心得数据集,它的预测结果y是什么样的。设定参数type=’response’,R可以输出形如P(Y=1|X)的概率。我们的预测边界就是0.5。如果P(Y=1|X)&0.5,y=1或y=0。记住,其它的决策边界可能是更好的选择。
fitted.results &- predict(model,newdata=subset(test,select=c(2,3,4,5,6,7,8)),type='response')
fitted.results &- ifelse(fitted.results & 0.5,1,0)
misClasificError &- mean(fitted.results != test$Survived)
print(paste('Accuracy',1-misClasificError))
"Accuracy 0.483"
精度为0.84,说明这个模型拟合效果不错。然而,你要记住,结果多少都依赖于手动改变数据集的数据。因此,如果你想要得到更精确的精度,你最好执行一些交叉检验,例如k次交叉检验。作为最后一步,我们会做ROC曲线并计算AUC(曲线下的面积),它常用于预测二元分类器的模型表现。ROC曲线是一种曲线,它可以通过设定各种极值来让正例律(TPR)来抵消反正例律(FPR),它就在ROC曲线之下。通常来说,一个预测能力强的模型应当能让ROC接近1(1是理想的)而不是0.5。
library(ROCR)
p &- predict(model, newdata=subset(test,select=c(2,3,4,5,6,7,8)), type="response")
pr &- prediction(p, test$Survived)
prf &- performance(pr, measure = "tpr", x.measure = "fpr")
auc &- performance(pr, measure = "auc")
auc &- auc@y.values[[1]]
下面是ROC图:
.cn/s/blog_vxwj.html&
.cn/s/blog_6934cecb0101a4qp.html
http://shujuren.org/article/143.html
阅读(...) 评论()

我要回帖

更多关于 r 删除变量 的文章

 

随机推荐