1"""Utilities for probabilistic error control at voxel- and 2cluster-level in brain imaging: cluster-level thresholding, false 3discovery rate control, false discovery proportion in clusters. 4 5Author: Bertrand Thirion, 2015 -- 2019 6 7""" 8import warnings 9 10import numpy as np 11from scipy.ndimage import label 12from scipy.stats import norm 13 14from nilearn.input_data import NiftiMasker 15from nilearn.image import get_data, math_img 16 17 18def _compute_hommel_value(z_vals, alpha, verbose=False): 19 """Compute the All-Resolution Inference hommel-value""" 20 if alpha < 0 or alpha > 1: 21 raise ValueError('alpha should be between 0 and 1') 22 z_vals_ = - np.sort(- z_vals) 23 p_vals = norm.sf(z_vals_) 24 n_samples = len(p_vals) 25 26 if len(p_vals) == 1: 27 return p_vals[0] > alpha 28 if p_vals[0] > alpha: 29 return n_samples 30 slopes = (alpha - p_vals[: - 1]) / np.arange(n_samples, 1, -1) 31 slope = np.max(slopes) 32 hommel_value = np.trunc(n_samples + (alpha - slope * n_samples) / slope) 33 if verbose: 34 try: 35 from matplotlib import pyplot as plt 36 except ImportError: 37 warnings.warn('"verbose" option requires the package Matplotlib.' 38 'Please install it using `pip install matplotlib`.') 39 else: 40 plt.figure() 41 plt.plot(p_vals, 'o') 42 plt.plot([n_samples - hommel_value, n_samples], [0, alpha]) 43 plt.plot([0, n_samples], [0, 0], 'k') 44 plt.show(block=False) 45 return np.minimum(hommel_value, n_samples) 46 47 48def _true_positive_fraction(z_vals, hommel_value, alpha): 49 """Given a bunch of z-avalues, return the true positive fraction 50 51 Parameters 52 ---------- 53 z_vals : array, 54 A set of z-variates from which the FDR is computed. 55 56 hommel_value: int 57 The Hommel value, used in the computations. 58 59 alpha : float 60 The desired FDR control. 61 62 Returns 63 ------- 64 threshold : float 65 Estimated true positive fraction in the set of values. 66 67 """ 68 z_vals_ = - np.sort(- z_vals) 69 p_vals = norm.sf(z_vals_) 70 n_samples = len(p_vals) 71 c = np.ceil((hommel_value * p_vals) / alpha) 72 unique_c, counts = np.unique(c, return_counts=True) 73 criterion = 1 - unique_c + np.cumsum(counts) 74 proportion_true_discoveries = np.maximum(0, criterion.max() / n_samples) 75 return proportion_true_discoveries 76 77 78def fdr_threshold(z_vals, alpha): 79 """Return the Benjamini-Hochberg FDR threshold for the input z_vals 80 81 Parameters 82 ---------- 83 z_vals : array 84 A set of z-variates from which the FDR is computed. 85 86 alpha : float 87 The desired FDR control. 88 89 Returns 90 ------- 91 threshold : float 92 FDR-controling threshold from the Benjamini-Hochberg procedure. 93 94 """ 95 if alpha < 0 or alpha > 1: 96 raise ValueError( 97 'alpha should be between 0 and 1. {} was provided'.format(alpha)) 98 z_vals_ = - np.sort(- z_vals) 99 p_vals = norm.sf(z_vals_) 100 n_samples = len(p_vals) 101 pos = p_vals < alpha * np.linspace( 102 .5 / n_samples, 1 - .5 / n_samples, n_samples) 103 if pos.any(): 104 return (z_vals_[pos][-1] - 1.e-12) 105 106 return np.infty 107 108 109def cluster_level_inference(stat_img, mask_img=None, 110 threshold=3., alpha=.05, verbose=False): 111 """ Report the proportion of active voxels for all clusters 112 defined by the input threshold. 113 114 This implements the method described in [1]_. 115 116 Parameters 117 ---------- 118 stat_img : Niimg-like object or None, optional 119 statistical image (presumably in z scale) 120 121 mask_img : Niimg-like object, optional, 122 mask image 123 124 threshold : list of floats, optional 125 Cluster-forming threshold in z-scale. Default=3.0. 126 127 alpha : float or list, optional 128 Level of control on the true positive rate, aka true dsicovery 129 proportion. Default=0.05. 130 131 verbose : bool, optional 132 Verbosity mode. Default=False. 133 134 Returns 135 ------- 136 proportion_true_discoveries_img : Nifti1Image 137 The statistical map that gives the true positive. 138 139 Notes 140 ----- 141 This function is experimental. 142 It may change in any future release of Nilearn. 143 144 References 145 ---------- 146 .. [1] Rosenblatt JD, Finos L, Weeda WD, Solari A, Goeman JJ. 147 All-Resolutions Inference for brain imaging. 148 Neuroimage. 2018 Nov 1;181:786-796. 149 doi: 10.1016/j.neuroimage.2018.07.060 150 151 """ 152 153 if not isinstance(threshold, list): 154 threshold = [threshold] 155 156 if mask_img is None: 157 masker = NiftiMasker(mask_strategy='background').fit(stat_img) 158 else: 159 masker = NiftiMasker(mask_img=mask_img).fit() 160 stats = np.ravel(masker.transform(stat_img)) 161 hommel_value = _compute_hommel_value(stats, alpha, verbose=verbose) 162 163 # embed it back to 3D grid 164 stat_map = get_data(masker.inverse_transform(stats)) 165 166 # Extract connected components above threshold 167 proportion_true_discoveries_img = math_img('0. * img', img=stat_img) 168 proportion_true_discoveries = masker.transform( 169 proportion_true_discoveries_img).ravel() 170 171 for threshold_ in sorted(threshold): 172 label_map, n_labels = label(stat_map > threshold_) 173 labels = label_map[get_data(masker.mask_img_) > 0] 174 for label_ in range(1, n_labels + 1): 175 # get the z-vals in the cluster 176 cluster_vals = stats[labels == label_] 177 proportion = _true_positive_fraction(cluster_vals, hommel_value, 178 alpha) 179 proportion_true_discoveries[labels == label_] = proportion 180 181 proportion_true_discoveries_img = masker.inverse_transform( 182 proportion_true_discoveries) 183 return proportion_true_discoveries_img 184 185 186def threshold_stats_img(stat_img=None, mask_img=None, alpha=.001, threshold=3., 187 height_control='fpr', cluster_threshold=0, 188 two_sided=True): 189 """ Compute the required threshold level and return the thresholded map 190 191 Parameters 192 ---------- 193 stat_img : Niimg-like object or None, optional 194 Statistical image (presumably in z scale) whenever height_control 195 is 'fpr' or None, stat_img=None is acceptable. 196 If it is 'fdr' or 'bonferroni', an error is raised if stat_img is None. 197 198 mask_img : Niimg-like object, optional, 199 Mask image 200 201 alpha : float or list, optional 202 Number controlling the thresholding (either a p-value or q-value). 203 Its actual meaning depends on the height_control parameter. 204 This function translates alpha to a z-scale threshold. 205 Default=0.001. 206 207 threshold : float, optional 208 Desired threshold in z-scale. 209 This is used only if height_control is None. Default=3.0. 210 211 height_control : string, or None optional 212 False positive control meaning of cluster forming 213 threshold: None|'fpr'|'fdr'|'bonferroni' 214 Default='fpr'. 215 216 cluster_threshold : float, optional 217 cluster size threshold. In the returned thresholded map, 218 sets of connected voxels (`clusters`) with size smaller 219 than this number will be removed. Default=0. 220 221 two_sided : Bool, optional 222 Whether the thresholding should yield both positive and negative 223 part of the maps. 224 In that case, alpha is corrected by a factor of 2. 225 Default=True. 226 227 Returns 228 ------- 229 thresholded_map : Nifti1Image, 230 The stat_map thresholded at the prescribed voxel- and cluster-level. 231 232 threshold : float 233 The voxel-level threshold used actually. 234 235 Notes 236 ----- 237 If the input image is not z-scaled (i.e. some z-transformed statistic) 238 the computed threshold is not rigorous and likely meaningless 239 240 This function is experimental. 241 It may change in any future release of Nilearn. 242 243 See also 244 -------- 245 nilearn.image.threshold_img 246 247 """ 248 height_control_methods = ['fpr', 'fdr', 'bonferroni', 249 'all-resolution-inference', None] 250 if height_control not in height_control_methods: 251 raise ValueError( 252 "height control should be one of {0}", height_control_methods) 253 254 # if two-sided, correct alpha by a factor of 2 255 alpha_ = alpha / 2 if two_sided else alpha 256 257 # if height_control is 'fpr' or None, we don't need to look at the data 258 # to compute the threshold 259 if height_control == 'fpr': 260 threshold = norm.isf(alpha_) 261 262 # In this case, and if stat_img is None, we return 263 if stat_img is None: 264 if height_control in ['fpr', None]: 265 return None, threshold 266 else: 267 raise ValueError('Map_threshold requires stat_img not to be None' 268 'when the height_control procedure ' 269 'is "bonferroni" or "fdr"') 270 271 if mask_img is None: 272 masker = NiftiMasker(mask_strategy='background').fit(stat_img) 273 else: 274 masker = NiftiMasker(mask_img=mask_img).fit() 275 stats = np.ravel(masker.transform(stat_img)) 276 n_voxels = np.size(stats) 277 278 # Thresholding 279 if two_sided: 280 # replace stats by their absolute value after storing the sign 281 sign = np.sign(stats) 282 stats = np.abs(stats) 283 284 if height_control == 'fdr': 285 threshold = fdr_threshold(stats, alpha_) 286 elif height_control == 'bonferroni': 287 threshold = norm.isf(alpha_ / n_voxels) 288 stats *= (stats > threshold) 289 if two_sided: 290 stats *= sign 291 292 # embed it back to 3D grid 293 stat_map = get_data(masker.inverse_transform(stats)) 294 295 # Extract connected components above threshold 296 label_map, n_labels = label(np.abs(stat_map) > threshold) 297 labels = label_map[get_data(masker.mask_img_) > 0] 298 299 for label_ in range(1, n_labels + 1): 300 if np.sum(labels == label_) < cluster_threshold: 301 stats[labels == label_] = 0 302 303 return masker.inverse_transform(stats), threshold 304