私は R で glm() を使用して、単一の二項描画を管理するロジット確率パラメーターの信頼区間を計算しています。
P <- 20 # Number of successes
D <- 1 # Number of failures
model1 <- glm(matrix(c(P,D), nrow=1) ~ 1, family="binomial") # Successes modeled as binomial draw from successes+failures
summary(model1)
confint(model1) # Confidence interval on binomial parameter p
成功または失敗の回数がゼロ (P=0またはD=0) の場合、返される信頼区間は無意味であることに気付きました。
次に、正規化された二項尤度を数値的に統合して、独自の信頼区間を計算しました。
binomial_fun <- function(p) {choose(N, K)*(p^K)*(1-p)^(N-K)} # A binomial function
logit_fun <- function(p) {log(p/(1-p))} # A logit function
f_upper <- function(a) {integrate(binomial_fun, 0, a )$value/integrate(binomial_fun, 0, 1 )$value - .975}
f_lower <- function(a) {integrate(binomial_fun, a, 1 )$value/integrate(binomial_fun, 0, 1 )$value - .975}
# These functions will take value zero when a corresponds to the 97.5% and 2.5% points
# of the normalized binomial likelihood given N and K.
N <- P+D # N is the number of trials
K <- P # K is the number of successes
uci2 <- logit_fun(uniroot( f_upper, c(0, 1) )$root) # Look for a solution in [0,1]
lci2 <- logit_fun(uniroot( f_lower, c(0, 1) )$root) # Look for a solution in [0,1]
P=0これにより、またはの意味のある信頼区間が得られましD=0た。ただし、 norがゼロでない場合でも、glm()+とは異なる信頼区間が得られます。confint()PD
confint(model1)
c(lci2, uci2)
glm()+と比較するとconfint()、私が計算した lci と uci はどちらもゼロに近い傾向があります。
間隔をどのようにconfint()計算していますか? より複雑な glms でのパフォーマンスに関しては、このようにする正当な理由があると確信していますが、この単純なケースでは奇妙な結果のように思えます。