编程语言及工具
R是用于统计分析、绘图的语言和操作环境。R是属于GNU系统的一个自由、免费、源代码开放的软件,它是一个用于统计计算和统计制图的优秀工具。
R语言是一个开源的数据分析环境,起初是由数位统计学家建立起来,以更好的进行统计计算和绘图,这篇wiki中包含了一些基本情况的介绍。由于R可以通过安装扩展包(Packages)而得到增强,所以其功能已经远远不限于统计分析。
R作为一种统计分析软件,是集统计分析与图形显示于一体的。它可以运行于UNIX,Windows和Macintosh的操作系统上,而且嵌入了一个非常方便实用的帮助系统,相比于其他统计分析软件,R还有以下特点:1.R是自由软件。这意味着它是完全免费,开放源代码的。可以在它的网站及其镜像中下载任何有关的安装程序、源代码、程序包及其源代码、文档资料。标准的安装文件身自身就带有许多模块和内嵌统计函数,安装好后可以直接实现许多常用的统计功能。
2.R是一种可编程的语言。作为一个开放的统计编程环境,语法通俗易懂,很容易学会和掌握语言的语法。而且学会之后,我们可以编制自己的函数来扩展现有的语言。这也就是为什么它的更新速度比一般统计软件,如,SPSS,SAS等快得多。大多数最新的统计方法和技术都可以在R中直接得到。
3. 所有R的函数和数据集是保存在程序包里面的。只有当一个包被载入时,它的内容才可以被访问。一些常用、基本的程序包已经被收入了标准安装文件中,随着新的统计分析方法的出现,标准安装文件中所包含的程序包也随着版本的更新而不断变化。在另外版安装文件中,已经包含的程序包有:base一R的基础模块、mle一极大似然估计模块、ts一时间序列分析模块、mva一多元统计分析模块、survival一生存分析模块等等。
4.R具有很强的互动性。除了图形输出是在另外的窗口处,它的输入输出窗口都是在同一个窗口进行的,输入语法中如果出现错误会马上在窗口口中得到提示,对以前输入过的命令有记忆功能,可以随时再现、编辑修改以满足用户的需要。输出的图形可以直接保存为JPG,BMP,PNG等图片格式,还可以直接保存为PDF文件。另外,和其他编程语言和数据库之间有很好的接口。[2] 5.如果加入R的帮助邮件列表一,每天都可能会收到几十份关于R的邮件资讯。可以和全球一流的统计计算方面的专家讨论各种问题,可以说是全世界最大、最前沿的统计学家思维的聚集地。
R是基于S语言的一个GNU项目,所以也可以当作S语言的一种实现,通常用S语言编写的代码都可以不作修改的在R环境下运行。 R的语法是来自Scheme。R的使用与S-PLUS有很多类似之处,这两种语言有一定的兼容性。S-PLUS的使用手册,只要稍加修改就可作为R的使用手册。所以有人说:R,是S-PLUS的一个“克隆”。但是请不要忘了:R是免费的(R is free)。R语言源代码托管在github,具体地址可以看参考资料。
R语言的下载可以通过CRAN的镜像来查找。
R语言有域名为.cn的下载地址,有六个,其中两个由Datagurn,由中国科学技术大学提供的。R语言Windows版,其中由两个下载地点是Datagurn和USTC提供的。
在继续学习本教程之前,您应该基本了解计算机编程术语。 对任何编程语言的基本理解将帮助您理解R语言编程概念,并在学习轨道上快速移动
本教程是为期待使用R编程开发统计软件的软件程序员,统计学家和数据挖掘者设计的。 如果你试图理解R编程语言作为一个初学者,本教程将给你足够的了解语言的几乎所有的概念,从那里你可以把自己的更高水平的专业知识。
可能你想说,“我已经学会了spss/sas/stata.。。,为什么还要去学习R呢?”如下几方面可能会吸引到你:
R是免费开源软件:现在很多学术期刊都对分析软件有版权要求,而免费的分析工具可以使你在这方面不会有什么担心。另一方面,如果学术界出现一种新的数据分析方法,那么要过很长一段时间才会出现在商业软件中。但开源软件的好处就在于,很快就会有人将这种方法编写成扩展包,或者你自己就可以做这件工作。
命令行工作方式:许多人喜欢类似SPSS菜单式的操作,这对于初学者来说很方便入门,但对于数据分析来说,命令行操作会更加的灵活,更容易进行编程和自动化处理。而且命令行操作会更容易耍酷,不是嘛,一般人看到你在狂敲一推代码后得到一个分析结果,对你投来的目光是会不一样的。
小巧而精悍:R语言的安装包更小,大约不到40M,相比其它几个大家伙它算是非常小巧精悍了。目前R语言非常受到专业人士欢迎,根据对数据挖掘大赛胜出者的调查可以发现,他们用的工具基本上都是R语言。此外,从最近几次R语言大会上可以了解到,咨询业、金融业、医药业都在大量的使用R语言,包括google/facebook的大公司都在用它。因此,学习R语言对你的职业发展一定是有帮助的。
R语言安装包可以在官方网站下载,windows版可直接点击这个连接
在ubuntu下面安装R则更容易,在终端里头运行如下命令即可
sudo apt-get update
sudo apt-get install r-base
此外,学习R语言时强烈推荐安装Rstudio做为R的图形界面,关于Rstudio之前的博文有过简单介绍,点这里可能转到它的官方网站。
学习R并不是一件非常轻松的事情,初学者需要记住的就是:
亲手键入代码并理解其意义
在笔记里记下一些重点或心得(个人推荐Evernote)
坚持练习,对手边的数据进行应用分析
理解背景知识,细节很重要。
1.官方网站 http://cran.csdb.cn/index.html (官方文献集中地)
2.统计之都论坛
3.人大经济论坛-R子论坛 (免费资料也不少)
4.http://library.nu/ 这是网上电子书最多的地方,其中有一个R语言专门书柜(也就是一个shelves)
5.关于R语言的教材小结
6.笔者在verycd上发的一个书单
7.一个国外著名的R语言群博 http://www.r-bloggers.com/
8.展示R语言的各类绘图 http://addictedtor.free.fr/graphiques/
本人博客里也有一些关于R语言的资料:xccds1977.blogspot.com (需FQ)
如果有一些简单的入门问题,也可以在推特上follow me twitter: @xccds
本系列入门的目的是为初学者提供最简洁清晰的资料,以迅速入门。所针对的读者人群是那些正在大学里学习初级统计学的同学。本系列计划包括内容有:基本命令,数据操作;描述统计和绘图;重要的R语言函数计算;统计推断和估计;非参数统计方法;方差分析;线性回归和一般线性模型。
#p##e#
对初学者来讲,面对一片空白的命令行窗口,第一道真正的难关也许就是数据的导入。数据导入有很多途径,例如从网页抓取、公共数据源获得、文本文件导入。为了快速入门,建议初学者采取R语言协同Excel电子表格的方法。也就是先用较为熟悉的Excel读取和整理你要处理的数据,然后“粘贴”到R中。
例如我们先从这个地址下载iris.csv演示数据,在Excel中打开,框选所有的样本然后“复制”。在R语言中输入如下命令:
data=read.table(‘clipboard’,T)
这的里read.table是R读取外部数据的常用命令,T表示第一行是表头信息,整个数据存在名为data的变量中。另一种更方便的导入方法是利用Rstudio的功能,在workspace菜单选择“import dataset”也是一样的。
在数据导入R语言后,会以数据框(dataframe)的形式储存。dataframe是一种R的数据格式,可以将它想象成类似统计表格,每一行都代表一个样本点,而每一列则代表了样本的不同属性或特征。初学者需要掌握的基本操作方法就是dataframe的编辑、抽取和运算。
尽管建议初学者在Excel中就把数据处理好,但有时候还是需要在R中对数据进行编辑,下面的命令可以让你有机会修改数据并存入到新的变量newdata中:
newdata=edit(data)
另一种情况就是我们可能只关注数据的一部分,例如从原数据中抽取第20到30号样本的Sepal.Width变量数据,因为Sepal.Width变量是第2个变量,所以此时键入下面的命令即可:
newdata=data[20:30,2]
如果需要抽取所有数据的Sepal.Width变量,那么下面两个命令是等价的:
newdata=data[,2]
newdata=data$Sepal.Width
第三种情况是需要对数据进行一些运算,例如需要将所有样本的Sepal.Width变量都放大10倍,我们先将原数据进行一个复制,再用$符号来提取运算对象即可:
newdata=data
newdata$Sepal.Width=newdata$Sepal.Width*10
描述统计是一种从大量数据中压缩提取信息的工具,最常用的就是summary命令,运行summary(data)得到结果如下:对于数值变量计算了五个分位点和均值,对于分类变量则计算了频数。
也可以单独计算Sepal.Width变量的平均值和标准差
mean(data$Sepal.Width)
sd(data$Sepal.Width)
计算分类数据Species变量的频数表和条形图
table(data$Species)
barplot(table(data$Species))
对于一元数值数据,绘制直方图和箱线图观察其分布是常用的方法:
hist(data$Sepal.Width)
boxplot(data$Sepal.Width)
对于二元数值数据,则可以通过散点图来观察规律
plot(data$Sepal.Width,Sepal.Length)
如果需要保存绘图结果,建议使用Rstudio中的plot菜单命令,选择save plot as image
在R语言中经常会用到函数,例如上节中讲到的求样本统计量就需要均值函数(mean)和标准差函数(sd)。对于二元数值数据还用到协方差(cov),对于二元分类数据则可以用交叉联列表函数(table)。下文讲述在初级统计学中最常用到的三类函数。
我们还是以R中自带的iris数据为例,输入head(iris)你可以获得数据的前6个样本及对应的5个变量。取出最后两列数据作为讲解的对象:Species表示花的种类,Petal.Width表示花瓣宽度
data=iris[,c(4,5)]
下一步我们想计算不同种类花瓣的平均宽度,可以使用tapply函数,在计算前先用attach命令将data这个数据框解包以方便直接操作其变量,而不需再用$符号。
attach(data)
tapply(X=Petal.Width,INDEX=Species,FUN=mean)
结果如下
setosa versicolor virginica
0.246 1.326 2.026
和tapply类似的还有sapply函数,在进一步讲解前初学者还需搞清楚两种数据表现方式,即stack(堆叠数据)和unstack(非堆叠数据),上面的data就是一个堆叠数据,每一行表示一个样本。而非堆叠数据可以根据unstack函数转换而来
data.unstack=unstack(data)
head(data.unstack)
你应该明白这二者之间的区别了,如果要对非堆叠数据计算不同种类花瓣的平均宽度,可以利用如下函数。
sapply(data.unstack,FUN=mean)
结果是一样的,也就是说tapply对应于stack数据,而sapply对应于unstack数据
如果给定一种概率分布,通常会有四类计算问题:
计算其概率密度density (d)
计算其概率分布probability(p)
计算其百分位数quantile (q)
随机数模拟random (r)
记住上面四类计算对应的英文首字母,再对照下表就很容易计算各种概率问题了。
举例来讲,我们求标准正态分布曲线下小于1的面积p(x《1),正态分布是norm,而分布函数是p,那么使用pnorm(1)就得出了结果0.84;若计算扔10次硬币实验中有3次正面向上的概率,类似的dbinom(x=3,size=10,prob=0.5)得出0.11
我们想从1到10中随机抽取5个数字,那么这样来做:首先产生一个序列,然后用sample函数进行无放回抽取。
x=1:10
sample(x,size=5)
有放回抽取则是
sample(x,size=5,replace=T)
sample函数在建模中经常用来对样本数据进行随机的划分,一部分作为训练数据,另一部分作为检验数据。
#p##e#
通常一个研究项目能够获得的数据是有限的,以有限的样本特征来推断总体特征就称为统计推断。推断又可细分为区间估计和假设检验,二者虽有区别,但却是一枚硬币的两面,之间有着紧密的关联。
假设我们从总体中抽得一个样本,希望根据样本均值判断总体均值的置信区间,如下例所示:
x=rnorm(50,mean=10,sd=5) #随机生成50个均值为10,标准差为5的随机数为作为研究对象
mean(x)-qt(0.975,49)*sd(x)/sqrt(50) #根据统计学区间估计公式,得到95%置信度下的区间下界
mean(x)+qt(0.975,49)*sd(x)/sqrt(50) #95%置信度下的区间上界
也可以直接利用R语言内置函数t.test
t.test(x,conf.level=0.95)
从如下结果可得95%置信区间为(9.56,12.36)
One Sample t-test
data: x
t = 15.7301, df = 49, p-value 《 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
9.563346 12.364729
sample estimates:
mean of x
10.96404
还是以上面的X数据作为对象,来检验总体均值是否为10
t.test(x,mu=10,alternative=‘two.sided’) #这里的原假设是总体均值(mu)为10,使用双侧检验,得到P值为0.17,可见P值不够小,不能够拒绝原假设。
T检验是极为常用的检验方法,除了单样本推断之外,t.test命令还可以实现两样本推断和配对样本推断。如果要对总体比率或总体方差进行推断,可以使用prop.test和var.test。
T检验的前提条件是总体服从正态分布,因此我们有必要先检验正态性。而且在评价回归模型时,对残差也需要检验正态性。检验正态性的函数是shapiro.test
shapiro.test(x)
结果所下:
Shapiro-Wilk normality test
data: x
W = 0.9863, p-value = 0.8265
该检验的原假设是服从正态分布,由P值为0.82可判断不能拒绝总体服从正态的假设
如果总体不服从正态分布,那么T检验就不再适用,此时我们可以利用非参数方法推断中位数。wilcoxon.test函数可实现符号秩检验。
wilcox.test(x,conf.int=T) #指定conf.int让函数返回中位数的置信区间
wilcox.test(x,mu=1) #指定mu让函数返回中位数为10的检验结果
卡方分布有一个重要应用就是根据样本数据来检验两个分类变量的独立性,我们以CO2数据为例来说明chisq.test函数的使用,help(CO2)可以了解更多信息。
data(CO2) #读入内置的数据包,其中Type和Treatmen是其中两个分类变量。
chisq.test(table(CO2$Type,CO2$Treatment)) #使用卡方检验函数来检验这两个因子之间是否独立
结果显示P值为0.82,因此可以认为两因子之间独立。在样本较小的情况下,还可以使用fisher精确检验,对应的函数是fisher.test。
线性回归可能是数据分析中最为常用的工具了,如果你认为手上的数据存在着线性定量关系,不妨先画个散点图观察一下,然后用线性回归加以分析。下面简单介绍一下如何在R中进行线性回归。
1 回归建模
我们利用R语言中内置的trees数据,其中包含了Volume(体积)、Girth(树围)、Height(树高)这三个变量,我们希望以体积为因变量,树围为自变量进行线性回归。
plot(Volume~Girth,data=trees,pch=16,col=‘red’)
model=lm(Volume~Girth,data=trees)
abline(model,lty=2)
summary(model)
首先绘制了两变量的散点图,然后用lm函数建立线性回归模型,并将回归直线加在原图上,最后用summary将模型结果进行了展示,从变量P值和F统计量可得回归模型是显著的。但截距项不应该为负数,所以也可以用下面方法将截距强制为0。
model2=lm(Volume~Girth-1,data=trees)
2 模型诊断
在模型建立后会利用各种方式来检验模型的正确性,对残差进行分析是常见的方法,下面我们来生成四种用于模型诊断的图形。
par(mfrow=c(2,2))
plot(model)
par(mfrow=c(1,1))
这里左上图是残差对拟合值作图,整体呈现出一种先下降后下升的模式,显示残差中可能还存在未提炼出来的影响因素。右上图残差QQ图,用以观察残差是否符合正态分布。左下图是标准化残差对拟合值,用于判断模型残差是否等方差。右下图是标准化残差对杠杆值,虚线表示的cooks距离等高线。我们发现31号样本有较大的影响。
3 变量变换
因为31号样本有着高影响力,为了降低其影响,一种方法就是将变量进行开方变换来改善回归结果,从残差标准误到残差图,各项观察都说明变换是有效的。
plot(sqrt(Volume)~Girth,data=trees,pch=16,col=‘red’)
model2=lm(sqrt(Volume)~Girth,data=trees)
abline(model2,lty=2)
summary(model2)
4 模型预测
下面根据上述模型计算预测值以及置信区间,predict函数可以获得模型的预测值,加入参数可以得到预测区间
plot(sqrt(Volume)~Girth,data=trees,pch=16,col=‘red’)
model2=lm(sqrt(Volume)~Girth,data=trees)
data.pre=data.frame(predict(model2,interval=‘prediction’))
lines(data.pre$lwr~trees$Girth,col=‘blue’,lty=2)
lines(data.pre$upr~trees$Girth,col=‘blue’,lty=2)
我们还可以将树围和树高都加入到模型中去,进行多元回归。如果要考虑的变量很多,可以用step函数进行变量筛选,它是以AIC作为评价指标来判断一个变量是否应该加入模型,建议使用这种自动判断函数时要谨慎。对于嵌套模型,还可以使用anova建立方差分析表来比较模型。对于变量变换的形式,则可以使用MASS扩展包中的boxcox函数来进行COX变换。
让我们用logistic回归来结束本系列的内容吧,本文用例来自于John Maindonald所著的《Data Analysis and Graphics Using R》一书,其中所用的数据集是anesthetic,数据集来自于一组医学数据,其中变量conc表示麻醉剂的用量,move则表示手术病人是否有所移动,而我们用nomove做为因变量,因为研究的重点在于conc的增加是否会使nomove的概率增加。
首先载入数据集并读取部分文件,为了观察两个变量之间关系,我们可以利cdplot函数来绘制条件密度图。
library(DAAG)
head(anesthetic)
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)
上面的方法是使用原始的0-1数据进行建模,即每一行数据均表示一个个体,另一种是使用汇总数据进行建模,先将原始数据按下面步骤进行汇总
anestot=aggregate(anesthetic[,c(‘move’,‘nomove’)],by=list(conc=anesthetic$conc),FUN=sum)
anestot$conc=as.numeric(as.character(anestot$conc))
anestot$total=apply(anestot[,c(‘move’,‘nomove’)],1,sum)
anestot$prop=anestot$nomove/anestot$total
得到汇总数据anestot如下所示
conc move nomove total prop
1 0.8 6 1 7 0.1428571
2 1.0 4 1 5 0.2000000
3 1.2 2 4 6 0.6666667
4 1.4 2 4 6 0.6666667
5 1.6 0 4 4 1.0000000
6 2.5 0 2 2 1.0000000
对于汇总数据,有两种方法可以得到同样的结果,一种是将两种结果的向量合并做为因变量,如anes2模型。另一种是将比率做为因变量,总量做为权重进行建模,如anes3模型。这两种建模结果是一样的。
anes2=glm(cbind(nomove,move)~conc,family=binomial(link=‘logit’),data=anestot)
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’)
全部0条评论
快来发表一下你的评论吧 !