Added README.md to .Rbuildignore so that R CMD check doesn't create a…

… NOTE

.Rbuildignore

@@ -3,5 +3,6 @@
 ^_pkgdown\.yml$
 ^docs$
 ^pkgdown$
+README.md
 ^README\.Rmd$
 ^\.github$

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: clinfun
 Title: Clinical Trial Design and Data Analysis Functions
-Version: 1.0.99.2
+Version: 1.1.0
 Depends: R (>= 2.0.0), graphics, stats
 Imports: mvtnorm
 Suggests: survival

NAMESPACE

@@ -2,10 +2,11 @@ useDynLib(clinfun)
 import(stats,graphics)
 importFrom(utils, combn)
 importFrom(mvtnorm, pmvnorm)
-export("aucVardiTest","calogrank","coxphQuantile","coxphCPE","coxphERR","deltaAUC","fe.power","fe.mdor","fe.ssize","CPS.ssize","mdrr","or2pcase","gsdesign.binomial","gsdesign.normal","gsdesign.survival","power.ladesign","permlogrank","jonckheere.test","ph2simon","oc.twostage.bdry","twostage.inference","ph2single","pselect","toxbdry","bdrycross.prob","roc.area.test","roc.curve","roc.perm.test")
+export("aucVardiTest","calogrank","coxphQuantile","coxphCPE","coxphERR","deltaAUC","fe.power","fe.mdor","fe.ssize","CPS.ssize","mdrr","or2pcase","gsdesign.binomial","gsdesign.normal","gsdesign.survival","power.ladesign","permlogrank","jonckheere.test","ph2simon","oc.twostage.bdry","twostage.inference","twostage.admissible","ph2single","pselect","toxbdry","futilbdry","bdrycross.prob","roc.area.test","roc.curve","roc.perm.test")
 S3method(print, ph2simon)
 S3method(plot, ph2simon)
 S3method(print, toxbdry)
+S3method(print, futilbdry)
 S3method(print, calogrank)
 S3method(print, permlogrank)
 S3method(print, gsdesign)

NEWS.md

@@ -1,9 +1,15 @@
+# clinfun 1.1.0 (02/22/2022)
+
+* Added added admissible two stage design function
+* Added added futilbdry function for sequential futility stopping
+
 # clinfun (development version)
 
 * Added a `NEWS.md` file to track changes to the package.
 * deposited development version on GitHub
 * created {pkgdown} website
 * added 1:r sample allocation in gsdesign functions
+* Fixed the warnings messages from deltaAUC.f
 
 # clinfun 1.0.15 (04/13/2018)
 

R/ph2simon.R

@@ -140,3 +140,35 @@ twostage.inference <- function(x, r1, n1, n, pu, alpha=0.05) {
   names(out) <- c("pumvue", "p.value", paste(100*alpha, "% LCL", sep=""), paste(100*(1-alpha), "% UCL", sep=""))
   out
 }
+
+# function to get all the admisssible 2 stage designs
+twostage.admissible  <- function(x) {
+  xout <- x$out
+  nmax <- x$nmax
+  nopt <- which(xout[,5]==min(xout[,5]))[1]
+  # admissible indicator for the rows
+  adm = rep(0, nopt)
+  adm[1] = 1 # minimax
+  adm[nopt] = 1 # optimal
+  mss = xout[1:nopt, 4] # max sample size
+  ess = xout[1:nopt, 5] # expected sample size
+  # search admissible points only if search space is not a singleton
+  istart = 1
+  while (nopt > istart) {
+    # slope of each point from the minimax
+    slopes = (ess[(istart+1):nopt] - ess[istart])/(mss[(istart+1):nopt]-mss[istart])
+    imin = which.min(slopes)[1]
+    adm[istart+imin] = 1
+    istart = istart+imin
+  }
+  xout = xout[1:nopt,][adm==1,]
+  # find weight probability at which each rule is admissible
+  dmss = diff(xout[,4])
+  dess = diff(xout[,5])
+  qwt = round(dess/(dess - dmss),3)
+  # admissible designs
+  xout = cbind(xout, c(qwt, 0), c(1, qwt))
+  colnames(xout)[7:8] = c("qLo", "qHi")
+  rownames(xout) = c("Minimax", rep("Admissible", nrow(xout)-2), "Optimal")
+  xout
+}

R/toxbdry.R

@@ -1,14 +1,14 @@
-toxbdry <- function(pLo, pHi, n, cP0=0.1, cP1=0.9, ngrid=6, niter=10, delta=0, priority=c("null","alt")) {
+# function to get high and low boundaries for parameters given
+getBoundary <- function(pLo, pHi, n, cP0=0.1, cP1=0.9, ngrid=6, niter=10, delta=0, priorityNull=TRUE) {
 # decide whether to prioritize the type I (null) or the type II (alt) error.
-  priority <- match.arg(priority)
   nlook <- length(n)
   ptox <- seq(pLo, pHi, length=ngrid)
 # boundary shape in the power family with (Pocock = 0 & O'Brien-Fleming = 0.5)
   if (delta < 0 | delta > 0.5) stop("delta should be between 0 and 0.5")
   ifrac <- n/n[nlook]
   bdryden <- ifrac^delta
 # error threshold is cP0 for null and cP1 for alt
-  ethresh <- ifelse(priority == "null", cP0, cP1)
+  ethresh <- ifelse(priorityNull, cP0, cP1)
   alpha0 <- cP0/nlook
   r0 <- qbinom(1-pnorm(qnorm(alpha0)/bdryden), n, pLo)
 #  print(r0)
@@ -18,14 +18,26 @@ toxbdry <- function(pLo, pHi, n, cP0=0.1, cP1=0.9, ngrid=6, niter=10, delta=0, p
 #  print(r1)
   pstop1 <- bdrycross.prob(n, r1, ptox)
 # depending on null or alt prob of interest is 1st or ngrid element
-  ii <- ifelse(priority == "null", 1, ngrid)
+  ii <- ifelse(priorityNull, 1, ngrid)
+# if priorityNull is false and cP1 is exceeded at alpha0 reduce alpha0
+  if (!priorityNull) {
+    while(pstop0[ii,2] > ethresh) {
+      alpha1 <- alpha0
+      pstop1 <- pstop0
+      alpha0 <- 0.5*alpha0
+      r0 <- qbinom(1 - pnorm(qnorm(alpha0)/bdryden), n, pLo)
+      pstop0 <- bdrycross.prob(n, r0, ptox)
+    }
+  }
+# increment alpha1 if cP1 is not met for alpha1
   while(pstop1[ii,2] < ethresh) {
     alpha0 <- alpha1
     pstop0 <- pstop1
     alpha1 <- 1.4*alpha1
     r1 <- qbinom(1 - pnorm(qnorm(alpha1)/bdryden), n, pLo)
     pstop1 <- bdrycross.prob(n, r1, ptox)
   }
+# now iterate to get boundaries that are above and bele ethresh values    
   iter <- 0
   while (sum(r0 - r1) > 1 & iter <= niter) {
     iter <- iter + 1
@@ -52,62 +64,137 @@ toxbdry <- function(pLo, pHi, n, cP0=0.1, cP1=0.9, ngrid=6, niter=10, delta=0, p
   bdry.oc <- cbind(pstop1, pstop0[,2:4])
   colnames(bdry.oc)[2:7] <- c("pcross.lo", "pstop.lo", "ess.lo", "pcross.hi", "pstop.hi", "ess.hi")
   if((pstop1[1,2] > cP0 & pstop0[1,2] > cP0) | (pstop1[ngrid,2] < cP1 & pstop0[ngrid,2] < cP1)) warning(paste("Max sample size", n[nlook], "may be small for the specified stopping probabilities\n"))
-  out <- list("looks"=n, "lo.bdry"=r1, "hi.bdry"=r0, "bdry.oc"=bdry.oc, "bdry.alpha"=c(alpha1, alpha0),"delta"=delta)
-  class(out) <- "toxbdry"
-  out
-}
-
-print.toxbdry <- function(x, ...) {
-# convert the boundary to human readable form
-  n <- x$looks
-  if (max(diff(n)) == 1) {
-    ii <- c(which(diff(x$lo.bdry) > 0), length(n))
-    bdry.lo <- paste(x$lo.bdry[ii], n[ii], sep="/")
-    ii <- c(which(diff(x$hi.bdry) > 0), length(n))
-    bdry.hi <- paste(x$hi.bdry[ii], n[ii], sep="/")
-  } else {
-    bdry.lo <- paste(x$lo.bdry, n, sep="/")
-    bdry.hi <- paste(x$hi.bdry, n, sep="/")
-  }
-  cat("\n Toxicity boundary based on repeated significance testing \n")
-  cat("    Boundary shape parameter delta =", x$delta, "\n\n")
-  cat(" ******************************************************************\n")
-  cat(" * Stop if the number of toxicities exceeds (i.e. >) the boundary *\n")
-  cat(" ******************************************************************\n")
-  cat("\n  Low boundary:",  bdry.lo, "\n", sep="   ")
-  cat(" High boundary:",  bdry.hi, "\n", sep="   ")
-  cat("\n Operating Characteristics: \n\n")
-  bdry.oc <- round(x$bdry.oc, digits=3)
-  bdry.oc[,c(4,7)] <- round(bdry.oc[,c(4,7)], digits=1)
-  print(bdry.oc)
+  list("looks"=n, "lo.bdry"=r1, "hi.bdry"=r0, "bdry.oc"=bdry.oc, "delta"=delta)
 }
 
+# function to estimate boundary crossing probabilities
 bdrycross.prob <- function(n, r, ptox) {
   nlook <- length(n)
   np <- length(ptox)
   dn <- diff(c(0,n))
+  # pcur is the probability that the path stops at current boundary 0...r[i] without stopping earlier
   pcur <- matrix(dbinom(rep(0:r[1],np), dn[1], rep(ptox, rep(r[1]+1, np))), r[1]+1, np)
+  # pstop is the cumulative probability of stopping at current look
   pstop <- 1-pbinom(r[1], dn[1], ptox)
   ess <- n[1]*pstop
   for(i in 2:nlook) {
     ll <- r[i-1] + 1
+    # pstop0 is the cumulative probability of stopping at current look (not having stopped before)
     pstop0 <- rep(0,np)
     for(j in (r[i]-r[i-1]):min(dn[i],r[i])) {
       pstop0 <- pstop0 + pcur[ll,]*(1-pbinom(j, dn[i], ptox))
       ll <- ll - 1
     }
     pstop <- pstop + pstop0
     ess <- ess + n[i]*pstop0
+    # pnext is probability is probability it doesn't stop at next look and is at 0...r[i+1]
     pnext <- matrix(0, r[i]+1, np)
     for(j in 0:r[i]) {
       for(k in max(0, j-dn[i]):min(j,r[i-1])) {
         pnext[j+1,] <- pnext[j+1,] + pcur[k+1,]*dbinom(j-k, dn[i], ptox)
       }
     }
+    # for the next look pnext becomes the new pcur
     pcur <- pnext
   }
   ess <- ess + n[nlook]*(1-pstop)
+  # pcross is probability of stopping (crossing) all the way through
   pcross <- pstop
+  # revised pstop is probability of stopping before the final look
   pstop <- pstop - pstop0
   cbind(ptox, pcross, pstop, ess)
 }
+
+# futility boundary is for response rates instead of toxicity
+# boundary is obtained because of non-response is equivalent to toxicity
+# desirable response rate same as acceptable non-response (toxicity) rate
+# baseline response rate same as high non-response (toxicity) rate
+# P(cross bdry | high non-response rate) = cP1 -> response baseline -> 1-size
+# P(cross bdry | low non-response rate) = cP0 -> response not promising -> 1-power
+# since size under null needs to be met priority choice should be "alt"
+futilbdry <- function(rLo, rHi, n, size=0.1, power=0.9, ngrid=3, niter=10, delta=0.5) {
+  # setup toxbdry parameters
+  pLo = 1-rHi
+  pHi = 1-rLo
+  cP0 = 1-power
+  cP1 = 1-size
+  # for response rate type 1 error (size) is prioritized; so priorityNull is FALSE
+  # now call the getBoundary function
+  out = getBoundary(pLo, pHi, n, cP0, cP1, ngrid, niter, delta, priorityNull=FALSE)
+  # modify the values to suit response rates
+  out$bdry.oc = out$bdry.oc[nrow(out$bdry.oc):1,] # reverse order i.e. increasing rates of response
+  out$bdry.oc[,1] = 1 - out$bdry.oc[,1] # change non-response to response rate
+  # convert excess non-response to max response
+  out$lo.bdry = out$looks - out$lo.bdry - 1
+  out$hi.bdry = out$looks - out$hi.bdry - 1
+  # convert pcross to p(effective)
+  out$bdry.oc[,"pcross.lo"] = 1 - out$bdry.oc[,"pcross.lo"]
+  out$bdry.oc[,"pcross.hi"] = 1 - out$bdry.oc[,"pcross.hi"]
+  # rename the columns
+  colnames(out$bdry.oc)[c(1,2,5)] = c("presp", "peffective.lo", "peffective.hi")
+  class(out) <- "futilbdry"
+  out  
+}
+
+# toxicity boundary
+toxbdry <- function(pLo, pHi, n, cP0=0.1, cP1=0.9, ngrid=6, niter=10, delta=0, priority=c("null","alt")) {
+# decide whether to prioritize the type I (null) or the type II (alt) error.
+  priority <- match.arg(priority)
+  priorityNull= priority=="null"
+  # now call the getBoundary function
+  out = getBoundary(pLo, pHi, n, cP0, cP1, ngrid, niter, delta, priorityNull)
+  class(out) <- "toxbdry"
+  out
+}
+
+# print function for toxicity boundary
+print.toxbdry <- function(x, ...) {
+# convert the boundary to human readable form
+  n <- x$looks
+  if (max(diff(n)) == 1) {
+    ii <- c(which(diff(x$lo.bdry) > 0), length(n))
+    bdry.lo <- paste(x$lo.bdry[ii], n[ii], sep="/")
+    ii <- c(which(diff(x$hi.bdry) > 0), length(n))
+    bdry.hi <- paste(x$hi.bdry[ii], n[ii], sep="/")
+  } else {
+    bdry.lo <- paste(x$lo.bdry, n, sep="/")
+    bdry.hi <- paste(x$hi.bdry, n, sep="/")
+  }
+  cat("\n Toxicity boundary based on repeated significance testing \n")
+  cat("    Boundary shape parameter delta =", x$delta, "\n\n")
+  cat(" ******************************************************************\n")
+  cat(" * Stop if the number of toxicities exceeds (i.e. >) the boundary *\n")
+  cat(" ******************************************************************\n")
+  cat("\n  Low boundary:",  bdry.lo, "\n", sep="   ")
+  cat(" High boundary:",  bdry.hi, "\n", sep="   ")
+  cat("\n Operating Characteristics: \n\n")
+  bdry.oc <- round(x$bdry.oc, digits=3)
+  bdry.oc[,c(4,7)] <- round(bdry.oc[,c(4,7)], digits=1)
+  print(bdry.oc)
+}
+
+# print function for futility boundary
+print.futilbdry <- function(x, ...) {
+# convert the boundary to human readable form
+  n <- x$looks
+  if (max(diff(n)) == 1) {
+    ii <- c(which(diff(x$lo.bdry) > 0), length(n))
+    bdry.lo <- paste(x$lo.bdry[ii], n[ii], sep="/")
+    ii <- c(which(diff(x$hi.bdry) > 0), length(n))
+    bdry.hi <- paste(x$hi.bdry[ii], n[ii], sep="/")
+  } else {
+    bdry.lo <- paste(x$lo.bdry, n, sep="/")
+    bdry.hi <- paste(x$hi.bdry, n, sep="/")
+  }
+  cat("\n Futility boundary based on repeated significance testing \n")
+  cat("    Boundary shape parameter delta =", x$delta, "\n\n")
+  cat(" ******************************************************************\n")
+  cat(" * Stop if the number of responses at most (i.e. <=) the boundary *\n")
+  cat(" ******************************************************************\n")
+  cat("\n  Low boundary:",  bdry.lo, "\n", sep="   ")
+  cat(" High boundary:",  bdry.hi, "\n", sep="   ")
+  cat("\n Operating Characteristics: \n\n")
+  bdry.oc <- round(x$bdry.oc, digits=3)
+  bdry.oc[,c(4,7)] <- round(bdry.oc[,c(4,7)], digits=1)
+  print(bdry.oc)
+}

man/clinfun-internal.Rd

@@ -10,6 +10,8 @@
 compareROC.paired(x,y,d,nperm=2500)
 compareROC.unpaired(x, dx, y, dy, nperm=2500, mp=NULL)
 exactNull(gsize, kw)
+getBoundary(pLo, pHi, n, cP0=0.1, cP1=0.9, ngrid=6, niter=10, delta=0,
+    priorityNull=TRUE) 
 gsd.bdryconstant(ifrac, eprob = 0.05, delta = 0.5,
     alternative = c("two.sided", "one.sided"), tol = 0.00001, ...)
 gsd.drift.efficacy(ifrac, delta.eb, sig.level=0.05, pow = 0.8,

man/toxbdry.Rd

@@ -1,8 +1,10 @@
 \name{toxbdry}
-\title{Stopping rule for toxicity monitoring}
+\title{Stopping rule for toxicity/futility monitoring}
 \alias{toxbdry}
+\alias{futilbdry}
 \alias{bdrycross.prob}
 \alias{print.toxbdry}
+\alias{print.futilbdry}
 \keyword{design}
 \description{
   Computes a stopping rule and its operating characteristics for
@@ -11,22 +13,29 @@
 \usage{
 toxbdry(pLo, pHi, n, cP0=0.1, cP1=0.9, ngrid=6, niter=10, delta=0,
           priority=c("null","alt"))
+futilbdry(rLo, rHi, n, size=0.1, power=0.9, ngrid=3, niter=10, delta=0.5)
 bdrycross.prob(n, r, ptox)
 \method{print}{toxbdry}(x, \dots)
+\method{print}{futilbdry}(x, \dots)
 }
 \arguments{
   \item{pLo}{the toxicity rate that is acceptable.}
+  \item{rLo}{baseline (null) response rate.}
   \item{pHi}{the toxicity rate that is too high and hence unacceptable.}
-  \item{n}{vector of times (sample size) when toxicty is moniroted.}
-  \item{r}{vector of maximum acceptable toxicities corresponding to n.}
+  \item{rHi}{desirable response rate. Stop when it is too unlikely.}
+  \item{n}{vector of times (sample size) when toxicty/response is monitored.}
+  \item{r}{vector of maximum acceptable toxicities (non-responders for
+    futility) corresponding to n.}
   \item{ptox}{the toxicity rates for which the operating characteristics
-    are calculated.}
+    are calculated. For futility this is the non-response rate.}
   \item{cP0}{boundary crossing probability under pLo i.e. type I error
     or the probability of declaring a treatment with toxicity rate pLo
     unacceptable.}
   \item{cP1}{boundary crossing probability under pHi i.e. power or the
     probability of declaring a treatment with toxicity rate pHi
     unacceptable.}
+  \item{size}{probability of calling drug effective if response rate is rLo.}
+  \item{power}{probability of calling drug effective if response rate is rHi.}
   \item{ngrid}{the number of toxicity rates from pLo to pHi for which
     the operating characteristics are computed.}
   \item{niter}{the number of iterations run to obtain the boundary.}

src/deltaAUC.f

@@ -144,6 +144,7 @@ subroutine daucmats(n, p, n0, x, xb, xb0, delta, dmat, wmat, valt)
          xb(i) = xb(i)*bwratio
  102  continue
 c     calculation of variance for the CI under non-null deltaauc
+      esq = 0.0d0
       do 120 i = n0+1, n
          i1 = i-n0
          do 110 j = 1, n0
@@ -185,7 +186,7 @@ subroutine rocauc(n, nn, marker, area)
       integer n, nn
       double precision marker(n), area
 
-      integer i, j, j0, k, nties0, nties1, nties, nd
+      integer i, nties0, nties1, nties, nd
       double precision rnd, disgt
 
 c     storage for x, loc used to store ith marker and order