So in attempting to replicate analysis pipeline from Tambini & Davachi, PNAS 2013, Persistence of hippocampal multivoxel patterns into postencoding rest is related to memory I'm hoping to use PCA as a noise reduction technique in MATLAB but I have never done this before and have some technical questions. I'm coding this myself and am very unsure I'm doing it correctly, which worries me.
Basically my input are individual NII files (128x128x33x240, here first three dimensions encode voxel coordinate in space and the last one encodes time) of masked CSF/WM (Cerebrospinal fluid & white matter) activations which I'm hoping to run PCA on - these are areas of no interest in that activations here do not represent 'real' neural data and so PCA is done on them so the most robust components can be removed via regression. Then I'd like to retain the top however-many components (e.g. top 8), create an output file consisting of just those components and use that output file in a regression via the AFNI program 3dTproject. Most of my questions are technical but they are numerous so I apologize.
Given that the input data sets are already masked - i.e. everything outside of these areas is set to zero so that in my original 128x128x33x240 matrix, I have mostly zeros and then values at the 3 dimensions in space (at all 240 time points) where the WM/CSF resides area equal to 'activation' values - in that everything outside of the WM/CSF (white matter & cerebrospinal fluid) is equal to zero, should I mask again once in MATLAB (with logical arrays or the like)? In other words, if I input the original 128x128x33x240 image with all of the zeros will all the zeroes affect the PCA in any way?
Okay this is a big one for me. I have done a lot of reading on this and think I have the answer but want to confirm because of how important it is. So if I ultimately use the MATLAB
pcafunction to run the PCA I need a 2D input. This input should be VOXELSxTIME and not TIMExVOXELS, correct? So e.g. will be 25654x240 and not 240x25654? I understand that the number of PCs will truncate at 240 regardless of the input matrix dimensions but obviously this needs to be correct in order to get the right output.
Okay so my current command for actually running the PCA is:
[COEFF, SCORE, LATENT] = pca(image);
COEFFoutputs a 240x240 matrix.
LATENTis 240x1. I know
LATENTrepresents the eigenvalues which I will be using to choose which components to keep.
COEFFis the eigenvectors?
SCOREis said (in the help) to be the representation of the input image into the PC space.
After figuring out how many PCs to keep I want to create my output file consisting of ONLY the top 8 (or whatever) PCs I'm retaining and I'm doing this by multiplying the original image (e.g. 25654x240) by the eigenvector matrix of the retained components (e.g. 240x8) so the result will be a 25654x8 matrix. This is my problem - I don't know what to do from here. I want to be able to use the PCs as a regressor and to have my output be a 3DxTime output. The exact line in the manuscript here is:
After selecting the number of PCs, the associated temporal projection or score of each selected PC was then removed from the BOLD data in a voxelwise fashion by linear regression.
The exact section I'm looking at is in the Supplementary Materials page 2/11 and is the bottom left paragraph (under heading, "Functional MRI Data Acquisition and Preprocessing" all the way to the end of the section).
If my understanding is correct, each column of my matrix represents the voxel-wide representation of the 1st, 2nd, etc. PC. So what does "temporal projection" mean in this case? And how do I use this in a regression?
I ask this because the regression programs in fMRI analysis softwares like AFNI's 3dTproject take input files of size TIMEx1 (so 240x1 in my case) which represents temporal PCA as opposed to a spatial map, which I have here. I'm assuming I have to reshape the matrix to the original size (e.g. 128x128x33x8) with the fourth dimension now being 8 instead of 240 so 128x128x33x1 is PC image 1, 128x128x33x2 is PC image 2, etc. Do I save all the 4D matrices separately? If I then use 3dTproject to regress these out, what option do I use as it is literally a single image and not a timeseries (240x1) file? Or does the 'temporal project' statement indicate that I DO need a time component to these inputs and if so how do I get that to be part of it here?
Is singular value decomposition something I should be looking into?
This is a big question and I would be so grateful for any help. I have been trying to work on this for awhile but I feel like my understanding is not where it needs to be and I definitely don't want to mess this up.