05-01 · Simulate before analysing
Your task
Your group will build a possible version of an experiment before anyone collects data. You will generate two datasets under the same assumptions and see what a GLM estimates from each.
Work in one of four groups of three. Rotate the driver, model translator and design auditor roles at 14:25 and 15:10. Open your assigned group-N-student.R file in Posit Cloud and run it section by section. No before-class video is required.
Files for the activity
Log in to Posit.cloud and open the 05-data-simulation project.
The project contains the worksheet, all four scripts and after-class materials. The instructor has already populated it. Open your assigned script and run it section by section. The links below are reference copies.
The study
A researcher proposes to compare control and high-fibre diets. There are four rooms, with four pens in each room and six pigs per pen. A shared feeder delivers diet to each pen. All pens have the same observation period.
Randomise two pens to each diet within each room. The pen is the experimental unit. Generate one pen-level response. For feed conversion ratio, assume every pen has positive gain and measured feed intake.
All numerical values below are invented teaching assumptions. The task is to understand and audit them.
Computing choices
We use tibble to construct the design table, dplyr for the familiar data manipulation verbs, and ggplot2 to compare the datasets visually. The stats package supplied with R contains the generators and glm(). The course Posit Cloud project should already have the software ready.
The starter contains four TODO expressions. Replace their NA placeholders and run the sections in order. An incomplete script will stop deliberately. Save your script regularly.
Task 1. The randomised design, 13:20–13:30
Complete TODO 1 with the blocking variable and TODO 2 with the operation that shuffles the four treatment labels without replacement.
Explain these operations to one another: rep(), set.seed(), group_by(), sample(), mutate() and ungroup().
Show the instructor:
- 16 rows and 16 distinct pen IDs
- two pens assigned to each diet in each room
- Control as the first level of the diet factor
What would go wrong if you sampled the labels with replacement? Why does a balanced count table alone not rule out duplicated pen IDs?
Your explanation:
____________________________________________________________________
Task 2. The assumption record, 13:30–13:40
Complete this for your assigned group before generating responses.
| Item | Your entry |
|---|---|
| Response definition and unit | |
| Possible response values | |
| Experimental unit and observation period | |
| Control expectation in Room 1 | |
| Diet effect and its scale | |
| Room effects and their scale | |
| Remaining variation | |
| One assumption needing evidence |
The room effects below are ordered Room 1, Room 2, Room 3, Room 4. Room 1 has zero effect, so the intercept describes the control response in Room 1 on the link scale. This matches the first factor level used by glm(). We treat these four room effects as fixed numbers.
Group 1 recipe: mean daily gain
One observation is the pen’s mean daily gain in grams/day. Use a Gaussian family with identity link.
| Quantity | Teaching assumption |
|---|---|
| Room 1 control expectation | 695 g/day |
| High-fibre diet effect | +45 g/day |
| Room effects | 0, +25, +45, +35 g/day |
| Conditional SD | 60 g/day |
The score and mean are the same: mu = eta. Draw observations with rnorm(n(), mean = mu, sd = 60).
The Gaussian model allows every real value. Daily gain can be negative, although extreme values would need biological scrutiny. The model describes variation among pen means.
Group 2 recipe: any pig treated
One observation is 1 if any pig in the pen requires veterinary treatment during the study, otherwise 0. It is not the number of pigs treated.
| Quantity | Teaching assumption |
|---|---|
| Room 1 control probability | 0.30 |
| High-fibre odds ratio | 0.60 |
| Room effects on log-odds scale | 0, +0.35, +0.65, +0.50 |
| Conditional variance | mu * (1 - mu) |
The script supplies beta_0 = qlogis(0.30) and beta_diet = log(0.60). Convert the score using mu = plogis(eta) and draw with rbinom(n(), size = 1, prob = mu).
The mean is a probability. The observation is 0 or 1. A binomial distribution with size 1 is Bernoulli. The odds ratio does not directly give the proportional change in probability.
Group 3 recipe: treatment events
One observation counts veterinary treatment events in the pen over the fixed follow-up period. Repeat events are possible, so the count is not bounded by six pigs.
| Quantity | Teaching assumption |
|---|---|
| Room 1 control expected count | 3 events |
| High-fibre expected-count ratio | 0.65 |
| Room effects on log scale | 0, +0.20, +0.40, +0.30 |
| Conditional variance | Equal to the expected count |
The script supplies beta_0 = log(3) and beta_diet = log(0.65). Convert the score using mu = exp(eta) and draw with rpois(n(), lambda = mu).
Observed counts are non-negative integers. The mean may be fractional. Equal follow-up and group sizes allow the count ratio to be read as a rate ratio here.
Group 4 recipe: feed conversion ratio
One observation is feed intake divided by positive mass gain for a pen, in kg feed/kg gain. Use a Gamma family with log link.
| Quantity | Teaching assumption |
|---|---|
| Room 1 control expectation | 2.80 kg/kg |
| High-fibre expected FCR ratio | 0.92 |
| Room effects on log scale | 0, −0.06, −0.10, −0.04 |
| Gamma shape | 25, giving a coefficient of variation of 0.20 |
The script supplies beta_0 = log(2.80) and beta_diet = log(0.92). Convert the score using mu = exp(eta) and draw with rgamma(n(), shape = 25, scale = mu / 25).
The mean and observations are positive. The scale conversion is supplied for today and explained in the after-class video.
Task 3. First dataset, 13:55–14:25
Complete TODO 3 with:
beta_0 + beta_diet * (diet == "High fibre") + room_effectExplain how the diet comparison acts as 0 or 1. Complete TODO 4 using your group’s inverse link.
The function shell is already supplied: generate_response(seed) runs the recipe again with the seed you give it. You do not need to write a loop.
Run the first draw with the supplied seed, 1501–1504. Inspect every calculated column. Show the instructor 16 responses, the response support check, and the range of the expected values.
Choose one pen and explain its score, mean and observed response in words.
Your pen and explanation:
____________________________________________________________________
Task 4. Another possible dataset, 14:25–14:55
Rotate roles. Before running the next draw, predict what will happen to diet, eta, mu and response.
Use the supplied second seed, 2501–2504. Keep the design and all model parameters fixed. Display the comparison table and two plots.
| Question | Your answer |
|---|---|
| Which columns stay identical, and why? | |
| Must every pen’s response change? | |
| Do both plots show the same apparent diet effect? | |
| What does this comparison tell us about observing one experiment? |
Retain both datasets, including unexpected patterns. Do not change seeds to obtain a preferred result.
Task 5. The planned GLM, 15:10–15:30
Rotate roles again. Run the supplied model for each dataset.
glm(response ~ diet + room, family = model_family, data = first)response is the observation, diet + room specifies the mean structure, and model_family supplies the family and link. Room appears because the generator allows systematic room differences.
The script extracts the high-fibre coefficient. For Group 1, report it in g/day. For Groups 2–4, also report its exponential as an odds ratio, expected-count ratio or expected FCR ratio.
| Quantity | First dataset | Second dataset |
|---|---|---|
| Generating treatment coefficient | ||
| Fitted treatment coefficient | ||
| Interpretation on the requested scale | ||
| Warnings or reasons for caution |
Small binary datasets can give unstable estimates, even with no warning. Ask the instructor about very large coefficients, standard errors, or fitted probabilities close to 0 or 1. A documented unstable fit is a valid result. Retain the data and explain the problem.
Task 6. Peer audit, 15:30–15:45
Pass the assumption record and results clockwise. Review the unit, support, mean, treatment interpretation and remaining variation.
Report one choice you can defend and one assumption worth changing in a sensitivity analysis. State why changing it might change the experiment’s ability to detect an effect.
Submission
Save your completed R script. Its final lines save both datasets into your group’s results folder.
Individual exit, 15:52–16:00
A single simulated dataset shows a clear diet difference. Why does that fail to establish adequate power? What would you keep fixed and what would you repeat?
Name one possible and one impossible value for your response. If your model excludes no finite real value, explain that.
After class
Follow 05-02 · GLM consolidation using the longer GLM deck and the recording when available. Bring your corrected simulation script and one sensitivity question to Week 6.