4
$\begingroup$

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?

$\endgroup$
  • 2
    $\begingroup$ 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). $\endgroup$ – Peter Flom - Reinstate Monica Jun 10 '13 at 11:28
  • $\begingroup$ 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. $\endgroup$ – Tim Jun 10 '13 at 11:53
  • 1
    $\begingroup$ 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$.) $\endgroup$ – whuber Jun 10 '13 at 14:31
  • $\begingroup$ 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. $\endgroup$ – Tim Jun 10 '13 at 14:59
4
+200
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Looks promising. While I'm absorbing this, please tell me whether there's a missing "-" in the line below "#Calculate the sum of p". $\endgroup$ – rolando2 Aug 31 '13 at 2:29
  • $\begingroup$ Good catch. That would make my confidence interval look more reasonable too :) $\endgroup$ – user25658 Aug 31 '13 at 2:32
  • $\begingroup$ @BabakP, correct me if I am wrong. Shouldn't the x.new be drawn from the same normal distribution rather than a binomial distribution? $\endgroup$ – Jerry Aug 31 '13 at 19:59
  • 1
    $\begingroup$ 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. $\endgroup$ – Scortchi - Reinstate Monica Sep 1 '13 at 20:22
  • 1
    $\begingroup$ @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. $\endgroup$ – user25658 Sep 1 '13 at 22:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.