『入門・社会統計学――2ステップで基礎から〔Rで〕学ぶ』
第11章 主成分分析と因子分析
SUGINO Isamu
2020-10-21 02:17:37, Build with R4.0.2
0 全体の章構成
1 主成分分析(Principal Component Analysis)
主成分分析(PCA; Principal Component Analysis)のモデルを図示すると以下の様になる。観測された変数を元に，主成分得点を合成するのである。
調査で実際に回答などの形で情報（数値）を得られている観測変数を var1 などと表記し，□で表す。それに対して，主成分(principal component)や因子(factor)など，直接観測されていない潜在変数を○で表す。
主成分分析のモデル
1-1 データが有する情報量の次元の縮約
3つの 主成分を合成するモデル
\[PC_{1}=h_{11}x_{1}+h_{12}x_{2}+...+h_{1j}x_{j}+...+h_{1p}x_{p}\] \[PC_{2}=h_{21}x_{1}+h_{22}x_{2}+...+h_{2j}x_{j}+...+h_{2p}x_{p}\] \[PC_{3}=h_{31}x_{1}+h_{32}x_{2}+...+h_{3j}x_{j}+...+h_{3p}x_{p}\] データフレイム名はd01とする。まずは使用する変数の度数分布から，欠損値処理が必要かどうかを確認し，必要なものは処理をする。また，頻度が高い方が数値が大きくなるよう，必要があれば変換する
d01 <- read.csv("practice.csv")# ここでは結果は省略する
table(d01$q0101, useNA="always") # 1~4
table(d01$q0102, useNA="always") # 1~5
table(d01$q0103, useNA="always") # 1~8
table(d01$q0200, useNA="always") # 1~6
table(d01$q0301, useNA="always") # 1~5
table(d01$q0302, useNA="always") # 1~5
table(d01$q0401, useNA="always") # 1~5
table(d01$q0402, useNA="always") # 1~5
table(d01$q0403, useNA="always") # 1~7# 欠損値処理と値の逆転
d01$q0103 <- c(8:1)[d01$q0103]
d01$q0200 <- c(1:6)[d01$q0200]
d01$q0301 <- c(5:1)[d01$q0301]
d01$q0302 <- c(5:1)[d01$q0302]
d01$q0402 <- c(5:1)[d01$q0402]次に，使用する変数だけ取り出してオブジェクトvarsに纏めておき，そこから完備ケース分析用にケース選択したvarsを作成する。
vars0 <- with(d01, cbind(q0101, q0102, q0103, q0200, q0301, q0302, q0401, q0402, q0403))
vars <- vars0[complete.cases(vars0),] # 完備ケース分析の準備こうすれば主成分分析は簡単に出来る。prcomp( )
(pr1 <- prcomp(vars, scale=T)) # 主成分分析を実行して結果を表示Standard deviations (1, .., p=9):
[1] 1.3480341 1.2762922 1.0299375 0.9838979 0.9363329 0.8974959 0.8434077
[8] 0.7614977 0.7427129
Rotation (n x k) = (9 x 9):
PC1 PC2 PC3 PC4 PC5 PC6
q0101 0.49802484 -0.05594720 0.04395227 -0.24386147 0.10064407 -0.60884535
q0102 0.55462152 -0.11620220 -0.06292039 -0.13546901 0.04559114 0.33674794
q0103 0.51380170 -0.14901883 -0.37208178 0.02327602 0.09483285 0.19862640
q0200 0.19839061 0.50629241 0.12871992 0.15505025 0.37457317 -0.35694017
q0301 0.15476394 0.45689311 0.30059386 0.17007928 -0.44947710 -0.09155721
q0302 0.25179900 -0.09923688 0.72421757 0.28591804 0.14384478 0.39301714
q0401 -0.14936189 0.45042646 -0.10294634 -0.20288696 0.64220223 0.30941826
q0402 0.15650630 0.37518387 -0.45674402 0.52435071 -0.24006680 0.16487432
q0403 -0.08280943 -0.37760050 -0.07095590 0.68615925 0.38172915 -0.25155723
PC7 PC8 PC9
q0101 0.01426529 0.5295311033 -0.16128354
q0102 0.18406205 -0.0005812592 0.71219681
q0103 0.09917056 -0.4259953369 -0.57884543
q0200 -0.34882468 -0.4825835005 0.20352623
q0301 0.64768373 -0.1063191634 -0.08247915
q0302 -0.24408349 0.1866299324 -0.22571622
q0401 0.33870309 0.2873705693 -0.13259326
q0402 -0.30042330 0.4194222320 0.01411737
q0403 0.38776191 0.0246223850 0.12143918scale=Tというオプションは元の変数を標準化するというオプションである。標準化の関数がscale( )であったことを思い出そう。prcomp( )関数のデフォルトはscale=Fであり，明示的に指定しないと標準化されない。標準化をしないと元の変数の単位や散布度の違いによって結果が影響される。元の変数の散布度の違いも重要であるならばそれでも良いが，ここでの分析例のようにそもそも単位がまったく違ったり散布度の違いが重要でない場合は，明示的にscale=Tとすることが必要である。scale=Fが分散共分散行列の固有値分解に，scale=Tが相関係数行列の固有値分解に対応する。
names(pr1) # 主成分分析の結果に含まれる情報名を確認[1] "sdev" "rotation" "center" "scale" "x" pr1$sdev^2 # 標準偏差を二乗して分散（＝固有値）を求める。[1] 1.8171959 1.6289217 1.0607713 0.9680552 0.8767194 0.8054988 0.7113366
[8] 0.5798788 0.5516225eigen(cor(vars)) # 相関係数行列を固有値分解すると，主成分分析の結果と一致する。eigen() decomposition
$values
[1] 1.8171959 1.6289217 1.0607713 0.9680552 0.8767194 0.8054988 0.7113366
[8] 0.5798788 0.5516225
$vectors
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] -0.49802484 0.05594720 0.04395227 -0.24386147 0.10064407 0.60884535
[2,] -0.55462152 0.11620220 -0.06292039 -0.13546901 0.04559114 -0.33674794
[3,] -0.51380170 0.14901883 -0.37208178 0.02327602 0.09483285 -0.19862640
[4,] -0.19839061 -0.50629241 0.12871992 0.15505025 0.37457317 0.35694017
[5,] -0.15476394 -0.45689311 0.30059386 0.17007928 -0.44947710 0.09155721
[6,] -0.25179900 0.09923688 0.72421757 0.28591804 0.14384478 -0.39301714
[7,] 0.14936189 -0.45042646 -0.10294634 -0.20288696 0.64220223 -0.30941826
[8,] -0.15650630 -0.37518387 -0.45674402 0.52435071 -0.24006680 -0.16487432
[9,] 0.08280943 0.37760050 -0.07095590 0.68615925 0.38172915 0.25155723
[,7] [,8] [,9]
[1,] 0.01426529 -0.5295311033 0.16128354
[2,] 0.18406205 0.0005812592 -0.71219681
[3,] 0.09917056 0.4259953369 0.57884543
[4,] -0.34882468 0.4825835005 -0.20352623
[5,] 0.64768373 0.1063191634 0.08247915
[6,] -0.24408349 -0.1866299324 0.22571622
[7,] 0.33870309 -0.2873705693 0.13259326
[8,] -0.30042330 -0.4194222320 -0.01411737
[9,] 0.38776191 -0.0246223850 -0.121439181-2 主成分の選出
主成分分析の結果pr1のsummary( )と，固有値のグラフ化を行う。本文よりも少しオプションを付けておこう。
summary(pr1) # 主成分分析の結果pr1の要約Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7
Standard deviation 1.3480 1.2763 1.0299 0.9839 0.93633 0.8975 0.84341
Proportion of Variance 0.2019 0.1810 0.1179 0.1076 0.09741 0.0895 0.07904
Cumulative Proportion 0.2019 0.3829 0.5008 0.6083 0.70574 0.7952 0.87428
PC8 PC9
Standard deviation 0.76150 0.74271
Proportion of Variance 0.06443 0.06129
Cumulative Proportion 0.93871 1.00000plot(pr1, main="主成分分析の固有値プロット", type="l")
1-3 主成分と元の変数の関係の解釈
主成分得点
prcomp( )の主成分得点
主成分分析の結果をオブジェクトpr1に格納してあるが，主成分得点はpr1$xで取り出せる。主成分得点は最大の9つまで算出されており，それが有効ケース数だけある。pr1$xをすべて表示させると長大になるので，head( )関数を使って最初の6ケース分だけ表示してみる。
head(pr1$x) # 主成分得点，先頭6ケースのみ PC1 PC2 PC3 PC4 PC5 PC6
[1,] -2.1456084 0.3557653 0.46480739 -0.37080957 0.12638217 0.8024176
[2,] 1.4349870 1.1018938 2.22131833 -0.03484862 -0.01421889 -0.5180879
[3,] -2.1404049 0.5573276 0.84895560 1.65980277 1.20047949 -0.4650972
[4,] -0.1859943 2.3187715 0.13240555 0.25299232 -0.02374193 0.6060723
[5,] 1.2554234 0.7640444 0.09494619 -0.53153289 0.11199489 0.9605717
[6,] 0.3320614 1.6194709 -0.98040416 -0.31037303 -0.34118294 -0.3479703
PC7 PC8 PC9
[1,] -0.02150042 -0.0187908 -1.0754498
[2,] 0.02444363 -0.1073828 0.3920544
[3,] 0.91611899 -0.5614812 -0.3386311
[4,] 0.55329350 -0.9033987 1.0548838
[5,] 1.42907857 -0.8930173 -0.8764615
[6,] 0.38068132 1.5110579 0.8498801prcomp( )の主成分ヴェクトル
主成分ヴェクトルpr1$rotationも表示させてみよう。
pr1$rotation # 主成分ヴェクトル PC1 PC2 PC3 PC4 PC5 PC6
q0101 0.49802484 -0.05594720 0.04395227 -0.24386147 0.10064407 -0.60884535
q0102 0.55462152 -0.11620220 -0.06292039 -0.13546901 0.04559114 0.33674794
q0103 0.51380170 -0.14901883 -0.37208178 0.02327602 0.09483285 0.19862640
q0200 0.19839061 0.50629241 0.12871992 0.15505025 0.37457317 -0.35694017
q0301 0.15476394 0.45689311 0.30059386 0.17007928 -0.44947710 -0.09155721
q0302 0.25179900 -0.09923688 0.72421757 0.28591804 0.14384478 0.39301714
q0401 -0.14936189 0.45042646 -0.10294634 -0.20288696 0.64220223 0.30941826
q0402 0.15650630 0.37518387 -0.45674402 0.52435071 -0.24006680 0.16487432
q0403 -0.08280943 -0.37760050 -0.07095590 0.68615925 0.38172915 -0.25155723
PC7 PC8 PC9
q0101 0.01426529 0.5295311033 -0.16128354
q0102 0.18406205 -0.0005812592 0.71219681
q0103 0.09917056 -0.4259953369 -0.57884543
q0200 -0.34882468 -0.4825835005 0.20352623
q0301 0.64768373 -0.1063191634 -0.08247915
q0302 -0.24408349 0.1866299324 -0.22571622
q0401 0.33870309 0.2873705693 -0.13259326
q0402 -0.30042330 0.4194222320 0.01411737
q0403 0.38776191 0.0246223850 0.12143918データとヴェクトルから得点合成
第1主成分得点はpr1$x[,1]で取り出せるが，元のデータ行列を標準化したものscale(vars)と第1主成分ヴェクトルpr1$rotation[,1]から主成分得点を手動で合成することもできる。以下のようにすると先頭10ケースについてその二つを並べて比較し，完全に一致していることが確認できる。%*%は行列の積の演算子である。
# 主成分得点を再現
score1 <- scale(vars) %*% pr1$rotation[,1]
cbind(pr1$x[,1], score1)[1:10, ] [,1] [,2]
[1,] -2.1456084 -2.1456084
[2,] 1.4349870 1.4349870
[3,] -2.1404049 -2.1404049
[4,] -0.1859943 -0.1859943
[5,] 1.2554234 1.2554234
[6,] 0.3320614 0.3320614
[7,] 0.4069525 0.4069525
[8,] -0.9670763 -0.9670763
[9,] -2.3684490 -2.3684490
[10,] 2.0735540 2.0735540主成分負荷量
積率相関としての計算
主成分負荷量（主成分得点と元の変数との積率相関係数）は以下の通り。行が元の9項目の観測変数である。
cor(scale(vars), pr1$x[, 1:2]) PC1 PC2
q0101 0.6713545 -0.07140497
q0102 0.7476487 -0.14830796
q0103 0.6926222 -0.19019156
q0200 0.2674373 0.64617703
q0301 0.2086271 0.58312909
q0302 0.3394336 -0.12665525
q0401 -0.2013449 0.57487576
q0402 0.2109758 0.47884423
q0403 -0.1116299 -0.48192855主成分負荷量の導出
主成分負荷量は，主成分分析の結果の情報から以下の様にも求められる。
t(pr1$rotation) * pr1$sdev q0101 q0102 q0103 q0200 q0301 q0302
PC1 0.67135446 0.7476487053 0.69262221 0.2674373 0.20862706 0.3394336
PC2 -0.07140497 -0.1483079572 -0.19019156 0.6461770 0.58312909 -0.1266552
PC3 0.04526809 -0.0648040732 -0.38322098 0.1325735 0.30959289 0.7458988
PC4 -0.23993480 -0.1332876775 0.02290123 0.1525536 0.16734065 0.2813142
PC5 0.09423635 0.0426884849 0.08879512 0.3507252 -0.42086021 0.1346866
PC6 -0.54643618 0.3022298795 0.17826637 -0.3203523 -0.08217222 0.3527313
PC7 0.01203146 0.1552393548 0.08364122 -0.2942014 0.54626145 -0.2058619
PC8 0.40323673 -0.0004426276 -0.32439448 -0.3674862 -0.08096180 0.1421183
PC9 -0.11978737 0.5289577776 -0.42991598 0.1511616 -0.06125833 -0.1676424
q0401 q0402 q0403
PC1 -0.20134492 0.21097582 -0.11162994
PC2 0.57487576 0.47884423 -0.48192855
PC3 -0.10602830 -0.47041780 -0.07308014
PC4 -0.19962007 0.51590758 0.67511067
PC5 0.60131510 -0.22478245 0.35742557
PC6 0.27770161 0.14797402 -0.22577157
PC7 0.28566479 -0.25337932 0.32704138
PC8 0.21883203 0.31938907 0.01874989
PC9 -0.09847873 0.01048515 0.09019445必要な部分の抽出
t( )は転置行列を与える関数である。
t((t(pr1$rotation) * pr1$sdev)[1:2,]) PC1 PC2
q0101 0.6713545 -0.07140497
q0102 0.7476487 -0.14830796
q0103 0.6926222 -0.19019156
q0200 0.2674373 0.64617703
q0301 0.2086271 0.58312909
q0302 0.3394336 -0.12665525
q0401 -0.2013449 0.57487576
q0402 0.2109758 0.47884423
q0403 -0.1116299 -0.48192855biplot
各ケースの位置関係と，元々の各変数の位置関係を同時に図示するのが，biplot( )関数である。
最も単純な図
第1主成分と第2主成分の関係を図示するには，biplot(pr1)だけでよい。
biplot(pr1)
見易くした図
いくつかのグラフオプションを指定して描いてみる。各変数は主成分ヴェクトルの定数倍（定数はケース数の平方根×標準偏差）でプロットされ（上目盛りと右目盛り），各ケースは主成分得点の定数倍（ケース数の平方根×標準偏差の逆数）でプロットされている（下目盛りと左目盛り）。ケース数が多いとケースのプロットはあまり役に立たないかも知れない。また，目盛りの値自体は気にせず，変数相互の位置関係を読み取るとよい。
# 作図の準備と作図
dimnames(pr1$rotation)[[1]] <- c("ネット通販", "スマホ", "SNS",
"蔵書", "図書館", "マンガ", "新聞", "奉仕", "喫煙")
biplot(pr1, cex=c(0.5, 0.8), col=c("#00000040", "red"),
font=2, family="serif", xlab="第1主成分", ylab="第2主成分")
2 探索的因子分析(Exploratory Factor Analysis)
（探索的）因子分析(EFA; Exploratory Factor Analysis)のモデルを図示すると以下の様になる。潜在変数と観測変数の間のパス＝矢印（規定関係）が主成分分析とは逆であるのに加え，観測変数にはそれ以外に誤差（独自因子）からのパスも引かれている。
因子分析においてはその他に独自因子（＝誤差; error）を e1 などと表記し，○も□もつけない事とする。
因みにこれらのグラフはGraphviz; Graph Visualization Softwareを用いて作成した。コマンドはWatanabe Shin, 「Graphvizとdot言語でグラフを描く方法のまとめ」(2020年10月20日閲覧)を見ながら書くとこの程度のものは何とか書ける様になるだろう。結果の保存の仕方はやや分かりにくいので注意する事。使用したコマンド（dot言語）のファイルも提供しておく： 主成分分析グラフ，（探索的）因子分析グラフ
2-1 潜在変数(latent variable)から観測変数(observed variable)への影響
例として，6変数から3因子を抽出するモデル
\[x_{1}=a_{11}f_{1}+a_{12}f_{2}+a_{13}f_{3}+\epsilon_{1}\] \[x_{2}=a_{21}f_{1}+a_{22}f_{2}+a_{23}f_{3}+\epsilon_{2}\] \[x_{3}=a_{31}f_{1}+a_{32}f_{2}+a_{33}f_{3}+\epsilon_{3}\] \[x_{4}=a_{41}f_{1}+a_{42}f_{2}+a_{43}f_{3}+\epsilon_{4}\] \[x_{5}=a_{51}f_{1}+a_{52}f_{2}+a_{53}f_{3}+\epsilon_{5}\] \[x_{6}=a_{61}f_{1}+a_{62}f_{2}+a_{63}f_{3}+\epsilon_{6}\]
2-2 因子の選出
変数群varsを固有値分解した結果の固有値のスクリープロットを，色々と装飾しながら描いてみる。赤が各因子の固有値の大きさ，青は累積割合である。
# 変数群varsの相関係数行列を固有値分解した固有値（分散）の値
eigen(cor(vars))$values[1] 1.8171959 1.6289217 1.0607713 0.9680552 0.8767194 0.8054988 0.7113366
[8] 0.5798788 0.5516225# その固有値の累積割合を計算してヴェクトルに格納
cum.prop <- cumsum(eigen(cor(vars))$values)/sum(eigen(cor(vars))$values)
cum.prop[1] 0.2019107 0.3829019 0.5007654 0.6083271 0.7057404 0.7952402 0.8742776
[8] 0.9387086 1.0000000# 各固有値のグラフ
plot(eigen(cor(vars))$values, type="b", family="serif", col="red",
xlab="因子番号", ylab="固有値", main="9変数の固有値分解の結果", yaxt="n")
axis(side=2, col="red", col.axis="red", family="serif")
abline(h=1, col="red") # 固有値1以上の基準線
par(new=T) # グラフの重ね描き
# 固有値の累積割合のグラフ
plot(cum.prop, type="b", col="blue", axes=F, xlab="", ylab="",
pch=4, lty=2)
axis(side=4, col="blue", col.axis="blue", family="serif")
abline(h=.6, col="blue", lty=2) # 累積割合60%の基準線
パッケイジ“psych”のfa.parallel( )関数による平行分析
適切な因子数を決定する為に，やや専門的な「平行分析」を行ってみる。パッケイジpsychが必要になるので，インストールしていなければ，install.packages(“psych”, repos=“http://cran.ism.ac.jp/”) などとしてインストールする必要がある。既にインストールしてある場合には，library(psych)としてパッケイジを有効化するか，或いは以下の様にダブルコロンを用いてパッケイジの中の特定の関数だけを呼び出す。この分析は乱数を発生させて，乱数の相関係数行列を固有値分解した結果と比較する為（清水裕士「因子分析における因子数選択のための基準」2012年5月23日），場合によっては何度か繰り返すうちに適切な因子数が1だけ変化する事が有り得る。この分析例では，2もしくは3の因子数が適切となった。 分析のたびにグラフを出力する。PCは主成分分析の場合，FAは因子分析の場合である。赤い点線で表されているシミュレイションの結果よりも固有値が大きければ採用すべきとなる。
# install.packages("psych", repos="http://cran.ism.ac.jp/")
library(psych)
psych::fa.parallel(vars, SMC=TRUE)
Parallel analysis suggests that the number of factors = 3 and the number of components = 2 2-3 因子負荷量(factor loadings)と寄与率(contribution)
Rで因子分析を行う関数も複数あるが，まずは標準の因子分析の関数factanal( )で分析する。出力結果を見やすくする為に，先に変数群のヴェクトルvarsに列名（変数ラベル）を付けておく。
colnames(vars) <- c("NetShop", "Devices", "SNS", "Books", "Library", "Comic", "NewsP", "Volunteer", "Smoking")因子数の指定
3因子解
独自因子（言い換えれば誤差）Uniquenesses，因子負荷量Loadings（線型方程式の係数），因子寄与（寄与度）SS loadings（因子負荷量の二乗和），分散割合Proportion Var，累積分散割合Cumulative Var，因子数についての検定結果を確認しよう。
factanal(vars, factors=3)
Call:
factanal(x = vars, factors = 3)
Uniquenesses:
NetShop Devices SNS Books Library Comic NewsP Volunteer
0.723 0.621 0.376 0.648 0.799 0.866 0.825 0.802
Smoking
0.912
Loadings:
Factor1 Factor2 Factor3
NetShop 0.364 0.367
Devices 0.516 0.335
SNS 0.787
Books 0.588
Library 0.438
Comic 0.359
NewsP -0.120 0.338 -0.216
Volunteer 0.188 0.345 -0.210
Smoking -0.295
Factor1 Factor2 Factor3
SS loadings 1.077 0.871 0.480
Proportion Var 0.120 0.097 0.053
Cumulative Var 0.120 0.216 0.270
Test of the hypothesis that 3 factors are sufficient.
The chi square statistic is 22.35 on 12 degrees of freedom.
The p-value is 0.0338 4因子解
因子数の検定結果は，4因子の場合も紹介しよう。「4因子で十分である」とのゼロ仮説は5%水準で棄却されなかった。
factanal(vars, factors=4)
Call:
factanal(x = vars, factors = 4)
Uniquenesses:
NetShop Devices SNS Books Library Comic NewsP Volunteer
0.722 0.579 0.590 0.663 0.819 0.005 0.836 0.005
Smoking
0.889
Loadings:
Factor1 Factor2 Factor3 Factor4
NetShop 0.508 0.119
Devices 0.640 0.103
SNS 0.622 0.115
Books 0.101 0.562
Library 0.103 0.404
Comic 0.121 0.989
NewsP -0.174 0.356
Volunteer 0.974 0.206
Smoking -0.328
Factor1 Factor2 Factor3 Factor4
SS loadings 1.111 1.013 0.989 0.779
Proportion Var 0.123 0.113 0.110 0.087
Cumulative Var 0.123 0.236 0.346 0.432
Test of the hypothesis that 4 factors are sufficient.
The chi square statistic is 12.49 on 6 degrees of freedom.
The p-value is 0.052 factanal( )の斜交解の因子間相関の誤り
標準のfactanal( )関数では，斜交promax解の因子間相関係数行列が正しく出力されない問題について，詳しくは説明しないが，一致しなければいけない筈の結果が食い違ってしまう事だけ示しておこう。以下のサイトも参考になる。
緒賀郷志「Rによる因子分析関連メモ&主成分分析とその回転メモ」(2020年10月21日閲覧)
裏 RjpWiki「因子間相関係数行列がへん?」(2020年10月21日閲覧)
factanal(vars, factors=3, rotation="promax") # 標準でのpromax解
Call:
factanal(x = vars, factors = 3, rotation = "promax")
Uniquenesses:
NetShop Devices SNS Books Library Comic NewsP Volunteer
0.723 0.621 0.376 0.648 0.799 0.866 0.825 0.802
Smoking
0.912
Loadings:
Factor1 Factor2 Factor3
NetShop 0.138 0.166 0.427
Devices 0.339 0.398
SNS 0.769
Books 0.603
Library 0.458 -0.111 0.126
Comic -0.111 0.398
NewsP 0.318 -0.227
Volunteer 0.318 0.243 -0.201
Smoking -0.303
Factor1 Factor2 Factor3
SS loadings 0.899 0.825 0.621
Proportion Var 0.100 0.092 0.069
Cumulative Var 0.100 0.192 0.261
Factor Correlations:
Factor1 Factor2 Factor3
Factor1 1.000 0.139 -0.378
Factor2 0.139 1.000 0.105
Factor3 -0.378 0.105 1.000
Test of the hypothesis that 3 factors are sufficient.
The chi square statistic is 22.35 on 12 degrees of freedom.
The p-value is 0.0338 # 標準で（回転させない）初期解を求めてpromax( )関数に与える
fa_n <- factanal(vars, factors=3, rotation="none"); promax(fa_n$loadings)$loadings
Loadings:
Factor1 Factor2 Factor3
NetShop 0.166 0.138 0.427
Devices 0.339 0.398
SNS 0.769
Books 0.603
Library -0.111 0.458 0.126
Comic -0.111 0.398
NewsP 0.318 -0.227
Volunteer 0.243 0.318 -0.201
Smoking -0.303
Factor1 Factor2 Factor3
SS loadings 0.825 0.899 0.621
Proportion Var 0.092 0.100 0.069
Cumulative Var 0.092 0.192 0.261
$rotmat
[,1] [,2] [,3]
[1,] 0.8072816 0.02995755 0.35516522
[2,] -0.1317796 1.01376498 0.02801236
[3,] -0.7370832 0.14863144 1.03697397# 因子間相関行列を求める
solve(t(promax(fa_n$loadings)$rotmat) %*% promax(fa_n$loadings)$rotmat) [,1] [,2] [,3]
[1,] 1.0000000 0.1388985 0.3780262
[2,] 0.1388985 1.0000000 -0.1050634
[3,] 0.3780262 -0.1050634 1.0000000psychパッケイジのfa( )関数では，結果がやや異なるものになる。
Package’psych’ October 5, 2020のp. 128辺りに関連の記述がある。
# psychパッケイジのfa( )関数を使用するので，インストールした事が無ければまずインストールする
# install.packages("psych", repos="https://cran.ism.ac.jp")
# library(psych)
fa_p2 <- psych::fa(vars, nfactors=3, rotate="promax", scores="tenBerge", fm="ml")Loading required namespace: GPArotationfa_p2; fa_p2$PhiFactor Analysis using method = ml
Call: psych::fa(r = vars, nfactors = 3, rotate = "promax", scores = "tenBerge",
fm = "ml")
Standardized loadings (pattern matrix) based upon correlation matrix
ML1 ML2 ML3 h2 u2 com
NetShop 0.30 0.11 0.34 0.277 0.72 2.2
Devices 0.47 0.01 0.27 0.379 0.62 1.6
SNS 0.83 -0.15 -0.11 0.624 0.38 1.1
Books 0.01 0.60 0.08 0.352 0.65 1.0
Library -0.06 0.46 0.13 0.201 0.80 1.2
Comic -0.01 0.06 0.37 0.134 0.87 1.1
NewsP -0.10 0.33 -0.19 0.175 0.82 1.8
Volunteer 0.22 0.30 -0.23 0.198 0.80 2.8
Smoking 0.04 -0.30 -0.05 0.088 0.91 1.1
ML1 ML2 ML3
SS loadings 1.08 0.87 0.48
Proportion Var 0.12 0.10 0.05
Cumulative Var 0.12 0.22 0.27
Proportion Explained 0.44 0.36 0.20
Cumulative Proportion 0.44 0.80 1.00
With factor correlations of
ML1 ML2 ML3
ML1 1.00 0.14 0.32
ML2 0.14 1.00 -0.14
ML3 0.32 -0.14 1.00
Mean item complexity = 1.5
Test of the hypothesis that 3 factors are sufficient.
The degrees of freedom for the null model are 36 and the objective function was 0.72 with Chi Square of 263.81
The degrees of freedom for the model are 12 and the objective function was 0.06
The root mean square of the residuals (RMSR) is 0.03
The df corrected root mean square of the residuals is 0.05
The harmonic number of observations is 373 with the empirical chi square 26.69 with prob < 0.0085
The total number of observations was 373 with Likelihood Chi Square = 22.35 with prob < 0.034
Tucker Lewis Index of factoring reliability = 0.863
RMSEA index = 0.048 and the 90 % confidence intervals are 0.013 0.079
BIC = -48.71
Fit based upon off diagonal values = 0.95
Measures of factor score adequacy
ML1 ML2 ML3
Correlation of (regression) scores with factors 0.84 0.74 0.64
Multiple R square of scores with factors 0.71 0.54 0.41
Minimum correlation of possible factor scores 0.41 0.08 -0.17 ML1 ML2 ML3
ML1 1.0000000 0.1424179 0.3238314
ML2 0.1424179 1.0000000 -0.1395861
ML3 0.3238314 -0.1395861 1.00000002-4 因子の解釈
プロマックス解の結果を出力する。
pro1 <- promax(fa_n$loadings) # 初期解のプロマックス回転の結果をpro1に保存
print(pro1$loadings, cutoff=0) # pro1の中のloadingsの情報を，省略せず表示
Loadings:
Factor1 Factor2 Factor3
NetShop 0.166 0.138 0.427
Devices 0.339 0.053 0.398
SNS 0.769 -0.078 0.065
Books -0.040 0.603 0.091
Library -0.111 0.458 0.126
Comic -0.111 0.064 0.398
NewsP -0.055 0.318 -0.227
Volunteer 0.243 0.318 -0.201
Smoking 0.061 -0.303 -0.046
Factor1 Factor2 Factor3
SS loadings 0.825 0.899 0.621
Proportion Var 0.092 0.100 0.069
Cumulative Var 0.092 0.192 0.261solve(t(pro1$rotmat) %*% pro1$rotmat) # 正しい因子間相関係数行列の計算 [,1] [,2] [,3]
[1,] 1.0000000 0.1388985 0.3780262
[2,] 0.1388985 1.0000000 -0.1050634
[3,] 0.3780262 -0.1050634 1.0000000グラフ
単純なbiplot( )
グラフを作図するには，以下の様に biplot( ) を使用する事も出来る。
fa_p <- factanal(vars, factors=3, rotation="promax", scores="regression")
biplot(fa_p$scores, fa_p$loadings)
オプション付き
biplot(fa_p$scores, fa_p$loadings,
cex=c(0.5, 0.8), col=c("#00000040", "red"), font=2, family="serif")
手動で作図
plot(pro1$loadings[,c(1,2)], xlim=c(-.15, .85), ylim=c(-.32, .65), xlab="第1因子", ylab="第2因子", pch=16, col="red")
abline(h=0, v=0, lty=2, col="gray") # 原点を通る座標を追加
text(pro1$loadings[,c(1,2)] + .04, colnames(vars), family="serif", col="firebrick")
パッケイジ“psych”のfa( )関数とalpha( )関数
psychパッケイジのfa( )関数と，クロンバックのα係数による内的整合性の検討の為のalpha( )関数の結果をフルサイズで掲載しておく。
fa(vars, nfactors=3, rotate="promax", scores="regression", fm="ml")Factor Analysis using method = ml
Call: fa(r = vars, nfactors = 3, rotate = "promax", scores = "regression",
fm = "ml")
Standardized loadings (pattern matrix) based upon correlation matrix
ML1 ML2 ML3 h2 u2 com
NetShop 0.30 0.11 0.34 0.277 0.72 2.2
Devices 0.47 0.01 0.27 0.379 0.62 1.6
SNS 0.83 -0.15 -0.11 0.624 0.38 1.1
Books 0.01 0.60 0.08 0.352 0.65 1.0
Library -0.06 0.46 0.13 0.201 0.80 1.2
Comic -0.01 0.06 0.37 0.134 0.87 1.1
NewsP -0.10 0.33 -0.19 0.175 0.82 1.8
Volunteer 0.22 0.30 -0.23 0.198 0.80 2.8
Smoking 0.04 -0.30 -0.05 0.088 0.91 1.1
ML1 ML2 ML3
SS loadings 1.08 0.87 0.48
Proportion Var 0.12 0.10 0.05
Cumulative Var 0.12 0.22 0.27
Proportion Explained 0.44 0.36 0.20
Cumulative Proportion 0.44 0.80 1.00
With factor correlations of
ML1 ML2 ML3
ML1 1.00 0.14 0.32
ML2 0.14 1.00 -0.14
ML3 0.32 -0.14 1.00
Mean item complexity = 1.5
Test of the hypothesis that 3 factors are sufficient.
The degrees of freedom for the null model are 36 and the objective function was 0.72 with Chi Square of 263.81
The degrees of freedom for the model are 12 and the objective function was 0.06
The root mean square of the residuals (RMSR) is 0.03
The df corrected root mean square of the residuals is 0.05
The harmonic number of observations is 373 with the empirical chi square 26.69 with prob < 0.0085
The total number of observations was 373 with Likelihood Chi Square = 22.35 with prob < 0.034
Tucker Lewis Index of factoring reliability = 0.863
RMSEA index = 0.048 and the 90 % confidence intervals are 0.013 0.079
BIC = -48.71
Fit based upon off diagonal values = 0.95
Measures of factor score adequacy
ML1 ML2 ML3
Correlation of (regression) scores with factors 0.84 0.74 0.64
Multiple R square of scores with factors 0.71 0.54 0.41
Minimum correlation of possible factor scores 0.41 0.08 -0.17alpha(vars[,1:3])
Reliability analysis
Call: alpha(x = vars[, 1:3])
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.54 0.61 0.51 0.34 1.5 0.032 3.2 1.4 0.32
lower alpha upper 95% confidence boundaries
0.48 0.54 0.6
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
NetShop 0.54 0.58 0.41 0.41 1.40 0.042 NA 0.41
Devices 0.30 0.45 0.29 0.29 0.81 0.043 NA 0.29
SNS 0.41 0.48 0.32 0.32 0.94 0.049 NA 0.32
Item statistics
n raw.r std.r r.cor r.drop mean sd
NetShop 373 0.54 0.72 0.46 0.35 2.2 0.89
Devices 373 0.75 0.77 0.59 0.45 3.9 1.71
SNS 373 0.88 0.76 0.56 0.44 3.4 2.66
Non missing response frequency for each item
1 2 3 4 5 6 7 8 miss
NetShop 0.25 0.44 0.23 0.09 0.00 0.00 0.00 0.0 0
Devices 0.24 0.02 0.02 0.05 0.68 0.00 0.00 0.0 0
SNS 0.24 0.39 0.05 0.02 0.05 0.01 0.04 0.2 0発展1 構造方程式モデリング(Structural Equation Modelling)
発展1-1 確証的因子分析(CFA)と構造方程式モデリング(SEM)
パス・ダイアグラム
パス・ダイアグラム2
発展1-2 構造方程式モデリングとパス解析(path analysis)
パス・ダイアグラム3