Getting parameters that define survival curve

We identified a comparable study that reported quintiles for the time-to-event outcome. Specifically, 20% of participants had a follow-up within 1.4 months, 40% by 4.7 months, 60% by 8.7 months, and 80% by just over 15 months. We used the survGetParams function from the simstudy package to estimate the Weibull distribution parameters—the intercept in the Weibull formula and the shape—that characterize this baseline survival curve.

q20 <- c(1.44, 4.68, 8.69, 15.32)

points <- list(c(q20[1], 0.80), c(q20[2], 0.60), c(q20[3], 0.40), c(q20[4], 0.20))
s <- survGetParams(points)

s
## [1] -1.868399  1.194869

We can visualize the idealized survival curve that will be generated using these parameters stored in the vector s:

survParamPlot(f = s[1], shape = s[2], points, limits = c(0, 20))

Visualizing the curve and assessing its properties

A plot of the survival curves from the two arms is shown below, with the control arm in yellow and the intervention arm in red:

Fitting a model

I fit a Cox proportional hazards model just to make sure I could recover the hazard ratio I used in generating the data:

fit <- coxph(Surv(time, event) ~  factor(A), data = dd)
tidy(fit, exponentiate = TRUE)
## # A tibble: 1 × 5
##   term       estimate std.error statistic p.value
##   <chr>         <dbl>     <dbl>     <dbl>   <dbl>
## 1 factor(A)1     1.99   0.00665      103.       0

Relationship of HR to median time-to-event

My collaborator is especially interested in how the interventions might shift the median time-to-event. This essentially raises the question of how the hazard ratio translates to a change in the median. To explore this, I generated a series of data sets using hazard ratios ranging from 1 to 2 and recorded the observed median for each. As expected, when the hazard ratio is 1, the median closely aligns with that of the baseline distribution.

getmedian <- function(eff = 0) {

  dd <- genData(100000, def)
  dd <- genSurv(dd, defS, timeName = "time", censorName = "censor")
  dd[A == 1, round(median(time), 1)]

}

dm <- rbindlist(lapply(
  log(seq(1, 2, by = .1)),
  function(x) data.table(HR = exp(x), median = getmedian(x))
))

Data definititions

While the model is semi-parametric (i.e., it does not assume a specific distribution for event times), the data generation process is fully parametric, based on the Weibull distribution. Despite this difference, the two are closely aligned: if all goes well, we should be able to recover the parameters used for data generation when fitting the semi-parametric model as we did in the simpler RCT case above.

Data generation occurs in three broad steps:

1. Clinic-level: generate the clinic-specific random effect (g).
2. Provider-level: generate the provider-specific random effect (b) and assign treatment.
3. Patient-level: generate individual time-to-event outcomes.

These steps are implemented using the following definitions:

defC <- defData(varname = "g", formula = 0, variance = "..s2_clinic")

defP <-
  defDataAdd(varname = "b", formula = 0, variance = "..s2_prov") |>
  defDataAdd(varname = "A", formula = "1;1;1", variance = "clinic", dist = "trtAssign")

defS <- defSurv(
  varname = "eventTime",
  formula = "..int + b + g + (A==2)*..eff_A + (A==3)*..eff_AB",
  shape = "..shape") 

Data generation

For this simulation, we assumed 16 clinics, each with 6 providers, and 48 patients per provider. A key element of the study is that recruitment occurs over 12 months, and patients are followed for up to 6 months after their recruitment period ends. Thus, follow-up duration varies depending on when a patient enters the study: patients recruited earlier have longer potential follow-up, while those recruited later are more likely to be censored.

This staggered follow-up is implemented in the final step of data generation:

nC <- 16               # number of clinics (clusters)
nP <- 6                # number of providers per clinic
nI <- 48               # number of patients per provider
s2_clinic <- 0.10      # variation across clinics (g)
s2_prov <- 0.25        # variation across providers (b)
eff_A <- log(c(1.4))   # log HR of intervention A (compared to control)
eff_AB <- log(c(1.6))  # log HR of combined A+B (compared to control)

ds <- genData(nC, defC, id = "clinic")
dp <- genCluster(ds, "clinic", nP, "provider")
dp <- addColumns(defP, dp)
dd <- genCluster(dp, "provider", nI, "id")
dd <- genSurv(dd, defS)

# assign a patient to a particular month - 4 per month

dd <- trtAssign(dd, nTrt = 12, strata = "provider", grpName = "month")

dd[, event := as.integer(eventTime <= 18 - month)]
dd[, obsTime := pmin(eventTime, 18 - month)]

Below is a Kaplan-Meier plot showing survival curves for each provider within each clinic, color-coded by study arm:

The mixed-effects Cox model recovers the variance components and coefficients used in data generation:

me.fit <- coxme(
  Surv(eventTime, event) ~ factor(A) + (1|provider) + (1|clinic),
  data = dd
)

summary(me.fit)
## Mixed effects coxme model
##  Formula: Surv(eventTime, event) ~ factor(A) + (1 | provider) + (1 | clinic)
##     Data: dd
##
##   events, n = 3411, 4608
##
## Random effects:
##      group  variable        sd   variance
## 1 provider Intercept 0.5437050 0.29561511
## 2   clinic Intercept 0.3051615 0.09312354
##                    Chisq    df p    AIC   BIC
## Integrated loglik  919.4  4.00 0  911.4 886.9
##  Penalized loglik 1251.0 86.81 0 1077.3 544.8
##
## Fixed effects:
##              coef exp(coef) se(coef)    z       p
## factor(A)2 0.3207    1.3781   0.1493 2.15 0.03173
## factor(A)3 0.4650    1.5921   0.1472 3.16 0.00158

Typically, I would use this data generation and model fitting code to estimate power or sample size requirements. While I did carry out those steps, I’ve left them out here so that you can try them yourself (though I’m happy to share my code if you’re interested). Beyond estimating sample size, simulation studies like this can also be used to evaluate the Type I error rates of the gatekeeping hypothesis testing framework.

References: