专业编程基础技术教程

网站首页 > 基础教程 正文

用R语言做数据分析——探索性因子分析

ccvgpt 2024-10-12 15:00:24 基础教程 73 ℃

EFA的目标是通过发掘隐藏在数据下的一组较少的、更为基本的无法观测的变量,来解释一组可观测变量的相关性。这些虚拟的、无法观测的变量称作因子。(每个因子被认为可解释多个观测变量间共有的方差,因此准确说,它们应该称作公共因子。)

用R语言做数据分析——探索性因子分析

模型的形式为:

其中Xi是第i个可观测变量(i=1...k),Fj是公共因子(j=1...p),并且p<k。Ui是Xi变量独有的部分(无法被公共因子解释)。ai可认为是每个因子对复合而成的可观测变量的贡献值。我们以Harman74.cor数据集为例,它包含了24个心理测验间的相互关系,受试对象为145个七年级或八年级的学生。我们认为每个个体在24个心理学测验上的观测得分,是根据四个潜在心理学因素的加权能力值组合而成。

虽然PCA和EFA存在差异,但是它们的许多分析步骤都是相似的。为阐述EFA的分析过程,我们用它来对六个心理学测验间的相关性进行分析。112个人参与了六个测验,包括非语言的普通测验(general)、画图测验(picture)、积木图案测验(blocks)、迷津测验(maze)、阅读测验(reading)和词汇测验(vocab)。我们如何用一组较少的、潜在的心理学因素来解释参与者的测验得分呢?

数据集ability.cov提供了变量的协方差矩阵,我们可用cov2cor()函数将其转化为相关系数矩阵。数据集没有缺失值。

> options(digits = 2)

> covariances <- ability.cov$cov

> correlations <- cov2cor(covariances)

> correlations

general picture blocks maze reading vocab

general 1.00 0.47 0.55 0.34 0.58 0.51

picture 0.47 1.00 0.57 0.19 0.26 0.24

blocks 0.55 0.57 1.00 0.45 0.35 0.36

maze 0.34 0.19 0.45 1.00 0.18 0.22

reading 0.58 0.26 0.35 0.18 1.00 0.79

vocab 0.51 0.24 0.36 0.22 0.79 1.00

因为要寻求用来解释数据的潜在结构,可使用EFA方法。与使用PCA相同,下一步工作为判断需要提取几个因子

判断需提取的公共因子数

可用fa.parallel()函数可判断需提取的因子数:

> covariances <- ability.cov$cov

> correlations <- cov2cor(covariances)

> fa.parallel(correlations, n.obs = 112, fa="both", n.iter = 100,main = "Scree plots with parallel analysis")

判断

上图是判断心理学测验需要保留的因子数。图中展示了PCA和EFA的结果。PCA结果建议提取一个或者两个成分,EFA建议提取两个因子

图形中有几个值得注意的地方。如果使用PCA方法,我们可能会选择一个成分(碎石检验和平行检验)或者两个成分(特征值大于1)。当摇摆不定时,高估因子数通常比低估因子数的结果好,因为高估因子数一般较少曲解“真实”情况。

观察EFA的结果,显然需提取两个因子。碎石检验的前两个特征值(三角形)都在拐角处之上,并且大于基于100次模拟数据矩阵额特征值的均值。对于EFA,Kaiser-Harris准则的特征值数大于0,而不是1.图形中该准则也建议选择两个因子。

提取公共因子

现在我们决定提取两个因子,可以使用fa()函数获得相应的结果。fa()函数的格式如下:

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

其中:

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

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

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

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

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

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

与PCA不同,提取公共因子的方法很多,包括极大似然法(ml)、主迭代法(pa)、加权最小二乘法(wls)、广义加权最小二乘法(gls)和最小残差法(minres)。统计学家青睐使用最大似然法,因为它有良好的统计性质。不过有时候最大似然法不会收敛,此时使用主迭代法效果会很好。下面的代码,将使用主迭代法(fm="pa")提取未旋转的因子。

> fa <- fa(correlations,nfactors = 2,rotate = "none", fm="pa")

> fa

Factor Analysis using method = pa

Call: fa(r = correlations, nfactors = 2, rotate = "none", fm = "pa")

Standardized loadings (pattern matrix) based upon correlation matrix

PA1 PA2 h2 u2 com

general 0.75 0.07 0.57 0.432 1.0

picture 0.52 0.32 0.38 0.623 1.7

blocks 0.75 0.52 0.83 0.166 1.8

maze 0.39 0.22 0.20 0.798 1.6

reading 0.81 -0.51 0.91 0.089 1.7

vocab 0.73 -0.39 0.69 0.313 1.5

PA1 PA2

SS loadings 2.75 0.83

Proportion Var 0.46 0.14

Cumulative Var 0.46 0.60

Proportion Explained 0.77 0.23

Cumulative Proportion 0.77 1.00

Mean item complexity = 1.5

Test of the hypothesis that 2 factors are sufficient.

The degrees of freedom for the null model are 15 and the objective function was 2.5

The degrees of freedom for the model are 4 and the objective function was 0.07

The root mean square of the residuals (RMSR) is 0.03

The df corrected root mean square of the residuals is 0.06

Fit based upon off diagonal values = 0.99

Measures of factor score adequacy

PA1 PA2

Correlation of scores with factors 0.96 0.92

Multiple R square of scores with factors 0.93 0.84

Minimum correlation of possible factor scores 0.86 0.68

可以看到,两个因子解释了六个心理学测验60%的方差。不过因子载荷阵的意义并不太好解释,此时使用因子旋转有助于因子的解释。

因子旋转

我们可以使用正交旋转或者斜交旋转来旋转两个因子的结果。首先使用正交旋转

> fa.varimax <- fa(correlations,nfactors = 2,rotate = "varimax", fm="pa")

> fa.varimax

Factor Analysis using method = pa

Call: fa(r = correlations, nfactors = 2, rotate = "varimax", fm = "pa")

Standardized loadings (pattern matrix) based upon correlation matrix

PA1 PA2 h2 u2 com

general 0.49 0.57 0.57 0.432 2.0

picture 0.16 0.59 0.38 0.623 1.1

blocks 0.18 0.89 0.83 0.166 1.1

maze 0.13 0.43 0.20 0.798 1.2

reading 0.93 0.20 0.91 0.089 1.1

vocab 0.80 0.23 0.69 0.313 1.2

PA1 PA2

SS loadings 1.83 1.75

Proportion Var 0.30 0.29

Cumulative Var 0.30 0.60

Proportion Explained 0.51 0.49

Cumulative Proportion 0.51 1.00

Mean item complexity = 1.3

Test of the hypothesis that 2 factors are sufficient.

The degrees of freedom for the null model are 15 and the objective function was 2.5

The degrees of freedom for the model are 4 and the objective function was 0.07

The root mean square of the residuals (RMSR) is 0.03

The df corrected root mean square of the residuals is 0.06

Fit based upon off diagonal values = 0.99

Measures of factor score adequacy

PA1 PA2

Correlation of scores with factors 0.96 0.92

Multiple R square of scores with factors 0.91 0.85

Minimum correlation of possible factor scores 0.82 0.71

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

使用正交旋转将人为地强制两个因子不相关。如果想允许两个因子相关该如何呢?这时可用斜交旋转法,代码如下:

> fa.promax <- fa(correlations,nfactors = 2,rotate = "promax", fm="pa")

Loading required namespace: GPArotation

Warning message:

In fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, :

A loading greater than abs(1) was detected. Examine the loadings carefully.

> fa.promax

Factor Analysis using method = pa

Call: fa(r = correlations, nfactors = 2, rotate = "promax", fm = "pa")

Warning: A Heywood case was detected.

Standardized loadings (pattern matrix) based upon correlation matrix

PA1 PA2 h2 u2 com

general 0.37 0.48 0.57 0.432 1.9

picture -0.03 0.63 0.38 0.623 1.0

blocks -0.10 0.97 0.83 0.166 1.0

maze 0.00 0.45 0.20 0.798 1.0

reading 1.00 -0.09 0.91 0.089 1.0

vocab 0.84 -0.01 0.69 0.313 1.0

PA1 PA2

SS loadings 1.83 1.75

Proportion Var 0.30 0.29

Cumulative Var 0.30 0.60

Proportion Explained 0.51 0.49

Cumulative Proportion 0.51 1.00

With factor correlations of

PA1 PA2

PA1 1.00 0.55

PA2 0.55 1.00

Mean item complexity = 1.2

Test of the hypothesis that 2 factors are sufficient.

The degrees of freedom for the null model are 15 and the objective function was 2.5

The degrees of freedom for the model are 4 and the objective function was 0.07

The root mean square of the residuals (RMSR) is 0.03

The df corrected root mean square of the residuals is 0.06

Fit based upon off diagonal values = 0.99

Measures of factor score adequacy

PA1 PA2

Correlation of scores with factors 0.97 0.94

Multiple R square of scores with factors 0.93 0.88

Minimum correlation of possible factor scores 0.86 0.77

根据以上结果,我们可以看出正交旋转和斜交旋转的不同之处。对于正交旋转,因子分析的重点在于因子结构矩阵(变量与因子的相关系数),而对于斜交旋转,因子分析会考虑三个矩阵:因子结构矩阵、因子模式矩阵和因子关联矩阵。

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

上面代码中,PA1和PA2栏中的值组成了因子模式矩阵。它们是标准化的回归系数,而不是相关系数。注意,矩阵的列仍用来对因子进行命名。

因子关联矩阵显示两个因子的相关系数为0.55,相关性很大。如果因子间的关联性很低,我们可能需要重新使用正交旋转来简化问题。

因子结构矩阵(或称因子载荷矩阵)没有被列出来,但我们可以使用公式F=P*Phi很轻松地得到它,其中F是因子载荷阵,P是因子模式矩阵,Phi是因子关联矩阵。下面的函数即可进行该乘法运算:

> fsm <- function(oblique){

+ if(class(oblique)[2]=="fa" & is.null(oblique$Phi)){

+ warning("Object dosen't look like oblique EFA")

+ } else {

+ P <- unclass(oblique$loading)

+ F <- P %*% oblique$Phi

+ colnames(F) <- c("PA1","PA2")

+ return (F)

+ }

+ }

> fsm(fa.promax)

PA1 PA2

general 0.64 0.69

picture 0.32 0.61

blocks 0.43 0.91

maze 0.25 0.45

reading 0.95 0.46

vocab 0.83 0.45

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

使用factor.plot()或fa.diagram()函数,我们可以绘制正交或斜交结果的图形,以下代码绘制了斜交结果的图形:

> factor.plot(fa.promax,labels = row.names(fa.promax$loadings))

上图是数据集ability.cov中心理学测验的两因子图形。词汇和阅读在第一个因子(PA1)上载荷较大。而积木图案、画图和迷宫在第二个因子(PA2)上载荷较大。普通智力测验在两个因子上较为平均。

>fa.diagram(fa.promax,simple=FALSE)

若使simple=TRUE,那么将仅显示每个因子下最大的载荷,以及因子间的相关系数。这类图形在有多个因子时十分实用。

当处理真实生活中的数据时,我们不可能只对这么少的变量进行因子分析。这里只是为了方便。如果我们像检测自己的能力,可以尝试Harman74.cor中的24个心理学测验进行因子分析,代码如下:

>library(psych)

>fa,24tests <- fa(Harman74.cor$cov, nfactors=4.rotate="promax")

因子得分

相比PCA,EFA并不那么关注计算因子得分。在fa()函数中添加score=TRUE选项便可轻松获得因子得分。另外还可以得到得分析系数(标准化的回归权重),它在返回对象的weights元素中。

对于ability.cov数据集,通过二因子斜交旋转法便可获得用来计算因子得分的权重:

> fa.promax$weights

PA1 PA2

general 0.078 0.211

picture 0.020 0.090

blocks 0.037 0.702

maze 0.027 0.035

reading 0.743 0.030

vocab 0.177 0.036

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

最近发表
标签列表