R语言CH7-统计回归分析课件.pptx

1、R语言基础与数据科学应用沈刚 主编人民邮电出版社内容导航C O N T E N T S数据的数值度量定性与定量数据概率分布与假设检验7.17.27.3回归分析7.4定性数据l 如果一个样本的取值属于一组已知的且互不重叠的类型，我们把这样的数据样本被称为定性数据，也称作分类数据。l iris数据是R自带的内置数据集之一，其中的鸢尾花分类信息就是一个定性数据的实际例子。str(iris)data.frame:150 obs.of 5 variables:$Sepal.Length:num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9.$Sepal.Width:num 3.

2、5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1.$Petal.Length:num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5.$Petal.Width:num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1.$Species :Factor w/3 levels setosa,versicolor,.:1 1 1 1 1 1 1 1 1 1.定性数据l 该数据集的最后一列是属性Species，包含了对鸢尾花品种的分类。品种的名称分别用setosa、versicolor、virginica等不同级别的因

3、子来表示，这就是一个定性数据的例子。iris$Species 1 setosa setosa setosa setosa setosa setosa 7 setosa setosa setosa setosa setosa setosa 145 virginica virginica virginica virginica virginica virginica Levels:setosa versicolor virginica定性数据l 数据变量的频数分布是对数据在一组不重叠的类别中出现次数的概括。species species.freq species.relfreq species.r

4、elfreqspecies setosa versicolor virginica 0.3333333 0.3333333 0.3333333定性数据l 使用柱状图可视化频数信息 barplot(species.freq)l 使用饼状图可视化频数信息 pie(species.freq)l 自定义饼图颜色 colors=c(red,yellow,blue)pie(species.freq,col=colors)定性数据 species s_species s_iris options(digits=3)mean(s_iris$Sepal.Length)1 5.01 tapply(iris$Sep

5、al.Length,iris$Species,mean)setosa versicolor virginica 5.01 5.94 6.59l 以setosa的花萼长度为例，可以分步骤得出其花萼长度的平均值。l option()函数指定R与用户的交互方式。l tapply()函数实现同样的效果。定量数据 head(faithful,3)eruptions waiting1 3.600 792 1.800 543 3.333 74l 计算喷发持续时间的频数分布步骤：使用range()函数得到数据上下界 将变量活动范围分成不重叠的区间 对变量进行区间分类 得到统计信息 l 预览faithful数据

6、集的前三行。定量数据l使用range()函数得到数据上下界 duration range(duration)1 1.6 5.1 l将变量活动范围分成不重叠的区间 breaks breaks 1 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5l 查看数据活动范围。l 以0.5为间隔的序列得到向量定量数据l对变量进行区间分类 duration.cut str(duration.cut)Factor w/8 levels 1.5,2),2,2.5),.:5 1 4 2 7 3 7 5 1 6.l得到统计信息 duration.freq duration.freq colors

7、 hist(duration,right=FALSE,col=colors,main=“老忠实喷发次数”,xlab=“持续时间（分钟),ylab=频数)定量数据l 累计频数分布 cumsum(1:10)1 1 3 6 10 15 21 28 36 45 55#breaks的长度比duration.freq多1，手动添加0。cumfreq0 plot(breaks,cumfreq0,main=“老忠实喷发持续时间,xlab=持续时间（分钟）,ylab=累积喷发频数)lines(breaks,cumfreq0)定量数据l 茎叶图 duration stem(duration)The decimal

8、 point is 1 digit(s)to the left of the|16|070355555588 18|000022233333335577777777888822335777888 20|00002223378800035778 48|00000022335800333 50|0370l 小数点位于“|”左侧一位，所以相邻两个树干的间隔为0.2。l 在树叶排列时按照升序排列，因此在最后一行“0370”中，第一个0表示的是5.0，而第二个0表示的则是5.1。l 第一行：1.60、1.67、1.70、1.73、1.75、1.75、1.75、1.75、1.75、1.75、1.78和1.

9、78分钟。定量数据 duration waiting head(cbind(duration,waiting),3)duration waiting1,3.600 792,1.800 543,3.333 74 plot(duration,waiting,xlab=喷发持续时间,ylab=等待时间)内容导航C O N T E N T S数据的数值度量定性与定量数据概率分布与假设检验7.17.27.3回归分析7.4均值l 观测样本的均值，等于数据值的和与数据个数的比值，是对数据值的中心位置的数值度量。如果数据的规模为n，其样本均值定义为：l 函数mean()返回的是一个数值向量的均值，可将na.r

10、m设置为TRUE排除向量中的缺失值。duration mean(duration)1 3.487783中位值、四分位数和百分位数l 中位值 median(duration)#调用median()函数1 4 l 四分位数 quantile(duration)#调用quantile()函数 0%25%50%75%100%1.60000 2.16275 4.00000 4.45425 5.10000l 百分位数 quantile(duration,c(.25,.5,.96)25%50%96%2.16275 4.00000 4.83572 四分位距l 四分位距是一个观测变量上下四分位数之差 IQR(d

11、uration)#调用IQR()函数1 2.2915 l 箱形图适用于直观地表示出四分位数和四分位距 par(mfrow=c(1,2)#绘图布局，用1行2列排列方式 boxplot(faithful$eruptions)boxplot(faithful$waiting)par(mfrow=c(1,1)summary()函数#R语言中提供的summary()函数可以直接对数据生成一些统计信息 summary(faithful)eruptions waiting Min.:1.600 Min.:43.0 1st Qu.:2.163 1st Qu.:58.0 Median:4.000 Median:

12、76.0 Mean :3.488 Mean :70.9 3rd Qu.:4.454 3rd Qu.:82.0 Max.:5.100 Max.:96.0 方差与标准差l 方差是对数据相对于其均值的分散程度的一个数值度量。l 函数var()计算给定数值向量的方差。观测变量的标准差是其方差的平方根，用函数sd()得到。var(duration)#调用var()函数1 1.3027 sd(duration)#调用sd()函数1 1.1414协方差l 数据集中两个变量x和y的协方差用于度量两者之间的线性相关度。如果协方差为正数，表明变量之间存在着正相关的关系，负的协方差表明两者之间是负相关关系。l co

13、v()函数可以计算出两个向量的协方差 cov(duration,waiting)#调用cov()函数1 13.978相关系数l 两个变量的相关系数等于它们的协方差除以各自标准差的乘积。l 当相关系数接近于1时，意味着变量之间线性正相关，在散点图上表现为散点几乎沿着一条斜率为正的直线分布。如果相关系数接近-1，则表示变量存在线性负相关关系，散点几乎落在一条斜率为负的直线附近。如果相关 系数为0，表明变量之间线性相关性很弱。l cor()函数就可以计算出它们的相关系数。cor(duration,waiting)#调用cor()函数计算相关系数1 0.901内容导航C O N T E N T S数据

14、的数值度量定性与定量数据概率分布与假设检验7.17.27.3回归分析7.4相关系数l R语言中提供了一组函数，分别以d、p、q和r开头，后面跟着概率分布的名称，用于返回一个给定参数的随机分布的概率密度、累积概率密度、分位数和按给定分布生成的（伪）随机数。以正态分布为例，其名称为norm，所对应的函数及参数默认值具有下列形式：dnorm(x,mean=0,sd=1,log=FALSE)pnorm(q,mean=0,sd=1,lower.tail=TRUE,log.p=FALSE)qnorm(p,mean=0,sd=1,lower.tail=TRUE,log.p=FALSE)rnorm(n,mea

15、n=0,sd=1)分布 R语言名称参数Beta分布betashape1,shape2二项式分布binomsize,prob柯西分布cauchylocation,scale分布chisqdf指数分布exprateF分布fdf1,df2Gamma分布gammashape,scale几何分布geomprob逻辑分布logislocation,scale对数正态分布lnormmeanlog,sdlog正态分布normmean,sd泊松分布poislambda学生t分布tdf均匀分布unifmin,max威布尔分布weibshape,scale二项分布l 二项式分布是一种离散概率分布，描述的是在n次独立

16、试验的最终结果。假定每次试验可以有两种结果，要么成功，要么失败。如果一次试验成功的概率为p，在n次独立试验中取得x次成功结果的概率如下，其中l 假设在一次考试中有10道多元选择题，每道题有4种可能的答案，其中只有一个答案是正确的。如果某个学生以随机方式回答了所有的问题，不及格（答对5道题及以下）的概率为：pbinom(5,size=10,prob=0.25)#返回P(x ppois(11,lambda=5)#默认计算下尾概率(x ppois(11,lambda=5,lower.tail=FALSE)#上尾(x11)1 0.0055连续均匀分布l 连续均匀分布是在从a到b的连续区间中随机选择数值

17、的概率分布。其概率密度函数的定义如下：l 例：例如，如果用户想在区间1,5中随机选取10个数，可以调用runif()函数返回10个随机值：runif(10,min=1,max=5)1 2.907 3.241 3.983 3.288 2.213 1.992 2.032 3.735 3.002 3.076指数分布l 指数分布表示了一系列随机重复发生的独立事件的到达时间的分布。假设到下一个事件发生的平均等待时间为，指数分布的概率密度函数形式如下：x dmu colors labels plot(x,dexp(x,1),type=l,lty=1,col=colors1,xlab=x,ylab=概率密度

18、,main=分布比较)for(i in 2:4)lines(x,dexp(x,1/dmui),lwd=2,lty=i,col=colorsi)legend(topright,inset=.05,title=均值,labels,lwd=1,lty=c(1,2,3,4),col=colors)正态分布l 正态分布用下列的概率密度函数定义，其中为均值而是方差：l 不同方差下均值为0的正态分布概率密度函数：l 假设在一次考试中全班的成绩符合正态分布，平均值是71分，标准差是12.8。计算一下，考试成绩在85分以上的概率是多少？前面已经使用过一些分布的累积概率函数，对于正态分布，这样的函数是pnorm(

19、).可以通过设置参数lower.tail=FALSE，也就是计算上尾概率求出所需结果。pnorm(85,mean=71,sd=12.8,lower.tail=FALSE)#85分以上，上尾累积概率1 0.137学生t分布统计假设检验统计假设检验l 举例来看，假设有一个制造商宣称所生产的灯泡平均寿命高达10000小时。检验时在30个灯泡样本中，发现其平均寿命只有9900小时。假定知道样本总体的标准差是120小时。给定5%的显著性水平，那么是否应该接受还是拒绝制造商的这一说法？。xbar mu0 sigma n z alpha z.alpha z.alpha#临界值 1-1.644854根据计算结

20、果，检验统计量-4.5644小于临界值-1.6449，因此，在5%的显著性水平上应该拒绝平均寿命超过10000小时的宣传。统计假设检验统计假设检验l 假设有一种食品商在曲奇包装上标签称每一块曲奇最多含有2克的饱和脂肪酸。在共计35块的样本中，发现平均的饱和脂肪酸含量为2.1克。假定已知总体标准差为0.25克，那么在5%的显著性水平上，是否应该拒绝或接受食品的标签？。xbar mu0 sigma n z alpha z.alpha z.alpha#临界值 1 1.644854检验统计量2.367大于临界值1.645，因此在5%的显著性水平，可以拒绝假设每一块曲奇只含有2克的饱和脂肪酸。统计假设检

21、验统计假设检验l 在南极洲发现了一群国王企鹅，去年它们的平均体重时15.4公斤。假设在今年的35只样本中，测量到的平均体重只有14.6公斤。如果已知总体标准差时2.5公斤。那么在5%的显著性水平下，能否拒绝企鹅平均体重与去年相同的假设？。xbar mu0 sigma n z alpha z.half.alpha c(z.half.alpha,z.half.alpha)1 1.9600 1.9600从计算结果可知，检验统计量-1.8931位于临界值-1.9600到1.9600之间。因此，在5%的显著性水平，无法拒绝这群企鹅平均体重与去年一样的假设。统计假设检验统计假设检验统计假设检验l 假设一个

22、制造商宣称一种灯泡的平均使用寿命超过10000小时。如果真实的平均使用寿命只有9900小时，总体方差为120小时。假设现在掌握了30个灯泡的样本，如何计算在5%的显著性水平，犯下第二类错误的概率是多少呢？。n-30;sigma sem alpha-.05;mu0 q mu pnorm(q,mean=mu,sd=sem,lower.tail=FALSE)#用SEM作为标准差，求上尾1 0.26196 内容导航C O N T E N T S数据的数值度量定性与定量数据概率分布与假设检验7.17.27.3回归分析7.4简单线性回归简单线性回归 eruption.lm coeffs coeffs#一元

23、线性回归的参数：截距、斜率(Intercept)waiting-1.87401599 0.07562795 plot(eruptions waiting,faithful,#绘图变量col=blue,#绘图参数main=老忠实线性回归结果,#标题xlab=等待时间,#x轴标签ylab=持续喷发时间)#y轴标签 fit abline(fit,col=red)#画出回归模型l lm参数：因变量eruptions，自变量waiting，数据集faithful。l 使用coefficients()来显示所得到的回归方程中的系数。简单线性回归 waiting duration duration(Inte

24、rcept)4.1762 newdata predict(eruption.lm,newdata)1 2 4.176220 1.907381 l 使用模型进行预测。简单线性回归简单线性回归简单线性回归简单线性回归l qq图可用于直观验证一组数据是否来自于某个给定的分布，或者验证两组数据是否来自同一分布。根据所讨论的分布计算出每个数据点的理论预期值，如果数据确实遵循假定的分布，那么在qq图上的点将大致散落在一条直线上。l 正态概率图就是一种把数据集与正态分布进行比较的图形化工具。例如，可以比较线性回归模型的标准化残差来检验残差是否真正地符合正态分布规律。qqnorm(eruption.stdre

25、s,ylab=标准化残差,xlab=正态得分,main=老忠实喷发持续时间)qqline(eruption.stdres)多元线性回归指标名说明crim按镇分布的人均犯罪率zn居住区域地块超过25,000 平方英尺的比例indus每个镇中非零售商业用地的比例chas与查尔斯河有关的哑数据（1表示河流范围，0为其他）nox一氧化氮浓度(parts per 10 million)rm每户平均房间数age1940年前修建的户主居住的单位数dis到5个波士顿就业中心的加权距离rad到达放射状高速公路方便程度的指数tax每万美元的全额房产税率ptratio每个镇的小学生师比 b1000(B-0.63)2

26、 其中B 是每个镇的黑人比例lstat低收入人口比例medv一千美元为单位的户主居住房屋的价格中位值多元线性回归l 把medv当作因变量，而把其余的指标作为自变量，可以建立起一个多元线性回归模型如下:l 使用lm()函数计算得出模型并保存在变量Boston.lm中 Boston.lm summary(Boston.lm)$r.squared1 0.7406427 多元线性回归逻辑回归逻辑回归l 封装测试数据到数据框newdata：newdata predict(am.glm,newdata,type=response)1 0.64181l 为了评价逻辑回归的分类效果，可以在mtcars数据集上生成混淆矩阵predict table(predict 0.5,mtcars$am)0 1 FALSE 18 1 TRUE 1 12逻辑回归l 用户也可以选择数据集中的一部分用于模型的选择，而把剩余的数据作为测试对象，检查逻辑回归模型的推广能力。例如，选择mtcars的前22条数据作为训练使用，而用后10条数据来测试。train test am.glm predict table(predict 0.5,test$am)0 1 FALSE 3 3 TRUE 0 4