# Constructed teaching example, not a random draw or a realistic effect estimate.
# stats supplies glm(); ggplot2 supplies the room-by-diet plot.
library("ggplot2")

separation_data <- expand.grid(
  pen_in_cell = 1:2,
  diet = c("Control", "High fibre"),
  room = paste("Room", 1:4)
)
separation_data$diet <- factor(separation_data$diet,
                              levels = c("Control", "High fibre"))
separation_data$room <- factor(separation_data$room)
separation_data$pen_id <- sprintf("P%02d", seq_len(nrow(separation_data)))
# This additive room + diet predictor strictly separates zeros from ones.
separation_data$eta_separator <- with(separation_data,
  1 - 2 * (diet == "High fibre") + c(0, -2, 2, 0)[as.integer(room)])
separation_data$response <- as.integer(separation_data$eta_separator > 0)
stopifnot(all(with(separation_data, (2 * response - 1) * eta_separator > 0)))
# Both outcomes occur within each diet overall, so diet alone cannot separate.
stopifnot(all(table(separation_data$diet, separation_data$response) > 0))

fit <- glm(response ~ diet + room, data = separation_data,
           family = binomial(link = "logit"))
print(summary(fit))
print(c(converged = fit$converged, iterations = fit$iter))
print(range(fitted(fit)))

# Horizontal offsets expose both pens; binary outcomes are never jittered vertically.
separation_data$plot_x <- as.integer(separation_data$diet) +
  ifelse(separation_data$pen_in_cell == 1, -.075, .075)
separation_plot <- ggplot(separation_data,
    aes(plot_x, response, colour = diet)) +
  geom_hline(yintercept = c(0, 1), colour = "grey85", linewidth = .4) +
  geom_point(size = 3.5) +
  facet_wrap(~room, ncol = 2) +
  scale_x_continuous(breaks = 1:2, labels = c("Control", "High fibre"),
                     limits = c(.7, 2.3)) +
  scale_y_continuous(breaks = c(0, 1), limits = c(-.15, 1.15)) +
  scale_colour_manual(values = c("Control" = "#003d73", "High fibre" = "#cc6600")) +
  labs(x = NULL, y = "Any pig treated (0/1)") +
  theme_minimal(base_size = 15) +
  theme(legend.position = "none", panel.grid = element_blank(),
        strip.text = element_text(face = "bold"),
        panel.spacing = grid::unit(.8, "lines"))
print(separation_plot)
# No finite unpenalised logistic MLE exists. The returned finite numbers depend
# on numerical stopping. A convergence flag or tiny deviance is not a cure.
# Source: https://stat.ethz.ch/CRAN/web/packages/detectseparation/vignettes/separation.html