Your study, the full rehearsal

Self-paced bonus · Module 7 · the county’s six steps, on your own numbers

What this page is

The M07 lecture ran a Colorado county’s survey through six steps before a single real interview: a sentence to report, a simulated county, one sample, the analysis, a thousand rehearsals, and the sample size the board’s target needs. In class you ran two of those steps on your own numbers. This page runs all six.

It is optional and ungraded. Budget about an hour to work through it with the county’s numbers, and more to adapt it to your own study: finding good planning numbers and calibrating a new outcome take time. If you are planning a thesis or a grant, the last section, on when a formula and a simulation disagree about sample size, is the part to read slowly.

You enter your numbers once, in the first chunk; every later chunk reads them by name. The county’s numbers are filled in to start, so every chunk runs before you change anything. Do Step 2’s calibration for the county first (two edits), and the rest of the page reproduces the lecture’s numbers, so you can check what you see against it before switching to your own study.

Your numbers

Two kinds of numbers go in here. The first four are your study plan’s numbers, the same ones you typed in class. The next two are your scale’s limits: the lowest and highest score a person can get. Benevolence runs from 1 to 5; a PHQ-9 score runs from 0 to 27; birth weight has no ceiling worth modelling, so use -Inf and Inf to switch the limits off.

The last two are dials. They are what you will feed to rnorm(), and they start out equal to your published mean and SD. Step 2 explains when to turn them. When you switch from the county to your own study, set the dials back to my_mean and my_sigma first, then calibrate again if your scale has limits.

The six steps, on your numbers

Step 1 — Write the sentence you’ll report

Before any data, write the finished sentence with blanks in it, in your outcome’s units. It gives every later step a target:

In our study of ___ (your population), the mean ___ (your outcome) was ___ (95% CI: ___, ___).

The requirement: the interval reaches no more than ±\(E\) either side of the estimate, where \(E\) is the my_E you entered above.

Step 2 — Build a simulated population, and check it against your source

The county’s scores could not follow a smooth, unlimited bell curve, because the scale stops at 1 and at 5. Yours may have the same feature. Build it in, then check that the simulated population reproduces your published mean and SD.

If your outcome has no limits, the mean and SD already match, and the two shares are 0. Move on.

If it has limits, the walls catch the people in the tails and pull the spread in, so the SD comes out smaller than your published one. That is what happened to the county: the published 3.14 and 1.04 produced an SD of about 0.99. The fix is to turn the dials a little past the target so that what comes out approximately matches. Raise dial_sd in the first chunk (the county went from 1.04 to 1.10), nudge dial_mean if the mean drifted (the county went from 3.14 to 3.15), rerun the first chunk and then this one, and repeat until the mean and SD come out close to your source. Close is enough: the county’s settle at 3.14 and about 1.03, against the trial’s 1.04. From here on, every chunk uses the dials, so you only turn them once.

Do it now for the county, before you enter your own numbers. In the first chunk, set dial_mean <- 3.15 and dial_sd <- 1.10, rerun it, then rerun this one. The mean and SD land on 3.14 and about 1.03, and every later chunk now reproduces the lecture: a coin flip at 67 interviews, about 8 in 10 at 75.

One thing to keep straight, because it matters in Step 6: the dials describe the latent column, before the limits apply. Whenever a calculation needs the spread of actual scores, as the sample-size formula does, use your published SD, my_sigma.

The simulation lets scores land anywhere between your limits. Many outcomes can’t: benevolence is the average of six items rated 1 to 5, so real scores move in steps of 1/6, and a sum score like the PHQ-9 moves in whole numbers. To build that in, change the walls line in the first chunk so it rounds after the limits:

walls <- function(x) round(scales::squish(x, range = c(my_min, my_max)) * 6) / 6   # an average of 6 items
walls <- function(x) round(scales::squish(x, range = c(my_min, my_max)))           # a sum of whole numbers

Use one of the two lines, not both. Rounding changes the model that generates the data, so after adding it, recalibrate the dials in this step and rerun the later steps. For the county it barely matters: the SD moves by about a thousandth, and the share of studies meeting ±0.25 at 75 interviews by less than a point.

Step 3 — One sample, and the checks the real data must pass

Simulate one sample of my_n people from the simulated population, and examine it the way you will examine the real data: a histogram first, then summary numbers.

Now write the checks. The simulated sample tells you what the real data must look like, so write that down as code. When the real data arrive, these lines run first. If a data-entry slip puts an impossible value in the file, the analysis stops instead of quietly averaging it in.

Starting point

You know what the real data must look like: as many people as the plan calls for, no missing scores, and every score inside the scale’s limits. You want R to check that automatically, every time, so a bad file stops the analysis instead of slipping through.

What the code does

n_planned <- my_n stores the planned sample size in one named place. pull() takes the y column out of my_sample as a plain list of numbers, saved as y. stopifnot() then runs each check in turn. Each check is a statement that should be TRUE, with a plain-language label in quotes before the =:

  • length(y) == n_planned counts the scores and asks whether there are exactly as many as the plan calls for. Because the check uses the name, not a number, a change of plan means changing one line. The double == asks whether two things are equal; a single = gives something a name, like the label in front of each check.
  • !anyNA(y): anyNA() asks “is any score missing?” and ! flips the answer, so the whole check reads “no score is missing.”
  • all(y >= my_min & y <= my_max): for each score, y >= my_min & y <= my_max asks “is it at least the minimum and at most the maximum?”, giving one TRUE or FALSE per score. all() asks whether every one of them is TRUE.

If every check is TRUE, stopifnot() stays silent and R moves on to the last line, which prints the success message. If a check is FALSE, R stops right there and shows that check’s label as an error.

Key functions

Function What it does
stopifnot(“label” = check, …) Stops with an error showing the label of the first check that is FALSE. Silent when every check is TRUE.
pull() Takes one column out of a data frame as a plain list of values.
anyNA(x) TRUE if any value of x is missing.
all(x) TRUE only if every value of x is TRUE.
==, !, & “Is equal to,” “not,” and “and.”

How to read output

If the data pass, you see the quoted success sentence. If a check fails, you see a red error such as Error: scores must lie within the scale's limits: the label of the check that failed. The label is only text for humans; what R actually tests is the code after the =.

Break one on purpose: in the third check, change y <= my_max to y <= my_mean and run it. Almost every sample has at least one score above the mean, so the check fails, and the red message names it. That red message is exactly what you want to see when data are wrong. Then change it back.

Step 4 — The interval, two ways

Write the analysis now, on simulated data, so it is ready and debugged when the real data arrive. First the parametric route, one step per line, then t_test() to confirm it:

Then the Module’s second route, the bootstrap: resample your my_n simulated people with replacement 1,000 times, and keep the middle 95% of the means.

The two routes give similar intervals in most samples. They need not match exactly, because the parametric and percentile-bootstrap methods make different approximations, and they can drift apart when the sample is small or the scores pile up against a limit. Fill in Step 1’s sentence with either interval, and label it simulated: these numbers describe your simulated population, not your real one. What you keep is the code. When the real data arrive, replace the chunk that creates my_sample with a line that reads them in, rerun the checks and this step, and fill in the same sentence with real numbers.

Step 5 — Rehearse your study 1,000 times

One simulated sample cannot tell you how precise your plan is. A thousand can. This is the lecture’s Step 5 on your numbers, and unlike the quick check you ran in class, it keeps your scale’s limits.

Read the sampling distribution first: its SD is the standard error at your my_n, and the dashed lines show how far a single study could land from the truth by chance alone. Then the widths: if most bars sit to the right of the gold line, your feasible sample is not enough for your target, and Step 6 is where you find out how many people are.

Step 6 — How many people: the formula, then the simulation

The formula first. Solve the half-width formula for \(n\) and you get an approximate starting sample size: roughly where a typical interval reaches your target:

\[n \approx \left( \frac{1.96 \cdot \sigma}{E} \right)^2\]

Then check it with a rehearsal. The formula aims at the typical sample, so at n_formula roughly half of the rehearsed studies meet the target and half miss (exactly how many depends on your outcome). This chunk rehearses your study at any size you name. Start with the formula’s answer, then raise n_try until about 8 in 10 studies meet your target (9 in 10 if your study is high-stakes).

Then the whole curve. Rather than trying sizes one at a time, rehearse every size from 10 up to twice the formula’s answer, and draw the result. This is the lecture’s precision curve, built by simulation on your numbers.

Starting point

You have code that rehearses your study at one sample size, and you want to run it at many sizes and collect the results in one table.

What the code does

rehearse_at() is the one-size rehearsal wrapped in a function: give it n, and it hands back a one-row table with that size’s typical margin, the 80th percentile of the margins (the margin 8 in 10 studies come in under), and the share of studies that met your target. sizes is the list of sample sizes to try. map() runs rehearse_at() once for each size and collects the one-row tables in a list; list_rbind() stacks them into one table, one row per size. The rest is the graph: a solid line for the typical study, a dashed line for the 80th percentile, the gold target, and two dotted lines at the formula’s \(n\) and at the first size where 8 in 10 studies met the target.

Key functions

Function What it does
function(n) { ... } Wraps a block of code so it can be run again with a different n. The last line of the block is what it hands back.
map(x, f) Runs the function f on each element of x and collects the results in a list.
list_rbind() Stacks a list of data frames into one data frame.
first(default = NA) The first value, or NA if there is none. Here: the first size that reached 8 in 10.
quantile(x, 0.80, names = FALSE) The value that 80% of the values fall below; names = FALSE drops the label so the column holds plain numbers.
linetype = "dashed" Draws a line dashed instead of solid, so the two lines are easy to tell apart.

How to read output

The line falls fast at first and then slowly, the \(\sqrt{n}\) pattern. Where the solid line crosses the gold target is the formula’s \(n\); where the dashed line crosses it is the size at which 8 in 10 studies meet the target. The last table gives both numbers.

Read the 8-in-10 \(n\) as approximate. It is the first size tried, and the 25 sizes are spaced several people apart, so the true minimum can sit a little lower. Each point also rests on only 500 rehearsals, so shares near 80% wobble by a few percentage points from run to run. The curve shows you the range to look in. Before a real decision, rerun Step 6’s one-size chunk at consecutive sizes around that \(n\), with more rehearsals (change replicate(1000, ...) to 5,000).

Formula or simulation? A closer look

The two dotted lines on your curve are the two answers to “how many people?” The formula’s \(n\) is where the typical study just meets your target. The rehearsal’s \(n\) is where you can be fairly sure that your one real study will. For the county those were 67 and about 75 interviews, and the gap between them is not a rounding error. Here is what each method assumes, and when the difference matters.

What the formula assumes.

  • The SD is known and fixed. It is neither. Your sample will have its own SD, a little above or below the published one, and its margin of error moves with it. The formula ignores that, so its \(n\) delivers the target in roughly the typical sample; in the county’s model, about half of samples miss it.
  • Nothing about the shape of the scores, or about the design. The formula needs the sampling distribution of the mean to be roughly Normal, not the scores themselves. But it can’t see a shape that makes each sample’s SD wobble more (skewed scores, scores piled against a limit or clumped at a few values), and it can’t see people clustered in clinics or dropping out.
  • \(z^{\star} = 1.96\) instead of \(t^{\star}\). For 95% confidence and a modest \(n\), \(t^{\star}\) is a little larger (2.02 at \(n = 40\)), so the formula runs slightly low. For planning, this one is minor.
  • One statistic. The formula is for a mean. A proportion, a difference between groups, a regression slope, or a change score each need their own formula, and some designs have none.

What the simulation adds.

  • Assurance. Because every rehearsed sample has its own SD, the rehearsal shows the whole spread of margins, not just the middle one, and lets you choose how sure you want to be: 8 in 10 (the county’s choice, borrowed from the usual target for power; there is no settled convention for an interval’s width), or 9 in 10 when falling short would be costly. For the county, that choice moves the answer from 67 interviews (a coin flip) to 75 (8 in 10) or 80 (9 in 10).
  • Your data’s real features. Limits, skew, clusters, dropout: if you can describe how the data will be generated, the rehearsal can include it. This page’s template builds in only the first, your scale’s limits (and, optionally, whole-number steps); skew, clustering, and dropout each need a few lines of your own code in the rehearsal. For the county, the 1-to-5 limits made each sample’s SD a little steadier, which is why the county’s 80% at 75 interviews is higher than the 72% a plain bell curve gives at the same size; a plain bell curve needs about 78 people to reach 8 in 10. That gap is small here. With a heavily skewed outcome, or people clustered in a handful of clinics, it can be large.
  • A test of the whole plan. The same simulated data already ran your checks and your analysis code before any real data existed.
  • It works when there is no formula. Any design you can simulate, you can size.
Question Formula Simulation
How many people for a typical interval of ±\(E\)? Yes, in one line Yes, and it agrees with the formula
How many for a dependable ±\(E\), 8 in 10 times? Not without an extension Yes: where the dashed line crosses your target
Does it respect my scale’s limits, skew, or clustering? No Yes, if you build them in
Does it work for a design with no textbook formula? No Yes
Will a reviewer recognise it? Yes, immediately Yes, if you show the code and the assumptions

The extension the table mentions. Statisticians have worked out formulas that add an assurance level to the precision target, under the heading accuracy in parameter estimation (Maxwell, Kelley, & Rausch, 2008, give an accessible review; the MBESS package in R implements the calculations for many designs). They are the right tool when a proposal wants a citable formula with assurance built in. Their assurance calculation assumes Normal scores, so for an outcome with limits, skew, or clusters, a rehearsal remains the check. G*Power, which many proposals cite, solves formulas for power (next week’s idea) rather than for an interval’s width, and it takes the same fixed-SD view.

What to do with this in a proposal. Report the formula’s \(n\) as the first, transparent number. Then report the rehearsal’s \(n\) as the one you are planning for, with two sentences on what the simulation assumed (the published mean and SD, the scale’s limits, the sample-to-sample variation in SD). Reviewers read that as care, not as hedging.

What this adds to your study plan

Three things to know about your study before next week:

  • The formula’s \(n\), for a typical interval of ±\(E\).
  • The rehearsal’s \(n\), for a dependable one, and the assurance you chose (8 in 10, or 9 in 10 for a high-stakes study).
  • What the rehearsal assumed that the real study might not honour: the borrowed mean and SD, the scale’s limits, and anything the simulation left out (nonresponse, clustering, dropout).

Next week the plan gains an intervention and a decision rule, and the same rnorm() line you rehearsed here becomes the first line of a power simulation.

Back to the M07 lecture → Creating, interpreting and utilizing confidence intervals