We demonstrate the basic simulation capabilities for time-to-event clinical trials that are implemented in the nphsim
package.
This includes:
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 to
rpwexp, it is designed to simulate enrollment with rates specified in
rate`.
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.
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())
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)))
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:
sim
)analysis
; e.g., a group sequential design would have multiple analysest
)NE
and NC
, respectively)DE
) and control (DC
) arms as well as the overall number of events at the analysis (D
)pval
), z
from the quantile function of the standard normal distribution based on the one-sided p-value), hr
); note that the sign for the above logrank z-value depends on this with negative z
corresponding to hr
< 1), and 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.
pval_FH(0,1)
, pvale_FH(1,1)
- Fleming-Harrington tests using the specified input values for gamma
and rho
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.