R 语言基础入门
这篇《R 语言基础入门》是我在别人的博客上看到的，完成时间不会晚于 2011 年。觉得总结得相当不错，于是搬运到我的博客上。需要提到，有些链接可能已经失效。
引言
什么是 R 语言
R 语言是一个开源的数据分析环境，起初是由数位统计学家建立起来，以更好的进行统计计算和绘图，这篇 wiki 中包含了一些基本情况的介绍。由于 R 可以通过安装扩展包（Packages）而得到增强，所以其功能已经远远不限于统计分析，如果感兴趣的话可以到官方网站了解关于其功能的更多信息。
至于 R 语言名称的由来则是根据两位主要作者的首字母 (Robert Gentleman and Ross Ihaka)，但过于简短的关键词也造成在搜索引擎中很不容易找到相关的资料。不过这个专门的搜索网站可以帮到你。
为什么要学习 R 语言
- R是免费开源软件： 现在很多学术期刊都对分析软件有版权要求，而免费的分析工具可以使你在这方面不会有什么担心。另一方面，如果学术界出现一种新的数据分析方法，那么要过很长一段时间才会出现在商业软件中。但开源软件的好处就在于，很快就会有人将这种方法编写成扩展包，或者你自己就可以做这件工作。
- 命令行工作方式： 许多人喜欢类似 SPSS 菜单式的操作，这对于初学者来说很方便入门，但对于数据分析来说，命令行操作会更加的灵活，更容易进行编程和自动化处理。而且命令行操作会更容易耍酷，不是嘛，一般人看到你在狂敲一推代码后得到一个分析结果，对你投来的目光是会不一样的。
- 小巧而精悍： R 语言的安装包更小，大约不到 40M，相比其它几个大家伙它算是非常小巧精悍了。目前 R 语言非常受到专业人士欢迎，根据对数据挖掘大赛胜出者的调查可以发现，他们用的工具基本上都是 R 语言。此外，从最近几次 R 语言大会上可以了解到，咨询业、金融业、医药业都在大量的使用 R 语言，包括 google / facebook 的大公司都在用它。因此，学习 R 语言对你的职业发展一定是有帮助的。
R 语言的下载和 GUI 界面
R 语言安装包可以在官方网站下载，windows 版可直接点击这个连接 在 ubuntu 下面安装 R 则更容易，在终端里头运行如下命令即可 sudo apt-get update sudo apt-get install r-base
此外，学习 R 语言时强烈推荐安装 Rstudio 做为 R 的图形界面，关于 Rstudio 之前的博文有过简单介绍，点这里可能转到它的官方网站。
R 语言的学习方法
学习 R 并不是一件非常轻松的事情，初学者需要记住的就是：
- 亲手键入代码并理解其意义
- 在笔记里记下一些重点或心得（个人推荐 Evernote）
- 坚持练习，对手边的数据进行应用分析
- 理解背景知识，细节很重要。
哪里可以得到参考资料
- 官方网站 http://cran.csdb.cn/index.html （官方文献集中地）
- 统计之都论坛
- 人大经济论坛－R 子论坛 （免费资料也不少）
- http://library.nu/ 这是网上电子书最多的地方，其中有一个 R 语言专门书柜（也就是一个 shelves）
- 关于 R 语言的教材小结
- 笔者在 verycd 上发的一个书单
- 一个国外著名的 R 语言群博 http://www.r-bloggers.com/
- 展示 R 语言的各类绘图 http://addictedtor.free.fr/graphiques/ 本人博客里也有一些关于 R 语言的资料：xccds1977.blogspot.com (需 FQ) 如果有一些简单的入门问题，也可以在推特上 follow me twitter: @xccds
本系列博文的目的
本系列入门的目的是为初学者提供最简洁清晰的资料，以迅速入门。所针对的读者人群是那些正在大学里学习初级统计学的同学。本系列计划包括内容有：基本命令，数据操作；描述统计和绘图；重要的 R 语言函数计算；统计推断和估计；非参数统计方法；方差分析；线性回归和一般线性模型。
数据导入和描述统计
数据导入
对初学者来讲，面对一片空白的命令行窗口，第一道真正的难关也许就是数据的导入。数据导入有很多途径，例如从网页抓取、公共数据源获得、文本文件导入。为了快速入门，建议初学者采取 R 语言协同 Excel 电子表格的方法。也就是先用较为熟悉的 Excel 读取和整理你要处理的数据，然后 「粘贴」 到 R 中。
例如我们先从这个地址下载 iris.csv 演示数据，在 Excel 中打开，框选所有的样本然后 「复制」。在 R 语言中输入如下命令：
data <- read.table('clipboard',T)
Data frame 操作
在数据导入 R 语言后，会以数据框 (data frame) 的形式储存。data frame 是一种 R 的数据格式，可以将它想象成类似统计表格，每一行都代表一个样本点，而每一列则代表了样本的不同属性或特征。初学者需要掌握的基本操作方法就是 data frame 的编辑、抽取和运算。
尽管建议初学者在 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 变量的频数表和条形图：
mean(data$Species)
barplot(table(data$Species))
对于一元数值数据，绘制直方图和箱线图观察其分布是常用的方法：
hist(data$Sepal.Width)
boxplot(data$Sepal.Width)
对于二元数值数据，则可以通过散点图来观察规律
plot(data$Sepal.Width, Sepal.Length)
常用统计函数运算
在 R 语言中经常会用到函数，例如上节中讲到的求样本统计量就需要均值函数 (mean) 和标准差函数 (sd)。对于二元数值数据还用到协方差 (cov)，对于二元分类数据则可以用交叉联列表函数 (table)。下文讲述在初级统计学中最常用到的三类函数。
数据汇总函数
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 函数在建模中经常用来对样本数据进行随机的划分，一部分作为训练数据，另一部分作为检验数据。
常用的统计推断
通常一个研究项目能够获得的数据是有限的，以有限的样本特征来推断总体特征就称为统计推断。推断又可细分为区间估计和假设检验，二者虽有区别，但却是一枚硬币的两面，之间有着紧密的关联。
对总体均值进行区间估计
假设我们从总体中抽得一个样本，希望根据样本均值判断总体均值的置信区间，如下例所示：
x <- rnorm(50, mean=10, sd=5)
mean(x)-qt(0.975, 49) * sd(x) / sqrt(50)
mean(x)+qt(0.975, 49) * sd(x) / sqrt(50)
也可以直接利用 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')
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)
wilcox.test(x, mu=1)
独立性检验（联列表检验）
卡方分布有一个重要应用就是根据样本数据来检验两个分类变量的独立性，我们以 CO2 数据为例来说明 chisq.test 函数的使用，help(CO2) 可以了解更多信息。
data(CO2)
chisq.test(table(CO2$Type, CO2$Treatment))
结果显示 P 值为 0.82，因此可以认为两因子之间独立。在样本较小的情况下，还可以使用 fisher 精确检验，对应的函数是 fisher.test。
简单线性回归
线性回归可能是数据分析中最为常用的工具了，如果你认为手上的数据存在着线性定量关系，不妨先画个散点图观察一下，然后用线性回归加以分析。下面简单介绍一下如何在 R 中进行线性回归。
回归建模
我们利用 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)
模型诊断
在模型建立后会利用各种方式来检验模型的正确性，对残差进行分析是常见的方法，下面我们来生成四种用于模型诊断的图形。
par(mfrow=c(2,2))
plot(model)
par(mfrow=c(1,1))
这里左上图是残差对拟合值作图，整体呈现出一种先下降后下升的模式，显示残差中可能还存在未提炼出来的影响因素。右上图残差 QQ 图，用以观察残差是否符合正态分布。左下图是标准化残差对拟合值，用于判断模型残差是否等方差。右下图是标准化残差对杠杆值，虚线表示的 cooks 距离等高线。我们发现 31 号样本有较大的影响。
变量变换
因为 31 号样本有着高影响力，为了降低其影响，一种方法就是将变量进行开方变换来改善回归结果，从残差标准误到残差图，各项观察都说明变换是有效的。
plot(sqrt(Volume) ~ Girth, data=trees, pch=16, col='red')
model2 <- lm(sqrt(Volume) ~ Girth, data=trees)
abline(model2, lty=2)
summary(model2)
模型预测
下面根据上述模型计算预测值以及置信区间，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 回归
让我们用 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')
