Hi!
I would like use fixelcfestats for a paired t-test. How can I specify the exchangability-block information for running the permutations? -Like what you do with the .grp file in FSL’s randomise where you “specify which data can be exchanged, that is, only condition A and B within the same subject (since the null hypothesis is that there is no difference between A and B, but each subject could have a different mean value and so we cannot exchange between subjects).”
Cheers,
Örjan
Hi Örjan,
Unfortunately it’s not currently possible to specify the exchangability block information with fixelcfestats or mrclusterstats. However, this is something we will need so I’ll implement it when I have time after ISMRM and get back to you.
Cheers,
Dave
Thank you for the quick response!
Is there a way then to calculate difference images and do a one-sample t-test? Maybe a workaround for now could be to mod the code a little somewhere around line 353 to do this subtraction just before calculating the beta coefficients?
Cheers,
Örjan
…I of course meant something that I perhaps could try myself in the meantime just to be able to get the results. So close… Any hint would be much appreciated.
Örjan
Hi Örjan,
Apologies for the slow reply. I’ve just committed a new command called fixelcalc. You can use this to generate your difference images before running fixelcfestats in the short term. Note that the command is very rudimentary, and does not have all the features of mrcalc. I have been holding back on implementing a fully featured fixelcalc, since we are still deciding on if it’s worth changing the fixel file format (in which case we may be able to just use mrcalc on fixel images).
Cheers,
Dave
Thank you for putting time on this! That function will be very handy. However, (I should probably have said) the paired fixel images are from MZ twins and based on FOD images warped to a group template, which means that the inter-subject fixel correspondence is going to be off slightly. It will not be possible to subtract one image from the other without taking some max angle threshold into account. Or should I use fixelcorrespondence to map all fixel images to the group template first? Sounds right.
Cheers,
Örjan
Hi Orjan,
Ah yes, good point. This would be an issue even for multiple images of the same subject.
As you are probably aware, fixel correspondence (i.e. matching the fixels within a single voxel between subjects and the template) is currently performed within fixelcfestats. However, if you look at the registration branch in Github (and the associated fixel-based analysis user documentation) you will see the fixel correspondence step is now performed as a separate pre-processing step using the creatively named command fixelcorrespondence. You can use this to assign all your subject fixels to your template fixel image/mask, prior to computing the difference image.
Performing fixelcorrespondence as a separate preprocessing step makes the whole analyse pipeline more general, so it can be used to investigate fixel-based measures derived from any DWI model by just changing two steps. We have a paper under review that describes this in more detail.
I’ve got a couple more tweaks to make, then I will merge the registration branch in the next couple of days.
Cheers,
Dave
Fantastic! I now tried the registration branch (and the corresponding workflow) and it works just fine with my previous FOD group template and warps. It seems that I would need numpy to run the population_template function. I might try installing it later (if it is fine to install it separately?).
Thanks again,
Örjan
Great. Yes you will need to install numpy separately.
…Just a little something for the wish list (if tweaks are going to be made 
) I would reeeally like if it was possible to output (and enter) the precomputed fixel-fixel connectivity matrix so one could skip that step (which is rather time consuming with 10M streamlines) when running fixelcfestats again - if I e.g. would want to add a covariate to the design matrix.
Cheers,
Örjan
Yes, I’ve thought the same. I guess the issue is that depending on your image resolution, number of tracks, and tractography dispersion, the fixel-fixel connectivity matrix can get quite large (>100GB).
Another way around it would be to speed this step up with multi-threading. Its tricky because the (sparse) matrix is stored as a std::map that is a shared resource between threads. I tried a while back using C++11 atomic types, but did not have much success. Though I admit I did not spend much time on it.
I’ll look into this when I get a sec.
Hi!
Referring to the original question and solution of making difference images - Can fixelcfestats actually handle the one-sample t-test? Will it recognize the situation (single column of ones in the design matrix) and do sign-flips instead of permuting the rows?
Cheers,
Örjan
Hi Örjan,
Unfortunately not. As you are finding out very quickly, the GLM code in MRtrix is still fairly basic. I’m happy to add this in, however at the moment I need to prioritise an overdue manuscript revision (also on fixel-based analysis).
In terms of the implementation, how would you envision this working? If the contrast contains a single 1 (e.g. [1 0 0 0], which corresponds to column of 1’s in the design matrix, then a one-sample t-test is ‘detected’ and performed?
I see FSL randomise requires that the user declares a one-sampled test with -1. I’d be happy with an automatic method providing it’s well documented in the help.
It should not be too hard to code up and test when I get time. However if you are capable and willing, we are always happy to have new contributors to MRtrix . I’ve just had a quick look and have a reasonable idea of how this should be added. I’d be happy to give you some pointers if you are interested? Otherwise, it’s likely to be a couple of weeks before I have the time to add this and the exchangeability block option - sorry.
Cheers,
Dave
Hi Dave,
Thank you for the quick reply. I would have loved to be able to contribute but my programming skills are unfortunately at the level of basic survival. I might be able to create some idiosyncratic hack… As for the implementation I was just secretly hoping for an automatic method since I couldn’t find the explicit option. I understand that there is no shortage of requests and I am obviously very grateful for any help I can get whenever it is possible. At least the feature should be relevant to others as well. Good luck with the revision and I look forward to reading that paper!
Cheers,
Örjan
Hi Dave,
Was the paired testing ever implemented?
Cheers,
Lee
Hi Lee,
Not yet sorry. I’ve been a little sidetracked by some changes to the fixel-based analysis pipeline required to make it more robust/automatic for the Stanford neuroimaging data processing platform. The good thing is these changes will hopefully make it easier for everyone to run fixel-based analysis on their own clusters too (including easy ‘installation’ and updates using Linux containers).
It’s high on my TODO list since we will need this for several in-house projects too. I will post back here when it’s available.
Cheers,
Dave
Hi Lee,
I implemented it, but the solution is a bit customized to our own data etc., since I just wanted to try it out. Unfortunately I had to work on other things once I got this up and running (and it’s not thoroughly tested). However since you asked, I will try to find some time to have another look at it.
Cheers,
Örjan
Hi Lee,
If you want to have a look at the work in progress - you can copy this branch. The fixelcfestatsOneSample.exe functions similar to the original. There you have my first attempts at programming in c++ and a non-optimal but from what I can tell a functional implementation of the one-sample version of fixelcfestats. The idea is to do sign-flipping instead of random permutations. The total number of unique permutations is calculated and run automatically based on the number of files you input (if total unique nperms < 5000 which is the default). If you want to run less than the total unique nperms, just use the -nperms flag as usual. However, since the sign-flips are (for now) generated in a specific way, running less than the total amount might have an influence on the outcome. If you look below under point 2.a., an example would be running only the first three rows. When it comes to the implementation I’ve simply copied and renamed three relevant files (cmd/fixelcfestats.cpp -> fixelcfestatsOneSample; lib/math/stats/glm.h -> glmOneSample.h; stats/permtest.h -> permtestOneSample.h). The idea is to later merge changes in these files with the original and use some onesamplettest-flag. Anyway, these are the current changes (in principle):
In fixelcfestatsOneSample I create a “OneSampleGLMTTest” object, which is a modified copy of the GLMTTest object in the original glm.h.
In permtestOneSample.h:
a. There is a class function PermsMatrix::fillPermsMatrix that creates a matrix with all unique
permutations/sign-flips (rows) of {1 -1} for a dataset of size n (columns). E.g.
for 3 subjects:
1 1 1
1 1 -1
1 -1 1
1 -1 -1
-1 1 1
-1 1 -1
-1 -1 1
-1 -1 -1
b. In precompute_default_permutation, the first row (permutation) is submitted as an argument to the stats_calculator instead of default_labelling.
c. In the execute() of PreProcessor and Processor, the “permsMatrix(index)” i.e the current permutation (i.e. sign-flips) is passed as an argument to process_permutation (and subsequently to stats_calculator).
In glmOneSample.h:
a. There is a modified copy of the GLMTTest class (OneSamleGLMTTest).
b. As indicated above, the operator() is given the current sign-flip-permutation as an argument instead of the perm_labelling.
c. The default design matrix is maintained.
d. Each row in tmp is multiplied (column-wise) with the current permutation.
Hi Orjan,
Thanks for sharing. I will definitely incorporate/add these features in the future when I get a spare few days. It’s not something I want to code up in a rush in between other tasks given the possible downstream implications of a mistake.
Cheers,
Dave
Thanks so much Orjan,
I’ll be finishing processing for this in the next couple of days and give it a crack next week.
I’ll let you know how I go