1"""Example of explicit fixed effects fMRI model fitting
2====================================================
3
4This example illustrates how to run a fixed effects model based on
5pre-computed statistics. This is helpful when the initial models
6have to be fit separately.
7
8For details on the data, please see:
9
10Dehaene-Lambertz G, Dehaene S, Anton JL, Campagne A, Ciuciu P, Dehaene
11G, Denghien I, Jobert A, LeBihan D, Sigman M, Pallier C, Poline
12JB. Functional segregation of cortical language areas by sentence
13repetition. Hum Brain Mapp. 2006: 27:360--371.
14http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=2653076#R11
15
16Please see `Simple example of two-session :term:`fMRI` model fitting
17<https://nistats.github.io/auto_examples/02_first_levels/plot_fiac_analysis.html>`_
18example for details.  The main difference is that
19the fixed-effects model is run explicitly here,
20after GLM fitting on two sessions.
21
22"""
23
24#########################################################################
25# Prepare data and analysis parameters
26# --------------------------------------
27#
28# Inspecting 'data', we note that there are two sessions
29
30from nilearn.datasets import func
31data = func.fetch_fiac_first_level()
32fmri_img = [data['func1'], data['func2']]
33
34#########################################################################
35# Create a mean image for plotting purpose
36from nilearn.image import mean_img
37mean_img_ = mean_img(fmri_img[0])
38
39#########################################################################
40# The design matrices were pre-computed, we simply put them in a list of
41# DataFrames
42design_files = [data['design_matrix1'], data['design_matrix2']]
43import pandas as pd
44import numpy as np
45design_matrices = [pd.DataFrame(np.load(df)['X']) for df in design_files]
46
47#########################################################################
48# GLM estimation
49# ----------------------------------
50# GLM specification. Note that the mask was provided in the dataset.
51# So we use it.
52
53from nilearn.glm.first_level import FirstLevelModel
54fmri_glm = FirstLevelModel(mask_img=data['mask'], smoothing_fwhm=5,
55                           minimize_memory=True)
56
57#########################################################################
58# Compute fixed effects of the two runs and compute related images
59# For this, we first define the contrasts as we would do for a single session
60n_columns = design_matrices[0].shape[1]
61contrast_val = np.hstack(([-1, -1, 1, 1], np.zeros(n_columns - 4)))
62
63#########################################################################
64# Statistics for the first session
65from nilearn import plotting
66cut_coords = [-129, -126, 49]
67contrast_id = 'DSt_minus_SSt'
68
69fmri_glm = fmri_glm.fit(fmri_img[0], design_matrices=design_matrices[0])
70summary_statistics_session1 = fmri_glm.compute_contrast(
71    contrast_val, output_type='all')
72plotting.plot_stat_map(
73    summary_statistics_session1['z_score'],
74    bg_img=mean_img_, threshold=3.0, cut_coords=cut_coords,
75    title='{0}, first session'.format(contrast_id))
76
77#########################################################################
78# Statistics for the second session
79
80fmri_glm = fmri_glm.fit(fmri_img[1], design_matrices=design_matrices[1])
81summary_statistics_session2 = fmri_glm.compute_contrast(
82    contrast_val, output_type='all')
83plotting.plot_stat_map(
84    summary_statistics_session2['z_score'],
85    bg_img=mean_img_, threshold=3.0, cut_coords=cut_coords,
86    title='{0}, second session'.format(contrast_id))
87
88#########################################################################
89# Fixed effects statistics
90from nilearn.glm.contrasts import compute_fixed_effects
91
92contrast_imgs = [summary_statistics_session1['effect_size'],
93                 summary_statistics_session2['effect_size']]
94variance_imgs = [summary_statistics_session1['effect_variance'],
95                 summary_statistics_session2['effect_variance']]
96
97fixed_fx_contrast, fixed_fx_variance, fixed_fx_stat = compute_fixed_effects(
98    contrast_imgs, variance_imgs, data['mask'])
99plotting.plot_stat_map(
100    fixed_fx_stat, bg_img=mean_img_, threshold=3.0, cut_coords=cut_coords,
101    title='{0}, fixed effects'.format(contrast_id))
102
103#########################################################################
104# Not unexpectedly, the fixed effects version displays higher peaks than the
105# input sessions. Computing fixed effects enhances the signal-to-noise ratio of
106# the resulting brain maps
107# Note however that, technically, the output maps of the fixed effects map is a
108# t statistic (not a z statistic)
109
110plotting.show()
111