R 数据分析个人总结 - 图文

Note of R

一、 Environment Setting ........................................................................................................................................ 2 二、 Data Processing ................................................................................................................................................ 2 三、 Decision Tree and Random Forest ................................................................................................................... 8 四、 Classification and Cluster ................................................................................................................................. 9 五、 Regression ...................................................................................................................................................... 11 六、 Principle Component Analysis ....................................................................................................................... 27 七、 SVM ................................................................................................................................................................ 38 八、 Neural Network ............................................................................................................................................. 39 九、 Web Mining ................................................................................................................................................. 39 十、 Time Serial Analysis ....................................................................................................................................... 39 十一、 Visualization .............................................................................................................................................. 40 十二、 Recommended Packages ........................................................................................................................... 44

1

一、Environment Setting

1. Installment

install.packages(\

“Error in loadNamespace”

install.packages(\install.packages(\install.packages(\警告: 无法将临时安装

‘C:\\Program Files\\R\\R-3.0.3\\library\\file30ec3b8b15ef\\Snowball 关闭百度防火墙和软件管家

Check all installed packages: packages(all.available=T) Installed.packages() getwd() Setwd() Ls(): list of

Rm(): delete variables

二、Data Processing

1. excel表格数据的读取

1. 1.1对于小文件的数据可以直接copy

打开excel表格->全选里面的内容->点击复制

data <- read.table(\

相信大家都是header=T或者F的区别，不知道可以自己尝试一下，越是自己试出来的记得越清楚，seq就是设定一个分隔符

使用这种方法的时候一定要注意复制！剪切板里面没有内容是无 1.2 read.csv

首先将excel文件另存为.csv格式然后直接使用read.csv(\

data<-read.csv(\write.csv(data,file=”directory.csv”, ); 默认header=T，可以自行改成F；

1.3 将xls另存为txt格式然后用read.table data<-read.table(\

1.4 安装RODBC包，使用里面的odbcConnectExcel()函数

Liabrary(RODBC);

Connection<-odbcConnect(dsn=”servername”, uid=”userid”,pwd=”****”); Query<-”SELECT ******** from”

myData<-sqlQuery(connection, query, errors=TRUE) odbcClose(connection)

2

a<-odbcConnectExcel(\ b<-sqlFetch(a,\ close(a)

1.5 xlsReadWrite包中的read.xls函数也可以读取excel文件

1.6 Save and load variables

Save(a,b, file=”/datafilename.Rdata”)# a, b -variables Rm(a,b): # delete variables

Load(“/datafilename.Rdata”) # load all variables into workplace

saveRDS(a, file=”/datafilename.rds”) # save variables into files; a2<- readRDS(“/datafilename.Rdata”) # load data into variables;

1.6 其他

R中导入数据的基本命令是scan()比较常用的命令有read.table()、read.csv()、read.fwf()等，这些命令其实在内部也是通过调用scan()实现的，它们可以看作是scan()的简化。读取数据到数据框的命令为read.table()。关于此文本的格式，有以下要求：行的分界符为换行/回车，列的分界符为空格或制表符(TAB)。特别注意，R把任意多个连续的空格或制表符只算作一个。由于R把连续空格或制表符只算作一个，如果数据表中有空缺的单元格，将会导致R读入时列发生错位。如果数据表中有空缺的单元格，一定要用NA来表示缺失的数据，否则将会导致读入数据错位。一般来说第一行用来存放变量名。如果第一行不是变量名，R用“V1, V2, ?”来命名变量。为了以后方便地跟踪数据，应当设置ID列。一般来说第一列用来存放ID。如果变量名比数据列少一个，R将假定第一列为ID列。

mydata <- read.table(\读入的数据进行小的修改

edit(x)或者fix(x)可以将数据打开成电子表格进行里面的一些修改 二进制文件的读取：对于二进制数据可以使用readBin()和writeBin()函数来读写或者write.foreign()来导出数据。 1.7 Attach

最近恰好看到“Symbols and Environments”这一章。发现LS诸位大大的分析其实均可以用“环境（environment）”来解释：R中的环境有很多级，最接近用户级的是global environment（全局环境？布吉岛这么翻译对不对）。而每个函数自己内部会创建一个local environment，同是函数有它的creating environment（创建环境）和calling environment（调用环境）。好，那么说一下attach干了什么，attach负责将data frame以及list中的所有元素添加到current environment（当前环境），而一般的current environment就是global environment。二楼“木友才”大大举的例子就很好的说明了这一点，数据data自己也有一个local environment，其中有x有y。但是如果不用attach的话，做线性回归的时候，lm函数的local environment中没有x和y，那怎么办？只能用通常从list中提取行or列的方法（按照名称提取），即data$x与data$y。但是attach可以将data中的object都加载进“当前环境”，因此调用lm函数的时候，data中的x和y都已经加载到了当前环境，因此可以直接写lm(y~x)。而问题也由此而生，假如有两个data，data1中有x、y，data2中也是x、y，那么把它们都加载到当前环境中就容易乱套。在调用时就不知道x和y指的是哪个data中的x和y。这就是为什么3L的“马甲1号”大大建议不要用attach，而是用with。用with的话就不会有这个问题。

1.8 SQL 设置

3

install.packages(\install.packages(\install.packages(\library(\library(\library(\

final <-sqldf(\ from a inner join b on a.id = b.id \

Error: could not find function \

SumAmount_by_year<- sqldf(\year, sum(GiftAmount) KeyIndicator='Individual' group by year\备注：select语句变量名不能含有”.”

from FundModified where

1.9.数据操作 1.7.1.变量变换

as.array(x),as.data.frame(x),as.numeric(x),as.logical(x),as.complex(x),as.character(x),...转换变量类型；使用如下命令可得到全部列表，methods(as) factor()：将一个向量转化为一个因子 1.7.2.变量信息

is.na(x),is.null(x),is.array(x),is.data.frame(x),is.numeric(x),is.complex(x),is.character (x),...检验变量的类型；使用如下命令得到全部列表，methods(is) length(x):x中元素的个数

dim(x):查看变量的维数；重新设置的维数，例如dim(x)=c(3,2) dimnames(x):重新设置对象的名称 nrow(x):行的个数 ncol(x):列的个数

class(x):得到或设置x的类；class(x)<-c(3,2) unclass(x):删除x的类

attr(x,which):得到或设置x的属性which attributes(obj):得到或设置obj的属性列表 1.7.2.fix,edit:对数据框数据进行表格形式的编辑 1.7.3.数据选取和操作

which.max(x):返回x中最大元素的指标 which.min(x):返回x中最小元素的指标 rev(x):翻转x中所有的元素

sort(x):升序排列x中的元素；降序排列使用：rev(sort(x))

cut(x,breaks):将x分割成为几段（或因子）；breaks为段数或分割点向量 match(x,y):返回一个和x长度相同且和y中元素相等的向量不等则返回NA which(x==a):如果比较操作为真（TRUE）,返回向量x的指针 choose(n,k):组合数的计算

na.omit(x):去除缺失值(NA)（去除相关行如果x为矩阵或数据框） na.fail(x):返回错误信息，如果x包含至少一个NA unique(x):如果x为向量或数据框，返回唯一值

table(x):返回一个由x不同值个数组成的表格(通常用于整数或因子), 即频数表

4

subset(x,...):根据条件（...选取x中元素，如x$V1<10）;如果x为数据框,选项select使用负号给出保留 或去除的变量 subset(x, subset, select,drop = FALSE, ...)

sample(x,size):不放回的随即在向量x中抽取size个元素,选项replace=TRUE允许放回抽取

prop.table(x,margin=):根据margin使用分数表示表格，wumargin时，所有元素和为1 *** seq() 产生一个向量序列

seq(from = 1, to = 1, by = ((to - from)/(length.out - 1)),length.out = NULL, along.with = NULL, ...) 其中length.out可简写为len。 *** rep() 重复一个对象

rep(x,times):x是要重复的对象（例如向量c(1,2,3)）,times为对象中每个元素重复的次数（如times=c(9,7,3)就是将x向量的1重复9次，2重复7次，3重复3次）。

除了上述主要的用法之外，还有一种特殊的用法：rep(x,times)重复x times次；使用each=来重复x元素each次；rep(c(1,2,3),2)得到1 2 3 1 2 3；rep(c(1,2,3),each=2)得到1 1 2 2 3 3。 *** 利用向量的索引实现定义分段函数

利用改变部分元素值的技术与逻辑值下标结合，可以定义向量的分段函数，例如，要定义

4.1.数组（包含矩阵、向量）创建

c(...) 常见的将一系列参数转化为向量的函数，通过recursive=TRUE降序排列列表并组合所有的元素为向量

from:to产生一个序列

\有较高的优先级；1:4+1得到\seq() 产生一个向量序列

seq(from = 1, to = 1, by = ((to - from)/(length.out - 1)),length.out = NULL, along.with = NULL, ...) 其中length.out可简写为len。 rep() 重复一个对象

rep(x,times):x是要重复的对象（例如向量c(1,2,3)）,times为对象中每个元素重复的次数（如times=c(9,7,3)就是将x向量的1重复9次，2重复7次，3重复3次）。

除了上述主要的用法之外，还有一种特殊的用法：rep(x,times)重复x times次；使用each=来重复x元素each次；rep(c(1,2,3),2)得到1 2 3 1 2 3；rep(c(1,2,3),each=2)得到1 1 2 2 3 3。 matrix(...)

创立矩阵 matrix(data = NA, nrow = 1, ncol = 1, byrow = FALSE,dimnames = NULL) array(x,dim=...)

产生由x组成的数组；使用类似dim=c(3,4,2)指定维数；如果x长度不够，则x自动循环 data.frame(...)

创建数据框，变量可能被命名或未被命名

data.frame(v=1:4,ch=c(\相对短的向量会被填充到最大向量长度。 注意数据框中各个属性的元素个数必须相同 list(...)

创建一个由变量组成的列表，变量可能被命名； list(a=c(1,2),b=\

factor(x,levels=) 把向量x编码称为因子

factor(x = character(), levels, labels = levels,exclude = NA, ordered = is.ordered(x)) gl() 产生因子变量

gl(n,k,lenth=n*k,labels=)

通过指定水平方式产生水平因子；k为水平的个数，n为重复的次数

5

rbind() cbind()

rbind(...)以行的形式组合矩阵，数据框，或其它 cbind(...)以列的方式组合，其他同rbind() 数据切割和分离 向量指标

x[n]:第n个元素

x[-n]:除了第n个元素的x x[1:n]:前n个元素

x[-(1:n)]:第n+1至最后的元素 x[c(1,4,2)]:指定元素

x[\名为\的元素 x[x>3]:所有大于3的元素 x[x>3 & x<5]:区间(3,5)的元素

x[x%in%c(\给定组中的元素 示例：

> x=c('and','a','b','c') > x

[1] \ \ \ > x[x%in%c(\[1] \

列表指标

x[n]:列表显示元素n

x［n］:列表的第n个元素 X[\名为\的元素 x$name :同上 矩阵指标

x[i,j]:下表为(i,j)的元素 x[i,]:第i行 x[,j]:第j列

x[,c(1,3)]:第1和3列

x[\名为\的行 x$name :同上 数据框指标

x［\］:列明为\的列 x$name :同上 4.2.R数学函数

4.2.1.泛用函数apply：

描述：让不同的函数作用数组array不同的边。 用法：apply(X, MARGIN, FUN, ...) 参数描述：

X为一数组array，当然包括矩阵。

MARGIN：为一向量，它指定函数要作用的下标。例如对一个矩阵来说，1代表行，2代表列，而c(1,2)代表行和列。

6

FUN作用于数组array上的函数。此函数可以是自己编写的。 注意：

当X时一个data.frame时，而不是array时，函数apply就会强制性的通过as.matrix或as.array将其转化为maxtrix或array。 例子：

## Compute row and column sums for a matrix: x <- cbind(x1 = 3, x2 = c(4:1, 2:5)) col.sums <- apply(x, 2, sum) row.sums <- apply(x, 1, sum)

4.2.2.基本符号：

+ :加 - :减 * :乘 / :除 ^ :幂 4.2.3.三角函数与其他基本数学函数

sin, cos, tan, asin, atan, atan2, log, log10, exp,max(x):

max(x):当x为一向量时，返回x中最大的元素。

当x为矩阵时，返回矩阵中所有元素的最大值。当要返回矩阵x每行或每列的最大值时，可以使用apply函数（apply(X, MARGIN, FUN, ...)），其中MARGIN表示要作用的方式，1代表行，2代表列，c(1,2)代表行和列；例如apply(data, MARGIN=1, max ) 当x为data.frame时，1仍然代表行，2仍然代表列。

当x为数组array时，1仍然代表行，2仍然代表列，其它更高维数用相应的编号数字来表示，如3就可以代表array的页。 其它详细情况参考apply函数 min(x):返回x中最小的元素，同max rev(x) 翻转x中的所有元素

sort(x) 升序排列x中的元素；降序排列使用rev(sort(x)),要得到排列的小标用order() choose(n,k) 组合数计算

which() 返回满足条件的下表

sample(x,size) 不放回的随机在向量x中抽取size个元素选项，replace+TRUE允许放回抽样 rank():计算向量的秩

range(x):返回c(min(x),max(x)) sum(x):x中各元素的加和 diff(x):向量x的差分

prod(x):x中元素的连乘积 mean(x):x的均值 median(x):x的中位数

quantile(x,probs=):满足给定概率的样本分位数（默认为0,0.25,0.75,1） weight.mean(x,w):加权平均，w即为weight,即权值。 rank(x):x中元素的秩

var(x):向量x的样本方差；如果x是矩阵或数据框，协方差矩阵将被计算 cor(x):如果x是矩阵或数据框，相关系数矩阵将被计算 sd(x):x的标准差；sd(x)=sqrt(var(x))

var(x,y) or cov(x,y):x和y间的协方差；如果x,y为矩阵或数据框，返回x和y各列的协方差 cor(x,y):x和y的线性相关系数；或者相关矩阵，如果x和y为矩阵或者数据框 round(x,n):x的约数，精确到n位

7

log(x,base):计算x以base为基的对数，默认基为exp(1)

scale(x):如果x是一个矩阵，则中心化和标准化数据；若只标准化数据，则使用选项center=FALSE，若只 中心化使用scale=FALSE(默认center=TRUE,scale=TRUE) pmin(x,y,...):x[i],y[i]相比较小者，组成新的向量 pmax(x,y,...):x[i],y[i]相比较大者，组成新的向量 Re(x):复数的实部 Im(x):复数的虚部 abs(x):绝对值

Arg(x):复数角度(in radians) Conj(x):共轭复数

fft(x):数组x的快速傅里叶变换

mvfft(x):矩阵x的每一列的傅里叶变换 factorial()：计算n!

cumsum():cumulative sums cumprod():cumulative products D(expression(exp(x^2)),\：求导 integrate(function(x) x^2,0,1)：积分

注意：大多数数学函数使用逻辑参数na.rm=FALSE来指定是否移除缺失值(NA)

4.3.R字符操作

paste(...):转化为字符后连接向量；seq=为分割界限（一个空格为默认）；选择collapse=可以分 割\结果

substr(x,start,stop):提取字符向量的子字段；同样可以赋值，使用substr(x,start,stop)<-value strsplit(x,split):在split的位置分割x,例如：

x<-strsplit(\

length(x［1］); #注意：对象被存放在列表x的第一项中，可以用x［1］提取 grep(pattern,x):搜索x中满足pattern条件；参见?regex

gsub(pattern,replacement,x):替换满足正则表达式的字段,sub()类似，但只替换第一个出现的字段 tolower(x):转化为lowercase toupper(x):转化为uppercase

match(x,table):table中匹配x元素位置组成的向量

x%in%table:table中匹配x元素位置组成的向量,返回值为逻辑值 pmatch(x,table):table中部分匹配x元素 nchar(x):字符的个数

三、Decision Tree and Random Forest

Examples:

dim(iris) names(iris) str(iris)

table(iris$Species) pie(table(iris$Species))

8

iris_ctree<-ctree(Species~Sepal.Length+Sepal.Width+Petal.Length+Petal.Width,data=iris) print(iris_ctree) plot(iris_ctree)

四、Classification and Cluster

聚类分析（Cluster Analysis） 是根据“物以类聚”的道理，对样品或指标进行分类的一种多元统计分析方法，它是在没有先验知识的情况下，对样本按各自的特性来进行合理的分类。

例子：

>X<-read.delim（\，header=T） >row.names（X）

<-c（\济南\，\青岛\，\淄博\，\枣庄\，\东营\，\烟台\，\潍坊\，\济宁\，\泰安\，\威海\，\日照\，\莱芜\，\临沂\，\德州\，\聊城\，\滨州\，\菏泽\） >d<-dist（scale（X）） >hc1<-hclust（d，\）#最短距离法 >hc2<-hclust（d，\）#最长距离法 >hc3<-hclust（d，\）#中间距离法 >hc4<-hclust（d，\）#Ward法 >opar<-par（mfrow=c（2，2））

>plot（hc1，hang=-1）；plot（hc2，hang=-1） >plot（hc3，hang=-1）；plot（hc4，hang=-1）

聚类分析被应用于很多方面，在商业上，聚类分析被用来发现不同的客户群，并且通过购买模式刻画不同的客户群的特征；在生物上，聚类分析被用来动植物分类和对基因进行分类，获取对种群固有结构的认识；在因特网应用上，聚类分析被用来在网上进行文档归类来修复信息。

聚类分析有两种主要计算方法，分别是凝聚层次聚类（Agglomerative hierarchical method）和K均值聚类（K-Means）。

1. 层次聚类

层次聚类又称为系统聚类，首先要定义样本之间的距离关系，距离较近的归为一类，较远的则属于不同的类。可用于定义“距离”的统计量包括了欧氏距离(euclidean)、马氏距离(manhattan)、 两项距离(binary)、明氏距离(minkowski)。还包括相关系数和夹角余弦。

9

层次聚类首先将每个样本单独作为一类，然后将不同类之间距离最近的进行合并，合并后重新计算类间距离。这个过程一直持续到将所有样本归为一类为止。在计算类间距离时则有六种不同的方法，分别是最短距离法、最长距离法、类平均法、重心法、中间距离法、离差平方和法。

下面我们用 iris 数据集来进行聚类分析，在R语言中所用到的函数为 hclust 。首先提取iris数据中的4个数值变量，然后计算其欧氏距离矩阵。然后将矩阵绘制热图，从图中可以看到颜色越深表示样本间距离越近，大致上可以区分出三到四个区块，其样本之间比较接近。 data=iris[,-5]

dist.e=dist(data,method='euclidean')

heatmap(as.matrix(dist.e),labRow = F, labCol = F)

然后使用hclust函数建立聚类模型，结果存在model1变量中，其中ward参数是将类间距离计算方法设置为离差平方和法。使用plot(model1)可以绘制出聚类树图。如果我们希望将类别设为3类，可以使用cutree函数提取每个样本所属的类别。 model1=hclust(dist.e,method='ward') result=cutree(model1,k=3)

为了显示聚类的效果，我们可以结合多维标度和聚类的结果。先将数据用MDS进行降维，然后以不同的的形状表示原本的分类，用不同的颜色来表示聚类的结果。可以看到setose品种聚类很成功，但有一些virginica品种的花被错误和virginica品种聚类到一起。 mds=cmdscale(dist.e,k=2,eig=T) x = mds$points[,1] y = mds$points[,2] library(ggplot2)

p=ggplot(data.frame(x,y),aes(x,y)) p+geom_point(size=3,alpha=0.8,

aes(colour=factor(result), shape=iris$Species))

10

2. K均值聚类

K均值聚类又称为动态聚类，它的计算方法较为简单，也不需要输入距离矩阵。首先要指定聚类的分类个数N，随机取N个样本作为初始类的中心，计算各样本与类中心的距离并进行归类，所有样本划分完成后重新计算类中心，重复这个过程直到类中心不再变化。

在R中使用 kmeans 函数进行K均值聚类，centers参数用来设置分类个数，nstart参数用来设置取随机初始中心的次数，其默认值为1，但取较多的次数可以改善聚类效果。model2$cluster可以用来提取每个样本所属的类别。

model2=kmeans(data,centers=3,nstart=10)

使用K均值聚类时需要注意，只有在类的平均值被定义的情况下才能使用，还要求事先给出分类个数。一种方法是先用层次聚类以决定个数，再用K均值聚类加以改进。或者以 轮廓系数 来判断分类个数。改善聚类的方法还包括对原始数据进行变换，如对数据进行降维后再实施聚类。

cluster扩展包中也有许多函数可用于聚类分析，如agnes函数可用于凝聚层次聚类，diana可用于划分层次聚类，pam可用于K均值聚类，fanny用于模糊聚类

五、Regression

11

R语言中的线性回归函数比较简单，就是lm()，比较复杂的是对线性模型的诊断和调整。这里结合Statistical Learning和杜克大学的Data Analysis and Statistical Inference的章节以及《R语言实战》的OLS(Ordinary Least Square)回归模型章节来总结一下，诊断多元线性回归模型的操作分析步骤。

1、选择预测变量

因变量比较容易确定，多元回归模型中难在自变量的选择。自变量选择主要可分为向前选择（逐次加使RSS最小的自变量），向后选择（逐次扔掉p值最大的变量）。个人倾向于向后选择法，一来p值比较直观，模型返回结果直接给出了各变量的p值，却没有直接给出RSS；二来当自变量比较多时，一个个加比较麻烦。

Call:

lm(formula=Sales~.+Income:Advertising+Age:Price,data=Carseats)

Residuals:

Min1QMedian3QMax

-2.9208-0.75030.01770.67543.3413

Coefficients:

EstimateStd.Errort valuePr(>|t|)

(Intercept)6.57556541.00874706.5192.22e-10*** CompPrice0.09293710.004118322.567<2e-16*** Income0.01089400.00260444.1833.57e-05*** Advertising0.07024620.02260913.1070.002030** Population0.00015920.00036790.4330.665330 Price-0.10080640.0074399-13.549<2e-16***

ShelveLocGood4.84867620.152837831.724<2e-16*** ShelveLocMedium1.95326200.125768215.531<2e-16*** Age-0.05794660.0159506-3.6330.000318*** Education-0.02085250.0196131-1.0630.288361 UrbanYes0.14015970.11240191.2470.213171 USYes-0.15755710.1489234-1.0580.290729

Income:Advertising0.00075100.00027842.6980.007290** Price:Age0.00010680.00013330.8010.423812 --- Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1

Residualstandard error:1.011on386degrees of freedom MultipleR-squared:0.8761,AdjustedR-squared:0.8719 F-statistic:210on13and386DF,p-value:<2.2e-16

构建一个回归模型后，先看F统计量的p值，这是对整个模型的假设检验，原假设是各系数都为0，如果连这个p值都不显著，无法证明至少有一个自变量对因变量有显著性影响，这个模型便不成立。然后看Adjusted R2，每调整一次模型，应该力使它变大；Adjusted R2越大说明模型中相关的自变量对因变量可解释的变异比例越大，模型的预测性就越好。

构建了线性模型后，如果是一元线性回归，可以画模型图初步判断一下线性关系（多元回归模型不好可视化）：

par(mfrow=c(1,1))

12

plot(medv~lstat,Boston)

fit1=lm(medv~lstat,data=Boston) abline(fit1,col=\

2、模型诊断

确定了回归模型的自变量并初步得到一个线性回归模型，并不是直接可以拿来用的，还要进行验证和诊断。诊断之前，先回顾多元线性回归模型的假设前提（by Data Analysis and Statistical Inference）：

（数值型）自变量要与因变量有线性关系； 残差基本呈正态分布；

残差方差基本不变（同方差性）； 残差（样本）间相关独立。

一个好的多元线性回归模型应当尽量满足这4点假设前提。

用lm()构造一个线性模型fit后，plot(fit)即可返回4张图（可以par(mfrow=c(2,2))一次展现），这4张

图可作初步检验：

左上图用来检验假设1，如果散点看不出什么规律，则表示线性关系良好，若有明显关系，则说明非线性关系明显。右上图用来检验假设2，若散点大致都集中在QQ图中的直线上，则说明残差正态性良好。左下图用来检验假设3，若点在曲线周围随机分布，则可认为假设3成立；若散点呈明显规律，比如方差随均值而增大，则越往右的散点上下间距会越大，方差差异就越明显。假设4的独立性无法通过这几张图来检验，只能通过数据本身的来源的意义去判断。 右下图是用来检验异常值。异常值与三个概念有关： 离群点：y远离散点主体区域的点

杠杆点：x远离散点主体区域的点，一般不影响回归直线的斜率

13

强影响点：影响回归直线的斜率，一般是高杠杆点。

对于多元线性回归，高杠杆点不一定就是极端点，有可能是各个变量的取值都正常，但仍然偏离散点主体。

对于异常值，可以谨慎地删除，看新的模型是否效果更好。 《R语言实战》里推荐了更好的诊断方法，总结如下。 （1）多元线性回归假设验证：

gvlma包的gvlma()函数可对拟合模型的假设作综合验证，并对峰度、偏度进行验证。 states<-as.data.frame(state.x77[,c(\\

fit<-lm(Murder~Population+Illiteracy+Income+ Frost,data=states) library(gvlma)

gvmodel<-gvlma(fit) summary(gvmodel) Call:

lm(formula=Murder~Population+Illiteracy+Income+Frost,data=states)

Residuals:

Min1QMedian3QMax

-4.7960-1.6495-0.08111.48157.6210

Coefficients:

EstimateStd.Errort valuePr(>|t|)

(Intercept)1.235e+003.866e+000.3190.7510 Population2.237e-049.052e-052.4710.0173* Illiteracy4.143e+008.744e-014.7382.19e-05*** Income6.442e-056.837e-040.0940.9253 Frost5.813e-041.005e-020.0580.9541 ---

Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1

Residualstandard error:2.535on45degrees of freedom MultipleR-squared:0.567,AdjustedR-squared:0.5285 F-statistic:14.73on4and45DF,p-value:9.133e-08

ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS USING THE GLOBAL TEST ON4DEGREES-OF-FREEDOM: LevelofSignificance=0.05 Call:

gvlma(x=fit)

14

Valuep-valueDecision

GlobalStat2.77280.5965Assumptionsacceptable. Skewness1.53740.2150Assumptionsacceptable. Kurtosis0.63760.4246Assumptionsacceptable. LinkFunction0.11540.7341Assumptionsacceptable.

Heteroscedasticity0.48240.4873Assumptionsacceptable.

最后的Global Stat是对4个假设条件进行综合验证，通过了即表示4个假设验证都通过了。最后的Heterosceasticity是进行异方差检测。注意这里假设检验的原假设都是假设成立，所以当p>0.05时，假设才能能过验证。

如果综合验证不通过，也有其他方法对4个假设条件分别验证： 线性假设 library(car) crPlots(fit)

返回的图是各个自变量与残差（因变量）的线性关系图，若存着明显的非线性关系，则需要对自变量作非线性转化。书中说这张图表明线性关系良好。

正态性 library(car)

qqPlot(fit,labels=row.names(states),id.method=\

qqPlot()可以生成交互式的qq图，选中异常点，就返回该点的名称。该图中除了Nevad点，其他点都在直线附近，可见正态性良好。

同方差性 library(car) ncvTest(fit)

15

Non-constantVarianceScoreTest Varianceformula:~fitted.values

Chisquare=1.746514Df=1p=0.1863156

p值大于0.05，可认为满足方差相同的假设。 独立性 library(car)

durbinWatsonTest(fit)

lagAutocorrelationD-WStatisticp-value 1-0.20069292.3176910.258 Alternativehypothesis:rho!=0

p值大于0.05，可认为误差之间相互独立。

除了以上4点基本假设，还有其他方面需要进行诊断。 多重共线性

理想中的线性模型各个自变量应该是线性无关的，若自变量间存在共线性，则会降低回归系数的准确性。一般用方差膨胀因子VIF(Variance Inflation Factor)来衡量共线性，《统计学习》中认为VIF超过5或10就存在共线性，《R语言实战》中认为VIF大于4则存在共线性。理想中的线性模型VIF=1，表完全不存在共线性。

library(car) vif(fit)

PopulationIlliteracyIncomeFrost

1.2452822.1658481.3458222.082547

可见这4个自变量VIF都比较小，可认为不存在多重共线性的问题。 异常值检验 离群点

离群点有三种判断方法：一是用qqPlot()画QQ图，落在置信区间（上图中两条虚线）外的即可认为是离群点，如上图中的Nevad点；一种是判断学生标准化残差值，绝对值大于2（《R语言实战》中认为2，《统计学习》中认为3）的可认为是离群点。

plot(x=fitted(fit),y=rstudent(fit)) abline(h=3,col=\

abline(h=-3,col=\

which(abs(rstudent(fit))>3) Nevada 28

还有一种方法是利用car包里的outlierTest()函数进行假设检验： library(car)

16

outlierTest(fit)

rstudent unadjusted p-valueBonferonnip Nevada3.5429290.000950880.047544

这个函数用来检验最大的标准化残差值，如果p>0.05，可以认为没有离群点；若p<0.05，则该点是离群点，但不能说明只有一个离群点，可以把这个点删除之后再作检验。第三种方法可以与第二种方法结合起来使用。

高杠杆点

高杠杆值观测点，即是与其他预测变量有关的离群点。换句话说，它们是由许多异常的预测变量值组合起来的，与响应变量值没有关系。《统计学习》中给出了一个杠杆统计量，《R语言实战》中给出了一种具体的操作方法。（两本书也稍有出入，《统计学习》中平均杠杆值为(p+1)/n，而在《R语言实战》中平均杠杆值为p/n；事实上在样本量n比较大时，几乎没有差别。）

hat.plot<-function(fit){ p<-length(coefficients(fit)) n<-length(fitted(fit))

plot(hatvalues(fit),main=\abline(h=c(2,3)*p/n,col=\

identify(1:n,hatvalues(fit),names(hatvalues(fit)))

#这句产生交互效果，选中某个点后，关闭后返回点的名称 }

hat.plot(fit)

超过2倍或3倍的平均杠杆值即可认为是高杠杆点，这里把Alaska和California作为高杠杆点。 强影响点

强影响点是那种若删除则模型的系数会产生明显的变化的点。一种方法是计算Cook距离，一般来说， Cook’s D值大于4/(n-k -1)，则表明它是强影响点，其中n 为样本量大小， k 是预测变量数目。

cutoff<-4/(nrow(states-length(fit$coefficients)-2)) #coefficients加上了截距项，因此要多减1 plot(fit,which=4,cook.levels=cutoff) abline(h=cutoff,lty=2,col=\

17

实际上这就是前面诊断的4张图之一，语句还是plot(fit)，which=4表示指定第4张图，cook.levels可设定标准值。红色虚线以上就返回了强影响点。

car包里的influencePlot()函数能一次性同时检查离群点、高杠杆点、强影响点。 library(car)

influencePlot(fit,id.method=\distance\

StudResHatCookD

Alaska1.7536920.432473190.4480510 Nevada3.5429290.095089770.2099157

纵坐标超过+2或小于-2的点可被认为是离群点，水平轴超过0.2或0.3的州有高杠杆值（通常为预测值的组合）。圆圈大小与影响成比例，圆圈很大的点可能是对模型参数的估计造成的不成比例影响的强影响点。

3、模型调整

到目前为止，《统计学习》中提到的多元线性回归模型潜在的问题，包括4个假设不成立、异常值、共线性的诊断方法在上面已经全部得到解决。这里总结、延伸《R语言实战》里提到的调整方法——

删除观测点

对于异常值一个最简单粗暴又高效的方法就是直接删除，不过有两点要注意。一是当数据量大的时候可以这么做，若数据量较小则应慎重；二是根据数据的意义判断，若明显就是错误就可以直接删除，否则需判断是否会隐藏着深层的现象。

另外删除观测点后要与删除之前的模型作比较，看模型是否变得更好。

18

变量变换

在进行非线性变换之前，先看看4个假设是否成立，如果成立可以不用变换；没必要追求更好的拟合效果而把模型搞得太复杂，这有可能出现过拟合现象。如果连假设检验都不通过，可以通过变量变换来调整模型。这里只讨论线性关系不佳的情况，其他情况遇到了再说。

（1）多项式回归

如果残差图中呈现明显的非线性关系，可以考虑对自变量进行多项式回归。举一个例子： library(MASS) library(ISLR)

fit1=lm(medv~lstat,data=Boston) plot(fit1,which=1)

可以看到这个一元线性回归模型的残差图中，散点的规律还是比较明显，说明线性关系较弱。

fit1_1<-lm(medv~poly(lstat,2),data=Boston) plot(fit1_1,which=1)

将自变量进行2次多项式回归后，发现现在的残差图好多了，散点基本无规律，线性关系较明显。再看看两个模型的整体效果——

summary(fit1)

19

summary(fit1_1)

可见多项式回归的模型Adjusted R2也增大了，模型的解释性也变强了。 多项式回归在《统计学习》后面的非线性模型中还会提到，到时候再讨论。

（2）Box-Tidwell变换

car包中的boxTidwell() 函数通过获得预测变量幂数的最大似然估计来改善线性关系。 library(car)

boxTidwell(Murder~Population+Illiteracy,data=states)

ScoreStatisticp-value MLE oflambda

Population-0.32280030.74684650.8693882 Illiteracy0.61938140.53566511.3581188

iterations=19

#这里看lambda，表示各个变量的幂次数

lmfit<-lm(Murder~Population+Illiteracy,data=states)

lmfit2<-lm(Murder~I(Population^0.87)+I(Illiteracy^1.36),data=states) plot(lmfit,which=1)

plot(lmfit2,which=1)

summary(lmfit)

20

summary(lmfit2)

可以发现残差图和Adjusted R2的提升都甚微，因此没有必要作非线性转换。 4、模型分析 （1）模型比较

前面只是简单得用Adjusted R2来比较模型，《R语言实战》里介绍了可以用方差分析来比较嵌套模型（即它的一些项完全包含在另一个模型中）有没有显著性差异。方差分析的思想是：如果线性模型y~x1+x2+x3与y~x1+x2没有显著性差异，若同时x3变量对模型也不显著，那就没必要加上变量x3。下面进行试验：

aovfit1<-lm(Murder~Population+Illiteracy+Income+Frost,data=states) aovfit2<-lm(Murder~Population+Illiteracy,data=states) anova(aovfit1,aovfit2)

AnalysisofVarianceTable

Model1:Murder~Population+Illiteracy+Income+Frost Model2:Murder~Population+Illiteracy Res.DfRSSDfSumofSqFPr(>F) 145289.17

247289.25-2-0.0785050.00610.9939

summary(aovfit1)

Coefficients:

EstimateStd.Errort valuePr(>|t|)

(Intercept)1.235e+003.866e+000.3190.7510 Population2.237e-049.052e-052.4710.0173* Illiteracy4.143e+008.744e-014.7382.19e-05*** Income6.442e-056.837e-040.0940.9253 Frost5.813e-041.005e-020.0580.9541

Residualstandard error:2.535on45degrees of freedom MultipleR-squared:0.567,AdjustedR-squared:0.5285 F-statistic:14.73on4and45DF,p-value:9.133e-08

summary(aovfit2)

Coefficients:

EstimateStd.Errort valuePr(>|t|)

21

(Intercept)1.652e+008.101e-012.0390.04713* Population2.242e-047.984e-052.8080.00724** Illiteracy4.081e+005.848e-016.9788.83e-09***

Residualstandard error:2.481on47degrees of freedom MultipleR-squared:0.5668,AdjustedR-squared:0.5484 F-statistic:30.75on2and47DF,p-value:2.893e-09

Income和Frost两个变量不显著，两个模型之间没有显著性差异，就可以不加这两个变量。删去这两个不显著的变量后，R2略微减少，Adjusted R2增大，这也符合二者的定义。

《R语言实战》里还介绍到了用AIC(Akaike Information Criterion，赤池信息准则)值来比较模型，AIC值越小的模型优先选择，原理不明。

aovfit1<-lm(Murder~Population+Illiteracy+Income+Frost,data=states) aovfit2<-lm(Murder~Population+Illiteracy,data=states) AIC(aovfit1,aovfit2)

df AIC

aovfit16241.6429 aovfit24237.6565

第二个模型AIC值更小，因此选第二个模型（真是简单粗暴）。注：ANOVA需限定嵌套模型，AIC则不需要。可见AIC是更简单也更实用的模型比较方法。

（2）变量选择

这里的变量选择与最开始的变量选择同也不同，虽然是一回事，但一开始是一个粗略的变量的选择，主要是为了构建模型；这里则要进行细致的变量选择来调整模型。

逐步回归

前面提到的向前或向后选择或者是同时向前向后选择变量都是逐步回归法。MASS包中的stepAIC() 函数可以实现逐步回归模型（向前、向后和向前向后），依据的是精确AIC准则。以下实例是向后回归法：

library(MASS)

aovfit1<-lm(Murder~Population+Illiteracy+Income+Frost,data=states) stepAIC(aovfit1,direction=\

# “forward”为向前选择，\为向后选择，\为混合选择

Start:AIC=97.75

Murder~Population+Illiteracy+Income+Frost

DfSumofSqRSS AIC

-Frost10.021289.1995.753 -Income10.057289.2295.759 289.1797.749

-Population139.238328.41102.111 -Illiteracy1144.264433.43115.986

Step:AIC=95.75

Murder~Population+Illiteracy+Income

22

DfSumofSqRSS AIC

-Income10.057289.2593.763 289.1995.753

-Population143.658332.85100.783 -Illiteracy1236.196525.38123.605

Step:AIC=93.76

Murder~Population+Illiteracy

DfSumofSqRSS AIC 289.2593.763

-Population148.517337.7699.516 -Illiteracy1299.646588.89127.311 Call:

lm(formula=Murder~Population+Illiteracy,data=states)

Coefficients:

(Intercept)PopulationIlliteracy 1.65154970.00022424.0807366

可见原本的4元回归模型向后退了两次，最终稳定成了2元回归模型，与前面模型比较的结果一致。 全子集回归

《R语言实战》里提到了逐步回归法的局限：不是每个模型都评价了，不能保证选择的是“最佳”模型。比如上例中，从Murder ~ Population + Illiteracy + Income + Frost到Murder ~ Population + Illiteracy + Income再到Murder~Population+Illiteracy虽然AIC值确实在减少，但Murder ~ Population + Illiteracy + Frost没评价，如果遇到变量很多的情况下，逐步回归只沿一个方向回归，就有可能错过最优的回归方向。

library(leaps)

leaps<-regsubsets(Murder~Population+Illiteracy+Income+Frost,data=states,nbest=4) plot(leaps,scale=\

横坐标是变量，纵坐标是Adjusted R2，可见除截距项以外，只选定Population和Illiteracy这两个变量，可以使线性模型有最大的Adjusted R2。

全子集回归比逐步回归范围更广，模型优化效果更好，但是一旦变量数多了之后，全子集回归迭代的次数就很多，就会很慢。

23

事实上，变量的选择不是机械式地只看那几个统计指标，更主要的是根据数据的实际意义，从业务角度上来选择合适的变量。

线性模型变量的选择在《统计学习》后面的第6章还会继续讲到，到时继续综合讨论。

（3）交互项

交互项《统计学习》中花了一定篇幅来描写，但在《R语言实战》是在方差分析章节中讨论。添加变量间的交互项有时可以改善线性关系，提高Adjusted R2。针对数据的实际意义，如果两个基本上是独立的，也很难产生交互、产生协同效应的变量，那就不必考虑交互项；只有从业务角度分析，有可能产生协同效应的变量间才考虑交互项。

涉及到交互项有一个原则：如果交互项是显著的，那么即使变量不显著，也要放在回归模型中；若变量和交互项都不显著，则可以都不放。

（4）交叉验证

Andrew Ng的Machine Learning中就提到了，模型对旧数据拟合得好不一定就对新数据预测得好。因此一个数据集应当被分两训练集和测试集两部分（或者训练集、交叉验证集、测试集三部分），训练好的模型还要在新数据中测试性能。

所谓交叉验证，即将一定比例的数据挑选出来作为训练样本，另外的样本作保留样本，先在训练样本上获取回归方程，然后在保留样本上做预测。由于保留样本不涉及模型参数的选择，该样本可获得比新数据更为精确的估计。

在k 重交叉验证中，样本被分为k个子样本，轮流将k-1个子样本组合作为训练集，另外1个子样本作为保留集。这样会获得k 个预测方程，记录k 个保留样本的预测表现结果，然后求其平均值。

bootstrap包中的crossval()函数可以实现k重交叉验证。 shrinkage<-function(fit,k=10){ require(bootstrap) # define functions

theta.fit<-function(x,y){ lsfit(x,y) }

theta.predict<-function(fit,x){ cbind(1,x)%*%fit$coef }

# matrix of predictors

x<-fit$model[,2:ncol(fit$model)] # vector of predicted values y<-fit$model[,1]

results<-crossval(x,y,theta.fit,theta.predict,ngroup=k) r2<-cor(y,fit$fitted.values)^2 r2cv<-cor(y,results$cv.fit)^2

cat(\

cat(k,\cat(\}

这个自定义的shrinkage()函数用来做k重交叉验证，比计算训练集和交叉验证集的R方差异。这个函数里涉及到一个概念：复相关系数。复相关系数实际上就是y和fitted(y)的简单相关系数。对于一元线性回归，

24

R2就是简单相关系数的平方；对于多元线性回归，R2是复相关系数的平方。这个我没有成功地从公式上推导证明成立，就记下吧。这个方法用到了自助法的思想，这个在统计学习后面会细致讲到。

fit<-lm(Murder~Population+Income+Illiteracy+ Frost,data=states) shrinkage(fit)

OriginalR-square=0.5669502

10FoldCross-ValidatedR-square=0.441954 Change=0.1249963

可见这个4元回归模型在交叉验证集中的R2下降了0.12之多。若换成前面分析的2元回归模型—— fit2<-lm(Murder~Population+Illiteracy,data=states) shrinkage(fit2)

OriginalR-square=0.5668327

10FoldCross-ValidatedR-square=0.517304 Change=0.04952868

这次R2下降只有约0.05。R2减少得越少，则预测得越准确。 5、模型应用 （1）预测

最重要的应用毫无疑问就是用建立的模型进行预测了。构建好模型后，可用predict()函数进行预测—— fit2<-lm(Murder~Population+Illiteracy,data=states) predict(fit2,

newdata=data.frame(Population=c(2000,3000),Illiteracy=c(1.7,2.2)), interval=\

fit lwr upr

19.0371748.00491110.06944 211.3017299.86685112.73661

这里newdata提供了两个全新的点供模型来预测。还可以用interval指定返回置信区间（confidence）或者预测区间（prediction），这也反映了统计与机器学习的一个差异——可解释性。注意置信区间考虑的是平均值，而预测区间考虑的是单个观测值，所以预测区间永远比置信区间广，因此预测区间考虑了单个观测值的不可约误差；而平均值同时也把不可约误差给抵消掉了。

（2）相对重要性

有的时候需要解释模型中各个自变量对因变量的重要程度，简单处理可以直接看系数即可，《R语言实战》里自定义了一个relweights()函数可以计算各个变量的权重：

relweights<-function(fit,...){ R<-cor(fit$model) nvar<-ncol(R)

rxx<-R[2:nvar,2:nvar] rxy<-R[2:nvar,1] svd<-eigen(rxx) evec<-svd$vectors ev<-svd$values

delta<-diag(sqrt(ev))

25

# correlations between original predictors and new orthogonal variables lambda<-evec%*Tlta%*%t(evec) lambdasq<-lambda^2

# regression coefficients of Y on orthogonal variables beta<-solve(lambda)%*%rxy rsquare<-colSums(beta^2) rawwgt<-lambdasq%*?ta^2 import<-(rawwgt/rsquare)*100 lbls<-names(fit$model[2:nvar]) rownames(import)<-lbls

colnames(import)<-\

# plot results

barplot(t(import),names.arg=lbls,ylab=\

xlab=\sub=paste(\...)

return(import) }

不要在意算法原理和代码逻辑这种细节，直接看结果： fit<-lm(Murder~Population+Illiteracy+Income+ Frost,data=states)

relweights(fit,col=\

Weights

Population14.723401 Illiteracy59.000195 Income5.488962 Frost20.787442

fit$coefficients

(Intercept)PopulationIlliteracyIncomeFrost

1.23456341120.00022367544.14283659030.00006442470.0005813055

26

barplot(fit$coefficients[2:5])

在本例中，相对权重与系数的排序结果一致。推荐用相对权重。

六、Principle Component Analysis

setwd(\

x1=c(149.3, 161.2, 171.5, 175.5, 180.8, 190.7,202.1, 212.4, 226.1, 231.9, 239.0) x2=c(4.2, 4.1, 3.1, 3.1, 1.1, 2.2, 2.1, 5.6, 5.0, 5.1, 0.7)

x3=c(108.1, 114.8, 123.2, 126.9, 132.1, 137.7,146.0, 154.1, 162.3, 164.3, 167.6) y=c(15.9, 16.4, 19.0, 19.1, 18.8, 20.4, 22.7,26.5, 28.1, 27.6, 26.3) p<-princomp(~x1+x2+x3,cor=T) summary(p,loadings=TRUE) loadings(p) pre<-predict(p) z1<-pre[,1] z2<-pre[,2]

r2<-lm(y~z1+z2) summary(r2)

主成分分析（PCA）是一种数据降维技巧，它能将大量相关变量转化为一组很少的不相关变量，这些无关变量称为主成分。

探索性因子分析（EFA）是一系列用来发现一组变量的潜在结构的方法，通过寻找一组更小的、潜在的或隐藏的结构来解释已观测到的、变量间的关系。

1.R中的主成分和因子分析R的基础安装包中提供了PCA和EFA的函数，分别为princomp （）和factanal（）psych包中有用的因子分析函数

函数 描述 principal（） 含多种可选的方差放置方法的主成分分析 fa（） 可用主轴、最小残差、加权最小平方或最大似然法估计的因子分析 fa.parallel（） 含平等分析的碎石图 factor.plot（） 绘制因子分析或主成分分析的结果 27

fa.diagram（） 绘制因子分析或主成分分析的载荷矩阵 scree（） 因子分析和主成分分析的碎石图

PCA/EFA 分析流程：（1）数据预处理；PCA和EFA都是根据观测变量间的相关性来推导结果。用户可以输入原始数据矩阵或相关系数矩阵列到principal()和fa（）函数中，若输出初始结果，相关系数矩阵将会被自动计算，在计算前请确保数据中没有缺失值；

（2）选择因子分析模型。判断是PCA（数据降维）还是EFA（发现潜在结构）更符合你的分析目标。若选择EFA方法时，还需要选择一种估计因子模型的方法（如最大似然估计）。 （3）判断要选择的主成分/因子数目； （4）选择主成分/因子； （5）旋转主成分/因子； （6）解释结果；

（7）计算主成分或因子得分。

2.主成分分析PCA的目标是用一组较少的不相关变量代替大量相关变量，同时尽可能保留初始变量的信息，这些推导所得的变量称为主成分，它们是观测变量的线性组合。如第一主成分为：

PC1=a1X1+a2X2+??+akXk 它是k个观测变量的加权组合，对初始变量集的方差解释性最大。

第二主成分是初始变量的线性组合，对方差的解释性排第二， 同时与第一主成分正交（不相关）。后面每一个主成分都最大化它对方差的解释程度，同时与之前所有的主成分都正交，但从实用的角度来看，都希望能用较少的主成分来近似全变量集。

（1）判断主成分的个数PCA中需要多少个主成分的准则： 根据先验经验和理论知识判断主成分数；

根据要解释变量方差的积累值的阈值来判断需要的主成分数； 通过检查变量间k*k的相关系数矩阵来判断保留的主成分数。

最常见的是基于特征值的方法，每个主成分都与相关系数矩阵的特征值 关联，第一主成分与最大的特征值相关联，第二主成分与第二大的特征值相关联，依此类推。

Kaiser-Harris准则建议保留特征值大于1的主成分，特征值小于1的成分所解释的方差比包含在单个变量中的方差更少。

Cattell碎石检验则绘制了特征值与主成分数的图形，这类图形可以展示图形弯曲状况，在图形变化最大处之上的主成分都保留。

最后，还可以进行模拟，依据与初始矩阵相同大小的随机数矩阵来判断要提取的特征值。若基于真实数据的某个特征值大于一组随机数据矩阵相应的平均特征值，那么该主成分可以保留。该方法称作平行分析。

利用fa.parallel（）函数，可同时对三种特征值判别准则进行评价。 [plain] view plaincopy library(psych)

fa.parallel(USJudgeRatings[,-1],fa=\plot with parallel analysis\

28

碎石头、特征值大于1准则和100次模拟的平行分析（虚线）都表明保留一个主成分即可保留数据集的大部分信息，下一步是使用principal（）函数挑选出相应的主成分。

（2）提取主成分principal（）函数可根据原始数据矩阵或相关系数矩阵做主成分分析 格式为：principal（的，nfactors=,rotate=,scores=） 其中：r是相关系数矩阵或原始数据矩阵； nfactors设定主成分数（默认为1）；

rotate指定旋转的方式[默认最大方差旋转（varimax）] scores设定是否需要计算主成分得分（默认不需要）。 [plain] view plaincopy 美国法官评分的主成分分析 library(psych)

pc<-principal(USJudgeRatings[,-1],nfactors=1) pc

此处，输入的是没有ONT变量的原始，并指定获取一个未旋转的主成分。由于PCA只对相关系数矩阵进行分析，在获取主成分前，原始数据将会被自动转换为相关系数矩阵。PC1栏包含了成分载荷，指观测变量与主成分的相关系数。如果提取不止一个主成分，则还将会有PC2、PC3等栏。成分载荷（component loadings）可用来解释主成分的含义。此处可看到，第一主成分（PC1）与每个变量都高度相关，也就是说，它是一个可用来进行一般性评价的维度。

29

h2柆指成分公因子方差-----主成分对每个变量的方差解释度。 u2栏指成分唯一性-------方差无法 被主成分解释的比例（1-h2）。

SS loadings行包含了主成分相关联的特征值，指的是与特定主成分相关联的标准化后的方差值。 Proportin Var行表示的是每个主成分对整个数据集的解释程度。

结果不止一个主成分的情况 [plain] view plaincopy library(psych)

fa.parallel(Harman23.cor$cov,n.obs=302,fa=\parallel analysis\

plot with

载荷阵解释了成分和因子的含义，第一成分与每个身体测量指标都正相关，看起来似乎是一个一般性的衡量因子；第二主成分与前四个变量负相关，与后四个变量正相关，因此它看起来似乎是一个长度容量因子。但理念上的东西都不容易构建，当提取了多个成分时，对它们进行旋转可使结果更具有解释性。

（3）主成分旋转旋转是一系列将成分载荷阵变得更容易解释的数学方法，它们尽可能地对成分去噪。旋转方法有两种：使选择的成分保持不相关（正效旋转），和让它们变得相关（斜交旋转）。旋转方法也会依据去噪定义的不同而不同。最流行的下次旋转是方差极大旋转，它试图对载荷阵的列进行去噪，使得每个成分只是由一组有限的变量来解释（即载荷阵每列只有少数几个很大的载荷，其他都是很小的载荷）。 [plain] view plaincopy

install.packages(\library(GPArotation)

rc<-principal(Harman23.cor$cov,nfactors=2,rotate=\rc

30

列

名从PC变成了RC，以表示成分被旋转

观察可以发现第一主成分主要由前四个变量来解释，第二主成分主要由变量5到变量8来解释。

注意两个主成分仍不相关，对变量的解释性不变，这是因为变量的群组没有发生变化。另外，两个主成分放置后的累积方差解释性没有变化，变的只是各个主成分对方差的解释（成分1从58%变为44%，成分2从22%变为37%）。各成分的方差解释度趋同，准确来说，此时应该称它们为成分而不是主成分。

（4）获取主成分得分利用principal（）函数，很容易获得每个调查对象在该主成分上的得分。 [plain] view plaincopy 从原始数据中获取成分得分

[plain] view plaincopy library(psych)

pc<-principal(USJudgeRatings[,-1],nfactors=1,score=TRUE) head(pc$scores)

当scores=TRUE时，主成

分得分存储在principal（）函数返回对象的scores元素中。

31

[plain] view plaincopy

cor(USJudgeRatings$CONT,PC$scores)

[plain] view plaincopy

获取主成分得分的系数

[plain] view plaincopy library(psych)

rc<-principal(Harman23.cor$cov,nfactor=2,rotate=\round(unclass(rc$weights),2)

得到主成分得分：

PC1=0.28*height+0.30*arm.span+0.30*forearm+0.29*lower.leg-0.06*weight-0.08*bitro.diameter-0.10*chest.girth-0.04*chest.width

PC2=-0.05*height-0.08*arm.span-0.09*forearm-0.06*lower.leg+0.33*weight+0.32*bitro.diameter+0.34*chest.girth+0.27*chest.width

3.探索性因子分析EFA的目标是通过发掘隐藏在数据下的一组较少的、更为基本的无法观测的变量，来解释一组可观测变量的相关性。这些虚拟的、无法观测的变量称作因子。（每个因子被认为可解释多个观测变量间共有的方差，也叫作公共因子） 模型的形式为：

Xi=a1F1+a2F2+??apFp+Ui

Xi是第i个可观测变量（i=1,2,??k） Fj是公共因子（j=1,2,??p） 并且p

[plain] view plaincopy options(digits=2)

covariances<-ability.cov$cov

correlations<-cov2cor(covariances) correlations

32

ability.cov提

供了变量的协方差矩阵

cov2cor（）函数将其转化为相关系数矩阵

（1）判断需提取的公共因子数[plain] view plaincopy library(psych)

convariances<-ability.cov$cov

correlations<-cov2cor(covariances)

fa.parallel(correlations,n.obs=112,fa=\

若使用PCA方法，可能会选择一个成分或两个成分。当摇摆不定时，高估因子数通常比低估因子数的结果好，因为高估因子数一般较少曲解“真实”情况。 （2）提取公共因子可使用fa（）函数来提取因子 fa（）函数的格式为：

fa（r，nfactors=，n.obs=，rotate=，scores=，fm） r是相关系数矩阵或原始数据矩阵；

nfactors设定提取的因子数（默认为1）；

n.obs是观测数（输入相关系数矩阵时需要填写）； rotate设定放置的方法（默认互变异数最小法）； scores设定是否计算因子得分（默认不计算）； fm设定因子化方法（默认极小残差法）。

与PCA不同，提取公共因子的方法很多，包括最大似然法（ml）、主轴迭代法（pa）、加权最小二乘法（wls）、广义加权最小二乘法（gls）和最小残差法（minres）。

[plain] view plaincopy

33

未旋转的主轴迭代因子法

[plain] view plaincopy

fa<-fa(correlations,nfactors=2,rotate=\fa

（3）因子旋转[plain] view plaincopy 用正交旋转提取因子

[plain] view plaincopy

fa.varimax<-fa(correlations,nfactors=2,rotate=\fa.varimax

[plain] view plaincopy

正交放置将人为地强制两个因子不相关

[plain] view plaincopy

正交旋转，因子分析的重点在于因子结构矩阵（变量与因子的相关系数）

34

用斜交旋转提取因子 [plain] view plaincopy

fa.promax<-fa(correlations,nfactors=2,rotate=\fa.promax

[plain] view plaincopy

对于斜交旋转，因子分析会考虑三个矩阵：因子结构矩阵、因子模式矩阵和因子关联矩阵

[plain] view plaincopy

因子模式矩阵即标准化的回归系数矩阵，它列出了因子的预测变量的权重；

[plain] view plaincopy

因子关联矩阵即因子相关系数矩阵；

[plain] view plaincopy

因子结构矩阵（或称因子载荷阵），可使用公式F=P*Phi来计算得到，其中F是载荷阵，P为因子模式矩阵，Phi为因子关联矩阵。

35

[plain] view plaincopy fsm<-function(oblique){

if(class(oblique)[2]==\ warning(\ }else{

P<-unclass(oblique$loading) F<-P%*%oblique$Phi

colnames(F)<-c(\ return (F) } }

fsm(fa.promax)

36

可以看到变量与因子间的相关系数。将它们与正交旋转所得因子载荷阵相比，发现该载荷阵列的噪音较大，这是因为之前允许潜在因子相关。虽然斜交方法更为复杂，但模型将更加符合真实数据。

使用factor.plot（）或fa.diagram（）函数，可绘制正交或斜交结果的图形 [plain] view plaincopy

factor.plot(fa.promax,labels=rownames(fa.promax$loadings))

[plain] view plaincopy

fa.diagram(fa.promax,simple=TRUE)

（4）因子得分

37

EFA并不十分关注因子得分，在fa（）函数中添加score=TRUE选项，便可轻松地得到因子得分。另外还可以得到得分系数（标准化的回归权重），它在返回对象的weights元素中。

[plain] view plaincopy fa.promax$weights

4.其他

（1）对因子分析非常有用的软件包，FactoMineR包不仅提供了PCA和EFA方法，还包含潜变量模型。FAiR包使用遗传算法来估计因子分析模型，增强了模型参数估计能力，能够处理不等式的约束条件；GPArotation包提供了许多因子旋转方法. nFactors包，提供了用来判断因子数目方法。

（2）其他潜变量模型

先验知识的模型：先从一些先验知识开始，比如变量背后有几个因子、变量在因子上的载荷是怎样的、因子间的相关性如何，然后通过收集数据检验这些先验知识。这种方法称作验证性因子分析（CFA）。做CFA的软件包：sem、openMx和lavaan等。

ltm包可以用来拟合测验和问卷中各项目的潜变量模型。

潜类别模型（潜在的因子被认为是类别型而非连续型）可通过FlexMix、lcmm、randomLCA和poLC包进行拟合。lcda包可做潜类别判别分析，而lsa可做潜在语义分析----一种自然语言处理中的方法。ca包提供了可做简单和多重对应分析的函数。

R中还包含了众多的多维标度法（MDS）计算工具。MDS即可用发现解释相似性和可测对象间距离的潜在维度。

cmdscale（）函数可做经典的MDS

MASS包中的isoMDS（）函数可做非线性MDS vagan包中则包含了两种MDS的函数

七、SVM

Install.packages(“e1071”,dependencies=T) data(iris) attach(iris)

subdata<-iris[iris$Species!='virginica',] subdata$Species<-factor(subdata$Species)

model<-svm(Species~Petal.Length+Petal.Width,data=subdata) summary(model)

38

八、Neural Network 九、Web Mining

十、Time Serial Analysis

library('forecast')

Amount_by_month<- sqldf(\Year, Month, sum(GiftAmount) from FundModified group by Year,Month order by Year,Month\qplot(,Amount_by_month[,3])

acf(Amount_by_month[,3])//自相关参数

pt<-ts(Amount_by_month[,3],start=2008,frequency=12) //配置数据 arima1<-auto.arima(pt,trace=T)//自动计算参数

airfore<-forecast(arima1,h=30,fan=T)//预测 plot(airfore) //画图

（1）根据趋势定差分

plot(lostjob,type=\查看图像总体趋势，确定如何差分 df1 = diff(lostjob) d=1阶差分

s4_df1=diff(df1,4) 对d=1阶差分结果进行k=4步（季节）差分

（2）根据所定差分检验平稳

adfTest(s4_df1,lag=6) 对差分结果进行平稳性检验

（3）ARIMA(p,d,q)中的pq定阶 acf(s4_df1) pacf(s4_df1)

（4）建立arima模型

ans=arima(lostjob,order=c(4,1,0),seasonal=list(order=c(1,0,1),period=4),include.mean=F,fixed=c(NA,0,0,NA,NA,NA))

（5）检验模型残差白噪声

//use natural log of T (the number ofobservations) which provides higher power (1 -Beta) Box.test(s4_df1,lag=5,type='Ljung')

Box.test(ans$residuals,lag=5,type='Ljung') 或者

tsdiag(ans)

（6）预测

39

predict(ans,10)

十一、Visualization

1)Plotting 1 - simples graph library(ggplot2)

data(mpg) # belongs to ggplot2 qplot(displ, hwy, data=mpg)

2)Plotting 2 - grouping by colour library(ggplot2)

qplot(displ,hwy, data = mpg, color=drv)

3)Plotting 3 - adding a geoms library(ggplot2)

qplot(displ,hwy, data = mpg, geom=c(\

40

4)Plotting 4 - adding a geoms

qplot(displ,hwy, data = mpg, geom=\

5)Plotting 5 - adding a geoms

qplot(displ, hwy, data = mpg, geom=\6)Plotting 6 - adding a geoms

qplot(hwy, data = mpg, geom=\

7)Plotting 7 - adding a geoms

qplot(hwy, data = mpg, geom=\

41

8)Histograms

library(ggplot2)qplot(hwy, data = mpg, fill=drv) 9)Withouth rows definition (1 row)

library(ggplot2)qplot(displ, hwy, data=mpg, facets = . ~ drv) 10)Withouth columns definition (1 column)

library(ggplot2)qplot(hwy, data=mpg, facets = drv ~ .) 11)Facets

#install.packages(\ library(UsingR) data(Galton) galton = Galton

##install.packages(\ library(ggplot2) library(reshape) long = melt(galton) ## Using as id variables

g = ggplot(long, aes(x=value, fill= variable))

g = g + geom_histogram(colour = \ g = g + facet_grid(. ~variable) g

12)Scatterplot (without color grouping) library(ggplot2)

qplot(displ, hwy, data = mpg, shape=drv) 13)Scatterplot (without color grouping) library(ggplot2)

qplot(displ, hwy, data = mpg, shape=drv, color=drv) 14)Scatterplot of multiple factors library(ggplot2)

qplot(displ, hwy, data = mpg, geom=c(\

15)Scatter plot of multiple factors with linear regression line in each one (same graph) library(ggplot2)

qplot(displ, hwy, data = mpg, geom=c(\ 16)Multiple scatterplot lines of best fit at same window grouped by one factor library(ggplot2)

qplot(displ, hwy, data = mpg, geom=c(\ 17)Multiple scatterplot linear regression lines at same window grouped by one factor library(ggplot2)

qplot(displ, hwy, data = mpg, geom=c(\ 18)Skeleton dataset tutorial geom_points()

skeleton = read.table(\library(ggplot2)

skeleton$Sex = factor(skeleton$Sex, levels=c(1,2))g = ggplot(skeleton, aes(Age, DGDifference))g = g + geom_point() # Adds a geom parameterg

42

Line of best fit

skeleton = read.table(\header=T)library(ggplot2)skeleton$Sex = factor(skeleton$Sex, levels=c(1,2))g = ggplot(skeleton, aes(Age, DGDifference))g = g + geom_point() # Adds a geom parameterg = g + geom_smooth() # Adds a line of best fitg Linear regression line

skeleton = read.table(\header=T)library(ggplot2)skeleton$Sex = factor(skeleton$Sex, levels=c(1,2))g = ggplot(skeleton, aes(Age, DGDifference))g = g + geom_point() # Adds a geom parameterg = g + geom_smooth(method=\Using facets (multple graphs)

skeleton = read.table(\header=T)library(ggplot2)skeleton$Sex = factor(skeleton$Sex, levels=c(1,2))g = ggplot(skeleton, aes(Age, DGDifference))g = g + geom_point() g = g + geom_smooth(method=\Using factets (colors)

library(ggplot2) skeleton = read.table(\header=T) qplot(Age, DGDifference,data=skeleton, facets=.~BMI,geom=c(\Facets by sex

library(ggplot2)skeleton = read.table(\header=T)skeleton$Sex = factor(skeleton$Sex, levels=c(1,2))qplot(Age, DGDifference,data=skeleton, facets=.~Sex,geom=c(\Controlling easthetics color, size, alpha

skeleton = read.table(\header=T)library(ggplot2)skeleton$Sex = factor(skeleton$Sex, levels=c(1,2))g = ggplot(skeleton, aes(Age, DGDifference))g = g + geom_point(color=\size=4, alpha=1/2) g

geom_point SIZE ALPHA AES

skeleton = read.table(\header=T)library(ggplot2)skeleton$Sex = factor(skeleton$Sex, levels=c(1,2))g = ggplot(skeleton, aes(Age, DGDifference))g = g + geom_point(size=4, alpha=1/2, aes(color=BMI)) g

Adding Labels, x axis and y axis text

skeleton = read.table(\header=T)library(ggplot2)skeleton$Sex = factor(skeleton$Sex, levels=c(1,2))g = ggplot(skeleton, aes(Age, DGDifference))g = g + geom_point(aes(color=BMI)) + labs(title=\Adding Labels, x axis and y axis text and changing the regression line

skeleton = read.table(\header=T)library(ggplot2)skeleton$Sex = factor(skeleton$Sex, levels=c(1,2))g = ggplot(skeleton, aes(Age, DGDifference))g = g + geom_point(size=4, alpha=1/2, aes(color=BMI))g = g + geom_smooth(color=\Applying a theme

skeleton = read.table(\header=T)library(ggplot2)skeleton$Sex = factor(skeleton$Sex, levels=c(1,2))g = ggplot(skeleton, aes(Age, DGDifference))g = g + geom_point(aes(color=BMI)) g = g + theme_bw() # simply changes the themeg Changing theme font

skeleton = read.table(\header=T)library(ggplot2)skeleton$Sex = factor(skeleton$Sex, levels=c(1,2))g = ggplot(skeleton, aes(Age, DGDifference))g = g + geom_point(aes(color=BMI)) g = g + theme_bw(base_family=\Cartesian Coordinates

testdat = data.frame(x=1:100, y = rnorm(100))testdat[50,2] = 100 ## Intencional outliersg = ggplot(testdat,

43

aes(x=x, y=y))g = g + geom_line() g = g + coord_cartesian(ylim= c(-3,3))g 2．Clustering

Dataframe, distance, hclust, plot dataFrame <- data.frame(x = x, y = y) distxy<-dist(dataFrame)

hClustering <- hclust(distxy) plot(hClustering)

十二、Recommended Packages

Many useful R function come in packages, free libraries of code written by R's active user community. To install an R package, open an R session and type at the command line

install.packages(\

R will download the package from CRAN, so you'll need to be connected to the internet. Once you have a package installed, you can make its contents available to use in your current R session by running

library(\

There are thousands of helpful R packages for you to use, but navigating them all can be a challenge. To help you out, we've compiled this guide to some of the best. We've used each of these, and found them to be outstanding – we've even written some of them. But you don't have to take our word for it, these packages are also some of the top most downloaded R packages.

To load data

RODBC, RMySQL, RPostgresSQL, RSQLite - If you'd like to read in data from a database, these packages are a good place to start. Choose the package that fits your type of database.

XLConnect, xlsx - These packages help you read and write Micorsoft Excel files from R. You can also just export your spreadsheets from Excel as .csv's.

foreign - Want to read a SAS data set into R? Or an SPSS data set? Foreign provides functions that help you load data files from other programs into R.

R can handle plain text files – no package required. Just use the functions read.csv, read.table, and read.fwf. If you have even more exotic data, consult the CRAN guide to data import and export.

44

To manipulate data

dplyr - Essential shortcuts for subsetting, summarizing, rearranging, and joining together data sets. dplyr is our go to package for fast data manipulation.

tidyr - Tools for changing the layout of your data sets. Use the gather and spread functions to convert your data into the tidy format, the layout R likes best.

stringr - Easy to learn tools for regular expressions and character strings.

lubridate - Tools that make working with dates and times easier.

To visualize data

ggplot2 - R's famous package for making beautiful graphics. ggplot2 lets you use the grammar of graphics to build layered, customizable plots.

ggvis - Interactive, web based graphics built with the grammar of graphics.

rgl - Interactive 3D visualizations with R

htmlwidgets - A fast way to build interactive (javascript based) visualizations with R. Packages that implement htmlwidgets include:

leaflet (maps)

dygraphs (time series) DT (tables)

diagrammeR (diagrams) network3D (network graphs)

threeJS (3D scatterplots and globes).

googleVis - Let's you use Google Chart tools to visualize data in R. Google Chart tools used to be called Gapminder, the graphing software Hans Rosling made famous in hie TED talk.

To model data

car - car's Anova function is popular for making type II and type III Anova tables.

mgcv - Generalized Additive Models

lme4/nlme - Linear and Non-linear mixed effects models

randomForest - Random forest methods from machine learning

45

multcomp - Tools for multiple comparison testing

vcd - Visualization tools and tests for categorical data

glmnet - Lasso and elastic-net regression methods with cross validation

survival - Tools for survival analysis

caret - Tools for training regression and classification models

To report results

shiny - Easily make interactive, web apps with R. A perfect way to explore data and share findings with non-programmers.

R Markdown - The perfect workflow for reproducible reporting. Write R code in your markdown reports. When you run render, R Markdown will replace the code with its results and then export your report as an HTML, pdf, or MS Word document, or a HTML or pdf slideshow. The result? Automated reporting. R Markdown is integrated straight into RStudio.

xtable - The xtable function takes an R object (like a data frame) and returns the latex or HTML code you need to paste a pretty version of the object into your documents. Copy and paste, or pair up with R Markdown.

For Spatial data

sp, maptools - Tools for loading and using spatial data including shapefiles.

maps - Easy to use map polygons for plots.

ggmap - Download street maps straight from Google maps and use them as a background in your ggplots.

For Time Series and Financial data

zoo - Provides the most popular format for saving time series objects in R.

xts - Very flexible tools for manipulating time series data sets.

quantmod - Tools for downloading financial data, plotting common charts, and doing technical analysis.

To write high performance R code

Rcpp - Write R functions that call C++ code for lightning fast speed.

data.table - An alternative way to organize data sets for very, very fast operations. Useful for big data.

46

parallel - Use parallel processing in R to speed up your code or to crunch large data sets.

To work with the web

XML - Read and create XML documents with R

jsonlite - Read and create JSON data tables with R

httr - A set of useful tools for working with http connections

To write your own R packages

devtools - An essential suite of tools for turning your code into an R package.

testthat - testthat provides an easy way to write unit tests for your code projects.

roxygen2 - A quick way to document your R packages. roxygen2 turns inline code comments into documentation pages and builds a package namespace.

You can also read about the entire package development process online in Hadley Wickham's R Packages book

47

本文来源：https://www.bwwdw.com/article/md28.html