Dylan_Gomes
Reputation: 2232
How do I manually plot SE/CI around a model estimate
I am trying to manually plot model estimates on top of data. My real problem is far more complicated than this, so I want to avoid using predict if I can, and would prefer to understand how these predictions would be calculated rather than relying on some package.
(data for a reproducible example at the bottom.)
So I first run a model, and grab the model estimates and standard errors:
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])
Then, I create a sequence of x values to predict from
x<-seq(from=min(Dat$DayL50),to=max(Dat$DayL50),length.out = length(Dat$DayL50))
Then I predict the y values with two different methods (using predict and writing the equation that should do the same thing)
ypred<-exp(dB_est*x+Int)
y<-predict(glmmLep,list(DayL50=x),type="response",se.fit = T)
We plot the two predicted lines (one as a smaller red line on top):
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))
We see that the equation I wrote works the same as the predict function. All good. However, when I go to add the SE / 95% CI ribbon around that line, I run into issues trying to recreate it (here I left as SE, since 95%CI leads to more of an unwieldy plot). I have played with the formula in many different ways, and can't seem to get it. For some reason, I cannot seem to find any posts about it, but perhaps I am not using the correct search terms. Can anyone explain to me what I am missing here. It seems as if I am missing quite a bit of complexity in my error ribbons (outlined in red).
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))
Or with 95%CI, like my predict ribbon, which is even further off:
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))
Upvotes: 7
Views: 900
Answers (1)
Allan Cameron
Reputation: 173858
This is more of a stats question than a programming question. The only thing "wrong" with your code is your method of calculating the standard errors for the fit of your model.
The se.fit of the prediction doesn't come from a single coefficient's standard error, but is calculated from a linear combination of the model's parameters (intercept and DayL50, or β₀ and β₁) using the model's variance-covariance matrix. So if I want to know the standard error of the fit with the given intercept i.e. (1 * β₀) and x = 50, i.e. (50 * β₁), I can do
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
So, provided you are able to get the covariance matrix for your model, you can calculate the standard errors to match se.fit as in the following worked example:
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))
Now we have our data frames, we can make the plot. Here I will show the equivalence of the predict and manual methods by overlapping blue (pred) and dashed red (manual) ribbons around the confidence intervals produced by both methods:
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))
Created on 2020-05-03 by the reprex package (v0.3.0)
Upvotes: 4
Related Questions
- Adding 95% confidence interval of prediction using ggplot2
- How to plot predicted data of a GLMM from the package glmmTMB?
- Plotting a 95% confidence interval band around a predicted regression line from a linear mixed model
- Plotting 95% confidence intervals of log-log linear model predicted values
- How to plot a CI only in the predicted line
- Adding confidence intervals from model predictions in mixed models in R -- ggplot2?
- why do ggplot2 95%CI and prediction 95%CI calculated manually differ?
- How to plot estimate values for a lmer regression model in R?
- How to plot mixed-effects model estimates in ggplot2 in R?
- How to visulaize linear model prediction in ggplot along with confidence interval?