Dear FreeSurfer users,
I am using R to apply a statistical method on surface data registered to fsaverage (for example, left hemisphere area smoothed at bandwidth=10mm). This method requires data be in a matrix Y with NS rows and V columns, where NS = #subjects and V=#vertices. The output produces a statistic Sj > 0 at each vertex j that, without loss of generality, rejects the null hypothesis for higher values with significance derived through permutation tests. Let's assume that Sj is a F-statistic with J & NS-NB-J degrees of freedom in a GLM with design matrix X (NS rows and NB columns) and contrast matrix CON (J rows and NB columns). While mri_glmfit can do the same for F-tests, my workflow uses matrix operations to produce the same F-statistics as mri_glmfit but considerably faster.
My purpose is to apply a statistics method not supported by FreeSurfer and also use GPUs for permutation tests (since these methods rely on matrix operations). I am asking for the best way to manage the data (observed and permutation null distributions) so I can perform spatial cluster corrections using tools from FreeSurfer.
My output will be a vector S = (S_1,...,S_V) containing observed statistics at the V vertices indexed in the same order as the input data matrix Y containing surface data for each subject. Performing Q permutation tests will produce matrix R with V rows that contain values of the statistic R_{j1}, ...,R_{jQ} at vertex j (j=1,...V) under the null distribution. P-value at a given vertex could be:
P_j := (1/(Q+1)) * \sum_{q=1}^Q I( S_j < R_{jq})
I can do all the computations described above (and any conversions, like p-values to -log10 scale).
Does FreeSurfer offer an efficient way of taking the vector S of observed statistics and matrix R of permutation distributions (or P-values on any scale) to identify significant clusters on the surface at some threshold? Thresholds could be applied equivalently to either P-values (P_j < P* := 0.001) or statistic value (S_j > S*). Significance of observed clusters would be evaluated under null distribution of the largest cluster from the permutation tests (at a given cluster defining threshold).
Please forgive my unfamiliarity with how spatial information is dealt with in mgh files.
I am also trying to make a cluster summary report (for example, the file cache.th30.abs.sig.cluster.summary when using cached simulations for a contrast in a GLM) and and also visualize significant clusters.
I can export or analyze this data using any matlab function in the ~/freesurfer/matlab folder ("save_mgh" for example). For exporting, I'm unsure about how to preserve or define appropriate header information with save_mgh.
To clarify, my workflow contains following steps:
1. mri_preproc: Produce DATA.mgh containing concatenated lh-area-fsaverage-fwhm10 for subjects 1,...,NS 2. Use "load_mgh" in Matlab to load DATA.mgh into matrix Y0 of dimension [V, 1, 1, NS], where V=#vertices. 3. Define matrix Y of dimension [NS, V], which is transpose of squeeze(Y0(:,1,1,:)). 4. Apply algorithm/test to matrix Y to obtain vector S of observed statistics at vertices and matrix R rows containing null distributions respective vertices. Indexing of vertices for entries in S and columns of R are in the same order as first dimension of matrix Y0 loaded using "load_mgh".
Thank you for your attention to this very long question!
Best,
Chintan