P.Mean: My entry into the Applications of R in Business Competition (created 2011-10-20, updated 2011-10-31).
News: Sign up for "The Monthly Mean," the newsletter that dares to call itself average, www.pmean.com/news. |

I recently heard about a contest sponsored by Revolution Analytics. *Revolution Analytics is offering $20,000 in prizes to the best examples of applying R to business problems. This competition is designed to grow the collection of on-line materials describing how to use R, and to spur adoption of R and Revolution R for business applications. The contest is open to all R users worldwide.* http://www.inside-r.org/howto/enter. I want to submit the R code that I developed for the accrual paper published by Byron Gajewski and me in 2008.

My submission is still an early draft, but you can look at it at

--> http://www.inside-r.org/howto/predict-duration-your-clinical-trial

The R code was originally published at

-->http://www.childrensmercy.org/stats/08/ExponentialAccrual.aspx

but I have tweaked it a bit here.

`duration.plot <- function(n,T,P,m,tm,B=1000,Tmax=2*T,sample.paths=0) {`

#

# n is the target sample size

# T is the target completion time

# P is the prior certainty (0 <= P <= 1)

# m is the sample observed to date

# tm is the time to date

# B is the number of simulated duration times

# Tmax is the upper limit on the time axis of the graph

#

# set P to zero for a non-informative prior

# set m and tm to zero if no data has been accumulated yet.

k <- n*P+m

V <- T*P+tm

#

# Draw B samples from the inverse gamma distribution

theta <- 1/rgamma(B,shape=k,rate=V)

#

# Create and fill a matrix of simulated trial durations
` `

simulated.duration <- matrix(NA,nrow=n-m,ncol=B)

for (i in 1:B) {

wait <- rexp(n-m,1/theta[i])

simulated.duration[,i] <- tm+cumsum(wait)

}

#

# Each column is a separate simulated trial of waiting

# times for the remaining n-m patients. Calculate quantiles

# across the rows to get confidence limits on the waiting

# times for patients m+1, m+2, ..., n.

lcl <- apply(simulated.duration,1,quantile,probs=0.025)

mid <- apply(simulated.duration,1,quantile,probs=0.500)

ucl <- apply(simulated.duration,1,quantile,probs=0.975)

#

# Graph the simulated accruals. Create a layout with
two graphs

# The top graph is a histogram of the total estimated accrual for

# the full trial. The bottom graph is a display of predicted

# waiting times for patients m+1, m+2, ..., n.

layout(matrix(c(1,2,2,2)))

#

# No margins on bottom, top, and right of histogram.

par(mar=c(0.1,4.1,0.1,0.1))

#

# create histogram with 40 bars.

duration.hist <- cut(simulated.duration[n-m,],

seq(0,Tmax,length=40))

barplot(table(duration.hist),horiz=FALSE,

axes=FALSE,xlab=" ",ylab=" ",space=0,

col="white",names.arg=rep(" ",39))

#

# No margins on boot, top, and right of graph of waiting

# times for patients m+1, m+2, ..., n.

par(mar=c(4.1,4.1,0.1,0.1))

#

# Layout graph and axes, but no data yet.

plot(c(0,Tmax),c(0,n),xlab="Time",

ylab="Number of patients",type="n")

#

# Draw a polygon in gray ranging from lower to upper

# quantile limits.

polygon(c(lcl,rev(ucl)),c((m+1):n,n:(m+1)),

density=-1,col="gray",border=NA)

#

# Draw a white line at the median.

lines(mid,(m+1):n,col="white")

#

# Draw a black line showing observed accrual through patient m.

segments(0,0,tm,m)

#

# Draw a samll number of simulated accrual pathways

while (sample.paths>0) {

lines(simulated.duration[,sample.paths],(m+1):n)

sample.paths <- sample.paths-1

}

#

# Calculate and return percentiles to estimated total accrual time.

pctiles <- matrix(NA,nrow=2,ncol=3)

dimnames(pctiles) <- list(c("Waiting time","Trial duration")

,c("2.5%","50%","97.5%"))

pctiles[1,] <- 1/qgamma(c(0.975,0.5,0.025),shape=k,rate=V)

pctiles[2,1] <- lcl[n-m]

pctiles[2,2] <- mid[n-m]

pctiles[2,3] <- ucl[n-m]

list(pct=pctiles,sd=simulated.duration)

}

The function calls that produced the three figures shown above are

`p0 <- duration.plot(n=350,T=3,P=0.5,m= 0,tm= 0,Tmax=10)`

` p1 <- duration.plot(n=350,T=3,P=0.5,m=41,tm=239/365,Tmax=10)`

` p2 <- duration.plot(n=350,T=3,P=0, m=41,tm=239/365,Tmax=10)`

The due date for a submission is October 31.

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 Category: R software.