fisher_test.R

#' Function to compute fisher z-transformation pvalues for the links
#'
#' @param score_matrix score matrix obtained from DPM or reg.DPM.
#' @param alpha a significance level
#' @param nsamp is the number of samples (number of rows of input matrix for (reg-)DPM ). Note that we need nsamp such that: nsamp*(nsamp+1)> 2*(nvar-5)
#' @param nvar is the number of variables (number of columns of input matrix for (reg-)DPM )
#'
#' @return a list containing
#' - the list of significant links as a data frame. Each line of the data frame corresponds to a link.
#' The first and second columns are the corresponding nodes of the link, the third column is the score of the link and the forth column is the corresponding pvalue.
#' - the threshold that all the links with scores more than it are significant at defined significance level.
#'
#' @examples
#' ## Load DREAM4 data (insilico_size10_1_knockouts)
#' data(dream4)
#'
#'## Run regularized DPM
#' res <- reg.dpm(data)
#'
#' ## Get the list of significant links at significance level alpha=0.05
#' fisher_test(res,alpha=0.05,10,10)$LinkList
#'
#' ##Get the threshold at significance level alpha=0.05
#' fisher_test(score_matrix=res,alpha=0.05,nsamp=10,nvar=10)$threshold
#'
#'
#' @export

###########################################################
### Function to compute fisher z transformation pvalues ###
###########################################################

fisher_test <- function(score_matrix, alpha, nsamp, nvar){
  nsamp <- (nsamp*(nsamp+1))/2
  if (nsamp-(nvar-2)-3 <= 0 ) stop("the test requires nsamp such that: nsamp*(nsamp+1)> 2*(nvar-5)")
  diag(score_matrix) <- NA
  score_matrix[upper.tri(score_matrix)] <- NA
  linkList <- reshape2::melt(score_matrix, na.rm=TRUE)
  colnames(linkList) <- c("Node1", "Node2", "Score")
  p <- linkList["Score"]
  z <- sqrt(nsamp-(nvar-2)-3) * abs(atan(p))
  pval <- 2 * pnorm(-as.matrix(z))
  colnames(pval) <- "Pvalue"
  linkList <- cbind(linkList,pval)
  linkList <- linkList[linkList$Pvalue<alpha,]
  temp <- exp( 2 * pval_tr / (sqrt(nsamp-(nvar-2)-3)) )
  threshold <- (temp -1)/(temp + 1)
  return(list(LinkList=linkList,threshold=threshold))
}