引き続きロジスティック回帰について調べている。
前回、ロジスティック回帰の最尤推定量にはバイアスがあることを調べた。
このバイアスは、サンプルサイズが大きい場合は無視できるが、入力変数の数に対してサンプルサイズが小さい場合には無視できないほど大きくなる。
今回はサンプルサイズをさらに小さくすると起こる問題について考える。
前回の最初の例を少しだけ変えて実行してみよう。
この例では、入力変数の数 として、パラメータ 300個の真の値を、最初の 100個は
、次の 100個は
、残りの 100個は
と設定した。
ここまでは同じだが、変更点として、前回はサンプルサイズ
だったのに対して、今回は
で実行する。
n <- 1300 p <- 300 # データの生成 set.seed(314) x <- rnorm(n * p, mean = 0, sd = sqrt(1/n)) X <- matrix(x, nrow = n, ncol = p) beta <- rep(c(10, -10, 0), each = p/3) prob <- plogis(X %*% beta) y <- rbinom(n, 1, prob) # ロジスティック回帰モデルの適用 fit <- glm(y ~ X, family = binomial, control = list(maxit = 50))
警告メッセージ: glm.fit: 数値的に 0 か 1 である確率が生じました
df <- data.frame(index = seq_len(p), coef = fit$coefficients[-1]) library(ggplot2) theme_set(theme_bw()) ggplot(df, aes(index, coef)) + geom_point(color = "blue") + annotate("segment", x = c(1, 101, 201), xend = c(100, 200, 300), y = c(10, -10, 0), yend = c(10, -10, 0), size = 1.5) + xlab("Index") + ylab("Coefficients (true and fitted)")
実行結果は見てのとおり、推定されたパラメータが 1e+16 (= 1016) などと、とんでもない値になっていることがわかる。 これは最尤推定量のバイアスとは異なる原因で起こっている。
何が起こっているかというと、実は、このデータに対して、ロジスティック回帰モデルには最尤推定量が存在しない。 ロジスティック回帰モデルには最尤推定量が存在しない場合があるのである。 しかし、R言語の最尤推定値を求めるアルゴリズムは、最尤推定量が存在しない場合にも、なんらかの値を返してしまう。
まとめると、入力変数の数に対してサンプルサイズが小さい場合、推定値には2つのケースが存在する。
- 最尤推定量が存在する場合、バイアスのある推定値
- 最尤推定量が存在しない場合、なんの意味もない値
この記事では、この2つのケースを判別する方法、すなわち、ロジスティック回帰に最尤推定量が存在するかどうかを判定する方法を紹介する。
判定材料
ロジスティック回帰に最尤推定量が存在しない場合、パラメータの推定値が異常な値を取るため、すぐにわかるのではないか？と思ってしまうが、そうでもない。
上記の例でサンプルサイズ を変えて推定されたパラメータの最初の100個 (
) がどのような値を取るかを見てみよう。
でだけ、推定値が異常な値を取ることがわかる。
の小さい範囲を拡大してみると、パラメータが大きめに推定されていることがわかるが、これは真の値が 10 であることを知っているため大きいとわかるのであって、真の値を知らない場合はこの推定値が大きいかどうかの判断は難しい。
判定方法
ロジスティック回帰に最尤推定量が存在するかどうかを判定する方法は、40年ほど前に研究された [Albert & Anderson, 1984][Santer & Duffy, 1986]。
文献 [Albert & Anderson, 1984] によれば、観測データが完全分離 (complete separation) または準完全分離 (quasi-complete separation) の場合に限り、最尤推定量は存在しない。
完全分離とは、入力変数が作る空間に、出力変数を完全に分離できる超平面が存在することを言う。
すなわち、入力変数 、出力変数
に対して、
\begin{cases} Xb > 0 & y = 1 \\ Xb < 0 & y = 0 \\ \end{cases}
を満たす が存在することを言う。ここで、
\begin{align} 2y - 1 = \begin{cases} 1 & y = 1 \\ -1 & y = 0 \\ \end{cases} \end{align}
であるので、
\begin{align} Xb \times (2y - 1) > 0 \end{align}
が成り立つ が存在すればよい。
これは、線形計画問題に実行可能解があるかどうかで判定できる。
ROI パッケージを使えば線形計画問題を簡単に解くことができる。
# ROI パッケージで線形計画問題を解く library(ROI) A <- cbind(X, 1) * (2 * y - 1) b <- rep(0.001, n) c <- rep(1, p + 1) lp <- ROI::OP(objective = ROI::L_objective(c), constraints = ROI::L_constraint(A, rep(">=", n), b), bounds = ROI::V_bound(lb = rep(-Inf, p + 1))) solution <- ROI::ROI_solve(lp) # 実行可能解 (feasible solution) があれば完全分離可能 status_code <- solution$status$msg[["code"]] STATUS_CODE_NO_FEASIBLE_SOLUTION_EXISTS <- 4L STATUS_CODE_SOLUTION_IS_OPTIMAL <- 5L STATUS_CODE_SOLUTION_IS_UNBOUNDED <- 6L is_separable <- if (status_code == STATUS_CODE_NO_FEASIBLE_SOLUTION_EXISTS) { FALSE } else if (status_code == STATUS_CODE_SOLUTION_IS_OPTIMAL) { TRUE } else if (status_code == STATUS_CODE_SOLUTION_IS_UNBOUNDED) { TRUE }
ここでは簡易的に完全分離を判定する方法を書いたが、完全分離、準完全分離、そしてオーバーラップ（最尤推定量が存在する）のいずれであるかを判別する線形計画問題の定式化が文献 [Santer & Duffy, 1986] に載っているので参照してほしい。
文献 [Konis, 2007] では、より洗練された方法が提案されている。 これを実装したのが detectseparation パッケージである。 このパッケージでは完全分離または準完全分離かどうかを簡単に判定できる。
library(detectseparation) result <- detect_separation(X, y, family = binomial()) is_separable <- result$outcome
完全分離または準完全分離であれば最尤推定量は存在しないので、これでロジスティック回帰に最尤推定量が存在するかどうかを判定できる。
参考文献
- Albert, Adelin & Anderson, J.A.. (1984). On the Existence of Maximum Likelihood Estimates in Logistic Regression Models.
- Santner, Thomas & Duffy, Diane. (1986). A Note on A. Albert and J. A. Anderson's Conditions for the Existence of Maximum Likelihood Estimates in Logistic Regression Models.
- Konis, K.P. (2007). Linear programming algorithms for detecting separated data in binary logistic regression models.
