I am facing the following problem: I have a regression similar to the following:

fit <- lm(I(hp > 100) ~ cyl + wt + disp + am,data = mtcars)

The idea is more or less to estimate the probability of having hp > 100. I then calculate the mean estimation for each category of am:


        0         1 
0.8421053 0.5384615 

This works for me. Now I would like to have confidence intervals for these estimates. I can do

predict(fit,interval ="confidence" )

         fit         lwr       upr
1  0.7848704  0.58406685 0.9856739
2  0.8131518  0.58796047 1.0383431
3  0.2948876  0.10013362 0.4896416
4  0.5771351  0.38357191 0.7706984
5  0.9762073  0.72139354 1.2310211
6  0.6588185  0.50067156 0.8169655
7  0.9906253  0.76233687 1.2189137
8  0.2156066 -0.04845223 0.4796655
9  0.2209163 -0.04090334 0.4827359
10 0.7514166  0.50829030 0.9945429
11 0.7514166  0.50829030 0.9945429
12 1.1851648  0.92950067 1.4408290
13 1.1474563  0.92047532 1.3744373
14 1.1530017  0.92392416 1.3820792
15 0.9919425  0.66900143 1.3148836

And I get the confidence of each estimation. How can I combine these confidence interval to get the confidence interval of the mean ?


Isabella Ghement pointed out that I should use glm instead of lm for binary variable.

I actually used lm because for the calculation I want to perform, both lm and glm are equivalent (I am not interested in individual estimation, but on the average probability on two subgroups):

model <- glm(I(hp > 100) ~ cyl + wt + disp + am, data = mtcars,
              family = binomial(link="logit"))
tapply(predict(model, type = "response"),mtcars$am,mean)

        0         1 
0.8421053 0.5384615 

But predict(model, type = "response",interval = "confidence") does not give the confidence interval of individual estimation, where it does for lm.

My question is to know if there is a way of calculating the CI of the two values 0.8421053 for am = 0 and 0.5384615 for am = 1 from the individual CI given by predict(fit,interval ="confidence" ), to avoid computer intensive calculation such as bootstraping.

  • $\begingroup$ If you want a 95% CI for the am coefficient (not the "mean"), does confint(fit) answer your question? Otherwise, you may have to define appropriate contrasts using, e.g., the rms or effects packages. $\endgroup$ – chl Nov 19 '20 at 11:45
  • $\begingroup$ Thanks for the tip. I actually want the CI for the calculated mean of the fit column. So my question is: how to I calculate that absed on the lwr and upr columns ? $\endgroup$ – denis Nov 19 '20 at 11:54
  • $\begingroup$ Do you mean you are interested in the 95% CI for tapply(predict(fit),mtcars$am,mean), when, e.g., am = 1? $\endgroup$ – chl Nov 19 '20 at 14:18
  • $\begingroup$ I want to compute the 95% CI for the two values given by tapply(predict(fit),mtcars$am,mean): 0.8421053 for am = 0 and 0.5384615 for am =1. I was wondering if there was a way to do so from the individual CI given by predict, I am trying to avoid bootstraping $\endgroup$ – denis Nov 20 '20 at 9:08
  • 1
    $\begingroup$ Do you mean this? library(rms); d <- datadist(mtcars); options(datadist = "d"); m <- lrm(I(hp > 100) ~ cyl + wt + disp + am, data = mtcars, x = TRUE); Predict(m, am = 0:1). (You can fix the other predictors to whichever values you like, and even perform joint estimation using conf.type="simultaneous".) $\endgroup$ – chl Nov 20 '20 at 11:05

Your response variable is binary - fitting a linear regression model, lm(), to a binary variable is not the best thing you can do. As you can see from your own response, a linear regression model can produce "probabilities" whose values are non-sensical (e.g., > 1). By definition, probabilities can only take values between 0 and 1.

Why not use a more appropriate model for your data, namely a binary logistic regression model? In R, you would fit this model with:

model <- glm(I(hp > 100) ~ cyl + wt + disp + am, data = mtcars,
    family = binomial(link="logit"))

The R command summary(model) will summarize the glm() model fit, showing how the log odds that hp > 100 depend on the values of the predictor variables included in the model.

The glm() model can be used as a basis for predicting the probability that hp > 100 as a function of the cyl, weight, disp and am values for the cars represented in your data:

prob <- predict(model, type = "response") 


Now, let's say you want to focus only on those cars with am = 0 represented in your data; for these cars, you can estimate the "average" probability that they have hp > 100. This involves subsetting prob, as computed above, so that it includes only probabilities for which am = 0, and then computing the average of the subsetted probabilities (across the combinations of cyl, weight and displ values represented in that subset).

You can get a 95% confidence interval for the "average" probability using bootstrapping. Essentially, repeat the subsetting process described above for each of B = 999 (say) bootstrap samples and save your "average" probability for each such sample. Then order the bootstrap "average" probabilities from smallest to largest and report the appropriate quantiles of their distribution as the endpoints of your confidence interval.

  • 1
    $\begingroup$ (+1) Oh, well done, I missed the fact that the response variable was binary! You may want to add a note about the effects package (or other related ones) for computing the marginal effect of am. $\endgroup$ – chl Nov 19 '20 at 21:05
  • $\begingroup$ Thank you, @chl! I’m not sure whether the effects package is helpful here, as it still requires the use of specific values for the other predictors in the model. For this reason, I suggested bootstrapping. $\endgroup$ – Isabella Ghement Nov 20 '20 at 3:02
  • $\begingroup$ Thank you. You are right glm is better, but it does not really change my problem, that was getting the 95% CI. As far as I understand, you are proposing bootstraping. I was wondering if there was a more direct way, i.e. a possible calculation from the individual CI of each estimation. $\endgroup$ – denis Nov 20 '20 at 8:56
  • $\begingroup$ That's exactly my problem: marginal effect would require specific values of the other predictors, and the proportion of the predictors are specific to my population. I was trying to avoid bootstraping $\endgroup$ – denis Nov 20 '20 at 9:03
  • $\begingroup$ @chl I eddited my question to be more specific $\endgroup$ – denis Nov 20 '20 at 9:19

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.