External Email - Use Caution        

Hi Tom,


regarding the first question:

I think that here the issue is that non-transformed and transformed p-values are mixed: the spval file should contain non-transformed values between p=0 and p=1, while the '1.3' (which is used as the threshold) should be the result of a -log10(p) transformation. 

Using the mri_surfcluster command with non-transformed values is possible, but a little bit tricky (and the tutorial is not too precise here): for example, one would need to set the range of values that are considered for clusters as [ 0 ... desired_p_threshold ] (keep all values below the threshold), as opposed to [ desired_transformed_p_threshold ... +Inf ] (keep all values above the threshold), which would be used for transformed p-values. Also, for non-transformed p-values the maximum statistics (Max and VtxMax) make little sense; one would rather be interested in the minima.This is, by the way, also how I interpret your mri_surfcluster output: it seems to me that there are a few very large clusters with a maximum value of p=1 - but you are not interested in those, I guess.

As a consequence, mri_surfcluster is (in my personal opinion) probably better used with -log10(p) - transformed p-values and the corresponding (transformed and FDR- or FDR2-corrected) threshold (sorry if my previous mail sent you in a different direction). 

In my eyes, you could use the 'sig.mgh' maps, which are also mentioned in the tutorial, for this purpose. This should work, because the lme_mass_FDR2 essentially just returns a new, adjusted threshold; it does not give 'adjusted p-values'. Therefore, thresholding the spval file with a p threshold (keeping all values below that threshold) should give the same clusters as the 'sig.mgh' file that is thresholded by a -log10(p) threshold (and keeping all values above that threshold). 

You could also compare the spval file to a map of uncorrected, non-transformed p-values that can also be written by the fs_write_fstats function (using the 'pval' argument instead of 'sig', type 'help fs_write_fstats' in Matlab for details). The pval map created by fs_write_fstats and the spval map should be identical; the only difference could be the sign that is attached to the values in the pval map (and in 'sig.mgh') to indicate the direction of the effect; the spval file, in contrast, should only contain absolute p-values.


regarding the second question:

That's difficult to tell, maybe changing the amount of smoothing could help, but I am speculating here.


Best regards,

Kersten

On Mi, 2019-07-24 at 21:16 +0100, tom parker wrote:

        External Email - Use Caution        

Dear Kersten,

Thank you so much for the quick reply! It was really helpful!!
I still have a couple of doubts regarding the visualization of the results if you don't mind.

After running the last command fs_write_Y(sided_pval,mri,'spval_lh.mgh'), I tried to visualize the results with tksurfer fsaverage lh inflated -overlay spval_lh.mgh -fthresh 1.3 (1.3 was the pth value that came out from lme_mass_FDR2) but there were no clusters on the surface.

If I run mri_surfcluster with spval_lh.mgh as input it seems there are 15 clusters, here is a sample of the first clusters that come out:
# ClusterNo  Max   VtxMax   Size(mm^2)  TalX   TalY   TalZ    NVtxs    WghtVtx
   1        1.000       8  22284.41     -2.2    8.6   -3.3  48988    46268.86
   2        1.000  140211   1542.44     -8.2  -16.7   58.6   3024     2633.36
   3        1.000   44065   1510.18    -16.3  -71.7   -4.6   2078     1810.43
  (...)

Am I doing something wrong?

I also tried the spatiotemporal mass-univariate approach that you suggested but matlab gets stuck with this command:
[lhRgs,lhRgMeans] = lme_mass_RgGrow(lhsphere,lhRe,lhTh0,lhcortex,2,95);

In the terminal it says:
"Starting region growing algorithm
Computing seeds ..."

Do you have any suggestions on how to solve this?

Thanks again for the help!


Hi Tom,

basically most of it looks OK to me, just a few remarks:

I assume that you deliberately chose the plain mass-univariate approach instead of the potentially more powerful spatiotemporal mass-univariate approach. 

For the interpretation of the intercepts, if desired, it may be advantageous to code sex as 0 and 1, and de-mean age. But that's not strictly necessary.

I suspect that it's just a typo, so just for the sake of completeness: for contrast estimation, I don't think that using CM.C1, CM.C2, etc. would work, as the program expects a 'C' substructure. So instead, one would use CM1.C, CM2.C etc for the different contrasts, if I am not mistaken.

About the mri_surfcluster command: if you want to do the FDR2 inference as suggested by the LME tutorial, just use the 'pth' value as returned by the lme_mass_FDR2 script for thresholding. I believe that you are mixing different kinds of inference (traditional FDR, clusterwise) otherwise. 

Finally, if you want to make corrected inferences for the whole cortex, it's probably better to compute one (FDR2) threshold for both hemispheres at once (see the very last part of the tutorial) instead of treating them separately.

And regarding your specific questions:

The order is fine.

You can visualize the 'spval_lh.mgh' by first loading e.g. the 'lh.inflated' surface of the fsaverage template, and then adding 'spval_lh.mgh' as an overlay. Same with e.g. 'rh.inflated' for the 'spval_rh.mgh' file.

About the Monte Carlo correction for this type of analysis others may know better, sorry.

Best regards,

Kersten

Missatge de tom parker <tomparker540@gmail.com> del dia dg., 21 de jul. 2019 a les 20:50:
Dear all,

I am currently running a Mass-Univariate LME analysis in Freesurfer. I have some doubts about the last steps and I was wondering whether someone could help me figure them out and also check what I have done so far (sorry for the long email).

My study has 2 time points (baseline, follow-up), 3 groups (controls, patients without cognitive problems, patients with cognitive problems) and 2 covariates of no interest (age, sex).
I am interested in looking at the interactions between group and time (mainly whether the 2 patient groups show more cortical thinning over time compared to controls).

These are the commands I have ran so far through terminal or in matlab:

1)
mris_preproc --qdec-long long.qdec.table.dat --target fsaverage --hemi lh --meas thickness --out lh.thickness.mgh
mris_preproc --qdec-long long.qdec.table.dat --target fsaverage --hemi rh --meas thickness --out rh.thickness.mgh

2)
mri_surf2surf --hemi lh --s fsaverage --sval lh.thickness.mgh --tval lh.thickness_sm15.mgh --fwhm-trg 15 --cortex  --noreshape
mri_surf2surf --hemi rh --s fsaverage --sval rh.thickness.mgh --tval rh.thickness_sm15.mgh --fwhm-trg 15 --cortex  --noreshape

3)
[Y,mri] = fs_read_Y('lh.thickness_sm15.mgh');
[Y,mri] = fs_read_Y('rh.thickness_sm15.mgh');

4)
lhsphere = fs_read_surf('fsaverage/surf/lh.sphere')
rhsphere = fs_read_surf('fsaverage/surf/rh.sphere')
lhcortex = fs_read_label('fsaverage/label/lh.cortex.label')
rhcortex = fs_read_label('fsaverage/label/rh.cortex.label')

5) Create X matrix in matlab, here is a sample:
intercept     time      group2     group2*time     group3     group3*time   age  sex
1                 0           0              0                      0              0                     41   1
1                 2.42      0              0                      0              0                     41   1
1                 0           1              0                      0              0                     53   1
1                 2.5        1              2.5                   0              0                     53   1
1                 0           0              0                      0              0                     54   2
1                 5.25      0              0                      0              0                     54   2

5)
lhstats = lme_mass_fit_vw(X, [1], Y, ni, lhcortex);
rhstats = lme_mass_fit_vw(X, [1], Y, ni, rhcortex);

6)
I entered 3 contrasts separately to check the interaction group*time.
CM.C1 = [0 0 0 1 0 0 0 0]; this contrast will show if group 1 and 2 show a different cortical thinning pattern over time
CM.C2 = [0 0 0 0 0 1 0 0]; this contrast will show if group 1 and 3 show a different cortical thinning pattern over time
CM.C3 = [0 0 0 -1 0 1 0 0]; this contrast will show if group 2 and 3 show a different cortical thinning pattern over time

7)
For each of the previous contrasts, i did the following:

F_lhstats = lme_mass_F(lhstats, CM);
F_rhstats = lme_mass_F(rhstats, CM);
fs_write_fstats(F_lhstats, mri,' sigL.mgh', 'sig');
fs_write_fstats(F_rhstats, mri,' sigR.mgh', 'sig');

mri.volsz(4) = 1
[detvtx,sided_pval,pth] = lme_mass_FDR2(F_lhstats.pval,F_lhstats.sgn,lhcortex,0.05,0);
[detvtx,sided_pval,pth] = lme_mass_FDR2(F_rhstats.pval,F_rhstats.sgn,rhcortex,0.05,0);

8)
fs_write_Y(sided_pval,mri,'spval_lh.mgh');
fs_write_Y(sided_pval,mri,'spval_rh.mgh');

mri_surfcluster --subject fsaverage --hemi lh --in spval_lh.mgh --cwpvalthresh 0.05 --fdr 0.05 --sign pos --o lh_time_cluster  --sum lh_time_cluster_sum
mri_surfcluster --subject fsaverage --hemi rh --in spval_rh.mgh --cwpvalthresh 0.05 --fdr 0.05 --sign pos --o rh_time_cluster  --sum rh_time_cluster_sum


My questions are:
- Is the order of the steps (1-8) correct?
- How can I visualize the significant FDR-corrected results?
- Is there any way I can apply a Monte Carlo correction using the generated files?

Thank you so much!


_______________________________________________
Freesurfer mailing list
Freesurfer@nmr.mgh.harvard.edu
https://mail.nmr.mgh.harvard.edu/mailman/listinfo/freesurfer