R语言基础应用

1 什么是R语言
R语言是一个开源的数据分析环境，起初是由数位统计学家创建起来，以更好的进行统计计算和绘图，这篇 wiki中包含了一些基本状况的介绍。因为R能够经过安装扩展包（Packages）而获得加强，因此其功能已经远远不限于统计分析，若是感兴趣的话能够到 官方网站了解关于其功能的更多信息。

至于R语言名称的由来则是根据两位主要做者的首字母(Robert Gentleman and Ross Ihaka)，但过于简短的关键词也形成在搜索引擎中很不容易找到相关的资料。不过这个专门的 搜索网站能够帮到你。


2 为何要学习R语言
可能你想说，“我已经学会了spss/sas/stata...，为何还要去学习R呢？”以下几方面可能会吸引到你：

• R是免费开源软件:如今不少学术期刊都对分析软件有版权要求，而免费的分析工具可使你在这方面不会有什么担忧。另外一方面，若是学术界出现一种新的数据分析方法，那么要过很长一段时间才会出如今商业软件中。但开源软件的好处就在于，很快就会有人将这种方法编写成扩展包，或者你本身就能够作这件工做。
• 命令行工做方式:许多人喜欢相似SPSS菜单式的操做，这对于初学者来讲很方便入门，但对于数据分析来讲，命令行操做会更加的灵活，更容易进行编程和自动化处理。并且命令行操做会更容易耍酷，不是嘛，通常人看到你在狂敲一推代码后获得一个分析结果，对你投来的目光是会不同的。
• 小巧而精悍:R语言的安装包更小，大约不到40M，相比其它几个你们伙它算是很是小巧精悍了。目前R语言很是受到专业人士欢迎，根据对数据挖掘大赛胜出者的调查能够发现，他们用的工具基本上都是R语言。此外，从最近几回R语言大会上能够了解到，咨询业、金融业、医药业都在大量的使用R语言，包括google/facebook的大公司都在用它。所以，学习R语言对你的职业发展必定是有帮助的。

3 R语言的下载和GUI界面
R语言安装包能够在 官方网站下载，windows版可直接点击这个 链接
在ubuntu下面安装R则更容易，在终端里头运行以下命令便可
sudo apt-get update
sudo apt-get install r-base

此外，学习R语言时强烈推荐安装Rstudio作为R的图形界面，关于Rstudio以前的 博文有过简单介绍，点 这里可能转到它的官方网站。

4 R语言的学习方法
学习R并非一件很是轻松的事情，初学者须要记住的就是：

• 亲手键入代码并理解其意义
• 在笔记里记下一些重点或心得（我的推荐Evernote）
• 坚持练习，对手边的数据进行应用分析
• 理解背景知识，细节很重要。

5 哪里能够获得参考资料
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 (需)
若是有一些简单的入门问题，也能够在推特上follow me twitter:  @xccds

6 本系列博文的目的
本系列入门的目的是为初学者提供最简洁清晰的资料，以迅速入门。所针对的读者人群是那些正在大学里学习 初级统计学的同窗。本系列计划包括内容有：基本命令，数据操做；描述统计和绘图；重要的R语言函数计算；统计推断和估计；非参数统计方法；方差分析；线性回归和通常线性模型。

R语言基础入门之二：数据导入和描述统计

 

1 数据导入
对初学者来说，面对一片空白的命令行窗口，第一道真正的难关也许就是数据的导入。数据导入有不少途径，例如从网页抓取、公共数据源得到、文本文件导入。为了快速入门，建议初学者采起R语言协同Excel电子表格的方法。也就是先用较为熟悉的Excel读取和整理你要处理的数据，而后“粘贴”到R中。

例如咱们先从 这个地址下载iris.csv演示数据，在Excel中打开，框选全部的样本而后“复制”。在R语言中输入以下命令：

data=read.table('clipboard',T)

这的里read.table是R读取外部数据的经常使用命令，T表示第一行是表头信息，整个数据存在名为data的变量中。另外一种更方便的导入方法是利用Rstudio的功能，在workspace菜单选择“import dataset”也是同样的。


2 Dataframe操做
在数据导入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

3 描述统计
描述统计是一种从大量数据中压缩提取信息的工具，最经常使用的就是 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语言基础入门之三：经常使用统计函数运算

 

在R语言中常常会用到函数，例如 上节中讲到的求样本统计量就须要均值函数(mean)和标准差函数(sd)。对于二元数值数据还用到协方差(cov)，对于二元分类数据则能够用交叉联列表函数(table)。下文讲述在初级统计学中最经常使用到的三类函数。

1、数据汇总函数
咱们仍是以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数据

2、几率计算函数
若是给定一种几率分布，一般会有四类计算问题：

• 计算其几率密度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

3、抽样函数
咱们想从1到10中随机抽取5个数字，那么这样来作：首先产生一个序列，而后用sample函数进行无放回抽取。

x=1:10
sample(x,size=5)

有放回抽取则是

sample(x,size=5,replace=T)

sample函数在建模中常常用来对样本数据进行随机的划分，一部分做为训练数据，另外一部分做为检验数据。

R语言基础入门之四：经常使用的统计推断

 

一般一个研究项目可以得到的数据是有限的，以有限的样本特征来推断整体特征就称为统计推断。推断又可细分为区间估计和假设检验，两者虽有区别，但倒是一枚硬币的两面，之间有着紧密的关联。

1 对整体均值进行区间估计
假设咱们从整体中抽得一个样本，但愿根据样本均值判断整体均值的置信区间，以下例所示：

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 

2 对整体均值进行假设检验
仍是以上面的X数据做为对象，来检验整体均值是否为10

t.test(x,mu=10,alternative='two.sided')  #这里的原假设是整体均值（mu）为10，使用双侧检验，获得P值为0.17，可见P值不够小，不可以拒绝原假设。

T检验是极为经常使用的检验方法，除了单样本推断以外，t.test命令还能够实现两样本推断和配对样本推断。若是要对整体比率或整体方差进行推断，可使用 prop.test和 var.test。

3 正态分布检验
T检验的前提条件是整体服从正态分布，所以咱们有必要先检验正态性。并且在评价回归模型时，对残差也须要检验正态性。检验正态性的函数是 shapiro.test

shapiro.test(x)

结果所下：

Shapiro-Wilk normality test
data:  x
W = 0.9863, p-value = 0.8265

该检验的原假设是服从正态分布，由P值为0.82可判断不能拒绝整体服从正态的假设

4 非参数检验
若是整体不服从正态分布，那么T检验就再也不适用，此时咱们能够利用非参数方法推断中位数。 wilcoxon.test函数可实现符号秩检验。

wilcox.test(x,conf.int=T)   #指定conf.int让函数返回中位数的置信区间
wilcox.test(x,mu=1）  #指定mu让函数返回中位数为10的检验结果

5 独立性检验（联列表检验）
卡方分布有一个重要应用就是根据样本数据来检验两个分类变量的独立性，咱们以CO2数据为例来讲明 chisq.test函数的使用，help(CO2)能够了解更多信息。

data(CO2)   #读入内置的数据包，其中Type和Treatmen是其中两个分类变量。
chisq.test(table(CO2$Type,CO2$Treatment))   #使用卡方检验函数来检验这两个因子之间是否独立

结果显示P值为0.82，所以能够认为两因子之间独立。在样本较小的状况下，还可使用fisher精确检验，对应的函数是 fisher.test。

R语言基础入门之五：简单线性回归

 

线性回归多是数据分析中最为经常使用的工具了，若是你认为手上的数据存在着线性定量关系，不妨先画个散点图观察一下，而后用线性回归加以分析。下面简单介绍一下如何在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变换。

R语言基础入门之六（完）：Logistic回归

 

让咱们用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')