R/sparseRegression.R

In neuroconductor-devel/ANTsR: ANTs in R: Quantification Tools for Biomedical Images

#' Sparse regression on input images.
#'
#' Compute a sparse, spatially coherent regression from a set of input images
#' (with mask) to an outcome variable.
#'
#' @param inmatrix Input data matrix, with dimension number of subjects by
#' number of voxels.
#' @param demog Input demographics data frame.  Contains outcome variable to
#' regress against.
#' @param outcome Name of column in \code{demog} to regress against.
#' @param mask Mask for reconstructing \code{inmatrix} in physical space.
#' @param sparseness Level of sparsity desired.  \code{0.05}, for example,
#' makes 5\% of the matrix be non-zero.
#' @param nvecs Number of eigenvectors to return.
#' @param its Number of cross-validation folds to run.
#' @param cthresh Cluster threshold.
#' @param statdir Where to put results.  If not provided, a temp directory is
#' created.
#' @param z Row (subject-wise) sparseness.
#' @param smooth Amount of smoothing.
#' @return A list of values:
#' \item{eigenanatomyimages}{Coefficient vector images.}
#' \item{umatrix}{Projections of input images on the sparse regression vectors.  Can be used for, e.g., subsequent classification/predictions.}
#' \item{projections}{Predicted values of outcome variable.}
#' @author Kandel BM, Avants BB.
#' @references Kandel B.M., D. Wolk, J. Gee, and B. Avants. Predicting
#' Cognitive Data from Medical Images Using Sparse Linear Regression.
#' Information Processing in Medical Imaging, 2013.  %% ~put references to the
#' literature/web site here ~
#' @examples
#'
#' nsubj <- 50
#' prop.train <- 1/2
#' subj.train <- sample(1:nsubj, prop.train*nsubj, replace = FALSE )
#' input <- t(replicate(nsubj, rnorm(125)))
#' outcome <- seq(1, 5, length.out=nsubj)
#' demog <- data.frame(outcome=outcome)
#' input[, 40:60] <- 30 + outcome + rnorm(length(input[, 40:60]), sd=2)
#' input.train <- input[subj.train, ]
#' input.test <- input[-subj.train, ]
#' demog.train <- data.frame(outcome=demog[subj.train, ])
#' demog.test <- data.frame(outcome=demog[-subj.train, ])
#' mymask <- as.antsImage(array(rep(1, 125), dim=c(5,5,5)))
#' myregression <- sparseRegression(input.train, demog.train, 'outcome', mymask,
#'     sparseness=0.05, nvecs=5, its=3, cthresh=250)
#' # visualization of results
#' sample <- rep(0, 125)
#' sample[40:60] <-1
#' signal.img <- as.antsImage(array(rep(0,125), dim=c(5, 5, 5)))
#' signal.img[signal.img >= 0 ] <- sample
#' # plot( signal.img, axis=2, slices='1x5x1') # actual source of signal
#' # compare against first learned regression vector
#' myimgs <- list()
#' for( i in 1:5){
#'   myarray <- as.array(myregression$eigenanatomyimages[[ i ]])
#'   myarray <- myarray / max(abs(myarray)) # normalize for visualization
#'   myimgs[[ i ]] <- antsImageClone(myregression$eigenanatomyimages[[ i ]])
#'   myimgs[[ i ]][mymask > 0] <- myarray
#' }
#' # plot(myimgs[[1]], axis=2, slices='1x5x1')
#' # use learned eigenvectors for prediction
#' result <- regressProjections(input.train, input.test, demog.train,
#'     demog.test, myregression$eigenanatomyimages, mymask, 'outcome')
#' # plot(result$outcome.comparison$real, result$outcome.comparison$predicted)
#'
#' @export sparseRegression
sparseRegression <- function(inmatrix, demog, outcome, mask = NULL, sparseness = 0.05,
  nvecs = 10, its = 5, cthresh = 250, statdir = NA, z = 0, smooth = 0) {
  if (missing(inmatrix))
    stop("Missing input image matrix.")
  if (missing(demog))
    stop("Missing demographics.")
  if (missing(outcome))
    stop("Missing outcome.")
  if (is.na(statdir))
    statdir <- paste(tempdir(), "/", sep = "")
  outfn <- paste(statdir, "spca.nii.gz", sep = "")
  decomp <- paste(statdir, "spcaprojectionsView1vec.csv", sep = "")
  matname <- paste(statdir, "spcamatrix.mha", sep = "")
  antsImageWrite(as.antsImage(inmatrix), matname)
  demog.name <- paste(statdir, "demog.csv", sep = "")
  write.csv(demog[, outcome], demog.name, row.names = FALSE)
  mfn <- NA
  if (!is.null(mask)) {
    mask = check_ants(mask)
    mfn <- paste(statdir, "spcamask.nii.gz", sep = "")
    antsImageWrite(mask, mfn)
  }
  args <- list("--svd", paste("network[", matname, ",", mfn, ",", sparseness, ",",
    demog.name, "]", sep = ""), "--l1", 1, "-i", its, "--PClusterThresh", cthresh,
    "-n", nvecs, "-o", outfn, "-z", z, "-s", smooth)
  .Call("sccanX", .int_antsProcessArguments(c(args)), PACKAGE = "ANTsR")
  mydecomp <- read.csv(decomp)
  if (!is.null(mask)) {
    glb <- paste("spca*View1vec*.nii.gz", sep = "")
    fnl <- list.files(path = statdir, pattern = glob2rx(glb), full.names = T,
      recursive = TRUE)[1:nvecs]
    fnll <- list()
    for (i in 1:length(fnl)) {
      img <- antsImageRead(fnl[i], length(dim(mask)))
      fnll <- lappend(fnll, img)
    }
    fnl <- fnll
  }
  if (is.null(mask)) {
    glb <- paste("spcaprojectionsView1vec.csv", sep = "")
    fnl <- list.files(path = statdir, pattern = glob2rx(glb), 
                      full.names = TRUE,
                      recursive = TRUE)
    fnl <- read.csv(fnl)
  }
  glb <- paste("spcaprojectionsView1vec.csv", sep = "")
  fnu <- list.files(path = statdir, pattern = glob2rx(glb), full.names = T, recursive = T)
  fnu <- read.csv(fnu)
  return(list(projections = mydecomp, eigenanatomyimages = fnl, umatrix = fnu))
}