1

我一直在 R 中使用 glm() 来计算控制单个二项式绘制的 logit 概率参数的置信区间。 

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。但是，即使两者都不为零，它也会给出与glm()+不同的置信区间。confint()PD

confint(model1)
c(lci2, uci2)

与glm()+相比confint()，我计算的 lci 和 uci 往往更接近于零。

如何confint()计算间隔？我敢肯定，就更复杂的 glms 的性能而言，这样做是有充分理由的，但在这种简单的情况下，这似乎是一个奇怪的结果。

4

0 回答 0