P.Mean: Fitting the homogenous accrual model in BUGS (created 2012-04-13).

News: Sign up for "The Monthly Mean," the newsletter that dares to call itself average, www.pmean.com/news.

Several years ago, I wrote some R code for the homegenous accrual model. This is the simplest case for accrual, with an inverse gamma prior on the waiting time between successive patients. I wanted to fit the same model in BUGS, because I want to look at some extensions and I wanted to start with something simple. I am not great at BUGS yet, but I got it to work in an hour. I'm using the R interface to Open BUGS (BRugs). Here is the code.

ac.model <- function() {
 k <- n*P
 V <- T*P
 lambda ~ dgamma(k,V)
 theta <- 1/lambda
 tm ~ dgamma(m,lambda)
 n.rem <- n-m
 t.pred ~ dgamma(n.rem,lambda)
}
writeModel(ac.model,"SimpleAccrualModel.txt")
ac.data <- list(k=175,V=547.5,tm=239,m=41,n=350)
bugsData(ac.data,"SimpleAccrualData.txt")
ac.init <- list(list(lambda=1/3,t.pred=1000))
bugsInits(ac.init,numChains=1,"SimpleAccrualInits.txt")
modelCheck("SimpleAccrualModel.txt")
modelData("SimpleAccrualData.txt")
modelCompile()
modelInits("SimpleAccrualInits.txt")
modelUpdate(1000)
samplesSet(c("theta","t.pred"))
modelUpdate(1000)
samplesStats("*")

Let's take this code apart line by line.

ac.model <- function() {

In BRugs, you need to create a "function" that looks like the BUGS program. You won't call this function, but rather you will save this later as a text file.

 k <- n*P
 V <- T*P

Suppose you are running a clinical trial that plans to recruit n patients. You want to monitor the accrual in the trial (how rapidly patients come in). To plan this trial, you need to specify a prior distribution for the waiting time between successive patients. The parameters for the prior distribution are defined by the two questions: How long do you think this trial will last, and on a scale of 1-10, how confident are you in this answer? Let T represent the answer to the first question. Take the answer to the second question and divide by 10 to get P. You may wish to adjust some of these parameters, but this is a starting point. Now define the two parameters of an inverse gamma as nP and TP, respectively. P represents the weight given to the prior distribution. If P=0.1, that means that the prior distribution has a weight equal to one-tenth the total sample size.

 lambda ~ dgamma(k,V)
 theta <- 1 / lambda

BUGS does not have an inverse gamma function. So specify a gamma prior on the rate (lambda) and set the mean (theta) equal to one over the rate.

 tm ~ dgamma(m,lambda)

In the accrual model, you observe a total of m patients in time tm. The waiting time between each individual patient is exponential, so the sum of the waiting times is gamma.

 n.rem <- n-m
 t.pred ~ dgamma(n.rem,lambda)
}

Finally, get a prediction for the sum of the remaining n-m waiting times (which is also gamma).

writeModel(ac.model,"SimpleAccrualModel.txt")

Now that you have the BUGS model, use the writeModel function to store it in a text file. You could skip all of the above steps and just have a text file waiting for BRugs. The advantage of specifying the model within an R function and then writing it to a text file is that you have everything (the data, the model, the initial values, etc.) in one place.

ac.data <- list(n=350,T=1095,P=0.5,tm=239,m=41)

This is the data that BUGS needs. The values of n and T are the planned sample size and planned duration of the clinical trial. P measures the relative strength of the prior distribution. As mentioned above, m is the number of patients observed so far, tm is the time required to get these patients, and n is the total number of patients required. For more complex problems, you will need the actual accrual times, but the homegenous model can do just fine with the summary statistics.

bugsData(ac.data,"SimpleAccrualData.txt")

Just like the model, the data needs to be written to a file.

ac.init <- list(list(lambda=1/3,t.pred=1000))

You don't really need to specify initial values. You could get BUGS to pick initial values using the bugsGenInits function. But if you want to look at generalizations of the homogenous accrual model, it probably makes sense to get in the habit of doing it now, because you'll probably need to specify initial values for many of the extensions. Notice that this is a list within a list. It's an odd structure, but if you need to look at multiple chains then you would have a list for each chain all nested inside one big list.

bugsInits(ac.init,numChains=1,"SimpleAccrualInits.txt")

You need a third text file for your inital values.

modelCheck("SimpleAccrualModel.txt")

Now just run the same steps you would in BUGS. Step one is checking the model. If you programmed this well, you will get the message "model is syntactically correct." If there's an error BUGS will try to tell you what it thinks the error is. Debugging BUGS is not easy or fun.

modelData("SimpleAccrualData.txt")

Step 2 is feeding the data to BUGS. If you do this properly, you will get the message "data loaded." If you left out some key data, BUGS will tell you.

modelCompile()

Step 3 is compiling the model. If you do this properly, you will get the message "model compiled." Sometimes there will be a bug in your program that does not appear until you compile the program.

modelInits("SimpleAccrualInits.txt")

Step 4 is feeding the initial values to BUGS. If this step goes well, you will see the message "Initializing chain 1: model is initialized."

modelUpdate(1000)

Step 5 is to run a "burn-in" period. The Markov Chain Monte Carlo sometimes takes a while to converge on the correct distribution. Normally you want to run things for about a thousand runs or so, but for complex problems, it may take longer to converge. There are diagnostics that you can use to check for convergence. This is a fairly simple problem, so I wasn't too worried about convergence.

samplesSet(c("theta","t.pred"))

Step 6 is to specify the parameters that you want to monitor. Here I am interested in the mean of the inverse gamma distribution and the prediction of the time for the remainder of the trial.

modelUpdate(1000)

Step 7 is to run some more of the Markov Chain Monte Carlo. Again, this is a simple problem so a thousand replications should work just fine. For more complex problems you may need a larger number here. Again, there are are diagnostics that can help you decide how many replications are needed. One simple check is to make sure that the simulation error (MC_error shown below) is smaller by one or two orders of magnitude than the standard deviation of your parameters (sd shown below).

samplesStats("*")

This produces the output you want. There are other options for additional statistics and graphs. Here is what the output looks like.

           mean      sd MC_error val2.5pc   median val97.5pc start sample
t.pred 1127.000 99.8100 3.147000  944.200 1122.000   1349.00  1001   1000
theta     3.656  0.2466 0.008574    3.197    3.642      4.17  1001   1000

You can check the calculations on the posterior distribution and the remaining predicted time using the statements

p <- c(0.025,0.5,0.975)
 c1 <- 350-41
c2 <- 350*0.5+41
c3 <- 1095*0.5+239
1/qgamma(rev(p),c2,c3)
bpct <- qbeta(p,c1,c2)
c3*bpct/(1-bpct)
c3*c1/c2*qf(p,2*c1,2*c2)

The output from the call to the qgamma function produces the list of the 2.5, 50, and 97.5 percentiles for the posterior distribution of theta.

[1] 3.200529 3.646830 4.180107

Notice that the these are fairly close to the percentiles that BUGS produced for theta. The call to the qbeta function and the call to the qf function both produce equivalent percentiles for the prediction of the remaining time for the clinical trial (t.pred).

[1] 946.7041 1125.6551 1340.8285

[1] 946.7041 1125.6551 1340.8285

You can find the mathematical derivation of the beta and F distribution results at
--> http://www.pmean.com/12/iowatalk.html

Since BUGS uses a sampling approach, its results will not agree perfectly (and should not agree perfectly) with the closed form solutions. The closed form solution is always preferable to a simulated result. It's worthwhile, though, to fit the simple homogenous accrual model in BUGS so that you can then look at some extensions where there is no closed form solution.

Creative Commons License 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 Accrual Problems.