R/meanConsTest.R

In MCARtest: Optimal Nonparametric Testing of Missing Completely at Random

#' Carry out a test of MCAR checking consistency of mean vectors.
#'
#' This is the implementation of Algorithm 2 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 mean vectors, as outlined in
#' Algorithm 2 in \insertCite{BB2024;textual}{MCARtest}.
#' @export
#'
#' @references \insertRef{BB2024}{MCARtest}
#'
#' @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(1,0), SigmaS[[1]])
#' X2 = mvrnorm(m, c(0,0), SigmaS[[2]])
#' X3 = mvrnorm(m, c(3,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)
#' 
#' meanConsTest(X, B)

meanConsTest = function(X, B){
  
  ####--------------------------------------------------------------------------
  #### compute M^0
  ####--------------------------------------------------------------------------
  
  result = get_SigmaS(X, min_diff = 10)
  d = result$ambient_dimension; N_S = result$n_S; n = sum(unlist(N_S))
  n_pattern = result$n_pattern; patterns = result$pattern
  
  mu_S = result$mu_S
  
  data_pattern = result$data_pattern
  
  ####--------------------------------------------------------------------------
  
  M_hat_0 = M(mu_S, patterns) 
  
  ####--------------------------------------------------------------------------
  #### rotate X, to make it look like it's from H0
  ####--------------------------------------------------------------------------
  
  avmu = av(mu_S, patterns)
  
  rot_data_pattern = list()
  for (i in 1:n_pattern){
    rot_data_pattern[[i]] = data_pattern[[i]] - 
      do.call(rbind, replicate(dim(data_pattern[[i]])[1], mu_S[[i]], simplify = FALSE)) +
      do.call(rbind, replicate(dim(data_pattern[[i]])[1], avmu[patterns[[i]]], simplify = FALSE))
  }
  
  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 M^b
    ####--------------------------------------------------------------------------
    
    result = get_SigmaS(X, min_diff = 10)
    mu_S_b = result$mu_S
    
    M_hat_b = M(mu_S_b, patterns)
    
    if (M_hat_b >= M_hat_0){
      sum_indicator = sum_indicator + 1
    }
  }
  
  p_hat = (1+sum_indicator)/(B+1)
  return(p_hat)
}