R实习指导模板.docVIP

第四章 R统计图形一图胜千言，优秀的图形可以将枯燥的数字以友好的形式展示在大家前面，从而让数据自己说话。统计图是统计描述的重要工具，它可以直观的反映出事物间的数量关系.在可视化技术如火如荼的发展之背景下，加上R软件在作图上的优越性,学习使用R统计作图已经成为必然趋势.本章主要向大家介绍R底层作图函数以及高水平作图函数，以及现今一些经典的R作图包.4.1 引入 图 4.1 mapped 2011想必大家看到这幅图，也许都想知道要怎么样才能画出来呢？是的，R就能完成，并且还能完成得非常漂亮.冰冻三尺非一日之寒，要想完成这么复杂但是美观的图案，还是得学习很多的东西，下面我们也就从R最基础的作图函数开始，达到你所见则能所画.或者想方法，让一般能实现的东西变得更加漂亮.4.2 R 统计描述作图4.2.1 直方图、经验分布图与QQ图 直方图的作图函数为hist(x),x为向量.其他参数详见帮助函数？hist. 与直方图相配套的核密度估计函数density(x)，其目的是用已知样本估计其密度.其他参数详见？density.经验分布函数的作图命令为ecdf(x).不论是直方图还是经验分布图，要从比较上鉴别样本是否近似于某种类型的分布较困难，而QQ图可以帮助我们鉴别样本的分布是否近似于正态分布.qqnorm()和qqline()提供了画正态QQ图和相应直线的方法，其他的参数详见帮助文件.例4.2.1 下面是15个学生的体重（单位:Kg）如下 75.0 64.0 47.4 66.9 62.2 62.2 58.7 63.5 66.6 64.0 57.0 69.0 56.9 50.0 72.0对此数据作出直方图、核密度估计图以及QQ图.命令为：whist(w,freq=F) #绘出密度直方图 lines(density(w),col=blue) #绘出蓝色的学生体重密度线 lines(44:76,dnorm(44:76,mean(w),sd(w),col=red)qqnorm(w);qqline(w) #绘出学生体重的QQ图我们得到下面的图：图4.2 学生体重的直方图、核密度估计图和QQ图由图4.2可以看出15个学生的体重是近似服从正态分布的.4.2.2 箱线图 箱线图的作图函数为boxplot()，参数详见？boxplot.在箱线图中，上（Q3）下（Q1）四分位数分别确定出中间箱体的顶部和底部，箱体中间的粗线是中位数所在的位置.由箱体向上下伸出的垂直部分成为“触须”，表示数据的散步范围，最远点为1.5倍四分位数间距，超出此范围的点称为异常指点，通常用“o”表示.在R软件中，InsectSprays是内部的数据框，它由两列数据构成，一列叫count,由数据构成，另一列叫spray,由因子构成，共有A-F六个水平，现在画出数据count在这六个水平下的箱线图，其命令为：boxplot(countspray,data=InsectSprays,notch=T,col=2:7)得到下面的图形：图4.3 InsectSpray的箱线图由上面可以看出在C、D两个水平上，存在极端异常值点（处于外限之外），都含有“o”标记.col=2:7代表箱线图的颜色向量从2（红色）开始到7（淡黄色），图中的C、E、F水平上四分位数离中位数的距离远比下四分位数距中位数的距离大，说明这些水平上，数据呈右偏分布，即数值较小的占多数；同理，B、D水平下，数据呈左偏分布，数值较大的占多数.当然boxplot()函数也还有其他很多的统计应用，希望您可以自己去学习.4.3 高水平作图函数4.3.1 高水平的函数作图简介 高水平作图函数有：plot()、pairs()、coplot()、dotchart()、contour()、image()、persp()等等(1) plot（）函数的使用方法： plot(x,y):x,y为向量，生成y关于x的散点图. plot(x):x是一时间序列，生成时序图；x是向量，生成x关于下标的散点图；x是复向量，生成复数的实部与虚部的散点图. plot(f,y):f是因子，y是数值向量,生成y关于f水平的箱线图. plot(df):df为数据框，生成df中所有变量的相对散点图.(2) pairs(x):x为矩阵，生成矩阵各列的散布图；x为数据框时，同plot(x)的功能.(3) coplot()函数的使用方法： coplot(ab|c)：在给定c值下，a关于b的散点图. coplot(ab|c+d):即按c,d划分下，a关于b的散点图.(4) dotchart(x):构造数据x的点图，点图中，y轴是数据x标记，x 轴是数据x的数值.(5) image(x,y,z,),contour(x,y,z,),persp(x,y,z,)其中x,y为数值型向量，z是与x和y对应的矩阵（z的行数是x的维数，z的列数是y的维数）.image()绘出三维图形的映像，contour()绘出三维图形的等值线，persp()绘出三维图形的表面曲线. 4.3.2 高水平作图例证例证1 R软件中，数据VADeaths给出了virginia(弗尼吉亚)州在1940年的人口死亡率.Rural Male Rural Female Urban Male Urban Female50-54 11.7 8.7 15.4 8.455-59 18.1 11.7 24.3 13.660-64 26.9 20.3 37.0 19.365-69 41.0 30.9 54.6 35.170-74 66.0 54.3 71.1 50.0画出该数据框的点图.R命令如下：dotchart(VADeaths,main=Death Rates in virginia -1940)dotchart(t(VADeaths),main=Death Rates in virginia -1940)图4.4 virginia州1940年的人口死亡率的点图例证2 对R中内部数据集iris中的前四列进行点图分析R命令如下：pairs(iris1:4, main = Andersons Iris Data - 3 species,+ pch = 21, bg = c(red, green3, blue)unclass(iris$Species)得到下面的点图：图4.5 Andersons iris的散点矩阵图例证 3 -2，2-2，2的正方形区域内绘函数的等值线图和三维曲面图. R程序如下： x-y f z contour(x,y,z,col=red,main=z=sinxsiny 的等值线图) persp(x,y,z,theta=30,phi=30,expand=.7,col=red) #增加了观 察图形的角度theta和图形的颜色.图4.6 的等值线图 画三维曲面图的时候，直接底层函数persp()画出的也许没那么美观，若选择rgl R包，名为：3D可视化装置，它便能画出优美的三维图，当然还有其他很多强大的功能:画3D图、电影动画、函数动画等等，这些都是基于OPENGL图形处理系统下的.那么我们用rgl画出上面这个函数的三维图时，命令为：library(rgl) zorder=rank(z) #z为函数f关系下的外积运算persp3d(x,y,z,col=rainbow(eger(max(zorder)zorder) # 加上了彩虹渐变颜色参数.得到下面的3D图：图4.7 的三维曲面图由3D图也可以很直观的看到，在上图达到尖刺的地方达到函数的最大或者最小值（上面的为最大值，下面的为最小值）.4.4. R 画图包 在cran官网上，在taskview一栏，分类列出了R的各种包，有一类别为Graphics即为图形包，现今共有32个，主要有：rgl、ggplot2、lattice、plotrix、vcd、RcolorBrewer、animation,gallery等等包，当然还有一些统计包中也会自带一些特殊的图形绘画函数：如corrplot绘制相关矩阵数据的可视化图. 我们主要通过一些例子，介绍下ggplot2和corrplot包的一些应用，其他详细的用法可以参见cran官网上的资料,也可以从google搜索到ggplot.pdf和corrplot.pdf.例 4.4.1 R内部数据quakes（Locations of Earthquakes off Fiji），包含了地震发生的经纬度，以及死亡人数、震级等等. 我们可以简单的画出地震发生代的级数与其经纬度的关系，画出散点图.命令为：require(ggplot2) #导入R包ggplot2head(quakes) #查看数据集quakes的前六行qpolt(lat,long,data=quakes,col=mag,size=exp(mag)得到下面的图形：图 4.8 地震震级散点图例4.4.2 给出我国外汇储备在1985-2011年的数据：包含GDP（国内生产总值）、FR（外汇储备量）、EM（出口额）、IM（进口额）、FDI（外商直接投资）、M2（广义货币供应量）、SD（短期债务）、LD（长期债务）等等.绘出外汇与其众多影响因素的相关矩阵可视化图形. R 命令：rtlibrary(corrplot)corrcorrplot(corr,method=pie,order=hclust,hclust.method=ward, addrect = 3,col=rev(rainbow(100)#采用饼图的方式，按照ward聚类(3类)排列.得到下面的图形:图4.9 外汇储备的影响因数相关矩阵可视化图图4.9鲜明活现的展示出了外汇储备与其他影响因素之间的相关程度.饼图中，白色部分越少的说明数值越接近于1，相关系数比较大，二者之间有密切的关系.并且按照ward聚类方法将相关关系较大的因素分为一组,形成了3类.与外汇fr影响比较大的因素为：GDP、EM、IM、M2、SD等.第五章 假设检验假设检验是统计推断中的一个重要内容，他是利用样本数据对某个事先做出的统计假设按照某种设计好的方法进行检验，判断此假设的正确与否.一般的统计理论都是基于统计假设之下的，所以懂得怎么做好统计运算的第一步就是要学会检验统计假设的正确性.在R中进行统计假设检验也有些特定的检验函数，当然也可以自定义一些函数. 假设检验分为参数检验和非参数检验，本章主要从这两方面来进行R的实际上机统计检验，为后面的一些统计建模的学习奠定基础.5.1 引例 例5.1 某种元件的寿命X（以h计）服从正态分布N（），其中均未知，现测得16只元件的寿命如下：159 280 101 212 224 379 179 264222 362 168 250 149 260 485 170问是否有理由认为元件的平均寿命大于225h? 例5.2 在平炉上进行一项试验以确定改变操作方法的建议是否会增加钢的得率，试验是在同一个平炉上进行的，每炼一炉钢时除操作方法外，其他条件都尽可能做到相同，先用标准方法炼一炉，然后用新方法炼一炉，以后交替进行，各炼了10炉，其得率分别为 标准方法 78.1 72.4 76.2 74.3 77.4 78.4 76.0 75.5 76.7 77.3 新方法 79.1 81.0 77.3 79.1 80.0 79.1 79.1 77.3 80.2 82.1设这两样本相互独立，且分别来自正态总体和，其中和未知.问新的操作能否提高得率？（取a=0.05）想必上面两个引例，大家在数理统计或者其他书籍中也见过，也学习过用理论的公式计算过，或者查表等等.下面我们就学习用R来解答这两个例题.5.2 t.test()函数和var.test()函数5.2.1 t.test()函数在R软件中，函数t.test()提供了t检验和相应的区间估计的功能，t.test()的使用格式如下：t.test(x,y=NULL,alternative=c(“two.sided”,“less”,“grater”),mu=0,paired=FALSE,var.equal=FALSE,conf.level=0. 95,)其中x,y是由数据构成的向量（如果只提供x,则作单个正态总体的均值检验，否则作两个总体的均值检验）；alternative表示备择假设；two.sided（默认）表示双边检验；greater表示单边检验；mu表示原假设；conf.level是置信水平，即1-a,通常是.2 var.test()函数 R软件中，var.test()函数提供作方差比的检验和相应的区间估计，该函数的使用格式：var.test(x,y,ratio=1,alternative=c(“two.sided”,“less”,“greater”),conf.level=0.95,)其中x,y是来自两样本数据构成的向量;ratio是方差比的原假设，默认值为1；alternative是备择假设；two.sided表示双边检验;less表示单边检验（）；greater表示单边检验（）5.2.3 应用实例解答 现在我们就可以利用所学的两个函数解答引例5.1的内容.解例5.1 按题意，需检验 此问题是单边检验问题.R程序如下： x t.test(x,alternative=greater,mu=225) One Sample t-testdata: x t = 0.6685, df = 15, p-value = 0.257alternative hypothesis: true mean is greater than 225 95 percent confidence interval: 198.2321 Inf sample estimates:mean of x 241.5解例5.2 按题意，需要假设 ，首先此题并未告诉标准方法与新方法得到的数据方差是否相等，因此应进行方差检验，用var.test()函数，如果证明了方差相等的假设的正确性，则选择t检验法，.R程序如下： x y var.test(x,y)F test to compare two variancesdata: x and y F = 1.4945, num df = 9, denom df = 9, p-value = 0.559alternative hypothesis: true ratio of variances is not equal to 1 95 percent confidence interval: 0.3712079 6.0167710 sample estimates:ratio of variances 1.494481上述结果表明在95%的置信水平下，p值为0.559大于0.05，所以接受原假设，因此我们可以采用t检验，用t.test()函数.认为这两个样本方差是相等的. t.test(x,y,var.equal=TRUE,alternative=less) #“less”在这里表达的含义如题中所假设的,与单样本均值检验中的含义不同. Two Sample t-testdata: x and y t = -4.2957, df = 18, p-value = 0.0002176alternative hypothesis:true difference in means is less than0 95 percent confidence interval: -Inf -1.908255 sample estimates:mean of x mean of y 76.23 79.43从计算结果可以看出，t.test()不但可以做单个总体的均值检验，还可以作两个总体的均值检验以及它们的均值区间估计，根据t.test（）函数的参数设置，我们也可得知道此函数还可以作成对数据（如二维数据）的t检验. x从而我们可以用t.test()函数完成均值检验与估计的全部工作.& 事实上，均值的区间估计与均值的假设检验本质上是对一个问题从不同角度讨论，有着内在的联系，这也是为什么t.test()将区间估计与假设检验放在一起的原因，可以使我们从多角度对问题进行判断，提高判断的准确性.5.3 一些非参数检验 在一些实际问题中，人们往往对总体的分布知之甚少，很难对总体的分布形式作出正确的假定，尽量从数据本身来获得所需要的信息的统计方法为非参数方法.对于非参数方法的检验问题称为非参数检验法，它涉及的范围很广，这里只介绍几种与R有关的、在应用上较为重要的检验法. 上一节讲到的参数检验目标是关于某个参数落在特定的范围内的假设是否正确，而现在所讲的非参数检验目标是针对分布的类型.5.3.1 chisq.test()函数Pearson拟合优度检验函数chisq.test()，其应用范围颇广，可以检验特定的数据是否服从某种分布，不过这些数据必须在检验前进行预处理进行分类，得到每个类别的频数数据，然后再进行检验.在学习chisq.test()这个函数之前必须知道两个跟此函数相关的函数cut()和table()函数. cut()函数是将变量的区域分成若干个区间，其使用方法为cut(x,breaks,labels=NULL,include.lowest=FALSE,right=TRUE,dig.lab=3,) x为数据构成的向量，breaks(简称br)是所分区间的端点值构成的向量. table()函数是计算因子合并后的个数，其使用方法为table(,exclude=c(NA,NaN),dnn=s(),deparse.level=1) chisq.test()函数的使用格式为chisq(x,y=NULL,correct=TRUE,p=rep(1/length(x),length(x),rescale.p=FALSE,simulate.p.value=FALSE,B=2000)其中x是由观测数据构成的向量或矩阵，y是数据向量（x为矩阵时，y无效）；correct是逻辑变量，表面是否用于连续修正，TRUE（默认值）表示修正；p是原假设落在小区间的理论概率，默认值表示均匀分布；rescale.p是逻辑向量，选择FALSE（默认值）表面要求输入的p值之和为1，选择TRUE时，并不要求这一点，程序将重新计算p值；simulate.p.value是逻辑变量（默认FALSE）,当为TRUE时，将用仿真的方法计算p值，此时B表示仿真的次数.5.3.2 Kolmogogrov-Smirnov检验 Kolmogogrov-Smirnov检验是通过经验分布与假设分布的上确界来构造统计量的，因此理论上可以检验任何分布，其用法为 ks.test(x, y, ., alternative = c(two.sided, less, greater),exact = NULL)其中参数详见帮助文件？ks.test5.3.3 fisher精确的独立检验在样本较小时（单元的期望频数小于4），需要用fisher精确检验来作独立性检验.fisher检验建立在超几何分布的基础上，对于单元频数小的表来说特别适合.其使用格式为fisher.test(x, y = NULL, workspace = 200000, hybrid = FALSE, control = list(), or = 1, alternative = two.sided, = TRUE, conf.level = 0.95, simulate.p.value = FALSE, B = 2000) 其中x是具有二维列联表形式的矩阵或是由因子构成的对象，y是由因子构成的对象，当x是矩阵时，y无效;workspace的输入值是一整数，其证书表示用于网络算法工作空间的大小；hybrid为逻辑变量，FALSE（默认值）表示精确计算概率，TRUE表示用混合算法计算概率；alternative表示备择，“two.sided”表示双边，“less”单边小于，“greater”单边大于；为逻辑变量，为TRUE给出区间估计；conf.level为置信水平，默认值为0.95，其余参数参见在线说明.5.3.4 McNemar检验McNemar检验虽然不是独立性检验，但它是关于列联表数据的检验，是在相同个体上的两次检验，检验二元数据的两个相关分布的频数变化的显著性.其使用格式为mcnemar.test(x, y = NULL, correct = TRUE)x,y的含义如同上面fisher.test()函数的参数一致，correct是逻辑向量，TRUE(默认值)表示计算检验统计量时连续修正，FALSE是不用修正.5.3.5 Wilcoxon秩和检验在实践中我们常常会遇到以下一些资料，如需比较患者和正常人的血铁蛋白、血铅值、不同药物的溶解时间、实验鼠发癌后的生存日数、护理效果评分等，这类数据特点：资料的总体分布类型未知或资料分布类型已知，但不符合正态分布或某些变量可能无法精确测量。对于此类数据，我们可以用非参数检验的Wilcoxon秩和检验.从两个非正态总体中所得到的两个样本之间的比较，其零假设为两个样本从同一总体中抽取的。这就是秩和检验,Wilcoxon秩和检验的使用格式为 wilcox.test(x, y = NULL, alternative = c(two.sided, less, greater), mu = 0, paired = FALSE, exact = NULL, correct = TRUE, = FALSE, conf.level = 0.95, .)其中x,y为观察数据构成的数据向量，alternative是备择假设；mu是待检参数，如中位数；paried是逻辑变量，说明x,y是否为成对数据；exact是逻辑变量，说明是否精确计算P值，当样本量较小时，此参数起作用，当样本量较大时，软件采用正态分布近似计算P值；correct是逻辑变量，说明是否对P值的计算采用连续性修正；是逻辑变量，说明是否给出相应的置信区间.5.4 应用非参数检验例5.4.1 为研究电话总机在某段时间内接到的呼叫次数是否服从Poisson分布，现收集了42个数据，如下所示，通过对数据的分析，问是否能确认在某段时间内接到的电话次数服从泊松分布（a=0.1）.表5.1 电话总机在某段时间内接到呼叫次数的频数接到的呼叫次数0 1 2 3 4 4 5 6 出现的频数7 10 12 8 3 2 0由题可知，零假设：呼叫次数的频数服从泊松分布；可以采用chisq.test()函数来检验，表格中后三组的频数分别为3,2,0，都小于5，但是Pearson 检验要求分组后，每组的频数至少要大于5，所以我们合并后三组.此时的频数为5，满足要求，下面是R程序：x-0:6;yqpzchisq.test(z,p=p)Chi-squared test for given probabilitiesdata: z X-squared = 0.5389, df = 4, p-value = 0.9696由上面的结果可以看出P0.1,因此能确认某段时间内接到的呼叫次数的频数服从泊松分布.例5.4.2 对一台设备进行寿命检验，记录10次无故障工作时间，并从小到大的次序排列如下（单位：h）420 500 920 1380 1510 1650 1760 2100 2300 2350选择合适的检验方法检验无故障工作时间是否服从的指数分布.解 由于事先未知，我们可以采用ks.test()函数来检验. xks.test(x,pexp,1/1500)One-sample Kolmogorov-Smirnov testdata: x D = 0.3015, p-value = 0.2654alternative hypothesis: two-sided由上面的结果显示p0.05,无法拒绝原假设，因此认为此设备无故障工作时间的分布服从的指数分布.例5.4.3 为了研究吸烟是否患肺癌有关，对63位肺癌患者及43名非肺癌患者（对照组）调查了其中的吸烟人数，得到22列联表，如图表5.2所示表5.2 肺癌列联表数据患肺癌 未患肺癌合计吸烟不吸烟60 323 119214合计63 43106解 列联表数据分析首先可以选择chisq.test()函数，得到吸烟跟患肺癌的相关性.R程序如下 x dim(x) chisq.test(x,correct=FALSE) Pearsons Chi-sq