気づかず使っているかも知れないプロファイル尤度信頼区間の計算

注意深いRユーザーはお気づきだと思いますが、glmでロジットモデルを推定したあと結果オブジェクトにconfintをかけると、waiting for profiling to be done...とメッセージが表示されます。

r_glm <- glm(ans ~ x + z, data = df02, family = binomial()) # df02の生成は後で説明
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.1783024

profiling? — 意味がわかりませんよね。?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)
    # X <- cbind(1, mvtnorm::rmvnorm(n, sigma = S))
    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)
# Wald型信頼区間（左右対象）
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){
    # パラメーターのi番目をvとし、i番目以外はcpの値をcpの並びで使う
    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){ # 最大ループ回数10は一般には小さい
        b0 <- b
        # グラディエントを一行変えて、ベクトル値関数の値にする
        f <- llfg(b)
        f[i] <- llf(b) - lwrupr_llf
        # ヘッシアンを一行変えて、ベクトル値関数のヤコビアンにする
        J <- hessian(b)
        J[i, ] <- llfg(b)
        b <- b - solve(J, f) # = 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回のラインサーチがいるわけですが、秒もかかりませんね。なお、ニュートン＝ラフソン法ではなく、最適化関数を複数組み合わせて計算しようとするともっさりとした動きになります。