R/gsBinomial.R

In gsDesign: Group Sequential Design

##################################################################################
#  Binomial functionality for the gsDesign package
#
#  Exported Functions:
#                   
#    ciBinomial
#    nBinomial
#    simBinomial
#    testBinomial
#
#  Hidden Functions:
#
#    bpdiff
#
#  Author(s): Keaven Anderson, PhD.
# 
#  Reviewer(s): REvolution Computing 19DEC2008 v.2.0 - William Constantine, Ph.D., Kellie Wills 
#
#  R Version: 2.7.2
#
##################################################################################

###
# Exported Functions
###

"ciBinomial" <- function(x1, x2, n1, n2, alpha=.05, adj=0, scale="Difference")
{
    # check input arguments
    checkScalar(n1, "integer", c(1, Inf))
    checkScalar(n2, "integer", c(1, Inf))
    checkScalar(x1, "integer", c(0, n1))
    checkScalar(x2, "integer", c(0, n2))
    checkScalar(alpha, "numeric", c(0, 1), c(FALSE, FALSE))    
    checkScalar(adj, "integer", c(0, 1))
    checkScalar(scale, "character")
    scale <- match.arg(tolower(scale), c("difference", "rr", "or"))
    checkLengths(x1, x2)
    
    if (scale == "difference")
    {    
        delta <- x1 / n1 - x2 / n2
        
        if (delta == -1)
        {
            lower <- -1
        }
        else if (testBinomial(delta0=-.9999, x1=x1, x2=x2, n1=n1, n2=n2, adj=adj)
                < qnorm(alpha / 2))
        {
            lower <- -1
        }
        else
        {
            lower <- uniroot(bpdiff, interval=c(-.9999, delta), x1=x1, x2=x2,
                    n1=n1, n2=n2, adj=adj, alpha=alpha)$root
        }
        
        if (delta == 1)
        {
            upper <- 1
        }
        else if (testBinomial(delta0=.9999, x1=x1, x2=x2, n1=n1, n2=n2, adj=adj)
                > -qnorm(alpha/2))
        {
            upper<-1
        }
        else
        {
            upper <- uniroot(bpdiff, interval=c(delta,.9999), lower.tail=TRUE,
                    x1=x1, x2=x2, n1=n1, n2=n2, adj=adj,
                    alpha=alpha)$root
        }
    }
    else if (scale=="rr")
    {
        if (x1 == 0) lower <- 0
        else
        {
            lower <- uniroot(bpdiff, interval=c(-20, 20), x1=x1, x2=x2,
                    n1=n1, n2=n2, adj=adj, scale="RR", alpha=alpha)$root
            lower <- exp(lower)
        }
        if (x2 == 0) upper <- Inf
        else
        {
            upper <- uniroot(bpdiff, interval=c(-20, 20), lower.tail=TRUE,
                    x1=x1, x2=x2, n1=n1, n2=n2, adj=adj, scale="RR",
                    alpha=alpha)$root
            upper <- exp(upper)
        }
    }
    else
    {   
        if (x1 == 0 || x2 == n2) lower <- -Inf
        else
        {
            lower <- uniroot(bpdiff, interval=c(-10, 10), x1=x1, x2=x2,
                    n1=n1, n2=n2, adj=adj, scale=scale, alpha=alpha)$root
            lower <- exp(lower)
        }
        if (x2 == 0 || x1 == n1) upper <- Inf
        else
        {
            upper <- uniroot(bpdiff, interval=c(-10, 10), lower.tail=TRUE,
                    x1=x1, x2=x2, n1=n1, n2=n2, adj=adj, scale=scale,
                    alpha=alpha)$root
            upper <- exp(upper)
        }
    }
    data.frame(lower=lower,upper=upper)
}

"nBinomial" <- function(p1, p2, alpha = 0.025, beta = 0.1, delta0 = 0, ratio = 1, 
                       sided = 1, outtype = 1, scale = "Difference", n = NULL) 
{
  checkVector(p1, "numeric", c(0, 1), c(FALSE, FALSE))
  checkVector(p2, "numeric", c(0, 1), c(FALSE, FALSE))
  checkScalar(sided, "integer", c(1, 2))
  checkScalar(alpha, "numeric", c(0, 1/sided), c(FALSE, FALSE))
  checkVector(beta, "numeric", c(0, 1 - alpha/sided), c(FALSE, 
                                                        FALSE))
  checkVector(delta0, "numeric")
  checkVector(ratio, "numeric", c(0, Inf), c(FALSE, FALSE))
  checkScalar(outtype, "integer", c(1, 3))
  checkScalar(scale, "character")
  if (!is.null(n)) 
    checkVector(n, "numeric")
  scale <- match.arg(tolower(scale), c("difference", "rr", 
                                       "or", "lnor"))
  if (is.null(n)) 
    checkLengths(p1, p2, beta, delta0, ratio, allowSingle = TRUE)
  else checkLengths(n, p1, p2, delta0, ratio, allowSingle = TRUE)
  len <- max(sapply(list(p1, p2, beta, delta0, ratio), length))
  if (len > 1) {
    if (length(p1) == 1) 
      p1 <- array(p1, len)
    if (length(p2) == 1) 
      p2 <- array(p2, len)
    if (length(alpha) == 1) 
      alpha <- array(alpha, len)
    if (length(beta) == 1) 
      beta <- array(beta, len)
    if (length(delta0) == 1) 
      delta0 <- array(delta0, len)
    if (length(ratio) == 1) 
      ratio <- array(ratio, len)
  }
  if (max(delta0 == 0) > 0 && max(p1[delta0 == 0] == p2[delta0 == 
                                                          0]) > 0) {
    stop("p1 may not equal p2 when delta0 is zero")
  }
  z.beta <- qnorm(1 - beta)
  sided[sided != 2] <- 1
  z.alpha <- qnorm(1 - alpha/sided)
  d0 <- (delta0 == 0)
  if (scale == "difference") {
    if (min(abs(p1 - p2 - delta0)) < 1e-11) {
      stop("p1 - p2 may not equal delta0 when scale is \"Difference\"")
    }
    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) + (v == 0)) * 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
    p10[d0] <- (p1[d0] + ratio[d0] * p2[d0])/(1 + ratio[d0])
    p20[d0] <- p10[d0]
    sigma0 <- sqrt((p10 * (1 - p10) + p20 * (1 - p20)/ratio) * 
                     (ratio + 1))
    sigma1 <- sqrt((p1 * (1 - p1) + p2 * (1 - p2)/ratio) * 
                     (ratio + 1))
    if (is.null(n)) {
      n <- ((z.alpha * sigma0 + z.beta * sigma1)/(p1 - p2 - delta0))^2
      if (outtype == 2) {
        return(data.frame(cbind(n1 = n/(ratio + 1), n2 = ratio * 
                                  n/(ratio + 1))))
      }
      else if (outtype == 3) {
        return(data.frame(cbind(n = n, n1 = n/(ratio + 1), 
                                n2 = ratio * n/(ratio + 1), alpha = alpha, 
                                sided = sided, beta = beta, Power = 1 - beta, 
                                sigma0 = sigma0, sigma1 = sigma1, p1 = p1, 
                                p2 = p2, delta0 = delta0, p10 = p10, p20 = p20)))
      }
      else return(n = n)
    }
    else {
      pwr <- pnorm(-(qnorm(1 - alpha/sided) - sqrt(n) * 
                       ((p1 - p2 - delta0)/sigma0)) * sigma0/sigma1)
      if (outtype == 2) {
        return(data.frame(cbind(n1 = n/(ratio + 1), n2 = ratio * 
                                  n/(ratio + 1), Power = pwr)))
      }
      else if (outtype == 3) {
        return(data.frame(cbind(n = n, n1 = n/(ratio + 1), 
                                n2 = ratio * n/(ratio + 1), alpha = alpha, 
                                sided = sided, beta = 1 - pwr, Power = pwr, 
                                sigma0 = sigma0, sigma1 = sigma1, p1 = p1, 
                                p2 = p2, delta0 = delta0, p10 = p10, p20 = p20)))
      }
      else return(Power = pwr)
    }
  }
  else if (scale == "rr") {
    RR <- exp(delta0)
    if (min(abs(p1/p2 - RR)) < 1e-07) {
      stop("p1/p2 may not equal exp(delta0) when scale=\"RR\"")
    }
    a <- (1 + ratio)
    b <- -(RR * (1 + ratio * p2) + ratio + p1)
    c <- RR * (p1 + ratio * p2)
    p10 <- (-b - sqrt(b^2 - 4 * a * c))/2/a
    p20 <- p10/RR
    p10[d0] <- (p1[d0] + ratio[d0] * p2[d0])/(1 + ratio[d0])
    p20[d0] <- p10[d0]
    sigma0 <- sqrt((ratio + 1) * (p10 * (1 - p10) + RR^2 * 
                                    p20 * (1 - p20)/ratio))
    sigma1 <- sqrt((ratio + 1) * (p1 * (1 - p1) + RR^2 * 
                                    p2 * (1 - p2)/ratio))
    if (is.null(n)) {
      n <- ((z.alpha * sigma0 + z.beta * sigma1)/(p1 - p2 * RR))^2
      if (outtype == 2) {
        return(data.frame(cbind(n1 = n/(ratio + 1),
                                n2 = ratio * n/(ratio + 1))))
      }
      else if (outtype == 3) {
        return(data.frame(cbind(n = n, n1 = n/(ratio + 1), 
                                n2 = ratio * n/(ratio + 1), alpha = alpha, 
                                sided = sided, beta = beta, Power = 1-beta, 
                                sigma0 = sigma0, sigma1 = sigma1, p1 = p1, 
                                p2 = p2, delta0 = delta0, p10 = p10, p20 = p20)))
      }
      else return(n = n)
    }
    else {
      pwr <- pnorm(-(qnorm(1 - alpha/sided) - sqrt(n) * 
                       ((p1 - p2 * RR)/sigma0)) * sigma0/sigma1)
      if (outtype == 2) {
        return(data.frame(cbind(n1 = n/(ratio + 1), n2 = ratio * 
                                  n/(ratio + 1), Power = pwr)))
      }
      else if (outtype == 3) {
        return(data.frame(cbind(n = n, n1 = n/(ratio + 1), 
                                n2 = ratio * n/(ratio + 1), alpha = alpha, 
                                sided = sided, beta = 1 - pwr, Power = pwr, 
                                sigma0 = sigma0, sigma1 = sigma1, p1 = p1, 
                                p2 = p2, delta0 = delta0, p10 = p10, p20 = p20)))
      }
      else return(Power = pwr)
    }
  }
  else {
    OR <- exp(-delta0)
    if (min(abs(p1/(1 - p1)/p2 * (1 - p2) * OR) - 1) < 1e-07) {
      stop("p1/(1-p1)/p2*(1-p2) may not equal exp(delta0) when scale=\"OR\"")
    }
    a <- OR - 1
    b <- 1 + ratio * OR + (1 - OR) * (ratio * p2 + p1)
    c <- -(ratio * p2 + p1)
    p10 <- (-b + sqrt(b^2 - 4 * a * c))/2/a
    p20 <- OR * p10/(1 + p10 * (OR - 1))
    p10[d0] <- (p1[d0] + ratio[d0] * p2[d0])/(1 + ratio[d0])
    p20[d0] <- p10[d0]
    sigma0 <- sqrt((ratio + 1) * (1/p10/(1 - p10) + 1/p20/(1 - p20)/ratio))
    sigma1 <- sqrt((ratio + 1) * (1/p1/(1 - p1) + 1/p2/(1 - p2)/ratio))
    if (is.null(n)) {
      n <- ((z.alpha * sigma0 + z.beta * sigma1)/log(OR/p2 * (1 - p2) * p1/(1 - p1)))^2
      if (outtype == 2) {
        return(data.frame(cbind(n1 = n/(ratio + 1), n2 = ratio * 
                                  n/(ratio + 1))))
      }
      else if (outtype == 3) {
        return(data.frame(cbind(n = n, n1 = n/(ratio + 1), 
                                n2 = ratio * n/(ratio + 1), alpha = alpha, 
                                sided = sided, beta = beta, Power = 1-beta, 
                                sigma0 = sigma0, sigma1 = sigma1, p1 = p1, 
                                p2 = p2, delta0 = delta0, p10 = p10, p20 = p20)))
      }
      else {
        return(n = n)
      }
    }
    else {
      pwr <- pnorm(-(qnorm(1 - alpha/sided) - sqrt(n) * 
                       (log(OR/p2 * (1 - p2) * p1/(1 - p1))/sigma0)) * 
                     sigma0/sigma1)
      if (outtype == 2) {
        return(data.frame(cbind(n1 = n/(ratio + 1), n2 = ratio * 
                                  n/(ratio + 1), Power = pwr)))
      }
      else if (outtype == 3) {
        return(data.frame(cbind(n = n, n1 = n/(ratio + 1), 
                                n2 = ratio * n/(ratio + 1), alpha = alpha, 
                                sided = sided, beta = 1 - pwr, Power = pwr, 
                                sigma0 = sigma0, sigma1 = sigma1, p1 = p1, 
                                p2 = p2, delta0 = delta0, p10 = p10, p20 = p20)))
      }
      else return(Power = pwr)
    }
  }
}
"simBinomial" <- function(p1, p2, n1, n2, delta0=0, nsim=10000, chisq=0, adj=0, 
        scale="Difference")
{
    # check input arguments 
    # Note: delta0, chisq, adj, and scale checked by testBinomial() function
    checkVector(p1, "numeric", c(0, 1), c(FALSE, FALSE))
    checkVector(p2, "numeric", c(0, 1), c(FALSE, FALSE))    
    checkScalar(n1, "integer", c(1, Inf))
    checkScalar(n2, "integer", c(1, Inf))
    checkScalar(nsim, "integer", c(1, Inf))
    checkLengths(p1, p2)
    
    x1 <- rbinom(prob=p1, size=n1, n=nsim)
    x2 <- rbinom(prob=p2, size=n2, n=nsim)
    scale <- match.arg(tolower(scale), c("difference", "rr", "or", "lnor"))
    
    testBinomial(x1=x1, x2=x2, n1=n1, n2=n2, delta0=delta0, adj=adj,
                    chisq=chisq, scale=scale)
}

"testBinomial" <- function(x1, x2, n1, n2, delta0=0, chisq=0, adj=0,
        scale="Difference", tol=.1e-10)
{
    # check input arguments
    checkVector(n1, "integer", c(1, Inf))
    checkVector(n2, "integer", c(1, Inf))
    checkVector(x1, "integer", c(0, Inf))
    checkVector(x2, "integer", c(0, Inf))
    checkVector(chisq, "integer", c(0, 1))
    checkVector(adj, "integer", c(0, 1))
    checkScalar(scale, "character")
    scale <- match.arg(tolower(scale), c("difference", "rr", "or", "lnor"))
    checkScalar(tol, "numeric", c(0, Inf), c(FALSE, TRUE))
    checkLengths(n1, n2, x1, x2, delta0, chisq, adj, allowSingle=TRUE)
    checkVector(n1 - x1, "integer", c(0, Inf))
    checkVector(n2 - x2, "integer", c(0, Inf))
 
    # make all vector arguments the same length (don't extend n1, n2)
    len<-max(sapply(list(n1, n2, x1, x2, chisq, adj, delta0),length))

    if (len > 1)
    {   if (length(x1) == 1) p1<-array(x1,len)
        if (length(x2) == 1) p2<-array(x2,len)
        if (length(chisq) == 1) alpha<-array(chisq,len)
        if (length(adj) == 1) beta<-array(adj,len)
        if (length(delta0) == 1) delta0<-array(delta0,len)
    }
   
    ntot <- n1 + n2
    xtot <- x1 + x2
    
    # risk difference test - from Miettinen and Nurminen eqn (9)
    if (scale == "difference")
    {   
        checkVector(delta0, "numeric", c(-1, 1), c(FALSE, FALSE))
        L2 <- (n1 + 2 * n2) * delta0 - ntot - xtot
        L1 <- (n2 * delta0 - ntot - 2 * x2) * delta0 + xtot
        L0 <- x2 * delta0 * (1 - delta0)
        q  <- (L2 / (3 * ntot)) ^ 3 - L1 * L2 / 6 / ntot ^ 2 + L0 / 2 / ntot
        p  <- (sign(q) + (q==0)) * sqrt((L2 / (3 * ntot)) ^ 2 - L1 / (3 * ntot))
        a  <- q / p ^ 3
        a[a > 1] <- 1
        a  <- (pi + acos(a)) / 3
        R0 <- 2 * p * cos(a) - L2 / 3 / ntot
        R1 <- R0 + delta0
        V  <- (R1 * (1 - R1) / n1 + R0 * (1 - R0) / n2)
        # V=0 only if no or all events and delta0=0
        # in which case test statistic will be 0 
        V[V <= tol] <- 1
        # compute z-statistic
        z <- x1 / n1 - x2 / n2 - delta0
    }
    # relative risk test - from Miettinen and Nurminen eqn (10)
    else if (scale == "rr")
    {   
        checkVector(delta0, "numeric", c(-Inf, Inf), c(FALSE, FALSE))
        delta0 <- exp(delta0) # value of 0 input represents equal rates
        A  <- delta0 * ntot
        B  <- -(n1 * delta0 + x1 + n2 + x2 * delta0)
        C  <- xtot
        R0 <- ( - B - sqrt(B ^ 2 - 4 * A * C)) / 2 / A
        R1 <- R0 * delta0
        V  <- R1 * (1 - R1) / n1 + delta0 ^ 2 * R0 * (1 - R0) / n2
        # V=0 only if no events or (all events and delta0=1)
        # in which case test statistic will be 0
        V[V <= 0 | is.na(V)] <- 1
        z <- x1 / n1 - x2 / n2 * delta0
    }
    # odds-ratio and log-odds-ratio
    else
    {   
        checkVector(delta0, "numeric", c(-Inf, Inf), c(FALSE, FALSE))
        delta0 <- exp(delta0) # change from log odds-ratio to odds-ratio
        A <- n2 * (delta0 - 1)
        B <- n1 * delta0 + n2 - xtot * (delta0 - 1)
        C <- -xtot
        R0 <- ( - B + sqrt(B ^ 2 - 4 * A * C)) / 2 / A
        R1 <- R0 * delta0 / (1 + R0 * (delta0 - 1))
        # next 2 lines deal with length(delta0)>length(xtot) or length(ntot)
        xtem <- (delta0 == delta0) * xtot
        ntem<-(delta0==delta0)*ntot
        R0[delta0 == 1] <- xtem[delta0 == 1] / ntem[delta0 == 1]
        R1[delta0 == 1] <- R0[delta0 == 1]
        
        # odds-ratio test - from Miettinen and Nurminen eqn (13)
        if (scale == "or")
        {   
            V <- 1 / (1 / n1 / R1 / (1-R1) + 1 / n2 / R0 / (1 - R0))
            V[xtot == 0 | xtot == ntot] <- 1
            z <- n1 * (x1 / n1 - R1)
        }
        # log-odds ratio - based on asymptotic distribution of log-odds
        # see vignette
        else if (scale == "lnor")
        {   
            V <- 1 / n1 / R1 / (1-R1) + 1 / n2 / R0 / (1-R0)
            V[xtot == 0 | xtot == ntot] <- 1
            z <- log(x1 / (n1 - x1) / x2 * (n2 - x2) / delta0)
            z[xtot == 0 | xtot == ntot] <- 0
            z[xtot > 0 & xtot < ntot] <- z[xtot > 0 & xtot < ntot]
        }   
    }
    
    # do continuity correction, where required
    one <- array(TRUE, max(length(V), length(adj)))
    adj <- (adj==1) & one
    V[adj] <- V[adj] * ntot / (ntot - 1)
    
    # square z where required to get chi-square
    one <- array(TRUE, max(length(z), length(chisq)))
    chisq <- chisq & one
    if (length(z) == 1) z <- z * one
    z[chisq] <- z[chisq] ^ 2 / V[chisq]
    z[ ! chisq] <- z[ ! chisq] / sqrt(V[ ! chisq])
    
    z
}

###
# Hidden Functions
###

"bpdiff" <- function(delta, x1, x2, n1, n2, alpha=.05, adj=0, scale="Difference",
        lower.tail=FALSE)
{
    pnorm(testBinomial(x1, x2, n1, n2, delta0=delta, adj=adj, scale=scale),
            lower.tail=lower.tail) - alpha / 2
}