Breakfast Foods: statistical model

set.seed(0)

n <- 100

tibble(
  Pancakes     = rnorm(n = n, mean =  0.8, sd = 1.0),
  Waffles      = rnorm(n = n, mean =  0.0, sd = 2.5),
  `French Toast` = rnorm(n = n, mean = -0.5, sd = 2.0),
) -> data

data |>
pivot_longer(
  cols = everything(),
  names_to = "Breakfast Food",
  values_to = "Quality"
  ) |>
  mutate(
    `Breakfast Food` = factor(`Breakfast Food`, levels = c("Pancakes", "French Toast", "Waffles"))
  ) -> data_long

data_long |>
ggplot(aes(x = `Breakfast Food`, y = Quality, color = `Breakfast Food`)) +
  geom_boxplot(show.legend = F) +
  scale_color_manual(values = c("skyblue3", "hotpink3", "seagreen4"))+
  theme_linedraw()

p <- list(
  pancakes = list(
    sd = sd(data$Pancakes),
    mean = mean(data$Pancakes)
  ),
  french_toast = list(
    sd = sd(data$`French Toast`),
    mean = mean(data$`French Toast`)
  ),
  waffles = list(
    sd = sd(data$Waffles),
    mean = mean(data$Waffles)
  )
)

xcrit = 3.0

p$pancakes$p <- round(1-pnorm(xcrit, mean = p$pancakes$mean, sd = p$pancakes$sd), 2)
p$french_toast$p <- round(1-pnorm(xcrit, mean = p$french_toast$mean, sd = p$french_toast$sd), 2)
p$waffles$p <- round(1-pnorm(xcrit, mean = p$waffles$mean, sd = p$waffles$sd), 2)

data |>
  ggplot() +
  geom_function(fun = dnorm, args = list(mean = p$pancakes$mean, sd = p$pancakes$sd), color = "skyblue3") +
  stat_function(fun = dnorm, args = list(mean = p$pancakes$mean, sd = p$pancakes$sd), geom = "area", fill = "skyblue", alpha = 0.2, xlim = c(xcrit, 10)) +
  geom_label(x = 6.0, y = 0.25, label = paste0("Pancakes : P = ", p$pancakes$p), hjust = "right", color = "skyblue3")+
  geom_function(fun = dnorm, args = list(mean = p$waffles$mean, sd = p$waffles$sd), color = "seagreen4") +
  stat_function(fun = dnorm, args = list(mean = p$waffles$mean, sd = p$waffles$sd), geom = "area", fill = "seagreen4", alpha = 0.2, xlim = c(xcrit, 10)) +
  geom_label(x = 6.0, y = 0.20, label = paste0("Waffles: P = ", p$waffles$p), hjust = "right", color = "seagreen4")+
  geom_function(fun = dnorm, args = list(mean = p$french_toast$mean, sd = p$french_toast$sd), color = "hotpink3") +
  stat_function(fun = dnorm, args = list(mean = p$french_toast$mean, sd = p$french_toast$sd), geom = "area", fill = "hotpink3", alpha = 0.2, xlim = c(xcrit, 10)) +
  geom_label(x = 6.0, y = 0.15, label = paste0("French Toast: P = ", p$french_toast$p), hjust = "right", color = "hotpink3")+
  geom_vline(xintercept = xcrit, linetype = "dotted") +
  labs(x = "Quality (sd)", y = "Frequency")+
  scale_x_continuous(limits = c(-4, 6), breaks = seq(-4.0, 6.0, 1))+
  theme_linedraw()