《R语言实战》自学笔记71-主成分和因子分析

Python028

《R语言实战》自学笔记71-主成分和因子分析,第1张

主成分分析

成分分析((Principal Component Analysis,PCA)是一种数据降维技巧,它能将大量相关变量转化为一组很少的不相关变量,这些无关变量称为主成分(原来变量的线性组合)。整体思想就是化繁为简,抓住问题关键,也就是降维思想。

主成分分析法是通过恰当的数学变换,使新变量——主成分成为原变量的线性组合,并选取少数几个在变差总信息量中比例较大的主成分来分析事物的一种方法。主成分在变差信息量中的比例越大,它在综合评价中的作用就越大。

因子分析

探索性因子分析法(Exploratory Factor Analysis,EFA)是一系列用来发现一组变量的潜在结构的方法。它通过寻找一组更小的、潜在的或隐藏的结构来解释已观测到的、显式的变量间的关系。

PCA与EFA模型间的区别

参见图14-1。主成分(PC1和PC2)是观测变量(X1到X5)的线性组合。形成线性组合的权重都是通过最大化各主成分所解释的方差来获得,同时还要保证个主成分间不相关。相反,因子(F1和F2)被当做是观测变量的结构基础或“原因”,而不是它们的线性组合。

R的基础安装包提供了PCA和EFA的函数,分别为princomp()和factanal()。

最常见的分析步骤

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

(2)选择因子模型。判断是PCA(数据降维)还是EFA(发现潜在结构)更符合你的研究目标。如果选择EFA方法,你还需要选择一种估计因子模型的方法(如最大似然估计)。

(3)判断要选择的主成分/因子数目。

(4)选择主成分/因子。

(5)旋转主成分/因子。

(6)解释结果。

(7)计算主成分或因子得分。

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

它是k个观测变量的加权组合,对初始变量集的方差解释性最大。第二主成分也是初始变量的线性组合,对方差的解释性排第二,同时与第一主成分正交(不相关)。后面每一个主成分都最大化它对方差的解释程度,同时与之前所有的主成分都正交。理论上来说,你可以选取与变量数相同的主成分,但从实用的角度来看,我们都希望能用较少的主成分来近似全变量集。

主成分与原始变量之间的关系

(1)主成分保留了原始变量绝大多数信息。

(2)主成分的个数大大少于原始变量的数目。

(3)各个主成分之间互不相关。

(4)每个主成分都是原始变量的线性组合。

数据集USJudgeRatings包含了律师对美国高等法院法官的评分。数据框包含43个观测,12个变量。

用来判断PCA中需要多少个主成分的准则:

根据先验经验和理论知识判断主成分数;

根据要解释变量方差的积累值的阈值来判断需要的主成分数;

通过检查变量间k × k的相关系数矩阵来判断保留的主成分数。

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

Kaiser-Harris准则建议保留特征值大于1的主成分,特征值小于1的成分所解释的方差比包含在单个变量中的方差更少。Cattell碎石检验则绘制了特征值与主成分数的图形。这类图形可以清晰地展示图形弯曲状况,在图形变化最大处之上的主成分都可保留。最后,你还可以进行模拟,依据与初始矩阵相同大小的随机数据矩阵来判断要提取的特征值。若基于真实数据的某个特征值大于一组随机数据矩阵相应的平均特征值,那么该主成分可以保留。该方法称作平行分析。

图形解读:线段和x符号组成的图(蓝色线):特征值曲线;

红色虚线:根据100个随机数据矩阵推导出来的平均特征值曲线;

绿色实线:特征值准则线(即:y=1的水平线)

判别标准:特征值大于平均特征值,且大于y=1的特征值准则线,被认为是可保留的主成分。根据判别标准,保留1个主成分即可。

fa.parallel函数学习

fa.parallel(data,n.obs=,fa=”pc”/”both”,n.iter=100,show.legend=T/F)

data:原始数据数据框;

n.obs:当data是相关系数矩阵时,给出原始数据(非原始变量)个数,data是原始数据矩阵时忽略此参数;

fa:“pc”为仅计算主成分,“fa”为因子分析,“both”为计算主成分及因子;

n.iter:模拟平行分析次数;

show.legend:显示图例。

principal(r, nfactors = , rotate = , scores = )

r:相关系数矩阵或原始数据矩阵;

nfactors:设定主成分数(默认为1);

rotate:指定旋转的方法,默认最大方差旋转(varimax)。

scores:设定是否需要计算主成分得分(默认不需要)。

PC1栏包含了成分载荷,指观测变量与主成分的相关系数。如果提取不止一个主成分,那么还将会有PC2、PC3等栏。成分载荷(component loadings)可用来解释主成分的含义,解释主成分与各变量的相关程度。

h2栏为成分公因子方差,即主成分对每个变量的方差解释度。

u2栏为成分唯一性,即方差无法被主成分解释的部分(1-h2)。

SS loadings包含了与主成分相关联的特征值,其含义是与特定主成分相关联的标准化后的方差值,即可以通过它来看90%的方差可以被多少个成分解释,从而选出主成分(即可使用nfactors=原始变量个数来把所有特征值查出,当然也可以直接通过eigen函数对它的相关矩阵进行查特征值)。

Proportion Var表示每个主成分对整个数据集的解释程度。

Cumulative Var表示各主成分解释程度之和。

Proportion Explained及Cumulative Proportion分别为按现有总解释方差百分比划分主成分及其累积百分比。

结果解读:第一主成分(PC1)与每个变量都高度相关,也就是说,它是一个可用来进行一般性评价的维度。ORAL变量99.1%的方差都可以被PC1来解释,仅仅有0.91%的方差不能被PC1解释。第一主成分解释了11个变量92%的方差。

结果解读:通过碎石图可以判定选择的主成分个数为2个。

结果解读:从结果Proportion Var: 0.58和0.22可以判定,第一主成分解释了身体测量指标58%的方差,而第二主成分解释了22%,两者总共解释了81%的方差。对于高度变量,两者则共解释了其88%的方差。

旋转是一系列将成分载荷阵变得更容易解释的数学方法,它们尽可能地对成分去噪。旋转方法有两种:使选择的成分保持不相关(正交旋转),和让它们变得相关(斜交旋转)。旋转方法也会依据去噪定义的不同而不同。最流行的正交旋转是方差极大旋转,它试图对载荷阵的列进行去噪,使得每个成分只是由一组有限的变量来解释(即载荷阵每列只有少数几个很大的载荷,其他都是很小的载荷)。 结果列表中列的名字都从PC变成了RC,以表示成分被旋转。

当scores = TRUE时,主成分得分存储在principal()函数返回对象的scores元素中。

如果你的目标是寻求可解释观测变量的潜在隐含变量,可使用因子分析。

EFA的目标是通过发掘隐藏在数据下的一组较少的、更为基本的无法观测的变量,来解释一

组可观测变量的相关性。这些虚拟的、无法观测的变量称作因子。(每个因子被认为可解释多个

观测变量间共有的方差,因此准确来说,它们应该称作公共因子。)

其中 是第i个可观测变量(i = 1…k), 是公共因子(j = 1…p),并且p<k。 是 变量独有的部分(无法被公共因子解释)。 可认为是每个因子对复合而成的可观测变量的贡献值。

碎石检验的前两个特征值(三角形)都在拐角处之上,并且大于基于100次模拟数据矩阵的特征值均值。对于EFA,Kaiser-Harris准则的特征值数大于0,而不是1。

结果解读:PCA结果建议提取一个或者两个成分,EFA建议提取两个因子。

fa(r, nfactors=, n.obs=, rotate=, scores=, fm=)

 r是相关系数矩阵或者原始数据矩阵;

 nfactors设定提取的因子数(默认为1);

 n.obs是观测数(输入相关系数矩阵时需要填写);

 rotate设定旋转的方法(默认互变异数最小法);

 scores设定是否计算因子得分(默认不计算);

 fm设定因子化方法(默认极小残差法)。

与PCA不同,提取公共因子的方法很多,包括最大似然法(ml)、主轴迭代法(pa)、加权最小二乘法(wls)、广义加权最小二乘法(gls)和最小残差法(minres)。统计学家青睐使用最大似然法,因为它有良好的统计性质。

结果解读:两个因子的Proportion Var分别为0.46和0.14,两个因子解释了六个心理学测试60%的方差。

结果解读:阅读和词汇在第一因子上载荷较大,画图、积木图案和迷宫在第二因子上载荷较大,非语言的普通智力测量在两个因子上载荷较为平均,这表明存在一个语言智力因子和一个非语言智力因子。

正交旋转和斜交旋转的不同之处。

对于正交旋转,因子分析的重点在于因子结构矩阵(变量与因子的相关系数),而对于斜交旋转,因子分析会考虑三个矩阵:因子结构矩阵、因子模式矩阵和因子关联矩阵。

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

图形解读:词汇和阅读在第一个因子(PA1)上载荷较大,而积木图案、画图和迷宫在第二个因子(PA2)上载荷较大。普通智力测验在两个因子上较为平均。

与可精确计算的主成分得分不同,因子得分只是估计得到的。它的估计方法有多种,fa()函数使用的是回归方法。

R包含了其他许多对因子分析非常有用的软件包。FactoMineR包不仅提供了PCA和EFA方法,还包含潜变量模型。它有许多此处我们并没考虑的参数选项,比如数值型变量和类别型变量的使用方法。FAiR包使用遗传算法来估计因子分析模型,它增强了模型参数估计能力,能够处理不等式的约束条件,GPArotation包则提供了许多因子旋转方法。最后,还有nFactors包,它提供了用来判断因子数目的许多复杂方法。

主成分分析

1.数据导入

数据结构:对10株玉米进行了生物学性状考察,考察指标有株高,穗位,茎粗,穗长,秃顶,穗粗,穗行数,行粒数。

结果解读:选择2个主成分即可保留样本大量信息。

3.提取主成分

结果解读:主成分1可解释44%的方差,主成分2解释了26%的方差,合计解释了70%的方差。

4.获取主成分得分

5.主成分方程

PC1 = 0.27 株高 - 0.04 穗位 + 0.29 茎粗 - 0.01 穗长 - 0.21 秃顶 - 0.13 穗粗 + 0.16 穗行数 + 0.24 行粒数

PC2 = -0.01 株高 + 0.36 穗位 - 0.10 茎粗 + 0.41 穗长 - 0.08 秃顶 + 0.43 穗粗 - 0.15 穗行数 + 0.01 行粒数

图形解读:此图反映了变量与主成分的关系,三个蓝点对应的RC2值较高,点上的标号2,4,6对应变量名穗位,穗长,穗粗,说明第2主成分主要解释了这些变量,与这些变量相关性强;黑点分别对应株高,茎粗,穗行数,行粒数,说明第一主成分与这些变量相关性强,第一主成分主要解释的也是这些变量,而5号点秃顶对于两个主成分均没有显示好的相关性。

因子分析

图解:可以看到需要提取4个因子。

2.提取因子

结果解读:因子1到4解释了80%的方差。

3.获取因子得分

图解:可以看出,因子1和因子2的相关系数为0.4,行粒数,株高,茎粗,秃顶在因子1的载荷较大,穗长,穗位在因子2上的载荷较大;因子3只有穗行数相关,因子4只有穗粗相关。

参考资料:

用R语言进行关联分析

关联是两个或多个变量取值之间存在的一类重要的可被发现的某种规律性。关联分析目的是寻找给定数据记录集中数据项之间隐藏的关联关系,描述数据之间的密切度。

几个基本概念

1. 项集

这是一个集合的概念,在一篮子商品中的一件消费品即为一项(Item),则若干项的集合为项集,如{啤酒,尿布}构成一个二元项集。

2. 关联规则

一般记为的形式,X为先决条件,Y为相应的关联结果,用于表示数据内隐含的关联性。如:,表示购买了尿布的消费者往往也会购买啤酒。

关联性强度如何,由三个概念——支持度、置信度、提升度来控制和评价。

例:有10000个消费者购买了商品,其中购买尿布1000个,购买啤酒2000个,购买面包500个,同时购买尿布和面包800个,同时购买尿布和面包100个。

3. 支持度(Support)

支持度是指在所有项集中{X, Y}出现的可能性,即项集中同时含有X和Y的概率:

该指标作为建立强关联规则的第一个门槛,衡量了所考察关联规则在“量”上的多少。通过设定最小阈值(minsup),剔除“出镜率”较低的无意义规则,保留出现较为频繁的项集所隐含的规则。

设定最小阈值为5%,由于{尿布,啤酒}的支持度为800/10000=8%,满足基本输了要求,成为频繁项集,保留规则;而{尿布,面包}的支持度为100/10000=1%,被剔除。

4. 置信度(Confidence)

置信度表示在先决条件X发生的条件下,关联结果Y发生的概率:

这是生成强关联规则的第二个门槛,衡量了所考察的关联规则在“质”上的可靠性。相似的,我们需要对置信度设定最小阈值(mincon)来实现进一步筛选。

具体的,当设定置信度的最小阈值为70%时,置信度为800/1000=80%,而的置信度为800/2000=40%,被剔除。

5. 提升度(lift)

提升度表示在含有X的条件下同时含有Y的可能性与没有X这个条件下项集中含有Y的可能性之比:

该指标与置信度同样衡量规则的可靠性,可以看作是置信度的一种互补指标。

R中Apriori算法

算法步骤:

1. 选出满足支持度最小阈值的所有项集,即频繁项集;

2. 从频繁项集中找出满足最小置信度的所有规则。

>library(arules) #加载arules包

>click_detail =read.transactions("click_detail.txt",format="basket",sep=",",cols=c(1)) #读取txt文档(文档编码为ANSI)

>rules <- apriori(click_detail, parameter =list(supp=0.01,conf=0.5,target="rules")) #调用apriori算法

>rules

set of419 rules

>inspect(rules[1:10]) #查看前十条规则

解释

1)library(arules):加载程序包arules,当然如果你前面没有下载过这个包,就要先install.packages(arules)

2)click_detail =read.transactions("click_detail.txt",format="basket",sep=",",cols=c(1)):读入数据

read.transactions(file, format =c("basket", "single"), sep = NULL,

cols = NULL, rm.duplicates =FALSE, encoding = "unknown")

file:文件名,对应click_detail中的“click_detail.txt”

format:文件格式,可以有两种,分别为“basket”,“single”,click_detail.txt中用的是basket。

basket: basket就是篮子,一个顾客买的东西都放到同一个篮子,所有顾客的transactions就是一个个篮子的组合结果。如下形式,每条交易都是独立的。

文件形式:

item1,item2

item1

item2,item3

读入后:

items

1 {item1,

item2}

2 {item1}

3 {item2,

item3}

single: single的意思,顾名思义,就是单独的交易,简单说,交易记录为:顾客1买了产品1, 顾客1买了产品2,顾客2买了产品3……(产品1,产品2,产品3中可以是单个产品,也可以是多个产品),如下形式:

trans1 item1

trans2 item1

trans2 item2

读入后:

items transactionID

1 {item1}trans1

2 {item1,

item2}trans2

sep:文件中数据是怎么被分隔的,默认为空格,click_detail里面用逗号分隔

cols:对basket, col=1,表示第一列是数据的transaction ids(交易号),如果col=NULL,则表示数据里面没有交易号这一列;对single,col=c(1,2)表示第一列是transaction ids,第二列是item ids

rm.duplicates:是否移除重复项,默认为FALSE

encoding:写到这里研究了encoding是什么意思,发现前面txt可以不是”ANSI”类型,如果TXT是“UTF-8”,写encoding=”UTF-8”,就OK了.

3)rules <- apriori(click_detail,parameter = list(supp=0.01,conf=0.5,target="rules")):apriori函数

apriori(data, parameter = NULL, appearance = NULL, control = NULL)

data:数据

parameter:设置参数,默认情况下parameter=list(supp=0.1,conf=0.8,maxlen=10,minlen=1,target=”rules”)

supp:支持度(support)

conf:置信度(confidence)

maxlen,minlen:每个项集所含项数的最大最小值

target:“rules”或“frequent itemsets”(输出关联规则/频繁项集)

apperence:对先决条件X(lhs),关联结果Y(rhs)中具体包含哪些项进行限制,如:设置lhs=beer,将仅输出lhs含有beer这一项的关联规则。默认情况下,所有项都将无限制出现。

control:控制函数性能,如可以设定对项集进行升序sort=1或降序sort=-1排序,是否向使用者报告进程(verbose=F/T)

补充

通过支持度控制:rules.sorted_sup = sort(rules, by=”support”)

通过置信度控制:rules.sorted_con = sort(rules, by=”confidence”)

通过提升度控制:rules.sorted_lift = sort(rules, by=”lift”)

Apriori算法

两步法:

1. 频繁项集的产生:找出所有满足最小支持度阈值的项集,称为频繁项集;

2. 规则的产生:对于每一个频繁项集l,找出其中所有的非空子集;然后,对于每一个这样的子集a,如果support(l)与support(a)的比值大于最小可信度,则存在规则a==>(l-a)。

频繁项集产生所需要的计算开销远大于规则产生所需的计算开销

频繁项集的产生

几个概念:

1, 一个包含K个项的数据集,可能产生2^k个候选集

2,先验原理:如果一个项集是频繁的,则它的所有子集也是频繁的(理解了频繁项集的意义,这句话很容易理解的);相反,如果一个项集是非频繁的,则它所有子集也一定是非频繁的。

3基于支持度(SUPPORT)度量的一个关键性质:一个项集的支持度不会超过它的子集的支持度(很好理解,支持度是共同发生的概率,假设项集{A,B,C},{A,B}是它的一个自己,A,B,C同时发生的概率肯定不会超过A,B同时发生的概率)。

上面这条规则就是Apriori中使用到的,如下图,当寻找频繁项集时,从上往下扫描,当遇到一个项集是非频繁项集(该项集支持度小于Minsup),那么它下面的项集肯定就是非频繁项集,这一部分就剪枝掉了。

一个例子(百度到的一个PPT上的):

当我在理解频繁项集的意义时,在R上简单的复现了这个例子,这里采用了eclat算法,跟apriori应该差不多:

代码:

item <- list(

c("bread","milk"),

c("bread","diaper","beer","eggs"),

c("milk","diaper","beer","coke"),

c("bread","milk","diaper","beer"),

c("bread","milk","diaper","coke")

)

names(item) <- paste("tr",c(1:5),sep = "")

item

trans <- as(item,"transactions") #将List转为transactions型

rules = eclat(trans,parameter = list(supp = 0.6,

target ="frequent itemsets"),control = list(sort=1))

inspect(rules) #查看频繁项集

运行后结果:

>inspect(rules)

items support

1{beer,

diaper}0.6

2{diaper,

milk} 0.6

3{bread,

diaper}0.6

4{bread,

milk} 0.6

5{beer} 0.6

6{milk} 0.8

7{bread} 0.8

8{diaper} 0.8

以上就是该例子的所有频繁项集,然后我发现少了{bread,milk,diaper}这个项集,回到例子一看,这个项集实际上只出现了两次,所以是没有这个项集的。

规则的产生

每个频繁k项集能产生最多2k-2个关联规则

将项集Y划分成两个非空的子集X和Y-X,使得X ->Y-X满足置信度阈值

定理:如果规则X->Y-X不满足置信度阈值,则X’->Y-X’的规则一定也不满足置信度阈值,其中X’是X的子集

Apriori按下图进行逐层计算,当发现一个不满足置信度的项集后,该项集所有子集的规则都可以剪枝掉了。