#' 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 T0s, A B by length(pow) T statistcs matrix based on B permutations
#' @return Pvalues of SPC and aSPC tests
#' @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 = 260 # sample size
#' p = 3000; q = 300;
#'
#' 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)
#' fit = aSPC_dcor_get_pvals(X, Y, pow = c(1:8, Inf), T0s = T_nulls)
#' fit
#' @export
#' @importFrom energy dcor
aSPC_dcor_get_pvals = function(df1, df2, pow = pow, T0s = T0s){
B = nrow(T0s)
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 )
# a = proc.time()
mat_obs = dcor_list2(ls_distCenterMat_X, ls_distCenterMat_Y,
ls_distvariance_X, ls_distvariance_Y)
# proc.time() - a
## obs statistics
T_obs = rep(NA,length(pow))
for(k in 1:length(pow)){
if(pow[k]<Inf) T_obs[k] = sum(mat_obs^pow[k]) else T_obs[k] = max(abs(mat_obs))
}
pPerm0 = rep(NA,length(pow))
for(j in 1:length(pow)){
pPerm0[j] = round((sum(abs(T_obs[j])<=abs(T0s[1:(B-1),j]))+1)/(B), digits=8)
P0s = (B-rank(abs(T0s[,j]))+1)/(B)
if (j==1) minp0=P0s else minp0[which(minp0>P0s)]=P0s[which(minp0>P0s)]
}
Paspu<-(sum(minp0<=min(pPerm0))+1)/(B+1)
pvs <- c(pPerm0, Paspu)
names(pvs) = c(paste0("SPC.",pow),"aSPC")
return(list(pvs = pvs, B = B))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.