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