【R实验.9】主成分和因子分析

小鱼儿 2022-10-16 15:19 45阅读 0赞

*解法并不单一，下列方法带有璇子个人的偏好，因此仅供参考。如有错误，欢迎在评论区斧正！*

9.1 用主成分方法探讨城市工业主体结构。表 9-4 是某市工业部门十三个行业，分别是冶金 （1）、电力（2）、煤炭（3）、化学（4）、机械（5）、建材（6）、森工（7）、食品（8）、 纺织（9）、缝纫（10）、皮革（11）、造纸（12）和文教艺术用品（13），八个指标，分别 是年未固定资产净值 X1（万元）、职工人数 X2（人）、工业总产值 X3（万元）、全员劳动生 产率 X4（元/人年）、百元固定原值实现产值 X5（元）、资金利税率 X6（%）、标准燃料消费 量 X7（吨）和能源利用效果 X8（万元/吨）的数据。  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70_pic_center]

（1）试用主成分分析方法确定 8 个指标的几个主成分，并对主成分进行解释；

#导入数据
    >data9.1 <- read.table("data9.1.txt",header=FALSE)
    > library(psych)
    > #碎石图-初步预计主成分个数
    > fa.parallel(data9.1[,], fa="pc", n.iter = 100, main = "碎石图")

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70]  
由kaiser-Harris准则，大致需要保留前两个主成分

#提取主成分
    > pc <- principal(data9.1[,], nfactors = 2, covar = FALSE, rotate = "none")> pc
    Principal Components Analysis
    Call: principal(r = data9.1[, ], nfactors = 2, rotate = "none", covar = FALSE)
    Standardized loadings (pattern matrix) based upon correlation matrix
         PC1   PC2   h2    u2 com
    V1  0.84  0.50 0.96 0.041 1.6
    V2  0.83  0.47 0.92 0.082 1.6
    V3  0.75  0.64 0.97 0.028 2.0
    V4 -0.38  0.77 0.73 0.269 1.5
    V5 -0.68  0.56 0.79 0.214 1.9
    V6 -0.62  0.69 0.86 0.144 2.0
    V7  0.38 -0.64 0.56 0.444 1.6
    V8  0.10  0.46 0.22 0.775 1.1
    
                           PC1  PC2
    SS loadings           3.10 2.90
    Proportion Var        0.39 0.36
    Cumulative Var        0.39 0.75
    Proportion Explained  0.52 0.48
    Cumulative Proportion 0.52 1.00
    
    Mean item complexity =  1.7
    Test of the hypothesis that 2 components are sufficient.
    
    The root mean square of the residuals (RMSR) is  0.09 
     with the empirical chi square  5.62  with prob <  0.96 
    Fit based upon off diagonal values = 0.96

从代码PC1和PC2栏可以看到，第一主成分解释了城市工业主体结构39%的方差，而第二主成分解释了36%，两者共解释了75%的方差，对于一个宏观经济问题算是不错了。

> #当提取了多个成分时，为使结果更具解释性，故进行如下旋转
    > pc <- principal(data9.1[,], nfactors = 2, covar = FALSE, rotate = "varimax")
    > pc
    Principal Components Analysis
    Call: principal(r = data9.1[, ], nfactors = 2, rotate = "varimax", 
        covar = FALSE)
    Standardized loadings (pattern matrix) based upon correlation matrix
         RC1   RC2   h2    u2 com
    V1  0.98 -0.09 0.96 0.041 1.0
    V2  0.95 -0.11 0.92 0.082 1.0
    V3  0.98  0.08 0.97 0.028 1.0
    V4  0.15  0.84 0.73 0.269 1.1
    V5 -0.22  0.86 0.79 0.214 1.1
    V6 -0.09  0.92 0.86 0.144 1.0
    V7 -0.08 -0.74 0.56 0.444 1.0
    V8  0.35  0.32 0.22 0.775 2.0
    
                           RC1  RC2
    SS loadings           3.03 2.97
    Proportion Var        0.38 0.37
    Cumulative Var        0.38 0.75
    Proportion Explained  0.51 0.49
    Cumulative Proportion 0.51 1.00

旋转后两成分的累积方差解释性没有变化（75%），变的只是各个主成分对方差的结解释度：第一主成分解释了城市工业主体结构38%的方差，而第二主成分解释了37%。

pc1 = 0.98x1 + 0.95x2 + 0.98x3 + 0.15x4 - 0.22x5 - 0.09x6 - 0.08x7 + 0.35x8

Pc2 = -0.09x1 - 0.11x2 + 0.08x3 + 0.84x4 + 0.86x5 + 0.92x6 - 0.74x7 + 0.32x8

由上，可解释主成分一为工业规模，主成分二为产出效益。

注：上述问题也可以选取三个主成分，因为题目并未明确信息损失率的最大限额。  
我在此处选用两个主成分，因为我认为70%以上的方差解释率对于一个宏观问题的研究已经足够了。

（3）利用主成分得分对 13 个行业进行排序和分类。  
利用上述公式代入数值即可得到主成分得分  
比如计算冶金行业第一主成分得分：

> pc1 = 0.98*90342 + 0.95*52455 + 0.98*101091 + 0.15*19272 - 0.22*82.0 - 0.09*16.1  - 0.08*197435 + 0.35*0.172
    > pc1
    [1]224513.2
    
    其他的方法类似，在此省略
    
    #获取主成分得分
    > pc <- principal(data9.1[,], nfactors = 2, scores = TRUE)
    > pc$scores
                  RC1        RC2
     [1,]  0.93850029 -0.1370301
     [2,] -0.67434595 -1.3940307
     [3,] -0.63978735 -1.8817249
     [4,]  0.62224308  0.4055288
     [5,]  2.85780141 -0.4520175
     [6,] -0.46651650 -0.9503237
     [7,] -0.61469825  0.2202598
     [8,] -0.22232063  1.8001172
     [9,] -0.07279029  0.7067208
    [10,] -0.64167788  1.0737326
    [11,] -0.58912765 -0.1177322
    [12,] -0.53868266  0.4183384
    [13,]  0.04140237  0.3081615

9.2 对某地区的某类消费品的销售量 Y 进行调查，它与下面 4 个变量有关：X1——居民可支配 收入，X2——该类消费品评价价格指数，X3——社会该消费品保有量，X4——其他消费品平均 价格指数，历史资料如表 9-5 所示。试利用主成分回归方法建立销售量与 4 个变量 X1，X2，X3，X4的回归方程。  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70_pic_center 1]

1.先找主成分

>data9.2 <- read.table("data9.2.txt")
    > library(psych)
    > #碎石图-初步预测主成分个数
    > fa.parallel(data9.2[,-5], fa="pc", n.iter = 100, main = "碎石图")

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70 1]  
初步预计，一个主成分即可（kaiser-harris准则）

pc <- principal(data9.2[,-5], nfactors = 1, covar = FALSE, rotate = "none")
    > pc
    Principal Components Analysis
    Call: principal(r = data9.2[, -5], nfactors = 1, rotate = "none", covar = FALSE)
    Standardized loadings (pattern matrix) based upon correlation matrix
        PC1   h2     u2 com
    V1 1.00 0.99 0.0078   1
    V2 0.99 0.99 0.0149   1
    V3 0.99 0.98 0.0221   1
    V4 0.99 0.99 0.0114   1
    
                    PC1
    SS loadings    3.94
    Proportion Var 0.99
    
    Proportion Var=0.99，预期正确，第一主成分表示为
    
    PC1 = x1 + 0.99(x2+x3+x4)，带入原始数据计算得向量z为 
    > x1 <- data9.2$V1
    > x2 <- data9.2$V2; x3 <- data9.2$V3; 
    > x4 <- data9.2$V4
    > y <- data9.2$V5
    > z <- x1 + 0.99*(x2+x3+x4)#R本身适用于向量运算
    
    2.主成分回归
    >fit <- lm(y ~ z)
    > summary(fit)
    Call:
    lm(formula = y ~ z)
    
    Residuals:
        Min      1Q  Median      3Q     Max 
    -0.6233 -0.3120 -0.0141  0.3308  0.5571 
    
    Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
    (Intercept) -11.049619   0.809174  -13.65 7.96e-07 ***
    z             0.068333   0.002176   31.41 1.15e-09 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Residual standard error: 0.4132 on 8 degrees of freedom
    Multiple R-squared:  0.992,	Adjusted R-squared:  0.9909 
    F-statistic: 986.4 on 1 and 8 DF,  p-value: 1.149e-09
    由上，解释变量z通过t检验，总体模型通过F检验，且调整后拟合优度>0.99，模型解释力度较好。表达式为：
    Y = -11.049619 + 0.068333

9.3 为考查学生的学习情况，学校随机抽取 12 名学生的 5 门课期末考试的成绩，其数据见表 9-6，试用因子分析的方法对这种数据进行分析。  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70_pic_center 2]

（1）找出 5 门课程的公共因子，并进行合理的解释；  
1.判断需要提取的公共因子数

data9.3 <- read.table("data9.3.txt",header=T)
    > library(psych)
    > #因子图形同时展示主成分和公共因子分析结果
    > fa.parallel(data9.3[,-1], n.obs = 112, fa="both", n.iter = 100)

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70 2]  
观察EFA结果，显然需要提取两个因子。碎石检验的前两个特征值（三角形)都在拐角处之上，并且大于基于100次模拟数据矩阵的特征值均值。（对EFA，kaiser-harris准则特征值数大于0，而不是1）

3.提取公共因子

#使用主轴迭代法（fm="pa）提取未旋转因子
    > fa <- fa(data9.3[,-1], nfactors = 2, rotate = "none", fm="pa")
    > fa
    Factor Analysis using method =  pa
    Call: fa(r = data9.3[, -1], nfactors = 2, rotate = "none", fm = "pa")
    Standardized loadings (pattern matrix) based upon correlation matrix
          PA1   PA2   h2     u2 com
    政治 0.82 -0.65 1.09 -0.092 1.9
    语文 0.84 -0.30 0.80  0.199 1.2
    外语 0.68  0.02 0.46  0.535 1.0
    数学 0.90  0.55 1.11 -0.108 1.7
    物理 0.62  0.44 0.58  0.419 1.8
    
                           PA1  PA2
    SS loadings           3.04 1.01
    Proportion Var        0.61 0.20
    Cumulative Var        0.61 0.81
    Proportion Explained  0.75 0.25
    Cumulative Proportion 0.75 1.00

可以看到，两个因子解释了期末成绩81%的方差。不过因子载荷的意义不好解释，遂使用因子旋转

fa.varimax <- fa(data9.3[,-1], nfactors = 2, rotate = "varimax", fm="pa")
    > fa.varimax
    Factor Analysis using method =  pa
    Call: fa(r = data9.3[, -1], nfactors = 2, rotate = "varimax", fm = "pa")
    Standardized loadings (pattern matrix) based upon correlation matrix
          PA1  PA2   h2     u2 com
    政治 1.04 0.10 1.09 -0.092 1.0
    语文 0.81 0.37 0.80  0.199 1.4
    外语 0.48 0.49 0.46  0.535 2.0
    数学 0.27 1.02 1.11 -0.108 1.1
    物理 0.14 0.75 0.58  0.419 1.1
    
                           PA1  PA2
    SS loadings           2.06 1.98
    Proportion Var        0.41 0.40
    Cumulative Var        0.41 0.81
    Proportion Explained  0.51 0.49
    Cumulative Proportion 0.51 1.00

结果显示：政治和语文在第一因子上载荷较大，物理和数学在第二因子上载荷较大，外语在两因子上载荷较为平均，这表明存在一个语言智力因子和一个逻辑智力因子。

（2）用回归方法或 Bartlett 方法计算样本的因子得分，画出因子得分的第 1、第 2 公共因子的散点图，通过这些散点图来分析这 12 名学生的学习情况。  
由上两公因子可表示为：  
Pa1 = 1.04*x1 + 0.81*x2 + 0.48*x3 + 0.27*x4 + 0.14*x5  
Pa2 = 0.1*x1 + 0.37*x2 + 0.49*x3 + 1.02*x4 + 0.75*x5

#计算因子得分
    x1<-data9.3[,2]; x2<-data9.3[,3]; x3 <- data9.3[,4]; x4 <- data9.3[,5]; x5 <- data9.3[,6]
    
    > Pa1 = 1.04*x1 + 0.81*x2 + 0.48*x3 + 0.27*x4 + 0.14*x5
    > Pa2 = 0.1*x1 + 0.37*x2 + 0.49*x3 + 1.02*x4 + 0.75*x5
    > par(mfrow=c(1,2))
    > plot(num <- data9.3[,1],Pa1)
    > plot(num <- data9.3[,1],Pa2)

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70 3]  
我们可以发现所抽取12名学生的学习情况在俩因子上是趋同的，即公因子1得分高的，公因子2也倾向于高得分。（有些类似于高一文理未分班时的情景）

#公共因子散点图
     factor.plot(fa.varimax,labels = rownames(fa.varimax$loadings))

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70 4]  
由图发现两种因子具有反向变动关系（此消彼长），从侧面也表现了文理分班教学的合理性。

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70_pic_center]: /images/20221014/b99533034d86478ebcb2b85ec902a159.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70]: /images/20221014/a5e8fb0921e04ea29af38e4c34597e22.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70_pic_center 1]: /images/20221014/d798a13ab4e0414e99a44e960e6a8418.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70 1]: /images/20221014/3b531c7c3e4049fbb6e29f4f6aac6d80.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70_pic_center 2]: /images/20221014/741b1c9f0301439d933452e9f793bdd1.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70 2]: /images/20221014/eb96b25161664a17b2b65b9716e07b3b.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70 3]: /images/20221014/190f1a401bba4055917726d8f85a56a9.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70 4]: /images/20221014/e246866b65104c50911406b7e0096e41.png