# Class 8 follow-up: analyse an individual animal's welfare-screen result.
# Run each section interactively from the supplied Class 8 project folder.
# Lines such as summary(binary) display their results in the Console.

# 1. Load the packages ------------------------------------------------------
# readr reads the CSV; dplyr prepares and counts the records.
# lme4 fits the model; performance checks it; marginaleffects calculates
# predicted probabilities and differences between diets.
library("readr")
library("dplyr", warn.conflicts = FALSE)
library("lme4", warn.conflicts = FALSE)
library("performance")
library("marginaleffects", warn.conflicts = FALSE)

# 2. Read and prepare the animal records ------------------------------------
# Each row is one animal. screened_positive is 1 for a positive welfare
# screen and 0 otherwise. These outcomes are invented for teaching.
# This is not the earlier pen-level outcome "any animal treated".
# Factors tell R that room, pen and diet identify named groups.
# Putting Control first makes it the reference diet for comparisons.
d <- read_csv("pen-animals.csv", show_col_types = FALSE) |>
  mutate(
    room = factor(room),
    pen = factor(pen),
    diet = factor(diet, levels = c("Control", "High fibre"))
  )
d |>
  count(diet, screened_positive)
# The table counts positive and negative animal screens under each diet.
# These animal counts are not counts of independent diet assignments:
# researchers assigned diets to whole pens, with six animals in each pen.

# 3. Fit a model for the yes/no outcome -------------------------------------
# glmer() fits a model with both named treatment differences and shared
# differences among groups. family = binomial() is for the 0/1 outcome.
# diet compares the two diets; room allows these four rooms to differ.
# (1 | pen) lets animals in the same pen share a pen difference.
binary <- lme4::glmer(
  screened_positive ~ diet + room + (1 | pen),
  family = binomial(),
  data = d
)
summary(binary)
# Read the full summary. The coefficients are on a log-odds scale, not a
# probability scale. The calls below express the diet results as probabilities.

# 4. Check the fit and the estimated pen spread -----------------------------
# TRUE from check_convergence() means the numerical fitting succeeded.
# This does not establish that the study assumptions are correct.
check_convergence(binary)

# TRUE from check_singularity() would flag an estimated group variation
# at or very near zero. FALSE means this check has not flagged the fit.
check_singularity(binary)

# VarCorr() displays the estimated spread of the shared pen differences.
# In this yes/no model, that pen spread is on the log-odds scale.
VarCorr(binary)
# The invented outcomes give a small non-zero fitted pen spread here.
# check_singularity() returns FALSE. Sixteen pens still give limited
# information about the spread of shared pen differences.

# 5. Predict the probability of a positive screen under each diet ------------
# newdata = "balanced" creates every combination of the model's categories.
# It gives each room equal weight, rather than weighting rooms by their
# number of recorded animals. This model has no numeric predictors.
# re.form = NA sets the fitted pen differences to zero in the predictions.
# by = "diet" returns one average probability for each diet.
avg_predictions(
  binary,
  newdata = "balanced",
  by = "diet",
  re.form = NA
)
# For a yes/no model, setting pen differences to zero generally gives a
# different answer from averaging probabilities across randomly varying pens.

# 6. Compare those probabilities -------------------------------------------
# variables = "diet" asks for High fibre minus Control. Both predictions
# keep the other settings the same; the function averages their differences
# over the balanced grid. A negative estimate means a lower predicted
# probability of a positive screen under High fibre.
avg_comparisons(
  binary,
  variables = "diet",
  newdata = "balanced",
  re.form = NA
)
# The estimate is a difference in probabilities, not an odds ratio.
# Multiply it by 100 to express the difference in percentage points.
# Std. Error and the approximate 95% interval describe uncertainty in the
# model's estimated comparison; they do not predict an individual animal's 0/1.