肉なし回鍋肉

なんでも。 主にNEW GAME!, 野球, 麻雀

【第2回】打者の「能力」を数値化してみる 【因子分析】

こんにちは。肉なし回鍋肉です。今回は「打者の「能力」を数値化してみる」の2回目です。実際に因子分析ってやつをやります。ちなみに今回は基本統計量、ヒストグラム、相関行列の作り方はやりませんが、やっとくと良い(らしい)です。

全体の目次

今回の目次

因子分析の解説

前回お話ししたように、因子分析とは、多次元のデータを分析する手法(多変量解析)の一つです。因子分析は、観測された変数の相関関係を使い、それらの背後に因子と呼ばれる共通の原因があると仮定して探し出す手法です。最低限これだけ分かってれば良いので、結果だけ知りたい人は「分析2 因子分析」まで飛ばしてもらって構いません。

因子分析の詳細

因子分析は、(教育)心理学とかマーケティングとかでよく使われる手法です。5教科の成績から文系能力と理系能力を定めるみたいな感じです。
用語とかを簡単に解説します。筆者も詳しく理解してるわけではないので多少間違いがあるかもしれません。

f:id:hui_guo:20200730204544p:plain
因子分析のイメージ

これは前回も使用した因子分析の簡単なイメージです。この図の「成績」は、私たちが実際に目に見える形で得られたデータです。このようなデータを「観測変数」と言います。図の中で「能力」としているものは、実際には目に見えない存在です。このように、観測変数に影響を与えている、目に見えない潜在的な要因を「因子」と言います。
因子には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()
f:id:hui_guo:20200731121628p:plain
スクリープロット1

第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です。ということで、

  1. 実行コード
  2. スクリープロット
  3. 出力結果
  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)
f:id:hui_guo:20200731123528p:plain
スクリープロット2
> 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)
f:id:hui_guo:20200731125220p:plain
スクリープロット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)
f:id:hui_guo:20200731125918p:plain
スクリープロット4
> 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)
f:id:hui_guo:20200731130600p:plain
スクリープロット5
> 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

今回の結果は、特にいちゃもんつける点がありません。ということで、これで分析パートは終了です。最後なので後に続くいろんな出力結果も載せておきます。
分かりやすくしたのはこちら↓
f:id:hui_guo:20200731132240p:plain

ここからは、各因子がどのような能力なのか分析するのですが・・・
すでに結構ボリュームあるので今回はここまでにします。

今回のまとめ

  • 因子数決める→計算→回転
  • 共通性や複雑性に気をつけながら、納得いくまで繰り返す
  • 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
この他にも色々あります。興味あったら探してみてください。