R/taskFMRI.R

Defines functions .arCorrection taskFMRI

Documented in taskFMRI

#' Simple taskFMRI function.
#'
#' Input 4D time series matrix. (Perform slice timing correction externally).
#' Estimate hemodynamicRF from block design. Compute brain mask on average bold
#' image.  Get nuisance variables : motion , compcor , globalsignal.
#' High-frequency filter the time series ( externally ). Correct for
#' autocorrelation using bullmore 1996 MRM and AR(2) model with parameters
#' derived from global residual signal. Estimate final glm.
#'
#'
#' @param mat input matrix
#' @param hrf input hrf
#' @param myvars output of getfMRInuisanceVariables
#' @param correctautocorr correction auto correlation boolean
#' @param residualizedesignmatrix boolean
#' @param myformula statistical equation to be assessed at each voxel
#' @return list of betas and other names entries is output
#' @author Avants BB
#' @examples
#'
#' \dontrun{
#' # read the fmri image in and maybe do slice timing correction
#'   fmri<-getANTsRData("pcasl")
#'   fmri<-antsImageRead( fmri )
#' #  fmri<-iMath(fmri,"SliceTimingCorrection","bspline") # optional
#'   myvars<-getfMRInuisanceVariables( fmri, moreaccurate = 0, maskThresh=100 )
#'   mat <- myvars$matrixTimeSeries
#'   mat<-frequencyFilterfMRI(mat, 2.5, freqLo=0.01, freqHi=0.1, opt="butt")
#'   blockfing = c(0, 36, 72 )
#'   hrf <- hemodynamicRF( scans=dim(fmri)[4] , onsets=blockfing ,
#'     durations=rep(  12,  length( blockfing ) ) ,  rt=2.5 )
#'   activationBeta<-taskFMRI(  mat , hrf , myvars )
#' }
#'
#' @export taskFMRI
taskFMRI <- function(mat, hrf, myvars,
  correctautocorr = FALSE, residualizedesignmatrix = FALSE,
  myformula = NA) {
  if (nargs() == 0) {
    print("Usage:  betas<-taskFMRI( x , hrf ) ")
    return(1)
  }
  if ( all(is.na(myformula) ) )
    myformula<-"globalsignal + motion1 +
      motion2 + motion3 + compcorr1 + compcorr2 + compcorr3"
  havemagic<-usePkg("magic")
  if ( ! havemagic )
    {
    print("Must have magic package to use this function")
    return(NA)
    }
  avg <- myvars$avgImage
  mask <- myvars$mask
  nuis <- (myvars$nuisancevariables)
  betas <- rep(NA, ncol(mat))
  desmat <- as.matrix(cbind(myvars$globalsignal, nuis))
  if (residualizedesignmatrix)
    desmat <- residuals(lm(desmat ~ hrf))
  desmat <- cbind(hrf, desmat)
  colnames(desmat) <- c(colnames(hrf), colnames(desmat)  )
  print(colnames(desmat))
  ####################################################### Statistical methods of estimation and inference for # functional MR image
  ####################################################### analysis Bullmore 1996 #
  amat <- mat
  adesmat <- desmat
  myformhrf <- as.formula(paste("amat ~  hrf + ", myformula))
  myformNOhrf <- as.formula(paste("amat ~ 1 + ", myformula))
  print(myformhrf)
  if (correctautocorr) {
    for (i in 1:11) {
      residmat <- residuals(lm(myformhrf, data = data.frame(adesmat)))
      residsig <- apply(residmat, FUN = mean, MARGIN = 1)
      arcoefs <- ar(residsig, FALSE, 2, method = "burg")$ar  # use something similar to  SPM's  autoregression estimate
      if (i == 1)
        initialarcoefs <- arcoefs
      if (abs(arcoefs[1]) > 0.05)
        {
          print(arcoefs)
          mat1 <- magic::ashift(amat, c(1, 0))
          mat1[1, ] <- amat[1, ]
          mat2 <- magic::ashift(amat, c(2, 0))
          mat2[1, ] <- amat[1, ]
          mat2[2, ] <- amat[2, ]
          if (length(arcoefs) == 2)
          amat <- amat - mat1 * arcoefs[1] - mat2 * arcoefs[2]
          if (length(arcoefs) == 1)
          amat <- amat - mat1 * arcoefs[1]
          mat1 <- magic::ashift(adesmat, c(1, 0))
          mat1[1, ] <- adesmat[1, ]
          mat2 <- magic::ashift(adesmat, c(2, 0))
          mat2[1, ] <- adesmat[1, ]
          mat2[2, ] <- adesmat[2, ]
          if (length(arcoefs) == 2)
          adesmat <- adesmat - mat1 * arcoefs[1] - mat2 * arcoefs[2]
          if (length(arcoefs) == 1)
          adesmat <- adesmat - mat1 * arcoefs[1]
        }  # if gt 0.05
    }
    # new regression
    amat <- residuals(lm(myformNOhrf, data = data.frame(adesmat)))
    residsig <- apply(amat, FUN = mean, MARGIN = 1)
    arcoefs <- ar(residsig, FALSE, length(arcoefs), method = "burg")$ar  # use something similar to  SPM's  autoregression estimate
    print(paste("final arcoefs", arcoefs))
    print(paste("initi arcoefs", initialarcoefs))
  }
  progress <- txtProgressBar(min = 0, max = ncol(mat), style = 3)
  for (i in 1:ncol(amat)) {
    vox <- amat[, i]
    # mdl<-lmrob( vox ~ hrf, data = data.frame( adesmat ) )
    mdl <- lm(vox ~ hrf, data = data.frame(desmat))
    betas[i] <- coefficients(summary(mdl))[2, 3]  # probably better way
    # betas[i]<-cor.test( vox , hrf )$est
    setTxtProgressBar(progress, i)
  }
  close(progress)
  betas[is.na(betas)] <- 0
  return(list(beta = betas, fmrimat = amat, adesmat = adesmat))
}



.arCorrection <- function(boldmatin, loops = 3, armethod = "burg") {
  boldmat <- boldmatin
  for (i in 1:loops) {
    arcoefs <- ar(boldmat, FALSE, 2, method = armethod)$ar  # use something similar to  SPM's  autoregression estimate
    if (i == 1)
      initialarcoefs <- arcoefs
    if (abs(arcoefs[1]) > 0.05) {
      mat1 <- magic::ashift(boldmat, c(1, 0))
      mat1[1, ] <- boldmat[1, ]
      mat2 <- magic::ashift(boldmat, c(2, 0))
      mat2[1, ] <- boldmat[1, ]
      mat2[2, ] <- boldmat[2, ]
      if (length(arcoefs) == 2)
        boldmat <- boldmat - mat1 * arcoefs[1] - mat2 * arcoefs[2]
      if (length(arcoefs) == 1)
        boldmat <- boldmat - mat1 * arcoefs[1]
    }
  }
  # new regression
  arcoefs <- ar(boldmat, FALSE, length(arcoefs), method = armethod)$ar  # use something similar to  SPM's  autoregression estimate
  if (mean(abs(sum(arcoefs))) > mean(abs(sum(initialarcoefs)))) {
    print(paste("final arcoefs", arcoefs))
    print(paste("initi arcoefs", initialarcoefs))
  }
  return(list(outmat = boldmat, initcoefs = initialarcoefs, arcoefs = arcoefs))
}
neuroconductor-devel/ANTsR documentation built on April 1, 2021, 1:02 p.m.