Monday, January 15, 2018

Equivalence Test for Two Means by TOST (Two One-Sided Test) Procedure


mean.equiv.tost.test <- function(x, y, margin=1, alpha=0.05) {
  ## Equivalence Test for Two Means by TOST (Two One-Sided Test) Procedure
  ## Assumption: Normality. 
  ## Reference: Understanding Equivalence and Noninferiority Testing (Walker and Nowacki)
  ## INPUT:
  ##    x: a vector of continuous scale.
  ##    y: another vector of continuous scale (to compare with x)
  ##    margin: Equivalence Margin (defaul to 1)
  ##    alpha: Significance Level.
  cl <- 1 - 2 * alpha
  a <- t.test(x, y, conf.level=cl)
  ci <- a$conf.int

  print(a)

  cat("Equivalence Test for Two Means by TOST Procedure:\n")
  cat("Equivalence Margin =", margin, "\n")

  if (ci[1] >= -margin && ci[2] <= margin) {
    cat("Equivalence Established! (Significance Level =", 100*alpha, "%)\n\n")
  } else cat("Equivalence NOT Established. (Significance Level =", 100*alpha, "%)\n\n")

  plot(0, 0, xlim=c(-2*margin, 2*margin), ylim=c(-1,1), axes=F, xlab="Difference in Efficacies", ylab="", type="n")
  axis(side=1, at=c(-2*margin, -margin, 0, margin, 2*margin), labels=c("", paste("-", margin), "0", margin, ""))
  abline(v=c(-margin, 0, margin), lty=c(2,1,2))
  lines(ci, c(0,0), lwd=3)
  lines(rep(ci[1], 2), c(-0.1, 0.1), lwd=2)
  lines(rep(ci[2], 2), c(-0.1, 0.1), lwd=2)
  lines(rep(-diff(a$estimate), 2), c(-0.1, 0.1), lwd=4)
}
if (F) {
  mean.equiv.tost.test(rnorm(20), rnorm(20, mean=0.3))
}

Equivalence Test for Proportions by TOST (Two One-Sided Test) Procedure


prop.equiv.tost.test <- function(x, n, margin=0.2, alpha=0.05) {
  ## Equivalence Test for Proportions by TOST (Two One-Sided Test) Procedure 
  ## Reference: Understanding Equivalence and Noninferiority Testing (Walker and Nowacki)
  ## INPUT:
  ##    x: a vector of two success counts,
  ##    n: a vector of two sample sizes
  ##    margin: Equivalence Margin (defaul to 20%)
  ##    alpha: Significance Level.
  cl <- 1 - 2 * alpha
  a <- prop.test(x, n, conf.level= cl)

  ## Using Score Interval for Difference of Proportions
  require(PropCIs)
  b <- diffscoreci(x[1], n[1], x[2], n[2], conf.level= cl)
  ci <- b$conf.int

  cat("\nEquivalence Test for Proportions by TOST Procedure:\n")
  cat("1st Trial =", x[1], "/", n[1], "=", round(x[1]/n[1]*100, 1), "% ,  2nd Trial =", x[2], "/", n[2], "=", round(x[2]/n[2]*100, 1), "%\n")
  pdiff <- -diff(a$estimate)            # p1 - p2
  cat("Difference in Efficacies =", 100 * round(pdiff, 3), "%\n")
  cat(100*cl, "% Confidence Interval = [", 100*round(ci[1], 3), "% ,", 100*round(ci[2], 3), "% ]\n")
  cat("Equivalence Margin =", 100*margin, "%\n")

  if (ci[1] >= -margin && ci[2] <= margin) {
    cat("Equivalence Established! (Significance Level =", 100*alpha, "%)\n\n")
  } else cat("Equivalence NOT Established. (Significance Level =", 100*alpha, "%)\n\n")

  cat("(Assumption: 1st Trial is New, 2nd one is Current (Reference), and Higher Proportion is Better.)\n")
  if (ci[1] >= -margin) {
    cat("Noninferiority Established! (Significance Level =", 100*alpha, "%)\n\n")
  } else cat("Noninferiority NOT Established. (Significance Level =", 100*alpha, "%)\n\n")


  ## Sample Size Determination (TOST) 
  p <- x/n
  c1 <- ((qnorm(0.95) + qnorm(0.8)))^2  # Type I Error = 5%, Type II Error = 20%, Power = 80%
  pd <- abs(pdiff)
  if (pd >= margin) N <- NA else N <- ceiling(c1 * (p[1]*(1-p[1])+p[2]*(1-p[2])) / (pd - margin)^2 )
  cat("Sample Size Required (using only input proportions and margin) =", N, "(Power = 80 % )\n")

  plot(0, 0, xlim=c(-2*margin, 2*margin), ylim=c(-1,1), axes=F, xlab="Difference in Efficacies", ylab="", type="n")
  axis(side=1, at=c(-2*margin, -margin, 0, margin, 2*margin), labels=c("", paste("-",
                                                                100*margin, "%"), "0",
                                                                paste(100*margin, "%"), ""))
  abline(v=c(-margin, 0, margin), lty=c(2,1,2))
  lines(ci, c(0,0), lwd=3)
  lines(rep(ci[1], 2), c(-0.1, 0.1), lwd=2)
  lines(rep(ci[2], 2), c(-0.1, 0.1), lwd=2)
  lines(rep(pdiff, 2), c(-0.1, 0.1), lwd=4)
}
if (F) {
    prop.equiv.tost.test(c(80, 90), c(100,100))
    prop.equiv.tost.test(c(143, 144), c(281,281), margin=0.12, alpha=0.025) 
}


Wednesday, January 3, 2018

Replace NA value with the last numeric value in each row


replace.NA.with.last.value.in.each.row <- function(dat)
{ 
    ## Purpose: Replace NA value with the last numeric value in each row. 
    ## Arguments:
    ##   dat: a data frame that may have NA in the tail of each row. 
    ## Return: an updated data frame with no NA values. 
    ## Author: Feiming Chen, Date:  3 Jan 2018, 11:34
    ## ________________________________________________

    last.value <- apply(dat, 1, function(x) as.numeric(tail(na.omit(x), 1))) # last numeric value in each row. 

    n <- is.na(dat)
    for (i in 1:nrow(dat)) {
        dat[i, n[i,]] <- last.value[i]
    }

    dat
}
if (F) {                                # Unit Test
    dat <- data.frame(ID=LETTERS[1:5], V1=c(1:4, NA), V2=c(6:7, NA, NA, NA), V3 = c(8, rep(NA, 4)))
    ##   ID V1 V2 V3
    ## 1  A  1  6  8
    ## 2  B  2  7 NA
    ## 3  C  3 NA NA
    ## 4  D  4 NA NA
    ## 5  E NA NA NA
    replace.NA.with.last.value.in.each.row(dat)
    ##   ID V1 V2 V3
    ## 1  A  1  6  8
    ## 2  B  2  7  7
    ## 3  C  3  3  3
    ## 4  D  4  4  4
    ## 5  E NA NA NA
}

Make a combined data frame with a source label using the list names


rbind.df.with.list.names <- function(obj)
{ 
    ## Purpose: Make a combined data frame with a source label using the list names. 
    ## Arguments:
    ##    obj: It is a list of data frames to be "rbind" together. 
    ## Return: a single combined data frame with all the sub data frames in "obj"
    ##         with an addtional column "SOURCE0" from the names of the list "obj",
    ##         and another addtional column "SOURCE1" that abbreviates "SOURCE0". 
    ## Author: Feiming Chen, Date:  2 Jan 2018, 11:00
    ## ________________________________________________

    o <- do.call(rbind, obj)
    o$SOURCE0 <- rep(names(obj), sapply(obj, nrow))
    o$SOURCE1 <- abbreviate(o$SOURCE0)
    o
}
if (F) {                                # Unit Test
    obj <- list(A = data.frame(x = 1:3), B = data.frame(x = 4:6))
    rbind.df.with.list.names(obj)
##     x SOURCE0 SOURCE1
## A.1 1       A       A
## A.2 2       A       A
## A.3 3       A       A
## B.1 4       B       B
## B.2 5       B       B
## B.3 6       B       B
}