How to use the nilearn.plotting.plot_stat_map function in nilearn

from nilearn.image import math_img
from nilearn.input_data import NiftiMasker
p_val = second_level_model.compute_contrast(second_level_contrast='age',
                                            output_type='p_value')
masker = NiftiMasker(mask_img=mask_img).fit(p_val)
n_voxel = np.size(masker.transform(p_val))
# Correcting the p-values for multiple testing and taking neg log
neg_log_pval = math_img("-np.log10(np.minimum(1, img * {}))"
                        .format(str(n_voxel)),
                        img=p_val)

###########################################################################
# Let us plot the second level contrast
from nilearn import plotting
cut_coords = [-4, 26]
display = plotting.plot_stat_map(
    neg_log_pval, colorbar=True, display_mode='z', cut_coords=cut_coords)
plotting.show()

##############################################################################
# Computing the (corrected) p-values with permutation test
from nistats.second_level_model import non_parametric_inference
neg_log_pvals_permuted_ols_unmasked = \
    non_parametric_inference(gray_matter_map_filenames,
                             design_matrix=design_matrix,
                             second_level_contrast='age',
                             model_intercept=True, n_perm=1000,
                             two_sided_test=False, mask=mask_img,
                             smoothing_fwhm=2.0, n_jobs=1)

###########################################################################
# Let us plot the second level contrast

# Two types of strategies can be used from this threshold function
# Type 1: strategy used will be based on scoreatpercentile
threshold_percentile_img = threshold_img(tmap_filename, threshold='97%', copy=False)


# Type 2: threshold strategy used will be based on image intensity
# Here, threshold value should be within the limits i.e. less than max value.
threshold_value_img = threshold_img(tmap_filename, threshold=3.0, copy=False)

################################################################################
# Visualization
# Showing thresholding results by importing plotting modules and its utilities
from nilearn import plotting

# Showing percentile threshold image
plotting.plot_stat_map(threshold_percentile_img, display_mode='z', cut_coords=5,
                       title='Threshold image with string percentile', colorbar=False)

# Showing intensity threshold image
plotting.plot_stat_map(threshold_value_img, display_mode='z', cut_coords=5,
                       title='Threshold image with intensity value', colorbar=False)

################################################################################
# Extracting the regions by importing connected regions function
from nilearn.regions import connected_regions

regions_percentile_img, index = connected_regions(threshold_percentile_img,
                                                  min_region_size=1500)

regions_value_img, index = connected_regions(threshold_value_img,
                                             min_region_size=1500)

# reverse feature selection
coef = feature_selection.inverse_transform(coef)
# reverse variance threshold
coef = variance_threshold.inverse_transform(coef)
# reverse masking
weight_img = nifti_masker.inverse_transform(coef)

# Create the figure
from nilearn.plotting import plot_stat_map, show
bg_filename = gray_matter_map_filenames[0]
z_slice = 0


fig = plt.figure(figsize=(5.5, 7.5), facecolor='k')
# Hard setting vmax to highlight weights more
display = plot_stat_map(weight_img, bg_img=bg_filename,
                        display_mode='z', cut_coords=[z_slice],
                        figure=fig, vmax=1)
display.title('SVM weights', y=1.2)

# Measure accuracy with cross validation
from sklearn.model_selection import cross_val_score
cv_scores = cross_val_score(anova_svr, gm_maps_masked, age)

# Return the corresponding mean prediction accuracy
prediction_accuracy = np.mean(cv_scores)
print("=== ANOVA ===")
print("Prediction accuracy: %f" % prediction_accuracy)
print("")

### Inference with massively univariate model #################################
print("Massively univariate model")

##########################################################################
# Estimating the contrast is very simple. We can just provide the column
# name of the design matrix.
z_map = second_level_model.compute_contrast(second_level_contrast=[1, 0, 0],
                                            output_type='z_score')

###########################################################################
# We threshold the second level contrast at uncorrected p < 0.001 and plot it.
from nistats.thresholding import map_threshold
from nilearn import plotting

_, threshold = map_threshold(
    z_map, alpha=.05, height_control='fdr')
print('The FDR=.05-corrected threshold is: %.3g' % threshold)

display = plotting.plot_stat_map(
    z_map, threshold=threshold, colorbar=True, display_mode='z',
    cut_coords=[-4, 26],
    title='age effect on grey matter density (FDR = .05)')
plotting.show()

###########################################################################
# We can also study the effect of sex by computing the contrast, thresholding it
# and plot the resulting map.

z_map = second_level_model.compute_contrast(second_level_contrast='sex',
                                            output_type='z_score')
_, threshold = map_threshold(
    z_map, alpha=.05, height_control='fdr')
plotting.plot_stat_map(
    z_map, threshold=threshold, colorbar=True,
    title='sex effect on grey matter density (FDR = .05)')

if not os.path.isdir(self.output_dir):
            os.makedirs(self.output_dir)

        if self.algorithm is 'lassopcr':
            coef = np.dot(self._pca.components_.T,self._lasso.coef_)
            coef_img = self.nifti_masker.inverse_transform(np.transpose(coef))
        elif self.algorithm is 'pcr':
            coef = np.dot(self._pca.components_.T,self._regress.coef_)
            coef_img = self.nifti_masker.inverse_transform(np.transpose(coef))
        else:
            coef_img = self.nifti_masker.inverse_transform(predicter.coef_)

        overlay_img = nib.load(os.path.join(get_resource_path(),'MNI152_T1_2mm_brain.nii.gz'))

        fig1 = plot_stat_map(coef_img, overlay_img, title=self.algorithm + " weights",
                            cut_coords=range(-40, 40, 10), display_mode='z')
        fig1.savefig(os.path.join(self.output_dir, self.algorithm + '_weightmap_axial.png'))

        if self.prediction_type == 'classification':
            if self.algorithm not in ['svm','ridgeClassifier','ridgeClassifierCV']:
                fig2 = probability_plot(self.stats_output)
                fig2.savefig(os.path.join(self.output_dir, self.algorithm + '_prob_plot.png'))
            else:
                fig2 = dist_from_hyperplane_plot(self.stats_output)
                fig2.savefig(os.path.join(self.output_dir, self.algorithm +
                            '_xVal_Distance_from_Hyperplane.png'))
                if self.algorithm is 'svm' and self.predicter.probability:
                    fig3 = probability_plot(self.stats_output)
                    fig3.savefig(os.path.join(self.output_dir, self.algorithm + '_prob_plot.png'))

        elif self.prediction_type == 'prediction':

motor_images = datasets.fetch_neurovault_motor_task()
motor_images.images

###############################################################################
# motor_images is a list of filenames. We need to take the first one
tmap_filename = motor_images.images[0]


###############################################################################
# Visualizing a 3D file
# ----------------------
#
# The file contains a 3D volume, we can easily visualize it as a
# statistical map:
from nilearn import plotting
plotting.plot_stat_map(tmap_filename)

###############################################################################
# Visualizing works better with a threshold
plotting.plot_stat_map(tmap_filename, threshold=3)


###############################################################################
# Visualizing one volume in a 4D file
# -----------------------------------
#
# We can download resting-state networks from the Smith 2009 study on
# correspondance between rest and task
rsn = datasets.fetch_atlas_smith_2009()['rsn10']
rsn

###############################################################################

def plot_double(img, img2, name, output_dir):
    import matplotlib
    matplotlib.use('agg')
    import matplotlib.pyplot as plt
    from nilearn.plotting import plot_stat_map, find_xyz_cut_coords
    vmax = np.abs(img.get_data()).max()
    cut_coords = find_xyz_cut_coords(img, activation_threshold=vmax / 3)
    fig, axes = plt.subplots(2, 1, figsize=(8, 8))
    plot_stat_map(img, title=name, figure=fig, axes=axes[0], threshold=0,
                  cut_coords=cut_coords)
    plot_stat_map(img2, figure=fig, axes=axes[1],
                  threshold=0,
                  cut_coords=cut_coords)
    plt.savefig(join(output_dir, '%s.png' % name))
    plt.close(fig)

# use a simple Bonferroni correction.
thresholded_map3, threshold3 = map_threshold(
    z_map, alpha=.05, height_control='bonferroni')
print('The p<.05 Bonferroni-corrected threshold is %.3g' % threshold3)

#########################################################################
# Visualize the results
# ---------------------
#
# First, the unthresholded map.
from nilearn import plotting
display = plotting.plot_stat_map(z_map, title='Raw z map')

#########################################################################
# Second, the p<.001 uncorrected-thresholded map (with only clusters > 10 voxels).
plotting.plot_stat_map(
    thresholded_map1, cut_coords=display.cut_coords, threshold=threshold1,
    title='Thresholded z map, fpr <.001, clusters > 10 voxels')

#########################################################################
# Third, the fdr-thresholded map.
plotting.plot_stat_map(thresholded_map2, cut_coords=display.cut_coords,
                       title='Thresholded z map, expected fdr = .05',
                       threshold=threshold2)

#########################################################################
# Fourth, the Bonferroni-thresholded map.
plotting.plot_stat_map(thresholded_map3, cut_coords=display.cut_coords,
                       title='Thresholded z map, expected fwer < .05',
                       threshold=threshold3)

#########################################################################

neg_log_pvals_anova = - np.log10(pvals_anova)
neg_log_pvals_anova_unmasked = nifti_masker.inverse_transform(
    neg_log_pvals_anova)

############################################################################
# Visualization
from nilearn.plotting import plot_stat_map, show

# Various plotting parameters
z_slice = 45  # plotted slice

threshold = - np.log10(0.1)  # 10% corrected

# Plot Anova p-values
fig = plt.figure(figsize=(5, 6), facecolor='w')
display = plot_stat_map(neg_log_pvals_anova_unmasked,
                        threshold=threshold,
                        display_mode='z', cut_coords=[z_slice],
                        figure=fig)

masked_pvals = np.ma.masked_less(neg_log_pvals_anova_unmasked.get_data(),
                                 threshold)

title = ('Negative $\log_{10}$ p-values'
         '\n(Parametric + Bonferroni correction)'
         '\n%d detections' % (~masked_pvals.mask).sum())

display.title(title, y=1.1, alpha=0.8)

show()

print(coef_img)

###########################################################################
# coef_img is now a NiftiImage.
#
# We can save the coefficients as a nii.gz file:
coef_img.to_filename('haxby_svc_weights.nii.gz')

###########################################################################
# Plotting the SVM weights
# .........................
#
# We can plot the weights, using the subject's anatomical as a background
from nilearn.plotting import plot_stat_map, show

plot_stat_map(coef_img, bg_img=haxby_dataset.anat[0],
              title="SVM weights", display_mode="yx")

show()