Overview for mand
2023-09-12 @Atsushi Kawaguchi
The mand package provides functions for multivariate
analysis for neuroimaging data.
library(mand)## Loading required package: msmaIntroduction
Template
One subject image as the template is available in the mand package. The coad to load it is as follows.
data(template)The image is compressed because of storage and computation time. The dimension is confirmed as follows.
dim(template)## [1] 30 36 30The image is plotted by the coat function.
coat(template)
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.
coat(x=template, plane="all")
Image Data Matrix
The 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
second argument 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, schange=TRUE, simscale=1/4, ROI=TRUE)Overlay
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
array format.
The overlay is implemented by the coat function.
data(diffimg)
coat(template, diffimg)
Atlas
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 = "aal3"
atlasdataset = atlasdatasets[[atlasname]]
head(atlasdataset)## ROIid ROIname
## 1 1 Left Precentral gyrus
## 2 2 Right Precentral gyrus
## 3 3 Left Superior frontal gyrus-dorsolateral
## 4 4 Right Superior frontal gyrus-dorsolateral
## 5 5 Left Middle frontal gyrus
## 6 6 Right Middle frontal gyrusIt is also neccesary to prepare the array data.
data(atlas)
tmpatlas = atlas[[atlasname]]
dim(tmpatlas)## [1] 30 36 30It has the ROIid as the element.
tmpatlas[,15:20,12]## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0
## [3,] 90 90 90 90 90 90
## [4,] 94 94 90 90 90 90
## [5,] 94 94 90 90 90 90
## [6,] 94 90 90 90 90 90
## [7,] 94 0 0 0 90 90
## [8,] 0 0 0 0 0 0
## [9,] 60 44 42 42 42 0
## [10,] 60 44 42 42 42 42
## [11,] 44 44 42 42 42 42
## [12,] 44 44 44 0 0 42
## [13,] 52 0 0 0 0 0
## [14,] 100 0 0 162 162 164
## [15,] 114 0 0 166 166 0
## [16,] 0 0 0 0 0 0
## [17,] 99 0 0 165 0 0
## [18,] 101 101 0 163 163 163
## [19,] 43 0 0 0 0 41
## [20,] 0 43 41 41 41 41
## [21,] 43 43 41 41 41 41
## [22,] 43 43 41 41 41 41
## [23,] 0 0 0 0 0 0
## [24,] 93 0 0 0 0 0
## [25,] 93 0 0 89 89 89
## [26,] 93 89 89 89 89 89
## [27,] 89 89 89 89 89 89
## [28,] 89 89 89 89 89 89
## [29,] 0 89 89 89 89 0
## [30,] 0 0 0 0 0 0The anatomical region can be plotted by the coat
function with regionplot=TRUE.
coat(template, tmpatlas, regionplot=TRUE,
atlasdataset=atlasdataset, ROIids = c(1:2, 41:44), regionlegend=TRUE)
The resulting map can be converted into the summary of the anatomical regions.
atlastable(tmpatlas, diffimg, atlasdataset, ROIids = c(1:2, 41:44))## ROIid ROIname sizepct sumvalue Min. Mean Max.
## 41 41 Left Hippocampus 0.591 -7.190 -1.000 0 0.448
## 42 42 Right Hippocampus 0.733 -12.110 -0.968 0 0.505
## 43 43 Left Parahippocampal gyrus 0.732 -5.495 -0.875 0 0.530
## 1 1 Left Precentral gyrus 0.617 8.300 -0.568 0 0.802
## 44 44 Right Parahippocampal gyrus 0.692 -3.147 -0.750 0 0.690
## 2 2 Right Precentral gyrus 0.632 -0.305 -0.735 0 0.750The outputs are the number of voxel in the sizenum
column, the percentage of the voxel in the sizepct column,
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 == 41) + (tmpatlas == 42)
template2 = template * hipmaskThe 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=12, paron=FALSE)
coat(hipmask, pseq=12, paron=FALSE)
coat(template2, pseq=12, 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## [1] 0.07090703Such a value is calculated for each region and a dataset with ROI values is created.
Principal Component Analysis
Generate Simulation Data
The mand package has a function to generate simulation
brain data from the base image, the difference image and the standard
deviation image. These basic images are loaded as follows.
data(baseimg)
data(diffimg)
data(mask)
data(sdevimg)The number of voxels in the original 3D image is as follows.
dim(baseimg)## [1] 30 36 30To understand the result easily, the difference region was restricted to Parahippocampus and Hippocampus.
diffimg2 = diffimg * (tmpatlas %in% 41:44)An image data matrix with subjects in the rows and voxels in the
columns was generated by using the simbrain function.
img1 = simbrain(baseimg = baseimg, diffimg = diffimg2,
sdevimg=sdevimg, mask=mask, n0=20, c1=0.01, sd1=0.05)The base image, the difference image and the standard deviation image were specified in the first three arguments. The out-of-brain region was specified by the mask argument, which was the binary image. The remaining arguments are the number of subjects per group, the coefficient multiplied by the difference image and the standard deviation for the noise.
The data matrix dimension was 40(subject) x 12620(voxel).
dim(img1$S)## [1] 40 12620The rec function creates the 3D image from the
vectorized data (the first subject).
coat(rec(img1$S[1,], img1$imagedim, mask=img1$brainpos))
The standard deviation image is created from the resulting data matrix.
sdimg = apply(img1$S, 2, sd)
coat(template, rec(sdimg, img1$imagedim, mask=img1$brainpos))
Principal Component Analysis
If the input is a matrix, a principal component analysis is
implemented by the msma function of the msma
package. Principal component analysis with the number of components of 2
for the image data matrix is performed as follows.
(fit111 = msma(img1$S, comp=2))## Call:
## msma.default(X = img1$S, comp = 2)
##
## Numbers of non-zeros for X block:
## comp1 comp2
## block1 9865 9865
##
## Numbers of non-zeros for X super:
## comp1 comp2
## comp1-1 1 1The scatter plots for the score vectors specified by the argument
v. The argument axes is specified by the two
length vectors on which components are displayed.
plot(fit111, v="score", axes = 1:2, plottype="scatter")
The weight (loading) vectors can be obtained and reconstructed as follows.
midx = 1 ## the index for the modality
vidx = 1 ## the index for the component
Q = fit111$wbX[[midx]][,vidx]
outstat1 = rec(Q, img1$imagedim, mask=img1$brainpos)The reconstructed loadings as image is overlayed on the template.
coat(template, outstat1)
The output is unclear; however, this will be improved later.
Two-steps dimension reduction
Basis Expansion
This is an example of the two-step dimension reduction.
Generate radial basis function.
B1 = rbfunc(imagedim=img1$imagedim, seppix=2, hispec=FALSE, mask=img1$brainpos)Multiplying the basis function to image data matrix.
SB1 = basisprod(img1$S, B1)The original dimension was 12620.
dim(img1$S)## [1] 40 12620The dimension was reduced to 1604.
dim(SB1)## [1] 40 1604Principal Component Analysis
The Principal Component Analysis (PCA) is applied to the dimension-reduced image.
(fit211 = msma(SB1, comp=2))## Call:
## msma.default(X = SB1, comp = 2)
##
## Numbers of non-zeros for X block:
## comp1 comp2
## block1 1604 1604
##
## Numbers of non-zeros for X super:
## comp1 comp2
## comp1-1 1 1The loading is reconstructed to the original space by using the
rec function.
Q = fit211$wbX[[1]][,1]
outstat1 = rec(Q, img1$imagedim, B=B1, mask=img1$brainpos)The plotted loading is smoother than the one without the dimension reduction by the basis function.
coat(template, outstat1)
Sparse PCA
If lambdaX (>0) is specified, a sparse principal
component analysis is implemented.
(fit112 = msma(SB1, comp=2, lambdaX=0.06))## Call:
## msma.default(X = SB1, comp = 2, lambdaX = 0.06)
##
## Numbers of non-zeros for X block:
## comp1 comp2
## block1 53 40
##
## Numbers of non-zeros for X super:
## comp1 comp2
## comp1-1 1 1The plotted loading is narrower than that from the PCA.
Q = fit112$wbX[[midx]][,vidx]
outstat1 = rec(Q, img1$imagedim, B=B1, mask=img1$brainpos)
outstat2 = outstat1
coat(template, outstat2)
The region is reported as follows to be compared with the next method.
atlastable(tmpatlas, outstat2, atlasdataset)## ROIid ROIname sizepct sumvalue Min. Mean Max.
## 154 154 Right Anterior cingulate cortex-pregenual 1.000 208.904 0 0.006 7.480
## 153 153 Left Anterior cingulate cortex-pregenual 1.000 133.103 0 0.004 6.944
## 156 156 Right Anterior cingulate cortex-supracallosal 1.000 136.198 0 0.004 6.880
## 155 155 Left Anterior cingulate cortex-supracallosal 1.000 106.948 0 0.003 6.449
## 19 19 Left Superior frontal gyrus-medial 1.000 237.574 0 0.007 6.166
## 20 20 Right Superior frontal gyrus-medial 1.000 226.415 0 0.007 6.114
## 38 38 Right Middle cingulate & paracingulate gyri 0.459 40.319 0 0.001 5.686
## 4 4 Right Superior frontal gyrus-dorsolateral 0.845 257.710 0 0.008 5.165
## 6 6 Right Middle frontal gyrus 0.976 127.759 0 0.004 3.423
## 152 152 Right Anterior cingulate cortex-subgenual 1.000 13.537 0 0.000 3.414Supervised Sparse PCA
The simbrain generates the synthetic brain image data
and the binary outcome. The outcome Z is obtained.
Z = img1$ZIf the outcome Z is specified in the msma function, a
supervised sparse principal component analysis is implemented.
(fit113 = msma(SB1, Z=Z, comp=2, lambdaX=0.06, muX=0.5))## Call:
## msma.default(X = SB1, Z = Z, comp = 2, lambdaX = 0.06, muX = 0.5)
##
## Numbers of non-zeros for X block:
## comp1 comp2
## block1 48 38
##
## Numbers of non-zeros for X super:
## comp1 comp2
## comp1-1 1 1The plotted loading is located differently from the sparse PCA.
Q = fit113$wbX[[1]][,1]
outstat1 = rec(Q, img1$imagedim, B=B1, mask=img1$brainpos)
outstat2 = -outstat1
coat(template, outstat2)
The region near the hippocampus, which differs from the sparse PCA (without supervision).
atlastable(tmpatlas, outstat2, atlasdataset)## ROIid ROIname sizepct sumvalue Min. Mean Max.
## 83 83 Left Heschls gyrus 1.000 63.252 0 0.002 6.878
## 41 41 Left Hippocampus 1.000 129.029 0 0.004 6.666
## 85 85 Left Superior temporal gyrus 1.000 297.309 0 0.009 6.196
## 145 145 Left Pulvinar medial 1.000 35.446 0 0.001 6.040
## 33 33 Left Insula 0.917 100.941 0 0.003 5.910
## 129 129 Left Ventral posterolateral 1.000 39.922 0 0.001 5.809
## 149 149 Left Pulvinar inferior 1.000 5.286 0 0.000 5.286
## 13 13 Left Rolandic operculum 1.000 125.476 0 0.004 5.222
## 43 43 Left Parahippocampal gyrus 1.000 71.732 0 0.002 5.095
## 51 51 Left Lingual gyrus 0.853 66.548 0 0.002 4.585The loading for the second component
Q = fit113$wbX[[1]][,2]
outstat1 = rec(Q, img1$imagedim, B=B1, mask=img1$brainpos)
outstat2 = -outstat1
coat(template, outstat2)
atlastable(tmpatlas, outstat2, atlasdataset)## ROIid ROIname sizepct sumvalue Min. Mean Max.
## 55 55 Left Middle occipital gyrus 0.988 186.365 0 0.006 7.968
## 53 53 Left Superior occipital gyrus 1.000 116.745 0 0.004 7.709
## 63 63 Left Superior parietal gyrus 1.000 322.454 0 0.010 7.646
## 71 71 Left Precuneus 1.000 307.242 0 0.009 7.467
## 65 65 Left Inferior parietal gyrus-excluding supramarginal and angular gyri 1.000 179.870 0 0.006 6.838
## 69 69 Left Angular gyrus 1.000 120.055 0 0.004 6.088
## 49 49 Left Cuneus 1.000 118.538 0 0.004 5.685
## 37 37 Left Middle cingulate & paracingulate gyri 0.894 35.244 0 0.001 5.265
## 39 39 Left Posterior cingulate gyrus 1.000 23.401 0 0.001 3.521
## 47 47 Left Calcarine fissure and surrounding cortex 1.000 45.335 0 0.001 2.301This is similar to the result from the sparse PCA (without supervision).
The following method can be used to plot the weights of several
components simultaneously. It is first reconstructed in three dimensions
with the multirec function and then plotted with the
multicompplot function. It is set to display four columns
per component.
ws = multirec(fit113, imagedim=img1$imagedim, B=B1, mask=img1$brainpos)
multicompplot(ws, template, col4comp=4)