#' 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)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.