R glm(): confusing confidence intervals for binomial regression
<p>I have been using glm() in R to compute confidence intervals for the logit probability parameter governing a single binomial draw.  </p>

<pre><code>P &lt;- 20 # Number of successes
D &lt;- 1  # Number of failures
model1 &lt;- 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
</code></pre>

<p>I notice that if either the number of successes or failures is zero (<code>P=0</code> or <code>D=0</code>), the confidence intervals returned are meaningless.</p>

<p>I then computed my own confidence intervals by numerically integrating the normalized binomial likelihood:</p>

<pre><code>binomial_fun &lt;- function(p) {choose(N, K)*(p^K)*(1-p)^(N-K)}  # A binomial function
logit_fun &lt;- function(p) {log(p/(1-p))}                       # A logit function
f_upper &lt;- function(a) {integrate(binomial_fun, 0, a )$value/integrate(binomial_fun, 0, 1 )$value - .975}
f_lower &lt;- 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 &lt;- P+D     # N is the number of trials
K &lt;- P       # K is the number of successes
uci2 &lt;- logit_fun(uniroot( f_upper, c(0, 1) )$root) # Look for a solution in [0,1]
lci2 &lt;- logit_fun(uniroot( f_lower, c(0, 1) )$root) # Look for a solution in [0,1]
</code></pre>

<p>This gave meaningful confidence intervals where <code>P=0</code> or <code>D=0</code>.
However, it gives different confidence intervals from <code>glm()</code> + <code>confint()</code> even when neither <code>P</code> nor <code>D</code> is zero.</p>

<pre><code>confint(model1)
c(lci2, uci2)
</code></pre>

<p>Compared to <code>glm()</code> + <code>confint()</code>, both the lci and uci that I compute tend to be closer to zero.</p>

<p>How is <code>confint()</code> computing the interval?  I'm sure there is a good reason for doing it this way in terms of performance on more complicated glms, but it seems like a weird result in this simple case.</p>