如何围绕模型估计手动绘制 SE/CI
Posted
技术标签:
【中文标题】如何围绕模型估计手动绘制 SE/CI【英文标题】：How do I manually plot SE/CI around a model estimate 【发布时间】：2020-07-20 09:49:30 【问题描述】：我正在尝试根据数据手动绘制模型估计值。我真正的问题远比这复杂，所以我想尽可能避免使用predict，并且更愿意了解这些预测是如何计算的，而不是依赖于某个包。
（底部可重复示例的数据。）
所以我首先运行一个模型，并获取模型估计值和标准误差：
library(glmmTMB)
glmmLep<-glmmTMB(Lepidoptera ~ DayL50,
data=Dat, family=nbinom2(link="log") )
dB_est<-(summary(glmmLep)$coeff$cond[2,1])
dB_SE<-(summary(glmmLep)$coeff$cond[2,2])
Int<-(summary(glmmLep)$coeff$cond[1,1])
Int_SE<-(summary(glmmLep)$coeff$cond[1,2])
然后，我创建一个 x 值序列来预测
x<-seq(from=min(Dat$DayL50),to=max(Dat$DayL50),length.out = length(Dat$DayL50))
然后我用两种不同的方法预测 y 值（使用 predict 并编写应该做同样事情的方程）
ypred<-exp(dB_est*x+Int)
y<-predict(glmmLep,list(DayL50=x),type="response",se.fit = T)
我们绘制两条预测线（一条为顶部较小的红线）：
ggplot(aes(x=DayL50,y=Lepidoptera),data=Dat)+
geom_point(size=2)+
geom_line(aes(y=y$fit,x=x),size=2)+
geom_ribbon(aes(ymax=y$fit+1.96*y$se.fit,ymin=y$fit-1.96*y$se.fit,x=x),alpha=0.2)+
geom_line(aes(y=ypred,x=x),size=1,color="red")+
# geom_ribbon(aes(ymax=ymax,ymin=ymin,x=x),alpha=0.2,color="red")+
coord_cartesian(ylim=c(0,1000))
我们看到我写的公式与predict 函数的工作原理相同。都好。但是，当我在该行周围添加 SE / 95% CI 功能区时，我在尝试重新创建它时遇到了问题（这里我保留为 SE，因为 95% CI 会导致更多笨拙的情节）。我以许多不同的方式玩过这个公式，但似乎无法理解。出于某种原因，我似乎找不到任何关于它的帖子，但也许我没有使用正确的搜索词。谁能向我解释我在这里缺少什么。似乎我的错误功能区中缺少相当多的复杂性（以红色标出）。
ymin<-exp((dB_est-dB_SE)*x+(Int))
ymax<-exp((dB_est+dB_SE)*x+(Int))
ggplot(aes(x=DayL50,y=Lepidoptera),data=Dat)+
geom_point(size=2)+
geom_line(aes(y=y$fit,x=x),size=2)+
geom_ribbon(aes(ymax=y$fit+1.96*y$se.fit,ymin=y$fit-1.96*y$se.fit,x=x),alpha=0.2)+
geom_line(aes(y=ypred,x=x),size=1,color="red")+
geom_ribbon(aes(ymax=ymax,ymin=ymin,x=x),alpha=0.2,color="red")+
coord_cartesian(ylim=c(0,1000))
或者使用 95%CI，比如我的 predict 功能区，它离得更远：
ymin<-exp((dB_est-1.96*dB_SE)*x+(Int))
ymax<-exp((dB_est+1.96*dB_SE)*x+(Int))
ggplot(aes(x=DayL50,y=Lepidoptera),data=Dat)+
geom_point(size=2)+
geom_line(aes(y=y$fit,x=x),size=2)+
geom_ribbon(aes(ymax=y$fit+1.96*y$se.fit,ymin=y$fit-1.96*y$se.fit,x=x),alpha=0.2)+
geom_line(aes(y=ypred,x=x),size=1,color="red")+
geom_ribbon(aes(ymax=ymax,ymin=ymin,x=x),alpha=0.2,color="red")+
coord_cartesian(ylim=c(0,1000))
Dat<-structure(list(Lepidoptera = c(0L, 0L, 1L, 0L, 1L, 1L, 807L,
103L, 6L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 63L, 0L, 0L, 3L, 1L, 94L,
0L, 0L, 0L, 0L, 27L, 0L, 0L, 117L, 0L, 0L, 95L, 0L, 0L, 0L, 11L,
0L, 0L, 0L, 0L, 0L, 0L, 2L, 11L, 0L, 0L, 0L, 5L, 26L, 0L, 0L,
0L, 0L, 0L, 76L, 0L, 610L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 56L, 0L,
1L, 119L, 0L, 14L, 0L, 0L, 302L, 0L, 0L, 113L, 312L, 0L, 0L,
0L, 1L, 323L, 53L, 0L, 0L, 0L, 0L, 0L, 0L, 2L, 2L, 720L, 0L,
2L, 0L, 2L, 152L, 0L, 1L, 0L, 2L, 172L, 0L, 0L, 55L, 0L, 136L,
0L, 5L, 0L, 108L, 0L, 0L, 912L, 34L, 0L, 1L, 6L, 1405L, 3L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 14L, 1236L, 7L, 8L, 11L, 231L, 1L, 0L,
163L, 531L, 7L, 2L, 155L, 3L, 0L, 16L, 69L, 2L, 1084L, 5L, 7L,
120L, 2L, 1L, 48L, 1L, 0L, 1303L, 107L, 0L, 0L, 0L, 463L, 13L,
36L, 2L, 0L, 0L, 2L, 0L, 77L, 0L, 0L, 374L, 0L, 0L, 18L, 1L,
0L, 0L, 158L, 269L, 0L, 0L, 0L, 1L, 16L, 6L, 0L, 1L, 258L, 0L,
8L, 0L, 22L, 2838L, 226L, 0L, 8L, 302L, 4196L, 16L, 1L, 0L, 0L,
1357L, 6L, 0L, 3L, 1L, 0L, 304L, 2257L, 0L, 0L, 2L, 34L, 142L,
0L, 0L, 2L, 0L, 402L, 154L, 480L, 461L, 1463L, 0L, 0L, 0L, 116L,
0L, 6L, 0L, 0L, 0L, 7L, 0L, 276L, 0L, 0L, 4L, 0L, 32L, 0L, 0L,
6L, 0L, 40L, 1L, 0L, 71L, 0L, 4L, 0L, 0L, 96L, 10L, 0L, 0L, 0L,
0L, 4L, 0L, 22L, 0L, 0L, 0L, 1L, 18L, 83L, 0L, 0L, 38L, 207L,
0L, 0L, 0L, 0L, 0L, 506L, 0L, 0L, 1L, 0L, 0L, 0L, 708L, 0L, 1L,
39L, 0L, 588L, 0L, 0L, 8L, 154L, 1L, 0L, 0L, 0L, 0L, 3L, 129L,
0L, 1L, 0L, 0L, 0L, 452L, 59L, 0L, 2L, 596L, 0L, 4L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 23L, 0L, 0L, 0L, 0L, 46L, 7L, 0L, 0L, 0L,
55L, 5L, 0L, 4L, 0L, 51L, 0L, 0L, 1L, 9L, 1L, 84L, 43L, 0L, 2L,
1L, 95L, 1L, 259L, 0L, 0L, 0L, 6L, 427L, 0L, 66L, 0L, 3L, 752L,
109L, 2L, 0L, 0L, 0L, 4L, 5L, 0L, 151L, 0L, 4L, 1L, 0L, 32L,
0L, 0L, 0L, 3L, 122L, 47L, 1L, 0L, 7L, 52L, 174L, 0L, 0L, 1L,
23L, 5L, 1L, 0L, 932L, 2L, 290L, 3L, 2078L, 48L, 0L, 3L, 0L,
0L, 37L, 0L, 169L, 0L, 0L, 142L, 2052L, 1L, 0L, 377L, 0L, 1L,
3857L, 19L, 220L, 2332L, 0L, 17L, 1L, 926L, 16L, 6815L, 39L,
0L, 6L, 289L, 626L, 1L, 1L, 0L, 1L, 0L, 30L, 0L, 0L, 395L, 0L,
450L, 1L, 679L, 0L, 0L, 17L, 817L, 4L, 10L, 300L, 41L, 1L, 1L,
164L), DayL50 = c(62.2, 45.4, 71.8, 60.4, 60.4, 60.4, 60.4, 60.4,
45.1, 45.1, 45.1, 45.1, 69.5, 71.3, 71.3, 71.3, 70.7, 74, 69.4,
69.4, 69.4, 69.4, 69.4, 67.3, 54.9, 71.5, 71.5, 71.5, 71.5, 71.5,
71.5, 74.1, 74.1, 74.1, 74.1, 66.5, 66.5, 66.5, 66.5, 66.5, 73.2,
55.8, 55.8, 70.3, 70.3, 70.3, 70.3, 68.2, 68.2, 68.2, 68.2, 68.2,
48.4, 50.6, 73.2, 73.2, 73.2, 73.2, 73.2, 52.2, 61.2, 66, 68.2,
58.1, 59.9, 59.9, 59.9, 59.9, 59.9, 54.8, 54.8, 54.8, 54.8, 54.8,
63.9, 63.9, 63.9, 63.9, 63.9, 69.8, 69.8, 69.8, 69.8, 69.8, 45.4,
47.2, 54.5, 48.8, 68.4, 39.7, 45.4, 45.4, 45.4, 45.4, 45.4, 46.8,
46.8, 46.8, 46.8, 46.8, 54.3, 54.3, 54.3, 54.3, 54.3, 49.2, 49.2,
49.2, 49.2, 49.2, 68.8, 68.8, 68.8, 68.8, 68.8, 39.6, 39.6, 39.6,
39.6, 39.6, 41.2, 70.7, 62.1, 44.5, 70.1, 49.8, 53.8, 72.5, 61.5,
61.5, 61.5, 61.5, 45.4, 45.4, 45.4, 45.4, 45.4, 69.5, 70.8, 70.8,
70.8, 70.8, 66.3, 73.2, 73.2, 73.2, 73.2, 73.2, 50.4, 50.4, 50.4,
50.4, 50.4, 54.1, 54.1, 54.1, 54.1, 54.1, 73.5, 67.9, 67.9, 67.9,
67.9, 67.9, 70.7, 74, 71.5, 74.1, 74.1, 74.1, 74.1, 74.1, 43.8,
71.5, 71.5, 71.5, 74.1, 74.1, 74.1, 74.1, 74.1, 48.7, 69, 69,
69, 69, 65.8, 45.4, 45.4, 45.4, 45.4, 47.9, 47.9, 47.9, 47.9,
39.9, 39.9, 39.9, 39.9, 39.9, 39.9, 67.7, 67.7, 67.7, 67.7, 70.9,
70.9, 70.9, 70.9, 70.9, 70.9, 57.3, 61.2, 59.9, 59.9, 59.9, 59.9,
63.9, 63.9, 63.9, 63.9, 63.9, 70, 70.4, 70.4, 63.6, 45.2, 45.2,
45.2, 45.2, 45.2, 55.1, 64.5, 64.1, 64.1, 64.1, 64.1, 54, 54,
54, 54, 54, 65, 65, 65, 65, 65, 61.9, 64.2, 62.3, 62.3, 62.3,
36.5, 64.2, 64.2, 64.2, 64.2, 64.2, 58.8, 38.3, 38.3, 38.3, 38.3,
38.3, 59.1, 59.1, 59.1, 59.1, 59.1, 58.6, 66.1, 66.1, 66.1, 66.1,
76.5, 76.5, 76.5, 76.5, 76.5, 76.5, 70.5, 72.7, 70.3, 70.3, 70.3,
70.3, 71.8, 71.8, 71.8, 71.8, 71.8, 45.4, 71, 37.2, 37.2, 37.2,
37.2, 61.2, 65, 69.8, 69.8, 69.8, 69.8, 69.8, 60.3, 60.3, 60.3,
60.3, 60.3, 64.9, 64.9, 64.9, 64.9, 64.9, 47.7, 54.3, 69.3, 54.5,
54.5, 54.5, 54.5, 54.5, 54.5, 47.8, 47.8, 47.8, 47.8, 47.8, 54.6,
54.6, 54.6, 54.6, 54.6, 69.1, 69.1, 69.1, 69.1, 69.1, 38.7, 57.1,
35.9, 35.9, 35.9, 35.9, 35.9, 56.7, 56.7, 56.7, 56.7, 56.7, 51.9,
61.8, 52.1, 52.1, 52.1, 52.1, 52.1, 63.2, 63.2, 63.2, 63.2, 63.2,
71.9, 74.7, 72, 72, 72, 72, 72, 74.6, 74.6, 74.6, 74.6, 74.6,
62, 69, 61.1, 61.1, 61.1, 61.1, 61.1, 68.4, 68.4, 68.4, 68.4,
68.4, 45.3, 58.6, 43.8, 43.8, 43.8, 43.8, 43.8, 60.3, 60.3, 60.3,
60.3, 60.3, 54, 54.4, 64.8, 55, 55, 55, 55, 55, 71, 71, 71, 71,
71, 52.8, 52.8, 52.8, 52.8, 52.8, 63.9, 63.9, 63.9, 63.9, 35.1,
35.1, 35.1, 35.1, 35.1, 35.1, 78.9, 78.9, 78.9, 78.9, 78.9, 48,
66.6, 54.2, 54.2, 54.2, 54.2, 54.2, 54.2, 49.5, 49.5, 49.5, 49.5,
49.5, 56.3, 56.3, 56.3, 56.3, 66.6, 66.6, 66.6, 66.6, 66.6)), class = "data.frame", row.names = c(1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L,
16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L,
29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 37L, 38L, 39L, 40L, 41L,
42L, 43L, 44L, 45L, 46L, 47L, 48L, 49L, 50L, 51L, 52L, 53L, 54L,
55L, 56L, 57L, 58L, 59L, 60L, 61L, 62L, 63L, 64L, 65L, 66L, 67L,
68L, 69L, 70L, 71L, 72L, 73L, 74L, 75L, 76L, 77L, 78L, 79L, 80L,
81L, 82L, 83L, 84L, 85L, 86L, 87L, 88L, 89L, 90L, 91L, 92L, 93L,
94L, 95L, 96L, 97L, 98L, 99L, 100L, 101L, 102L, 103L, 104L, 105L,
106L, 107L, 108L, 109L, 110L, 111L, 112L, 113L, 114L, 115L, 116L,
117L, 118L, 119L, 120L, 121L, 122L, 123L, 124L, 125L, 126L, 127L,
128L, 129L, 130L, 131L, 132L, 133L, 134L, 135L, 136L, 137L, 138L,
139L, 140L, 141L, 142L, 143L, 144L, 145L, 146L, 147L, 148L, 149L,
150L, 151L, 152L, 153L, 154L, 155L, 156L, 157L, 158L, 159L, 160L,
161L, 162L, 163L, 164L, 165L, 166L, 167L, 168L, 169L, 170L, 171L,
172L, 173L, 175L, 176L, 177L, 178L, 179L, 180L, 181L, 182L, 183L,
184L, 185L, 186L, 187L, 188L, 189L, 190L, 191L, 192L, 193L, 194L,
195L, 196L, 197L, 198L, 199L, 200L, 201L, 202L, 203L, 204L, 205L,
206L, 207L, 208L, 209L, 210L, 211L, 212L, 213L, 214L, 215L, 216L,
217L, 218L, 219L, 220L, 221L, 222L, 223L, 224L, 225L, 226L, 227L,
228L, 229L, 230L, 231L, 232L, 233L, 234L, 235L, 236L, 237L, 238L,
239L, 240L, 241L, 242L, 243L, 244L, 245L, 246L, 247L, 248L, 249L,
250L, 251L, 252L, 253L, 254L, 255L, 256L, 257L, 258L, 259L, 260L,
262L, 263L, 264L, 265L, 266L, 267L, 268L, 269L, 270L, 271L, 272L,
273L, 274L, 275L, 276L, 277L, 278L, 279L, 280L, 281L, 282L, 283L,
284L, 285L, 286L, 287L, 288L, 289L, 290L, 291L, 292L, 293L, 294L,
295L, 296L, 297L, 298L, 299L, 300L, 301L, 302L, 303L, 304L, 305L,
306L, 307L, 308L, 309L, 310L, 311L, 312L, 313L, 314L, 315L, 316L,
317L, 318L, 319L, 320L, 321L, 322L, 323L, 324L, 325L, 326L, 327L,
328L, 329L, 330L, 331L, 332L, 333L, 334L, 335L, 336L, 337L, 338L,
339L, 340L, 341L, 342L, 343L, 344L, 345L, 346L, 347L, 348L, 349L,
350L, 351L, 352L, 353L, 354L, 355L, 356L, 357L, 358L, 359L, 360L,
361L, 362L, 363L, 364L, 365L, 366L, 367L, 368L, 369L, 370L, 371L,
372L, 373L, 374L, 375L, 376L, 377L, 378L, 379L, 380L, 381L, 382L,
383L, 384L, 385L, 386L, 387L, 388L, 389L, 390L, 391L, 392L, 393L,
394L, 395L, 396L, 397L, 398L, 399L, 400L, 401L, 402L, 403L, 404L,
405L, 406L, 407L, 408L, 409L, 410L, 411L, 412L, 413L, 414L, 415L,
416L, 417L, 418L, 419L, 420L, 421L, 422L, 423L, 424L, 425L, 426L,
427L, 428L, 429L, 430L, 431L, 432L, 433L, 434L, 435L, 436L, 437L,
438L, 439L, 440L, 441L, 442L, 443L, 444L, 445L, 446L, 447L, 448L,
449L, 450L, 451L, 452L, 453L, 454L, 455L))
【问题讨论】：
我在这里标记了lme4，因为glmmTMB 没有标记，我假设这两个包的问题都是一样的。希望绘制这些效果的人可能正在使用任何一个包，所以我认为这是一个相关的标签。如果这是对标签的不当使用，请纠正我
【参考方案1】：
这更像是一个统计问题，而不是编程问题。您的代码唯一“错误”的是您计算模型拟合标准误差的方法。
预测的se.fit 不是来自单个系数的标准误差，而是使用模型的方差-协方差矩阵根据模型参数（截距和 DayL50，或 β₀ 和 β₁）的线性组合计算得出的。因此，如果我想知道给定截距的拟合标准误差，即 (1 * β₀) 和 x = 50，即 (50 * β₁)，我可以这样做
fit_at <- c(1, 50)
covariance_matrix <- vcov(mod)$cond
se_50 <- sqrt(t(fit_at) %*% covariance_matrix %*% fit_at)
se_50
#> [,1]
#> [1,] 0.03696078
因此，只要您能够获得模型的协方差矩阵，您就可以计算标准误差以匹配 se.fit，如下面的示例所示：
library(glmmTMB)
library(ggplot2)
# Generate our model
mod <- glmmTMB(Lepidoptera ~ DayL50, data = Dat, family = nbinom2())
# Create sequence for x variable in predictions
x <- with(Dat, seq(min(DayL50), max(DayL50), length.out = length(DayL50)))
# Generate a data frame of prediction with upper and lower SEM using predict()
pred <- predict(mod, list(DayL50 = x), type = "link", se.fit = TRUE)
pred <- data.frame(x = x, y = exp(pred$fit),
ymin = exp(pred$fit - 1.96 * pred$se.fit),
ymax = exp(pred$fit + 1.96 * pred$se.fit))
# Generate a data frame of prediction with upper and lower SEM manually
Int <- summary(mod)$coefficients$cond[1]
dB <- summary(mod)$coefficients$cond[2]
se <- sapply(x, function(y) sqrt(t(c(1, y)) %*% vcov(mod)$cond %*% c(1, y)))
manual <- data.frame(x = x, y = exp(x * dB + Int),
ymax = exp(x * dB + Int + 1.96 * se),
ymin = exp(x * dB + Int - 1.96 * se))
现在我们有了数据框，我们可以制作绘图了。在这里，我将通过在两种方法产生的置信区间周围重叠蓝色 (pred) 和红色虚线 (manual) 色带来展示 predict 和手动方法的等效性：
ggplot() +
geom_point(aes(x = DayL50, y = Lepidoptera), data = Dat, size = 2) +
geom_line(aes(x = x, y = y), data = pred, size = 2) +
geom_ribbon(aes(x = x, ymin = ymin, ymax = ymax), data = pred,
colour = "blue", size = 2, alpha = 0.2) +
geom_ribbon(aes(x = x, ymin = ymin, ymax = ymax), data = manual,
colour = "red", size = 2, linetype = 2, alpha = 0) +
coord_cartesian(ylim = c(0, 1000))
由reprex package (v0.3.0) 于 2020 年 5 月 3 日创建
【讨论】：
以上是关于如何围绕模型估计手动绘制 SE/CI的主要内容，如果未能解决你的问题，请参考以下文章
如何使用 Google MAps API v2 围绕引脚绘制圆圈