266 lines
8.7 KiB
Plaintext
266 lines
8.7 KiB
Plaintext
---
|
||
title: "Censoring and Inverse Probability of Censoring Weighting"
|
||
format:
|
||
html:
|
||
embed-resources: true
|
||
docx: default
|
||
execute:
|
||
echo: true
|
||
warning: false
|
||
message: false
|
||
---
|
||
|
||
## Goal
|
||
|
||
This notebook adds **inverse probability of censoring weighting (IPCW)** to the cloned dataset from notebook 05.
|
||
|
||
The key idea is that clones are censored for a reason: sicker patients might die before the 2-hour window closes, or clinicians might deviate from the protocol based on severity.
|
||
|
||
IPCW weights the uncensored clones so they represent the full population that started in each trial arm.
|
||
|
||
## Setup
|
||
|
||
```{r}
|
||
suppressPackageStartupMessages({
|
||
library(dplyr)
|
||
library(gt)
|
||
library(tibble)
|
||
library(tidyr)
|
||
})
|
||
|
||
source("../R/simulate_icu_cohort.R")
|
||
source("../R/simulate_icu_cohort_longitudinal.R")
|
||
source("../R/clone_trial_arms.R")
|
||
source("../R/estimate_iptw_and_ipcw_effect.R")
|
||
source("../R/estimate_iptw_vasopressor_mortality_effect.R")
|
||
```
|
||
|
||
## Simulate, Clone, and Weight
|
||
|
||
```{r}
|
||
longitudinal_data <- simulate_icu_cohort_longitudinal(n_patients = 1000, seed = 20260531)
|
||
clones <- clone_trial_arms(longitudinal_data)
|
||
weighted_effects <- estimate_iptw_and_ipcw_effect(clones$clone_baseline)
|
||
```
|
||
|
||
## Propensity Score Model (IPTW)
|
||
|
||
The IPTW model predicts the probability of receiving early vasopressors given baseline confounders.
|
||
|
||
It is identical in purpose to the model in notebook 04, but it is now applied to the **cloned** dataset.
|
||
|
||
Because both clones of the same patient share the same observed treatment, we fit the model once per patient, then assign the propensity score to both clones.
|
||
|
||
```{r}
|
||
propensity_coefficients <- summary(
|
||
weighted_effects$propensity_score_model
|
||
)$coefficients |>
|
||
as.data.frame() |>
|
||
rownames_to_column("model_term") |>
|
||
as_tibble()
|
||
|
||
names(propensity_coefficients) <- c(
|
||
"model_term",
|
||
"log_odds_estimate",
|
||
"standard_error",
|
||
"z_statistic",
|
||
"p_value"
|
||
)
|
||
|
||
propensity_coefficients |>
|
||
mutate(odds_ratio = exp(log_odds_estimate)) |>
|
||
select(model_term, log_odds_estimate, odds_ratio, standard_error, p_value) |>
|
||
gt() |>
|
||
tab_header(title = "Propensity Score Model (Treatment Assignment)") |>
|
||
cols_label(
|
||
model_term = "Model term",
|
||
log_odds_estimate = "Log-odds estimate",
|
||
odds_ratio = "Odds ratio",
|
||
standard_error = "Standard error",
|
||
p_value = "P-value"
|
||
) |>
|
||
fmt_number(
|
||
columns = c(log_odds_estimate, odds_ratio, standard_error, p_value),
|
||
decimals = 3
|
||
)
|
||
```
|
||
|
||
Higher SOFA score and lactate are associated with a higher probability of early vasopressors.
|
||
|
||
This is the same confounding pattern we adjusted for in notebook 04.
|
||
|
||
## Censoring Model (IPCW)
|
||
|
||
The IPCW model predicts the probability that a clone is censored at the 2-hour mark.
|
||
|
||
Censoring can happen for two reasons:
|
||
|
||
1. **Protocol deviation**: the patient did not follow the assigned clone strategy.
|
||
2. **Death before treatment**: the patient died before the 2-hour window closed.
|
||
|
||
Both reasons may depend on baseline severity, so we model censoring using the same confounders plus the clone strategy.
|
||
|
||
```{r}
|
||
censoring_coefficients <- summary(
|
||
weighted_effects$censoring_model
|
||
)$coefficients |>
|
||
as.data.frame() |>
|
||
rownames_to_column("model_term") |>
|
||
as_tibble()
|
||
|
||
names(censoring_coefficients) <- c(
|
||
"model_term",
|
||
"log_odds_estimate",
|
||
"standard_error",
|
||
"z_statistic",
|
||
"p_value"
|
||
)
|
||
|
||
censoring_coefficients |>
|
||
mutate(odds_ratio = exp(log_odds_estimate)) |>
|
||
select(model_term, log_odds_estimate, odds_ratio, standard_error, p_value) |>
|
||
gt() |>
|
||
tab_header(title = "Censoring Model (Probability of Being Censored at 2 Hours)") |>
|
||
cols_label(
|
||
model_term = "Model term",
|
||
log_odds_estimate = "Log-odds estimate",
|
||
odds_ratio = "Odds ratio",
|
||
standard_error = "Standard error",
|
||
p_value = "P-value"
|
||
) |>
|
||
fmt_number(
|
||
columns = c(log_odds_estimate, odds_ratio, standard_error, p_value),
|
||
decimals = 3
|
||
)
|
||
```
|
||
|
||
A positive coefficient means the factor is associated with a **higher** probability of being censored.
|
||
|
||
If sicker patients (higher SOFA, higher lactate) are more likely to be censored, the IPCW weight will give more weight to the uncensored sicker patients so the analysis represents the full population.
|
||
|
||
## Combined Weight Diagnostics
|
||
|
||
The combined weight for each clone is:
|
||
|
||
`combined_weight = IPTW_weight × IPCW_weight`
|
||
|
||
- IPTW balances baseline confounders across treatment arms.
|
||
- IPCW accounts for informative censoring.
|
||
|
||
Extreme weights indicate either a positivity problem (some patients have near-zero probability of treatment or near-zero probability of remaining uncensored) or model misspecification.
|
||
|
||
```{r}
|
||
weighted_effects$weight_diagnostics |>
|
||
gt() |>
|
||
tab_header(title = "Combined IPTW + IPCW Weight Diagnostics") |>
|
||
cols_label(
|
||
clone_strategy = "Clone strategy",
|
||
n_clones = "Clones",
|
||
min_weight = "Minimum weight",
|
||
max_weight = "Maximum weight",
|
||
mean_weight = "Mean weight",
|
||
median_weight = "Median weight"
|
||
) |>
|
||
fmt_integer(columns = n_clones) |>
|
||
fmt_number(columns = c(min_weight, max_weight, mean_weight, median_weight), decimals = 3)
|
||
```
|
||
|
||
If the maximum weights are very large, we would consider trimming them (capping at the 99th percentile) or stabilizing them.
|
||
|
||
For this teaching example the weights are moderate, so we proceed.
|
||
|
||
## Weighted Mortality Risks
|
||
|
||
Using the combined weights, we compute the weighted 28-day mortality risk in each clone arm.
|
||
|
||
Only uncensored clones are included in the numerator, but their weights account for the censored clones who would have had similar outcomes.
|
||
|
||
```{r}
|
||
weighted_effects$weighted_effect_estimates |>
|
||
head(2) |>
|
||
gt() |>
|
||
tab_header(title = "IPTW + IPCW 28-Day Mortality Risks") |>
|
||
cols_label(
|
||
estimate = "Estimate",
|
||
value = "Value"
|
||
) |>
|
||
fmt_number(columns = value, decimals = 3)
|
||
```
|
||
|
||
These are **per-protocol** risks: the risk under each strategy for patients who were able to follow the strategy through the 2-hour window.
|
||
|
||
## IPTW + IPCW Mortality Effect Estimates
|
||
|
||
```{r}
|
||
weighted_effects$weighted_effect_estimates |>
|
||
tail(2) |>
|
||
gt() |>
|
||
tab_header(title = "IPTW + IPCW Mortality Effect Estimates") |>
|
||
cols_label(
|
||
estimate = "Estimate",
|
||
value = "Value"
|
||
) |>
|
||
fmt_number(columns = value, decimals = 3)
|
||
```
|
||
|
||
## Comparison with Cross-Sectional IPTW
|
||
|
||
How does the longitudinal IPTW + IPCW estimate compare with the cross-sectional IPTW estimate from notebook 04?
|
||
|
||
```{r}
|
||
# Pull the cross-sectional IPTW estimate from notebook 04
|
||
icu_data <- simulate_icu_cohort(n_patients = 1000, seed = 20260531)
|
||
iptw_analysis <- estimate_iptw_vasopressor_mortality_effect(icu_data)
|
||
|
||
cross_sectional_rd <- iptw_analysis$iptw_mortality_effect_estimates |>
|
||
filter(estimate == "IPTW risk difference") |>
|
||
pull(value)
|
||
|
||
cross_sectional_rr <- iptw_analysis$iptw_mortality_effect_estimates |>
|
||
filter(estimate == "IPTW risk ratio") |>
|
||
pull(value)
|
||
|
||
longitudinal_rd <- weighted_effects$weighted_effect_estimates |>
|
||
filter(estimate == "IPTW + IPCW risk difference") |>
|
||
pull(value)
|
||
|
||
longitudinal_rr <- weighted_effects$weighted_effect_estimates |>
|
||
filter(estimate == "IPTW + IPCW risk ratio") |>
|
||
pull(value)
|
||
|
||
comparison <- tibble(
|
||
method = c("Cross-sectional IPTW", "Longitudinal IPTW + IPCW"),
|
||
risk_difference = c(cross_sectional_rd, longitudinal_rd),
|
||
risk_ratio = c(cross_sectional_rr, longitudinal_rr)
|
||
)
|
||
|
||
comparison |>
|
||
gt() |>
|
||
tab_header(title = "Cross-Sectional Versus Longitudinal IPTW Estimates") |>
|
||
cols_label(
|
||
method = "Method",
|
||
risk_difference = "Risk difference",
|
||
risk_ratio = "Risk ratio"
|
||
) |>
|
||
fmt_number(columns = c(risk_difference, risk_ratio), decimals = 3)
|
||
```
|
||
|
||
### What to notice
|
||
|
||
- The **cross-sectional IPTW** estimate treats the treatment as assigned at baseline, with no treatment window.
|
||
- The **longitudinal IPTW + IPCW** estimate respects the grace period, censors protocol deviations, and reweights for informative censoring.
|
||
- The two estimates may differ because the longitudinal design removes patients who could not adhere to the assigned strategy, and the IPCW adjustment accounts for the severity of those who were censored.
|
||
|
||
Both are valid approaches, but they answer slightly different questions:
|
||
|
||
- Cross-sectional IPTW answers: "What is the effect of receiving early vasopressors versus not, among all eligible patients?"
|
||
- Longitudinal IPTW + IPCW answers: "What is the effect of the early vasopressor strategy versus the no early vasopressor strategy, among patients who can adhere to the strategy through the 2-hour window?"
|
||
|
||
The second is closer to the **per-protocol effect** of the target trial.
|
||
|
||
## Next Step
|
||
|
||
The next tutorial step is to put the full pipeline together in one place: simulate, clone, weight, and estimate.
|
||
|
||
We will also compare all four methods we have learned so far: naive, standardized, cross-sectional IPTW, and longitudinal IPTW + IPCW.
|