工数の分析にあたり、画面、帳票、バッチの重み付けをシステマティックに行いたいという積年の課題があり、単純に考えると重回帰分析を行えばという話なのだけど、実際にやってみるとまったく予想外の結果がでる。
df= pd.read_csv('./IPA_2014-2019.csv') model = smf.ols('MH ~ VC + RC + BC - 1', df) result = model.fit() result.summary()
coef std err t P>|t| [0.025 0.975] VC 276.7189 23.494 11.778 0.000 230.532 322.906 RC 173.3967 49.852 3.478 0.001 75.391 271.403 BC -15.2673 9.259 -1.649 0.100 -33.469 2.934
なんとバッチの係数がマイナス！ これを信じて工数計算を行うとバッチ本数を大量に増やすことでゼロカロリーならぬゼロ工数を実現できてしまう。*1
色々考えていたところ、uncorrelated 氏に「多重共線性が出ているので、何か説明変数を落とすか、主成分分析で変数を合成するなりした方が良さそう」との助言を頂いた。
ようは変数間の相関が無視できないということだ。この事自体は前から気付いていたのだけど、説明変数を落とすのは工数を決めるための式を求めるそもそもの目的にそぐわないし、主成分分析は軸を変えてフィッティングするという程度しか知識がなく画面、帳票、バッチの数値からどうやって軸を変えて工数を算出すればいいのかよくわからなかったので、手を出すのを躊躇していたのであった。
助言をもらったからには重い腰を上げてやってみようということで、今回のエントリはその作業メモ。
まず、主成分分析だけど基本的には、元の説明変数の軸を変換してよりフィットしやすい説明変数を作り出そうというもの。単に説明変数を作り出すだけだとあまり意味はないけれども、作り出された説明変数を落とすことで、元の説明変数すべての影響を残すことができる。一般的には累積寄与率が7、8割を超えるところまで残すのがよいとされているようだ。
Python で主成分分析をする方法を調べてみると、statsmodels を使う方法と scikit-learn を使う方法があることがわかった。ただ、statsmodels の方はあまりこなれていない印象で、scikit-learn を使う方がよさそうだ。
scikit-learn で主成分分析
from sklearn.decomposition import PCA df = pd.read_csv('./IPA_2014-2019.csv') X = df[["VC", "RC", "BC"]].values pca = PCA() pca.fit(X)
結果は、pca オブジェクト内に保持される。まず、寄与率を調べる。寄与率は pca.explained_variance_ratio_ で取得できる。
array([0.85622803, 0.11874323, 0.02502874])
結果を見る限り、最初の要素だけで8割超えているので、第1主成分だけでほとんど説明できると解釈できる（本当か？）。そこで、第一主成分のみで回帰を行う。
df2 = pd.DataFrame(columns=["C1", "C2", "C3"], data=pca.transform(X)) df2["MH"] = df["MH"] model = smf.ols('MH ~ C1 - 1', df2) # C2、C3 は捨てる ols = model.fit()
回帰分析の結果は、ols.params で取得でき、C1の係数 21.47954 だった。
これを画面、帳票、バッチ数からなる式にする必要がある。transform は、(元データの行列（X） - 元データの列ごとの平均（np.mean(X, axis=0)）) × 主成分行列（pca.components_） という計算を行っているので、同じことを行って主成分に変換すればいい。
今回のデータの場合、最終的には次の結果が得られる。
人時工数 = 2.50 * 画面数 + 0.62 * 帳票数 + 21.32 * バッチ数 - 1459.15
見れば明らかだけれど、この結果も妥当とは言いがたい。係数がマイナスなのもそうだが、画面、帳票、バッチの比率も明らかにおかしい。これは線形回帰、主成分分析に共通する特徴として、外れ値に弱いことが原因と考えられる。
この問題に対応するためにロバストPCAという方法があるものの、scikit-learn にも statsmodels にも組み込まれておらず、野良ライブラリの dganguli/robust-pca を使う必要があるようだ（現段階では調査中）。
[追記] robust-pca を試してみたが、ほとんど外れ値として扱われてしまい元の行列とは似てもつかないものになってしまった。
以下の結果の S が外れ値、Lが外れ値を考慮した行列
import libs.r_pca as r_pca pca3 = r_pca.R_pca(D=df[["VC", "RC", "BC"]]) L, S = pca3.fit()
statsmodels で主成分分析
statsmodels では次の方法で実行できるが、なぜか scikit-learn での結果とまったく一致せず断念してしまった。オプションの違いなのだろうか……
import statsmodels.multivariate.pca as pca df = pd.read_csv('./IPA_2014-2019.csv') result = pca.PCA(df[["VC", "RC", "BC"]])
