統計検定１級の問題をRを用いて再現する。
標準偏差と相関係数が分かれば分散共分散行列が分かる。
平均と分散共分散行列が分かればデータを再現することができる。
Rによる分析を追体験しよう。

1 2013_問3

library(MASS)
mu <- c(49.81,44.81,58.97,40.17,58.39,43.75)#平均
R <- matrix(c(1,0.79,0.83,0.52,0.46,0.56,#相関行列
              0.79,1,0.8,0.41,0.41,0.53,
              0.83,0.8,1,0.48,0.42,0.54,
              0.52,0.41,0.48,1,0.66,0.68,
              0.46,0.41,0.42,0.66,1,0.66,
              0.56,0.53,0.54,0.68,0.66,1), 6, 6)
A <- matrix(c(15.21,0,0,0,0,0,
              0,10.69,0,0,0,0,
              0,0,14.79,0,0,0,
              0,0,0,16.59,0,0,
              0,0,0,0,18.84,0,
              0,0,0,0,0,16.05), 6, 6)
Sigma<- A %*% R %*% A #分散共分散行列
data <- mvrnorm(200, mu, Sigma)#200個の乱数を生成
d <- data.frame(X1=data[,1],X2=data[,2],X3=data[,3],X4=data[,4],X5=data[,5],X6=data[,6])#データフレーム化
round(cov(d), digits = 2)#分散共分散行列
##        X1     X2     X3     X4     X5     X6
## X1 216.83 117.82 186.36 131.82 126.82 142.48
## X2 117.82 107.77 126.42  72.20  79.01  94.71
## X3 186.36 126.42 225.78 128.97 125.88 150.58
## X4 131.82  72.20 128.97 302.74 221.93 207.83
## X5 126.82  79.01 125.88 221.93 357.15 225.36
## X6 142.48  94.71 150.58 207.83 225.36 280.01
round(cor(d), digits = 2)#相関行列
##      X1   X2   X3   X4   X5   X6
## X1 1.00 0.77 0.84 0.51 0.46 0.58
## X2 0.77 1.00 0.81 0.40 0.40 0.55
## X3 0.84 0.81 1.00 0.49 0.44 0.60
## X4 0.51 0.40 0.49 1.00 0.67 0.71
## X5 0.46 0.40 0.44 0.67 1.00 0.71
## X6 0.58 0.55 0.60 0.71 0.71 1.00

相関行列の固有値

eigen_r <- eigen(cor(d))#相関行列の固有値
round(eigen_r$vectors,digits=2)#固有ベクトル
##      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]
## [1,] 0.43 -0.34  0.18  0.34 -0.49  0.56
## [2,] 0.41 -0.44 -0.26 -0.16  0.71  0.22
## [3,] 0.43 -0.38  0.07  0.04 -0.20 -0.79
## [4,] 0.38  0.45  0.71  0.14  0.36 -0.02
## [5,] 0.37  0.50 -0.61  0.48  0.00 -0.06
## [6,] 0.42  0.31 -0.12 -0.78 -0.30  0.11
eigen_r$values#固有値
## [1] 3.9965884 1.0462283 0.3363842 0.2632765 0.2088208 0.1487018

分散共分散行列の固有値

#attach(d)
eigen_v <- eigen(cov(d))#分散共分散行列の固有値
round(eigen_v$vectors,digits=2)#固有ベクトル
##      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]
## [1,] 0.37 -0.49 -0.02 -0.23  0.76  0.01
## [2,] 0.23 -0.35 -0.12  0.02 -0.32 -0.84
## [3,] 0.37 -0.52 -0.05 -0.05 -0.55  0.53
## [4,] 0.46  0.31  0.75 -0.33 -0.10 -0.08
## [5,] 0.50  0.50 -0.64 -0.31 -0.03  0.02
## [6,] 0.47  0.14  0.05  0.86  0.13  0.04
eigen_v$values#固有値
## [1] 1000.91916  246.98496  105.53091   74.92314   35.25651   26.65175

理論上の固有値

eigen_s <- eigen(Sigma)#理論上の固有値
round(eigen_s$vectors,digits=2)
##       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]
## [1,] -0.40 -0.48 -0.07 -0.15  0.73 -0.23
## [2,] -0.25 -0.34 -0.12  0.07 -0.09  0.89
## [3,] -0.37 -0.51 -0.08 -0.08 -0.67 -0.37
## [4,] -0.44  0.30  0.66 -0.52 -0.05  0.11
## [5,] -0.50  0.53 -0.67 -0.13 -0.02 -0.02
## [6,] -0.44  0.16  0.30  0.83  0.04 -0.08
eigen_s$values
## [1] 949.35979 247.30422 105.94505  83.65828  37.95727  27.91589

これは問題文のものと一致している。

2 主成分分析

2.1 分散共分散行列による主成分分析

library("psych")
#(describe(d))#基本統計情報
par(family = "HiraKakuProN-W3") #日本語フォントの指定
result <- prcomp(d, scale=F) #scale=Tは相関行列から分析
summary(result)
## Importance of components:
##                            PC1     PC2      PC3     PC4     PC5     PC6
## Standard deviation     31.6373 15.7158 10.27282 8.65582 5.93772 5.16253
## Proportion of Variance  0.6716  0.1657  0.07081 0.05027 0.02366 0.01788
## Cumulative Proportion   0.6716  0.8374  0.90818 0.95846 0.98212 1.00000
biplot(result)

biplot からは、変数が２つのカテゴリーに分かれることが読み取れる。
第１主成分得点
基準化されていることがわかる。

describe(result$x[,1])
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 200 0 31.63731 2.579019 0.9356176 28.87393 -84.80238 87.12533 171.9277 -0.223019 -0.0086125 2.237095 
#hist(result$x[,1])
sd(result$x[,1])*0.84
## [1] 26.57534

上位２０％での基準化された得点は、ほぼ理論値に近い。

2.2 相関行列による主成分分析

result <- prcomp(d, scale=T) 
eigen_r <- eigen(cor(d))#相関行列の固有値
round(eigen_r$vectors,digits=2)
##      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]
## [1,] 0.43 -0.34  0.18  0.34 -0.49  0.56
## [2,] 0.41 -0.44 -0.26 -0.16  0.71  0.22
## [3,] 0.43 -0.38  0.07  0.04 -0.20 -0.79
## [4,] 0.38  0.45  0.71  0.14  0.36 -0.02
## [5,] 0.37  0.50 -0.61  0.48  0.00 -0.06
## [6,] 0.42  0.31 -0.12 -0.78 -0.30  0.11
eigen_r$values
## [1] 3.9965884 1.0462283 0.3363842 0.2632765 0.2088208 0.1487018

標準化したデータを固有ベクトル方向に固有値倍に拡大するので、
第１主成分の分散は第１固有値に等しい。

describe(result$x[,1])
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 200 0 1.999147 0.2777244 0.0813763 1.981412 -5.396052 5.398751 10.7948 -0.3198059 -0.1837341 0.141361 
#hist(result$x[,1])

主成分分析を因子分析に進めるには、主成分得点を標準化する必要がある。固有ベクトル（rotation）に標準偏差（sdev）をかけたものが主成分負荷量になる。

sum(eigen_r$vectors[,1]^2)
## [1] 1
fc.l <- sweep(result$rotation, MARGIN=2, result$sdev, FUN="*")
round(fc.l,digits=3)
##      PC1    PC2    PC3    PC4    PC5    PC6
## X1 0.859 -0.353  0.107 -0.173 -0.222 -0.216
## X2 0.811 -0.449 -0.151  0.081  0.325 -0.083
## X3 0.866 -0.384  0.041 -0.021 -0.093  0.305
## X4 0.767  0.459  0.412 -0.070  0.163  0.009
## X5 0.741  0.512 -0.355 -0.248  0.000  0.025
## X6 0.845  0.316 -0.069  0.400 -0.138 -0.043

2.3 principalによる因子分析

因子分析をprincipalを用いて行えば、主成分負荷量は自動で計算される。

library(psych)
dat <- d
fit <- principal(dat, nfactors=2, rotate=FALSE)#回転ない
## Specified rotation not found, rotate='none' used
fit
## Principal Components Analysis
## Call: principal(r = dat, nfactors = 2, rotate = FALSE)
## Standardized loadings (pattern matrix) based upon correlation matrix
##     PC1   PC2   h2   u2 com
## X1 0.86 -0.35 0.86 0.14 1.3
## X2 0.81 -0.45 0.86 0.14 1.6
## X3 0.87 -0.38 0.90 0.10 1.4
## X4 0.77  0.46 0.80 0.20 1.6
## X5 0.74  0.51 0.81 0.19 1.8
## X6 0.85  0.32 0.81 0.19 1.3
## 
##                        PC1  PC2
## SS loadings           4.00 1.05
## Proportion Var        0.67 0.17
## Cumulative Var        0.67 0.84
## Proportion Explained  0.79 0.21
## Cumulative Proportion 0.79 1.00
## 
## Mean item complexity =  1.5
## Test of the hypothesis that 2 components are sufficient.
## 
## The root mean square of the residuals (RMSR) is  0.05 
##  with the empirical chi square  17.65  with prob <  0.0014 
## 
## Fit based upon off diagonal values = 0.99
summary(fit)
## 
## Factor analysis with Call: principal(r = dat, nfactors = 2, rotate = FALSE)
## 
## Test of the hypothesis that 2 factors are sufficient.
## The degrees of freedom for the model is 4  and the objective function was  0.44 
## The number of observations was  200  with Chi Square =  84.97  with prob <  1.5e-17 
## 
## The root mean square of the residuals (RMSA) is  0.05
biplot(fit)

因子負荷量は主成分負荷量とほぼ一致している。

3 farによる因子分析

3.1 none回転

fit <- fa(r=dat, nfactors=2 ,rotate="", fm="ml", scores=T)
## Specified rotation not found, rotate='none' used
print(fit, digits=5)# 因子得点を出すためにはscores=T
## Factor Analysis using method =  ml
## Call: fa(r = dat, nfactors = 2, rotate = "", scores = T, fm = "ml")
## Standardized loadings (pattern matrix) based upon correlation matrix
##        ML1      ML2      h2      u2    com
## X1 0.87684 -0.17941 0.80104 0.19896 1.0836
## X2 0.82819 -0.24151 0.74423 0.25577 1.1689
## X3 0.91188 -0.23452 0.88653 0.11347 1.1317
## X4 0.66309  0.48637 0.67624 0.32376 1.8345
## X5 0.62935  0.53617 0.68356 0.31644 1.9508
## X6 0.76442  0.42651 0.76625 0.23375 1.5676
## 
##                           ML1     ML2
## SS loadings           3.70639 0.85146
## Proportion Var        0.61773 0.14191
## Cumulative Var        0.61773 0.75964
## Proportion Explained  0.81319 0.18681
## Cumulative Proportion 0.81319 1.00000
## 
## Mean item complexity =  1.5
## Test of the hypothesis that 2 factors are sufficient.
## 
## df null model =  15  with the objective function =  4.46551 with Chi Square =  875.9842
## df of  the model are 4  and the objective function was  0.03092 
## 
## The root mean square of the residuals (RMSR) is  0.01178 
## The df corrected root mean square of the residuals is  0.02282 
## 
## The harmonic n.obs is  200 with the empirical chi square  0.83319  with prob <  0.93394 
## The total n.obs was  200  with Likelihood Chi Square =  6.02441  with prob <  0.19733 
## 
## Tucker Lewis Index of factoring reliability =  0.991121
## RMSEA index =  0.050054  and the 90 % confidence intervals are  0 0.127096
## BIC =  -15.16885
## Fit based upon off diagonal values = 0.99963
## Measures of factor score adequacy             
##                                                       ML1     ML2
## Correlation of (regression) scores with factors   0.97466 0.87579
## Multiple R square of scores with factors          0.94996 0.76700
## Minimum correlation of possible factor scores     0.89992 0.53401
biplot(fit)

fa.diagram( fit )

3.2 promax回転

fit.ml.promax <- fa(r=dat, nfactors=2 ,rotate="promax", fm="ml", scores=T)# 因子得点を出すためにはscores=T
##  要求されたパッケージ GPArotation をロード中です
print(fit.ml.promax, digits=5)
## Factor Analysis using method =  ml
## Call: fa(r = dat, nfactors = 2, rotate = "promax", scores = T, fm = "ml")
## Standardized loadings (pattern matrix) based upon correlation matrix
##         ML1      ML2      h2      u2    com
## X1  0.85571  0.06089 0.80104 0.19896 1.0101
## X2  0.88219 -0.03179 0.74423 0.25577 1.0026
## X3  0.93911  0.00391 0.88653 0.11347 1.0000
## X4  0.00870  0.81687 0.67624 0.32376 1.0002
## X5 -0.06825  0.86769 0.68356 0.31644 1.0124
## X6  0.14775  0.77543 0.76625 0.23375 1.0725
## 
##                           ML1     ML2
## SS loadings           2.47533 2.08253
## Proportion Var        0.41255 0.34709
## Cumulative Var        0.41255 0.75964
## Proportion Explained  0.54309 0.45691
## Cumulative Proportion 0.54309 1.00000
## 
##  With factor correlations of 
##         ML1     ML2
## ML1 1.00000 0.62467
## ML2 0.62467 1.00000
## 
## Mean item complexity =  1
## Test of the hypothesis that 2 factors are sufficient.
## 
## df null model =  15  with the objective function =  4.46551 with Chi Square =  875.9842
## df of  the model are 4  and the objective function was  0.03092 
## 
## The root mean square of the residuals (RMSR) is  0.01178 
## The df corrected root mean square of the residuals is  0.02282 
## 
## The harmonic n.obs is  200 with the empirical chi square  0.83319  with prob <  0.93394 
## The total n.obs was  200  with Likelihood Chi Square =  6.02441  with prob <  0.19733 
## 
## Tucker Lewis Index of factoring reliability =  0.991121
## RMSEA index =  0.050054  and the 90 % confidence intervals are  0 0.127096
## BIC =  -15.16885
## Fit based upon off diagonal values = 0.99963
## Measures of factor score adequacy             
##                                                       ML1     ML2
## Correlation of (regression) scores with factors   0.96867 0.94049
## Multiple R square of scores with factors          0.93832 0.88451
## Minimum correlation of possible factor scores     0.87664 0.76903
biplot(fit.ml.promax)

fa.diagram( fit.ml.promax)

3.3 varimax回転

fit.ml.varimax <- fa(r=dat, nfactors=2 ,rotate="varimax", fm="ml", scores=T)
print(fit.ml.varimax, digits=5)
## Factor Analysis using method =  ml
## Call: fa(r = dat, nfactors = 2, rotate = "varimax", scores = T, fm = "ml")
## Standardized loadings (pattern matrix) based upon correlation matrix
##        ML1     ML2      h2      u2    com
## X1 0.83023 0.33431 0.80104 0.19896 1.3160
## X2 0.82393 0.25567 0.74423 0.25577 1.1908
## X3 0.88986 0.30769 0.88653 0.11347 1.2358
## X4 0.28454 0.77155 0.67624 0.32376 1.2671
## X5 0.22892 0.79445 0.68356 0.31644 1.1649
## X6 0.40207 0.77755 0.76625 0.23375 1.4991
## 
##                           ML1     ML2
## SS loadings           2.45502 2.10284
## Proportion Var        0.40917 0.35047
## Cumulative Var        0.40917 0.75964
## Proportion Explained  0.53863 0.46137
## Cumulative Proportion 0.53863 1.00000
## 
## Mean item complexity =  1.3
## Test of the hypothesis that 2 factors are sufficient.
## 
## df null model =  15  with the objective function =  4.46551 with Chi Square =  875.9842
## df of  the model are 4  and the objective function was  0.03092 
## 
## The root mean square of the residuals (RMSR) is  0.01178 
## The df corrected root mean square of the residuals is  0.02282 
## 
## The harmonic n.obs is  200 with the empirical chi square  0.83319  with prob <  0.93394 
## The total n.obs was  200  with Likelihood Chi Square =  6.02441  with prob <  0.19733 
## 
## Tucker Lewis Index of factoring reliability =  0.991121
## RMSEA index =  0.050054  and the 90 % confidence intervals are  0 0.127096
## BIC =  -15.16885
## Fit based upon off diagonal values = 0.99963
## Measures of factor score adequacy             
##                                                       ML1     ML2
## Correlation of (regression) scores with factors   0.94564 0.90704
## Multiple R square of scores with factors          0.89423 0.82273
## Minimum correlation of possible factor scores     0.78847 0.64546
biplot(fit.ml.varimax)

fa.diagram( fit.ml.varimax)

3.4 最小残差法minres

rotate = “oblimin”

result = fa( dat, nfactors = 2, fm = "minres", rotate = "oblimin", use = "complete.obs" )
print( result$loadings, digits = 2, cutoff = 0.3 )
## 
## Loadings:
##    MR1   MR2  
## X1  0.85      
## X2  0.89      
## X3  0.94      
## X4        0.83
## X5        0.87
## X6        0.78
## 
##                 MR1  MR2
## SS loadings    2.43 2.06
## Proportion Var 0.40 0.34
## Cumulative Var 0.40 0.75
biplot(result)

fa.diagram( result)