[Gsdesign-commits] r352 - in pkg/gsDesign: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed May 29 15:34:09 CEST 2013
Author: keaven
Date: 2013-05-29 15:34:09 +0200 (Wed, 29 May 2013)
New Revision: 352
Added:
pkg/gsDesign/R/ssrCP.R
pkg/gsDesign/R/varBinomial.R
pkg/gsDesign/R/xtable.R
pkg/gsDesign/man/ssrCP.Rd
Log:
Added varBinomial.R and sscrCP source code required with recent updates.
Added: pkg/gsDesign/R/ssrCP.R
===================================================================
--- pkg/gsDesign/R/ssrCP.R (rev 0)
+++ pkg/gsDesign/R/ssrCP.R 2013-05-29 13:34:09 UTC (rev 352)
@@ -0,0 +1,29 @@
+ssrCP <- function(z, theta=NULL, maxinc=2, overrun=0, beta = 0.1, cpadj=c(.5,1-beta), x=gsDesign(k=2, timing=.5, beta=beta)){
+ if (class(x)!="gsDesign") stop("x must be passed as an object of class gsDesign")
+ if (2 != x$k) stop("input group sequential design must have 2 stages (k=2)")
+ w <- x$timing[1]
+ if (is.null(theta)) theta <- z / sqrt(x$n.I[1])
+ CP <- pnorm(theta*sqrt(x$n.I[2]*(1-w))-(x$upper$bound[2]-z*sqrt(w))/sqrt(1-w))
+ n2 <- array(x$n.I[1]+overrun,length(z))
+ indx <- ((z>x$lower$bound[1])&(z<x$upper$bound[1]))
+ indx2 <- indx & ((CP < cpadj[1])|(CP>cpadj[2]))
+ n2[indx2] <- x$n.I[2]
+ indx <- indx & !indx2
+ n2[indx] <- (((x$upper$bound[2] - z[indx] * sqrt(w)) / sqrt(1 - w) - qnorm(beta))/theta[indx])^2 + x$n.I[1]
+ n2[n2>maxinc*x$n.I[2]]<-maxinc*x$n.I[2]
+ return(n2)
+}
+# Type I error if sufficient statistic used instead of
+# combination test
+
+#unadjTypeIErr <- function(maxinc=2, beta = 0.1, cpadj=c(.5,1-beta),
+# x=gsDesign(k=2, timing=.5, beta=beta)){
+# z <- normalGrid() # grid for interim z
+# B <- sqrt(x$timing[1])*z # interim B-value
+# thetahat <- z/sqrt(n.I[1]) # interim standardized effect size
+# cp <-
+ # stage 2 sample size
+# n2 <- ssrCP(z=z,maxinc=maxinc,beta=beta,cpadj=cpadj,x=x)
+ # Type I error for non-promising zone
+ # assuming non-binding futility
+
\ No newline at end of file
Added: pkg/gsDesign/R/varBinomial.R
===================================================================
--- pkg/gsDesign/R/varBinomial.R (rev 0)
+++ pkg/gsDesign/R/varBinomial.R 2013-05-29 13:34:09 UTC (rev 352)
@@ -0,0 +1,55 @@
+# blinded estimate of variance for 2-sample binomial
+varBinomial<-function(x,n,delta0=0,ratio=1,scale="Difference")
+{ # check input arguments
+ checkVector(x, "integer", c(1, Inf))
+ checkScalar(n, "integer", c(1, Inf))
+ checkScalar(ratio, "numeric", c(0, Inf), c(FALSE,FALSE))
+ scale <- match.arg(tolower(scale), c("difference", "rr", "or"))
+ # risk difference test - from Miettinen and Nurminen eqn (9)
+ p<-x/n
+ phi<-array(0,max(length(delta0),length(x),length(ratio),length(n)))
+ if (scale=="difference")
+ { checkScalar(delta0, "numeric", c(-1, 1),c(FALSE,FALSE))
+ p1<-p+ratio*delta0/(ratio+1)
+ p2<-p1-delta0
+ a<-1+ratio
+ b<- -(a+p1+ratio*p2-delta0*(ratio+2))
+ c<- delta0^2-delta0*(2*p1+a)+p1+ratio*p2
+ d<- p1*delta0*(1-delta0)
+ v<-(b/(3*a))^3-b*c/6/a^2+d/2/a
+ u<-sign(v)*sqrt((b/3/a)^2-c/3/a)
+ w<-(pi+acos(v/u^3))/3
+ p10<-2*u*cos(w)-b/3/a
+ p20<-p10+delta0
+ phi <- (p10*(1-p10)+p20*(1-p20)/ratio)*(ratio+1)
+ phi[delta0==0]<-p*(1-p)/ratio*(1+ratio)^2
+ }
+ else if (scale=="rr") # log(p2/p1)
+ { checkScalar(delta0, "numeric", c(-Inf, Inf),c(FALSE,FALSE))
+ RR<-exp(delta0)
+ if (delta0==0) phi<-(1-p)/p/ratio*(1+ratio)^2
+ else
+ { p1<-p*(ratio+1)/(ratio*RR+1)
+ p2<-RR*p1
+ a<-(1+ratio)*RR
+ b<- -(RR*ratio+p2*ratio+1+p1*RR)
+ c<- ratio*p2+p1
+ p10<-(-b-sqrt(b^2-4*a*c))/2/a
+ p20<-RR*p10
+ phi<-(ratio+1)*((1-p10)/p10+(1-p20)/ratio/p20)
+ } }
+ # log-odds ratio - based on asymptotic distribution of log-odds
+ # see vignette
+ else
+ { checkScalar(delta0, "numeric", c(-Inf, Inf),c(FALSE,FALSE))
+ OR<-exp(delta0)
+ a<-OR-1
+ c<- -p*(ratio+1)
+ b<- 1+ratio*OR+(OR-1)*c
+ p10<-(-b+sqrt(b^2-4*a*c))/2/a
+ p20<-OR*p10/(1+p10*(OR-1))
+ phi<-(ratio+1)*(1/p10/(1-p10)+1/p20/(1-p20)/ratio)
+ phi[delta0==0]<-1/p/(1-p)*(1+1/ratio)*(1+ratio)
+ }
+ return(phi/n)
+}
Added: pkg/gsDesign/R/xtable.R
===================================================================
--- pkg/gsDesign/R/xtable.R (rev 0)
+++ pkg/gsDesign/R/xtable.R 2013-05-29 13:34:09 UTC (rev 352)
@@ -0,0 +1,63 @@
+xtable.gsDesign <- function (x, caption = NULL, label=NULL, align=NULL, digits=NULL, display=NULL,
+ footnote = NULL, fnwid = "9cm", deltaname="delta",
+ Nname="N",logdelta=FALSE, ...)
+{
+ k <- x$k
+ deltafutility <- gsDelta(x=x,i=1:x$k,z=x$lower$bound[1:x$k])
+ deltaefficacy <- gsDelta(x=x,i=1:x$k,z=x$upper$bound[1:x$k])
+ deltavals <- c(x$delta0,x$delta1)
+ if (logdelta){
+ deltafutility <- exp(deltafutility)
+ deltaefficacy <- exp(deltaefficacy)
+ deltavals <- exp(deltavals)
+ }
+ stat <- c("Z-value", "p (1-sided)",
+ paste(deltaname,"at bound"),
+ paste("P\\{Cross\\} if ", deltaname,"=",
+ deltavals[1], sep = ""),
+ paste("P\\{Cross\\} if ", deltaname,"=",
+ deltavals[2],sep = ""))
+ st <- stat
+ for (i in 2:k) stat <- c(stat, st)
+ an <- array(" ", 5 * k)
+ tim <- an
+ enrol <- an
+ fut <- an
+ eff <- an
+ an[5 * (0:(k - 1)) + 1] <- c(paste("IA ",
+ as.character(1:(k - 1)), ": ", as.character(round(100 * x$timing[1:(k - 1)],1)),
+ "\\%", sep = ""), "Final analysis")
+ an[5 * (1:(k - 1)) + 1] <- paste("\\hline", an[5 * (1:(k -1)) + 1])
+ an[5 * (0:(k - 1)) + 2] <- paste(Nname,":", ceiling(x$n.I[1:k]))
+ fut[5 * (0:(k - 1)) + 1] <- as.character(round(x$lower$bound,2))
+ eff[5 * (0:(k - 1)) + 1] <- as.character(round(x$upper$bound,2))
+ asp <- as.character(round(pnorm(-x$upper$bound), 4))
+ asp[asp == "0"] <- "$< 0.0001$"
+ eff[5 * (0:(k - 1)) + 2] <- asp
+ bsp <- as.character(round(pnorm(-x$lower$bound), 4))
+ bsp[bsp == "0"] <- " $< 0.0001$"
+ fut[5 * (0:(k - 1)) + 2] <- bsp
+ asp <- as.character(round(deltafutility[1:x$k], 4))
+ fut[5 * (0:(k - 1)) + 3] <- asp
+ bsp <- as.character(round(deltaefficacy[1:x$k], 4))
+ eff[5 * (0:(k - 1)) + 3] <- bsp
+ asp <- as.character(round(cumsum(x$upper$prob[, 1]), 4))
+ asp[asp == "0"] <- "$< 0.0001$"
+ eff[5 * (0:(k - 1)) + 4] <- asp
+ bsp <- as.character(round(cumsum(x$lower$prob[, 1]), 5))
+ bsp[bsp == "0"] <- "$< 0.0001$"
+ fut[5 * (0:(k - 1)) + 4] <- bsp
+ asp <- as.character(round(cumsum(x$upper$prob[, 2]), 4))
+ asp[asp == "0"] <- "$< 0.0001$"
+ eff[5 * (0:(k - 1)) + 5] <- asp
+ bsp <- as.character(round(cumsum(x$lower$prob[, 2]), 4))
+ bsp[bsp == "0"] <- "$< 0.0001$"
+ fut[5 * (0:(k - 1)) + 5] <- bsp
+ neff <- length(eff)
+ if (!is.null(footnote))
+ eff[neff] <- paste(eff[neff], "\\\\ \\hline \\multicolumn{4}{p{",
+ fnwid, "}}{\\footnotesize", footnote, "}")
+ x <- data.frame(cbind(an, stat, fut, eff))
+ colnames(x) <- c("Analysis", "Value", "Futility", "Efficacy")
+ xtable(x,caption = caption, label=label, align=align, digits=digits, display=display, ...)
+}
Added: pkg/gsDesign/man/ssrCP.Rd
===================================================================
--- pkg/gsDesign/man/ssrCP.Rd (rev 0)
+++ pkg/gsDesign/man/ssrCP.Rd 2013-05-29 13:34:09 UTC (rev 352)
@@ -0,0 +1,91 @@
+\name{ssrCP}
+\alias{ssrCP}
+\title{Sample size re-estimation based on conditional power}
+\description{
+\code{ssrCP()} is used with 2-stage designs to update sample size at an interim analysis based on conditional power.
+Either the estimated treatment effect at the interim analysis or any chosen effect size can be used for sample size re-estimation.
+If not done carefully, these designs can be very inefficient. It is probably a good idea to compare any design to a simpler group sequential design; see, for example, Jennison and Turnbull, Statistics in Medicine, 2003.
+Method assumes a small Type I error is included with the interim analysis and that the design is an adaptation from a 2-stage group sequential design (Lehmacher and Wassmer, Biometrics, 1999).}
+\usage{
+ssrCP(z, theta = NULL, maxinc = 2, overrun = 0, beta = 0.1, cpadj = c(0.5, 1 - beta), x = gsDesign(k = 2, timing = 0.5, beta = beta))
+}
+\arguments{
+ \item{z}{Scalar or vector with interim standardized Z-value(s). Input of multiple values makes it easy to plot the revised sample size as a function of the interim test statistic.}
+ \item{theta}{If \code{NULL} (default), conditional power calculation will be based on estimated interim treatment effect. Otherwise, \code{theta} is the standardized effect size used for conditional power calculation. Using the alternate hypothesis treatment effect can be more efficient than the estimated effect size; see Liu and Chi, Biometrics, 2001.
+There is a large literature for designs with conditional power; two commonly discussed methods allowing use of sufficient statistics for testing rather than combination tests have been proposed by Chen, DeMets and Lan, 2004, and Mehta and Pocock, 200x
+}
+ \item{maxinc}{Maximum-fold increase in sample size from planned sample size.}
+ \item{overrun}{Number of patients enrolled but not analyzed at interim analysis.}
+ \item{beta}{Targeted Type II error (1 - targeted conditional power)}
+ \item{cpadj}{Interval of conditional power values for which sample size is to be increased.}
+ \item{x}{2-stage group sequential design generated using \code{gsDesign().}}
+}
+\details{See references and examples.}
+\value{The initial version of \code{ssrCP()} returns a vector of sample sizes corresponding to the interim test statistics supplied in \code{z}.}
+\references{
+Chen, YHJ, DeMets, DL and Lan, KKG. Increasing the sample size when the unblinded interim result is promising, Statistics in Medicine,
+23:1023-1038, 2004.
+
+Jennison, C and Turnbull, BW. Mid-course sample size modification in clinical trials based on the observed treatment effect. Statistics in Medicine, 22:971-993", 2003.
+
+Lehmacher, W and Wassmer, G. Adaptive sample size calculations in group sequential trials, Biometrics, 55:1286-1290, 1999.
+
+Liu, Q and Chi, GY., On sample size and inference for two-stage adaptive designs, Biometrics, 57:172-177, 2001.
+
+Mehta, C and Pocock, S. Adaptive increase in sample size when interim results are promising: A practical guide with examples, Statistics in Medicine, 30:3267-3284, 2011.
+}
+\seealso{\code{\link{gsDesign}}}
+
+\author{Keaven Anderson \email{keaven\_anderson at merck.}}
+
+\examples{
+# timing of interim analysis
+timing <- .5
+# upper spending function
+sfu <- sfHSD
+# upper spending function paramater
+sfupar <- -12
+# maximum sample size inflation
+maxinflation <- 2
+# assumed enrollment overrrun at IA
+overrun <- 25
+# interim z-values for plotting
+z <- seq(0,4,.025)
+# Type I error (1-sided)
+alpha <- .025
+# Type II error for design
+beta <- .1
+# Fixed design sample size
+n.fix <- 100
+# conditional power interval where sample
+# size is to be adjusted
+cpadj <- c(.3,.9)
+# targeted Type II error when adapting sample size
+betastar <- beta
+
+# generate a 2-stage group sequential design with
+# desired planned sample size (first 2 lines set planned sample size to that from above; normally)
+x<-gsDesign(k=2,n.fix=n.fix,timing=timing,sfu=sfu,sfupar=sfupar,alpha=alpha,beta=beta)
+nalt <- maxinflation*(x$n.I[2])
+par(mar=c(7, 4, 4, 4)+.1)
+plot(z,ssrCP(x=x,z=z,overrun=overrun,beta=betastar,cpadj=cpadj),type="l",xlim=c(0,3.5),axes=FALSE,xlab="",ylab="")
+lines(z,80+240*dnorm(z,mean=0),col=2)
+lines(z,80+240*dnorm(z,mean=sqrt(x$n.I[1])*x$theta[2]),col=3)
+lines(z,80+240*dnorm(z,mean=sqrt(x$n.I[1])*x$theta[2]/2),col=4)
+lines(z,80+240*dnorm(z,mean=sqrt(x$n.I[1])*x$theta[2]*.75),col=5)
+axis(side=2, at=75+25*(0:5), labels=as.character(75+25*(0:5)))
+mtext(side=2,"Sample size",line=2)
+axis(side=4, at=80+240*seq(0,.4,.1), labels=as.character(seq(0,.4,.1)))
+mtext(side=4,expression(paste("Density for ",z[1])),line=2)
+lines(x=c(-3.5,3.5),y=c(nalt,nalt),lty=2)
+w <- x$timing[1]
+theta <- seq(.5,3.5,.5) / sqrt(x$n.I[1])
+CP <- pnorm(theta*sqrt(x$n.I[2]*(1-w))-(x$upper$bound[2]-seq(.5,3.5,.5)*sqrt(w))/sqrt(1-w))
+axis(side=1,line=3,at=seq(.5,3.5,.5),labels=as.character(round(CP,2)))
+mtext(side=1,"CP",line=3.5,at=.25)
+axis(side=1,line=1,at=seq(0,3.5,.5),labels=as.character(seq(0,3.5,.5)))
+mtext(side=1,expression(z[1]),line=.75,at=-.15)
+axis(side=1,line=5,at=seq(.5,3.5,.5),labels=as.character(round(theta/x$theta[2],2)))
+mtext(side=1,expression(hat(theta)/theta[1]),line=5.5,at=.25)
+}
+\keyword{design}
More information about the Gsdesign-commits
mailing list