R/ms_trc.R

In trc: Truncated Rank Correlation

# Truncated Rank Correlation (TRC) tau for two vector observations
# Arguments
# X: first observed data vector 
# Y: second observed data vector
# m: threshold rank for the calculation of TRC tau
# Values
# tau: calculated TRC tau value

trc_tau <- function(X,Y,m=5)
{
   n <- as.integer(length(X))
   X <- as.double(X)
   Y <- as.double(Y)
   m <- as.integer(m)
   tau <- as.double(0)
   
   #dyn.load("ms_trc.dll")
   res <- .C("trc_tau",n,X,Y, m, tau=tau)
   #dyn.unload("ms_trc.dll")
   tau = res$tau
   return(tau)
}

# Kendall's tau for two vector observations (for the purpose of checking inner calculation)
# Arguments
# X: first observed data vector 
# Y: second observed data vector
# Values
# tau: calculated Kendall's tau value
k_tau <- function(X,Y)
{
  n <- as.integer(length(X))
  X <- as.double(X)
  Y <- as.double(Y)
  tau <- as.double(0)
  
  #dyn.load("ms_trc.dll")
  res <- .C("k_tau",n,X,Y, tau=tau)
  #dyn.unload("ms_trc.dll")
  tau = res$tau
  return(tau)
}

# Pearson's correlation for two vector observations (for the purpose of checking inner calculation)
# Arguments
# X: first observed data vector 
# Y: second observed data vector
# Values
# rho: calculated Pearson's correlation value
rho <- function(X,Y)
{
  n <- as.integer(length(X))
  X <- as.double(X)
  Y <- as.double(Y)
  rho <- as.double(0)
  
  #dyn.load("ms_trc.dll")
  res <- .C("rho",n,X,Y, rho=rho)
  #dyn.unload("ms_trc.dll")
  rho = res$rho
  return(rho)
}

# Procedure for the choice of m for the TRC tau
# Arguments
# X: first observed data vector 
# Y: second observed data vector
# start: lower bound of a search region for the threshold rank m
# range_m: proportion of length of X for specifying the end of the search region for m 
# span: the parameter alpha which controls the degree of smoothing in loess function
# Values
# tau: calculated TRC tau value with the chosen m value (chs_m)
# chs_m: the chosen m value
# km_tau_vec: calculated TRC tau values for the given values of m [start, floor(range_m*n)]
# km_tau_loess:  fitted values by the local regression with loess function for tau_vec 

trc_m_search <- function(X,Y,start=3,range_m=0.8,span=0.3)
{

   n <- as.integer(length(X))
   X <- as.double(X)
   Y <- as.double(Y)
   range_m <- as.double(range_m)
   mlen =  floor(range_m*n)
   hist_km_tau = as.double(rep(0,mlen))
   #dyn.load("ms_trc.dll")
   res <- .C("km_trc_search",n,X,Y,st=as.integer(start), range_m,hist_km_tau=hist_km_tau)
   #dyn.unload("ms_trc.dll")
   
   km_tau_vec = res$hist_km_tau[start:mlen]
   idx = start:mlen
   n_idx = length(idx)
   fit = loess(km_tau_vec~idx,span=span)
   sm_km_tau = fitted(fit)
   chs_m = start
   for(k in 1:(mlen-start-1))
   {
      if(sm_km_tau[k+1] - sm_km_tau[k] > 0 & sm_km_tau[k+2] - sm_km_tau[k+1] <=0)
      {
          chs_m = k + start
          break
      }
      else if(k==(mlen-start-1))
          chs_m = mlen
   } 

   tau = trc_tau(X,Y,m=chs_m)
   names(km_tau_vec) = idx
   names(sm_km_tau) = idx
   return(list(tau=tau, chs_m = chs_m, km_tau_vec=km_tau_vec,km_tau_loess=sm_km_tau))
}


# Procedure for estimating the null distribution of the TRC tau with a given m0 value
# Arguments
# X: first observed data vector 
# Y: second observed data vector
# nperm: the number of permutations to estimate the null distribution
# m: A rank threshold for the calculation of TRC tau
# seed: an initial seed for the permutation
# Values
# perm_tau: TRC tau values from the permuted samples

null_perm_m0<- function(X,Y,nperm=1000,m=5,seed=21)
{

   n <- as.integer(length(X))
   X <- as.double(X)
   Y <- as.double(Y)
   nperm <- as.integer(nperm)
   m <- as.integer(m)
   perm_tau <- as.double(rep(0,nperm))
   
   set.seed(seed)
   
   #dyn.load("ms_trc.dll")
   res <- .C("null_perm_m0",n,X,Y,nperm, m, perm=perm_tau)
   #dyn.unload("ms_trc.dll")
   perm_tau = res$perm
   return(perm_tau)
}


# Procedure for estimating the null distribution of the TRC tau with the m value chosen by the proposed rule
# Arguments
# X: first observed data vector 
# Y: second observed data vector
# nperm: the number of permutations to estimate the null distribution
# start: lower bound of a search region for the threshold rank m
# range_m: proportion of length of X for specifying the end of the search region for m 
# span: the parameter alpha which controls the degree of smoothing in loess function
# seed: a seed number for permutation
# all_m: a logical flag for returning permuted TRC tau values for all m values
# Values
# perm_trc: TRC tau values from the permuted samples with the m value chosen by the proposed rule
# hist_m: a vector of the chosen m values for permutations 
# perm_ktau: Kendall's tau values from the permuted samples
# perm_rho: Pearson's correlation values from the permuted samples
# perm_trc_all_m: a matrix of permuted TRC tau values for all m values,
#                 in which each column stores the permuted TRC tau values for corresponding m value.
null_perm <- function(X,Y,nperm=1000,start=3,range_m=0.5,span=0.5,seed=21,all_m=FALSE)
{
   n <- as.integer(length(X))
   X <- as.double(X)
   Y <- as.double(Y)
   nperm <- as.integer(nperm)
   range_m <- as.double(range_m)
   mlen = floor(range_m*n)
   perm_km_tau <- as.double(rep(0,nperm*mlen))
   perm_taum <- as.double(rep(0,nperm*mlen))
   perm_tau  = rep(0,nperm)
   perm_ktau  = rep(0,nperm)
   perm_rho  = rep(0,nperm)
   hist_m <- rep(0,nperm)
   set.seed(seed)
   #dyn.load("ms_trc.dll")
   res <- .C("null_perm",n,X,Y,nperm, as.integer(start), range_m, perm_km_tau=perm_km_tau, perm_taum=perm_taum,
             perm_ktau=perm_ktau, perm_rho=perm_rho)
   #dyn.unload("ms_trc.dll")
   
   perm_km_tau = matrix(res$perm_km_tau,nperm,mlen)
   perm_taum = matrix(res$perm_taum,nperm,mlen)
   colnames(perm_taum) = 1:mlen
   idx = start:mlen
   n_idx = length(idx)
   for(k in 1:nperm)
   {
     km_tau_vec = perm_km_tau[k,start:mlen]
     fit = loess(km_tau_vec~idx,span=span)
     sm_km_tau = fitted(fit)
     hist_m[k] = start
     for(s in 1:(mlen-start-1))
     {
       if(sm_km_tau[s+1] - sm_km_tau[s] > 0 & sm_km_tau[s+2] - sm_km_tau[s+1] <=0)
       {
         hist_m[k] = s + start
         break
       }
       else if(s==(mlen-start-1))
         hist_m[k] = mlen
     } 
     perm_tau[k] = perm_taum[k,hist_m[k]]
   }
   
   if(all_m == TRUE)
     return(list(perm_trc=perm_tau, hist_m = hist_m, perm_ktau=res$perm_ktau,
                 perm_rho=res$perm_rho, perm_trc_all_m = perm_taum[,idx]))
   
   return(list(perm_trc=perm_tau, hist_m = hist_m,
               perm_ktau=res$perm_ktau, perm_rho=res$perm_rho))
}

# Procedure for calculating p-values of Pearson's rho, Kendall's tau, TRC tau 
# for two-sided test for the null hypothesis correaltion is equal to 0
# based on the estimated null distribution by permutation
# Arguments
# X: first observed data vector 
# Y: second observed data vector
# nperm: the number of permutations to estimate the null distribution
# start: lower bound of a search region for the threshold rank m
# range_m: proportion of length of X for specifying the end of the search region for m 
# span: the parameter alpha which controls the degree of smoothing in loess function
# seed: a seed number for permutation
# m0: a specific m0 value for p-value of the TRC tau with m0 (defalut: NULL (not reported))
# Values
# measure:
# a vector of calculated Pearson's rho, Kendall's tau, and TRC tau with m chosen by the proposed rule if m0 = NULL
# a vector of calculated Pearson's rho, Kendall's tau, TRC tau with m0, TRC tau with m chosen by the proposed rule
# if m0 is specified.
# p_val: 
# a vector of p-values for Pearson's rho, Kendall's tau, and TRC tau with m chosen by the proposed rule if m0 = NULL
# a vector of p-values for Pearson's rho, Kendall's tau, TRC tau with m0, TRC tau with m chosen by the proposed rule
# if m0 is specified.
# chs_m: the chosen m value by the proposed procedure
# mean_perm_trc: a mean value of the estimated null distribution of TRC tau by permutation
trc_cor_test <- function(X,Y, nperm=10000,start=3,range_m=0.8,span=0.5,seed=21, m0=NULL)
{

  rho_val = rho(X,Y)
  ktau_val = k_tau(X,Y)
  temp = trc_m_search(X,Y,start,range_m,span)
  trc_val = temp$tau
  chs_m = temp$chs_m
  
  measure = c(rho=rho_val, ktau=ktau_val, trc=trc_val)
  
  if(is.null(m0))
  {
    out = null_perm(X,Y,nperm,start,range_m, span, seed, all_m=FALSE)
  } else {
    trc_m0 = trc_tau(X,Y,m=m0)
    measure = c(rho=rho_val, ktau=ktau_val, trc_m0=trc_m0, trc=trc_val)
    out = null_perm(X,Y,nperm,start,range_m, span, seed, all_m=TRUE)
  }
  
  p_rho = mean(abs(out$perm_rho)>abs(rho_val))
  p_ktau = mean(abs(out$perm_ktau)>abs(ktau_val))
  mean_trc_perm = mean(out$perm_trc)
  p_trc = mean(abs(out$perm_trc-mean_trc_perm)>abs(trc_val-mean_trc_perm))

  p_val = c(rho=p_rho, ktau=p_ktau, trc=p_trc)
  if(!is.null(m0))
  {
    indx = which(colnames(out$perm_trc_all_m) == as.character(m0))
    perm_trc_m0 = out$perm_trc_all_m[,indx]
    p_trc_m0 = mean(abs(perm_trc_m0)>abs(trc_m0))
    p_val = c(rho=p_rho, ktau=p_ktau, trc_m0=p_trc_m0, trc=p_trc)
  }
  
  return(list(measure=measure, p_val=p_val, chs_m=chs_m, mean_perm_trc = mean_trc_perm))
}