R/ffdEvidenceFETS.R

In BayesDLMfMRI: Statistical Analysis for Task-Based Fmri Data

#' @import pbapply
#' @importFrom  Rcpp evalCpp

#' @name ffdEvidenceFETS
#' @title ffdEvidenceFETS
#' @description
#' \loadmathjax
#' This function can be used to build activation maps for task-based fMRI data. 
#' @references
#' \insertRef{CARDONAJIMENEZ2021107297}{BayesDLMfMRI}
#' 
#' \insertRef{cardona2021bayesdlmfmri}{BayesDLMfMRI}
#' @details
#' Every voxel from the 4D array image is clustered with its nearest neighbors. There are as many clusters as voxels in the image. Then, activation maps are obtained by fitting a multivariate dynamic linear model on every cluster of voxels. 
#' The resulting activation evidence measure for every voxel is obtained using the Forward Estimated Trajectories Sampler (FETS) algorithm. To deeply understand the method implemented in this package, a reading of \insertCite{CARDONAJIMENEZ2021107297}{BayesDLMfMRI} and \insertCite{cardona2021bayesdlmfmri}{BayesDLMfMRI} is mandatory.
#' @param ffdc  a 4D array (\code{ffdc[i,j,k,t]}) that contains the sequence of MRI images that are meant to be analyzed. \code{(i,j,k)} define the position of the observed voxel at time t.
#' @param covariates a data frame or matrix whose columns contain the covariates related to the expected BOLD response obtained from the experimental setup.
#' @param m0 the constant prior mean value for the covariates parameters and common to all voxels within every neighborhood at \code{t=0} (\code{m=0} is the default value when no prior information is available). For the case of available prior information, \code{m0} can be defined as a \mjseqn{p\times q} matrix, where \mjseqn{p} is the number of columns in the covariates object and \mjseqn{q} is the cluster size.
#' @param Cova a positive constant that defines the prior variances for the covariates parameters at \code{t=0} (\code{Cova=100} is the default value when no prior information is available). For the case of available prior information, \code{Cova} can be defined as a \mjseqn{p\times p} matrix, where \mjseqn{p} is the number of columns in the covariates object.
#' @param delta a discount factor related to the evolution variances. Recommended values between \code{0.85<delta<1}. \code{delta=1} will yield results similar to the classical general linear model.
#' @param S0 prior covariance structure among voxels within every cluster at \code{t=0}. \code{S0=1} is the default value when no prior information is available and defines an \mjseqn{q\times q} identity matrix. For the case of available prior information, \code{S0} can be defined as an \mjseqn{q\times q} matrix, where \mjseqn{q} is the common number of voxels in every cluster.
#' @param n0 a positive hyperparameter of the prior distribution for the covariance matrix \code{S0} at \code{t=0} (\code{n0=1} is the default value when no prior information is available). For the case of available prior information, \code{n0} can be set as \code{n0=np}, where \code{np} is the number of MRI images in the pilot sample.
#' @param N1 is the number of images (\code{2<N1<T}) from the \code{ffdc} array employed in the model fitting. \code{N1=NULL} (or equivalently \code{N1=T}) is its default value, taking all the images in the \code{ffdc} array for the fitting process.
#' @param Nsimu1 is the number of simulated on-line trajectories related to the state parameters. These simulated curves are later employed to compute the posterior probability of voxel activation.
#' @param Cutpos1 a cutpoint time from where the on-line trajectories begin. This parameter value is related to an approximation from a t-student distribution to a normal distribution. Values equal to or greater than 30 are recommended (\code{30<Cutpos1<T}).  
#' @param r1 a positive integer number that defines the distance from every voxel with its most distant neighbor. This value determines the size of the cluster. The users can set a range of different \code{r1} values: \mjseqn{r1 = 0, 1, 2, 3, 4}, which leads to \mjseqn{q = 1, 7, 19, 27, 33}, where \mjseqn{q} is the size of the cluster.
#' @param perVol helps to define a threshold for the voxels considered in the analysis. For example, \code{Min.vol = 0.10} means that all the voxels with values
#' below to \code{max(ffdc)*perVol} can be considered irrelevant and discarded from the analysis.
#' @param Test test type either \code{"LTT"} (Average cluster effect) or \code{"JointTest"} (Joint effect).
#' @param Ncores a positive integer indicating the number of threads or cores to be used in the computation of the activation maps.
#' @param seed random seed.
#' @return It returns a list of the type \code{res[[p]][x,y,z]}, where \code{p} represents the column position in
#' the covariates matrix and \code{[x,y,z]} represent the voxel position in the brain image.
#' @examples
#'\dontrun{
#' fMRI.data  <- get_example_fMRI_data()
#' data("covariates", package="BayesDLMfMRI")
#' res <- ffdEvidenceFETS(ffdc = fMRI.data, covariates = Covariates,
#'                        m0 = 0, Cova = 100, delta = 0.95, S0 = 1, 
#'                        n0 = 1, Nsimu1 = 100, Cutpos1 = 30, 
#'                        r1 = 2, Test = "JointTest",  Ncores = 1)
#' str(res)
#' }
#' @export
ffdEvidenceFETS = function(ffdc, covariates, m0=0, Cova=100,
                           delta=0.95, S0=1, n0=1, N1=FALSE, 
                           Nsimu1 = 100, Cutpos1=30, r1 = 1, 
                           perVol = 0.10, Test = "LTT", Ncores = NULL, seed=NULL){
  
  if(is.logical(N1)) {
    if(N1==FALSE){N1=dim(ffdc)[4]}
  }
  
  # validation
  Ncores  <- .get_n_cores(Ncores)

  .validate_input(
    N1=N1, Test=Test, 
    ffdc=ffdc, covariates=covariates,
    n0=n0, Nsimu1=Nsimu1, perVol=perVol,
    Cutpos1=Cutpos1, r1=r1, delta=delta
  )


  #TAKING THE POSITIONS FROM THE 4D IMAGE WITH NON-NULL VALUES 
  posiffd1 <- which(ffdc[,,,1] != 0, arr.ind = TRUE)
  
  
  if(Test == "LTT"){
    #COMPUTING THE EVIDENCE FOR BRAIN ACTIVATION: VOXEL-WISE ANALYSIS
    set.seed(seed)
    ffd.out = pbapply::pbapply(posiffd1, 1, .ffdSingleVoxelFETS, covariates, ffdc, m0, Cova,
                               delta, S0, n0, N1, Nsimu1, Cutpos1, Min.vol = perVol*max(ffdc), r1, Test, cl = Ncores)
    
    
    vol.evidence <- list()
    
    
    for(k in 1:(dim(covariates)[2])){
      vol.evidence[[k]] <- array(0, dim(ffdc)[1:3])
    }
    
    
    for(j in 1:dim(covariates)[2]){
      for(ii in 1:dim(posiffd1)[1]){
        vol.evidence[[j]][posiffd1[ii,1], posiffd1[ii,2], posiffd1[ii,3]] <- ffd.out[j, ii]
      }
    }

    attr(vol.evidence, "class") <- "fMRI_single_evidence"
    return(vol.evidence) 
  }
  if(Test == "JointTest"){
    #COMPUTING THE EVIDENCE FOR BRAIN ACTIVATION: VOXEL-WISE ANALYSIS
    set.seed(seed)
    ffd.out = pbapply::pbapply(posiffd1, 1, .ffdSingleVoxelFETS, covariates, ffdc, m0, Cova,
                               delta, S0, n0, N1, Nsimu1, Cutpos1, Min.vol = perVol*max(ffdc), r1, Test, cl = Ncores)
    #number of tests from the output of ffdIndividualVoxelLTT  (Joint and marginal)
    Ntest <- 2
    pAux <- dim(covariates)[2]
    vol.evidence <- list()
    
    
    for(k in 1:(dim(covariates)[2]*Ntest)){
      vol.evidence[[k]] <- array(0, dim(ffdc)[1:3])
    }
    
    
    for(j in 1:dim(covariates)[2]){
      for(ii in 1:dim(posiffd1)[1]){
        vol.evidence[[j]][posiffd1[ii,1], posiffd1[ii,2], posiffd1[ii,3]] <- ffd.out[[ii]]$EvidenceJoint[j]
        vol.evidence[[Ntest+j]][posiffd1[ii,1], posiffd1[ii,2], posiffd1[ii,3]] <- ffd.out[[ii]]$EvidenceMargin[j]
      }
    }
    attr(vol.evidence, "class") <- "fMRI_single_evidence"
    return(vol.evidence)
  }
  
  
  
}