How to deal with the Poisson-gamma model to forecast patients' recruitment in clinical trials when there are pauses in recruitment dynamic?

Recruiting patients is a crucial step of a clinical trial. Estimation of the trial duration is a question of paramount interest. Most techniques are based on deterministic models and various ad hoc methods neglecting the variability in the recruitment process. To overpass this difficulty the so-called Poisson-gamma model has been introduced involving, for each centre, a recruitment process modelled by a Poisson process whose rate is assumed constant in time and gamma-distributed. The relevancy of this model has been widely investigated. In practice, rates are rarely constant in time, there are breaks in recruitment (for instance week-ends or holidays). Such information can be collected and included in a model considering piecewise constant rate functions yielding to an inhomogeneous Cox model. The estimation of the trial duration is much more difficult. Three strategies of computation of the expected trial duration are proposed considering all the breaks, considering only large breaks and without considering breaks. The bias of these estimations procedure are assessed by means of simulation studies considering three scenarios of breaks simulation. These strategies yield to estimations with a very small bias. Moreover, the strategy with the best performances in terms of prediction and with the smallest bias is the one which does not take into account of breaks. This result is important as, in practice, collecting breaks data is pretty hard to manage.

Introduction

In order to get marketing authorization, a new product has to succeed in clinical trials. A clinical trial is based on statistical considerations in order to show the product efficacy, taking into account the variability of the environment. It is a well known fact that the power of this test is linked to the number of patients we deal with. If an inadequate number of enrolled patients is used, then the study may fail to reject the null hypothesis due to lack of power. In many cases the goal number of patients to include is thus a fixed parameter of the trial. There has been much effort in computing the sample size for clinical trials. Its computation is now standard in trial protocols (see ICH E9 guidance) and mandatory for most of the publications (see Consort Group works Schulz et al. [20]). On the other side relatively little attention is focused on improving the prediction of the recruitment process. Indeed, till now the most of techniques used by pharmaceutical companies are based on deterministic models and various ad hoc techniques. Rojavin [19] says “Patient recruitment and retention remains until now more of an art rather than a science”.

The problem of predicting patients recruitment and evaluating the recruitment time in clinical trials is of paramount interest for planning trials because of scientific concern, economic and ethical reasons.

Ethical concern, because it is not satisfactory to continue a study in vain.

Economical concern, a clinical trial is an expensive study in itself and, as the duration of the trials is included in the duration of the exclusive right to exploit the drug (20 years) [22], a delay generates an enormous loss of income. Moreover, an improvement of the planning and monitoring of a trial may reduce costs.

Scientific concern, because new drugs are increasingly developed and approved by regulatory agencies. When accrual rates are too low, there may be new information available during the enrolment period such as the results of other trials or a change in the understanding of the underlying biology.

For these reasons, stopping or continuing a trial is a decision with huge consequences and it will be useful to have some objective tools based on scientific criteria to take it.

Few authors have considered the problem of patients' recruitment. The reader can refer to Barnard et al. [11] for a systematic review of the existing models for recruitment and Anisimov [7], Bakhshi et al. [10], Gajewski et al. [13], Heitjan et al. [14] for a presentation and/or a discussion of methods. As far as we know, the pioneer work is the one of Morgan [18] where an estimation of the total study duration is proposed as a function of inclusion duration and based on data from previous clinical trials. Let us cite Lee [15] for a model of recruitment by Poisson processes. In Senn [21] a model for multicentric trials based on Poisson process is introduced. Poisson process appears as a natural assumption in literature [21]. However, Poisson processes depends on only one parameter, which is the rate of enrolment in our setting and Carter [12] have noticed that the use of the historic mean is a too simple model and highlighted the necessity to take into consideration the variability of the rates.

Inspired by queueing theory, Anisimov and Fedorov [9] proposed to use a doubly stochastic Poisson process to take into consideration the variation in recruitment rates between different centres. This model, called as a Poisson-gamma model, assumes that the patients arrive at different centres according to Poisson processes with the rates viewed as independent gamma distributed random variables. In Anisimov and Fedorov [9] the procedure of parameters estimation from data collected up to an interim time using empirical Bayesian technique have been suggested. The model has been validated using data from a large number of real trials [3]. The Poisson-gamma model was developed further for predicting recruitment process at initial and interim stages [1], to account for the situations when the centres opening dates may not be known and assumed to be uniformly distributed in some intervals [2], [17] or when some centres can be closed or opened in the future [6]. Sensitivity analysis to errors in parameters' estimation has been investigated in Mijoule et al. [17]. Poisson-gamma model allows to develop techniques for analysing the effects of unstratified and centre-stratified randomization [4], for predictive event modelling [5], for predict randomization process [6] and for modelling clinical trials' cost [16].

The Poisson-gamma approach on the modelling of patients' recruitment is more and more popular. However, the assumption on the rates, which are assumed to be constant in time, has to be discussed. Indeed, in practice, there are breaks in recruitment process. A break is defined as a period during which a centre does not recruit any patient (holidays, week-end,..). This information is observed and can be collected but, in practice, it is a huge and complicated process. The Poisson-gamma model can be enriched considering piecewise constant rate function yielding to inhomogeneous Cox models. In this case the estimation of the trial duration is much more difficult. Three strategies of estimation of recruitment duration are proposed in this paper: one taking account of breaks, one not taking account of breaks and a third taking into account of only large breaks. These estimations are biased. The bias is assessed by means of simulation study and the strategies are compared in terms of bias in estimation and in terms of predictive performances. Finally, these investigations allow us to deal with a question of paramount interest: Is it really useful to enrich the model in order to take into account breaks in recruitments?

The paper is organized as follows. Section 2 describes the Poisson-gamma model. The procedure of parameters estimation and the computation of the expected duration of the trial are given. Section 3 presents strategies which can be used to take into account of breaks in the recruitment model. In Section 4 authors explain the data generation procedure which allows to perform the simulation studies. These studies investigate the bias and the predictive performances of the strategies introduced in Section 3 together with the comparison of these strategies. The results are presented in Section 5. Finally the paper ends with a concluding Section 6.

The Poisson-gamma model without breaks in recruitment

Notations

Consider a multicentric clinical trial where C centres are involved to recruit n patients. Denote u the opening date of the i-th centre which is assumed to be observed. The recruitment process of centre i is denoted and is modelled by a Poisson process with intensity λ. The global inclusion process is . The parameter of paramount interest is the stopping time:

In the sequel we consider, for any i, λ as a random variable which is gamma-distributed with parameters (α,β) whose probability density function is:where .

Estimation

In most settings, parameter (α,β) is unknown. To estimate this parameter, an empirical Bayesian strategy is used. Fix an interim time of analysis t1. Data collected on [0,t1] are used to calibrate the model. For any centre i, denote the duration of activity up to t1 of centre i and the number of patients recruited by centre i up to t1. Notice that {(τ1,,k1,),i = 1, …, C} are observed data.

(Anisimov et al. [8]). Maximum likelihood estimation of the parameter (α,β) is obtained by maximization of the function:

Prediction

If centre i is initiated at time 0 (u = 0), the so-called forward rate knowing is Gamma-distributed with parameters . If centre i is not initiated at time 0, the forward instantaneous intensity is time dependent and expressed as:where is a predictive posterior rate in center i with distribution . The prediction process is now written as:

(Anisimov et al. [8]). Given, , the predictive recruitment process N expressed, for any t > t1 as:and is a Cox process whose predictive rate is -distributed starting at time t1.

Expected duration

Consider the number of patients remaining to recruit after t1 and the remaining inclusion time. Denote .

(Anisimov [6]). • Assume that all centres are initiated at the same time (u = u > 0, for any i) and denote . Then is -distributed where PVI (n, a, b) denotes the Pearson VI distribution whose probability density function is:where is the beta function.

In practice, the τ1,'s may be different. The distribution of can be approximated for large n by a distribution withwhere

Recall that for X a random variable P(n, a, b)-distributed with a > 2, the expectation and the variance are:

As a consequence, by Theorem 2 the expression of the expected duration of the trial is:

The Poisson-gamma model with breaks in recruitment

Assume that the recruitment process for centre i stops at some time denoted b for periods denoted d. As in previous section, fix an interim time of analysis t1. Data collected on [0,t1] will be used to calibrate the model. For centre i, the data collected are:

the number of recruited patients up to t1 denoted k1,,

the number of breaks up to t1:

(b1,,j = 1, …, ,j1,) (resp. (d1,,j = 1, …, ,j1,)) the breaks times (resp. durations) up to t1.

The duration of activity for centre i up to t1 is thus . The recruitment dynamic for centre i, still denoted as N, is a non-homogeneous Cox process with intensity governed by λ which is Γ(α,β) distributed. Indeed, the instantaneous intensity for centre i is time dependent and expressed as:with .

The following theorem insures that (α,β) are estimated following the same strategy as in Section 2. Only the definition of τ1, differ.

Maximum likelihood estimation of the parameter (α, β) is obtained by maximization of the function:

When considering breaks, the difficulty comes from the Cox process N modelling the recruitment which is non-homogeneous because of potential breaks in the dynamic. To overpass this difficulty, consider a homogeneous Cox process starting at t1 with intensity . Recall the remaining recruitment process defined in Theorem 2. Processes and have the same intensities but the first one allows breaks while the second one does not (see Fig. 1).

Fig. 1

Processes involved in the prediction of T.

Consider the “true” duration of the clinical trial and (respectively ) the remaining time after t1 considering there are (respectively there are no) breaks in the recruitment. Obviously . The duration denoted by (B) on Fig. 1 is thus interpreted as the cumulated breaks duration.

The expected duration cannot be estimated directly but, as , an estimation can be proposed. Indeed, on the one hand, Theorem 3 gives us an estimation of since is a homogeneous Cox process and, on the other hand, assuming that breaks behave after t1 as before t1 and that the cumulated breaks duration is proportional to the duration of the follow up, can be estimated bysince is the cumulated breaks duration on [0,t1] and is the cumulated duration of activity of the centres on [0,t1]. Replacing , which cannot be computed, by yield to the following estimation of .

The simulation study of Section 4 aims to quantify the bias when is estimated by .

Another strategy

In practice, it may be useful to consider only sufficiently large breaks (greater than a predefined value denoted dmax). Consider t1 an interim time of analysis of the recruitment process. For the center i, the data collected on [0,t1] are:

the number of recruited patients up to t1 denoted k1,,

the number of breaks which duration is greater than dmax up to t1

the breaks times and durations up to t1 greater than dmax denoted (b1,,j∈J) and (d1,,j∈J), with

The duration of activity up to t1 expressed as . The rate of recruitment for centre i is governed by λ which is Γ(α,β)-distributed. The instantaneous intensities are still time dependent and expressed, for any centre i, as:with . Considering these new notations and definitions, the estimation of the parameters of the model, the predictive process and expected duration of the trial can be computed as explained in sections 3.1, 3.2, 3.3.

Data generation procedure

Consider a multicentric trial involving C = 60 centres. We aim to recruit n = 720 patients in 365 days.

The whole numerical values introduced in this section to perform the data generation procedure are chosen in such a way that the recruitment dynamic is consistent with what is observed on real trial. To do so, authors have taken the parameters of the real case study of Anisimov and Fedorov [9].

In order to investigate different approaches of the breaks dynamic, different scenarios are proposed. Scenarios differ by the breaks generation procedure (times of breaks and durations of breaks). The scenarios are:

Scenario 1: Exponential generation. The instants and durations of breaks are generated according to exponential distributions. The breaks times are exponentially distributed with intensity 1/60 this means a break appears on average every 60 days. The breaks durations are exponentially distributed with intensity 1/14 this means the average break duration is 14 days.

Scenario 2: Multinomial generation. The instants of breaks are generated according to an exponential distribution with parameter 1/60 meanwhile the durations are generated according to a multinomial probability. Five levels of duration (2, 4, 8, 16 and 32 days) are involved. The corresponding probability vector is built in such a way that it ensures that the total durations for each level are the same, (the breaks of 2 days happen twice more than the one of 4 days).

Scenario 3: Deterministic generation. The instants of the breaks are the Saturdays and Sundays and a complete week every two months for holidays.

Recruitment dynamic is generated involving breaks themselves generated according to the scenarios defined above. The whole of the dynamic is known thus the true duration denoted T0 of the trial is known. The study consists in considering the data collected on [0,t1] and to make use of the results of sections 2, 3 to estimate the duration of the trial. The three strategies of estimation of the trial duration introduced in sections 2.4, 3.3, 3.4 are considered:

Strategy 1: not taking into account of breaks. The breaks times and durations are not collected. The parameters of the Poisson-gamma model for recruitment are estimated following results of Section 2. The expected duration is denoted by T1 and is computed by means of Theorem 3.

Strategy 2: Taking into account of all the breaks. The breaks times and durations are collected. This allow us to make use of the estimation of the trial duration as explained in Section 3.3. The expected duration is denoted T2.

Strategy 3: Taking into account of only large breaks. Only the breaks that have a duration greater than a predefined value of 14 days are collected and used as explained in Section 3.4. This strategy leads to the expected duration denoted T3.

For a sake of simplicity, all centres are initiated at t = 0 (u = 0 for all i). The data generation procedure splits in two steps:

Step 1: Generate R = 1000 recruitment processes according to scenario 1, 2 or 3. Consider, for any 1 ≤ r ≤ R, , where denotes the first time verifying N(t) = 720. The generation procedure is the following one:

Generate the breaks according to the scenario 1, 2, and 3 considered for a period of 730 days. The duration of 730 days is here arbitrary and chosen in order to be sure to catch the true duration of the trial.

Generate the rates according to a Γ(2, 60.8) distribution.

Consider the modified rate function as defined in Equation (1).

Generate the recruitment process up to 730 days.

Identify and shrink the recruitment process to .

Step 2: Consider an interim time at t1 = 182 days. For each simulation run r = 1, …, R and each strategy s = 1, 2, 3,

Estimate parameters of the gamma distribution applying Theorem 1 for s = 1 or Theorem 4 for s = 2 and s = 3 from data collected on [0,t1].

Compute the expected duration of the trial through the application of Theorem 2 for s = 1 or following strategies explained in Section 3.3 for s = 2 and Section 3.4 for s = 3.

The performances of the model at interim time t1 are measured by means of the absolute error defined by:

The value chosen for the minimal break duration in analysis strategy 3 (14 days) is not coherent with the values chosen in scenario 3 (2 and 7 days). There is thus no difference between the first and the third scenario in this setting.

Results and discussion

Results on the bias

The first point investigated by the simulation study is the assessment of the bias embedded by the strategy of estimation of the expected duration for each scenario. The results are collected in Table 1. For each scenario, the mean duration (over the simulation runs) of the simulated recruitment dynamic together with its 95% confidence interval (which corresponds to the 25-th and 975-th values of the sorted samples) are identified. For Strategy p (p = 1,2,3), the average expected trial duration and its 95% confidence interval are computed. The bias is assessed by the average (over r) of the absolute errors made between the expected durations and the true value of the trial duration denoted E0,. Table 1 presents the value of E0,, its 95% confidence interval and the parameter denoted S0, which measure the proportion of over-estimation defined as the proportion of runs for which the expected duration computed by means of Strategy p is greater than the true duration .

Table 1

Average duration of the trial (in days) as a function of p, the strategy, together with its 95% confidence Interval. Absolute error between expected duration and true duration (E0,) as a function of the strategy together with its 95% confidence interval and the proportion of over-estimation (S0,).

		Real duration	Strategy 1	Strategy 2	Strategy 3
Scenario 1	Mean	384.37	381.92	374.66	377.73
	CI (95%)	[351,415]	[340,431]	[336,421]	[337,426]
	E0,p	–	13.81	15.42	14.41
	CI (95%)	–	[2,13]	[1,40]	[0,38]
	S0,p	–	0.41	0.27	0.32
Scenario 2	Mean	337.92	337.8	336.41	337.16
	CI (95%)	[303,373]	[304,374]	[303,372]	[353,406]
	E0,p	–	9.78	9.88	9.82
	CI (95%)	–	[0,29]	[1,29]	[0,29]
	S0,p	–	0.44	0.40	0.42
Scenario 3	Mean	369.17	367.95	356.94	–
	CI (95%)	[343,398]	[330,412]	[321,398]	–
	E0,p	–	11.04	14.88	–
	CI (95%)	–	[1,32]	[1,38]	–
	S0,p	–	0.44	0.16	–

These results are enriched by Fig. 2 (resp. 3 and 4) which refers to Scenario 1 (resp. 2, 3). Each figure is composed of the plots of the histograms of the (R = 1000) expected durations for each strategy together with the plot of the empirical density curves for each strategy completed by the one of the true durations. Keep in mind that, in accordance with Remark 2, Strategy 3 is not investigated for Scenario 3.

Fig. 2

Scenario 1. Histograms of expected duration for all strategies and comparison of the densities (solid line: T0, dashed line: T1, dotted line: T2, longdash line: T3).

First, the computation strategies are efficient. Indeed, on Table 1, the mean value of the real duration is close to the mean durations estimated by the different strategies. This is confirmed by the values of E0,. Moreover, the width of the confidence intervals of E0, is not large whatever the setting.

Second, for any scenario, the three strategies under-evaluate the trial duration. This is not surprising since the estimation (3) of involves instead of and . The values of S0, on Table 1 and Fig. 2, Fig. 3, Fig. 4 confirmed this results. Indeed, the histograms are uni-modal, not asymmetric and shift to the left comparing with the one of the real duration. The strategies are thus moderately biased.

Fig. 3

Scenario 2. Histograms of expected duration for all strategies and comparison of the densities (solid line: T0, dashed line: T1, dotted line: T2, longdash line: T3).

Fig. 4

Scenario 3. Histograms of expected duration for all strategies and comparison of the densities (solid line: T0, dashed line: T1, dotted line: T2).

Comparison of the strategies

Table 1 informs us that the error generated by Strategy 1 is minimal whatever the setting and that the confidence intervals for that strategy are the smallest. Strategy 1 is the less biased strategy with values of S0, close to 0.5. Fig. 2, Fig. 3, Fig. 4 highlight this result.

To specify the comparisons of the three strategies, results analogous to those of Table 1 are collected in Table 2 for Scenario 1, Table 3 for Scenario 2 and Table 4 for Scenario 3. Each block of the table corresponds to a comparison of two strategies. These comparisons are given in terms of absolute error together with its 95% confidence interval (which corresponds to the 25-th and 975-th values of the sorted samples) and the proportion of over-estimations, this means the proportion of runs for which the expected duration computed by strategy s′ is greater than the one computed by strategy s.

Table 2

Results for Scenario 1. Comparison of the strategies of analysis in terms of absolute error together with its 95% confidence interval and the proportion of over-estimation.

		Strategy 1	Strategy 2
Strategy 2	E2,p	7.27	–
	CI (95%)	[2,13]	–
	S2,p	0.99	–
Strategy 3	E3,p	4.25	3.07
	CI (95%)	[0,10]	[1,5]
	S3,p	0.94	0.01

Table 3

Results for scenario 2. Comparison of the strategies of analysis in terms of absolute error together with its 95% confidence interval and the proportion of over-estimation.

		Strategy 1	Strategy 2
Strategy 2	E2,p	1.44	–
	CI (95%)	[0,3]	–
	S2,p	0.82	–
Strategy 3	E3,p	0.82	0.83
	CI (95%)	[0,2]	[0,2]
	S3,p	0.55	0.04

Table 4

Results for scenario 3. Comparison of the strategies of analysis in terms of absolute error together with its 95% confidence interval and the proportion of over-estimation.

		Strategy 1
Strategy 2	E2,1	11.01
	CI (95%)	[9,14]
	S2,1	1

Whatever the scenario, Table 2, Table 3, Table 4 indicate that the results obtained by the different strategies are more or less the same. Notice that, for Scenario 1, the largest difference is around 7 days for a study of 365 days which is very small. The confidence intervals are not large and the width is more or less the same whatever the strategy of computation. Finally the proportion of over-estimation S indicates that results from Strategy 1 over-estimate those from Strategy 3 which over-estimate those from Strategy 2.

Predictive performance of the strategies

On contrary to bias analysis, predictive behaviour analysis of a strategy is an individual concept. To investigate the predictive performances of a strategy, Fig. 5 related to Scenario 1, Fig. 6 to Scenario 2 and Fig. 7 to Scenario 3 are considered. Each figure splits in sub-plots which are the plots of the 1000 values of the real duration in abscissa and the expected duration in ordinate for each strategy of computation. The plots are enriched with the regression line and the straight line y = x.

Fig. 5

Scenario 1. Linear Regression of the expected duration as a function of the real duration for each strategy. Solid line: y = x line, dashed line: regression line.

Fig. 6

Scenario 2. Linear Regression of the expected duration as a function of the real duration for each strategy. Solid line: y = x line, dashed line: regression line.

Fig. 7

Scenario 3. Linear Regression of the expected duration as a function of the real duration for each strategy. Solid line: y = x line, dashed line: regression line.

Whatever the scenario and whatever the strategy, the regression analysis yields to the same conclusions: the quality of the regression model is good (quality of the fitting, normality, homoscedasticity and autocorrelation of the residuals values). It is thus possible to compare the slope to 1 and the intercept to 0 by means of Wald's tests and allow us to conclude that the differences observed on the different plots are not significant. The predictive performances are thus very good.

Differences between strategies and between scenarios are not significant but exist. These differences are connected to results on bias and illustrate the over-estimation issue already discussed. For instance, for Scenario 1, Fig. 5 shows that Strategy 1 is closest to the line y = x and thus appears to be the best strategy (the slope are 0.94 for the first strategy, 0.87 for the second and 0.89 for the third). The under-estimation of the strategies regardless to the real duration is observable, the regression line is under the line y = x. It is easily seen on the plots that the difference between the regression line and the y = x line increases with the true duration. The bias in the estimation is thus time dependent unlike what is observed for Scenario 3 with Strategy 2 (subplot on the right of Fig. 7).

Conclusions and recommendations

To conclude, results of Section 5 allow us to claim that the procedures proposed in this paper to estimate the expected duration of the recruitment accounting of breaks in the recruitment process are moderately biased and have very good predictive performances. Notice that the assumption that the breaks behave after the interim time as they behave before the interim time is reasonable.

The main result states that the strategy with the best results is the one ignoring breaks. This is not surprising since, for this strategy, estimations (2) and (3) are the same and the induced bias is null. That strategy applies directly the Poisson-gamma model on the data collected. The consequence of the breaks is balanced by the estimation of the parameters of the gamma distribution.

Thus, finally, as a recommendation one suggests to not worry with breaks and to make use of a standard Poisson-gamma model. Indeed, to collect breaks data during recruitment process is a difficult process in practice and the impact on the recruitment process modelling is really moderate.