月趋势的负二项回归

Question

我读了一篇关于负二项回归的论文：“我们使用负二项回归(合并过度分散)对每月大肠杆菌血流感染和大肠杆菌尿路感染的数量进行建模，假设基础人群相同(没有偏移量)。”图如下所示

​

我也有一组数据，想要用月份/年来计算感染情况，我该怎么做呢？非常感谢

df <- structure(list(Year = c(2013, 2013, 2013, 2013, 2013, 2013, 2013, 
2013, 2013, 2013, 2013, 2013, 2014, 2014, 2014, 2014, 2014, 2014, 
2014, 2014, 2014, 2014, 2014, 2014, 2015, 2015, 2015, 2015, 2015, 
2015, 2015, 2015, 2015, 2015, 2015, 2015), Month = c(1, 2, 3, 
4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 
11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12), Incidence = c(2.25538216197745, 
3.49502862307924, 2.76311704439615, 2.9836483329794, 3.09375, 
3.0368028900429, 3.82920688208141, 3.9154960734432, 3.33517393705135, 
3.54593329432417, 3.27586206896552, 3.25655281969817, 3.35912052117264, 
3.21672101986362, 2.78237182605312, 2.58435732397113, 2.72516428295323, 
3.1227603153476, 2.6300688599847, 2.66324718879463, 2.62653374233129, 
2.45256358498183, 2.39520958083832, 3.58683926645092, 3.41995942421022, 
3.61001317523057, 2.62718158187895, 2.86944045911047, 2.77978993118435, 
2.89282762420792, 2.69410829432029, 3.22232223222322, 3.39818882811799, 
3.36725958337297, 2.90030211480363, 3.20789124668435), Inpatient = c(8779, 
6638, 9663, 9418, 9600, 8858, 9532, 9041, 9055, 8545, 9280, 10072, 
9824, 6746, 10279, 10254, 10348, 9767, 10456, 10138, 10432, 9908, 
9853, 11124, 10351, 7590, 10772, 11152, 11044, 10889, 11321, 
11110, 11153, 10513, 11585, 12064), infection = c(198, 232, 267, 
281, 297, 269, 365, 354, 302, 303, 304, 328, 330, 217, 286, 265, 
282, 305, 275, 270, 274, 243, 236, 399, 354, 274, 283, 320, 307, 
315, 305, 358, 379, 354, 336, 387)), row.names = c(NA, -36L), class = c("tbl_df", 
"tbl", "data.frame"))

参考文献: Vihta K D，Stoesser N，Llewelyn M J，et al.1998年至2016年英国牛津郡大肠杆菌血液感染、尿路感染和抗生素敏感性的趋势:电子健康recordsJ的研究。柳叶刀传染病，2018，18(10)：1138-1149。

Answer 1

使用上面的数据，可以做到以下几点：

library(MASS) # for function glm.nb
library(ggplot2)
library(broom) # for tidy model outputs

创建日期，以简化绘图

df$t <- as.Date(paste("01", df$Month, df$Year, sep = "-"), format = "%d-%m-%Y")

绘制数据。geom_smooth使用日期作为预测值，添加趋势线和置信区间。

p <- ggplot(data = df, aes(x = t, y = infection)) + 
  geom_point() + 
  geom_smooth(method = "glm.nb")
p

要执行回归，请将感染计数设置为因变量，并将第n个月设置为自变量，在month_as_integer之下。

df$month_as_integer <- seq_along(df$Month)
m1 <- glm.nb(infection ~ month_as_integer, data = df)

使用broom包中的tidy，可以获得作为数据帧的估计和置信区间。

out1 <- as.data.frame(tidy(m1, exponentiate = TRUE, conf.int = TRUE) )

out1
              term  estimate   std.error  statistic     p.value   conf.low  conf.high
1      (Intercept) 264.44399 0.048006493 116.184897 0.000000000 240.943378 290.556355
2 month_as_integer   1.00697 0.002250993   3.085763 0.002030303   1.002569   1.011394