--- title: "Fit and explain a GLM" subtitle: "Week 7: compare diets, explain results and check answers" order: 18 html-math-method: katex format: aarhus-revealjs: embed-resources: true css: class-5-simulation.css transition: none incremental: false code-line-numbers: false execute: echo: true warning: true message: false freeze: false eval: true author: - name: Gavin Simpson affiliations: Aarhus University date: "2026-10-21" title-slide-attributes: data-background-color: "#003d73" presenter: name: Gavin Simpson institute: Department of Animal & Veterinary Sciences --- ## Does the diet difference change with housing? This is a **separate study** with 48 pens in four rooms. Pens receive one of two diets and one of two housing conditions, Standard or Enriched. Every room contains all four combinations, with three pens per combination. Both treatments are assigned to pens. We ask whether the High fibre minus Control difference in average gain is the same in both housing conditions. ::: notes 13:55–13:58. Recreate the four combinations with labelled Lego pens. Explain that a factorial design is one in which treatment factors are combined; it does not require one particular analysis model. Use the same example throughout this section. Do not add housing to the diet-only group datasets. ::: ## Read the data for the diet-and-housing study The supplied `read_factorial_data()` function reads `data/factorial-pens.csv` into `f`. Its reference categories are Control diet, Standard housing and Room 1. ```{r} #| label: load-factorial source("glm-helpers.R") f <- read_factorial_data() head(f[c("diet", "housing", "room", "gain")], 2) ``` `gain` is the pen's average daily gain in g/day. `f` is different from the group datasets we analyse later. ## Two ways to describe the diet difference An **additive model** requires the same High fibre minus Control gain difference in Standard and Enriched housing. A model with an **interaction** allows that difference to change with housing. In R, `diet * housing` includes both variables and their interaction. ```{r} #| label: factorial-formula additive_formula <- gain ~ diet + housing + room interaction_formula <- gain ~ diet * housing + room additive_formula interaction_formula ``` ::: notes 13:58–14:03. Explain the two assumptions before asking students to interpret the plot. Both models retain room. In this Gaussian identity-link example, the comparison is directly in g/day. Detailed interpretations under other links are not required now. ::: ## Fit the two models to the same observations We use the same gain data, rooms and Gaussian identity family/link for both models. Only the treatment terms change. ```{r} #| label: fit-factorial additive <- glm(additive_formula, data = f, family = gaussian("identity")) interaction <- glm(interaction_formula, data = f, family = gaussian("identity")) # The two diet-related coefficients; other terms are omitted here. coef(interaction)[ c("dietHigh fibre", "dietHigh fibre:housingEnriched")] ``` In the interaction model, `dietHigh fibre` describes the diet difference in Standard housing. The `dietHigh fibre:housingEnriched` coefficient tells us how much larger or smaller that difference is in Enriched housing. ## List the four combinations to compare To compare diets fairly, keep room fixed. We choose Room 1 and list both diets under both housing conditions. ```{r} #| label: factorial-grid factorial_grid <- expand.grid( diet = levels(f$diet), housing = levels(f$housing), room = "Room 1") factorial_grid ``` The four rows describe the conditions for predictions, not four new measurements. ## Calculate the predicted average gains `predict_means()` uses each fitted model and the four combinations above. We label the results so the plot can show both models. The printed table shows the interaction model's expected gains, in g/day. ```{r} #| label: factorial-predictions pred_add <- predict_means(additive, factorial_grid) pred_add$model <- "Additive" pred_int <- predict_means(interaction, factorial_grid) pred_int$model <- "Interaction" factorial_predictions <- rbind(pred_add, pred_int) pred_int[c("diet", "housing", "expected")] ``` ## Read the two plots In each panel, compare Control and High fibre within Standard housing, then make the same comparison within Enriched housing. ```{r} #| label: factorial-plot #| fig-width: 8 #| fig-height: 3.7 #| out-width: "85%" plot_factorial_predictions(factorial_predictions) ``` Which model gives the same increase from Control to High fibre in both housing conditions? ::: notes Give students time to point to the two comparisons. The additive panel forces the same increase. The interaction panel allows different increases. The bars show uncertainty in the estimated averages; a difference between drawn lines is not, by itself, evidence that a practical diet recommendation is justified. ::: ## Calculate how the diet difference changes First subtract Control from High fibre within each housing condition. Then subtract the Standard-housing difference from the Enriched-housing difference. ```{r} #| label: factorial-contrasts diet_difference <- with(pred_int, expected[diet == "High fibre"] - expected[diet == "Control"]) names(diet_difference) <- levels(f$housing) diet_difference c(difference_of_differences = unname(diff(diet_difference))) ``` The final number is in g/day. A positive value means the estimated diet difference is larger in Enriched housing. ::: notes 14:03–14:10. The fitted differences are about 43.0 and 78.7 g/day; the increase in the diet difference is about 35.7 g/day. Ask students to say what was subtracted from what. Do not require memorisation of the phrase difference of differences. Update Study C together; optional factorial-student.R reproduces these calculations. ::: ## Now open your group's script Use `group-1-student.R`, `group-2-student.R`, `group-3-student.R` or `group-4-student.R` in the supplied Posit Cloud project. One person runs R, a second explains each result, and a third checks the explanation against the study description. Swap jobs at the next break. Run sections 1 and 2 first. Check that there are eight pens on each diet in every room. All groups fit `response ~ diet + room`; their measured outcomes and families differ. ## Worked example: treatment-event counts Here we demonstrate Group 3. `response_spec(3)` supplies its family, link and labels; `spec$family` selects the family/link settings from that list. `read_pen_data(3)` reads the counts into `d`. ```{r} #| label: fit-counts spec <- response_spec(3) # worked example: Group 3 # These d rows count treatment events per pen. d <- read_pen_data(3) fit <- glm(response ~ diet + room, family = spec$family, data = d) coef(summary(fit))[1:2, c("Estimate", "Std. Error")] ``` With a log link, the printed diet estimate is a log ratio. The next step converts it to an easier-to-read ratio and shows expected counts. ::: notes 14:10–14:20. Other groups run their own scripts rather than changing to Group 3. Counts refer to events per pen during 28 days. A pig can contribute more than one event. The Gaussian group reads a difference directly; Group 2 uses probabilities plus an odds ratio; Group 4 compares expected feed conversion ratios. ::: ## Expected counts under the two diets The first table below gives expected treatment-event counts for pens in Room 1 over 28 days. The second gives the High fibre / Control expected-count ratio and its approximate 95% confidence interval. ```{r} #| label: count-predictions grid <- prediction_grid(d) predictions <- predict_means(fit, grid) subset(predictions, room == "Room 1", select = c(diet, expected, lower, upper)) diet_comparison(fit) ``` `estimate` in the second table is the High fibre / Control expected-count ratio; `lower` and `upper` are its confidence limits. Which diet has the lower expected count? Explain what the Control count is being multiplied by. ## Write an explanation of your own results Use Room 1 for your first explanation. Write: 1. What you measured, including units, and that diet was assigned to pens. 2. Which response, explanatory variables, family and link you used. 3. The estimated outcome for each diet, and the difference or ratio between them. 4. The interval showing uncertainty, and one assumption that still needs checking. Use the worksheet's examples and recap if you need help with the wording. ::: notes 14:20–14:55. Groups run predictions, write and discuss explanations, then compare with another group. Group 2 must say probability of at least one pig treated in a pen, not individual pig probability. Confidence intervals concern estimated averages/probabilities under model assumptions. A successful fit does not show that those assumptions are appropriate. ::: ## Apply the same steps to a calf study There are 40 pens in five barns, with four calves in every pen. In each barn, four pens receive Standard feed and four receive Supplemented feed, assigned at random. The outcome `any_treated` is 1 if at least one calf in a pen needed treatment during 21 days, and 0 otherwise. Use the worksheet to describe a model comparing the two feeds while allowing for differences between barns. No new dataset or R fit is needed for this question. ::: notes 15:10–15:20. Read the complete study description before asking for a formula. Students need to identify the pen as the treatment unit and barn as a predictor. They have already been given the binary-outcome family/link explanation. Do not ask for an animal-level model. ::: ## Check another group's answer Exchange answers. For each point below, explain whether you agree and why: - Does one row describe a pen, with feed assigned to that pen? - Is the outcome correctly described as 0 or 1 for any calf treated? - Does the formula include both feed and barn? - Are the family and link stated? - Do the proposed predictions compare both feeds in the same barn? Give the group one specific suggestion, then revise your own answer. ::: notes 15:20–15:30. Five minutes checking and five revising. This replaces the old 35-minute instructions that did not fit the 20-minute slot. Help students connect an identified error to its consequence; do not ask them merely to label an answer supported or unsupported. ::: ## Discuss your group project Bring out the short project description from the week 42 task on the course website. Tell us the question you want to answer, what you will measure, and which groups or conditions you will compare. Which dataset will you use, and can you access it? Write down any information you still need about how the data were collected. If your preferred data are unavailable, choose one of the provided studies on the course website to work with. ::: notes 15:30–15:45. Allow three minutes per group and three minutes to summarise next steps. Do not require a sophisticated fitted model today. Record who will obtain missing information about treatment assignment or sampling before 28 October. ::: ## A short question before you leave For the pen study, a student writes only `response ~ diet + room` as their description of the model. What else must they state about the measured outcome, the family and the link so that someone else could fit the same model? Then use your own output to write one sentence comparing the two diets in Room 1. Include units, and mention uncertainty. ::: notes 15:45–16:00 includes discussion and after-class instructions. Accept clear explanations in everyday words. Check that students know which outcome their group measured, why the pen is the unit assigned diet, and what the model predicts. Introduce the animals-within-pens drawing for next week without requiring mixed-model notation. :::