2023-09-12 @Atsushi Kawaguchi
In this vignette, the output is omitted. Please refer to the following book for the output.
Kawaguchi A. (2021). Multivariate Analysis for Neuroimaging Data. CRC Press.
The first installation of the package should be performed from the Comprehensive R Archive Network (CRAN) using the following code.
Once the package is installed, this code is no longer required. The package is loaded with the following command line code.
library(oro.nifti) library(oro.dicom) library("imager")
This code must be executed every time to start the R. Other packages are also available. Some depend on the mand package being preinstalled.
One subject image as the template is available in the mand package. The coad to load it is as follows.
The image is compressed because of storage and computation time. The dimension is confirmed as follows.
The image is plotted by the
Other options with the plane argument (such as “axial,” “coronal,”
“sagittal,” and “all”) are available. The default setting is “axial”. If
the argument is specified as
plane="all", three directional
slices at a certain coordinate are plotted.
imgdatamat function reads image files saved in the
nifti format and creates data matrix with subjects in row and voxel in
column (this example does not work).
fnames1 = c("data1.nii", "data2.nii") imgmat = imgdatamat(fnames1, simscale=1/4)
The first argument is file names with the length equaling the number
of subjects (the number of rows in the resulting data matrix). The
simscale is the image resize scale. In this
example, the all sizes (number of voxel) for three direction was reduced
into its quarter size. The output is the list form where the “S” is data
matrix, the “brainpos” is a binary image indicating brain region, and
the “imagedim” is image dimension. The ROI (Region Of Interest) volume
is computed in the “roi” if the ROI argument is TRUE.
imgmat = imgdatamat(fnames1, simscale=1/4, ROI=TRUE)
The resulting map from the statistical analysis such as the t
statistics map from the SPM is represented with overlapping on the
template. For example, the
mand package has the average
difference assuming Alzheimer’s disease and healthy subjects with the
The overlay is implemented by the
data(diffimg) coat(template, diffimg)
Anatomical brain segmentation region is useful for the plot and the interpretation. For example, the Automated Anatomical Labeling (AAL) atlas is used. The data.frame has two columns (“ROIid” and “ROIname”) format.
data(atlasdatasets) atlasname = "aal" atlasdataset = atlasdatasets[[atlasname]] head(atlasdataset)
It is also neccesary to prepare the array data.
data(atlas) tmpatlas = atlas[[atlasname]] dim(tmpatlas)
It has the ROIid as the element.
The anatomical region can be plotted by the
function with regionplot=TRUE.
coat(template, tmpatlas, regionplot=TRUE, atlasdataset=atlasdataset, ROIids = c(1:2, 37:40), regionlegend=TRUE)
The resulting map can be converted into the summary of the anatomical regions.
atlastable(tmpatlas, diffimg, atlasdataset, ROIids = c(1:2, 37:40))
The outputs are the number of voxel in the
column, the percentage of the voxel in the
and the minimum, mean, and maximum valued of the region of the
overlaying map. The order of the table row is in the larger absolute
value of the minimum or maximum values.
The brain image corresponding to the region of interest can be extracted as follows. First, we create a mask image in which the hippocampal region is represented by 1 and others by 0. Then the product of the template and the mask image is taken for each voxel.
hipmask = (tmpatlas == 37) + (tmpatlas == 38) template2 = template * hipmask
The images generated by these processes are plotted from left to right in one slice.
par(mfrow=c(1,3), mar=rep(1,4)) coat(template, pseq=11, paron=FALSE) coat(hipmask, pseq=11, paron=FALSE) coat(template2, pseq=11, paron=FALSE)
The template image (left) and the mask image (middle) are multiplied voxel by voxel to obtain the only hippocampus region image (right).
The sum of the voxel values in the region is calculated as follows.
sum(template[which(hipmask==1, arr.ind = TRUE)])/1000
Such a value is calculated for each region and a dataset with ROI values is created.