[Gsdesign-commits] r222 - in pkg/gsDesign: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Dec 13 04:27:42 CET 2009
Author: keaven
Date: 2009-12-13 04:27:39 +0100 (Sun, 13 Dec 2009)
New Revision: 222
Added:
pkg/gsDesign/R/gsqplot.R
pkg/gsDesign/R/nEvents.R
pkg/gsDesign/man/gsDensity.Rd
Modified:
pkg/gsDesign/DESCRIPTION
pkg/gsDesign/NAMESPACE
pkg/gsDesign/R/gsCP.R
pkg/gsDesign/R/gsDesign.R
pkg/gsDesign/R/gsMethods.R
pkg/gsDesign/R/gsNormalGrid.R
pkg/gsDesign/man/gsCP.Rd
pkg/gsDesign/man/normalGrid.Rd
pkg/gsDesign/man/plot.gsDesign.Rd
Log:
Update to gsDesign 2.2 - (ggplot2, gsPP...)
Modified: pkg/gsDesign/DESCRIPTION
===================================================================
--- pkg/gsDesign/DESCRIPTION 2009-12-13 00:03:31 UTC (rev 221)
+++ pkg/gsDesign/DESCRIPTION 2009-12-13 03:27:39 UTC (rev 222)
@@ -1,8 +1,9 @@
Package: gsDesign
-Version: 2.1-0
+Version: 2.2-0
Title: Group Sequential Design
Author: Keaven Anderson
Maintainer: Keaven Anderson <keaven_anderson at merck.com>
+Depends: ggplot2
Description: gsDesign is a package that derives group sequential designs and describes their properties.
License: GPL (version 2 or later)
-Packaged: Mon Jan 14 15:39:21 2008; Anderkea
+Packaged: Thu Aug 6 17:28:16 2009; anderkea
Modified: pkg/gsDesign/NAMESPACE
===================================================================
--- pkg/gsDesign/NAMESPACE 2009-12-13 00:03:31 UTC (rev 221)
+++ pkg/gsDesign/NAMESPACE 2009-12-13 03:27:39 UTC (rev 222)
@@ -1,5 +1,5 @@
useDynLib(gsDesign)
-export(gsBoundCP, gsCP)
+export(gsBoundCP, gsCP, gsPP, gsPOS, gsCPOS)
export(gsBound, gsBound1, gsDesign, gsProbability)
export(ciBinomial, nBinomial, simBinomial, testBinomial, gsBinomialExact)
export(nSurvival)
Modified: pkg/gsDesign/R/gsCP.R
===================================================================
--- pkg/gsDesign/R/gsCP.R 2009-12-13 00:03:31 UTC (rev 221)
+++ pkg/gsDesign/R/gsCP.R 2009-12-13 03:27:39 UTC (rev 222)
@@ -4,6 +4,7 @@
# Exported Functions:
#
# gsCP
+# gsPP
# gsBoundCP
#
# Hidden Functions:
@@ -70,6 +71,63 @@
gsProbability(k=knew, theta=theta, n.I=Inew, a=anew, b=bnew, r=r)
}
+gsPP <- function(x, i=1, zi=0, theta=c(0,3), wgts=c(.5,.5), r=18, total=TRUE)
+{ if (!(is(x, "gsProbability") || is(x, "gsDesign")))
+ {
+ stop("gsPP must be called with class of x either gsProbability or gsDesign")
+ }
+
+ if (i < 1 || i >= x$k)
+ {
+ stop("gsPP must be called with i from 1 to x$k-1")
+ }
+
+ test.type <- ifelse(is(x, "gsProbability"), 3, x$test.type)
+
+ if (zi > x$upper$bound[i])
+ {
+ stop("gsPP must have x$lower$bound[i] <= zi <= x$upper$bound[i]")
+ }
+ else if (test.type > 1 && zi < x$lower$bound[i])
+ {
+ stop("gsPP must have x$lower$bound[i]<=zi<=x$upper$bound[i]")
+ }
+
+ if (!is.real(theta) || is.na(theta))
+ {
+ theta <- c(zi/sqrt(x$n.I[i]), x$theta)
+ }
+
+ if (i < 1 || i >= x$k)
+ {
+ stop("gsPP must be called with i from 1 to x$k-1")
+ }
+
+ test.type <- ifelse(is(x, "gsProbability"), 3, x$test.type)
+
+ if (zi > x$upper$bound[i])
+ {
+ stop("gsCP must have x$lower$bound[i] <= zi <= x$upper$bound[i]")
+ }
+ else if (test.type > 1 && zi < x$lower$bound[i])
+ {
+ stop("gsCP must have x$lower$bound[i]<=zi<=x$upper$bound[i]")
+ }
+
+ if (!is.real(theta) || is.na(theta))
+ {
+ theta <- c(zi/sqrt(x$n.I[i]), x$theta)
+ }
+
+ cp <- gsCP(x = x, i = i, theta = theta, zi = zi, r = r)
+ cp$upper$prob %*% wgts
+ gsDen <- gsDensity(x, theta=theta, i=i, zi=zi, r=r)
+ if (!total) return(cp$upper$prob %*% t(gsDen$density))
+ one <- array(1, x$k-i)
+ cpmarg <- one %*% cp$upper$prob
+ sum(gsDen$density * cpmarg)
+}
+
"gsBoundCP" <- function(x, theta="thetahat", r=18)
{
len <- x$k-1
@@ -78,10 +136,9 @@
if (is(x, "gsDesign") || theta != "thetahat")
{
- thetahi <- x$delta
+ thetahi <- array(x$delta, len)
if (test.type > 1) thetalow <- thetahi
- }
- else
+ }else
{
if (test.type>1) thetalow <- x$lower$bound[1:len]/sqrt(x$n.I[1:len])
thetahi <- x$upper$bound[1:len]/sqrt(x$n.I[1:len])
Modified: pkg/gsDesign/R/gsDesign.R
===================================================================
--- pkg/gsDesign/R/gsDesign.R 2009-12-13 00:03:31 UTC (rev 221)
+++ pkg/gsDesign/R/gsDesign.R 2009-12-13 03:27:39 UTC (rev 222)
@@ -7,6 +7,9 @@
# gsBound1
# gsDesign
# gsProbability
+# gsDensity
+# gsPOS
+# gsCPOS
#
# Hidden Functions:
#
@@ -247,6 +250,63 @@
x
}
+"gsPOS" <- function(x, theta, wgts)
+{
+ if (class(x) != "gsProbability" && class(x) != "gsDesign")
+ stop("x must have class gsProbability or gsDesign")
+ checkVector(theta, "numeric")
+ checkVector(wgts, "numeric")
+ checkLengths(theta, wgts)
+ x <- gsProbability(theta = theta, d=x)
+ one <- array(1, x$k)
+ as.real(one %*% x$upper$prob %*% wgts)
+}
+
+"gsCPOS" <- function(i, x, theta, wgts)
+{
+ if (class(x) != "gsProbability" && class(x) != "gsDesign")
+ stop("x must have class gsProbability or gsDesign")
+ checkScalar(i, "integer", c(1, x$k), c(TRUE, FALSE))
+ checkVector(theta, "numeric")
+ checkVector(wgts, "numeric")
+ checkLengths(theta, wgts)
+ x <- gsProbability(theta = theta, d=x)
+ v <- c(array(1, i), array(0, (x$k - i)))
+ pAi <- 1 - as.real(v %*% (x$upper$prob + x$lower$prob) %*% wgts)
+ v <- 1 - v
+ pAiB <- as.real(v %*% x$upper$prob %*% wgts)
+ pAiB / pAi
+}
+
+"gsDensity" <- function(x, theta=0, i=1, zi=0, r=18)
+{ if (class(x) != "gsDesign" && class(x) != "gsProbability")
+ stop("x must have class gsDesign or gsProbability.")
+ checkVector(theta, "numeric")
+ checkScalar(i, "integer", c(0,x$k), c(FALSE, TRUE))
+ checkVector(zi, "numeric")
+ checkScalar(r, "integer", c(1, 80))
+ den <- array(0, length(theta) * length(zi))
+ xx <- .C("gsdensity", den, as.integer(i), length(theta),
+ as.double(theta), as.double(x$n.I),
+ as.double(x$lower$bound), as.double(x$upper$bound),
+ as.double(zi), length(zi), as.integer(r))
+ list(zi=zi, theta=theta, density=matrix(xx[[1]], nrow=length(zi), ncol=length(theta)))
+}
+
+#"gsPosterior" <- function(x, theta=NULL, wgts=NULL, i=1, zi=0, r=18)
+#{ if (class(x) != "gsDesign" && class(x) != "gsProbability")
+# stop("x must have class gsDesign or gsProbability.")
+# checkVector(theta, "numeric")
+# checkVector(wgts, "numeric")
+# checkLengths(
+# checkScalar(i, "integer", c(0,x$k), c(FALSE, TRUE))
+# checkVector(zi, "numeric")
+# checkScalar(r, "integer", c(1, 80))
+
+
+#}
+
+
###
# Hidden Functions
###
Modified: pkg/gsDesign/R/gsMethods.R
===================================================================
--- pkg/gsDesign/R/gsMethods.R 2009-12-13 00:03:31 UTC (rev 221)
+++ pkg/gsDesign/R/gsMethods.R 2009-12-13 03:27:39 UTC (rev 222)
@@ -3,8 +3,6 @@
#
# Exported Functions:
#
-# plot.gsDesign
-# plot.gsProbability
# print.gsDesign
# print.gsProbability
#
@@ -33,18 +31,6 @@
# Exported Functions
###
-"plot.gsDesign" <- function(x, plottype=1, ...)
-{
-# checkScalar(plottype, "integer", c(1, 7))
- invisible(do.call(gsPlotName(plottype), list(x, ...)))
-}
-
-"plot.gsProbability" <- function(x, plottype=2, ...)
-{
-# checkScalar(plottype, "integer", c(1, 7))
- invisible(do.call(gsPlotName(plottype), list(x, ...)))
-}
-
"print.gsProbability" <- function(x, ...)
{
ntxt <- "N "
@@ -264,384 +250,6 @@
c(expression(paste("Reject ", H[0])), "Continue", expression(paste("Reject ", H[1]))))
}
-"gsPlotName" <- function(plottype)
-{
- # define plots and associated valid plot types
- plots <- list(plotgsZ=c("1","z"),
- plotgsPower=c("2","power"),
- plotreleffect=c("3","xbar","thetahat","theta"),
- plotgsCP=c("4","cp"),
- sfplot=c("5","sf"),
- plotASN=c("6","asn", "e{n}","n"),
- plotBval=c("7","b","b-val","b-value"))
-
- # 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, main="Group sequential stopping boundaries", ylab="Normal critical value",
- xlab=ifelse(x$n.I[x$k] < 3, "Sample size relative to fixed design", "N"),
- lty=1, col=1, pch=22, textcex=1, legtext=gsLegendText(test.type), ...)
-{
- test.type <- ifelse(is(x,"gsProbability"), 3, x$test.type)
- ymax <- ceiling(4 * max(x$upper$bound)) / 4
- ymin <- ifelse(test.type == 1, floor(4 * min(x$upper$bound)) / 4, floor(4 * min(x$lower$bound)) / 4)
-
- if (ymin == ymax)
- {
- ymin <- ymin - .25
- }
-
- ymid <- (x$upper$bound[2] + x$lower$bound[2]) / 2
-
- if (test.type == 1)
- {
- plot(x$n.I, x$upper$bound, xlab=xlab, ylab=ylab, main = main,
- ylim=c(ymin, ymax), type="l", ...)
- points(x$n.I, x$upper$bound, ...)
- ymid <- ymin
- }
- else
- {
- matplot(x$n.I, cbind(x$lower$bound, x$upper$bound), xlab=xlab, ylab=ylab, main = main,
- lty=lty, col=col, ylim=c(ymin, ymax), type="l", ...)
- matpoints(x$n.I, cbind(x$lower$bound, x$upper$bound), pch=pch, col=col, ...)
- text(x$n.I[2], ymin, legtext[3], cex=textcex)
- }
-
- text(x$n.I[2], ymid, legtext[2], cex=textcex)
- text(x$n.I[2], ymax, legtext[1], cex=textcex)
-
- invisible(x)
-}
-
-"plotBval" <- function(x, main="B-values at stopping boundaries", ylab="B-value",
- xlab=ifelse(x$n.I[x$k] < 3, "Sample size relative to fixed design", "N"),
- lty=1, col=1, pch=22, textcex=1, legtext=gsLegendText(test.type), ...)
-{
- test.type <- ifelse(is(x,"gsProbability"), 3, x$test.type)
-
- sqrtt <- sqrt(x$n.I / x$n.I[x$k])
- ymax <- ceiling(4 * max(sqrtt * x$upper$bound)) / 4
- ymin <- ifelse(test.type == 1, 0, min(0, floor(4 * min(sqrtt * x$lower$bound)) / 4))
-
- if (ymin == ymax)
- {
- ymin <- ymin - .25
- }
-
- ymid <- sqrtt[2] * (x$upper$bound[2] + x$lower$bound[2]) / 2
-
- if (test.type == 1)
- {
- plot(x$n.I, sqrtt * x$upper$bound, xlab=xlab, ylab=ylab, main = main,
- ylim=c(ymin, ymax), xlim=c(0, x$n.I[x$k]), type="l", ...)
- points(x$n.I, sqrtt * x$upper$bound, ...)
- ymid <- ymin
- }
- else
- {
- matplot(x$n.I, cbind(sqrtt * x$lower$bound, sqrtt * x$upper$bound), xlab=xlab, ylab=ylab, main = main,
- lty=lty, col=col, ylim=c(ymin, ymax), xlim=c(0, x$n.I[x$k]), type="l", ...)
- matpoints(x$n.I, cbind(sqrtt * x$lower$bound, sqrtt * x$upper$bound), pch=pch, col=col, ...)
- text(x$n.I[2], ymin, legtext[3], cex=textcex)
- }
- text(x$n.I[2], ymid, legtext[2], cex=textcex)
- text(x$n.I[2], ymax, legtext[1], cex=textcex)
-
- invisible(x)
-}
-
-"plotreleffect" <- function(x, main="Observed treatment effect at stopping boundaries",
- ylab=ifelse(ses && is(x, "gsDesign"), expression(paste("Observed effect (", hat(theta) / delta, ")")),
- expression(paste("Observed effect (", hat(theta), ")"))),
- xlab=ifelse(x$n.I[x$k] < 3, "Sample size relative to fixed design", "N"),
- lty=1, col=1, pch=22, textcex=1, legtext=gsLegendText(test.type), ses=TRUE, ...)
-{
- test.type <- ifelse(is(x,"gsDesign"), x$test.type, 3)
- psi <- if (ses && is(x, "gsDesign")) 1 / x$delta else 1
- thathi <- x$upper$bound / sqrt(x$n.I) * psi
- ymax <- ceiling(4 * max(thathi)) / 4
-
- if (test.type == 1)
- {
- ymid <- floor(4 * min(thathi)) / 4
- ymin <- ymid
- }
- else
- {
- thatlow <- x$lower$bound / sqrt(x$n.I) * psi
- ymin <- floor(4 * min(thatlow)) / 4
- ymid <- (thathi[2] + thatlow[2]) / 2
-
- if (ymin == ymax)
- {
- ymid <- ymid-.25
- }
- }
-
- if (test.type == 1)
- {
- plot(x$n.I, thathi, xlab=xlab, ylab=ylab, main = main,
- ylim=c(ymin, ymax), type="l", ...)
- points(x$n.I, thathi, ...)
- }
- else
- {
- matplot(x$n.I, cbind(thatlow, thathi), xlab=xlab, ylab=ylab, main = main,
- lty=lty, col=col, ylim=c(ymin, ymax), type="l", ...)
- matpoints(x$n.I, cbind(thatlow, thathi), pch=pch, col=col, ...)
- text(x$n.I[2], ymin, legtext[3], cex=textcex)
- }
-
- text(x$n.I[2], ymid, legtext[2], cex=textcex)
- text(x$n.I[2], ymax, legtext[1], cex=textcex)
-
- invisible(list(n.I=x$n.I, lower=thatlow, upper=thathi))
-}
-
-"plotgsCP" <- function(x, theta="thetahat", main="Conditional power at interim stopping boundaries",
- ylab=ifelse(theta == "thetahat", expression(paste("Conditional power at", " ", theta, "=", hat(theta))),
- expression(paste("Conditional power at", " ", theta, "=", delta))),
- xlab=ifelse(x$n.I[x$k] < 3, "Sample size relative to fixed design", "N"),
- lty=1, col=1, pch=22, textcex=1, legtext=gsLegendText(test.type), ...)
-{
- test.type <- ifelse(is(x,"gsProbability"), 3, x$test.type)
- y <- gsBoundCP(x, theta=theta)
- ymax <- 1.05
- ymin <- - 0.05
-
- if (x$k > 3)
- {
- xtext <- x$n.I[2]
- }
- else if (x$k == 3)
- {
- xtext <- (x$n.I[2] + x$n.I[1]) / 2
- }
- else
- {
- xtext <- x$n.I[1]
- }
-
- if (test.type > 1)
- {
- ymid <- (y[2, 2] + y[2, 1]) / 2
- }
-
- if (test.type == 1)
- {
- plot(x$n.I[1:(x$k-1)], y, xlab=xlab, ylab=ylab, main = main,
- ylim=c(ymin, ymax), type="l", ...)
- points(x$n.I[1:(x$k-1)], y, ...)
- ymid <- ymin
- }
- else
- {
- matplot(x$n.I[1:(x$k-1)], y, xlab=xlab, ylab=ylab, main = main,
- lty=lty, col=col, ylim=c(ymin, ymax), type="l", ...)
- matpoints(x$n.I[1:(x$k-1)], y, pch=pch, col=col, ...)
- text(xtext, ymin, legtext[3], cex=textcex)
- }
- text(xtext, ymid, legtext[2], cex=textcex)
- text(xtext, 1.03, legtext[1], cex=textcex)
-
- invisible(x)
-}
-
-"sfplot" <- function(x,
- xlab=ifelse(max(x$n.I) > 3, "N", "Proportion of total sample size"),
- ylab=expression(paste(alpha, "-spending")),
- ylab2=expression(paste(beta, "-spending")), oma=c(2, 2, 2, 2),
- legtext=if (x$test.type > 4) c("Upper bound", "Lower bound") else c(expression(paste(alpha, "-spending")),
- expression(paste(beta, "-spending"))),
- mai=c(.85, .75, .5, .5), lty=1:2, xmax=1, ...)
-{
- # K. Wills (GSD-31)
- if (is(x, "gsProbability"))
- {
- stop("Spending function plot not available for gsProbability object")
- }
-
- # K. Wills (GSD-30)
- if (x$upper$name %in% c("WT", "OF", "Pocock") || x$upper$parname=="Points")
- {
- stop("Spending function plot not available for pointwise spending functions")
- }
-
- t <- 0:40 / 40 * xmax
-
- if (x$test.type > 2)
- {
- par(mai=mai, oma=oma)
- }
-
- plot(t, x$upper$sf(x$alpha, t, x$upper$param)$spend, type="l", ylab=ylab, xlab=xlab, ...)
-
- if (x$test.type > 2)
- {
-# if (is.character(legtext) && legtext != "")
-# if (legtext != "")
- legend(x=c(.0, .43), y=x$alpha * c(.85, 1), lty=lty, legend=legtext)
- par(new=TRUE)
- plot(t, x$lower$sf(x$beta, t, x$lower$param)$spend, ylim=c(0, x$beta), type="l", ylab="", main="", yaxt="n", xlab=xlab, lty=2)
- axis(4)
- mtext(text=ylab2, side = 4, outer=TRUE)
- }
-}
-
-"plotASN" <- function(x, xlab=NULL, ylab=NULL, main=NULL, theta=NULL, xval=NULL, type="l", ...)
-{
- if (is(x, "gsDesign") && x$n.fix == 1)
- {
- if (is.null(ylab))
- {
- ylab <- "E{N} relative to fixed design"
- }
-
- if (is.character(main) && main == "Default")
- {
- main <- "Expected sample size relative to fixed design"
- }
- }
-
- if (is.null(ylab))
- {
- ylab <- "E{N} relative to fixed design"
- }
-
- if (is.null(main))
- {
- main <- "Expected sample size by treatment difference"
- }
-
- if (is.null(theta))
- {
- if (is(x,"gsDesign"))
- {
- theta <- seq(0, 2, .05) * x$delta
- }
- else
- {
- theta <- x$theta
- }
- }
-
- if (is.null(xval))
- {
- if (is(x, "gsDesign") && is.null(xlab))
- {
- xval <- theta / x$delta
-
- if (is.null(xlab))
- {
- xlab <- expression(theta / delta)
- }
- }
- else
- {
- xval <- theta
-
- if (is.null(xlab))
- {
- xlab <- expression(theta)
- }
- }
- }
-
- if (is.null(xlab))
- {
- xlab <- ""
- }
-
- x <- if (is(x, "gsDesign")) gsProbability(d=x, theta=theta) else
- gsProbability(k=x$k, a=x$lower$bound, b=x$upper$bound, n.I=x$n.I, theta=theta)
-
- plot(xval, x$en, type=type, ylab=ylab, xlab=xlab, main=main, ...)
-
- invisible(x)
-}
-
-"plotgsPower" <- function(x, main="Group sequential power plot", ylab="Boundary Crossing Probalities", xlab=NULL,
- title="Boundary", legtext=c("Upper", "Lower"), lty=c(1, 2), col=c(1, 2), lwd=1,
- theta=if (is(x, "gsDesign")) seq(0, 2, .05) * x$delta else x$theta, xval=NULL, ...)
-{
- if (is.null(xval))
- {
- if (is(x, "gsDesign") && is.null(xlab))
- {
- xval <- theta / x$delta
- xlab <- expression(theta/delta)
- }
- else
- {
- xval <- theta
- if (is.null(xlab)) xlab <- expression(theta)
- }
- }
-
- if (is.null(xlab)) xlab <- ""
-
- x <- if (is(x, "gsDesign")) gsProbability(d=x, theta=theta) else
- gsProbability(k=x$k, a=x$lower$bound, b=x$upper$bound, n.I=x$n.I, theta=theta)
-
- col2 <- ifelse(length(col) > 1, col[2], col)
- lwd2 <- ifelse(length(lwd) > 1, lwd[2], lwd)
- lty2 <- ifelse(length(lty) > 1, lty[2], lty)
-
- ylim <- if (is(x, "gsDesign") && x$test.type<=2) c(0, 1) else c(0, 1.25)
-
- plot(xval, x$upper$prob[1, ], xlab=xlab, main=main, ylab=ylab,
- ylim=ylim, type="l", col=col[1], lty=lty[1], lwd=lwd[1], yaxt = "n")
-
- if (is(x, "gsDesign") && x$test.type <= 2)
- {
- axis(2, seq(0, 1, 0.1))
- axis(4, seq(0, 1, 0.1))
- }
- else
- {
- axis(4, seq(0, 1, by =0.1))
- axis(2, seq(0, 1, .1), labels=1 - seq(0, 1, .1), col.axis=col2, col=col2)
- }
-
- if (x$k == 1)
- {
- return(invisible(x))
- }
-
- if ((is(x, "gsDesign") && x$test.type > 2) || !is(x, "gsDesign"))
- {
- lines(xval, 1-x$lower$prob[1, ], lty=lty2, col=col2, lwd=lwd2)
- plo <- x$lower$prob[1, ]
-
- for (i in 2:x$k)
- {
- plo <- plo + x$lower$prob[i, ]
- lines(xval, 1 - plo, lty=lty2, col=col2, lwd=lwd2)
- }
-
- temp <- legend("topleft", legend = c(" ", " "), col=col,
- text.width = strwidth("Lower"), lwd=lwd,
- lty = lty, xjust = 1, yjust = 1,
- title = title)
-
- text(temp$rect$left + temp$rect$w, temp$text$y,
- legtext, col=col, pos=2)
- }
-
- phi <- x$upper$prob[1, ]
-
- for (i in 2:x$k)
- {
- phi <- phi + x$upper$prob[i, ]
- lines(xval, phi, col=col[1], lwd=lwd[1], lty=lty[1])
- }
-
- invisible(x)
-}
-
"sfprint" <- function(x)
{
# print spending function information
Modified: pkg/gsDesign/R/gsNormalGrid.R
===================================================================
--- pkg/gsDesign/R/gsNormalGrid.R 2009-12-13 00:03:31 UTC (rev 221)
+++ pkg/gsDesign/R/gsNormalGrid.R 2009-12-13 03:27:39 UTC (rev 222)
@@ -53,7 +53,7 @@
xx <- .C("stdnorpts", r, b, z, w)
len <- sum(xx[[3]] <= b[2])
z <- xx[[3]][1:len] * sigma + mu
- w <- xx[[4]][1:len] * sigma * dnorm(z, mean=mu, sd=sigma)
-
- list(z=z, wgts=w)
+ w <- xx[[4]][1:len] * sigma
+ d <- dnorm(z, mean=mu, sd=sigma)
+ list(z=z, density=d, gridwgts=w, wgts=d*w)
}
Added: pkg/gsDesign/R/gsqplot.R
===================================================================
--- pkg/gsDesign/R/gsqplot.R (rev 0)
+++ pkg/gsDesign/R/gsqplot.R 2009-12-13 03:27:39 UTC (rev 222)
@@ -0,0 +1,555 @@
+##################################################################################
+# S3 methods for the gsDesign package
+#
+# Exported Functions:
+#
+# plot.gsDesign
+# plot.gsProbability
+#
+# Hidden Functions:
+#
+# gsLegendText
+# gsPlotName
+# plotgsZ
+# plotBval
+# plotreleffect
+# plotgsCP
+# sfplot
+# plotASN
+# plotgsPower
+# qplotit
+#
+# Author(s): Keaven Anderson, PhD.
+#
+# Reviewer(s):
+#
+# R Version: 2.9.1
+#
+##################################################################################
+
+###
+# Exported Functions
+###
+
+"plot.gsDesign" <- function(x, plottype=1, base=FALSE,...)
+{
+# checkScalar(plottype, "integer", c(1, 7))
+ # if (base) invisible(do.call(gsPlotName(plottype), list(x, ...)))
+ do.call(gsPlotName(plottype), list(x, base=base,...))
+}
+
+"plot.gsProbability" <- function(x, plottype=2, base=FALSE,...)
+{
+# checkScalar(plottype, "integer", c(1, 7))
+ invisible(do.call(gsPlotName(plottype), list(x, base=base,...)))
+}
+
+###
+# Hidden Functions
+###
+"gsPlotName" <- function(plottype)
+{
+ # define plots and associated valid plot types
+ plots <- list(plotgsZ=c("1","z"),
+ plotgsPower=c("2","power"),
+ plotreleffect=c("3","xbar","thetahat","theta"),
+ plotgsCP=c("4","copp"),
+ plotsf=c("5","sf"),
+ 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",...){qplotit(x,ylab=ylab,fn=function(z,...){z},...)}
+"plotBval" <- function(x, ylab="B-value",...){qplotit(x, fn=gsBvalue, ylab=ylab,...)}
+"plotreleffect" <- function(x, ylab=expression(hat(theta)/theta[1]), delta=1, delta0=0,...){qplotit(x, fn=gsDeltaHat, ylab=ylab, delta=delta, delta0=delta0,...)}
+"plotHR" <- function(x, ylab="Estimated hazard ratio",...){qplotit(x, fn=gsHRHat, ylab=ylab,...)}
+gsBvalue <- function(z,i,x,...)
+{ Bval <- z * sqrt(x$timing[i])
+ Bval
+}
+gsThetaHat <- function(z, i, x,...)
+{ thetaHat <- z / sqrt(x$n.I[i])/x$delta
+ thetaHat
+}
+gsDeltaHat <- function(z, i, x, delta, delta0=0,...)
+{ deltaHat <- z / sqrt(x$n.I[i]) * delta / x$delta - delta0
+ deltaHat
+}
+gsHRHat <- function(z, i, x, ratio,...)
+{ c <- 1 / (1 + ratio)
+ psi <- c * (1 - c)
+ hrHat <- exp(-z / sqrt(x$n.I[i] * psi))
+}
+gsCPfn <- function(z, i, x, theta, ...)
+{ 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)
+ }
+ cp
+}
+# qplots for z-values and transforms of z-values
+"qplotit" <- function(x, xlim=NULL, ylim=NULL, geom=c("line", "text"), zround=2, lty=c(1,1), col=c(1,1),
+ lwd=c(1,1), nlabel="TRUE", xlab="", ylab="", fn=function(z,i,x,...){z},
+ ratio=1, delta0=0, delta=.05, cex=1, base=FALSE,...)
+{ if(x$n.I[x$k] < 3)
+ { nround <- 3
+ ntx <- "r="
+ if (xlab=="") xlab <- "Information relative to fixed sample design"
+ } else
+ { nround <- 0
+ ntx <- "n="
+ if (xlab=="") xlab <- "N"
+ }
+ 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)
+ 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=as.character(round(z, zround)))
+ if (!is.numeric(ylim))
+ { ylim <- range(y$Z)
+ ylim[1] <- ylim[1] -.1 * (ylim[2] - ylim[1])
+ }
+ if (!is.numeric(xlim))
+ { xlim <- range(x$n.I)
+ xlim <- xlim + c(-.05,.05) * (xlim[2] - xlim[1])
+ }
+ if (base==TRUE)
+ { plot(x=y$N[y$Bound=="Upper"], y=y$Z[y$Bound=="Upper"], type="l", xlim=xlim, ylim=ylim,
+ lty=lty[1], col=col[1], lwd=lwd[1], xlab=xlab, ylab=ylab,...)
+ lines(x=y$N[y$Bound=="Lower"], y=y$Z[y$Bound=="Lower"], lty=lty[2], col=col[2], lwd=lwd[2])
+ }else
+ { p <- qplot(x=as.numeric(N), y=as.numeric(Z), data=y,
+ group=Bound, geom=geom, label=Ztxt,
+ xlab=xlab, ylab=ylab,
+ ylim=ylim, xlim=xlim,...)
+ }
+ if (nlabel==TRUE)
+ { y2 <- data.frame(
+ N=as.numeric(x$n.I),
+ Z=as.numeric(array(ylim[1], x$k)),
+ Bound=array("Ntxt", x$k),
+ Ztxt=as.character(round(x$n.I,nround)))
+ if (base)
+ { text(x=y2$N, y=y$Z, y$Ztxt, cex=cex)
+ }
+ if (max(x$n.I) < 3)
+ { if (base)
+ { text(x=y2$N, y=y2$Z, paste(array("r=",x$k), y2$Ztxt, sep=""), cex=cex)
+ }else
+ { p <- p + geom_text(data=y2, label=paste(array("r=",x$k), y2$Ztxt, sep=""))
+ }
+ }else
+ { if(base)
+ { text(x=y2$N, y=y2$Z, paste(array("N=",x$k), y2$Ztxt, sep=""), cex=cex)
+ }else
+ { p <- p + geom_text(data=y2, label=paste(array("N=",x$k), y2$Ztxt, sep=""))
+ }
+ } }
+ if (base)
+ { invisible(x)
+ }else
+ { return(p)
+ }
+}
+
+"plotgsCP" <- function(x, theta="thetahat", main="Conditional power at interim stopping boundaries",
+ ylab=ifelse(theta == "thetahat",
+ expression(paste("Conditional power at",
+ theta, "=", hat(theta),sep=" ")),
+ expression(paste("Conditional power at", theta, "=", delta))),
+ geom="line",
+ xlab=ifelse(x$n.I[x$k] < 3, "Sample size relative to fixed design", "N"), xlim=NULL,
+ lty=1, col=1, pch=22, textcex=1, legtext=gsLegendText(test.type), zround=3, nlabel=TRUE, base=FALSE,...)
+{
+ if (!is.numeric(xlim))
+ { xlim <- range(x$n.I[1:(x$k-1)])
+ xlim <- xlim + c(-.05,.05) * (xlim[2] - xlim[1])
+ }
+ if(x$n.I[x$k] < 3)
+ { nround <- 3
+ ntx <- "r="
+ if (xlab=="") xlab <- "Information relative to fixed sample design"
+ }else
+ { nround <- 0
+ ntx <- "n="
+ if (xlab=="") xlab <- "N"
+ }
+ test.type <- ifelse(is(x,"gsProbability"), 3, x$test.type)
+ y <- gsBoundCP(x, theta=theta)
+ ymax <- 1.05
+ ymin <- - 0.1
+
+ if (x$k > 3)
+ {
+ xtext <- x$n.I[2]
+ }else if (x$k == 3)
+ {
+ xtext <- (x$n.I[2] + x$n.I[1]) / 2
+ }else
+ {
+ xtext <- x$n.I[1]
+ }
+
+ if (test.type > 1)
+ {
+ ymid <- (y[2, 2] + y[2, 1]) / 2
+ }
+
+ if (base)
+ { 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, type="l", ...)
+ points(x$n.I[1:(x$k-1)], y, ...)
+ ymid <- ymin
+ }
+ else
+ {
+ matplot(x$n.I[1:(x$k-1)], y, xlab=xlab, ylab=ylab, main = main,
+ lty=lty, col=col, ylim=c(ymin, ymax), xlim=xlim, type="l", ...)
+ matpoints(x$n.I[1:(x$k-1)], y, pch=pch, col=col, ...)
+ text(xtext, ymin, legtext[3], cex=textcex)
+ }
+ text(xtext, ymid, legtext[2], cex=textcex)
+ text(xtext, 1.03, legtext[1], cex=textcex)
+ invisible(x)
+ }else
+ { N <- as.numeric(x$n.I[1:(x$k-1)])
+ CP <- y[,2]
+ Bound <- array("Upper", x$k-1)
+ Ztxt <- as.character(round(CP[1:(x$k-1)], zround))
+ if (test.type > 1)
+ { N <- c(N, N)
+ CP <- c(CP, y[,1])
+ Bound <- c(Bound, array("Lower", x$k-1))
+ Ztxt <- as.character(c(Ztxt ,round(y[,1],zround)))
+ }
+ y <- data.frame(N=N, CP=CP, Bound=Bound, Ztxt=Ztxt)
+ p <- qplot(x=as.numeric(N), y=as.numeric(CP), data=y,
+ group=Bound, geom=c("line","text"), label=Ztxt,
+ xlab=xlab, ylab=ylab, ylim=c(ymin, ymax), xlim=xlim,...)
+ }
+ if (nlabel==TRUE)
+ { y2 <- data.frame(
+ N=x$n.I[1:(x$k-1)],
+ CP=array(ymin/2, x$k-1),
+ Bound=array("Ntxt", x$k-1),
+ Ztxt=as.character(round(x$n.I[1:(x$k-1)],nround)))
+ if (base)
+ { text(x=y2$N, y=y$CP, y2$Ztxt, cex=textcex)
+ }
+ if (max(x$n.I) < 3)
+ { if (base)
+ { text(x=y2$N, y=y2$CP, paste(array("r=",x$k), y2$Ztxt, sep=""), cex=textcex)
+ }else
+ { p <- p + geom_text(data=y2,label=paste(array("r=",x$k), y2$Ztxt, sep=""))
+ }
+ }else
+ { if(base)
+ { text(x=y2$N, y=y2$Z, paste(array("N=",x$k), y2$Ztxt, sep=""), cex=textcex)
+ }else
+ { p <- p + geom_text(data=y2,label=paste(array("N=",x$k-1), y2$Ztxt, sep=""))
+ }
+ } }
+ if (base)
+ { invisible(x)
+ }else
+ { return(p)
+ }
+}
+"plotsf" <- function(x,
+ xlab="Proportion of total sample size",
+ ylab=expression(paste(alpha, "-spending")),
+ ylab2=expression(paste(beta, "-spending")), oma=c(2, 2, 2, 2),
+ legtext=if (x$test.type > 4) c("Upper bound", "Lower bound") else c(expression(paste(alpha, "-spending")),
+ expression(paste(beta, "-spending"))),
+ col=c(1,1), lwd=c(.5,.5), lty=c(1,2),
+ mai=c(.85, .75, .5, .5), xmax=1, base=FALSE,...)
+{
+ # K. Wills (GSD-31)
+ if (is(x, "gsProbability"))
+ {
+ stop("Spending function plot not available for gsProbability object")
+ }
+
+ # K. Wills (GSD-30)
+ if (x$upper$name %in% c("WT", "OF", "Pocock"))
+ {
+ stop("Spending function plot not available for boundary families")
+ }
+ if (x$upper$parname == "Points"){x$sfupar <- sfLinear}
+
+ t <- 0:40 / 40 * xmax
+
+ if (x$test.type > 2 && base)
+ {
+ par(mai=mai, oma=oma)
+ }
+ if (base)
+ { plot(t, x$upper$sf(x$alpha, t, x$upper$param)$spend, type="l", ylab=ylab, xlab=xlab, lty=lty[1],
+ lwd=lwd[1], col=col[1],...)
+ }
+ else if (x$test.type < 3)
+ { spend <- x$upper$sf(x$alpha, t, x$upper$param)$spend
+ q <- data.frame(t=t, spend=spend)
+ p <- qplot(x=t, y=spend, data=q, geom="line", ylab=ylab, xlab=xlab,...)
+ return(p)
+ }
+
+ if (x$test.type > 2)
+ {
+ if (base)
+ { legend(x=c(.0, .43), y=x$alpha * c(.85, 1), lty=lty, col=col, lwd=lwd, legend=legtext)
+ par(new=TRUE)
+ plot(t, x$lower$sf(x$beta, t, x$lower$param)$spend,
+ ylim=c(0, x$beta), type="l", ylab="", main="",
+ yaxt="n", xlab=xlab, lty=lty[2], lwd=lwd[2], col=col[2],...)
+ axis(4)
+ mtext(text=ylab2, side = 4, outer=TRUE)
+ }
+ else
+ { spenda <- x$upper$sf(x$alpha, t, x$upper$param)$spend/x$alpha
+ 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
+ }
+ group <- array(1, length(t))
+ q <- data.frame(t=c(t,t), spend=c(spenda, spendb), group=c(group,2*group))
+ ylab <- "Proportion of spending"
+ p <- qplot(x=t, y=spend, data=q, geom="line", ylab=ylab, xlab=xlab, group=factor(group),
+ linetype=factor(group),
+ colour=factor(group)) +
+ scale_colour_manual(name="Spending",values=col, labels=c(expression(alpha),expression(beta))) +
+ scale_linetype_manual(name="Spending",values=lty, labels=c(expression(alpha),expression(beta)))
+ return(p)
+ }
+ }
+}
+
+
+"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)
+ {
+ if (is.null(ylab))
+ {
+ ylab <- "E{N} relative to fixed design"
+ }
+
+ if (is.character(main) && main == "Default")
+ {
+ main <- "Expected sample size relative to fixed design"
+ }
+ }
+
+
+ if (is.null(main))
+ {
+ main <- "Expected sample size by treatment difference"
+ }
+
+ if (is.null(theta))
+ {
+ if (is(x,"gsDesign"))
+ {
+ theta <- seq(0, 2, .05) * x$delta
+ }
+ else
+ {
+ theta <- x$theta
+ }
+ }
+
+ if (is.null(xval))
+ {
+ if (is(x, "gsDesign") && is.null(xlab))
+ {
+ xval <- theta / x$delta
+
+ if (is.null(xlab))
+ {
+ xlab <- expression(theta / theta[1])
+ }
+ }
+ else
+ {
+ xval <- theta
+
+ if (is.null(xlab))
+ {
+ xlab <- expression(theta)
+ }
+ }
+ }
+
+ if (is.null(xlab))
+ {
+ xlab <- ""
+ }
+
+ x <- if (is(x, "gsDesign")) gsProbability(d=x, theta=theta) else
+ gsProbability(k=x$k, a=x$lower$bound, b=x$upper$bound, n.I=x$n.I, theta=theta)
+
+ if (is.null(ylab))
+ {
+ if (max(x$n.I) < 3) ylab <- "E{N} relative to fixed design"
+ else ylab <- "Expected sample size"
+ }
+ if (base)
+ { plot(xval, x$en, type=type, ylab=ylab, xlab=xlab, main=main, ...)
+ return(invisible(x))
+ }
+ else
+ { q <- data.frame(x=xval, y=x$en)
+ p <- qplot(x=x, y=y, data=q, geom="line", ylab=ylab, xlab=xlab, main=main,...)
+ return(p)
+ }
+}
+"plotgsPower" <- function(x, main="Group sequential power plot",
+ ylab="Cumulative Boundary Crossing Probality",
+ xlab=NULL, title="Boundary", legtext=c("Upper", "Lower"), lty=c(1, 2), col=c(1, 2), lwd=1, cex=1,
+ theta=if (is(x, "gsDesign")) seq(0, 2, .05) * x$delta else x$theta, xval=NULL, base=FALSE,...)
+{
+ if (is.null(xval))
+ {
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/gsdesign -r 222
More information about the Gsdesign-commits
mailing list