R/aSPC_dcor_getTstat.R

In jasonzyx/aSPC: An Adaptive Sum of Powered Correlation Test (aSPC) to Test Global Association Between Two Groups of Variables

#' An Adaptive Sum of Powered Correlation Test (aSPC) with dcor
#'
#' @param df1, first matrix
#' @param df2, second matrix
#' @param pow, power integer candidates, default c(1:8, Inf)
#' @param B, number of permutations to calculate a P-value
#' @return A B by length(pow) T statistcs matrix based on B permutations
#' @references Xu Z., Pan W. An adaptive and powerful test for two groups of variables with high dimension
#' @examples
#' library(mvtnorm)
#' sigma = diag(0.6, 10) + 0.4
#' n = 200 # sample size
#' p = 5; q = 5;
#'
#' Z = rmvnorm(n=n, mean=rep(0,10), sigma=sigma)
#' X = rmvnorm(n=n, mean=rep(0,p), sigma=diag(1, p))
#' Y = rmvnorm(n=n, mean=rep(0,q), sigma=diag(1, q))
#' X = as.data.frame(cbind(Z[,1:5], X))
#' Y = as.data.frame(cbind(Z[,6:10], Y))
#' dim(X)
#' dim(Y)
#' set.seed(123) # to ensure we can replicate the permutation P-value
#' T_nulls = aSPC_dcor_get_Tstat(X, Y, pow = c(1:8, Inf), B = 100, show_b = T)
#'
#'
#'
#' @export
#' @importFrom energy dcor


aSPC_dcor_get_Tstat = function(df1, df2, pow = pow, B = B, show_b = FALSE){
  # test parameter
  # df1 = X; df2 = Y; pow = c(1:8, Inf); B=2


  X = scale(df1)
  Y = scale(df2)

  n = nrow(X) ## number of subjects
  ### X and Y has to be standardize before input
  ls_X = lapply(1:ncol(X), function(x) get_doublyCenterDist(X[,x]))
  ls_distCenterMat_X = lapply(1:ncol(X), function(x) ls_X[[x]]$mat )
  ls_distvariance_X = lapply(1:ncol(X), function(x) ls_X[[x]]$dVar )

  ls_Y = lapply(1:ncol(Y), function(x) get_doublyCenterDist(Y[,x]))
  ls_distCenterMat_Y = lapply(1:ncol(Y), function(x) ls_Y[[x]]$mat )
  ls_distvariance_Y = lapply(1:ncol(Y), function(x) ls_Y[[x]]$dVar )



  T0s = matrix(nrow = B, ncol = length(pow))
  for(i in 1: B){
    index = sample(nrow(X))
    ls_distCenterMat_X_sample = lapply(1:length(ls_distCenterMat_X), function(x) ls_distCenterMat_X[[x]][index, index])
    mat0 = dcor_list2(ls_distCenterMat_X_sample, ls_distCenterMat_Y,
                      ls_distvariance_X, ls_distvariance_Y)
    for(j in 1:length(pow)){
      if(pow[j] < Inf) T0s[i,j] = sum(mat0^pow[j]) else T0s[i,j] = max(abs(mat0))
    }
    if(show_b == T) print(i)
  }


  return(T0s)


}