[Gsdesign-commits] r265 - in pkg/gsDesign: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue May 25 12:14:12 CEST 2010
Author: keaven
Date: 2010-05-25 12:14:12 +0200 (Tue, 25 May 2010)
New Revision: 265
Modified:
pkg/gsDesign/DESCRIPTION
pkg/gsDesign/NAMESPACE
pkg/gsDesign/R/gsMethods.R
pkg/gsDesign/R/gsqplot.R
Log:
V 2.3 adds gsBoundSummary and xtable to summarize boundary properties; also includes several plot fixes for ggplot versions of plots.
Modified: pkg/gsDesign/DESCRIPTION
===================================================================
--- pkg/gsDesign/DESCRIPTION 2010-03-11 17:08:19 UTC (rev 264)
+++ pkg/gsDesign/DESCRIPTION 2010-05-25 10:14:12 UTC (rev 265)
@@ -1,9 +1,10 @@
Package: gsDesign
-Version: 2.2-11
+Version: 2.3
Title: Group Sequential Design
Author: Keaven Anderson
Maintainer: Keaven Anderson <keaven_anderson at merck.com>
-Depends: ggplot2
+Depends: ggplot2, xtable
Description: gsDesign is a package that derives group sequential designs and describes their properties.
License: GPL (>= 2)
+Copyright: Copyright 2010, Merck Research Laboratories
Packaged: Thu Aug 6 17:28:16 2009; anderkea
Modified: pkg/gsDesign/NAMESPACE
===================================================================
--- pkg/gsDesign/NAMESPACE 2010-03-11 17:08:19 UTC (rev 264)
+++ pkg/gsDesign/NAMESPACE 2010-05-25 10:14:12 UTC (rev 265)
@@ -1,10 +1,12 @@
useDynLib(gsDesign)
export(gsBoundCP, gsCP, gsPP, gsPOS, gsCPOS, gsDensity)
+export(gsCPz, gsHR, gsDelta, gsBValue)
export(gsBound, gsBound1, gsDesign, gsProbability)
export(ciBinomial, nBinomial, simBinomial, testBinomial, gsBinomialExact)
export(nSurvival, nEvents)
export(normalGrid)
-export(plot.gsDesign, plot.gsProbability, print.gsProbability, print.gsDesign, print.nSurvival)
+export(plot.gsDesign, plot.gsProbability, print.gsProbability, print.gsDesign)
+export(print.nSurvival, xtable.gsDesign, gsBoundSummary)
export(sfPower, sfHSD, sfExponential, sfBetaDist, sfLDOF, sfLDPocock, sfPoints, sfLogistic, sfExtremeValue, sfExtremeValue2, sfLinear)
export(sfCauchy, sfNormal, sfTDist, spendingFunction)
export(checkScalar, checkVector, checkRange, checkLengths, isInteger)
Modified: pkg/gsDesign/R/gsMethods.R
===================================================================
--- pkg/gsDesign/R/gsMethods.R 2010-03-11 17:08:19 UTC (rev 264)
+++ pkg/gsDesign/R/gsMethods.R 2010-05-25 10:14:12 UTC (rev 265)
@@ -5,6 +5,9 @@
#
# print.gsDesign
# print.gsProbability
+# print.nSurvival
+# gsBoundSummary
+# xtable.gsDesign
#
# Hidden Functions:
#
@@ -316,3 +319,67 @@
cat("\n")
}
+gsBoundSummary <- function(x, upper=TRUE, ratio=1)
+{ if (upper)
+ { sn <- -1
+ bnd <- x$upper
+ }else
+ { sn <- 1
+ bnd <- x$lower
+ }
+ p0 <- cumsum(bnd$prob[,1])
+ p1 <- cumsum(bnd$prob[,x$theta==x$delta])
+ CP <- p0
+ CP[x$k] <- NA
+ CPhat <- CP
+ for(i in 1:x$k)
+ { if (i < x$k)
+ { CPtem <- gsCP(zi=bnd$bound[i], i=i, x=x, theta=NULL)$upper$prob
+ CP[i] <- sum(CPtem[,3])
+ CPhat[i] <- sum(CPtem[,1])
+ }
+ }
+ if (is.null(x$endpoint) || tolower(x$endpoint)!="survival")
+ { effect <- gsDelta(bnd$bound, i=1:x$k, x=x)
+ ename <- "delta"
+ }else
+ { effect <- gsHR(z=bnd$bound, i=1:x$k, x=x, ratio=1)
+ ename <- "HR"
+ }
+ nval <- x$n.I
+ if (x$n.fix != 1) nval <- ceiling(nval/2)*2
+ tab <- cbind(100*x$timing, nval, bnd$bound, pnorm(sn * bnd$bound), p0, p1, CP, CPhat,
+ effect,
+ gsBValue(bnd$bound, i=1:x$k, x=x))
+ rnames <- c(paste("IA",1:(x$k-1)), "Final")
+ cnames <- c("Timing (%)", ifelse(x$n.fix==1,"r","N"), "Z", "Nominal p", "H0 spend", "H1 spend",
+ "CP theta1", "CP thetahat", ename , "B-value")
+ dimnames(tab) <- list(rnames, cnames)
+ class(tab) <- "gsBoundSummary"
+ tab
+}
+xtable.gsDesign <- function(x, caption=NULL, label=NULL, align=NULL, digits=c(0,0,3,4,4,4,3,3,3,3),
+ display=NULL, upper=TRUE, rnames=NULL, cnames=NULL, ratio=1,
+ sanitize.text.function=function(x){x},
+ sanitize.rownames.function=function(x){x},...)
+{ bnd <- round(t(gsBoundSummary(x, upper, ratio)), digits)
+ if (is.null(cnames)) cnames <- colnames(bnd)
+ if (is.null(rnames))
+ { rnames <- rownames(bnd)
+ if (rnames[9] == "delta") ename <- "$\\hat{\\delta}$"
+ else if (rnames[9] == "HR") ename <- "HR"
+
+ if (x$n.fix==1) nname <- "r"
+ else if (!is.null(x$endpoint) && tolower(x$endpoint)=="survival") nname <- "Events"
+ else nname <- "N"
+ rnames <- c("Timing (\\%)", nname, "Z", "Nominal p", "H$_0$-spend", "$\\beta$-spend",
+ "$\\hbox{CP}(\\theta_1,\\hbox{Z})$", "$\\hbox{CP}(\\hat{\\theta}, \\hbox{Z})$",
+ ename, "B-value")
+ }
+ b <- matrix(as.character(bnd), nrow=nrow(bnd), ncol=ncol(bnd),
+ dimnames=list(rnames,cnames))
+ b
+ print(xtable(b, caption=caption, label=label, display=display, align=align,
+ sanitize.text.function=sanitize.text.function,...),
+ sanitize.rownames.function=sanitize.rownames.function)
+}
Modified: pkg/gsDesign/R/gsqplot.R
===================================================================
--- pkg/gsDesign/R/gsqplot.R 2010-03-11 17:08:19 UTC (rev 264)
+++ pkg/gsDesign/R/gsqplot.R 2010-05-25 10:14:12 UTC (rev 265)
@@ -5,6 +5,10 @@
#
# plot.gsDesign
# plot.gsProbability
+# gsCPz
+# gsHR
+# gsDelta
+# gsBValue
#
# Hidden Functions:
#
@@ -40,8 +44,14 @@
"plot.gsProbability" <- function(x, plottype=2, base=FALSE,...)
{
-# checkScalar(plottype, "integer", c(1, 7))
- invisible(do.call(gsPlotName(plottype), list(x, base=base,...)))
+# checkScalar(plottype, "integer", c(1, 7))
+ y <- x
+ if (max(x$n.I)>3) y$n.fix<-max(x$n.I)
+ else y$n.fix<-1
+ y$nFixSurv <- 0
+ y$test.type <- 3
+ y$timing <- y$n.I/max(y$n.I)
+ do.call(gsPlotName2(plottype), list(y, base=base,...))
}
###
@@ -63,32 +73,49 @@
plottype <- match.arg(tolower(as.character(plottype)), as.vector(unlist(plots)))
names(plots)[which(unlist(lapply(plots, function(x, type) is.element(type, x), type=plottype)))]
}
+"gsPlotName2" <- function(plottype)
+{
+ # define plots and associated valid plot types
+ plots <- list(plotgsZ=c("1","z"),
+ plotgsPower=c("2","power"),
+ plotgsCP=c("4", "cp", "copp"),
+ plotASN=c("6","asn", "e{n}","n"),
+ plotBval=c("7","b","b-val","b-value"),
+ plotHR=c("8","hr","hazard"))
+
+ # perform partial matching on plot type and return name
+ plottype <- match.arg(tolower(as.character(plottype)), as.vector(unlist(plots)))
+ names(plots)[which(unlist(lapply(plots, function(x, type) is.element(type, x), type=plottype)))]
+}
"plotgsZ" <- function(x, ylab="Normal critical value",main="Normal test statistics at bounds",...){qplotit(x=x,ylab=ylab,main=main,fn=function(z,...){z},...)}
-"plotBval" <- function(x, ylab="B-value",main="B-values at bounds",...){qplotit(x=x, fn=gsBvalue, ylab=ylab,main=main,...)}
+"plotBval" <- function(x, ylab="B-value", main="B-values at bounds",...){qplotit(x=x, fn=gsBValue, ylab=ylab, main=main, ...)}
"plotreleffect" <- function(x=x, ylab=NULL, main="Treatment effect at bounds",...){
- qplotit(x, fn=gsDeltaHat, main=main, ylab=ifelse(!is.null(ylab), ylab,
+ qplotit(x, fn=gsDelta, main=main, ylab=ifelse(!is.null(ylab), ylab,
ifelse(tolower(x$endpoint)=="binomial",
expression(hat(p)[C]-hat(p)[E]),
expression(hat(theta)/theta[1]))),...)
}
-"plotHR" <- function(x=x, ylab="Estimated hazard ratio",main="Hazard ratio at bounds",...){qplotit(x, fn=gsHRHat, ylab=ylab,main=main,...)}
-gsBvalue <- function(z,i,x,ylab="B-value",...)
+"plotHR" <- function(x=x, ylab="Estimated hazard ratio",main="Hazard ratio at bounds",...){qplotit(x, fn=gsHR, ylab=ylab, main=main, ratio=1, ...)}
+gsBValue <- function(z,i,x,ylab="B-value",...)
{ Bval <- z * sqrt(x$timing[i])
Bval
}
-gsDeltaHat <- function(z, i, x, ylab=NULL,...)
+gsDelta <- function(z, i, x, ylab=NULL,...)
{ deltaHat <- z / sqrt(x$n.I[i]) * (x$delta1-x$delta0) / x$delta + x$delta0
deltaHat
}
-gsHRHat <- function(z, i, x, ratio, ylab="Estimated hazard ratio",...)
+gsHR <- function(z, i, x, ratio=1, ylab="Estimated hazard ratio", ...)
{ c <- 1 / (1 + ratio)
psi <- c * (1 - c)
hrHat <- exp(-z / sqrt(x$n.I[i] * psi))
+ hrHat
}
-gsCPfn <- function(z, i, x, theta, ...)
+gsCPz <- function(z, i, x, theta=NULL, ylab=NULL, ...)
{ cp <- array(0,length(z))
- for(j in 1:length[z])
- { cp[i] <- sum(gsCP(x=x, theta=theta, i=j, zi=z[j])$upper$prob)
+ if (is.null(theta)) theta <- z / sqrt(x$n.I[i])
+ if (length(theta)==1 && length(z)>1) theta <- array(theta,length(z))
+ for(j in 1:length(z))
+ { cp[j] <- sum(gsCP(x=x, theta=theta[j], i=i[j], zi=z[j])$upper$prob)
}
cp
}
@@ -114,15 +141,26 @@
ntx <- "n="
if (is.null(xlab)) xlab <- "Sample size"
}
- z <- fn(z=c(x$upper$bound,x$lower$bound), i=c(1:x$k, 1:x$k), x=x,
+ if (x$test.type > 1)
+ { z <- fn(z=c(x$upper$bound,x$lower$bound), i=c(1:x$k, 1:x$k), x=x,
ratio=ratio, delta0=delta0, delta=delta)
- Ztxt <- as.character(c(round(z[1:(x$k-1)],dgt[2]), round(z[x$k],max(dgt)),
+ Ztxt <- as.character(c(round(z[1:(x$k-1)],dgt[2]), round(z[x$k],max(dgt)),
round(z[(x$k+1):(2*x$k-1)], dgt[1]), round(z[2*x$k],max(dgt))))
- y <- data.frame(
- N=as.numeric(c(x$n.I,x$n.I)),
- Z=as.numeric(z),
- Bound=c(array("Upper", x$k), array("Lower", x$k)),
- Ztxt=Ztxt)
+ y <- data.frame(
+ N=as.numeric(c(x$n.I,x$n.I)),
+ Z=as.numeric(z),
+ Bound=c(array("Upper", x$k), array("Lower", x$k)),
+ Ztxt=Ztxt)
+ }else{
+ z <- fn(z=x$upper$bound, i=1:x$k, x=x,
+ ratio=ratio, delta0=delta0, delta=delta)
+ Ztxt <- as.character(round(z,dgt[1]))
+ y <- data.frame(
+ N=as.numeric(x$n.I),
+ Z=as.numeric(z),
+ Bound=array("Upper", x$k),
+ Ztxt=Ztxt)
+ }
if (!is.numeric(ylim))
{ ylim <- range(y$Z)
ylim[1] <- ylim[1] -.1 * (ylim[2] - ylim[1])
@@ -138,12 +176,22 @@
}else
{ GeomText$guide_geom <- function(.) "blank"
lbls <- c("Lower", "Upper")
- p <- qplot(x=as.numeric(N), y=as.numeric(Z), data=y,
- group=factor(Bound), colour=factor(Bound), geom=geom, label=Ztxt,
- xlab=xlab, ylab=ylab, main=main,
- ylim=ylim, xlim=xlim,...) + aes(lty=factor(Bound)) +
- scale_colour_manual(name= "Bound", values=col, labels=lbls, breaks=lbls) +
- scale_linetype_manual(name= "Bound", values=lty, labels=lbls, breaks=lbls)
+ if (x$test.type > 1)
+ { p <- qplot(x=as.numeric(N), y=as.numeric(Z), data=y,
+ group=factor(Bound), colour=factor(Bound), geom=geom, label=Ztxt,
+ xlab=xlab, ylab=ylab, main=main,
+ ylim=ylim, xlim=xlim,...) + aes(lty=factor(Bound)) +
+ scale_colour_manual(name= "Bound", values=col, labels=lbls, breaks=lbls) +
+ scale_linetype_manual(name= "Bound", values=lty, labels=lbls, breaks=lbls)
+ }else{
+ p <- qplot(x=as.numeric(N), y=as.numeric(Z), data=y,
+ aes(colour=col[1], lty=lty[1], lwd=lwd[1]), geom=geom, label=Ztxt,
+ xlab=xlab, ylab=ylab, main=main,
+ ylim=ylim, xlim=xlim,...)
+ p <- ggplot(aes(x=as.numeric(N), y=as.numeric(Z), label=Ztxt), data=y) +
+ geom_line(colour=col[1], lty=lty[1], lwd=lwd[1]) +
+ geom_text() + xlab(xlab) + ylab(ylab) + opts(title=main)
+ }
}
if (nlabel==TRUE)
{ y2 <- data.frame(
@@ -151,7 +199,6 @@
Z=as.numeric(array(ylim[1], x$k)),
Bound=array("Lower", x$k),
Ztxt=ifelse(array(nround,x$k) > 0, as.character(round(x$n.I, nround)), ceiling(x$n.I)))
-#browser()
if (base)
{ text(x=y2$N, y=y$Z, y$Ztxt, cex=cex)
}
@@ -234,7 +281,7 @@
{ if (test.type == 1)
{
plot(x$n.I[1:(x$k-1)], y, xlab=xlab, ylab=ylab, main = main,
- ylim=c(ymin, ymax), xlim=xlim, col=col[2], lwd=lwd[2], lty=lty[2], type="l", ...)
+ ylim=c(ymin, ymax), xlim=xlim, col=col[1], lwd=lwd[1], lty=lty[1], type="l", ...)
points(x$n.I[1:(x$k-1)], y, ...)
text(x$n.I[1:(x$k-1)], y, as.character(round(y,dgt[2])), cex=textcex)
ymid <- ymin
@@ -274,7 +321,7 @@
}else
{ p <- qplot(x=as.numeric(N), y=as.numeric(CP), data=y, main=main,
label=Ztxt, geom="text",
- xlab=xlab, ylab=ylab, ylim=c(ymin, ymax), xlim=xlim) + geom_line(colour=col[2])
+ xlab=xlab, ylab=ylab, ylim=c(ymin, ymax), xlim=xlim) + geom_line(colour=col[1],lty=lty[1],lwd=lwd[1])
} }
if (nlabel==TRUE)
{ y2 <- data.frame(
@@ -371,21 +418,23 @@
if (x$test.type < 5)
{ spendb <- x$lower$sf(x$beta, t, x$lower$param)$spend/x$beta
}else
- { spendb <- x$lower$sf(x$alphastar, t, x$lower$param)$spend/x$alphastar
+ { spendb <- x$lower$sf(x$astar, t, x$lower$param)$spend/x$astar
}
group <- array(1, length(t))
q <- data.frame(t=c(t,t), spend=c(spenda, spendb), group=c(group,2*group))
p <- qplot(x=t, y=spend, data=q, geom="line", ylab=ylab, xlab=xlab, main=main,
- group=factor(group), linetype=factor(group), colour=factor(group)) +
- scale_colour_manual(name="Spending",values=col, labels=c(expression(alpha),expression(beta)), breaks=1:2) +
- scale_linetype_manual(name="Spending",values=lty, labels=c(expression(alpha),expression(beta)), breaks=1:2)
+ group=factor(group), linetype=factor(group), colour=factor(group))
+ p <- p +
+ scale_colour_manual(name="Spending",values=col,
+ labels=c(expression(alpha),ifelse(x$test.type<5,expression(beta),expression(1-alpha))), breaks=1:2) +
+ scale_linetype_manual(name="Spending",values=lty,
+ labels=c(expression(alpha),ifelse(x$test.type<5,expression(beta),expression(1-alpha))), breaks=1:2)
return(p)
}
}
}
-
-
-"plotASN" <- function(x, xlab=NULL, ylab=NULL, main=NULL, theta=NULL, xval=NULL, type="l", base=FALSE,...)
+"plotASN" <- function(x, xlab=NULL, ylab=NULL, main=NULL, theta=NULL, xval=NULL, type="l",
+ base=FALSE,...)
{
if (is(x, "gsDesign") && x$n.fix == 1)
{
@@ -452,12 +501,11 @@
else ylab <- "Expected sample size"
}
if (base)
- { plot(xval, x$en, type=type, ylab=ylab, xlab=xlab, main=main, ...)
+ { plot(xval, x$en, type=type, ylab=ylab, xlab=xlab, main=main,...)
return(invisible(x))
}
else
{ q <- data.frame(x=xval, y=x$en)
-#browser()
p <- qplot(x=x, y=y, data=q, geom="line", ylab=ylab, xlab=xlab, main=main)
return(p)
}
@@ -469,6 +517,7 @@
{ if (is.null(main)) main <- "Boundary crossing probabilities by effect size"
if (length(col==1)) col=array(col,2)
if (length(lty==1)) lty=array(lty,2)
+ if (length(lwd==1)) lwd=array(lwd,2)
if (is.null(xval))
{
if (is(x, "gsDesign") && is.null(xlab))
@@ -505,7 +554,7 @@
itxt <- array("Interim",x$k-1)
itxt <- paste(itxt,1:(x$k-1),sep=" ")
- if ((is(x, "gsDesign") && test.type > 1) || !is(x, "gsDesign"))
+ if (is(x, "gsProbability") || (is(x, "gsDesign") && test.type > 1))
{
itxt <- c(itxt,"Final",itxt)
boundprob <- array(1, length(xval))
@@ -537,7 +586,6 @@
interim <- c(interim, 1:(x$k-1))
}
yt <- data.frame(theta=xv, interim=interim, bound=bound, prob=yval, itxt=itxt)
-#browser()
if (base)
{ col2 <- ifelse(length(col) > 1, col[2], col)
lwd2 <- ifelse(length(lwd) > 1, lwd[2], lwd)
@@ -597,15 +645,26 @@
invisible(x)
}
else
- { p <- qplot(x=theta, y=prob, data=subset(y,interim==1), main=main,
- colour=factor(bound), geom="line", xlab = xlab, ylab = ylab, ylim=c(0,1),
- group=factor(bound)) + aes(lty=factor(bound))
- p <- p + scale_colour_manual(name= "Probability", values=col, breaks=1:2, labels=c("Upper bound","1-Lower bound")) +
- scale_linetype_manual(name="Probability", values=lty, breaks=1:2, labels=c("Upper bound","1-Lower bound"))
+ { GeomText$guide_geom <- function(.) "blank"
+ p <- qplot(x=theta, y=prob, data=subset(y,interim==1), main=main,
+ colour=factor(bound), lty=lty, geom="line", xlab = xlab, ylab = ylab, ylim=c(0,1),
+ group=factor(bound)) + aes(lty=factor(bound))
+ if(test.type == 1)
+ { p <- p + scale_colour_manual(name= "Probability", values=col, breaks=1,
+ labels="Upper bound") +
+ scale_linetype_manual(name="Probability", values=lty[1], breaks=1,
+ labels="Upper bound")
+ }else{
+ p <- p + scale_colour_manual(name= "Probability", values=col, breaks=1:2,
+ labels=c("Upper bound","1-Lower bound")) +
+ scale_linetype_manual(name="Probability", values=lty, breaks=1:2,
+ labels=c("Upper bound","1-Lower bound"))
+ }
p <- p + geom_text(data=yt, aes(theta, prob, colour=factor(bound), group=1, label=itxt))
- for(i in 1:x$k) p <- p + geom_line(data=subset(y,interim==i&bound==1), colour=col[1], lty=lty[1])
+ for(i in 1:x$k) p <- p + geom_line(data=subset(y,interim==i&bound==1),
+ colour=col[1], lty=lty[1], lwd=lwd[1])
if (test.type > 2) for(i in 1:(x$k-1)) {
- p <- p + geom_line(data=subset(y,interim==i&bound==2), colour=col[2], lty=lty[2])
+ p <- p + geom_line(data=subset(y,interim==i&bound==2), colour=col[2], lty=lty[2], lwd=lwd[2])
}
return(p)
}
More information about the Gsdesign-commits
mailing list