R/varConsTest.R

In MCARtest: Optimal Nonparametric Testing of Missing Completely at Random

#' Carry out a test of MCAR checking consistency of variance vectors.
#'
#' This is the implementation of Algorithm 3 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 variance vectors, as outlined in
#' Algorithm 3 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(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)
#' 
#' varConsTest(X, B)

varConsTest = function(X, B){
  
  ####--------------------------------------------------------------------------
  #### rescale the data
  ####--------------------------------------------------------------------------
  
  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
  avsigma = av(result$sigma_squared_S, patterns)
  
  for (j in 1:d){
    X[,j] = X[,j]/sqrt(avsigma[j])
  }
  
  ####--------------------------------------------------------------------------
  #### compute V^0
  ####--------------------------------------------------------------------------
  result = get_SigmaS(X, min_diff = 10)
  sigma_squared_S = result$sigma_squared_S
  
  data_pattern = result$data_pattern
  
  V_hat_0 = V(sigma_squared_S, patterns) 
  
  ####--------------------------------------------------------------------------
  #### rotate X, to make it look like it's from H0
  ####--------------------------------------------------------------------------
  
  avsigma = av(sigma_squared_S, patterns)
  
  rot_data_pattern = list()
  for (i in 1:n_pattern){
    rot_data_pattern[[i]] = data_pattern[[i]]%*%
      diag(1/sqrt(sigma_squared_S[[i]]), ncol = length(sigma_squared_S[[i]]))%*% 
      diag(sqrt(avsigma[patterns[[i]]]), ncol = length(avsigma[patterns[[i]]]))
  }
  
  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,])
    
    ####--------------------------------------------------------------------------
    #### rescale data and compute V^b
    ####--------------------------------------------------------------------------
    
    result = get_SigmaS(X, min_diff = 10)
    avsigma = av(result$sigma_squared_S, result$patterns)
    
    for (j in 1:d){
      X[,j] = X[,j]/sqrt(avsigma[j])
    }
    
    result = get_SigmaS(X, min_diff = 10)
    sigma_squared_S_b = result$sigma_squared_S
    
    V_hat_b = V(sigma_squared_S_b, patterns)
    
    if (V_hat_b >= V_hat_0){
      sum_indicator = sum_indicator + 1
    }
  }
  
  p_hat = (1+sum_indicator)/(B+1)
  return(p_hat)
}