CODEWITHSUNDEEP
r
PPG
PPG

Reputation: 23

why I can't get a confidence interval using predict function in R

I am trying to get a confidence interval for my response in a poisson regression model. Here is my data:

X <- c(1,0,2,0,3,1,0,1,2,0)
Y <- c(16,9,17,12,22,13,8,15,19,11)

What I've done so far: (i) read my data (ii) fit a Y by poisson regression with X as a covariate

model <- glm(Y ~ X, family = "poisson", data = mydata)

(iii) use predict()

predict(model,newdata=data.frame(X=4),se.fit=TRUE,interval="confidence",level=0.95, type = "response")

I was expecting to get "fit, lwr, upr" for my response but I got the following instead:

$fit
       1 
30.21439 

$se.fit
       1 
6.984273 

$residual.scale
[1] 1

Could anyone offer some suggestions? I am new to R and struggling with this problem for a long time. Thank you very much.

Upvotes: 2

Views: 3257

Answers (1)

user11599
user11599

Reputation: 326

First, the function predict() that you are using is the method predict.glm(). If you look at its help file, it does not even have arguments 'interval' or 'level'. It doesn't flag them as erroneous because predict.glm() has the (in)famous ... argument, that absorbs all 'extra' arguments. You can write confidence=34.2 and interval="woohoo" and it still gives the same answer. It only produces the estimate and the standard error.

Second, one COULD then take the fit +/- 2*se to get an approximate 95 percent confidence interval. However, without getting into the weeds of confidence intervals, pivotal statistics, non-normality in the response scale, etc., this doesn't give very satisfying intervals because, for instance, they often include impossible negative values.

So, I think a better approach is to form an interval in the link scale, then transform it (this is still an approximation, but probably better):

X <- c(1,0,2,0,3,1,0,1,2,0)
Y <- c(16,9,17,12,22,13,8,15,19,11)
model <- glm(Y ~ X, family = "poisson")
tmp <- predict(model, newdata=data.frame(X=4),se.fit=TRUE, type = "link")
exp(tmp$fit - 2*tmp$se.fit)
       1 
19.02976
exp(tmp$fit + 2*tmp$se.fit)
   1 
47.97273

Upvotes: 2

Related Questions