文章来源于:sci666 国际惯例,前面我们谈到了概率与分布,那在这基础上我们就可以分析数据的集中趋势和离散趋势。详情点击:R语言系列第二期(番外篇):R先生教你统计概率与分布 大家应该都知道统计分析分为统计描述和统计推断,我们这部分就来给大家唠叨一下统计描述。在医学科学研究中,一般需要根据研究目的先对收集的资料进行统计描述,查看一下资料的分布类型,平均水平和离散程度。往往这部分的内容可以帮助我们来确定重点研究因素和发现异常值等等。因此简单的特征分析是灰常有用的。 A. 首先我们来学习下单组数据或者说不分组数据的统计汇总。 ①下面是R计算均值、标准差、方差和中位数的程序。使用R来处理这部分数据是非常简单的。 > x<-rnorm(100) > mean(x) [1] 0.07446497 > sd(x) [1] 1.000719 > var(x) [1] 1.001438 > median(x) [1] 0.05286808 #Tips:上述例子是100个来自正态分布的数据,关于rnorm()的使用见前一个部分。因为这是一个随机数,当你再次运行这个例子的时候,得到的结果会完全不同。 常用的分位数可以用quantile()获得: > quantile(x) 0% 25% 50% 75% 100% -2.98480720 -0.62049495 0.05286808 0.67416785 2.67264258 上述数据给出的是最大值、最小值、中位数以及第一个和第三个四分位数,而后两个数值的差又称为四分位数间距(IQR)是离散程度的稳健估计,尤其适用于非正态分布。 当然有的时候你需要得到自己的定义的分位数比如十分位数: > pvec<-seq(0,1,0.1) > quantile(x,pvec) 0% 10% 20% 30% 40% -2.98480720 -1.07898557 -0.69323246 -0.53447332 -0.19438670 50% 60% 70% 80% 90% 0.05286808 0.26937874 0.59325224 0.88904176 1.34013940 100% 2.67264258 这样,不论是连续型数据还是离散型数据,我们都已经了解了他们的集中趋势和离散趋势的描述方式。 ② 不过当出现缺失值,情况就不同了。 我们以ISwR包中的juul数据集为例子, > library(ISwR) > head(juul) age menarche sex igf1 tanner testvol 1 NA NA NA 90 NA NA 2 NA NA NA 88 NA NA 3 NA NA NA 164 NA NA 4 NA NA NA 166 NA NA 5 NA NA NA 131 NA NA 6 0.17 NA 1 101 1 NA #Tips: 如果在加载ISwR的过程中提示:Error in library(ISwR) : 不存在叫‘ISwR’这个名字的程辑包。需要安装ISwR包。> install.packages(“ISwR”) Juul数据集报告的是健康人群,主要是儿童血清IGF-1(类胰岛素生长因子)的水平。里面含有6个变量,这里我们暂时只用到IGF-1的值。 然而当我们试着计算igf1均值的时候发现了一个问题: > attach(juul) > mean(igf1) [1] NA R不能跳过缺失值,除非明确要求跳过。有缺失值的向量的均值也是缺失的。那要如何解决呢?可以用na.rm(not available,remove)参数来要求移除缺失值: > mean(igf1,na.rm=T) [1] 340.168 #Tips:这里计算均值是先排除掉缺失值,在计算均值,所以分母是有效数字的个数不包括缺失值的数量。 不过,有一个不好的事情是我们没法用length()函数来算有多少个有效数值,因为它不包含na.rm的参数,只有一个参数。所以在这里我们也就不能用它计算igf1中非缺失观测的个数。 > length(igf1) [1] 1339 > length(igf1,na.rm=T) Error in length(igf1, na.rm = T) : 2个参数给‘length’,但它只需要1个 但是我们可以用以下方法查看有多少个非缺失观测数: > sum(!is.na(igf1)) [1] 1018 #Tips:上述结果利用了逻辑向量作为数字计算,即TRUE转化成1,FALSE转化成0,那么如果出现缺失值,则会判断is.na为TRUE,而!则把它转化成FALSE,因此,最后所有的缺失值判断结果都是FALSE,所有的有效值都转化成TRUE,因此最后的合计就是所有TRUE的数量,也就是有效数值的数量。 ③ 同时我们可以用汇总函数查看数据情况。 summary()函数就是对数据进行汇总报告的一个函数: > summary(igf1) Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s 25.0 202.2 313.5 340.2 462.8 915.0 321 对整个数据集进行汇总: > summary(juul) 同时我们从原始数据中可以发现menarche(月经初潮),sex,tanner(青春期分期)这三个变量在原始数据中被编成数值变量,实际上他们是分类变量。那么我们可以改进一下: > juul$sex<-factor(juul$sex,labels=c(“M”,”F”)) > juul$menarche<-factor(juul$menarche,labels=c(“No”,”Yes”)) > juul$tanner<-factor(juul$tanner,labels=c(“I”,”II”,”III”,”IV”,”V”)) > juul > summary(juul) #Tips:我们可以观察下因子变量的变化,另外summary()对因子的总结就不再是原来的各个分位点的值了,而是各水平的频数(数量)。 另外,我们有其他方法,把这些用数值表示的分类变量转换成以适当水平填充的因子变量,然后转化后的变量重新放回数据框中,代替原来的三个变量。就是我们之前《R语言系列2(1)》接触的transform()函数(或within()函数): > juul<-transform(juul, +sex=factor(sex,labels=c(“M”,”F”)), +menarche=factor(menarche,labels=c(“No”,”Yes”)), +tanner=factor(tanner,labels=c(“I”,”II”,”III”,”IV”,”V”))) B. 分布图形展示 01 直方图 我们可以通过绘制直方图得到对数据分布的一个适当的印象: > hist(x) 通过在hist()函数中指定breaks=n,你就可以在直方图中得到n个柱,算法会自动为我们创建合适的分割点。我们也可以通过指定breaks成一个向量来实现对区间分割的完全控制。比如我们对年龄分组:0-4岁,5-9岁,10-15岁,10-15岁,16岁,17岁,18-19岁,20-24岁,25-59岁和60-79岁: > mid.age<-c(2.5,7.5,13,16.5,17.5,19,22.5,44.5,70.5) > acc.count<-c(28,46,58,20,31,64,149,316,103) > age.acc<-rep(mid.age,acc.count) > brk<-c(0,5,10,16,17,18,20,25,60,80) > hist(age.acc,breaks=brk) #Tips:这里前三行不是原始数据,是通过原始数据分组的中位值和其数量做成的一组对应的数据,即0-4岁有28个,5-9岁有46个人。另外,每个柱子的面积与组内观测数成正比。Y轴是密度单位,直方图的总面积是1.,而如果想要直方图的列高度代表每个区间的频数而不是密度,那么你需要用参数freq=T指定。 02 累积经验分布 累积经验分布函数定义为小于等于x的数据占总数据的比例。 > n<-length(x) > plot(sort(x),(1:n)/n,type=”s”,ylim=c(0,1)) #Tips:x是前一部分的100个随机数组成的向量。type=”s”画出来的是阶梯函数。 03 Q-Q图 我们画累计经验分布函数的目的之一就是观察数据是否呈正态分布。但是我们有更加简便的方法,那就是qqnorm()函数,你只需要写: > qqnorm(x) #Tips:相对于上一个图形,这里的Q-Q图需要注意的是观测的结果画在了y轴上,x轴表示理论分位数,y轴表示样本分位数。对于Q-Q图而言,如果像如图这种呈一种斜率接近1的直线的话,证明数据正态性很好。而如果两侧比中间陡峭,那分布就是重尾的。数据分布两端偏多。 04 箱式图 箱式图或者叫箱须图也是一个分布的图形概括。我们拿ISwR里的IgM数据集作为作图的例子。 #Tips:中间的箱是‘中枢’(上下边框分别是上下四分位数,中间加粗的黑线是中位数),上下的线则表示落在离最近箱边1.5倍箱子尺寸距离的最大或最小值。其它更远的观测则被视为极端值全部以点的形式被画出来。 #Tips:par()是个很重要的高级画图函数,这里两个并列图的布局使用mfrow参数指定,它的全称是multiframe,rowwise,1*2 layout。单个图被组织在一行两列上。那么mfcol参数就是用来绘制按列排的图形。在2*2的布局里,他们俩的区别就是第二幅图在右上角还是左下角。记住,最后一定要布局参数调回到c(1,1),除非你还想这样画出并列在一起的图。关于par()函数的应用知识,之后会为大家揭晓。 这部分就是不分组数据的描述汇总,下个部分是分组数据的汇总统计,敬请期待。 参考资料: End 往期推文 R统计分析与绘图系列 R语言系列第一期(番外篇 ):R的6种对象—向量、矩阵、数组、因子、列表、数据框 生信神器系列 【神器分享】自从用了这个神器,大规模RNA-seq数据挖掘我也可以 DAVID&Metascape:专注于基因功能注释和富集通路分析的网站 SCI写作系列 玩转SCI:短平直快查询 IF & JCR分区&审稿周期&国产发文 GEO数据挖掘系列 后台回复“入群”,百味科研交流群等你加入 多点好看,少点脱发 |
|