######################################################################
# Computer notes                 Statistics for Laboratory Scientists 
# Lecture 18                                 Johns Hopkins University 
######################################################################
# These notes provide R code related to the course lectures, to
# further illustrate the topics in the lectures and to assist the
# students to learn R.
#
# Lines beginning with the symbol '#' are comments in R.  All other
# lines contain code.
#
# In R for Windows, you may wish to open this file from the menu bar
# (File:Display file); you can then easily copy commands into the
# command window.  (Use the mouse to highlight one or more lines; then
# right-click and select "Paste to console".)
######################################################################

#######################################
# Functions for doing permutation tests
#######################################
# Utility function
#     returns binary representation of 1:(2^n)
binary.v <-
function(n)
{
  x <- 1:(2^n)
  mx <- max(x)
  digits <- floor(log2(mx))
  ans <- 0:(digits-1); lx <- length(x)
  x <- matrix(rep(x,rep(digits, lx)),ncol=lx)
  x <- (x %/% 2^ans) %% 2
}

# Function to perform a paired permutation test
#     Input: differences, d
#            no. permutations
# Use paired.perm.test(d, n.perm=NULL) to 
#     do the exact test
paired.perm.test <-
function(d,n.perm=1000)
{
  n <- length(d)
  if(is.null(n.perm)) { # do exact test
    ind <- binary.v(n)
    allt <- apply(ind,2,function(x,y)
                  t.test((2*x-1)*y)$statistic,d)
  }
  else { # do n.perm samples
    allt <- 1:n.perm
    for(i in 1:n.perm) 
      allt[i] <- t.test(d*sample(c(-1,1),n,repl=TRUE))$statistic
  }
  allt
}


# Function to perform permutation test
#     x, y = the two samples
#     n.perm = number of permutations
#     var.equal = passed to the function t.test
# Use perm.test(x, y, n.perm=NULL)
#     to get exact P-value
perm.test <-
function(x, y, n.perm=1000, var.equal=TRUE)
{
  # number of data points
  kx <- length(x)
  ky <- length(y)
  n <- kx + ky

  # Data re-compiled
  X <- c(x,y)
  z <- rep(1:0,c(kx,ky))

  if(is.null(n.perm)) { # do exact permutation test
    o <- binary.v(n)  # indicator of all possible samples
    o <- o[,apply(o,2,sum)==kx]  
    nc <- choose(n,kx)
    allt <- 1:nc
    for(i in 1:nc) {
      xn <- X[o[,i]==1]
      yn <- X[o[,i]==0]
      allt[i] <- t.test(xn,yn,var.equal=var.equal)$statistic
    }
  }
  else { # do 1000 permutations of the data
    allt <- 1:n.perm
    for(i in 1:n.perm) {
      z <- sample(z)
      xn <- X[z==1]
      yn <- X[z==0]
      allt[i] <- t.test(xn,yn,var.equal=var.equal)$statistic
    }
  }
  allt
}

##############################
# paired example
##############################
# the data 
d <- c(28.6, -5.3, 13.5, -12.9, 37.3, 25.0,
       5.1, 34.6, -12.1, 9.0, 39.4)
x <- c(117.3, 100.1, 94.5, 135.5, 92.9, 118.9,
       144.8, 103.9, 103.8, 153.6, 163.1)
y <- x+d

# T-test
tobs <- t.test(d)$statistic

# Sign test
2*pbinom(sum(d < 0), length(d), 0.5)

# Signed rank test
wilcox.test(d)

# Permutation test: all possible permutations
tall <- paired.perm.test(d, n.perm=NULL)
# p-value
mean(abs(tall) >= abs(tobs))

# Permutation test: 1000 permutations
tsamp <- paired.perm.test(d, n.perm=1000)
# p-value
mean(abs(tsamp) >= abs(tobs))


##############################
# Two-sample example
##############################
x2 <- c(43.3, 57.1, 35.0, 50.0, 38.2, 61.2)
y2 <- c(51.9, 95.1, 90.0, 49.7, 101.5, 74.1, 
        84.5, 46.8, 75.1)

# t-test
tobs2 <- t.test(x2,y2, var.equal=TRUE)$statistic

# Permutation test: all possible permutations
tall2 <- perm.test(x2, y2, n.perm=NULL)
# p-value
mean(abs(tall2) >= abs(tobs2))

# Permutation test: 1000 permutations
tsamp2 <- perm.test(x2, y2, n.perm=1000)
# p-value
mean(abs(tsamp2) >= abs(tobs2))

# Wilcoxon rank-sum test
wilcox.test(x2,y2)

##################
# End of comp18.R
##################