Hello Spharm experts, 

I am currently running shapeAnalysisMANCOVA on my hippocampal data. Everything is going along smoothly in my examination of a statistically significant main effect of group (diagnosis) and I have just a few questions on how to move forward.

My input .csv file has the following variables/headings:
diagnosis,scaling,inputfile,sex,totalbrainvolume,age

...and I am running the following command to obtain a significant main effect of group, while accounting for sex, totalbrainvolume and age:

shapeAnalysisMANCOVA input.csv --infileColumn 2 --numGroupTypes 2 --numIndependent 2 --columnGroupTypes 0,3 --columnIndependent 4,5

My questions are as follows:

1. how can I explicitly model the interaction of diagnosis and sex, while controlling for totalbrainvolume and age?

2. is there a singular overall p value that is produced by shapeAnalysisMANCOVA that would be similar to the significance of an omnibus test? I'm a bit confused as to how correction for multiple comparisons is performed. I see that FDRP and Bonferroni text files are produced. But can one report that a statistically significant difference exists if the clusters in the FDRP and Bonferroni files are exceptionally small?

3. In visualizing the magnitude of structural differences between two groups, it's my understanding that I should load and view the *_meanAll.vtk file which happens to have various attributes already embedded within it (including "DiffMagnitude" and "DiffMagnitude_original"). The DiffMagnitude_original attribute has values ranging from -3.22 to 3.52, whereas the "DiffMagnitude" has a range of 0 to 100. Can you tell me what these values indicate (including their units)?

4. do I need to mean center my continuous covariates prior to running shapeAnalysisMANCOVA? I see that this question was previously raised by another user, but the answer was not clear.

Many thanks in advance for the help!
Lara

Hello again, 

Reposting to see if users had suggestions/answers to my questions, especially question #1:

how can I explicitly model the interaction of diagnosis and sex, while controlling for total brain volume and age?

Many thanks for any advice,
Lara

Hi Lara,

From my understanding:

1. I'm not certain about your interaction test, but it should be possible (from this post https://www.nitrc.org/forum/message.php?...). I believe you could set up the interaction test with the additional flags --interactionTest --simpleCorrs, but I think you'll have the interaction between brainvolume and age together. You could perhaps use brainvolume as the ICV correction factor, which in that case shouldn't get included in the interaction test. You'll definitely get an output, I'm just not certain if it's the model you're after! Beatriz or Martin will need to confirm.

2. No overall single value is produced in the MANCOVA tool, although there was in the previous one (StatNonParam, but this didn't allow covariates). But a single value isn't strictly necessary.

The MANCOVA occurs between each point (The X, Y and Z coordinates form the multiple dependant variables), and this is repeated for all points (1002 times if you're using the standard settings). So the correction for multiple comparisons is done over all those points, because with 1002 comparisons 50 of those may have a p value less than .05 by chance alone. Correction for multiple comparisons takes all these into account. The Bonferroni is the standard .05/1002 which is very conservative. FDR is less so. In my laymans terms, I think of the areas that remain significant after correction (despite how small) are the areas where the significant difference between groups is "so significant" that it is beyond a chance occurrence.

3. DiffMagnitude_original should be in millimeters (but I think this may change if your original MRI wasn't in 1x1x1mm spacing). pvalue_original is the actual p_value. The ones which range from 0 to 100 are just re-scaled. I use the *_original for clarity's sake.

4. I haven't mean centered. From my understanding this is necessary if interested in parameter estimates around the intercept, which I'm not (mainly interested in the significance maps, which centering doesn't change).

I'm very keen to hear Beatriz or Martin's take as well.

Regards,
David.

Hi Lara,



Question 1:

shapeAnalysisMANCOVA allows for including variables that ICV correct any measurements, given that you know total brain volume.
(in the wizard, in the test parameters tab, under the "reinterprets the scaling column values etc...") is where you would click tthis option to have your data ICV corrected. For correcting for age you will have to list them as independent variable. For testing diagnosis and sex, probably you could model them using a set of binary variables similar to what we propose here (http://www.insight-journal.org/browse/pu...)
Examples of commandlines are also provided in the aforementioned paper.



Question 2:

I think it would be safe to report significant shape changes, even if the FDR or Bonferroni (especially with Bonferroni, since it is a much more strict correction test) result into small clusters of significance. I would actually look at how those clusters look in the raw P maps as well.



Question 3:

The differences between DiffMagnitude_original attribute and "DiffMagnitude" (range 0-100) are sadly a result of having no maintenance support for the module. The released and discontinued version of SAM that you are using, generates MRML 3DSlicer scenes for visualization of statistical results of this program in Slicer. Slicer back in the years when we developed this (2010) did not have a good visualization scheme for continuous variables (such as structural differences). Slicer allowed though visualizing color labels (that they had capped at 100) so in order to use that infrastructure we had to rescale our continuous variable maps from 0-100. DiffMagnitude is that scaled map. Your real map and the one you want to report results with is DiffMagnitude_original. Sorry about the confusion. We are working on gathering funding for a new shape statistical package and we will keep you posted, as well as have all these bugs in mind for correction.



I hope that helps!
Thanks,


Beatriz

Thank you, Beatriz and David, for your help with this!

I will follow your suggestion of creating 4 binary, mutually exclusive variables to model the interaction of diagnosis and sex. This would change the my input .csv file from having the following variables/headings:

diagnosis,scaling,inputfile,sex,totalbrainvolume,age

to having the following variable headings instead:

control_male, control_female, patient_male, patient_female, inputfile, totalbrainvolume, age

So in setting up the command line for usage of shapeAnalysisMANCOVA, would the following be appropriate (with no --interactionTest flag)?

shapeAnalysisMANCOVA input.csv --infileColumn 4 --numGroupTypes 4 --numIndependent 2 --columnGroupTypes 0,1,2,3 --columnIndependent 5,6

This would then produce significance maps that visualize areas having a significant interaction of diagnosis and sex, correct? My plan would then be to decompose significant interactions using pairwise comparisons , through screening that show both significant effects in the interaction and in the pairwise comparison. Please know if this reasoning is correct and thanks so much again for the help!

Lara

Hi again, 

In playing with different ways of setting up a GLM in SPSS, I see that another (more straightforward?) way of modeling the interaction is to code a single grouping variable has having one of four values (e.g., 0 if male control, 1, if female control, 2 if male patient, 3 if female patient), then testing this main effect of group to get at the interaction of sex and diagnosis. Is this correct?

Many thanks again for the help,
Lara

Hi Lara,
 
I think the SPHARM MANOVA tool is only built to handle 2 groups within a single variable, so having 4 groups in the same variable will ignore the last 2.
 
David.

Indeed it will only be able use binary group variables. Not sure what it will do if you provide more
Martin

Okay. So then Beatriz's suggestion (of creating 4 binary, mutually exclusive variables to model diagnosis and sex) is not possible then either?

In this case, my input .csv file now has the following column headers:
control_male, control_female, patient_male, patient_female, inputfile, totalbrainvolume, age

I have entered the following command:
shapeAnalysisMANCOVA input.csv --infileColumn 4 --numGroupTypes 2 --numIndependent 2 --columnGroupTypes 0,1,2,3 --columnIndependent 5,6 --testColumn 0,1,2,3

...but get the error: "PARSE ERROR: Argument: -t (--testColumn). Couldn't read argument value from string '0,1,2,3'

Thanks very much,
Lara

Hi Lara

if you want a gender * group interaction, you will have to exactly model that, and thus you would need 3 group variables:
a) gender: (e.g. male = 0, female = 1)
b) group: (e.g. control = 0, group = 1)
c) group * gender: simply multiply a * b

you can still have your other independent effects.

a) captures the pure gender effect, b) the pure group effect and c) their interaction.

You can then test for whichever you are interested in. If you are interested to find out if there is a significant group by gender interaction, then test for column (c)

Martin

Hi Martin, 

Thanks so much. That worked!

Some additional minor questions for you:

1. In viewing the FDR corrected results in ShapePopulationViewer, should one set the upper threshold limit to be equal to the value that is printed in the terminal window at the end of shapeAnalysisMANCOVA? For example, I see the following:
...
No statistic type selected, defaulting to Hotelling.
Now running nonparametric tests:
Performing MANCOVA permutation testing.
Using 10000 permutations.
Completed MANCOVA permutation testing.
fdr thresh (MANCOVA) is = 0.0477
Runtime was 14 minute(s) and 2 second(s).
... 
would I therefore set the upper limit to 0.0477 in SPV? What about for a bonferroni correction?

2. if the *_mancovaFDRP.txt file contains all 1's, does that mean that there are simply no significant effects following FDR correction? 

3. Should brain volume be accounted for in the model using the scaling input/flags (--scale and --scaleColumn)? or by entering this as an independent variable (with input/flags --numIndependent and --columnIndependent)? I've modeled my data both ways and see that my results vary substantially depending on whether TBV is entered as a scale factor or as a continuous covariate.

4. Lastly, should I mean center my continuous covariates? Again, my results vary quite a bit depending on whether I demean.

Many thanks for the continued help. I'm very excited to have gotten this far and couldn't have made it through my analysis without the support of your group!

Lara

Hi Lara,


I am not sure about 1). I would think this is not how you want to show your results though, and display in SPV the user inputed value for significance instead, maybe Martin can provide more information.

Yes, that is correct for 2)

I believe that if you want to ICV correct your data is just going to be introducing a scaling factor to your data, if you model them as covariates, my understanding is that they will be used in the GLM trace to compute significance. Therefore the processing is different.

Maybe centering your continuous variables makes sense if the magnitude of each different variable is not comparable to each other. What are your continuous variables? What are the ranges for those?


Lara, we are happy to be of help. Let us know how it goes.
Best,


Beatriz

Hi Lara
Re 1) the fdr maps are corrected, so if you want a 5% threshold then set the colormaps in SPV such that everything above 5% shows a non-significance color. The output on the terminal is not relevant for any settings in SPV.

Re 3) Just like Bea says the processing is different. If you scale with ICV/TBV, each individual measurement is corrected with its ICV perfectly accounting for individual head/skull size or brain volume (and thus highly non-linear as viewed across the population). Using it as a covariate, the GLM fits a linear interaction and thus you "correct" your shape data linearly for head/brain size. It all depends on what you want the interaction between shape and brain size to be. 

The norm in the field is linear correction, but whether that is the "best/correct" thing to do is unclear. In our own studies, we currently normalize linearly via the GLM.

Re centering) shapeAnalysisMANCOVA currently does not center/normalize the covariate data. Best practice would be to center the data, particularly if you have multiple covariates that have vastly different ranges. 

Martin

Hi Martin and Beatriz,

Thank you very much for the clear and helpful feedback. Everything is coming along smoothly!

I am now visualizing group differences in shape alterations (magnitude of differences) in ShapePopulationViewer. I see in your 2009 Insight Journal publication (http://www.insight-journal.org/browse/pu...) that the difference maps show negative values if group A has shrinkage relative to group B, and positive values if group A has expansion relative to group B.

My question is: how can I confirm which group is group A in my study, and which is group B?

I see that, in the terminal in which I'm running shapeAnalysisMANCOVA, there is a printout that reads:
"data in group_type 0 has been relabeled: 1 --> group A = -1 ; 0 --> group B = 1
#(A)= 142; #(B)= 64
data in group_type 1 has been relabeled: 0 --> group A = -1 ; 1 --> group B = 1
#(A)= 250; #(B)= 162"
...which offers some hints. But given my analysis setup, I'm a bit confused as to what this means. Recall that my inputfile.csv file has the following headers:

diagnosis,sex,diagnosis_sex_interaction,Scaling,File,tbv,tbv_dm,age,age_dm

...To inspect a main effect of sex, controlling for tbv_dm and age_dm, with sex, diagnosis and the interaction of sex and diagnosis in the model, my shape command is as follows:

shapeAnalysisMANCOVA inputfile.csv --infileColumn 4 --numGroupTypes 2 --numIndependent 2 --columnGroupTypes 0,1,2 --columnIndependent 6,8 --testColumn 1

From the printout in the terminal, I'm inferring that re: "diagnosis", group A = individuals with a "1" in the diagnosis column, and group B = individuals with a "0" in the diagnosis column. But what about sex? Is there a way I can tell whether individuals with a "1" in the sex column are in group A or group B? Likewise for the "diagnosis_sex_interaction" column.

Many thanks for the continued help!
Lara

Hi Lara,



Happy to hear there has been progress.
A is the group that you labeled as "1"for grouptype 0 and grouptype 1 value "0" in your CSV and B is the group you labeled as "0" in your CSV for grouptype=0 and "1" for group type 1. SAM will use A for the first label found in the CSV (you can verify your first row in the CSV is a "1" case). The *_DiffMagnitde.txt file will show signed differences from A to B (+ value means B is bigger than A, and - value means B is smaller than A). You can use any vtk visualization tool (Paraview, Slicer) to make sure this is correct.


Could you please share few rows of your CSV with us in order I can think about this a little further?

Thanks,
Beatriz

If you load *_meanA and *_meanB you should see the shape differences captured by *_DiffMagnitude.txt.

Thanks so much for the quick response. Here are the first few rows:

diagnosis,sex,diagnosis_sex_interaction,Scaling,File,tbv,tbv_dm,age,age_dm
1,0,0,1,/media/sf_Spharm/14753.1/spharm/14753.1_Lhipp_surfSPHARM_ellalign.vtk,0.1618064,0.007095606,5.61,-1.384223301
0,0,0,1,/media/sf_Spharm/15397.1/spharm/15397.1_Lhipp_surfSPHARM_ellalign.vtk,0.1607335,0.006022706,6.75,-0.244223301
1,0,0,1,/media/sf_Spharm/15403.1/spharm/15403.1_Lhipp_surfSPHARM_ellalign.vtk,0.1522853,-0.002425494,7.82,0.825776699
0,0,0,1,/media/sf_Spharm/15407.1/spharm/15407.1_Lhipp_surfSPHARM_ellalign.vtk,0.1629996,0.008288806,4.86,-2.134223301
1,0,0,1,/media/sf_Spharm/15433.1/spharm/15433.1_Lhipp_surfSPHARM_ellalign.vtk,0.1516102,-0.003100594,5.15,-1.844223301
0,0,0,1,/media/sf_Spharm/15461.1/spharm/15461.1_Lhipp_surfSPHARM_ellalign.vtk,0.1483783,-0.006332494,7.37,0.375776699
1,0,0,1,/media/sf_Spharm/15462.1/spharm/15462.1_Lhipp_surfSPHARM_ellalign.vtk,0.1542314,-0.000479394,6.91,-0.084223301
1,0,0,1,/media/sf_Spharm/15467.1/spharm/15467.1_Lhipp_surfSPHARM_ellalign.vtk,0.1723395,0.017628706,6.26,-0.734223301
1,0,0,1,/media/sf_Spharm/15473.1/spharm/15473.1_Lhipp_surfSPHARM_ellalign.vtk,0.1689724,0.014261606,6.11,-0.884223301

Thanks again,
Lara

Hi Beatriz and Martin, 

Following up on my last post - can you please confirm, with respect to my input .csv file formatted as specified in my last post, that the correct command to model the interaction of diagnosis and sex, controlling for age and total brain volume is the following:

shapeAnalysisMANCOVA inputfile.csv --infileColumn 4 --numGroupTypes 2 --numIndependent 2 --columnGroupTypes 0,1,2 --columnIndependent 6,8 --testColumn 2

Specifically, I would like to know whether the value for the --numGroupTypes is correct if I have 3 separate grouping columns (diagnosis, sex, diagnosis*sex). I am assuming that --columnGroupTypes takes care of the number of grouping columns, and --numGroupTypes is where one indicates how many groups are referred to within each of the columns of c--columnGroupTypes. Is this correct?

Many thanks in advance,
Lara

Hello Spharmers,

I hope everyone had a nice holiday. Just checking in to see if I could get some feedback on my last post?

Many thanks!
Lara

Hi Lara,



Sorry it has taken so long. I actually had to dig this information back from the code.

My personal way to figure out the direction of the differences of the P-value maps is by looking at meanA and meanB combined morphologies (in Slicer, for example) and making comparisons with the P-value maps.,Mean A and mean B are created by averaging each files of each group for the first group type (in your case "diagnosis"), and I dont think the program has any special way to create meanA and meanB. You could use MeshMath distributed with SPHARM to average your files according to group_type 0 (diagnosis) and group_type 1 (gender) and see what your actual meanA and meanB look like compared to those. I am sure that computing p-values is done with the GLM including all variables for testing.

Another way to check directionality of differences, is the DiffMagnitude file, but even so I am not so sure what that files means in your case, since it is just an arithmetic subtraction of meanB to meanA so I dont think its something you can use to figure out the directionality of the changes.


I am sorry I am still unsure about this question. I will check with Martin and get back to you asap.


Thank you,

Beatriz

Hi Lara,



Martin and I discussed your problem today. 

First let us highlight that the problems you are facing to interpret your results are one of the main limitations of this software and the main reason why we have decided to discontinue its development and rather work on an alternative that can provide more meaningful results. This alternative solution would provide beta maps that could make interpretation easier (among other features).

Second, after careful evaluation of the information we have about your experiment, we think you are doing the things right (or as right as possible given the known limitations of the program).

Although imperfect, one can interpret the raw mean maps (MeanDiff) as I explained before, to understand the directionality of your p-values, where negative differences between meanA and meanB will mean meanB is smaller and positive differences will mean meanB is larger. This does not work in your case because you dont have 2 groups but 4 (G1: diag=0 and sex=0, G2: diag=0 and sex=1, G3: diag=1 sex=0, G4: diag=1 and sex=1). In the way the program is encoded only one of the meanA and meanB will be correct, since one will encode diag=1 and sex=1 meaning diag*sex=1 and the other will encode all other groups or diag*sex=0.

Our suggestion is that you generate means and shape difference maps (your own DiffMagnitude maps) that can help you interpret the directionality of your p-value maps better for G3: diag=1 sex=0, G4: diag=1 and sex=1 yourself. After analyzing how all the different Betas in the GLM will interact with each other in the experiment you describe in your commandline, we believe this would be the best thing to do.


SPHARM-PDM includes some commands that will help you do that. Here is the recipe:

1. Generate mean for G3 and G4. For that you will use MeshMath (included in the SPHARM-PDM distribution) 

Commandline: MeshMath subject1_G3.vtk average_output_G3.vtk -avgMesh subject2_G3.vtk ... subjectN_G3.vtk (repeat this for all subjects in G4)

2. Compute distances between avgG3 and avgG4. For that we need to first convert the format of your newly created average_output_G3.vtk and average_output_G4.vtk to Meta format. VTK2Meta is also part of SPHARM.

Commandline: VTK2Meta average_output_G4.vtk average_output_G4.meta (repeat for G3)

Now you compute differences between your two averages using MeshMath:

Commandline: MeshMath average_output_G4.meta -subtract average_output_G3.meta (this will do G3 - G4, meaning if positive differences are computed G4 is bigger than G3 and if negative differences are computed G4 is smaller than G3, you can change the direction by switching the order of your average files in your commandline) 

You will have to rename your output vectors file (mv ./-subtract vectors_G3_to_G4.txt) and also map them into your mean (I usually chose the reference average shape, for example if A-B I will map into A) for visualization in ShapePopulationViewer (MeshMath average_output_G3.vtk average_output_G3_with_attributes.vtk -KWtoPolyData vectors_G3_to_G4.txt vectorMap).

Next step in computation of differences is to compute the signed differences using also MeshMath.

Commandline: MeshMath average_output_G4.meta -magNormDir vectors_G3_to_G4.txt 

Again, you will have to rename your signed distances file (mv ./-magNormDir diffmagnitude_G3_to_G4.txt) and also map them into your mean for visualization in ShapePopulationViewer (MeshMath average_output_G3_with_attributes.vtk average_output_G3_with_attributes.vtk -KWtoPolyData diffmagnitude_G3_to_G4.txt DiffMagnitude).


Load average_output_G3_with_attributes.vtk in ShapePopulationViewer to see a) the vector map named vectorMap, and b) the signed differences in the map named DiffMagnitude. You can choose to change the name of your files or maps to whatever makes sense to you. Finally, let me highlight it is very important to keep the order consistent if you switch the order of the groups to visualize different signs in your maps.



Please, let me know if this makes sense as it comprises a lot of information. I tried to simplify the explanations to make it coherent for you, but if you want me to I can further explain anything that is not clear enough.



I hope that helps.


Best regards,
Beatriz

Thank you so much for your reply. As I read it now, it makes perfect sense. I will give this a go and let you know if I run into problems.

Many thanks again!
Lara

Hello everyone,

I am going to reply to a question Lara submitted in this forum, since it contains information that might be interesting for other users:
I think we are talking the same thing, but with two different approaches. I understand your point of view, but let me explain mine a little better.
Lets say we have the following 2 structures (figure 1), from which green is meanA and semitransparent red is meanB, and I want to see how meanB changes with respect to meanA (distances computed from A to B, so order in MeshMath is MeshMath meanB -subtract meanA).

Figure 1: https://drive.google.com/open?id=0B06eo8efbhotZTl4STBzeWRObHc

After we do -subtract operation in MeshMath in that orders our vectors have the following form (figure 2):

Figure 2: https://drive.google.com/open?id=0B06eo8efbhotRklJSVZtS25ab2s

That makes sense, right? our vectors will be pointing from meanA to meanB. They are a 3D data map that should be displayed on top of meanA. This is what your VectorMaps corresponds to.

Now, to signed distances. The way MeshMath calculates the sign of the distances (all distances are by definition positive) is calculating if the distances are coming out or into the reference mesh (meanA). MeshMath uses the normal vectors of meanA ([url=https://en.wikipedia.org/wiki/Normal_(geometry)]https://en.wikipedia.org/wiki/Normal_(geometry)[/url]). If the geometry of meanA is right, normal vectors should be pointing outwards the shape in each of its triangles. So we compare our difference vectors with the normal vectors of meanA and if they match the distance will be positive (both pointing outwards) if they don't the distance will be negative (diff vector points inwards, normal points outwards).
For the example above, the signed distances are (figure 3):

Figure 3: https://drive.google.com/open?id=0B06eo8efbhotX3pKMV9YQjUwRjQ

Where blue is positive and red is negative. Now compare with the first picture. Areas in which meanB (red) is bigger than meanA (green) have positive sigdistances. Areas in which meanB is smaller than meanA have negative distances. This is what DiffMagnitude should be showing in your case.

Now, for how VectorMaps and DiffMagnitude look like, of course they should look different, they are different things: VectorMaps are 3D maps, DiffMagnitude are 1D maps. In order to understand your problem a little better, could you please tell me what program are you using to display them?

I hope all the information above makes sense. Let me know if it doesn't, I am happy to clarify.
Thanks,


Beatriz