Search code examples
rpredictionglmconfidence-interval

Confidence interval of prediction from glm model


I am trying to obtain the 95% confidence interval of the prediction from a glm. Below is my model

library(caret)
library(ISLR)
data(Smarket)
glm.fit <- glm(Direction ~ Lag1 + Lag2, data = Smarket, family = binomial)

Now I tried with below code to estimate the 95% confidence interval for prediction

> head(predict(glm.fit, interval="confidence", level = 0.025))
          1           2           3           4           5           6 
 0.05554775 -0.01128128 -0.04222017  0.07288089  0.05806364  0.03169772 
> head(predict(glm.fit, interval="confidence", level = 0.975))
          1           2           3           4           5           6 
 0.05554775 -0.01128128 -0.04222017  0.07288089  0.05806364  0.03169772 

Clearly this is not giving me right interval as in the both cases, the values are same.

Could you please help me to get the right interval?


Solution

  • predict.glm() doesn't take the same arguments as predict.lm (see ?predict.glm): you have to do this by hand (or find a package with helper functions). The following code constructs the lower and upper 95% Wald confidence limits on the logit (log-odds) scale and then uses plogis() to back-transform to the probability scale ...

    pp <- predict(glm.fit, se.fit = TRUE)
    ci_lwr <- with(pp, plogis(fit + qnorm(0.025)*se.fit))
    ci_upr <- with(pp, plogis(fit + qnorm(0.975)*se.fit))
    
    > head(ci_lwr)
            1         2         3         4         5         6 
    0.4842931 0.4596593 0.4451171 0.4780052 0.4796479 0.4759596 
    > head(ci_upr)
            1         2         3         4         5         6 
    0.5433766 0.5347319 0.5339426 0.5581846 0.5492351 0.5398233