Прогнозирование пуассоновского регрессии

Question

В настоящее время я работаю над набором данных с моделью

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

Теперь мне нужно, чтобы сделать прогноз количества падений в один год для женщины, кто находится в контрольной группе.

Мне нужно сделать функцию predict, но я не уверен, как это сделать. Я пытался сделать несколько вещей и в последний раз пробовал это:

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') 

Но это дает мне только странный перспективный сюжет. Что мне нужно сделать?

---

Редактировать данные: По запросу воспроизводимости

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)) 

С наилучшими пожеланиями.

---

Edit: вопрос о доверительном интервале

У меня есть еще один вопрос. Я создал доверительные интервалы следующим образом:

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 

Но возможно ли построить на этом же графике?

Answer 1

Когда вы используете predict, вам необходимо установить newdata. Просто позвонив predict без newdata, просто верните установленные значения. Таким образом, ваш звонок predict по существу доставит вам glm1$fitted.values.

Посмотрите, вы хотите прогноз на SEX == "FEMALE" от GRP == "CONTROL" с FU == 12. Используйте

## 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") 

Первоначально я предложил:

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

, но это плохая идея. Это даст вам пустой фрейм данных, поскольку в вашем dat нет соответствующих столбцов с критериями выбора.

---

Последующая деятельность (на самом деле больше, чем стоит отвечать выше)

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?

Ваш доверительный интервал не правильно вычислить. Ответ обычно не распространяется, поэтому вы не можете использовать 1.96. Линейный предиктор является асимптотически нормальным, поэтому вам нужно создать доверительный диапазон для линейного предиктора, а затем преобразовать его в шкалу ответа с использованием функции обратной ссылки.

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]]) 

Чтобы построить их на том же участке, вы можете использовать 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)) 

Или, вы можете использовать matplot:

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