library("mgcViz")
library("DHARMa")
library("see")Case 1 · Pen event counts
Read the output from two count models and decide what the evidence supports.
All four cases · Diagnostic report
Use this page for the shared count example or revisit it on your own. Complete each formula in your R console or write it down before revealing the solution. Point to a number, label or pattern in the displayed output when you answer an interpretation question. To run the analysis, start from the supplied Posit Cloud project root and work through the R blocks in order. If you have yet to fit a model yourself, run the fit in its solution callout before using it in the next block.
Load the diagnostic packages
Run this block in Posit Cloud before starting the case. If RStudio prompts you to install a missing package, install it and run the block again. You can also install the three packages from the Console with install.packages(c("mgcViz", "DHARMa", "see")).
What did the researchers count?
Researchers assigned Control or High fibre diets to 64 pens in four rooms. Each room contained eight pens per diet. The researchers recorded one count of treatment events per pen over the same follow-up time. They want to know how expected event counts compare between the two diets after allowing rooms to differ.
The data are invented for teaching. Set Control and R1 as the reference categories, then check how many rows each room and diet combination contributes.
library("dplyr")
d <- readr::read_csv("data/case-1.csv")
d <- d |>
mutate(
diet = factor(diet, levels = c("Control", "High fibre")),
room = factor(room, levels = c("R1", "R2", "R3", "R4"))
)
head(d, 3)# A tibble: 3 × 4
room pen diet events
<fct> <chr> <fct> <dbl>
1 R1 R1-P01 Control 1
2 R1 R1-P02 Control 1
3 R1 R1-P03 Control 2d |> count(room, diet)# A tibble: 8 × 3
room diet n
<fct> <fct> <int>
1 R1 Control 8
2 R1 High fibre 8
3 R2 Control 8
4 R2 High fibre 8
5 R3 Control 8
6 R3 High fibre 8
7 R4 Control 8
8 R4 High fibre 8
NoteQuestion 1. Read the study design
What does one row measure? Which unit did researchers assign a diet to? Use the table to explain whether any room contributes more pens than the others.
TipReveal answer 1
A row gives the number of treatment events recorded in one pen. Researchers assigned a diet to each pen. Every room contributes eight Control pens and eight High fibre pens, so averaging over the recorded pens gives each room the same weight.
Fit a Poisson model
The Poisson model assumes that the variance equals the mean at each diet and room combination. Its default log link puts the coefficients on a log scale.
Read the formula below. Which column records the count? Which terms compare diets while allowing the expected counts to differ across rooms?
poisson_fit <- glm(
events ~ diet + room,
family = poisson(),
data = d
)
summary(poisson_fit)
Call:
glm(formula = events ~ diet + room, family = poisson(), data = d)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.9165 0.1821 5.032 4.86e-07 ***
dietHigh fibre -0.4313 0.1361 -3.168 0.00153 **
roomR2 0.9223 0.2058 4.482 7.40e-06 ***
roomR3 0.6306 0.2155 2.927 0.00343 **
roomR4 0.3747 0.2261 1.657 0.09752 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 250.77 on 63 degrees of freedom
Residual deviance: 216.28 on 59 degrees of freedom
AIC: 380.22
Number of Fisher Scoring iterations: 5
NoteQuestion 2. Read a model summary
Find the row for High fibre. What is the estimated coefficient and its standard error? Is the coefficient a difference in event counts, or is it on another scale? What does its sign tell you?
TipReveal answer 2
The estimate is about −0.431, with a standard error of 0.136. The coefficient is the log of the expected count ratio for High fibre versus Control. A negative log ratio corresponds to a count ratio below one, so the model estimates fewer events under High fibre. The standard error is on the same log scale as the coefficient.
Check variation and the residual pattern
check_overdispersion() checks whether counts vary more than the model allows. A dispersion ratio near one indicates that the amount of variation is close to what the check expects.
For each pen, simulate_residuals() uses DHARMa to simulate counts from the fitted model and locate the observed count among them. Those locations are called randomised quantile residuals. If the model fits, the residuals should be spread evenly from zero to one. The QQ plot compares their quantiles with uniform quantiles on that range, so points should lie roughly along the diagonal. The seed argument makes the simulation reproducible.
library("performance")
check_overdispersion(poisson_fit)# Overdispersion test
dispersion ratio = 3.975
Pearson's Chi-Squared = 234.526
p-value = < 0.001plot(simulate_residuals(poisson_fit, seed = 9401))
NoteQuestion 3. Diagnose the Poisson fit
Use the dispersion ratio and the QQ plot to name one concern about the fit. Suggest two possible explanations. What additional records would help distinguish them?
TipReveal answer 3
The dispersion ratio is about 4.0, and the QQ points flatten near the bottom and top of the residual range. Both point to a mismatch between observed and expected count variation. Pens could have different underlying event rates within the same diet and room. Shared management conditions within rooms could also matter. Records such as pen-level health, stocking or management information would help investigate these explanations. Equal follow-up makes unequal time at risk a less plausible explanation in this example.
Allow extra variation and check again
A negative binomial model allows extra variation among pens with the same expected count. Keep the outcome and the diet and room terms fixed so the two models answer the same question. method = "ML" fits the model by maximum likelihood. The dispersion check below uses 500 simulations.
NoteWrite the negative binomial formula
Complete the model formula. Which argument changes how much variation the model can allow? Predict how the dispersion ratio might change.
Formula template. Fill in the blanks before running it.
nb_fit <- mgcv::gam(___, family = mgcv::nb(), data = d, method = "ML")
TipCheck the negative binomial formula
family = mgcv::nb() allows the variance to exceed the mean. Both formulas use the same terms for the expected count.
nb_fit <- mgcv::gam(
events ~ diet + room,
family = mgcv::nb(),
data = d,
method = "ML"
)
summary(nb_fit)
Family: Negative Binomial(1.252)
Link function: log
Formula:
events ~ diet + room
Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.9007 0.3097 2.908 0.00364 **
dietHigh fibre -0.3967 0.2642 -1.501 0.13331
roomR2 0.9071 0.3790 2.394 0.01668 *
roomR3 0.6083 0.3847 1.581 0.11382
roomR4 0.4219 0.3891 1.084 0.27834
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.0796 Deviance explained = 10.7%
-ML = 149.04 Scale est. = 1 n = 64check_convergence(nb_fit)[1] TRUEset.seed(9401)
check_overdispersion(nb_fit, iterations = 500)# Overdispersion test (using simulated residuals)
dispersion ratio = 0.868
p-value = 0.78No overdispersion detected.plot(simulate_residuals(nb_fit, seed = 9401))
NoteQuestion 4. Compare the checks
What does the numerical fitting check report? How have the dispersion ratio and the QQ pattern changed? Name one question about the pens that these checks leave open.
TipReveal answer 4
The numerical fitting check returns TRUE. The dispersion ratio drops to about 0.87, closer to one, and the QQ points follow the diagonal more closely. The model better describes the amount of count variation in this example. We would still need information about the pens to know why their event rates differ.
Compare expected counts under the two diets
newdata = "balanced" supplies one prediction row for each diet and room combination. Each room contributes equally to the mean prediction for each diet. comparison = "lnratioavg" takes the log of the ratio of the two mean predictions. transform = "exp" returns the estimate and interval to the count-ratio scale. A ratio of one means equal expected counts.
library("marginaleffects")
avg_comparisons(
poisson_fit,
variables = "diet",
newdata = "balanced",
comparison = "lnratioavg",
transform = "exp"
)
Estimate Pr(>|z|) S 2.5 % 97.5 %
0.65 0.00153 9.3 0.497 0.848
Term: diet
Type: response
Comparison: ln(mean(High fibre) / mean(Control))avg_comparisons(
nb_fit,
variables = "diet",
newdata = "balanced",
comparison = "lnratioavg",
transform = "exp"
)
Estimate Pr(>|z|) S 2.5 % 97.5 %
0.673 0.133 2.9 0.401 1.13
Term: diet
Type: response
Comparison: ln(mean(High fibre) / mean(Control))
NoteQuestion 5. Interpret both estimates and intervals
For each model, read the estimate and the two 95% limits. Translate each estimate into a percentage change in expected events. Which interval includes one? Explain how the change in model affects the strength of your conclusion.
TipReveal answer 5
The Poisson ratio is about 0.65, or 35% fewer expected events, with an approximate 95% interval from 0.50 to 0.85. The negative binomial ratio is about 0.67, or 33% fewer expected events, with an interval from 0.40 to 1.13. The second interval includes one. The estimated reduction is similar across the two fits, but once the model allows extra variation, its interval also includes equal expected counts and some increases.
Assess proposed reductions
NoteDecide whether the evidence supports a target
Suppose the researchers would recommend High fibre only if the evidence supported at least a 20% reduction in expected events. A ratio of 0.80 marks that threshold. For each fit, find whether its entire 95% interval is below 0.80, above it, or crosses it. Then repeat for a 30% reduction, using 0.70 and the same fits and room weights. Does either target change your recommendation?
TipCheck the target decisions
Both intervals cross both thresholds. The two estimates suggest reductions greater than 30%, but the uncertainty leaves smaller reductions plausible under either model. Neither model’s interval supports recommending High fibre on the basis of either target.
NoteQuestion 6. Write a conclusion you could defend
In two sentences, say what you would report about expected event counts under High fibre. Name the comparison, its uncertainty and one limitation. Use the checks to explain which model you would report.
TipReveal answer 6
Across the four study rooms, I would report that the negative binomial model estimates about 33% fewer expected treatment events under High fibre than Control, with an approximate 95% ratio interval of 0.40 to 1.13. Its dispersion and QQ checks describe the count variation better than the Poisson model, but its interval still allows a reduction, equal expected counts or an increase.