注意深いRユーザーはお気づきだと思いますが、glmでロジットモデルを推定したあと結果オブジェクトにconfintをかけると、waiting for profiling to be done...とメッセージが表示されます。
r_glm <- glm(ans ~ x + z, data = df02, family = binomial())
confint(r_glm)
Waiting for profiling to be done...
2.5 % 97.5 %
(Intercept) -0.851407 0.4332838
x 1.263595 2.9298628
z -5.110815 -2.1783024profiling? — 意味がわかりませんよね。?confint.glmとしても、詳しい説明はありません。
左右非対称であったりします。
matrix(c(coef(r_glm), apply(confint(r_glm), 1, mean)),
2, length(coef(r_glm)), byrow = TRUE,
dimnames = list(c("coef", "mean(confint)"), names(coef(r_glm))))
Waiting for profiling to be done...
(Intercept) x z
coef -0.1865190 1.984496 -3.448283
mean(confint) -0.2090616 2.096729 -3.644559これ、プロファイル尤度信頼区間を計算していると言うメッセージです。教科書で説明してあるのは、だいたい標準誤差を元にした左右対称となるワルド信頼区間の計算方法です。
伊藤 (2019)に詳しい説明があるのでそちらを参照して頂きたいのですが、プロファイル尤度信頼区間は尤度比検定に対応した信頼区間です。プロファイル信頼区間を外れる値を帰無仮説に尤度比検定を行うと、帰無仮説は棄却されます。正規分布への漸近性を仮定しないで計算されるので、性質が良いとされています。安心して、デフォルトの値を用いましょう。
話の本題はこれで終わりですが、最尤法でワルド信頼区間とプロファイル尤度信頼区間の計算をしてみます。
DGPからデータセットを作成します。
set.seed(1234)
df02 <- with(new.env(), {
n <- 100
S <- matrix(c(2, 0.5, 0.5, 1), 2, 2)
C <- chol(S)
X <- cbind(1, t(t(C) %*% matrix(rnorm(ncol(C) * n), ncol(C), n)))
beta <- c(0.5, 2, -3)
y <- X %*% beta
p <- 1/(1 + exp(-y))
ans <- 1*(runif(n) < p)
data.frame(ans, x = X[,2], z = X[,3])
})
最尤法に用いる関数
今回は後で使うのでグラディエントも用意します。
mf <- model.frame(ans ~ x + z, df02)
X <- model.matrix(mf, mf)
ans <- model.response(mf)
llf <- function(beta){
y <- X %*% beta
p <- 1/(1 + exp(-y))
min <- 1e-16
r <- sum(ans*log(pmax(p, min)) + (1 - ans)*log(pmax(1 - p, min)))
browser(expr = is.na(r) || !is.finite(r))
r
}
llfg <- function(beta){
sxp <- X %*% beta
apply(sapply(1:nrow(X), \(i) X[i, ] * (ans[i] - 1 + exp(-sxp[i])/(1 + exp(-sxp[i])))), 1, sum)
}
ワルド信頼区間
まず、最尤法をかけます。
r_ml <- optim(rep(0, ncol(X)), llf, llfg, method = "L-BFGS-B", control = list(fnscale = -1), hessian = TRUE)
beta <- r_ml$par
se <- sqrt(diag(-solve(r_ml$hessian)))
次に、係数の標準誤差から信頼区間を求めます。
a <- 0.05
lwrupr <- c(a/2, 1 - a/2)
confint.wald <- t(sapply(1:length(beta), \(i) qnorm(lwrupr, mean = beta[i], sd = se[i])))
rownames(confint.wald) <- names(beta) <- names(se) <- colnames(X)
colnames(confint.wald) <- sprintf("%.1f%%", lwrupr*100)
print(confint.wald)
2.5% 97.5%
(Intercept) -0.8198741 0.4468386
x 1.1614836 2.8075079
z -4.9013764 -1.9951880
プロファイル尤度信頼区間と微妙に異なる値が出ました。
制約付き最尤法に用いる関数
多次元の尤度関数を一次元にしたものがプロファイル尤度関数なのですが、尤度比検定に用いる制約付き最尤法に使う関数が必要になります。
pparam <- function(cp, i, v){
p <- numeric(length(cp) + 1)
if(1 < i) p[1:(i-1)] <- cp[1:(i-1)]
p[i] <- v
if(i < length(p)) p[(i + 1):length(p)] <- cp[i:length(cp)]
p
}
pllf <- function(cp, i, v){
llf(pparam(cp, i, v))
}
pllfg <- function(cp, i, v){
llfg(pparam(cp, i, v))[-i]
}
さらに、ヘッシアンが必要になるので計算する関数を用意します。
hessian <- function(beta){
H <- matrix(NA, length(beta), length(beta))
d <- 1e-12
for(i in 1:length(beta)){
b0 <- b1 <- beta
b0[i] <- beta[i] - d
b1[i] <- beta[i] + d
H[i, ] <- (llfg(b1) - llfg(b0))/d/2
}
H
}
プロファイル尤度信頼区間の計算
まず、最尤法の結果から、プロファイル尤度信頼区間の端点における対数尤度を求めます。
lwrupr_q <- qchisq(1 - a, df = 1, lower.tail = TRUE)
lwrupr_llf <- r_ml$value - lwrupr_q/2
ベイズ統計学で出てくる最高事後密度信用区間（HPD: highest posterior density）と同じような特性になります。
プロファイル尤度信頼区間の端点にはもう一つ条件があり、信頼区間を計算するパラメーター以外のパラメーターのグラディエントはゼロです。
パラメーターを、対数尤度と信頼区間を計算する変数以外の変数のグラディエントに移すベクトル値関数を考え、これをニュートン・ラフソン法で推定することになります。信頼区間の上限と下限だけに解は二つあるのですが、最尤推定量のグラディエントはゼロで領域が分割されるため、初期値を最尤推定量の左側にするか右側にするかで、上限と下限をそれぞれ推定できます。
lsearch <- function(i, init.p){
r_optim <- optim(rep(0, ncol(X) - 1), pllf, pllfg, i = i, v = init.p, method = "BFGS", control = list(fnscale = -1))
b <- pparam(r_optim$par, i, init.p)
for(j in 1:10){
b0 <- b
f <- llfg(b)
f[i] <- llf(b) - lwrupr_llf
J <- hessian(b)
J[i, ] <- llfg(b)
b <- b - solve(J, f)
if(any(is.na(b))) stop("b has NA!")
if(all(b - b0 < 1e-6)) break
}
b[i]
}
confint.lr <- matrix(NA, length(beta), 2)
for(i in 1:length(beta)){
confint.lr[i, ] <- c(lsearch(i, beta[i] - 2*se[i]), lsearch(i, beta[i] + 2*se[i]))
}
print(confint.lr)
[,1] [,2]
[1,] -0.851403 0.4338161
[2,] 1.263573 2.9297657
[3,] -5.110678 -2.1783299
微妙に違いますが、概ねconfint.glmと同じ値が出てきました。
パラメーターの数×2回のラインサーチがいるわけですが、秒もかかりませんね。なお、ニュートン＝ラフソン法ではなく、最適化関数を複数組み合わせて計算しようとするともっさりとした動きになります。
