R/I2C2.R

In neuroconductor-devel/I2C2: Image Intraclass Correlation Coefficient

#' @rdname I2C2
#' @title Image Intraclass Correlation Coefficient
#' @description Calculate image intraclass correlation
#'  coefficient (I2C2) of balanced/unbalanced data using the trace method
#'
#' @param y An n by p data matrix containing n vectorized image data with p voxels.
#' Each row contains one observed image data at a particular visit for one subject.
#' Each column contains image values for all subjects and visits at a particular voxel.
#'
#' The rows are organized by subjects and then visits, EX)
#' (Y11, Y12, Y21, Y22, ... , YI1 , YI2)
#' @param id Vector of IDs, EX) c(1, 1, 2, 2, 3, 3, 4, 4, ... , I, I)
#' @param visit Vector of visits, EX) (1, 2, 1, 2, 1, 2, ... , 1, 2)
#' @param demean if TRUE, include the demean step and
#' output the demeaned dataset
#' @param twoway a logical argument indicating whether a oneway or twoway
#' mean subtraction is more appropriate for the problem. If FALSE, only the overall sample
#' mean will be removed only; if TRUE, it will also remove visit specific means to
#' avoid scanner or batch effects
#' @param symmetric if FALSE then the function uses the
#' method of moments estimator formula;
#' if TRUE, pairwise symmetric sum formula, default is FALSE
#' @param truncate if TRUE, set negative I2C2 to zero
#' @param return_demean return de-meaned matrix.  Set to \code{FALSE}
#' @param ... not used
#'
#' @return The output of the function is a list that contains the
#' following elements.
#' lambda:       estimated I2C2
#' Kx:           the trace of between-cluster variance operator
#' Ku:           the trace of within-cluster variance operator
#' demean_y:     if demean == TRUE, output the demeaned dataset
#'
#' @author Haochang Shou, Ani Eloyan, Seonjoo Lee, Vadim Zipunnikov, Adina N. Crainiceanu,
#' Mary Beth Nebel, Brian Caffo, Martin Lindquist, Ciprian M. Crainiceanu
#' @references
#' Shou, H, Eloyan, A, Lee, S, Zipunnikov, V, Crainiceanu, AN, Nebel, NB, Caffo, B, Lindquist, MA, Crainiceanu, CM (2013).
#' Quantifying the reliability of image replication studies: the image intraclass correlation coefficient (I2C2).
#' Cogn Affect Behav Neurosci, 13, 4:714-24.
#'
#' @export
#' @examples
#' id = c(1:10, 10:1)
#' visit = rep(1:2, each = 10)
#' visit = as.character(visit)
#' n = length(id)
#' p = 100
#' y = matrix(rnorm(n * p), nrow = n, ncol = p)
#' res = I2C2(y, id = id, visit = visit)
#' res$lambda
I2C2 <- function(
  y,
  id,
  visit,
  symmetric = FALSE,
  truncate = FALSE,
  twoway = TRUE,
  demean = TRUE,
  return_demean = TRUE,
  ...){


  L = check_id_visit(y = y, id = id, visit = visit)
  rm(y); gc()
  n = L$n
  y = L$y
  id = L$id
  I = L$I
  visit = L$visit
  rm(L);
  gc(); gc();
  # p = L$p

  n_I0 = as.numeric(table(id))  # visit number for each id cluster
  k2 = sum(n_I0 ^ 2) # \sum J_i^2

  ### If demean == TRUE, we calculate the overall mean function and subtract
  ### the mean function from the data
  ### If twoway mean subtraction is needed ("twoway==TRUE"),  the visit specific
  ### mean function also computed and removed from the raw data.
  demean = as.logical(demean)
  tol = 0

  if (demean) {
    # need resd as in output
    resd = demean_matrix(y = y, visit = visit, twoway = twoway, tol = tol)
    rm(y); gc()
    W = resd
    if (!return_demean) {
      rm(resd)
    }
  } else {
    W <- y
    rm(y); gc()
  }

  # population average for the demeaned dataset W
  Wdd = colMeans(W)

  # subject-specific sum for the demeaned dataset W
  Ni = !is.na(W)
  class(Ni) = "numeric"
  Ni = rowsum(Ni, group = id)

  Si = rowsum(W, group = id)
  Si = Si / Ni

  # repeat the rows in order for the ids
  Wi = Si[id,]
  ### If symmetric is FALSE, use the method of moments estimator
  ### formula from the manuscript; otherwise, use pairwise symmetric sum estimator

  if (!symmetric) {
    trKu <- sum((W - Wi) ^ 2) / (n - I)
    trKw <- sum((t(W) - Wdd) ^ 2) / (n - 1)
    trKx <- (trKw - trKu) #/ (1 + (1 - k2 / n) / (n - 1))  #removed constant in the denominator
  } else {
    Si = Si * Ni
    trKu <- (sum(W ^ 2 * n_I0[id]) - sum(Si ^ 2)) / (k2 - n)
    trKw <-
      (sum(W ^ 2) * n - sum((n * Wdd) ^ 2) - trKu * (k2 - n)) / (n ^ 2 - k2)
    trKx <-  trKw - trKu
  }

  ## estimated I2C2 values
  lambda <-  trKx / (trKx + trKu)
  if (truncate) {
    lambda[ lambda <= 0] = 0
    ## If trun==TRUE, truncate negative lambdas to 0
  }

  ###  Return the results from I2C2 calculation as a list, with 'lambda' as I2C2 value,
  ###  Kx and Ku being the trace of between and within cluster variance operators;
  ###  If demean == TRUE, also return demeaned data
  L = list(
    lambda = lambda,
    Kx = trKx,
    Ku = trKu
  )
  if (demean & return_demean) {
    class(resd)  = c(class(resd), "demeaned_matrix")
    L$demean_y = resd
  }

  return(L)

}

#' @rdname I2C2
#' @export
#' @aliases i2c2
i2c2 = function(...){
  res = I2C2(...)
  return(res)
}