Overview

In tackling this problem, there were four key elements that I needed to work out before actually conducting the power simulations. First, I needed to determine the hypothetical survival curve in the context of a single (control) arm and simulate data to confirm that I could replicate the desired curve. Second, I wanted to generate cluster-level variation so that I could assess the implications of the variance assumptions (still in a single-arm context). Third, I generated two intervention arms (without any clustering) to assess effect size assumptions. And lastly, I generated a full data set that included clustering, randomization, and censoring, and then fit a model that would be the basis for the power analysis to ensure that everything was working as expected. Once this was all completed, I was confident that I could move on to generating estimates of power under a range of sample size and variability assumptions. I apologize in advance for a post that is a bit long, but the agenda is quite packed and includes a lot of code.

Defining shape of survival curve

Defining the shape of the survival curve is made relatively easy using the function survGetParams in the simstudy package. All we need to do is specify some (at least two) coordinates along the curve and the function will return the parameters for the mean and shape of a Weibull function that best fit the points. These parameters are used in the data generation process. In this case, the study’s investigators provided me with a couple of points, indicating that approximately 10% of the sample would have an event by day 30, and half would have an event at day 365. Since the study is following patients at most for 365 days, we will consider anything beyond that to be censored (more on censoring later).

To get things started, here are the libraries needed for all the code that follows:

library(simstudy)
library(data.table)
library(survival)
library(survminer)
library(GGally)
library(coxme)
library(parallel)

Now, we get the parameters that define the survival curve:

points <- list(c(30, 0.90), c(365, .50))
r <- survGetParams(points)
r
## [1] -4.8  1.3

To simulate data from this curve, the time-to-event variable tte is defined using these parameters generating by survGetParams. The observed time is the minimum of one year and the time-to-event, and we create an event indicator if the time-to-even is less than one year.

defs <- defSurv(varname = "tte", formula = r[1], shape = r[2])

defa <- defDataAdd(varname = "time", formula = "min(365, tte)", dist = "nonrandom")
defa <- defDataAdd(defa, "event", "1*(tte <= 365)", dist = "nonrandom")

Generating the data is quite simple in this case:

set.seed(589823)

dd <- genData(1000)
dd <- genSurv(dd, defs, digits = 0)
dd <- addColumns(defa, dd)

dd
##         id  tte time event
##    1:    1  234  234     1
##    2:    2 1190  365     0
##    3:    3  395  365     0
##    4:    4  178  178     1
##    5:    5   72   72     1
##   ---                     
##  996:  996  818  365     0
##  997:  997  446  365     0
##  998:  998  118  118     1
##  999:  999  232  232     1
## 1000: 1000 5308  365     0

The plots below show the source function determined by the parameters on the left and the actual data generated on the right. It appears that the generated data matches the data generation process:

splot <- survParamPlot(r[1], r[2], points = points, n = 1000, limits = c(0, 365) )

fit <- survfit(Surv(time, event) ~ 1, data = dd)

j <- ggsurv(fit, CI = FALSE, surv.col = "#ed7c67", size.est = 0.8) + 
  theme(panel.grid = element_blank(),
        axis.text = element_text(size = 7.5),
        axis.title = element_text(size = 8, face = "bold"),
        plot.title = element_blank()) +
  ylim(0, 1) +
  xlab("time") + ylab("probability of survival")

ggarrange(splot, j, ncol = 2, nrow = 1)

Evaluating cluster variation

Cluster variation in the context of survival curves implies that there is a cluster-specific survival curve. This variation is induced with a random effect in in the data generation process. In this case, I am assuming a normally distributed random effect with mean 0 and some variance (distributions other than a normal distribution can be used). The variance assumption is a key one (which will ultimately impact the estimates of power), and I explore that a bit more in the second part of this section.

defc <- defData(varname = "b", formula = 0, variance = 0.1)

defs <- defSurv(varname = "tte", formula = "r[1] + b", shape = r[2])

defa <- defDataAdd(varname = "time", formula = "min(365, tte)", dist = "nonrandom")
defa <- defDataAdd(defa, "event", "1*(tte <= 365)", dist = "nonrandom")
dc <- genData(20, defc, id = "pcp")
dd <- genCluster(dc, "pcp", numIndsVar = 1000, "id")
dd <- genSurv(dd, defs, digits = 0)
dd.10 <- addColumns(defa, dd)
dd.10
##        pcp    b    id  tte time event
##     1:   1 0.25     1  788  365     0
##     2:   1 0.25     2   93   93     1
##     3:   1 0.25     3  723  365     0
##     4:   1 0.25     4  151  151     1
##     5:   1 0.25     5 1367  365     0
##    ---                               
## 19996:  20 0.22 19996  424  365     0
## 19997:  20 0.22 19997  207  207     1
## 19998:  20 0.22 19998   70   70     1
## 19999:  20 0.22 19999  669  365     0
## 20000:  20 0.22 20000  994  365     0

defc <- defData(varname = "b", formula = 0, variance = 0.100)

dc <- genData(2000, defc, id = "pcp")
dd <- genCluster(dc, "pcp", numIndsVar = 500, "id")
dd <- genSurv(dd, defs, digits = 0)
dd <- addColumns(defa, dd)

ds <- dd[, .(p = mean(event)), keyby = pcp]
ds
##        pcp    p
##    1:    1 0.41
##    2:    2 0.43
##    3:    3 0.36
##    4:    4 0.51
##    5:    5 0.50
##   ---          
## 1996: 1996 0.82
## 1997: 1997 0.53
## 1998: 1998 0.88
## 1999: 1999 0.42
## 2000: 2000 0.63

ds[, .(mu = mean(p), s = sd(p), cv = sd(p)/mean(p))]
##     mu    s   cv
## 1: 0.5 0.11 0.22
ds[, .(quantile(p, probs = c(0.025, .975)))]
##      V1
## 1: 0.31
## 2: 0.72

Evaluating the effect size

Next, we generate data that includes treatment assignment (but excludes cluster variability and censoring before one year). The treatment effect is expressed as a log hazard ratio, which in this case 0.4 (equal to a hazard ratio of just about 1.5).

The data generation starts with treatment assignment, adds the time-to-event survival data, and then adds the one-year censoring data, as before:

defa <- defData(varname = "rx", formula = "1;1", dist = "trtAssign")

defs <- defSurv(varname = "tte", formula = "r[1] + 0.4 * rx", shape = r[2])

defe <- defDataAdd(varname = "time", formula = "min(365, tte)", dist = "nonrandom")
defe <- defDataAdd(defe, "event", "1*(tte <= 365)", dist = "nonrandom")

dd <- genData(1000, defa)
dd <- genSurv(dd, defs, digits = 0)
dd <- addColumns(defe, dd)

dd
##         id rx  tte time event
##    1:    1  0 1902  365     0
##    2:    2  0  286  286     1
##    3:    3  0  480  365     0
##    4:    4  0   32   32     1
##    5:    5  1   12   12     1
##   ---                        
##  996:  996  0  663  365     0
##  997:  997  1  962  365     0
##  998:  998  1   19   19     1
##  999:  999  1  502  365     0
## 1000: 1000  0   85   85     1

The plot of the survival curves by treatment arms provides a visualization of the treatment effect:

Again, since we have no censoring, we can estimate the probability of an event within 365 days for each arm:

dd[, .(p = mean(event)), keyby = rx]
##         rx    p
## 1: control 0.51
## 2:      rx 0.63

An effect size of 0.4 on the log hazard scale translates to an odds ratio of about 1.64, a risk ratio of 1.24, and risk difference of 12 percentage points. In the absence of any external data on the potential effect size, we can use an effect size that is minimally clinically meaningful based on any or all these effect size measurements.

Complete data generation and model estimation

With the pieces in place, we are ready to put it all together and add censoring to the mix to finalize the full data generating process. We will fit a mixed-effects Cox proportional hazards model to see if we can recover the parameters that we have used to generate the data, and, if that goes well, we will be ready to estimate power.

We start in defc by defining the cluster-level random effect variation and treatment assignment design (in this case 1 to 1, treatment to control). We add a censoring process in defa. This assumes that we will be enrolling patients for six months, spread out across this time period. The study will last exactly one year, so every patient will be followed for at least six months, and only some will be followed for a full year (i.e., those who join the study on the first day).

Finally, defs defines the data generation process for the survival outcome, which we’ve seen above, though now we have both a treatment effect and a random effect, in addition to the baseline parameters in the vector r.

defc <- defData(varname = "b", formula = 0, variance = 0.05)
defc <- defData(defc, varname = "rx", formula = "1;1", dist = "trtAssign")

defa <- defDataAdd(varname = "start_day", formula = "1;182", dist = "uniformInt")
defa <- defDataAdd(defa, varname = "censor", 
  formula = "365 - start_day ", dist = "nonrandom")

defs <- defSurv(varname = "tte", formula = "r[1] + 0.4 * rx + b", shape = r[2])

The data generation is the same as before, though there is an additional censoring process, which is done with the function addCompRisk:

dc <- genData(500, defc, id = "pcp")
dd <- genCluster(dc, "pcp", numIndsVar = 200, "id")
dd <- addColumns(defa, dd)
dd <- genSurv(dd, defs, digits = 0)
dd <- addCompRisk(dd, events = c("tte", "censor"), 
  timeName = "time", censorName = "censor", keepEvents = TRUE)

dd
##         pcp       b rx     id start_day censor  tte time event   type
##      1:   1 -0.0158  1      1        65    300  307  300     0 censor
##      2:   1 -0.0158  1      2       113    252  907  252     0 censor
##      3:   1 -0.0158  1      3       121    244  110  110     1    tte
##      4:   1 -0.0158  1      4        60    305  170  170     1    tte
##      5:   1 -0.0158  1      5        95    270 2673  270     0 censor
##     ---                                                              
##  99996: 500 -0.0015  0  99996       160    205  291  205     0 censor
##  99997: 500 -0.0015  0  99997        79    286  193  193     1    tte
##  99998: 500 -0.0015  0  99998        64    301   28   28     1    tte
##  99999: 500 -0.0015  0  99999       153    212  459  212     0 censor
## 100000: 500 -0.0015  0 100000        19    346  298  298     1    tte

Since we have generated a rather large data set, we should be able to recover the parameters pretty closely if we are using the correct model. We are going to fit a mixed effects survival model (also known as a frailty model) to see how well we did.

fit_coxme <- coxme(Surv(time, event) ~ rx + (1 | pcp), data = dd)
summary(fit_coxme)
## Cox mixed-effects model fit by maximum likelihood
##   Data: dd
##   events, n = 49955, 1e+05
##   Iterations= 11 48 
##                   NULL Integrated  Fitted
## Log-likelihood -555429    -553708 -553058
## 
##                   Chisq  df p  AIC  BIC
## Integrated loglik  3442   2 0 3438 3421
##  Penalized loglik  4743 414 0 3914  259
## 
## Model:  Surv(time, event) ~ rx + (1 | pcp) 
## Fixed coefficients
##    coef exp(coef) se(coef)  z p
## rx 0.39       1.5    0.022 18 0
## 
## Random effects
##  Group Variable  Std Dev Variance
##  pcp   Intercept 0.22    0.05

Pretty good! The estimated HR of 0.39 (95% CI: 0.35 - 0.43) is on target (we used 0.40 in the data generation process), and the estimated variance for the PCP random effect was 0.05, also on the mark. I’d say we are ready to proceed to the final step.

Power estimation

To conduct the power estimation, I’ve essentially wrapped the data generation and model estimation code in a collection of functions that can be called repeatedly to generate multiple data sets and model estimates The goal is to calculate the proportion of data sets with a statistically significant result for a particular set of assumptions (i.e., the estimate of power for the assumed effect size, variation, and sample sizes). I’ve provided the code below in the addendum in case you haven’t grown weary of all this detail. I described a general framework for using simulation to estimate sample size/power, and I’m largely following that process here.

I’ve written a little function scenario_list (which I’m now thinking I should add to simstudy) to create different parameter combinations that will determine the power estimation. In this case, the parameters I am interested in are the number of PCPs that should be randomized and the variance assumption. The number of patients per PCP (cluster size) is also important to vary, but for illustration purposes here I am keeping it constant.

Here is the simplified scenario list with four possible combinations:

scenario_list <- function(...) {
  argmat <- expand.grid(...)
  return(asplit(argmat, MARGIN = 1))
}

npcps <- c(20, 30)
npats <- c(15)
s2 <- c(0.03, 0.04)

scenarios <- scenario_list(npcps = npcps, npats = npats, s2 = s2)

scenarios
## [[1]]
## npcps npats    s2 
## 20.00 15.00  0.03 
## 
## [[2]]
## npcps npats    s2 
## 30.00 15.00  0.03 
## 
## [[3]]
## npcps npats    s2 
## 20.00 15.00  0.04 
## 
## [[4]]
## npcps npats    s2 
## 30.00 15.00  0.04

I use the mclapply function in the parallel package to generate three iterations for each scenario:

model.ests <- mclapply(scenarios, function(a) s_scenarios(a, nrep = 3))

model.ests
## [[1]]
##    npcps npats   s2 est_s re.var_s   p_s
## 1:    20    15 0.03  0.41   0.1596 0.100
## 2:    20    15 0.03  0.32   0.0306 0.065
## 3:    20    15 0.03  0.32   0.0004 0.041
## 
## [[2]]
##    npcps npats   s2 est_s re.var_s     p_s
## 1:    30    15 0.03  0.13    0.017 3.5e-01
## 2:    30    15 0.03  0.62    0.029 7.3e-05
## 3:    30    15 0.03  0.45    0.040 3.8e-03
## 
## [[3]]
##    npcps npats   s2 est_s re.var_s    p_s
## 1:    20    15 0.04  0.55  7.9e-02 0.0084
## 2:    20    15 0.04  0.07  4.0e-04 0.6700
## 3:    20    15 0.04  0.24  8.2e-05 0.1500
## 
## [[4]]
##    npcps npats   s2 est_s re.var_s     p_s
## 1:    30    15 0.04  0.42  8.3e-05 2.7e-03
## 2:    30    15 0.04  0.57  4.0e-04 2.6e-05
## 3:    30    15 0.04  0.36  1.0e-01 4.2e-02

In the actual power calculations, which are reported below, I used 60 scenarios defined by these data generation parameters:

npcps <- c(18, 24, 30, 36)
npats <- c(15, 18, 21)
s2 <- c(0.01, 0.02, 0.03, 0.04, 0.05)

For each of these scenarios, I generated 5000 data sets and estimated models for each (i.e., a total of 300,000 data sets and model fits). For each of the 60 scenarios, I estimated the proportion of the 5000 model fits that yielded a p-value < 0.05 for the estimated log hazard ratio. I had the benefit of using a high performance computer, because running this on my laptop would have taken well over 10 hours (only about 10 minutes on the HPC).

At the end, we have a plot of “power curves” that shows estimated power for each of the scenarios. If we assume that we can expect at least 18 patients per PCP and that the between-PCP variance will be around 0,03 or 0.04, we should be OK randomizing 30 PCPs (15 in each arm), though it might more prudent to go with 36, just to be safe:

Addendum

Here is the code I used to generate the data for the power curve plot. It is based on the framework I mentioned earlier. There is one extra function here, extract_coxme_table, which I pulled from stackoverflow, because there is currently no obvious way to extract data from the coxme model fit.

extract_coxme_table <- function (mod) {
  beta <- mod$coefficients 
  nvar <- length(beta)
  nfrail <- nrow(mod$var) - nvar
  se <- sqrt(diag(mod$var)[nfrail + 1:nvar])
  z <- round(beta/se, 2)
  p <- signif(1 - pchisq((beta/se)^2, 1), 2)
  table=data.table(beta = beta, se = se, z = z, p = p)
  return(table)
}

s_def <- function() {
  
  defc <- defData(varname = "b", formula = 0, variance = "..s2")
  defc <- defData(defc, varname = "rx", formula = "1;1", dist = "trtAssign")
  
  defa <- defDataAdd(varname = "start_day", formula = "1;182", dist = "uniformInt")
  defa <- defDataAdd(defa, varname = "censor", 
                     formula = "365 - start_day ", dist = "nonrandom")
  
  defs <- defSurv(varname = "tte", formula = "-4.815 + 0.4 * rx + b", shape = 1.326)
  
  defa2 <- defDataAdd(varname = "event6", 
                      formula = "1*(tte <= 182)", dist = "nonrandom")
  
  return(list(defc = defc, defa = defa, defs = defs, defa2 = defa2))
  
}

s_generate <- function(argsvec, list_of_defs) {
  
  list2env(list_of_defs, envir = environment())
  list2env(as.list(argsvec), envir = environment())
  
  dc <- genData(npcps, defc, id = "pcp")
  dd <- genCluster(dc, "pcp", npats, "id")
  dd <- addColumns(defa, dd)
  dd <- genSurv(dd, defs, digits = 0)
  dx <- addCompRisk(dd, events = c("tte", "censor"), 
                    timeName = "time", censorName = "censor", keepEvents = TRUE)
  dx <- addColumns(defa2, dx)
  
  dx[]
  
}

s_replicate <- function(argsvec, list_of_defs) {
  
  dx <- s_generate(argsvec, list_of_defs)
  
  coxfitm <-coxme(Surv(time, event) ~ rx + (1 | pcp), data = dx)
  
  list2env(as.list(argsvec), envir = environment())
  
  return(data.table(
    npcps = npcps,
    npats = npats,
    s2 = s2,
    est_s = fixef(coxfitm), 
    re.var_s = VarCorr(coxfitm)$pcp,
    p_s = extract_coxme_table(coxfitm)$p
  ))
  
}

s_scenarios <- function(argsvec, nreps) {
  
  list_of_defs <- s_def()
  
  rbindlist(
    parallel::mclapply(
      X = 1 : nreps, 
      FUN = function(x) s_replicate(argsvec, list_of_defs), 
      mc.cores = 4)
  )
  
}