# Logistic regression, “sum” confidence interval

I have an existing fitted logit regression model.

Model:

$\hat{p}(x)=\frac{1}{1+e^{-\hat{\beta}x}}$

With parameter estimates $\hat{\beta}$, and observation $x$

Given a new set of datapoints $x_1,\ldots,x_n$ and $o_1,\ldots,o_n$, I want to construct a confidence interval for $E=\sum_i \hat{p}(x_i)$. To goal is to check if $O=\sum_i o_i$ is (not) significantly different.

For a single observation, you can construct a confidence interval for the probability quite easily. Given $\hat{\beta}$ is normally distributed, and the confidence interval for these estimates is ($\hat{\beta}_{0.025}, \hat{\beta}_{0.975}$), the confidence interval for the predictor is: ($\frac{1}{1 + e^{-\hat{\beta}_{0.025}x}}$, $\frac{1}{1 + e^{-\hat{\beta}_{0.975}x}}$).

Is there a good way to get the confidence interval for the "sum" of the predicted probabilities?

• Ummm, $\hat{\beta}$ are not variables, they are parameter estimates. And $x$ is usually regarded as a matrix of independent variables (each observation includes both the DV and these IVs). – Peter Flom - Reinstate Monica Jun 10 '13 at 11:28
• You're right. I think this edit makes things clear enough. In my case, $x$ is relatively simple, just the intercept and some observations as a vector. – Tim Jun 10 '13 at 11:53
• It looks like you should instead be asking for a prediction interval for $O$, assuming your $o_i$ are the (binary) observations corresponding to the $x_i$. ($E$ itself has little direct relevance to your question; in particular, a CI for $E$ does not translate into an interval for $O$.) – whuber Jun 10 '13 at 14:31
• Yes, the $o_i$ are binary observations. But given that you know the distribution of $\hat{p}$, and because the model parameters $\hat{B}$ are estimated from a much larger sample (in this case), wouldn't it make more sense to estimate an interval around E? Why doesn't it translate? I am tempted to just do a bootstrap interval and be done with it. – Tim Jun 10 '13 at 14:59

One possible solution is to use a boot strapping approach, given the new set of data points, to construct a boot strap estimate and confidence interval for $\sum_i\hat p(x_i)$

set.seed(44)

#Pseudo Data
prob = .2
x = sort(rnorm(100))
y = rbinom(100,1,prob)

#Fit the logistic regression model
model = glm(y~x,family="binomial")

#More data arrives
x.new = rnorm(10000)

#The number of bootstrap samples
B = 10000

sum.p = rep(NA,)
for(i in 1:B){
#Create a bootstrap sample
x.boot = sample(x.new,length(x.new)*.1,replace=TRUE)

#Calculate the sum of p
sum.p[i] = sum(1/(1+exp(-(model$coef[1]+model$coef[2]*x.boot))))
}

#Get the 2.5% and 97.5% quantile from the bootstrap estimator
lower = quantile(sum.p,prob=.025)
upper = quantile(sum.p,prob=.975)

#Construct a 95% confidence interval
ci = c(lower,upper)


An alternative method would be to take a Bayesian approach and recaculate $\sum_i\hat p(x_i)$ for every sample of $\beta$ at every step of an MCMC type algorithm. Then at the end of our MCMC we would have a sample of $\sum_i\hat p(x_i)^{(j)}$ for every $j^{th}$ step of the MCMC that we could take the quantiles of in order to obtain our 95% confidence interval for $\sum_i\hat p(x_i)$ .

# Using Scortchi's Suggestion here is the revised code:

#Scortchi's suggestion
set.seed(44)
prob = .2

#More data arrives
x.new = rnorm(10000)
y.new = rbinom(10000,1,prob)

#The number of bootstrap samples
B = 10000

sum.p = rep(NA,B)
for(i in 1:B){
#Create a bootstrap sample
index = sample(1:length(x.new),length(x.new)*.1,replace=TRUE)
x.boot = x.new[index]
y.boot = y.new[index]

model = glm(y.boot~x.boot,family="binomial")

#Calculate the sum of p
sum.p[i] = sum(1/(1+exp(-(model$coef[1]+model$coef[2]*x.boot))))
}

#Get the 2.5% and 97.5% quantile from the bootstrap estimator
lower = quantile(sum.p,prob=.025)
upper = quantile(sum.p,prob=.975)

#Construct a 95% confidence interval
ci = c(lower,upper)


Now interestingly, the confidence interval from using Scortchi's suggestion results in

> ci
2.5% 97.5%
174   223


where as using my original code we obtain the following:

> ci
2.5%    97.5%
242.4727 247.0230


So there is clearly a difference between the two methods.

• Looks promising. While I'm absorbing this, please tell me whether there's a missing "-" in the line below "#Calculate the sum of p". – rolando2 Aug 31 '13 at 2:29
• Good catch. That would make my confidence interval look more reasonable too :) – user25658 Aug 31 '13 at 2:32
• @BabakP, correct me if I am wrong. Shouldn't the x.new be drawn from the same normal distribution rather than a binomial distribution? – Jerry Aug 31 '13 at 19:59
• This seems to be answering the question "What's the distribution of the sum of predicted probabilities given a randomly-drawn predictor vector & known coefficients?". I think the OP wanted to know how to transform uncertainty in estimated coefficients to uncertainty in the sum of predicted probabilities given a known predictor vector. Though @whuber's surely right in saying that what the OP ought to want to know is a prediction interval. – Scortchi - Reinstate Monica Sep 1 '13 at 20:22
• @Scortchi, I guess I do not get that interpretation when I read the question since they state that they have parameter estimates $\hat \beta$ and $x$. Never the less, if you want to change it to what you are saying then just make the simple adjustment of re-estimating the model parameters each time inside of the for-loop with each bootstrap sample. Not sure why we need to fuss with this anyway though. – user25658 Sep 1 '13 at 22:43