# Note you can now also use the p.adjust function in R for doing FDR and other multiple comparisons mentioned in the talk.
#
# Here's p.adjust for working out FDR which is in R (at least in version 2.14.2 and above)
p <- c(0.021,0.001,0.017,0.041,0.005,0.036,0.042,0.023,0.07,0.1)
p.adjust(p,method="fdr",n=10)


#
# FDR implementation, adapted from Peter Watson's Stats wiki entry by Adam Wagner
# http://imaging.mrc-cbu.cam.ac.uk/statswiki/FAQ/FDR

#  fdrPValue
#    Inputs
#      pValues.....the P-values that we wish to apply the FDR correction 
#                  multiple correction to; can be a vector or matrix
#      alpha.....the size of alpha we wish to use; default=0.05
#    Outputs.....the largest P-value that is still significant with FDR applied -or-
#                halts execution if no P-values are significant
fdrPValue <- function(pValues, alpha=0.05){
  sortedPValues <- sort(pValues)  
  # Error catching - stop if there are no p-values smaller than alpha
  if (sum(sortedPValues<alpha)==0) 
    stop(paste("There no p-values smaller than alpha at", alpha, "- FDR pointless"))
  
  # Calculate the "prototypical" pvalues
  protoP <- 1:length(sortedPValues)*(alpha/length(sortedPValues))
  diff = sortedPValues-protoP

  # FDR P value corresponds to the largest observed P-value which is smaller than or equal
  # to its prototype
  # Error checking to deal with the case when no P-values are still significant
  if (sum(diff<0)) fdrP <- sortedPValues[which(diff==max(diff[diff<0]))]
  else stop("No P-Values remain significant with the application of FDR")
  
  fdrP
}

testMatrix <- matrix(c(0.021,0.001,0.017,0.041,0.005,0.036,0.042,0.023,0.07,0.1), nrow=2)
fdrPValue(testMatrix)
fdrPValue(testMatrix, alpha=0.009)

# NB: The P-values do *not* need to be in matrix form - 
#     that is simple used here for testing



