【第2回】打者の「能力」を数値化してみる 【因子分析】
こんにちは。肉なし回鍋肉です。今回は「打者の「能力」を数値化してみる」の2回目です。実際に因子分析ってやつをやります。ちなみに今回は基本統計量、ヒストグラム、相関行列の作り方はやりませんが、やっとくと良い(らしい)です。
全体の目次
今回の目次
因子分析の解説
前回お話ししたように、因子分析とは、多次元のデータを分析する手法(多変量解析)の一つです。因子分析は、観測された変数の相関関係を使い、それらの背後に因子と呼ばれる共通の原因があると仮定して探し出す手法です。最低限これだけ分かってれば良いので、結果だけ知りたい人は「分析2 因子分析」まで飛ばしてもらって構いません。
因子分析の詳細
因子分析は、(教育)心理学とかマーケティングとかでよく使われる手法です。5教科の成績から文系能力と理系能力を定めるみたいな感じです。
用語とかを簡単に解説します。筆者も詳しく理解してるわけではないので多少間違いがあるかもしれません。
これは前回も使用した因子分析の簡単なイメージです。この図の「成績」は、私たちが実際に目に見える形で得られたデータです。このようなデータを「観測変数」と言います。図の中で「能力」としているものは、実際には目に見えない存在です。このように、観測変数に影響を与えている、目に見えない潜在的な要因を「因子」と言います。
因子には2種類あります。一つはこの図の「能力」のように、複数の観測変数に影響する因子です。これは「共通因子」と呼ばれます。もう一つは、この図には書かれていないのですが1つの観測変数にのみ影響を与えている因子です。例としては、「打点」におけるランナーの存在みたいな感じです。ランナーがいても(だいたい)打率やホームラン数には影響しないですが、打点はランナーの存在で大きく変わりますよね。イメージはこんな感じです。
各因子が、それぞれの説明変数を同じ程度説明しているわけではなく、そこには強弱があります。その強弱の値を「因子負荷量」って言います。詳しい資料は記事の最後の「参考」に載せます。
因子分析の流れ
ここでは因子分析をどういう流れで行うか説明します。
1.観測変数の選出 : 前回紹介した29変数すべてが使えるわけではありません。分析すると不要な変数が出てきたりします。ということで、まずは変数選びからです。
2.共通因子数の決定 : そもそも因子が何個あるとして分析するかを決めます。少ない共通因子数で多くの観測変数を説明できると天晴れですので、多すぎず少なすぎずな因子数を決めることは大切なんです。適当に決めてはまずいのです。具体的な方法はいくつかありますが、今回はスクリー法によって決定します。スクリー法とは、固有値をプロットした図であるスクリープロットにおいて、固有値の減少がなだらかになる因子の直前の因子までを採用する方法です。固有値とは、各因子が観測変数をどれだけ説明しているかを表す数値です。例えばある因子の固有値が3.7ならば、その因子は観測変数3.7個分を説明していることになります。理論上、因子は最大で観測変数の数だけ存在して、全因子の固有値の和は観測変数の総数になります。(なんでなだらかになる直前の因子までを採用するとかは知りません)
3.因子の抽出 : 色々計算するパートです。計算方法にも色々あるのですが、今回は一般化最小二乗法を採用します。最尤法ってのもポピュラーなんですけど、これはデータの正規性を仮定した手法なんです。今回のデータは正規分布じゃないデータがいっぱいあるので今回は一般化最小二乗法です。
4.因子を回転させる : 突然回転とか言い出して横転。抽出したままの結果は、そのままだと意味不明です。そこで、因子軸を回転させても元の観測変数の関係は不変である性質を利用して、抽出された因子の解釈を簡単にするために因子軸を回転させます。ここの説明は詳しくできるほど理解してないので「結果を分かりやすくする」程度の認識です。
回転方法には大きく分けると直交回転と斜交回転があります。直交回転は、因子同士が無相関であると仮定した回転法で、斜交回転は因子間に相関があると仮定した回転法です。今回は因子間の相関を仮定して、斜交回転の一つであるプロマックス回転を採用します。
これを満足するまで繰り返します。
分析2 因子分析
ようやく、ここから本番です。
1回目
まずは因子数を決定します。
#因子分析1 #分析に使う観測変数を決める statsNames1 <- c("打率", "試合", "打席", "打数", "得点", "安打", "二塁打", "三塁打", "本塁打", "塁打", "打点", "盗塁", "盗塁刺", "犠打", "犠飛", "四球", "故意四球", "死球", "三振", "併殺打", "長打率", "出塁率", "OPS", "BABIP", "BB.", "K.", "HR.", "BB.K", "ISO") zstats1 <- zdatas[statsNames1] #スクリープロットの出力 png("screePlot1.png", width = 500, height = 500) VSS.scree(zstats1) dev.off()
第6因子からなだらかですね。ということで因子数を5にします。
f1 <- fa(zstats1, nfactors = 5, fm = "gls", rotate = "promax", scores = T) print(f1, sort = T, digit = 3)
結果は以下の通り。
> f1 <- fa(zstats1, nfactors = 5, fm = "gls", + rotate = "promax", scores = T) solve.default(S) でエラー: システムは数値的に特異です: 条件数の逆数 = 1.13542e-18 追加情報: 警告メッセージ: 1: cor.smooth(R) で: Matrix was not positive definite, smoothing was done 2: cor.smooth(R) で: Matrix was not positive definite, smoothing was done > print(f1, sort = T, digit = 3) print(f1, sort = T, digit = 3) でエラー: オブジェクト 'f1' がありません
エラー出ちゃった!原因としては、観測変数同士で相関係数が高い組み合わせがあったことです。計算の過程でおかしくなるらしい。ということで、相関係数の絶対値が0.9を超える組み合わせが以下の通り。
相関係数の高い組み合わせ:打席と打数、本塁打とHR%、本塁打とISO、四球とBB%、三振とK%、長打率とOPS、長打率とISO、HR%とISO
これらの組み合わせと解釈性を考慮して、試合、打席、打数、本塁打、四球、故意四球、三振、長打率、OPSの9変数を削除して2回目の因子分析を行います。
2回目
ネタバレ(?)ですが、ここからは最後まで因子数は4です。ということで、
- 実行コード
- スクリープロット
- 出力結果
- コメント
の流れで進めます。
#因子分析2 statsNames2 <- c("打率", "得点", "安打", "二塁打", "三塁打", "塁打", "打点", "盗塁", "盗塁刺", "犠打", "犠飛", "死球", "併殺打", "出塁率", "BABIP", "BB.", "K.", "HR.", "BB.K", "ISO") zstats2 <- zdatas[statsNames2] png("screePlot2.png", width = 500, height = 500) VSS.scree(zstats2) dev.off() f2 <- fa(zstats2, nfactors = 4, fm = "gls", rotate = "promax", scores = T) print(f2, sort = T, digit = 3)
> f2 <- fa(zstats2, nfactors = 4, fm = "gls", + rotate = "promax", scores = T) Loading required namespace: GPArotation 警告メッセージ: fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, で: A loading greater than abs(1) was detected. Examine the loadings carefully. > print(f2, sort = T, digit = 3) Factor Analysis using method = gls Call: fa(r = zstats2, nfactors = 4, rotate = "promax", scores = T, fm = "gls") Warning: A Heywood case was detected. Standardized loadings (pattern matrix) based upon correlation matrix item GLS4 GLS1 GLS2 GLS3 h2 u2 com 安打 3 1.015 -0.111 0.054 -0.133 0.8817 0.1183 1.06 打率 1 0.929 -0.266 -0.031 0.127 0.8308 0.1692 1.21 BABIP 15 0.721 -0.099 0.342 0.036 0.6058 0.3942 1.48 二塁打 4 0.699 0.108 -0.222 -0.162 0.5735 0.4265 1.37 塁打 6 0.632 0.526 -0.135 -0.095 0.9516 0.0484 2.09 得点 2 0.545 0.440 0.368 0.085 0.7192 0.2808 2.78 K. 17 -0.450 1.005 0.211 -0.191 0.7786 0.2214 1.57 ISO 20 0.010 0.826 -0.235 0.043 0.9103 0.0897 1.17 HR. 18 -0.074 0.822 -0.300 0.013 0.9110 0.0890 1.28 打点 7 0.262 0.571 -0.409 -0.016 0.8574 0.1426 2.27 死球 12 0.047 0.210 0.063 0.024 0.0512 0.9488 1.32 盗塁 8 0.194 0.044 0.804 -0.022 0.6476 0.3524 1.12 盗塁刺 9 0.187 -0.037 0.751 -0.042 0.5996 0.4004 1.14 三塁打 5 0.202 0.058 0.709 -0.061 0.5072 0.4928 1.19 併殺打 13 0.084 0.007 -0.555 -0.077 0.3194 0.6806 1.09 犠飛 11 0.203 -0.012 -0.408 -0.014 0.2095 0.7905 1.47 犠打 10 -0.092 -0.302 0.379 -0.151 0.4159 0.5841 2.39 BB. 16 -0.324 0.377 0.052 0.927 0.9359 0.0641 1.60 BB.K 19 0.043 -0.374 -0.122 0.914 0.8315 0.1685 1.37 出塁率 14 0.423 0.080 0.070 0.698 0.9354 0.0646 1.70
本当はもっと色々出てくるのですが、あまりにも長いので省略させていただきます。。簡単に解説すると、GLS1~4は第1因子から第4因子。各項目の数字は因子負荷量です。「h2」は共通性といって、その項目がどれだけ共通因子の影響を受けているかを表しています。小さいと因子とは関係ないってことです。「u2」は独自性です。1-共通性で求められ、どれだけ独自因子の影響が強いかです。「com」は複雑性と言います。これが1に近いほど、一つの因子からのみ影響を受けていると言えます。
h2の小さい観測変数は、各因子から大きな影響を受けていないとみなせます。これらを外すとさらに良くなるので、0.4以下を基準とし、死球、併殺打、犠飛を除いて3回目の因子分析を行います。
3回目
#因子分析3 statsNames3 <- c("打率", "得点", "安打", "二塁打", "三塁打", "塁打", "打点", "盗塁", "盗塁刺", "犠打", "出塁率", "BABIP", "BB.", "K.", "HR.", "BB.K", "ISO") zstats3 <- zdatas[statsNames3] png("screePlot3.png", width = 500, height = 500) VSS.scree(zstats3) dev.off() f3 <- fa(zstats3, nfactors = 4, fm = "gls", rotate = "promax", scores = T) print(f3, sort = T, digit = 3)
> print(f3, sort = T, digit = 3) Factor Analysis using method = gls Call: fa(r = zstats3, nfactors = 4, rotate = "promax", scores = T, fm = "gls") Standardized loadings (pattern matrix) based upon correlation matrix item GLS4 GLS1 GLS3 GLS2 h2 u2 com 安打 3 0.968 -0.040 -0.134 0.119 0.869 0.1312 1.07 打率 1 0.951 -0.245 0.087 -0.071 0.864 0.1362 1.16 BABIP 12 0.724 -0.160 0.010 0.262 0.608 0.3916 1.37 二塁打 4 0.704 0.200 -0.174 -0.149 0.563 0.4371 1.39 塁打 6 0.619 0.616 -0.083 -0.009 0.951 0.0489 2.04 K. 14 -0.444 0.954 -0.156 0.268 0.743 0.2575 1.66 HR. 15 -0.030 0.880 0.023 -0.197 0.919 0.0806 1.10 ISO 17 0.044 0.879 0.055 -0.131 0.918 0.0817 1.06 打点 7 0.270 0.688 -0.002 -0.244 0.822 0.1784 1.58 犠打 10 -0.125 -0.376 -0.145 0.291 0.402 0.5979 2.48 BB. 13 -0.352 0.360 0.952 0.118 0.926 0.0742 1.61 BB.K 16 0.014 -0.353 0.920 -0.110 0.835 0.1650 1.32 出塁率 11 0.415 0.063 0.692 0.072 0.921 0.0788 1.68 盗塁 8 0.061 -0.049 0.035 0.805 0.685 0.3155 1.02 盗塁刺 9 0.064 -0.125 0.009 0.742 0.632 0.3679 1.07 三塁打 5 0.085 -0.016 -0.014 0.709 0.528 0.4720 1.03 得点 2 0.457 0.445 0.129 0.458 0.730 0.2697 3.15
犠打が怪しいですけど、いずれの観測変数も共通性が0.4を超えましたね。しかし、これだとまだ各因子がどういう能力なのか分かりにくいですね。共通性の高い観測変数を除くことで分かりやすくなりますので、2以上を基準に犠打、塁打、得点を取り除いて再分析します。
4回目
#因子分析4 statsNames4 <- c("打率", "安打", "二塁打", "三塁打", "打点", "盗塁", "盗塁刺", "出塁率", "BABIP", "BB.", "K.", "HR.", "BB.K", "ISO") zstats4 <- zdatas[statsNames4] png("screePlot4.png", width = 500, height = 500) VSS.scree(zstats4) dev.off() f4 <- fa(zstats4, nfactors = 4, fm = "gls", rotate = "promax", scores = T) print(f4, sort = T, digit = 3)
> print(f4, sort = T, digit = 3) Factor Analysis using method = gls Call: fa(r = zstats4, nfactors = 4, rotate = "promax", scores = T, fm = "gls") Warning: A Heywood case was detected. Standardized loadings (pattern matrix) based upon correlation matrix item GLS2 GLS1 GLS3 GLS4 h2 u2 com 打率 1 0.927 -0.172 0.089 0.020 0.890 0.1097 1.09 安打 2 0.912 -0.072 -0.097 0.104 0.782 0.2183 1.06 BABIP 9 0.723 0.019 -0.002 0.427 0.720 0.2801 1.62 二塁打 3 0.718 0.154 -0.149 -0.167 0.561 0.4386 1.30 K. 11 -0.326 1.014 -0.139 0.323 0.841 0.1593 1.47 ISO 14 0.142 0.814 0.097 -0.165 0.914 0.0860 1.18 HR. 12 0.062 0.798 0.065 -0.241 0.907 0.0929 1.21 打点 5 0.330 0.594 0.041 -0.288 0.785 0.2146 2.08 BB. 10 -0.278 0.359 0.946 0.119 0.924 0.0757 1.51 BB.K 13 -0.006 -0.394 0.904 -0.143 0.863 0.1371 1.42 出塁率 8 0.447 0.119 0.689 0.134 0.946 0.0543 1.88 盗塁 6 0.057 -0.002 0.031 0.787 0.624 0.3764 1.01 盗塁刺 7 0.048 -0.072 0.006 0.743 0.605 0.3946 1.03 三塁打 4 0.083 0.051 -0.009 0.732 0.520 0.4800 1.04
打点の複雑性が2を超えているので、取り除きます。
ちなみに、K%の第1因子の因子負荷量、1を超えていることに気がつきましたか?本来は1を超えてはいけないのですが、Heywood caseって問題が発生しているためこのような結果になってます(結果の最初の方に警告があります)。因子負荷量が1を超えたからって、K%を取り除けばいいわけではないので注意です。むしろ1に近いなら重要な観測変数と言えます。
5回目
#因子分析5 statsNames5 <- c("打率", "安打", "二塁打", "三塁打", "盗塁", "盗塁刺", "出塁率", "BABIP", "BB.", "K.", "HR.", "BB.K", "ISO") zstats5 <- zdatas[statsNames5] png("screePlot5.png", width = 500, height = 500) VSS.scree(zstats5) dev.off() f5 <- fa(zstats5, nfactors = 4, fm = "gls", rotate = "promax", scores = T) print(f5, sort = T, digit = 3)
> print(f5, sort = T, digit = 3) Factor Analysis using method = gls Call: fa(r = zstats5, nfactors = 4, rotate = "promax", scores = T, fm = "gls") Standardized loadings (pattern matrix) based upon correlation matrix item GLS1 GLS4 GLS3 GLS2 h2 u2 com 打率 1 0.926 -0.172 0.083 -0.013 0.899 0.1008 1.09 安打 2 0.888 -0.092 -0.086 0.084 0.761 0.2391 1.06 BABIP 8 0.740 0.041 -0.016 0.384 0.734 0.2662 1.51 二塁打 3 0.729 0.133 -0.138 -0.207 0.552 0.4477 1.31 K. 10 -0.249 0.995 -0.129 0.259 0.858 0.1424 1.30 ISO 13 0.202 0.762 0.118 -0.227 0.896 0.1040 1.38 HR. 11 0.121 0.746 0.086 -0.299 0.885 0.1148 1.40 BB. 9 -0.242 0.328 0.948 0.104 0.924 0.0764 1.41 BB.K 12 -0.026 -0.407 0.898 -0.114 0.869 0.1308 1.43 出塁率 7 0.471 0.101 0.686 0.101 0.951 0.0491 1.88 盗塁 5 0.030 0.011 0.033 0.799 0.634 0.3661 1.01 盗塁刺 6 0.018 -0.056 0.006 0.757 0.611 0.3892 1.01 三塁打 4 0.057 0.059 -0.003 0.743 0.533 0.4671 1.02 GLS1 GLS4 GLS3 GLS2 SS loadings 3.132 2.480 2.258 2.236 Proportion Var 0.241 0.191 0.174 0.172 Cumulative Var 0.241 0.432 0.605 0.777 Proportion Explained 0.310 0.245 0.223 0.221 Cumulative Proportion 0.310 0.555 0.779 1.000 With factor correlations of GLS1 GLS4 GLS3 GLS2 GLS1 1.000 0.127 0.344 0.086 GLS4 0.127 1.000 0.117 -0.391 GLS3 0.344 0.117 1.000 -0.087 GLS2 0.086 -0.391 -0.087 1.000 Mean item complexity = 1.3 Test of the hypothesis that 4 factors are sufficient. The degrees of freedom for the null model are 78 and the objective function was 20.011 with Chi Square of 15004.7 The degrees of freedom for the model are 32 and the objective function was 9.761 The root mean square of the residuals (RMSR) is 0.066 The df corrected root mean square of the residuals is 0.104 The harmonic number of observations is 756 with the empirical chi square 518.525 with prob < 3.31e-89 The total number of observations was 756 with Likelihood Chi Square = 7293.031 with prob < 0 Tucker Lewis Index of factoring reliability = -0.19 RMSEA index = 0.5511 and the 90 % confidence intervals are 0.5377 0.5588 BIC = 7080.934 Fit based upon off diagonal values = 0.966 Measures of factor score adequacy GLS1 GLS4 Correlation of (regression) scores with factors 0.996 0.999 Multiple R square of scores with factors 0.992 0.997 Minimum correlation of possible factor scores 0.983 0.995 GLS3 GLS2 Correlation of (regression) scores with factors 0.989 0.967 Multiple R square of scores with factors 0.978 0.934 Minimum correlation of possible factor scores 0.957 0.869
今回の結果は、特にいちゃもんつける点がありません。ということで、これで分析パートは終了です。最後なので後に続くいろんな出力結果も載せておきます。
分かりやすくしたのはこちら↓
ここからは、各因子がどのような能力なのか分析するのですが・・・
すでに結構ボリュームあるので今回はここまでにします。
今回のまとめ
- 因子数決める→計算→回転
- 共通性や複雑性に気をつけながら、納得いくまで繰り返す
- Heywood caseに注意
今回はここまでです。次回は因子の命名と下位尺度の構成、因子得点の予測をします。
参考
当時参考にしてたページを貼っておきます。高校の図書室には参考になる本が無かったので(そりゃあそう)、ネットでいろんな大学のpdfとか読んで0から勉強しました。
授業関係
Rの練習帳
この二つは特に重宝してました。気が狂うほど読みました。
http://www.team1mile.com/asarin/hus/ip-hus/materials/factor.pdf
http://www.u.tsukuba.ac.jp/~hirai.akiyo.ft/forstudents/eigokyouikugaku7.files/M4_9(N.K.).pdf
http://cogpsy.educ.kyoto-u.ac.jp/personal/Kusumi/datasem13/masuda.pdf
https://koumurayama.com/koujapanese/factor.htm
http://psy.isc.chubu.ac.jp/~oshiolab/teaching_folder/datakaiseki_folder/08_folder/da08_01.html
この他にも色々あります。興味あったら探してみてください。