R/sparseDecom2.R

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

#' Convenience wrapper for 2-view eigenanatomy decomposition.
#'
#' Decomposes two matrices into paired sparse eigenevectors to maximize
#' canonical correlation. Note: we do not scale the matrices internally.
#' We leave scaling choices to the user.
#'
#' @param inmatrix input as inmatrix=list(mat1,mat2). n by p input matrix and n
#' by q input matrix , spatial variable lies along columns.
#' @param inmask optional pair of antsImage masks
#' @param sparseness a c(.,.) pair of values e.g c(0.01,0.1) enforces an
#' unsigned 99 percent and 90 percent sparse solution for each respective view
#' @param nvecs number of eigenvector pairs
#' @param its number of iterations, 10 or 20 usually sufficient
#' @param cthresh cluster threshold pair
#' @param statdir temporary directory if you want to look at full output
#' @param perms number of permutations. settings permutations greater than 0
#' will estimate significance per vector empirically. For small datasets, these
#' may be conservative.  p-values depend on how one scales the input matrices.
#' @param uselong enforce solutions of both views to be the same - requires
#'  matrices to be the same size
#' @param z subject space (low-dimensional space) sparseness value
#' @param smooth smooth the data (only available when mask is used)
#' @param robust rank transform input matrices
#' @param mycoption enforce 1 - spatial orthogonality, 2 - low-dimensional
#' orthogonality or 0 - both
#' @param initializationList initialization for first view
#' @param initializationList2 initialization for 2nd view
#' @param ell1 gradient descent parameter, if negative then l0 otherwise use l1
#' @param priorWeight Scalar value weight on prior between 0 (prior is weak)
#' and 1 (prior is strong).  Only engaged if initialization is used
#' @param verbose activates verbose output to screen
#' @param rejector rejects small correlation solutions
#' @param maxBased boolean that chooses max-based thresholding
#' @return outputs a decomposition of a pair of matrices
#' @author Avants BB
#' @examples
#'
#' mat<-replicate(100, rnorm(20))
#' mat2<-replicate(100, rnorm(20))
#' mat<-scale(mat)
#' mat2<-scale(mat2)
#' mydecom<-sparseDecom2(
#' inmatrix = list(mat,mat2),
#' sparseness=c(0.1,0.3),
#' nvecs=3, its=3, perms=0)
#' wt<-0.666
#' mat3<-mat*wt+mat2*(1-wt)
#' mydecom<-sparseDecom2( inmatrix=list(mat,mat3),
#'   sparseness=c(0.2,0.2), nvecs=5, its=10, perms=5 )
#'
#' \dontrun{
#' # a masked example
#' im<-antsImageRead( getANTsRData("r64"))
#' mask<-thresholdImage( im, 250, Inf )
#' dd = sum( mask == 1 )
#' mat1<-matrix( rnorm(dd*10) , nrow=10 )
#' mat2<-matrix( rnorm(dd*10) , nrow=10 )
#' initlist<-list()
#' for ( nvecs in 1:2 ) {
#'   init1<-antsImageClone( mask )
#'   init1[ mask == 1 ]<-rnorm( dd )
#'   initlist<-lappend( initlist, init1 )
#' }
#' ff<-sparseDecom2( inmatrix=list(mat1,mat2), inmask=list(mask,mask),
#'   sparseness=c(0.1,0.1) ,nvecs=length(initlist) , smooth=1,
#'   cthresh=c(0,0), initializationList = initlist ,ell1 = 11 )
#' ### now SNPs ###
#' rf<-usePkg('randomForest')
#' bg<-usePkg('BGLR')
#' if ( bg & rf ) {
#' data(mice)
#' snps<-mice.X
#' numericalpheno<-as.matrix( mice.pheno[,c(4,5,13,15) ] )
#' numericalpheno<-residuals( lm( numericalpheno ~
#'   as.factor(mice.pheno$Litter) ) )
#' nfolds<-6
#' train<-sample( rep( c(1:nfolds), 1800/nfolds ) )
#' train<-( train < 4 )
#' snpd<-sparseDecom2( inmatrix=list( ( as.matrix(snps[train,]) ),
#'   numericalpheno[train,] ), nvecs=20, sparseness=c( 0.001, -0.5 ),
#'   its=3, ell1=0.1 , z=-1 )
#' for ( j in 3:3) {
#' traindf<-data.frame( bmi=numericalpheno[ train,j] ,
#'    snpse=as.matrix( snps[train, ] ) %*% as.matrix( snpd$eig1 ) )
#' testdf <-data.frame( bmi=numericalpheno[!train,j] ,
#'    snpse=as.matrix( snps[!train,] ) %*% as.matrix( snpd$eig1 ) )
#' myrf<-randomForest( bmi ~ . , data=traindf )
#' preddf<-predict(myrf, newdata=testdf )
#' print( cor.test(preddf, testdf$bmi ) )
#' plot( preddf, testdf$bmi )
#' }
#' } # check bg and rf
#' }
#'
#' @export sparseDecom2
sparseDecom2 <- function(
  inmatrix,
  inmask = list(NULL, NULL),
  sparseness = c(0.01, 0.01),
  nvecs = 3,
  its = 20,
  cthresh = c(0, 0),
  statdir = NA,
  perms = 0,
  uselong = 0,
  z = 0,
  smooth = 0,
  robust = 0,
  mycoption = 0,
  initializationList = list(),
  initializationList2 = list(),
  ell1 = 10,
  priorWeight = 0,
  verbose = FALSE,
  rejector=0,
  maxBased=FALSE  ) {
  idim=3
  # safety 1 & 2
  if ( ! is.null( inmask[[1]] ) ) {
    if (!is.antsImage(inmask[[1]]) && !is.na(inmask[[1]])) {
      if ( ncol( inmatrix[[1]] )  !=  sum( inmask[[1]] ) ) {
        stop("Number of columns in matrix X must equal the size of the maskx")
      }
    }
  }
  if ( ! is.null( inmask[[2]] ) ) {
    if (!is.antsImage(inmask[[2]]) && !is.na(inmask[[2]])) {    
      if ( ncol( inmatrix[[2]] )  !=  sum( inmask[[2]] ) ) {
        stop("Number of columns in matrix Y must equal the size of the masky")
      }
    }
  }
  if (is.antsImage(inmask[[1]])) idim=inmask[[1]]@dimension
  if (is.antsImage(inmask[[2]])) idim=inmask[[2]]@dimension
  if (!is.antsImage(inmask[[1]]))
    maskx = new("antsImage", "float",idim) else maskx = antsImageClone( inmask[[1]] )
    if (!is.antsImage(inmask[[2]]))
      masky = new("antsImage", "float",idim) else masky = antsImageClone( inmask[[2]] )
      if ( nrow(inmatrix[[1]]) != nrow(inmatrix[[2]]) )
        stop("Matrices must have same number of rows")
      inmask = list(maskx, masky )
      verbose = as.numeric( verbose )
      if ( robust > 0 )
      {
        inputMatrices = list(
          robustMatrixTransform( inmatrix[[1]] ),
          robustMatrixTransform( inmatrix[[2]] )
        )
      } else inputMatrices = inmatrix
      # helper function allows easier R-based permutation
      sccaner=.sparseDecom2helper2(
        inputMatrices,
        inmask,
        sparseness,
        nvecs,
        its,
        cthresh,
        statdir,
        uselong,
        z,
        smooth,
        robust,
        mycoption,
        initializationList,
        initializationList2,
        ell1,
        priorWeight,
        verbose,
        maxBased
      )
      ccasummary = data.frame(
        corrs = sccaner$corrs,
        pvalues = rep(NA,nvecs)
      )
      if ( perms >  0 )
      {
        ccasummary$pvalues = rep(0 , nvecs )
        nsubs = nrow( inputMatrices[[1]] )
        for ( permer in 1:perms )
        {
          m1 = data.matrix( inputMatrices[[1]][ sample( 1:nsubs ) ,  ] )
          m2 = data.matrix( inputMatrices[[2]][ sample( 1:nsubs ) ,  ] )
          permmatrix = list( m1, m2 )
          sccanerp=.sparseDecom2helper2(
            permmatrix,
            inmask,
            sparseness,
            nvecs,
            its,
            cthresh,
            statdir,
            uselong,
            z,
            smooth,
            robust,
            mycoption,
            initializationList,
            initializationList2,
            ell1,
            priorWeight,
            verbose,
            maxBased
          )
          counter = as.numeric( abs(ccasummary$corrs) < abs(sccanerp$corrs)   )
          ccasummary$pvalues = ccasummary$pvalues + counter
        }
        ccasummary$pvalues = ccasummary$pvalues / perms
      }
      mynames = paste(0:(nvecs-1),sep='')
      if ( length(mynames) <= 10 ) mynames = paste( "00", mynames, sep='')
      if ( length(mynames) >10 ) {
        mynames[1:10] = paste( "00", mynames[1:10], sep='')
        tt = nvecs
        if ( tt > 99 ) tt = 99
        mynames[11:tt] = paste( "0", mynames[11:tt], sep='')
      }
      mynames=paste("Variate",mynames,sep='')
      if ( nvecs > 1 )
      {
        colnames( sccaner$eig1 ) = mynames[1:nvecs]
        colnames( sccaner$eig2 ) = mynames[1:nvecs]
        colnames( sccaner$projections ) = mynames[1:nvecs]
        colnames( sccaner$projections2 ) = mynames[1:nvecs]
      }
      if ( rejector > 0 )
      {
        selector=abs(ccasummary$corrs) > rejector
        if ( sum( selector ) > 1 )
        {
          sccaner$eig1 = sccaner$eig1[ , selector ]
          sccaner$eig2 = sccaner$eig2[ , selector ]
          ccasummary = ccasummary[ selector, ]
        }
      }
      return(
        list(
          projections = sccaner$projections,
          projections2 = sccaner$projections2,
          eig1 = sccaner$eig1,
          eig2 = sccaner$eig2,
          ccasummary = ccasummary,
          sparseness = sparseness
        )
      )
}


.sparseDecom2helper2 <- function(
  inputMatrices,
  inmask,
  sparseness,
  nvecs,
  its,
  cthresh,
  statdir,
  uselong,
  z,
  smooth,
  robust,
  mycoption,
  initializationList,
  initializationList2,
  ell1,
  priorWeight,
  verbose,
  maxBased ) {
  outval = .Call( "sccanCpp",
                  inputMatrices[[1]],
                  inputMatrices[[2]],
                  inmask[[1]],
                  inmask[[2]],
                  sparseness[1],
                  sparseness[2],
                  nvecs,
                  its,
                  cthresh[1],
                  cthresh[2],
                  z,
                  smooth,
                  initializationList,
                  initializationList2,
                  mycoption,
                  ell1,
                  verbose,
                  priorWeight,
                  maxBased,
                  PACKAGE="ANTsR" )
  p1 = inputMatrices[[1]] %*% t(outval$eig1)
  p2 = inputMatrices[[2]] %*% t(outval$eig2)
  outcorrs = diag( cor( p1 , p2  ) )
  if ( priorWeight < 1.e-10 )
  {
    myord = rev( order( abs( outcorrs ) ) )
    outcorrs = outcorrs[ myord ]
    p1 = p1[ , myord ]
    p2 = p2[ , myord ]
    outval$eig1 = outval$eig1[ myord, ]
    outval$eig2 = outval$eig2[ myord, ]
  }
  return(
    list(
      projections = p1,
      projections2 = p2,
      eig1 = t(outval$eig1),
      eig2 = t(outval$eig2),
      corrs = outcorrs
    )
  )
}