# Display every row when checking the number of pigs in each pen.
options(tibble.print_max = Inf)
# Class 8. Run from this project folder; all data are invented teaching
# examples.
library("readr")
library("dplyr", warn.conflicts = FALSE)
library("tibble")
library("tidyr")
library("ggplot2")
library("lme4")
library("performance")
library("marginaleffects")
dir.create("outputs", showWarnings = FALSE)

# 1. Researchers assign one diet per pen and measure six pigs in each pen.
d <- readr::read_csv("pen-animals.csv", show_col_types = FALSE)
d |>
  select(room, pen, diet, gain_g_day) |>
  slice_head(n = 2)
d <- d |>
  mutate(
    room = factor(room, levels = paste0("R", 1:4)),
    pen = factor(pen),
    diet = factor(diet, levels = c("Control", "High fibre"))
  )
pen_design <- d |>
  distinct(room, pen, diet)
pen_design |> count(diet, name = "pens")
pen_design |> count(room, diet, name = "pens")
d |>
  count(pen, name = "animals")

# Show each pig before fitting: jitter changes horizontal positions only.
pig_plot <- ggplot(d, aes(pen, gain_g_day, colour = diet)) +
  geom_point(
    position = position_jitter(width = 0.15, height = 0, seed = 828),
    size = 2,
    alpha = 0.8
  ) +
  facet_wrap(~room, scales = "free_x", nrow = 2) +
  labs(x = "Pen", y = "Daily gain (g/day)", colour = "Diet") +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom")
pig_plot
ggsave("outputs/pig-gains-by-pen.png", pig_plot, width = 10, height = 4.2)

independent <- lm(
  gain_g_day ~ diet + room,
  data = d
)
shared_pen <- lme4::lmer(
  gain_g_day ~ diet + room + (1 | pen),
  data = d,
  REML = TRUE
)
summary(shared_pen)
# Formula: diet and room are fixed effects. (1 | pen) gives each pen a shared
# intercept difference from its diet-and-room mean. 1 means intercept, | means
# "for each group identified by", and pen is the grouping column.
# Pen differences have a common normal distribution centred on zero, with
# spread estimated from all pens. This is different from randomising diets.
pen_differences <- ranef(shared_pen)$pen |>
  rownames_to_column("pen")
pen_differences |>
  slice_head(n = 4)
# Positive means above that pen's diet-and-room expectation, in g/day.
# These are partly pooled fitted differences, not raw pen averages.

VarCorr(shared_pen)
check_singularity(shared_pen)
check_convergence(shared_pen)
avg_comparisons(independent, variables = "diet")
avg_comparisons(shared_pen, variables = "diet", re.form = NA)
# The approximate interval describes the estimated difference in mean gain,
# not the spread of individual animal gains. Only 16 pens were allocated diet.

# The same mean and varying-start idea is illustrated in random-effects-plots.R.
# Run that file separately for forty calf lines and the random-slope extension.

# 2. Split plot: diet goes to whole pens; bedding goes to sections within pens.
s <- read_csv("split-plot.csv", show_col_types = FALSE) |>
  mutate(
    room = factor(room, levels = paste0("R", 1:4)),
    pen = factor(pen),
    diet = factor(diet, levels = c("Control", "High fibre")),
    bedding = factor(bedding, levels = c("Standard", "Deep"))
  )
s |>
  select(pen, section, diet, bedding, gain_g_day) |>
  slice_head(n = 4)
s |>
  count(diet, bedding, name = "sections", .drop = FALSE)
split_fit <- lme4::lmer(
  gain_g_day ~ room + diet * bedding + (1 | pen),
  data = s,
  REML = TRUE
)
summary(split_fit)
VarCorr(split_fit)
# A * B includes the two separate changes and whether one changes the other.
# The balanced shortcut gives each room equal weight within each bedding.

avg_comparisons(
  split_fit,
  variables = "diet",
  newdata = "balanced",
  by = "bedding",
  re.form = NA
)

# Alternative: omit newdata to use every observed section with equal weight.
# Predict each row under both diets, keeping its room and bedding unchanged;
# subtract Control from High fibre, then average separately within bedding.
avg_comparisons(
  split_fit,
  variables = "diet",
  by = "bedding",
  re.form = NA
)
# No sections are missing: every room contributes four records per bedding,
# so both methods give each room equal weight and produce the same results.
# The model also gives the same diet difference in every room within bedding.

avg_comparisons(
  split_fit,
  variables = "bedding",
  newdata = "balanced",
  by = "diet",
  re.form = NA
)


avg_predictions(
  split_fit,
  newdata = "balanced",
  by = c("diet", "bedding"),
  re.form = NA
)

check_convergence(split_fit)
check_singularity(split_fit)
png("outputs/split-checks.png", width = 1200, height = 600, res = 130)
plot(check_model(split_fit, check = c("linearity", "qq")))
dev.off()
# Look for a curve or a changing spread in the differences between observations
# and fitted values. A reassuring picture does not prove the design assumptions.

# 3. Latin square: four feeds in sixteen broiler pens, four rows and columns.
# Each pen receives one feed throughout; one record is its mean daily gain
# in grams per bird per day. This is a spatial trial, not a crossover.
l <- read_csv("broiler-latin-square.csv", show_col_types = FALSE) |>
  mutate(
    row = factor(row, levels = paste0("Row", 1:4)),
    column = factor(column, levels = paste0("Col", 1:4)),
    feed = factor(feed, levels = c("A", "B", "C", "D"))
  )
l |>
  select(row, column, feed) |>
  pivot_wider(names_from = column, values_from = feed)
l |>
  count(row, feed, name = "pens", .drop = FALSE)
l |>
  count(column, feed, name = "pens", .drop = FALSE)
latin_fit <- lm(
  gain_g_day ~ feed + row + column,
  data = l
)
summary(latin_fit)

avg_comparisons(
  latin_fit,
  variables = list(feed = c("A", "C")),
  df = df.residual(latin_fit)
)

png("outputs/latin-checks.png", width = 1200, height = 600, res = 130)
plot(check_model(latin_fit, check = c("linearity", "qq")))
dev.off()
# Each pen counts equally; the additive contrast is the same at every
# position. No random effects. The t interval uses six error degrees of
# freedom and assumes independent normal errors with common variance.
# One square cannot estimate every feed-by-row and feed-by-column change.
# The greenhouse example and latin-square.csv remain in the backup slides.

# 4. Unequal group sizes after a measuring-equipment fault, unrelated to growth.
u <- read_csv("unequal-split-plot.csv", show_col_types = FALSE) |>
  mutate(
    room = factor(room, levels = levels(s$room)),
    pen = factor(pen),
    diet = factor(diet, levels = levels(s$diet)),
    bedding = factor(bedding, levels = levels(s$bedding))
  )
u |>
  count(diet, bedding, name = "sections", .drop = FALSE)
u |>
  distinct(room, pen, diet) |>
  count(room, diet, name = "pens", .drop = FALSE)
unequal_fit <- lme4::lmer(
  gain_g_day ~ room + diet * bedding + (1 | pen),
  data = u,
  REML = TRUE
)
u |>
  group_by(diet, bedding) |>
  summarise(mean_gain_g_day = mean(gain_g_day), .groups = "drop")
# Use the balanced shortcut again: each room keeps equal weight after
# observations are lost. re.form = NA excludes fitted pen differences.

avg_predictions(
  unequal_fit,
  newdata = "balanced",
  by = c("diet", "bedding"),
  re.form = NA
)


avg_comparisons(
  unequal_fit,
  variables = "diet",
  newdata = "balanced",
  by = "bedding",
  re.form = NA
)

check_convergence(unequal_fit)
# The model can handle unequal counts. It cannot fix losses caused by unrecorded
# poor health, or recover a treatment comparison absent from the design.

# Save the output you will interpret in the worksheet.
dir.create("outputs", showWarnings = FALSE)
avg_comparisons(
  split_fit,
  variables = "diet",
  newdata = "balanced",
  by = "bedding",
  re.form = NA
) |>
  write_csv("outputs/split-diet-differences.csv")
avg_predictions(
  unequal_fit,
  newdata = "balanced",
  by = c("diet", "bedding"),
  re.form = NA
) |>
  write_csv("outputs/unequal-adjusted-means.csv")

# 5. After the second break: three checks of the complete split-plot model.
# A residual is observed section gain minus fitted section gain, in g/day.
# Describe a sustained pattern in the green line, rather than expecting a
# perfectly flat line from 32 records.
png(
  "outputs/split-residuals-vs-fitted.png",
  width = 1000,
  height = 800,
  res = 140
)
plot(check_model(split_fit, check = "linearity"), n_columns = 1)
dev.off()

# These are remaining section residuals, not the calf random effects.
# Compare the middle and both ends with the normal reference line.
png("outputs/split-residual-qq.png", width = 1000, height = 800, res = 140)
plot(
  check_model(split_fit, check = "qq", detrend = FALSE),
  n_columns = 1
)
dev.off()

# Compare observed gains with 50 datasets simulated under the fitted model.
# re_formula = NA draws new random pen differences; section errors are also
# simulated. The room and treatment assignments stay the same.
set.seed(82850)
png("outputs/split-predictive-check.png", width = 1000, height = 800, res = 140)
plot(check_predictions(split_fit, iterations = 50, re_formula = NA))
dev.off()
# Describe agreement or differences in centre, spread and tails. These checks
# cannot establish randomisation, independent pens or no bedding spillover.