R/fsl_generate_fsf_ev_syntax.R

In UNCDEPENdLab/fmri.pipeline: Analysis of fMRI Data on HPC clusters

#' This function generates syntax for FSL Feat .fsf files for the EVs tab of a higher-level fMRI analysis.
#' It accepts a numeric design matrix whose colum names correspond to individual EVs in the model
#'
#' @details This allows you to generate a design matrix dynamically based on a numeric matrix that has been
#' setup in R using a function such as lm() to handle more complex designs appropriately. The syntax
#' generated by this function can be combined with a .fsf template file to implement the entire LVL2 or LVL3
#' analysis setup.
#'
#' @param inputs A vector of inputs (usually lower-level .feat directories) corresponding to rows of \code{dmat}.
#' @param dmat A numeric matrix whose rows specify individual inputs and columns specify EVs. Columns should be named by EV!
#' @param group_membership An optional vector specifying group membership of each row in \code{dmat} in order to
#'      estimate separate variances for multiple groups. If not specified, variances will be modeled as one group.
#'
#' @return A character vector containing .fsf syntax for the contrasts portion of a Feat analysis based on \code{cmat}.
#'
#' @examples
#'   dmat <- matrix(rnorm(1000), ncol=10)
#'   colnames(dmat) <- letters[1:10]
#'
#'   inputs <- replicate(100, tempfile())
#'
#'   result <- fsl_generate_fsf_ev_syntax(inputs, dmat)
#'
#' @author Michael Hallquist
#' @export
fsl_generate_fsf_ev_syntax <- function(inputs, dmat, group_membership=NULL) {

  stopifnot(length(inputs) == nrow(dmat))

  #overall number of inputs
  fsf_syntax <- c(
    "# Total volumes",
    paste0("set fmri(npts) ", length(inputs)),
    "",
    "# Number of first-level analyses",
    paste0("set fmri(multiple) ", length(inputs)),
    ""
  )

  #locations of inputs
  for (i in seq_along(inputs)) {
    fsf_syntax <- c(fsf_syntax,
      paste0("# 4D AVW data or FEAT directory (", i, ")"),
      paste0("set feat_files(", i, ") \"", inputs[i], "\""),
      ""
    )
  }

  #number of EVs
  fsf_syntax <- c(fsf_syntax,
    "# Number of EVs",
    paste0("set fmri(evs_orig) ", ncol(dmat)),
    paste0("set fmri(evs_real) ", ncol(dmat)),
    "set fmri(evs_vox) 0"
  )

  #loop over EV columns in the design
  for (j in seq_len(ncol(dmat))) {
    fsf_syntax <- c(fsf_syntax,
      paste0("# EV ", j, " title"),
      paste0("set fmri(evtitle", j, ") \"", colnames(dmat)[j], "\""),
      "",
      paste0("# Basic waveform shape (EV ", j, ")"),
      "# 0 : Square",
      "# 1 : Sinusoid",
      "# 2 : Custom (1 entry per volume)",
      "# 3 : Custom (3 column format)",
      "# 4 : Interaction",
      "# 10 : Empty (all zeros)",
      paste0("set fmri(shape", j, ") 2"),
      "",
      paste0("# Convolution (EV ", j, ")"),
      "# 0 : None",
      "# 1 : Gaussian",
      "# 2 : Gamma",
      "# 3 : Double-Gamma HRF",
      "# 4 : Gamma basis functions",
      "# 5 : Sine basis functions",
      "# 6 : FIR basis functions",
      "# 8 : Alternate Double-Gamma",
      paste0("set fmri(convolve", j, ") 0"),
      "",
      paste0("# Convolve phase (EV ", j, ")"),
      paste0("set fmri(convolve_phase", j, ") 0"),
      "",
      paste0("# Apply temporal filtering (EV ", j, ")"),
      paste0("set fmri(tempfilt_yn", j, ") 0"),
      "",
      paste0("# Add temporal derivative (EV ", j, ")"),
      paste0("set fmri(deriv_yn", j, ") 0"),
      "",
      paste0("# Custom EV file (EV ", j, ")"),
      paste0("set fmri(custom", j, ") \"dummy\""),
      ""
    )

    #setup EV orthogonalizations (none supported at present)
    #also, not sure why you need a section for orthogonalization wrt EV 0, but keeping it for consistency with Feat
    for (w in 0:ncol(dmat)) {
      fsf_syntax <- c(fsf_syntax,
        paste0("# Orthogonalise EV ", j, " wrt EV ", w),
        paste0("set fmri(ortho", j, ".", w, ") 0"),
        "")
    }

    #setup regressor values for this EV
    for (i in seq_len(nrow(dmat))) {
      fsf_syntax <- c(fsf_syntax,
        paste0("# Higher-level EV value for EV ", j, " and input ", i),
        paste0("set fmri(evg", i, ".", j, ") ", round(dmat[i, j], 5)),
        ""
      )
    }
  }

  #setup group membership
  if (is.null(group_membership)) {
    group_membership <- rep(1, length(inputs)) #one group
  } else if (length(group_membership) != length(inputs)) {
    stop("group_membership vector has different length than inputs")
  }

  for (i in seq_along(group_membership)) {
    fsf_syntax <- c(fsf_syntax,
      paste0("# Group membership for input ", i),
      paste0("set fmri(groupmem.", i, ") ", group_membership[i]),
      ""
    )

  }

  return(fsf_syntax)

}