GLMs and modelling

After Week 5: explanation, consolidation and preparation for power

Gavin Simpson

gavin@anivet.au.dk

Aarhus University

Wednesday, 21 October 2026

How to use this recording

This is the follow-up to the simulation workshop. Keep your group’s script and the after-class worksheet beside you.

Chapter	Focus	Recording time
1	The model behind the code	About 10 minutes
2	Four data generating processes	About 15 minutes
3	Fitting and repeated experiments	About 10 minutes

Pause to answer the questions. The worksheet adds about 25 minutes of practice.

The model behind the code

Chapter 1: about 10 minutes

The class model

We kept the treatment allocation and numerical assumptions fixed, then drew two sets of responses.

The GLM described the expected response for each pen and the variation remaining among comparable pens.

The fitted coefficients changed because each dataset contained different information about the same generating model.

The generative sequence

Sketch connecting an experimental design to a model calculation and a random observation.

A conditional model

For a pen with a particular diet and room, write:

E(Y_i \mid \text{diet}_i,\text{room}_i)=\mu_i

The family describes responses given those predictors. Pooling pens with different means can produce a distribution with two peaks.

Histogram of 2,000 simulated responses pooled from two overlapping Gaussian groups, with peaks near 640 and 800 g per day.

This histogram pools one simulated dataset: 1,000 independent Gaussian observations with mean 640 g/day and 1,000 with mean 800 g/day, both with SD 55 g/day. These illustrative parameters differ from the workshop parameters. The two groups overlap, but pooling them produces two broad peaks. The fixed seed makes the example reproducible, and the histogram uses 25 g/day bins. Ask students whether this raw response histogram rules out a Gaussian conditional model. Each group’s generating distribution is Gaussian even though the pooled observations have a different shape. Pooling does not always produce multiple peaks: the separation, spread and mixture weights matter, and histogram appearance also depends on sample size and bin width. A Gaussian GLM concerns conditional responses or errors around fitted means, not normality of the raw response histogram. Our scenario also assumes independent responses conditional on diet and the specified room effects. Fixed room effects describe mean differences, but do not automatically model residual dependence within rooms.

Components of a GLM

\eta_i=\boldsymbol{x}_i^\mathsf{T}\boldsymbol{\beta}, \qquad g(\mu_i)=\eta_i, \qquad Y_i \mid \boldsymbol{x}_i \sim F_i

The linear predictor combines treatment and room terms. The link connects the expected response to that predictor. The family determines response support and a mean–variance relationship.

For simulation, we must also specify any remaining scale or dispersion parameters.

\eta — “eta”, the values on the linear predictor
\mu — “mu”, the expected value or mean of the response on the response scale
*g — the link function

Link and inverse link

Family used in class	Link, from mean to score	Inverse link, from score to mean
Gaussian	\eta=\mu	\mu=\eta
Binomial	\eta=\log\{p/(1-p)\}	p=1/(1+\exp(-\eta))
Poisson and Gamma	\eta=\log(\mu)	\mu=\exp(\eta)

R functions for the logit and inverse logit

R uses

qlogis(p) for the logit, i.e. \eta=\log\{p/(1-p)\}, and
plogis(eta) for its inverse, i.e. p=1/(1+\exp(-\eta))

These are the links chosen for our scenarios, but other links are available.

Notice that we never transform the data themselves!

Room effects and the intercept

Our generator uses:

\eta_i=b_0+b_\mathrm{diet}I(\text{high fibre})+\gamma_{\mathrm{room}(i)}

Room 1 has \gamma_1=0, so b_0 describes its control mean on the link scale.

Room 1 is also the first factor level in our fitted model. With treatment contrasts, the fitted intercept estimates b_0, and each room coefficient compares that room with Room 1.

Fixed room means and independence

In the workshop, each room effect is a fixed number. We draw responses independently after setting those means.

Shared room conditions could create additional dependence, and diet effects could differ among rooms. Either would require a richer generator and an appropriate analysis.

Simulation forces us to state which variation we include and which questions remain open.

Pause: model anatomy

In your script, locate:

the allocation
the numerical model parameters
the linear predictor and inverse link
the random response generator

Explain which would change if you changed only the response seed. Then explain what would change if you increased the treatment effect.

Four data generating processes

Chapter 2: about 15 minutes

Gaussian pen means

Y_i\mid\boldsymbol{x}_i \sim N(\mu_i,\sigma^2), \qquad \mu_i=\eta_i

In our scenario, \sigma=60 g/day and the diet changes expected daily gain by +45 g/day.

The Gaussian model assigns probability across the real line. Negative daily gain can make sense, but extreme tails may still be biologically implausible.

The assumed SD describes variation among pen means after accounting for diet and room.

Binary pens and baseline probability

Y_i\sim\operatorname{Bernoulli}(p_i), \qquad \operatorname{Var}(Y_i\mid\boldsymbol{x}_i)=p_i(1-p_i)

The response records whether any pig in a pen required treatment. There is one binary trial per pen, even though each pen contains six pigs.

A fixed odds ratio can imply different absolute probability changes at different baseline probabilities.

Odds and probabilities

In Room 1, control probability is 0.30, giving odds 0.30/0.70=0.4286.

Multiplying those odds by 0.60 gives high-fibre odds of 0.2571. Converting back gives:

p_\mathrm{high}=\frac{0.2571}{1+0.2571}\approx0.2045

This is about a 9.55 percentage-point decrease, not a probability of 0.18 (0.3 / 0.6). The odds ratio and probability ratio are different quantities.

The event here is that at least one pig in a pen requires veterinary treatment during the study. The probability is for that pen-level event, not for an individual pig.

Probability versus odds. A probability compares events with all possible outcomes. A control probability of 0.30 means that, across many comparable pens, we expect about 30 in every 100 to have the event and 70 not to. Odds compare events with non-events: 30/70 = 0.30/(1 - 0.30) = 0.4286. This is about 3 to 7, not a probability of 42.86%. Probabilities lie between 0 and 1, whereas odds can range from 0 to infinity.

Apply the treatment effect on the correct scale. Our binomial GLM uses a logit link, so eta = log(p / (1 - p)). The diet coefficient is log(0.60), approximately -0.5108. Adding this coefficient to the control log-odds multiplies the odds by exp(log(0.60)) = 0.60. Thus the high-fibre odds are 0.4286 * 0.60 = 0.2571. An odds ratio of 0.60 means 40% lower odds, comparing diets within the same room.

Convert back to probability. If o = p / (1 - p), then o(1 - p) = p, so p = o / (1 + o). Hence 0.2571 / 1.2571 is approximately 0.2045. Use unrounded values in calculations. Across many comparable high-fibre pens, we would expect about 20.45 in every 100 to have the event. This is an expectation, not a claim that the small class dataset will have that exact proportion.

Three summaries of the same contrast. The odds ratio is 0.60, a 40% reduction in odds. The probability ratio is 0.2045 / 0.30 = 0.6818, approximately a 31.8% relative reduction in probability. The absolute probability change is 0.2045 - 0.30 = -0.09545, a decrease of about 9.55 percentage points. Distinguish “percentage points” from a relative percentage reduction. Multiplying 0.30 by 0.60 to get 0.18 would treat 0.60 as a probability ratio, which is not our model assumption.

Connect to the R code. Calculate plogis(qlogis(0.30) + log(0.60)). qlogis() converts the control probability to log-odds, adding log(0.60) applies the diet effect, and plogis() converts the result back to probability. Equivalently, for control probability p and odds ratio OR, the high-fibre probability is OR * p / (1 - p + OR * p).

Both the generator and fitted model use Room 1 as the reference. Its control probability is 0.30, supplied by beta_0 = qlogis(0.30).

Transition to the next slide. Ask: “If the control probability changes but the odds ratio stays at 0.60, will the probability fall by the same number of percentage points?” No. For example, a control probability of 0.15 gives plogis(qlogis(0.15) + log(0.60)) = 0.09574, a decrease of about 5.43 percentage points. The next slide shows this dependence on the starting probability.

A constant odds ratio

The curve uses odds ratio 0.60 throughout. The dashed line marks equal probabilities.

Orient the audience. The horizontal axis is the control probability of at least one pig in a pen requiring treatment. The vertical axis is the corresponding high-fibre probability for comparable pens in the same room. Each point on the blue curve is a different possible starting probability, with the diet odds ratio held at 0.60. This is a calculated relationship, not simulated observations, fitted data or a power curve.

Read the reference line. The dashed diagonal is y = x: equal probabilities under both diets, equivalent to an odds ratio of 1 for probabilities strictly between 0 and 1. The blue curve lies below it because an odds ratio of 0.60 gives a lower high-fibre probability. At a chosen control probability, move vertically to the blue curve and read across to the vertical axis.

Explain the pink point. It marks Room 1’s teaching assumption: control probability 0.30 and high-fibre probability about 0.2045. The vertical distance from the blue curve to the dashed line at x = 0.30 is 0.30 - 0.2045 = 0.09545, or about 9.55 percentage points. Room 1 is the reference level in the generator and fitted model.

Why the curve bends. A constant odds ratio means a constant difference on the log-odds scale: logit(p_high) = logit(p_control) + log(0.60). The inverse logit is nonlinear, so that fixed difference does not translate into a fixed probability difference or a fixed probability ratio. In the code, qlogis(p) computes control log-odds, log(.60) is the diet effect, and plogis() converts their sum back to probability.

Walk through contrasting examples. With control probability 0.15, the high-fibre probability is approximately 0.0957, a decrease of 5.43 percentage points. At 0.30, it is 0.2045, a decrease of 9.55 points. At 0.50, it is 0.375, a decrease of 12.5 points. At 0.90, it is approximately 0.8438, a decrease of 5.63 points. The absolute decrease therefore does not simply keep increasing with baseline probability: it is small near either endpoint and larger in the middle. All these examples use exactly the same odds ratio.

Distinguish alternative effect assumptions. A constant probability ratio of 0.60 would give the straight line p_high = 0.60 * p_control. That is not the blue curve. A constant absolute decrease would give a line parallel to the diagonal over its valid range. Our model assumes neither of those. In a logit model without a diet-by-room interaction, the odds ratio is common across rooms even though their probability differences can differ.

Check understanding and connect to simulation. Ask: “If the control probability is 0.50, should the high-fibre probability be 0.30?” No: the control odds are 1, the high-fibre odds are 0.60, and the probability is 0.60 / 1.60 = 0.375. Then ask what should remain fixed when exploring a different baseline in the simulator: keep beta_diet = log(0.60), change the intercept, and recompute the probabilities. Baseline probability also changes Bernoulli variation p(1 - p), which is one reason Week 6 power simulations must specify baseline assumptions as well as the treatment effect. The curve alone does not establish which scenario has higher power.

The plot uses stats and graphics supplied with R; no additional package is required.

Poisson counts and exposure

Y_i\sim\operatorname{Poisson}(\mu_i), \qquad \operatorname{Var}(Y_i\mid\boldsymbol{x}_i)=\mu_i

Our expected-count ratio is 0.65. It is also a rate ratio because observation time and pen size are equal.

If exposure differed, a count model could use \log(\text{exposure}_i) as an offset when the scientific assumptions justify a rate model.

Fractional means, integer observations

The dashed line marks the mean, 1.95. The response places probability only on integer counts.

Extra variation in counts

Pens with the same predicted mean may vary more than a Poisson model allows.

Possible causes include unmodelled heterogeneity, dependent events or a different biological mechanism. A negative-binomial generator is one possible sensitivity scenario.

Choosing a more variable generator also requires checking the planned analysis. A Poisson fit may understate uncertainty under extra variation.

Gamma means and spread

R’s rgamma() uses a shape k and scale \theta parameterisation for a Gamma draw:

E(Y)=k\theta,\qquad \operatorname{Var}(Y)=k\theta^2

Setting \theta=\mu/k gives mean \mu, variance \mu^2/k and coefficient of variation 1/\sqrt{k}.

Our shape is k=25, so the coefficient of variation is 0.20: 1 / \sqrt{25}. The log-link treatment coefficient gives a ratio of expected feed-conversion-ratio values.

Introduce the reason for the comparison before doing the arithmetic: “Suppose we are uncertain about the assumed variability in FCR. As a sensitivity check, we could increase its coefficient of variation from 20% to about 33%, while keeping the expected FCR unchanged. This corresponds to reducing the Gamma shape from 25 to 9. How must we adjust scale?”

Shape 9 is an illustrative sensitivity scenario from the Group 4 after-class exercise, not a required value or an estimate supported by data. We choose it because CV = 1 / sqrt(k): shape 25 gives CV = 1/5, while shape 9 gives CV = 1/3. These make the comparison easy to explain. Another smaller positive shape could illustrate greater relative variation at the same mean. The aim is to distinguish uncertainty about the expected response from uncertainty about variation around that expectation.

The Gamma mean is the product of shape and scale: E(Y) = k * theta. Shape therefore does not control only spread when scale is held fixed. If we change k without adjusting theta, we also change the mean.

To isolate a change in variation, first keep the model’s expected response mu fixed, then calculate scale = mu / shape. Reducing shape requires increasing scale so that their product remains mu. With this adjustment, Var(Y) = mu^2 / k, SD(Y) = mu / sqrt(k), and CV = SD / mean = 1 / sqrt(k). Smaller shape gives more relative variation and greater right-skewness.

For the Room 1 control example, mu = 2.80 kg feed/kg gain. With shape 25, scale is 2.80 / 25 = 0.112, giving mean 25 * 0.112 = 2.80 and SD 2.80 / sqrt(25) = 0.56. If we reduce shape to 9 but leave scale at 0.112, the mean becomes 9 * 0.112 = 1.008. That would change expected FCR as well as its variation, so it would not answer our intended sensitivity question.

To implement this hypothetical increase in variability while preserving the mean, use shape 9 together with scale 2.80 / 9, approximately 0.3111. The mean stays at 2.80, while the SD rises to 2.80 / 3, approximately 0.9333. The CV increases from 0.20 to one third. The next slide shows these two distributions with the same mean.

In the activity code, change rgamma(n(), shape = 25, scale = mu / 25) to rgamma(n(), shape = 9, scale = mu / 9). A useful way to avoid mismatching the two values is to define gamma_shape <- 9 and use rgamma(n(), shape = gamma_shape, scale = mu / gamma_shape). Each pen has its own mu from diet and room, so calculate its scale from that mu. Keeping the mean fixed refers to each conditional generating mean, not to the realised sample mean, which still varies across simulated datasets.

Ask: “If we halve shape and want the same mean, what must happen to scale?” It must double. The conditional variance then doubles, and the SD increases by a factor of sqrt(2).

FCR means feed conversion ratio: pen feed intake divided by pen mass gain, in kg/kg. This scenario assumes positive intake and gain so FCR is positive and defined. R also accepts a rate parameter, equal to 1 / scale; our code explicitly uses scale throughout. R Gamma documentation: https://stat.ethz.ch/R-manual/R-patched/library/stats/html/GammaDist.html.

Same Gamma mean, different spread

Both means are 2.80. Reducing shape increases relative variation.

Conditional variation summary

Response	Conditional variance	Extra simulation parameter
Gaussian	\sigma^2	SD \sigma=60
Binary	p(1-p)	None beyond p
Poisson	\mu	None beyond \mu
Gamma	\mu^2/k	Shape k=25

For the last three models, changing the mean changes the variance. Increasing sample size and changing the response variation are different design questions.

Pause: assumption audit

Which assumption in your model most needs external evidence?

Explain how you could change it while preserving valid response values. State what else must change in the code to keep the intended mean.

Use the sensitivity task in the after-class worksheet.

Students examining the assumptions behind a simulated animal experiment.

Pause the recording before giving these suggested answers. There is no single correct choice of the assumption most in need of evidence. A defensible answer identifies the assumption, a relevant source of evidence, and a coherent change to the generator. Pilot data or previous studies should match the response definition, experimental unit and follow-up period.

Group 1: Gaussian pen mean daily gain. A reasonable concern is whether the conditional SD of 60 g/day is plausible for pen means after accounting for diet and room. Seek comparable pen-level data, rather than substituting the SD among individual pigs. For the worksheet sensitivity check, change the supplied SD from 60 to 90 in rnorm(), leaving eta, mu, the diet coefficient and room effects unchanged. Each pen’s generating mean stays fixed, while conditional variance increases from 3600 to 8100. Gaussian draws remain real-valued; negative daily gain is possible under this model. Wider tails may nevertheless challenge biological plausibility and should be inspected. No compensating change to mu is required.

Group 2: binary treatment indicator. A reasonable concern is the assumed baseline probability that at least one pig in a pen needs treatment. Evidence must concern this pen-level event over a comparable follow-up period, not simply the individual-pig treatment probability. The worksheet changes the Room 1 control probability from 0.30 to 0.15. Set beta_0 <- qlogis(0.15), retain beta_diet <- log(0.60) and the room contrasts, and recompute eta and mu = plogis(eta). Continue using rbinom(n(), size = 1, prob = mu), so responses remain 0 or 1. This task deliberately changes the means; the high-fibre probability in Room 1 becomes about 0.0957. It preserves the odds ratio, not the absolute probability difference. Bernoulli variance is mu * (1 - mu), so there is no independent dispersion parameter to change while holding mu fixed.

Group 3: Poisson treatment-event count. A reasonable concern is the expected event count over the specified follow-up period, or the assumption that conditional variance equals the mean. Seek comparable pen-level event counts, allowing repeat events, with exposure recorded. The worksheet changes the Room 1 control mean from 3 to 1.5. Set beta_0 <- log(1.5), retain beta_diet <- log(0.65) and the room contrasts, and recompute eta and mu = exp(eta). Continue using rpois(n(), lambda = mu), so observations remain non-negative integers. This deliberately changes the generating means and hence Poisson variances. The Room 1 high-fibre mean becomes 1.5 * 0.65 = 0.975. There is no separate Poisson variance setting that preserves mu. If the concern is extra variation at a fixed mean, a different generator, such as a negative-binomial model, would be an extension beyond this worksheet and would require an explicit dispersion assumption and a review of the analysis.

Group 4: Gamma feed conversion ratio. A reasonable concern is the assumed coefficient of variation of 0.20, corresponding to shape 25. Seek pen-level feed conversion ratios with positive intake and gain under comparable conditions. For the worksheet sensitivity check, keep eta and mu = exp(eta) fixed, change shape from 25 to 9, and also change scale from mu / 25 to mu / 9 in rgamma(). Since mean = shape * scale, each pen’s generating mean stays at mu. Variance increases from mu^2 / 25 to mu^2 / 9, and CV increases from 0.20 to one third. Responses remain positive. Changing shape alone while leaving scale at mu / 25 would change the mean to 9 * mu / 25, so it would confound a change in spread with a change in mean.

Resolve the wording about keeping the mean. For Groups 1 and 4, the worksheet asks for a spread change at a fixed conditional mean. For Groups 2 and 3, it explicitly asks for a new baseline mean. Here “the intended mean” means the new target, not preservation of the original mean. In all groups, changing the generator does not force a finite simulated dataset’s sample mean to equal its generating mean.

Common checks. Keep the randomised design and treatment effect fixed for these sensitivity comparisons. Use the same response seed as instructed, verify the support and the recalculated mu values, and record what changed. One pair of datasets cannot establish a power difference. Independence between pens conditional on the model is another defensible assumption to question; if it is doubtful, both the dependence structure in the generator and the planned analysis need attention.

Fitting and repeated experiments

Chapter 3: about 10 minutes

Generating and estimating

During simulation, we choose the parameters and draw possible responses.

During fitting, we give R the responses and model structure, and it estimates the unknown coefficients. It does not receive our generating values.

Likelihood compares how well different parameter values account for the observed data. glm() carries out the numerical fitting.

The model call

fit <- glm(response ~ diet + room,
           family = poisson(link = "log"),
           data = first)
coef(fit)["dietHigh fibre"]

dietHigh fibre 
    -0.4700036

The formula allows a diet contrast and room differences. The family call states both the data generating process and link.

The fitted diet coefficient estimates a log expected-count ratio. Its exponential estimates that ratio under this model.

Predictions on the response scale

new_pens <- data.frame(
  diet = factor(c("Control", "High fibre"),
                levels = levels(first$diet)),
  room = factor(rep("Room 1", 2),
                levels = levels(first$room))
)
predict(fit, newdata = new_pens, type = "response")

       1        2 
2.769231 1.730769

These are estimated mean event counts for Control and High fibre pens in Room 1, respectively. They do not include the random variation of a new observation.

Coefficients and uncertainty

A coefficient estimate summarises one dataset. Its standard error describes estimated sampling uncertainty under the model.ß

A confidence interval adds that uncertainty to the report. For a log-link coefficient, exponentiating interval endpoints gives a ratio interval.

With small or sparse datasets, familiar large-sample intervals can behave poorly.

Sparse binary data

Some simulated rooms may contain only zeros. Diet and room patterns can create separation or near-separation, leaving ordinary logistic estimates poorly determined.

Inspect warnings, fitted probabilities and coefficient standard errors. A convergence flag alone does not establish reliable estimation.

Retain these outcomes. Week 6 needs a predeclared strategy for unreliable fits and must report how often they occur.

What power asks

For a stated nonzero effect, model, design and analysis, power is the probability that the chosen test rejects its null hypothesis.

Simulation estimates this by drawing many possible experiments and applying the same analysis and rejection rule each time.

A single estimate close to the generating effect, or a single small p-value, cannot establish that probability.

A Week 6 simulation run

At each candidate sample size:

Construct the proposed design and specify the generating parameters.
Draw a response dataset.
Fit the planned model and apply a predeclared rejection rule.
Repeat and summarise rejection and analysis-failure rates.

We will first condition on a balanced allocation. Re-randomising within each repetition answers a related question averaged over treatment assignments.

Monte Carlo uncertainty

Suppose 800 of 1,000 simulated experiments reject the null under an assumed effect. The estimated power is 0.80.

Its approximate Monte Carlo standard error is:

\sqrt{\frac{0.80(1-0.80)}{1000}}\approx0.013

More repetitions reduce simulation noise. They do not make unrealistic biological assumptions more credible.

Sample size refers to independent pens

Pens per room	Rooms	Total pens	Pens per diet
4	4	16	8
8	4	32	16
12	4	48	24
16	4	64	32

This candidate grid increases independent pens while retaining four fixed room effects and balanced diets. Feasibility requires separate consideration.

Sensitivity and false positives

A useful power curve should state its effect, baseline, variation, allocation and analysis assumptions.

Compare it with at least one plausible alternative scenario. Also set the treatment effect to its null value and check the rejection rate.

An inflated false-positive rate can make an apparent gain in power misleading.

Your Week 6 preparation

Complete the after-class worksheet and bring:

a corrected and reproducible group script
one explained effect on the response scale
one assumption you would vary
a draft detection rule and a question about it

You will learn to repeat the full simulation and analysis in Week 6. You do not need to write the power loop beforehand.

Sources and further reference

McCullagh and Nelder (1989). Generalized Linear Models, 2nd ed.
R Core Team. Families, GLMs, predictions, and Gamma distributions.
Colegrave and Ruxton (2020). Power Analysis: An Introduction for the Life Sciences. Oxford University Press.