泊松回归预测
Prediction of poisson regression
我目前正在使用模型处理数据集
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))
亲切的问候。
编辑：关于置信区间的问题
我还有一个问题。我创建了这样的置信区间：
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
但是否可以将其绘制在同一张图中？
使用predict时，需要设置newdata。简单地调用 predict 而不调用 newdata 只会 return 拟合值。因此，您的 predict 呼叫实质上是让您 glm1$fitted.values.
看，您想要从 GRP == "CONTROL" 和 FU == 12 预测 SEX == "FEMALE"。使用
## 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.
中没有匹配的列与选择标准
Follow-up（其实比上面更值得回答）
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。线性预测变量是渐近正态的，因此您需要为线性预测变量生成置信带，然后使用反 link 函数将其转换为响应尺度。
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))
我目前正在使用模型处理数据集
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))
亲切的问候。
编辑：关于置信区间的问题
我还有一个问题。我创建了这样的置信区间：
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
但是否可以将其绘制在同一张图中？
使用predict时，需要设置newdata。简单地调用 predict 而不调用 newdata 只会 return 拟合值。因此，您的 predict 呼叫实质上是让您 glm1$fitted.values.
看，您想要从 GRP == "CONTROL" 和 FU == 12 预测 SEX == "FEMALE"。使用
## 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.
Follow-up（其实比上面更值得回答）
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。线性预测变量是渐近正态的，因此您需要为线性预测变量生成置信带，然后使用反 link 函数将其转换为响应尺度。
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))