In elong0527/qte: Quantile Treatment Effects under Censoring and Truncation by Death

library(dplyr)
devtools::load_all()
load("../simulation/simu_qte_v3.Rdata")

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Simulation Setup

We conducted simulation studies to verify the validity of the proposed IPW and AIPW estimators.

Simulation set up

• Total sample size: n = 400
• Analysis time $L = 1$
• Covariates $X$ is generated from $\mathcal{N}(0.5, 1)$
• Treatment group $A$ is generated from $Bernoulli(p)$, where $p=logit^{-1}(-0.25 + 0.5 * x)$ such that the randomization ratio is close to 1:1.
• Outcome $Y$ is generated from $\mathcal{N}(\eta, 1)$, where $\eta = 0.4 + 0.8 * A + 0.3 * X$.
• Yilong: add dependency with $X$
• Death time: $T \sim \exp(\lambda)$, where $\lambda = exp(-2 + 0.5 * X)$ (for a simple Cox PH model)
• Censoring time: $C \sim \exp(-2 + 0.4 * x)$
• Yilong: update censoring distribution depends on $X$
• Observed time: $U = T \wedge C$

• Event indicator: $\delta = T < C$

• Missing data: y_obs in simulated data

• $Y$ is missing if $U < L$
• we summarize the missing data proportion based on the event indicator.
• Actual data: the underline true value of y_actual in simulated data
• y_actual = $\min(y) - 1$ if $T < L$

Example of simulated data

head(db) %>% select(a, x, y, y_obs, y_actual, event_time, censor_time, obs_time, status, analysis_time) %>% 
             knitr::kable(digits = 3)

Method to calculate true value

The true value is estimated by y_actual without using propensity score and censoring adjustment.

The sample size to calculate true value is 10,000 with replication of 1,000 times.

truth_est

Model Estimation

The calculation is within each treatment group

1. Propensity score: Fit a logistic model: glm(a ~ x, family = "binomial", data = db)
2. IPW weight $\pi$: Coxph model: coxph(Surv(obs_time, 1 - status) ~ x, data = db) and calculate the probability at the minimal of obs_time and analysis_time
3. $F(y\mid x; \theta)$
4. Fit a linear model based on observed outcome lm(y_obs ~ x, data = db)
5. Calculate the scaled CDF for f_q

y_std <- (q - db$y_mean) / db$y_sigma
rho + (1 - rho) * pnorm(y_std)

1. $R$: (obs_time > analysis_time) | (obs_time <= analysis_time & status = TRUE)

Simulation Setup

• a: treatment group
• xi the quantile of interest
• n: sample size in each group for true value estimation
• pct_missing: Percent of subject with missing value at analysis time $L$
• pct_missing_death: Percent of subject dead on or before analysis time $L$
• pct_missing_censor: Percent of subject censored on or before analysis time $L$

truth_summary %>% select(-q, -diff) %>% knitr::kable(digits = 1)

Simulation Results

• xi the quantile of interest
• use_propensity: whether to use propensity score is the estimator equation.
• use_km: whether to use KM estimator to adjust censoring
• km_rho: method to estimate the probability of death at $L$
• rho: raw estimator (observed death before $L$) / number of subjects
• rho_km: estimated from Cox model coxph(Surv(obs_time, 1 - status) ~ x, data = db) at time $L$.
• rho_km_simple: estimated from KM estimator at time $L$.
• type: using IPW (formula 5) or AIPW (formula 6).
• x0, x0_true: estimated or true value for the percentile in control group
• x1, x1_true: estimated or true value for the percentile in treatment group
• diff, diff_true: estimated or true value for the percentile difference between two group

Results with observing censoring

res_summary %>% ungroup() %>% 
  subset(censoring) %>% 
  mutate_if(is.numeric, round, 3)  %>% show_db()

Results without observing censoring

res_summary %>% ungroup() %>% 
  subset( (!censoring) & n == 400) %>% 
  mutate_if(is.numeric, round, 3) %>% show_db()

elong0527/qte documentation built on Jan. 19, 2021, 10:24 p.m.