概要

Pythonの統計解析statsmodelsで重回帰分析を行っていたら，あるサンプル数以上から回帰係数が正しく出力されないという現象に出会ったので，その原因と対策について共有します．
結論をざっくりいうと，原因は（たぶん）因子であるカテゴリ変数をOne-Hotエンコーディングをしていたことで，対策は因子ごとのダミー変数を1つずつ抜くことです．

起こった問題

A, B, C, D, 4つの因子に対してそれぞれa, b, cという水準(値)を取りうるデータとその評価点（1から7の数値）があるデータに対し, 各因子の値を説明変数，評価点を目的変数としてPythonのstatsmodeslで重回帰分析をすることで，各因子の値が評価点にどのように影響を与えるかを分析していました．（抽象的でスイマセン...．要は得点評価法を用いたコンジョイント分析です）

下記のような感じでカテゴリ変数である因子A, B, C, Dを0または1の値を取るダミー変数(A_a, A_b,...)に置き換えました．いわゆるOne-Hotエンコーディングというものです[参考2]．

df_x2 = pd.get_dummies(df_x[['A','B','C','D']])

置き換えたデータに対して，次のように最小二乗法による重回帰分析を実行した結果, あるサンプル数まではそれっぽい結果を出すのですが，それを超えると不可思議な結果を出すようになりました．

model= sm.OLS(df_y, sm.add_constant(df_x2))
result = model.fit()
result.summary()

サンプル数が842まではこんな感じなのですが...

サンプル数が843を超えると...

こんな感じで, 1サンプル増えただけなのに回帰係数が大きく変わって，しかも同じ因子から作ったダミー変数の回帰係数がまったく同じ値になるという謎現象が起こりました（涙）．この謎の閾値はデータの値によって違いましたが，データ数が増えると似たような現象が起こることは同様でした．

原因と対策

どうにも問題が解決しなかったので，Rの方で同様の重回帰分析を行った結果，次のようになりました．

"NA"という値が含まれていました．これはある説明変数が他の説明変数で完全に説明できることを意味しています[参考3]．重回帰分析では説明変数間に相関があると多重共線性により正しい結果が出ないため，説明変数から除く必要があります（...そうだった）．
※ これとは別に，各説明変数の多重共線性をチェックする必要があります[参考4]. いっぺんにやってくれてもいいのに．

このRが自動でやってくれている処理を手動でやってあげる（データから各因子ごとにダミー変数を1つずつ抜いてあげる）とstatsmodelsでもRと同じ結果を出力するようになりました．

回帰係数が最初の図とずいぶんと変わったように感じるかもしれませんが，これは各説明変数の意味が変わるからです．

今回の教訓は"One-Hotエンコーディングしたデータを，そのまま重回帰分析に喰わすな"ということでした．

おまけ

statsmodesのGithubのissue，なんと2.5k（2500件！？）ありました@@;．Tensorflowをはじめとした最新の機械学習のライブラリを使うならPythonが超絶便利ですが，クラシックな統計解析は，昔からあるRなどの統計解析ソフトを使う方が動作が安定していてよいのかな～と思いました．どなたかのご参考になれば幸いです．

今回はコンジョイント分析のやり方自体は説明しませんでしたが，次回はコンジョイント分析のやり方についてまとめて書いてみようと思います♪

参考

以下のサイトを参考にさせて頂きました．ありがとうございます．

参考1: GitHub - statsmodels/statsmodels: Statsmodels: statistical modeling and econometrics in Python
- https://github.com/statsmodels/statsmodels
参考2: ダミー変数（One-Hotエンコーディング）とは？実装コードを交えて徹底解説
-　https://www.codexa.net/get_dummies/
参考3: r - Rで重回帰分析を行った際のNA値について - スタック・オーバーフロー
- https://ja.stackoverflow.com/questions/35536/r%E3%81%A7%E9%87%8D%E5%9B%9E%E5%B8%B0%E5%88%86%E6%9E%90%E3%82%92%E8%A1%8C%E3%81%A3%E3%81%9F%E9%9A%9B%E3%81%AEna%E5%80%A4%E3%81%AB%E3%81%A4%E3%81%84%E3%81%A6
参考4: 多重共線性をRでチェックする ― RでVIFを計算するには？ - 統計ER
- https://toukeier.hatenablog.com/entry/how-to-calculate-vif-by-r/