Advanced models for predicting event occurrence in event-driven clinical trials accounting for patient dropout, cure and ongoing recruitment

08/20/2021
by   Vladimir Anisimov, et al.
0

We consider event-driven clinical trials, where the analysis is performed once a pre-determined number of clinical events has been reached. For example, these events could be progression in oncology or a stroke in cardiovascular trials. At the interim stage, one of the main tasks is predicting the number of events over time and the time to reach specific milestones, where we need to account for events that may occur not only in patients already recruited and are followed-up but also in patients yet to be recruited. Therefore, in such trials we need to model patient recruitment and event counts together. In the paper we develop a new analytic approach which accounts for the opportunity of patients to be cured, as well as for them to dropout and be lost to follow-up. Recruitment is modelled using a Poisson-gamma model developed in previous publications. When considering the occurrence of events, we assume that the time to the main event and the time to dropout are independent random variables, and we have developed a few advanced models with cure using exponential, Weibull and log-normal distributions. This technique is supported by well developed, tested and documented software. The results are illustrated using simulation and a real dataset with reference to the developed software.

READ FULL TEXT VIEW PDF

Authors

page 1

page 2

page 3

page 4

02/14/2022

Modelling and forecasting patient recruitment in clinical trials with patients' dropout

This paper focuses on statistical modelling and prediction of patient re...
11/02/2021

Comparison of Time-to-First-Event and Recurrent Event Methods in Multiple Sclerosis Trials

Suppression of disability progression is an important goal in the treatm...
06/01/2018

Adversity Index for Clinical Trials: An Inclusive Approach for Analysis of Safety Data

This paper proposes a new method for analysis of adverse event data in c...
09/06/2019

A simulation study of methods for handling disease progression in dose-finding clinical trials

In traditional dose-finding studies, dose-limiting toxicity (DLT) is det...
12/02/2018

In-silico Risk Analysis of Personalized Artificial Pancreas Controllers via Rare-event Simulation

Modern treatments for Type 1 diabetes (T1D) use devices known as artific...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

An important aspect of event-driven trials is the operational design at the initial and interim stages, i.e. predicting the event counts over time and the time to reach specific milestones, accounting for events that may occur not only in patients already recruited and are followed-up but also in patients yet to be recruited. Therefore, in event-driven trials we need to model patient recruitment and event counts together.

There are different techniques for recruitment modelling described in the literature and one of the main directions is using mixed Poisson models. This direction has a long history, with several papers devoted to the use of Poisson processes with fixed recruitment rates to describe the recruitment process (Carter et al., [16]; Senn [17, 18]). However, in real clinical trials the recruitment rates in different centres vary. Therefore, to model this variation, Anisimov and Fedorov[1]

introduced a Poisson-gamma model, where the variation in rates in different centres is modelled using a gamma distribution (see also Anisimov

[3]). Some applications to real trials are considered in Anisimov et al.[2]. This technique was developed further in several publications for predicting and interim re-forecasting the recruitment process under various conditions Anisimov[4, 8]. Other approaches to recruitment modelling primarily deal with global recruitment. These approaches use different techniques and we refer interested readers to survey papers by Barnard et al.[11]; Heitjan et al.[15]) and Gkioni et al.[14], and also to a discussion paper Anisimov[7] on using Poisson models with random parameters with other references therein.

A larger number of clinical trials are event-driven where the number of clinical events is required to be large enough to allow for reliable statistical conclusions about the parameters of patient responses. For such trials, one of the main tasks is predicting not only the required number of recruited patients, but also the number of events that may occur and the time to reach particular milestones. A useful review of different approaches for event-driven trials is provided in Heitjan et al.[15]. However, for predicting the number of events over time and the time to stop the trial, the authors of papers cited primarily use a Monte Carlo simulation technique, e.g. Bagiella and Heitjan[10]. Therefore, Anisimov[5] developed an analytic methodology for predictive modelling the event counts together with patient recruitment in ongoing event-driven trials accounting also for patient dropout. This methodology is developed further to forecasting multiple events at start-up and interim stages under exponential assumptions, Anisimov[8], and predicting some operational characteristics during follow-up times, Anisimov[6].

It is also of interest in event-driven trials that a number of patients under treatment will not experience the event within their exposure time (i.e. time from randomisation to a particular milestone) as the therapy of different diseases is improving. Therefore, an interesting direction is to investigate the opportunity of cure. Here we note the paper by Chen[13], but he also uses a simulation technique for predicting event timing and does not consider the case of patient dropout.

Therefore, in the paper we developed a new analytic approach to this problem which accounts for patient dropout and also for the opportunity for patients to be cured with some probability.

We assume that the patient recruitment is modelled using a Poisson-gamma model developed in Anisimov and Fedorov[1], Anisimov[4]. We consider non-repeated events and assume that the times to the main event and to dropout are independent random variables, and there is an opportunity of cure. Several new models have been developed using exponential, Weibull and log-normal distributions. The focus is on the interim stage, where the parameters of different models are estimated using maximum likelihood technique. The predictive distributions of the number of future events for all considered models are derived in the closed forms, thus, Monte Carlo simulation is not required.

The developed technique and R-tools allow for forecasting the event counts over time and also the time to stop the trial with mean and predictive bounds. The results are illustrated in the paper using Monte Carlo simulation and the real dataset.

The paper is organised as follows:

Section 2, the basic models for the process of event occurrence; Section 3, predicting event counts for patients at risk; Section 4, predicting event counts accounting for ongoing recruitment; Section 5, testing of the Weibull model with cure using Monte Carlo simulation; Section 6, software development and an R-package; Section 7, implementation to a real clinical trial, and; Section 8, fitting models to real data.

2 Modelling the process of event occurrence

Consider a trial at some interim time and assume that there is one type of non-repeated events:  the main event of interest , and the patients also can be lost to follow-up (call it dropout).

Then all patients that were recruited in the trial until a given interim time can be divided into three groups:

1) group : patients experienced event . Denote by the total number of patients in this group and by the lengths of follow-up periods from randomisation date until the event;

2) group : each patient is censored at interim time, thus the patients have neither experienced an event nor are they lost to follow-up. Denote by the total number of patients, and by the lengths of follow-up periods from randomisation date until interim time;

3) group : patients are lost to follow-up. Denote by the total number of patients, and by the lengths of follow-up periods until censoring by dropout.

Consider now the following cure model describing the process of event occurrence for every patient.

Assume that after randomisation, the patient can be either cured with some probability , or with probability can experience event after some random time . If a patient is cured, then event cannot occur.

Denote time to dropout , and if event doesn’t occur before , the patient experiences dropout regardless of whether this patient is cured or not (in this case event cannot occur).

Assume that the events for different patients occur independently and the times and

are also independent random variables with cumulative distribution functions (CDF)

and , respectively. Suppose that , and probability of cure are the same for all patients, though potentially we can consider different treatment groups with different parameters. Assume also that these functions are continuously differentiable and denote by and

the corresponding probability density functions (pdf).

2.1 Estimating parameters of the model

Consider maximum likelihood method. Denote for convenience, and .

For a patient in group with exposure time , the probability that event and dropout will not occur is  .

For a patient in group with exposure time , the probability that event occurs in a small interval before dropout is .

For a patient in group with exposure time , the probability that dropout occurs in a small interval before event is .

Given data, the maximum likelihood function has the form

Correspondingly, the log-likelihood function is

For different types of distributions this expression will have a different form.

2.1.1 Exponential with cure model

This model assumes that the variables and

are exponentially distributed with rates

and respectively. This is a three parameter model: . The log-likelihood function:

where and .

Consider equating the partial derivatives of the log-likelihood function to zero to find a relationship between parameters. Partial derivatives are:

Equating the last derivative to zero, we get that

Substituting into relation (2.1.1) we get a simpler relation depending only on two variables :

To find the estimators, optimisation is carried out by maximising the log-likelihood function. Initial values are set as and taken to be some range of values in . In optimisation, new variables are considered:

After optimisation, the variables are transformed back to the original parameters.

2.1.2 Weibull with cure model

By definition, the pdf and CDF of a Weibull distribution are

where are shape and scale parameters. For ease of notation, we use the parametrisation . Then pdf and CDF have the form

Weibull with cure model assumes that the variables and have Weibull distribution with parameters ( and respectively. This is a five parameter model: . The log-likelihood function:

Optimisation is carried out in the same way as for the exponential model, with initial values: taken to be some range of values in . The new variables are:

Similar relations can be written for the combination of the distributions, e.g. Weibull distribution for time to event and exponential distribution for time to dropout , and vice versa.

Note that the Weibull model is in some sense a generalisation of the exponential model. Indeed, if in particular in the relations above we fix the value , then we get the combined exponential-Weibull with cure model (time to event has an exponential distribution). By setting both values, and , we get the exponential with cure model.

In a similar way the log-likelihood function can be derived also for a log-normal with cure model.

3 Predicting event counts for patients at risk

Let us introduce for convenience the time of the occurrence of event , , so . Note that if , then is an improper random variable as .

Consider a conditional probability for a patient in group to experience an event in the future time interval given that the follow-up period until the interim time is :

For the exponential model can be calculated in a closed form:

(3)

where .

Note that for the exponential model, if , , so this expression does not depend on and we have a memoryless property. However, for , the memoryless property is lost.

For the Weibull with cure model has the following form:

(4)

where

(5)

To compute this function in applications we can use a numerical integration.

Similar relations can be written for the combination of the distributions, and also for a log-normal with cure model.

3.1 Global prediction

Assume now that the recruitment of new patients is already completed, thus, the events in the future may occur only in patients at risk in group . Denote by the total predictive number of events that may occur in future time interval   for patients in group where are the times of exposure. Let be a Bernoulli random variable, .

Lemma 3.1

The process can be represented in the form:

(6)

where the variables are independent and the probability is defined above in Section 3 and depends on the type of the distributions used in the event model.

For a rather large number of patients in group , , we can apply a normal approximation for the process

using simple formulae for the mean and the variance:

(7)

Then and -predictive interval at time is , where is an

-quantile of a standard normal distribution.

For a not so large number of patients, a distribution of can be calculated numerically as a convolution of the sum of Bernoulli variables.

Let us evaluate the predictive distribution for the time to reach a given target for the total planned number of events in the study.

Recall that in previous notation denotes the total number of events that occurred prior to interim time (size of group ). The remaining number of events that are left to achieve is .

Let be the remaining time to reach events after the interim time . Then the following relation holds: for any ,

(8)

As the distribution of can be evaluated for any time , this relation allows us to calculate also the distribution of .

Consider the calculation of PoS (probability to complete study before a planned time ). Denote it as . From (8) we get

(9)

If we use a normal approximation for the process , then

(10)

where is the CDF of a standard normal distribution.

4 Predicting event counts accounting for ongoing recruitment

Consider now the situation when at the interim time the planned number of patients to be recruited is not reached yet, that means, the recruitment is still ongoing. In this case we need also to predict the future recruitment and how many events may occur for patients to be recruited in the future.

4.1 Modelling and predicting patient recruitment

Assume that patients arrive at clinical centres according to Poisson processes with some rates . To model the variation in the rates among different centres we assume that are jointly independent gamma distributed random variables with parameters (shape and rate) and pdf

(11)

where is a gamma function.

This model is called a Poisson-gamma (PG) recruitment model and was developed in Anisimov & Fedorov[1] and further extended in Anisimov[3, 4, 8].

Denote by a standard Poisson process with rate and by

a random variable which has a Poisson distribution with parameter

. Then a mixed Poisson process where the rate is gamma distributed with parameters is a PG process (Bernardo and Smith[12]) with parameters :

(12)

Note that for a mixed Poisson process with random rate ,

(13)

Assume now that some centre is active only in time interval . Denote by the duration of recruitment window (duration of active recruitment) in a centre up to time :

(14)

Assume that the recruitment rate in this centre is which is gamma distributed with some parameters. Then the recruitment process in this centre for any can be represented as a PG process with a cumulative rate . That means, the number of patients recruited in interval has a mixed Poisson distribution with the rate .

Consider now predicting the remaining recruitment at some interim time . Assume for simplicity that all centres are active and in every centre the following data are available:   - the duration of active recruitment (recruitment window) and the number of patients recruited.

In Anisimov and Fedorov[1] (see also Anisimov[4]) it was developed a maximum likelihood technique for estimating parameters of a PG model assuming that in all active centres the rates have a gamma distribution with the same parameters. In [4], the Bayesian technique was also developed for predicting future recruitment using the property that the posterior rate in a centre , , which is adjusted to the data in this centre, also has a gamma distribution with parameters .

Consider now a given interim time . Let be some active centre. Denote by the parameters of a PG model estimated using data in all active centres as noted above. Then the future recruitment process in centre can be modelled as a PG process with posterior recruitment rate . Assume that the recruitment in this centre can be closed due to some operational reasons at some time . Then for any the recruitment process in centre in time interval can be represented as a PG process with a cumulative rate .

Assume now that is some new centre that is planned to be initiated at time and let be the closing date of recruitment in this centre. Denote by the recruitment rate in this centre. Note that the rates in the new centres can be provided by clinical teams using expert estimates or evaluated using historical data from similar trials.

Then centre will be active only in time interval . Thus, for any , the recruitment process in time interval can be represented as a PG process with a cumulative rate .

Consider the prediction of the remaining global recruitment.

Denote by a set of active centres with posterior rates . Assume also that it can be some set of new centres that are planned to be initiated after interim time at times . Denote by the rates in the new centres. Then the predictive total number of patients to be recruited in the time interval can be represented as

(15)

This means, has a mixed Poisson distribution with a cumulative rate

(16)

For a rather large number of centres, the predictive bounds for can be evaluated using a normal approximation, as the mean and the variance of can be easily calculated using the property (13) and relations   . In particular, the mean predicted time to reach a required remaining number of patients can be numerically calculated as the point when the line hits level .

Note that for a not so large number of centres, for predicting one can use a PG approximation developed in Anisimov[8], Anisimov and Austin[9].

4.2 Predicting event counts

Consider now predicting event counts accounting for ongoing recruitment.

Denote by the time it takes until event occurs first (before dropout), and let , , be its CDF.

For cure model with dropout defined in Section 2, in previous notation,

(17)

In particular, for the exponential with cure model, using notation ,

(18)

For Weibull model, using parametrisation ,

(19)

where

(20)

Consider now one clinical centre. Assume that patients arrive according to a mixed Poisson process with possibly random rate . Assume also that the centre is active only in a fixed time interval . In Anisimov[5, 8] the following result is proved.

Lemma 4.1

The predicted number of events in interval that occur in the newly recruited patients in this centre has a mixed Poisson distribution with rate , where

(21)

For the exponential model, the function can be easy calculated. Consider the duration of recruitment window in a centre at time defined in (14). Then, using parameters ,

(22)

For Weibull Model with parameters ,

(23)

where is defined in (20). This function can be numerically calculated.

Similar relations in the integral form can be written for the combination of the distributions, and also for a log-normal with cure model.

These results form the basis for creating predictions of the event counts in any active centre and globally.

4.3 Global forecasting event counts at interim stage

Consider now forecasting the total number of events at some interim time .

Denote the times of initiation of new centres (if any) by and the times of closure for all centres by . In general it is assumed that centres will be closed for recruitment at the time when recruitment hits the recruitment target. Thus, in applications, we usually assume that where is the predicted mean remaining time to reach the recruitment target.

Theorem 4.2

The predictive total number of new events , , that may occur in future time interval , can be represented as a convolution of two independent random variables:

(24)

where according to (16),

(25)

and is the predictive number of events in group defined in (6), Section 3.1.

Here the function is defined by (21) (for exponential and Weibull models we have the expressions (22) and (23), respectively). The first sum in (25) is taken across all active centres and are the posterior rates defined in Section 4.1. and the second sum is taken across new centres.

Correspondingly, the probability to complete trial in time is

(26)

where is the planned remaining time to complete the trial and is the remaining number of events left to achieve.

The proof follows from results of Lemmas 3.1, 4.1 and Section 4.1.

Note that the mean and the variance of the process can be calculated explicitly in terms of functions and parameters of the rates.

As typically in real trials the number of centres is rather large, to create predictive bounds for one can use a normal approximation. This technique is realised in R package (EventPrediction), see Section 6.

5 Monte Carlo simulation

Monte Carlo simulation was used to test each model’s performance. We considered 1000 patients assuming uniform distribution of centre initiation over 6 months, and took the target number of events 550. At a specified cut-off date the model parameters were estimated using maximum likelihood estimation, see Section

2.1. Using these estimators, predictions of the future occurrence of events were created. For the Weibull model two different cases for the initial parameters were considered, and , see Fig 1 and Fig  2 respectively.

Figure 1: Plot of number of events against time (in days), following the timeline of a simulated trial. The simulated trajectory of events A is marked by the black solid line, the initial parameters: , , , and . An interim analysis was taken at 7 months, the estimated parameters: , , , and . Predictions on future event counts were created using the estimated parameters; the mean trajectory is shown by the blue dashed line, the 90% confidence bounds by the red dotted lines.
Figure 2: Plot of number of events against time (in days), following the timeline of a simulated trial. The simulated trajectory of events A is marked by the black solid line, the initial parameters: , , , and . An interim analysis was taken at 7 months, the estimated parameters: , , , and . Predictions on future event counts were created using the estimated parameters; the mean trajectory is shown by the blue dashed line, the 90% confidence bounds by the red dotted lines.

In both cases, the model successfully predicts the trajectory of the number of events with the real trajectory falling within the 90% predicted bounds. Furthermore, the parameters estimated at the cut-off time using maximum likelihood technique are close to the initial parameters showing an appropriate estimation.

6 Software development

In order to expose the event and recruitment prediction models (as detailed in the previous sections) to a large number of key stakeholders, an R package (EventPrediction) has been developed, tested and deployed to a centralised R server. The EventPrediction package allows a user to easily pass the data required (i.e. subject event data, centre level data and configuration) and return back key parameter estimates and predictions with bounds, for both events and recruitment.

6.1 R package design

R was chosen over other programming languages, and an R Package was developed over standalone R scripts, for a number of reasons, including: R has an easy to use and to setup testing framework; R ships with easy to use code coverage tools; R has a Comprehensive R Archive Network ("CRAN") set of packages that are easily accessible; R seamlessly integrates with GitLab (and other source control software); R allows for an Object Oriented ("OO") approach (i.e. S3, S4 and R6); and also, primarily, it is simply straightforward to develop, test, document and centrally deploy an R Package for key stakeholders to use.

R’s S3 lightweight OO solution was a key design feature of the EventPrediction package, as using such an OO approach yields four main benefits: first, S3’s simple to use OO benefit of polymorphism (aka in R as method dispatch); secondly, S3 gives the OO benefit of inheritance; thirdly, S3 is ubiquitously used by R contributors, easy to use and for others to comment; and fourthly, S3 is in accordance with the functional programming paradigm, when an object is passed into an S3 function it is not going to change (unlike full OO approaches like R6).

6.1.1 Good software engineering principles

Another major benefit of developing an R Package (and utilising R’s OO approach) is to ensure adherence to good Software Engineering principles, with code that is at a minimum: reliable; easy to use; efficient; well tested, with tests traceable to requirements and/or design; well documented; and (importantly) easy to maintain. The EventPrediction package conforms to each of these key programming elements, not only because these are simply good Software Engineering practices, but also as the biotechnology sector is highly regulated and there is a requirement to document a number of Software Development Life Cycle ("SDLC") tasks in accordance with departmental, company and regulatory policies.

6.1.2 R package SDLC and platform architecture

Before the design, development and/or testing of any code was initiated, two further key platform architectural design decisions were made: first, GitLab was used for source control, continuous integration, documentation, vignettes, readme files and also as part of the full deployment process; and secondly, R Studio Server Pro was used for development and testing of code, a Docker Image with a physical R server on AWS.

6.1.3 Further R package design: function layers

With a large number of complex R scripts and source papers another design choice (primarily, to make the code easier to use and easier to maintain) was grouping the code into four layers using R’s S3 OO approach. The four programming layers are as follows:

Layer One: Highest Level: Main Exposed Application Programming Interface (API).

This level is exposed to the user and contains: S3 Classes (functions) that allow instantiation of the objects that contain the input data required and configuration; functions to predict events and recruitment; plotting and printing functionality; and key getter functions.

Layer Two: Second Level Functions.

This level is not exposed to the user and is simply used to dispatch to the third level functions based on the S3 configuration objects instantiated in Layer One.

Layer Three: Third Level Functions.

This level is not exposed to the user and contains the main set of controller code and does all of the hard work of the package.

Layer Four: Lowest Level Functions.

This level is not exposed to the user and contains a large number of complex R scripts/algorithms that have been developed and tested using Monte Carlo Simulation, as detailed in the previous section.

6.2 R package input data required

The following set of input data is required by the EventPrediction package to predict events and recruitment (if recruitment is ongoing), with each set of data instantiated using R’s S3 approach (as detailed in the previous sections).

6.2.1 Event data

analysis_time_days  censor_flag  drop_out_flag  randomisation_date
28 0 0 YYYY-MM-DD
33 0 0 YYYY-MM-DD
87 1 0 YYYY-MM-DD
42 0 0 YYYY-MM-DD
77 1 1 YYYY-MM-DD

This data is in accordance with how the key stakeholders produce their data, it is transformed into the values as described in the previous sections, such that:

  • analysis_time_days is the number of days from randomisation to either the event date or censoring date (i.e. the dropout date for subjects that have dropped out or the date used to censor at the cut off if a subject has not dropped out).

  • censor_flag == 0 represents group A, a subject experienced event A.

  • censor_flag == 1 & drop_out_flag == 0 represents group O, a subject did not experience an event nor dropout.

  • drop_out_flag == 1 represents group L, a subject dropped out before the interim time.

6.2.2 Site recruitment data

study_centre_id  centre_actual_enrol  centre_recruitment_window_days
xx001 0 140
xx002 1 224
xx003 2 238
xx004 1 221
xx005 0 201
  • centre_actual_enrol represents the number of subjects recruited at the unique centre ID.

  • centre_recruitment_window_days represents the actual duration of recruitment at the unique centre ID (not including any screening period). The centre is active only during this interval .

6.2.3 New Sites

A vector of days for new centres to be initiated {

}, e.g. .

6.2.4 Configuration

The following key pieces of information are accepted by the EventPrediction package (with appropriate defaults) which are used to select the appropriate algorithms and to provide key modelling values:

  • distributions_to_use: A list detailing the distributions to model the dropouts and events: e.g. list(events = "Exponential", drop_outs = "Exponential")

  • target_number_of_events: The target number of events for the analysis to be predicted

  • sample_size: The number of patients planned to recruit

  • confidence_level: The confidence probability for the upper and lower bounds

7 R package and implementation in a clinical trial

7.1 Introduction

In order to help the key stakeholders with the operational planning of a clinical trial and to test the quality of the prediction, the EventPrediction package was used on several historic studies. The following is one such case study in a historical oncology clinical trial, using the data at a given interim time when recruitment had not completed. The task was to predict the future recruitment and event counts with bounds and compare the results with the real trajectory of the recruitment and the events that have already occurred in the past.

The event and centre data was provided in accordance with the package API’s (as detailed in the previous section) along with a target number of events of 250 and patients sample size of 405. At the interim cut-off time the data for the study had the following recruitment and event status:

  • 152 Events (i.e. censor_flag == 1)

  • 155 At Risk (i.e. censor_flag == 1 and drop_out_flag == 0)

  • 13 Drop Outs( i.e. drop_out_flag == 1)

  • 85 patients left to recruit.

7.1.1 Key predictions

Given the above input and implementing the developed model, the EventPrediction package predicted:

Recruitment: Predicted number of days until target number of patients is reached with 90% bounds, (mean, lower bound, upper bound): 151, 120, 191.

The estimated parameters of a PG model are: , and the prediction is constructed according to (15) where it was used some schedule of closing centres.

Events: Predicted number of days until target number of events reached, with 90% bounds, (mean, lower bound, upper bound):
exponential model: 227, 181, 322
Weibull model: 241, 188, 423

Figure 3: Prediction of the remaining recruitment against time from cut-off (in days). Real trajectory of patient recruitment is shown by black line. Mean prediction and 90% bounds are shown by the blue dashed and red dotted lines.

7.1.2 Plots and parameter estimates

Further, the EventPrediction package produced the following three key plots, along with key parameter estimates, for the key stakeholders to consume:
1) prediction of the remaining recruitment
2) prediction of the remaining number of events using exponential model
3) prediction of the remaining number of events using Weibull model

Figure 4: Prediction of the remaining number of events against time from cut-off (in days). Real trajectory of events is shown by black line. Exponential model, mean prediction and bounds depicted by the blue dashed and red dotted lines.
Figure 5: Prediction of the remaining number of events against time from cut-off (in days). Real trajectory of events is shown by black line. Weibull model, mean prediction and bounds depicted by the blue dashed and red dotted lines.

Fig 3 shows very good fit of the predictive area of recruitment where the real trajectory of the historical recruitment falls into the predictive area.

As one can see from Fig 4 and Fig 5, the predictions for both types of models, exponential and Weibull, are rather close, with the following estimated parameters for each:
Exponential model: , and .
Weibull model: , , , and .

The actual number of days when the target number of patients was reached in this trial was 152 and the actual number of days when the planned number of events occurred is 244. As seen in the figures, the predictions were indeed very close to the actuals.

7.2 Predicting events when recruitment complete

We also considered the same case study as above, but at a later interim time when recruitment had completed. Therefore the task here was to predict the future event counts only.

As detailed in the previous section, this study has a target number of events of 250 and patients sample size of 405. At the interim cut-off time the data for the study had the following recruitment and event status:

  • 220 Events (i.e. censor_flag == 1)

  • 163 At Risk (i.e. censor_flag == 1 and drop_out_flag == 0)

  • 22 Drop Outs( i.e. drop_out_flag == 1)

7.2.1 Key predictions

Given the above input and implementing the developed model, the EventPrediction package predicted:

Events: Predicted number of days until target number of events reached, with 90% bounds, (mean, lower bound, upper bound):
exponential model: 84, 60, 118
Weibull model: 83, 59, 116

7.2.2 Plots and parameter estimates

Further, the EventPrediction package produced the following two key plots, along with key parameter estimates, for the key stakeholders to consume:
1) prediction of the remaining number of events using exponential model
2) prediction of the remaining number of events using Weibull model

Figure 6: Prediction of the remaining number of events against time from cut-off (in days). Real trajectory of events is shown by black line. Exponential model, mean prediction and bounds depicted by the blue dashed and red dotted lines.
Figure 7: Prediction of the remaining number of events against time from cut-off (in days). Real trajectory of events is shown by black line. Weibull model, mean prediction and bounds depicted by the blue dashed and red dotted lines.

As one can see from Fig 6 and Fig 7, the predictions for both types of models, exponential and Weibull, are rather close, with the following estimated parameters for each:
Exponential model: , and .
Weibull model: , , , and .

The actual number of days when the planned number of events occurred is 91. As seen in the figures, the predictions were indeed very close to the actuals.

8 Fitting models to real data

8.1 Kaplan-Meier plots

To assess the model fit, we looked at the Kaplan-Meier (KM) curve for each interim dataset and compared these to the predicted survival functions for exponential and Weibull distributions. This provides a visualisation for model fit: the best fit model for the interim data will be the model which distribution best maps the KM curve. This is for the occurrence of events and so will only inform of the best distribution for modelling events . However, similar curves can be created to test the fit of dropout distribution.

The survival functions for exponential and Weibull models with cure are calculated using their respective estimated parameters, with formulae:

Exponential:
Weibull:

In Fig 8, one can see that the predicted survival functions for exponential and Weibull models are very close and both map the KM curve well.

Figure 8:

Kaplan-Meier plot showing survival function using interim data with associated confidence intervals (blue lines), alongside survival functions for exponential (violet solid line) and Weibull (green dashed line) cure models.

8.2 AIC and BIC criteria

We also used Akaike information criterion (AIC) and Bayesian information criterion (BIC) to test which model best fits the real dataset. These are calculated by:

AIC

BIC

where = number of parameters in the model, = total number of patients at interim time, = the log-likelihood from parameter estimation.

For the stopped recruitment case:

Model AIC BIC
Exponential model, no cure 3333.6 3341.6
Exponential cure model 3319.0 3327.1
Weibull (A) and exponential (L) cure model 3321.0 3333.0
Exponential (A) and Weibull (L) cure model 3283.1 3299.1
Weibull cure model 3285.0 3305.1

For the ongoing recruitment case:

Model AIC BIC
Exponential model, no cure 2222.9 2230.4
Exponential cure model 2210.9 2218.4
Weibull (A) and exponential (L) cure model 2209.3 2220.6
Exponential (A) and Weibull (L) cure model 2183.7 2198.7
Weibull cure model 2182.1 2200.9

The lower the AIC or BIC, the better the model fit. From the two tables above, one can see that the cure models show an improvement on the ”no cure” model. The exponential (A) and Weibull (L) cure model is the best fit for the stopped recruitment dataset. By small margins, the values of criteria for ongoing recruitment suggest either the Weibull cure model or the exponential (A) and Weibull (L) cure model is best suited for the data.

Conclusions

We have developed a new analytic approach for the prediction of event counts in event-driven trials when recruitment is complete or ongoing. We use the exponential and Weibull models and account for patient dropout and opportunity of cure. Not only can we predict the future occurrence of events, but we can also predict any remaining recruitment, with mean and bounds, using the Poisson-gamma recruitment model. The developed results can be easily extended to the combined cure models using the combination of exponential and Weibull distributions for the time to event and the time to dropout and also to log-normal with cure model.

Using these novel advanced models and with access to real subject level and centre level data we are now able to address key business use cases for a number of key stakeholders in order to better forecast the operational design of event-driven clinical trials. Furthermore, by centralising an exposed R Package EventPrediction, utilising good software engineering principles, each of our key stakeholders have access to the package and can obtain plots, parameter estimates and predictions with bounds, without contacting the mathematical modellers nor the R package developer.

We have many opportunities for future improvements to the mathematical modelling and EventPrediction package. One of the major priorities is to evaluate the predictions against real clinical trial operational data to ensure the existing and any future models are as accurate as we found in testing on historical data. We are also looking to incorporate other statistical distributions for modelling time to both the main event and dropout.

References

  • [1] V. Anisimov and V. Fedorov. Modeling, prediction and adaptive adjustment of recruitment in multicentre trials, Statistics in Medicine, 26, 27, 4958–4975, 2007.
  • [2] V. Anisimov, D. Downing and V. Fedorov. Recruitment in multicentre trials: prediction and adjustment, mODa 8 - Advances in Model-Oriented Design and Analysis, 1–8, 2007.
  • [3] V. Anisimov. Using mixed Poisson models in patient recruitment, Proc. of the World Congress on Engineering, II, 1046–1049, 2008.
  • [4] V. Anisimov. Statistical modeling of clinical trials (recruitment and randomization), Communications in Statistics - Theory and Methods, 40, 19-20, 3684–3699, 2011.
  • [5] V. Anisimov. Predictive event modelling in multicentre clinical trials with waiting time to response, Pharmaceutical Statistics, 10, 6, 517–522, 2011.
  • [6] V. Anisimov. Predictive hierarchic modelling of operational characteristics in clinical trials, Communications in Statistics - Simulation and Computation, 45, 5, 1477–1488, 2016.
  • [7] V. Anisimov. Discussion on the paper "Real-time prediction of clinical trial enrollment and event counts: a review" by D.F. Heitjan et al. Contemporary Clinical Trials, 40, 7–10, 2016.
  • [8] V. Anisimov. Modern analytic techniques for predictive modelling of clinical trial operations, Quantitative Methods in Pharmaceutical Research and Development: Concepts and Applications, Springer International Publ., 361–408, 2020.
  • [9] V. Anisimov and M. Austin. Centralized statistical monitoring of clinical trial enrollment performance, Communications in Statistics - Case Studies and Data Analysis, 6, 4, 2020, 392–410.
  • [10] E. Bagiella and D.F. Heitjan. Predicting analysis times in randomized clinical trials, Statistics in Medicine, 20, 2055–2063, 2001.
  • [11] K.D. Barnard, L. Dent and A. Cook. A systematic review of models to predict recruitment to multicentre clinical trials, BMC Medical Research Methodology, 10, 63, 2010.
  • [12] J.M. Bernardo and A.F.M. Smith. Bayesian Theory, John Wiley & Sons: Hoboken, NJ, USA, 2004.
  • [13] T.T. Chen. Predicting analysis times in randomized clinical trials with cancer immunotherapy, BMC medical research methodology, 16(1), 1–10, 2016.
  • [14] E. Gkioni, R. Riusd, S. Dodda and C. Gamblea. A systematic review describes models for recruitment prediction at the design stage of a clinical trial, Journal of Clinical Epidemiology, 115:141–149, 2019.
  • [15] D.F. Heitjan, Z. Ge and G.S. Ying. Real-time prediction of clinical trial enrollment and event counts: a review, Contemporary Clinical Trials, 45, part A, 26–33, 2015.
  • [16] R.E. Carter, S.C. Sonne and K.T. Brady. Practical considerations for estimating clinical trial accrual periods: Application to a multi-center effectiveness study, BMC Medical Research Methodology 5:11–15, 2005.
  • [17] S. Senn. Statistical Issues in Drug Development. Wiley: Chichester, 1997.
  • [18] S. Senn. Some controversies in planning and analysis multi-center trials. Statistics in Medicine, 17, 1753–1756, 1998.