ロジスティック回帰の最尤推定量にはバイアスがある

ロジスティック回帰について調べている。

ロジスティック回帰モデルのパラメータの最尤推定量は、不偏推定量ではなく、バイアスがある。

例として、サンプルサイズ $n = 1500$ 、入力変数の数 $p = 300$ のときを考える。パラメータ 300個の真の値を、最初の 100個は $\beta = 10$ 、次の 100個は $\beta = -10$ 、残りの 100個は $\beta = 0$ に設定して推定してみよう。

n <- 1500
p <- 300

# データの生成
set.seed(314)
x <- rnorm(n * p, mean = 0, sd = sqrt(1/n))
X <- matrix(x, nrow = n, ncol = p)
beta <- matrix(rep(c(10, -10, 0), each = p/3))
prob <- plogis(X %*% beta)
y <- rbinom(n, 1, prob)

# ロジスティック回帰モデルの適用
fit <- glm(y ~ X, family = binomial)

# パラメータの最尤推定値の抽出
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)")

最初の 100個のパラメータは、真の値 $\beta = 10$ よりも大きめに推定されている。次の 100個のパラメータは、真の値 $\beta = -10$ よりも小さめに推定されている。最後の 100個のパラメータは真の値 ( $\beta = 0$ ) の周りに均等にばらついている。

つまり、ロジスティック回帰のパラメータの最尤推定量には、絶対値を大きめに推定するようなバイアスが乗っている。

最尤推定量には一致性と漸近正規性があるので、サンプルサイズ $n$ が十分大きければ、このバイアスは無視できるようになる。上記の例でサンプルサイズを変えて同じシミュレーションをしてみよう。

入力変数の数 $p = 300$ に対して $n = 6000$ まで増やせば、最尤推定量のバイアスは無視できるほど小さくなるようだ。

しかし、逆に言うと、入力変数の数に対してサンプルサイズが小さい場合は無視できないほどのバイアスが乗ってしまう。

このバイアスは、統計解析において次の問題を引き起こす。

効果量を過大に見積もってしまう
出力の確率を極端に見積もってしまう

上記の例はパラメータの真の値が作為的すぎるので、もう少し現実的な例で説明しよう。

真のパラメータ $\beta$ を正規分布 $N(3, 4)$ から抽出し、パラメータの真の値と最尤推定値 (MLE) の散布図を描いてみる。

n <- 4000
p <- 800

# データの生成
set.seed(314)
x <- rnorm(n * p, mean = 0, sd = sqrt(1/n))
X <- matrix(x, nrow = n, ncol = p)
beta <- rnorm(p, mean = 3, sd = 4)
prob <- plogis(X %*% beta)
y <- rbinom(n, 1, prob)

# ロジスティック回帰モデルの適用
fit <- glm(y ~ X, family = binomial)

# パラメータの最尤推定値の抽出
df <- data.frame(true = beta, mle = fit$coefficients[-1])

ggplot(df, aes(true, mle)) +
  geom_point(color = "blue") +
  geom_abline(slope = 1, size = 1) +
  xlab("True signal") + ylab("MLE")

黒線は傾き 1 の直線であり、もしバイアスがなければ推定値はこの線の周りに分布するはずである。しかし、明らかに分布がずれていることから、この例でも最尤推定量にバイアスがあることがわかる。

また、散布図の傾きが大きくなっていることから、絶対値が大きくなる方向にバイアスが乗っていることがわかる。

絶対値が大きくなる方向にバイアスの乗った推定値を、入力変数の効果量として解釈すると、その入力変数は過大評価され、誤った結論を導いてしまう恐れがある。

もう一つの問題は、このようなバイアスの乗った推定値で算出された出力の確率は、極端な値になりやすいということである。

出力 = 1 となる確率をパラメータの推定値から算出してみる。

x <- X %*% beta
pred <- predict(fit, type = "response")

df <- data.frame(x = x, prob = prob, pred = pred, pred2 = pred2)
ggplot(df, aes(x, pred)) +
  geom_point(color = "blue", alpha = 0.1) +
  geom_point(y = prob) +
  scale_x_continuous(breaks = NULL, limits = c(-6, 6)) +
  scale_y_continuous(breaks = 0:5 * 1/5, minor_breaks = NULL) +
  xlab(NULL) + ylab("Probabilities (True and predicted)")