Simulate exact continuous-time enrollment — enrollment • goldilocks

Simulates enrollment times from a Poisson process with a piecewise-constant rate.

Usage

enrollment(lambda = 1, N_total, lambda_time = NULL)

Arguments

lambda: A numeric vector of finite, positive enrollment rates per unit of calendar time. Supply one rate for each interval defined by lambda_time, so length(lambda) must equal length(lambda_time) + 1. The default is 1.
N_total: A required positive integer giving the total sample size.
lambda_time: NULL (the default), or a numeric vector of finite, positive, strictly increasing interior times at which the enrollment rate changes. The initial boundary at time zero is implicit and must not be supplied. Use NULL for a constant enrollment rate.

Value

A non-decreasing numeric vector of N_total continuous enrollment times, measured from first patient in and expressed in the same time unit used for lambda_time. The first value is zero.

Details

Major behavior change from goldilocks 0.5.0 and earlier. Versions through 0.5.0 generated Poisson counts in unit-time bins. enrollment() returned rebased integer bin times, after which sim_comp_data() added independent uniform jitter and sorted the result. From version 0.6.0, enrollment times are drawn directly from the exact continuous-time piecewise-constant Poisson process. Consequently:

seeded simulations do not reproduce enrollment or downstream trial results obtained with version 0.5.0 or earlier;
enrollment times and operating-characteristic estimates can change, particularly when rates are low or a rate change is not an integer time;
the lambda_time argument now contains internal change times only: change lambda_time = 0 to lambda_time = NULL for a constant rate, and change lambda_time = c(0, t1, t2) to lambda_time = c(t1, t2) for a piecewise rate.

enrollment() treats time zero as first patient in, matching the time origin used throughout goldilocks. The first returned enrollment time is therefore fixed at zero. The remaining N_total - 1 arrivals form a continuous-time non-homogeneous Poisson process whose rate is constant between the supplied internal knots. Thus, lambda is measured in enrollments per unit of lambda_time; for example, when time is measured in months, lambda = 5 means five enrollments per month on average.

Write the internal knots as $0 < \tau_1 < \cdots < \tau_K$, with $\tau_0 = 0$ implicit. The enrollment intensity is

$$ \lambda(t) = \lambda_j, \qquad \tau_{j-1} \le t < \tau_j, $$

for $j = 1, \ldots, K$, and $\lambda(t) = \lambda_{K+1}$ after the final knot. The final rate therefore continues for as long as needed to reach N_total; there is no finite accrual horizon in this function. The value of the intensity at an isolated knot does not change the Poisson process. When assigning a realized enrollment time to an interval, goldilocks follows the survival counting-process convention: intervals are open on the left and closed on the right, so an arrival exactly at $\tau_j$ belongs to the interval ending at $\tau_j$. The first patient at time zero is a fixed origin and is handled separately.

Arrivals are generated exactly by the time-rescaling theorem. If $E_2, \ldots, E_N$ are independent unit-rate exponential variables and $S_i = \sum_{k=2}^i E_k$, then the enrollment times after first patient in are

$$T_1 = 0, \qquad T_i = \Lambda^{-1}(S_i),$$

where $\Lambda(t) = \int_0^t \lambda(u)\,du$ is the cumulative enrollment intensity. This construction gives independent Poisson counts on disjoint calendar intervals, with expected count $\int_a^b \lambda(u)\,du$ over $(a,b]$. For a constant rate, successive enrollment gaps are independent Exponential(lambda) variables.

The rate-change times are measured from first patient in, not from an earlier site-opening or trial-activation date. If operational delays before first patient in are important, they must be modelled separately before using the returned relative times.

For example, lambda = c(0.3, 0.7, 0.9, 1.2) with lambda_time = c(5, 10, 15) specifies average enrollment rates of 0.3 over positive times in $(0,5]$, 0.7 over $(5,10]$, 0.9 over $(10,15]$, and 1.2 after time 15. Fractional knots such as lambda_time = 2.5 are handled exactly; no unit-time binning or post-hoc jitter is used.

Examples

# Constant enrollment: first patient at zero, then exponential gaps.
enrollment(lambda = 0.7, N_total = 10)
#>  [1]  0.000000  3.387777  4.343014  4.630902  6.979419 14.372719 15.590170
#>  [8] 16.255547 17.034180 17.068086

# Three internal rate changes define four enrollment intervals.
enrollment(
  lambda = c(0.3, 0.7, 0.9, 1.2),
  N_total = 50,
  lambda_time = c(5, 10, 15)
)
#>  [1]  0.000000  3.248967  6.110215  8.898290 10.595670 12.926809 14.313422
#>  [8] 14.536687 16.280805 16.310522 16.636895 16.951155 19.969558 20.313615
#> [15] 20.417782 20.718054 21.351142 21.800600 22.025401 22.820496 24.306914
#> [22] 24.359537 25.055098 26.824180 27.749159 28.291158 28.948175 29.064916
#> [29] 29.928269 30.088707 31.687406 32.000755 32.673702 33.267533 33.856746
#> [36] 34.857605 34.957459 35.179734 35.347463 35.962991 36.674114 37.483202
#> [43] 38.672304 38.758390 40.229620 43.388860 43.912279 45.002938 46.399839
#> [50] 46.412589

# Fractional change times are supported exactly.
enrollment(
  lambda = c(0.25, 1),
  N_total = 20,
  lambda_time = 2.5
)
#>  [1]  0.000000  1.282874  1.377205  2.890406  3.307968  3.315925  3.564742
#>  [8]  3.832695  4.008059  4.180333  5.485005  6.024266  8.217778  8.544191
#> [15]  8.927070  9.081201  9.456097 10.873318 10.899277 11.086384