ロジットを用いて、複数の説明変数から二値データを分類するモデルのパラメータを決定する。
RのGLMなんかに解かせると、各説明変数についてP値が算出され、この変数がどの程度有意か、というのが把握できる。
ここでのP値とは何か。ロジスティック回帰によって各 $x_i$ の係数 $\beta_i$ が推定される。するとP値は「もし $x_i$ が実は結果に全く影響なかったとして、$x_i$ の分散度合いなどから $\beta_i$ が“たまたま”今回の値以上となってしまう確率」を表す。変数 $x_i$ が結果を左右する重要な値なのであれば、$\beta_i$ のP値は低くなる。(逆にP値が高いからといって影響が低いとは言えず、この辺が統計的な解釈として難しいところだが、あくまで影響があるという保証が得られなかった、といえるに留まる)
分析内容によるが、一般的な指標としてP値は1%, 5%, 10%以下あたりが基準とされる。ただ機械学習の場合は20%程度でも採用することもあるという意見もあるし、業界の習慣と試行錯誤によるのだろう。P値が基準値以上となった説明変数は、組み込むと過学習を招くため、取り除くことを検討する。(ただし、取り除くことで他の変数への影響が出るため、絶対的な基準とはならない)
で、本題。
sklearn.linear_model.LogisticRegressionには、P値を自動で算出してくれる機能は無い。クラスを拡張して、自前で計算する方法が以下に紹介されている。ただしコメント欄にて微修正が行われているため、コピペした以下のコードのままで動くか未確認。
from sklearn import linear_model
from scipy import stats
import numpy as np
class LinearRegression(linear_model.LinearRegression):
"""
LinearRegression class after sklearn's, but calculate t-statistics
and p-values for model coefficients (betas).
Additional attributes available after .fit()
are `t` and `p` which are of the shape (y.shape[1], X.shape[1])
which is (n_features, n_coefs)
This class sets the intercept to 0 by default, since usually we include it
in X.
"""
def __init__(self, *args, **kwargs):
if not "fit_intercept" in kwargs:
kwargs['fit_intercept'] = False
super(LinearRegression, self)\
.__init__(*args, **kwargs)
def fit(self, X, y, n_jobs=1):
self = super(LinearRegression, self).fit(X, y, n_jobs)
sse = np.sum((self.predict(X) - y) ** 2, axis=0) / float(X.shape[0] - X.shape[1])
se = np.array([
np.sqrt(np.diagonal(sse[i] * np.linalg.inv(np.dot(X.T, X))))
for i in range(sse.shape[0])
])
self.t = self.coef_ / se
self.p = 2 * (1 - stats.t.cdf(np.abs(self.t), y.shape[0] - X.shape[1]))
return self
回帰分析では説明変数は独立であることが理想とされており、互いに相関の高い2変数を同時に組み込むと結果がおかしくなりやすい。いずれかを除くか、統合して1つの指標値とすることが望ましい。
これを検出するには、一般的な回帰分析ではVIFという統計量を見る手法が有効だが、ロジスティック回帰では有効に働かないらしい。また、他にも決定打となる検出方法は確立されていないとのこと。
ロジスティック回帰分析における多重共線性の評価に関して|薬剤師のためのEBMお悩み相談所-基礎から実践まで
とりあえず、ダミー変数の内1つは削除しとくのと、各説明変数毎に相関取ってみて0.9とかだとあやしいよね、くらいと説明されている。