R/pcor.R

In YilinZhang97/PC: Projection correlation

#' @param Z A numeric matrix
#' @param r the index of row
#' @return A numeric matrix
Subtractr=function(Z,r){
  n=nrow(Z);p=ncol(Z);
  Zctr =Z - matrix(1,n,1)%*%matrix(Z[r,],1,p);                        #  Z(i,:) - Z(r,:)
  Zctr=Zctr[-r,];


  ZrNorm =matrix(sqrt(rowSums(Zctr^2)),n-1,1);                          #norm|Z(i,:) - Z(r,:)|
  ZrStd =Zctr/(ZrNorm%*%matrix(1,1,p));                                 #Z(i,:) - Z(r,:)/|Z(i,:) - Z(r,:)|
  A=ZrStd%*%t(ZrStd);
  A[A>1]=1
  A[A<(-1)]=-1
  A=acos(A)
  A[is.nan(A)]=0;
  return(A)
}


#' Projection Correlatoin
#'
#' Calculate the projection correlation of two random vectors. Two vectors can be with different dimensions, but with equal sample sizes.
#' @param X A numeric matrix, n*p, each row of the matrix is i.i.d generated from one vector
#' @param Y A numeric matrix, n*q, each row of the matrix is i.i.d generated from the other vector
#' @return The projection correlation of X and Y
#' @examples X = matrix(rnorm(100*34,1),100,34)
#' @examples Y = matrix(rnorm(100*62,1),100,62)
#' @examples pcor(X,Y)
#' @seealso \code{\link{pcor.test}}
#' @references L.Zhu, K.Xu, R.Li, W.Zhong(2017). Projection correlation between two random vectors. Biometrika, Volume 104, Pages 829–843. https://doi.org/10.1093/biomet/asx043
#' @export

pcor=function(X,Y){

  n=nrow(X);p=ncol(X); q=ncol(Y); ##nrow(X)=nrow(Y)
  IN = rep(0,n);
  SA = rep(0,n);

  for (r in 1:n) {

    MX=Subtractr(X,r);
    MY=Subtractr(Y,r);
    XYn2=MX%*%MY

    diagMX<-rep(diag(MX),each=n-2)
    diagMY<-rep(diag(MY),each=n-2)
    offdiagM<-as.logical(lower.tri(MX)+upper.tri(MX))
    VecOffdiagMX<-matrix(MX,(n-1)^2,1)[offdiagM]
    VecOffdiagMY<-matrix(MY,(n-1)^2,1)[offdiagM]

    Sr3=(sum(XYn2)-sum(diag(XYn2))-sum(diagMX*VecOffdiagMY)-sum(diagMY*VecOffdiagMX))/((n-1)*(n-2)*(n-3))

    IN[r] =(sum(MX*MY)-sum(diag(MX*MY)))/((n-1)*(n-2))+((sum(MX)-sum(diag(MX)))/((n-1)*(n-2)))*((sum(MY)-sum(diag(MY)))/((n-1)*(n-2)))-2*Sr3
    SA[r]=((sum(MX)-sum(diag(MX)))/((n-1)*(n-2)))*((sum(MY)-sum(diag(MY)))/((n-1)*(n-2)));
  }
  #Indep_Index= mean(IN);
  Indep_Index= mean(IN)/(pi^2-mean(SA));
  return(Indep_Index)
}

#' Projection Correlatoin Permutation Test
#'
#' Return the test result of projection correlation test. Test whether two vectors are independnet. Two vectors can be with different dimensions, but with equal sample sizes.
#' @param X A numeric matrix, n*p, each row of the matrix is i.i.d generated from one vector
#' @param Y A numeric matrix, n*q, each row of the matrix is i.i.d generated from the other vector
#' @return \code{pcor.value} projection correlation of X and Y
#' @return \code{p.value} the p-value under the null hypothesis two vectors are independent
#' @examples X = matrix(rnorm(100*34,1),100,34)
#' @examples Y = matrix(rnorm(100*62,1),100,62)
#' @examples pcor.test(X,Y)
#' @seealso \code{\link{pcor}}
#' @references L.Zhu, K.Xu, R.Li, W.Zhong(2017). Projection correlation between two random vectors. Biometrika, Volume 104, Pages 829–843. https://doi.org/10.1093/biomet/asx043
#' @export
pcor.test=function(X,Y){
    pcor.value <- pcor(X,Y)
    n <- nrow(X)
    times <-  2000
    pcor.permu <- replicate(times, pcor(x[sample(1:n,n),],Y), simplify = TRUE)
    p.value <- 1-mean(pcor.value>pcor.permu)
    dname <- paste(deparse(substitute(X)), "and", deparse(substitute(Y)))
    method <- "Projection Correlation Permutation Test of Independence"
    rval <- list(pcor.value = pcor.value, p.value = p.value, method = method, data.name = dname)
    return(rval)
}