Please skip down to the analysis part.

1 Setup

1.1 Paths/Libraries

library(plyr)
suppressMessages(library(niftir))
## 
## Loading niftir version < 1.0.
source("/mnt/nfs/psych4/haxby01/scripts/information/50_info_measures/info_functions.R")

base       <- "/mnt/nfs/psych4/haxby01"
roidir     <- file.path(base, "rois/highres2func")
scriptdir  <- file.path(base, "scripts/information")
taskdir    <- file.path(base, "analysis/task_activity")
prepdir    <- file.path(base, "analysis/preprocessed")

subjects  <- as.character(read.table(file.path(scriptdir, "sublist.txt"))[,1])
roi.names <- c("rh_ventrolateral", "rh_ventromedial")
#categories<- c("face", "house", "chair", "scrambled")
categories<- c("face", "house") # keep it 'simple' to start
runs      <- c("even", "odd")
categories.title<- c("Face", "House")#, "Chair", "Scrambled")
runs.title      <- c("Even", "Odd")


s <- sprintf

1.2 Load the Data

NOT USING THIS….LOAD DATA.

This will load the data as a fairly complicated structure. Here’s an example, if I want to access data for one ROI and one subject, here’s how:

res <- betas.ave\(rh_ventrolateral\)sub001

If you have a script, then the above might look like:

roi <- “rh_ventrolateral” subject <- “sub001” res <- betas.ave[[roi]][[subject]]

The res output above is an array of voxels x categories x runs.

Note that you may not actually need this data….so it’s commented out for now.

# all.betas <- llply(roi.names, function(roi.name) {
#   cat("\n", roi.name, "\n")
#   ret <- llply(subjects, function(subject) {
#     cat("-", subject, "\n")
#     
#     # Load the ROI
#     roifn <- file.path(roidir, subject, s("%s.nii.gz", roi.name))
#     hdr   <- read.nifti.header(roifn)
#     roi   <- read.mask(roifn)
#     
#     # Get data in ROI across runs and categories
#     # output (ret) will be voxel x category x run
#     ret <- laply(runs, function(run) {
#       laply(categories, function(category) {
#         fn <- file.path(taskdir, subject, s("%s_runs.reml", run), "stats", 
#                         s("coef_%s.nii.gz", category))
#         img <- read.nifti.image(fn)[roi]
#         img
#       })
#     })
#     dimnames(ret) <- list(run=runs.title, category=categories.title, voxel=NULL)
#     ret <- aperm(ret, c(3,2,1))
#     
#     ret
#   })
#   names(ret) <- subjects
#   ret
# })
# names(all.betas) <- roi.names

1.3 Get input file paths

What you will need is the ROI and data file paths to feed into the info measures function below.

The output below will be a list of lists of lists etc.

all.filepaths\(rh_ventrolateral\)sub001\(even\)face

all.filepaths <- llply(roi.names, function(roi.name) {
  ret <- llply(subjects, function(subject) {
    ret <- llply(runs, function(run) {
      ret <- llply(categories, function(category) {
        roifn <- file.path(roidir, subject, s("%s.nii.gz", roi.name))
        datfn <- file.path(taskdir, subject, s("%s_runs.reml", run), "stats", 
                           s("coef_%s.nii.gz", category))
        list(roi=roifn, dat=datfn)
      })
      names(ret) <- categories
      ret
    })
    names(ret) <- runs
    ret
  })
  names(ret) <- subjects
  ret
})
names(all.filepaths) <- roi.names

2 Analysis

Let’s now calculate each of our measures to summarize the pattern of activity in each of our ROIs. This includes:

  1. Mean
  2. Variance
  3. Skewness
  4. Kurtosis
  5. Non-Gaussianity (higher values mean more non-gaussian)
  6. Differential Entropy
  7. Binned Entropy
  8. Compression (using NIFTI-GZ and output is in bytes)
  9. Number of local maxima (peaks)
  10. Smoothness of data in ROI (FWHM in mm)

2.1 Function

First we will make a function to run everything for one spatial pattern.

spatial.pattern.analysis <- function(dat.fname, roi.fname) {
  library(moments)
  library(nortest)
  
  # Load the ROI
  roi   <- read.mask(roi.fname)
  # Load the data
  vec   <- read.nifti.image(dat.fname)[roi]
  
  # Calculate Measures
  m <- c()
  
  # Mean
  m$mean <- mean(vec)
  
  # Standard-Deviation
  m$sd <- sd(vec)
  
  # Skewness
  m$skew <- skewness(vec)
  
  # Kurtosis
  m$kurtosis <- kurtosis(vec)
  
  # Gaussianity of Data (higher values give more non-gaussian)
  m$nongaussian <- as.numeric(ad.test(vec)$statistic)
  
  # Entropy
  ## Use R's KNN approach (differential)
  m$entropy.knn <- FNN::entropy(vec, k=3)[3]
  ## R's binned approach
  m$entropy.bin <- entropy::entropy(entropy::discretize(vec, numBins=30)) # arbitrary binning
  
  # Compression (convert to bytes)
  tmp <- compressed_nifti_sizes(as.matrix(vec), digits=0, odt='short')[2]
  m$compressed <- as.numeric(tmp)*1024
  rm(tmp)
  
  # Number of local maxima (peaks)
  tmp <- roi.local.maxima(dat.fname, roi.fname) # gives the peak table
  m$peaks <- nrow(tmp)
  
  # 3D Smoothness (FWHM)
  m$fwhm <- roi.smoothness(dat.fname, roi.fname)
  
  as.data.frame(m)
}

2.2 Call

Sample call for one subject and even runs for the R ventrolateral ROI. Note the output seems to suggest that the houses are really where the action is at here for this subject/runs.

paths <- all.filepaths$rh_ventrolateral$sub001$even$face
res1 <- spatial.pattern.analysis(paths$dat, paths$roi)
paths <- all.filepaths$rh_ventrolateral$sub001$even$house
res2 <- spatial.pattern.analysis(paths$dat, paths$roi)
res <- rbind(face=res1, house=res2)
print(res, digits=3)
##       mean   sd skew kurtosis nongaussian entropy.knn entropy.bin
## face  4.11 12.4 0.66     6.68        2.50        3.92        2.51
## house 7.09 12.9 1.43     8.17        8.48        3.78        2.38
##       compressed peaks fwhm
## face         655    15 4.43
## house        645    18 4.37