在通常的情况下，我习惯写一些比较有意思的感悟，这样喜欢看的受众也广泛一些。但这次恐怕要让部分朋友失望了，我想总结下统计软件R的学习精要。R软件做为使用最为广泛的开源统计软件，有着其独特的设计和吸引人的魅力。但这次与大家谈的，更多为R软件的入门精要，期望对统计感兴趣或者工作中用到统计的朋友能够快速上手，体会到用R来完成数据分析工作的乐趣。
整篇文章分为两个部分。在基础部分会向大家介绍R的一些基本概念，包括：R是什么、相关资源和安装、数据结构、实用函数、输入输出、作图、脚本编程。而在实战部分会向大家介绍如何用R来进行一些常用的数据分析工作：相关系数、回归分析、分类判别、聚类分析、时间序列、主成分分析、T校验。
基础部分
- R是什么？英文的解释为：The R Project for Statistical Computing，它既是一种统计语言，也是一种统计软件。
- 为什么用R？总结了如下四点原因：
a) 多领域的统计资源：目前在R网站上约有2400个程序包，涵盖了基础统计学、社会学、经济学、生态学、空间分析、系统发育分析、生物信息学等诸多方面。
b) 跨平台：R可在多种操作系统下运行，如Windows、MacOS、多种Linux和UNIX等。
c) 命令行 和 脚本编程:R即时解释，输入命令，即可获得相应的结果。也支持相对复杂功能的脚本编程，处理输入输出文件等。
d) 开源，庞大的社区支持:可以免费、持续的获取最新的算法支持。
3. R的资源
a) 官方网站：http://www.r-project.org/
b) 中文网站：http://www.r-user.org/
c) 统计之都： http://cos.name/
4. 安装流程：共分为四步
a) 下载源代码版本（Source Code for all Platforms）
b) 配置安装路径：./configure –prefix=/home/work/tools/R/
c) 编译：make
d) 安装：make install
5. 数据结构：R与常见的编程语言一样，具有数值、字符、复数、逻辑四种基本类型。同时，为了方便的应对统计和分析处理，提供了矩阵、数据框等复合类型来表示二维表格。常用的复合类型如下，可以用函数mode()查看一个变量的类型。
a) 向量：一个变量
b) 因子：一个分类变量
c) 数组：你懂的
d) 矩阵：相当于2维数组
e) 数据框：相当于含有多种数据类型元素的矩阵，而矩阵元素必须是同一类型
f) 时间序列：附加时间、频率等额外的数据
g) 列表：相当于结构体，里面什么都可以有。R很多函数的返回值都是这个类型，可以用 列表变量名$内部变量名，来取出该变量内容。
6. R的一些实用函数
a) ？ lm：查看该函数的说明文档，用于已知函数名时的精准查询
b) apropos(“lm”)： 模糊查询所有相关函数
c) attributes(x)：查看列表变量x所含有的变量名，相应变量结果可以用‘x$变量名’得到。通常用来查看函数返回结果中包含哪些可用内容。
d) summary(x)：查看列表变量x的汇总内容，通常用于查看函数的返回结果
e) setwd()、getwd()：设置工作路径
f) library(MASS)：加载一些非默认的函数库。R将所有的函数库分为三个层次，最为通用+常用的自动加载到内存，可以直接在程序中调用。而次要不太常用的放在硬盘上，需要使用的时候需要用library函数现加载到内存中。最专用、只能满足小众需求的，比如DNA生物库，是不在默认安装包里的，需要使用时可以用install.packages()来下载安装。
7. 输入输出：最常用的为读入输出变量，和读入输出表格。
a) 直接写变量名：打印到屏幕
b) 从文件输入
i. X=scan(“file.txt”,nrow=3,ncol=5,byrow=T) #将文件内容读入变量
ii. X=read.table(“file.txt”, header=T, sep=“\t”) #读入表格，数据有表头，以”\t”作为列分隔符
iii. X=read.delim(“file.txt”, header=F, sep=”\t”) #读入表格，但每行的列数可以不同
c) 写入文件
i. write(X，file=“output”,append=F)
ii. write.table(X，file=“output”，append=T，row.names=F，col.names=F，sep=”\t”)
8. 作图：R将作图分为高级作图函数和低级作图函数，其中高级函数提供的是标准、成型的图形服务，如散点图、柱状图等。而低级函数提供基本的绘画功能，高级作图函数是调用低级函数来实现的。如果喜欢，也可以用低级函数编写自己想要的图形。
9. 最后再进入实战之前，介绍下写程序所需的最后一部分知识，即程序的控制和执行。R软件的条件语句和循环控制语句与C语言一致，使用很方便。
a) R的条件判断：If（condition） statement1 else statment2
b) R的循环语句：for（i in 1:n） { statement } while（condition） { statement }
c) 保存脚本文件为 program.R，执行脚本： R -f “program.R”
实战部分
上一部分为R的基本使用知识，但只有这些我们并没有必要去学习R，任何一种编程语言均能实现这些基本的运算，在大部分方面甚至可以实现的更好。那么R吸引我们的魅力就要在这一部分向大家展现，即实用的数学、统计、算法函数。
1. 回归分析
回归是寻找自变量和因变量之间的逻辑关系的方法。R提供的常用回归分析有：线性回归lm()、非线性回归nls()、Logistic回归glm()。线性回归和非线性回归，需要给出一个设计好的公式，而Logistic回归则要求回归变量为分类变量。下面以一个实例来阐述R软件的处理流程：
a) 读入数据文件，为每行有6列数据的表格
b) 对每一行数据进行线性回归，取得结果中的斜率项
c) 将这行数字、数字之和、线性拟合的数据输出到结果文件中
| #共六个数据项
x<-c(1:6) #读入数据文件 data = read.delim(“input”,header=FALSE,sep=”\t”) for(i in 1:length(data[,1])) { #取出一列数据 y<-as.numeric(data[i,]) cl<-data.frame(x,y) #最简单的线性公式 result<-lm(y~x, data=cl) #复杂自设计的非线性公式 #result1<-nls(y~a+(0.49-a)*exp(-b*(x-8)), data=cl, start=list(a=0.1, b=0.01)) #Logistic回归，回归结果变量为定性变量 #result2=glm(Opinion~Age+Income, family=”binomial”, data=cl) #取出结果中的线性斜率的系数项 rate=as.matrix(coef(result)) #打印结果，附加求和项和系数项 result=c(data[i,], sum(data[i,]),rate[2,1]) write.table(result,file=”output”, append=TRUE,row.names=FALSE,col.names=FALSE,sep=”\t”) } |
2. 相关系数
线性相关系数为确定两个变量是否数据相关的产用方法，R实现了三种算法：pearson、kendall、spearman。输入文件为两列，第一列的title是first，第二列的title是second。
| x=read.table(“input”, head=T)
#三种线性相关系数均支持, 默认为Pearson系数 result=cor(x$first, x$second) #result=cor(x$first, x$second,method=”kendall”) #result=cor(x$first, x$second,method=”spearman”)
Result |
3. 聚类分析
聚类是自动根据数据特征进行分类的方法，其特点时事先并不知道分类依据，完全根据自变量数据表现分类。常用方法除了大家都知道的KMeans聚类之外，还有分层聚类：将每步都将聚类最近的点进行合并，直到合并成只有一个节点，根据需要取某一步为最终结果。
| x=read.table(“cities1.txt”, head=T, sep=”,”)
#得到数量为4的聚类结果 result=kmeans(x[,-1], 4, iter.max = 10) attributes(result) result$cluster result$centers #分层聚类的用法，计算距离的方案有: 离差平方和”ward”, 最短距离法”single”, 最长距离法”complete”, 类平均法”average”, “mcquitty”, 中间距离法”median” or 重心法”centroid”. hc=hclust(dist(x[,-1]),”ave”) #在linux下执行，会将画图结果保存到当前目录下，为PDF文件 plot(hc, label=x[,1]) |
4. 判别分析
分类的方案有很多，这里仅仅取一个用线性模型判别的方法。需要注意，分类与聚类不同，需要事先给出一个训练集，同时需要预留出测试集来对分类效果进行测试。
| #该数据分析函数不在R默认包中，需要引入
library(MASS) #数据集使用内置数据集 iris iris #设定随机种子 set.seed(5) #从150个数据记录中取一半作为训练集，一半作为测试集 trainset = c(sample(1:50, 25), sample(51:100,25), sample(101:150, 25)) #使用线性判别模型进行判别, 将Species的列作为判别结果，其余列作为输入 model = lda(Species~., data=iris, subset=trainset) #使用训练的模型对测试集进行判别 result=predict(model, iris[-trainset,]) #打印判别的结果 result$class |
5. 时间序列分析
对时间序列分析有两个常用的目标：第一个是拆解，将数据中的趋势变量，周期变量拆解出来，以便对数据的趋势有更清晰的认识。第二个是预测，预估曲线未来的发展趋势。
R中的曲线拆解将曲线分为了三个部分：趋势分量、周期分量、剩余部分。而我的个人观点，将任何一条时间序列曲线拆成如下四个部分趋势+季节性+例外（异常）+随机波动更加合理一些。另：ts函数是将普通的数字串组织成时间序列的格式。
| x=scan(“input”)
#转换为时间序列的格式 #start参数或者为数字，或者为年月的二元组 charge=ts(x,frequency=7,start=1) #时间序列的分解函数 a=stl(charge, “period”) result=a$time.series #将分解结果的三部分，输出到文件中 第一列为季节分量，第二列为趋势分量，第三列为剩下的部分 write.table(result,file=”output”, append=FALSE,row.names=FALSE,col.names=FALSE,sep=”\t”) #ts.plot(a$time.series) |
曲线趋势预估，最常用的方法有移动平均法、指数平滑法、ARIMA方法，其中移动平均法非常简单，这里仅提供R软件关于指数平滑法和ARIMA方法的使用程序。
| x=scan(“input”)
charge=ts(x,frequency=7,start=1) b = HoltWinter(charge, beta=0) #输出拟合的曲线 b$fitted #对未来三十天作预测 predict(b, n.ahead=30) |
ARIMA预测模型，需要自行根据自相关图和偏相关图来确定P与Q，函数会自动训练系数值。
| x=scan(“input”)
charge=ts(x,frequency=7,start=1) #查看曲线的自相关图和偏自相关图 #acf(charge) pacf(charge) #arima模型拟合 #第二项参数为p,d,q 分别为自回归模型项数,差分等级,移动平均项数,seasonal的三个参数为季节性分量的配置 w=arima(charge, c(0,1,1),seasonal=list(order=(1,2,1), period=7) #模型的拟合系数 w$coef #使用模型预测 result=predict(w, n.ahead=30) write.table(result,file=”output”, append=FALSE,row.names=FALSE,col.names=FALSE,sep=”\t”) |
6. 主成分分析
确定一个多维数据的主成分，将多维数据进行降维又能最大化的保留原始信息。
| #使用R软件自带的iris数据集
iris #进行主成分分析 #算法使用相关矩阵还是协方差矩阵来计算 prin_result=princomp(iris[c(2,3,4)], cor=T) prin_result$sdev prin_result$loadings #因子分析，只选择一个因子 factor_result=factanal(iris[c(2,3,4)], fact=1) #载荷结果 factor_result$loadings |
主成分分析结果，有两部分信息是最为重要的。第一是每个新的主成分所占原数据信息的比重，如下面comp.1的比重是1.49是最大的；comp.2的比重是0.86，也还可以；而最后comp.3则没有多大分量，可以忽略。在这个数据结果下，我们可以取前两个合成变量做为数据的主成分。Loading载荷表示每个合成变量如何由原始变量来生成，如comp.1=0.418*Sepal.width+petal.length*-0.648+petal.width*-0.636。
| Comp.1 Comp.2 Comp.3
1.4904641 0.8624409 0.1863127 Loadings: Comp.1 Comp.2 Comp.3 Sepal.Width 0.418 0.907 Petal.Length -0.648 0.256 -0.717 Petal.Width -0.636 0.335 0.695 |
7. T校验
最后给大家介绍下统计中最简单最基本的T校验（即通过统计实验来判断一个结论是否是真的），让大家轻松下吧~
| #读入数据，厂商声明每袋糖的平均重量是500克，实际观测的50袋重量倾向于小于500克
x=scan(“sugar.txt”) #假设整体数据分布为正态分布，且中心为500，得到这样观测的概率p #由于实际观测的数据倾向于小于500克，所以用“less”的单边校验 result=t.test(x, m=500,alt=”less”) #对两个样本总体均值差异的校验，可用于验证产品试验效果 #result=t.test(data1，data2) result$p.value #如果p小于1%，即说明这种情况是小概率事件，有理由认为厂商的声明是错误的 if (result$p.value<0.01) { “the assumption is wrong!” } |
最近我面试了好多有意向来百度的同学，我深切的感觉到：学校中学习知识的思路太刻板了，同学们都绞尽脑汁的去回忆各种算法或概念的公式。这恰恰是忽略了事物的本质，学习一个知识或算法，就要站在这个知识或算法的发明者的角度去思考，究竟是为了解决怎样的问题？而这个方案中运用怎样的原理解决问题的？或许这就是工作后的学习更加务实了吧，所以期望这篇日志能对日常用统计来处理数据的同学有一定的帮助。第二部分的程序都是直接可运行来查看结果的，有需要相关实验数据的可以单独联系我。
ys2887 