cwatson / brainGraph

Graph theory analysis of brain MRI data
166 stars 48 forks source link

Anova tests between three or more groups of graph descriptors #11

Closed giova1588 closed 5 years ago

giova1588 commented 5 years ago

Hi everyone, i have a question about structural covariance matrices analysis. Following the manual i built adjacency matrices and graphs for my structural data. At this point I would like performing statistical analysis between groups at graph and vertex level, but i have three groups (i.e 1st patient group, 2nd patient group, 3 control group). How could i perform Statistics analysis in this situation (ANOVA)? The permutation analysis I read in the manual does a comparison only between two groups. An other question is if it is possible choosing the groups to compare in the permutation analysis. Thanks

cwatson commented 5 years ago

Yes, if you only have a single graph per group (as is common with structural covariance networks), then you will not be able to do something like an ANOVA. This is because you only would have 3 values for a network metric, instead of n values, 1 for each subject. To select only 2 of your 3 groups, you could subset the residuals object in the call to brainGraph_permute; for example, brainGraph_permute(densities, myResids[, c('Control', 'Patient1')], ...)

replacing with the correct group names and the correct name of your residuals object.

giova1588 commented 5 years ago

Ok, thanks. So I could perform the following analysis in the way that you describe: Patient1 vs Patient2, Patient1 vs Control, Patient2 vs Control, and there are no possibility to include all three groups (Patient1 vs Patient2 vs Control) to evaluate the differences between all groups like ANOVA. But statistically doing this type of analysis is right? Because in this way is like doing three t-test increasing the probability of type I error. How could i control for it? I should apply a multiple comparison correction? Addendum: i have tried doing the analysis as you suggested but it doesn't work. The three output (one for each comparison) are exactly the same.

cwatson commented 5 years ago

Can you post the exact commands you used?

giova1588 commented 5 years ago

This is an example: perms.all.graph_VHvsNVH <- brainGraph_permute(densities, all.dat.resids[, c("DLB_VH", "DLB_NVH")], perms = myPerms, level = "graph", atlas = atlas) all.dat.resid is a brainGraph_resids object (i have done this one through the get.resid function and it contains a variable named Groups with my three group labels).

cwatson commented 5 years ago

Ok, please send me some more information:

  1. What version of brainGraph are you using?
  2. What is the output of str(all.dat.resids, max=1)
  3. What is the output of str(all.dat.resids[, c('DLB_VH', 'DLB_NVH')]
  4. What is the output of all.dat.resids$resids.all[, levels(Group)]
giova1588 commented 5 years ago

1) version 2.7 2) str(all.dat.resids, max = 1) List of 6 $ X :List of 3 $ method : chr "sep.groups" $ use.mean : logi FALSE $ all.dat.tidy:Classes ‘data.table’ and 'data.frame': 3700 obs. of 7 variables: ..- attr(, ".internal.selfref")= ..- attr(, "index")= int(0) .. ..- attr(, "__Study.ID")= int [1:3700] 1 14 27 40 53 66 79 92 105 118 ... $ resids.all :Classes ‘data.table’ and 'data.frame': 37 obs. of 102 variables: ..- attr(, ".internal.selfref")= ..- attr(*, "sorted")= chr [1:2] "Group" "Study.ID" $ groups : chr [1:3] "DLB_NVH" "DLB_VH" "FENO"

3) str(all.dat.resids[, c('DLB_VH', 'DLB_NVH')]) List of 6 $ X :List of 3 ..$ DLB_NVH: num [1:13, 1:3] 1 1 1 1 1 1 1 1 1 1 ... .. ..- attr(, "dimnames")=List of 2 .. .. ..$ : NULL .. .. ..$ : chr [1:3] "Intercept" "Age" "SexM" ..$ DLB_VH : num [1:10, 1:3] 1 1 1 1 1 1 1 1 1 1 ... .. ..- attr(, "dimnames")=List of 2 .. .. ..$ : NULL .. .. ..$ : chr [1:3] "Intercept" "Age" "SexM" ..$ FENO : num [1:14, 1:3] 1 1 1 1 1 1 1 1 1 1 ... .. ..- attr(, "dimnames")=List of 2 .. .. ..$ : NULL .. .. ..$ : chr [1:3] "Intercept" "Age" "SexM" $ method : chr "sep.groups" $ use.mean : logi FALSE $ all.dat.tidy:Classes ‘data.table’ and 'data.frame': 2300 obs. of 7 variables: ..$ Study.ID: chr [1:2300] "subj001" "subj003" "subj009" "subj011" ... ..$ Group : Factor w/ 3 levels "DLB_NVH","DLB_VH",..: 1 1 1 1 1 1 1 1 1 1 ... ..$ Sex : Factor w/ 2 levels "F","M": 2 1 2 1 1 1 2 2 1 2 ... ..$ Age : int [1:2300] 73 50 63 77 69 71 69 76 77 68 ... ..$ region : Factor w/ 100 levels "lL_VisCent_ExStr_1",..: 1 1 1 1 1 1 1 1 1 1 ... ..$ value : num [1:2300] 1.1 1.6 1.53 1.89 1.94 ... ..$ resids : num [1:2300] -2.632 -0.577 0.227 1.141 1.472 ... ..- attr(, ".internal.selfref")= $ resids.all :Classes ‘data.table’ and 'data.frame': 23 obs. of 102 variables: ..$ Study.ID : chr [1:23] "subj001" "subj003" "subj009" "subj011" ... ..$ Group : Factor w/ 2 levels "DLB_NVH","DLB_VH": 1 1 1 1 1 1 1 1 1 1 ... ..$ lL_VisCent_ExStr_1 : num [1:23] -2.632 -0.577 0.227 1.141 1.472 ... ..$ lL_VisCent_ExStr_2 : num [1:23] 0.783 -0.345 0.244 0.114 2.152 ... ..$ lL_VisCent_ExStr_3 : num [1:23] -0.5684 0.0272 1.4868 2.3602 0.0583 ... ..$ lL_VisCent_ExStr_4 : num [1:23] -1.239 -1.009 0.777 0.564 2.181 ... ..$ lL_VisPeri_ExStrInf_1 : num [1:23] -1.0219 -0.0903 0.75 1.7024 0.9992 ... ..$ lL_VisPeri_ExStrInf_2 : num [1:23] 0.1378 0.7123 0.0309 2.5843 0.1806 ... ..$ lL_VisPeri_ExStrInf_3 : num [1:23] -1.037 -0.812 0.785 1.264 1.733 ... ..$ lL_SomMotA_1 : num [1:23] 0.63433 0.18363 0.00142 1.28885 1.16187 ... ..$ lL_SomMotA_2 : num [1:23] 0.292 0.115 0.123 2.621 0.367 ... ..$ lL_SomMotB_Aud_1 : num [1:23] 1.2726 0.0683 0.2522 0.8098 1.2448 ... ..$ lL_SomMotB_Aud_2 : num [1:23] -0.741 0.551 -0.068 0.598 0.685 ... ..$ lL_SomMotB_Aud_3 : num [1:23] -0.342 -0.742 0.593 0.664 2.684 ... ..$ lL_SomMotB_Aud_4 : num [1:23] 0.5315 -0.0423 -0.7707 0.1696 2.1371 ... ..$ lL_DorsAttnA_SPL_1 : num [1:23] -2.17 -1.827 1.687 0.765 0.87 ... .. [list output truncated] ..- attr(, ".internal.selfref")= ..- attr(, "sorted")= chr [1:2] "Group" "Study.ID" $ groups : chr [1:3] "DLB_NVH" "DLB_VH" "FENO"

4) all.dat.resids$resids.all[, levels(Group)] [1] "DLB_NVH" "DLB_VH" "FENO"

cwatson commented 5 years ago

I see. Please try running the permutation analysis, but this time, when you subset the residuals add g=; i.e.,

brainGraph_permute(densities, all.dat.resids[, g=c('DLB_VH', 'DLB_NVH')], ...)

giova1588 commented 5 years ago

ok, i think that the permutation works in this way, because the output are different for the three comparison. Now how i could visualize them? Because if i run print(plot_global(dt.G.tidy, vline = densities[15], level.names = level.names, exclude = NULL, perms = perms.all.graph_VHvsNVH, g = g)) it return "Error in obs[, i] : index out of limits" I thought that it is because the variable "dt.G.tidy" contains all three groups and so i exclude the group FENO with print(plot_global(dt.G.tidy[!Group == "FENO"], vline = densities[15], level.names = level.names, exclude = NULL, perms = perms.all.graph_VHvsNVH, g = g)), but there is still the same error. I have the same error with summary(perms.all.graph_VHvsNVH, measure = 'Lp', alt = 'two.sided')

cwatson commented 5 years ago

Hi, plot_global is a pretty old function but if you try subsetting g to contain only the groups in the permutation object, then it may work; e.g., g=g[1:2] or whatever the group indices would be. (In addition to excluding the relevant group from the data.table, as you did in the second command)

RE the summary error: you receive the same error?

giova1588 commented 5 years ago

With summary i have the same error and I don't know how I could visualize the summary of the permutation statistics...

cwatson commented 5 years ago

If you email to me the permutation results in the form of ".rda" files, along with any other variables needed to inspect the data, I will take a look.

cwatson commented 5 years ago

I will close this issue, as it was solved via email.

For posterity, this was due to a bug in groups vector of the permutation results object was incorrect. In this case, it included group names that were not part of the analysis. This bug will be fixed in an upcoming release.