library(dplyr)
library(gt)
library(gtsummary)
library(readr)
library(skimr)
library(tibble)Target Trial Basics: Early Vasopressors in Septic Shock
Goal
This notebook introduces the target trial we want to emulate.
The clinical question is:
Among ICU patients with suspected septic shock at ICU admission, what is the effect of starting vasopressors within 2 hours compared with not starting vasopressors within 2 hours on 28-day mortality?
Setup
These packages make the routine analysis code easier to read.
dplyr handles data manipulation, readr reads rectangular files, skimr gives quick data summaries, gt builds display tables, and gtsummary builds analytic summary tables.
Why A Target Trial?
Observational ICU data are not randomized.
Sicker patients are often treated earlier, so a simple comparison between patients who received early vasopressors and patients who did not can be biased.
The target trial framework asks us to describe the randomized trial we wish we had run, then emulate it as closely as possible using observational data.
Target Trial Protocol
target_trial_protocol <- tribble(
~component, ~definition,
"Eligibility criteria", "ICU admission, suspected sepsis, hypotension, elevated lactate",
"Time zero", "ICU admission",
"Treatment strategy 1", "Start vasopressors within 2 hours",
"Treatment strategy 2", "Do not start vasopressors within 2 hours",
"Outcome", "Death within 28 days",
"Causal contrast", "Risk difference and risk ratio"
)
target_trial_protocol |>
gt() |>
tab_header(title = "Target Trial Protocol") |>
cols_label(
component = "Component",
definition = "Definition"
)| Target Trial Protocol | |
| Component | Definition |
|---|---|
| Eligibility criteria | ICU admission, suspected sepsis, hypotension, elevated lactate |
| Time zero | ICU admission |
| Treatment strategy 1 | Start vasopressors within 2 hours |
| Treatment strategy 2 | Do not start vasopressors within 2 hours |
| Outcome | Death within 28 days |
| Causal contrast | Risk difference and risk ratio |
tribble() creates a small tibble by typing the rows directly.
gt() turns that tibble into a clearer presentation table.
Load The Simulated Cohort
This first version uses a CSV generated by scripts/01_simulate_icu_data_base_r.R.
icu_data <- read_csv("../data/icu_septic_shock_base_r.csv", show_col_types = FALSE)read_csv() reads a rectangular CSV file into R.
The object icu_data is a tibble, where each row is one ICU patient.
Inspect The Data
icu_data |>
slice_head(n = 6)# A tibble: 6 × 13
patient_id age sex sofa_score lactate map suspected_sepsis hypotension_at_baseline elevated_lactate_at_baseline
<dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 58 female 4 1.9 65 1 0 0
2 2 69 female 7 1.3 59 1 1 0
3 3 70 male 7 6.4 61 1 1 1
4 4 50 female 10 3.7 55 1 1 1
5 5 41 female 6 3.1 45 1 1 1
6 6 58 male 3 2.3 49 1 1 1
# ℹ 4 more variables: eligible <dbl>, early_vasopressor <dbl>, time_to_vasopressor_hours <dbl>, death_28d <dbl>slice_head() prints the first few rows so we can inspect the structure before analyzing anything.
skim(icu_data)| Name | icu_data |
| Number of rows | 1000 |
| Number of columns | 13 |
| _______________________ | |
| Column type frequency: | |
| character | 1 |
| numeric | 12 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| sex | 0 | 1 | 4 | 6 | 0 | 2 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| patient_id | 0 | 1 | 500.50 | 288.82 | 1.0 | 250.75 | 500.50 | 750.25 | 1000.00 | ▇▇▇▇▇ |
| age | 0 | 1 | 64.98 | 13.92 | 18.0 | 56.00 | 65.00 | 75.00 | 95.00 | ▁▂▇▇▃ |
| sofa_score | 0 | 1 | 7.01 | 2.57 | 1.0 | 5.00 | 7.00 | 9.00 | 17.00 | ▃▇▆▂▁ |
| lactate | 0 | 1 | 3.50 | 1.92 | 0.5 | 2.10 | 3.05 | 4.30 | 13.10 | ▇▆▂▁▁ |
| map | 0 | 1 | 62.16 | 10.18 | 35.0 | 55.00 | 62.00 | 70.00 | 95.00 | ▂▆▇▃▁ |
| suspected_sepsis | 0 | 1 | 0.88 | 0.33 | 0.0 | 1.00 | 1.00 | 1.00 | 1.00 | ▁▁▁▁▇ |
| hypotension_at_baseline | 0 | 1 | 0.59 | 0.49 | 0.0 | 0.00 | 1.00 | 1.00 | 1.00 | ▆▁▁▁▇ |
| elevated_lactate_at_baseline | 0 | 1 | 0.81 | 0.40 | 0.0 | 1.00 | 1.00 | 1.00 | 1.00 | ▂▁▁▁▇ |
| eligible | 0 | 1 | 0.42 | 0.49 | 0.0 | 0.00 | 0.00 | 1.00 | 1.00 | ▇▁▁▁▆ |
| early_vasopressor | 0 | 1 | 0.46 | 0.50 | 0.0 | 0.00 | 0.00 | 1.00 | 1.00 | ▇▁▁▁▇ |
| time_to_vasopressor_hours | 0 | 1 | 7.48 | 7.48 | 0.0 | 1.17 | 3.83 | 13.77 | 23.99 | ▇▂▂▂▂ |
| death_28d | 0 | 1 | 0.24 | 0.43 | 0.0 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▂ |
skim() gives a quick summary of variable types, missingness, and distributions.
Apply Eligibility Criteria
eligible_data <- icu_data |>
filter(eligible == 1)filter() keeps rows that satisfy a condition.
Here, we keep only patients who satisfy the simulated eligibility criteria.
Count Treatment Groups
eligible_data |>
count(early_vasopressor)# A tibble: 2 × 2
early_vasopressor n
<dbl> <int>
1 0 181
2 1 235count() counts how many eligible patients were observed under each treatment group.
In this first simplified dataset:
1means vasopressors started within 2 hours.0means vasopressors were not started within 2 hours.
Estimate Naive Mortality Risks
naive_risks <- eligible_data |>
group_by(early_vasopressor) |>
summarize(
n_patients = n(),
risk_death_28d = mean(death_28d),
.groups = "drop"
)
naive_risks |>
gt() |>
tab_header(title = "Naive 28-Day Mortality Risk") |>
cols_label(
early_vasopressor = "Early vasopressor",
n_patients = "Patients",
risk_death_28d = "28-day mortality risk"
) |>
fmt_number(columns = risk_death_28d, decimals = 3)| Naive 28-Day Mortality Risk | ||
| Early vasopressor | Patients | 28-day mortality risk |
|---|---|---|
| 0 | 181 | 0.260 |
| 1 | 235 | 0.357 |
Because death_28d is coded as 0 or 1, its mean is the proportion who died.
This is a naive comparison because it does not yet adjust for the fact that treatment decisions depend on patient severity.
Estimate Naive Contrasts
risk_early <- naive_risks |>
filter(early_vasopressor == 1) |>
pull(risk_death_28d)
risk_not_early <- naive_risks |>
filter(early_vasopressor == 0) |>
pull(risk_death_28d)
naive_contrasts <- tibble(
measure = c("Risk difference", "Risk ratio"),
value = c(
risk_early - risk_not_early,
risk_early / risk_not_early
)
)
naive_contrasts |>
gt() |>
tab_header(title = "Naive Treatment Contrast") |>
cols_label(
measure = "Measure",
value = "Value"
) |>
fmt_number(columns = value, decimals = 3)| Naive Treatment Contrast | |
| Measure | Value |
|---|---|
| Risk difference | 0.098 |
| Risk ratio | 1.377 |
The risk difference is an absolute difference in 28-day mortality risk.
The risk ratio is a relative comparison of 28-day mortality risk.
Check Confounding By Severity
eligible_data |>
mutate(
early_vasopressor = factor(
early_vasopressor,
levels = c(0, 1),
labels = c("No early vasopressor", "Early vasopressor")
)
) |>
select(early_vasopressor, age, sex, sofa_score, lactate, map, death_28d) |>
tbl_summary(
by = early_vasopressor,
statistic = list(
all_continuous() ~ "{mean} ({sd})",
all_categorical() ~ "{n} ({p}%)"
),
missing = "no"
) |>
add_overall()| Characteristic | Overall N = 4161 |
No early vasopressor N = 1811 |
Early vasopressor N = 2351 |
|---|---|---|---|
| age | 66 (14) | 63 (14) | 67 (14) |
| sex | |||
| female | 201 (48%) | 91 (50%) | 110 (47%) |
| male | 215 (52%) | 90 (50%) | 125 (53%) |
| sofa_score | 7.00 (2.60) | 6.29 (2.50) | 7.55 (2.55) |
| lactate | 3.98 (1.86) | 3.68 (1.60) | 4.22 (2.01) |
| map | 55.4 (6.3) | 56.4 (6.0) | 54.6 (6.4) |
| death_28d | 131 (31%) | 47 (26%) | 84 (36%) |
| 1 Mean (SD); n (%) | |||
This compares baseline severity between the two treatment groups.
If early vasopressor patients have higher SOFA scores, higher lactate, or lower MAP, then the naive comparison mixes the treatment effect with baseline severity differences.
That problem is one reason we need target trial emulation methods rather than a simple treated-versus-untreated comparison.
Next Step
The next lesson should walk through R/simulate_icu_cohort.R line by line.
After that, we can use scripts/02_naive_analysis_base_r.R to compute the first naive association and then discuss why it is not yet a causal estimate.