Estimating R

How Epidemics Spread

Chapter 8 · Estimating R

Every model so far began with a number handed to us: an R of 2, a contact rate, a transmission chance. Nobody hands those to a health officer. What arrives each morning is a list of positive tests, a few hospital admissions, and some interviews with people who fell ill. In 1948 the Framingham Heart Study began following thousands of ordinary residents of one Massachusetts town for years, to see who developed heart disease and what they had in common (per search summaries, around 5,200 people were recruited). In 1950 Richard Doll and Austin Bradford Hill published a very different design in the British Medical Journal: they asked hospital patients with lung cancer, and matched patients with other illnesses, about their smoking, and compared the two groups. The first design starts with the exposure and waits for the illness; the second starts with the illness and looks back at the exposure. This chapter uses simulations to ask how we get from messy counts to a number like R, which studies answer which questions, and how the way we count can mislead us before any model is applied.

What was run, and what to trust

Every Python snippet in this chapter and in the three exercise solutions was executed, and the output in the comments is what Python printed (random seeds are fixed, so you should get the same numbers with the same Python random generator). The history and methods come from search-result summaries, which are often truncated or model-written, so check the originals before quoting: Doll and Hill, “Smoking and carcinoma of the lung: preliminary report”, British Medical Journal 1950 (the summaries I saw gave its design as hospital patients with lung cancer compared with matched patients with other diseases, but I did not see its tables); the Framingham Heart Study began in 1948 as a prospective study; the link between a growth rate and R through the shape of the generation interval is Wallinga and Lipsitch, Proceedings of the Royal Society B 274 (2007); and the correction of test results for sensitivity and specificity is Rogan and Gladen, American Journal of Epidemiology 107 (1978). The search summaries gave me the general form R = 1 / M(−r), where M is the moment generating function of the generation interval. The special case for a gamma-shaped interval, R = (1 + r × scale) ^ shape, is my own derivation from that form (one search summary said the same), and it is checked here only against my simulations: treat it as an inference, not a quotation. I did not check the methods named in passing in this chapter’s text (the Cori and Wallinga-Teunis ways of estimating R over time). Every population size, risk, sensitivity, specificity, test count, serial interval and true R is assumed or illustrative. No calculus is used; the only maths beyond arithmetic is logarithms and a least-squares line.

Estimating R from a Growth Curve

In the early days of an outbreak, with nobody immune and nothing changed, cases grow exponentially, as in How Epidemics Spread 2. If we plot the logarithm of daily cases against the day, the points fall near a straight line, and the slope of that line is the growth rate r per day. A growth rate is easy to measure. R is what we want, and it needs one more ingredient: how long it takes for one case to infect the next. The serial interval is the time between the start of symptoms in one case and the start of symptoms in the person they infected. It is measured from interviews. The generation interval is the time between the two infections, which usually cannot be seen. The two have the same mean in simple models but different spreads, and in practice the serial interval is used as a stand-in (a point I am relying on from general knowledge, not from a source I checked this session).

A shortcut is R = exp(r × mean interval). It is exact only if every infection takes exactly the mean time. When the interval varies, as it always does, the correct conversion is smaller. For an interval with a gamma shape (a flexible curve with a mean and a spread), my derivation gives R = (1 + r × scale) ^ shape. The code simulates outbreaks with a known R, fits the growth rate on days 10 to 20, and compares the two conversions. Each case has a random number of children (Poisson, mean R) and each child is infected a random whole number of days later. The serial interval is assumed gamma-shaped with a mean of 5 days and a standard deviation of 2.5 days. The four helper functions are reused in the exercises.

# Estimating R from early growth. Serial interval, window, sizes and true R are all assumed.
import math
import random
import statistics

# ---- helpers (the same four functions are reused in Exercises 1 and 2) ----
def poisson(lam):
    """Poisson random number (Knuth's method; fine for small lam)."""
    limit, k, p = math.exp(-lam), 0, 1.0
    while True:
        p *= random.random()
        if p <= limit:
            return k
        k += 1

def simulate(R, days, shape, scale, start=20):
    """New cases per day. Each case has Poisson(R) children; each child is infected a whole number
    of days later, drawn from a gamma 'serial interval' (shape, scale) and rounded, minimum 1 day."""
    daily = [0] * (days + 1)
    daily[0] = start
    for day in range(days + 1):
        for _ in range(daily[day]):
            for _ in range(poisson(R)):
                lag = max(1, round(random.gammavariate(shape, scale)))
                if day + lag <= days:
                    daily[day + lag] += 1
    return daily

def growth_rate(daily, first, last):
    """Slope of log(cases) against day, by least squares over days first..last (None if a day has no cases)."""
    if min(daily[first:last + 1]) == 0:
        return None
    xs = list(range(first, last + 1))
    ys = [math.log(daily[d]) for d in xs]
    mx, my = sum(xs) / len(xs), sum(ys) / len(ys)
    return sum((x - mx) * (y - my) for x, y in zip(xs, ys)) / sum((x - mx) ** 2 for x in xs)

def r_to_R(r, shape, scale):
    """Growth rate to R for a gamma serial interval: R = (1 + r * scale) ** shape."""
    return (1 + r * scale) ** shape

random.seed(3)
shape, scale = 4, 1.25            # assumed serial interval: gamma, mean 5 days, sd 2.5 days
print("true R   mean estimate   sd     naive exp(r * 5)   died out")
for true_R in [1.2, 1.5, 2.0, 3.0]:
    est, naive, died = [], [], 0
    for _ in range(100):
        r = growth_rate(simulate(true_R, 30, shape, scale), 10, 20)
        if r is None:                 # the outbreak died out before day 20: nothing to fit
            died += 1
            continue
        est.append(r_to_R(r, shape, scale))
        naive.append(math.exp(r * shape * scale))
    print(f"{true_R:6.1f} {statistics.mean(est):14.2f} {statistics.stdev(est):6.2f} {statistics.mean(naive):16.2f} {died:10d}")

print()
print("true R = 2.0, but we assume a different mean serial interval (same shape):")
daily_runs = [simulate(2.0, 30, shape, scale) for _ in range(100)]
rates = [r for r in (growth_rate(d, 10, 20) for d in daily_runs) if r is not None]
for assumed_mean in [3, 4, 5, 6, 7]:
    est = [r_to_R(r, shape, assumed_mean / shape) for r in rates]
    print(f"assumed mean {assumed_mean} days: mean R estimate {statistics.mean(est):5.2f}")
# true R   mean estimate   sd     naive exp(r * 5)   died out
#    1.2           1.24   0.22             1.25          2
#    1.5           1.56   0.19             1.61          0
#    2.0           2.04   0.12             2.18          0
#    3.0           3.01   0.07             3.56          0
#
# true R = 2.0, but we assume a different mean serial interval (same shape):
# assumed mean 3 days: mean R estimate  1.56
# assumed mean 4 days: mean R estimate  1.78
# assumed mean 5 days: mean R estimate  2.04
# assumed mean 6 days: mean R estimate  2.32
# assumed mean 7 days: mean R estimate  2.62
    

The gamma conversion recovers the true R well: 2.04 for a true 2.0 and 3.01 for a true 3.0. The shortcut overshoots, by a little at R = 1.2 (1.25) and by a lot at R = 3 (3.56), because it pretends that nobody infects anyone sooner than the average. The small upward drift at the low values (1.24 for 1.2 and 1.56 for 1.5) has two likely causes that I did not separate: rounding of lags to whole days, and, at R = 1.2, the two runs that died out and were left out, which removes the unluckiest outbreaks (both are inferences, untested). The second table shows the weak point. The simulation knows the serial interval exactly; a real analyst does not. If the true mean is 5 days but we assume 3, the estimate for a true R of 2.0 falls to 1.56, and an assumed 7 gives 2.62. An estimate of R is only as good as the interval it was converted with.

	Growth rate r	Serial interval	R
What it is	slope of log daily cases	onset of symptoms to onset in the next case	new cases per case
Observed directly?	yes, from the counts	from interviews (a sample)	no, always inferred
Main risk	noisy when counts are small	wrong mean gives wrong R	inherits both errors
Result in the simulation (true R = 2.0)	fitted on days 10 to 20	mean 5 days, sd 2.5 assumed	2.04 (gamma), 2.18 (shortcut)

Cohort and Case-Control Studies

Counts tell us how fast an epidemic grows, but not why. To find what puts people at risk (a food, a gathering, a type of work) we compare groups. A cohort study takes a group of people, records who was exposed, follows everyone, and compares the risk of illness in the exposed with the risk in the unexposed. Their ratio is the risk ratio (RR). It is the cleaner design but needs everyone followed, which is expensive when the illness is rare. A case-control study goes the other way: collect the ill people (cases), pick some well people (controls), and ask both groups about the exposure. It cannot give risks, because the numbers of cases and controls were chosen by the investigator, but it gives the odds ratio (OR): the odds of exposure among cases divided by the odds of exposure among controls. The code builds a made-up town of 20,000 people with 30 per cent exposed, then runs both designs on it. Two illnesses are tried with the same true risk ratio of 5: a rare one (5 per cent of the exposed ill against 1 per cent of the unexposed) and a common one (50 per cent against 10 per cent). The intervals use the standard log method and are approximate.

# Cohort versus case-control study in a made-up population. All risks and sizes are assumed.
import math
import random

def make_population(n, exposed_share, risk_exposed, risk_unexposed):
    """List of (exposed, ill) pairs."""
    people = []
    for _ in range(n):
        exposed = random.random() < exposed_share
        ill = random.random() < (risk_exposed if exposed else risk_unexposed)
        people.append((exposed, ill))
    return people

def cohort(people):
    """Follow everyone: risk ratio with an approximate 95% interval (log method)."""
    a = sum(1 for e, i in people if e and i)          # exposed, ill
    b = sum(1 for e, i in people if e and not i)      # exposed, well
    c = sum(1 for e, i in people if not e and i)      # unexposed, ill
    d = sum(1 for e, i in people if not e and not i)  # unexposed, well
    rr = (a / (a + b)) / (c / (c + d))
    se = math.sqrt(1 / a - 1 / (a + b) + 1 / c - 1 / (c + d))
    return rr, rr * math.exp(-1.96 * se), rr * math.exp(1.96 * se)

def case_control(people, n_cases, n_controls):
    """Sample some ill people and some well people: odds ratio with an approximate 95% interval."""
    cases = random.sample([p for p in people if p[1]], n_cases)
    controls = random.sample([p for p in people if not p[1]], n_controls)
    a = sum(1 for e, i in cases if e)
    b = n_cases - a
    c = sum(1 for e, i in controls if e)
    d = n_controls - c
    orr = (a * d) / (b * c)
    se = math.sqrt(1 / a + 1 / b + 1 / c + 1 / d)
    return orr, orr * math.exp(-1.96 * se), orr * math.exp(1.96 * se)

random.seed(11)
print("scenario                     study size  true RR  cohort RR (95% interval)  true OR  case-control OR (95% interval)")
scenarios = [("rare illness (5% v 1%)", 0.05, 0.01, 200), ("common illness (50% v 10%)", 0.50, 0.10, 1000)]
for label, risk_e, risk_u, size in scenarios:
    people = make_population(20000, 0.3, risk_e, risk_u)
    rr, rr_lo, rr_hi = cohort(people)
    orr, or_lo, or_hi = case_control(people, size, size)
    true_or = (risk_e / (1 - risk_e)) / (risk_u / (1 - risk_u))
    print(f"{label:27s} {size:5d}+{size:<5d} {risk_e / risk_u:6.1f} {rr:8.2f} ({rr_lo:.2f} to {rr_hi:.2f}) {true_or:8.2f} {orr:10.2f} ({or_lo:.2f} to {or_hi:.2f})")
# scenario                     study size  true RR  cohort RR (95% interval)  true OR  case-control OR (95% interval)
# rare illness (5% v 1%)        200+200      5.0     4.60 (3.77 to 5.61)     5.21       4.53 (2.97 to 6.90)
# common illness (50% v 10%)   1000+1000     5.0     4.96 (4.69 to 5.24)     9.00       8.44 (6.88 to 10.35)
    

In the rare case the two designs agree with the truth within their intervals: the cohort risk ratio is 4.60 (3.77 to 5.61), against a true 5.0, and the case-control odds ratio is 4.53 (2.97 to 6.90), against a true odds ratio of 5.21 calculated by hand. The case-control interval is much wider, even though the cohort followed 20,000 people and the case-control study interviewed only 400. In the common case the story changes. The cohort still recovers the risk ratio (4.96), but the case-control study estimates 8.44 (6.88 to 10.35), close to the true odds ratio of 9.00 and well away from the risk ratio of 5. The odds ratio is a different quantity from the risk ratio. It is close to the risk ratio only when the illness is rare, a standard result in epidemiology that this run illustrates, not proves. The practical lesson: a case-control odds ratio from a common illness should not be quoted as “times more likely”.

	Cohort study	Case-control study
Starts from	exposure	illness
Measures	risk, so risk ratio	odds of exposure, so odds ratio
People needed	everyone at risk (20,000 here)	cases plus a few controls (400 or 2,000 here)
Good for	common exposures, several outcomes	rare illnesses, quick outbreak investigations
Main risks	cost, loss to follow-up	choosing controls, recall of past exposure
Rare illness (RR 5 true)	4.60	4.53 (true OR 5.21)
Common illness (RR 5 true)	4.96	8.44 (true OR 9.00)

Bias in Testing

So far the counts were correct. Real counts are not. A sensitive test (sensitivity: the share of infected people who test positive) misses few cases; a specific test (specificity: the share of uninfected people who test negative) raises few false alarms. Even a decent test can mislead when the infection is uncommon. Worse, the people tested are rarely a random sample: those with symptoms, or contacts of a known case, or people with the time to queue, are tested more. The code makes a town of 100,000 with 2 per cent truly infected, a test with 80 per cent sensitivity and 98 per cent specificity, and two plans. Plan A tests 5,000 random people. Plan B tests only people who come forward, which we assume 40 per cent of infected people do against 3 per cent of the rest. Each plan reports the share positive and a corrected figure from the Rogan-Gladen formula, (observed + specificity − 1) / (sensitivity + specificity − 1), which removes the effect of test errors only.

# Testing bias: random testing versus people who choose to be tested. Every number is assumed.
import random

random.seed(5)
n, prevalence = 100000, 0.02          # assumed true share infected
sens, spec = 0.80, 0.98               # assumed test sensitivity and specificity
infected = [random.random() < prevalence for _ in range(n)]

def positive(is_infected):
    """Test result: infected people test positive with chance sens, others with chance 1 - spec."""
    return random.random() < (sens if is_infected else 1 - spec)

def rogan_gladen(observed):
    """Correct for test errors only: (observed + spec - 1) / (sens + spec - 1)."""
    return (observed + spec - 1) / (sens + spec - 1)

# Plan A: test 5,000 people chosen at random.
sample = random.sample(range(n), 5000)
observed_a = sum(positive(infected[i]) for i in sample) / len(sample)

# Plan B: test only people who come forward (40% of infected people do, 3% of the rest).
tested = [i for i in range(n) if random.random() < (0.40 if infected[i] else 0.03)]
observed_b = sum(positive(infected[i]) for i in tested) / len(tested)

print(f"true prevalence                       {sum(infected) / n:6.2%}")
print(f"A random sample ({len(sample)} tests):   positive {observed_a:6.2%}  corrected {rogan_gladen(observed_a):6.2%}")
print(f"B volunteers ({len(tested)} tests):     positive {observed_b:6.2%}  corrected {rogan_gladen(observed_b):6.2%}")
# true prevalence                        1.93%
# A random sample (5000 tests):   positive  3.24%  corrected  1.59%
# B volunteers (3716 tests):     positive 18.25%  corrected 20.83%
    

The true share infected is 1.93 per cent. Plan A sees 3.24 per cent positive, about 1.7 times the truth, because at this low prevalence the false alarms from the uninfected outnumber the real finds: by my arithmetic, using the expected values, roughly 56 per cent of positive results here are false. The corrected figure, 1.59 per cent, is close but not exact, and the gap is plausibly sampling noise (about 160 positives out of 5,000 tests, so a small change in the count moves the result by a few tenths of a point; my reasoning, not tested here). Plan B is badly off: 18.25 per cent positive, and 20.83 per cent after correction, more than ten times the truth. The formula repaired the test errors but could not repair the selection, since infected people were far more likely to be tested. A bigger sample would not fix this either: Plan B’s sample is large (3,716 tests) and still wrong. The common remedy is a random survey, even a small one; the size and the randomness are separate questions, and randomness is the one that matters here.

Source of error	What goes wrong	Does it fix itself with more tests?	Remedy
Test errors	false positives and missed cases distort the share positive	no	correct with sensitivity and specificity
Who gets tested	the tested are more often the infected (Plan B)	no	random sampling; report who was tested
Changing test effort	more tests find more cases with no change in infections (Exercise 2)	makes it worse	report the share positive as well as counts
Result in this run	truth 1.93%	A 3.24%, 1.59% corrected	B 18.25%, 20.83% corrected

Common mistakes

Using exp(r × mean interval) as if it were exact (it overestimates when the interval varies); assuming the serial interval is known when it is itself an estimate from a handful of interviews; fitting a growth rate on days when interventions or immunity have already changed the epidemic, so that the answer is no longer the R you wanted; reading a case-control odds ratio as a risk ratio when the illness is common; choosing controls who differ from the cases in ways other than the exposure; trusting a share-positive figure without asking who was tested and how good the test is; treating a rise in reported cases as a rise in infections when testing has also risen; correcting for test errors and believing that also corrects for who got tested; reporting a single run instead of a spread over many; and using the illustrative numbers here (5-day interval, 2 per cent prevalence, 80 and 98 per cent, 20,000 people) as estimates for a real disease.

Exercise 1

Simulate 100 outbreaks with a true R of 2.0 and the 5-day gamma serial interval (assumed), over 30 days. Estimate R from the growth rate in five windows (days 0 to 10, 5 to 15, 10 to 20, 15 to 25, 20 to 30), then, for the window 10 to 20, using assumed serial-interval standard deviations of 1.0, 2.5 and 4.0 days with the mean fixed at 5. Which window is least noisy here, and what would make a later window misleading in a real epidemic?

📄 View solution

Exercise 2

Keep the true share infected fixed at 1 per cent (assumed) but let the number of daily tests double every 14 days, from 2,000 tests on day 0, over 28 days, with 80 per cent sensitivity and 98 per cent specificity (assumed). Estimate the growth rate and the apparent R from the raw positive counts, then from the share of tests that are positive. Which is misleading, and by how much?

📄 View solution

Exercise 3

In a town of 20,000 (30 per cent exposed, 5 per cent of the exposed and 1 per cent of the unexposed ill, assumed), repeat a case-control study with 100 cases 500 times, using 1, 2, 4 and 10 controls per case. Print the mean odds ratio and the spread of its logarithm for each. How much does adding controls help, and what limits it?

📄 View solution

Quick Reference

Growth rate r: slope of the log of daily cases, fitted over a window where nothing else changes
Serial interval: symptom onset to symptom onset between infector and infected; stands in for the generation interval (infection to infection)
From r to R: R = exp(r × mean) only if the interval never varies; for a gamma interval R = (1 + r × scale) ^ shape (my derivation, checked by simulation)
Cohort study: follow exposed and unexposed, compare risks; risk ratio = risk if exposed / risk if unexposed
Case-control study: compare exposure in cases and controls; odds ratio approximates the risk ratio only if the illness is rare
Sensitivity / specificity: share of infected who test positive / share of uninfected who test negative
Rogan-Gladen: true share = (observed + specificity − 1) / (sensitivity + specificity − 1); fixes test errors, not selection
Testing bias: volunteers, symptomatic people and rising test numbers distort counts; a random sample and the share positive help
Assumes: constant R during the fitting window, known test accuracy, illustrative numbers

Coming next

How Epidemics Spread 9 is about controls and interventions: quarantine, contact tracing, closures, timing effects, comparing strategies by simulation.