生成碎石图，参考主成分个数

library(foreign)
library(psych)

## Warning: package 'psych' was built under R version 3.6.1

data<-read.spss("./例5-2.sav")
data

## $地区
##  [1] "北 京               " "天 津               " "河 北               "
##  [4] "山 西               " "内蒙古              " "辽 宁               "
##  [7] "吉 林               " "黑龙江              " "上 海               "
## [10] "江 苏               " "浙 江               " "安 徽               "
## [13] "福 建               " "江 西               " "山 东               "
## [16] "河 南               " "湖 北               " "湖 南               "
## [19] "广 东               " "广 西               " "海 南               "
## [22] "重 庆               " "四 川               " "贵 州               "
## [25] "云 南               " "西 藏               " "陕 西               "
## [28] "甘 肃               " "青 海               " "宁 夏               "
## [31] "新 疆               "
## 
## $X1
##  [1]  17837.50  26654.45  46906.78  12757.28  19884.38  21035.90  23412.38
##  [8]  11103.56  31056.80 155820.09  66628.47  42329.72  43309.15  32928.80
## [15] 148872.26  79404.82  47295.43  39319.29 129840.69  23406.97   1765.10
## [22]  23497.42  42103.39  11550.45  10080.39    163.73  21788.47   6527.42
## [29]   2663.51   3899.68   8105.79
## 
## $X2
##  [1]  43093.68  25075.09  44562.88  33621.95  30900.83  36106.92  18969.47
##  [8]  14951.92  39838.24 114536.32  69468.91  33563.37  32081.30  21811.92
## [15] 105046.32  60454.73  37942.33  25518.07 105604.17  16023.46   2764.18
## [22]  20214.63  41514.58  14319.98  19474.18   1110.65  30828.91  12263.36
## [29]   6143.77   8521.18  19538.65
## 
## $X3
##  [1] 19798.13 15385.02 24449.56 25579.36 19445.75 23272.90  9932.67
##  [8]  8399.97 19588.27 59466.56 38304.18 19039.87 16779.94 10549.25
## [15] 56837.87 28805.88 20355.90 13343.81 59318.72  9825.40  1538.21
## [22] 12374.58 24234.79  9074.84 12431.16   550.52 17380.70  8076.14
## [29]  4203.12  5773.39 12525.08
## 
## $X4
##  [1] 23272.45 10095.20 19977.44  8041.22 11370.20 12286.17  9022.79
##  [8]  6542.24 20004.96 54939.31 30863.45 14403.25 15169.60 11192.47
## [15] 47850.60 31347.06 17502.70 12171.52 45694.31  6185.34  1225.97
## [22]  7709.01 17167.51  5238.94  7031.04   554.21 13368.89  4187.87
## [29]  1937.98  2745.60  7005.82
## 
## $X5
##  [1]  19746.96  25888.20  47318.60  14226.45  20056.67  22038.95  23431.37
##  [8]  11347.77  34315.15 156591.04  65453.88  42190.46  42537.24  35961.32
## [15] 150641.21  79657.15  45850.64  39134.64 129151.31  22231.30   1668.92
## [22]  23467.03  41529.25  11172.44  10149.03    171.82  21027.90   7850.29
## [29]   2244.47   3646.10   8300.96
## 
## $X6
##  [1]  1608.26  2046.69  2815.11   294.78  1344.41   575.39  1268.49
##  [8]   295.54  2913.91 10574.40  4469.42  2242.26  2889.26  2443.93
## [15]  8820.02  5240.61  2713.46  2028.59  8383.04  1393.35   101.87
## [22]  1648.36  2339.82   847.02   334.98    16.94  1589.00    72.68
## [29]    80.02   143.23   386.59
## 
## $X7
##  [1]  635.55   -6.52 -229.49   92.52   29.95   77.04   95.45    7.67
##  [9]  604.70  215.01  326.89   45.46   40.51   31.24  189.58   87.93
## [17]  185.24   16.73  463.62 -163.42   10.52   56.99  -56.51   18.63
## [25]   66.62    0.56   70.87  -46.23    1.95    5.53   32.02
## 
## $X8
##  [1]  104.45  146.98  367.32  190.72  120.71  228.05  143.47  117.40
##  [9]  215.34 1111.84  690.30  330.49  421.66  269.13  905.76  724.19
## [17]  343.86  336.31 1435.86  174.62   10.87  194.95  338.15  103.41
## [25]   90.18    2.00  175.01   59.71   20.13   31.12   71.57
## 
## attr(,"label.table")
## attr(,"label.table")$地区
## NULL
## 
## attr(,"label.table")$X1
## NULL
## 
## attr(,"label.table")$X2
## NULL
## 
## attr(,"label.table")$X3
## NULL
## 
## attr(,"label.table")$X4
## NULL
## 
## attr(,"label.table")$X5
## NULL
## 
## attr(,"label.table")$X6
## NULL
## 
## attr(,"label.table")$X7
## NULL
## 
## attr(,"label.table")$X8
## NULL
## 
## attr(,"codepage")
## [1] 936

df <-as.data.frame(data)
df <- df[,2:9]
df

##           X1        X2       X3       X4        X5       X6      X7
## 1   17837.50  43093.68 19798.13 23272.45  19746.96  1608.26  635.55
## 2   26654.45  25075.09 15385.02 10095.20  25888.20  2046.69   -6.52
## 3   46906.78  44562.88 24449.56 19977.44  47318.60  2815.11 -229.49
## 4   12757.28  33621.95 25579.36  8041.22  14226.45   294.78   92.52
## 5   19884.38  30900.83 19445.75 11370.20  20056.67  1344.41   29.95
## 6   21035.90  36106.92 23272.90 12286.17  22038.95   575.39   77.04
## 7   23412.38  18969.47  9932.67  9022.79  23431.37  1268.49   95.45
## 8   11103.56  14951.92  8399.97  6542.24  11347.77   295.54    7.67
## 9   31056.80  39838.24 19588.27 20004.96  34315.15  2913.91  604.70
## 10 155820.09 114536.32 59466.56 54939.31 156591.04 10574.40  215.01
## 11  66628.47  69468.91 38304.18 30863.45  65453.88  4469.42  326.89
## 12  42329.72  33563.37 19039.87 14403.25  42190.46  2242.26   45.46
## 13  43309.15  32081.30 16779.94 15169.60  42537.24  2889.26   40.51
## 14  32928.80  21811.92 10549.25 11192.47  35961.32  2443.93   31.24
## 15 148872.26 105046.32 56837.87 47850.60 150641.21  8820.02  189.58
## 16  79404.82  60454.73 28805.88 31347.06  79657.15  5240.61   87.93
## 17  47295.43  37942.33 20355.90 17502.70  45850.64  2713.46  185.24
## 18  39319.29  25518.07 13343.81 12171.52  39134.64  2028.59   16.73
## 19 129840.69 105604.17 59318.72 45694.31 129151.31  8383.04  463.62
## 20  23406.97  16023.46  9825.40  6185.34  22231.30  1393.35 -163.42
## 21   1765.10   2764.18  1538.21  1225.97   1668.92   101.87   10.52
## 22  23497.42  20214.63 12374.58  7709.01  23467.03  1648.36   56.99
## 23  42103.39  41514.58 24234.79 17167.51  41529.25  2339.82  -56.51
## 24  11550.45  14319.98  9074.84  5238.94  11172.44   847.02   18.63
## 25  10080.39  19474.18 12431.16  7031.04  10149.03   334.98   66.62
## 26    163.73   1110.65   550.52   554.21    171.82    16.94    0.56
## 27  21788.47  30828.91 17380.70 13368.89  21027.90  1589.00   70.87
## 28   6527.42  12263.36  8076.14  4187.87   7850.29    72.68  -46.23
## 29   2663.51   6143.77  4203.12  1937.98   2244.47    80.02    1.95
## 30   3899.68   8521.18  5773.39  2745.60   3646.10   143.23    5.53
## 31   8105.79  19538.65 12525.08  7005.82   8300.96   386.59   32.02
##         X8
## 1   104.45
## 2   146.98
## 3   367.32
## 4   190.72
## 5   120.71
## 6   228.05
## 7   143.47
## 8   117.40
## 9   215.34
## 10 1111.84
## 11  690.30
## 12  330.49
## 13  421.66
## 14  269.13
## 15  905.76
## 16  724.19
## 17  343.86
## 18  336.31
## 19 1435.86
## 20  174.62
## 21   10.87
## 22  194.95
## 23  338.15
## 24  103.41
## 25   90.18
## 26    2.00
## 27  175.01
## 28   59.71
## 29   20.13
## 30   31.12
## 31   71.57

df_zscale <-scale(df)
df_zscale_cor <- cor(df_zscale )
fa.parallel(df_zscale_cor, n.obs = 112, fa = "both", n.iter = 100)

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs
## = np.obs, : The estimated weights for the factor scores are probably
## incorrect. Try a different factor extraction method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs
## = np.obs, : The estimated weights for the factor scores are probably
## incorrect. Try a different factor extraction method.

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : A loading greater than abs(1) was detected. Examine the loadings
## carefully.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs
## = np.obs, : The estimated weights for the factor scores are probably
## incorrect. Try a different factor extraction method.

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : An ultra-Heywood case was detected. Examine the results carefully

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : A loading greater than abs(1) was detected. Examine the loadings
## carefully.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs
## = np.obs, : The estimated weights for the factor scores are probably
## incorrect. Try a different factor extraction method.

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : An ultra-Heywood case was detected. Examine the results carefully

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs
## = np.obs, : The estimated weights for the factor scores are probably
## incorrect. Try a different factor extraction method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs
## = np.obs, : The estimated weights for the factor scores are probably
## incorrect. Try a different factor extraction method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs
## = np.obs, : The estimated weights for the factor scores are probably
## incorrect. Try a different factor extraction method.

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : A loading greater than abs(1) was detected. Examine the loadings
## carefully.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs
## = np.obs, : The estimated weights for the factor scores are probably
## incorrect. Try a different factor extraction method.

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : An ultra-Heywood case was detected. Examine the results carefully

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : A loading greater than abs(1) was detected. Examine the loadings
## carefully.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs
## = np.obs, : The estimated weights for the factor scores are probably
## incorrect. Try a different factor extraction method.

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : An ultra-Heywood case was detected. Examine the results carefully

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : A loading greater than abs(1) was detected. Examine the loadings
## carefully.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs
## = np.obs, : The estimated weights for the factor scores are probably
## incorrect. Try a different factor extraction method.

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : An ultra-Heywood case was detected. Examine the results carefully

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : A loading greater than abs(1) was detected. Examine the loadings
## carefully.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs
## = np.obs, : The estimated weights for the factor scores are probably
## incorrect. Try a different factor extraction method.

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : An ultra-Heywood case was detected. Examine the results carefully

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : A loading greater than abs(1) was detected. Examine the loadings
## carefully.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs
## = np.obs, : The estimated weights for the factor scores are probably
## incorrect. Try a different factor extraction method.

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : An ultra-Heywood case was detected. Examine the results carefully

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs
## = np.obs, : The estimated weights for the factor scores are probably
## incorrect. Try a different factor extraction method.

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : A loading greater than abs(1) was detected. Examine the loadings
## carefully.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs
## = np.obs, : The estimated weights for the factor scores are probably
## incorrect. Try a different factor extraction method.

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : An ultra-Heywood case was detected. Examine the results carefully

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : A loading greater than abs(1) was detected. Examine the loadings
## carefully.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs
## = np.obs, : The estimated weights for the factor scores are probably
## incorrect. Try a different factor extraction method.

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : An ultra-Heywood case was detected. Examine the results carefully

## Parallel analysis suggests that the number of factors =  1  and the number of components =  1

因子分析

fa_model4 <- fa(df_zscale_cor, nfactors = 2,rotate = "none", fm = "ml")#fa函数进行主成分分析，fa.parallel函数生成碎石图
fa_model4

## Factor Analysis using method =  ml
## Call: fa(r = df_zscale_cor, nfactors = 2, rotate = "none", fm = "ml")
## Standardized loadings (pattern matrix) based upon correlation matrix
##     ML1   ML2   h2     u2 com
## X1 0.99 -0.11 1.00 0.0027 1.0
## X2 0.98  0.17 1.00 0.0033 1.1
## X3 0.97  0.22 0.98 0.0179 1.1
## X4 0.98  0.10 0.98 0.0232 1.0
## X5 0.99 -0.11 1.00 0.0027 1.0
## X6 0.99 -0.10 0.98 0.0166 1.0
## X7 0.41  0.47 0.40 0.6033 2.0
## X8 0.96 -0.02 0.91 0.0864 1.0
## 
##                        ML1  ML2
## SS loadings           6.90 0.35
## Proportion Var        0.86 0.04
## Cumulative Var        0.86 0.91
## Proportion Explained  0.95 0.05
## Cumulative Proportion 0.95 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  28  and the objective function was  31.7
## The degrees of freedom for the model are 13  and the objective function was  8.43 
## 
## The root mean square of the residuals (RMSR) is  0.02 
## The df corrected root mean square of the residuals is  0.03 
## 
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy             
##                                                   ML1  ML2
## Correlation of (regression) scores with factors     1 0.97
## Multiple R square of scores with factors            1 0.94
## Minimum correlation of possible factor scores       1 0.87

# fa(r, nfactors=, n.obs=, rotate=, scores=, fm=) r:相关系数矩阵或原始数据矩阵， 
#   nfactors:设定主提取的因子数（默认为1） n.obs:观测数（输入相关系数矩阵时需要填写） 
#   rotate:设定旋转的方法（默认互变异数最小法） scores:设定是否需要计算因子得分（默认不需要） 
#   fm:设定因子化方法（默认极小残差法）
#   
#   提取公因子的方法（fm），方法包括： ml：最大似然法 pa：主轴迭代法 wls：加权最小二乘法 gls：广义加权最小二乘法 

factor.plot(fa_model4)

fa.diagram(fa_model4, simple = FALSE)#看出提取两个因子

因子旋转：正交旋转

#使用正交旋转将人为地强制两个因子不相关，使用斜交旋转法则允许两个因子相关。
fa_model2 <- fa(df_zscale_cor, nfactors = 2, rotate = "varimax", fm = "pa")# varimax表示旋转方式为正交因子旋转

## maximum iteration exceeded

## Warning in cor.smooth(model): Matrix was not positive definite, smoothing
## was done

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs
## = np.obs, : The estimated weights for the factor scores are probably
## incorrect. Try a different factor extraction method.

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : An ultra-Heywood case was detected. Examine the results carefully

fa_model2#系数都变了

## Factor Analysis using method =  pa
## Call: fa(r = df_zscale_cor, nfactors = 2, rotate = "varimax", fm = "pa")
## Standardized loadings (pattern matrix) based upon correlation matrix
##     PA1  PA2   h2      u2 com
## X1 0.98 0.22 1.00 -0.0044 1.1
## X2 0.88 0.47 0.99  0.0086 1.5
## X3 0.87 0.42 0.93  0.0658 1.4
## X4 0.87 0.50 1.00  0.0012 1.6
## X5 0.97 0.23 1.00  0.0013 1.1
## X6 0.94 0.31 0.97  0.0309 1.2
## X7 0.20 0.69 0.52  0.4847 1.2
## X8 0.91 0.30 0.91  0.0852 1.2
## 
##                        PA1  PA2
## SS loadings           5.92 1.41
## Proportion Var        0.74 0.18
## Cumulative Var        0.74 0.92
## Proportion Explained  0.81 0.19
## Cumulative Proportion 0.81 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  28  and the objective function was  31.7
## The degrees of freedom for the model are 13  and the objective function was  2934395 
## 
## The root mean square of the residuals (RMSR) is  0.01 
## The df corrected root mean square of the residuals is  0.01 
## 
## Fit based upon off diagonal values = 1

 factor.plot(fa_model2)#现在即可直观地看到变量和因子间的相关系数（即因子载荷矩阵）

 fa.diagram(fa_model4, simple = FALSE)#使用fa.diagram()函数绘制，如果simple=TRUE，那么将仅显示每个因子下的最大载荷，以及因子间的相关系数。

因子旋转：斜交旋转

fa.promax <- fa(df_zscale_cor, nfactors=1, rotate="promax", fm="pa")#这里nfactors=2就报错了

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs
## = np.obs, : The estimated weights for the factor scores are probably
## incorrect. Try a different factor extraction method.

fa.promax#对于正交旋转，因子分析的重点在于因子载荷矩阵（变量和因子的相关系数）；而对于斜交旋转，因子分析会考虑三个矩阵，因子载荷矩阵、因子模式矩阵和因子关联矩阵。因子模式矩阵即PA1那的值，是标准化的回归系数（而不是载荷矩阵中的相关系数）

## Factor Analysis using method =  pa
## Call: fa(r = df_zscale_cor, nfactors = 1, rotate = "promax", fm = "pa")
## Standardized loadings (pattern matrix) based upon correlation matrix
##     PA1   h2    u2 com
## X1 0.98 0.96 0.036   1
## X2 0.99 0.99 0.015   1
## X3 0.97 0.94 0.065   1
## X4 0.99 0.98 0.017   1
## X5 0.98 0.97 0.034   1
## X6 0.98 0.96 0.036   1
## X7 0.43 0.19 0.814   1
## X8 0.95 0.91 0.090   1
## 
##                 PA1
## SS loadings    6.89
## Proportion Var 0.86
## 
## Mean item complexity =  1
## Test of the hypothesis that 1 factor is sufficient.
## 
## The degrees of freedom for the null model are  28  and the objective function was  31.7
## The degrees of freedom for the model are 20  and the objective function was  13.94 
## 
## The root mean square of the residuals (RMSR) is  0.03 
## The df corrected root mean square of the residuals is  0.04 
## 
## Fit based upon off diagonal values = 1