P.Mean: Accrual with a slow start-up (created 2012-10-13). News: Sign up for "The Monthly Mean," the newsletter that dares to call itself average, www.pmean.com/news.

Whenever we get feedback about our accrual models, we seem to get a comment along the lines that accrual is typically slower at the start of the study. It makes intuitive sense. There are inefficiencies at the start because the research team is getting used to the new protocol and it takes time for word to get out about the study. I wanted to examine an accrual model where the accrual rate increases linearly to a point and then levels off. It requires the use of WinBUGS software, but, surprisingly, is not all that hard to fit. Here are the details.

First, you need to understand a flat or constant accrual model. You specify a planned sample size (N) and a planned time duration (T), but then you specify your degree of certainty using a prior distribution. In Gajeweski 2008, we suggest that you get an initial feel for the strength of the prior by asking the researcher, "On a scale of 1-10, how certain are you that you will recruit N patients in T time?" Define P as 0.1 times the answer that you get to this question. Then consider a gamma(NP, NT) distribution for your prior distribution on the accrual rate.

You can define accrual in terms of an exponential waiting time between successive patients or you can define accrual in terms of a Poisson count of the number of patients that arrive on a particular day. While they are theoretically equivalent, there are pragmatic distinctions for these models. For example, if your trial has a fixed sample size and you will continue the trial for as long as it takes, then you will get a simpler answer using an exponential waiting time model. If your trial has a fixed duration, and you will stop on a specific date, no matter what, then you will get a simpler answer using a Poisson model.

Assume that you are partly through the trial, and you have an observed sample size (n) over an observed time period (t). The posterior distribution of the accrual rate is a gamma (NP+n, TP+t) distribution. The posterior mean accrual rate is (NP+n) / (TP+t) which can be shown to be a weighted average of the prior mean (N/T) and the mean accrual rate observed in the data (n/t). This posterior distribution is the same whether you use a model with exponential waiting times or Poisson counts.

You can get predictions for accrual for the remainder of the trial, and here there is a difference between the exponential and Poisson model.

Exponential model. For a trial with a fixed sample size, you can easily simulate the predicted duration of the remainder of the trial. Randomly draw a large number (1,000 or so) of values from the posterior gamma (NP+n, TP+t) distribution. The waiting time of the remaining N-n observations is going to be a sum of N-n exponential random variables or more simply a gamma distribution with shape parameter equal to N-n. Generate a large number of these gammas, each the same N-n value for the shape parameter, but with a different scale parameter drawn from the posterior distribution. Then calculate the percentiles of the preducted duration using simple percentiles from the large number of gamma variables that you just generated. To get the toal estimated time, add the observed waiting time so far (t).

It's a bit confusing, perhaps because you have gamma random variables from the posterior distribution plugging into the gamma distribution representing the estimated remaining waiting time. There's also some confusion because there are two equivalent forms of the exponential distribution.

Here's some simple code in R to do this simulation.

compute.exponentialsimulation <- function(n.current, t.current, n.planned, t.planned, P) {
  k <- n.planned*P+n.current
  V <- t.planned*P+t.current
  theta <- 1/rgamma(1000, shape=k, rate=V)
  lambda <- 1/theta
  t.projected <- t.current + rgamma(1000, shape=n.planned-n.current, rate=lambda)
  cbind(t.projected, theta, lambda)
}

The simulation here produces estimates of theta and lambda, which are inversely related, just for ease in debugging.

There is also a way to get this same simulation done in BUGS, and it is instructive because it shows how to generalize to more complex accrual patterns. Here is the BUGS code for the exponential simulation.

model {
  k <- n.planned * P
  V <- t.planned * P
  alpha[1] ~ dgamma(k, V)
  for (i in 1:n.current) {
    lambda[i] <- alpha[1]
    wait[i] ~ dexp(lambda[i])
  }
  for (j in 1:(n.planned - n.current)) {
    pred.lambda[j] <- alpha[1]
    pred.wait[j] ~ dexp(pred.lambda[j])
  }
  proj.time <- sum(wait[1:nobs]) + sum(pred.wait[])
}

There is some deliberate redundancy in this code to allow a smooth path for more complex accrual models. This code can be run with R using the BRugs package.

There is also a closed form solution that you can use. It involves the use of the beta prime distribution, a distibution closely related to the beta distribution, but one that involves a ratio of two gamma distributions. You can get percentiles of the distribution of the estimated remaining time for the study using percentiles from the beta distribution or percentiles from an F distribution. Here is the R code to get these percentiles.

compute.exponentialclosedform <- function(n.current, t.current, n.planned, t.planned, P) {
  p <- c(0.025,0.5,0.975)
  c1 <- n.planned-n.current
  c2 <- n.planned*P+n.current
  c3 <- t.planned*P+t.current
  bpct <- qbeta(p,c1,c2)
  solution.1 <- c3*bpct/(1-bpct)
  solution.2 <- c3*c1/c2*qf(p,2*c1,2*c2)
  t.current + rbind(solution.1,solution.2)
}

The derivation of this closed form solution appears elsewhere on this website.

Poisson model. Suppose instead that the trial has a fixed time duration and that it will end on a specific day, not matter how good or bad the recruited sample size is. In this setting, a Poisson model will allow you to predict the number of patients accrued before the end of the study. Randomly draw a large number (1,000 or so) of values from the posterior gamma (NP+n, TP+t) distribution. The number of subjects recruited for the remaining T-t days of the study is going to be a sum of T-t Poisson random variables, with rate parameter determined by the posterior distribution. This simplifies to a single Poisson random variable with a rate that is (T-t) times the rate drawn from the posterior distribution. Generate a large number of these Poissons. This gives an estimate of the remaining number of patients that you would expect to see. To get the toal estimated time, add the observed patient count so far (n).

Here's some simple R code to run this simulation.

compute.poissonsimulation <- function(n.current, t.current, n.planned, t.planned, P) {
  k <- n.planned*P+n.current
  V <- t.planned*P+t.current
  theta <- 1/rgamma(1000, shape=k, rate=V)
  lambda <- 1/theta
  n.projected <- n.current + rpois(1000, lambda=lambda*(t.planned-t.current))
  cbind(n.projected, theta, lambda)
}

As above, you can also get this simulation run in BUGS.

model {
  k <- nexp * P
  V <- texp * P
  alpha[1] ~ dgamma(k, V)
  for (i in 1:tobs) {
    lambda[i] <- alpha[1]
    coun[i] ~ dpois(lambda[i])
  }
  for (j in 1:(texp - tobs)) {
    pred.lambda[j] <- alpha[1]
    pred.coun[j] ~ dpois(pred.lambda[j])
  }
  proj.coun <- sum(coun[1:tobs]) + sum(pred.coun[])
}

Again, there is a bit of built-in redundancy.

There is a closed form solution based on the negative binomial distribution.

compute.poissonclosedform <- function(n.current, t.current, n.planned, t.planned, P) {
  p <- c(0.025,0.5,0.975)
  pr <- 1 - (t.planned-t.current) / (t.planned*P+t.planned)
  r <- n.planned*P+n.current
  n.current + qnbinom(p,r,pr)
}

The derivation appears elsewhere.

Model for a slow start-up rate. You can use a piece-wise linear function to model a slow start-up rate, though there are several other models that might also be reasonable. The piece-wise linear function starts out with an accrual rate equal to beta1, increases linearly with a slope beta2, and levels off at a certain point, beta3. The rate for the rest of the study is beta1+beta2*beta3. The following graph illustrates this piece-wise linear function.

You should, however, consider a re-parameterization of this model. The re-parameterization may produce more stable estimates, estimates with simpler interpretation, and estimates that are easier to specify prior distributions for.

Let alpha1 represent the final accrual rate, the accrual rate after the piece-wise linear function has leveled off. Then set alpha2 equal to the fraction of the accrual rate that you start at. Finally set alpha3 equal to the fraction of the total planned accrual time at which the piece-wise linear function transitions from a linear increase to a level accrual rate.

Here is a graph illustrating the reparameterization.

With this reparameterization, you can continue to use a gamma prior distribution for the final accrual rate, alpha1. Since alpha2 and alpha3 are constrained to be between 0 and 1, a beta distribution is a natural choice for the prior distribution for these parameters.You can also add additional constraints (e.g., the accrual rate has to level off before half of the planned accrual time).

Eliciting prior distributions for alpha2 and alpha3 would require that you ask the researcher how much slower the accrual rate is likely to be at the start of the study and how long into the study will it take before accrual levels off. Then you can inquire about the level of certainty that the researcher has for each of these estimates.

These graphs show accrual rate making a transition in a specified time frame. It is also possible to specify a transition in a specified number of subjects. You might, for example, hypothesize that the first ten subjects will accrue at a slower rate, but then the rate levels off for the remaining subjects.

Piece-wise linear exponential model. The code for this model in WinBUGS is

model {
  k <- nexp * P
  V <- texp * P
  alpha[1] ~ dgamma(k, V)
  alpha[2] ~ dbeta(1,1)
  alpha[3] ~ dbeta(1,1)
  beta1 <- alpha[1]*alpha[2]
  beta2 <- alpha[1]*(1-alpha[2])/alpha[3]
  beta3 <- alpha[3]*nexp
  for (i in 1:nobs) {
    lambda[i] <- beta1+beta2*min(i,beta3)
    wait[i] ~ dexp(lambda[i])
  }
  for (j in 1:(nexp - nobs)) {
    pred.lambda[j] <- beta1+beta2*min(j,beta3)
    pred.wait[j] ~ dexp(pred.lambda[j])
  }
  proj.time <- sum(wait[1:nobs]) + sum(pred.wait[])
}

This model uses a flat prior (beta(1,1) is a uniform distribution) for both alpha2 and alpha3. You can easily modify this to incorporate informative priors.

Piece-wise linear Poisson model. The code for this model in WinBUGS is

model {
  k <- nexp * P
  V <- texp * P
  alpha[1] ~ dgamma(k, V)
  alpha[2] ~ dbeta(1,1)
  alpha[3] ~ dunif(0,1)
  beta1 <- alpha[1]*alpha[2]
  beta2 <- alpha[1]*(1-alpha[2])/alpha[3]
  beta3 <- alpha[3]*nexp
  for (i in 1:tobs) {
    lambda[i] <- beta1+beta2*min(i,beta3)
    coun[i] ~ dpois(lambda[i])
  }
  for (j in 1:(texp - tobs)) {
    pred.lambda[j] <- beta1+beta2*min(j,beta3)
    pred.coun[j] ~ dpois(pred.lambda[j])
  }
  proj.coun <- sum(coun[1:tobs]) + sum(pred.coun[])
}

Again, it uses flat priors for both alpha2 and alpha3.

Here is a graph illustrating a fit to the piece-wise linear exponential model. The density on the left side of the graph represents estimates of accrual at the very start of the study (beta1 or alpha1*alpha2). The density on the right side of the graph represents estimates of accrual once things have leveled off (alpha1). The density at the top represents estimates of the transition between linearly increasing accrual and constant accrual (beta3 or alpha3*nexp). The red dots represent the "elbow" of the model, the (alpha2*nexp, alpha1) pair where the regression function bends from a linear increase to a flat line.

This graph uses a prior distribution which is uniform, but uniform across a restricted range (0 to 100). I want to experiment with different real and simulated accrual patterns and with different prior distributions. I'll report those results when I get them on these webpages.

This page was written by Steve Simon and is licensed under the Creative Commons Attribution 3.0 United States License. Need more information? I have a page with general help resources. You can also browse for pages similar to this one at Incomplete pages.