About the data Every record and result in this guide is simulated. The examples are reproducible and contain no identifiable health information. Simulated associations are teaching devices, not evidence about a real population.
Each chapter follows the same rhythm: why the idea matters, the core concept, a public health example, R code, and an interpretation check. Read the explanation first, run or adapt the code, and then open the answer panels. Code is visible by default and can be collapsed with the Code buttons.
You do not need to memorize every formula. Focus on three habits:
By the end of the module, you should be able to:
Biostatistics applies statistical reasoning to questions about health, disease, health services, and populations. It is not only a collection of tests. It is a disciplined path from a question to evidence and then to an appropriately cautious decision.
What health events, exposures, and patterns are present?
How large are differences between meaningful groups or time periods?
What can a sample tell us about a target population?
A useful workflow is:
Question → design → measurement → data check → description → estimation → uncertainty → interpretation → action
Statistics cannot repair a vague question, biased sampling, or poor measurement. Those decisions happen before a p-value is calculated.
Suppose a city surveys 520 adults to estimate influenza vaccination coverage.
| Term | Meaning | Example |
|---|---|---|
| Target population | The full group to which the question refers | All non-institutionalized adults living in the city this season |
| Sample | The people actually observed | The 520 survey participants |
| Parameter | A fixed but usually unknown population quantity | The true vaccination proportion among all target adults |
| Statistic / estimate | A quantity calculated from the sample | The observed vaccination proportion in the survey |
| Estimand | A precise description of what is to be estimated | Vaccination coverage during this season in the target population |
The unit of observation is also essential. A row might represent a person, household, clinic, neighborhood, or person-day. The method must respect that unit and any clustering.
Before opening R, write down:
Worked question Among adults living in the city at the start of winter (population), what is the 12-week risk of laboratory-confirmed respiratory infection (outcome) among vaccinated adults compared with unvaccinated adults (comparison), expressed as a risk ratio (estimand)?
In a random sample of 800 high-school students, 18% report current vaping. What is the parameter?
Answer: The parameter is the true proportion of students in the defined target population who currently vape. The observed 18% is a sample statistic. Whether it estimates the parameter well depends on sampling, measurement, and nonresponse.A variable’s role depends on the question.
| Variable type | What it records | Public health example | Common summary |
|---|---|---|---|
| Nominal categorical | Unordered labels | Neighborhood, clinic | Count and proportion |
| Ordinal categorical | Ordered labels with unknown spacing | Self-rated health: poor to excellent | Count, proportion, median category |
| Binary | Two categories | Vaccinated: yes/no | Count and proportion |
| Discrete count | Non-negative event count | Emergency visits in one year | Mean, variance, rate; often a count model |
| Continuous | Measured numeric value | Blood pressure, age, BMI | Mean/SD or median/IQR |
| Time-to-event | Time until an event or censoring | Time to relapse | Survival probability, hazard, median survival |
Do not choose a summary from the storage format alone. A postal code stored as a number is still categorical, and a five-point rating is usually ordinal rather than truly continuous.
| Design | How it starts | Useful for | Main caution |
|---|---|---|---|
| Cross-sectional survey | Sample people at one time or short period | Prevalence and current patterns | Exposure and outcome timing may be unclear |
| Cohort study | Group people by exposure and follow outcomes | Incidence, risk, rates, temporal order | Loss to follow-up and confounding |
| Case-control study | Sample cases and controls, then assess prior exposure | Rare outcomes or long latency | Selection/recall bias; risk is not directly estimated from the sampled table |
| Randomized trial | Randomly assign an intervention | Intervention effects under the trial conditions | Adherence, loss to follow-up, ethics, and generalizability |
| Cluster randomized trial | Randomize schools, clinics, or communities | Population-level programs | Outcomes within a cluster are correlated |
| Surveillance system | Continuously collect defined events | Trends, signals, program monitoring | Changes in case definition, testing, and completeness |
Design determines interpretation A cross-sectional association cannot establish which variable came first. A case-control odds ratio is valid under its sampling design, but the sampled case proportion is not disease prevalence. Clustered data need methods that account for within-cluster similarity.
These are different problems:
| Issue | Plain-language meaning | Example | Typical response |
|---|---|---|---|
| Chance | Sample results vary from one random sample to another | Coverage is 68% in one sample and 71% in another | Quantify with standard errors and confidence intervals |
| Selection bias | Inclusion is related to variables important to the question | An online survey misses residents without internet access | Improve sampling/recruitment; assess nonresponse |
| Information bias | Exposure or outcome is measured differently or inaccurately | Cases remember a past exposure more completely than controls | Use valid, standardized, blinded measurement when possible |
| Confounding | A common cause creates or distorts an exposure–outcome association | Age affects both vaccination uptake and infection risk | Use design and analysis informed by subject matter and causal reasoning |
Confounding is not simply “any third variable” and is not diagnosed solely by a small p-value. Adjustment helps only when the covariates, model, measurements, and causal assumptions are appropriate.
Researchers select 300 people with lung cancer and 300 people without lung cancer and ask about occupational asbestos exposure 30 years earlier.
Answer: This is a case-control study. Differential recall or reconstruction of old exposure histories could cause information bias; selection of controls could also create selection bias. The odds ratio is the natural association measure from the sampled 2×2 table.The first rows of the simulated teaching dataset are shown below. Never begin with a test before checking variable names, units, categories, impossible values, duplicates, and missingness.
## [1] 520 10
head(
ph_data[c(
"participant_id", "age", "neighborhood", "smoking_status",
"physical_activity_min_week", "systolic_bp"
)],
6
)| participant_id | age | neighborhood | smoking_status | physical_activity_min_week | systolic_bp |
|---|---|---|---|---|---|
| P001 | 20 | North | Not current | 80 | 114 |
| P002 | 27 | South | Current | 81 | 110 |
| P003 | 58 | Central | Not current | 76 | 130 |
| P004 | 28 | Central | Not current | 6 | 116 |
| P005 | 35 | Central | Not current | 284 | 98 |
| P006 | 64 | West | Not current | 68 | 124 |
## participant_id age neighborhood
## 0 0 0
## smoking_status physical_activity_min_week bmi
## 0 18 0
## systolic_bp access_to_care vaccinated
## 0 0 0
## respiratory_infection
## 0
The missing physical-activity values are coded as NA,
not zero. Zero means a measured value of no activity; NA
means the value is unknown.
Use mean with SD for a reasonably symmetric distribution when those quantities answer the question. Use median with IQR for a strongly skewed distribution or when a typical rank is more meaningful. Always look at the distribution.
summary_table <- data.frame(
Variable = c("Age (years)", "Systolic blood pressure (mmHg)",
"Physical activity (min/week)"),
N_observed = c(sum(!is.na(ph_data$age)),
sum(!is.na(ph_data$systolic_bp)),
sum(!is.na(ph_data$physical_activity_min_week))),
Mean = c(mean(ph_data$age), mean(ph_data$systolic_bp),
mean(ph_data$physical_activity_min_week, na.rm = TRUE)),
SD = c(sd(ph_data$age), sd(ph_data$systolic_bp),
sd(ph_data$physical_activity_min_week, na.rm = TRUE)),
Median = c(median(ph_data$age), median(ph_data$systolic_bp),
median(ph_data$physical_activity_min_week, na.rm = TRUE)),
IQR = c(IQR(ph_data$age), IQR(ph_data$systolic_bp),
IQR(ph_data$physical_activity_min_week, na.rm = TRUE))
)
knitr::kable(summary_table, digits = 1,
caption = "Descriptive statistics for simulated participant data")| Variable | N_observed | Mean | SD | Median | IQR |
|---|---|---|---|---|---|
| Age (years) | 520 | 44 | 15.5 | 44 | 23 |
| Systolic blood pressure (mmHg) | 520 | 124 | 14.4 | 123 | 19 |
| Physical activity (min/week) | 502 | 119 | 97.3 | 97 | 112 |
old_par <- par(mfrow = c(1, 2), mar = c(4.2, 4.2, 2.2, 1), las = 1)
hist(
ph_data$systolic_bp,
breaks = "FD",
col = palette_ph["sky"],
border = "white",
main = "Distribution of systolic BP",
xlab = "Systolic blood pressure (mmHg)",
ylab = "Participants"
)
boxplot(
systolic_bp ~ smoking_status,
data = ph_data,
col = c(palette_ph["teal"], palette_ph["orange"]),
border = palette_ph["navy"],
main = "Systolic BP by smoking status",
xlab = "Smoking status",
ylab = "Systolic blood pressure (mmHg)"
)Systolic blood pressure is roughly unimodal, and its distribution differs by current smoking status in this simulated sample.
Graphical integrity Label units, show the denominator behind percentages, preserve time order, and explain exclusions. A truncated vertical axis can exaggerate small differences. Three-dimensional effects and rainbow palettes add visual noise and can reduce accessibility.
Always state the population, time window, numerator, and denominator. “The rate was 12%” is ambiguous because percentages usually describe proportions, while rates use person-time units.
vaccination_counts <- table(ph_data$vaccinated)
vaccination_summary <- data.frame(
Vaccinated = names(vaccination_counts),
Count = as.vector(vaccination_counts),
Proportion = as.vector(vaccination_counts) / sum(vaccination_counts)
)
knitr::kable(vaccination_summary, digits = 3,
caption = "Vaccination status in the simulated sample")| Vaccinated | Count | Proportion |
|---|---|---|
| No | 151 | 0.29 |
| Yes | 369 | 0.71 |
Emergency-department waiting time is strongly right-skewed because a few patients wait many hours. Which pair is most informative: mean/SD or median/IQR?
Answer: Median and IQR are usually more representative of the center and spread of a strongly skewed distribution. Reporting selected percentiles or the mean as an additional policy-relevant measure can still be useful if clearly labeled.A probability ranges from 0 (impossible) to 1 (certain). For an event :
A conditional probability, , is the probability of among units for which is true. This denominator matters. Sensitivity, for example, is the probability of a positive test among people who truly have the condition.
Events and are independent when learning that one occurred does not change the probability of the other:
Independence is an assumption to justify, not a default. Repeated observations from one person and residents of the same household are often correlated.
| Distribution | What it can represent | Key conditions |
|---|---|---|
| Binomial | Number infected among a fixed number of similarly exposed people | Fixed trials, two outcomes, common probability, independent trials |
| Normal | Symmetric continuous measurements or sampling distributions | Bell-shaped approximation with suitable context |
| Poisson | Event counts in a specified amount of time or space | Events occur at an approximately stable rate; mean and variance relationship may be restrictive |
Real data may violate these conditions. For example, infectious cases cluster and can show more variability than a simple Poisson model allows.
old_par <- par(mfrow = c(1, 3), mar = c(4, 3.8, 2.4, 0.8), las = 1)
x_bin <- 0:12
plot(x_bin, dbinom(x_bin, size = 12, prob = 0.20), type = "h", lwd = 6,
lend = 1, col = palette_ph["blue"],
xlab = "Cases among 12", ylab = "Probability", main = "Binomial")
x_norm <- seq(-3.5, 3.5, length.out = 300)
plot(x_norm, dnorm(x_norm), type = "l", lwd = 3,
col = palette_ph["teal"],
xlab = "Standardized value", ylab = "Density", main = "Normal")
x_pois <- 0:14
plot(x_pois, dpois(x_pois, lambda = 4), type = "h", lwd = 6,
lend = 1, col = palette_ph["vermillion"],
xlab = "Events per day", ylab = "Probability", main = "Poisson")Examples of binomial, normal, and Poisson probability models. Each model answers a different kind of question.
If we repeatedly draw random samples and calculate an estimate, the estimates form a sampling distribution. Its standard deviation is the estimate’s standard error (SE). For independent observations from a common population with finite variance, an estimated SE of the sample mean is:
The SD describes variability among observations. The SE describes uncertainty in an estimate. They are not interchangeable.
The sampling distribution is unknown because we have only one sample. A nonparametric bootstrap approximates it by repeatedly resampling the observed data, with replacement, using the original sample size. The empirical sample temporarily stands in for the unknown population; this approximation does not correct a biased or unrepresentative sample.
set.seed(410)
bootstrap_means <- replicate(
1000,
mean(sample(ph_data$systolic_bp, size = nrow(ph_data), replace = TRUE))
)
hist(bootstrap_means, breaks = 24, col = palette_ph["sky"], border = "white",
main = "Bootstrap distribution of the mean",
xlab = "Mean systolic BP in bootstrap samples",
ylab = "Bootstrap samples")
abline(v = mean(ph_data$systolic_bp), col = palette_ph["vermillion"],
lwd = 3, lty = 2)
legend("topright", legend = "Observed sample mean", lty = 2, lwd = 3,
col = palette_ph["vermillion"], bty = "n")The bootstrap distribution of the sample mean is centered near the observed mean. The dashed line marks the observed sample mean.
Because SE usually decreases in proportion to , quadrupling the sample size approximately halves the SE—not quarters it. The simple formula does not apply unchanged to clustered, paired, weighted, or otherwise dependent observations; the design must be reflected in the SE.
se_example <- data.frame(
Sample_size = c(25, 100, 400),
Assumed_SD = 16,
Standard_error = 16 / sqrt(c(25, 100, 400))
)
knitr::kable(se_example, digits = 1,
caption = "How sample size changes the standard error when SD = 16")| Sample_size | Assumed_SD | Standard_error |
|---|---|---|
| 25 | 16 | 3.2 |
| 100 | 16 | 1.6 |
| 400 | 16 | 0.8 |
A report says, “Average systolic blood pressure was 128 mmHg (SD 17).” What does 17 describe?
Answer: It describes variability among participant blood-pressure values around their sample mean. It does not directly describe uncertainty in the estimated population mean; that role belongs to the standard error.A point estimate is a single best estimate from the data. An interval estimate communicates sampling uncertainty. A 95% confidence interval (CI) combines an estimate with its standard error and an appropriate reference distribution.
For a mean under common conditions:
# A t-based confidence interval for the mean systolic blood pressure.
n_bp <- sum(!is.na(ph_data$systolic_bp))
mean_bp <- mean(ph_data$systolic_bp, na.rm = TRUE)
se_bp <- sd(ph_data$systolic_bp, na.rm = TRUE) / sqrt(n_bp)
ci_mean_bp <- mean_bp + c(-1, 1) * qt(0.975, df = n_bp - 1) * se_bp
# A score-based interval for the vaccination proportion.
x_vax <- sum(ph_data$vaccinated == "Yes")
n_vax <- sum(!is.na(ph_data$vaccinated))
vax_ci_object <- prop.test(x_vax, n_vax, correct = FALSE)
vax_estimate <- unname(vax_ci_object$estimate)
vax_ci <- vax_ci_object$conf.int
ci_table <- data.frame(
Quantity = c("Mean systolic BP (mmHg)", "Vaccination proportion"),
Estimate = c(mean_bp, vax_estimate),
Lower_95_CI = c(ci_mean_bp[1], vax_ci[1]),
Upper_95_CI = c(ci_mean_bp[2], vax_ci[2])
)
knitr::kable(ci_table, digits = 3,
caption = "Point and interval estimates from the simulated sample")| Quantity | Estimate | Lower_95_CI | Upper_95_CI |
|---|---|---|---|
| Mean systolic BP (mmHg) | 124.30 | 123.060 | 125.536 |
| Vaccination proportion | 0.71 | 0.669 | 0.747 |
The estimated mean systolic blood pressure is 124.3 mmHg (95% CI 123.1 to 125.5). Under the model and sampling assumptions, the interval describes uncertainty in the population mean—not the range containing 95% of individual blood-pressure values.
If we repeatedly used the same valid sampling and interval procedure, about 95% of the resulting intervals would contain the true parameter. After one interval is computed, the parameter is fixed; in the standard frequentist interpretation we do not say there is a 95% probability that this particular interval contains it.
CI width reflects random uncertainty under a model. It does not automatically include selection bias, measurement error, uncontrolled confounding, model misspecification, or data-processing mistakes.
Intervals are not pass/fail devices A narrow interval around a trivial effect can be unimportant; a wide interval may contain both meaningful benefit and harm. Interpret the estimate, interval limits, units, study design, and public health context together.
A risk difference is −4 percentage points with a 95% CI from −9 to +1 points. What is a careful interpretation?
Answer: The data are compatible with effects ranging from a meaningful reduction to a small increase, under the assumptions. Because the interval includes zero, the data do not reject a zero difference at the two-sided 5% level, but “no statistically significant difference” is not proof of no effect.A hypothesis test begins with a null hypothesis , often no difference or no association, and an alternative . A p-value is:
the probability, assuming the null hypothesis and model assumptions are true, of obtaining data at least as incompatible with the null as the observed data.
A p-value is not the probability that the null hypothesis is true, the probability that results occurred “by chance,” or the size of an effect.
Power depends on sample size, variability, effect size, outcome frequency, design, missingness, and the chosen significance threshold. It should be planned before data collection when possible.
Welch’s t-test is a sensible default for comparing two independent means because it does not assume equal group variances. Here it compares mean systolic blood pressure by current smoking status.
bp_test <- t.test(systolic_bp ~ smoking_status, data = ph_data)
bp_means <- aggregate(systolic_bp ~ smoking_status, data = ph_data, mean)
bp_means| smoking_status | systolic_bp |
|---|---|
| Not current | 123 |
| Current | 129 |
##
## Welch Two Sample t-test
##
## data: systolic_bp by smoking_status
## t = -3, df = 137, p-value = 0.001
## alternative hypothesis: true difference in means between group Not current and group Current is not equal to 0
## 95 percent confidence interval:
## -8.91 -2.28
## sample estimates:
## mean in group Not current mean in group Current
## 123 129
The estimated mean difference is 5.6 mmHg when expressed as current minus not-current smoking. The two-sided p-value is 0.0011. This simulated, observational comparison is an association; it does not show that smoking caused the difference, because other characteristics may differ between groups.
| Question | Common introductory method | Important check |
|---|---|---|
| Compare two independent means | Welch two-sample t-test | Independent units; distribution and influential observations |
| Compare paired before/after measurements | Paired t-test on within-person differences | Pairing is correct; distribution of differences |
| Compare two categorical variables | Chi-square test of independence | Expected cell counts; use Fisher’s exact test when sparse |
| Compare ordered or highly skewed outcomes | Rank-based method when appropriate | A rank-sum test is not automatically a test of medians |
| Estimate an adjusted continuous association | Linear regression | Functional form, residual patterns, influential observations |
| Estimate an adjusted binary association | Logistic regression | Correct specification, sparse data, odds-ratio interpretation |
Statistical significance is not public health importance Large samples can make small effects yield small p-values. Small studies can miss important effects. Report the effect estimate, confidence interval, units, absolute risk when possible, assumptions, and consequences—not only whether p < 0.05.
Multiple testing also matters. If many hypotheses are tested, false-positive findings become more likely. Pre-specify primary questions, limit opportunistic testing, and use multiplicity methods when the scientific setting calls for them.
Which is correct: “There is a 3% probability that the null hypothesis is true,” or “If the null and model assumptions were true, results at least this incompatible with the null would occur about 3% of the time”?
Answer: The second statement. The first incorrectly treats the p-value as . A p-value does not supply that probability.Suppose 240 workers are followed through one respiratory-outbreak period.
| Ill | Not ill | Total | |
|---|---|---|---|
| Exposed | 120 | ||
| Unexposed | 120 |
The group risks are:
Three common association measures answer different questions:
a <- 48; b <- 72; c <- 24; d <- 96
risk_exposed <- a / (a + b)
risk_unexposed <- c / (c + d)
risk_difference <- risk_exposed - risk_unexposed
risk_ratio <- risk_exposed / risk_unexposed
odds_ratio <- (a * d) / (b * c)
# Large-sample confidence intervals for teaching purposes.
se_rd <- sqrt(risk_exposed * (1 - risk_exposed) / (a + b) +
risk_unexposed * (1 - risk_unexposed) / (c + d))
ci_rd <- risk_difference + c(-1, 1) * 1.96 * se_rd
se_log_rr <- sqrt(1 / a - 1 / (a + b) + 1 / c - 1 / (c + d))
ci_rr <- exp(log(risk_ratio) + c(-1, 1) * 1.96 * se_log_rr)
se_log_or <- sqrt(1 / a + 1 / b + 1 / c + 1 / d)
ci_or <- exp(log(odds_ratio) + c(-1, 1) * 1.96 * se_log_or)
effect_table <- data.frame(
Measure = c("Risk difference", "Risk ratio", "Odds ratio"),
Estimate = c(risk_difference, risk_ratio, odds_ratio),
Lower_95_CI = c(ci_rd[1], ci_rr[1], ci_or[1]),
Upper_95_CI = c(ci_rd[2], ci_rr[2], ci_or[2])
)
knitr::kable(effect_table, digits = 2,
caption = "Unadjusted association measures for the outbreak example")| Measure | Estimate | Lower_95_CI | Upper_95_CI |
|---|---|---|---|
| Risk difference | 0.20 | 0.09 | 0.31 |
| Risk ratio | 2.00 | 1.31 | 3.04 |
| Odds ratio | 2.67 | 1.50 | 4.75 |
Illness risk was 40% in exposed workers and 20% in unexposed workers. The estimated risk difference was 20 percentage points, or 20 additional cases per 100 workers over this outbreak period. The risk ratio was 2.0: exposed workers had twice the observed risk. The odds ratio was 2.67 and should not be mislabeled as a risk ratio.
These calculations are unadjusted. Confounding, clustering, censoring, unequal follow-up, and complex sampling can require other methods.
An intervention reduces risk from 2% to 1%. What are the risk difference and risk ratio?
Answer: The risk difference is −1 percentage point (1% − 2%); the risk ratio is 0.50 (1% / 2%), a 50% relative reduction. Both are correct but communicate different information.Consider a screening program for 1,000 people. A reference standard identifies 80 people with the condition. The screening test gives the following results:
| Condition present | Condition absent | Total | |
|---|---|---|---|
| Test positive | TP = 68 | FP = 92 | 160 |
| Test negative | FN = 12 | TN = 828 | 840 |
| Total | 80 | 920 | 1,000 |
tp <- 68; fp <- 92; fn <- 12; tn <- 828
screening <- data.frame(
Measure = c("Sensitivity", "Specificity", "Positive predictive value",
"Negative predictive value"),
Estimate = c(tp / (tp + fn), tn / (tn + fp),
tp / (tp + fp), tn / (tn + fn))
)
knitr::kable(screening, digits = 3,
caption = "Screening performance in the hypothetical program")| Measure | Estimate |
|---|---|
| Sensitivity | 0.850 |
| Specificity | 0.900 |
| Positive predictive value | 0.425 |
| Negative predictive value | 0.986 |
Sensitivity is 85% and specificity is 90%. Among those who test positive, however, only 42.5% have the condition (PPV). The difference occurs because the condition prevalence in this screened population is 8% and false positives accumulate among the much larger condition-absent group.
For fixed sensitivity and specificity, PPV generally rises as prevalence rises, while NPV generally falls.
prevalence_grid <- seq(0.01, 0.50, by = 0.01)
sens <- 0.85
spec <- 0.90
ppv_grid <- sens * prevalence_grid /
(sens * prevalence_grid + (1 - spec) * (1 - prevalence_grid))
npv_grid <- spec * (1 - prevalence_grid) /
((1 - sens) * prevalence_grid + spec * (1 - prevalence_grid))
plot(prevalence_grid, ppv_grid, type = "l", lwd = 3,
col = palette_ph["vermillion"], ylim = c(0, 1),
xlab = "Condition prevalence", ylab = "Predictive value",
main = "Predictive values depend on prevalence")
lines(prevalence_grid, npv_grid, lwd = 3, lty = 2,
col = palette_ph["blue"])
legend("right", legend = c("PPV", "NPV"), lwd = 3, lty = c(1, 2),
col = c(palette_ph["vermillion"], palette_ph["blue"]), bty = "n")Even with fixed sensitivity of 85% and specificity of 90%, predictive values change as condition prevalence changes.
A threshold that increases sensitivity usually reduces specificity, and vice versa. The preferred trade-off depends on the consequences of missed cases, false alarms, follow-up resources, and equity. Sensitivity and specificity can also vary across populations because of disease spectrum and verification processes.
To calculate sensitivity, should the denominator be all positive tests or all people who truly have the condition?
Answer: All people who truly have the condition, (TP+FN). The proportion of positive tests that are true positives is PPV, which uses (TP+FP) as its denominator.Pearson correlation summarizes the strength of a linear relationship between two numeric variables. Spearman correlation summarizes a monotonic rank relationship. Both can be altered by outliers, restricted ranges, mixed subgroups, and nonlinear patterns.
plot(ph_data$age, ph_data$systolic_bp,
pch = 16, cex = 0.65, col = rgb(0, 114/255, 178/255, 0.38),
xlab = "Age (years)", ylab = "Systolic blood pressure (mmHg)",
main = "Age and systolic blood pressure")
abline(lm(systolic_bp ~ age, data = ph_data),
col = palette_ph["vermillion"], lwd = 3)Systolic blood pressure tends to increase with age in the simulated data, but individuals vary substantially.
The Pearson correlation is 0.62. That value does not prove that changing age would change blood pressure by any particular amount. Causal interpretation requires a defensible design and assumptions beyond correlation.
A multiple linear model can describe mean systolic blood pressure as a function of several covariates:
bp_model <- lm(systolic_bp ~ age + bmi + smoking_status, data = ph_data)
coef_table <- summary(bp_model)$coefficients
knitr::kable(coef_table, digits = 3,
caption = "Linear regression of systolic blood pressure")| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 90.533 | 3.220 | 28.12 | 0.000 |
| age | 0.569 | 0.032 | 17.99 | 0.000 |
| bmi | 0.308 | 0.108 | 2.85 | 0.005 |
| smoking_statusCurrent | 3.532 | 1.265 | 2.79 | 0.005 |
The age coefficient is 0.57 mmHg per year. In this model, participants who differ by one year in age but have the same modeled BMI and smoking status differ in mean systolic blood pressure by about 0.57 mmHg. The coefficient is an adjusted association, not automatically a causal effect.
Checks should include residual patterns, linearity, influential observations, dependence, and whether a linear mean model makes scientific sense. A large sample does not rescue the wrong functional form.
Logistic regression models the log odds of a binary outcome. Exponentiating a coefficient gives an odds ratio, not a risk ratio.
infection_model <- glm(
respiratory_infection ~ vaccinated + age + smoking_status + neighborhood,
data = ph_data,
family = binomial()
)
model_coef <- summary(infection_model)$coefficients
non_intercept <- rownames(model_coef) != "(Intercept)"
or_table <- data.frame(
Term = rownames(model_coef)[non_intercept],
Odds_ratio = exp(model_coef[non_intercept, "Estimate"]),
Lower_95_CI = exp(model_coef[non_intercept, "Estimate"] -
1.96 * model_coef[non_intercept, "Std. Error"]),
Upper_95_CI = exp(model_coef[non_intercept, "Estimate"] +
1.96 * model_coef[non_intercept, "Std. Error"]),
P_value = model_coef[non_intercept, "Pr(>|z|)"],
row.names = NULL
)
knitr::kable(or_table, digits = 3,
caption = "Adjusted odds ratios for respiratory infection (intercept omitted)")| Term | Odds_ratio | Lower_95_CI | Upper_95_CI | P_value |
|---|---|---|---|---|
| vaccinatedYes | 0.510 | 0.324 | 0.802 | 0.004 |
| age | 1.000 | 0.986 | 1.014 | 0.979 |
| smoking_statusCurrent | 1.554 | 0.931 | 2.593 | 0.091 |
| neighborhoodNorth | 0.980 | 0.544 | 1.767 | 0.947 |
| neighborhoodSouth | 0.893 | 0.491 | 1.623 | 0.710 |
| neighborhoodWest | 1.014 | 0.546 | 1.883 | 0.966 |
The adjusted odds ratio comparing vaccinated with unvaccinated participants is 0.51 (95% CI 0.32 to 0.8). It describes a conditional association in this simulated model. It is not an adjusted risk ratio and does not by itself establish vaccine effectiveness.
Categorical terms compare the displayed factor level with its reference level; the age odds ratio is for a one-year difference. The exponentiated intercept is a baseline odds, not an odds ratio, so it is intentionally omitted from the table.
Adjustment is not magic Do not select covariates only because their univariable p-values are small. Pre-specify variables using subject-matter knowledge and a causal framework. Avoid adjusting for consequences of the exposure when estimating a total effect, and remember that unmeasured or poorly measured confounding can remain.
An association can genuinely differ across subgroups or contexts. This is effect modification or interaction, not necessarily bias. Assess it on a scientifically meaningful scale, pre-specify key subgroup questions, and report subgroup estimates with uncertainty rather than comparing which subgroup p-value crosses 0.05.
A logistic model gives an adjusted odds ratio of 0.70 for an intervention versus comparison. Can you say risk is 30% lower?
Answer: Not without additional information or a method that estimates risks. The odds are 30% lower because . The risk ratio may be similar when the outcome is uncommon, but odds and risks can differ substantially for common outcomes.A defensible analysis should let another analyst understand what was done and reproduce the result.
# Compact checks that should precede modeling.
data_quality <- data.frame(
Check = c("Rows", "Duplicate participant IDs", "Missing activity values",
"Age outside 18–85", "Systolic BP outside 70–250"),
Result = c(
nrow(ph_data),
sum(duplicated(ph_data$participant_id)),
sum(is.na(ph_data$physical_activity_min_week)),
sum(ph_data$age < 18 | ph_data$age > 85),
sum(ph_data$systolic_bp < 70 | ph_data$systolic_bp > 250)
)
)
knitr::kable(data_quality, caption = "Selected reproducible data-quality checks")| Check | Result |
|---|---|
| Rows | 520 |
| Duplicate participant IDs | 0 |
| Missing activity values | 18 |
| Age outside 18–85 | 0 |
| Systolic BP outside 70–250 | 0 |
Ask why values are missing, whether missingness differs across groups, and what assumptions an analysis makes. Complete-case analysis can reduce precision and introduce bias. Never silently convert missing values to zero or automatically delete rows without reporting the effect.
More advanced options include multiple imputation, inverse-probability weighting, likelihood-based methods, and sensitivity analysis. Their validity depends on assumptions that should be stated and examined.
Health data represent people and systems Protect privacy, avoid stigmatizing language, explain how categories were defined, and include communities in decisions about collection and reporting. A technically correct model can still cause harm when the question, labels, or use of results are inappropriate.
Race and ethnicity variables usually reflect social, political, historical, and measurement processes; they should not be treated as simple biological causes. When disparities are observed, investigate structural conditions, access, discrimination, environment, and measurement practices. Avoid publishing small cells that could identify people, especially when geography and rare conditions are combined.
Survey weights, stratification, clustering, repeated measurements, multilevel structures, time-to-event outcomes, spatial correlation, and causal estimands often require methods beyond this introductory guide. The correct next step is to recognize the structure and seek an appropriate method—not to ignore it.
A map shows exact home locations for six people with a rare infection. What is the primary concern?
Answer: Re-identification and privacy risk. Aggregate or mask geography, apply disclosure-control rules, and involve data stewards and affected communities before release. Public benefit does not remove the obligation to protect participants.Question: In the simulated community sample, is vaccination associated with a lower observed 12-week risk of respiratory infection?
For teaching, we treat the sample as a cohort observed over a common period. The analysis is unadjusted first, then interpreted alongside the adjusted logistic model above. Because the data are simulated and observational, this is not an estimate of real-world effectiveness.
vax_table <- with(
ph_data,
table(Vaccinated = vaccinated, Infection = respiratory_infection)
)
vax_table## Infection
## Vaccinated No Yes
## No 106 45
## Yes 303 66
risk_vaccinated <- vax_table["Yes", "Yes"] / sum(vax_table["Yes", ])
risk_unvaccinated <- vax_table["No", "Yes"] / sum(vax_table["No", ])
rd_vax <- risk_vaccinated - risk_unvaccinated
rr_vax <- risk_vaccinated / risk_unvaccinated
a_v <- unname(vax_table["Yes", "Yes"])
b_v <- unname(vax_table["Yes", "No"])
c_v <- unname(vax_table["No", "Yes"])
d_v <- unname(vax_table["No", "No"])
se_log_rr_vax <- sqrt(1 / a_v - 1 / (a_v + b_v) +
1 / c_v - 1 / (c_v + d_v))
ci_rr_vax <- exp(log(rr_vax) + c(-1, 1) * 1.96 * se_log_rr_vax)
se_rd_vax <- sqrt(
risk_vaccinated * (1 - risk_vaccinated) / (a_v + b_v) +
risk_unvaccinated * (1 - risk_unvaccinated) / (c_v + d_v)
)
ci_rd_vax <- rd_vax + c(-1, 1) * 1.96 * se_rd_vax
ci_risk_vaccinated <- prop.test(a_v, a_v + b_v, correct = FALSE)$conf.int
ci_risk_unvaccinated <- prop.test(c_v, c_v + d_v, correct = FALSE)$conf.int
mini_results <- data.frame(
Measure = c("Risk among vaccinated", "Risk among unvaccinated",
"Risk difference", "Risk ratio"),
Estimate = c(risk_vaccinated, risk_unvaccinated, rd_vax, rr_vax),
Lower_95_CI = c(ci_risk_vaccinated[1], ci_risk_unvaccinated[1],
ci_rd_vax[1], ci_rr_vax[1]),
Upper_95_CI = c(ci_risk_vaccinated[2], ci_risk_unvaccinated[2],
ci_rd_vax[2], ci_rr_vax[2])
)
knitr::kable(mini_results, digits = 3,
caption = "Unadjusted estimates and 95% confidence intervals in the simulated mini case study")| Measure | Estimate | Lower_95_CI | Upper_95_CI |
|---|---|---|---|
| Risk among vaccinated | 0.179 | 0.143 | 0.221 |
| Risk among unvaccinated | 0.298 | 0.231 | 0.375 |
| Risk difference | -0.119 | -0.202 | -0.036 |
| Risk ratio | 0.600 | 0.432 | 0.833 |
The two risks use score intervals. The risk-difference interval uses a simple large-sample normal approximation, and the risk-ratio interval uses a large-sample log approximation; these teaching methods may be unsuitable for sparse data.
risks <- c(Unvaccinated = risk_unvaccinated, Vaccinated = risk_vaccinated)
bar_positions <- barplot(
risks,
ylim = c(0, max(risks) * 1.30),
col = c(palette_ph["orange"], palette_ph["teal"]),
border = NA,
ylab = "Observed 12-week infection risk",
main = "Respiratory infection risk by vaccination status",
las = 1
)
text(bar_positions, risks, labels = pct(risks), pos = 3, font = 2)Observed infection risk was lower among vaccinated participants in the simulated dataset.
Among 369 vaccinated participants, the observed infection risk was 17.9%, compared with 29.8% among 151 unvaccinated participants. The unadjusted risk ratio was 0.6 (95% CI 0.43 to 0.83), and the risk difference was -11.9 percentage points (95% CI -20.2 to -3.6). In this simulated sample, vaccination was associated with lower infection risk. Because vaccination was not randomized, differences in access, age, neighborhood, health behavior, or other factors may affect the association; this result should not be interpreted as proof of causation.
Among [population], the estimated [outcome] was [estimate] in group 1 and [estimate] in group 0. The estimated [risk difference/risk ratio/mean difference] was [value] with a 95% CI of [lower, upper]. This suggests [plain-language interpretation with units and time]. Because the data were [design], [specific limitation] may affect the estimate and causation should not be assumed.
Why might the unadjusted risk ratio and adjusted odds ratio differ?
Answer: They are different measures (risk ratio versus odds ratio), and the adjusted model conditions on age, smoking status, and neighborhood. Differences can therefore reflect both the measure’s scale and covariate adjustment. Model specification, sparse subgroups, and sampling variation can also contribute.| Quantity | Formula | Interpretation |
|---|---|---|
| Prevalence | Existing cases / defined population | Proportion with a condition at a specified time or period |
| Incidence proportion | New cases / population initially at risk | Risk over a defined period |
| Incidence rate | New cases / person-time at risk | Speed of new-event occurrence |
| Risk difference | Absolute excess or reduction | |
| Risk ratio | Relative risk | |
| Odds ratio | in a 2×2 table | Ratio of group odds |
| Estimated SE of a mean | Estimated sampling variability of a mean | |
| Approximate SE of a proportion | Simple large-sample uncertainty approximation | |
| Sensitivity | Positive test among those with the condition | |
| Specificity | Negative test among those without the condition | |
| PPV | Condition present among positive tests | |
| NPV | Condition absent among negative tests |
For small samples or proportions close to 0 or 1, basic Wald intervals can perform poorly. Score/Wilson or exact methods are often preferable, depending on the goal.
Before choosing a method, ask:
| Term | Meaning |
|---|---|
| Bias | A systematic tendency for an estimate to differ from the target quantity |
| Confidence interval | A range produced by a procedure designed to quantify sampling uncertainty |
| Confounder | A variable that can create or distort an association because of the causal structure |
| Estimate | A sample-based approximation of a target parameter |
| Estimand | The precisely defined quantity the analysis aims to estimate |
| Generalizability | How well findings apply to a target population or setting |
| Precision | How little random variability an estimate has; often reflected by CI width |
| p-value | A measure of data compatibility with a null hypothesis under model assumptions |
| Standard deviation | Variation among observed values |
| Standard error | Estimated sampling variability of a statistic |
Next topics for public health training include sample-size planning, survey weights and design effects, generalized linear models, longitudinal and multilevel data, survival analysis, spatial methods, causal inference, Bayesian methods, missing-data methods, and transparent reporting guidelines.