# Week 5, Group 4: Feed conversion ratio (kg feed/kg gain)
# Run section by section in Posit Cloud. Complete the four TODO expressions.
# Packages: tibble constructs the design table; dplyr supplies the familiar
# data manipulation verbs; ggplot2 compares the simulated datasets visually.
# R's supplied stats package provides the generators and glm().
library("tibble")
library("dplyr")
library("ggplot2")

# 1. Design: one row per pen, six pigs per pen, equal follow-up.
pens <- tibble(
  pen_id = sprintf("P%02d", 1:16),
  room = rep(paste0("Room ", 1:4), each = 4),
  position = rep(1:4, times = 4)
)
labels <- rep(c("Control", "High fibre"), each = 2)
set.seed(505)
design <- pens |>
  group_by(NA) |> # TODO 1: replace NA with the blocking variable
  mutate(diet = NA_character_) |> # TODO 2: shuffle all labels without replacement
  ungroup()
if (anyNA(design$diet)) stop("Complete the treatment shuffle in TODO 2.")
if (!"room" %in% names(design)) stop("Keep the room column.")
design <- design |>
  select(pen_id, room, position, diet) |>
  mutate(diet = factor(diet, levels = c("Control", "High fibre")),
         room = factor(room, levels = paste0("Room ", 1:4)))
allocation_check <- design |> count(room, diet, .drop = FALSE)
print(allocation_check)
stopifnot(nrow(design) == 16, n_distinct(design$pen_id) == 16,
          nrow(allocation_check) == 8, all(allocation_check$n == 2))
# Explain why passing these checks matters:
# ...

# 2. Assumptions supplied for your scenario, not estimated from real data.
# The intercept refers to CONTROL IN ROOM 1 (the first factor level).
# Room effects are contrasts with Room 1, on the link scale.
beta_0 <- log(2.80)
beta_diet <- log(0.92)
room_effects <- c(0, -0.06, -0.10, -0.04)
model_family <- Gamma(link = "log")
# Response support:
# Experimental unit:
# Meaning of beta_diet, and its response-scale interpretation:
# What variation remains after accounting for diet and room?
# One assumption to question:

# 3. This supplied function repeats your recipe when given a response seed.
# Complete its two model lines, then rerun the function definition.
generate_response <- function(response_seed) {
  set.seed(response_seed)
  simulated <- design |>
    mutate(
      room_effect = room_effects[as.integer(room)],
      eta = NA_real_, # TODO 3: complete the model score
      mu = NA_real_ # TODO 4: use your recipe's inverse link
    )
  if (anyNA(simulated$eta) || anyNA(simulated$mu)) {
    stop("Complete eta and mu in TODOs 3 and 4, then rerun the function.")
  }
  stopifnot(all(is.finite(simulated$eta)), all(is.finite(simulated$mu)))
  simulated |>
    mutate(response = rgamma(n(), shape = 25, scale = mu / 25))
}
first <- generate_response(1504)
print(first |> select(pen_id, room, diet, eta, mu, response))
stopifnot(nrow(first) == 16, all(first$response > 0))
range(first$mu)
range(first$response)

# 4. Predict which columns change. Rotate roles before running this section.
second <- generate_response(2504)
stopifnot(identical(first$pen_id, second$pen_id),
          identical(first$diet, second$diet),
          identical(first$eta, second$eta), identical(first$mu, second$mu))
comparison <- first |>
  select(pen_id, room, diet, mu, first_response = response) |>
  mutate(second_response = second$response)
print(comparison)
plot_data <- bind_rows(first = first, second = second, .id = "dataset")
ggplot(plot_data, aes(diet, response, colour = room)) +
  geom_point(position = position_jitter(width = 0.08, height = 0, seed = 505),
             size = 3) +
  facet_wrap(~dataset) +
  labs(x = NULL, y = "Feed conversion ratio (kg feed/kg gain)") +
  theme_minimal(base_size = 14)
# Which columns stayed fixed? Did every response have to change?
# ...

# 5. Supplied analysis. Rotate roles again.
# Keep warnings and extreme estimates in the record, especially for Group 2.
fit_first <- glm(response ~ diet + room, family = model_family, data = first)
fit_second <- glm(response ~ diet + room, family = model_family, data = second)
summary(fit_first)
fit_first$converged
fit_first$boundary
b <- c(first = unname(coef(fit_first)["dietHigh fibre"]),
       second = unname(coef(fit_second)["dietHigh fibre"]))
print(b)
print(exp(b)) # ratio of expected feed conversion ratios
# Reference-room predictions, on the scale of the response:
new_pens <- data.frame(
  diet = factor(c("Control", "High fibre"), levels = levels(design$diet)),
  room = factor(rep("Room 1", 2), levels = levels(design$room))
)
predict(fit_first, newdata = new_pens, type = "response")
# Record generating effect, fitted effects, interpretation, warnings:
# ...

# 6. Save your two datasets and script. Use your group folder to avoid clashes.
dir.create("group-4-results", showWarnings = FALSE)
write.csv(first, "group-4-results/first.csv", row.names = FALSE)
write.csv(second, "group-4-results/second.csv", row.names = FALSE)
# Peer audit: one defended assumption and one sensitivity question:
# ...
# Individual exit answer belongs in the worksheet.