Prediction of poisson regression

Question

I'm currently working on a data set with the model

glm1 <- glm(FALL ~ GRP + AGE + SEX + offset(log(FU)), family=poisson, data=dat)

Now I need to make a prediction of the amount of falls in one year for a female who's in the control group.

I need to do the predict function, but I'm not sure how. I tried to do several things and last tried this:

levels(dat$GRP)
levels(dat$SEX)
SEX="FEMALE"
GRP="CONTROL"
FU="12"
y<- predict(glm1, type = 'response')
plot(x=dat$AGE[order(dat$AGE)],y=y[order(dat$FALL)],type='l')

But this gives me only a weird looking plot. What do I need to do?

---

Edit: data added on request for reproducibility

dat <- structure(list(FALL = c(0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 2L, 1L, 
2L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 
3L, 0L, 1L, 1L, 0L, 0L, 2L, 3L, 0L, 0L, 3L, 1L, 0L, 0L, 2L, 1L, 
2L, 2L, 1L, 1L, 0L, 0L, 0L, 4L, 1L, 0L, 0L, 0L, 0L, 2L, 3L, 1L, 
0L, 1L, 2L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 
3L, 4L, 0L, 1L, 0L, 0L, 1L, 1L, 2L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 
1L, 0L, 1L, 0L, 0L, 3L, 0L, 0L, 2L, 0L, 0L, 2L, 0L, 3L, 1L, 0L, 
0L, 1L, 1L, 2L, 1L, 0L, 0L, 0L, 0L, 1L, 0L), GRP = structure(c(1L, 
2L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 
2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 
2L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
1L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 
1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 
2L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 
2L, 2L, 2L, 1L), .Label = c("CONTROL", "TAI CHI"), class = "factor"), 
FU = c(18, 12, 17, 4, 23, 16, 22, 24, 23, 11, 22, 9, 23, 
8, 20, 17, 23, 17, 15, 17, 19, 21, 22, 16, 14, 21, 20, 21, 
7, 22, 19, 12, 15, 21, 24, 11, 23, 21, 10, 15, 19, 19, 16, 
24, 17, 23, 16, 17, 18, 18, 20, 8, 21, 16, 15, 19, 23, 14, 
13, 6, 16, 18, 9, 7, 16, 14, 16, 18, 13, 12, 15, 22, 17, 
17, 20, 21, 11, 24, 9, 13, 24, 12, 21, 20, 19, 17, 21, 15, 
17, 11, 24, 10, 18, 9, 16, 19, 6, 13, 22, 18, 10, 15, 14, 
21, 21, 5, 24, 21, 11, 23, 21, 16, 22, 6, 24, 18, 21), AGE = c(71, 
81, 71, 79, 77, 79, 76, 86, 75, 75, 76, 83, 71, 80, 77, 79, 
77, 74, 83, 81, 83, 79, 74, 79, 78, 85, 82, 71, 81, 78, 82, 
74, 73, 75, 83, 78, 83, 83, 65, 75, 75, 75, 75, 78, 80, 69, 
80, 73, 74, 79, 76, 78, 70, 77, 77, 76, 84, 71, 73, 76, 80, 
77, 74, 78, 68, 76, 77, 76, 72, 72, 76, 82, 72, 80, 78, 83, 
80, 73, 79, 75, 79, 75, 80, 77, 81, 78, 74, 79, 78, 74, 79, 
77, 77, 85, 79, 73, 78, 73, 70, 68, 74, 82, 75, 77, 77, 73, 
73, 83, 74, 87, 76, 81, 77, 78, 66, 79, 82), SEX = structure(c(1L, 
1L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 
2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 
2L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 
1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 
1L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 2L, 1L), .Label = c("FEMALE", 
"MALE"), class = "factor")), .Names = c("FALL", "GRP", "FU", 
"AGE", "SEX"), class = "data.frame", row.names = c(NA, -117L))

Kind regards.

---

Edit: question on confidence interval

I have one more question. I created the confidence intervals like this:

prs <- predict(glm1, newdata = newdat, type = "response", se.fit=TRUE)
newdat$pred <- prs[[1]]
newdat$se <- prs[[2]]
newdat$lo <- newdat$pred - 1.96 * newdat$se 
newdat$up <- newdat$pred + 1.96 * newdat$se

But is it possible to plot this in the same graph?

Answer 1

When you use predict, you need to set newdata. Simply calling predict without newdata will just return fitted values. So your predict call is essentially getting you glm1$fitted.values.

Look, you want prediction for SEX == "FEMALE" from GRP == "CONTROL" with FU == 12. Use

## I use `AGE = 65:87` because this is what `range(dat$AGE)` gives
## we must provide all covariates used in model formula to make `predict` work
## recycling rule is applied here.
## `GRP`, `SEX` and `FU` are given a single value, while `AGE` has length 23
## they will be recycled 23 times
newdat <- data.frame(AGE = 65:87, GRP = "CONTROL", SEX = "FEMALE", FU = 12)
pred <- predict(glm1, newdata = newdat, type = "response")
plot(newdat$AGE, pred, type = "l")

Initially I suggested:

newdat <- subset(dat, GRP == "CONTROL" & SEX == "FEMALE" & FU == 12)

but this is a bad idea. It will give you an empty data frame, since there are no matching columns with selection criteria in your dat.

---

Follow-up (actually more worth answering than above)

I have one more question. I created the confidence intervals like this:

prs <- predict(glm1, newdata = newdat, type = "response", se.fit=TRUE)
newdat$pred <- prs[[1]]
newdat$se <- prs[[2]]
newdat$lo <- newdat$pred - 1.96 * newdat$se 
newdat$up <- newdat$pred + 1.96 * newdat$se

But is it possible to plot this in the same graph?

Your confidence interval is not correctly computed. Response is not normally distributed, so you can't use 1.96. Linear predictor is asymptotically normal, so you need produce confidence band for linear predictor, then transform it to response scale using inverse link function.

ginv <- glm1$family$linkinv  ## inverse link function
prs <- predict(glm1, newdata = newdat, type = "link", se.fit=TRUE)
newdat$pred <- ginv(prs[[1]])
newdat$lo <- ginv(prs[[1]] - 1.96 * prs[[2]])
newdat$up <- ginv(prs[[1]] + 1.96 * prs[[2]])

To plot them on the same plot, you can use plot + lines:

with(newdat, plot(AGE, pred, type = "l", ylim = c(min(lo), max(up)) ))
with(newdat, lines(AGE, lo, lty = 2))
with(newdat, lines(AGE, up, lty = 2))

Or, you may use matplot:

matplot(newdat[c("pred", "lo", "up")], type = "l", col = 1, lty = c(1, 2, 2))