Plot Simultaneous and Pointwise Confidence Bands for Linear Regression

library(dplyr)

setosa <- iris %>% filter(Species == "setosa")

                %>% select(Sepal.Length,

                           Sepal.Width, 

                           Species)

virginica <- iris %>% filter(Species == "virginica")

                    %>% select(Sepal.Length, 

                               Sepal.Width,

                               Species)

setosa <- setosa %>% arrange(Sepal.Length)

virginica <- virginica %>% arrange(Sepal.Length)

 compute linear models

Model <- as.formula(Sepal.Width ~ poly(Sepal.Length,2))

fit1  <- lm(Model, data = setosa)

fit2  <- lm(Model, data = virginica)

X1  <- model.matrix(fit1)

X2  <- model.matrix(fit2)

cht1 <- multcomp::glht(fit1, linfct = X1) 

cht2 <- multcomp::glht(fit2, linfct = X2) 

cc1 <- confint(cht1); cc1 <- as.data.frame(cc1$confint)

cc2 <- confint(cht2); cc2 <- as.data.frame(cc2$confint)

setosa$LowerBound <- cc1$lwr

setosa$UpperBound <- cc1$upr

virginica$LowerBound <- cc2$lwr

virginica$UpperBound <- cc1$upr

mydata <- rbind(setosa, virginica)

library(ggplot2)

ggplot(data = mydata, aes(x = Sepal.Length,

                          y = Sepal.Width, 

                          color = Species)) + 

  geom_point() +

  geom_smooth(method ="lm", formula = y ~ poly(x,2)) +

  geom_ribbon(aes(ymin = LowerBound,

                  ymax = UpperBound),

              alpha = 0.5, 

              fill = "grey70")

library(dplyr)

setosa <- iris %>% filter(Species == "setosa")

                %>% select(Sepal.Length,

                           Sepal.Width, 

                           Species)

virginica <- iris %>% filter(Species == "virginica")

                    %>% select(Sepal.Length, 

                               Sepal.Width,

                               Species)

setosa <- setosa %>% arrange(Sepal.Length)

virginica <- virginica %>% arrange(Sepal.Length)

 compute linear models

Model <- as.formula(Sepal.Width ~ poly(Sepal.Length,2))

fit1  <- lm(Model, data = setosa)

fit2  <- lm(Model, data = virginica)

X1  <- model.matrix(fit1)

X2  <- model.matrix(fit2)

cht1 <- multcomp::glht(fit1, linfct = X1) 

cht2 <- multcomp::glht(fit2, linfct = X2) 

cc1 <- confint(cht1); cc1 <- as.data.frame(cc1$confint)

cc2 <- confint(cht2); cc2 <- as.data.frame(cc2$confint)

setosa$LowerBound <- cc1$lwr

setosa$UpperBound <- cc1$upr

virginica$LowerBound <- cc2$lwr

virginica$UpperBound <- cc1$upr

mydata <- rbind(setosa, virginica)

library(ggplot2)

ggplot(data = mydata, aes(x = Sepal.Length,

                          y = Sepal.Width, 

                          color = Species)) + 

  geom_point() +

  geom_smooth(method ="lm", formula = y ~ poly(x,2)) +

  geom_ribbon(aes(ymin = LowerBound,

                  ymax = UpperBound),

              alpha = 0.5, 

              fill = "grey70")