'Residual mean squares' is property of model not the contrast

[cmaumet]

Discussed with @nicholst, @TomMaullin and @AlexBowring, the computation of residual mean square at the second-level needs to be fixed.

https://github.com/incf-nidash/nidmresults-fsl/blob/master/nidmfsl/fsl_exporter/fsl_exporter.py#L806-L809

RMS is property of the model not of a contrast. Actually the contrast number is hard-coded to '1'.

[cmaumet]

Waiting on @nicholst's feedback on the best way to compute this.

[nicholst]

Background:

1. FEAT first level fMRI produces a stats/sigmasquared.nii.gz; these are the whitened residual MSE (whitening happens around film_gls.cc:229, and writing out around film_gls.cc:229 via the glimGls object defined in glimGls.cc). They are an attribute of the entire model (at each voxel).
2. FEAT second level fMRI does not produce any analgous file. It produces a copeX.feat/stats/mean_random_effects_var1.nii.gz; it appears (but I hope @AlexBowring can confirm) that this file is produced whether using "Mixed effects: FLAME1" or "Mixed effects: Simple OLS"; I don't know about "Fixed Effects". However this file (files, for multiple contrasts) are contrast-specific, and not a property of the model over all.
3. Variance groups. Confusingly, the "1" in mean_random_effects_var1.nii.gz does not refer to the contrast number but rather to variance group.
4. The variance at the 2nd level varies with scan; the for each scan is the sum of the copeX.feat/var_filtered_func_data and the mean_random_effects_var?.nii.gz where ? is the variance group to which that scan belongs.

If we must create a single residual mean squared variance, it seems that the most direct thing is to take the average of the scan-specific variance. For a single variance group, this is just the sum of mean_random_effects_var1.nii.gz and the ("temporal") mean of the 4D copeX.feat/var_filtered_func_data. For arbitrary K groups. the K mean_random_effects_var?.nii.gz images would need to be combined with a weighted mean, where the weights are determined by the number of scans in that variance group.

So! The question: At this point, do we want to support more than one variance group?

[nicholst]

Notes from further reflection:

The appropriate sigma squared should be easy to calculate for the multi-group case

Where \sigma^2_i is the i-th varcope from copeX.feat/var_filtered_func_data, \sigma^2_g(i) is the between subject variance for the i-th subject (who is in group g(i)) found in copeX.feat/stats/mean_random_effects_var?.nii.gz (?=g), and N_g is the number of subjects in the g-th group.

It remains to be seen if we want to support multiple groups right now.