1. Background and
Research Question
This example demonstrates a typical epidemiologic
effect-estimation workflow using a cohort-study framework.
Research question:
Among adults initially free of cardiovascular disease, is baseline
cigarette smoking associated with a higher risk of incident
cardiovascular disease during follow-up?
The causal question can be written conceptually as:
\[
Smoking \rightarrow Cardiovascular\ Disease
\]
Because smoking is not randomly assigned, the exposed and unexposed
groups may differ in important baseline characteristics. Therefore,
confounding must be considered carefully.
library(dplyr)
library(survival)
library(tableone)
library(MatchIt)
library(cobalt)
set.seed(2026)
2. Study Design and
Data Source
A common epidemiologic design for this question is a
prospective or retrospective cohort study.
In a real epidemiologic study, data might come from:
- a population-based cohort,
- an epidemiologic survey linked to outcomes,
- EHR or claims data,
- a disease registry,
- or linked administrative databases.
For teaching purposes, we simulate a cohort with the same basic
structure.
n <- 6000
epi <- data.frame(
id = 1:n,
age = round(rnorm(n, 55, 10)),
male = rbinom(n, 1, 0.48),
bmi = round(rnorm(n, 28, 5), 1),
systolic_bp = round(rnorm(n, 128, 16)),
diabetes = rbinom(n, 1, 0.14),
hypertension = rbinom(n, 1, 0.38),
high_cholesterol = rbinom(n, 1, 0.32),
physical_inactivity = rbinom(n, 1, 0.35)
)
# Smoking probability depends on several baseline characteristics.
# This intentionally creates confounding.
lp_smoking <-
-1.1 +
0.012 * (epi$age - 55) +
0.25 * epi$male +
0.18 * epi$physical_inactivity +
0.15 * epi$hypertension -
0.08 * (epi$bmi - 28)
p_smoking <- plogis(lp_smoking)
epi$smoking <- rbinom(n, 1, p_smoking)
epi$smoking <- factor(
epi$smoking,
levels = c(0, 1),
labels = c("Non-smoker", "Smoker")
)
table(epi$smoking)
##
## Non-smoker Smoker
## 4183 1817
3. Define the Target
Population and Study Cohort
The cohort should be defined before outcome analysis.
Example eligibility criteria:
- adults aged 35 to 80 years,
- no cardiovascular disease at baseline,
- baseline smoking status available,
- baseline covariates available,
- eligible for follow-up.
For this simulated example, we restrict the cohort to plausible
baseline values.
cohort <- epi %>%
filter(
age >= 35,
age <= 80,
bmi >= 15,
bmi <= 55,
systolic_bp >= 80,
systolic_bp <= 220
)
nrow(cohort)
## [1] 5825
table(cohort$smoking)
##
## Non-smoker Smoker
## 4057 1768
4. Define Exposure,
Outcome, and Time Zero
Exposure
The exposure is baseline cigarette smoking:
Outcome
The outcome is incident cardiovascular disease (CVD)
during follow-up.
Time zero
Time zero is the baseline study visit or cohort entry date.
This is important because exposure and baseline covariates should be
defined before follow-up begins.
5. Identify Potential
Confounders
Potential confounders should be selected using:
- subject-matter knowledge,
- previous epidemiologic evidence,
- temporal ordering,
- and causal reasoning / DAGs.
For illustration, we consider:
- age,
- sex,
- BMI,
- systolic blood pressure,
- diabetes,
- hypertension,
- high cholesterol,
- physical inactivity.
A confounder is generally a variable that is associated with both the
exposure and the outcome and is not a downstream consequence of the
exposure for the causal effect being estimated.
The goal is not to include every available variable
automatically.
6. Generate Follow-up
and Outcome
We now simulate time to incident CVD.
Smoking truly increases the event hazard in this simulated dataset.
Several baseline covariates also affect CVD risk.
lp_event <-
-4.2 +
0.45 * (cohort$smoking == "Smoker") +
0.035 * (cohort$age - 55) +
0.30 * cohort$male +
0.025 * (cohort$bmi - 28) +
0.012 * (cohort$systolic_bp - 128) +
0.55 * cohort$diabetes +
0.40 * cohort$hypertension +
0.30 * cohort$high_cholesterol +
0.18 * cohort$physical_inactivity
event_rate <- exp(lp_event)
true_event_time <- rexp(
nrow(cohort),
rate = event_rate
)
censor_time <- runif(
nrow(cohort),
min = 2,
max = 10
)
cohort$followup_time <- pmin(true_event_time, censor_time)
cohort$cvd_event <- as.integer(true_event_time <= censor_time)
table(cohort$cvd_event)
##
## 0 1
## 4790 1035
7. Data Quality
Assessment
Before estimating an exposure effect, assess the quality of the
epidemiologic data.
Important checks include:
- duplicate participants,
- missing data,
- implausible values,
- exposure misclassification,
- outcome misclassification,
- temporal consistency,
- and completeness of follow-up.
# Duplicate participant IDs
sum(duplicated(cohort$id))
## [1] 0
# Missingness
colSums(is.na(cohort))
## id age male bmi
## 0 0 0 0
## systolic_bp diabetes hypertension high_cholesterol
## 0 0 0 0
## physical_inactivity smoking followup_time cvd_event
## 0 0 0 0
# Basic distributions
summary(
cohort[, c(
"age",
"bmi",
"systolic_bp",
"followup_time"
)]
)
## age bmi systolic_bp followup_time
## Min. :35.00 Min. :15.00 Min. : 80.0 Min. :0.002586
## 1st Qu.:49.00 1st Qu.:24.70 1st Qu.:118.0 1st Qu.:3.212607
## Median :55.00 Median :28.00 Median :129.0 Median :5.135016
## Mean :55.31 Mean :28.07 Mean :128.5 Mean :5.292894
## 3rd Qu.:62.00 3rd Qu.:31.40 3rd Qu.:139.0 3rd Qu.:7.335437
## Max. :80.00 Max. :50.20 Max. :192.0 Max. :9.999102
8. Describe the Cohort
and Compare Baseline Characteristics
Before modeling the outcome, describe the exposed and unexposed
groups.
vars <- c(
"age",
"male",
"bmi",
"systolic_bp",
"diabetes",
"hypertension",
"high_cholesterol",
"physical_inactivity"
)
table1 <- CreateTableOne(
vars = vars,
strata = "smoking",
data = cohort,
test = FALSE
)
print(
table1,
smd = TRUE
)
## Stratified by smoking
## Non-smoker Smoker SMD
## n 4057 1768
## age (mean (SD)) 55.09 (9.31) 55.82 (9.32) 0.078
## male (mean (SD)) 0.46 (0.50) 0.51 (0.50) 0.108
## bmi (mean (SD)) 28.60 (4.88) 26.87 (4.76) 0.358
## systolic_bp (mean (SD)) 128.57 (16.03) 128.38 (15.88) 0.012
## diabetes (mean (SD)) 0.15 (0.36) 0.14 (0.35) 0.040
## hypertension (mean (SD)) 0.36 (0.48) 0.43 (0.50) 0.156
## high_cholesterol (mean (SD)) 0.32 (0.47) 0.33 (0.47) 0.018
## physical_inactivity (mean (SD)) 0.36 (0.48) 0.39 (0.49) 0.066
The standardized mean difference (SMD) helps describe baseline
differences between smokers and non-smokers.
However, baseline imbalance alone does not determine whether a
variable is a confounder. Confounder selection should primarily follow
causal reasoning rather than statistical significance testing.
9. Estimate the Crude
Association
First estimate the unadjusted association between smoking and
incident CVD.
cox_crude <- coxph(
Surv(followup_time, cvd_event) ~ smoking,
data = cohort
)
summary(cox_crude)
## Call:
## coxph(formula = Surv(followup_time, cvd_event) ~ smoking, data = cohort)
##
## n= 5825, number of events= 1035
##
## coef exp(coef) se(coef) z Pr(>|z|)
## smokingSmoker 0.47251 1.60402 0.06357 7.433 1.06e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## smokingSmoker 1.604 0.6234 1.416 1.817
##
## Concordance= 0.553 (se = 0.008 )
## Likelihood ratio test= 53.04 on 1 df, p=3e-13
## Wald test = 55.25 on 1 df, p=1e-13
## Score (logrank) test = 56.28 on 1 df, p=6e-14
exp(
cbind(
HR = coef(cox_crude),
confint(cox_crude)
)
)
## HR 2.5 % 97.5 %
## smokingSmoker 1.604018 1.416115 1.816854
The crude hazard ratio compares the observed CVD hazard in smokers
versus non-smokers without controlling for baseline differences.
10. Estimate the
Adjusted Exposure Effect
Next fit a multivariable Cox model adjusting for prespecified
baseline confounders.
cox_adjusted <- coxph(
Surv(followup_time, cvd_event) ~
smoking +
age +
male +
bmi +
systolic_bp +
diabetes +
hypertension +
high_cholesterol +
physical_inactivity,
data = cohort
)
summary(cox_adjusted)
## Call:
## coxph(formula = Surv(followup_time, cvd_event) ~ smoking + age +
## male + bmi + systolic_bp + diabetes + hypertension + high_cholesterol +
## physical_inactivity, data = cohort)
##
## n= 5825, number of events= 1035
##
## coef exp(coef) se(coef) z Pr(>|z|)
## smokingSmoker 0.444877 1.560299 0.064752 6.871 6.40e-12 ***
## age 0.037288 1.037992 0.003326 11.210 < 2e-16 ***
## male 0.354852 1.425969 0.062696 5.660 1.52e-08 ***
## bmi 0.012722 1.012803 0.006445 1.974 0.048389 *
## systolic_bp 0.012370 1.012447 0.001934 6.397 1.58e-10 ***
## diabetes 0.662160 1.938975 0.073602 8.996 < 2e-16 ***
## hypertension 0.402836 1.496061 0.062632 6.432 1.26e-10 ***
## high_cholesterol 0.240616 1.272032 0.064203 3.748 0.000178 ***
## physical_inactivity 0.169128 1.184271 0.063382 2.668 0.007622 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## smokingSmoker 1.560 0.6409 1.374 1.771
## age 1.038 0.9634 1.031 1.045
## male 1.426 0.7013 1.261 1.612
## bmi 1.013 0.9874 1.000 1.026
## systolic_bp 1.012 0.9877 1.009 1.016
## diabetes 1.939 0.5157 1.679 2.240
## hypertension 1.496 0.6684 1.323 1.691
## high_cholesterol 1.272 0.7861 1.122 1.443
## physical_inactivity 1.184 0.8444 1.046 1.341
##
## Concordance= 0.66 (se = 0.009 )
## Likelihood ratio test= 366.8 on 9 df, p=<2e-16
## Wald test = 371.2 on 9 df, p=<2e-16
## Score (logrank) test = 379.7 on 9 df, p=<2e-16
exp(
cbind(
HR = coef(cox_adjusted),
confint(cox_adjusted)
)
)
## HR 2.5 % 97.5 %
## smokingSmoker 1.560299 1.374330 1.771432
## age 1.037992 1.031247 1.044781
## male 1.425969 1.261081 1.612417
## bmi 1.012803 1.000090 1.025678
## systolic_bp 1.012447 1.008617 1.016291
## diabetes 1.938975 1.678503 2.239867
## hypertension 1.496061 1.323234 1.691460
## high_cholesterol 1.272032 1.121626 1.442607
## physical_inactivity 1.184271 1.045924 1.340919
The coefficient for smokingSmoker is the
covariate-adjusted hazard ratio comparing smokers with non-smokers.
Under the required causal assumptions, this adjusted estimate may be
interpreted as an estimate of the exposure effect. Without those
assumptions, it should be interpreted more cautiously as an adjusted
association.
11. Compare Crude and
Adjusted Estimates
A useful epidemiologic step is to compare the crude and adjusted
exposure estimates.
extract_hr <- function(model, term, label) {
est <- coef(model)[term]
ci <- confint(model)[term, ]
data.frame(
Model = label,
HR = exp(est),
Lower95 = exp(ci[1]),
Upper95 = exp(ci[2])
)
}
comparison <- bind_rows(
extract_hr(cox_crude, "smokingSmoker", "Crude"),
extract_hr(cox_adjusted, "smokingSmoker", "Adjusted")
)
comparison
## Model HR Lower95 Upper95
## smokingSmoker...1 Crude 1.604018 1.416115 1.816854
## smokingSmoker...2 Adjusted 1.560299 1.374330 1.771432
A difference between crude and adjusted estimates may indicate the
impact of covariate adjustment, but it does not by
itself prove which variables are true confounders.
13. Sensitivity
Analysis 1: Propensity Score Weighting
Propensity-score methods can also be used in epidemiologic exposure
studies when appropriate.
Here the propensity score is the probability of being a smoker
conditional on measured baseline covariates.
ps_model <- glm(
I(smoking == "Smoker") ~
age +
male +
bmi +
systolic_bp +
diabetes +
hypertension +
high_cholesterol +
physical_inactivity,
family = binomial(),
data = cohort
)
cohort$ps <- predict(
ps_model,
type = "response"
)
cohort <- cohort %>%
mutate(
exposed = as.integer(smoking == "Smoker"),
iptw = ifelse(
exposed == 1,
1 / ps,
1 / (1 - ps)
)
)
summary(cohort$ps)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0547 0.2393 0.2966 0.3035 0.3604 0.6534
summary(cohort$iptw)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.058 1.343 1.514 2.000 2.422 9.494
14. Check Covariate
Balance After Weighting
bal.tab(
smoking ~
age +
male +
bmi +
systolic_bp +
diabetes +
hypertension +
high_cholesterol +
physical_inactivity,
data = cohort,
weights = cohort$iptw,
method = "weighting",
un = TRUE,
thresholds = c(m = 0.10)
)
## Note: `s.d.denom` not specified; assuming "pooled".
## Balance Measures
## Type Diff.Un Diff.Adj M.Threshold
## age Contin. 0.0779 0.0016 Balanced, <0.1
## male Binary 0.0540 -0.0026 Balanced, <0.1
## bmi Contin. -0.3584 -0.0065 Balanced, <0.1
## systolic_bp Contin. -0.0120 0.0117 Balanced, <0.1
## diabetes Binary -0.0141 -0.0011 Balanced, <0.1
## hypertension Binary 0.0761 0.0003 Balanced, <0.1
## high_cholesterol Binary 0.0086 -0.0011 Balanced, <0.1
## physical_inactivity Binary 0.0320 -0.0030 Balanced, <0.1
##
## Balance tally for mean differences
## count
## Balanced, <0.1 8
## Not Balanced, >0.1 0
##
## Variable with the greatest mean difference
## Variable Diff.Adj M.Threshold
## systolic_bp 0.0117 Balanced, <0.1
##
## Effective sample sizes
## Non-smoker Smoker
## Unadjusted 4057. 1768.
## Adjusted 3988.07 1616.19
An absolute SMD below 0.10 is commonly used as a practical indicator
of acceptable balance for a measured covariate.
Good measured balance does not establish that unmeasured confounding
is absent.
15. Sensitivity
Analysis 1: IPTW Outcome Model
# Truncate extreme weights at the 1st and 99th percentiles
limits <- quantile(
cohort$iptw,
probs = c(0.01, 0.99)
)
cohort <- cohort %>%
mutate(
iptw_trim = pmin(
pmax(iptw, limits[1]),
limits[2]
)
)
cox_iptw <- coxph(
Surv(followup_time, cvd_event) ~ smoking,
data = cohort,
weights = iptw_trim,
robust = TRUE
)
summary(cox_iptw)
## Call:
## coxph(formula = Surv(followup_time, cvd_event) ~ smoking, data = cohort,
## weights = iptw_trim, robust = TRUE)
##
## n= 5825, number of events= 1035
##
## coef exp(coef) se(coef) robust se z Pr(>|z|)
## smokingSmoker 0.41449 1.51360 0.04352 0.06550 6.328 2.48e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## smokingSmoker 1.514 0.6607 1.331 1.721
##
## Concordance= 0.549 (se = 0.008 )
## Likelihood ratio test= 92.29 on 1 df, p=<2e-16
## Wald test = 40.05 on 1 df, p=2e-10
## Score (logrank) test = 92.02 on 1 df, p=<2e-16, Robust = 34.93 p=3e-09
##
## (Note: the likelihood ratio and score tests assume independence of
## observations within a cluster, the Wald and robust score tests do not).
exp(
cbind(
HR = coef(cox_iptw),
confint(cox_iptw)
)
)
## HR 2.5 % 97.5 %
## smokingSmoker 1.513602 1.331253 1.720929
16. Sensitivity
Analysis 2: Alternative Covariate Specification
Suppose hypertension has uncertain temporal ordering relative to
smoking history or disease development. We can examine how strongly the
smoking estimate depends on adjustment for hypertension.
cox_without_htn <- coxph(
Surv(followup_time, cvd_event) ~
smoking +
age +
male +
bmi +
systolic_bp +
diabetes +
high_cholesterol +
physical_inactivity,
data = cohort
)
bind_rows(
extract_hr(cox_adjusted, "smokingSmoker", "Adjusted including hypertension"),
extract_hr(cox_without_htn, "smokingSmoker", "Adjusted excluding hypertension")
)
## Model HR Lower95 Upper95
## smokingSmoker...1 Adjusted including hypertension 1.560299 1.374330 1.771432
## smokingSmoker...2 Adjusted excluding hypertension 1.599692 1.409354 1.815736
If the results are similar, the exposure estimate is relatively
robust to this particular modeling choice. If the results differ
substantially, the conclusion is sensitive to assumptions about
hypertension.
This sensitivity analysis does not determine the true causal role of
hypertension; it only shows how much the estimate depends on the
assumption.
17. Compare Main and
Sensitivity Analyses
results <- bind_rows(
extract_hr(cox_crude, "smokingSmoker", "Crude Cox"),
extract_hr(cox_adjusted, "smokingSmoker", "Multivariable-adjusted Cox"),
extract_hr(cox_iptw, "smokingSmoker", "IPTW Cox"),
extract_hr(cox_without_htn, "smokingSmoker", "Adjusted Cox without hypertension")
)
results
## Model HR Lower95 Upper95
## smokingSmoker...1 Crude Cox 1.604018 1.416115 1.816854
## smokingSmoker...2 Multivariable-adjusted Cox 1.560299 1.374330 1.771432
## smokingSmoker...3 IPTW Cox 1.513602 1.331253 1.720929
## smokingSmoker...4 Adjusted Cox without hypertension 1.599692 1.409354 1.815736
18. Interpretation and
Limitations
A typical epidemiologic effect-estimation workflow is:
- Define the epidemiologic / causal research question.
- Choose an appropriate study design and data source.
- Define the target population and analytic cohort.
- Define exposure, outcome, time zero, and follow-up.
- Identify potential confounders using temporal and causal
reasoning.
- Assess data quality and missingness.
- Describe exposed and unexposed groups.
- Estimate the crude exposure-outcome association.
- Estimate the confounder-adjusted effect / association.
- Check model assumptions.
- Evaluate covariate balance when propensity-score methods are
used.
- Perform sensitivity analyses under alternative reasonable
assumptions.
- Interpret the estimate in light of residual confounding, measurement
error, selection bias, and temporality uncertainty.
The key epidemiologic idea is:
The objective is not simply to make two groups statistically similar.
The objective is to estimate the exposure-outcome effect as validly as
possible by using an appropriate design, defining time correctly,
selecting covariates based on causal reasoning, and evaluating the
robustness of the result.
19. Main Limitations
of Observational Epidemiologic Effect Studies
Even after careful adjustment, important limitations may remain:
- unmeasured confounding,
- residual confounding,
- exposure misclassification,
- outcome misclassification,
- measurement error,
- uncertain temporal ordering,
- selection bias,
- informative censoring,
- missing data,
- positivity / overlap problems,
- model misspecification,
- and limited generalizability.
Therefore, an adjusted observational estimate should not
automatically be interpreted as proof of causality.
