We recommend that you install a complete scientific Python
- distribution like 64 bit Anaconda
- . Since it meets all the requirements of nistats, it will save
- you time and trouble. You could also check PythonXY
- as an alternative.
+ distribution like 64 bit
+ Anaconda . Since it meets all the requirements of nistats,
+ it will save you time and trouble. You could also check
+ PythonXY as an alternative.
Nistats requires a Python installation and the following
dependencies: ipython, scikit-learn, matplotlib and nibabel.
diff --git a/doc/introduction.rst b/doc/introduction.rst
index 6a921b56..47dcbe40 100644
--- a/doc/introduction.rst
+++ b/doc/introduction.rst
@@ -1,5 +1,5 @@
=====================================
-Introduction: nistats in a nutshell
+Introduction: Nistats in a nutshell
=====================================
.. contents:: **Contents**
@@ -7,19 +7,129 @@ Introduction: nistats in a nutshell
:depth: 1
-What is nistats: Bold-fMRI analysis
-===========================================================================
+What is Nistats?
+================
-.. topic:: **Why use nistats?**
+.. topic:: **What is Nistats?**
- Nistats is a python module to perform voxel-wise univariate analyses of functional magnetic resonance images (fMRI), using general linear models. It provides functions to create design matrices, at the subject and group levels, to estimate them from images series and to compute statistical maps (contrasts)
+ Nistats is a Python module to perform voxel-wise analyses of functional magnetic resonance images (fMRI) using linear models. It provides functions to create design matrices, at the subject and group levels, to estimate them from images series and to compute statistical maps (contrasts). It allows to perform the same statistical analyses as `SPM`_ or `FSL`_ (but it does not provide tools for preprocessing stages (realignment, spatial normalization, etc.); for this, see `NiPype`_.
- For tutorials, please check out the `Examples `_ gallery
+.. _SPM: https://www.fil.ion.ucl.ac.uk/spm/
+
+.. _FSL: https://www.fmrib.ox.ac.uk/fsl
+
+.. _NiPype: https://nipype.readthedocs.io/en/latest/
+
+
+
+A primer on BOLD-fMRI data analysis
+===================================
+
+What is fMRI ?
+--------------
+
+Functional magnetic resonance imaging (fMRI) is based on the fact that when local neural activity increases, increases in metabolism and blood flow lead to fluctuations of the relative concentrations of oxyhaemoglobin (the red cells in the blood that carry oxygen) and deoxyhaemoglobin (the same red cells after they have delivered the oxygen). Oxy-haemoglobin and deoxy-haemoglobin have different magnetic properties (diamagnetic and paramagnetic, respectively), and they affect the local magnetic field in different ways. The signal picked up by the MRI scanner is sensitive to these modifications of the local magnetic field. To record cerebral activity, during functional sessions, the scanner is tuned to detect this "Blood Oxygen Level Dependent" (BOLD) signal.
+
+Brain activity is measured in sessions that span several minutes, while the participant performs some a cognitive task and the scanner acquires brain images, typically every 2 or 3 seconds (the time between two successive image acquisition is called the Repetition time, or TR).
+
+A cerebral MR image provides a 3D image of the brain that can be decomposed into `voxels`_ (the equivalent of pixels, but in 3 dimensions). The series of images acquired during a functional session provides, in each voxel, a time series of positive real number representing the MRI signal, sampled at the TR.
+
+.. _voxels: https://en.wikipedia.org/wiki/Voxel
+
+.. note:: Before fMRI images can be used to do meaningful comparisons, they must be processed to ensure that the voxels that are being compared represent the same brain regions, irrespective of the variability in size and shape of the brain and its microarchitecture across different subjects in the experiment. The process is called spatial registration or spatial normalization. During this procedure, the voxels of all the brain images are 'registered' to correspond to the same region of the brain. Usually, the images (their voxels) are registered to a standard 'template' brain image (its voxels). One often used standard template is the MNI152 template from the Montreal Neurological Institute. Once this is done, the coordinates of a voxel are in the same space as the template and can be used to estimate its brain location using brain atlases based on that same template. As already mentioned, the Nistats package is not meant to perform spatial preprocessing, but only statistical analyses on the voxel times series, regardless of the coordinate system.
+
+fMRI data modelling
+-------------------
+
+One way to analyze times series consists in comparing them to a *model* built from our knowledge of the events that occurred during the functional session. Events can correspond to actions of the participant (e.g. button presses), presentations of sensory stimui (e.g. sound, images), or hypothesized internal processes (e.g. memorization of a stimulus), ...
+
+
+.. figure:: images/stimulation-time-diagram.png
+
+
+One expects that a brain region involved in the processing of a certain type of event (e.g. the auditory cortex for sounds), would show a time course of activation that correlates with the time-diagram of these events. If the fMRI signal directly showed neural activity and did not contain any noise, we could just look at it in various voxels and detect those that conform to the time-diagrams.
+
+Yet, we know, from previous measurements, that the BOLD signal does not follow the exact time course of stimulus processing and the underlying neural activity. The BOLD response reflects changes in blood flow and concentrations in oxy-deoxy haemoglobin, all together forming a `haemodynamic response`_ which is sluggish and long-lasting, as can be seen on the following figure showing the response to an impulsive event (for example, an auditory click played to the participants).
+
+.. figure:: images/spm_iHRF.png
+
+From the knowledge of the impulse haemodynamic response, we can build a predicted time course from the time-diagram of a type of event (The operation is known as `convolution`_. Simply stated, how the shape of one function's plot would affect the shape of another function's plot. **Remark:** it assumes linearity of the BOLD response, an assumption that may be wrong in some scenarios). It is this predicted time course, also known as a *predictor*, that is compared to the actual fMRI signal. If the correlation between the predictor and the signal is higher than expected by chance, the voxel is said to exhibit a significant response to the event type.
+
+
+.. _haemodynamic response: https://en.wikipedia.org/wiki/Haemodynamic_response
+.. _convolution: https://en.wikipedia.org/wiki/Convolution
+
+
+.. figure:: images/time-course-and-model-fit-in-a-voxel.png
+
+Correlations are computed separately at each voxel and a correlation map can be produced displaying the values of correlations (real numbers between -1 and +1) at each voxel. Generally, however, the maps presented in the papers report the significance of the correlations at each voxel, in forms of T, Z or p values for the null hypothesis test of no correlation (see below). For example, the following figure displays a Z-map showing voxels responding to auditory events. Large (positive or negative) Z values are unlikely to be due to chance alone. The map is tresholded so that only voxels with a p-value less than 1/1000 are coloured.
+
+.. note::
+ In this approach, hypothesis tests are conducted in parallel at many voxels, increasing the liklelihood of False Positives. This is known as the Problem of `Multiple Comparisons`_. Some common strategies for dealing with this are discussed later in this page. This issue can also be addressed in Nistats by using random permutations tests.
+
+.. figure:: images/example-spmZ_map.png
+
+
+In most fMRI experiments, several predictors are needed to fully describe the events occuring during the session -- for example, the experimenter may want to distinguish brain activities linked to the perception of auditory stimuli or to button presses. To find the effect specific to each predictor, a multiple `linear regression`_ approach is typically used: all predictors are entered as columns in a *design-matrix* and the software finds the linear combination of these columns that best fits the signal. The weights assigned to each predictor by this linear combination are estimates of the contribution of this predictor to the response in the voxel. One can plot this using effect size maps or, maps showing their statistical significance (how unlikely they are under the null hypothesis of no effect).
+
+
+.. _linear regression: https://en.wikipedia.org/wiki/Linear_regression
+
+In brief, the analysis of fMRI images involves:
+
+1. describing the paradigm in terms of events of various types occuring at certain times and having some durations.
+2. from this description, creating predictors for each type of event, typically using a convolution by the impulse haemodynamic response.
+3. assembling these predictors in a design-matrix, providing a *linear model*
+4. estimate the parameters of the model, that is, the weights associated with each predictors at each voxel, using linear regression.
+5. display the coefficients, or linear combination of them, and/or their statistical significance.
+
+fMRI statistical analysis
+-------------------------
+
+As explained in the previous section, the basic statistical analysis of fMRI is conceptually a correlation analysis, where one seeks whether a certain combination (contrast) of columns of the design matrix fits a significant proportion of the fMRI signal at a given location.
+
+It can be shown that this is equivalent to studying whether the estimated contrast effect is large with respect to the uncertainty about its exact value. Concretely, we compute the effect size estimate and the uncertainty about its value and divide the two. The resulting number has no physical dimension, it is a statistic --- a Student or t-statistic, which we will denote `t`.
+Next, based on `t`, we want to decide whether the true effect was indeed greater than zero or not.
+
+If the true effect were zero, `t` would not necessarily be 0: by chance, the noise in the data may be partly explained by the contrast of interest.
+However, if we assume that the noise is Gaussian, and that the model is correctly specified, then we know that `t` should follow a Student distribution with `dof` degrees of freedom, where q is the number of free parameters in the model: in practice, the number of observations (i.e. the number of time points), `n_scans` minus the number of effects modelled (i.e. the number of columns `n_columns`) of the design matrix:
+
+:math: `dof = n_scans - n_columns`
+
+With this we can do statistical inference: Given a pre-defined error rate :math:`\alpha`, we compare the observed `t` to the :math:`(1-\alpha)` quantile of the Student distribution with `dof` degrees of freedom. If t is greater than this number, we can reject the null hypothesis with a *p-value* :math:`\alpha`, meaning: if there were no effect, the probability of oberving an effect as large as `t` would be less than `\alpha`.
+
+.. figure:: images/student.png
+
+.. note:: A frequent misconception consists in interpreting :math:`1-\alpha` as the probability that there is indeed an effect: this is not true ! Here we rely on a frequentist approach, that does not support Bayesian interpretation. See e.g. https://en.wikipedia.org/wiki/Frequentist_inference
+
+.. note:: It is cumbersome to work with Student distributions, since those always require to specify the number `dof` of degrees of freedom. To avoid this, we can transform `t` to another variable `z` such that comparing `t` to the Student distribution with `dof` degrees of freedom is equivalent to comparing `z` to a standard normal distribution. We call this a z-transform of `t`. We call the :math:`(1-\alpha)` quantile of the normal distribution the *threshold*, since we use this value to declare voxels active or not.
+
+Multiple Comparisons
+--------------------
+
+A well-known issue that arrives then is that of multiple comparisons:
+ when a statistical tests is repeated a large number times, say one for each voxel, i.e. `n_voxels` times, then one can expect that, in the absence of any effect, the number of detections ---false detections since there is no effect--- will be roughly :math:`n\_voxels \alpha`. Then, take :math:`\alpha=.001` and :math:`n=10^5`, the number of false detections will be about 100. The danger is that one may no longer trust the detections, i.e. values of `z` larger than the :math:`(1-\alpha)`-quantile of the standard normal distribution.
+
+The first idea that one might think of is to take a much smaller :math:`\alpha`: for instance, if we take, :math:`\alpha=\frac{0.05}{n\_voxels}` then the expected number of false discoveries is only about 0.05, meaning that there is a 5% chance to declare active a truly inactive voxel. This correction on the signifiance is known as Bonferroni procedure. It is fairly accurate when the different tests are independent or close to independent, and becomes conservative otherwise.
+The problem with his approach is that truly activated voxel may not surpass the corresponding threshold, which is typically very high, because `n\_voxels` is large.
+
+A second possibility is to choose a threshold so that the proportion of true discoveries among the discoveries reaches a certain proportion `0`_ gallery, especially those of the Tutorial section.
.. _installation:
-Installing nistats
-====================
+Installing Nistats
+==================
.. raw:: html
:file: install_doc_component.html
diff --git a/doc/modules/reference.rst b/doc/modules/reference.rst
index e49d5f57..efbc480d 100644
--- a/doc/modules/reference.rst
+++ b/doc/modules/reference.rst
@@ -78,9 +78,9 @@ uses.
:toctree: generated/
:template: function.rst
- make_design_matrix
+ make_first_level_design_matrix
check_design_matrix
- create_second_level_design
+ make_second_level_design_matrix
.. _experimental_paradigm_ref:
@@ -99,8 +99,7 @@ uses.
:toctree: generated/
:template: function.rst
- check_paradigm
- paradigm_from_csv
+ check_events
.. _model_ref:
@@ -243,7 +242,6 @@ uses.
fdr_threshold
map_threshold
- get_clusters_table
.. _reporting_ref:
@@ -265,6 +263,7 @@ uses.
compare_niimgs
plot_design_matrix
plot_contrast_matrix
+ get_clusters_table
.. _utils_ref:
@@ -284,9 +283,9 @@ uses.
:template: function.rst
z_score
- multiple_fast_inv
+ multiple_fast_inverse
multiple_mahalanobis
full_rank
- pos_recipr
+ positive_reciprocal
get_bids_files
parse_bids_filename
diff --git a/doc/themes/nistats/layout.html b/doc/themes/nistats/layout.html
index bcb354b8..cfacc5a8 100644
--- a/doc/themes/nistats/layout.html
+++ b/doc/themes/nistats/layout.html
@@ -152,7 +152,7 @@
First Level
diff --git a/doc/user_guide.rst b/doc/user_guide.rst
index aca2549f..ccd62904 100644
--- a/doc/user_guide.rst
+++ b/doc/user_guide.rst
@@ -19,5 +19,5 @@ User guide: table of contents
:numbered:
introduction.rst
- modules/reference.rst
auto_examples/index.rst
+ modules/reference.rst
diff --git a/doc/whats_new.rst b/doc/whats_new.rst
index 17cb6969..5849e2f2 100644
--- a/doc/whats_new.rst
+++ b/doc/whats_new.rst
@@ -1,3 +1,60 @@
+0.0.1b
+=======
+
+Changelog
+---------
+
+* Updated the minimum versions of the dependencies
+ * Numpy >= 1.11
+ * SciPy >= 0.17
+ * Nibabel >= 2.0.2
+ * Nilearn >= 0.4.0
+ * Pandas >= 0.18.0
+ * Sklearn >= 0.18.0
+
+* Added comprehensive tutorial
+
+* Second-level model accepts 4D images as input.
+
+* Changes in function parameters
+ * third argument of map_threshold is now called ``level``.
+ * Changed the defaut oversampling value for the hemodynamic response
+ to 50 and exposed this parameter.
+ * changed the term ``paradigm`` to ``events`` and made it
+ BIDS-compliant. Set the event file to be tab-separated
+
+* Certain functions and methods have been renamed for clarity
+ * ``nistats.design_matrix``
+ * ``make_design_matrix() -> make_first_level_design_matrix()``
+ * ``create_second_level_design() -> make_second_level_design_matrix()``
+ * ``nistats.utils``
+ * ``pos_recipr() -> positive_reciprocal()``
+ * ``multiple_fast_inv() -> multiple_fast_inverse()``
+
+* Python2 Deprecation:
+ Python 2 is now deprecated and will not be supported in a future version.
+ A DeprecationWarning is displayed in Python 2 environments with a suggestion to move to Python 3.
+
+
+Contributors
+------------
+
+The following people contributed to this release::
+
+ 45 Bertrand Thirion
+ 70 Kshitij Chawla
+ 16 Taylor Salo
+ 6 KamalakerDadi
+ 5 chrplr
+ 5 hcherkaoui
+ 5 rschmaelzle
+ 4 mannalytics
+ 3 Martin Perez-Guevara
+ 2 Christopher J. Markiewicz
+ 1 Loïc Estève
+
+
+
0.0.1a
=======
diff --git a/examples/01_tutorials/auditory_block_paradigm.csv b/examples/01_tutorials/auditory_block_paradigm.csv
new file mode 100644
index 00000000..9511381d
--- /dev/null
+++ b/examples/01_tutorials/auditory_block_paradigm.csv
@@ -0,0 +1,17 @@
+duration,onset,trial_type
+42.0,0.0,rest
+42.0,42.0,active
+42.0,84.0,rest
+42.0,126.0,active
+42.0,168.0,rest
+42.0,210.0,active
+42.0,252.0,rest
+42.0,294.0,active
+42.0,336.0,rest
+42.0,378.0,active
+42.0,420.0,rest
+42.0,462.0,active
+42.0,504.0,rest
+42.0,546.0,active
+42.0,588.0,rest
+42.0,630.0,active
diff --git a/examples/01_tutorials/plot_first_level_model_details.py b/examples/01_tutorials/plot_first_level_model_details.py
new file mode 100644
index 00000000..64927f8c
--- /dev/null
+++ b/examples/01_tutorials/plot_first_level_model_details.py
@@ -0,0 +1,470 @@
+"""Studying first-level-model details in a trials-and-error fashion
+===================================================================
+
+In this tutorial, we study the parametrization of the first-level
+model used for fMRI data analysis and clarify their impact on the
+results of the analysis.
+
+We use an exploratory approach, in which we incrementally include some
+new features in the analysis and look at the outcome, i.e. the
+resulting brain maps.
+
+Readers without prior experience in fMRI data analysis should first
+run the `Analysis of a single session, single subject fMRI dataset`_
+tutorial to get a bit more
+familiar with the base concepts, and only then run this tutorial example.
+
+To run this example, you must launch IPython via ``ipython
+--matplotlib`` in a terminal, or use ``jupyter-notebook``.
+
+.. contents:: **Contents**
+ :local:
+ :depth: 1
+
+.. _Analysis of a single session, single subject fMRI dataset: plot_single_subject_single_run.html
+
+"""
+
+###############################################################################
+# Retrieving the data
+# -------------------
+#
+# We use a so-called localizer dataset, which consists in a 5-minutes
+# acquisition of a fast event-related dataset.
+#
+from nistats import datasets
+data = datasets.fetch_localizer_first_level()
+fmri_img = data.epi_img
+
+###############################################################################
+# Define the paradigm that will be used. Here, we just need to get the provided file.
+#
+# This task, described in Pinel et al., BMC neuroscience 2007 probes
+# basic functions, such as button with the left or right hand, viewing
+# horizontal and vertical checkerboards, reading and listening to
+# short sentences, and mental computations (subractions).
+#
+# Visual stimuli were displayed in four 250-ms epochs, separated by
+# 100ms intervals (i.e., 1.3s in total). Auditory stimuli were drawn
+# from a recorded male voice (i.e., a total of 1.6s for motor
+# instructions, 1.2-1.7s for sentences, and 1.2-1.3s for
+# subtraction). The auditory or visual stimuli were shown to the
+# participants for passive viewing or button response in
+# event-related paradigms. Post-scan questions verified that the
+# experimental tasks were understood and followed correctly.
+#
+# This task comprises 10 conditions:
+#
+# * clicGaudio: Left-hand three-times button press, indicated by visual instruction
+# * clicDaudio: Right-hand three-times button press, indicated by visual instruction
+# * clicGvideo: Left-hand three-times button press, indicated by auditory instruction
+# * clicDvideo: Right-hand three-times button press, indicated by auditory instruction
+# * damier_H: Visualization of flashing horizontal checkerboards
+# * damier_V: Visualization of flashing vertical checkerboards
+# * phraseaudio: Listen to narrative sentences
+# * phrasevideo: Read narrative sentences
+# * calculaudio: Mental subtraction, indicated by auditory instruction
+# * calculvideo: Mental subtraction, indicated by visual instruction
+#
+
+t_r = 2.4
+events_file = data['events']
+import pandas as pd
+events= pd.read_table(events_file)
+
+###############################################################################
+# Running a basic model
+# ---------------------
+#
+# First specify a linear model.
+# the fit() model creates the design matrix and the beta maps.
+#
+from nistats.first_level_model import FirstLevelModel
+first_level_model = FirstLevelModel(t_r)
+first_level_model = first_level_model.fit(fmri_img, events=events)
+design_matrix = first_level_model.design_matrices_[0]
+
+#########################################################################
+# Let us take a look at the design matrix: it has 10 main columns corresponding to 10 experimental conditions, followed by 3 columns describing low-frequency signals (drifts) and a constant regressor.
+from nistats.reporting import plot_design_matrix
+plot_design_matrix(design_matrix)
+import matplotlib.pyplot as plt
+plt.show()
+
+#########################################################################
+# Specification of the contrasts.
+#
+# For this, let's create a function that, given the design matrix,
+# generates the corresponding contrasts. This will be useful to
+# repeat contrast specification when we change the design matrix.
+import numpy as np
+
+def make_localizer_contrasts(design_matrix):
+ """ returns a dictionary of four contrasts, given the design matrix"""
+
+ # first generate canonical contrasts
+ contrast_matrix = np.eye(design_matrix.shape[1])
+ contrasts = dict([(column, contrast_matrix[i])
+ for i, column in enumerate(design_matrix.columns)])
+
+ # Add more complex contrasts
+ contrasts['audio'] = (contrasts['clicDaudio']
+ + contrasts['clicGaudio']
+ + contrasts['calculaudio']
+ + contrasts['phraseaudio']
+ )
+ contrasts['video'] = (contrasts['clicDvideo']
+ + contrasts['clicGvideo']
+ + contrasts['calculvideo']
+ + contrasts['phrasevideo']
+ )
+ contrasts['computation'] = contrasts['calculaudio'] + contrasts['calculvideo']
+ contrasts['sentences'] = contrasts['phraseaudio'] + contrasts['phrasevideo']
+
+ # Short dictionary of more relevant contrasts
+ contrasts = {
+ 'left-right': (contrasts['clicGaudio']
+ + contrasts['clicGvideo']
+ - contrasts['clicDaudio']
+ - contrasts['clicDvideo']
+ ),
+ 'H-V': contrasts['damier_H'] - contrasts['damier_V'],
+ 'audio-video': contrasts['audio'] - contrasts['video'],
+ 'computation-sentences': (contrasts['computation'] -
+ contrasts['sentences']
+ ),
+ }
+ return contrasts
+
+#########################################################################
+# So let's look at these computed contrasts
+#
+# * "left - right button press" probes motor activity in left versus right button presses
+# * 'H-V': probes the differential activity in viewing a horizontal vs vertical checkerboard
+# * "audio - video" probes the difference of activity between listening to some content or reading the same type of content (instructions, stories)
+# * "computation - sentences" looks at the activity when performing a mental comptation task versus simply reading sentences.
+#
+contrasts = make_localizer_contrasts(design_matrix)
+plt.figure(figsize=(5, 9))
+from nistats.reporting import plot_contrast_matrix
+for i, (key, values) in enumerate(contrasts.items()):
+ ax = plt.subplot(len(contrasts) + 1, 1, i + 1)
+ plot_contrast_matrix(values, design_matrix=design_matrix, ax=ax)
+
+plt.show()
+
+#########################################################################
+# Contrast estimation and plotting
+#
+# Since this script will be repeated several times, for the sake of readability,
+# we encapsulate it in a function that we call when needed.
+#
+from nilearn import plotting
+
+def plot_contrast(first_level_model):
+ """ Given a first model, specify, enstimate and plot the main contrasts"""
+ design_matrix = first_level_model.design_matrices_[0]
+ # Call the contrast specification within the function
+ contrasts = make_localizer_contrasts(design_matrix)
+ fig = plt.figure(figsize=(11, 3))
+ # compute the per-contrast z-map
+ for index, (contrast_id, contrast_val) in enumerate(contrasts.items()):
+ ax = plt.subplot(1, len(contrasts), 1 + index)
+ z_map = first_level_model.compute_contrast(
+ contrast_val, output_type='z_score')
+ plotting.plot_stat_map(
+ z_map, display_mode='z', threshold=3.0, title=contrast_id, axes=ax,
+ cut_coords=1)
+
+#########################################################################
+# Let's run the model and look at the outcome.
+
+plot_contrast(first_level_model)
+plt.show()
+
+#########################################################################
+# Changing the drift model
+# ------------------------
+#
+# By default the drift model is a set of slow oscillating functions (Discrete Cosine transform), with a cutoff at frequency 1/128 hz.
+# We can change this cut-off, e.g. to 1/64Hz.
+# This is done by setting period_cut=64(s)
+
+first_level_model = FirstLevelModel(t_r, period_cut=64)
+first_level_model = first_level_model.fit(fmri_img, events=events)
+design_matrix = first_level_model.design_matrices_[0]
+plot_design_matrix(design_matrix)
+
+#########################################################################
+# Does the model perform worse or better ?
+
+plot_contrast(first_level_model)
+plt.show()
+
+#########################################################################
+# Note that the design matrix has more columns to model drifts in the data.
+#
+# We notice however that this model performs rather poorly.
+#
+# Another solution is to remove these drift terms. Maybe they're simply useless.
+# this is done by setting drift_model to None.
+
+first_level_model = FirstLevelModel(t_r, drift_model=None)
+first_level_model = first_level_model.fit(fmri_img, events=events)
+design_matrix = first_level_model.design_matrices_[0]
+plot_design_matrix(design_matrix)
+plot_contrast(first_level_model)
+plt.show()
+
+#########################################################################
+# Is it better than the original ? No !
+#
+# Note that the design matrix has changed with no drift columns.
+# the event columns, on the other hand, haven't changed.
+# Another alternative to get a drift model is to specify a set of polynomials
+# Let's take a basis of 5 polynomials
+
+first_level_model = FirstLevelModel(t_r, drift_model='polynomial',
+ drift_order=5)
+first_level_model = first_level_model.fit(fmri_img, events=events)
+design_matrix = first_level_model.design_matrices_[0]
+plot_design_matrix(design_matrix)
+plot_contrast(first_level_model)
+plt.show()
+
+#########################################################################
+# Is it good ? No better, no worse. Let's turn to another parameter.
+
+#########################################################################
+# Changing the hemodynamic response model
+# ---------------------------------------
+#
+# This is the filter used to convert the event sequence into a
+# reference BOLD signal for the design matrix.
+#
+# The first thing that we can do is to change the default model (the
+# so-called Glover hrf) for the so-called canonical model of SPM
+# --which has slightly weaker undershoot component.
+
+first_level_model = FirstLevelModel(t_r, hrf_model='spm')
+first_level_model = first_level_model.fit(fmri_img, events=events)
+design_matrix = first_level_model.design_matrices_[0]
+plot_design_matrix(design_matrix)
+plot_contrast(first_level_model)
+plt.show()
+
+#########################################################################
+# No strong --positive or negative-- effect.
+#
+# We could try to go one step further: using not only the so-called
+# canonical hrf, but also its time derivative. Note that in that case,
+# we still perform the contrasts and obtain statistical significance
+# for the main effect ---not the time derivative. This means that the
+# inclusion of time derivative in the design matrix has the sole
+# effect of discounting timing misspecification from the error term,
+# which vould decrease the estimated variance and enhance the
+# statistical significance of the effect. Is it the case ?
+
+first_level_model = FirstLevelModel(t_r, hrf_model='spm + derivative')
+first_level_model = first_level_model.fit(fmri_img, events=events)
+design_matrix = first_level_model.design_matrices_[0]
+plot_design_matrix(design_matrix)
+plot_contrast(first_level_model)
+plt.show()
+
+#########################################################################
+# Not a huge effect, but rather positive overall. We could keep that one.
+#
+# Bzw, a benefit of this approach is that we can test which voxels are
+# well explined by the derivative term, hinting at misfit regions, a
+# possibly valuable information This is implemented by an F test on
+# the time derivative regressors.
+
+contrast_val = np.eye(design_matrix.shape[1])[1:21:2]
+plot_contrast_matrix(contrast_val, design_matrix)
+plt.show()
+
+z_map = first_level_model.compute_contrast(
+ contrast_val, output_type='z_score')
+plotting.plot_stat_map(
+ z_map, display_mode='z', threshold=3.0, title='effect of time derivatives')
+plt.show()
+
+#########################################################################
+# Well, there seems to be something here. Maybe we could adjust the
+# timing, by increasing the slice_time_ref parameter: 0 to 0.5 now the
+# reference for model sampling is not the beginning of the volume
+# acquisition, but the middle of it.
+first_level_model = FirstLevelModel(t_r, hrf_model='spm + derivative',
+ slice_time_ref=0.5)
+first_level_model = first_level_model.fit(fmri_img, events=events)
+z_map = first_level_model.compute_contrast(
+ contrast_val, output_type='z_score')
+plotting.plot_stat_map(
+ z_map, display_mode='z', threshold=3.0,
+ title='effect of time derivatives after model shift')
+plt.show()
+#########################################################################
+# The time derivatives regressors capture less signal: it's better so.
+
+#########################################################################
+# We can also consider adding the so-called dispersion derivative to
+# capture some mis-specification in the shape of the hrf.
+#
+# This is done by specifying `hrf_model='spm + derivative + dispersion'`
+#
+first_level_model = FirstLevelModel(t_r,slice_time_ref=0.5,
+ hrf_model='spm + derivative + dispersion')
+first_level_model = first_level_model.fit(fmri_img, events=events)
+design_matrix = first_level_model.design_matrices_[0]
+plot_design_matrix(design_matrix)
+plot_contrast(first_level_model)
+plt.show()
+
+#########################################################################
+# Not a huge effect. For the sake of simplicity and readibility, we
+# can drop that one.
+
+#########################################################################
+# The noise model ar(1) or ols ?
+# ------------------------------
+#
+# So far,we have implicitly used a lag-1 autoregressive model ---aka
+# ar(1)--- for the temporal structure of the noise. An alternative
+# choice is to use an ordinaly least squares model (ols) that assumes
+# no temporal structure (time-independent noise)
+
+first_level_model = FirstLevelModel(t_r, slice_time_ref=0.5,
+ hrf_model='spm + derivative',
+ noise_model='ols')
+first_level_model = first_level_model.fit(fmri_img, events=events)
+design_matrix = first_level_model.design_matrices_[0]
+plot_design_matrix(design_matrix)
+plot_contrast(first_level_model)
+plt.show()
+
+#########################################################################
+# While the difference is not obvious you should rather stick to the
+# ar(1) model, which is arguably more accurate.
+
+#########################################################################
+# Removing confounds
+# ------------------
+#
+# A problematic feature of fMRI is the presence of unconctrolled
+# confounds in the data, sue to scanner instabilities (spikes) or
+# physiological phenomena, such as motion, heart and
+# respiration-related blood oxygenation flucturations. Side
+# measurements are sometimes acquired to charcterise these
+# effects. Here we don't have access to those. What we can do instead
+# is to estimate confounding effects from the data themselves, using
+# the compcorr approach, and take those into account in the model.
+#
+# For this we rely on the so-called `high_variance_confounds`_
+# routine of Nilearn.
+#
+# .. _high_variance_confounds: https://nilearn.github.io/modules/generated/nilearn.image.high_variance_confounds.html
+
+
+from nilearn.image import high_variance_confounds
+confounds = pd.DataFrame(high_variance_confounds(fmri_img, percentile=1))
+first_level_model = FirstLevelModel(t_r, hrf_model='spm + derivative',
+ slice_time_ref=0.5)
+first_level_model = first_level_model.fit(fmri_img, events=events,
+ confounds=confounds)
+design_matrix = first_level_model.design_matrices_[0]
+plot_design_matrix(design_matrix)
+plot_contrast(first_level_model)
+plt.show()
+
+#########################################################################
+# Note the five additional columns in the design matrix
+#
+# The effect on activation maps is complex: auditory/visual effects
+# are killed, probably because they were somewhat colinear to the
+# confounds. On the other hand, some of the maps become cleaner (H-V,
+# computation) after this addition.
+
+
+#########################################################################
+# Smoothing
+# ----------
+#
+# Smoothing is a regularization of the model. It has two benefits:
+# decrease the noise level in images, and reduce the discrepancy
+# between individuals. The drawback is that it biases the shape and
+# position of activation. We simply illustrate here the statistical
+# gains. We use a mild smoothing of 5mm full-width at half maximum
+# (fwhm).
+
+first_level_model = FirstLevelModel(
+ t_r, hrf_model='spm + derivative', smoothing_fwhm=5,
+ slice_time_ref=0.5).fit(fmri_img, events=events, confounds=confounds)
+design_matrix = first_level_model.design_matrices_[0]
+plot_design_matrix(design_matrix)
+plot_contrast(first_level_model)
+plt.show()
+
+#########################################################################
+# The design is unchanged but the maps are smoother and more contrasted
+#
+
+#########################################################################
+# Masking
+# --------
+#
+# Masking consists in selecting the region of the image on which the
+# model is run: it is useless to run it outside of the brain.
+#
+# The approach taken by FirstLeveModel is to estimate it from the fMRI
+# data themselves when no mask is explicitly provided. Since the data
+# have been resampled into MNI space, we can use instead a mask of the
+# grey matter in MNI space. The benefit is that it makes voxel-level
+# comparisons easier across subjects and datasets, and removed
+# non-grey matter regions, in which no BOLD signal is expected. The
+# downside is that the mask may not fit very well these particular
+# data.
+
+data_mask = first_level_model.masker_.mask_img_
+from nilearn.datasets import fetch_icbm152_brain_gm_mask
+icbm_mask = fetch_icbm152_brain_gm_mask()
+
+from nilearn.plotting import plot_roi
+plt.figure(figsize=(16, 4))
+ax = plt.subplot(121)
+plot_roi(icbm_mask, title='ICBM mask', axes=ax)
+ax = plt.subplot(122)
+plot_roi(data_mask, title='Data-driven mask', axes=ax)
+plt.show()
+
+#########################################################################
+# For the sake of time saving, we reample icbm_mask to our data
+# For this we call the resample_to_img routine of Nilearn.
+# We use interpolation = 'nearest' to keep the mask a binary image
+from nilearn.image import resample_to_img
+resampled_icbm_mask = resample_to_img(icbm_mask, data_mask,
+ interpolation='nearest')
+
+#########################################################################
+# Impact on the first-level model
+first_level_model = FirstLevelModel(
+ t_r, hrf_model='spm + derivative', smoothing_fwhm=5, slice_time_ref=0.5,
+ mask=resampled_icbm_mask).fit(
+ fmri_img, events=events, confounds=confounds)
+design_matrix = first_level_model.design_matrices_[0]
+plot_design_matrix(design_matrix)
+plot_contrast(first_level_model)
+plt.show()
+
+#########################################################################
+# Note that it removed suprious spots in the white matter.
+
+#########################################################################
+# Conclusion
+# ----------
+#
+# Interestingly, the model used here seems quite resilient to
+# manipulation of modeling parameters: this is reassuring. It shows
+# that Nistats defaults ('cosine' drift, cutoff=128s, 'glover' hrf,
+# ar(1) model) are actually reasonable. Note that these conclusions
+# are specific to this dataset and may vary with other ones.
diff --git a/examples/01_tutorials/plot_single_subject_single_run.py b/examples/01_tutorials/plot_single_subject_single_run.py
new file mode 100644
index 00000000..73c22489
--- /dev/null
+++ b/examples/01_tutorials/plot_single_subject_single_run.py
@@ -0,0 +1,316 @@
+"""
+Analysis of a single session, single subject fMRI dataset
+=========================================================
+
+In this tutorial, we compare the fMRI signal during periods of auditory
+stimulation versus periods of rest, using a General Linear Model (GLM).
+
+The dataset comes from an experiment conducted at the FIL by Geriant Rees
+under the direction of Karl Friston. It is provided by FIL methods
+group which develops the SPM software.
+
+According to SPM documentation, 96 scans were acquired (repetition time TR=7s) in one session. The paradigm consisted of alternating periods of stimulation and rest, lasting 42s each (that is, for 6 scans). The sesssion started with a rest block.
+Auditory stimulation consisted of bi-syllabic words presented binaurally at a
+rate of 60 per minute. The functional data starts at scan number 4, that is the
+image file ``fM00223_004``.
+
+The whole brain BOLD/EPI images were acquired on a 2T Siemens
+MAGNETOM Vision system. Each scan consisted of 64 contiguous
+slices (64x64x64 3mm x 3mm x 3mm voxels). Acquisition of one scan took 6.05s, with the scan to scan repeat time (TR) set arbitrarily to 7s.
+
+The analyse described here is performed in the native space, directly on the
+original EPI scans without any spatial or temporal preprocessing.
+(More sensitive results would likely be obtained on the corrected,
+spatially normalized and smoothed images).
+
+
+To run this example, you must launch IPython via ``ipython
+--matplotlib`` in a terminal, or use ``jupyter-notebook``.
+
+.. contents:: **Contents**
+ :local:
+ :depth: 1
+
+"""
+
+###############################################################################
+# Retrieving the data
+# -------------------
+#
+# .. note:: In this tutorial, we load the data using a data downloading
+# function. To input your own data, you will need to provide
+# a list of paths to your own files in the ``subject_data`` variable.
+# These should abide to the Brain Imaging Data Structure (BIDS)
+# organization.
+
+from nistats.datasets import fetch_spm_auditory
+subject_data = fetch_spm_auditory()
+print(subject_data.func) # print the list of names of functional images
+
+###############################################################################
+# We can display the first functional image and the subject's anatomy:
+from nilearn.plotting import plot_stat_map, plot_anat, plot_img, show
+plot_img(subject_data.func[0])
+plot_anat(subject_data.anat)
+
+###############################################################################
+# Next, we concatenate all the 3D EPI image into a single 4D image,
+# then we average them in order to create a background
+# image that will be used to display the activations:
+
+from nilearn.image import concat_imgs, mean_img
+fmri_img = concat_imgs(subject_data.func)
+mean_img = mean_img(fmri_img)
+
+###############################################################################
+# Specifying the experimental paradigm
+# ------------------------------------
+#
+# We must now provide a description of the experiment, that is, define the
+# timing of the auditory stimulation and rest periods. This is typically
+# provided in an events.tsv file. The path of this file is
+# provided in the dataset.
+import pandas as pd
+events = pd.read_table(subject_data['events'])
+print(events)
+
+###############################################################################
+# Performing the GLM analysis
+# ---------------------------
+#
+# It is now time to create and estimate a ``FirstLevelModel`` object, that will generate the *design matrix* using the information provided by the ``events`` object.
+
+from nistats.first_level_model import FirstLevelModel
+
+###############################################################################
+# Parameters of the first-level model
+#
+# * t_r=7(s) is the time of repetition of acquisitions
+# * noise_model='ar1' specifies the noise covariance model: a lag-1 dependence
+# * standardize=False means that we do not want to rescale the time series to mean 0, variance 1
+# * hrf_model='spm' means that we rely on the SPM "canonical hrf" model (without time or dispersion derivatives)
+# * drift_model='cosine' means that we model the signal drifts as slow oscillating time functions
+# * period_cut=160(s) defines the cutoff frequency (its inverse actually).
+fmri_glm = FirstLevelModel(t_r=7,
+ noise_model='ar1',
+ standardize=False,
+ hrf_model='spm',
+ drift_model='cosine',
+ period_cut=160)
+
+###############################################################################
+# Now that we have specified the model, we can run it on the fMRI image
+fmri_glm = fmri_glm.fit(fmri_img, events)
+
+###############################################################################
+# One can inspect the design matrix (rows represent time, and
+# columns contain the predictors).
+design_matrix = fmri_glm.design_matrices_[0]
+
+###############################################################################
+# Formally, we have taken the first design matrix, because the model is
+# implictily meant to for multiple runs.
+from nistats.reporting import plot_design_matrix
+plot_design_matrix(design_matrix)
+import matplotlib.pyplot as plt
+plt.show()
+
+###############################################################################
+# Save the design matrix image to disk
+# first create a directory where you want to write the images
+
+import os
+outdir = 'results'
+if not os.path.exists(outdir):
+ os.mkdir(outdir)
+
+from os.path import join
+plot_design_matrix(design_matrix, output_file=join(outdir, 'design_matrix.png'))
+
+###############################################################################
+# The first column contains the expected reponse profile of regions which are
+# sensitive to the auditory stimulation.
+# Let's plot this first column
+
+plt.plot(design_matrix['active'])
+plt.xlabel('scan')
+plt.title('Expected Auditory Response')
+plt.show()
+
+###############################################################################
+# Detecting voxels with significant effects
+# -----------------------------------------
+#
+# To access the estimated coefficients (Betas of the GLM model), we
+# created constrast with a single '1' in each of the columns: The role
+# of the contrast is to select some columns of the model --and
+# potentially weight them-- to study the associated statistics. So in
+# a nutshell, a contrast is a weigted combination of the estimated
+# effects. Here we can define canonical contrasts that just consider
+# the two condition in isolation ---let's call them "conditions"---
+# then a contrast that makes the difference between these conditions.
+
+from numpy import array
+conditions = {
+ 'active': array([1., 0., 0., 0., 0., 0., 0., 0., 0., 0.]),
+ 'rest': array([0., 1., 0., 0., 0., 0., 0., 0., 0., 0.]),
+}
+
+###############################################################################
+# We can then compare the two conditions 'active' and 'rest' by
+# defining the corresponding contrast:
+
+active_minus_rest = conditions['active'] - conditions['rest']
+
+###############################################################################
+# Let's look at it: plot the coefficients of the contrast, indexed by
+# the names of the columns of the design matrix.
+
+from nistats.reporting import plot_contrast_matrix
+plot_contrast_matrix(active_minus_rest, design_matrix=design_matrix)
+
+###############################################################################
+# Below, we compute the estimated effect. It is in BOLD signal unit,
+# but has no statistical guarantees, because it does not take into
+# account the associated variance.
+
+eff_map = fmri_glm.compute_contrast(active_minus_rest,
+ output_type='effect_size')
+
+###############################################################################
+# In order to get statistical significance, we form a t-statistic, and
+# directly convert is into z-scale. The z-scale means that the values
+# are scaled to match a standard Gaussian distribution (mean=0,
+# variance=1), across voxels, if there were now effects in the data.
+
+z_map = fmri_glm.compute_contrast(active_minus_rest,
+ output_type='z_score')
+
+###############################################################################
+# Plot thresholded z scores map.
+#
+# We display it on top of the average
+# functional image of the series (could be the anatomical image of the
+# subject). We use arbitrarily a threshold of 3.0 in z-scale. We'll
+# see later how to use corrected thresholds. we show to display 3
+# axial views: display_mode='z', cut_coords=3
+
+plot_stat_map(z_map, bg_img=mean_img, threshold=3.0,
+ display_mode='z', cut_coords=3, black_bg=True,
+ title='Active minus Rest (Z>3)')
+plt.show()
+
+###############################################################################
+# Statistical signifiance testing. One should worry about the
+# statistical validity of the procedure: here we used an arbitrary
+# threshold of 3.0 but the threshold should provide some guarantees on
+# the risk of false detections (aka type-1 errors in statistics). One
+# first suggestion is to control the false positive rate (fpr) at a
+# certain level, e.g. 0.001: this means that there is.1% chance of
+# declaring active an inactive voxel.
+
+from nistats.thresholding import map_threshold
+_, threshold = map_threshold(z_map, level=.001, height_control='fpr')
+print('Uncorrected p<0.001 threshold: %.3f' % threshold)
+plot_stat_map(z_map, bg_img=mean_img, threshold=threshold,
+ display_mode='z', cut_coords=3, black_bg=True,
+ title='Active minus Rest (p<0.001)')
+plt.show()
+
+###############################################################################
+# The problem is that with this you expect 0.001 * n_voxels to show up
+# while they're not active --- tens to hundreds of voxels. A more
+# conservative solution is to control the family wise errro rate,
+# i.e. the probability of making ony one false detection, say at
+# 5%. For that we use the so-called Bonferroni correction
+
+_, threshold = map_threshold(z_map, level=.05, height_control='bonferroni')
+print('Bonferroni-corrected, p<0.05 threshold: %.3f' % threshold)
+plot_stat_map(z_map, bg_img=mean_img, threshold=threshold,
+ display_mode='z', cut_coords=3, black_bg=True,
+ title='Active minus Rest (p<0.05, corrected)')
+plt.show()
+
+###############################################################################
+# This is quite conservative indeed ! A popular alternative is to
+# control the false discovery rate, i.e. the expected proportion of
+# false discoveries among detections. This is called the false
+# disovery rate
+
+_, threshold = map_threshold(z_map, level=.05, height_control='fdr')
+print('False Discovery rate = 0.05 threshold: %.3f' % threshold)
+plot_stat_map(z_map, bg_img=mean_img, threshold=threshold,
+ display_mode='z', cut_coords=3, black_bg=True,
+ title='Active minus Rest (fdr=0.05)')
+plt.show()
+
+###############################################################################
+# Finally people like to discard isolated voxels (aka "small
+# clusters") from these images. It is possible to generate a
+# thresholded map with small clusters removed by providing a
+# cluster_threshold argument. here clusters smaller than 10 voxels
+# will be discarded.
+
+clean_map, threshold = map_threshold(
+ z_map, level=.05, height_control='fdr', cluster_threshold=10)
+plot_stat_map(clean_map, bg_img=mean_img, threshold=threshold,
+ display_mode='z', cut_coords=3, black_bg=True,
+ title='Active minus Rest (fdr=0.05), clusters > 10 voxels')
+plt.show()
+
+
+
+###############################################################################
+# We can save the effect and zscore maps to the disk
+z_map.to_filename(join(outdir, 'active_vs_rest_z_map.nii.gz'))
+eff_map.to_filename(join(outdir, 'active_vs_rest_eff_map.nii.gz'))
+
+###############################################################################
+# Report the found positions in a table
+
+from nistats.reporting import get_clusters_table
+table = get_clusters_table(z_map, stat_threshold=threshold,
+ cluster_threshold=20)
+print(table)
+
+###############################################################################
+# the table can be saved for future use
+
+table.to_csv(join(outdir, 'table.csv'))
+
+###############################################################################
+# Performing an F-test
+#
+# "active vs rest" is a typical t test: condition versus
+# baseline. Another popular type of test is an F test in which one
+# seeks whether a certain combination of conditions (possibly two-,
+# three- or higher-dimensional) explains a significant proportion of
+# the signal. Here one might for instance test which voxels are well
+# explained by combination of the active and rest condition.
+import numpy as np
+effects_of_interest = np.vstack((conditions['active'], conditions['rest']))
+plot_contrast_matrix(effects_of_interest, design_matrix)
+plt.show()
+
+###############################################################################
+# Specify the contrast and compute the correspoding map. Actually, the
+# contrast specification is done exactly the same way as for t
+# contrasts.
+
+z_map = fmri_glm.compute_contrast(effects_of_interest,
+ output_type='z_score')
+plt.show()
+
+###############################################################################
+# Note that the statistic has been converted to a z-variable, which
+# makes it easier to represent it.
+
+clean_map, threshold = map_threshold(
+ z_map, level=.05, height_control='fdr', cluster_threshold=10)
+plot_stat_map(clean_map, bg_img=mean_img, threshold=threshold,
+ display_mode='z', cut_coords=3, black_bg=True,
+ title='Effects of interest (fdr=0.05), clusters > 10 voxels')
+plt.show()
+
+###############################################################################
+# Oops, there is a lot of non-neural signal in there (ventricles, arteries)...
diff --git a/examples/01_tutorials/plot_spm_auditory.py b/examples/01_tutorials/plot_spm_auditory.py
deleted file mode 100644
index 678115a3..00000000
--- a/examples/01_tutorials/plot_spm_auditory.py
+++ /dev/null
@@ -1,253 +0,0 @@
-"""
-Univariate analysis of block design, one condition versus rest, single subject
-==============================================================================
-
-In this tutorial, we compare the fMRI signal during periods of auditory
-stimulation versus periods of rest, using a General Linear Model (GLM). We will
-use a univariate approach in which independent tests are performed at
-each single-voxel.
-
-The dataset comes from experiment conducted at the FIL by Geriant Rees
-under the direction of Karl Friston. It is provided by FIL methods
-group which develops the SPM software.
-
-According to SPM documentation, 96 acquisitions were made (RT=7s), in
-blocks of 6, giving 16 42s blocks. The condition for successive blocks
-alternated between rest and auditory stimulation, starting with rest.
-Auditory stimulation was bi-syllabic words presented binaurally at a
-rate of 60 per minute. The functional data starts at acquisiton 4,
-image fM00223_004.
-
-The whole brain BOLD/EPI images were acquired on a modified 2T Siemens
-MAGNETOM Vision system. Each acquisition consisted of 64 contiguous
-slices (64x64x64 3mm x 3mm x 3mm voxels). Acquisition took 6.05s, with
-the scan to scan repeat time (RT) set arbitrarily to 7s.
-
-
-This analyse described here is performed in the native space, on the
-original EPI scans without any spatial or temporal preprocessing.
-(More sensitive results would likely be obtained on the corrected,
-spatially normalized and smoothed images).
-
-
-To run this example, you must launch IPython via ``ipython
---matplotlib`` in a terminal, or use the Jupyter notebook.
-
-.. contents:: **Contents**
- :local:
- :depth: 1
-
-"""
-
-
-###############################################################################
-# Retrieving the data
-# -------------------
-#
-# .. note:: In this tutorial, we load the data using a data downloading
-# function.To input your own data, you will need to pass
-# a list of paths to your own files.
-
-from nistats.datasets import fetch_spm_auditory
-subject_data = fetch_spm_auditory()
-
-
-###############################################################################
-# We can list the filenames of the functional images
-print(subject_data.func)
-
-###############################################################################
-# Display the first functional image:
-from nilearn.plotting import plot_stat_map, plot_anat, plot_img
-plot_img(subject_data.func[0])
-
-###############################################################################
-# Display the subject's anatomical image:
-plot_anat(subject_data.anat)
-
-
-###############################################################################
-# Next, we concatenate all the 3D EPI image into a single 4D image:
-
-from nilearn.image import concat_imgs
-fmri_img = concat_imgs(subject_data.func)
-
-###############################################################################
-# And we average all the EPI images in order to create a background
-# image that will be used to display the activations:
-
-from nilearn import image
-mean_img = image.mean_img(fmri_img)
-
-###############################################################################
-# Specifying the experimental paradigm
-# ------------------------------------
-#
-# We must provide a description of the experiment, that is, define the
-# timing of the auditory stimulation and rest periods. According to
-# the documentation of the dataset, there were 16 42s blocks --- in
-# which 6 scans were acquired --- alternating between rest and
-# auditory stimulation, starting with rest. We use standard python
-# functions to create a pandas.DataFrame object that specifies the
-# timings:
-
-import numpy as np
-tr = 7.
-slice_time_ref = 0.
-n_scans = 96
-epoch_duration = 6 * tr # duration in seconds
-conditions = ['rest', 'active'] * 8
-n_blocks = len(conditions)
-duration = epoch_duration * np.ones(n_blocks)
-onset = np.linspace(0, (n_blocks - 1) * epoch_duration, n_blocks)
-
-import pandas as pd
-events = pd.DataFrame(
- {'onset': onset, 'duration': duration, 'trial_type': conditions})
-
-###############################################################################
-# The ``events`` object contains the information for the design:
-print(events)
-
-
-###############################################################################
-# Performing the GLM analysis
-# ---------------------------
-#
-# We need to construct a *design matrix* using the timing information
-# provided by the ``events`` object. The design matrix contains
-# regressors of interest as well as regressors of non-interest
-# modeling temporal drifts:
-
-frame_times = np.linspace(0, (n_scans - 1) * tr, n_scans)
-drift_model = 'Cosine'
-period_cut = 4. * epoch_duration
-hrf_model = 'glover + derivative'
-
-###############################################################################
-# It is now time to create a ``FirstLevelModel`` object
-# and fit it to the 4D dataset:
-
-from nistats.first_level_model import FirstLevelModel
-
-fmri_glm = FirstLevelModel(tr, slice_time_ref, noise_model='ar1',
- standardize=False, hrf_model=hrf_model,
- drift_model=drift_model, period_cut=period_cut)
-fmri_glm = fmri_glm.fit(fmri_img, events)
-
-###############################################################################
-# One can inspect the design matrix (rows represent time, and
-# columns contain the predictors):
-
-from nistats.reporting import plot_design_matrix
-design_matrix = fmri_glm.design_matrices_[0]
-plot_design_matrix(design_matrix)
-
-
-###############################################################################
-# The first column contains the expected reponse profile of regions which are
-# sensitive to the auditory stimulation.
-
-import matplotlib.pyplot as plt
-plt.plot(design_matrix['active'])
-plt.xlabel('scan')
-plt.title('Expected Auditory Response')
-plt.show()
-
-
-###############################################################################
-# Detecting voxels with significant effects
-# -----------------------------------------
-#
-# To access the estimated coefficients (Betas of the GLM model), we
-# created constrasts with a single '1' in each of the columns:
-
-contrast_matrix = np.eye(design_matrix.shape[1])
-contrasts = dict([(column, contrast_matrix[i])
- for i, column in enumerate(design_matrix.columns)])
-
-"""
-contrasts::
-
- {
- 'active': array([ 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]),
- 'active_derivative': array([ 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]),
- 'constant': array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.]),
- 'drift_1': array([ 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.]),
- 'drift_2': array([ 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.]),
- 'drift_3': array([ 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.]),
- 'drift_4': array([ 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.]),
- 'drift_5': array([ 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.]),
- 'drift_6': array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.]),
- 'drift_7': array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.]),
- 'rest': array([ 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.]),
- 'rest_derivative': array([ 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.])}
-"""
-
-###############################################################################
-# We can then compare the two conditions 'active' and 'rest' by
-# generating the relevant contrast:
-
-active_minus_rest = contrasts['active'] - contrasts['rest']
-
-eff_map = fmri_glm.compute_contrast(active_minus_rest,
- output_type='effect_size')
-
-z_map = fmri_glm.compute_contrast(active_minus_rest,
- output_type='z_score')
-
-###############################################################################
-# Plot thresholded z scores map
-
-plot_stat_map(z_map, bg_img=mean_img, threshold=3.0,
- display_mode='z', cut_coords=3, black_bg=True,
- title='Active minus Rest (Z>3)')
-
-
-###############################################################################
-# We can use ``nibabel.save`` to save the effect and zscore maps to the disk
-
-import os
-outdir = 'results'
-if not os.path.exists(outdir):
- os.mkdir(outdir)
-
-import nibabel
-from os.path import join
-nibabel.save(z_map, join('results', 'active_vs_rest_z_map.nii'))
-nibabel.save(eff_map, join('results', 'active_vs_rest_eff_map.nii'))
-
-###############################################################################
-# Extract the signal from a voxels
-# --------------------------------
-#
-# We search for the voxel with the larger z-score and plot the signal
-# (warning: double dipping!)
-
-
-# Find the coordinates of the peak
-
-from nibabel.affines import apply_affine
-values = z_map.get_data()
-coord_peaks = np.dstack(np.unravel_index(np.argsort(values.ravel()),
- values.shape))[0, 0, :]
-coord_mm = apply_affine(z_map.affine, coord_peaks)
-
-###############################################################################
-# We create a masker for the voxel (allowing us to detrend the signal)
-# and extract the time course
-
-from nilearn.input_data import NiftiSpheresMasker
-mask = NiftiSpheresMasker([coord_mm], radius=3,
- detrend=True, standardize=True,
- high_pass=None, low_pass=None, t_r=7.)
-sig = mask.fit_transform(fmri_img)
-
-##########################################################
-# Let's plot the signal and the theoretical response
-
-plt.plot(frame_times, sig, label='voxel %d %d %d' % tuple(coord_mm))
-plt.plot(design_matrix['active'], color='red', label='model')
-plt.xlabel('scan')
-plt.legend()
-plt.show()
diff --git a/examples/02_first_level_models/plot_adhd_dmn.py b/examples/02_first_level_models/plot_adhd_dmn.py
index e6fbef62..54cb9b44 100644
--- a/examples/02_first_level_models/plot_adhd_dmn.py
+++ b/examples/02_first_level_models/plot_adhd_dmn.py
@@ -9,9 +9,8 @@
1. A sequence of fMRI volumes are loaded
2. A design matrix with the Posterior Cingulate Cortex seed is defined
-4. A GLM is applied to the dataset (effect/covariance,
- then contrast estimation)
-5. The Default Mode Network is displayed
+3. A GLM is applied to the dataset (effect/covariance, then contrast estimation)
+4. The Default Mode Network is displayed
"""
import numpy as np
@@ -20,7 +19,7 @@
from nilearn.input_data import NiftiSpheresMasker
from nistats.first_level_model import FirstLevelModel
-from nistats.design_matrix import make_design_matrix
+from nistats.design_matrix import make_first_level_design_matrix
#########################################################################
# Prepare data and analysis parameters
@@ -47,9 +46,9 @@
memory_level=1, verbose=0)
seed_time_series = seed_masker.fit_transform(adhd_dataset.func[0])
frametimes = np.linspace(0, (n_scans - 1) * t_r, n_scans)
-design_matrix = make_design_matrix(frametimes, hrf_model='spm',
- add_regs=seed_time_series,
- add_reg_names=["pcc_seed"])
+design_matrix = make_first_level_design_matrix(frametimes, hrf_model='spm',
+ add_regs=seed_time_series,
+ add_reg_names=["pcc_seed"])
dmn_contrast = np.array([1] + [0]*(design_matrix.shape[1]-1))
contrasts = {'seed_based_glm': dmn_contrast}
diff --git a/examples/01_tutorials/plot_bids_features.py b/examples/02_first_level_models/plot_bids_features.py
similarity index 99%
rename from examples/01_tutorials/plot_bids_features.py
rename to examples/02_first_level_models/plot_bids_features.py
index fef976bd..40ad117c 100644
--- a/examples/01_tutorials/plot_bids_features.py
+++ b/examples/02_first_level_models/plot_bids_features.py
@@ -141,5 +141,5 @@
###############################################################################
# We can get a latex table from a Pandas Dataframe for display and publication
-from nistats.thresholding import get_clusters_table
+from nistats.reporting import get_clusters_table
print(get_clusters_table(z_map, norm.isf(0.001), 10).to_latex())
diff --git a/examples/02_first_level_models/plot_fiac_analysis.py b/examples/02_first_level_models/plot_fiac_analysis.py
index fa0a4fdd..2f4e6c2b 100644
--- a/examples/02_first_level_models/plot_fiac_analysis.py
+++ b/examples/02_first_level_models/plot_fiac_analysis.py
@@ -1,9 +1,9 @@
-"""
-Simple example of GLM fitting in fMRI
-======================================
+"""Simple example of two-session fMRI model fitting
+================================================
Full step-by-step example of fitting a GLM to experimental data and visualizing
the results. This is done on two runs of one subject of the FIAC dataset.
+
For details on the data, please see:
Dehaene-Lambertz G, Dehaene S, Anton JL, Campagne A, Ciuciu P, Dehaene
@@ -20,20 +20,17 @@
4. A GLM is applied to the dataset (effect/covariance,
then contrast estimation)
-"""
-from os import mkdir, path, getcwd
-
-import numpy as np
-import pandas as pd
-
-from nilearn import plotting
-from nilearn.image import mean_img
+Technically, this example shows how to handle two sessions that
+contain the same experimental conditions. The model directly returns a
+fixed effect of the statistics across the two sessions.
-from nistats.first_level_model import FirstLevelModel
-from nistats import datasets
+"""
-# write directory
+###############################################################################
+# Create a write directory to work
+# it will be a 'results' subdirectory of the current directory.
+from os import mkdir, path, getcwd
write_dir = path.join(getcwd(), 'results')
if not path.exists(write_dir):
mkdir(write_dir)
@@ -41,16 +38,31 @@
#########################################################################
# Prepare data and analysis parameters
# --------------------------------------
+#
+# Note that there are two sessions
+
+from nistats import datasets
data = datasets.fetch_fiac_first_level()
fmri_img = [data['func1'], data['func2']]
+
+#########################################################################
+# Create a mean image for plotting purpose
+from nilearn.image import mean_img
mean_img_ = mean_img(fmri_img[0])
+
+#########################################################################
+# The design matrices were pre-computed, we simply put them in a list of DataFrames
design_files = [data['design_matrix1'], data['design_matrix2']]
+import pandas as pd
+import numpy as np
design_matrices = [pd.DataFrame(np.load(df)['X']) for df in design_files]
#########################################################################
# GLM estimation
# ----------------------------------
-# GLM specification
+# GLM specification. Note that the mask was provided in the dataset. So we use it.
+
+from nistats.first_level_model import FirstLevelModel
fmri_glm = FirstLevelModel(mask=data['mask'], minimize_memory=True)
#########################################################################
@@ -58,13 +70,16 @@
fmri_glm = fmri_glm.fit(fmri_img, design_matrices=design_matrices)
#########################################################################
-# compute fixed effects of the two runs and compute related images
+# Compute fixed effects of the two runs and compute related images
+# For this, we first define the contrasts as we would do for a single session
n_columns = design_matrices[0].shape[1]
-
def pad_vector(contrast_, n_columns):
+ """A small routine to append zeros in contrast vectors"""
return np.hstack((contrast_, np.zeros(n_columns - len(contrast_))))
+#########################################################################
+# Contrast specification
contrasts = {'SStSSp_minus_DStDSp': pad_vector([1, 0, 0, -1], n_columns),
'DStDSp_minus_SStSSp': pad_vector([-1, 0, 0, 1], n_columns),
@@ -75,20 +90,60 @@ def pad_vector(contrast_, n_columns):
'Deactivation': pad_vector([-1, -1, -1, -1, 4], n_columns),
'Effects_of_interest': np.eye(n_columns)[:5]}
+#########################################################################
+# Compute and plot statistics
+
+from nilearn import plotting
print('Computing contrasts...')
for index, (contrast_id, contrast_val) in enumerate(contrasts.items()):
print(' Contrast % 2i out of %i: %s' % (
index + 1, len(contrasts), contrast_id))
- z_image_path = path.join(write_dir, '%s_z_map.nii' % contrast_id)
+ # estimate the contasts
+ # note that the model implictly compute a fixed effects across the two sessions
z_map = fmri_glm.compute_contrast(
contrast_val, output_type='z_score')
+
+ # Write the resulting stat images to file
+ z_image_path = path.join(write_dir, '%s_z_map.nii.gz' % contrast_id)
z_map.to_filename(z_image_path)
- # make a snapshot of the contrast activation
- if contrast_id == 'Effects_of_interest':
- display = plotting.plot_stat_map(
- z_map, bg_img=mean_img_, threshold=2.5, title=contrast_id)
- display.savefig(path.join(write_dir, '%s_z_map.png' % contrast_id))
+#########################################################################
+# Comparing session-specific and fixed effects.
+# Here, we compare the activation mas produced from each separately then, the fixed effects version
+
+contrast_id = 'Effects_of_interest'
+
+#########################################################################
+# Statistics for the first session
+
+fmri_glm = fmri_glm.fit(fmri_img[0], design_matrices=design_matrices[0])
+z_map = fmri_glm.compute_contrast(
+ contrasts[contrast_id], output_type='z_score')
+plotting.plot_stat_map(
+ z_map, bg_img=mean_img_, threshold=3.0,
+ title='%s, first session' % contrast_id)
+
+#########################################################################
+# Statistics for the second session
+
+fmri_glm = fmri_glm.fit(fmri_img[1], design_matrices=design_matrices[1])
+z_map = fmri_glm.compute_contrast(
+ contrasts[contrast_id], output_type='z_score')
+plotting.plot_stat_map(
+ z_map, bg_img=mean_img_, threshold=3.0,
+ title='%s, second session' % contrast_id)
+
+#########################################################################
+# Fixed effects statistics
+
+fmri_glm = fmri_glm.fit(fmri_img, design_matrices=design_matrices)
+z_map = fmri_glm.compute_contrast(
+ contrasts[contrast_id], output_type='z_score')
+plotting.plot_stat_map(
+ z_map, bg_img=mean_img_, threshold=3.0,
+ title='%s, fixed effects' % contrast_id)
+
+#########################################################################
+# Not unexpectedly, the fixed effects version looks displays higher peaks than the input sessions. Computing fixed effects enhances the signal-to-noise ratio of the resulting brain maps
-print('All the results were witten in %s' % write_dir)
plotting.show()
diff --git a/examples/02_first_level_models/plot_localizer_analysis.py b/examples/02_first_level_models/plot_localizer_analysis.py
deleted file mode 100644
index 647bb2d0..00000000
--- a/examples/02_first_level_models/plot_localizer_analysis.py
+++ /dev/null
@@ -1,89 +0,0 @@
-"""
-First level analysis of localizer dataset
-=========================================
-
-Full step-by-step example of fitting a GLM to experimental data and visualizing
-the results.
-
-More specifically:
-
-1. A sequence of fMRI volumes are loaded
-2. A design matrix describing all the effects related to the data is computed
-3. a mask of the useful brain volume is computed
-4. A GLM is applied to the dataset (effect/covariance,
- then contrast estimation)
-
-"""
-from os import mkdir, path
-
-import numpy as np
-import pandas as pd
-from nilearn import plotting
-
-from nistats.first_level_model import FirstLevelModel
-from nistats import datasets
-
-#########################################################################
-# Prepare data and analysis parameters
-# -------------------------------------
-# Prepare timing
-t_r = 2.4
-slice_time_ref = 0.5
-
-# Prepare data
-data = datasets.fetch_localizer_first_level()
-paradigm_file = data.paradigm
-paradigm = pd.read_csv(paradigm_file, sep=' ', header=None, index_col=None)
-paradigm.columns = ['session', 'trial_type', 'onset']
-fmri_img = data.epi_img
-
-#########################################################################
-# Perform first level analysis
-# ----------------------------
-# Setup and fit GLM
-first_level_model = FirstLevelModel(t_r, slice_time_ref,
- hrf_model='glover + derivative')
-first_level_model = first_level_model.fit(fmri_img, paradigm)
-
-#########################################################################
-# Estimate contrasts
-# ------------------
-# Specify the contrasts
-design_matrix = first_level_model.design_matrices_[0]
-contrast_matrix = np.eye(design_matrix.shape[1])
-contrasts = dict([(column, contrast_matrix[i])
- for i, column in enumerate(design_matrix.columns)])
-
-contrasts["audio"] = contrasts["clicDaudio"] + contrasts["clicGaudio"] +\
- contrasts["calculaudio"] + contrasts["phraseaudio"]
-contrasts["video"] = contrasts["clicDvideo"] + contrasts["clicGvideo"] + \
- contrasts["calculvideo"] + contrasts["phrasevideo"]
-contrasts["computation"] = contrasts["calculaudio"] + contrasts["calculvideo"]
-contrasts["sentences"] = contrasts["phraseaudio"] + contrasts["phrasevideo"]
-
-#########################################################################
-# Short list of more relevant contrasts
-contrasts = {
- "left-right": (contrasts["clicGaudio"] + contrasts["clicGvideo"]
- - contrasts["clicDaudio"] - contrasts["clicDvideo"]),
- "H-V": contrasts["damier_H"] - contrasts["damier_V"],
- "audio-video": contrasts["audio"] - contrasts["video"],
- "video-audio": -contrasts["audio"] + contrasts["video"],
- "computation-sentences": (contrasts["computation"] -
- contrasts["sentences"]),
- "reading-visual": contrasts["phrasevideo"] - contrasts["damier_H"]
- }
-
-#########################################################################
-# contrast estimation
-for index, (contrast_id, contrast_val) in enumerate(contrasts.items()):
- print(' Contrast % 2i out of %i: %s' %
- (index + 1, len(contrasts), contrast_id))
- z_map = first_level_model.compute_contrast(contrast_val,
- output_type='z_score')
-
- # Create snapshots of the contrasts
- display = plotting.plot_stat_map(z_map, display_mode='z',
- threshold=3.0, title=contrast_id)
-
-plotting.show()
diff --git a/examples/02_first_level_models/plot_localizer_surface_analysis.py b/examples/02_first_level_models/plot_localizer_surface_analysis.py
new file mode 100644
index 00000000..bdabe0f3
--- /dev/null
+++ b/examples/02_first_level_models/plot_localizer_surface_analysis.py
@@ -0,0 +1,205 @@
+"""Example of surface-based first-level analysis
+=============================================
+
+Full step-by-step example of fitting a GLM to experimental data
+sampled on the cortical surface and visualizing the results.
+
+More specifically:
+
+1. A sequence of fMRI volumes are loaded
+2. fMRI data are projected onto a reference cortical surface (the freesurfer template, fsaverage)
+3. A design matrix describing all the effects related to the data is computed
+4. A GLM is applied to the dataset (effect/covariance, then contrast estimation)
+
+The result of the analysis are statistical maps that are defined on
+the brain mesh. We display them using Nilearn capabilities.
+
+The projection of fMRI data onto a given brain mesh requires that both
+are initially defined in the same space.
+
+* The functional data should be coregistered to the anatomy from which the mesh was obtained.
+
+* Another possibility, used here, is to project the normalized fMRI data to an MNI-coregistered mesh, such as fsaverage.
+
+The advantage of this second approach is that it makes it easy to run
+second-level analyses on the surface. On the other hand, it is
+obviously less accurate than using a subject-tailored mesh.
+
+"""
+
+#########################################################################
+# Prepare data and analysis parameters
+# -------------------------------------
+# Prepare timing parameters
+t_r = 2.4
+slice_time_ref = 0.5
+
+#########################################################################
+# Prepare data
+# First the volume-based fMRI data.
+from nistats.datasets import fetch_localizer_first_level
+data = fetch_localizer_first_level()
+fmri_img = data.epi_img
+
+#########################################################################
+# Second the experimental paradigm.
+events_file = data['events']
+import pandas as pd
+events = pd.read_table(events_file)
+
+#########################################################################
+# Project the fMRI image to the surface
+# -------------------------------------
+#
+# For this we need to get a mesh representing the geometry of the
+# surface. we could use an individual mesh, but we first resort to a
+# standard mesh, the so-called fsaverage5 template from the Freesurfer
+# software.
+
+import nilearn
+fsaverage = nilearn.datasets.fetch_surf_fsaverage5()
+
+#########################################################################
+# The projection function simply takes the fMRI data and the mesh.
+# Note that those correspond spatially, are they are bothin MNI space.
+from nilearn import surface
+texture = surface.vol_to_surf(fmri_img, fsaverage.pial_right)
+
+#########################################################################
+# Perform first level analysis
+# ----------------------------
+#
+# This involves computing the design matrix and fitting the model.
+# We start by specifying the timing of fMRI frames
+
+import numpy as np
+n_scans = texture.shape[1]
+frame_times = t_r * (np.arange(n_scans) + .5)
+
+#########################################################################
+# Create the design matrix
+#
+# We specify an hrf model containing Glover model and its time derivative
+# the drift model is implicitly a cosine basis with period cutoff 128s.
+from nistats.design_matrix import make_first_level_design_matrix
+design_matrix = make_first_level_design_matrix(frame_times,
+ events=events,
+ hrf_model='glover + derivative'
+ )
+
+#########################################################################
+# Setup and fit GLM.
+# Note that the output consists in 2 variables: `labels` and `fit`
+# `labels` tags voxels according to noise autocorrelation.
+# `estimates` contains the parameter estimates.
+# We keep them for later contrast computation.
+
+from nistats.first_level_model import run_glm
+labels, estimates = run_glm(texture.T, design_matrix.values)
+
+#########################################################################
+# Estimate contrasts
+# ------------------
+# Specify the contrasts
+# For practical purpose, we first generate an identity matrix whose size is
+# the number of columns of the design matrix
+contrast_matrix = np.eye(design_matrix.shape[1])
+
+#########################################################################
+# first create basic contrasts
+basic_contrasts = dict([(column, contrast_matrix[i])
+ for i, column in enumerate(design_matrix.columns)])
+
+#########################################################################
+# add some intermediate contrasts
+basic_contrasts["audio"] = (basic_contrasts["clicDaudio"]
+ + basic_contrasts["clicGaudio"]
+ + basic_contrasts["calculaudio"]
+ + basic_contrasts["phraseaudio"]
+ )
+basic_contrasts["video"] = (basic_contrasts["clicDvideo"]
+ + basic_contrasts["clicGvideo"]
+ + basic_contrasts["calculvideo"]
+ + basic_contrasts["phrasevideo"]
+ )
+basic_contrasts["computation"] = (basic_contrasts["calculaudio"]
+ + basic_contrasts["calculvideo"]
+ )
+basic_contrasts["sentences"] = (basic_contrasts["phraseaudio"]
+ + basic_contrasts["phrasevideo"]
+ )
+
+#########################################################################
+# Finally make a dictionary of more relevant contrasts
+#
+# * "left - right button press" probes motor activity in left versus right button presses
+# * "audio - video" probes the difference of activity between listening to some content or reading the same type of content (instructions, stories)
+# * "computation - sentences" looks at the activity when performing a mental comptation task versus simply reading sentences.
+#
+# Of course, we could define other contrasts, but we keep only 3 for simplicity.
+
+contrasts = {
+ "left - right button press": (basic_contrasts["clicGaudio"]
+ + basic_contrasts["clicGvideo"]
+ - basic_contrasts["clicDaudio"]
+ - basic_contrasts["clicDvideo"]
+ ),
+ "audio - video": basic_contrasts["audio"] - basic_contrasts["video"],
+ "computation - sentences": (basic_contrasts["computation"] -
+ basic_contrasts["sentences"]
+ )
+ }
+
+#########################################################################
+# contrast estimation
+from nistats.contrasts import compute_contrast
+from nilearn import plotting
+
+#########################################################################
+# iterate over contrasts
+for index, (contrast_id, contrast_val) in enumerate(contrasts.items()):
+ print(' Contrast % i out of %i: %s, right hemisphere' %
+ (index + 1, len(contrasts), contrast_id))
+ # compute contrast-related statistics
+ contrast = compute_contrast(labels, estimates, contrast_val,
+ contrast_type='t')
+ # we present the Z-transform of the t map
+ z_score = contrast.z_score()
+ # we plot it on the surface, on the inflated fsaverage mesh,
+ # together with a suitable background to give an impression
+ # of the cortex folding.
+ plotting.plot_surf_stat_map(
+ fsaverage.infl_right, z_score, hemi='right',
+ title=contrast_id, colorbar=True,
+ threshold=3., bg_map=fsaverage.sulc_right)
+
+#########################################################################
+# Analysing the left hemisphere
+# -----------------------------
+#
+# Note that it requires little additional code!
+
+#########################################################################
+# Project the fMRI data to the mesh
+texture = surface.vol_to_surf(fmri_img, fsaverage.pial_left)
+
+#########################################################################
+# Estimate the General Linear Model
+labels, estimates = run_glm(texture.T, design_matrix.values)
+
+#########################################################################
+# Create contrast-specific maps
+for index, (contrast_id, contrast_val) in enumerate(contrasts.items()):
+ print(' Contrast % i out of %i: %s, left hemisphere' %
+ (index + 1, len(contrasts), contrast_id))
+ # compute contrasts
+ contrast = compute_contrast(labels, estimates, contrast_val,
+ contrast_type='t')
+ z_score = contrast.z_score()
+ # Plot the result
+ plotting.plot_surf_stat_map(
+ fsaverage.infl_left, z_score, hemi='left',
+ title=contrast_id, colorbar=True,
+ threshold=3., bg_map=fsaverage.sulc_left)
+
+plotting.show()
diff --git a/examples/02_first_level_models/plot_spm_multimodal_faces.py b/examples/02_first_level_models/plot_spm_multimodal_faces.py
index 22e94f44..8b9de1b2 100644
--- a/examples/02_first_level_models/plot_spm_multimodal_faces.py
+++ b/examples/02_first_level_models/plot_spm_multimodal_faces.py
@@ -1,6 +1,13 @@
-"""
-Minimal script for preprocessing single-subject data (two session)
-==================================================================
+"""Single-subject data (two sessions) in native space
+==================================================
+
+The example shows the analysis of an SPM dataset studying face perception.
+The anaylsis is performed in native spave. Realignment parameters are provided
+with the input images, but those have not been resampled to a common space.
+
+The experimental paradigm is simple, with two conditions; viewing a
+face image or a scrambled face image, supposedly with the same
+low-level statistical properties, to find face-specific responses.
For details on the data, please see:
Henson, R.N., Goshen-Gottstein, Y., Ganel, T., Otten, L.J., Quayle, A.,
@@ -8,37 +15,32 @@
perception, recognition and priming. Cereb Cortex. 2003 Jul;13(7):793-805.
http://www.dx.doi.org/10.1093/cercor/13.7.793
-Note: this example takes a lot of time because the input are lists of 3D images
+This example takes a lot of time because the input are lists of 3D images
sampled in different position (encoded by different) affine functions.
+
"""
print(__doc__)
-# Standard imports
-import numpy as np
-from scipy.io import loadmat
-import pandas as pd
-
-# Imports for GLM
-from nilearn.image import concat_imgs, resample_img, mean_img
-from nistats.design_matrix import make_design_matrix
-from nistats.first_level_model import FirstLevelModel
-from nistats.datasets import fetch_spm_multimodal_fmri
#########################################################################
# Fetch spm multimodal_faces data
+from nistats.datasets import fetch_spm_multimodal_fmri
subject_data = fetch_spm_multimodal_fmri()
#########################################################################
-# Experimental paradigm specification
-tr = 2.
-slice_time_ref = 0.
-drift_model = 'Cosine'
-hrf_model = 'spm + derivative'
-period_cut = 128.
+# Timing and design matrix parameter specification
+tr = 2. # repetition time, in seconds
+slice_time_ref = 0. # we will sample the design matrix at the beggining of each acquisition
+drift_model = 'Cosine' # We use a discrete cosin transform to model signal drifts.
+period_cut = 128. # The cutoff for the drift model is 1/128 Hz.
+hrf_model = 'spm + derivative' # The hemodunamic response finction is the SPM canonical one
#########################################################################
-# Resample the images
+# Resample the images.
+#
+# This is achieved by the concat_imgs function of Nilearn.
+from nilearn.image import concat_imgs, resample_img, mean_img
fmri_img = [concat_imgs(subject_data.func1, auto_resample=True),
concat_imgs(subject_data.func2, auto_resample=True)]
affine, shape = fmri_img[0].affine, fmri_img[0].shape
@@ -51,59 +53,84 @@
#########################################################################
# Make design matrices
+import numpy as np
+import pandas as pd
+from nistats.design_matrix import make_first_level_design_matrix
design_matrices = []
-for idx in range(len(fmri_img)):
- # Build paradigm
- n_scans = fmri_img[idx].shape[-1]
- timing = loadmat(getattr(subject_data, "trials_ses%i" % (idx + 1)),
- squeeze_me=True, struct_as_record=False)
-
- faces_onsets = timing['onsets'][0].ravel()
- scrambled_onsets = timing['onsets'][1].ravel()
- onsets = np.hstack((faces_onsets, scrambled_onsets))
- onsets *= tr # because onsets were reporting in 'scans' units
- conditions = (['faces'] * len(faces_onsets) +
- ['scrambled'] * len(scrambled_onsets))
- paradigm = pd.DataFrame({'trial_type': conditions, 'onset': onsets})
-
- # Build design matrix
+
+#########################################################################
+# loop over the two sessions
+for idx, img in enumerate(fmri_img, start=1):
+ # Build experimental paradigm
+ n_scans = img.shape[-1]
+ events = pd.read_table(subject_data['events{}'.format(idx)])
+ # Define the sampling times for the design matrix
frame_times = np.arange(n_scans) * tr
- design_matrix = make_design_matrix(
- frame_times, paradigm, hrf_model=hrf_model, drift_model=drift_model,
- period_cut=period_cut)
+ # Build design matrix with the reviously defined parameters
+ design_matrix = make_first_level_design_matrix(
+ frame_times,
+ events,
+ hrf_model=hrf_model,
+ drift_model=drift_model,
+ period_cut=period_cut,
+ )
+
+ # put the design matrices in a list
design_matrices.append(design_matrix)
#########################################################################
-# We can specify basic contrasts (To get beta maps)
+# We can specify basic contrasts (to get beta maps).
+# We start by specifying canonical contrast that isolate design matrix columns
contrast_matrix = np.eye(design_matrix.shape[1])
-contrasts = dict([(column, contrast_matrix[i])
+basic_contrasts = dict([(column, contrast_matrix[i])
for i, column in enumerate(design_matrix.columns)])
+
#########################################################################
-# Instead in this example we define more interesting contrasts
+# We actually want more interesting contrasts. The simplest contrast
+# just makes the difference between the two main conditions. We
+# define the two opposite versions to run one-tail t-tests. We also
+# define the effects of interest contrast, a 2-dimensional contrasts
+# spanning the two conditions.
+
contrasts = {
- 'faces-scrambled': contrasts['faces'] - contrasts['scrambled'],
- 'scrambled-faces': -contrasts['faces'] + contrasts['scrambled'],
- 'effects_of_interest': np.vstack((contrasts['faces'],
- contrasts['scrambled']))
+ 'faces-scrambled': basic_contrasts['faces'] - basic_contrasts['scrambled'],
+ 'scrambled-faces': -basic_contrasts['faces'] + basic_contrasts['scrambled'],
+ 'effects_of_interest': np.vstack((basic_contrasts['faces'],
+ basic_contrasts['scrambled']))
}
#########################################################################
-# Fit GLM
+# Fit the GLM -- 2 sessions.
+# Imports for GLM, the sepcify, then fit.
+from nistats.first_level_model import FirstLevelModel
print('Fitting a GLM')
-fmri_glm = FirstLevelModel(tr, slice_time_ref)
+fmri_glm = FirstLevelModel()
fmri_glm = fmri_glm.fit(fmri_img, design_matrices=design_matrices)
#########################################################################
-# Compute contrast maps
+# Compute contrast-related statistical maps (in z-scale), and plot them
print('Computing contrasts')
from nilearn import plotting
+# Iterate on contrasts
for contrast_id, contrast_val in contrasts.items():
print("\tcontrast id: %s" % contrast_id)
+ # compute the contrasts
z_map = fmri_glm.compute_contrast(
contrast_val, output_type='z_score')
+ # plot the contrasts as soon as they're generated
+ # the display is overlayed on the mean fMRI image
+ # a threshold of 3.0 is used. More sophisticated choices are possible.
plotting.plot_stat_map(
z_map, bg_img=mean_image, threshold=3.0, display_mode='z',
cut_coords=3, black_bg=True, title=contrast_id)
+#########################################################################
+# Show the resulting maps: We observe that the analysis results in
+# wide activity for the 'effects of interest' contrast, showing the
+# implications of large portions of the visual cortex in the
+# conditions. By contrast, the differential effect between "faces" and
+# "scambled" involves sparser, more anterior and lateral regions. It
+# displays also some responses in the frontal lobe.
+
plotting.show()
diff --git a/examples/03_second_level_models/plot_oasis.py b/examples/03_second_level_models/plot_oasis.py
new file mode 100644
index 00000000..15c2b5ef
--- /dev/null
+++ b/examples/03_second_level_models/plot_oasis.py
@@ -0,0 +1,127 @@
+"""Voxel-Based Morphometry on Oasis dataset
+========================================
+
+This example uses Voxel-Based Morphometry (VBM) to study the relationship
+between aging, sex and gray matter density.
+
+The data come from the `OASIS `_ project.
+If you use it, you need to agree with the data usage agreement available
+on the website.
+
+It has been run through a standard VBM pipeline (using SPM8 and
+NewSegment) to create VBM maps, which we study here.
+
+VBM analysis of aging
+---------------------
+
+We run a standard GLM analysis to study the association between age
+and gray matter density from the VBM data. We use only 100 subjects
+from the OASIS dataset to limit the memory usage.
+
+Note that more power would be obtained from using a larger sample of subjects.
+"""
+# Authors: Bertrand Thirion, , July 2018
+# Elvis Dhomatob, , Apr. 2014
+# Virgile Fritsch, , Apr 2014
+# Gael Varoquaux, Apr 2014
+
+
+n_subjects = 100 # more subjects requires more memory
+
+############################################################################
+# Load Oasis dataset
+# ------------------
+
+from nilearn import datasets
+oasis_dataset = datasets.fetch_oasis_vbm(n_subjects=n_subjects)
+gray_matter_map_filenames = oasis_dataset.gray_matter_maps
+age = oasis_dataset.ext_vars['age'].astype(float)
+
+###############################################################################
+# Sex is encoded as 'M' or 'F'. make it a binary variable
+sex = oasis_dataset.ext_vars['mf'] == b'F'
+
+###############################################################################
+# Print basic information on the dataset
+print('First gray-matter anatomy image (3D) is located at: %s' %
+ oasis_dataset.gray_matter_maps[0]) # 3D data
+print('First white-matter anatomy image (3D) is located at: %s' %
+ oasis_dataset.white_matter_maps[0]) # 3D data
+
+###############################################################################
+# Get a mask image: A mask of the cortex of the ICBM template
+gm_mask = datasets.fetch_icbm152_brain_gm_mask()
+
+###############################################################################
+# Resample the images, since this mask has a different resolution
+from nilearn.image import resample_to_img
+mask_img = resample_to_img(
+ gm_mask, gray_matter_map_filenames[0], interpolation='nearest')
+
+#############################################################################
+# Analyse data
+# ------------
+#
+# First create an adequate design matrix with three columns: 'age',
+# 'sex', 'intercept'.
+import pandas as pd
+import numpy as np
+intercept = np.ones(n_subjects)
+design_matrix = pd.DataFrame(np.vstack((age, sex, intercept)).T,
+ columns=['age', 'sex', 'intercept'])
+
+#############################################################################
+# Plot the design matrix
+from nistats.reporting import plot_design_matrix
+ax = plot_design_matrix(design_matrix)
+ax.set_title('Second level design matrix', fontsize=12)
+ax.set_ylabel('maps')
+
+##########################################################################
+# Specify and fit the second-level model when loading the data, we
+# smooth a little bit to improve statistical behavior
+
+from nistats.second_level_model import SecondLevelModel
+second_level_model = SecondLevelModel(smoothing_fwhm=2.0, mask=mask_img)
+second_level_model.fit(gray_matter_map_filenames,
+ design_matrix=design_matrix)
+
+##########################################################################
+# Estimate 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.
+# First compute the threshold.
+from nistats.thresholding import map_threshold
+_, threshold = map_threshold(
+ z_map, level=.05, height_control='fdr')
+print('The FDR=.05-corrected threshold is: %.3g' % threshold)
+
+###########################################################################
+# The plot it
+from nilearn import plotting
+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()
+
+###########################################################################
+# Can also study the effect of sex: compute the stat, compute the
+# threshold, plot the map
+
+z_map = second_level_model.compute_contrast(second_level_contrast='sex',
+ output_type='z_score')
+_, threshold = map_threshold(
+ z_map, level=.05, height_control='fdr')
+plotting.plot_stat_map(
+ z_map, threshold=threshold, colorbar=True,
+ title='sex effect on grey matter density (FDR = .05)')
+
+###########################################################################
+# Note that there does not seem to be any significant effect of sex on
+# grey matter density on that dataset.
+
diff --git a/examples/03_second_level_models/plot_second_level_association_test.py b/examples/03_second_level_models/plot_second_level_association_test.py
new file mode 100644
index 00000000..d862b437
--- /dev/null
+++ b/examples/03_second_level_models/plot_second_level_association_test.py
@@ -0,0 +1,78 @@
+"""
+Example of generic design in second-level models
+================================================
+
+This example shows the results obtained in a group analysis using a more
+complex contrast than a one- or two-sample t test.
+We use the [left button press (auditory cue)] task from the Localizer
+dataset and seek association between the contrast values and a variate
+that measures the speed of pseudo-word reading. No confounding variate
+is included in the model.
+
+
+"""
+# Author: Virgile Fritsch, Bertrand Thirion, 2014 -- 2018
+
+
+##############################################################################
+# Load Localizer contrast
+from nilearn import datasets
+n_samples = 94
+localizer_dataset = datasets.fetch_localizer_contrasts(
+ ['left button press (auditory cue)'], n_subjects=n_samples)
+
+##############################################################################
+# print basic information on the dataset
+print('First contrast nifti image (3D) is located at: %s' %
+ localizer_dataset.cmaps[0])
+
+##############################################################################
+# Load the behavioral variable
+tested_var = localizer_dataset.ext_vars['pseudo']
+print(tested_var)
+
+##############################################################################
+# Quality check / Remove subjects with bad tested variate
+import numpy as np
+mask_quality_check = np.where(tested_var != b'None')[0]
+n_samples = mask_quality_check.size
+contrast_map_filenames = [localizer_dataset.cmaps[i]
+ for i in mask_quality_check]
+tested_var = tested_var[mask_quality_check].astype(float).reshape((-1, 1))
+print("Actual number of subjects after quality check: %d" % n_samples)
+
+############################################################################
+# Estimate second level model
+# ---------------------------
+# We define the input maps and the design matrix for the second level model
+# and fit it.
+import pandas as pd
+design_matrix = pd.DataFrame(
+ np.hstack((tested_var, np.ones_like(tested_var))),
+ columns=['fluency', 'intercept'])
+
+###########################################################################
+# Fit of the second-level model
+from nistats.second_level_model import SecondLevelModel
+model = SecondLevelModel(smoothing_fwhm=5.0)
+model.fit(contrast_map_filenames, design_matrix=design_matrix)
+
+##########################################################################
+# To estimate the contrast is very simple. We can just provide the column
+# name of the design matrix.
+z_map = model.compute_contrast('fluency', output_type='z_score')
+
+###########################################################################
+# We compute the fdr-corrected p = 0.05 threshold for these data
+from nistats.thresholding import map_threshold
+_, threshold = map_threshold(z_map, level=.05, height_control='fdr')
+
+###########################################################################
+#Let us plot the second level contrast at the computed thresholds
+from nilearn import plotting
+plotting.plot_stat_map(
+ z_map, threshold=threshold, colorbar=True,
+ title='Group-level association between motor activity \n'
+ 'and reading fluency (fdr<0.05')
+
+plotting.show()
diff --git a/examples/03_second_level_models/plot_second_level_button_press.py b/examples/03_second_level_models/plot_second_level_one_sample_test.py
similarity index 73%
rename from examples/03_second_level_models/plot_second_level_button_press.py
rename to examples/03_second_level_models/plot_second_level_one_sample_test.py
index d4cf61d5..30b29d41 100644
--- a/examples/03_second_level_models/plot_second_level_button_press.py
+++ b/examples/03_second_level_models/plot_second_level_one_sample_test.py
@@ -1,9 +1,9 @@
"""
-GLM fitting in second level fMRI
-================================
+Second-level fMRI model: one sample test
+========================================
-Full step-by-step example of fitting a GLM to perform a second level analysis
-in experimental data and visualizing the results.
+Full step-by-step example of fitting a GLM to perform a second-level analysis (one-sample test)
+and visualizing the results.
More specifically:
@@ -11,23 +11,16 @@
2. a mask of the useful brain volume is computed
3. A one-sample t-test is applied to the brain maps
-(as fixed effects, then contrast estimation)
+We focus on a given contrast of the localizer dataset: the motor response to left versus right button press. Both at the ndividual and group level, this is expected to elicit activity in the motor cortex (positive in the right hemisphere, negative in the left hemisphere).
"""
-import pandas as pd
-from nilearn import plotting
-from scipy.stats import norm
-import matplotlib.pyplot as plt
-
-from nilearn.datasets import fetch_localizer_contrasts
-from nistats.second_level_model import SecondLevelModel
-
#########################################################################
# Fetch dataset
# --------------
# We download a list of left vs right button press contrasts from a
-# localizer dataset.
+# localizer dataset. Note that we fetc individual t-maps that represent the Bold activity estimate divided by the uncertainty about this estimate.
+from nilearn.datasets import fetch_localizer_contrasts
n_subjects = 16
data = fetch_localizer_contrasts(["left vs right button press"], n_subjects,
get_tmaps=True)
@@ -38,6 +31,8 @@
# We plot a grid with all the subjects t-maps thresholded at t = 2 for
# simple visualization purposes. The button press effect is visible among
# all subjects
+from nilearn import plotting
+import matplotlib.pyplot as plt
subjects = [subject_data[0] for subject_data in data['ext_vars']]
fig, axes = plt.subplots(nrows=4, ncols=4)
for cidx, tmap in enumerate(data['tmaps']):
@@ -53,10 +48,14 @@
# ---------------------------
# We define the input maps and the design matrix for the second level model
# and fit it.
+import pandas as pd
second_level_input = data['cmaps']
design_matrix = pd.DataFrame([1] * len(second_level_input),
columns=['intercept'])
+############################################################################
+# Model specification and fit
+from nistats.second_level_model import SecondLevelModel
second_level_model = SecondLevelModel(smoothing_fwhm=8.0)
second_level_model = second_level_model.fit(second_level_input,
design_matrix=design_matrix)
@@ -68,10 +67,13 @@
###########################################################################
# We threshold the second level contrast at uncorrected p < 0.001 and plot
+from scipy.stats import norm
p_val = 0.001
-z_th = norm.isf(p_val)
+p001_unc = norm.isf(p_val)
display = plotting.plot_glass_brain(
- z_map, threshold=z_th, colorbar=True, plot_abs=False, display_mode='z',
+ z_map, threshold=p001_unc, colorbar=True, display_mode='z', plot_abs=False,
title='group left-right button press (unc p<0.001')
+###########################################################################
+# As expected, we find the motor cortex
plotting.show()
diff --git a/examples/03_second_level_models/plot_second_level_two_sample_test.py b/examples/03_second_level_models/plot_second_level_two_sample_test.py
new file mode 100644
index 00000000..7b6ccc22
--- /dev/null
+++ b/examples/03_second_level_models/plot_second_level_two_sample_test.py
@@ -0,0 +1,89 @@
+"""Second-level fMRI model: a two-sample test
+==========================================
+
+Full step-by-step example of fitting a GLM to perform a second level analysis
+in experimental data and visualizing the results.
+
+More specifically:
+
+1. A sample of n=16 visual activity fMRIs are downloaded.
+2. A two-sample t-test is applied to the brain maps in order to see the effect of the contrast difference across subjects.
+
+The contrast is between responses to vertical versus horizontal
+checkerboards than are retinotopically distinct. At the individual
+level, these stimuli are sometimes used to map the borders of primary
+visual areas. At the group level, such a mapping is not possible. Yet,
+we may observe some significant effects in these areas.
+
+"""
+
+import pandas as pd
+from nilearn import plotting
+from nilearn.datasets import fetch_localizer_contrasts
+
+#########################################################################
+# Fetch dataset
+# --------------
+# We download a list of left vs right button press contrasts from a
+# localizer dataset.
+n_subjects = 16
+sample_vertical = fetch_localizer_contrasts(
+ ["vertical checkerboard"], n_subjects, get_tmaps=True)
+sample_horizontal = fetch_localizer_contrasts(
+ ["horizontal checkerboard"], n_subjects, get_tmaps=True)
+
+# What remains implicit here is that there is a one-to-one
+# correspondence between the two samples: the first image of both
+# samples comes from subject S1, the second from subject S2 etc.
+
+############################################################################
+# Estimate second level model
+# ---------------------------
+# We define the input maps and the design matrix for the second level model
+# and fit it.
+second_level_input = sample_vertical['cmaps'] + sample_horizontal['cmaps']
+
+############################################################################
+# model the effect of conditions (sample 1 vs sample 2)
+import numpy as np
+condition_effect = np.hstack(([1] * n_subjects, [- 1] * n_subjects))
+
+############################################################################
+# model the subject effect: each subject is observed in sample 1 and sample 2
+subject_effect = np.vstack((np.eye(n_subjects), np.eye(n_subjects)))
+subjects = ['S%02d' % i for i in range(1, n_subjects + 1)]
+
+############################################################################
+# Assemble those in a design matrix
+design_matrix = pd.DataFrame(
+ np.hstack((condition_effect[:, np.newaxis], subject_effect)),
+ columns=['vertical vs horizontal'] + subjects)
+
+############################################################################
+# plot the design_matrix
+from nistats.reporting import plot_design_matrix
+plot_design_matrix(design_matrix)
+
+############################################################################
+# formally specify the analysis model and fit it
+from nistats.second_level_model import SecondLevelModel
+second_level_model = SecondLevelModel().fit(
+ second_level_input, design_matrix=design_matrix)
+
+##########################################################################
+# Estimating the contrast is very simple. We can just provide the column
+# name of the design matrix.
+z_map = second_level_model.compute_contrast('vertical vs horizontal',
+ output_type='z_score')
+
+###########################################################################
+# We threshold the second level contrast and plot it
+threshold = 3.1 # correponds to p < .001, uncorrected
+display = plotting.plot_glass_brain(
+ z_map, threshold=threshold, colorbar=True, plot_abs=False,
+ title='vertical vs horizontal checkerboard (unc p<0.001')
+
+###########################################################################
+# Unsurprisingly, we see activity in the primary visual cortex, both positive and negative.
+
+plotting.show()
diff --git a/examples/03_second_level_models/plot_thresholding.py b/examples/03_second_level_models/plot_thresholding.py
index a956798d..5ee3446f 100644
--- a/examples/03_second_level_models/plot_thresholding.py
+++ b/examples/03_second_level_models/plot_thresholding.py
@@ -1,41 +1,45 @@
"""
-Example of simple second level analysis
-=======================================
+Statistical testing of a second-level analysis
+==============================================
-Perform a one-sample t-test on a bunch of images
-(a.k.a. second-level analyis in fMRI) and threshold a statistical image.
-This is based on the so-called localizer dataset.
+Perform a one-sample t-test on a bunch of images (a.k.a. second-level analyis in fMRI) and threshold the resulting statistical map.
+
+This example is based on the so-called localizer dataset.
It shows activation related to a mental computation task,
as opposed to narrative sentence reading/listening.
"""
-from nilearn import datasets
-from nilearn.input_data import NiftiMasker
#########################################################################
# Prepare some images for a simple t test
# ----------------------------------------
# This is a simple manually performed second level analysis
+from nilearn import datasets
n_samples = 20
localizer_dataset = datasets.fetch_localizer_calculation_task(
n_subjects=n_samples)
-# mask data
-nifti_masker = NiftiMasker(
- smoothing_fwhm=5,
- memory='nilearn_cache', memory_level=1) # cache options
+#########################################################################
+# Get the set of individual statstical maps (contrast estimates)
cmap_filenames = localizer_dataset.cmaps
#########################################################################
# Perform the second level analysis
# ----------------------------------
-# perform a one-sample test on these values
+#
+# First define a design matrix for the model. As the model is trivial (one-sample test), the design matrix is just one column with ones.
import pandas as pd
design_matrix = pd.DataFrame([1] * n_samples, columns=['intercept'])
+#########################################################################
+# Specify and estimate the model
from nistats.second_level_model import SecondLevelModel
second_level_model = SecondLevelModel().fit(
cmap_filenames, design_matrix=design_matrix)
+
+#########################################################################
+# Compute the only possible contrast: the one-sample test. Since there
+# is only one possible contrast, we don't need to specify it in detail
z_map = second_level_model.compute_contrast(output_type='z_score')
#########################################################################
@@ -43,22 +47,54 @@
# false positive rate < .001, cluster size > 10 voxels
from nistats.thresholding import map_threshold
thresholded_map1, threshold1 = map_threshold(
- z_map, threshold=.001, height_control='fpr', cluster_threshold=10)
+ z_map, level=.001, height_control='fpr', cluster_threshold=10)
#########################################################################
-# Now use FDR <.05, no cluster-level threshold
+# Now use FDR <.05, (False Discovery Rate) no cluster-level threshold
thresholded_map2, threshold2 = map_threshold(
- z_map, threshold=.05, height_control='fdr')
+ z_map, level=.05, height_control='fdr')
+print('The FDR=.05 threshold is %.3g' % threshold2)
+
+#########################################################################
+# Now use FWER <.05, (Familywise Error Rate) no cluster-level
+# threshold. As the data have not been intensively smoothed, we can
+# use a simple Bonferroni correction
+thresholded_map3, threshold3 = map_threshold(
+ z_map, level=.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)
+
+#########################################################################
+# These different thresholds correspond to different statistical
+# guarantees: in the FWER corrected image there is only a
+# probability<.05 of observing any false positive voxel. In the
+# FDR-corrected image, 5% of the voxels found are likely to be false
+# positive. In the uncorrected image, one expects a few tens of false
+# positive voxels.
+
plotting.show()
diff --git a/examples/04_low_level_functions/plot_design_matrix.py b/examples/04_low_level_functions/plot_design_matrix.py
index dfe811c5..b6db7373 100644
--- a/examples/04_low_level_functions/plot_design_matrix.py
+++ b/examples/04_low_level_functions/plot_design_matrix.py
@@ -9,65 +9,76 @@
Requires matplotlib
"""
-import numpy as np
try:
import matplotlib.pyplot as plt
except ImportError:
raise RuntimeError("This script needs the matplotlib library")
-from nistats.design_matrix import make_design_matrix
-from nistats.reporting import plot_design_matrix
-import pandas as pd
-
-
#########################################################################
# Define parameters
# ----------------------------------
# first we define parameters related to the images acquisition
-tr = 1.0
-n_scans = 128
-frame_times = np.arange(n_scans) * tr
+import numpy as np
+tr = 1.0 # repetition time is 1 second
+n_scans = 128 # the acquisition comprises 128 scans
+frame_times = np.arange(n_scans) * tr # here are the correspoding frame times
#########################################################################
# then we define parameters related to the experimental design
+# these are the types of the different trials
conditions = ['c0', 'c0', 'c0', 'c1', 'c1', 'c1', 'c3', 'c3', 'c3']
+duration = [1., 1., 1., 1., 1., 1., 1., 1., 1.]
+# these are the corresponding onset times
onsets = [30., 70., 100., 10., 30., 90., 30., 40., 60.]
-motion = np.cumsum(np.random.randn(n_scans, 6), 0) # simulate motion
+# Next, we simulate 6 motion parameters jointly observed with fMRI acquisitions
+motion = np.cumsum(np.random.randn(n_scans, 6), 0)
+# The 6 parameters correspond to three translations and three
+# rotations describing rigid body motion
add_reg_names = ['tx', 'ty', 'tz', 'rx', 'ry', 'rz']
#########################################################################
# Create design matrices
# -------------------------------------
# The same parameters allow us to obtain a variety of design matrices
-# First we compute an event-related design matrix
-paradigm = pd.DataFrame({'trial_type': conditions, 'onset': onsets})
+# We first create an events object
+import pandas as pd
+events = pd.DataFrame({'trial_type': conditions, 'onset': onsets,
+ 'duration': duration})
+#########################################################################
+# We sample the events into a design matrix, also including additional regressors
hrf_model = 'glover'
-X1 = make_design_matrix(
- frame_times, paradigm, drift_model='polynomial', drift_order=3,
+from nistats.design_matrix import make_first_level_design_matrix
+X1 = make_first_level_design_matrix(
+ frame_times, events, drift_model='polynomial', drift_order=3,
add_regs=motion, add_reg_names=add_reg_names, hrf_model=hrf_model)
#########################################################################
# Now we compute a block design matrix. We add duration to create the blocks.
+# For this we first define an event structure that includes the duration parameter
duration = 7. * np.ones(len(conditions))
-paradigm = pd.DataFrame({'trial_type': conditions, 'onset': onsets,
+events = pd.DataFrame({'trial_type': conditions, 'onset': onsets,
'duration': duration})
-X2 = make_design_matrix(frame_times, paradigm, drift_model='polynomial',
- drift_order=3, hrf_model=hrf_model)
+#########################################################################
+# Then we sample the design matrix
+X2 = make_first_level_design_matrix(frame_times, events, drift_model='polynomial',
+ drift_order=3, hrf_model=hrf_model)
#########################################################################
# Finally we compute a FIR model
-paradigm = pd.DataFrame({'trial_type': conditions, 'onset': onsets})
+events = pd.DataFrame({'trial_type': conditions, 'onset': onsets,
+ 'duration': duration})
hrf_model = 'FIR'
-X3 = make_design_matrix(frame_times, paradigm, hrf_model='fir',
- drift_model='polynomial', drift_order=3,
- fir_delays=np.arange(1, 6))
+X3 = make_first_level_design_matrix(frame_times, events, hrf_model='fir',
+ drift_model='polynomial', drift_order=3,
+ fir_delays=np.arange(1, 6))
#########################################################################
# Here the three designs side by side
+from nistats.reporting import plot_design_matrix
fig, (ax1, ax2, ax3) = plt.subplots(figsize=(10, 6), nrows=1, ncols=3)
plot_design_matrix(X1, ax=ax1)
ax1.set_title('Event-related design matrix', fontsize=12)
@@ -75,5 +86,8 @@
ax2.set_title('Block design matrix', fontsize=12)
plot_design_matrix(X3, ax=ax3)
ax3.set_title('FIR design matrix', fontsize=12)
+
+#########################################################################
+# Improve the layout and show the result
plt.subplots_adjust(left=0.08, top=0.9, bottom=0.21, right=0.96, wspace=0.3)
plt.show()
diff --git a/examples/04_low_level_functions/plot_hrf.py b/examples/04_low_level_functions/plot_hrf.py
index 59f6a2db..6b1a96a6 100644
--- a/examples/04_low_level_functions/plot_hrf.py
+++ b/examples/04_low_level_functions/plot_hrf.py
@@ -1,38 +1,56 @@
-"""
-Example of hemodynamic reponse functions.
+"""Example of hemodynamic reponse functions.
=========================================
-Plot the hrf model in SPM together with the hrf shape proposed by
-G.Glover, as well as their time and dispersion derivatives.
+Plot the hemodynamic reponse function (hrf) model in SPM together with
+the hrf shape proposed by G.Glover, as well as their time and
+dispersion derivatives.
-Requires matplotlib
+Requires matplotlib.
-"""
+The hrf is the filter that couples neural responses to the
+metabolic-related changes in the MRI signal. hrf models are simply
+phenomenological.
-import numpy as np
-import matplotlib.pyplot as plt
-from nistats import hemodynamic_models
+In current analysis frameworks, the choice of hrf model is essentially
+left to the user. Fortunately, using spm or Glover model does not make
+a huge difference. Adding derivatives should be considered whenever
+timing information has some degree of uncertainty. It is actually
+useful to detect timing issues.
+"""
#########################################################################
-# A first step: looking at our data
-# ----------------------------------
+# Set up some parameters for model display
+# ----------------------------------------
#
-# Let's quickly plot this file:
+# To get an impulse reponse, we simulate a single event occurring at time t=0, with duration 1s.
+import numpy as np
frame_times = np.linspace(0, 30, 61)
onset, amplitude, duration = 0., 1., 1.
+exp_condition = np.array((onset, duration, amplitude)).reshape(3, 1)
+
+#########################################################################
+# Sample this on a fris for display
stim = np.zeros_like(frame_times)
stim[(frame_times > onset) * (frame_times <= onset + duration)] = amplitude
-exp_condition = np.array((onset, duration, amplitude)).reshape(3, 1)
+
+#########################################################################
+# Define the candidate hrf models
hrf_models = [None, 'glover + derivative', 'glover + derivative + dispersion']
#########################################################################
# sample the hrf
+# --------------
+from nistats import hemodynamic_models
+import matplotlib.pyplot as plt
+
fig = plt.figure(figsize=(9, 4))
for i, hrf_model in enumerate(hrf_models):
+ # obtain the signal of interest by convolution
signal, name = hemodynamic_models.compute_regressor(
exp_condition, hrf_model, frame_times, con_id='main',
oversampling=16)
+ # plot this
plt.subplot(1, 3, i + 1)
plt.fill(frame_times, stim, 'k', alpha=.5, label='stimulus')
for j in range(signal.shape[1]):
@@ -41,5 +59,6 @@
plt.legend(loc=1)
plt.title(hrf_model)
+# adjust the plot
plt.subplots_adjust(bottom=.12)
plt.show()
diff --git a/examples/04_low_level_functions/plot_second_level_design_matrix.py b/examples/04_low_level_functions/plot_second_level_design_matrix.py
index b3205e75..38b00a84 100644
--- a/examples/04_low_level_functions/plot_second_level_design_matrix.py
+++ b/examples/04_low_level_functions/plot_second_level_design_matrix.py
@@ -1,7 +1,14 @@
-"""
-Example of second level design matrix
+"""Example of second level design matrix
=====================================
+The shows how a second-level design matrix is specified: assuming that
+the data refer to a group of individuals, with one image per subject,
+the design matrix typically holds the characteristics of each
+individual.
+
+This is used in a second-level analysis to assess the impact of these
+characteristics on brain signals.
+
Requires matplotlib.
"""
@@ -11,10 +18,6 @@
except ImportError:
raise RuntimeError("This script needs the matplotlib library")
-from nistats.design_matrix import create_second_level_design
-from nistats.reporting import plot_design_matrix
-import pandas as pd
-
#########################################################################
# Create a simple experimental paradigm
# --------------------------------------
@@ -25,18 +28,20 @@
##############################################################################
# Specify extra information about the subjects to create confounders
# Without confounders the design matrix would correspond to a one sample test
+import pandas as pd
extra_info_subjects = pd.DataFrame({'subject_label': subjects_label,
'age': range(15, 15 + n_subjects),
'sex': [0, 1] * int(n_subjects / 2)})
-
#########################################################################
# Create a second level design matrix
# -----------------------------------
+from nistats.design_matrix import make_second_level_design_matrix
+design_matrix = make_second_level_design_matrix(subjects_label, extra_info_subjects)
-design_matrix = create_second_level_design(subjects_label, extra_info_subjects)
-
+#########################################################################
# plot the results
+from nistats.reporting import plot_design_matrix
ax = plot_design_matrix(design_matrix)
ax.set_title('Second level design matrix', fontsize=12)
ax.set_ylabel('maps')
diff --git a/examples/04_low_level_functions/write_events_file.py b/examples/04_low_level_functions/write_events_file.py
new file mode 100644
index 00000000..65ba3248
--- /dev/null
+++ b/examples/04_low_level_functions/write_events_file.py
@@ -0,0 +1,65 @@
+"""Example of a events.tsv file generation: the neurospin/localizer events.
+=============================================================================
+
+The protocol described is the so-called "archi standard" localizer
+event sequence. See Pinel et al., BMC neuroscience 2007 for reference
+"""
+
+print(__doc__)
+
+#########################################################################
+# Define the onset times in seconds. Those are typically extracted
+# from the stimulation software used.
+import numpy as np
+onset = np.array([
+ 0., 2.4, 8.7, 11.4, 15., 18., 20.7, 23.7, 26.7, 29.7, 33., 35.4, 39.,
+ 41.7, 44.7, 48., 56.4, 59.7, 62.4, 69., 71.4, 75., 83.4, 87., 89.7,
+ 96., 108., 116.7, 119.4, 122.7, 125.4, 131.4, 135., 137.7, 140.4,
+ 143.4, 146.7, 149.4, 153., 156., 159., 162., 164.4, 167.7, 170.4,
+ 173.7, 176.7, 188.4, 191.7, 195., 198., 201., 203.7, 207., 210.,
+ 212.7, 215.7, 218.7, 221.4, 224.7, 227.7, 230.7, 234., 236.7, 246.,
+ 248.4, 251.7, 254.7, 257.4, 260.4, 264., 266.7, 269.7, 275.4, 278.4,
+ 284.4, 288., 291., 293.4, 296.7])
+
+#########################################################################
+# Associated trial types: these are numbered between 0 and 9, hence
+# correspond to 10 different conditions
+trial_idx = np.array(
+ [7, 7, 0, 2, 9, 4, 9, 3, 5, 9, 1, 6, 8, 8, 6, 6, 8, 0, 3, 4, 5, 8, 6,
+ 2, 9, 1, 6, 5, 9, 1, 7, 8, 6, 6, 1, 2, 9, 0, 7, 1, 8, 2, 7, 8, 3, 6,
+ 0, 0, 6, 8, 7, 7, 1, 1, 1, 5, 5, 0, 7, 0, 4, 2, 7, 9, 8, 0, 6, 3, 3,
+ 7, 1, 0, 0, 4, 1, 9, 8, 4, 9, 9])
+
+#########################################################################
+# We may want to map these indices to explicit condition names.
+# For that, we define a list of 10 strings.
+condition_ids = ['horizontal checkerboard',
+ 'vertical checkerboard',
+ 'right button press, auditory instructions',
+ 'left button press, auditory instructions',
+ 'right button press, visual instructions',
+ 'left button press, visual instructions',
+ 'mental computation, auditory instructions',
+ 'mental computation, visual instructions',
+ 'visual sentence',
+ 'auditory sentence']
+
+trial_type = np.array([condition_ids[i] for i in trial_idx])
+
+#########################################################################
+# We also define a duration (required by BIDS conventions)
+duration = np.ones_like(onset)
+
+
+#########################################################################
+# Form an event dataframe from these information
+import pandas as pd
+events = pd.DataFrame({'trial_type': trial_type,
+ 'onset': onset,
+ 'duration': duration})
+
+#########################################################################
+# Export them to a tsv file
+tsvfile = 'localizer_events.tsv'
+events.to_csv(tsvfile, sep='\t', index=False)
+print("Created the events file in %s " % tsvfile)
diff --git a/examples/04_low_level_functions/write_paradigm_file.py b/examples/04_low_level_functions/write_paradigm_file.py
deleted file mode 100644
index fda58570..00000000
--- a/examples/04_low_level_functions/write_paradigm_file.py
+++ /dev/null
@@ -1,42 +0,0 @@
-"""
-Example of a paradigm .csv file generation: the neurospin/localizer paradigm.
-=============================================================================
-
-See Pinel et al., BMC neuroscience 2007 for reference
-"""
-
-print(__doc__)
-
-import numpy as np
-import pandas as pd
-
-#########################################################################
-# onset times in milliseconds
-time = np.array([
- 0., 2.4, 8.7, 11.4, 15., 18., 20.7, 23.7, 26.7, 29.7, 33., 35.4, 39.,
- 41.7, 44.7, 48., 56.4, 59.7, 62.4, 69., 71.4, 75., 83.4, 87., 89.7,
- 96., 108., 116.7, 119.4, 122.7, 125.4, 131.4, 135., 137.7, 140.4,
- 143.4, 146.7, 149.4, 153., 156., 159., 162., 164.4, 167.7, 170.4,
- 173.7, 176.7, 188.4, 191.7, 195., 198., 201., 203.7, 207., 210.,
- 212.7, 215.7, 218.7, 221.4, 224.7, 227.7, 230.7, 234., 236.7, 246.,
- 248.4, 251.7, 254.7, 257.4, 260.4, 264., 266.7, 269.7, 275.4, 278.4,
- 284.4, 288., 291., 293.4, 296.7])
-
-#########################################################################
-# corresponding onset types
-trial_idx = np.array(
- [7, 7, 0, 2, 9, 4, 9, 3, 5, 9, 1, 6, 8, 8, 6, 6, 8, 0, 3, 4, 5, 8, 6,
- 2, 9, 1, 6, 5, 9, 1, 7, 8, 6, 6, 1, 2, 9, 0, 7, 1, 8, 2, 7, 8, 3, 6,
- 0, 0, 6, 8, 7, 7, 1, 1, 1, 5, 5, 0, 7, 0, 4, 2, 7, 9, 8, 0, 6, 3, 3,
- 7, 1, 0, 0, 4, 1, 9, 8, 4, 9, 9])
-
-condition_ids = ['damier_H', 'damier_V', 'clicDaudio', 'clicGaudio',
- 'clicDvideo', 'clicGvideo', 'calculaudio', 'calculvideo',
- 'phrasevideo', 'phraseaudio']
-
-trial_type = np.array([condition_ids[i] for i in trial_idx])
-events = pd.DataFrame({'trial_type': trial_type, 'onset': time})
-csvfile = 'localizer_paradigm.csv'
-events.to_csv(csvfile)
-
-print("Created the paradigm file in %s " % csvfile)
diff --git a/examples/05_complete_examples/README.txt b/examples/05_complete_examples/README.txt
new file mode 100644
index 00000000..398942b8
--- /dev/null
+++ b/examples/05_complete_examples/README.txt
@@ -0,0 +1,4 @@
+Complete examples
+-----------------
+
+Complete step-by-step examples, including First and Second Level Models.
diff --git a/examples/01_tutorials/plot_bids_analysis.py b/examples/05_complete_examples/plot_bids_analysis.py
similarity index 84%
rename from examples/01_tutorials/plot_bids_analysis.py
rename to examples/05_complete_examples/plot_bids_analysis.py
index 96ffe519..784fdf85 100644
--- a/examples/01_tutorials/plot_bids_analysis.py
+++ b/examples/05_complete_examples/plot_bids_analysis.py
@@ -15,8 +15,6 @@
in this case the preprocessed bold images were already normalized to the
same MNI space.
-
-
To run this example, you must launch IPython via ``ipython
--matplotlib`` in a terminal, or use the Jupyter notebook.
@@ -37,6 +35,10 @@
from nistats.datasets import fetch_bids_langloc_dataset
data_dir, _ = fetch_bids_langloc_dataset()
+##############################################################################
+# Here is the location of the dataset on disk
+print(data_dir)
+
##############################################################################
# Obtain automatically FirstLevelModel objects and fit arguments
# --------------------------------------------------------------
@@ -82,23 +84,32 @@
# ----------------------------
# Now we simply fit each first level model and plot for each subject the
# contrast that reveals the language network (language - string). Notice that
-# we can define a contrast using the names of the conditions especified in the
+# we can define a contrast using the names of the conditions specified in the
# events dataframe. Sum, substraction and scalar multiplication are allowed.
+
+############################################################################
+# Set the threshold as the z-variate with an uncorrected p-value of 0.001
+from scipy.stats import norm
+p001_unc = norm.isf(0.001)
+
+############################################################################
+# Prepare figure for concurrent plot of individual maps
from nilearn import plotting
import matplotlib.pyplot as plt
-from scipy.stats import norm
-fig, axes = plt.subplots(nrows=2, ncols=5)
+fig, axes = plt.subplots(nrows=2, ncols=5, figsize=(8, 4.5))
model_and_args = zip(models, models_run_imgs, models_events, models_confounds)
for midx, (model, imgs, events, confounds) in enumerate(model_and_args):
+ # fit the GLM
model.fit(imgs, events, confounds)
+ # compute the contrast of interest
zmap = model.compute_contrast('language-string')
- plotting.plot_glass_brain(zmap, colorbar=False, threshold=norm.isf(0.001),
+ plotting.plot_glass_brain(zmap, colorbar=False, threshold=p001_unc,
title=('sub-' + model.subject_label),
axes=axes[int(midx / 5), int(midx % 5)],
plot_abs=False, display_mode='x')
fig.suptitle('subjects z_map language network (unc p<0.001)')
-plt.show()
+plotting.show()
#########################################################################
# Second level model estimation
@@ -109,6 +120,9 @@
# column names)
from nistats.second_level_model import SecondLevelModel
second_level_input = models
+
+#########################################################################
+# Note that we apply a smoothing of 8mm.
second_level_model = SecondLevelModel(smoothing_fwhm=8.0)
second_level_model = second_level_model.fit(second_level_input)
@@ -117,12 +131,14 @@
# Since we are not providing confounders we are performing an one-sample test
# at the second level with the images determined by the specified first level
# contrast.
-zmap = second_level_model.compute_contrast(first_level_contrast='language-string')
+zmap = second_level_model.compute_contrast(
+ first_level_contrast='language-string')
#########################################################################
# The group level contrast reveals a left lateralized fronto-temporal
# language network
-plotting.plot_glass_brain(zmap, colorbar=True, threshold=norm.isf(0.001),
+plotting.plot_glass_brain(zmap, colorbar=True, threshold=p001_unc,
title='Group language network (unc p<0.001)',
plot_abs=False, display_mode='x')
plotting.show()
+
diff --git a/examples/05_complete_examples/plot_surface_bids_analysis.py b/examples/05_complete_examples/plot_surface_bids_analysis.py
new file mode 100644
index 00000000..9f5f2b82
--- /dev/null
+++ b/examples/05_complete_examples/plot_surface_bids_analysis.py
@@ -0,0 +1,161 @@
+"""Surface-based dataset first and second level analysis of a dataset
+==================================================================
+
+
+Full step-by-step example of fitting a GLM (first and second level
+analysis) in a 10-subjects dataset and visualizing the results.
+
+More specifically:
+
+1. Download an fMRI BIDS dataset with two language conditions to contrast.
+2. Project the data to a standard mesh, fsaverage5, aka the Freesurfer template mesh downsampled to about 10k nodes per hemisphere.
+3. Run the first level model objects
+4. Fit a second level model on the fitted first level models.
+
+Notice that in this case the preprocessed bold images were already
+ normalized to the same MNI space.
+
+To run this example, you must launch IPython via ``ipython
+--matplotlib`` in a terminal, or use the Jupyter notebook.
+
+.. contents:: **Contents**
+ :local:
+ :depth: 1
+
+"""
+
+##############################################################################
+# Fetch example BIDS dataset
+# --------------------------
+# We download an simplified BIDS dataset made available for illustrative
+# purposes. It contains only the necessary
+# information to run a statistical analysis using Nistats. The raw data
+# subject folders only contain bold.json and events.tsv files, while the
+# derivatives folder with preprocessed files contain preproc.nii and
+# confounds.tsv files.
+from nistats.datasets import fetch_bids_langloc_dataset
+data_dir, _ = fetch_bids_langloc_dataset()
+
+##############################################################################
+# Here is the location of the dataset on disk
+print(data_dir)
+
+##############################################################################
+# Obtain automatically FirstLevelModel objects and fit arguments
+# --------------------------------------------------------------
+# From the dataset directory we obtain automatically FirstLevelModel objects
+# with their subject_id filled from the BIDS dataset. Moreover we obtain
+# for each model a dictionary with run_imgs, events and confounder regressors
+# since in this case a confounds.tsv file is available in the BIDS dataset.
+# To get the first level models we only have to specify the dataset directory
+# and the task_label as specified in the file names.
+from nistats.first_level_model import first_level_models_from_bids
+task_label = 'languagelocalizer'
+space_label = 'MNI152nonlin2009aAsym'
+_, models_run_imgs, models_events, models_confounds = \
+ first_level_models_from_bids(
+ data_dir, task_label, space_label,
+ img_filters=[('variant', 'smoothResamp')])
+
+#############################################################################
+# We also need to get the TR information. For that we use a json file
+# of the dataset
+import os
+json_file = os.path.join(data_dir, 'sub-01', 'ses-02', 'func',
+ 'sub-01_ses-02_task-languagelocalizer_bold.json')
+import json
+with open(json_file, 'r') as f:
+ t_r = json.load(f)['RepetitionTime']
+
+#############################################################################
+# Project fMRI data to the surface: First get fsaverage5
+from nilearn.datasets import fetch_surf_fsaverage
+fsaverage = fetch_surf_fsaverage(mesh='fsaverage5')
+
+#########################################################################
+# The projection function simply takes the fMRI data and the mesh.
+# Note that those correspond spatially, are they are bothin MNI space.
+import numpy as np
+from nilearn import surface
+from nistats.design_matrix import make_first_level_design_matrix
+from nistats.first_level_model import run_glm
+from nistats.contrasts import compute_contrast
+
+#########################################################################
+# Empty lists in which we are going to store activation values.
+z_scores_right = []
+z_scores_left = []
+for (fmri_img, confound, events) in zip(
+ models_run_imgs, models_confounds, models_events):
+ texture = surface.vol_to_surf(fmri_img[0], fsaverage.pial_right)
+ n_scans = texture.shape[1]
+ frame_times = t_r * (np.arange(n_scans) + .5)
+
+ # Create the design matrix
+ #
+ # We specify an hrf model containing Glover model and its time derivative
+ # the drift model is implicitly a cosine basis with period cutoff 128s.
+ design_matrix = make_first_level_design_matrix(
+ frame_times, events=events[0], hrf_model='glover + derivative',
+ add_regs=confound[0])
+
+ # contrast_specification
+ contrast_values = (design_matrix.columns == 'language') * 1.0 -\
+ (design_matrix.columns == 'string')
+
+ # Setup and fit GLM.
+ # Note that the output consists in 2 variables: `labels` and `fit`
+ # `labels` tags voxels according to noise autocorrelation.
+ # `estimates` contains the parameter estimates.
+ # We input them for contrast computation.
+ labels, estimates = run_glm(texture.T, design_matrix.values)
+ contrast = compute_contrast(labels, estimates, contrast_values,
+ contrast_type='t')
+ # we present the Z-transform of the t map
+ z_score = contrast.z_score()
+ z_scores_right.append(z_score)
+
+ # Do the left hemipshere exactly in the same way
+ texture = surface.vol_to_surf(fmri_img, fsaverage.pial_left)
+ labels, estimates = run_glm(texture.T, design_matrix.values)
+ contrast = compute_contrast(labels, estimates, contrast_values,
+ contrast_type='t')
+ z_scores_left.append(contrast.z_score())
+
+############################################################################
+# Individual activation maps have been accumulated in the z_score_left
+# and az_scores_right lists respectively. We can now use them in a
+# group study (one -sample study)
+
+############################################################################
+# Group study
+# -----------
+#
+# Prepare figure for concurrent plot of individual maps
+# compute population-level maps for left and right hemisphere
+# we directetly do that on the values arrays
+from scipy.stats import ttest_1samp, norm
+t_left, pval_left = ttest_1samp(np.array(z_scores_left), 0)
+t_right, pval_right = ttest_1samp(np.array(z_scores_right), 0)
+
+############################################################################
+# What we have so far are p-values: we convert them to z-values for plotting
+z_val_left = norm.isf(pval_left)
+z_val_right = norm.isf(pval_right)
+
+############################################################################
+# Plot the resulting maps.
+# Left hemipshere
+from nilearn import plotting
+plotting.plot_surf_stat_map(
+ fsaverage.infl_left, z_val_left, hemi='left',
+ title="language-string, left hemisphere", colorbar=True,
+ threshold=3., bg_map=fsaverage.sulc_left)
+############################################################################
+# Right hemisphere
+plotting.plot_surf_stat_map(
+ fsaverage.infl_right, z_val_left, hemi='right',
+ title="language-string, right hemisphere", colorbar=True,
+ threshold=3., bg_map=fsaverage.sulc_right)
+
+plotting.show()
diff --git a/examples/README.txt b/examples/README.txt
index 03e887c8..54853efa 100644
--- a/examples/README.txt
+++ b/examples/README.txt
@@ -12,8 +12,3 @@ Nistats usage examples
.. contents:: **Contents**
:local:
:depth: 1
-
-General examples
-----------------
-
-General-purpose and introductory examples for nistats.
diff --git a/nistats/__init__.py b/nistats/__init__.py
index d22ad687..c584ed57 100644
--- a/nistats/__init__.py
+++ b/nistats/__init__.py
@@ -26,9 +26,24 @@
"""
import gzip
+import six
+import warnings
from .version import _check_module_dependencies, __version__
+
+def py2_deprecation_warning():
+ warnings.simplefilter('once')
+ py2_warning = ('Python2 support is deprecated and will be removed in '
+ 'a future release. Consider switching to Python3.')
+ if six.PY2:
+ warnings.warn(message=py2_warning,
+ category=DeprecationWarning,
+ stacklevel=4,
+ )
+
_check_module_dependencies()
+
__all__ = ['__version__', 'datasets', 'design_matrix']
+py2_deprecation_warning()
diff --git a/nistats/contrasts.py b/nistats/contrasts.py
index 7bc4f275..fb5f80cd 100644
--- a/nistats/contrasts.py
+++ b/nistats/contrasts.py
@@ -54,23 +54,31 @@ def compute_contrast(labels, regression_result, con_val, contrast_type=None):
'"{0}" is not a known contrast type. Allowed types are {1}'.
format(contrast_type, acceptable_contrast_types))
- effect_ = np.zeros((dim, labels.size))
- var_ = np.zeros((dim, dim, labels.size))
if contrast_type == 't':
+ effect_ = np.zeros((1, labels.size))
+ var_ = np.zeros(labels.size)
for label_ in regression_result:
label_mask = labels == label_
resl = regression_result[label_].Tcontrast(con_val)
effect_[:, label_mask] = resl.effect.T
- var_[:, :, label_mask] = (resl.sd ** 2).T
- else:
+ var_[label_mask] = (resl.sd ** 2).T
+ elif contrast_type == 'F':
+ from scipy.linalg import sqrtm
+ effect_ = np.zeros((dim, labels.size))
+ var_ = np.zeros(labels.size)
for label_ in regression_result:
label_mask = labels == label_
- resl = regression_result[label_].Fcontrast(con_val)
- effect_[:, label_mask] = resl.effect
- var_[:, :, label_mask] = resl.covariance
+ reg = regression_result[label_]
+ cbeta = np.atleast_2d(np.dot(con_val, reg.theta))
+ invcov = np.linalg.inv(np.atleast_2d(
+ reg.vcov(matrix=con_val, dispersion=1.0)))
+ wcbeta = np.dot(sqrtm(invcov), cbeta)
+ rss = reg.dispersion
+ effect_[:, label_mask] = wcbeta
+ var_[label_mask] = rss
dof_ = regression_result[label_].df_resid
- return Contrast(effect=effect_, variance=var_, dof=dof_,
+ return Contrast(effect=effect_, variance=var_, dim=dim, dof=dof_,
contrast_type=contrast_type)
@@ -105,7 +113,7 @@ class Contrast(object):
(high-dimensional F constrasts may lead to memory breakage).
"""
- def __init__(self, effect, variance, dof=DEF_DOFMAX, contrast_type='t',
+ def __init__(self, effect, variance, dim=None, dof=DEF_DOFMAX, contrast_type='t',
tiny=DEF_TINY, dofmax=DEF_DOFMAX):
"""
Parameters
@@ -113,29 +121,30 @@ def __init__(self, effect, variance, dof=DEF_DOFMAX, contrast_type='t',
effect : array of shape (contrast_dim, n_voxels)
the effects related to the contrast
- variance : array of shape (contrast_dim, contrast_dim, n_voxels)
+ variance : array of shape (n_voxels)
the associated variance estimate
+ dim: int or None,
+ the dimension of the contrast
+
dof : scalar
- the degrees of freedom of the resiudals
+ the degrees of freedom of the residuals
contrast_type: {'t', 'F'}
specification of the contrast type
"""
- if variance.ndim != 3:
- raise ValueError('Variance array should have 3 dimensions')
+ if variance.ndim != 1:
+ raise ValueError('Variance array should have 1 dimension')
if effect.ndim != 2:
raise ValueError('Effect array should have 2 dimensions')
- if variance.shape[0] != variance.shape[1]:
- raise ValueError('Inconsistent shape for the variance estimate')
- if ((variance.shape[1] != effect.shape[0]) or
- (variance.shape[2] != effect.shape[1])):
- raise ValueError('Effect and variance have inconsistent shape')
self.effect = effect
self.variance = variance
self.dof = float(dof)
- self.dim = effect.shape[0]
+ if dim is None:
+ self.dim = effect.shape[0]
+ else:
+ self.dim = dim
if self.dim > 1 and contrast_type is 't':
print('Automatically converted multi-dimensional t to F contrast')
contrast_type = 'F'
@@ -154,7 +163,7 @@ def effect_size(self):
def effect_variance(self):
"""Make access to summary statistics more straightforward when
computing contrasts"""
- return self.variance[0, 0, :]
+ return self.variance
def stat(self, baseline=0.0):
""" Return the decision statistic associated with the test of the
@@ -173,22 +182,13 @@ def stat(self, baseline=0.0):
self.baseline = baseline
# Case: one-dimensional contrast ==> t or t**2
- if self.dim == 1:
+ if self.contrast_type == 'F':
+ stat = np.sum((self.effect - baseline) ** 2, 0) / self.dim /\
+ np.maximum(self.variance, self.tiny)
+ elif self.contrast_type == 't':
# avoids division by zero
stat = (self.effect - baseline) / np.sqrt(
np.maximum(self.variance, self.tiny))
- if self.contrast_type == 'F':
- stat = stat ** 2
- # Case: F contrast
- elif self.contrast_type == 'F':
- # F = |t|^2/q , |t|^2 = e^t inv(v) e
- if self.effect.ndim == 1:
- self.effect = self.effect[np.newaxis]
- if self.variance.ndim == 1:
- self.variance = self.variance[np.newaxis, np.newaxis]
- stat = (multiple_mahalanobis(
- self.effect - baseline, self.variance) / self.dim)
- # Unknwon stat
else:
raise ValueError('Unknown statistic type')
self.stat_ = stat
@@ -252,11 +252,13 @@ def __add__(self, other):
if self.dim != other.dim:
raise ValueError(
'The two contrasts do not have compatible dimensions')
+ dof_ = self.dof + other.dof
+ if self.contrast_type == 'F':
+ warn('Running approximate fixed effects on F statistics.')
effect_ = self.effect + other.effect
variance_ = self.variance + other.variance
- dof_ = self.dof + other.dof
- return Contrast(effect=effect_, variance=variance_, dof=dof_,
- contrast_type=self.contrast_type)
+ return Contrast(effect=effect_, variance=variance_, dim=self.dim,
+ dof=dof_, contrast_type=self.contrast_type)
def __rmul__(self, scalar):
"""Multiplication of the contrast by a scalar"""
diff --git a/nistats/datasets.py b/nistats/datasets.py
index 4767cb64..721e37bb 100644
--- a/nistats/datasets.py
+++ b/nistats/datasets.py
@@ -3,18 +3,27 @@
Author: Gael Varoquaux
"""
-import os
import fnmatch
-import re
import glob
import json
+import os
+import re
+import warnings
+
import nibabel as nib
+import numpy as np
+import pandas as pd
+from nilearn.datasets.utils import (_fetch_file,
+ _fetch_files,
+ _get_dataset_dir,
+ _uncompress_file,
+ )
+from scipy.io import loadmat
+from scipy.io.matlab.miobase import MatReadError
from sklearn.datasets.base import Bunch
-from nilearn.datasets.utils import (
- _get_dataset_dir, _fetch_files, _fetch_file, _uncompress_file)
+from nistats.utils import _verify_events_file_uses_tab_separators
-from botocore.handlers import disable_signing
SPM_AUDITORY_DATA_FILES = ["fM00223/fM00223_%03i.img" % index
for index in range(4, 100)]
@@ -102,6 +111,7 @@ def fetch_openneuro_dataset_index(
urls: list of string
Sorted list of dataset directories
"""
+ from botocore.handlers import disable_signing
boto3 = _check_import_boto3("boto3")
data_prefix = '{}/{}/uncompressed'.format(
dataset_version.split('_')[0], dataset_version)
@@ -273,6 +283,31 @@ def fetch_openneuro_dataset(
return data_dir, sorted(downloaded)
+def _make_events_file_localizer_first_level(events_file):
+ """ Makes the first-level localizer fMRI dataset events file
+ BIDS compliant. Overwrites the original file.
+ Adds headers in first row.
+ Removes first column (spurious data).
+ Uses Tab character as value separator.
+
+ Parameters
+ ----------
+ events_file: string
+ path to the localizer_first_level dataset's events file.
+
+ Returns
+ -------
+ None
+ """
+ events = pd.read_csv(events_file, sep=' ', header=None, index_col=None,
+ names=['session', 'trial_type', 'onset'],
+ )
+ events.drop(labels='session', axis=1, inplace=True)
+ # duration is required in BIDS specification
+ events['duration'] = np.ones_like(events.onset)
+ events.to_csv(events_file, sep='\t', index=False)
+
+
def fetch_localizer_first_level(data_dir=None, verbose=1):
""" Download a first-level localizer fMRI dataset
@@ -284,15 +319,15 @@ def fetch_localizer_first_level(data_dir=None, verbose=1):
Returns
-------
data: sklearn.datasets.base.Bunch
- dictionary-like object, keys are:
+ dictionary-like object, with the keys:
epi_img: the input 4D image
- paradigm: a csv file describing the paardigm
+ events: a csv file describing the paardigm
"""
url = 'ftp://ftp.cea.fr/pub/dsv/madic/download/nipy'
dataset_name = "localizer_first_level"
files = dict(epi_img="s12069_swaloc1_corr.nii.gz",
- paradigm="localizer_paradigm.csv")
+ events="localizer_paradigm.csv")
# The options needed for _fetch_files
options = [(filename, os.path.join(url, filename), {})
for _, filename in sorted(files.items())]
@@ -303,10 +338,128 @@ def fetch_localizer_first_level(data_dir=None, verbose=1):
verbose=verbose)
params = dict(zip(sorted(files.keys()), sub_files))
-
+ try:
+ _verify_events_file_uses_tab_separators(params['events'])
+ except ValueError:
+ _make_events_file_localizer_first_level(events_file=
+ params['events']
+ )
+
return Bunch(**params)
+def _download_spm_auditory_data(data_dir, subject_dir, subject_id):
+ print("Data absent, downloading...")
+ url = ("http://www.fil.ion.ucl.ac.uk/spm/download/data/MoAEpilot/"
+ "MoAEpilot.zip")
+ archive_path = os.path.join(subject_dir, os.path.basename(url))
+ _fetch_file(url, subject_dir)
+ try:
+ _uncompress_file(archive_path)
+ except:
+ print("Archive corrupted, trying to download it again.")
+ return fetch_spm_auditory(data_dir=data_dir, data_name="",
+ subject_id=subject_id)
+
+
+def _prepare_downloaded_spm_auditory_data(subject_dir):
+ """ Uncompresses downloaded spm_auditory dataset and organizes
+ the data into apprpriate directories.
+
+ Parameters
+ ----------
+ subject_dir: string
+ Path to subject's data directory.
+
+ Returns
+ -------
+ _subject_data: skl.Bunch object
+ Scikit-Learn Bunch object containing data of a single subject
+ from the SPM Auditory dataset.
+
+ """
+ subject_data = {}
+ for file_name in SPM_AUDITORY_DATA_FILES:
+ file_path = os.path.join(subject_dir, file_name)
+ if os.path.exists(file_path):
+ subject_data[file_name] = file_path
+ else:
+ print("%s missing from filelist!" % file_name)
+ return None
+
+ _subject_data = {}
+ _subject_data["func"] = sorted(
+ [subject_data[x] for x in subject_data.keys()
+ if re.match("^fM00223_0\d\d\.img$", os.path.basename(x))])
+
+ # volumes for this dataset of shape (64, 64, 64, 1); let's fix this
+ for x in _subject_data["func"]:
+ vol = nib.load(x)
+ if len(vol.shape) == 4:
+ vol = nib.Nifti1Image(vol.get_data()[:, :, :, 0],
+ vol.affine)
+ nib.save(vol, x)
+
+ _subject_data["anat"] = [subject_data[x] for x in subject_data.keys()
+ if re.match("^sM00223_002\.img$",
+ os.path.basename(x))][0]
+
+ # ... same thing for anat
+ vol = nib.load(_subject_data["anat"])
+ if len(vol.shape) == 4:
+ vol = nib.Nifti1Image(vol.get_data()[:, :, :, 0],
+ vol.affine)
+ nib.save(vol, _subject_data["anat"])
+
+ return Bunch(**_subject_data)
+
+
+def _make_path_events_file_spm_auditory_data(spm_auditory_data):
+ """
+ Accepts data for spm_auditory dataset as Bunch
+ and constructs the filepath for its events descriptor file.
+ Parameters
+ ----------
+ spm_auditory_data: Bunch
+
+ Returns
+ -------
+ events_filepath: string
+ Full path to the events.tsv file for spm_auditory dataset.
+ """
+ events_file_location = os.path.dirname(spm_auditory_data['func'][0])
+ events_filename = os.path.basename(events_file_location) + '_events.tsv'
+ events_filepath = os.path.join(events_file_location, events_filename)
+ return events_filepath
+
+
+def _make_events_file_spm_auditory_data(events_filepath):
+ """
+ Accepts destination filepath including filename and
+ creates the events.tsv file for the spm_auditory dataset.
+
+ Parameters
+ ----------
+ events_filepath: string
+ The path where the events file will be created;
+
+ Returns
+ -------
+ None
+
+ """
+ tr = 7.
+ epoch_duration = 6 * tr # duration in seconds
+ conditions = ['rest', 'active'] * 8
+ n_blocks = len(conditions)
+ duration = epoch_duration * np.ones(n_blocks)
+ onset = np.linspace(0, (n_blocks - 1) * epoch_duration, n_blocks)
+ events = pd.DataFrame(
+ {'onset': onset, 'duration': duration, 'trial_type': conditions})
+ events.to_csv(events_filepath, sep='\t', index=False,
+ columns=['onset', 'duration', 'trial_type'])
+
+
def fetch_spm_auditory(data_dir=None, data_name='spm_auditory',
subject_id="sub001", verbose=1):
"""Function to fetch SPM auditory single-subject data.
@@ -334,69 +487,139 @@ def fetch_spm_auditory(data_dir=None, data_name='spm_auditory',
data_dir = _get_dataset_dir(data_name, data_dir=data_dir,
verbose=verbose)
subject_dir = os.path.join(data_dir, subject_id)
-
- def _glob_spm_auditory_data():
- """glob data from subject_dir.
-
- """
-
- if not os.path.exists(subject_dir):
+ if not os.path.exists(subject_dir):
+ _download_spm_auditory_data(data_dir, subject_dir, subject_id)
+ spm_auditory_data = _prepare_downloaded_spm_auditory_data(subject_dir)
+ try:
+ spm_auditory_data['events']
+ except KeyError:
+ events_filepath = _make_path_events_file_spm_auditory_data(
+ spm_auditory_data)
+ if not os.path.isfile(events_filepath):
+ _make_events_file_spm_auditory_data(events_filepath)
+ spm_auditory_data['events'] = events_filepath
+ return spm_auditory_data
+
+
+def _get_func_data_spm_multimodal(subject_dir, session, _subject_data):
+ session_func = sorted(glob.glob(
+ os.path.join(
+ subject_dir,
+ ("fMRI/Session%i/fMETHODS-000%i-*-01.img" % (
+ session, session + 4)
+ )
+ )
+ ))
+ if len(session_func) < 390:
+ print("Missing %i functional scans for session %i." % (
+ 390 - len(session_func), session))
+ return None
+
+ _subject_data['func%i' % (session)] = session_func
+ return _subject_data
+
+
+def _get_session_trials_spm_multimodal(subject_dir, session, _subject_data):
+ sess_trials = os.path.join(
+ subject_dir,
+ "fMRI/trials_ses%i.mat" % (session))
+ if not os.path.isfile(sess_trials):
+ print("Missing session file: %s" % sess_trials)
+ return None
+
+ _subject_data['trials_ses%i' % (session)] = sess_trials
+ return _subject_data
+
+
+def _get_anatomical_data_spm_multimodal(subject_dir, _subject_data):
+ anat = os.path.join(subject_dir, "sMRI/smri.img")
+ if not os.path.isfile(anat):
+ print("Missing structural image.")
+ return None
+
+ _subject_data["anat"] = anat
+ return _subject_data
+
+
+def _glob_spm_multimodal_fmri_data(subject_dir):
+ """glob data from subject_dir."""
+ _subject_data = {'slice_order': 'descending'}
+
+ for session in range(1, 3):
+ # glob func data for session
+ _subject_data = _get_func_data_spm_multimodal(subject_dir, session, _subject_data)
+ if not _subject_data:
return None
+ # glob trials .mat file
+ _subject_data = _get_session_trials_spm_multimodal(subject_dir, session, _subject_data)
+ if not _subject_data:
+ return None
+ try:
+ events = _make_events_file_spm_multimodal_fmri(_subject_data, session)
+ except MatReadError as mat_err:
+ warnings.warn('{}. An events.tsv file cannot be generated'.format(str(mat_err)))
+ else:
+ events_filepath = _make_events_filepath_spm_multimodal_fmri(_subject_data, session)
+ events.to_csv(events_filepath, sep='\t', index=False)
+ _subject_data['events{}'.format(session)] = events_filepath
+
+
+ # glob for anat data
+ _subject_data = _get_anatomical_data_spm_multimodal(subject_dir, _subject_data)
+ if not _subject_data:
+ return None
+
+ return Bunch(**_subject_data)
+
+
+def _download_data_spm_multimodal(data_dir, subject_dir, subject_id):
+ print("Data absent, downloading...")
+ urls = [
+ # fmri
+ ("http://www.fil.ion.ucl.ac.uk/spm/download/data/mmfaces/"
+ "multimodal_fmri.zip"),
- subject_data = {}
- for file_name in SPM_AUDITORY_DATA_FILES:
- file_path = os.path.join(subject_dir, file_name)
- if os.path.exists(file_path):
- subject_data[file_name] = file_path
- else:
- print("%s missing from filelist!" % file_name)
- return None
-
- _subject_data = {}
- _subject_data["func"] = sorted(
- [subject_data[x] for x in subject_data.keys()
- if re.match("^fM00223_0\d\d\.img$", os.path.basename(x))])
-
- # volumes for this dataset of shape (64, 64, 64, 1); let's fix this
- for x in _subject_data["func"]:
- vol = nib.load(x)
- if len(vol.shape) == 4:
- vol = nib.Nifti1Image(vol.get_data()[:, :, :, 0],
- vol.affine)
- nib.save(vol, x)
+ # structural
+ ("http://www.fil.ion.ucl.ac.uk/spm/download/data/mmfaces/"
+ "multimodal_smri.zip")
+ ]
- _subject_data["anat"] = [subject_data[x] for x in subject_data.keys()
- if re.match("^sM00223_002\.img$",
- os.path.basename(x))][0]
+ for url in urls:
+ archive_path = os.path.join(subject_dir, os.path.basename(url))
+ _fetch_file(url, subject_dir)
+ try:
+ _uncompress_file(archive_path)
+ except:
+ print("Archive corrupted, trying to download it again.")
+ return fetch_spm_multimodal_fmri(data_dir=data_dir,
+ data_name="",
+ subject_id=subject_id)
- # ... same thing for anat
- vol = nib.load(_subject_data["anat"])
- if len(vol.shape) == 4:
- vol = nib.Nifti1Image(vol.get_data()[:, :, :, 0],
- vol.affine)
- nib.save(vol, _subject_data["anat"])
+ return _glob_spm_multimodal_fmri_data(subject_dir)
- return Bunch(**_subject_data)
- # maybe data_dir already contains the data ?
- data = _glob_spm_auditory_data()
- if data is not None:
- return data
+def _make_events_filepath_spm_multimodal_fmri(_subject_data, session):
+ key = 'trials_ses{}'.format(session)
+ events_file_location = os.path.dirname(_subject_data[key])
+ events_filename = 'session{}_events.tsv'.format(session)
+ events_filepath = os.path.join(events_file_location, events_filename)
+ return events_filepath
- # No. Download the data
- print("Data absent, downloading...")
- url = ("http://www.fil.ion.ucl.ac.uk/spm/download/data/MoAEpilot/"
- "MoAEpilot.zip")
- archive_path = os.path.join(subject_dir, os.path.basename(url))
- _fetch_file(url, subject_dir)
- try:
- _uncompress_file(archive_path)
- except:
- print("Archive corrupted, trying to download it again.")
- return fetch_spm_auditory(data_dir=data_dir, data_name="",
- subject_id=subject_id)
- return _glob_spm_auditory_data()
+def _make_events_file_spm_multimodal_fmri(_subject_data, session):
+ tr = 2.
+ timing = loadmat(_subject_data["trials_ses%i" % (session)],
+ squeeze_me=True, struct_as_record=False)
+ faces_onsets = timing['onsets'][0].ravel()
+ scrambled_onsets = timing['onsets'][1].ravel()
+ onsets = np.hstack((faces_onsets, scrambled_onsets))
+ onsets *= tr # because onsets were reporting in 'scans' units
+ conditions = (['faces'] * len(faces_onsets) +
+ ['scrambled'] * len(scrambled_onsets))
+ duration = np.ones_like(onsets)
+ events = pd.DataFrame({'trial_type': conditions, 'onset': onsets,
+ 'duration': duration})
+ return events
def fetch_spm_multimodal_fmri(data_dir=None, data_name="spm_multimodal_fmri",
@@ -427,77 +650,16 @@ def fetch_spm_multimodal_fmri(data_dir=None, data_name="spm_multimodal_fmri",
"""
- data_dir = _get_dataset_dir(data_name, data_dir=data_dir,
- verbose=verbose)
+ data_dir = _get_dataset_dir(data_name, data_dir=data_dir, verbose=verbose)
subject_dir = os.path.join(data_dir, subject_id)
- def _glob_spm_multimodal_fmri_data():
- """glob data from subject_dir."""
- _subject_data = {'slice_order': 'descending'}
-
- for session in range(2):
- # glob func data for session s + 1
- session_func = sorted(glob.glob(
- os.path.join(
- subject_dir,
- ("fMRI/Session%i/fMETHODS-000%i-*-01.img" % (
- session + 1, session + 5)))))
- if len(session_func) < 390:
- print("Missing %i functional scans for session %i." % (
- 390 - len(session_func), session))
- return None
-
- _subject_data['func%i' % (session + 1)] = session_func
-
- # glob trials .mat file
- sess_trials = os.path.join(
- subject_dir,
- "fMRI/trials_ses%i.mat" % (session + 1))
- if not os.path.isfile(sess_trials):
- print("Missing session file: %s" % sess_trials)
- return None
-
- _subject_data['trials_ses%i' % (session + 1)] = sess_trials
-
- # glob for anat data
- anat = os.path.join(subject_dir, "sMRI/smri.img")
- if not os.path.isfile(anat):
- print("Missing structural image.")
- return None
-
- _subject_data["anat"] = anat
-
- return Bunch(**_subject_data)
-
# maybe data_dir already contains the data ?
- data = _glob_spm_multimodal_fmri_data()
+ data = _glob_spm_multimodal_fmri_data(subject_dir)
if data is not None:
return data
# No. Download the data
- print("Data absent, downloading...")
- urls = [
- # fmri
- ("http://www.fil.ion.ucl.ac.uk/spm/download/data/mmfaces/"
- "multimodal_fmri.zip"),
-
- # structural
- ("http://www.fil.ion.ucl.ac.uk/spm/download/data/mmfaces/"
- "multimodal_smri.zip")
- ]
-
- for url in urls:
- archive_path = os.path.join(subject_dir, os.path.basename(url))
- _fetch_file(url, subject_dir)
- try:
- _uncompress_file(archive_path)
- except:
- print("Archive corrupted, trying to download it again.")
- return fetch_spm_multimodal_fmri(data_dir=data_dir,
- data_name="",
- subject_id=subject_id)
-
- return _glob_spm_multimodal_fmri_data()
+ return _download_data_spm_multimodal(data_dir, subject_dir, subject_id)
def fetch_fiac_first_level(data_dir=None, verbose=1):
@@ -515,7 +677,7 @@ def _glob_fiac_data():
_subject_data = {}
subject_dir = os.path.join(data_dir, 'nipy-data-0.2/data/fiac/fiac0')
for session in [1, 2]:
- # glob func data for session session + 1
+ # glob func data for session
session_func = os.path.join(subject_dir, 'run%i.nii.gz' % session)
if not os.path.isfile(session_func):
print('Missing functional scan for session %i.' % session)
diff --git a/nistats/design_matrix.py b/nistats/design_matrix.py
index 696acdbf..7f461b22 100644
--- a/nistats/design_matrix.py
+++ b/nistats/design_matrix.py
@@ -5,7 +5,7 @@
Design matrices are represented by Pandas DataFrames
Computations of the different parts of the design matrix are confined
-to the make_design_matrix function, that create a DataFrame
+to the make_first_level_design_matrix function, that create a DataFrame
All the others are ancillary functions.
Design matrices contain three different types of regressors:
@@ -13,14 +13,15 @@
1. Task-related regressors, that result from the convolution
of the experimental paradigm regressors with hemodynamic models
A hemodynamic model is one of:
- 'spm' : linear filter used in the SPM software
- 'glover' : linear filter estimated by G.Glover
- 'spm + derivative', 'glover + derivative': the same linear models,
+
+ - 'spm' : linear filter used in the SPM software
+ - 'glover' : linear filter estimated by G.Glover
+ - 'spm + derivative', 'glover + derivative': the same linear models,
plus their time derivative (2 regressors per condition)
- 'spm + derivative + dispersion', 'glover + derivative + dispersion':
+ - 'spm + derivative + dispersion', 'glover + derivative + dispersion':
idem plus the derivative wrt the dispersion parameter of the hrf
(3 regressors per condition)
- 'fir' : finite impulse response model, generic linear filter
+ - 'fir' : finite impulse response model, generic linear filter
2. User-specified regressors, that represent information available on
the data, e.g. motion parameters, physiological data resampled at
@@ -32,6 +33,7 @@
estimates.
Author: Bertrand Thirion, 2009-2015
+
"""
from __future__ import with_statement
from warnings import warn
@@ -41,7 +43,7 @@
import pandas as pd
from .hemodynamic_models import compute_regressor, _orthogonalize
-from .experimental_paradigm import check_paradigm
+from .experimental_paradigm import check_events
from .utils import full_rank, _basestring
######################################################################
@@ -163,15 +165,15 @@ def _make_drift(drift_model, frame_times, order=1, period_cut=128.):
return drift, names
-def _convolve_regressors(paradigm, hrf_model, frame_times, fir_delays=[0],
- min_onset=-24):
+def _convolve_regressors(events, hrf_model, frame_times, fir_delays=[0],
+ min_onset=-24, oversampling=50):
""" Creation of a matrix that comprises
the convolution of the conditions onset with a certain hrf model
Parameters
----------
- paradigm : DataFrame instance,
- Descriptor of an experimental paradigm
+ events : DataFrame instance,
+ Events data describing the experimental paradigm
see nistats.experimental_paradigm to check the specification
for these to be valid paradigm descriptors
@@ -191,6 +193,10 @@ def _convolve_regressors(paradigm, hrf_model, frame_times, fir_delays=[0],
Minimal onset relative to frame_times[0] (in seconds) events
that start before frame_times[0] + min_onset are not considered.
+ oversampling: int or None, optional, default:50,
+ Oversampling factor used in temporal convolutions.
+ Should be 1 whenever hrf_model is 'fir'.
+
Returns
-------
regressor_matrix : array of shape (n_scans, n_regressors),
@@ -210,11 +216,14 @@ def _convolve_regressors(paradigm, hrf_model, frame_times, fir_delays=[0],
regressor_names = []
regressor_matrix = None
if hrf_model == 'fir':
+ if oversampling not in [1, None]:
+ warn('Forcing oversampling factor to 1 for a finite'
+ 'impulse response hrf model')
oversampling = 1
- else:
- oversampling = 16
+ elif oversampling is None:
+ oversampling = 50
- trial_type, onset, duration, modulation = check_paradigm(paradigm)
+ trial_type, onset, duration, modulation = check_events(events)
for condition in np.unique(trial_type):
condition_mask = (trial_type == condition)
exp_condition = (onset[condition_mask],
@@ -271,10 +280,10 @@ def _full_rank(X, cmax=1e15):
######################################################################
-def make_design_matrix(
- frame_times, paradigm=None, hrf_model='glover',
+def make_first_level_design_matrix(
+ frame_times, events=None, hrf_model='glover',
drift_model='cosine', period_cut=128, drift_order=1, fir_delays=[0],
- add_regs=None, add_reg_names=None, min_onset=-24):
+ add_regs=None, add_reg_names=None, min_onset=-24, oversampling=50):
"""Generate a design matrix from the input parameters
Parameters
@@ -282,11 +291,11 @@ def make_design_matrix(
frame_times : array of shape (n_frames,)
The timing of acquisition of the scans in seconds.
- paradigm : DataFrame instance, optional
- Description of the experimental paradigm. The DataFrame instance might
- have those keys:
+ events : DataFrame instance, optional
+ Events data that describes the experimental paradigm.
+ The DataFrame instance might have these keys:
'onset': column to specify the start time of each events in
- seconds. An exception is raised if this key is missing.
+ seconds. An error is raised if this key is missing.
'trial_type': column to specify per-event experimental conditions
identifier. If missing each event are labelled
'dummy' and considered to form a unique condition.
@@ -296,11 +305,13 @@ def make_design_matrix(
'modulation': column to specify the amplitude of each
events. If missing the default is set to
ones(n_events).
- A paradigm is considered as valid whenever it has an 'onset' key, if
- this key is missing an exception will be thrown. For the others keys
- only a simple warning will be displayed. A particular attention should
- be given to the 'trial_type' key which defines the different conditions
- in the paradigm.
+
+ An experimental paradigm is valid if it has an 'onset' key
+ and a 'duration' key.
+ If these keys are missing an error will be raised.
+ For the others keys a warning will be displayed.
+ Particular attention should be given to the 'trial_type' key
+ which defines the different conditions in the experimental paradigm.
hrf_model : {'spm', 'spm + derivative', 'spm + derivative + dispersion',
'glover', 'glover + derivative', 'glover + derivative + dispersion',
@@ -332,6 +343,10 @@ def make_design_matrix(
Minimal onset relative to frame_times[0] (in seconds)
events that start before frame_times[0] + min_onset are not considered.
+ oversampling: int or None, optional,
+ Oversampling factor used in temporal convolutions.
+ Should be 1 whenever hrf_model is 'fir'.
+
Returns
-------
design_matrix : DataFrame instance,
@@ -363,13 +378,14 @@ def make_design_matrix(
names = []
matrix = None
- # step 1: paradigm-related regressors
- if paradigm is not None:
+ # step 1: events-related regressors
+ if events is not None:
# create the condition-related regressors
if isinstance(hrf_model, _basestring):
hrf_model = hrf_model.lower()
matrix, names = _convolve_regressors(
- paradigm, hrf_model, frame_times, fir_delays, min_onset)
+ events, hrf_model, frame_times, fir_delays, min_onset,
+ oversampling)
# step 2: additional regressors
if add_regs is not None:
@@ -428,7 +444,7 @@ def check_design_matrix(design_matrix):
return frame_times, matrix, names
-def create_second_level_design(subjects_label, confounds=None):
+def make_second_level_design_matrix(subjects_label, confounds=None):
"""Sets up a second level design.
Construct a design matrix with an intercept and subject specific confounds.
@@ -478,7 +494,7 @@ def create_second_level_design(subjects_label, confounds=None):
# check design matrix is not singular
epsilon = sys.float_info.epsilon
- if np.linalg.cond(design_matrix.as_matrix()) > design_matrix.size:
+ if np.linalg.cond(design_matrix.values) > design_matrix.size:
warn('Attention: Design matrix is singular. Aberrant estimates '
'are expected.')
return design_matrix
diff --git a/nistats/experimental_paradigm.py b/nistats/experimental_paradigm.py
index afe2b003..6f8e8c51 100644
--- a/nistats/experimental_paradigm.py
+++ b/nistats/experimental_paradigm.py
@@ -4,29 +4,32 @@
An experimental protocol is handled as a pandas DataFrame
that includes an 'onset' field.
-This yields the onset time of the events in the paradigm.
+This yields the onset time of the events in the experimental paradigm.
It can also contain:
+
* a 'trial_type' field that yields the condition identifier.
* a 'duration' field that yields event duration (for so-called block
paradigms).
* a 'modulation' field that associated a scalar value to each event.
Author: Bertrand Thirion, 2015
+
"""
from __future__ import with_statement
-import warnings
import numpy as np
+import pandas
+import warnings
-def check_paradigm(paradigm):
- """Test that the DataFrame is describes a valid experimental paradigm
+def check_events(events):
+ """Test that the events data describes a valid experimental paradigm
- A DataFrame is considered as valid whenever it has an 'onset' key.
+ It is valid if the events data has an 'onset' key.
Parameters
----------
- paradigm : pandas DataFrame
- Describes a functional paradigm.
+ events : pandas DataFrame
+ Events data that describes a functional experimental paradigm.
Returns
-------
@@ -45,40 +48,20 @@ def check_paradigm(paradigm):
Per-event modulation, (in seconds)
defaults to ones(n_events) when no duration is provided
"""
- if 'onset' not in paradigm.keys():
- raise ValueError('The provided paradigm has no onset key')
+ if 'onset' not in events.keys():
+ raise ValueError('The provided events data has no onset column.')
+ if 'duration' not in events.keys():
+ raise ValueError('The provided events data has no duration column.')
- onset = np.array(paradigm['onset'])
+ onset = np.array(events['onset'])
+ duration = np.array(events['duration']).astype(np.float)
n_events = len(onset)
- trial_type = np.repeat('dummy', n_events)
- duration = np.zeros(n_events)
+ trial_type = np.array(events['trial_type'])
modulation = np.ones(n_events)
- if 'trial_type' in paradigm.keys():
- warnings.warn("'trial_type' key not found in the given paradigm.")
- trial_type = np.array(paradigm['trial_type'])
- if 'duration' in paradigm.keys():
- warnings.warn("'duration' key not found in the given paradigm.")
- duration = np.array(paradigm['duration']).astype(np.float)
- if 'modulation' in paradigm.keys():
- warnings.warn("'modulation' key not found in the given paradigm.")
- modulation = np.array(paradigm['modulation']).astype(np.float)
+ if 'trial_type' not in events.keys():
+ warnings.warn("'trial_type' column not found in the given events data.")
+ trial_type = np.repeat('dummy', n_events)
+ if 'modulation' in events.keys():
+ warnings.warn("'modulation' column found in the given events data.")
+ modulation = np.array(events['modulation']).astype(np.float)
return trial_type, onset, duration, modulation
-
-
-def paradigm_from_csv(csv_file):
- """Utility function to directly read the paradigm from a csv file
-
- This is simply meant to avoid explicitly import pandas everywhere.
-
- Parameters
- ----------
- csv_file : string,
- Path to a csv file.
-
- Returns
- -------
- paradigm : pandas DataFrame,
- Holding the paradigm information.
- """
- import pandas
- return pandas.read_csv(csv_file)
diff --git a/nistats/first_level_model.py b/nistats/first_level_model.py
index d4e9902d..c92350fc 100644
--- a/nistats/first_level_model.py
+++ b/nistats/first_level_model.py
@@ -20,7 +20,7 @@
import numpy as np
import pandas as pd
-from nibabel import Nifti1Image, AnalyzeImage
+from nibabel import Nifti1Image
from sklearn.base import BaseEstimator, TransformerMixin, clone
from sklearn.externals.joblib import Memory
@@ -31,10 +31,11 @@
from patsy import DesignInfo
from .regression import OLSModel, ARModel, SimpleRegressionResults
-from .design_matrix import make_design_matrix
+from .design_matrix import make_first_level_design_matrix
from .contrasts import _fixed_effect_contrast
-from .utils import (_basestring, _check_run_tables, get_bids_files,
- parse_bids_filename)
+from .utils import (_basestring, _check_run_tables,
+ _verify_events_file_uses_tab_separators,
+ get_bids_files, parse_bids_filename)
def mean_scaling(Y, axis=0):
@@ -161,10 +162,10 @@ class FirstLevelModel(BaseEstimator, TransformerMixin, CacheMixin):
expressed as a percentage of the t_r (time repetition), so it can have
values between 0. and 1.
- hrf_model : string, optional
- This parameter specifies the hemodynamic response function (HRF) for
- the design matrices. It can be 'canonical', 'canonical with derivative'
- or 'fir'.
+ hrf_model : {'spm', 'spm + derivative', 'spm + derivative + dispersion',
+ 'glover', 'glover + derivative', 'glover + derivative + dispersion',
+ 'fir', None}
+ String that specifies the hemodynamic response function. Defaults to 'glover'.
drift_model : string, optional
This parameter specifies the desired drift model for the design
@@ -286,6 +287,7 @@ def __init__(self, t_r=None, slice_time_ref=0., hrf_model='glover',
self.memory = memory
self.memory_level = memory_level
self.standardize = standardize
+
if signal_scaling is False:
self.signal_scaling = signal_scaling
elif signal_scaling in [0, 1, (0, 1)]:
@@ -295,6 +297,7 @@ def __init__(self, t_r=None, slice_time_ref=0., hrf_model='glover',
else:
raise ValueError('signal_scaling must be "False", "0", "1"'
' or "(0, 1)"')
+
self.noise_model = noise_model
self.verbose = verbose
self.n_jobs = n_jobs
@@ -316,18 +319,20 @@ def fit(self, run_imgs, events=None, confounds=None,
Parameters
----------
run_imgs: Niimg-like object or list of Niimg-like objects,
- See http://nilearn.github.io/building_blocks/manipulating_mr_images.html#niimg.
+ See http://nilearn.github.io/manipulating_images/input_output.html#inputing-data-file-names-or-image-objects
Data on which the GLM will be fitted. If this is a list,
the affine is considered the same for all.
events: pandas Dataframe or string or list of pandas DataFrames or
- strings,
+ strings
+
fMRI events used to build design matrices. One events object
expected per run_img. Ignored in case designs is not None.
If string, then a path to a csv file is expected.
confounds: pandas Dataframe or string or list of pandas DataFrames or
- strings,
+ strings
+
Each column in a DataFrame corresponds to a confound variable
to be included in the regression model of the respective run_img.
The number of rows must match the number of volumes in the
@@ -341,6 +346,9 @@ def fit(self, run_imgs, events=None, confounds=None,
"""
# Check arguments
# Check imgs type
+ if events is not None:
+ _verify_events_file_uses_tab_separators(
+ events_files=events)
if not isinstance(run_imgs, (list, tuple)):
run_imgs = [run_imgs]
if design_matrices is None:
@@ -356,7 +364,6 @@ def fit(self, run_imgs, events=None, confounds=None,
# Also check that events and confound files can be loaded as DataFrame
if events is not None:
events = _check_run_tables(run_imgs, events, 'events')
-
if confounds is not None:
confounds = _check_run_tables(run_imgs, confounds, 'confounds')
@@ -431,11 +438,11 @@ def fit(self, run_imgs, events=None, confounds=None,
start_time = self.slice_time_ref * self.t_r
end_time = (n_scans - 1 + self.slice_time_ref) * self.t_r
frame_times = np.linspace(start_time, end_time, n_scans)
- design = make_design_matrix(frame_times, events[run_idx],
- self.hrf_model, self.drift_model,
- self.period_cut, self.drift_order,
- self.fir_delays, confounds_matrix,
- confounds_names, self.min_onset)
+ design = make_first_level_design_matrix(frame_times, events[run_idx],
+ self.hrf_model, self.drift_model,
+ self.period_cut, self.drift_order,
+ self.fir_delays, confounds_matrix,
+ confounds_names, self.min_onset)
else:
design = design_matrices[run_idx]
self.design_matrices_.append(design)
@@ -449,11 +456,11 @@ def fit(self, run_imgs, events=None, confounds=None,
if self.verbose > 1:
t_masking = time.time() - t_masking
- sys.stderr.write('Masker took %d seconds \n' % t_masking)
+ sys.stderr.write('Masker took %d seconds \n' % t_masking)
if self.signal_scaling:
Y, _ = mean_scaling(Y, self.scaling_axis)
- if self.memory is not None:
+ if self.memory:
mem_glm = self.memory.cache(run_glm, ignore=['n_jobs'])
else:
mem_glm = run_glm
@@ -462,7 +469,7 @@ def fit(self, run_imgs, events=None, confounds=None,
if self.verbose > 1:
t_glm = time.time()
sys.stderr.write('Performing GLM computation\r')
- labels, results = mem_glm(Y, design.as_matrix(),
+ labels, results = mem_glm(Y, design.values,
noise_model=self.noise_model,
bins=100, n_jobs=self.n_jobs)
if self.verbose > 1:
@@ -494,13 +501,14 @@ def compute_contrast(self, contrast_def, stat_type=None,
----------
contrast_def : str or array of shape (n_col) or list of (string or
array of shape (n_col))
+
where ``n_col`` is the number of columns of the design matrix,
(one array per run). If only one array is provided when there
are several runs, it will be assumed that the same contrast is
desired for all runs. The string can be a formula compatible with
the linear constraint of the Patsy library. Basically one can use
the name of the conditions as they appear in the design matrix of
- the fitted model combined with operators /*+- and numbers.
+ the fitted model combined with operators /\*+- and numbers.
Please checks the patsy documentation for formula examples:
http://patsy.readthedocs.io/en/latest/API-reference.html#patsy.DesignInfo.linear_constraint
@@ -561,7 +569,7 @@ def compute_contrast(self, contrast_def, stat_type=None,
def first_level_models_from_bids(
- dataset_path, task_label, space_label, img_filters=[],
+ dataset_path, task_label, space_label, img_filters=None,
t_r=None, slice_time_ref=0., hrf_model='glover', drift_model='cosine',
period_cut=128, drift_order=1, fir_delays=[0], min_onset=-24,
mask=None, target_affine=None, target_shape=None, smoothing_fwhm=None,
@@ -587,7 +595,7 @@ def first_level_models_from_bids(
Specifies the space label of the preproc.nii images.
As they are specified in the file names like _space-_.
- img_filters: list of tuples (str, str), optional (default: [])
+ img_filters: list of tuples (str, str), optional (default: None)
Filters are of the form (field, label). Only one filter per field
allowed. A file that does not match a filter will be discarded.
Possible filters are 'acq', 'rec', 'run', 'res' and 'variant'.
@@ -619,6 +627,7 @@ def first_level_models_from_bids(
Items for the FirstLevelModel fit function of their respective model.
"""
# check arguments
+ img_filters = img_filters if img_filters else []
if not isinstance(dataset_path, str):
raise TypeError('dataset_path must be a string, instead %s was given' %
type(task_label))
diff --git a/nistats/hemodynamic_models.py b/nistats/hemodynamic_models.py
index 55842a45..22c0b931 100644
--- a/nistats/hemodynamic_models.py
+++ b/nistats/hemodynamic_models.py
@@ -3,7 +3,7 @@
Here we provide for SPM, Glover hrfs and finite timpulse response (FIR) models.
This module closely follows SPM implementation
-Author: Bertrand Thirion, 2011--2015
+Author: Bertrand Thirion, 2011--2018
"""
import warnings
@@ -11,7 +11,7 @@
from scipy.stats import gamma
-def _gamma_difference_hrf(tr, oversampling=16, time_length=32., onset=0.,
+def _gamma_difference_hrf(tr, oversampling=50, time_length=32., onset=0.,
delay=6, undershoot=16., dispersion=1.,
u_dispersion=1., ratio=0.167):
""" Compute an hrf as the difference of two gamma functions
@@ -22,7 +22,7 @@ def _gamma_difference_hrf(tr, oversampling=16, time_length=32., onset=0.,
tr : float
scan repeat time, in seconds
- oversampling : int, optional (default=16)
+ oversampling : int, optional (default=50)
temporal oversampling factor
time_length : float, optional (default=32)
@@ -52,7 +52,7 @@ def _gamma_difference_hrf(tr, oversampling=16, time_length=32., onset=0.,
hrf sampling on the oversampled time grid
"""
dt = tr / oversampling
- time_stamps = np.linspace(0, time_length, float(time_length) / dt)
+ time_stamps = np.linspace(0, time_length, np.rint(float(time_length) / dt).astype(np.int))
time_stamps -= onset
hrf = gamma.pdf(time_stamps, delay / dispersion, dt / dispersion) -\
ratio * gamma.pdf(
@@ -61,7 +61,7 @@ def _gamma_difference_hrf(tr, oversampling=16, time_length=32., onset=0.,
return hrf
-def spm_hrf(tr, oversampling=16, time_length=32., onset=0.):
+def spm_hrf(tr, oversampling=50, time_length=32., onset=0.):
""" Implementation of the SPM hrf model
Parameters
@@ -86,7 +86,7 @@ def spm_hrf(tr, oversampling=16, time_length=32., onset=0.):
return _gamma_difference_hrf(tr, oversampling, time_length, onset)
-def glover_hrf(tr, oversampling=16, time_length=32., onset=0.):
+def glover_hrf(tr, oversampling=50, time_length=32., onset=0.):
""" Implementation of the Glover hrf model
Parameters
@@ -113,7 +113,7 @@ def glover_hrf(tr, oversampling=16, time_length=32., onset=0.):
u_dispersion=.9, ratio=.35)
-def spm_time_derivative(tr, oversampling=16, time_length=32., onset=0.):
+def spm_time_derivative(tr, oversampling=50, time_length=32., onset=0.):
"""Implementation of the SPM time derivative hrf (dhrf) model
Parameters
@@ -141,7 +141,7 @@ def spm_time_derivative(tr, oversampling=16, time_length=32., onset=0.):
return dhrf
-def glover_time_derivative(tr, oversampling=16, time_length=32., onset=0.):
+def glover_time_derivative(tr, oversampling=50, time_length=32., onset=0.):
"""Implementation of the Glover time derivative hrf (dhrf) model
Parameters
@@ -166,7 +166,7 @@ def glover_time_derivative(tr, oversampling=16, time_length=32., onset=0.):
return dhrf
-def spm_dispersion_derivative(tr, oversampling=16, time_length=32., onset=0.):
+def spm_dispersion_derivative(tr, oversampling=50, time_length=32., onset=0.):
"""Implementation of the SPM dispersion derivative hrf model
Parameters
@@ -196,7 +196,7 @@ def spm_dispersion_derivative(tr, oversampling=16, time_length=32., onset=0.):
return dhrf
-def glover_dispersion_derivative(tr, oversampling=16, time_length=32.,
+def glover_dispersion_derivative(tr, oversampling=50, time_length=32.,
onset=0.):
"""Implementation of the Glover dispersion derivative hrf model
@@ -230,7 +230,7 @@ def glover_dispersion_derivative(tr, oversampling=16, time_length=32.,
return dhrf
-def _sample_condition(exp_condition, frame_times, oversampling=16,
+def _sample_condition(exp_condition, frame_times, oversampling=50,
min_onset=-24):
"""Make a possibly oversampled event regressor from condition information.
@@ -265,12 +265,13 @@ def _sample_condition(exp_condition, frame_times, oversampling=16,
min_onset) * oversampling) + 1
hr_frame_times = np.linspace(frame_times.min() + min_onset,
- frame_times.max() * (1 + 1. / (n - 1)), n_hr)
+ frame_times.max() * (1 + 1. / (n - 1)),
+ np.rint(n_hr).astype(np.int))
# Get the condition information
onsets, durations, values = tuple(map(np.asanyarray, exp_condition))
if (onsets < frame_times[0] + min_onset).any():
- warnings.warn(('Some stimulus onsets are earlier than %d in the' +
+ warnings.warn(('Some stimulus onsets are earlier than %s in the'
' experiment and are thus not considered in the model'
% (frame_times[0] + min_onset)), UserWarning)
@@ -374,7 +375,7 @@ def _regressor_names(con_name, hrf_model, fir_delays=None):
return [con_name + "_delay_%d" % i for i in fir_delays]
-def _hrf_kernel(hrf_model, tr, oversampling=16, fir_delays=None):
+def _hrf_kernel(hrf_model, tr, oversampling=50, fir_delays=None):
""" Given the specification of the hemodynamic model and time parameters,
return the list of matching kernels
@@ -432,7 +433,7 @@ def _hrf_kernel(hrf_model, tr, oversampling=16, fir_delays=None):
def compute_regressor(exp_condition, hrf_model, frame_times, con_id='cond',
- oversampling=16, fir_delays=None, min_onset=-24):
+ oversampling=50, fir_delays=None, min_onset=-24):
""" This is the main function to convolve regressors with hrf model
Parameters
@@ -472,16 +473,17 @@ def compute_regressor(exp_condition, hrf_model, frame_times, con_id='cond',
Notes
-----
The different hemodynamic models can be understood as follows:
- 'spm': this is the hrf model used in SPM
- 'spm + derivative': SPM model plus its time derivative (2 regressors)
- 'spm + time + dispersion': idem, plus dispersion derivative (3 regressors)
- 'glover': this one corresponds to the Glover hrf
- 'glover + derivative': the Glover hrf + time derivative (2 regressors)
- 'glover + derivative + dispersion': idem + dispersion derivative
+ - 'spm': this is the hrf model used in SPM
+ - 'spm + derivative': SPM model plus its time derivative (2 regressors)
+ - 'spm + time + dispersion': idem, plus dispersion derivative (3 regressors)
+ - 'glover': this one corresponds to the Glover hrf
+ - 'glover + derivative': the Glover hrf + time derivative (2 regressors)
+ - 'glover + derivative + dispersion': idem + dispersion derivative
(3 regressors)
'fir': finite impulse response basis, a set of delayed dirac models
with arbitrary length. This one currently assumes regularly spaced
frame times (i.e. fixed time of repetition).
+
It is expected that spm standard and Glover model would not yield
large differences in most cases.
diff --git a/nistats/model.py b/nistats/model.py
index 007eead0..a8bbfd79 100644
--- a/nistats/model.py
+++ b/nistats/model.py
@@ -12,7 +12,7 @@
from nibabel.onetime import setattr_on_read
-from .utils import pos_recipr
+from .utils import positive_reciprocal
# Inverse t cumulative distribution
inv_t_cdf = t_distribution.ppf
@@ -75,7 +75,7 @@ def __init__(self, theta, Y, model, cov=None, dispersion=1., nuisance=None,
self.df_model = model.df_model
# put this as a parameter of LikelihoodModel
self.df_resid = self.df_total - self.df_model
-
+
@setattr_on_read
def logL(self):
"""
@@ -98,7 +98,7 @@ def t(self, column=None):
_cov = self.vcov(column=column)
if _cov.ndim == 2:
_cov = np.diag(_cov)
- _t = _theta * pos_recipr(np.sqrt(_cov))
+ _t = _theta * positive_reciprocal(np.sqrt(_cov))
return _t
def vcov(self, matrix=None, column=None, dispersion=None, other=None):
@@ -200,7 +200,7 @@ def Tcontrast(self, matrix, store=('t', 'effect', 'sd'), dispersion=None):
if 'sd' in store:
st_sd = np.squeeze(sd)
if 't' in store:
- st_t = np.squeeze(effect * pos_recipr(sd))
+ st_t = np.squeeze(effect * positive_reciprocal(sd))
return TContrastResults(effect=st_effect, t=st_t, sd=st_sd,
df_den=self.df_resid)
@@ -258,8 +258,8 @@ def Fcontrast(self, matrix, dispersion=None, invcov=None):
q = matrix.shape[0]
if invcov is None:
invcov = inv(self.vcov(matrix=matrix, dispersion=1.0))
- F = np.add.reduce(np.dot(invcov, ctheta) * ctheta, 0) *\
- pos_recipr((q * dispersion))
+ F = np.add.reduce(np.dot(invcov, ctheta) * ctheta, 0) * \
+ positive_reciprocal((q * dispersion))
F = np.squeeze(F)
return FContrastResults(
effect=ctheta, covariance=self.vcov(
@@ -300,10 +300,12 @@ def conf_int(self, alpha=.05, cols=None, dispersion=None):
Notes
-----
+
Confidence intervals are two-tailed.
- TODO:
+
tails : string, optional
- `tails` can be "two", "upper", or "lower"
+ Possible values: 'two' | 'upper' | 'lower'
+
'''
if cols is None:
lower = self.theta - inv_t_cdf(1 - alpha / 2, self.df_resid) *\
diff --git a/nistats/regression.py b/nistats/regression.py
index 081fdffc..4635d606 100644
--- a/nistats/regression.py
+++ b/nistats/regression.py
@@ -29,7 +29,7 @@
from nibabel.onetime import setattr_on_read
-from .utils import pos_recipr
+from .utils import positive_reciprocal
from .model import LikelihoodModelResults
@@ -282,15 +282,16 @@ def __init__(self, theta, Y, model, wY, wresid, cov=None, dispersion=1.,
dispersion, nuisance)
self.wY = wY
self.wresid = wresid
+ self.wdesign = model.wdesign
- @setattr_on_read
+ @property
def resid(self):
"""
Residuals from the fit.
"""
return self.Y - self.predicted
- @setattr_on_read
+ @property
def norm_resid(self):
"""
Residuals, normalized to have unit length.
@@ -308,30 +309,30 @@ def norm_resid(self):
See: Montgomery and Peck 3.2.1 p. 68
Davidson and MacKinnon 15.2 p 662
"""
- return self.resid * pos_recipr(np.sqrt(self.dispersion))
+ return self.resid * positive_reciprocal(np.sqrt(self.dispersion))
- @setattr_on_read
+ @property
def predicted(self):
""" Return linear predictor values from a design matrix.
"""
beta = self.theta
# the LikelihoodModelResults has parameters named 'theta'
- X = self.model.design
+ X = self.wdesign
return np.dot(X, beta)
- @setattr_on_read
+ @property
def SSE(self):
"""Error sum of squares. If not from an OLS model this is "pseudo"-SSE.
"""
return (self.wresid ** 2).sum(0)
- @setattr_on_read
+ @property
def MSE(self):
""" Mean square (error) """
return self.SSE / self.df_resid
-class SimpleRegressionResults(LikelihoodModelResults):
+class SimpleRegressionResults(RegressionResults):
"""This class contains only information of the model fit necessary
for contast computation.
@@ -353,19 +354,23 @@ def __init__(self, results):
# put this as a parameter of LikelihoodModel
self.df_resid = self.df_total - self.df_model
+ self.wdesign = results.model.wdesign
+
def logL(self, Y):
"""
The maximized log-likelihood
"""
- raise ValueError('can not use this method for simple results')
+ raise NotImplementedError('logL not implemented for '
+ 'SimpleRegressionsResults. '
+ 'Use RegressionResults')
- def resid(self, Y):
+ def resid(self, Y, X):
"""
Residuals from the fit.
"""
- return Y - self.predicted
+ return Y - self.predicted(X)
- def norm_resid(self, Y):
+ def norm_resid(self, Y, X):
"""
Residuals, normalized to have unit length.
@@ -382,12 +387,11 @@ def norm_resid(self, Y):
See: Montgomery and Peck 3.2.1 p. 68
Davidson and MacKinnon 15.2 p 662
"""
- return self.resid(Y) * pos_recipr(np.sqrt(self.dispersion))
+ return self.resid(Y, X) * positive_reciprocal(np.sqrt(self.dispersion))
- def predicted(self):
+ def predicted(self, X):
""" Return linear predictor values from a design matrix.
"""
beta = self.theta
# the LikelihoodModelResults has parameters named 'theta'
- X = self.model.design
return np.dot(X, beta)
diff --git a/nistats/reporting.py b/nistats/reporting.py
index d831e264..3f0e5dff 100644
--- a/nistats/reporting.py
+++ b/nistats/reporting.py
@@ -6,15 +6,180 @@
import os
import warnings
+from string import ascii_lowercase
+
import numpy as np
+import pandas as pd
+import nibabel as nib
from scipy import stats
+from scipy import ndimage
import nilearn.plotting # overrides the backend on headless servers
-import matplotlib
+from nilearn.image.resampling import coord_transform
import matplotlib.pyplot as plt
from patsy import DesignInfo
+
from .design_matrix import check_design_matrix
+def _local_max(data, affine, min_distance):
+ """Find all local maxima of the array, separated by at least min_distance.
+ Adapted from https://stackoverflow.com/a/22631583/2589328
+
+ Parameters
+ ----------
+ data : array_like
+ 3D array of with masked values for cluster.
+
+ min_distance : `int`
+ Minimum distance between local maxima in ``data``, in terms of mm.
+
+ Returns
+ -------
+ ijk : `numpy.ndarray`
+ (n_foci, 3) array of local maxima indices for cluster.
+
+ vals : `numpy.ndarray`
+ (n_foci,) array of values from data at ijk.
+ """
+ # Initial identification of subpeaks with minimal minimum distance
+ data_max = ndimage.filters.maximum_filter(data, 3)
+ maxima = (data == data_max)
+ data_min = ndimage.filters.minimum_filter(data, 3)
+ diff = ((data_max - data_min) > 0)
+ maxima[diff == 0] = 0
+
+ labeled, n_subpeaks = ndimage.label(maxima)
+ ijk = np.array(ndimage.center_of_mass(data, labeled,
+ range(1, n_subpeaks + 1)))
+ ijk = np.round(ijk).astype(int)
+
+ vals = np.apply_along_axis(arr=ijk, axis=1, func1d=_get_val,
+ input_arr=data)
+
+ # Sort subpeaks in cluster in descending order of stat value
+ order = (-vals).argsort()
+ vals = vals[order]
+ ijk = ijk[order, :]
+ xyz = nib.affines.apply_affine(affine, ijk) # Convert to xyz in mm
+
+ # Reduce list of subpeaks based on distance
+ keep_idx = np.ones(xyz.shape[0]).astype(bool)
+ for i in range(xyz.shape[0]):
+ for j in range(i + 1, xyz.shape[0]):
+ if keep_idx[i] == 1:
+ dist = np.linalg.norm(xyz[i, :] - xyz[j, :])
+ keep_idx[j] = dist > min_distance
+ ijk = ijk[keep_idx, :]
+ vals = vals[keep_idx]
+ return ijk, vals
+
+
+def _get_val(row, input_arr):
+ """Small function for extracting values from array based on index.
+ """
+ i, j, k = row
+ return input_arr[i, j, k]
+
+
+def get_clusters_table(stat_img, stat_threshold, cluster_threshold=None,
+ min_distance=8.):
+ """Creates pandas dataframe with img cluster statistics.
+
+ Parameters
+ ----------
+ stat_img : Niimg-like object,
+ Statistical image (presumably in z- or p-scale).
+
+ stat_threshold: `float`
+ Cluster forming threshold in same scale as `stat_img` (either a
+ p-value or z-scale value).
+
+ cluster_threshold : `int` or `None`, optional
+ Cluster size threshold, in voxels.
+
+ min_distance: `float`, optional
+ Minimum distance between subpeaks in mm. Default is 8 mm.
+
+ Returns
+ -------
+ df : `pandas.DataFrame`
+ Table with peaks and subpeaks from thresholded `stat_img`. For binary
+ clusters (clusters with >1 voxel containing only one value), the table
+ reports the center of mass of the cluster, rather than any peaks/subpeaks.
+ """
+ cols = ['Cluster ID', 'X', 'Y', 'Z', 'Peak Stat', 'Cluster Size (mm3)']
+ stat_map = stat_img.get_data()
+ conn_mat = np.zeros((3, 3, 3), int) # 6-connectivity, aka NN1 or "faces"
+ conn_mat[1, 1, :] = 1
+ conn_mat[1, :, 1] = 1
+ conn_mat[:, 1, 1] = 1
+ voxel_size = np.prod(stat_img.header.get_zooms())
+
+ # Binarize using CDT
+ binarized = stat_map > stat_threshold
+ binarized = binarized.astype(int)
+
+ # If the stat threshold is too high simply return an empty dataframe
+ if np.sum(binarized) == 0:
+ warnings.warn('Attention: No clusters with stat higher than %f' %
+ stat_threshold)
+ return pd.DataFrame(columns=cols)
+
+ # Extract connected components above cluster size threshold
+ label_map = ndimage.measurements.label(binarized, conn_mat)[0]
+ clust_ids = sorted(list(np.unique(label_map)[1:]))
+ for c_val in clust_ids:
+ if cluster_threshold is not None and np.sum(label_map == c_val) < cluster_threshold:
+ stat_map[label_map == c_val] = 0
+ binarized[label_map == c_val] = 0
+
+ # If the cluster threshold is too high simply return an empty dataframe
+ # this checks for stats higher than threshold after small clusters
+ # were removed from stat_map
+ if np.sum(stat_map > stat_threshold) == 0:
+ warnings.warn('Attention: No clusters with more than %d voxels' %
+ cluster_threshold)
+ return pd.DataFrame(columns=cols)
+
+ # Now re-label and create table
+ label_map = ndimage.measurements.label(binarized, conn_mat)[0]
+ clust_ids = sorted(list(np.unique(label_map)[1:]))
+ peak_vals = np.array([np.max(stat_map * (label_map == c)) for c in clust_ids])
+ clust_ids = [clust_ids[c] for c in (-peak_vals).argsort()] # Sort by descending max value
+
+ rows = []
+ for c_id, c_val in enumerate(clust_ids):
+ cluster_mask = label_map == c_val
+ masked_data = stat_map * cluster_mask
+
+ cluster_size_mm = int(np.sum(cluster_mask) * voxel_size)
+
+ # Get peaks, subpeaks and associated statistics
+ subpeak_ijk, subpeak_vals = _local_max(masked_data, stat_img.affine,
+ min_distance=min_distance)
+ subpeak_xyz = np.asarray(coord_transform(subpeak_ijk[:, 0],
+ subpeak_ijk[:, 1],
+ subpeak_ijk[:, 2],
+ stat_img.affine)).tolist()
+ subpeak_xyz = np.array(subpeak_xyz).T
+
+ # Only report peak and, at most, top 3 subpeaks.
+ n_subpeaks = np.min((len(subpeak_vals), 4))
+ for subpeak in range(n_subpeaks):
+ if subpeak == 0:
+ row = [c_id + 1, subpeak_xyz[subpeak, 0],
+ subpeak_xyz[subpeak, 1], subpeak_xyz[subpeak, 2],
+ subpeak_vals[subpeak], cluster_size_mm]
+ else:
+ # Subpeak naming convention is cluster num + letter (1a, 1b, etc.)
+ sp_id = '{0}{1}'.format(c_id + 1, ascii_lowercase[subpeak - 1])
+ row = [sp_id, subpeak_xyz[subpeak, 0], subpeak_xyz[subpeak, 1],
+ subpeak_xyz[subpeak, 2], subpeak_vals[subpeak], '']
+ rows += [row]
+ df = pd.DataFrame(columns=cols, data=rows)
+ return df
+
+
def compare_niimgs(ref_imgs, src_imgs, masker, plot_hist=True, log=True,
ref_label="image set 1", src_label="image set 2",
output_dir=None, axes=None):
@@ -102,7 +267,7 @@ def compare_niimgs(ref_imgs, src_imgs, masker, plot_hist=True, log=True,
return corrs
-def plot_design_matrix(design_matrix, rescale=True, ax=None):
+def plot_design_matrix(design_matrix, rescale=True, ax=None, output_file=None):
"""Plot a design matrix provided as a DataFrame
Parameters
@@ -116,6 +281,11 @@ def plot_design_matrix(design_matrix, rescale=True, ax=None):
ax : axis handle, optional
Handle to axis onto which we will draw design matrix.
+ output_file: string or None, optional,
+ The name of an image file to export the plot to. Valid extensions
+ are .png, .pdf, .svg. If output_file is not None, the plot
+ is saved to a file, and the display is closed.
+
Returns
-------
ax: axis handle
@@ -126,12 +296,11 @@ def plot_design_matrix(design_matrix, rescale=True, ax=None):
from nilearn.plotting import _set_mpl_backend
# avoid unhappy pyflakes
_set_mpl_backend
- import matplotlib.pyplot as plt
# normalize the values per column for better visualization
_, X, names = check_design_matrix(design_matrix)
if rescale:
- X = X / np.maximum(1.e-12, np.sqrt(np.sum(X ** 2, 0)))
+ X = X / np.maximum(1.e-12, np.sqrt(np.sum(X ** 2, 0))) # pylint: disable=no-member
if ax is None:
plt.figure()
ax = plt.subplot(1, 1, 1)
@@ -144,24 +313,29 @@ def plot_design_matrix(design_matrix, rescale=True, ax=None):
ax.set_xticklabels(names, rotation=60, ha='right')
plt.tight_layout()
-
+ if output_file is not None:
+ plt.savefig(output_file)
+ plt.close()
+ ax = None
return ax
-def plot_contrast_matrix(contrast_def, design_matrix, colorbar=False, ax=None):
+def plot_contrast_matrix(contrast_def, design_matrix, colorbar=False, ax=None,
+ output_file=None):
"""Creates plot for contrast definition.
Parameters
----------
contrast_def : str or array of shape (n_col) or list of (string or
array of shape (n_col))
+
where ``n_col`` is the number of columns of the design matrix,
(one array per run). If only one array is provided when there
are several runs, it will be assumed that the same contrast is
desired for all runs. The string can be a formula compatible with
the linear constraint of the Patsy library. Basically one can use
the name of the conditions as they appear in the design matrix of
- the fitted model combined with operators /*+- and numbers.
+ the fitted model combined with operators /\*+- and numbers.
Please checks the patsy documentation for formula examples:
http://patsy.readthedocs.io/en/latest/API-reference.html#patsy.DesignInfo.linear_constraint
@@ -172,39 +346,50 @@ def plot_contrast_matrix(contrast_def, design_matrix, colorbar=False, ax=None):
ax: matplotlib Axes object, optional (default None)
Directory where plotted figures will be stored.
+
+ output_file: string or None, optional,
+ The name of an image file to export the plot to. Valid extensions
+ are .png, .pdf, .svg. If output_file is not None, the plot
+ is saved to a file, and the display is closed.
+
Returns
-------
Plot Axes object
+
"""
design_column_names = design_matrix.columns.tolist()
if isinstance(contrast_def, str):
- di = DesignInfo(design_column_names)
- contrast_def = di.linear_constraint(contrast_def).coefs
+ design_info = DesignInfo(design_column_names)
+ contrast_def = design_info.linear_constraint(contrast_def).coefs
+
+ maxval = np.max(np.abs(contrast_def))
+ con_matrix = np.asmatrix(contrast_def)
if ax is None:
plt.figure(figsize=(8, 4))
ax = plt.gca()
- maxval = np.max(np.abs(contrast_def))
+ mat = ax.matshow(con_matrix, aspect='equal',
+ extent=[0, con_matrix.shape[1], 0, con_matrix.shape[0]],
+ cmap='gray', vmin=-maxval, vmax=maxval)
- con_mx = np.asmatrix(contrast_def)
- mat = ax.matshow(con_mx, aspect='equal', extent=[0, con_mx.shape[1],
- 0, con_mx.shape[0]], cmap='gray', vmin=-maxval,
- vmax=maxval)
ax.set_label('conditions')
ax.set_ylabel('')
ax.set_yticklabels(['' for x in ax.get_yticklabels()])
# Shift ticks to be at 0.5, 1.5, etc
- ax.xaxis.set(ticks=np.arange(1.0, len(design_column_names) + 1.0),
- ticklabels=design_column_names)
- ax.set_xticklabels(design_column_names, rotation=90, ha='right')
+ ax.xaxis.set(ticks=np.arange(len(design_column_names)))
+ ax.set_xticklabels(design_column_names, rotation=60, ha='left')
if colorbar:
plt.colorbar(mat, fraction=0.025, pad=0.04)
plt.tight_layout()
+ if output_file is not None:
+ plt.savefig(output_file)
+ plt.close()
+ ax = None
return ax
diff --git a/nistats/second_level_model.py b/nistats/second_level_model.py
index e2bda531..ede1fec8 100644
--- a/nistats/second_level_model.py
+++ b/nistats/second_level_model.py
@@ -16,9 +16,11 @@
from nibabel import Nifti1Image
from sklearn.base import BaseEstimator, TransformerMixin, clone
+from sklearn.externals.joblib import Memory
from nilearn._utils.niimg_conversions import check_niimg
from nilearn._utils import CacheMixin
from nilearn.input_data import NiftiMasker
+from nilearn.image import mean_img
from patsy import DesignInfo
from .first_level_model import FirstLevelModel
@@ -26,7 +28,7 @@
from .regression import SimpleRegressionResults
from .contrasts import compute_contrast
from .utils import _basestring
-from .design_matrix import create_second_level_design
+from .design_matrix import make_second_level_design_matrix
def _infer_effect_maps(second_level_input, contrast_def):
@@ -100,11 +102,14 @@ class SecondLevelModel(BaseEstimator, TransformerMixin, CacheMixin):
"""
def __init__(self, mask=None, smoothing_fwhm=None,
- memory=None, memory_level=1, verbose=0,
+ memory=Memory(None), memory_level=1, verbose=0,
n_jobs=1, minimize_memory=True):
self.mask = mask
self.smoothing_fwhm = smoothing_fwhm
- self.memory = memory
+ if isinstance(memory, _basestring):
+ self.memory = Memory(memory)
+ else:
+ self.memory = memory
self.memory_level = memory_level
self.verbose = verbose
self.n_jobs = n_jobs
@@ -123,6 +128,7 @@ def fit(self, second_level_input, confounds=None, design_matrix=None):
----------
second_level_input: list of `FirstLevelModel` objects or pandas
DataFrame or list of Niimg-like objects.
+
Giving FirstLevelModel objects will allow to easily compute
the second level contast of arbitrary first level contrasts thanks
to the first_level_contrast argument of the compute_contrast
@@ -155,7 +161,7 @@ def fit(self, second_level_input, confounds=None, design_matrix=None):
from second_level_input.
Ensure that the order of maps given by a second_level_input
list of Niimgs matches the order of the rows in the design matrix.
- Must contain a column of 1s with column name 'intercept'.
+
"""
# Check parameters
# check first level input
@@ -210,6 +216,12 @@ def fit(self, second_level_input, confounds=None, design_matrix=None):
if not isinstance(labels_dtype, np.object):
raise ValueError('subject_label column must be of dtype '
'object instead of dtype %s' % labels_dtype)
+ elif isinstance(second_level_input, (str, Nifti1Image)):
+ if design_matrix is None:
+ raise ValueError('List of niimgs as second_level_input'
+ ' require a design matrix to be provided')
+ second_level_input = check_niimg(niimg=second_level_input,
+ ensure_ndim=4)
else:
raise ValueError('second_level_input must be a list of'
' `FirstLevelModel` objects, a pandas DataFrame'
@@ -238,18 +250,14 @@ def fit(self, second_level_input, confounds=None, design_matrix=None):
if design_matrix is not None:
if not isinstance(design_matrix, pd.DataFrame):
raise ValueError('design matrix must be a pandas DataFrame')
- if 'intercept' not in design_matrix.columns:
- raise ValueError('design matrix must contain "intercept"')
# sort a pandas dataframe by subject_label to avoid inconsistencies
# with the design matrix row order when automatically extracting maps
if isinstance(second_level_input, pd.DataFrame):
- # Avoid pandas df.sort_value to keep compatibility with numpy 1.8
- # also pandas df.sort since it is completely deprecated.
columns = second_level_input.columns.tolist()
column_index = columns.index('subject_label')
sorted_matrix = sorted(
- second_level_input.as_matrix(), key=lambda x: x[column_index])
+ second_level_input.values, key=lambda x: x[column_index])
sorted_input = pd.DataFrame(sorted_matrix, columns=columns)
second_level_input = sorted_input
@@ -267,6 +275,8 @@ def fit(self, second_level_input, confounds=None, design_matrix=None):
sample_map = second_level_input['effects_map_path'][0]
labels = second_level_input['subject_label']
subjects_label = labels.values.tolist()
+ elif isinstance(second_level_input, Nifti1Image):
+ sample_map = mean_img(second_level_input)
elif isinstance(second_level_input[0], FirstLevelModel):
sample_model = second_level_input[0]
sample_condition = sample_model.design_matrices_[0].columns[0]
@@ -276,12 +286,12 @@ def fit(self, second_level_input, confounds=None, design_matrix=None):
subjects_label = labels
else:
# In this case design matrix had to be provided
- sample_map = second_level_input[0]
+ sample_map = mean_img(second_level_input)
# Create and set design matrix, if not given
if design_matrix is None:
- design_matrix = create_second_level_design(subjects_label,
- confounds)
+ design_matrix = make_second_level_design_matrix(subjects_label,
+ confounds)
self.design_matrix_ = design_matrix
# Learn the mask. Assume the first level imgs have been masked.
@@ -309,7 +319,7 @@ def fit(self, second_level_input, confounds=None, design_matrix=None):
return self
def compute_contrast(
- self, second_level_contrast='intercept', first_level_contrast=None,
+ self, second_level_contrast=None, first_level_contrast=None,
second_level_stat_type=None, output_type='z_score'):
"""Generate different outputs corresponding to
the contrasts provided e.g. z_map, t_map, effects and variance.
@@ -321,17 +331,17 @@ def compute_contrast(
The string can be a formula compatible with the linear constraint
of the Patsy library. Basically one can use the name of the
conditions as they appear in the design matrix of
- the fitted model combined with operators /*+- and numbers.
+ the fitted model combined with operators /\*+- and numbers.
Please check the patsy documentation for formula examples:
http://patsy.readthedocs.io/en/latest/API-reference.html#patsy.DesignInfo.linear_constraint
-
- VERY IMPORTANT: The 'intercept' corresponds to the second level
- effect after taking confounders in consideration. If there are
- no confounders then this will be equivalent to a simple t test.
- By default we compute the 'intercept' second level contrast.
+ The default (None) is accepted if the design matrix has a single
+ column, in which case the only possible contrast array([1]) is
+ applied; when the design matrix has multiple columns, an error is
+ raised.
first_level_contrast: str or array of shape (n_col) with respect to
FirstLevelModel, optional
+
In case a list of FirstLevelModel was provided as
second_level_input, we have to provide a contrast to apply to
the first level models to get the corresponding list of images
@@ -366,6 +376,11 @@ def compute_contrast(
'compute_contrast method of FirstLevelModel')
# check contrast definition
+ if second_level_contrast is None:
+ if self.design_matrix_.shape[1] == 1:
+ second_level_contrast = np.ones([1])
+ else:
+ raise ValueError('No second-level contrast is specified.')
if isinstance(second_level_contrast, np.ndarray):
con_val = second_level_contrast
if np.all(con_val == 0):
@@ -388,21 +403,20 @@ def compute_contrast(
# Get effect_maps appropriate for chosen contrast
effect_maps = _infer_effect_maps(self.second_level_input_,
first_level_contrast)
- # check design matrix X and effect maps Y agree on number of rows
+ # Check design matrix X and effect maps Y agree on number of rows
if len(effect_maps) != self.design_matrix_.shape[0]:
raise ValueError(
'design_matrix does not match the number of maps considered. '
'%i rows in design matrix do not match with %i maps' %
(self.design_matrix_.shape[0], len(effect_maps)))
- # Fit an OLS regression for parametric statistics
+ # Fit an Ordinary Least Squares regression for parametric statistics
Y = self.masker_.transform(effect_maps)
- if self.memory is not None:
- arg_ignore = ['n_jobs']
- mem_glm = self.memory.cache(run_glm, ignore=arg_ignore)
+ if self.memory:
+ mem_glm = self.memory.cache(run_glm, ignore=['n_jobs'])
else:
mem_glm = run_glm
- labels, results = mem_glm(Y, self.design_matrix_.as_matrix(),
+ labels, results = mem_glm(Y, self.design_matrix_.values,
n_jobs=self.n_jobs, noise_model='ols')
# We save memory if inspecting model details is not necessary
if self.minimize_memory:
@@ -412,7 +426,7 @@ def compute_contrast(
self.results_ = results
# We compute contrast object
- if self.memory is not None:
+ if self.memory:
mem_contrast = self.memory.cache(compute_contrast)
else:
mem_contrast = compute_contrast
diff --git a/nistats/tests/test_datasets.py b/nistats/tests/test_datasets.py
index 9c50dd7e..842cb712 100644
--- a/nistats/tests/test_datasets.py
+++ b/nistats/tests/test_datasets.py
@@ -1,8 +1,7 @@
import os
import json
+from tempfile import NamedTemporaryFile
import zipfile
-import numpy as np
-from nose.tools import assert_true, assert_false, assert_equal, assert_raises
import nibabel
from nilearn._utils.testing import (mock_request, wrap_chunk_read_,
@@ -12,6 +11,9 @@
from nilearn.datasets.utils import _get_dataset_dir
from nilearn._utils.compat import _basestring
from nose import with_setup
+from nose.tools import assert_true, assert_false, assert_equal, assert_raises
+import numpy as np
+import pandas as pd
from nistats import datasets
@@ -148,13 +150,92 @@ def test_fetch_openneuro_dataset():
assert_true(len(dl_files) == 9)
+def test_make_events_file_localizer_first_level():
+ def _input_data_for_test_file():
+ file_data = [
+ [0, 'calculvideo', 0.0],
+ [0, 'calculvideo', 2.400000095],
+ [0, 'damier_H', 8.699999809],
+ [0, 'clicDaudio', 11.39999961],
+ ]
+ return pd.DataFrame(file_data)
+
+ def _expected_output_data_from_test_file():
+ file_data = [
+ ['calculvideo', 0.0, 1.0],
+ ['calculvideo', 2.400000095, 1.0],
+ ['damier_H', 8.699999809, 1.0],
+ ['clicDaudio', 11.39999961, 1.0],
+ ]
+ file_data = pd.DataFrame(file_data)
+ file_data.columns = ['trial_type', 'onset', 'duration']
+ return file_data
+
+ def run_test():
+ data_for_tests = _input_data_for_test_file()
+ expected_data_from_test_file = _expected_output_data_from_test_file()
+ with NamedTemporaryFile(mode='w',
+ dir=os.getcwd(),
+ suffix='.csv') as temp_csv_obj:
+ data_for_tests.to_csv(temp_csv_obj.name,
+ index=False,
+ header=False,
+ sep=' ',
+ )
+ datasets._make_events_file_localizer_first_level(
+ temp_csv_obj.name
+ )
+ data_from_test_file_post_mod = pd.read_csv(temp_csv_obj.name,
+ sep='\t')
+ assert_true(all(
+ expected_data_from_test_file == data_from_test_file_post_mod
+ )
+ )
+
+ run_test()
+
+
@with_setup(setup_mock, teardown_mock)
def test_fetch_localizer():
dataset = datasets.fetch_localizer_first_level()
- assert_true(isinstance(dataset.paradigm, _basestring))
+ assert_true(isinstance(dataset['events'], _basestring))
assert_true(isinstance(dataset.epi_img, _basestring))
+def test_make_spm_auditory_events_file():
+ def create_expected_data():
+ expected_events_data = {
+ 'onset': [factor * 42.0 for factor in range(0, 16)],
+ 'duration': [42.0] * 16,
+ 'trial_type': ['rest', 'active'] * 8,
+ }
+ expected_events_data = pd.DataFrame(expected_events_data)
+ expected_events_data_string = expected_events_data.to_csv(
+ sep='\t',
+ index=0,
+ columns=['onset', 'duration', 'trial_type'],
+ )
+ return expected_events_data_string
+
+ def create_actual_data():
+ events_filepath = os.path.join(os.getcwd(), 'tests_events.tsv')
+ datasets._make_events_file_spm_auditory_data(
+ events_filepath=events_filepath)
+ with open(events_filepath, 'r') as actual_events_file_obj:
+ actual_events_data_string = actual_events_file_obj.read()
+ return actual_events_data_string, events_filepath
+
+ def run_test():
+ try:
+ actual_events_data_string, events_filepath = create_actual_data()
+ finally:
+ os.remove(events_filepath)
+ expected_events_data_string = create_expected_data()
+ assert_equal(actual_events_data_string, expected_events_data_string)
+
+ run_test()
+
+
@with_setup(setup_mock, teardown_mock)
@with_setup(tst.setup_tmpdata, tst.teardown_tmpdata)
def test_fetch_spm_auditory():
diff --git a/nistats/tests/test_dmtx.py b/nistats/tests/test_dmtx.py
index 1d067173..9133182e 100644
--- a/nistats/tests/test_dmtx.py
+++ b/nistats/tests/test_dmtx.py
@@ -13,9 +13,9 @@
from nilearn._utils.testing import assert_raises_regex
from nistats.design_matrix import (
- _convolve_regressors, make_design_matrix,
+ _convolve_regressors, make_first_level_design_matrix,
_cosine_drift, check_design_matrix,
- create_second_level_design)
+ make_second_level_design_matrix)
from nibabel.tmpdirs import InTemporaryDirectory
@@ -30,14 +30,14 @@
def design_matrix_light(
- frame_times, paradigm=None, hrf_model='glover',
+ frame_times, events=None, hrf_model='glover',
drift_model='cosine', period_cut=128, drift_order=1, fir_delays=[0],
add_regs=None, add_reg_names=None, min_onset=-24, path=None):
- """ Idem make_design_matrix, but only returns the computed matrix
+ """ Idem make_first_level_design_matrix, but only returns the computed matrix
and associated names """
- dmtx = make_design_matrix(frame_times, paradigm, hrf_model,
- drift_model, period_cut, drift_order, fir_delays,
- add_regs, add_reg_names, min_onset)
+ dmtx = make_first_level_design_matrix(frame_times, events, hrf_model,
+ drift_model, period_cut, drift_order, fir_delays,
+ add_regs, add_reg_names, min_onset)
_, matrix, names = check_design_matrix(dmtx)
return matrix, names
@@ -45,41 +45,45 @@ def design_matrix_light(
def basic_paradigm():
conditions = ['c0', 'c0', 'c0', 'c1', 'c1', 'c1', 'c2', 'c2', 'c2']
onsets = [30, 70, 100, 10, 30, 90, 30, 40, 60]
- paradigm = pd.DataFrame({'trial_type': conditions,
- 'onset': onsets})
- return paradigm
+ durations = 1 * np.ones(9)
+ events = pd.DataFrame({'trial_type': conditions,
+ 'onset': onsets,
+ 'duration': durations})
+ return events
def modulated_block_paradigm():
conditions = ['c0', 'c0', 'c0', 'c1', 'c1', 'c1', 'c2', 'c2', 'c2']
onsets = [30, 70, 100, 10, 30, 90, 30, 40, 60]
- duration = 5 + 5 * np.random.rand(len(onsets))
+ durations = 5 + 5 * np.random.rand(len(onsets))
values = 1 + np.random.rand(len(onsets))
- paradigm = pd.DataFrame({'trial_type': conditions,
+ events = pd.DataFrame({'trial_type': conditions,
'onset': onsets,
- 'duration': duration,
+ 'duration': durations,
'modulation': values})
- return paradigm
+ return events
def modulated_event_paradigm():
conditions = ['c0', 'c0', 'c0', 'c1', 'c1', 'c1', 'c2', 'c2', 'c2']
onsets = [30, 70, 100, 10, 30, 90, 30, 40, 60]
+ durations = 1 * np.ones(9)
values = 1 + np.random.rand(len(onsets))
- paradigm = pd.DataFrame({'trial_type': conditions,
+ events = pd.DataFrame({'trial_type': conditions,
'onset': onsets,
+ 'duration': durations,
'modulation': values})
- return paradigm
+ return events
def block_paradigm():
conditions = ['c0', 'c0', 'c0', 'c1', 'c1', 'c1', 'c2', 'c2', 'c2']
onsets = [30, 70, 100, 10, 30, 90, 30, 40, 60]
- duration = 5 * np.ones(9)
- paradigm = pd.DataFrame({'trial_type': conditions,
+ durations = 5 * np.ones(9)
+ events = pd.DataFrame({'trial_type': conditions,
'onset': onsets,
- 'duration': duration})
- return paradigm
+ 'duration': durations})
+ return events
def test_cosine_drift():
@@ -94,10 +98,10 @@ def test_cosine_drift():
def test_design_matrix0():
- # Test design matrix creation when no paradigm is provided
+ # Test design matrix creation when no experimental paradigm is provided
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128)
- _, X, names = check_design_matrix(make_design_matrix(
+ _, X, names = check_design_matrix(make_first_level_design_matrix(
frame_times, drift_model='polynomial', drift_order=3))
assert_equal(len(names), 4)
x = np.linspace(- 0.5, .5, 128)
@@ -109,7 +113,7 @@ def test_design_matrix0c():
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128)
ax = np.random.randn(128, 4)
- _, X, names = check_design_matrix(make_design_matrix(
+ _, X, names = check_design_matrix(make_first_level_design_matrix(
frame_times, drift_model='polynomial',
drift_order=3, add_regs=ax))
assert_almost_equal(X[:, 0], ax[:, 0])
@@ -117,12 +121,12 @@ def test_design_matrix0c():
assert_raises_regex(
AssertionError,
"Incorrect specification of additional regressors:.",
- make_design_matrix, frame_times, add_regs=ax)
+ make_first_level_design_matrix, frame_times, add_regs=ax)
ax = np.random.randn(128, 4)
assert_raises_regex(
ValueError,
"Incorrect number of additional regressor names.",
- make_design_matrix, frame_times, add_regs=ax, add_reg_names='')
+ make_first_level_design_matrix, frame_times, add_regs=ax, add_reg_names='')
def test_design_matrix0d():
@@ -130,7 +134,7 @@ def test_design_matrix0d():
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128)
ax = np.random.randn(128, 4)
- _, X, names = check_design_matrix(make_design_matrix(
+ _, X, names = check_design_matrix(make_first_level_design_matrix(
frame_times, drift_model='polynomial', drift_order=3, add_regs=ax))
assert_equal(len(names), 8)
assert_equal(X.shape[1], 8)
@@ -140,11 +144,12 @@ def test_convolve_regressors():
# tests for convolve_regressors helper function
conditions = ['c0', 'c1']
onsets = [20, 40]
- paradigm = pd.DataFrame({'trial_type': conditions,
- 'onset': onsets})
+ duration = [1, 1]
+ events = pd.DataFrame(
+ {'trial_type': conditions, 'onset': onsets, 'duration': duration})
# names not passed -> default names
frame_times = np.arange(100)
- f, names = _convolve_regressors(paradigm, 'glover', frame_times)
+ f, names = _convolve_regressors(events, 'glover', frame_times)
assert_equal(names, ['c0', 'c1'])
@@ -152,9 +157,9 @@ def test_design_matrix1():
# basic test based on basic_paradigm and glover hrf
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128)
- paradigm = basic_paradigm()
+ events = basic_paradigm()
hrf_model = 'glover'
- X, names = design_matrix_light(frame_times, paradigm, hrf_model=hrf_model,
+ X, names = design_matrix_light(frame_times, events, hrf_model=hrf_model,
drift_model='polynomial', drift_order=3)
assert_equal(len(names), 7)
assert_equal(X.shape, (128, 7))
@@ -166,9 +171,9 @@ def test_design_matrix2():
# idem test_design_matrix1 with a different drift term
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128)
- paradigm = basic_paradigm()
+ events = basic_paradigm()
hrf_model = 'glover'
- X, names = design_matrix_light(frame_times, paradigm, hrf_model=hrf_model,
+ X, names = design_matrix_light(frame_times, events, hrf_model=hrf_model,
drift_model='cosine', period_cut=63)
assert_equal(len(names), 7) # was 8 with old cosine
@@ -177,9 +182,9 @@ def test_design_matrix3():
# idem test_design_matrix1 with a different drift term
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128)
- paradigm = basic_paradigm()
+ events = basic_paradigm()
hrf_model = 'glover'
- X, names = design_matrix_light(frame_times, paradigm, hrf_model=hrf_model,
+ X, names = design_matrix_light(frame_times, events, hrf_model=hrf_model,
drift_model=None)
assert_equal(len(names), 4)
@@ -188,46 +193,48 @@ def test_design_matrix4():
# idem test_design_matrix1 with a different hrf model
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128)
- paradigm = basic_paradigm()
+ events = basic_paradigm()
hrf_model = 'glover + derivative'
- X, names = design_matrix_light(frame_times, paradigm, hrf_model=hrf_model,
+ X, names = design_matrix_light(frame_times, events, hrf_model=hrf_model,
drift_model='polynomial', drift_order=3)
assert_equal(len(names), 10)
def test_design_matrix5():
- # idem test_design_matrix1 with a block paradigm
+ # idem test_design_matrix1 with a block experimental paradigm
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128)
- paradigm = block_paradigm()
+ events = block_paradigm()
hrf_model = 'glover'
- X, names = design_matrix_light(frame_times, paradigm, hrf_model=hrf_model,
+ X, names = design_matrix_light(frame_times, events, hrf_model=hrf_model,
drift_model='polynomial', drift_order=3)
assert_equal(len(names), 7)
def test_design_matrix6():
- # idem test_design_matrix1 with a block paradigm and the hrf derivative
+ # idem test_design_matrix1 with a block experimental paradigm and the hrf derivative
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128)
- paradigm = block_paradigm()
+ events = block_paradigm()
hrf_model = 'glover + derivative'
- X, names = design_matrix_light(frame_times, paradigm, hrf_model=hrf_model,
+ X, names = design_matrix_light(frame_times, events, hrf_model=hrf_model,
drift_model='polynomial', drift_order=3)
assert_equal(len(names), 10)
def test_design_matrix7():
- # idem test_design_matrix1, but odd paradigm
+ # idem test_design_matrix1, but odd experimental paradigm
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128)
conditions = [0, 0, 0, 1, 1, 1, 3, 3, 3]
+ durations = 1 * np.ones(9)
# no condition 'c2'
onsets = [30, 70, 100, 10, 30, 90, 30, 40, 60]
- paradigm = pd.DataFrame({'trial_type': conditions,
- 'onset': onsets})
+ events = pd.DataFrame({'trial_type': conditions,
+ 'onset': onsets,
+ 'duration':durations})
hrf_model = 'glover'
- X, names = design_matrix_light(frame_times, paradigm, hrf_model=hrf_model,
+ X, names = design_matrix_light(frame_times, events, hrf_model=hrf_model,
drift_model='polynomial', drift_order=3)
assert_equal(len(names), 7)
@@ -236,9 +243,9 @@ def test_design_matrix8():
# basic test based on basic_paradigm and FIR
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128)
- paradigm = basic_paradigm()
+ events = basic_paradigm()
hrf_model = 'FIR'
- X, names = design_matrix_light(frame_times, paradigm, hrf_model=hrf_model,
+ X, names = design_matrix_light(frame_times, events, hrf_model=hrf_model,
drift_model='polynomial', drift_order=3)
assert_equal(len(names), 7)
@@ -247,9 +254,9 @@ def test_design_matrix9():
# basic test based on basic_paradigm and FIR
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128)
- paradigm = basic_paradigm()
+ events = basic_paradigm()
hrf_model = 'FIR'
- X, names = design_matrix_light(frame_times, paradigm, hrf_model=hrf_model,
+ X, names = design_matrix_light(frame_times, events, hrf_model=hrf_model,
drift_model='polynomial', drift_order=3,
fir_delays=range(1, 5))
assert_equal(len(names), 16)
@@ -259,12 +266,12 @@ def test_design_matrix10():
# Check that the first column o FIR design matrix is OK
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128)
- paradigm = basic_paradigm()
+ events = basic_paradigm()
hrf_model = 'FIR'
- X, names = design_matrix_light(frame_times, paradigm, hrf_model=hrf_model,
+ X, names = design_matrix_light(frame_times, events, hrf_model=hrf_model,
drift_model='polynomial', drift_order=3,
fir_delays=range(1, 5))
- onset = paradigm.onset[paradigm.trial_type == 'c0'].astype(np.int)
+ onset = events.onset[events.trial_type == 'c0'].astype(np.int)
assert_true(np.all((X[onset + 1, 0] == 1)))
@@ -272,12 +279,12 @@ def test_design_matrix11():
# check that the second column of the FIR design matrix is OK indeed
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128)
- paradigm = basic_paradigm()
+ events = basic_paradigm()
hrf_model = 'FIR'
- X, names = design_matrix_light(frame_times, paradigm, hrf_model=hrf_model,
+ X, names = design_matrix_light(frame_times, events, hrf_model=hrf_model,
drift_model='polynomial', drift_order=3,
fir_delays=range(1, 5))
- onset = paradigm.onset[paradigm.trial_type == 'c0'].astype(np.int)
+ onset = events.onset[events.trial_type == 'c0'].astype(np.int)
assert_true(np.all(X[onset + 3, 2] == 1))
@@ -285,12 +292,12 @@ def test_design_matrix12():
# check that the 11th column of a FIR design matrix is indeed OK
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128)
- paradigm = basic_paradigm()
+ events = basic_paradigm()
hrf_model = 'FIR'
- X, names = design_matrix_light(frame_times, paradigm, hrf_model=hrf_model,
+ X, names = design_matrix_light(frame_times, events, hrf_model=hrf_model,
drift_model='polynomial', drift_order=3,
fir_delays=range(1, 5))
- onset = paradigm.onset[paradigm.trial_type == 'c2'].astype(np.int)
+ onset = events.onset[events.trial_type == 'c2'].astype(np.int)
assert_true(np.all(X[onset + 4, 11] == 1))
@@ -298,12 +305,12 @@ def test_design_matrix13():
# Check that the fir_duration is well taken into account
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128)
- paradigm = basic_paradigm()
+ events = basic_paradigm()
hrf_model = 'FIR'
- X, names = design_matrix_light(frame_times, paradigm, hrf_model=hrf_model,
+ X, names = design_matrix_light(frame_times, events, hrf_model=hrf_model,
drift_model='polynomial', drift_order=3,
fir_delays=range(1, 5))
- onset = paradigm.onset[paradigm.trial_type == 'c0'].astype(np.int)
+ onset = events.onset[events.trial_type == 'c0'].astype(np.int)
assert_true(np.all(X[onset + 1, 0] == 1))
@@ -312,12 +319,12 @@ def test_design_matrix14():
# time shift
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128) + tr / 2
- paradigm = basic_paradigm()
+ events = basic_paradigm()
hrf_model = 'FIR'
- X, names = design_matrix_light(frame_times, paradigm, hrf_model=hrf_model,
+ X, names = design_matrix_light(frame_times, events, hrf_model=hrf_model,
drift_model='polynomial', drift_order=3,
fir_delays=range(1, 5))
- onset = paradigm.onset[paradigm.trial_type == 'c0'].astype(np.int)
+ onset = events.onset[events.trial_type == 'c0'].astype(np.int)
assert_true(np.all(X[onset + 1, 0] > .9))
@@ -325,10 +332,10 @@ def test_design_matrix15():
# basic test based on basic_paradigm, plus user supplied regressors
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128)
- paradigm = basic_paradigm()
+ events = basic_paradigm()
hrf_model = 'glover'
ax = np.random.randn(128, 4)
- X, names = design_matrix_light(frame_times, paradigm, hrf_model=hrf_model,
+ X, names = design_matrix_light(frame_times, events, hrf_model=hrf_model,
drift_model='polynomial', drift_order=3, add_regs=ax)
assert_equal(len(names), 11)
assert_equal(X.shape[1], 11)
@@ -338,10 +345,10 @@ def test_design_matrix16():
# Check that additional regressors are put at the right place
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128)
- paradigm = basic_paradigm()
+ events = basic_paradigm()
hrf_model = 'glover'
ax = np.random.randn(128, 4)
- X, names = design_matrix_light(frame_times, paradigm, hrf_model=hrf_model,
+ X, names = design_matrix_light(frame_times, events, hrf_model=hrf_model,
drift_model='polynomial', drift_order=3, add_regs=ax)
assert_almost_equal(X[:, 3: 7], ax)
@@ -350,11 +357,11 @@ def test_design_matrix17():
# Test the effect of scaling on the events
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128)
- paradigm = modulated_event_paradigm()
+ events = modulated_event_paradigm()
hrf_model = 'glover'
- X, names = design_matrix_light(frame_times, paradigm, hrf_model=hrf_model,
+ X, names = design_matrix_light(frame_times, events, hrf_model=hrf_model,
drift_model='polynomial', drift_order=3)
- ct = paradigm.onset[paradigm.trial_type == 'c0'].astype(np.int) + 1
+ ct = events.onset[events.trial_type == 'c0'].astype(np.int) + 1
assert_true((X[ct, 0] > 0).all())
@@ -362,11 +369,11 @@ def test_design_matrix18():
# Test the effect of scaling on the blocks
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128)
- paradigm = modulated_block_paradigm()
+ events = modulated_block_paradigm()
hrf_model = 'glover'
- X, names = design_matrix_light(frame_times, paradigm, hrf_model=hrf_model,
+ X, names = design_matrix_light(frame_times, events, hrf_model=hrf_model,
drift_model='polynomial', drift_order=3)
- ct = paradigm.onset[paradigm.trial_type == 'c0'].astype(np.int) + 3
+ ct = events.onset[events.trial_type == 'c0'].astype(np.int) + 3
assert_true((X[ct, 0] > 0).all())
@@ -374,21 +381,21 @@ def test_design_matrix19():
# Test the effect of scaling on a FIR model
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128)
- paradigm = modulated_event_paradigm()
+ events = modulated_event_paradigm()
hrf_model = 'FIR'
- X, names = design_matrix_light(frame_times, paradigm, hrf_model=hrf_model,
+ X, names = design_matrix_light(frame_times, events, hrf_model=hrf_model,
drift_model='polynomial', drift_order=3,
fir_delays=range(1, 5))
- idx = paradigm.onset[paradigm.trial_type == 0].astype(np.int)
+ idx = events.onset[events.trial_type == 0].astype(np.int)
assert_array_equal(X[idx + 1, 0], X[idx + 2, 1])
def test_design_matrix20():
# Test for commit 10662f7
frame_times = np.arange(0, 128) # was 127 in old version of _cosine_drift
- paradigm = modulated_event_paradigm()
+ events = modulated_event_paradigm()
X, names = design_matrix_light(
- frame_times, paradigm, hrf_model='glover', drift_model='cosine')
+ frame_times, events, hrf_model='glover', drift_model='cosine')
# check that the drifts are not constant
assert_true(np.all(np.diff(X[:, -2]) != 0))
@@ -398,10 +405,10 @@ def test_design_matrix21():
# basic test on repeated names of user supplied regressors
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128)
- paradigm = basic_paradigm()
+ events = basic_paradigm()
hrf_model = 'glover'
ax = np.random.randn(128, 4)
- assert_raises(ValueError, design_matrix_light, frame_times, paradigm,
+ assert_raises(ValueError, design_matrix_light, frame_times, events,
hrf_model=hrf_model, drift_model='polynomial', drift_order=3,
add_regs=ax, add_reg_names=['aha'] * ax.shape[1])
@@ -421,14 +428,37 @@ def test_fir_block():
assert_true((X[idx + 2, 6] == 1).all())
assert_true((X[idx + 3, 7] == 1).all())
+def test_oversampling():
+ events = basic_paradigm()
+ frame_times = np.linspace(0, 127, 128)
+ X1 = make_first_level_design_matrix(
+ frame_times, events, drift_model=None)
+ X2 = make_first_level_design_matrix(
+ frame_times, events, drift_model=None, oversampling=50)
+ X3 = make_first_level_design_matrix(
+ frame_times, events, drift_model=None, oversampling=10)
+
+ # oversampling = 16 is the default so X2 = X1, X3 \neq X1, X3 close to X2
+ assert_almost_equal(X1.values, X2.values)
+ assert_almost_equal(X2.values, X3.values, 0)
+ assert_true(np.linalg.norm(X2.values - X3.values) / np.linalg.norm(X2.values) > 1.e-4)
+
+ # fir model, oversampling is forced to 1
+ X4 = make_first_level_design_matrix(
+ frame_times, events, hrf_model='fir', drift_model=None,
+ fir_delays=range(0, 4), oversampling=1)
+ X5 = make_first_level_design_matrix(
+ frame_times, events, hrf_model='fir', drift_model=None,
+ fir_delays=range(0, 4), oversampling=3)
+ assert_almost_equal(X4.values, X5.values)
def test_csv_io():
# test the csv io on design matrices
tr = 1.0
frame_times = np.linspace(0, 127 * tr, 128)
- paradigm = modulated_event_paradigm()
- DM = make_design_matrix(frame_times, paradigm, hrf_model='glover',
- drift_model='polynomial', drift_order=3)
+ events = modulated_event_paradigm()
+ DM = make_first_level_design_matrix(frame_times, events, hrf_model='glover',
+ drift_model='polynomial', drift_order=3)
path = 'design_matrix.csv'
with InTemporaryDirectory():
DM.to_csv(path)
@@ -446,9 +476,11 @@ def test_spm_1():
frame_times = np.linspace(0, 99, 100)
conditions = ['c0', 'c0', 'c0', 'c1', 'c1', 'c1', 'c2', 'c2', 'c2']
onsets = [30, 50, 70, 10, 30, 80, 30, 40, 60]
- paradigm = pd.DataFrame({'trial_type': conditions,
- 'onset': onsets})
- X1 = make_design_matrix(frame_times, paradigm, drift_model=None)
+ durations = 1 * np.ones(9)
+ events = pd.DataFrame({'trial_type': conditions,
+ 'onset': onsets,
+ 'duration': durations})
+ X1 = make_first_level_design_matrix(frame_times, events, drift_model=None)
_, matrix, _ = check_design_matrix(X1)
spm_design_matrix = DESIGN_MATRIX['arr_0']
assert_true(((spm_design_matrix - matrix) ** 2).sum() /
@@ -461,11 +493,11 @@ def test_spm_2():
frame_times = np.linspace(0, 99, 100)
conditions = ['c0', 'c0', 'c0', 'c1', 'c1', 'c1', 'c2', 'c2', 'c2']
onsets = [30, 50, 70, 10, 30, 80, 30, 40, 60]
- duration = 10 * np.ones(9)
- paradigm = pd.DataFrame({'trial_type': conditions,
+ durations = 10 * np.ones(9)
+ events = pd.DataFrame({'trial_type': conditions,
'onset': onsets,
- 'duration': duration})
- X1 = make_design_matrix(frame_times, paradigm, drift_model=None)
+ 'duration': durations})
+ X1 = make_first_level_design_matrix(frame_times, events, drift_model=None)
spm_design_matrix = DESIGN_MATRIX['arr_1']
_, matrix, _ = check_design_matrix(X1)
assert_true(((spm_design_matrix - matrix) ** 2).sum() /
@@ -487,7 +519,7 @@ def test_create_second_level_design():
subjects_label = ['02', '01'] # change order to test right output order
regressors = [['01', 0.1], ['02', 0.75]]
regressors = pd.DataFrame(regressors, columns=['subject_label', 'f1'])
- design = create_second_level_design(subjects_label, regressors)
+ design = make_second_level_design_matrix(subjects_label, regressors)
expected_design = np.array([[0.75, 1], [0.1, 1]])
assert_array_equal(design, expected_design)
assert_true(len(design.columns) == 2)
diff --git a/nistats/tests/test_first_level_model.py b/nistats/tests/test_first_level_model.py
index 6a4b0827..91e97aa3 100644
--- a/nistats/tests/test_first_level_model.py
+++ b/nistats/tests/test_first_level_model.py
@@ -13,7 +13,7 @@
from nistats.first_level_model import (mean_scaling, run_glm, FirstLevelModel,
first_level_models_from_bids)
-from nistats.design_matrix import check_design_matrix, make_design_matrix
+from nistats.design_matrix import check_design_matrix, make_first_level_design_matrix
from nose.tools import assert_true, assert_equal, assert_raises
from numpy.testing import (assert_almost_equal, assert_array_equal)
@@ -76,32 +76,32 @@ def test_high_level_glm_one_session():
def test_high_level_glm_with_data():
# New API
- shapes, rk = ((7, 8, 7, 15), (7, 8, 7, 16)), 3
- mask, fmri_data, design_matrices = write_fake_fmri_data(shapes, rk)
-
- multi_session_model = FirstLevelModel(mask=None).fit(
- fmri_data, design_matrices=design_matrices)
- n_voxels = multi_session_model.masker_.mask_img_.get_data().sum()
- z_image = multi_session_model.compute_contrast(np.eye(rk)[1])
- assert_equal(np.sum(z_image.get_data() != 0), n_voxels)
- assert_true(z_image.get_data().std() < 3.)
-
- # with mask
- multi_session_model = FirstLevelModel(mask=mask).fit(
- fmri_data, design_matrices=design_matrices)
- z_image = multi_session_model.compute_contrast(
- np.eye(rk)[:2], output_type='z_score')
- p_value = multi_session_model.compute_contrast(
- np.eye(rk)[:2], output_type='p_value')
- stat_image = multi_session_model.compute_contrast(
- np.eye(rk)[:2], output_type='stat')
- effect_image = multi_session_model.compute_contrast(
- np.eye(rk)[:2], output_type='effect_size')
- variance_image = multi_session_model.compute_contrast(
- np.eye(rk)[:2], output_type='effect_variance')
- assert_array_equal(z_image.get_data() == 0., load(mask).get_data() == 0.)
- assert_true(
- (variance_image.get_data()[load(mask).get_data() > 0] > .001).all())
+ with InTemporaryDirectory():
+ shapes, rk = ((7, 8, 7, 15), (7, 8, 7, 16)), 3
+ mask, fmri_data, design_matrices = write_fake_fmri_data(shapes, rk)
+ multi_session_model = FirstLevelModel(mask=None).fit(
+ fmri_data, design_matrices=design_matrices)
+ n_voxels = multi_session_model.masker_.mask_img_.get_data().sum()
+ z_image = multi_session_model.compute_contrast(np.eye(rk)[1])
+ assert_equal(np.sum(z_image.get_data() != 0), n_voxels)
+ assert_true(z_image.get_data().std() < 3.)
+
+ # with mask
+ multi_session_model = FirstLevelModel(mask=mask).fit(
+ fmri_data, design_matrices=design_matrices)
+ z_image = multi_session_model.compute_contrast(
+ np.eye(rk)[:2], output_type='z_score')
+ p_value = multi_session_model.compute_contrast(
+ np.eye(rk)[:2], output_type='p_value')
+ stat_image = multi_session_model.compute_contrast(
+ np.eye(rk)[:2], output_type='stat')
+ effect_image = multi_session_model.compute_contrast(
+ np.eye(rk)[:2], output_type='effect_size')
+ variance_image = multi_session_model.compute_contrast(
+ np.eye(rk)[:2], output_type='effect_variance')
+ assert_array_equal(z_image.get_data() == 0., load(mask).get_data() == 0.)
+ assert_true(
+ (variance_image.get_data()[load(mask).get_data() > 0] > .001).all())
def test_high_level_glm_with_paths():
@@ -142,7 +142,7 @@ def test_run_glm():
n, p, q = 100, 80, 10
X, Y = np.random.randn(p, q), np.random.randn(p, n)
- # ols case
+ # Ordinary Least Squares case
labels, results = run_glm(Y, X, 'ols')
assert_array_equal(labels, np.zeros(n))
assert_equal(list(results.keys()), [0.0])
@@ -217,9 +217,11 @@ def test_fmri_inputs():
def basic_paradigm():
conditions = ['c0', 'c0', 'c0', 'c1', 'c1', 'c1', 'c2', 'c2', 'c2']
onsets = [30, 70, 100, 10, 30, 90, 30, 40, 60]
- paradigm = pd.DataFrame({'trial_type': conditions,
- 'onset': onsets})
- return paradigm
+ durations = 1 * np.ones(9)
+ events = pd.DataFrame({'trial_type': conditions,
+ 'onset': onsets,
+ 'duration': durations})
+ return events
def test_first_level_model_design_creation():
@@ -230,20 +232,20 @@ def test_first_level_model_design_creation():
FUNCFILE = FUNCFILE[0]
func_img = load(FUNCFILE)
# basic test based on basic_paradigm and glover hrf
- t_r = 1.0
+ t_r = 10.0
slice_time_ref = 0.
- paradigm = basic_paradigm()
+ events = basic_paradigm()
model = FirstLevelModel(t_r, slice_time_ref, mask=mask,
drift_model='polynomial', drift_order=3)
- model = model.fit(func_img, paradigm)
+ model = model.fit(func_img, events)
frame1, X1, names1 = check_design_matrix(model.design_matrices_[0])
# check design computation is identical
n_scans = func_img.get_data().shape[3]
start_time = slice_time_ref * t_r
end_time = (n_scans - 1 + slice_time_ref) * t_r
frame_times = np.linspace(start_time, end_time, n_scans)
- design = make_design_matrix(frame_times, paradigm,
- drift_model='polynomial', drift_order=3)
+ design = make_first_level_design_matrix(frame_times, events,
+ drift_model='polynomial', drift_order=3)
frame2, X2, names2 = check_design_matrix(design)
assert_array_equal(frame1, frame2)
assert_array_equal(X1, X2)
@@ -257,24 +259,42 @@ def test_first_level_model_glm_computation():
FUNCFILE = FUNCFILE[0]
func_img = load(FUNCFILE)
# basic test based on basic_paradigm and glover hrf
- t_r = 1.0
+ t_r = 10.0
slice_time_ref = 0.
- paradigm = basic_paradigm()
- # ols case
+ events = basic_paradigm()
+ # Ordinary Least Squares case
model = FirstLevelModel(t_r, slice_time_ref, mask=mask,
drift_model='polynomial', drift_order=3,
minimize_memory=False)
- model = model.fit(func_img, paradigm)
+ model = model.fit(func_img, events)
labels1 = model.labels_[0]
results1 = model.results_[0]
labels2, results2 = run_glm(
model.masker_.transform(func_img),
- model.design_matrices_[0].as_matrix(), 'ar1')
+ model.design_matrices_[0].values, 'ar1')
# ar not giving consistent results in python 3.4
# assert_almost_equal(labels1, labels2, decimal=2) ####FIX
# assert_equal(len(results1), len(results2)) ####FIX
+def test_first_level_glm_computation_with_memory_caching():
+ with InTemporaryDirectory():
+ shapes = ((7, 8, 9, 10),)
+ mask, FUNCFILE, _ = write_fake_fmri_data(shapes)
+ FUNCFILE = FUNCFILE[0]
+ func_img = load(FUNCFILE)
+ # initialize FirstLevelModel with memory option enabled
+ t_r = 10.0
+ slice_time_ref = 0.
+ events = basic_paradigm()
+ # Ordinary Least Squares case
+ model = FirstLevelModel(t_r, slice_time_ref, mask=mask,
+ drift_model='polynomial', drift_order=3,
+ memory='nilearn_cache', memory_level=1,
+ minimize_memory=False)
+ model.fit(func_img, events)
+
+
def test_first_level_model_contrast_computation():
with InTemporaryDirectory():
shapes = ((7, 8, 9, 10),)
@@ -282,10 +302,10 @@ def test_first_level_model_contrast_computation():
FUNCFILE = FUNCFILE[0]
func_img = load(FUNCFILE)
# basic test based on basic_paradigm and glover hrf
- t_r = 1.0
+ t_r = 10.0
slice_time_ref = 0.
- paradigm = basic_paradigm()
- # ols case
+ events = basic_paradigm()
+ # Ordinary Least Squares case
model = FirstLevelModel(t_r, slice_time_ref, mask=mask,
drift_model='polynomial', drift_order=3,
minimize_memory=False)
@@ -293,7 +313,7 @@ def test_first_level_model_contrast_computation():
# asking for contrast before model fit gives error
assert_raises(ValueError, model.compute_contrast, c1)
# fit model
- model = model.fit([func_img, func_img], [paradigm, paradigm])
+ model = model.fit([func_img, func_img], [events, events])
# smoke test for different contrasts in fixed effects
model.compute_contrast([c1, c2])
# smoke test for same contrast in fixed effects
diff --git a/nistats/tests/test_hemodynamic_models.py b/nistats/tests/test_hemodynamic_models.py
index ba333b7f..d84f1067 100644
--- a/nistats/tests/test_hemodynamic_models.py
+++ b/nistats/tests/test_hemodynamic_models.py
@@ -17,7 +17,7 @@ def test_spm_hrf():
"""
h = spm_hrf(2.0)
assert_almost_equal(h.sum(), 1)
- assert_equal(len(h), 256)
+ assert_equal(len(h), 800)
def test_spm_hrf_derivative():
@@ -25,10 +25,10 @@ def test_spm_hrf_derivative():
"""
h = spm_time_derivative(2.0)
assert_almost_equal(h.sum(), 0)
- assert_equal(len(h), 256)
+ assert_equal(len(h), 800)
h = spm_dispersion_derivative(2.0)
assert_almost_equal(h.sum(), 0)
- assert_equal(len(h), 256)
+ assert_equal(len(h), 800)
def test_glover_hrf():
@@ -36,10 +36,10 @@ def test_glover_hrf():
"""
h = glover_hrf(2.0)
assert_almost_equal(h.sum(), 1)
- assert_equal(len(h), 256)
+ assert_equal(len(h), 800)
h = glover_dispersion_derivative(2.0)
assert_almost_equal(h.sum(), 0)
- assert_equal(len(h), 256)
+ assert_equal(len(h), 800)
@@ -48,7 +48,7 @@ def test_glover_time_derivative():
"""
h = glover_time_derivative(2.0)
assert_almost_equal(h.sum(), 0)
- assert_equal(len(h), 256)
+ assert_equal(len(h), 800)
def test_resample_regressor():
@@ -195,11 +195,11 @@ def test_hkernel():
h = _hrf_kernel('fir', tr, fir_delays=np.arange(4))
assert_equal(len(h), 4)
for dh in h:
- assert_equal(dh.sum(), 16.)
+ assert_equal(dh.sum(), 50.)
#
h = _hrf_kernel(None, tr)
assert_equal(len(h), 1)
- assert_almost_equal(h[0], np.hstack((1, np.zeros(15))))
+ assert_almost_equal(h[0], np.hstack((1, np.zeros(49))))
def test_make_regressor_1():
@@ -221,7 +221,7 @@ def test_make_regressor_2():
frame_times = np.linspace(0, 69, 70)
hrf_model = 'spm'
reg, reg_names = compute_regressor(condition, hrf_model, frame_times)
- assert_almost_equal(reg.sum() * 16, 3, 1)
+ assert_almost_equal(reg.sum() * 50, 3, 1)
assert_equal(reg_names[0], 'cond')
diff --git a/nistats/tests/test_model.py b/nistats/tests/test_model.py
index bf1ff7e2..d08aa3ea 100644
--- a/nistats/tests/test_model.py
+++ b/nistats/tests/test_model.py
@@ -59,11 +59,7 @@ def test_model():
assert_array_almost_equal(RESULTS.theta[1], np.mean(Y))
# Check we get the same as R
assert_array_almost_equal(RESULTS.theta, [1.773, 2.5], 3)
- try:
- percentile = np.percentile
- except AttributeError:
- # Numpy <=1.4.1 does not have percentile function
- raise SkipTest('Numpy does not have percentile function')
+ percentile = np.percentile
pcts = percentile(RESULTS.resid, [0, 25, 50, 75, 100])
assert_array_almost_equal(pcts, [-1.6970, -0.6667, 0, 0.6667, 1.6970], 4)
diff --git a/nistats/tests/test_paradigm.py b/nistats/tests/test_paradigm.py
index 15368742..f9ef45de 100644
--- a/nistats/tests/test_paradigm.py
+++ b/nistats/tests/test_paradigm.py
@@ -9,67 +9,70 @@
import os
import pandas as pd
-from nistats.experimental_paradigm import paradigm_from_csv
-
from nose.tools import assert_true
def basic_paradigm():
conditions = ['c0', 'c0', 'c0', 'c1', 'c1', 'c1', 'c2', 'c2', 'c2']
onsets = [30, 70, 100, 10, 30, 90, 30, 40, 60]
- paradigm = pd.DataFrame({'name': conditions, 'onset': onsets})
- return paradigm
+ durations = 1 * np.ones(9)
+ events = pd.DataFrame({'name': conditions,
+ 'onset': onsets,
+ 'duration': durations})
+ return events
def modulated_block_paradigm():
conditions = ['c0', 'c0', 'c0', 'c1', 'c1', 'c1', 'c2', 'c2', 'c2']
onsets = [30, 70, 100, 10, 30, 90, 30, 40, 60]
- duration = 5 + 5 * np.random.rand(len(onsets))
+ durations = 5 + 5 * np.random.rand(len(onsets))
values = np.random.rand(len(onsets))
- paradigm = pd.DataFrame({'name': conditions,
+ events = pd.DataFrame({'name': conditions,
'onset': onsets,
- 'duration': duration,
+ 'duration': durations,
'modulation': values})
- return paradigm
+ return events
def modulated_event_paradigm():
conditions = ['c0', 'c0', 'c0', 'c1', 'c1', 'c1', 'c2', 'c2', 'c2']
onsets = [30, 70, 100, 10, 30, 90, 30, 40, 60]
+ durations = 1 * np.ones(9)
values = np.random.rand(len(onsets))
- paradigm = pd.DataFrame({'name': conditions,
+ events = pd.DataFrame({'name': conditions,
'onset': onsets,
+ 'durations': durations,
'amplitude': values})
- return paradigm
+ return events
def block_paradigm():
conditions = ['c0', 'c0', 'c0', 'c1', 'c1', 'c1', 'c2', 'c2', 'c2']
onsets = [30, 70, 100, 10, 30, 90, 30, 40, 60]
- duration = 5 * np.ones(9)
- paradigm = pd.DataFrame({'name': conditions,
+ durations = 5 * np.ones(9)
+ events = pd.DataFrame({'name': conditions,
'onset': onsets,
- 'duration': duration})
- return paradigm
+ 'duration': durations})
+ return events
-def write_paradigm(paradigm, tmpdir):
- """Function to write a paradigm to a file and return the address
+def write_events(events, tmpdir):
+ """Function to write events of an experimental paradigm to a file and return the address.
"""
- csvfile = os.path.join(tmpdir, 'paradigm.csv')
- paradigm.to_csv(csvfile)
- return csvfile
+ tsvfile = os.path.join(tmpdir, 'events.tsv')
+ events.to_csv(tsvfile, sep='\t')
+ return tsvfile
-def test_read_paradigm():
- """ test that a paradigm is correctly read
+def test_read_events():
+ """ test that a events for an experimental paradigm are correctly read.
"""
import tempfile
tmpdir = tempfile.mkdtemp()
- for paradigm in (block_paradigm(),
+ for events in (block_paradigm(),
modulated_event_paradigm(),
modulated_block_paradigm(),
basic_paradigm()):
- csvfile = write_paradigm(paradigm, tmpdir)
- read_paradigm = paradigm_from_csv(csvfile)
- assert_true((read_paradigm['onset'] == paradigm['onset']).all())
+ csvfile = write_events(events, tmpdir)
+ read_paradigm = pd.read_table(csvfile)
+ assert_true((read_paradigm['onset'] == events['onset']).all())
diff --git a/nistats/tests/test_regression.py b/nistats/tests/test_regression.py
index 45842da2..de810615 100644
--- a/nistats/tests/test_regression.py
+++ b/nistats/tests/test_regression.py
@@ -4,7 +4,7 @@
import numpy as np
-from nistats.regression import OLSModel, ARModel
+from nistats.regression import OLSModel, ARModel, SimpleRegressionResults
from nose.tools import assert_equal, assert_true
from numpy.testing import assert_array_almost_equal, assert_array_equal
@@ -41,3 +41,13 @@ def test_AR_degenerate():
model = ARModel(design=Xd, rho=0.9)
results = model.fit(Y)
assert_equal(results.df_resid, 31)
+
+def test_simple_results():
+ model = OLSModel(X)
+ results = model.fit(Y)
+ residfull = results.resid
+ predictedfull = results.predicted
+
+ simple_results = SimpleRegressionResults(results)
+ assert_array_equal(results.predicted, simple_results.predicted(X))
+ assert_array_equal(results.resid, simple_results.resid(Y, X))
diff --git a/nistats/tests/test_reporting.py b/nistats/tests/test_reporting.py
index 3bf7c421..b0f970a6 100644
--- a/nistats/tests/test_reporting.py
+++ b/nistats/tests/test_reporting.py
@@ -1,7 +1,13 @@
-from nistats.design_matrix import make_design_matrix
-from nistats.reporting import plot_design_matrix
+import os
+import nibabel as nib
import numpy as np
from numpy.testing import dec
+from nose.tools import assert_true
+from nibabel.tmpdirs import InTemporaryDirectory
+
+from nistats.design_matrix import make_first_level_design_matrix
+from nistats.reporting import (plot_design_matrix, get_clusters_table,
+ _local_max, plot_contrast_matrix)
# Set the backend to avoid having DISPLAY problems
from nilearn.plotting import _set_mpl_backend
@@ -21,7 +27,66 @@
def test_show_design_matrix():
# test that the show code indeed (formally) runs
frame_times = np.linspace(0, 127 * 1., 128)
- DM = make_design_matrix(
+ dmtx = make_first_level_design_matrix(
+ frame_times, drift_model='polynomial', drift_order=3)
+ ax = plot_design_matrix(dmtx)
+ assert (ax is not None)
+ with InTemporaryDirectory():
+ ax = plot_design_matrix(dmtx, output_file='dmtx.png')
+ assert os.path.exists('dmtx.png')
+ assert (ax is None)
+ plot_design_matrix(dmtx, output_file='dmtx.pdf')
+ assert os.path.exists('dmtx.pdf')
+
+@dec.skipif(not have_mpl)
+def test_show_contrast_matrix():
+ # test that the show code indeed (formally) runs
+ frame_times = np.linspace(0, 127 * 1., 128)
+ dmtx = make_first_level_design_matrix(
frame_times, drift_model='polynomial', drift_order=3)
- ax = plot_design_matrix(DM)
+ contrast = np.ones(4)
+ ax = plot_contrast_matrix(contrast, dmtx)
assert (ax is not None)
+ with InTemporaryDirectory():
+ ax = plot_contrast_matrix(contrast, dmtx, output_file='contrast.png')
+ assert os.path.exists('contrast.png')
+ assert (ax is None)
+ plot_contrast_matrix(contrast, dmtx, output_file='contrast.pdf')
+ assert os.path.exists('contrast.pdf')
+
+
+def test_local_max():
+ shape = (9, 10, 11)
+ data = np.zeros(shape)
+ # Two maxima (one global, one local), 10 voxels apart.
+ data[4, 5, :] = [4, 3, 2, 1, 1, 1, 1, 1, 2, 3, 4]
+ data[5, 5, :] = [5, 4, 3, 2, 1, 1, 1, 2, 3, 4, 6]
+ data[6, 5, :] = [4, 3, 2, 1, 1, 1, 1, 1, 2, 3, 4]
+ affine = np.eye(4)
+
+ ijk, vals = _local_max(data, affine, min_distance=9)
+ assert_true(np.array_equal(ijk, np.array([[5., 5., 10.], [5., 5., 0.]])))
+ assert_true(np.array_equal(vals, np.array([6, 5])))
+
+ ijk, vals = _local_max(data, affine, min_distance=11)
+ assert_true(np.array_equal(ijk, np.array([[5., 5., 10.]])))
+ assert_true(np.array_equal(vals, np.array([6])))
+
+
+def test_get_clusters_table():
+ shape = (9, 10, 11)
+ data = np.zeros(shape)
+ data[2:4, 5:7, 6:8] = 5.
+ stat_img = nib.Nifti1Image(data, np.eye(4))
+
+ # test one cluster extracted
+ cluster_table = get_clusters_table(stat_img, 4, 0)
+ assert_true(len(cluster_table) == 1)
+
+ # test empty table on high stat threshold
+ cluster_table = get_clusters_table(stat_img, 6, 0)
+ assert_true(len(cluster_table) == 0)
+
+ # test empty table on high cluster threshold
+ cluster_table = get_clusters_table(stat_img, 4, 9)
+ assert_true(len(cluster_table) == 0)
diff --git a/nistats/tests/test_second_level_model.py b/nistats/tests/test_second_level_model.py
index 4faf5a1f..ad02ded4 100644
--- a/nistats/tests/test_second_level_model.py
+++ b/nistats/tests/test_second_level_model.py
@@ -18,7 +18,7 @@
from numpy.testing import (assert_almost_equal, assert_array_equal)
from nibabel.tmpdirs import InTemporaryDirectory
import pandas as pd
-
+from nilearn.image import concat_imgs
# This directory path
BASEDIR = os.path.dirname(os.path.abspath(__file__))
@@ -46,7 +46,7 @@ def test_high_level_glm_with_paths():
mask, FUNCFILE, _ = write_fake_fmri_data(shapes)
FUNCFILE = FUNCFILE[0]
func_img = load(FUNCFILE)
- # ols case
+ # Ordinary Least Squares case
model = SecondLevelModel(mask=mask)
# asking for contrast before model fit gives error
assert_raises(ValueError, model.compute_contrast, [])
@@ -92,6 +92,7 @@ def test_fmri_inputs():
['03', 'a', FUNCFILE]]
niidf = pd.DataFrame(dfrows, columns=dfcols)
niimgs = [FUNCFILE, FUNCFILE, FUNCFILE]
+ niimg_4d = concat_imgs(niimgs)
confounds = pd.DataFrame([['01', 1], ['02', 2], ['03', 3]],
columns=['subject_label', 'conf1'])
sdes = pd.DataFrame(X[:3, :3], columns=['intercept', 'b', 'c'])
@@ -110,6 +111,9 @@ def test_fmri_inputs():
SecondLevelModel().fit(niidf, None, sdes)
# niimgs as input
SecondLevelModel().fit(niimgs, None, sdes)
+ # 4d niimg as input
+ SecondLevelModel().fit(niimg_4d, None, sdes)
+
# test wrong input errors
# test first level model requirements
assert_raises(ValueError, SecondLevelModel().fit, flm)
@@ -146,7 +150,7 @@ def test_second_level_model_glm_computation():
mask, FUNCFILE, _ = write_fake_fmri_data(shapes)
FUNCFILE = FUNCFILE[0]
func_img = load(FUNCFILE)
- # ols case
+ # Ordinary Least Squares case
model = SecondLevelModel(mask=mask)
Y = [func_img] * 4
X = pd.DataFrame([[1]] * 4, columns=['intercept'])
@@ -157,7 +161,7 @@ def test_second_level_model_glm_computation():
results1 = model.results_
labels2, results2 = run_glm(
- model.masker_.transform(Y), X.as_matrix(), 'ols')
+ model.masker_.transform(Y), X.values, 'ols')
assert_almost_equal(labels1, labels2, decimal=1)
assert_equal(len(results1), len(results2))
@@ -168,7 +172,7 @@ def test_second_level_model_contrast_computation():
mask, FUNCFILE, _ = write_fake_fmri_data(shapes)
FUNCFILE = FUNCFILE[0]
func_img = load(FUNCFILE)
- # ols case
+ # Ordinary Least Squares case
model = SecondLevelModel(mask=mask)
# asking for contrast before model fit gives error
assert_raises(ValueError, model.compute_contrast, 'intercept')
@@ -196,3 +200,28 @@ def test_second_level_model_contrast_computation():
assert_raises(ValueError, model.compute_contrast, c1, None, '')
assert_raises(ValueError, model.compute_contrast, c1, None, [])
assert_raises(ValueError, model.compute_contrast, c1, None, None, '')
+ # check that passing no explicit contrast when the dsign
+ # matrix has morr than one columns raises an error
+ X = pd.DataFrame(np.random.rand(4, 2), columns=['r1', 'r2'])
+ model = model.fit(Y, design_matrix=X)
+ assert_raises(ValueError, model.compute_contrast, None)
+
+
+def test_second_level_model_contrast_computation_with_memory_caching():
+ with InTemporaryDirectory():
+ shapes = ((7, 8, 9, 1),)
+ mask, FUNCFILE, _ = write_fake_fmri_data(shapes)
+ FUNCFILE = FUNCFILE[0]
+ func_img = load(FUNCFILE)
+ # Ordinary Least Squares case
+ model = SecondLevelModel(mask=mask, memory='nilearn_cache')
+ # fit model
+ Y = [func_img] * 4
+ X = pd.DataFrame([[1]] * 4, columns=['intercept'])
+ model = model.fit(Y, design_matrix=X)
+ ncol = len(model.design_matrix_.columns)
+ c1 = np.eye(ncol)[0, :]
+ # test memory caching for compute_contrast
+ model.compute_contrast(c1, output_type='z_score')
+ # or simply pass nothing
+ model.compute_contrast()
diff --git a/nistats/tests/test_thresholding.py b/nistats/tests/test_thresholding.py
index 2a76945c..537d2690 100644
--- a/nistats/tests/test_thresholding.py
+++ b/nistats/tests/test_thresholding.py
@@ -2,11 +2,10 @@
"""
import numpy as np
from scipy.stats import norm
-from nose.tools import assert_true
+from nose.tools import assert_true, assert_raises
from numpy.testing import assert_almost_equal, assert_equal
import nibabel as nib
-from nistats.thresholding import (fdr_threshold, map_threshold,
- get_clusters_table)
+from nistats.thresholding import fdr_threshold, map_threshold
def test_fdr():
@@ -17,20 +16,22 @@ def test_fdr():
np.random.shuffle(x)
assert_almost_equal(fdr_threshold(x, .1), norm.isf(.0005))
assert_true(fdr_threshold(x, .001) == np.infty)
+ assert_raises(ValueError, fdr_threshold, x, -.1)
+ assert_raises(ValueError, fdr_threshold, x, 1.5)
def test_map_threshold():
shape = (9, 10, 11)
p = np.prod(shape)
data = norm.isf(np.linspace(1. / p, 1. - 1. / p, p)).reshape(shape)
- threshold = .001
+ alpha = .001
data[2:4, 5:7, 6:8] = 5.
stat_img = nib.Nifti1Image(data, np.eye(4))
mask_img = nib.Nifti1Image(np.ones(shape), np.eye(4))
# test 1
th_map, _ = map_threshold(
- stat_img, mask_img, threshold, height_control='fpr',
+ stat_img, mask_img, alpha, height_control='fpr',
cluster_threshold=0)
vals = th_map.get_data()
assert_equal(np.sum(vals > 0), 8)
@@ -44,7 +45,7 @@ def test_map_threshold():
# test 3:excessive size threshold
th_map, z_th = map_threshold(
- stat_img, mask_img, threshold, height_control='fpr',
+ stat_img, mask_img, alpha, height_control='fpr',
cluster_threshold=10)
vals = th_map.get_data()
assert_true(np.sum(vals > 0) == 0)
@@ -72,21 +73,25 @@ def test_map_threshold():
vals = th_map.get_data()
assert_equal(np.sum(vals > 0), 8)
+ # test 7 without a map
+ th_map, threshold = map_threshold(
+ None, None, 3.0, height_control=None,
+ cluster_threshold=0)
+ assert_equal(threshold, 3.0)
+ assert_equal(th_map, None)
-def test_get_clusters_table():
- shape = (9, 10, 11)
- data = np.zeros(shape)
- data[2:4, 5:7, 6:8] = 5.
- stat_img = nib.Nifti1Image(data, np.eye(4))
-
- # test one cluster extracted
- cluster_table = get_clusters_table(stat_img, 4, 0)
- assert_true(len(cluster_table) == 1)
+ th_map, threshold = map_threshold(
+ None, None, level=0.05, height_control='fpr',
+ cluster_threshold=0)
+ assert (threshold > 1.64)
+ assert_equal(th_map, None)
- # test empty table on high stat threshold
- cluster_table = get_clusters_table(stat_img, 6, 0)
- assert_true(len(cluster_table) == 0)
+ assert_raises(ValueError, map_threshold, None, None, level=0.05,
+ height_control='fdr')
+ assert_raises(ValueError, map_threshold, None, None, level=0.05,
+ height_control='bonferroni')
- # test empty table on high cluster threshold
- cluster_table = get_clusters_table(stat_img, 4, 9)
- assert_true(len(cluster_table) == 0)
+ # test 8 wrong procedure
+ assert_raises(ValueError, map_threshold, None, None, level=0.05,
+ height_control='plop')
+
diff --git a/nistats/tests/test_utils.py b/nistats/tests/test_utils.py
index 655408e4..fa86b1b6 100644
--- a/nistats/tests/test_utils.py
+++ b/nistats/tests/test_utils.py
@@ -11,8 +11,8 @@
from nibabel.tmpdirs import InTemporaryDirectory
from nose import with_setup
-from nistats.utils import (multiple_mahalanobis, z_score, multiple_fast_inv,
- pos_recipr, full_rank, _check_run_tables,
+from nistats.utils import (multiple_mahalanobis, z_score, multiple_fast_inverse,
+ positive_reciprocal, full_rank, _check_run_tables,
_check_and_load_tables, _check_list_length_match,
get_bids_files, parse_bids_filename,
get_design_from_fslmat)
@@ -69,26 +69,26 @@ def test_multiple_fast_inv():
for i in range(shape[0]):
X[i] = np.dot(X[i], X[i].T)
X_inv_ref[i] = spl.inv(X[i])
- X_inv = multiple_fast_inv(X)
+ X_inv = multiple_fast_inverse(X)
assert_almost_equal(X_inv_ref, X_inv)
def test_pos_recipr():
X = np.array([2, 1, -1, 0], dtype=np.int8)
eX = np.array([0.5, 1, 0, 0])
- Y = pos_recipr(X)
+ Y = positive_reciprocal(X)
yield assert_array_almost_equal, Y, eX
yield assert_equal, Y.dtype.type, np.float64
X2 = X.reshape((2, 2))
- Y2 = pos_recipr(X2)
+ Y2 = positive_reciprocal(X2)
yield assert_array_almost_equal, Y2, eX.reshape((2, 2))
# check that lists have arrived
XL = [0, 1, -1]
- yield assert_array_almost_equal, pos_recipr(XL), [0, 1, 0]
+ yield assert_array_almost_equal, positive_reciprocal(XL), [0, 1, 0]
# scalars
- yield assert_equal, pos_recipr(-1), 0
- yield assert_equal, pos_recipr(0), 0
- yield assert_equal, pos_recipr(2), 0.5
+ yield assert_equal, positive_reciprocal(-1), 0
+ yield assert_equal, positive_reciprocal(0), 0
+ yield assert_equal, positive_reciprocal(2), 0.5
def test_img_table_checks():
@@ -119,9 +119,9 @@ def write_fake_bold_img(file_path, shape, rk=3, affine=np.eye(4)):
def basic_paradigm():
conditions = ['c0', 'c0', 'c0', 'c1', 'c1', 'c1', 'c2', 'c2', 'c2']
onsets = [30, 70, 100, 10, 30, 90, 30, 40, 60]
- paradigm = pd.DataFrame({'trial_type': conditions,
+ events = pd.DataFrame({'trial_type': conditions,
'onset': onsets})
- return paradigm
+ return events
def basic_confounds(length):
diff --git a/nistats/tests/test_verify_events_file_uses_tab_separators.py b/nistats/tests/test_verify_events_file_uses_tab_separators.py
new file mode 100644
index 00000000..1c7be885
--- /dev/null
+++ b/nistats/tests/test_verify_events_file_uses_tab_separators.py
@@ -0,0 +1,137 @@
+import csv
+import os
+from tempfile import NamedTemporaryFile
+
+import pandas as pd
+from nose.tools import (assert_raises,
+ assert_true,
+ )
+
+from nistats.utils import _verify_events_file_uses_tab_separators
+
+
+def make_data_for_test_runs():
+ data_for_temp_datafile = [
+ ['csf', 'constant', 'linearTrend', 'wm'],
+ [13343.032102491035, 1.0, 0.0, 9486.199545677482],
+ [13329.224068063204, 1.0, 1.0, 9497.003324892803],
+ [13291.755627241291, 1.0, 2.0, 9484.012965365506],
+ ]
+
+ delimiters = {
+ 'tab': '\t',
+ 'comma': ',',
+ 'space': ' ',
+ 'semicolon': ';',
+ 'hyphen': '-',
+ }
+
+ return data_for_temp_datafile, delimiters
+
+
+def _create_test_file(temp_csv, test_data, delimiter):
+ csv_writer = csv.writer(temp_csv, delimiter=delimiter)
+ for row in test_data:
+ csv_writer.writerow(row)
+ temp_csv.flush()
+
+
+def _run_test_for_invalid_separator(filepath, delimiter_name):
+ if delimiter_name not in ('tab', 'comma'):
+ with assert_raises(ValueError):
+ _verify_events_file_uses_tab_separators(events_files=filepath)
+ else:
+ result = _verify_events_file_uses_tab_separators(events_files=filepath)
+ assert_true(result is None)
+
+
+def test_for_invalid_separator():
+ data_for_temp_datafile, delimiters = make_data_for_test_runs()
+ for delimiter_name, delimiter_char in delimiters.items():
+ tmp_file_prefix, temp_file_suffix = (
+ 'tmp ', ' ' + delimiter_name + '.csv')
+ with NamedTemporaryFile(mode='w', dir=os.getcwd(),
+ prefix=tmp_file_prefix,
+ suffix=temp_file_suffix) as temp_csv_obj:
+ _create_test_file(temp_csv=temp_csv_obj,
+ test_data=data_for_temp_datafile,
+ delimiter=delimiter_char)
+ _run_test_for_invalid_separator(filepath=temp_csv_obj.name,
+ delimiter_name=delimiter_name)
+
+
+def test_with_2D_dataframe():
+ data_for_pandas_dataframe, _ = make_data_for_test_runs()
+ events_pandas_dataframe = pd.DataFrame(data_for_pandas_dataframe)
+ result = _verify_events_file_uses_tab_separators(
+ events_files=events_pandas_dataframe)
+ assert_true(result is None)
+
+
+def test_with_1D_dataframe():
+ data_for_pandas_dataframe, _ = make_data_for_test_runs()
+ for dataframe_ in data_for_pandas_dataframe:
+ events_pandas_dataframe = pd.DataFrame(dataframe_)
+ result = _verify_events_file_uses_tab_separators(
+ events_files=events_pandas_dataframe)
+ assert_true(result is None)
+
+def test_for_invalid_filepath():
+ filepath = 'junk_file_path.csv'
+ result = _verify_events_file_uses_tab_separators(events_files=filepath)
+ assert_true(result is None)
+
+
+def test_for_pandas_dataframe():
+ events_pandas_dataframe = pd.DataFrame([['a', 'b', 'c'], [0, 1, 2]])
+ result = _verify_events_file_uses_tab_separators(
+ events_files=events_pandas_dataframe)
+ assert_true(result is None)
+
+
+def test_binary_opening_an_image():
+ img_data = bytearray(
+ b'GIF87a\x01\x00\x01\x00\xe7*\x00\x00\x00\x00\x01\x01\x01\x02\x02'
+ b'\x07\x08\x08\x08\t\t\t\n\n\n\x0b\x0b\x0b\x0c\x0c\x0c\r;')
+ with NamedTemporaryFile(mode='wb', suffix='.gif',
+ dir=os.getcwd()) as temp_img_obj:
+ temp_img_obj.write(img_data)
+ with assert_raises(ValueError):
+ _verify_events_file_uses_tab_separators(
+ events_files=temp_img_obj.name)
+
+
+def test_binary_bytearray_of_ints_data():
+ temp_data_bytearray_from_ints = bytearray([0, 1, 0, 11, 10])
+ with NamedTemporaryFile(mode='wb', dir=os.getcwd(),
+ suffix='.bin') as temp_bin_obj:
+ temp_bin_obj.write(temp_data_bytearray_from_ints)
+ with assert_raises(ValueError):
+ result = _verify_events_file_uses_tab_separators(
+ events_files=temp_bin_obj.name)
+
+
+if __name__ == '__main__':
+
+ def _run_tests_print_test_messages(test_func):
+ from pprint import pprint
+ pprint(['Running', test_func.__name__])
+ test_func()
+ pprint('... complete')
+
+
+ def run_test_suite():
+ tests = [
+ test_for_invalid_filepath,
+ test_with_2D_dataframe,
+ test_with_1D_dataframe,
+ test_for_invalid_filepath,
+ test_for_pandas_dataframe,
+ test_binary_opening_an_image,
+ test_binary_bytearray_of_ints_data,
+ ]
+ for test_ in tests:
+ _run_tests_print_test_messages(test_func=test_)
+
+
+ run_test_suite()
diff --git a/nistats/thresholding.py b/nistats/thresholding.py
index 4bdae008..cb4c8df8 100644
--- a/nistats/thresholding.py
+++ b/nistats/thresholding.py
@@ -7,54 +7,101 @@
from scipy.ndimage import label
from scipy.stats import norm
from nilearn.input_data import NiftiMasker
-from nilearn.image.resampling import coord_transform
-import pandas as pd
-from warnings import warn
def fdr_threshold(z_vals, alpha):
- """ return the BH fdr for the input z_vals"""
+ """ return the Benjamini-Hochberg FDR threshold for the input z_vals
+
+ Parameters
+ ----------
+ z_vals: array,
+ a set of z-variates from which the FDR is computed
+ alpha: float,
+ desired FDR control
+
+ Returns
+ -------
+ threshold: float,
+ FDR-controling threshold from the Benjamini-Hochberg procedure
+ """
+ if alpha < 0 or alpha > 1:
+ raise ValueError('alpha should be between 0 and 1')
z_vals_ = - np.sort(- z_vals)
p_vals = norm.sf(z_vals_)
n_samples = len(p_vals)
pos = p_vals < alpha * np.linspace(
.5 / n_samples, 1 - .5 / n_samples, n_samples)
if pos.any():
- return (z_vals_[pos][-1] - 1.e-8)
+ return (z_vals_[pos][-1] - 1.e-12)
else:
return np.infty
-def map_threshold(stat_img, mask_img=None, threshold=.001,
+def map_threshold(stat_img=None, mask_img=None, level=.001,
height_control='fpr', cluster_threshold=0):
- """ Threshold the provided map
+ """ Compute the required threshold level and return the thresholded map
Parameters
----------
- stat_img : Niimg-like object,
+ stat_img : Niimg-like object or None, optional
statistical image (presumably in z scale)
+ whenever height_control is 'fpr' or None,
+ stat_img=None is acceptable.
+ If it is 'fdr' or 'bonferroni', an error is raised if stat_img is None.
mask_img : Niimg-like object, optional,
mask image
- threshold: float, optional
- cluster forming threshold (either a p-value or z-scale value)
+ level: float, optional
+ number controling the thresholding (either a p-value or z-scale value).
+ Not to be confused with the z-scale threshold: level can be a p-values,
+ e.g. "0.05" or another type of number depending on the
+ height_control parameter. The z-scale threshold is actually returned by
+ the function.
- height_control: string, optional
+ height_control: string, or None optional
false positive control meaning of cluster forming
- threshold: 'fpr'|'fdr'|'bonferroni'|'none'
+ threshold: 'fpr'|'fdr'|'bonferroni'\|None
cluster_threshold : float, optional
- cluster size threshold
+ cluster size threshold. In the returned thresholded map,
+ sets of connected voxels (`clusters`) with size smaller
+ than this number will be removed.
Returns
-------
thresholded_map : Nifti1Image,
- the stat_map theresholded at the prescribed voxel- and cluster-level
+ the stat_map thresholded at the prescribed voxel- and cluster-level
threshold: float,
the voxel-level threshold used actually
+
+ Note
+ ----
+ If the input image is not z-scaled (i.e. some z-transformed statistic)
+ the computed threshold is not rigorous and likely meaningless
"""
+ # Check that height_control is correctly specified
+ if height_control not in ['fpr', 'fdr', 'bonferroni', None]:
+ raise ValueError(
+ "height control should be one of ['fpr', 'fdr', 'bonferroni', None]")
+
+ # if height_control is 'fpr' or None, we don't need to look at the data
+ # to compute the threhsold
+ if height_control == 'fpr':
+ threshold = norm.isf(level)
+ elif height_control is None:
+ threshold = level
+
+ # In this case, and is stat_img is None, we return
+ if stat_img is None:
+ if height_control in ['fpr', None]:
+ return None, threshold
+ else:
+ raise ValueError(
+ 'Map_threshold requires stat_img not to be None'
+ 'when the heigh_control procedure is bonferroni or fdr')
+
# Masking
if mask_img is None:
masker = NiftiMasker(mask_strategy='background').fit(stat_img)
@@ -64,105 +111,21 @@ def map_threshold(stat_img, mask_img=None, threshold=.001,
n_voxels = np.size(stats)
# Thresholding
- if height_control == 'fpr':
- z_th = norm.isf(threshold)
- elif height_control == 'fdr':
- z_th = fdr_threshold(stats, threshold)
+ if height_control == 'fdr':
+ threshold = fdr_threshold(stats, level)
elif height_control == 'bonferroni':
- z_th = norm.isf(threshold / n_voxels)
- else: # Brute-force thresholding
- z_th = threshold
- stats *= (stats > z_th)
+ threshold = norm.isf(level / n_voxels)
+ stats *= (stats > threshold)
# embed it back to 3D grid
stat_map = masker.inverse_transform(stats).get_data()
# Extract connected components above threshold
- label_map, n_labels = label(stat_map > z_th)
+ label_map, n_labels = label(stat_map > threshold)
labels = label_map[masker.mask_img_.get_data() > 0]
for label_ in range(1, n_labels + 1):
if np.sum(labels == label_) < cluster_threshold:
stats[labels == label_] = 0
- return masker.inverse_transform(stats), z_th
-
-
-def get_clusters_table(stat_img, stat_threshold, cluster_threshold):
- """Creates pandas dataframe with img cluster statistics.
-
- Parameters
- ----------
- stat_img : Niimg-like object,
- statistical image (presumably in z scale)
-
- stat_threshold: float, optional
- cluster forming threshold (either a p-value or z-scale value)
-
- cluster_threshold : int, optional
- cluster size threshold
-
- Returns
- -------
- Pandas dataframe with img clusters
- """
-
- stat_map = stat_img.get_data()
-
- # If the stat threshold is too high simply return an empty dataframe
- if np.sum(stat_map > stat_threshold) == 0:
- warn('Attention: No clusters with stat higher than %f' %
- stat_threshold)
- return pd.DataFrame()
-
- # Extract connected components above threshold
- label_map, n_labels = label(stat_map > stat_threshold)
-
- for label_ in range(1, n_labels + 1):
- if np.sum(label_map == label_) < cluster_threshold:
- stat_map[label_map == label_] = 0
-
- # If the cluster threshold is too high simply return an empty dataframe
- # this checks for stats higher than threshold after small clusters
- # were removed from stat_map
- if np.sum(stat_map > stat_threshold) == 0:
- warn('Attention: No clusters with more than %d voxels' %
- cluster_threshold)
- return pd.DataFrame()
-
- label_map, n_labels = label(stat_map > stat_threshold)
- label_map = np.ravel(label_map)
- stat_map = np.ravel(stat_map)
-
- peaks = []
- max_stat = []
- clusters_size = []
- coords = []
- for label_ in range(1, n_labels + 1):
- cluster = stat_map.copy()
- cluster[label_map != label_] = 0
-
- peak = np.unravel_index(np.argmax(cluster),
- stat_img.get_data().shape)
- peaks.append(peak)
-
- max_stat.append(np.max(cluster))
-
- clusters_size.append(np.sum(label_map == label_))
-
- x_map, y_map, z_map = peak
- mni_coords = np.asarray(
- coord_transform(
- x_map, y_map, z_map, stat_img.get_affine())).tolist()
- mni_coords = [round(x) for x in mni_coords]
- coords.append(mni_coords)
-
- vx, vy, vz = zip(*peaks)
- x, y, z = zip(*coords)
-
- columns = ['Vx', 'Vy', 'Vz', 'X', 'Y', 'Z', 'Peak stat', 'Cluster size']
- clusters_table = pd.DataFrame(
- list(zip(vx, vy, vz, x, y, z, max_stat, clusters_size)),
- columns=columns)
-
- return clusters_table
+ return masker.inverse_transform(stats), threshold
diff --git a/nistats/utils.py b/nistats/utils.py
index 6668d2a3..b6ff672d 100644
--- a/nistats/utils.py
+++ b/nistats/utils.py
@@ -2,6 +2,7 @@
Authors: Bertrand Thirion, Matthew Brett, 2015
"""
+import csv
import sys
import scipy.linalg as spl
import numpy as np
@@ -22,15 +23,39 @@ def _check_list_length_match(list_1, list_2, var_name_1, var_name_2):
% (str(var_name_1), len(list_1), str(var_name_2), len(list_2)))
+def _read_events_table(table):
+ """
+ Accepts the path to en event.tsv file and loads it as a Pandas Dataframe.
+ Raises an error if loading fails.
+ Parameters
+ ----------
+ table: string
+ Accepts the path to an events file
+
+ Returns
+ -------
+ loaded: pandas.Dataframe object
+ Pandas Dataframe witht e events data.
+ """
+ try:
+ # kept for historical reasons, a lot of tests use csv with index column
+ loaded = pd.read_csv(table, index_col=0)
+ except:
+ raise ValueError('table path %s could not be loaded' % table)
+ if loaded.empty:
+ try:
+ loaded = pd.read_table(table)
+ except:
+ raise ValueError('table path %s could not be loaded' % table)
+ return loaded
+
+
def _check_and_load_tables(tables_, var_name):
"""Check tables can be loaded in DataFrame to raise error if necessary"""
tables = []
for table_idx, table in enumerate(tables_):
if isinstance(table, _basestring):
- try:
- loaded = pd.read_csv(table, index_col=0)
- except:
- raise ValueError('table path %s could not be loaded' % table)
+ loaded = _read_events_table(table)
tables.append(loaded)
elif isinstance(table, pd.DataFrame):
tables.append(table)
@@ -40,6 +65,61 @@ def _check_and_load_tables(tables_, var_name):
(var_name, type(table), table_idx))
return tables
+
+def _verify_events_file_uses_tab_separators(events_files):
+ """
+ Raises a ValueError if provided list of text based data files
+ (.csv, .tsv, etc) do not enforce the BIDS convention of using Tabs
+ as separators.
+
+ Only scans their first row.
+ Does nothing if:
+ If the separator used is BIDS compliant.
+ Paths are invalid.
+ File(s) are not text files.
+
+ Does not flag comma-separated-values-files for compatibility reasons;
+ this may change in future as commas are not BIDS compliant.
+
+ parameters
+ ----------
+ events_files: str, List/Tuple[str]
+ A single file's path or a collection of filepaths.
+ Files are expected to be text files.
+ Non-text files will raise ValueError.
+
+ Returns
+ -------
+ None
+
+ Raises
+ ------
+ ValueError:
+ If value separators are not Tabs (or commas)
+ """
+ valid_separators = [',', '\t']
+ events_files = [events_files] if not isinstance(events_files, (list, tuple)) else events_files
+ for events_file_ in events_files:
+ try:
+ with open(events_file_, 'r') as events_file_obj:
+ events_file_sample = events_file_obj.readline()
+ except TypeError as type_err: # events is Pandas dataframe.
+ pass
+ except UnicodeDecodeError as unicode_err: # py3:if binary file
+ pass
+ except IOError as io_err: # if invalid filepath.
+ pass
+ else:
+ try:
+ csv.Sniffer().sniff(sample=events_file_sample,
+ delimiters=valid_separators,
+ )
+ except csv.Error:
+ raise ValueError(
+ 'The values in the events file are not separated by tabs; '
+ 'please enforce BIDS conventions',
+ events_file_)
+
def _check_run_tables(run_imgs, tables_, tables_name):
"""Check fMRI runs and corresponding tables to raise error if necessary"""
@@ -57,7 +137,7 @@ def z_score(pvalue):
return norm.isf(pvalue)
-def multiple_fast_inv(a):
+def multiple_fast_inverse(a):
"""Compute the inverse of a set of arrays.
Parameters
@@ -146,7 +226,7 @@ def multiple_mahalanobis(effect, covariance):
Xt, Kt = np.ascontiguousarray(effect.T), np.ascontiguousarray(covariance.T)
# compute the inverse of the covariances
- Kt = multiple_fast_inv(Kt)
+ Kt = multiple_fast_inverse(Kt)
# derive the squared Mahalanobis distances
sqd = np.sum(np.sum(Xt[:, :, np.newaxis] * Xt[:, np.newaxis] * Kt, 1), 1)
@@ -185,7 +265,7 @@ def full_rank(X, cmax=1e15):
return X, cmax
-def pos_recipr(X):
+def positive_reciprocal(X):
""" Return element-wise reciprocal of array, setting `X`>=0 to 0
Return the reciprocal of an array, setting all entries less than or
@@ -314,13 +394,14 @@ def parse_bids_filename(img_path):
'file_tag', 'file_type' and 'file_fields'.
The 'file_tag' field refers to the last part of the file under the
- BIDS convention that is of the form *_tag.type. Contrary to the rest
+ BIDS convention that is of the form \*_tag.type. Contrary to the rest
of the file name it is not a key-value pair. This notion should be
revised in the case we are handling derivatives since so far the
convention will keep the tag prepended to any fields added in the
case of preprocessed files that also end with another tag. This parser
will consider any tag in the middle of the file name as a key with no
value and will be included in the 'file_fields' key.
+
"""
reference = {}
reference['file_path'] = img_path
diff --git a/nistats/version.py b/nistats/version.py
index d15f977b..5b8d3c44 100644
--- a/nistats/version.py
+++ b/nistats/version.py
@@ -21,7 +21,7 @@
# Dev branch marker is: 'X.Y.dev' or 'X.Y.devN' where N is an integer.
# 'X.Y.dev0' is the canonical version of 'X.Y.dev'
#
-__version__ = '0.0.1a'
+__version__ = '0.0.1b'
_NISTATS_INSTALL_MSG = 'See %s for installation information.' % (
'http://nistats.github.io/introduction.html#installation')
@@ -30,32 +30,32 @@
# in some meaningful order (more => less 'core').
REQUIRED_MODULE_METADATA = (
('numpy', {
- 'min_version': '1.8.2',
+ 'min_version': '1.11',
'install_info': _NISTATS_INSTALL_MSG}),
('scipy', {
- 'min_version': '0.14',
+ 'min_version': '0.17',
'install_info': _NISTATS_INSTALL_MSG}),
('sklearn', {
'pypi_name': 'scikit-learn',
- 'min_version': '0.15.0',
+ 'min_version': '0.18',
'install_info': _NISTATS_INSTALL_MSG}),
('nibabel', {
'min_version': '2.0.2',
'required_at_installation': True,
'install_info': _NISTATS_INSTALL_MSG}),
('nilearn', {
- 'min_version': '0.2.0',
+ 'min_version': '0.4',
'install_info': _NISTATS_INSTALL_MSG}),
('pandas', {
- 'min_version': '0.13.0',
+ 'min_version': '0.18.0',
'install_info': _NISTATS_INSTALL_MSG}),
('patsy', {
- 'min_version': '0.2.0',
+ 'min_version': '0.4.1',
'install_info': _NISTATS_INSTALL_MSG}),
)
-OPTIONAL_MATPLOTLIB_MIN_VERSION = '1.3.1'
-OPTIONAL_BOTO3_MIN_VERSION = '1.0.0'
+OPTIONAL_MATPLOTLIB_MIN_VERSION = '1.5.1'
+OPTIONAL_BOTO3_MIN_VERSION = '1.4.0'
def _import_module_with_version_check(