20
$\begingroup$

Ok, I have a logistic regression and have used the predict() function to develop a probability curve based on my estimates.

## LOGIT MODEL:
library(car)
mod1 = glm(factor(won) ~ as.numeric(bid), data=mydat, family=binomial(link="logit"))

## PROBABILITY CURVE:
all.x <- expand.grid(won=unique(won), bid=unique(bid))
y.hat.new <- predict(mod1, newdata=all.x, type="response")
plot(bid<-000:1000,predict(mod1,newdata=data.frame(bid<-c(000:1000)),type="response"), lwd=5, col="blue", type="l")

This is great but I'm curious about plotting the confidence intervals for the probabilities. I've tried plot.ci() but had no luck. Can anyone point me to some ways to get this done, preferably with the car package or base R.

$\endgroup$
  • 4
    $\begingroup$ (+1) In response to the votes to close as off topic: Apparently the basis for those votes is that the question appears to ask a purely software-related question ("how to plot such-and-such in R"), a question that indeed ought to appear on SO. Note, however, that buried in the current reply are statistical formulas to create the plotting points. This suggests there is statistical interest to the question, so I am reluctant to vote for migration. A good reply here would highlight and explain this statistical point. $\endgroup$ – whuber May 24 '12 at 14:13
26
$\begingroup$

The code you used estimates a logistic regression model using the glm function. You didn't include data, so I'll just make some up.

set.seed(1234)
mydat <- data.frame(
    won=as.factor(sample(c(0, 1), 250, replace=TRUE)), 
    bid=runif(250, min=0, max=1000)
)
mod1 <- glm(won~bid, data=mydat, family=binomial(link="logit"))

A logistic regression model models the relationship between a binary response variable and, in this case, one continuous predictor. The result is a logit-transformed probability as a linear relation to the predictor. In your case, the outcome is a binary response corresponding to winning or not winning at gambling and it is being predicted by the value of the wager. The coefficients from mod1 are given in logged odds (which are difficult to interpret), according to:

$$\text{logit}(p)=\log\left(\frac{p}{(1-p)}\right)=\beta_{0}+\beta_{1}x_{1}$$

To convert logged odds to probabilities, we can translate the above to

$$p=\frac{\exp(\beta_{0}+\beta_{1}x_{1})}{(1+\exp(\beta_{0}+\beta_{1}x_{1}))}$$

You can use this information to set up the plot. First, you need a range of the predictor variable:

plotdat <- data.frame(bid=(0:1000))

Then using predict, you can obtain predictions based on your model

preddat <- predict(mod1, newdata=plotdat, se.fit=TRUE)

Note that the fitted values can also be obtained via

mod1$fitted

By specifying se.fit=TRUE, you also get the standard error associated with each fitted value. The resulting data.frame is a matrix with the following components: the fitted predictions (fit), the estimated standard errors (se.fit), and a scalar giving the square root of the dispersion used to compute the standard errors (residual.scale). In the case of a binomial logit, the value will be 1 (which you can see by entering preddat$residual.scale in R). If you want to see an example of what you've calculated so far, you can type head(data.frame(preddat)).

The next step is to set up the plot. I like to set up a blank plotting area with the parameters first:

with(mydat, plot(bid, won, type="n", 
    ylim=c(0, 1), ylab="Probability of winning", xlab="Bid"))

Now you can see where it is important to know how to calculate the fitted probabilities. You can draw the line corresponding to the fitted probabilities following the second formula above. Using the preddat data.frame you can convert the fitted values to probabilities and use that to plot a line against the values of your predictor variable.

with(preddat, lines(0:1000, exp(fit)/(1+exp(fit)), col="blue"))

Finally, answer your question, the confidence intervals can be added to the plot by calculating the probability for the fitted values +/- 1.96 times the standard error:

with(preddat, lines(0:1000, exp(fit+1.96*se.fit)/(1+exp(fit+1.96*se.fit)), lty=2))
with(preddat, lines(0:1000, exp(fit-1.96*se.fit)/(1+exp(fit-1.96*se.fit)), lty=2))

The resulting plot (from the randomly generated data) should look something like this:

enter image description here

For expediency's sake, here's all the code in one chunk:

set.seed(1234)
mydat <- data.frame(
    won=as.factor(sample(c(0, 1), 250, replace=TRUE)), 
    bid=runif(250, min=0, max=1000)
)
mod1 <- glm(won~bid, data=mydat, family=binomial(link="logit"))
plotdat <- data.frame(bid=(0:1000))
preddat <- predict(mod1, newdata=plotdat, se.fit=TRUE)
with(mydat, plot(bid, won, type="n", 
    ylim=c(0, 1), ylab="Probability of winning", xlab="Bid"))
with(preddat, lines(0:1000, exp(fit)/(1+exp(fit)), col="blue"))
with(preddat, lines(0:1000, exp(fit+1.96*se.fit)/(1+exp(fit+1.96*se.fit)), lty=2))
with(preddat, lines(0:1000, exp(fit-1.96*se.fit)/(1+exp(fit-1.96*se.fit)), lty=2))

(Note: This is a heavily edited answer in an attempt to make it more relevant to stats.stackexchange.)

$\endgroup$
  • $\begingroup$ where is the variable se.fit defined? $\endgroup$ – Macro May 24 '12 at 16:43
  • $\begingroup$ In predict(..., se.fit=TRUE). $\endgroup$ – smillig May 24 '12 at 17:09
  • $\begingroup$ (-1) These CIs are for each for individual cases? If so, for a binary outcome, the only sensible CI for a predicted probability is [0,1]. Even though this may be a technically proficient answer. $\endgroup$ – rolando2 Mar 23 '13 at 17:19
  • $\begingroup$ Per @whuber's comment, I think a good answer should include a formula for how the SE is calculated. Could someone perhaps edit and improve the answer? $\endgroup$ – Heisenberg Apr 5 '15 at 5:57
  • 1
    $\begingroup$ Your answer seems to give just the 'mean prediction interval'. How would I add the 'point prediction interval'? $\endgroup$ – Bob Hopez Aug 5 '16 at 17:16
0
$\begingroup$

Here's a modification of @smillig's solution. I use tidyverse tools here, and also use the linkinv function that is a part of the GLM model object mod1. That way, you don't have to manually invert the logistic function, and this approach will work regardless of what specific GLM you fit.

library(tidyverse)
library(magrittr)


set.seed(1234)

# create fake data on gambling. Does prob win depend on bid size? 
mydat <- data.frame(
  won=as.factor(sample(c(0, 1), 250, replace=TRUE)), 
  bid=runif(250, min=0, max=1000)
)

# logistic regression model: 
mod1 <- glm(won~bid, data=mydat, family=binomial(link="logit"))

# new predictor values to use for prediction: 
plotdat <- data.frame(bid=(0:1000))

# df with predictions, lower and upper limits of CIs: 
preddat <- predict(mod1,
               type = "link",
               newdata=plotdat,
               se.fit=TRUE) %>% 
  as.data.frame() %>% 
  mutate(bid = (0:1000), 

         # model object mod1 has a component called linkinv that 
         # is a function that inverts the link function of the GLM:
         lower = mod1$family$linkinv(fit - 1.96*se.fit), 
         point.estimate = mod1$family$linkinv(fit), 
         upper = mod1$family$linkinv(fit + 1.96*se.fit)) 


# plotting with ggplot: 
preddat %>% ggplot(aes(x = bid, 
                   y = point.estimate)) + 
  geom_line(colour = "blue") + 
  geom_ribbon(aes(ymin = lower,
                  ymax = upper), 
              alpha = 0.5) + 
  scale_y_continuous(limits = c(0,1))
$\endgroup$
  • 3
    $\begingroup$ Although implementation is often mixed with substantive content in questions, we are supposed to be a site for providing information about statistics, machine learning, etc., not code. It can be good to provide code as well, but please elaborate your substantive answer in text for people who don't read this language well enough to recognize & extract the answer from the code. $\endgroup$ – gung - Reinstate Monica Oct 5 '18 at 1:56

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.