[Gsdesign-commits] r351 - in pkg/gsDesign: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon May 27 20:30:09 CEST 2013
Author: keaven
Date: 2013-05-27 20:30:09 +0200 (Mon, 27 May 2013)
New Revision: 351
Modified:
pkg/gsDesign/NAMESPACE
pkg/gsDesign/R/gsBinomialExact.R
pkg/gsDesign/man/gsBinomialExact.Rd
pkg/gsDesign/man/gsDesign.Rd
Log:
binomialSPRT: new function to produce truncated sequential probability ratio test (SPRT) for a single sample binomial trial. Also included plotting for this design.
gsBinomialExact: updated plotting capabilities.
Modified: pkg/gsDesign/NAMESPACE
===================================================================
--- pkg/gsDesign/NAMESPACE 2013-05-27 11:27:22 UTC (rev 350)
+++ pkg/gsDesign/NAMESPACE 2013-05-27 18:30:09 UTC (rev 351)
@@ -3,7 +3,7 @@
export(gsCPz, gsHR, gsDelta, gsBValue, gsRR, hrn2z, hrz2n, zn2hr)
export(gsBound, gsBound1, gsDesign, gsProbability)
export(ciBinomial, nBinomial, simBinomial, testBinomial, gsBinomialExact, varBinomial)
-export(nSurvival, nEvents, nNormal)
+export(binomialSPRT, plot.binomialSPRT, plot.gsBinomialExact, nSurvival, nEvents, nNormal)
export(normalGrid,ssrCP)
export(plot.gsDesign, plot.gsProbability, print.gsProbability, print.gsDesign)
export(print.nSurvival, xtable.gsDesign, gsBoundSummary, xtable.gsSurv)
Modified: pkg/gsDesign/R/gsBinomialExact.R
===================================================================
--- pkg/gsDesign/R/gsBinomialExact.R 2013-05-27 11:27:22 UTC (rev 350)
+++ pkg/gsDesign/R/gsBinomialExact.R 2013-05-27 18:30:09 UTC (rev 351)
@@ -60,3 +60,82 @@
class(x) <- c("gsBinomialExact", "gsProbability")
return(x)
}
+# see http://theriac.org/DeskReference/viewDocument.php?id=65&SectionsList=3
+binomialSPRT<-function(p0=.05,p1=.25,alpha=.1,beta=.15,minn=10,maxn=35){
+ lnA <- log((1-beta)/alpha)
+ lnB <- log(beta/(1-alpha))
+ a <- log((1-p1)/(1-p0))
+ b <- log(p1/p0)-a
+ slope <- -a / b
+ intercept <- c(lnA,lnB)/b
+ upper <- ceiling(slope*(minn:maxn)+intercept[1])
+ lower <- floor(slope*(minn:maxn)+intercept[2])
+ lower[lower< -1] <- -1
+ indx <- (minn:maxn >= upper)|(lower>=0)
+ # compute exact boundary crossing probabilities
+ y <- gsBinomialExact(k=sum(indx),n.I=(minn:maxn)[indx],
+ theta=c(p0,p1),a=lower[indx],b=upper[indx])
+ y$intercept <- intercept
+ y$slope <- slope
+ y$alpha <- alpha
+ y$beta <- beta
+ y$p0 <- p0
+ y$p1 <- p1
+ class(y) <- c("binomialSPRT","gsBinomialExact","gsProbability")
+ return(y)
+}
+plot.gsBinomialExact <- function(x,plottype=1,...){
+ if (plottype==6){
+ theta<-(max(x$theta)-min(x$theta))*(0:50)/50+min(x$theta)
+ y <- gsBinomialExact(k=x$k,theta=theta,n.I=x$n.I,a=x$lower$bound,b=x$upper$bound)
+ xx <- data.frame(p=theta,EN=y$en)
+ p<-ggplot(data=xx,aes(x=p,y=EN)) + geom_line() + ylab("Expected sample size")
+ }else if(plottype==3){
+ xx <- data.frame(N=x$n.I,p=x$upper$bound/x$n.I,Bound="Upper")
+ xx <- rbind(xx, data.frame(N=x$n.I,p=x$lower$bound/x$n.I,Bound="Lower"))
+ p<-ggplot(data=xx,aes(x=N,y=p,group=Bound))+
+ geom_point() +
+ ylab("Rate at bound")
+ }else if (plottype==2){
+ theta<-(max(x$theta)-min(x$theta))*(0:50)/50+min(x$theta)
+ # compute exact boundary crossing probabilities
+ y <- gsBinomialExact(k=x$k,n.I=x$n.I,
+ theta=theta,a=x$lower$bound,b=x$upper$bound)
+ # compute probability of crossing upper bound for each theta
+ Power <- data.frame(rr=theta,
+ Percent=100*as.vector(matrix(1,ncol=length(y$n.I),nrow=1)%*%
+ y$upper$prob),
+ Outcome="Reject H0")
+ # compute probability of crossing lower bound
+ futility <- data.frame(rr=theta,
+ Percent=100*as.vector(matrix(1,ncol=length(y$n.I),nrow=1)%*%
+ y$lower$prob),
+ Outcome="Reject H1")
+ # probability of no boundary crossing
+ indeterminate <- data.frame(rr=theta,Percent=100-Power$Percent-futility$Percent,
+ Outcome="Indeterminate")
+ #combine and plot
+ outcome <- rbind(Power,futility,indeterminate)
+ p <- ggplot(data=outcome,aes(x=rr,y=Percent,lty=Outcome))+
+ geom_line()+
+ xlab("Underlying response rate")
+ }else{
+ xx <- data.frame(N=x$n.I,x=x$upper$bound,Bound="Upper")
+ xx <- rbind(xx, data.frame(N=x$n.I,x=x$lower$bound,Bound="Lower"))
+ p<-ggplot(data=xx,aes(x=N,y=x,group=Bound))+
+ geom_point() +
+ ylab("Number of responses")
+ }
+ return(p)
+}
+plot.binomialSPRT <- function(x,plottype=1,...){
+ p <- plot.gsBinomialExact(x,plottype=plottype,...)
+ if (plottype==1){
+ p <- p + geom_abline(intercept=x$intercept[1],
+ slope=x$slope) +
+ geom_abline(intercept=x$intercept[2],
+ slope=x$slope)
+ }
+ return(p)
+}
+
Modified: pkg/gsDesign/man/gsBinomialExact.Rd
===================================================================
--- pkg/gsDesign/man/gsBinomialExact.Rd 2013-05-27 11:27:22 UTC (rev 350)
+++ pkg/gsDesign/man/gsBinomialExact.Rd 2013-05-27 18:30:09 UTC (rev 351)
@@ -1,14 +1,27 @@
\name{gsBinomialExact}
\alias{gsBinomialExact}
\alias{print.gsBinomialExact}
+\alias{plot.gsBinomialExact}
+\alias{binomialSPRT}
+\alias{plot.binomialSPRT}
\title{3.4: One-Sample Exact Binomial Boundary Crossing Probabilities}
-\description{Computes power/Type I error and expected sample size for a group sequential design
+\description{
+
+\code{gsBinomialExact} computes power/Type I error and expected sample size for a group sequential design
in a single-arm trial with a binary outcome.
+This can also be used to compare event rates in two-arm studies.
The print function has been extended using \code{print.gsBinomialExact} to print \code{gsBinomialExact} objects; see examples.
+Similarly, a plot function has been extended using \code{plot.gsBinomialExact} to plot \code{gsBinomialExact} objects; see examples.
+
+\code{binomialSPRT} computes a truncated binomial sequential probability ratio test (SPRT) which is a specific instance of an exact binomial group sequential design for a single arm trial with a binary outcome. See, for example, \url{http://theriac.org/DeskReference/viewDocument.php?id=65&SectionsList=3}.
}
\usage{
gsBinomialExact(k=2, theta=c(.1, .2), n.I=c(50, 100), a=c(3, 7), b=c(20,30))
+binomialSPRT(p0,p1,alpha,beta,minn,maxn)
+\method{plot}{gsBinomialExact}(x,plottype=1,\dots)
+\method{plot}{binomialSPRT}(x,plottype=1,\dots)
+
}
\arguments{
\item{k}{Number of analyses planned, including interim and final.}
@@ -16,10 +29,18 @@
\item{n.I}{Sample size at analyses (increasing positive integers); vector of length k.}
\item{a}{Number of "successes" required to cross lower bound cutoffs for futility or harm at each analysis; vector of length k; -1 means no lower bound.}
\item{b}{Number of "successes" required to cross upper bound cutoffs for futility or harm at each analysis; vector of length k.}
-
+ \item{p0}{Lower of the two response (event) rates hypothesized.}
+ \item{p1}{Higher of the two response (event) rates hypothesized.}
+ \item{alpha}{Nominal probability of rejecting response (event) rate \code{p0} when it is true (Type I error rate).}
+ \item{beta}{Nominal probability of rejecting response (event) rate \code{p1} when it is trun (Type II error rate).}
+ \item{minn}{Minimum sample size at which sequential testing begins.}
+ \item{maxn}{Maximum sample size.}
+ \item{x}{Item of class \code{gsBinomialExact} or \code{binomialSPRT} for \code{print.gsBinomialExact}. Item of class \code{gsBinomialExact} for \code{plot.gsBinomialExact}. Item of class \code{binomialSPRT} for item of class \code{binomialSPRT}.}
+ \item{plottype}{1 produces a plot with counts of response at bounds (for \code{binomialSPRT}, also produces linear SPRT bounds); 2 produces a plot with power to reject null and alternate response rates as well as the probability of not crossing a bound by the maximum sample size; 3 produces a plot with the response rate at the boundary as a function of sample size when the boundary is crossed; 6 produces a plot of the expected sample size by the underlying event rate (this assumes there is no enrollment beyond the sample size where the boundary is crossed).}
+\item{\dots}{arguments passed through to \code{ggplot}}
}
\details{
-Based on the book "Group Sequential Methods with Applications to Clinical Trials,"
+\code{gsBinomialExact} is based on the book "Group Sequential Methods with Applications to Clinical Trials,"
Christopher Jennison and Bruce W. Turnbull, Chapter 12, Section 12.1.2 Exact Calculations for Binary Data.
This computation is often used as an approximation for the distribution of the number of events in one treatment group out of all events when the probability of an event is small and sample size is large.
@@ -27,7 +48,13 @@
On output, the values of \code{theta} input to \code{gsBinomialExact} will be the parameter values for which the boundary crossing probabilities and expected sample sizes are computed.
Note that a[1] equal to -1 lower bound at n.I[1] means 0 successes continues at interim 1; a[2]==0 at interim 2 means 0 successes stops trial for futility at 2nd analysis.
-For final analysis, set a[k] equal to b[k]-1 to incorporate all possibilities into non-positive trial; see example.
+For final analysis, set a[k] equal to b[k]-1 to incorporate all possibilities into non-positive trial; see example.
+
+The sequential probability ratio test (SPRT) is a sequential testing scheme allowing testing after each observation. This likelihood ratio is used to determine upper and lower cutoffs which are linear and parallel in the number of responses as a function of sample size.
+\code{binomialSPRT} produces a variation the the SPRT that tests only within a range of sample sizes.
+While the linear SPRT bounds are continuous, actual bounds are the integer number of response at or beyond each linear bound for each sample size where testing is performed. Because of the truncation and discretization of the bounds, power and Type I error achieve will be lower than the nominal levels specified by \code{alpha} and \code{beta} which can be altered to produce desired values that are achieved by the planned sample size. See also example that shows computation of Type I error when futility bound is considered non-binding.
+
+Plots produced include boundary plots, expected sample size, response rate at the boundary and power.
}
\value{
@@ -41,11 +68,19 @@
\item{upper}{A list of the same form as \code{lower} containing the upper bound and upper boundary crossing probabilities.}
\item{en}{A vector of the same length as \code{theta} containing expected sample sizes for the trial design
corresponding to each value in the vector \code{theta}.}
+
+\code{binomialSPRT} produces an object of class \code{binomialSPRT} that is an extension of the \code{gsBinomialExact} class. The values returned in addition to those returned by \code{gsBinomialExact} are:
+ \item{intercept}{A vector of length 2 with the intercepts for the two SPRT bounds.}
+ \item{slope}{A scalar with the common slope of the SPRT bounds.}
+ \item{alpha}{As input. Note that this will exceed the actual Type I error achieved by the design returned.}
+ \item{beta}{As input. Note that this will exceed the actual Type II error achieved by the design returned.}
+ \item{p0}{As input.}
+ \item{p1}{As input.}
}
\seealso{\code{\link{gsProbability}}}
\note{The manual is not linked to this help file, but is available in library/gsdesign/doc/gsDesignManual.pdf
in the directory where R is installed.}
-\author{Jon Hartzel with modifications for gsDesign package by Yevgen Tymofyeyev and Keaven Anderson \email{keaven\_anderson at merck.}}
+\author{Jon Hartzel, Yevgen Tymofyeyev and Keaven Anderson \email{keaven\_anderson at merck.}}
\references{
Jennison C and Turnbull BW (2000), \emph{Group Sequential Methods with Applications to Clinical Trials}.
Boca Raton: Chapman and Hall.
@@ -62,5 +97,19 @@
# because of "gsProbability" class above, following is equivalent to
# print.gsProbability(zz)
zz
+# also plot (see also plots below for \code{binomialSPRT})
+# add lines using geom_line()
+plot(zz) + geom_line()
+
+# now for SPRT examples
+x <- binomialSPRT(p0=.05,p1=.25,alpha=.1,beta=.2)
+# boundary plot
+plot(x)
+# power plot
+plot(x,plottype=2)
+# Response (event) rate at boundary
+plot(x,plottype=3)
+# Expect sample size at boundary crossing or end of trial
+plot(x,plottype=6)
}
\keyword{design}
Modified: pkg/gsDesign/man/gsDesign.Rd
===================================================================
--- pkg/gsDesign/man/gsDesign.Rd 2013-05-27 11:27:22 UTC (rev 350)
+++ pkg/gsDesign/man/gsDesign.Rd 2013-05-27 18:30:09 UTC (rev 351)
@@ -12,7 +12,8 @@
\method{print}{gsDesign}(x,\dots)
\method{xtable}{gsDesign}(x, caption = NULL, label=NULL, align=NULL, digits=NULL, display=NULL,
footnote = NULL, fnwid = "9cm", deltaname="delta",
- Nname="N",logdelta=FALSE, \dots)}
+ Nname="N",logdelta=FALSE, \dots)
+}
\arguments{
\item{k}{Number of analyses planned, including interim and final.}
\item{test.type}{\code{1=}one-sided \cr
More information about the Gsdesign-commits
mailing list