R/test_significance.R

In thom82/whassocr: An Implementation of Whitehead's Association Analysis

drill<-function(x,row,col) {
  l<-length(x)
  results<-vector()
  for(i in c(1:l)) {
    results<-append(results,x[[i]][row,col])
  }
  return(results)
}

#' Calculate significance statistics
#'
#' Calculate significance statistics based on a list returned
#' by \code{\link{randomGraphs}}
#'
#' @param x a n x n association matrix returned by \code{\link{makeAssociation}}
#' @param rgraphs a list object returned by \code{\link{randomGraphs}}
#'
#' @return A list of significance statistics:
#'    \item{pl}{A n x n matrix containing the percentage for each cell
#'      of random association graphs that have values below the value in \code{x}.
#'      This is essentially a p-value for whether \code{x_ij} is lower than expected}
#'    \item{pu}{A n x n matrix containing the percentage for each cell
#'      of random association graphs that have values above the value in \code{x}.
#'      This is essentially a p-value for whether \code{x_ij} is higher than expected}
#'    \item{S}{The mean square deviation of the association matrix \code{x}}
#'    \item{Sp}{The percentage of random association matrices with mean square deviations
#'      greater than \code{S}. This is essentially a p-value for whether the mean square
#'      deviation of \code{S} is significantly larger than random, and thus that associations
#'      overall are occuring significantly more often than random.}
#' @export


getSig<-function(x,rgraphs) {
  nr<-dim(x)[1]
  nc<-dim(x)[2]


  pu.matrices<-lapply(rgraphs$a,FUN=function(x,e){ (x<=e)*1 },e=x)
  pu<-Reduce("+",pu.matrices)/length(pu.matrices)

  pl.matrices<-lapply(rgraphs$a,FUN=function(x,e){ (x>=e)*1 },e=x)
  pl<-Reduce("+",pl.matrices)/length(pl.matrices)

  std<-apply(simplify2array(rgraphs$a),c(1,2),sd)
  S<-msdAssoc(x,e=rgraphs$expected)
  ses<-sapply(rgraphs$a,msdAssoc,e=rgraphs$expected)
  Sp<-length(ses[which(ses>S)])/length(rgraphs$a)
  results<-list(pl=pl,pu=pu,sd=std,S=S,Sp=Sp)
  return(results)
}

#' Calculate the mean square deviation of an association matrix
#'
#' This calculates the mean square deviation of an association
#' matrix returned by \code{\link{makeAssociation}}. Used internally by
#' \code{\link{getSig}}.
#'
#' @param g an association matrix generated by \code{\link{makeAssociation}}
#' @param e the expected association matrix returned by \code{\link{randomGraphs}}
#'  (The \code{expected} object in the returned list)
#'
#' @return A numeric value equal to the mean square deviation of \code{g}

msdAssoc<-function(g,e) {
  g[is.nan(g)]<-0
  o<-g-e
  o<-o^2
  s<-sum(o)
  D<-dim(g)[1]
  D<-D^2
  S<-s/D
  return(S)
}

#' Significance test for network statistics of the association matrix
#'
#' Test whether a network statistic of an association matrix is
#' significantly different than the mean of the network statistic
#' for the random graphs.
#'
#' @param assoc A n X n association matrix
#' @param rgraphs A list object returned by \code{\link{randomGraphs}}
#' @param FUN a function to calculate a network statistic. The function should return a singular
#' numeric value for which a mean and standard deviation can be calculated.
#' @param ... Additional arguments to pass to \code{FUN}
#'
#' @return A list containing the following values:
#'    \item{value}{the value of the test statistic for \code{assoc}}
#'    \item{random.value}{the mean value of the test statistic for all random graphs in \code{rgraphs}}
#'    \item{sd.value }{the standard deviation of \code{random.value}}
#'    \item{p.up}{the percentage of random graphs with values higher than \code{value}}
#'    \item{p.down}{the percentage of random graphs with values lower than \code{value}}
#' @export

test.assoc<-function(assoc,rgraphs,FUN=net.density,...) {
  xvalue<-FUN(assoc,...)
  xrand<-sapply(rgraphs$a,FUN=FUN,...)
  pscore.up<-length(which(xrand>=xvalue))/length(xrand)
  pscore.down<-length(which(xrand<=xvalue))/length(xrand)
  mean.xrand<-mean(xrand)
  sd.xrand<-sd(xrand)
  return(list(value=xvalue,random.value=mean.xrand,sd.value=sd.xrand,p.up=pscore.up,p.down=pscore.down))
}