--- title: "From experimental designs to GLMs" subtitle: "Week 7: choosing a model and explaining its results" order: 17 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 --- ## What we will do today We will use the experiment from Weeks 5–6 to choose a statistical model, fit it in R, and explain what it says about the two diets. | Time | Work | |---|---| | 12:30–13:10 | Review the experiment; learn how to describe and fit a GLM | | 13:10–13:40 | Choose models for three short studies | | 13:40–13:55 | Break | | 13:55–14:55 | Compare diet effects in two housing conditions; analyse your group's data | | 14:55–15:10 | Break | | 15:10–16:00 | Apply the steps to a calf study; discuss projects; finish with a short question | ::: notes No new preparation is required. Explain the experiment again before asking students to recall terminology. All scripts and data are already in the Posit Cloud project. Use this deck until the first break. ::: ## The experiment we are analysing There are **64 pens in four rooms**, with six pigs in each pen. Within each room, eight pens receive Control feed and eight receive High fibre feed. Diet is randomly assigned to whole pens. We follow every pen for **28 days** and record one outcome per pen. The four groups analyse different outcomes from this same study design. This is a simulated teaching dataset, not evidence about the effects of a real diet. ::: notes The smaller experiment shown in earlier illustrations had 16 pens. State visibly that this dataset has 64 pens; do not ask students to infer the number from the old drawing. The pen is the independently treated unit because all six pigs share its assigned diet. ::: ## Opening recall: discuss the experiment Use the description on the previous slide. Discuss these questions with your group for two minutes: 1. Is diet assigned separately to each pig, or to the whole pen? 2. If the data contain one row per pen, how many rows should there be? 3. Why is it useful to have both diets represented in every room? Be ready to explain one answer in a sentence. ::: notes 12:30–12:40 including the experiment description. Expected answers: whole pens; 64 rows; comparing diets within rooms helps account for differences between rooms. Explain that the experimental unit is the unit independently assigned a treatment. Avoid asking students to name a family or link before these have been reviewed. ::: ## What does each group measure? The **response** is the outcome we want the model to explain. | Group | One measurement for each pen | |---|---| | 1 | Average daily weight gain, in grams per day | | 2 | Whether any pig needed treatment: 0 means none, 1 means at least one | | 3 | Number of treatment events during the 28 days; a pig may be treated more than once | | 4 | Feed conversion ratio: kg of feed consumed per kg of weight gained | For Group 2, one row describes the whole pen. It does not tell us how many individual pigs were treated. ## Four things to state when choosing a GLM A **generalised linear model (GLM)** describes how an expected outcome changes with explanatory variables, such as diet and room. State: 1. **Response:** what is measured, and what does one row describe? 2. **Formula:** which explanatory variables should be included? 3. **Family:** which probability distribution describes variation in the response? 4. **Link:** how is the expected response related to the sum of the model terms? The next slides explain these choices using daily weight gain. ## Reading the formula For our study, the formula is `response ~ diet + room`. - `response` is the outcome column in the data. - `~` separates the outcome from the explanatory variables. - `diet` lets the expected outcome differ between Control and High fibre. - `room` allows for differences between the four rooms. Including `room` lets us compare diets while accounting for where the pens were kept. ::: notes 12:40 onwards. This is a formula, not an instruction to add the raw data columns. The family and link are separate arguments to glm(). Use the same formula for all four groups but do not imply that it alone fully describes a GLM. ::: ## The linear predictor and the link The **linear predictor** is the intercept plus the contributions of diet and room. An intercept is the starting value for the reference categories: here, Control diet in Room 1. For daily gain we use an **identity link**. This means the linear predictor is already the expected gain in grams per day; no conversion is needed. Other links are useful for other outcomes. With a log link, use `exp()` to convert the linear predictor back to a positive expected value. With a logit link, use `plogis()` to convert it back to a probability. ::: notes Use the standard term linear predictor, not an invented replacement. Explain its meaning before any symbols. More detailed link calculations are in the after-class material. The identity-link example supplies everything needed for the worked Gaussian fit. ::: ## The families and links supplied for this activity | Outcome | Family | Link | |---|---|---| | Average daily gain, a continuous measurement | Gaussian, also called normal | Identity | | Whether any pig was treated, recorded as 0 or 1 | Binomial | Logit | | Number of treatment events | Poisson | Log | | Positive feed conversion ratio | Gamma | Log | These are starting assumptions for the teaching examples. A suitable outcome type does not, by itself, establish that a model fits the data well. ## Open the daily-gain example in Posit Cloud Open `group-1-student.R` in the supplied Week 7 project. Run the first section. `source()` reads the supplied file of functions. `read_pen_data(1)` reads Group 1's CSV into a data frame called `d` and sets Control and Room 1 as the reference categories. ```{r} #| label: load-pen-data source("glm-helpers.R") d <- read_pen_data(1) head(d[c("pen_id", "diet", "room", "response")], 2) ``` Each displayed row describes one pen; `response` is that pen's average daily gain. ## Fit the model for daily gain `glm()` estimates the model from `d`. We show the first two coefficient rows and only their estimates and standard errors. ```{r} #| label: fit-gain fit <- glm(response ~ diet + room, family = gaussian("identity"), data = d) # Intercept and diet estimate/SE; remaining rows are room terms. coef(summary(fit))[1:2, c("Estimate", "Std. Error")] ``` `(Intercept)` is the estimated average gain for Control pens in Room 1. `dietHigh fibre` is the estimated High fibre minus Control difference, in grams per day. ## How to read the coefficient table The **Estimate** column contains the fitted coefficient values. A positive diet estimate means the model predicts higher average gain with High fibre than with Control, within the same room. The **Std. Error** column describes uncertainty in each coefficient estimate. It is not the spread of individual pen measurements. Look back at the diet row. Does its estimate describe an increase or a decrease in expected gain? What are its units? ::: notes Allow students to identify the positive estimate and g/day before moving on. The fitted difference is about 57.6 g/day. Keep the distinction between the assumed effect used to simulate data and the estimate fitted to this dataset. Do not use this synthetic result as empirical evidence about diets. ::: ## Choose the conditions for our predictions We want to predict the expected gain for both diets **in the same room**. The supplied `prediction_grid(d)` function makes a table containing both diets in each room. It does not contain new observations: it lists the conditions for which we want predictions. ```{r} #| label: prediction-grid grid <- prediction_grid(d) head(grid, 2) ``` The first two rows ask for Control and High fibre predictions in Room 1. ## Read the predictions for Room 1 `predict_means(fit, grid)` calculates the expected outcomes and approximate 95% confidence intervals. ```{r} #| label: gain-predictions predictions <- predict_means(fit, grid) subset(predictions, room == "Room 1", select = c(diet, room, expected, lower, upper)) ``` `expected` is estimated average gain, in g/day. `lower` and `upper` show uncertainty in that average, assuming the model is appropriate. They do not describe the range containing 95% of individual pens. ## Compare the predictions in all four rooms `response_spec(1)` supplies Group 1's labels and settings, stored in `spec`. The plotting function uses these to label gain in g/day. Its points show estimated averages; its bars show their confidence intervals. ```{r} #| label: gain-prediction-plot #| fig-width: 10 #| fig-height: 3.4 #| out-width: "90%" spec <- response_spec(1) plot_predictions(predictions, spec) ``` In Room 1, which diet has the higher estimated average gain? Check that your answer agrees with the diet coefficient. ::: notes The supplied response_spec(1) function provides Group 1 labels and family/link settings. The plot uses its units. Compare the two points within one room, then show that an additive identity-link model uses the same gain difference across rooms. Do not ask students to interpret intervals as model checks. ::: ## Describe the diet comparison in the correct units | Group | What the model's diet coefficient tells us | |---|---| | 1 | Difference in expected gain, in g/day | | 3 and 4 | Exponentiating the coefficient gives the High fibre / Control ratio of expected outcomes | | 2 | Exponentiating the coefficient gives a ratio of odds; this needs a separate explanation | For example, an expected-count ratio of 0.60 means the expected count under High fibre is 60% of the Control count in the same room. ## For the binary outcome: probability and odds A probability of 0.30 means that about 30 out of 100 comparable pens would have at least one pig treated. **Odds** are probability divided by one minus probability: here, `0.30 / 0.70`. An **odds ratio** compares those odds between the two diets. An odds ratio of 0.60 does not mean the probability is multiplied by 0.60. Group 2 should report the model's predicted probability for each diet; the after-class example shows the odds calculation step by step. ## Choose models for the short studies Open the [class worksheet](student-activity.html). For Studies A and B, write down what is measured, which unit receives diet, the formula, and the supplied family and link. For Study C, identify the four diet-and-housing combinations and the averages we need to compare. We will explain the interaction formula after the break. Use the response table and the worked daily-gain example as references while answering. ::: notes 13:10–13:40. Give 10 minutes each to A and B, then 7 minutes identifying the comparisons for C and 3 minutes to compare answers. Do not require an unexplained interaction formula before teaching it. Students may ask for the recap; participation does not depend on completing earlier after-class work. ::: ## Before the break Keep your written answers. If your group is unsure about a choice, write the question beside it so we can discuss it. After the break we will look at a study with both diet and housing treatments, then analyse each group's pen data.