06-01 · Simulation-based power analysis

Build power curves and assess the assumptions behind a design recommendation.

Start here

No new preparation is required. Work in your Week 5 group. Open group-N-student.R in the supplied Posit Cloud project and run sections with Ctrl/Cmd+Enter. The script is complete: your tasks are to predict, inspect and explain the baseline results, then predict one assumption change. Save your answers in a document alongside the script.

Driver: runs code. Model translator: explains the mean, response and effect. Design auditor: checks units, assumptions and the result. Rotate at 13:55 and 15:10.

Files for the activity

Open the instructor-supplied Posit Cloud project for 06-power-analysis. Its scripts, data and reference files are already provided; open your assigned script there. Workshop slides

The supplied Posit Cloud project contains the shared helper and reference results required by the scripts. Individual script copies for reference:

Shared simulation helper

Your model reminder

Four rooms, initially four pens per room, six pigs per pen, and equal follow-up. Each room has half its pens on each diet. The pen receives diet independently; one row is one pen. Room 1 has zero room effect in the data generating process and is also the fitted reference room. Control is the reference diet in every room.

The score is beta_0 + beta_diet * (diet == "High fibre") + room_effect. Convert it to a mean, then draw an observation:

Group Room 1 control mean Mean from score Observation Diet effect
1: daily gain 695 g/day eta rnorm(..., sd = 60) +45 g/day
2: any pig treated 0.30 plogis(eta) rbinom(..., size = 1) Odds ratio 0.60
3: treatment events 3 exp(eta) rpois(...) Count ratio 0.65
4: FCR 2.80 kg/kg exp(eta) rgamma(..., shape = 25, scale = mu/25) Expected FCR ratio 0.92

The room effects are contrasts with Room 1, on the link scale, matching Week 5:

Group Room 1 Room 2 Room 3 Room 4
Gaussian 0 25 45 35
Binary 0 0.35 0.65 0.50
Poisson 0 0.20 0.40 0.30
Gamma 0 −0.06 −0.10 −0.04

For Groups 2–4 the diet coefficient is the logarithm of the displayed ratio. These are illustrative assumptions, not observed diet effects. The family’s support constrains observations; the inverse link constrains the mean.

1. Guided reconstruction — 12:30–13:00

Watch the shared Gaussian example in simulation-recap.R; no new coding is required. Work through these questions before using the power helpers:

  1. What does one row represent? Which unit receives a diet?
  2. We assume Room 1 control mean 695 g/day, a +45 diet effect and SD 60. What is the high-fibre mean? Are these assumed values or estimates?
  3. Why can two pens with the same expected response have different observations?
  4. Fit the model to the generated observations. Why need the estimated diet effect not equal +45?
  5. Generate and analyse a second dataset. Which of design, expected response, observation and fitted effect stay fixed?

Connect the example to your group’s response family. The stages are design, model score (\(\eta\)), expected response (\(\mu\)), observation, fitted model and decision. For Gaussian identity-link responses, \(\eta = \mu\). Explain why 200 simulated experiments of 16 pens do not create a real experiment with 3,200 independent pens.

2. Declare the question — 13:00

Write your response, effect scale and relevant effect. We use the supplied full and reduced models:

full:    response ~ diet + room
reduced: response ~ room

The null hypothesis is no diet effect. Detection means a comparison with p < 0.05 that passes the declared numerical screens; it does not additionally require the estimated sign to be correct. Record direction separately. Gaussian and Gamma comparisons use an F test; binomial and Poisson use a likelihood-ratio chi-squared approximation. The details are supplied in analyse_one(). These screens concern errors, warnings and numerical fit problems; they do not require QQ plots or other model-adequacy diagnostics to pass.

Before running: does a statistically interesting result establish biological importance? Does a result with p ≥ 0.05 establish no effect?

3. One repetition — 13:15

Run sections 1–2 of your script. Inspect the design counts and the columns eta, mu, and response. Inspect the returned estimate, p-value, convergence, boundary flag, warning text and unresolved.

Record the generating effect, fitted effect, test result and decision. If unresolved, explain why the declared procedure cannot count it as a detection. Do not change seeds until you obtain a preferred result.

4. Repetition — 13:55–14:15

Rotate roles. Run section 3. Count TRUE in detected, divide by all attempts, and check the summary. Record unresolved attempts separately. With a nonzero generating effect the question concerns power; with a zero effect it concerns false positives. Return to the supplied null reference during interpretation at 14:40; no new null run is required.

Prediction: what changes if we increase repetitions while keeping the experimental design and parameters fixed?

5. Baseline curve — 14:15–14:40

By 14:15, run section 4 with 200 repetitions at 16 pens. Then run section 5 across the four sizes. Reuse the 16-pen results in the curve.

Total pens Pens per room Pens per diet Pigs per pen
16 4 8 6
32 8 16 6
48 12 24 6
64 16 32 6

Only pens per room increase. Adding rooms is a different planning question. The design is held fixed within each sample-size batch, while responses are drawn afresh.

Stop computing at 14:40 and interpret until 14:55. If calculation is slow, use your group’s Baseline rows in the supplied reference and say that these are supplied results. Record total pens, attempts, detections, unresolved proportion and the 95% Monte Carlo interval. Compare your curve with reference/group-N-curve.png and reference/power-summary.csv (2,000 repetitions).

  1. What does each point estimate?
  2. Why can the exploratory curve have small irregularities?
  3. Do its intervals cover uncertainty about the assumed effect or only simulation noise?
  4. Does any tested size support an 80% target? If none does, say so.
  5. Is the proposed increase in room capacity feasible?

The vertical axis is operational detection probability: unresolved fits count as no detection and remain in the denominator. Always show their proportion beside the curve. Do not reinterpret the result as power conditional on a successful fit.

6. Reports — 15:10–15:25

Rotate roles. Give a two-minute report: response and assumed effect, one baseline curve finding and one limitation. Different response types answer different scientific questions; do not choose an outcome solely because its curve is higher.

7. Recommend and predict — 15:25–15:40

Write:

“Under [data generating process and parameter assumptions], with [allocation and model], the estimated probability of [testing event] at [total pens; pens per diet] is [estimate and interval], from [B] simulations with [unresolved proportion]. We recommend [design or further investigation], conditional on [feasibility and uncertain assumptions].”

Add one feasible design or measurement improvement. Predict the direction of your assigned change and explain what stays fixed. Do not run the sensitivity curve in class.

Group Sensitivity change
Gaussian SD 60 to 90; mean unchanged
Binary Room 1 control probability 0.30 to 0.15; odds ratio unchanged
Poisson Room 1 control mean 3 to 1.5; count ratio unchanged
Gamma Shape 25 to 9; mean unchanged, scale = mean/shape

Save baseline results using section 6 of the script. The later sensitivity section is optional code; supplied results support the after-class comparison without running it.

Project launch — 15:40–15:55

Discuss the group project and autumn-break scoping task. No project model fitting is required today.

8. Individual exit — 15:55

“A group reports that 32 pens give 80% power and proposes increasing the simulation count to improve the experiment. What information is missing, and what would increasing the simulation count actually improve?”

Keep your work

Save your annotated script, baseline curve, summary, recommendation and sensitivity prediction. Complete the week 42 revision and project-scoping task, capped at 150 minutes per student including revision. Compare the supplied sensitivity results within its 30-minute revision block. The consolidation worksheet provides optional depth, not an extra required hour. Bring your project brief, model map and one question on 21 October.