In keaven/nphsim: Non proportional hazards sample size and simulation

We demonstrate the basic simulation capabilities for time-to-event clinical trials that are implemented in the nphsim package. This includes:

Simulating enrollment with piecewise constant enrollment rates
Simulating times-to-events with a piecewise exponential distribution
Simulating a clinical trial (fixed design) by combining the above capabilities
Performing various statistical tests on a simulated clinical trial

Simulating Enrollment

We use piecewise constant enrollment and piecewise exponential failure rates as methods to approximate arbitrary patterns of enrollment and dropouts. We specify an example non-proportional hazards scenario.

We approximate an arbitrary enrollment pattern by assuming constant enrollment rates within fixed time periods. The length of enrollment periods prior to the final stable enrollment (2, 4, respectively) are specified here; the final duration of 12 in the third period may be extended or shortened to achieve enrollment.

enrollIntervals <- c(2, 4)

Thus, we assume enrollment is constant in months 1-2, 3-6, and from month 7 onward. Next we specify relative rates of enrollment in these time periods:

enrollRates <- c(.5, 3, 16)

Let's say we wish to enroll 200 patients, starting with 7 patients per month and doubling after month 2 and month 6 as specified above.

We generate random enrollment times as follows:

library(nphsim)
enrollTimes <- rpwexp(n = 300, rate = 7 * enrollRates, intervals = enrollIntervals, cumulative = TRUE)

When 'cumulative = TRUEis input torpwexp, it is designed to simulate enrollment with rates specified inrate`.

library(ggplot2)
qplot(x = c(0, enrollTimes), y = 0:length(enrollTimes), geom="step", 
      ylab="Number enrolled", xlab="Time") +
      scale_x_continuous(breaks=c(0,6,12,18))

You can see the enrollment rate increases from the first 2 months to months 3-6 and again in the period after month 6. The last patient in this simulation is enrolled at max(enrollTimes) = r round(max(enrollTimes),1) months.

Generating Piecewise Exponential Failure Times

We can also use the rpwexp function to generate piecewise exponential failure times. We have a different set of intervals to define the periods of time for which different failure rates apply:

# Failure rates for piecewise exponential time periods
failRates <- c(.3, .6, .5)
# Interval duration(s) before final stable rate
# Note that length is 1 fewer than for failRates
# and should be NULL if there is only 1 failure rate
failIntervals <- c(1,4)

Now we generate the piecewise exponential failure times corresponding to the patient enrollment previously simulated, again using rpwexp. We do not need to specify the default value cumulative = FALSE which distinguishes this from the call generating enrollment times.

n <- length(enrollTimes)
y <- rpwexp(n = n, rate = failRates, intervals = failIntervals)

The longest duration among these is max(y) = r round(max(y),1). Although we will eventually be interested in generating censored survival times, this is not done here; there is a value for event time each simulated patient. We plot the duration of time until an event for all patients ordered by decreasing time which reproduces the shape of the Kaplan-Meier curve for the points; the theme function removes the y-axis as well as its labels and tick-marks. Note that at random there are some patients with very long simulated times.

id <- 1:n
# sort patient time-to-event and create a 0 starting point for each patient
dta <- data.frame(N=factor(c(id,id)),Time=c(array(0,n), sort(y, decreasing = TRUE)))
ggplot(dta, aes(x = Time, y = N, grp = N)) + geom_line() +
       xlab("Time-to-event") + ylab("Patients ordered by time-to-event") +
  scale_x_continuous(breaks=(0:4)*6) +
  theme(axis.text.y=element_blank(),axis.ticks.y=element_blank())

Simulating a 2-Arm Trial Instance with Non-Proportional Hazards

We now proceed to simulating a clinical trial with the control group rate specified above and with an experimental group having piecewise exponential rates generated according to a non-proportional hazards assumption. The parameters here appear somewhat different than above in rpwexp examples as they follow variable names used in the gsDesign R package. We make the sample size large so that the resulting estimated survival should approximate the underlying distributions closely.

# Hazard ratio corresponding to control group event rates above
hr <- c(1, .6, .3)
# Sample size of control and experimental arms
ssC <- 5000
ssE <- 5000
# We assume dropouts at a rate of .1 per time unit, increasing to .3 at end
# dropout rates specified for same intervals as failure rates
eta <- c(.1,.1,.3)
# Simulate a single trial instance
# We add an arbitrary enrollment interval for indefinite enrollment duration
trial <- nphsim(nsim = 1, lambdaC = failRates, lambdaE = failRates*hr,
                intervals = failIntervals, ssC = ssC, ssE = ssE,
                gamma = 2 * enrollRates, R = enrollIntervals, fixEnrollTime = FALSE,
                eta = eta, etaE = eta)
# show a few lines
head(trial$simd, n=5)

In the above, we see the simulation ID in sim, treatment group in treatment, time of study entry in enterT, duration followed for survival in survival and censoring value (0 = event, 1 = censored) in cnsr. We examine censoring by treatment group:

table(trial$simd$treatment, trial$simd$cnsr)

Summarizing times of events relative to start of enrollment, we have

summary(trial$simd$enterT+trial$simd$survival)

We compare the distributions of times to events by computing Kaplan-Meier curves. The curves appear to separate more as time goes on and the hazard ratio becomes more extreme. We would expect no separation over the first r enrollIntervals[1] where the hazard ratio is 1, with increasing separation over time.

library(survival)
plot(with(trial$simd, survfit(Surv(survival, 1-cnsr) ~ treatment)))

Performing an analysis

Now suppose we wish to test after 900 events. First, we test with a logrank statistic.

testOut <- simtest(x = trial, anaD=900, method='LR')
testOut$result

The values after that you see above are:

ID number (sequential) for each simulation (sim)
analysis number for the analysis (analysis; e.g., a group sequential design would have multiple analyses
the time since start of randomization until the analysis (t)
the simulation sample size in the experimental and control arms (NE and NC, respectively)
the number of events in the expermental (DE) and control (DC) arms as well as the overall number of events at the analysis (D)
the logrank test p-value (pval),
a Z-value corresponding to the logrank test (z from the quantile function of the standard normal distribution based on the one-sided p-value),
the hazard ratio from a Cox model for experimental compared to control treatment (hr); note that the sign for the above logrank z-value depends on this with negative z corresponding to hr < 1), and
the standard error for the hazard ratio estimate (sehr).

Now we consider several weighted logrank test statistics at the same time. For this we apply method=wlr.Stat below; we will explain the output after showing the code.

testOut <- simtest(x = trial, anaD=900, method=wlr.Stat, fparam=list(rho=c(0,1), gamma=c(1,1)))
testOut$result

The first columns through D are as explained above. Following this are one-sided p-values for Fleming-Harrington Test with sequence of $\rho$ and $\gamma$ parameters specified.