# How Epidemics Spread 8, Exercise 3
# Exercise 3: how many controls per case? Repeat a case-control study 500 times (all numbers assumed).
import math
import random
import statistics

random.seed(51)
n = 20000
people = []
for _ in range(n):
    exposed = random.random() < 0.3
    ill = random.random() < (0.05 if exposed else 0.01)
    people.append((exposed, ill))
ill_people = [p for p in people if p[1]]
well_people = [p for p in people if not p[1]]

def odds_ratio(cases, controls):
    a = sum(1 for e, i in cases if e)
    b = len(cases) - a
    c = sum(1 for e, i in controls if e)
    d = len(controls) - c
    return (a * d) / (b * c)

whole = odds_ratio(ill_people, well_people)
print(f"odds ratio using everyone: {whole:.2f}  ({len(ill_people)} ill people)")
print("100 cases; controls per case   mean OR   sd of log OR   studies with OR below 2 or above 12")
for per_case in [1, 2, 4, 10]:
    logs = []
    for _ in range(500):
        cases = random.sample(ill_people, 100)
        controls = random.sample(well_people, 100 * per_case)
        logs.append(math.log(odds_ratio(cases, controls)))
    odd = sum(1 for x in logs if math.exp(x) < 2 or math.exp(x) > 12)
    print(f"{per_case:28d} {statistics.mean(math.exp(x) for x in logs):10.2f} {statistics.stdev(logs):14.3f} {odd:36d}")

# Output when run:
# odds ratio using everyone: 4.90  (437 ill people)
# 100 cases; controls per case   mean OR   sd of log OR   studies with OR below 2 or above 12
#                            1       5.29          0.301                                    4
#                            2       5.16          0.247                                    0
#                            4       5.09          0.221                                    0
#                           10       4.99          0.198                                    0

# Reading the output: using everyone, the odds ratio is 4.90. With 100 cases and a single control per case the studies scatter widely (sd of the log odds ratio 0.301, and 4 of 500 studies landed below 2 or above 12). Adding controls helps, but with diminishing returns: 0.301, 0.247, 0.221 and 0.198 for 1, 2, 4 and 10 controls per case. The floor is set by the number of cases (only 100), so the better lever is more cases, not unlimited controls. The mean odds ratio drifts from 5.29 (1 control) to 4.99 (10 controls); a small-sample effect is my guess and was not tested. All sizes and risks are assumed.