R/corrCompTest.R

In MCARtest: Optimal Nonparametric Testing of Missing Completely at Random

#' Carry out a test of MCAR checking compatibility of correlation matrices.
#'
#' This is the implementation of Algorithm 1 in \insertCite{BB2024;textual}{MCARtest}.
#'
#' @param X The dataset with incomplete data.
#' @param B The bootstrap sample \eqn{B} for the bootstrap test.
#'
#' @return The p-value of the test of MCAR based on correlation matrices, as outlined in
#' Algorithm 1 in \insertCite{BB2024;textual}{MCARtest}.
#' @export
#'
#' @references \insertRef{BB2024}{MCARtest}
#'
#' @importFrom pracma sqrtm
#'
#' @examples
#' library(MASS)
#' alpha = 0.05
#' B = 20
#' m = 500
#' 
#' SigmaS=list() #Random 2x2 correlation matrices (necessarily consistent)
#' for(j in 1:3){
#' x=runif(2,min=-1,max=1); y=runif(2,min=-1,max=1)
#' SigmaS[[j]]=cov2cor(x%*%t(x) + y%*%t(y))
#' }
#' 
#' X1 = mvrnorm(m, c(0,0), SigmaS[[1]])
#' X2 = mvrnorm(m, c(0,0), SigmaS[[2]])
#' X3 = mvrnorm(m, c(0,0), SigmaS[[3]])
#' columns = c("X1","X2","X3")
#' X = data.frame(matrix(nrow = 3*m, ncol = 3))
#' X[1:m, c("X1", "X2")] = X1
#' X[(m+1):(2*m), c("X2", "X3")] = X2
#' X[(2*m+1):(3*m), c("X1", "X3")] = X3
#' X = as.matrix(X)
#' 
#' corrCompTest(X, B)

corrCompTest = function(X, B){
  
  ####--------------------------------------------------------------------------
  #### compute R^0
  ####--------------------------------------------------------------------------
  
  result = get_SigmaS(X, min_diff = 1)
  d = result$ambient_dimension; N_S = result$n_S; n = sum(unlist(N_S))
  n_pattern = result$n_pattern; patterns = result$pattern
  
  SigmaS = result$SigmaS 
  
  data_pattern = result$data_pattern
  
  ####--------------------------------------------------------------------------
  #### compute c
  ####--------------------------------------------------------------------------
  c = 1
  for(i in 1:n_pattern){
    cand = min(eigen(SigmaS[[i]])$values)
    if (cand < c){
      c = cand
    }
  }
  
  ####--------------------------------------------------------------------------
  
  tmp = computeR(patterns, SigmaS)
  Q_hat = tmp$Sigma/(1-tmp$R)
  
  QS_hat = list()
  for (i in 1:n_pattern){
    QS_hat[[i]] = Q_hat[patterns[[i]], patterns[[i]]]
  }
  
  R_hat_0 = tmp$R
  
  ####--------------------------------------------------------------------------
  # Early return if R_hat_0 exceeds threshold
  ####--------------------------------------------------------------------------
  
  if (R_hat_0 >= 3 / 4) {
    return(as.numeric(0))
  }
  
  ####--------------------------------------------------------------------------
  #### rotate X, to make it look like it's from H0
  ####--------------------------------------------------------------------------
  rot_data_pattern = list()
  for (i in 1:n_pattern){
    rot_data_pattern[[i]] = scale(data_pattern[[i]]) %*%
      solve(sqrtm(as.matrix(SigmaS[[i]]))$B) %*%
      sqrtm(as.matrix(QS_hat[[i]]))$B
  }
  
  sum_indicator = 0
  for (b in 1:B){
    r_ind = 0
    X = data.frame(matrix(nrow = d*n, ncol = d))
    for (i in 1:n_pattern){
      n_S = dim(rot_data_pattern[[i]])[1]
      tmp_data = as.matrix(rot_data_pattern[[i]][sample(1:n_S, n_S, replace = T),])
      while(dim(unique(tmp_data))[1] <= dim(tmp_data)[2]){
        tmp_data = as.matrix(rot_data_pattern[[i]][sample(1:n_S, n_S, replace = T),])
      }
      X[(1+r_ind):(n_S+r_ind), patterns[[i]]] = tmp_data
      r_ind = r_ind + n_S
    }
    X = as.matrix(X[1:r_ind,])
    
    ####--------------------------------------------------------------------------
    #### compute R^b
    ####--------------------------------------------------------------------------
    result = get_SigmaS(X, min_diff = 1)
    SigmaS_b = result$SigmaS
    
    R_hat_b = computeR.reg(patterns, SigmaS_b, 2 / c)$R
    
    if (R_hat_b >= R_hat_0){
      sum_indicator = sum_indicator + 1
    }
  }
  
  p_hat = (1+sum_indicator)/(B+1)
  return(p_hat)
  
}