R/el.test2.R

In Compositional: Compositional Data Analysis

el.test2 <- function(y1, y2, R = 0, ncores = 1, graph = FALSE) {
  ## y1 and y2 are the two matrices containing the multivariate (or univariate data)
  ## R is the type of calibration
  ## If R = 0, the chi-square distribution is used
  ## If R = 1, the James corrected chi-square distribution is used
  ## If R = 2, the MNV F distribution is used
  ## If R > 2, bootstrap calibration is implemented
  ## the next function is for the EL
  elpa <- function(mu, y1, y2, n) {
    t1 <- emplik::el.test(y1, mu, maxit = 1000)
    t2 <- emplik::el.test(y2, mu, maxit = 1000)
    g1 <- as.numeric( t1$"-2LLR" )
    g2 <- as.numeric( t2$"-2LLR" )
    if ( round( sum(t1$wts) ) + round( sum(t2$wts) ) == n ) {
      g <- g1 + g2
    } else  g <- 1000
    g
  }
  d <- dim(y1)[2]  ## number of variables
  n1 <- dim(y1)[1]   ;   n2 <- dim(y2)[1]  ## sample sizes
  n <- n1 + n2  ## total sample size
  s1 <- (n1 - 1) / n1 * Rfast::cova(y1)   ;   s2 <- (n2 - 1) / n2 * Rfast::cova(y2)
  m1 <- Rfast::colmeans(y1)
  m2 <- Rfast::colmeans(y2)  ## mean vectors
  ## mu is the estimate of the common mean vector
  v1 <- solve(s1)    ;   v2 <- solve(s2)
  a1 <- n1 * v1    ;    a2 <-  n2 * v2
  muc <- solve( a1 + a2, a1 %*% m1 + a2 %*% m2 )
  runtime <- proc.time()
  apot <- nlm( elpa, muc, y1 = y1, y2 = y2, n = n )
  test <- apot$minimum
  mu <- apot$estimate
  runtime <- proc.time() - runtime
  if ( R == 0 ) {
    pvalue <- pchisq(test, d, lower.tail = FALSE)
    result <- list( test = test, dof = d, pvalue = pvalue, mu = mu, runtime = runtime,
                    note = paste("Chi-square approximation") )
  } else if ( R == 1 ) {
    delta <- Compositional::james(y1, y2, R = 1)$info[3]
    stat <- as.numeric( test / delta )
    pvalue <- as.numeric( pchisq(test / delta, d, lower.tail = FALSE) )
    result <- list( test = test, modif.test = stat, dof = d, pvalue = pvalue, mu = mu,
                    runtime = runtime, note = paste("James corrected chi-square approximation") )
  } else if ( R == 2 ) {
    dof <- Compositional::james(y1, y2, R = 2)$info[5]
    v <- dof + d - 1
    stat <- as.numeric( ( dof / (v * d) ) * test )
    pvalue <- as.numeric( pf(stat, d, dof, lower.tail = FALSE) )
    dof <- c(d, v - d + 1)
    names(dof) <- c("numer df", "denom df")
    result <- list( test = test, modif.test = stat, dof = dof, pvalue = pvalue,
                    mu = mu, runtime = runtime, note = paste("F approximation") )
    ## else bootstrap calibration is implemented
  } else if (R > 2) {
    ybar1 <- Rfast::colmeans(y1)
    ybar2 <- Rfast::colmeans(y2)
    z1 <- y1 - rep( ybar1 - muc, rep(n1, d) )
    z2 <- y2 - rep( ybar2 - muc, rep(n2, d) )
    if ( ncores == 1 ) {
      runtime <- proc.time()
      tb <- numeric(R)
      for ( i in 1:R ) {
        b1 <- Rfast2::Sample.int(n1, n1, replace = TRUE)
        b2 <- Rfast2::Sample.int(n2, n2, replace = TRUE)
        y1 <- z1[b1, ]    ;    y2 <- z2[b2, ]
        tb[i] <- nlm( elpa, muc, y1 = y1, y2 = y2, n = n )$minimum
      }
      runtime <- proc.time() - runtime
    } else {
      runtime <- proc.time()
      
      cl <- parallel::makeCluster(ncores)
      # Load required packages on workers
      parallel::clusterEvalQ(cl, {
        library(Rfast2)
        library(emplik)
      })
      # Export only what workers need
      parallel::clusterExport(cl, 
                             varlist = c("elpa", "muc", "z1", "z2", "n1", "n2", "n"), 
                             envir = environment())
      
      tb <- parallel::parSapply(cl, 1:R, function(i) {
        b1 <- Rfast2::Sample.int(n1, n1, replace = TRUE)
        b2 <- Rfast2::Sample.int(n2, n2, replace = TRUE)
        nlm( elpa, muc, y1 = z1[b1, ], y2 = z2[b2, ], n = n )$minimum
      })
      
      parallel::stopCluster(cl)    
      runtime <- proc.time() - runtime
    }
    pvalue <- ( sum(tb > test) + 1 ) / (R + 1)
    result <- list( test = test, pvalue = pvalue, mu = mu, runtime = runtime,
                    note = paste("Bootstrap calibration") )
    if ( graph ) {
      hist(tb, xlab = "Bootstrapped test statistic", main = " ")
      abline(v = test, lty = 2, lwd = 2)  ## The line is the test statistic
    }
  }
  result
}