ピープルアナリティクス・ケーススタディ：決定木分析で離職リスク分析に挑戦！
決定木分析を利用するメリット
前回までロジスティック回帰で離職リスク分析を行いましたが、今回は同じデータを決定木分析を使って実施したいと思います。
決定木分析は、ツリー構造の条件分岐でデータを説明するモデルです。機械学習の代表的な手法であり、モデル内容が理解しやすいため、実務でも多用されているということです。
また、決定木分析の場合、実行前に独立変数の分布の偏りなどの調整をする必要がないといった点もメリットだそうです。この退職者分析はRを用いて実施したケーススタディを基にPythonで再実施しているのですが、ロジスティック回帰を行う前に年齢分布をより正規分布に近づけるために指数化を行っています。
こういうことを決定木分析ではしなくてもいいということです。
ちなみに、何度か紹介していますが、元のケーススタディはこちらから見ることができます。日本語でこのようなケースがWebで紹介されるとより知識は深まるのですが。。
データを分割するルールをつ決定木分析を実行！
Pythonで決定木を実行するには、scikit-learnにあるtree モジュールのDecisioTreeClassifierクラスを利用します。決定木の深さはこちらで指定しないといけませんが、今回は2(max_depth=2) で指定しています。
from sklearn.tree import DecisionTreeClassifier
tree=DecisionTreeClassifier(max_depth=2)
tree.fit(x_train,y_train)
次は可視化をmatplotlibとscikit-learnにあるtree モジュールのplot_treeで行います。注意点としては①max_depth=の数は、上記の決定木より少ないとすべてが表示されない（例えば上記で3と指定しなおして、そのまま可視化のmax_depthを2のままにしておくと、3番目の要素が表示されません）、②proportion＝TrueをFalseにすると比率でなく観測値の絶対値が表示されます。
import matplotlib.pyplot as plt
from sklearn.tree import plot_tree
plt.figure(figsize=(20,20))
plot_tree(
tree,
max_depth=2,
proportion=True,
filled=True,
rotate=False,
precision=2,
fontsize=12,
feature_names=x_train.columns
)
あと、決定木の中で独立変数名をどうやったら表示されるのか最初はわからなかったのですが、"feature_names= x_train.columns)で表示されました。
結果はこんな感じです。
この図の解釈ですが、上から見ていきましょう。最初はSales（営業）かどうかで別れます。Salesが0.5以下の場合（この場合Salesだと1のフラグがたつのでSales以外の職種ということですね）は左側（オレンジ）、Salesの場合は右側（青）の場合分けになります。
2行目にあるEntropy（エントロピー）というのは不純度の指標です。決定木で分岐をする基準は「不純度を減らしていく」という考え方に基づいているからです。
サンプル数は最初ですので100％になっています。そして全体としては0である（＝辞めずに在籍している）人が62％、１である（＝退職した人）が38％ということになります。
次の階層を見ていきましょう。
営業以外の場合の次の分岐はDfrommedian（給与の中央値）からの差が‐19.56より小さければ右、大きければ左に行きます。給与がやめる原因ということですね。それに対し営業では男性か女性かが次の分岐点で、女性（値は0）であれば左、男性（値は1）であれば右に行きます。
この2の段階で見ると、離職率は、営業以外が30％（Value値の右側）であるのに対し、営業だと56％と営業が高いことがわかります。
3回目の階層です。
営業以外の人は中央値からの給与が‐19.56より小さい（すなわち‐19.56よりさらに低い）場合の離職率は42％とそうでない人（18%）に比べて格段に高くなります。
一方営業職は女性だと66％が離職するという衝撃的な事実です。営業はただでさえ離職率が高いのに女性は男性（45％）よりはるかに高く何らかの対策が必要と考えられます。
決定木分析のメリットと注意点
ロジスティック回帰の結果を見ても、Salesが有意でプラスになっていること、一方性別や給与はマイナスになっていることを考えると同様の傾向が出ています。ただし、営業での離職要因が性別が強く、営業以外だと給与だというような分類はロジスティック回帰ではできないですね。
なにより、決定木での分析結果はわかりやすいというメリットがあると思います。わかりやすいというのは説明がしやすいということにつながるので、ロジスティック回帰の数値を見せるよりこの決定木の結果を見せた方が、マネジメント層の説得がしやすいのではないでしょうか？
実際、ロジスティック回帰は係数の数値が影響度を直接表しているわけではない（ロジスティック回帰の傾きが一定でないため）ため、どの要素が退職に関わっていますよと説明できてもそれ以上の説明が難しいということが言えます。
ちなみに、このモデルの精度をテストデータを使って前回と同様①分類精度と②予測確率の正確さ（ROC曲線、AUC)を出してみましょう。（方法は前回ロジスティック回帰で行った方法と全く同じです）
全体の正解率は68%、退職に限った適合率は67％とそんなに悪くはないですね。
ただし、予測確率の正確さ（ROC曲線、AUC)はロジスティック回帰に比べてかなり悪いです。これは決定木のモデルでの分類要因が少ないことを考えると当然の帰結なのかもしれませんが。。説明のしやすさと予測モデルの正確性の両立はなかなか難しいのかもしれません。（ちなみに決定木の階層を2から4にしたところ、AUCは0.63まではアップしました）。
ピープルアナリティクス・ケーススタディ：ロジスティック回帰分析で離職リスク分析に挑戦！②
検証データでモデルの確からしさを計測する方法
前回のブログでロジスティック回帰を用いて離職リスク分析を行いました。
ただ、実際にこの回帰分析がどれだけ正確に離職者を導き出してくれるのでしょうか？
Pythonでは作成したモデルに対する検証を行うことができます！
今回は前回のブログでご紹介した、Scikit-learnのLogistic Regression関数で作成したモデルの評価を1)カテゴリーの分類精度、2）予測確率の正確さという2点の観点で行いたいと思います。
モデル評価までの下準備はこんな感じです。
1.データを学習データセットとテストデータセットに分類＝＞Scikit-learnのmodel_selectionモジュールのtrain_test_split関数を利用します。式はこんな感じですね。（注意：xは独立変数(Numpy配列もしくはDataFrame), yは従属変数を表すNumpy配列が入っています
from sklearn.model_selection import train_test_split x_train, x_test, y_train, y_test=train_test_split(x, y, train_size=0.67)
2. 学習データセットを用いて分類モデル（ロジスティック回帰分s系）を構築=＞これは前回ご紹介したScikit-learnでロジスティック回帰分析を行うと同じです。
from sklearn.linear_model import LogisticRegression
log_model=LogisticRegression()
result= log_model.fit(x_train, y_train)
3. 予測＝＞テストデータセットで指定した独立変数(x_test)から従属変数を予測従属変数の予測の出し方には、そのものずばりを予測するpredict関数と確率を予測するpredict_proba関数があります。前者は従属変数と同様に0（このケースだと会社にとどまる）か1（退職する）で予測データが算出されるのに対し、後者は0と1それぞれの確率を算出します。ですので、Y_probaの方は1の方の確率を出すように[:,1]をつけることが必要になります。
Y_pred =log_model.predict(x_test)
Y_proba=log_model.predict_proba(x_test)[:,1]
こうした分析を行うにあたってぶち当たるといろいろブログなどをググるほかに、「Python 3 エンジニア認定データ分析試験」の教科書として指定されている「Pythonによる新しいデータ分析の教科書」をいつも参考にしています。
今回の評価でもお世話になったのでご紹介させていただきますね。
カテゴリーの分類精度を抽出！
こちらは上記の予測で出したY_predを利用します。Scikit-learnのmetricsモジュールのclassification_reportという関数を使うと一気に出ます。
from sklearn.metrics import classification_report
print(classification_report(y_test, Y_pred))
結果はこんな感じです。
正解率（正例（このケースでは退職した：記号では1）、負例（このケースでは退職していない：記号では0）どちらにおいても予測が正しかった）は69％、正例（退職者のみに限った正答率）は64％ということですね。
辞める人の3分の2くらいが推測できるというのは、全体の離職率が38％くらいだと考えるとそこそこ予測できているということなのでしょうか。
予測確率の正確さ（ROC曲線、AUC)を出す
こちらは予測確率ですのでY_probaの方を利用します。ROC（Receiver Operating Characteristic）曲線や、そこから算出するAUC(Area Under the Curve)は、各データが正例（今回のケースだと退職する）確率を算出し、確率の大きい順にデータを並べた時にその順序がどのくらい性格であるかを定量化するというものです。
ROC, AUCの詳細を説明するのがこのブログの目的ではないので、とりあえずPythonでどう出すかをご紹介しますね。太字がそのための数式で、それ以降は可視化するためのものになります。
import matplotlib.pyplot as plt
from sklearn.metrics import roc_curve, auc
fpr,tpr, threshholds= roc_curve(y_true=y_test, y_score=Y_proba)
plt.plot(fpr, tpr, label='roc curve (area = %0.3f)' % auc(fpr, tpr))
plt.plot([0, 1], [0, 1], linestyle='--', label='random')
plt.plot([0, 0, 1], [0, 1, 1], linestyle='--', label='ideal')
plt.legend()
plt.xlabel('false positive rate')
plt.ylabel('true positive rate')
plt.show()
結果はこんな感じです。本来予測が100％的中しているなら緑、でたらめに予測すればオレンジ、今回の予測モデルは青のカーブで表されています。
AUCスコアは凡例に書いてある通り0.72です。このスコアについては退職者分析の解析で紹介されていた参照資料によると以下の通りなので、「まあそこそこ」という感じでしょうか。
- 90-1 = excellent
- 80-.90 = good
- 70-.80 = fair
- 60-.70 = poor
- 50-.60 = fail
次は同じ退職者分析を決定木を使って分析を行いたいと思います。
ピープルアナリティクス・ケーススタディ：ロジスティック回帰分析で離職リスク分析に挑戦！①
ロジスティック回帰分析とは
ピープルアナリティクスで離職リスク分析（どういう人が退職のリスクが高いかを過去のデータから導き出す。さらに予測モデルを作るものも）はポピュラーなものの一つです。
組織にとって社員の退職は大きな痛手ですし、事前にどういう人が退職しそうかということがわかれば、少なからず手を打つことができるかもしれません。
ということで、以前Plotly Expressを使ったPythonでのグラフの可視化でご紹介した、離職リスク分析(Turnover Analysis) のケーススタディを用いて実際の離職リスク分析を行っていきたいと思います。
データ分析のための元データもあり、ケーススタディとしてはもってこいですね。
こちらのケーススタディで利用されている分析手法が①ロジスティック回帰、②決定木分析です。
理由は扱いが簡単で、どの要因（統計的に言うと独立変数（従属変数）にあたる項目ですね）が退職に影響を与えているのかを説明しやすいからということでした。
実は、ロジスティック回帰は2年前ほどにPythonを使ったデータ分析講座で習った時にいまいちピンときませんでした。
単回帰や重回帰分析のようにＸ， Y で示される計算式だと何となく想像つくのですが、「0か1かで判断する回帰分析手法です」といわれてもピンとこないですよね。。
この分析を通じてやっとその便利さをわかったような気がします！
0か1かで判断するということは「2つの値に分ける」ということです。
例えば
＊男性か女性か
＊（ウィルスが）陰性か陽性か ＊今ホットトピックですが。。
＊退職したか在籍中か
ちょうど前回のブログで前処理で0か1かにpythonで数値化する方法を書きましたが、「ダミー変数を考慮した分析ですと書いてあった本」＊1もありました
ここではロジスティック回帰の根本的な原理については説明しませんが、以下のような前提を満たすことが必要です。
- 従属変数は2つの値であること（いわゆる0、1ですね）
- それぞれの観察データは独立していること（影響を与えないこと）
- 従属変数同士の値に相関性がないこと
- 独立変数とロジット変換した確率には直線的な関係があること
- サンプル数が多く必要：最低限必要なイベント発生（退職分析の場合退職者）＝従属変数の数 × 10
最後の前提を考えると、変数を絞っていくということも重要ですね。実際ケーススタディの説明を見ていても、一度関係がありそうな従属変数群で分析を行って、さらに絞るという手法を使っています。
ただし、単回帰・重回帰分析と比べて、従属変数・独立変数が正規性に従わなくてもいい、残差のばらつきを考慮しなくてもいいなど、自由な部分もあります。
それにしてもそれぞれの分析手法での注意点はきちんとまとめていかないと分からなくなりますね。。😅
＊1参考本
Pythonでのロジスティック回帰はやり方がいろいろあるが、一長一短。。
Pythonでのロジスティック回帰はいろいろな方法があります。おそらくもっとあるのでしょうが、statsmodelを使った事例とScikit-learnを使った事例を紹介します。
１．statsmodel.apiのロジット関数を使う
オリジナル文献：statsmodels.formula.api.logit — statsmodels
以下のような式で結果を出してくれます。
import statsmodels.formula.api as smf
import statsmodels.api as sm
logistic_model = sm.Logit(y, sm.add_constant(x)) *sm.add_constant(x)はxに定数項を追加するために必要
result=logistic_model.fit(disp=0)
print(result.summary())
*xは独立変数、yは従属変数です。以下同様
結果サンプル
利点
- それぞれの従属変数のp値も出してくれるため、どの変数が有効かどうかが一目でわかる＝＞次のロジスティック回帰分析をするために有効な数値
- 説明変数の列を式で明記しなくてもきちんと表示してくれる。
難点
- 予測モデルを作って検証することが難しそう。。できないことはないのでしょうが、次に紹介する、Sckit-learnのLogisticRegression関数よりずっとめんどくさそうです。（少なくとも今の私の時点ではできませんでした。。）
２．Scikit-learnのLogistic Regression関数を利用
オリジナル文献：sklearn.linear_model.LogisticRegression — scikit-learn 0.24.2 documentation
参考文献：Lec76-79_ロジスティック回帰 (tsjshg.info)
以下のような結果で出ます。
from sklearn.linear_model import LogisticRegression
logistic_model2=LogisticRegression()
result2= logistic_model2.fit(x, y)
summary2 =pd.DataFrame([x.columns, logistic_model2.coef_[0]]).T
summary2
statsmodelと数値の大小について傾向は似てると言えますが（特にp値的に有意なもの）、値が違いますね。。n参考文献にあるように機械学習だと考えると仕方がないのかしら。。
利点
- Accuracyスコアや予測モデルを作成していくことが簡単（作成したモデルのフィットがわかりやすい）
難点
- p値がでない。。
ロジスティック回帰の理解があいまいな中でこんなことをいうのはおこがましいかもしれませんが、各変数の有意性をstatsmodelで確かめ、変数を絞った後、Scikit-learnでモデルの正確性を確かめるというのが適切なのでしょうか。。
修業は続きます。
Pythonでデータ分析に必要な前処理：数値化
統計分析で一番時間がかかるのが前処理。。
Pythonでデータ分析を行うときに一番時間がかかるのが前処理です。
多くの分析の下となるデータはデータ分析がしやすいように集められていないことが普通です。
まず、データ分析でとりわけ回帰分析を行うにおいて必要なのがデータを数値化や数値の重みを統一化すること。でも人事関連のデータには数値でないものがたくさん含まれていますし、数値の値もバラバラだったりします。
それをどう数値化するか？今受講中の講座でいただいたヒントを基にPythonでの数値化の方法と正規化について自分も備忘録もかねてまとめたいと思います。
カテゴリー変数を0か1かに変換（map変数を活用）
Pythonにはpandasにダミー変数化できる関数(get_dummies)があり、こちらを使うことが多いかもしれませんが、この関数を使うと1）変数が増える、2）また2択といった相関性が明らかなものの場合回帰分析には使えなくなります。
まあ、後で削除すればいいだけといえばそうなのですが、別のやり方としてはmap関数を用いて、同じ列である場合は0、別の場合は1という風に0か1に変換します。
例えば採用区分に以下のようなデータが入っていた場合、新卒を0、中途を1に変換するということですね。
こちらはmap関数を使えば2行で簡単にできます。
hiring_map={"新卒": 0, "中途": 1}
df["採用区分"] = df["採用区分"].map(hiring_map)
注意：dfはデータが入っているDataFrameの変数名です。以下同様
私は列数を増やしたくなかったので、採用区分の列を上書きして0、1に変換しました。違う列にしたい場合は列名を変更してくださいね。
この方法ですが、もとのデータが3つ以上の分類のものを0、1の2つの分類に変えたいときにも使えます。
例えば、役職名データを管理職（課長クラス、次長クラス）であれば1、非管理職であれば0というような場合ですね。
その場合は以下のように辞書の定義（以下の場合ですと、mgr_mapという変数に入っている辞書です）をそのように設定すればいいだけです。
mgr_map={"主担当": 0,"副担当": 0, "課長クラス": 1,"次長クラス": 1,"部長クラス": 1}
df["タイトル（管理職有無）"] = df["タイトル"].map(mgr_map)
条件式を使って0か1に変換
0か1の変更は条件式を使って行いたい場合もあるかと思います。
例えば、転職回数が3回以上の人には1というフラグを、それ以下の人には0というフラグを立てる。
その場合はapply関数とlambdaを利用して以下のように書きます。
change=lambda x:1 if x>=3 else 0
df["転職回数3回以上"] = df["転職回数"].apply(change)
さすがに内容が異なるので違う列名で作成しました！
これを応用して2つの列の比較から新たな0、1のフラグを立てる列を作るということもできます。
以下のように前年と直近の評価を比べて、直近の評価が前年評価より低い場合（直近評価－前年評価が0より小さい場合）は1、そうでない場合は0というような感じです。
これは以下のように書けますね。
df["前年からの評価低下_P"]=df["直近評価"]-df["前年評価"]
p_change=lambda x:1 if x<0 else 0
df["前年からの評価低下"] = df["前年からの評価低下_P"].apply(p_change)
（最小最大）正規化で連続数値の重みを統一化
（最小最大）正規化というのは、数値データからなる項目を最小値が0、最大値が1となるように特徴量を正規化する処理です。
一般的には単に正規化といわれているのですが、私がよく参照させてもらっている”Pythonによるデータ分析の教科書"では、分散正規化（平均を0、標準偏差が1となるように変換する処理。一般的に標準化といわれています）と区別してそのように書いていたのでカッコ書きで使うことにしました（自分も標準化とよく間違うので😋)
数式で書くと以下のような式になります。
X’＝（X-Xmin）/ X max -Xmin
が、実際の例で見てもらった方がよくわかると思います。
人事関係のデータではこのようにアンケートであったり在籍年数、時間外、有給取得など、例えば退職者のプロフィールの予測の回帰分析を作ったりするのに幅の異なる連続数値を使うことがあります。
そういう時に便利な方法ですね。
PythonではScikit-Learnでpreprocessingモジュールというのを使って行います。
データフレームの一部の列を変更したい場合はminmax_scaleが便利ですね。ただし、数値は浮動小数点でなければならないので、整数の場合浮動小数点に変更する処理が必要となります。
from sklearn import preprocessing
a=df['在籍年数'].astype(float)
df['在籍年数']=preprocessing.minmax_scale(a)
ただし、このminmax_scaleを使う方法はその後データの精度を検証するために学習データとテストデータに分類する前に行うのには向かないそうです。。(私にはまったく理由がわからなかったのですが、Scikit-Learnの説明をご覧ください）
MinMaxScaler()を使いつつ一番めんどくさくない方法は。。ということでいろいろググると以下のような形だと比較的楽かなと思いました。
from sklearn.preprocessing import MinMaxScaler
mmsc= MinMaxScaler()
df['在籍年数']=mmsc.fit_transform(df['在籍年数'].values.reshape(-1,1))
for文とかを使えば複数の列も同じようにできるのかもしれませんが、相変わらず苦手で😥..次回の課題にしたいと思います。
Plotly Expressを用いてグラフを可視化
Pythonの可視化で何を使うか？
データ分析を行うにあたって、重要なことの一つがデータのグラフ化するということです。
データを分析するにあたってすぐ分析ツールを使いたくなるのですが、そもそもどんなデータなのか、それを理解するには、なかなか表では理解が難しい場合がありますよね。
グラフにすることで一目で「こういう傾向があるんだ」ということがわかり、例えばどの変数を分析対象にするかや加工が必要なのかがわかったりします。
Pythonでグラフの描画として一般的なのはMatplotlibというライブラリです。
ですが、あまりMatplotlibを使うのが得意ではなく毎回ググったり本を見ながら試行錯誤をしていました。。
ある講座で「Plotly」をお勧めされたので、こちらを今回の取り組みでは使ってみようと決めました。
実は、今回はAIHRで参考ページとして挙げられていたページを参考にTurnover analysis（離職リスク分析）を行おうとしています。
分析に行うにあたってのデータもついており、練習するにはいいのですが、説明がRを使った説明になっているのでこちらをPythonでやってみようと思います。
Plotlyの情報は少ない。。
インストール方法を含めPythonでPlotlyを利用する情報はこちらを見ていただいたらいいかと思います。
ただ、何故か最初はうまくいかなかったのですが、いろいろググって以下を入れることで解決しました！
from plotly.offline import download_plotlyjs,init_notebook_mode
init_notebook_mode(connected=True)
ちなみに私はAnacondaからJupyter Notebookを利用してPythonを使っています。
Plotlyはやはりマイナーなのか、特に日本語での情報は少ないです。。なので、今回のまとめもむしろ自分の備忘録を兼ねてまとめようと思いました。
Plotly Express で描画してみる！
Jupyter Notebookで描画するにあたり、できるだけシンプルにコードを書きたかったので何となく楽そうな？”Plotly Express”を今回は用いました。
いくつかグラフを書いてみましたが、パターンが簡単なので使いやすかったです。
ただ、色を変えたりはできないような気がしますが。。
ヒストグラム
>|python|
import plotly.express as px
graph1=px.histogram(df, x="perf",title='Distribution of Performance Rating')
graph1.update_xaxes(title_text="performance")
graph1.update_yaxes(title_text="# employees")
graph1.update_layout(bargap=0.2)
graph1.show()
||<
結果はこんな感じです。
graph1=px.histogram(df, x="perf",title='Distribution of Performance Rating')
- pxの後に描画したいグラフの種類（今回はヒストグラム）
- ()の中は1）データの変数名、2）X軸に入るdataframeの列名、3)タイトル名
graph1.update_xaxes(title_text="performance")
graph1.update_yaxes(title_text="# employees")
- X軸と、Y軸の名称を更新（そのままDataframeの列名を作るなら不要）
graph1.update_layout(bargap=0.2)
- ヒストグラム同士に間隔をあけるための操作
今回はこれが肝でした。。デフォルトだと以下のようにもっさりとしたグラフになります。。特にパフォーマンスは１と２と３しかないのに1.5や2.5にも値がありそうな誤解を生みますよね。
graph1.show()
グラフの表示
一見するとPlotlyのメリットはわかりにくいのですが、以下のようにそのグラフにカーソルを置くとそのデータ情報を見ることができます。
グラフで情報量が多いとこういうインタラクティブさはとてもいいですよね。
横並び棒グラフ
「どんなグラフ？？」て思われたかもしれないのでまずは完成形から。
カテゴリーごとで（今回の場合は性別）棒グラフを横並びにしたいというものですね。
こちらはこんな感じです。。（plotly expressはImportされているので省略されています）
>|Python|
graph5=px.bar(pro_leave_ba_sex, x="area", y=0, color="sex",barmode="group",title='Voluntary Termination Rate by Business Area')
graph5.update_xaxes(title_text="Business Area")
graph5.update_yaxes(title_text="Proportion Leaving")
graph5.show()
||<graph5=px.bar(pro_leave_ba_sex, x="area", y=0, color="sex",barmode="group",title='Voluntary Termination Rate by Business Area')
- px.bar(棒グラフなので今回はbarになっています）
- ｙはgroupbyしたnumpy arrayをDataframe化した結果、元の列名が0になっているからです。。見苦しくてすみません。
- color: 横並びで表示したいDataframeの列名です。今回は性別なので”sex”を選んでいます。
- barmode="group": このグラフではここに苦労しました！これがないと横並びでなく積み上げ棒グラフになります。このグラフでは意味がないですよね。。
barmode="group"がない場合は積み上げグラフになる
箱ひげ図
こちらはそんなに難しくないですね。pxの後をboxにすると簡単に描けます。
>|python|
graph7=px.box(df, x="role", y="age",title=)
graph7.update_xaxes(title_text="role")
graph7.update_yaxes(title_text="age")
graph7.show()
||<
またはPlotlyで箱ひげ図を書くと、四分位数や最小値、最大値などもマウスを動かすだけで見れるのはいいですね。matplotlibではこうはいきません。
にしても使わないとすぐ忘れるので、ちょこちょこ使っていくことが重要ですね。。
ピープル・アナリティクスケーススタディ（補足②）：Pythonで偏相関係数を計算する
2つの変数に本当に相関がある？
以前のブログで相関係数をPythonで計算しましたが、この時どうしても拭い去れなかった思いが、「この係数は2者間の純粋な係数を算出しているのか」ということです。
人事のデータを扱う際には、いくつもの変数の相関性を一気に出したくなる場面があります。ただ、それぞれの変数が完全に独立しているというようなことがはっきり言えない場合が多いと思うので、出てきた相関係数を信じていいのか？という疑いがどうしても残ります。
その疑いをはらしてくれそう？なのが偏相関係数です。
以下のように第三者の影響を除いた「見かけの相関」を取り除いて相関係数を出してくれるというものですね。
Pythonで偏相関係数を出してみる！
実装はいつものとおりPython頼みです！
Scipy.Statsあたりに関数あるかな～と思ったのですが、ないみたいですね。。
ということで、またあまたあるブログよりお知恵を拝借して作成しました。
今回参考にさせていただいたのは以下のブログです。ありがとうございます。
まず、普通にピアソンの相関係数を出してみます。data_2に相関を出したいデータがすべて入っていると想定してください。
>|Python|
import pandas as pd
cor_2= data_2.corr()
print(cor_2)
||<cor_2の結果はこちらです。
array([[ 2.06213096, 0.23729669, -0.28963555, -0.34506658, -0.30974387,
-1.01414826, -0.14092344],
[ 0.23729669, 1.29826588, 0.29162403, 0.34677453, 0.2800417 ,
-0.06692941, -0.02397585],
[-0.28963555, 0.29162403, 1.38245585, 0.81203609, 0.00587289,
-0.05940651, -0.01374919],
[-0.34506658, 0.34677453, 0.81203609, 1.96962014, -0.25722581,
-0.36393852, -0.36610405],
[-0.30974387, 0.2800417 , 0.00587289, -0.25722581, 1.28210144,
0.00363923, 0.21239925],
[-1.01414826, -0.06692941, -0.05940651, -0.36393852, 0.00363923,
1.92073679, -0.29606956],
[-0.14092344, -0.02397585, -0.01374919, -0.36610405, 0.21239925,
-0.29606956, 1.29148533]])
そのあと、偏相関係数行列を出す関数を作成し（cor2pcor)、pcorに入れてその値を返します。正直あまり計算式はよくわかっていません。。
>|python|
import scipy as sp
from scipy import linalg
def cor2pcor(cor_2):
inv_cor = linalg.inv(cor_2)
rows = inv_cor.shape[0]
regu_1 = 1 / sp.sqrt(sp.diag(inv_cor))
regu_2 = sp.repeat(regu_1, rows).reshape(rows, rows)
pcor = (-inv_cor) * regu_1 * regu_2
sp.fill_diagonal(pcor, 1)
return pcor
pcor = cor2pcor(cor_2)
pcor
||<出てきた結果はこちらです。
array([[ 1. , -0.14502796, 0.17154107, 0.17121962, 0.19049492,
0.50957651, 0.08635359],
[-0.14502796, 1. , -0.21767866, -0.21685734, -0.21705995,
0.04238393, 0.01851601],
[ 0.17154107, -0.21767866, 1. , -0.49210577, -0.00441129,
0.03645646, 0.0102898 ],
[ 0.17121962, -0.21685734, -0.49210577, 1. , 0.16186841,
0.18711253, 0.22954538],
[ 0.19049492, -0.21705995, -0.00441129, 0.16186841, 1. ,
-0.00231907, -0.16506197],
[ 0.50957651, 0.04238393, 0.03645646, 0.18711253, -0.00231907,
1. , 0.1879815 ],
[ 0.08635359, 0.01851601, 0.0102898 , 0.22954538, -0.16506197,
0.1879815 , 1. ]])
これだけだとよくわからないので、相関係数及び偏相関係数をヒートマップで表してみましょう。
>|Python|
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
sns.heatmap(cor_2, annot=True, square=True, vmin=-1, vmax=1, fmt=".2f", cmap="jet")
plt.savefig("cor_2.png")
plt.show()
||<ピアソンの相関係数のヒートマップはこんな感じです。
同じように偏相関係数をヒートマップで表すと。。
>|python|
sns.heatmap(pcor, annot=True, square=True, vmin=-1, vmax=1, fmt=".2f", cmap="jet")
plt.savefig("pcor.png")
plt.show()
||<タイトルが上手く表せていませんが、ピアソンの相関係数とはだいぶん違うことがお分かりいただけるかと思います。
数学の知識不足で完全に理解はできていませんが、偏相関係数という概念が理解できたことで少しすっきりしました！
ピープル・アナリティクスケーススタディ（補足①）：Pythonで正規性を確認する
正規性を前提とした分析手法はいくつもある！
これまで、Pythonを使ってピープル・アナリティクスのケーススタディをいろいろな分析手法を使ってやってみました。
その方法としての前提で多いのが正規性に従う（＝観測値の分布が正規分布である）というものです。
正規分布とは、平均値が一番値が高く、そこからなだらかに減少していく以下のような分布をしています。
平均から標準偏差の１倍まで離れた範囲に全体のデータの約68.3％、2倍までの範囲に95.4％、3倍までの範囲に99.7％までが含まれます。
ここでは深く立ち入りませんが、この正規分布に従うという前提で行われている分析や統計モデルはいくつもあります。
私がAIHRの中にある”Statistics in HR”という講座で習った分析手法や統計モデルで正規分布を前提しているものとしては以下のようなものがありました。
- ピアソンの相関係数
- T-テスト（2つのグループ間における平均の違い）
- ANOVA（分散分析）の説明変数
- 単回帰分析
- 重回帰分析の説明変数・被説明変数の誤差
以前のブログで「正規性の確認はヒストグラムと歪度でいい」ということを述べましたが、より厳密に正規性を確認する方法をGW中に学びましたのでこちらをご紹介したいと思います。
この記事の参考にさせていただいた本及びブログは以下のものです。
以下の本ではモデリングで気を付けるべきことや以前から疑問に思っていたことをわかりやすく解説されていました。正規性についてはT-テストで触れられてのみでしたが、その注釈で検定手法が書かれていたのが正規性の実装につながる一歩となりました。
相関係数の検定方法や、偏相関係数についても学ぶことができたので、こちらはまた別の機会に書きたいと思います。
以下のブログはPythonでの実装に参考（というよりそのままコピペ！）にさせていただきました。ありがとうございます！
ヒストグラム以外にも確認方法がある！
「データ解釈学入門」の注釈にあったのが、「シャピロウィルク検定」と「コルモゴロフ・スミルノフ検定（K-S検定）」です。
こちらは両方とも「帰無仮説が有意水準で正規分布に従う」ということを前提としているため、P値が設定した有意水準（例えば0.05）より大きければ帰無仮説を棄却できない＝正規分布に従うということができるということだそうです。
また「シャピロウィルク検定」と「コルモゴロフ・スミルノフ検定（K-S検定）」の違いは、前者がデータ数が小さい（1000以下）に利用し、後者が大きいデータ数に利用するというような状況ということです。
以前のブログでヒストグラムと歪度で正規性を確認したイノベーションの評価についての正規性について出してみましょう。
Pythonでの出し方はとても簡単です。 シャピロウィルク検定についてはscipy.statsでshapiroメソッドを出せば簡単に出ます。
>|Python|
import scipy.stats as stats
stats.shapiro(df["MEANINNOVATION"])
||<結果は以下の通り。前者が「W」と略される検定統計量と、後者がp値になります。
(0.9798566699028015, 0.20555362105369568)
P値は0.2 ＞ 0.05ですので、有意水準を0.05とした場合正規性を棄却できない＝正規分布であるといえるということになります。
続いて コルモゴロフ・スミルノフ検定もご紹介しておきます。今回データ数はそこまで多くないのでシャピロ検定でいいかと思いますが、出し方は同じくkstestメソッドを使ってとても簡単に出ます。
>|python|
stats.kstest(df["MEANINNOVATION"], "norm")
||結果は以下のような感じです。6.668048395734443の10の‐98乗だからこちらは棄却されちゃいますね。。
KstestResult(statistic=0.9286802359624794, pvalue=6.668048395734443e-98)
もう一つ、ご紹介したブログにあったのがQ_Qプロットです。
こちらは前回のブログの残差診断の右上にあった方法と同じです。
すなわち赤い線に従うほど正規性が高いというものです。
前回は4つ一度に出す方法だったので、今回はQ_Qプロットだけ出す方法で見てみたいと思います。
>|Python|
import seaborn as sns
import matplotlib.pyplot as plt
sns.set(font="Hiragino Maru Gothic Pro",context="talk")
fig = plt.subplots(figsize=(8,8))
stats.probplot(df["MEANINNOVATION"], dist="norm", plot=plt)
plt.show()
||<結果はこんな感じです。ちょっと上の方が外れているのが気になりますがほぼ正規性に従うといっていいのかな？？
終わりに
以前正規性の確認で出したのが以下のヒストグラムです。このヒストグラムをみて「正規性がある」と断定するのは迷いがありますが、特にシャピロウィルク検定での結果は「大丈夫！」と念押しをしてくれるいいツールだと思いました。
以前のブログでも書いた通り、データの分布が「正規分布でないことも多い！」ということも念頭に置いておくことが重要なのでこうした検定は、手間であったとしてもやったうえで分析することが必要だと思いました。統計手法はデータをできる限り正確に分析するための手段であり目的ではないですからね😁
ピープル・アナリティクスケーススタディ：重回帰分析に挑戦！③残差診断
残差診断のためのコマンドlmdiagをインストール
いよいよ、重回帰分析を使ったピープル・アナリティクスも最終回です。
今回は、最初の回で重回帰分析を使うための前提の3つ目である「残差」について取り上げたいと思います。
3つの前提は以下のようになっていました。
1．説明変数と被説明変数が線形関係にあること＝＞散布図を描いて確認
２．多重共線性がないこと＝＞説明変数同士で相関係数0.6以上の相関関係がない
３．残差にばらつきがないこと
３について私が今回の課題を解いたAIHRではRを用いて4つの指標をグラフ化していました。
Pythonで同じことを行う場合、”lmdiag”というコマンドをインストールしないといけません。
以下のようにすぐにインストールできるかなと思ったのですが
"linearmodels 4.24 has requirement statsmodels>=0.11, but you'll have statsmodels 0.10.0 which is incompatible"
というメッセージが出てきました。
どうもstatsmodelが古かったようですね。
ただ、そのままstats modelをアップデートすればいいのではなく、pip自体をまずアップデートしなければなりませんでした。
なのでまずpipをアップグレード
statsmodelを再(?)インストール
で再度lmdiagをインストール
でうまくいきました！！
なれている方にはたいしたことではないんでしょうが、私のような万年ひよっこPython利用者にはそれなりに大変でした😅。
4つのグラフはすぐに出ました！
AIHRの解説例で出ていたRを用いた4つのグラフは以下のような形ですぐ出ました！
>|Python|
import lmdiag
plt.figure(figsize=(8,7))
lmdiag.plot(res_4)
pass
||<注:res_4は重回帰分析を行った結果を格納している変数です
左上：残差に非線形のパターンがないか （単独で出す場合の関数は resid_fit())
こちらが上向きやU字型になっていると非線形の可能性が示唆されるそうですが問題なさそうですね。理想は示されている赤線が横軸に平行になっている状況です。
右上：残差が正規分布に従っているかどうか (単独で出す場合の関数は q_q())
正規分布に従っていれば直線状ということですがこちらもほぼに直線であることが確認できました。
左下：残差が均一分散かどうか (単独で出す場合の関数は scale_loc())
縦の散らばりがおおむね均等であれば均一分散といえるそうです。示されている赤線が横軸に平行になるのが理想だそうです。
右下：大きい観測値を特定し、推定値が少数の観測値に大きく左右されていないか (単独で出す場合の関数はresid_lev())
大きな外れ値によって影響を受けていないかですね。特に右上にちらっと見える点線の外に出ているような外れ値がある場合は問題があるようです。
この分析にあたっては以下のブログに大変お世話になりました。正直考えの背景まで理解できているのかは不明ですが今回の課題取り組みを通じて何に気をつけないといけないかということが理解できてよかったと思っています。
これからも実践を通じていろいろ試したいと思います！
ピープル・アナリティクスケーススタディ：重回帰分析に挑戦！②ステップワイズ法で説明変数を絞る！
重回帰分析で出した変数を調整していく
前回のブログのお題で重回帰分析を行いました。
しかし、ここで終わりではありません。前回は学術研究で関連性のあると思われる指標をすべて説明変数に加えて分析を行いましたが、各変数のp値を見ても5%以下となっているものは”MEANTI”のみです。
では、どう調整するのか？AIHRの”Statistics in HR”のこの課題の解説ビデオで紹介されていたのが、AICという指標を参考にする方法です。
AICというのは（Akaike’ s Information Criterion：赤池情報量基準）のことで、元統計数理研究所所長の赤池弘次さんという日本人の方が考案したものだそうです。
Wikipediaで調べると、「モデルの複雑さとデータの適合度のバランス」を取るためにモデルの良さを評価するための指標だということです。
変数（や次数）を増やすと、分析したデータには適合するけれどもじゃあ、同じような状況にこのモデルが使えるかというとかが適合（Overfitting）が起こり、使えなくなる。
どれくらいの数に抑えるかということの一つの解を与えてくれるのがAICで最小のモデルであれば多くの場合良いモデルが選択できるそうです。
また変数の調整方法ですが、ステップワイズ法を用いていました。ステップワイズ法には、増減法（forward-backward:最も単純なモデルから変数を増やしていって探していく）、減増法(backward-foward:すべての変数を投入してそこから減らしていき最適な変数を探していく）がありますが、今回は学術的に考えうる変数をすべて使って行った重回帰分析の結果を使うので、減増法で行います。
ステップワイズ法をPythonで実践するのは大変。。
しかし、Pythonには利用できるパッケージ等を含めてもステップワイズ法を実装できる関数がありません（見つかったら教えてください）。。
ということでこちらのブログを参考にさせていただき、自作関数を作りました。
私のPythonのレベルからですとゼロからこんな関数は作れないので、本当に感謝、感謝です🙇♂️。
===============================================
desc_list = ['CLGENDER', 'CLWORKEXP', 'CSMnetworking', 'CSMmobility', 'MEANTI', 'MEANPERSONALINIT']
def stepAIC(descs_l):
import copy
descriptors = descs_l
f_model = ols(formula='MEANINNOVATION~ ' + ' + '.join(descriptors), data=data_2).fit()
best_aic = f_model.aic
best_model = f_model
while descriptors:
desc_selected = ''
flag = 0
for desk in descriptors:
used_desks = copy.deepcopy(descriptors)
used_desks.remove(desk)
formula = 'MEANINNOVATION ~ ' + ' + '.join(used_desks)
model = ols(formula=formula, data=data_2).fit()
if model.aic < best_aic:
best_aic = model.aic
best_model = model
desc_selected = desk
flag = 1
if flag:
descriptors.remove(desc_selected)
else:
break
return best_model
model_a = stepAIC(desc_list)
print(model_a.summary())
===============================================
注意：
この間の続きですのですでに”from statsmodels.formula.api import ols”をimport済みです。
でてきた結果がこちらです。Rを使って分析したAIHRの解説ビデオとも合致しました!
変数は6つなったものから4つに絞られましたね。
どの値もP値が改善され、一番大きいCLWORKEXPでも5.6％となりました（なぜかサマリーの表が古めかしくなりましたが。。）。
ただ、AICについてはもともと191.5だったので、189.9への変更がどこまで大きいといえるのかは難しいところですが。。（まだまだ勉強が必要です）
次回は重回帰分析で分析できる前提の最後にあった”残差のばらつき”についてご紹介したいと思います。こちらもPythonでは苦労しました。。
ピープル・アナリティクスケーススタディ：重回帰分析に挑戦！①
いよいよ？重回帰分析
前回は、人事の課題をPythonで単回帰分析を行い、分析する事例を紹介しました。
今回は、いよいよ重回帰分析に挑戦したいと思います。
といっても重回帰分析と単回帰分析の大きな違いは説明変数が1つではなく2つ以上であるということのみ。
とはいうものの、普通ある出来事を予測する際に一つの変数のみで説明することができないのが当然なので、実際には重回帰分析で行うことが多いのかなと思います。
単回帰分析と同様、こちらも線形モデルばかりではないのですが、今回は線形モデルを作りたいと思います。
お題は以下の通り。
イノベーションを起こす力を持っている人を予測できるモデルを作ってください。過去の学術研究によると、性別、年齢、勤務年数、社歴、ネットワーキング活動、組織に関係なく人をつなぐ活動、自主性、異動の多さなどが影響があるといわれています。
お題及びそれに基づくデータはあくまでも私が受講しているAIHRの”Statistics HR”講座にあった課題の中で与えられたものですので、今回紹介する分析結果をうのみにしないでくださいね😅。
人事の事例としてこういう分析方法があるんだという参考になれば幸いです（というより自分の備忘録に近いのですが。。）。
重回帰分析の前提条件（の一部）を確認！
重回帰分析での最初のチャレンジは、どの要素を説明変数として取り入れるかということです。
今回は素直に、過去の学術研究に従った変数のみでまずは分析を行うことにします。
以下のような形で、全体のデータが入ったdfから被説明変数である”MEANINNOVATION”と今回説明変数として取り上げる6つの要素（"CLGENDER", "CLWORKEXP", "CSMnetworking", "CSMmobility", "MEANTI", "MEANPERSONALINIT"）に絞って新しいデータフレーム「df_2」を作成します。
>|python|
data_2=df.loc[:,["MEANINNOVATION", "CLGENDER", "CLWORKEXP", "CSMnetworking", "CSMmobility", "MEANTI", "MEANPERSONALINIT",]]
||<
講座の解説ビデオによると、確認する前提としては3つ挙げられていました。
1．説明変数と被説明変数が線形関係にあること
２．多重共線性がないこと
３．残差にばらつきがないこと
3は結果を出してからのみしか調べられないのでまた次回以降に説明します。
まず１．説明変数と被説明変数が線形関係にあることは散布図を描いて確かめるのが手っ取り早いです。
Seaboornのpairplotを用いて散布図を取り出します
>|python|
import seaborn as sns
sns.pairplot(data_2)
||<すると以下のように出てきます。一番左はイノベーション行動のデータの分布ですが、他はイノベーション行動との相関性ですね。
左から2番目は性別なので０か1かに偏っています。ただ、グラフを見る限りU字型になっているような散布図はないので、とりあえず線形モデルで試していいということです。
２．多重共感性がないことですが、解説ビデオでは、説明変数同士で0.6以上の強い相関性がある場合は排除をした方がいいという程度の緩い？基準を満たしていればいいしたのでどちらにせよ相関係数を見てみることにしましょう。
>|python|
import pandas as pd
cor_2= data_2.corr()
print(cor_2)
||<
被説明変数との相関係数も出ていますが、いずれも±0.6以上の相関性はないので重回帰分析をすすめましょう！
Statsmodelでの重回帰分析はとても簡単！
重回帰分析は前回の単回帰分析と同じstatsモデルで行います。
式も"被説明変数 ~ 説明変数（それぞれを +でつなげる）"で簡単に設定することができます。
>|python|
from scipy import stats
from statsmodels.formula.api import ols
formula_2 = "MEANINNOVATION ~ CLGENDER + CLWORKEXP + CSMnetworking + CSMmobility+ MEANTI + MEANPERSONALINIT"
res_2 = ols(formula_2,data = data_2).fit()
res_2.params
||<
重回帰分析の式となる係数をparamsで出していますので、結果は以下の通りとなります。
Intercept -0.130093
CLGENDER -0.227898
CLWORKEXP 0.014701
CSMnetworking 0.266372
CSMmobility 0.171613
MEANTI 0.436122
MEANPERSONALINIT 0.160258
詳細の分析結果は単回帰分析と同様.summary()で出します。
>|python|
res_2.summary()
||<
結果はこんな感じです
決定係数は、重回帰解析の場合調整済み決定係数（Addj. R-squared）を見るのがいいそうですので、0.478ということですね。
各説明変数のP値を見ると0．1以下は２つ（CSMmobility(異動の多さ）、MEANTI(組織に関係なく人をつなぐ行動））のみですね。。
ではここからどう調整するかを次回観ていきたいと思います。
