The display of results for voxelwise WB does not seem to use WB p-values for uncorrected thresholds

[BryanGuillaume]

For a voxelwise WB, the display of results seems to threshold the non-parametric p-value images when an FDR or an FWER correction is used (see xCon(Ic).VspmFDRP and xCon(Ic).VspmFWEP in the code below). That is what I would expect as behaviour. However, when the treshold is an uncorrected p-value, it seems that the threshold is done on the equivalent parametric Z-score (variable Z in the code below)and not on the non-parametric uncorrected p-values. I would have expected that would be done on the non-parametric uncorrected p-values. Either this current behaviour is expected and we should document this better or we refactor the code to threshold the non-parametric uncorrected p-value instead. We could also offer the possibility to threshold both in the GUI.

          switch thresDesc

              case 'FWE' % Family-wise false positive rate
                  % This is performed on the voxelwise FWE P value map
                  %--------------------------------------------------------
                  try
                      pu = xSwE.u;
                  catch
                      pu = spm_input('p value (FWE)','+0','r',0.05,1,[0,1]);
                  end
                  thresDesc = ['p<' num2str(pu) ' (' thresDesc ')'];

                  FWE_ps = 10.^-spm_get_data(xCon(Ic).VspmFWEP,XYZum);

                  Q      = find(FWE_ps  < pu);

                  % Obtain statistic threshold (this will be the maximum
                  % statistic value that did not pass thresholding when Q
                  % was applied to the statistic).
                  if ~isempty(Q)
                      u = max(Z(Z<min(Z(Q))));
                  else
                      u = Inf;
                  end

              case 'FDR' % False discovery rate
                  % This is performed on the FDR P value map
                  %--------------------------------------------------------
                  try
                      pu = xSwE.u;
                  catch
                      pu = spm_input('p value (FDR)','+0','r',0.05,1,[0,1]);
                  end
                  thresDesc = ['p<' num2str(pu) ' (' thresDesc ')'];

                  FDR_ps = 10.^-spm_get_data(xCon(Ic).VspmFDRP,XYZum);

                  Q      = find(FDR_ps  < pu);

                  % Obtain statistic threshold
                  switch STAT
                      case 'T'
                         u = spm_uc_FDR(pu,Inf,'Z',n,VspmSv,0); 
                      case 'F'
                         u = spm_uc_FDR(pu,[1 1],'X',n,VspmSv,0); 
                  end

              case 'none'  % No adjustment: p for conjunctions is p of the conjunction SwE
                  % This is performed in the normal manor on the Z map.
                  %--------------------------------------------------------
                  try
                      u = xSwE.u;
                  catch
                      u = spm_input(['threshold {',eSTAT,' or p value}'],'+0','r',0.001,1);
                  end
                  if u <= 1
                      thresDesc = ['p<' num2str(u) ' (unc.)'];
                      switch STAT
                          case 'T'
                              u  = swe_invNcdf(1-u^(1/n));
                          case 'F'
                              u  = spm_invXcdf(1-u^(1/n),1);
                      end
                  else
                      thresDesc = '';
                  end

                  Q      = find(Z > u);

              otherwise
                  %--------------------------------------------------------------
                  fprintf('\n');                                              %-#
                  error('Unknown control method "%s".',thresDesc);
          end

[BryanGuillaume]

The same type of behaviour seems to exist for clusterwise inference for WB analyses which does not use an FWER correction (option none). Indeed, the clusters are formed based on the parametric Z-scores. I am not sure that this behaviour is obvious for the users because we are actually within a non-parametric WB analysis. I believe there is a need to clarify what the toolbox is doing or simply offer the option to the user to threshold either the parametric or the non-parametric uncorrected p-values.

              case 'none'  % No adjustment: p for conjunctions is p of the conjunction SwE
                  % This is performed in the normal manor on the Z map.
                  %--------------------------------------------------------
                  % Record what type of clusterwise inference we are doing.
                  clustWise = 'Uncorr';

                  % Cluster-forming threshold.
                  try
                      u = xSwE.u;
                  catch
                      u = spm_input(['threshold {',eSTAT,' or p value}'],'+1','r',0.001,1);
                  end
                  if u <= 1
                      thresDesc = ['p<' num2str(u) ' (unc.)'];
                      switch STAT
                          case 'T'
                              u  = swe_invNcdf(1-u^(1/n));
                          case 'F'
                              u  = spm_invXcdf(1-u^(1/n),1);
                      end
                  else
                      thresDesc = '';
                  end

                  up  = NaN;
                  Pp  = NaN;
                  uc  = NaN;
                  ue  = NaN;
                  Pc  = [];
                  uu = [];

                  Q      = find(Z > u);

[nicholst]

I agree we need to improve the clarity of what's done. As discussed on https://github.com/NISOx-BDI/SwE-toolbox/issues/124#issuecomment-502129965 for FDR I suggested a user-definable action.

But we don't need that level of flexibility over all. Let's review the possible actions for WB and what I would have thought would be the natural derivation:

Inference type / threshold	Image / method
Voxelwise data visualised in MIP results display	Z or X (??)
Voxelwise uncorrected statistic threshold	Z or X
Voxelwise uncorrected p-value threshold	Nonparam P
Voxelwise FDR corrected p-value threshold	BH applied to Nonparam P
Voxelwise FWE corrected p-value threshold	Max dist from Z or X
Cluster FWE corrected p-value threshold	Clusters formed on Z or X
TFCE FWE corrected p-value threshold	TFCE formed on Z or X

... before I write more... does this make sense?

[nicholst]

Also, going back to your original comment,

For a voxelwise WB, the display of results seems to threshold the non-parametric p-value images when an FDR or an FWER correction is used (see xCon(Ic).VspmFDRP and xCon(Ic).VspmFWEP

I'm very surprised that we use non-parametric uncorrected p-value images to determine the voxelwise FWE correction. The non-parametric uncorrected p-values are highly discrete and as a result we usually don't use them with FWE. Can I confirm this is the case?

[BryanGuillaume]

I will review the code and make a summary (e.g., in a table like you did above) of what the toolbox is actually doing now for each option. That will be easier to understand what is globally going on and decide what to do.

[BryanGuillaume]

This issue has been handled by PR #131. Now, when a WB is used to estimate the model, only WB p-values are used.