R/score.R

In miheerdew/BMD-metrics: Metrics to compare cbce output.

#' Assign score to a bi-community.
#'
#' Given X and Y set of a community, compute a "p-value"
#' (warning: this should be thought of as a score as it hasn't been corrected for multiple comparisons)
#' for the community being a fixed point.
#'
#' @param X,Y the X and Y matrices corresponding to the entries in the Bimodule.
#' @param method The score evalulation method to use. Must be one from c('chisq','normal')
#' @export
score <- function(X, Y, method='chisq', matrix_only=FALSE) {
    stopifnot(method %in% c('chisq', 'normal'))

    dx <- ncol(X)
    dy <- ncol(Y)
    n <- nrow(X)
    N <- n - 1

    X <- scale(X)
    Y <- scale(Y)


    X2 <- X^2
    Y2 <- Y^2

    XYs <- do.call("cbind", rlist::list.map(1:dy, X*Y[,.]))
    X2s <- matrix(rep(X2, times=dy), nrow=n)
    Y2s <- do.call("cbind", rep(rlist::list.map(1:dy, Y2[,.]), each=dx))

    #The indexing is (i, j) -> i + (j-1)*dx
    Rs <- as.vector(stats::cor(X, Y))

    S1 <- crossprod(XYs)/N
    S2 <- tcrossprod(Rs) * crossprod(X2s + Y2s)/N
    S3 <- Rs*crossprod(X2s + Y2s, XYs)/N

    #See 5.1 Steiger-Hakstian

    S <- S1 + S2/4 - (S3 + t(S3))/2

    if(matrix_only) {
      return(list(S1=S1, S2=S2, S3=S3, S=S))
    }

    if(method == "chisq"){
      Tstat <- n*sum(Rs^2)
      #Use satterthwite approx
      f1 <- sum(S^2) # trace(StS)
      f2 <- sum(diag(S))

      a <- f1/f2
      b <- f2^2/f1

      # Tstat ~ zSz ~ a * chisq(df=b)
      stats::pchisq(Tstat / a, df = b, lower.tail = FALSE)
    } else {
      Tstat <- sqrt(n)*sum(Rs)
      stats::pnorm(Tstat, sd=sqrt(sum(S)), lower.tail = FALSE)
    }
}