1. Background and
Research Question
This educational example mirrors the design of a published real-world
retrospective cohort study among patients with type 2 diabetes who were
initially treated with metformin and subsequently added either insulin
or a sulfonylurea.
Research question:
Among eligible patients with type 2 diabetes receiving metformin, is
adding insulin associated with a different risk of myocardial
infarction, stroke, or all-cause death compared with adding a
sulfonylurea?
Because treatment is not randomized, baseline confounding must be
addressed.
library(dplyr)
library(MatchIt)
library(cobalt)
library(survival)
library(tableone)
set.seed(2026)
2. Create / Import
Real-World Data
In a real study, this step would import EHR, claims, registry,
pharmacy, laboratory, or linked administrative data.
Here we simulate a dataset with the same basic structure.
n <- 5000
rwd <- data.frame(
id = 1:n,
age = round(rnorm(n, 62, 9)),
male = rbinom(n, 1, 0.94),
bmi = round(rnorm(n, 31, 5), 1),
hba1c = round(rnorm(n, 8.1, 1.3), 1),
creatinine = round(pmax(rnorm(n, 1.1, 0.35), 0.4), 2),
systolic_bp = round(rnorm(n, 132, 15)),
smoking = rbinom(n, 1, 0.25),
prior_cvd = rbinom(n, 1, 0.25),
hypertension = rbinom(n, 1, 0.65),
ckd = rbinom(n, 1, 0.18)
)
# Treatment selection depends on baseline patient characteristics.
# This intentionally creates confounding by indication.
lp_treat <-
-3.0 +
0.02 * (rwd$age - 60) +
0.55 * (rwd$hba1c - 8) +
0.35 * rwd$prior_cvd +
0.45 * rwd$ckd +
0.02 * (rwd$bmi - 30)
p_insulin <- plogis(lp_treat)
rwd$treatment <- ifelse(
rbinom(n, 1, p_insulin) == 1,
"Insulin",
"Sulfonylurea"
)
rwd$treatment <- factor(
rwd$treatment,
levels = c("Sulfonylurea", "Insulin")
)
table(rwd$treatment)
##
## Sulfonylurea Insulin
## 4579 421
3. Define the Study
Cohort
A real-world study should prespecify eligibility criteria, index
date, baseline period, treatment strategy, follow-up, and outcome.
For demonstration, we restrict the cohort to adults with plausible
baseline values.
cohort <- rwd %>%
filter(
age >= 40,
age <= 85,
hba1c >= 5,
hba1c <= 14,
creatinine > 0
)
nrow(cohort)
## [1] 4912
table(cohort$treatment)
##
## Sulfonylurea Insulin
## 4492 420
4. Define Exposure,
Outcome, and Baseline Covariates
The exposure is treatment intensification:
Potential baseline confounders include age, sex, BMI, HbA1c,
creatinine, blood pressure, smoking, prior cardiovascular disease,
hypertension, and chronic kidney disease.
We now simulate a time-to-event outcome. In a real study, outcome
events would be identified from diagnosis/procedure codes, death
records, or adjudication.
lp_event <-
-4.3 +
0.22 * (cohort$treatment == "Insulin") +
0.025 * (cohort$age - 60) +
0.30 * cohort$prior_cvd +
0.20 * cohort$ckd +
0.15 * cohort$smoking +
0.08 * (cohort$hba1c - 8)
event_rate <- exp(lp_event)
true_event_time <- rexp(
nrow(cohort),
rate = event_rate
)
censor_time <- runif(
nrow(cohort),
min = 1,
max = 5
)
cohort$followup_time <- pmin(true_event_time, censor_time)
cohort$event <- as.integer(
true_event_time <= censor_time
)
table(cohort$event)
##
## 0 1
## 4654 258
5. Data Quality
Checks
Before analysis, check missingness, implausible values, duplicates,
coding consistency, and follow-up.
# Duplicate patient IDs
sum(duplicated(cohort$id))
## [1] 0
# Missingness
colSums(is.na(cohort))
## id age male bmi hba1c
## 0 0 0 0 0
## creatinine systolic_bp smoking prior_cvd hypertension
## 0 0 0 0 0
## ckd treatment followup_time event
## 0 0 0 0
# Basic distributions
summary(
cohort[, c(
"age",
"bmi",
"hba1c",
"creatinine",
"systolic_bp",
"followup_time"
)]
)
## age bmi hba1c creatinine
## Min. :40.00 Min. :11.80 Min. : 5.000 Min. :0.40
## 1st Qu.:56.00 1st Qu.:27.70 1st Qu.: 7.300 1st Qu.:0.87
## Median :62.00 Median :31.00 Median : 8.100 Median :1.12
## Mean :61.97 Mean :31.04 Mean : 8.152 Mean :1.11
## 3rd Qu.:68.00 3rd Qu.:34.50 3rd Qu.: 9.000 3rd Qu.:1.35
## Max. :85.00 Max. :53.20 Max. :12.700 Max. :2.51
## systolic_bp followup_time
## Min. : 81.0 Min. :0.002528
## 1st Qu.:122.0 1st Qu.:1.898389
## Median :132.0 Median :2.904539
## Mean :132.1 Mean :2.918759
## 3rd Qu.:143.0 3rd Qu.:3.941082
## Max. :186.0 Max. :4.998977
6. Examine Baseline
Imbalance Before Adjustment
Because treatment was not randomized, the treatment groups may differ
substantially at baseline.
vars <- c(
"age",
"male",
"bmi",
"hba1c",
"creatinine",
"systolic_bp",
"smoking",
"prior_cvd",
"hypertension",
"ckd"
)
table1_before <- CreateTableOne(
vars = vars,
strata = "treatment",
data = cohort,
test = FALSE
)
print(
table1_before,
smd = TRUE
)
## Stratified by treatment
## Sulfonylurea Insulin SMD
## n 4492 420
## age (mean (SD)) 61.77 (8.60) 64.10 (8.38) 0.275
## male (mean (SD)) 0.93 (0.25) 0.95 (0.21) 0.083
## bmi (mean (SD)) 31.00 (5.00) 31.54 (5.13) 0.108
## hba1c (mean (SD)) 8.08 (1.25) 8.93 (1.26) 0.679
## creatinine (mean (SD)) 1.11 (0.34) 1.13 (0.36) 0.064
## systolic_bp (mean (SD)) 132.13 (15.12) 131.92 (14.79) 0.014
## smoking (mean (SD)) 0.24 (0.43) 0.27 (0.44) 0.056
## prior_cvd (mean (SD)) 0.25 (0.43) 0.30 (0.46) 0.115
## hypertension (mean (SD)) 0.64 (0.48) 0.63 (0.48) 0.027
## ckd (mean (SD)) 0.17 (0.38) 0.24 (0.43) 0.168
7. Estimate the
Propensity Score
The propensity score estimates each patient’s probability of
receiving insulin conditional on measured baseline covariates.
ps_model <- glm(
I(treatment == "Insulin") ~
age +
male +
bmi +
hba1c +
creatinine +
systolic_bp +
smoking +
prior_cvd +
hypertension +
ckd,
family = binomial(),
data = cohort
)
cohort$ps <- predict(
ps_model,
type = "response"
)
summary(cohort$ps)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.006097 0.040404 0.067552 0.085505 0.110650 0.586436
8. Propensity Score
Matching
Perform 1:1 nearest-neighbor propensity-score matching.
m.out <- matchit(
treatment ~
age +
male +
bmi +
hba1c +
creatinine +
systolic_bp +
smoking +
prior_cvd +
hypertension +
ckd,
data = cohort,
method = "nearest",
distance = "glm",
ratio = 1,
caliper = 0.2
)
summary(m.out)
##
## Call:
## matchit(formula = treatment ~ age + male + bmi + hba1c + creatinine +
## systolic_bp + smoking + prior_cvd + hypertension + ckd, data = cohort,
## method = "nearest", distance = "glm", caliper = 0.2, ratio = 1)
##
## Summary of Balance for All Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance 0.1334 0.0810 0.6047 2.0004 0.2155
## age 64.1024 61.7674 0.2786 0.9489 0.0508
## male 0.9524 0.9330 0.0910 . 0.0194
## bmi 31.5438 30.9956 0.1068 1.0556 0.0193
## hba1c 8.9317 8.0789 0.6771 1.0094 0.1108
## creatinine 1.1312 1.1085 0.0630 1.0926 0.0146
## systolic_bp 131.9167 132.1273 -0.0142 0.9563 0.0069
## smoking 0.2690 0.2444 0.0555 . 0.0246
## prior_cvd 0.2976 0.2464 0.1119 . 0.0512
## hypertension 0.6286 0.6418 -0.0274 . 0.0132
## ckd 0.2381 0.1705 0.1586 . 0.0676
## eCDF Max
## distance 0.3328
## age 0.1217
## male 0.0194
## bmi 0.0571
## hba1c 0.2911
## creatinine 0.0500
## systolic_bp 0.0346
## smoking 0.0246
## prior_cvd 0.0512
## hypertension 0.0132
## ckd 0.0676
##
## Summary of Balance for Matched Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance 0.1327 0.1326 0.0011 1.0057 0.0002
## age 64.0740 63.9165 0.0188 0.9976 0.0126
## male 0.9523 0.9523 0.0000 . 0.0000
## bmi 31.5527 31.9504 -0.0775 1.1270 0.0160
## hba1c 8.9253 8.9103 0.0119 1.0273 0.0093
## creatinine 1.1314 1.1252 0.0174 1.0714 0.0128
## systolic_bp 131.9069 130.5012 0.0951 0.9802 0.0144
## smoking 0.2697 0.2840 -0.0323 . 0.0143
## prior_cvd 0.2959 0.3031 -0.0157 . 0.0072
## hypertension 0.6301 0.6492 -0.0395 . 0.0191
## ckd 0.2387 0.2387 0.0000 . 0.0000
## eCDF Max Std. Pair Dist.
## distance 0.0072 0.0027
## age 0.0430 1.0433
## male 0.0000 0.0811
## bmi 0.0597 1.0654
## hba1c 0.0358 0.5695
## creatinine 0.0477 1.0836
## systolic_bp 0.0668 1.1619
## smoking 0.0143 0.8719
## prior_cvd 0.0072 0.8300
## hypertension 0.0191 0.9978
## ckd 0.0000 0.3532
##
## Sample Sizes:
## Control Treated
## All 4492 420
## Matched 419 419
## Unmatched 4073 1
## Discarded 0 0
matched <- match.data(m.out)
table(matched$treatment)
##
## Sulfonylurea Insulin
## 419 419
9. Check Covariate
Balance
A common rule of thumb is an absolute standardized mean difference
(SMD) below 0.10.
bal.tab(
m.out,
un = TRUE,
thresholds = c(m = 0.10)
)
## Balance Measures
## Type Diff.Un Diff.Adj M.Threshold
## distance Distance 0.6047 0.0011 Balanced, <0.1
## age Contin. 0.2786 0.0188 Balanced, <0.1
## male Binary 0.0194 0.0000 Balanced, <0.1
## bmi Contin. 0.1068 -0.0775 Balanced, <0.1
## hba1c Contin. 0.6771 0.0119 Balanced, <0.1
## creatinine Contin. 0.0630 0.0174 Balanced, <0.1
## systolic_bp Contin. -0.0142 0.0951 Balanced, <0.1
## smoking Binary 0.0246 -0.0143 Balanced, <0.1
## prior_cvd Binary 0.0512 -0.0072 Balanced, <0.1
## hypertension Binary -0.0132 -0.0191 Balanced, <0.1
## ckd Binary 0.0676 0.0000 Balanced, <0.1
##
## Balance tally for mean differences
## count
## Balanced, <0.1 11
## Not Balanced, >0.1 0
##
## Variable with the greatest mean difference
## Variable Diff.Adj M.Threshold
## systolic_bp 0.0951 Balanced, <0.1
##
## Sample sizes
## Control Treated
## All 4492 420
## Matched 419 419
## Unmatched 4073 1
love.plot(
m.out,
threshold = 0.10,
abs = TRUE
)
## Warning: Standardized mean differences and raw mean differences are present in
## the same plot. Use the `stars` argument to distinguish between them and
## appropriately label the x-axis. See `?love.plot` for details.
Good balance among measured covariates does not
prove that unmeasured confounding has been eliminated.
10. Primary Outcome
Analysis
Use a Cox proportional hazards model to compare time to the composite
event between treatment groups.
cox_fit <- coxph(
Surv(followup_time, event) ~ treatment,
data = matched,
weights = weights,
robust = TRUE
)
summary(cox_fit)
## Call:
## coxph(formula = Surv(followup_time, event) ~ treatment, data = matched,
## weights = weights, robust = TRUE)
##
## n= 838, number of events= 57
##
## coef exp(coef) se(coef) robust se z Pr(>|z|)
## treatmentInsulin -0.1788 0.8363 0.2660 0.2645 -0.676 0.499
##
## exp(coef) exp(-coef) lower .95 upper .95
## treatmentInsulin 0.8363 1.196 0.498 1.404
##
## Concordance= 0.519 (se = 0.035 )
## Likelihood ratio test= 0.45 on 1 df, p=0.5
## Wald test = 0.46 on 1 df, p=0.5
## Score (logrank) test = 0.45 on 1 df, p=0.5, Robust = 0.46 p=0.5
##
## (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_fit),
confint(cox_fit)
)
)
## HR 2.5 % 97.5 %
## treatmentInsulin 0.8362955 0.4979732 1.404473
Interpretation:
- HR > 1: higher event hazard in the insulin group.
- HR < 1: lower event hazard in the insulin group.
- A 95% confidence interval containing 1 indicates that the difference
is not statistically significant at the conventional 0.05 level.
Because this is an observational study, the estimate should be
interpreted as an adjusted association / causal estimate under the
required causal assumptions rather than proof of causality.
11. Sensitivity
Analysis 1: IPTW
Instead of matching, estimate the average treatment effect using
inverse probability of treatment weighting.
cohort <- cohort %>%
mutate(
treat_num = as.integer(treatment == "Insulin"),
weight_iptw = ifelse(
treat_num == 1,
1 / ps,
1 / (1 - ps)
)
)
# Truncate extreme weights at the 1st and 99th percentiles
limits <- quantile(
cohort$weight_iptw,
probs = c(0.01, 0.99)
)
cohort <- cohort %>%
mutate(
weight_trim = pmin(
pmax(weight_iptw, limits[1]),
limits[2]
)
)
summary(cohort$weight_trim)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.013 1.043 1.076 1.867 1.140 19.961
cox_iptw <- coxph(
Surv(followup_time, event) ~ treatment,
data = cohort,
weights = weight_trim,
robust = TRUE
)
summary(cox_iptw)
## Call:
## coxph(formula = Surv(followup_time, event) ~ treatment, data = cohort,
## weights = weight_trim, robust = TRUE)
##
## n= 4912, number of events= 258
##
## coef exp(coef) se(coef) robust se z Pr(>|z|)
## treatmentInsulin -0.03438 0.96621 0.09311 0.23578 -0.146 0.884
##
## exp(coef) exp(-coef) lower .95 upper .95
## treatmentInsulin 0.9662 1.035 0.6087 1.534
##
## Concordance= 0.52 (se = 0.03 )
## Likelihood ratio test= 0.14 on 1 df, p=0.7
## Wald test = 0.02 on 1 df, p=0.9
## Score (logrank) test = 0.14 on 1 df, p=0.7, Robust = 0.02 p=0.9
##
## (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 %
## treatmentInsulin 0.9662085 0.608663 1.533786
12. Sensitivity
Analysis 2: Conventional Multivariable Adjustment
cox_adjusted <- coxph(
Surv(followup_time, event) ~
treatment +
age +
male +
bmi +
hba1c +
creatinine +
systolic_bp +
smoking +
prior_cvd +
hypertension +
ckd,
data = cohort
)
summary(cox_adjusted)
## Call:
## coxph(formula = Surv(followup_time, event) ~ treatment + age +
## male + bmi + hba1c + creatinine + systolic_bp + smoking +
## prior_cvd + hypertension + ckd, data = cohort)
##
## n= 4912, number of events= 258
##
## coef exp(coef) se(coef) z Pr(>|z|)
## treatmentInsulin -0.0134323 0.9866575 0.2123178 -0.063 0.949555
## age 0.0250074 1.0253227 0.0072740 3.438 0.000586 ***
## male 0.5347248 1.7069784 0.3228067 1.656 0.097623 .
## bmi 0.0020259 1.0020279 0.0124538 0.163 0.870778
## hba1c 0.1531915 1.1655482 0.0493613 3.103 0.001913 **
## creatinine -0.1953313 0.8225621 0.1813973 -1.077 0.281563
## systolic_bp 0.0000788 1.0000788 0.0041133 0.019 0.984715
## smoking 0.1121143 1.1186407 0.1407215 0.797 0.425619
## prior_cvd 0.3020455 1.3526228 0.1348161 2.240 0.025063 *
## hypertension 0.1516527 1.1637560 0.1328983 1.141 0.253821
## ckd -0.0242161 0.9760748 0.1643188 -0.147 0.882838
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## treatmentInsulin 0.9867 1.0135 0.6508 1.496
## age 1.0253 0.9753 1.0108 1.040
## male 1.7070 0.5858 0.9067 3.214
## bmi 1.0020 0.9980 0.9779 1.027
## hba1c 1.1655 0.8580 1.0581 1.284
## creatinine 0.8226 1.2157 0.5765 1.174
## systolic_bp 1.0001 0.9999 0.9920 1.008
## smoking 1.1186 0.8939 0.8490 1.474
## prior_cvd 1.3526 0.7393 1.0385 1.762
## hypertension 1.1638 0.8593 0.8969 1.510
## ckd 0.9761 1.0245 0.7073 1.347
##
## Concordance= 0.587 (se = 0.018 )
## Likelihood ratio test= 33.14 on 11 df, p=5e-04
## Wald test = 32.78 on 11 df, p=6e-04
## Score (logrank) test = 32.91 on 11 df, p=5e-04
exp(
cbind(
HR = coef(cox_adjusted),
confint(cox_adjusted)
)
)
## HR 2.5 % 97.5 %
## treatmentInsulin 0.9866575 0.6507905 1.495862
## age 1.0253227 1.0108087 1.040045
## male 1.7069784 0.9066808 3.213673
## bmi 1.0020279 0.9778655 1.026787
## hba1c 1.1655482 1.0580685 1.283946
## creatinine 0.8225621 0.5764518 1.173747
## systolic_bp 1.0000788 0.9920486 1.008174
## smoking 1.1186407 0.8490011 1.473917
## prior_cvd 1.3526228 1.0385347 1.761702
## hypertension 1.1637560 0.8968890 1.510029
## ckd 0.9760748 0.7073179 1.346950
13. Compare Results
Across Methods
extract_hr <- function(model, label) {
est <- coef(model)["treatmentInsulin"]
ci <- confint(model)["treatmentInsulin", ]
data.frame(
Method = label,
HR = exp(est),
Lower95 = exp(ci[1]),
Upper95 = exp(ci[2])
)
}
results <- bind_rows(
extract_hr(cox_fit, "Propensity-score matching"),
extract_hr(cox_iptw, "IPTW"),
extract_hr(cox_adjusted, "Multivariable Cox")
)
results
## Method HR Lower95 Upper95
## treatmentInsulin...1 Propensity-score matching 0.8362955 0.4979732 1.404473
## treatmentInsulin...2 IPTW 0.9662085 0.6086630 1.533786
## treatmentInsulin...3 Multivariable Cox 0.9866575 0.6507905 1.495862
If the estimates are reasonably similar across multiple prespecified
approaches, confidence in robustness is increased. However, agreement
does not eliminate the possibility of unmeasured confounding or other
biases.
14. Interpretation and
Limitations
The real-world evidence workflow is:
- Define the causal question.
- Identify the real-world data source.
- Define eligibility criteria and the index date.
- Define treatment strategies.
- Define outcomes and follow-up.
- Prespecify baseline confounders using clinical knowledge and causal
reasoning.
- Assess data quality.
- Evaluate baseline imbalance.
- Adjust measured confounding using matching, weighting, or
regression.
- Check covariate balance and positivity.
- Estimate the treatment effect.
- Perform sensitivity analyses.
- Interpret the findings while acknowledging residual and unmeasured
confounding.