创建和绘制置信区间
Creating and plotting confidence intervals
我已经为我的数据拟合了一个高斯 GLM 模型，我现在希望创建 95% CIs 并将它们拟合到我的数据中。我在绘图时遇到了几个问题，因为我无法让它们捕获我的数据，它们似乎只是绘制了与模型相同的线而没有捕获数据点。此外，我也不确定我是否以正确的方式在此处创建了 CIs 的均值。如果有人知道如何解决这个问题，我在下面输入了我的数据和代码
使用的数据
aids
cases quarter date
1 2 1 83.00
2 6 2 83.25
3 10 3 83.50
4 8 4 83.75
5 12 1 84.00
6 9 2 84.25
7 28 3 84.50
8 28 4 84.75
9 36 1 85.00
10 32 2 85.25
11 46 3 85.50
12 47 4 85.75
13 50 1 86.00
14 61 2 86.25
15 99 3 86.50
16 95 4 86.75
17 150 1 87.00
18 143 2 87.25
19 197 3 87.50
20 159 4 87.75
21 204 1 88.00
22 168 2 88.25
23 196 3 88.50
24 194 4 88.75
25 210 1 89.00
26 180 2 89.25
27 277 3 89.50
28 181 4 89.75
29 327 1 90.00
30 276 2 90.25
31 365 3 90.50
32 300 4 90.75
33 356 1 91.00
34 304 2 91.25
35 307 3 91.50
36 386 4 91.75
37 331 1 92.00
38 368 2 92.25
39 416 3 92.50
40 374 4 92.75
41 412 1 93.00
42 358 2 93.25
43 416 3 93.50
44 414 4 93.75
45 496 1 94.00
我的代码用于在绘图之前创建模型和间隔
#creating the model
model3 = glm(cases ~ date,
data = aids,
family = poisson(link='log'))
#now to add approx. 95% confidence envelope around this line
#predict again but at the linear predictor level along with standard errors
my_preds <- predict(model3, newdata=data.frame(aids), se.fit=T, type="link")
#calculate CI limit since linear predictor is approx. Gaussian
upper <- my_preds$fit+1.96*my_preds$se.fit #this might be logit not log
lower <- my_preds$fit-1.96*my_preds$se.fit
#transform the CI limit to get one at the level of the mean
upper <- exp(upper)/(1+exp(upper))
lower <- exp(lower)/(1+exp(lower))
#plotting data
plot(aids$date, aids$cases,
xlab = 'Date', ylab = 'Cases', pch = 20)
#adding CI lines
plot(aids$date, exp(my_preds$fit), type = "link",
xlab = 'Date', ylab = 'Cases') #add title
lines(aids$date,exp(my_preds$fit+1.96*my_preds$se.fit),lwd=2,lty=2)
lines(aids$date,exp(my_preds$fit-1.96*my_preds$se.fit),lwd=2,lty=2)
我目前没有数据点的结果，这里的模型是正确的，但是 CI 不是因为我没有数据点，所以 CI 是错误的，我认为在某个地方
编辑：对OP提供完整数据集的回应。
这最初是关于在同一图表上绘制数据和模型的问题，但已经发生了很大变化。您似乎对原始问题有答案。下面是解决其余问题的一种方法。
看看你的（和我的）图，泊松 glm 显然不是一个好的模型。换句话说，案例数量可能会随日期而变化，但也会受到模型中不存在的其他因素（外部回归变量）的影响。
仅绘制您的数据强烈表明您至少有两种甚至更多的制度：病例增长遵循不同模型的时间框架。
ggplot(aids, aes(x=date)) + geom_point(aes(y=cases))
这表明 segmented regression。与 R 中的大多数东西一样，有一个包（实际上不止一个）。下面的代码使用 segmented 包使用 1 个断点（两个区域）构建连续的泊松 glm。
library(data.table)
library(ggplot2)
library(segmented)
setDT(aids) # convert aids to a data.table
aids[, pred:=
predict(
segmented(glm(cases~date, .SD, family = poisson), seg.Z = ~date, npsi=1),
type='response', se.fit=TRUE)$fit]
ggplot(aids, aes(x=date))+ geom_line(aes(y=pred))+ geom_point(aes(y=cases))
请注意，我们需要告诉 segmented 断点的数量，而不是断点的位置 - 算法会为您计算出来。所以在这里，我们看到 3Q87 之前的状态使用泊松 glm 很好地建模，而之后的状态则不是。这是一种奇特的说法，即 87 年 3 季度左右“发生了一些事情”，它改变了疾病的进程（至少在这个数据中是这样）。
下面的代码做同样的事情，但有 1 到 4 个断点。
get.pred <- \(p.n, p.DT) {
fit <- glm(cases~date, p.DT, family=poisson)
seg.fit <- segmented(fit, seg.Z = ~date, npsi=p.n)
predict(seg.fit, type='response', se.fit=TRUE)[c('fit', 'se.fit')]
}
gg.dt <- rbindlist(lapply(1:4, \(x) { copy(aids)[, c('pred', 'se'):=get.pred(x, .SD)][, npsi:=x] } ))
ggplot(gg.dt, aes(x=date))+
geom_ribbon(aes(ymin=pred-1.96*se, ymax=pred+1.96*se), fill='grey80')+
geom_line(aes(y=pred))+
geom_point(aes(y=cases))+
facet_wrap(~npsi)
请注意，第一个断点的位置似乎没有改变，而且，尽管使用了泊松 glm，但除第一个区域外，所有区域的增长都是线性的。
包文档中描述了 goodness-of-fit 指标，可以帮助您确定有多少断点与您的数据最一致。
最后，还有 mcp 包，它更强大一些，但使用起来也更复杂一些。
原始响应：这是构建模型预测和标准的一种方法。 data.table 中的错误，然后使用 ggplot.
绘图
library(data.table)
library(ggplot2)
setDT(aids) # convert aids to a data.table
aids[, c('pred', 'se', 'resid.scale'):=predict(glm(cases~date, data=.SD, family=poisson), type='response', se.fit=TRUE)]
ggplot(aids, aes(x=date))+
geom_ribbon(aes(ymin=pred-1.96*se, ymax=pred+1.96*se), fill='grey80')+
geom_line(aes(y=pred))+
geom_point(aes(y=cases))
或者，您可以让 ggplot 为您完成所有工作。
ggplot(aids, aes(x=date, y=cases))+
stat_smooth(method = glm, method.args=list(family=poisson))+
geom_point()
我已经为我的数据拟合了一个高斯 GLM 模型，我现在希望创建 95% CIs 并将它们拟合到我的数据中。我在绘图时遇到了几个问题，因为我无法让它们捕获我的数据，它们似乎只是绘制了与模型相同的线而没有捕获数据点。此外，我也不确定我是否以正确的方式在此处创建了 CIs 的均值。如果有人知道如何解决这个问题，我在下面输入了我的数据和代码
使用的数据
aids
cases quarter date
1 2 1 83.00
2 6 2 83.25
3 10 3 83.50
4 8 4 83.75
5 12 1 84.00
6 9 2 84.25
7 28 3 84.50
8 28 4 84.75
9 36 1 85.00
10 32 2 85.25
11 46 3 85.50
12 47 4 85.75
13 50 1 86.00
14 61 2 86.25
15 99 3 86.50
16 95 4 86.75
17 150 1 87.00
18 143 2 87.25
19 197 3 87.50
20 159 4 87.75
21 204 1 88.00
22 168 2 88.25
23 196 3 88.50
24 194 4 88.75
25 210 1 89.00
26 180 2 89.25
27 277 3 89.50
28 181 4 89.75
29 327 1 90.00
30 276 2 90.25
31 365 3 90.50
32 300 4 90.75
33 356 1 91.00
34 304 2 91.25
35 307 3 91.50
36 386 4 91.75
37 331 1 92.00
38 368 2 92.25
39 416 3 92.50
40 374 4 92.75
41 412 1 93.00
42 358 2 93.25
43 416 3 93.50
44 414 4 93.75
45 496 1 94.00
我的代码用于在绘图之前创建模型和间隔
#creating the model
model3 = glm(cases ~ date,
data = aids,
family = poisson(link='log'))
#now to add approx. 95% confidence envelope around this line
#predict again but at the linear predictor level along with standard errors
my_preds <- predict(model3, newdata=data.frame(aids), se.fit=T, type="link")
#calculate CI limit since linear predictor is approx. Gaussian
upper <- my_preds$fit+1.96*my_preds$se.fit #this might be logit not log
lower <- my_preds$fit-1.96*my_preds$se.fit
#transform the CI limit to get one at the level of the mean
upper <- exp(upper)/(1+exp(upper))
lower <- exp(lower)/(1+exp(lower))
#plotting data
plot(aids$date, aids$cases,
xlab = 'Date', ylab = 'Cases', pch = 20)
#adding CI lines
plot(aids$date, exp(my_preds$fit), type = "link",
xlab = 'Date', ylab = 'Cases') #add title
lines(aids$date,exp(my_preds$fit+1.96*my_preds$se.fit),lwd=2,lty=2)
lines(aids$date,exp(my_preds$fit-1.96*my_preds$se.fit),lwd=2,lty=2)
我目前没有数据点的结果，这里的模型是正确的，但是 CI 不是因为我没有数据点，所以 CI 是错误的，我认为在某个地方
编辑：对OP提供完整数据集的回应。
这最初是关于在同一图表上绘制数据和模型的问题，但已经发生了很大变化。您似乎对原始问题有答案。下面是解决其余问题的一种方法。
看看你的（和我的）图，泊松 glm 显然不是一个好的模型。换句话说，案例数量可能会随日期而变化，但也会受到模型中不存在的其他因素（外部回归变量）的影响。
仅绘制您的数据强烈表明您至少有两种甚至更多的制度：病例增长遵循不同模型的时间框架。
ggplot(aids, aes(x=date)) + geom_point(aes(y=cases))
这表明 segmented regression。与 R 中的大多数东西一样，有一个包（实际上不止一个）。下面的代码使用 segmented 包使用 1 个断点（两个区域）构建连续的泊松 glm。
library(data.table)
library(ggplot2)
library(segmented)
setDT(aids) # convert aids to a data.table
aids[, pred:=
predict(
segmented(glm(cases~date, .SD, family = poisson), seg.Z = ~date, npsi=1),
type='response', se.fit=TRUE)$fit]
ggplot(aids, aes(x=date))+ geom_line(aes(y=pred))+ geom_point(aes(y=cases))
请注意，我们需要告诉 segmented 断点的数量，而不是断点的位置 - 算法会为您计算出来。所以在这里，我们看到 3Q87 之前的状态使用泊松 glm 很好地建模，而之后的状态则不是。这是一种奇特的说法，即 87 年 3 季度左右“发生了一些事情”，它改变了疾病的进程（至少在这个数据中是这样）。
下面的代码做同样的事情，但有 1 到 4 个断点。
get.pred <- \(p.n, p.DT) {
fit <- glm(cases~date, p.DT, family=poisson)
seg.fit <- segmented(fit, seg.Z = ~date, npsi=p.n)
predict(seg.fit, type='response', se.fit=TRUE)[c('fit', 'se.fit')]
}
gg.dt <- rbindlist(lapply(1:4, \(x) { copy(aids)[, c('pred', 'se'):=get.pred(x, .SD)][, npsi:=x] } ))
ggplot(gg.dt, aes(x=date))+
geom_ribbon(aes(ymin=pred-1.96*se, ymax=pred+1.96*se), fill='grey80')+
geom_line(aes(y=pred))+
geom_point(aes(y=cases))+
facet_wrap(~npsi)
请注意，第一个断点的位置似乎没有改变，而且，尽管使用了泊松 glm，但除第一个区域外，所有区域的增长都是线性的。
包文档中描述了 goodness-of-fit 指标，可以帮助您确定有多少断点与您的数据最一致。
最后，还有 mcp 包，它更强大一些，但使用起来也更复杂一些。
原始响应：这是构建模型预测和标准的一种方法。 data.table 中的错误，然后使用 ggplot.
library(data.table)
library(ggplot2)
setDT(aids) # convert aids to a data.table
aids[, c('pred', 'se', 'resid.scale'):=predict(glm(cases~date, data=.SD, family=poisson), type='response', se.fit=TRUE)]
ggplot(aids, aes(x=date))+
geom_ribbon(aes(ymin=pred-1.96*se, ymax=pred+1.96*se), fill='grey80')+
geom_line(aes(y=pred))+
geom_point(aes(y=cases))
或者，您可以让 ggplot 为您完成所有工作。
ggplot(aids, aes(x=date, y=cases))+
stat_smooth(method = glm, method.args=list(family=poisson))+
geom_point()