Friday, February 22, 2013

Using SPM.mat to Stay on Track (II)

Previously I discoursed at length upon the fatalistic nature of FMRI analysis, and how it leads to disorientation, depression, anxiety, and, finally, insanity. Those who are engaged in such a hopeless enterprise are easily identified by the voids that swirl like endless whorls within their eyes. I once knew a man driven out of his wits by just such an interminable cycle; he was found five days later on the sandy deserts of southern Indiana, ragged and filthy and babbling in alien tongues, scraping with pottery shards the boils that stippled his skin as though they had been painted on. Curse SPM, and die.

For those willing to prolong their agony, however, it is useful to know exactly what is contained within the SPM structure. Although not well-documented (if at all) on the SPM website, there are vast riches of information you can assay within the depths of the scripts themselves, which can be found in your spm8/spm5 library. (If they are in your path, you should be able to open them from anywhere; e.g., using a command such as "open spm_getSPM".) I have copied and pasted below the information from three such scripts, as well as a brief description about what they do. Good fortune to you.



================================================




1) spm_getSPM: Provides information about the xSPM structure which is updated each time Results are loaded and viewed


function [SPM,xSPM] = spm_getSPM(varargin)
% computes a specified and thresholded SPM/PPM following parameter estimation
% FORMAT [SPM,xSPM] = spm_getSPM;
% Query SPM in interactive mode.
%
% FORMAT [SPM,xSPM] = spm_getSPM(xSPM);
% Query SPM in batch mode. See below for a description of fields that may
% be present in xSPM input. Values for missing fields will be queried
% interactively.
%
% xSPM      - structure containing SPM, distribution & filtering details
% .swd      - SPM working directory - directory containing current SPM.mat
% .title    - title for comparison (string)
% .Z        - minimum of Statistics {filtered on u and k}
% .n        - conjunction number <= number of contrasts       
% .STAT     - distribution {Z, T, X, F or P}    
% .df       - degrees of freedom [df{interest}, df{residual}]
% .STATstr  - description string    
% .Ic       - indices of contrasts (in SPM.xCon)
% .Im       - indices of masking contrasts (in xCon)
% .pm       - p-value for masking (uncorrected)
% .Ex       - flag for exclusive or inclusive masking
% .u        - height threshold
% .k        - extent threshold {voxels}
% .XYZ      - location of voxels {voxel coords}
% .XYZmm    - location of voxels {mm}
% .S        - search Volume {voxels}
% .R        - search Volume {resels}
% .FWHM     - smoothness {voxels}    
% .M        - voxels -> mm matrix
% .iM       - mm -> voxels matrix
% .VOX      - voxel dimensions {mm} - column vector
% .DIM      - image dimensions {voxels} - column vector
% .Vspm     - Mapped statistic image(s)
% .Ps       - list of P values for voxels at SPM.xVol.XYZ (used by FDR)
% .thresDesc - description of height threshold (string)
%
% Required fields of SPM
%
% xVol   - structure containing details of volume analysed
%
% xX     - Design Matrix structure
%        - (see spm_spm.m for structure)
%
% xCon   - Contrast definitions structure array
%        - (see also spm_FcUtil.m for structure, rules & handling)
% .name  - Contrast name
% .STAT  - Statistic indicator character ('T', 'F' or 'P')
% .c     - Contrast weights (column vector contrasts)
% .X0    - Reduced design matrix data (spans design space under Ho)
%          Stored as coordinates in the orthogonal basis of xX.X from spm_sp
%          (Matrix in SPM99b)  Extract using X0 = spm_FcUtil('X0',...
% .iX0   - Indicates how contrast was specified:
%          If by columns for reduced design matrix then iX0 contains the
%          column indices. Otherwise, it's a string containing the
%          spm_FcUtil 'Set' action: Usually one of {'c','c+','X0'}
% .X1o   - Remaining design space data (X1o is orthogonal to X0)
%          Stored as coordinates in the orthogonal basis of xX.X from spm_sp
%          (Matrix in SPM99b)  Extract using X1o = spm_FcUtil('X1o',...
% .eidf  - Effective interest degrees of freedom (numerator df)
%        - Or effect-size threshold for Posterior probability
% .Vcon  - Name of contrast (for 'T's) or ESS (for 'F's) image
% .Vspm  - Name of SPM image
%
% Evaluated fields in xSPM (input)
%
% xSPM      - structure containing SPM, distribution & filtering details
% .swd      - SPM working directory - directory containing current SPM.mat
% .title    - title for comparison (string)
% .Ic       - indices of contrasts (in SPM.xCon)
% .Im       - indices of masking contrasts (in xCon)
% .pm       - p-value for masking (uncorrected)
% .Ex       - flag for exclusive or inclusive masking
% .u        - height threshold
% .k        - extent threshold {voxels}
% .thresDesc - description of height threshold (string)
%
% In addition, the xCon structure is updated. For newly evaluated
% contrasts, SPM images (spmT_????.{img,hdr}) are written, along with
% contrast (con_????.{img,hdr}) images for SPM{T}'s, or Extra
% Sum-of-Squares images (ess_????.{img,hdr}) for SPM{F}'s.



2) spm_spm: Command to estimate the general linear model for a dataset; contains information common to both 1st and 2nd-level SPM.mat files

function [SPM] = spm_spm(SPM)
% [Re]ML Estimation of a General Linear Model
% FORMAT [SPM] = spm_spm(SPM)
%
% required fields of SPM:
%
% xY.VY - nScan x 1 struct array of mapped image volumes
%         Images must have the same orientation, voxel size and data type
%       - Any scaling should have already been applied via the image handle
%         scalefactors.
%
% xX    - Structure containing design matrix information
%       - Required fields are:
%         xX.X      - Design matrix (raw, not temporally smoothed)
%         xX.name   - cellstr of parameter names corresponding to columns
%                     of design matrix
%       - Optional fields are:
%         xX.K      - cell of session-specific structures (see spm_filter)
%                   - Design & data are pre-multiplied by K
%                     (K*Y = K*X*beta + K*e)
%                   - Note that K should not smooth across block boundaries
%                   - defaults to speye(size(xX.X,1))
%         xX.W      - Optional whitening/weighting matrix used to give
%                     weighted least squares estimates (WLS). If not specified
%                     spm_spm will set this to whiten the data and render
%                     the OLS estimates maximum likelihood
%                     i.e. W*W' = inv(xVi.V).
%
% xVi   - structure describing intrinsic temporal non-sphericity
%       - required fields are:
%         xVi.Vi    - array of non-sphericity components
%                   - defaults to {speye(size(xX.X,1))} - i.i.d.
%                   - specifying a cell array of constraints (Qi)
%                     These constraints invoke spm_reml to estimate
%                     hyperparameters assuming V is constant over voxels.
%                     that provide a high precise estimate of xX.V
%       - Optional fields are:
%         xX.V      - Optional non-sphericity matrix.  Cov(e) = sigma^2*V
%                     If not specified spm_spm will compute this using
%                     a 1st pass to identify significant voxels over which
%                     to estimate V.  A 2nd pass is then used to re-estimate
%                     the parameters with WLS and save the ML estimates
%                     (unless xX.W is already specified)
%
% xM    - Structure containing masking information, or a simple column vector
%         of thresholds corresponding to the images in VY.
%       - If a structure, the required fields are:
%         xM.TH - nVar x nScan matrix of analysis thresholds, one per image
%         xM.I  - Implicit masking (0=>none, 1 => implicit zero/NaN mask)
%         xM.VM - struct array of mapped explicit mask image volumes
%       - (empty if no explicit masks)
%               - Explicit mask images are >0 for valid voxels to assess.
%               - Mask images can have any orientation, voxel size or data
%                 type. They are interpolated using nearest neighbour
%                 interpolation to the voxel locations of the data Y.
%       - Note that voxels with constant data (i.e. the same value across
%         scans) are also automatically masked out.
....
% The following SPM.fields are set by spm_spm (unless specified)
%
%     xVi.V      - estimated non-sphericity trace(V) = rank(V)
%     xVi.h      - hyperparameters  xVi.V = xVi.h(1)*xVi.Vi{1} + ...
%     xVi.Cy     - spatially whitened (used by ReML to estimate h)
%     xVi.CY     - <(Y - )*(Y - )'>    (used by spm_spm_Bayes)
%
%                           ----------------
%
%     Vbeta     - struct array of beta image handles (relative)
%     VResMS    - file struct of ResMS image handle  (relative)
%     VM        - file struct of Mask  image handle  (relative)
%
%                           ----------------
%     xX.W      - if not specified W*W' = inv(x.Vi.V)
%     xX.V      - V matrix (K*W*Vi*W'*K') = correlations after K*W is applied
%     xX.xKXs   - space structure for K*W*X, the 'filtered and whitened'
%                 design matrix
%               - given as spm_sp('Set',xX.K*xX.W*xX.X) - see spm_sp
%     xX.pKX    - pseudoinverse of K*W*X, computed by spm_sp
%     xX.Bcov   - xX.pKX*xX.V*xX.pKX - variance-covariance matrix of
%                 parameter estimates
%         (when multiplied by the voxel-specific hyperparameter ResMS)
%                 (of the parameter estimates. (ResSS/xX.trRV = ResMS)    )
%     xX.trRV   - trace of R*V, computed efficiently by spm_SpUtil
%     xX.trRVRV - trace of RVRV
%     xX.erdf   - effective residual degrees of freedom (trRV^2/trRVRV)
%     xX.nKX    - design matrix (xX.xKXs.X) scaled for display
%                 (see spm_DesMtx('sca',... for details)
%
%                           ----------------
%
%     xVol.M    - 4x4 voxel->mm transformation matrix
%     xVol.iM   - 4x4 mm->voxel transformation matrix
%     xVol.DIM  - image dimensions - column vector (in voxels)
%     xVol.XYZ  - 3 x S vector of in-mask voxel coordinates
%     xVol.S    - Lebesgue measure or volume       (in voxels)
%     xVol.R    - vector of resel counts           (in resels)
%     xVol.FWHM - Smoothness of components - FWHM, (in voxels)
%
%                           ----------------
%
% xCon  - See Christensen for details of F-contrasts.  These are specified
%         at the end of spm_spm after the non-sphericity V has been defined
%         or estimated. The fist contrast tests for all effects of interest
%         assuming xX.iB and xX.iG index confounding or nuisance effects.
%
%     xCon      - Contrast structure (created by spm_FcUtil.m)
%     xCon.name - Name of contrast
%     xCon.STAT - 'F', 'T' or 'P' - for F/T-contrast ('P' for PPMs)
%     xCon.c    - (F) Contrast weights
%     xCon.X0   - Reduced design matrix (spans design space under Ho)
%                 It is in the form of a matrix (spm99b) or the
%                 coordinates of this matrix in the orthogonal basis
%                 of xX.X defined in spm_sp.
%     xCon.iX0  - Indicates how contrast was specified:
%                 If by columns for reduced design matrix then iX0 contains the
%                 column indices. Otherwise, it's a string containing the
%                 spm_FcUtil 'Set' action: Usually one of {'c','c+','X0'}
%                 (Usually this is the input argument F_iX0.)
%     xCon.X1o  - Remaining design space (orthogonal to X0).
%                 It is in the form of a matrix (spm99b) or the
%                 coordinates of this matrix in the orthogonal basis
%                 of xX.X defined in spm_sp.
%     xCon.eidf - Effective interest degrees of freedom (numerator df)
%     xCon.Vcon - ...for handle of contrast/ESS image (empty at this stage)
%     xCon.Vspm - ...for handle of SPM image (empty at this stage)
%
%                           ----------------
%
%
% The following images are written to file
%
% mask.{img,hdr}                                   - analysis mask image
% 8-bit (uint8) image of zero-s & one's indicating which voxels were
% included in the analysis. This mask image is the intersection of the
% explicit, implicit and threshold masks specified in the xM argument.
% The XYZ matrix contains the voxel coordinates of all voxels in the
% analysis mask. The mask image is included for reference, but is not
% explicitly used by the results section.
%
%                           ----------------
%
% beta_????.{img,hdr}                                 - parameter images
% These are 16-bit (float) images of the parameter estimates. The image
% files are numbered according to the corresponding column of the
% design matrix. Voxels outside the analysis mask (mask.img) are given
% value NaN.
%
%                           ----------------
%
% ResMS.{img,hdr}                    - estimated residual variance image
% This is a 32-bit (double) image of the residual variance estimate.
% Voxels outside the analysis mask are given value NaN.
%
%                           ----------------
%
% RPV.{img,hdr}                      - estimated resels per voxel image
% This is a 32-bit (double) image of the RESELs per voxel estimate.
% Voxels outside the analysis mask are given value 0.  These images
% reflect the nonstationary aspects the spatial autocorrelations.




3) spm_fMRI_design: Contains information for 1st-level SPM.mat files (e.g., the basis function and session information for that GLM).

function [SPM] = spm_fMRI_design(SPM,save_SPM)
% Assembles a design for fMRI studies
% FORMAT [SPM] = spm_fMRI_design(SPM)
%
% 1st level
%---------------------------------------------------------------------------
% SPM.
%
%       xY: [1x1 struct] - data structure
%    nscan: [1xs double] - nscan(s) = number of scans in session s
%      xBF: [1x1 struct] - Basis function structure
%     Sess: [1xs struct] - Session structure array
%       xX: [1x1 struct] - Design matrix structure
%
%
%    2nd level
%    -----------------------------------------------------------------------
%    SPM.xY
%           RT: - repetition time {seconds)
%
%    SPM.xBF
%            T: - number of time bins per scan
%           T0: - first time bin (see slice timing)
%        UNITS: - 'scans'|'secs' - units in which onsets are specified
%     Volterra: - 1|2 - order of [Volterra] convolution
%           dt: - length of time bin {seconds}
%         name: - name of basis set
%       length: - support of basis set {seconds}
%        order: - order of basis set
%           bf: - basis set matrix
%
%    SPM.Sess(s)
%            U: - Input structure array
%            C: - User specified covariate structure
%          row: - scan   indices for session s
%          col: - effect indices for session s
%           Fc: - F Contrast information for input-specific effects
%
%    SPM.xX
%            X: - design matrix
%           iH: - vector of H partition (indicator variables) indices
%           iC: - vector of C partition (covariates)          indices
%           iB: - vector of B partition (block effects)          indices
%           iG: - vector of G partition (nuisance variables)  indices
%         name: - cellstr of names for design matrix columns
%
%
%        3rd level
%        -------------------------------------------------------------------
%        SPM.Sess(s).U
%               dt: - time bin length {seconds}
%             name: - {1 x j} cell of names for each input or cause
%              ons: - (q x 1) onsets for q  trials {in UNITS}
%              dur: - (q x 1) durations for trials {in UNITS}
%                P: - Parameter stucture
%                u: - (t x j) inputs or stimulus function matrix
%              pst: - (1 x k) peristimulus times (seconds)
%
%
%        SPM.Sess(s).C    
%
%                C: - [kx1 double] of user specified regressors
%             name: - {1xk} cellstr of regressor names
%
%
%        SPM.Sess(s).Fc    
%
%                i: - F Contrast colums for input-specific effects
%             name: - F Contrast names  for input-specific effects
%
%
%            4th level
%            ---------------------------------------------------------------
%            SPM.Sess(s).U(i).P(p)
%
%                 name: - parameter name
%                    P: - (q x 1) parameter matrix
%                    h: - order of polynomial expansion (0 = none)
%                    i: - sub-indices of U(i).u for plotting


Thanks to Will Moore. When he's right, he's right; when he's wrong, he's dead wrong. But when he's right...he's really right.

10 comments:

  1. Life saving post!

    ReplyDelete
  2. Hi, how do I create a SPM.mat file?
    I have the time series already from a DCM file, but I want to create other model.
    Thank you

    ReplyDelete
  3. Hi Andy!
    I really need your help. Once I get the SPM.mat how can I extract the fitted time series of the BOLD signal for ALL the voxels in the activation map? I guess in theory I can do it at a ROI level (don't know how) and I know I can do it for a single voxel by plotting the 'fitted response' and selecting the contrast I choose and this would generate a 'y' variable which corresponds to the fitting response. How do I get this fitted response for all voxels? Also, if a voxel is active in different contrasts, why is the fitting function different?
    Thank you in advance for your help!

    ReplyDelete
    Replies
    1. Hi Jess,

      Someone recently asked me a similar question, and I recommended using AFNI's 3dmaskdump tool. I've reprinted the email below:


      I don't know how to do that in SPM, but you can do it in AFNI using the command 3dmaskdump.

      First, concatenate all of your .hdr files into a single dataset using 3dTcat; e.g., 3dTcat -prefix myRun.nii swar*.hdr

      Then, use your mask to extract the timecourse from each voxel:

      3dmaskdump -mask mask.nii -noijk -xyz myRun.nii > allTimeCourses.txt

      This will print the timecourses to the file allTimeCourses.txt. Each row will be a timecourse from each voxel, and the number of columns will be the number of volumes in your dataset. You can import this into Excel for further manipulation or statistics. Check the help of 3dmaskdump to learn more about the options -noijk and -xyz.


      =========

      The only problem is that I don't know where SPM keeps its fitted timecourse; unlike AFNI, which has the option of generating a separate dataset with the fitted timecourse at each voxel, SPM does not create a separate file. How it calculates the fitted timecourse in the GUI, I do not know.

      -Andy

      Delete
    2. Thank you for your help! The problem is that I don't want the initial fMRI time series which are encoded in the fMRI images, I would like to have the fitted responses (after preprocessing and after determining the contrast and pvalue) which as you said are not encoded in any separate file so I wouldn't know which .hdr file to concatenate since the images generated are all 3D. I assume that since you can extract the information for a single voxel, in theory you could extract the same information for all voxels (maybe writing a matlab script). If you have any ideas please let me know; my final goal is to calculate the pair correlations of the timeseries of the activated voxels, this is why I need the fitted responses instead of the raw signal.
      Thank you again
      Jess

      Delete
    3. I don't know how SPM generates the fitted time series. I'm sure you could find out by looking at how SPM's .m files create it, but unfortunately I don't know where to direct you to.

      Good luck!

      -Andy

      Delete
    4. Thank you again for your help! I will hopefully figure it out -

      Delete
  4. Hi Andy,
    Thanks a lot for your always helpful tutorials on fMRI!
    I have done all the preprocessing and 3Ddeconvole using AFNI. Now I need do DCM using SPM. So pb05.sub.scale+tlrc and roi.mask+tlrc will be transformed to nifti format using 3dAFNItoNIFTI and the voi.mat will be generated later using SPM functions. The time series of first scaled eigenvariates of VOI will be inputted to DCM functions. But now I have two questions:
    1.The DCM in SPM also needs the input 'SPM.Sess.U.dt' and 'SPM.Sess.U.u'.
    Do you know how to create SPM.Sess.U.dt and SPM.Sess.U(1).u by myself without using GLM function in SPM?
    2.Timing slices has already been done in AFNI preprocessing and all slices have already been aligned to the first slice. And I have not too much time for repeating the preprocessing which omits ‘slice-timing correction’. So I will set the timing slice input 'DCM.delays' to zeros because the paper (which has been attached following) says that the DCM models whether including timing slice or not doesn't matter too much for very short TR ( TR 2s in my experiment). Do you think I will do something wrong with the timing slices for DCM?
    Your answers will be very appreciated!
    Po
    (Kiebel, S. J., Klöppel, S., Weiskopf, N., & Friston, K. J. (2007). Dynamic causal modeling: a generative model of slice timing in fMRI. Neuroimage, 34(4), 1487-1496.)

    ReplyDelete
    Replies
    1. Hi Po,

      I haven't done DCM analysis before, so I don't know all the details. That said, here's my best guess about the answers to your questions:

      1. You could try creating them manually by creating them as structures within SPM and assigning variables to the fields that you need. My feeling is that SPM may look for the entire structure, not just the fields you specified, so you may need to populate the rest of the fields with null values.

      2. I can't say how much the slice timing will affect your DCM analysis. If you set the delays to zero, which sounds as though it will prevent DCM from doing any additional slice timing, then I don't see how that will be problem. It could be that SPM looks for the DCM.delays field later on and uses it for the actual model comparison, so I would look in the scripts to see if that's the case.

      -Andy

      Delete
    2. Hello Andy,
      Thanks a lot for your patient answers! I am very sorry to trouble you but I have another general question for beginners like me:
      How can we know which specific functions involved in a preprocessing step in SPM?
      Although we can define structures and conduct these structures using the function spm_jobman to preprocess data of many subjects. We still don't know the specific Matlab function behind a preprocessing step. For example, we all know the structure jobs{1}.stats{1}.fmri_spec is in charge of GLM in SPM but we don't know which functions are in charge of GLM. Now I want to know how to find out the specific Matlab functions used in a preprocessing step.
      Thank you in advance!
      Po

      Delete