From experimental designs to GLMs

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.

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.

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.

What does each group measure?

The response is the outcome we want the model to explain.

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.

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.

The families and links supplied for this activity

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.

source("glm-helpers.R")
d <- read_pen_data(1)
head(d[c("pen_id", "diet", "room", "response")], 2)

  pen_id    diet   room response
1   P001 Control Room 1 695.1326
2   P002 Control Room 1 688.3944

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.

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")]

                Estimate Std. Error
(Intercept)    689.98514   16.46264
dietHigh fibre  57.56809   14.72463

(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?

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.

grid <- prediction_grid(d)
head(grid, 2)

        diet   room
1    Control Room 1
2 High fibre Room 1

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.

predictions <- predict_means(fit, grid)
subset(predictions, room == "Room 1",
       select = c(diet, room, expected, lower, upper))

        diet   room expected    lower    upper
1    Control Room 1 689.9851 657.0435 722.9268
2 High fibre Room 1 747.5532 714.6116 780.4949

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.

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.

Describe the diet comparison in the correct units

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. 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.

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.