From 2b7df92bddbf022f736e0c0fb452bda0a50b994d Mon Sep 17 00:00:00 2001 From: Gilles de Hollander Date: Fri, 7 Dec 2018 17:40:05 +0100 Subject: [PATCH 01/32] Allow for storing residuals in First Level Models --- nistats/first_level_model.py | 44 ++++++++++++++++++++++ nistats/regression.py | 49 +++++++++---------------- nistats/tests/test_first_level_model.py | 38 +++++++++++++++++++ nistats/tests/test_regression.py | 31 +++++++++++++++- 4 files changed, 129 insertions(+), 33 deletions(-) diff --git a/nistats/first_level_model.py b/nistats/first_level_model.py index cc47160a..1482e75a 100644 --- a/nistats/first_level_model.py +++ b/nistats/first_level_model.py @@ -18,6 +18,7 @@ import numpy as np import pandas as pd from nibabel import Nifti1Image +from nibabel.onetime import setattr_on_read from sklearn.base import (BaseEstimator, clone, @@ -309,6 +310,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 @@ -583,6 +585,48 @@ def compute_contrast(self, contrast_def, stat_type=None, return outputs if output_type == 'all' else output + def get_voxelwise_model_attribute_(self, attribute, timeseries=True): + + if self.minimize_memory: + raise ValueError('minimize_memory should be set to False to make residuals or predictions.') + + if self.labels_ is None or self.results_ is None: + raise ValueError('The model has not been fit yet') + + output = [] + + for design_matrix, labels, results in zip(self.design_matrices_, self.labels_, self.results_): + + if timeseries: + voxelwise_attribute = np.zeros((design_matrix.shape[0], len(labels))) + else: + voxelwise_attribute = np.zeros((1, len(labels))) + + for label_ in results: + label_mask = labels == label_ + voxelwise_attribute[:, label_mask] = getattr(results[label_], attribute) + + + output.append(self.masker_.inverse_transform(voxelwise_attribute)) + + if len(output) == 1: + return output[0] + else: + return output + + @setattr_on_read + def residuals(self): + return self.get_voxelwise_model_attribute_('resid') + + @setattr_on_read + def predicted(self): + return self.get_voxelwise_model_attribute_('predicted') + + @setattr_on_read + def rsq(self): + return self.get_voxelwise_model_attribute_('rsq', timeseries=False) + + @replace_parameters({'mask': 'mask_img'}, end_version='next') def first_level_models_from_bids( diff --git a/nistats/regression.py b/nistats/regression.py index 670e7676..39554051 100644 --- a/nistats/regression.py +++ b/nistats/regression.py @@ -61,7 +61,6 @@ class OLSModel(object): df_resid : scalar Degrees of freedom of the residuals. Number of observations less the rank of the design. - df_model : scalar Degrees of freedome of the model. The rank of the design. """ @@ -276,6 +275,7 @@ 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 def resid(self): @@ -310,7 +310,7 @@ def predicted(self): """ beta = self.theta # the LikelihoodModelResults has parameters named 'theta' - X = self.model.design + X = self.wdesign return np.dot(X, beta) @setattr_on_read @@ -319,13 +319,20 @@ def SSE(self): """ return (self.wresid ** 2).sum(0) + @setattr_on_read + def rsq(self): + """Proportion of explained variance. + If not from an OLS model this is "pseudo"-R2. + """ + return np.var(self.predicted, 0) / np.var(self.wY, 0) + @setattr_on_read 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. @@ -347,41 +354,19 @@ 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 ValueError('minimize_memory should be set to False to make residuals or predictions.') - def resid(self, Y): + def resid(self): """ Residuals from the fit. """ - return Y - self.predicted - - def norm_resid(self, Y): - """ - Residuals, normalized to have unit length. - - Notes - ----- - Is this supposed to return "stanardized residuals," - residuals standardized - to have mean zero and approximately unit variance? + raise ValueError('minimize_memory should be set to False to make residuals or predictions.') - d_i = e_i / sqrt(MS_E) - - Where MS_E = SSE / (n - k) - - See: Montgomery and Peck 3.2.1 p. 68 - Davidson and MacKinnon 15.2 p 662 - """ - return self.resid(Y) * positive_reciprocal(np.sqrt(self.dispersion)) - - def predicted(self): - """ 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) + def norm_resid(self): + raise ValueError('minimize_memory should be set to False to make residuals or predictions.') diff --git a/nistats/tests/test_first_level_model.py b/nistats/tests/test_first_level_model.py index 7a12176b..2dfb3b35 100644 --- a/nistats/tests/test_first_level_model.py +++ b/nistats/tests/test_first_level_model.py @@ -523,3 +523,41 @@ def test_param_mask_deprecation_first_level_models_from_bids(): for param_warning_ in raised_param_deprecation_warnings: assert str(param_warning_.message) == deprecation_msg assert param_warning_.category is DeprecationWarning + + +def test_first_level_model_residuals(): + shapes = [(10, 10, 10, 100)] * 3 + mask, fmri_data, design_matrices = generate_fake_fmri_data(shapes) + + for i in range(len(design_matrices)): + design_matrices[i].iloc[:, 0] = 1 + + model = FirstLevelModel(mask=mask, minimize_memory=False, + noise_model='ols') + model.fit(fmri_data, design_matrices=design_matrices) + + for resid in model.residuals: + mean_resids = model.masker_.transform(resid).mean(0) + assert_array_almost_equal(mean_resids, 0) + + +def test_first_level_model_predictions_rsq(): + shapes = [(10, 10, 10, 25)] * 3 + mask, fmri_data, design_matrices = generate_fake_fmri_data(shapes, 25) + + model = FirstLevelModel(mask=mask, + signal_scaling=False, + minimize_memory=False, + noise_model='ols') + model.fit(fmri_data, design_matrices=design_matrices) + + for pred, data, rsq in zip(model.predicted, fmri_data, model.rsq): + y_predicted = model.masker_.transform(pred) + y_measured = model.masker_.transform(data) + + assert_array_almost_equal(y_predicted, y_measured) + + rsq = model.masker_.transform(rsq) + assert_array_almost_equal(rsq, 1.0) + + diff --git a/nistats/tests/test_regression.py b/nistats/tests/test_regression.py index ff25e289..75e05241 100644 --- a/nistats/tests/test_regression.py +++ b/nistats/tests/test_regression.py @@ -6,7 +6,8 @@ from nose.tools import assert_equal -from nistats.regression import OLSModel, ARModel +from nose.tools import assert_equal, assert_true, assert_almost_equal +from numpy.testing import assert_array_almost_equal, assert_array_equal RNG = np.random.RandomState(20110902) @@ -18,12 +19,40 @@ def test_OLS(): model = OLSModel(design=X) results = model.fit(Y) assert_equal(results.df_resid, 30) + assert_equal(results.resid.shape[0], 40) + assert_equal(results.predicted.shape[0], 40) def test_AR(): model = ARModel(design=X, rho=0.4) results = model.fit(Y) assert_equal(results.df_resid, 30) + assert_equal(results.resid.shape[0], 40) + assert_equal(results.predicted.shape[0], 40) + +def test_residuals(): + Xintercept = X.copy() + + # If design matrix contains an intercept, the + # mean of the residuals should be 0 (short of + # some numerical rounding errors) + X[:, 0] = 1 + model = OLSModel(design=X) + results = model.fit(Y) + assert_almost_equal(results.resid.mean(), 0) + +def test_predicted_rsq(): + Xshort = X.copy()[:10, :] + Yshort = Y.copy()[:10] + + # Signal of 10 elements should be completely + # predicted by 10 predictors (short of some numerical + # rounding errors) + model = OLSModel(design=Xshort) + results = model.fit(Yshort) + assert_almost_equal(results.resid.sum(), 0) + assert_array_almost_equal(results.predicted, Yshort) + assert_almost_equal(results.rsq, 1.0) def test_OLS_degenerate(): From 0663e7ae507900294e3f7a8f5fa1287b143d06b2 Mon Sep 17 00:00:00 2001 From: James Kent Date: Tue, 19 Nov 2019 21:22:56 -0600 Subject: [PATCH 02/32] import missing modules and update calls to functions --- nistats/tests/test_first_level_model.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/nistats/tests/test_first_level_model.py b/nistats/tests/test_first_level_model.py index 2dfb3b35..f3c4a8f5 100644 --- a/nistats/tests/test_first_level_model.py +++ b/nistats/tests/test_first_level_model.py @@ -18,6 +18,7 @@ assert_true, ) from numpy.testing import (assert_almost_equal, + assert_array_almost_equal, assert_array_equal, ) from nibabel.tmpdirs import InTemporaryDirectory @@ -526,8 +527,8 @@ def test_param_mask_deprecation_first_level_models_from_bids(): def test_first_level_model_residuals(): - shapes = [(10, 10, 10, 100)] * 3 - mask, fmri_data, design_matrices = generate_fake_fmri_data(shapes) + shapes, rk = [(10, 10, 10, 100)], 3 + mask, fmri_data, design_matrices = _generate_fake_fmri_data(shapes, rk) for i in range(len(design_matrices)): design_matrices[i].iloc[:, 0] = 1 @@ -542,8 +543,8 @@ def test_first_level_model_residuals(): def test_first_level_model_predictions_rsq(): - shapes = [(10, 10, 10, 25)] * 3 - mask, fmri_data, design_matrices = generate_fake_fmri_data(shapes, 25) + shapes, rk = [(10, 10, 10, 25)], 3 + mask, fmri_data, design_matrices = _generate_fake_fmri_data(shapes, rk) model = FirstLevelModel(mask=mask, signal_scaling=False, From 1b58ebbe93aa360cb1439c54d29e49cc8843dd6f Mon Sep 17 00:00:00 2001 From: James Kent Date: Tue, 19 Nov 2019 21:30:08 -0600 Subject: [PATCH 03/32] fix flake8 errors --- nistats/first_level_model.py | 25 ++++++++++++------------- nistats/regression.py | 2 +- nistats/tests/test_first_level_model.py | 6 ++---- nistats/tests/test_regression.py | 7 +++++-- 4 files changed, 20 insertions(+), 20 deletions(-) diff --git a/nistats/first_level_model.py b/nistats/first_level_model.py index 1482e75a..733deb55 100644 --- a/nistats/first_level_model.py +++ b/nistats/first_level_model.py @@ -118,14 +118,14 @@ def run_glm(Y, X, noise_model='ar1', bins=100, n_jobs=1, verbose=0): acceptable_noise_models = ['ar1', 'ols'] if noise_model not in acceptable_noise_models: raise ValueError( - "Acceptable noise models are {0}. You provided 'noise_model={1}'".\ - format(acceptable_noise_models, noise_model)) + "Acceptable noise models are {0}. You provided 'noise_model={1}'". + format(acceptable_noise_models, noise_model)) if Y.shape[0] != X.shape[0]: raise ValueError( 'The number of rows of Y should match the number of rows of X.' - ' You provided X with shape {0} and Y with shape {1}'.\ - format(X.shape, Y.shape)) + ' You provided X with shape {0} and Y with shape {1}'. + format(X.shape, Y.shape)) # Create the model ols_result = OLSModel(X).fit(Y) @@ -586,27 +586,26 @@ def compute_contrast(self, contrast_def, stat_type=None, return outputs if output_type == 'all' else output def get_voxelwise_model_attribute_(self, attribute, timeseries=True): - + if self.minimize_memory: raise ValueError('minimize_memory should be set to False to make residuals or predictions.') - + if self.labels_ is None or self.results_ is None: raise ValueError('The model has not been fit yet') - + output = [] - + for design_matrix, labels, results in zip(self.design_matrices_, self.labels_, self.results_): - + if timeseries: voxelwise_attribute = np.zeros((design_matrix.shape[0], len(labels))) else: voxelwise_attribute = np.zeros((1, len(labels))) - + for label_ in results: - label_mask = labels == label_ + label_mask = labels == label_ voxelwise_attribute[:, label_mask] = getattr(results[label_], attribute) - - + output.append(self.masker_.inverse_transform(voxelwise_attribute)) if len(output) == 1: diff --git a/nistats/regression.py b/nistats/regression.py index 39554051..a881a8c1 100644 --- a/nistats/regression.py +++ b/nistats/regression.py @@ -321,7 +321,7 @@ def SSE(self): @setattr_on_read def rsq(self): - """Proportion of explained variance. + """Proportion of explained variance. If not from an OLS model this is "pseudo"-R2. """ return np.var(self.predicted, 0) / np.var(self.wY, 0) diff --git a/nistats/tests/test_first_level_model.py b/nistats/tests/test_first_level_model.py index f3c4a8f5..cec5e947 100644 --- a/nistats/tests/test_first_level_model.py +++ b/nistats/tests/test_first_level_model.py @@ -528,7 +528,7 @@ def test_param_mask_deprecation_first_level_models_from_bids(): def test_first_level_model_residuals(): shapes, rk = [(10, 10, 10, 100)], 3 - mask, fmri_data, design_matrices = _generate_fake_fmri_data(shapes, rk) + mask, fmri_data, design_matrices = _generate_fake_fmri_data(shapes, rk) for i in range(len(design_matrices)): design_matrices[i].iloc[:, 0] = 1 @@ -544,7 +544,7 @@ def test_first_level_model_residuals(): def test_first_level_model_predictions_rsq(): shapes, rk = [(10, 10, 10, 25)], 3 - mask, fmri_data, design_matrices = _generate_fake_fmri_data(shapes, rk) + mask, fmri_data, design_matrices = _generate_fake_fmri_data(shapes, rk) model = FirstLevelModel(mask=mask, signal_scaling=False, @@ -560,5 +560,3 @@ def test_first_level_model_predictions_rsq(): rsq = model.masker_.transform(rsq) assert_array_almost_equal(rsq, 1.0) - - diff --git a/nistats/tests/test_regression.py b/nistats/tests/test_regression.py index 75e05241..fcc38c9d 100644 --- a/nistats/tests/test_regression.py +++ b/nistats/tests/test_regression.py @@ -9,6 +9,7 @@ from nose.tools import assert_equal, assert_true, assert_almost_equal from numpy.testing import assert_array_almost_equal, assert_array_equal +from nistats.regression import OLSModel, ARModel RNG = np.random.RandomState(20110902) X = RNG.standard_normal((40, 10)) @@ -30,21 +31,23 @@ def test_AR(): assert_equal(results.resid.shape[0], 40) assert_equal(results.predicted.shape[0], 40) + def test_residuals(): Xintercept = X.copy() # If design matrix contains an intercept, the # mean of the residuals should be 0 (short of # some numerical rounding errors) - X[:, 0] = 1 + Xintercept[:, 0] = 1 model = OLSModel(design=X) results = model.fit(Y) assert_almost_equal(results.resid.mean(), 0) + def test_predicted_rsq(): Xshort = X.copy()[:10, :] Yshort = Y.copy()[:10] - + # Signal of 10 elements should be completely # predicted by 10 predictors (short of some numerical # rounding errors) From 4909febe191871658290a943e21441ca66df5c52 Mon Sep 17 00:00:00 2001 From: James Kent Date: Tue, 19 Nov 2019 21:33:05 -0600 Subject: [PATCH 04/32] remove @set_attr_on_read --- nistats/first_level_model.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/nistats/first_level_model.py b/nistats/first_level_model.py index 733deb55..fe89f803 100644 --- a/nistats/first_level_model.py +++ b/nistats/first_level_model.py @@ -613,15 +613,12 @@ def get_voxelwise_model_attribute_(self, attribute, timeseries=True): else: return output - @setattr_on_read def residuals(self): return self.get_voxelwise_model_attribute_('resid') - @setattr_on_read def predicted(self): return self.get_voxelwise_model_attribute_('predicted') - @setattr_on_read def rsq(self): return self.get_voxelwise_model_attribute_('rsq', timeseries=False) From acac09b594cb70af407ca837a0da8a19b23b6966 Mon Sep 17 00:00:00 2001 From: James Kent Date: Wed, 20 Nov 2019 00:37:41 -0600 Subject: [PATCH 05/32] modify tests to reflect updated code base --- nistats/tests/test_first_level_model.py | 22 ++++++++++++++-------- nistats/tests/test_regression.py | 2 +- 2 files changed, 15 insertions(+), 9 deletions(-) diff --git a/nistats/tests/test_first_level_model.py b/nistats/tests/test_first_level_model.py index cec5e947..415ca630 100644 --- a/nistats/tests/test_first_level_model.py +++ b/nistats/tests/test_first_level_model.py @@ -18,8 +18,8 @@ assert_true, ) from numpy.testing import (assert_almost_equal, - assert_array_almost_equal, assert_array_equal, + assert_array_less, ) from nibabel.tmpdirs import InTemporaryDirectory @@ -537,7 +537,7 @@ def test_first_level_model_residuals(): noise_model='ols') model.fit(fmri_data, design_matrices=design_matrices) - for resid in model.residuals: + for resid in model.residuals(): mean_resids = model.masker_.transform(resid).mean(0) assert_array_almost_equal(mean_resids, 0) @@ -546,17 +546,23 @@ def test_first_level_model_predictions_rsq(): shapes, rk = [(10, 10, 10, 25)], 3 mask, fmri_data, design_matrices = _generate_fake_fmri_data(shapes, rk) + for i in range(len(design_matrices)): + design_matrices[i].iloc[:, 0] = 1 + model = FirstLevelModel(mask=mask, signal_scaling=False, minimize_memory=False, noise_model='ols') model.fit(fmri_data, design_matrices=design_matrices) - for pred, data, rsq in zip(model.predicted, fmri_data, model.rsq): - y_predicted = model.masker_.transform(pred) - y_measured = model.masker_.transform(data) + pred = model.predicted() + data = fmri_data[0] + rsq = model.rsq() + + y_predicted = model.masker_.transform(pred) + y_measured = model.masker_.transform(data) - assert_array_almost_equal(y_predicted, y_measured) + assert_almost_equal(np.mean(y_predicted - y_measured), 0) - rsq = model.masker_.transform(rsq) - assert_array_almost_equal(rsq, 1.0) + rsq = model.masker_.transform(rsq) + assert_array_lessl(0., rsq) diff --git a/nistats/tests/test_regression.py b/nistats/tests/test_regression.py index fcc38c9d..9721f325 100644 --- a/nistats/tests/test_regression.py +++ b/nistats/tests/test_regression.py @@ -39,7 +39,7 @@ def test_residuals(): # mean of the residuals should be 0 (short of # some numerical rounding errors) Xintercept[:, 0] = 1 - model = OLSModel(design=X) + model = OLSModel(design=Xintercept) results = model.fit(Y) assert_almost_equal(results.resid.mean(), 0) From 4e7c897696e635ddf3e6515b6a78bc7191bdeddd Mon Sep 17 00:00:00 2001 From: James Kent Date: Wed, 20 Nov 2019 00:45:26 -0600 Subject: [PATCH 06/32] fix typo and simplifiy loop --- nistats/tests/test_first_level_model.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/nistats/tests/test_first_level_model.py b/nistats/tests/test_first_level_model.py index 415ca630..31d8476e 100644 --- a/nistats/tests/test_first_level_model.py +++ b/nistats/tests/test_first_level_model.py @@ -18,6 +18,7 @@ assert_true, ) from numpy.testing import (assert_almost_equal, + assert_array_almost_equal, assert_array_equal, assert_array_less, ) @@ -537,9 +538,9 @@ def test_first_level_model_residuals(): noise_model='ols') model.fit(fmri_data, design_matrices=design_matrices) - for resid in model.residuals(): - mean_resids = model.masker_.transform(resid).mean(0) - assert_array_almost_equal(mean_resids, 0) + resid = model.residuals() + mean_resids = model.masker_.transform(resid).mean(0) + assert_array_almost_equal(mean_resids, 0) def test_first_level_model_predictions_rsq(): @@ -565,4 +566,4 @@ def test_first_level_model_predictions_rsq(): assert_almost_equal(np.mean(y_predicted - y_measured), 0) rsq = model.masker_.transform(rsq) - assert_array_lessl(0., rsq) + assert_array_less(0., rsq) From 0b15ec66a878406f2775f5830513eeeb649009fa Mon Sep 17 00:00:00 2001 From: James Kent Date: Mon, 2 Dec 2019 12:15:58 -0600 Subject: [PATCH 07/32] respond to review comments: - add docstrings - re-add @setattr_on_read - change rsq to r_square --- nistats/first_level_model.py | 51 ++++++++++++++++++++++++++++++++++-- 1 file changed, 49 insertions(+), 2 deletions(-) diff --git a/nistats/first_level_model.py b/nistats/first_level_model.py index fe89f803..db3e0183 100644 --- a/nistats/first_level_model.py +++ b/nistats/first_level_model.py @@ -586,7 +586,25 @@ def compute_contrast(self, contrast_def, stat_type=None, return outputs if output_type == 'all' else output def get_voxelwise_model_attribute_(self, attribute, timeseries=True): + """Transform RegressionResults instances within a dictionary + (whose keys represent the autoregressive coefficient under the 'ar1' + noise model or only 0.0 under 'ols' noise_model and values are the + RegressionResults instances) into input nifti space. + Parameters + ---------- + attribute : str + an attribute of a RegressionResults instance + timeseries : bool, optional + whether the RegressionResult attribute has a value + per timepoint of the input nifti image. + + Returns + ------- + output : list or Nifti1Image + a list of Nifti1Images if FirstLevelModel is fit with + a list of Nifti1Images, or a single Nifti1Image otherwise. + """ if self.minimize_memory: raise ValueError('minimize_memory should be set to False to make residuals or predictions.') @@ -613,15 +631,44 @@ def get_voxelwise_model_attribute_(self, attribute, timeseries=True): else: return output + @setattr_on_read def residuals(self): + """Transform voxelwise residuals to the same shape + as the input Nifti1Image(s) + + Returns + ------- + output : list or Nifti1Image + a list of Nifti1Images if FirstLevelModel is fit with + a list of Nifti1Images, or a single Nifti1Image otherwise. + """ return self.get_voxelwise_model_attribute_('resid') + @setattr_on_read def predicted(self): + """Transform voxelwise predicted values to the same shape + as the input Nifti1Image(s) + + Returns + ------- + output : list or Nifti1Image + a list of Nifti1Images if FirstLevelModel is fit with + a list of Nifti1Images, or a single Nifti1Image otherwise. + """ return self.get_voxelwise_model_attribute_('predicted') - def rsq(self): - return self.get_voxelwise_model_attribute_('rsq', timeseries=False) + @setattr_on_read + def r_square(self): + """Transform voxelwise r-squared values to the same shape + as the input Nifti1Image(s) + Returns + ------- + output : list or Nifti1Image + a list of Nifti1Images if FirstLevelModel is fit with + a list of Nifti1Images, or a single Nifti1Image otherwise. + """ + return self.get_voxelwise_model_attribute_('r_square', timeseries=False) @replace_parameters({'mask': 'mask_img'}, end_version='next') From 5f260c9c0530b57271f8e4de82a98110f351c7a2 Mon Sep 17 00:00:00 2001 From: James Kent Date: Mon, 2 Dec 2019 12:16:24 -0600 Subject: [PATCH 08/32] change rsq to r_square --- nistats/regression.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nistats/regression.py b/nistats/regression.py index a881a8c1..ed3e57d3 100644 --- a/nistats/regression.py +++ b/nistats/regression.py @@ -320,7 +320,7 @@ def SSE(self): return (self.wresid ** 2).sum(0) @setattr_on_read - def rsq(self): + def r_square(self): """Proportion of explained variance. If not from an OLS model this is "pseudo"-R2. """ From adb504cff1653c22f4b0f67c902e8f428078bd35 Mon Sep 17 00:00:00 2001 From: James Kent Date: Mon, 2 Dec 2019 12:16:58 -0600 Subject: [PATCH 09/32] change rsq to r_square in tests --- nistats/tests/test_first_level_model.py | 7 +++---- nistats/tests/test_regression.py | 4 ++-- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/nistats/tests/test_first_level_model.py b/nistats/tests/test_first_level_model.py index 31d8476e..8d16b0a2 100644 --- a/nistats/tests/test_first_level_model.py +++ b/nistats/tests/test_first_level_model.py @@ -543,7 +543,7 @@ def test_first_level_model_residuals(): assert_array_almost_equal(mean_resids, 0) -def test_first_level_model_predictions_rsq(): +def test_first_level_model_predictions_r_square(): shapes, rk = [(10, 10, 10, 25)], 3 mask, fmri_data, design_matrices = _generate_fake_fmri_data(shapes, rk) @@ -558,12 +558,11 @@ def test_first_level_model_predictions_rsq(): pred = model.predicted() data = fmri_data[0] - rsq = model.rsq() + r_square = model.r_square() y_predicted = model.masker_.transform(pred) y_measured = model.masker_.transform(data) assert_almost_equal(np.mean(y_predicted - y_measured), 0) - rsq = model.masker_.transform(rsq) - assert_array_less(0., rsq) + assert_array_less(0., r_square) diff --git a/nistats/tests/test_regression.py b/nistats/tests/test_regression.py index 9721f325..81801446 100644 --- a/nistats/tests/test_regression.py +++ b/nistats/tests/test_regression.py @@ -44,7 +44,7 @@ def test_residuals(): assert_almost_equal(results.resid.mean(), 0) -def test_predicted_rsq(): +def test_predicted_r_square(): Xshort = X.copy()[:10, :] Yshort = Y.copy()[:10] @@ -55,7 +55,7 @@ def test_predicted_rsq(): results = model.fit(Yshort) assert_almost_equal(results.resid.sum(), 0) assert_array_almost_equal(results.predicted, Yshort) - assert_almost_equal(results.rsq, 1.0) + assert_almost_equal(results.r_square, 1.0) def test_OLS_degenerate(): From 9b434bc0814ff441d4a619ad93b3eb6ef582ea66 Mon Sep 17 00:00:00 2001 From: James Kent Date: Mon, 2 Dec 2019 12:27:07 -0600 Subject: [PATCH 10/32] fix function calls --- nistats/tests/test_first_level_model.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/nistats/tests/test_first_level_model.py b/nistats/tests/test_first_level_model.py index 8d16b0a2..25657729 100644 --- a/nistats/tests/test_first_level_model.py +++ b/nistats/tests/test_first_level_model.py @@ -538,7 +538,7 @@ def test_first_level_model_residuals(): noise_model='ols') model.fit(fmri_data, design_matrices=design_matrices) - resid = model.residuals() + resid = model.residuals mean_resids = model.masker_.transform(resid).mean(0) assert_array_almost_equal(mean_resids, 0) @@ -556,13 +556,14 @@ def test_first_level_model_predictions_r_square(): noise_model='ols') model.fit(fmri_data, design_matrices=design_matrices) - pred = model.predicted() + pred = model.predicted data = fmri_data[0] - r_square = model.r_square() + r_square_3d = model.r_square y_predicted = model.masker_.transform(pred) y_measured = model.masker_.transform(data) assert_almost_equal(np.mean(y_predicted - y_measured), 0) - assert_array_less(0., r_square) + r_square_2d = model.masker_.transform(r_square_3d) + assert_array_less(0., r_square_2d) From d416d0bfb515a7a2c2216b8586b4fb45ee64ddb7 Mon Sep 17 00:00:00 2001 From: Gilles de Hollander Date: Wed, 18 Dec 2019 15:13:04 +0100 Subject: [PATCH 11/32] Example of how to use for . --- .../plot_predictions_residuals.py | 110 ++++++++++++++++++ 1 file changed, 110 insertions(+) create mode 100644 examples/02_first_level_models/plot_predictions_residuals.py diff --git a/examples/02_first_level_models/plot_predictions_residuals.py b/examples/02_first_level_models/plot_predictions_residuals.py new file mode 100644 index 00000000..500f25c8 --- /dev/null +++ b/examples/02_first_level_models/plot_predictions_residuals.py @@ -0,0 +1,110 @@ +""" +Predicted time series and residuals +=================================== + +Here we fit a First Level GLM with the `minimize_memory`-argument set to `False`. +By doing so, the `FirstLevelModel`-object stores the residuals, which we can then inspect. +Also, the predicted time series can be extracted, which are useful to assess the quality of the model fit. +""" + + +######################################################################### +# Import modules +# ------------------------------------- +from nistats.datasets import fetch_spm_auditory +from nilearn.plotting import plot_stat_map +from nilearn import input_data, image +import matplotlib.pyplot as plt +from nilearn import plotting, masking +from nistats.reporting import get_clusters_table +from nistats.first_level_model import FirstLevelModel +import pandas as pd + +# load fMRI data +subject_data = fetch_spm_auditory() +fmri_img = image.concat_imgs(subject_data.func) + +# Make an average +mean_img = image.mean_img(fmri_img) +mask = masking.compute_epi_mask(mean_img) + +# Clean (demean) and smooth data +fmri_img = image.clean_img(fmri_img) +fmri_img = image.smooth_img(fmri_img, 5.) + +# load events +events = pd.read_table(subject_data['events']) + + +######################################################################### +# Fit model +# ------------------------------------- +# Note that `minimize_memory` is set to `False` so that `FirstLevelModel` +# stores the residuals. +# `signal_scaling` is set to False, so we keep the same scaling as the +# original data in `fmri_img`. + +fmri_glm = FirstLevelModel(t_r=7, + drift_model='cosine', + signal_scaling=False, + mask_img=mask, + minimize_memory=False) + +fmri_glm = fmri_glm.fit(fmri_img, events) + + +######################################################################### +# Calculate and plot contrast +# ------------------------------------- +z_map = fmri_glm.compute_contrast('active - rest') + +plotting.plot_stat_map(z_map, bg_img=mean_img, threshold=3.1) + +######################################################################### +# Extract the largest clusters +# ------------------------------------- + +table = get_clusters_table(z_map, stat_threshold=3.1, + cluster_threshold=20).set_index('Cluster ID', drop=True) +table.head() + +masker = input_data.NiftiSpheresMasker(table.loc[range(1,7), ['X', 'Y', 'Z']].values) + +real_timeseries = masker.fit_transform(fmri_img) +predicted_timeseries = masker.fit_transform(fmri_glm.predicted) + + +######################################################################### +# Plot predicted and actual time series for 6 most significant clusters +# ------------------------------------- +for i in range(1, 7): + plt.subplot(2, 3, i) + plt.title('Cluster peak {}\n'.format(table.loc[i, ['X', 'Y', 'Z']].tolist())) + plt.plot(real_timeseries[:, i-1], c='k', lw=2) + plt.plot(predicted_timeseries[:, i-1], c='r', ls='--', lw=2) + plt.xlabel('Time') + plt.ylabel('Signal intensity') + +plt.gcf().set_size_inches(12, 7) +plt.tight_layout() + +######################################################################### +# Get residuals +# ------------------------------------- +resid = masker.fit_transform(fmri_glm.residuals) + + +######################################################################### +# Plot distribution of residuals +# ------------------------------------- +# Note that residuals are not really distributed normally. + + +for i in range(1, 7): + plt.subplot(2, 3, i) + plt.title('Cluster peak {}\n'.format(table.loc[i, ['X', 'Y', 'Z']].tolist())) + plt.hist(resid[:, i-1]) + print('Mean residuals: {}'.format(resid[:, i-1].mean())) + +plt.gcf().set_size_inches(12, 7) +plt.tight_layout() From 65010cfaf9f3d46462a91bf19a0a4b207b3d62ef Mon Sep 17 00:00:00 2001 From: Gilles de Hollander Date: Wed, 18 Dec 2019 15:26:11 +0100 Subject: [PATCH 12/32] Made ValueError for storing model attributes more verbose. --- nistats/first_level_model.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/nistats/first_level_model.py b/nistats/first_level_model.py index db3e0183..2bb410cc 100644 --- a/nistats/first_level_model.py +++ b/nistats/first_level_model.py @@ -606,7 +606,10 @@ def get_voxelwise_model_attribute_(self, attribute, timeseries=True): a list of Nifti1Images, or a single Nifti1Image otherwise. """ if self.minimize_memory: - raise ValueError('minimize_memory should be set to False to make residuals or predictions.') + raise ValueError('To access voxelwise attributes like R-squared, residuals, ' + 'and predictions, the `FirstLevelModel`-object needs to store ' + 'there attributes. To do so, set `minimize_memory` to `False` ' + 'when initializing the `FirstLevelModel`-object.') if self.labels_ is None or self.results_ is None: raise ValueError('The model has not been fit yet') From 0b5a170317fbabcc4e3c24dc0b0c9fcc9be76193 Mon Sep 17 00:00:00 2001 From: Gilles de Hollander Date: Wed, 18 Dec 2019 15:30:25 +0100 Subject: [PATCH 13/32] Also include R-squared --- .../plot_predictions_residuals.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/examples/02_first_level_models/plot_predictions_residuals.py b/examples/02_first_level_models/plot_predictions_residuals.py index 500f25c8..94e66343 100644 --- a/examples/02_first_level_models/plot_predictions_residuals.py +++ b/examples/02_first_level_models/plot_predictions_residuals.py @@ -108,3 +108,14 @@ plt.gcf().set_size_inches(12, 7) plt.tight_layout() + + +######################################################################### +# Plot R-squared +# ------------------------------------- +# Because we stored the residuals, we can plot the R-squared: the proportion +# of explained variance of the GLM as a whole. Note that the R-squared is markedly +# lower deep down the brain, where there is more physiological noise and we +# are further away from the receive coils. +plotting.plot_stat_map(fmri_glm.r_square, + bg_img=mean_img, threshold=.1, display_mode='z', cut_coords=7) From 5b61be4378ebe8f0deeba317f2e51367a35f0618 Mon Sep 17 00:00:00 2001 From: James Kent Date: Wed, 18 Dec 2019 12:47:51 -0600 Subject: [PATCH 14/32] fix heading underlines in example --- .../plot_predictions_residuals.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/examples/02_first_level_models/plot_predictions_residuals.py b/examples/02_first_level_models/plot_predictions_residuals.py index 94e66343..ad55097b 100644 --- a/examples/02_first_level_models/plot_predictions_residuals.py +++ b/examples/02_first_level_models/plot_predictions_residuals.py @@ -10,7 +10,7 @@ ######################################################################### # Import modules -# ------------------------------------- +# -------------- from nistats.datasets import fetch_spm_auditory from nilearn.plotting import plot_stat_map from nilearn import input_data, image @@ -38,7 +38,7 @@ ######################################################################### # Fit model -# ------------------------------------- +# --------- # Note that `minimize_memory` is set to `False` so that `FirstLevelModel` # stores the residuals. # `signal_scaling` is set to False, so we keep the same scaling as the @@ -55,14 +55,14 @@ ######################################################################### # Calculate and plot contrast -# ------------------------------------- +# --------------------------- z_map = fmri_glm.compute_contrast('active - rest') plotting.plot_stat_map(z_map, bg_img=mean_img, threshold=3.1) ######################################################################### # Extract the largest clusters -# ------------------------------------- +# ---------------------------- table = get_clusters_table(z_map, stat_threshold=3.1, cluster_threshold=20).set_index('Cluster ID', drop=True) @@ -76,7 +76,7 @@ ######################################################################### # Plot predicted and actual time series for 6 most significant clusters -# ------------------------------------- +# --------------------------------------------------------------------- for i in range(1, 7): plt.subplot(2, 3, i) plt.title('Cluster peak {}\n'.format(table.loc[i, ['X', 'Y', 'Z']].tolist())) @@ -90,13 +90,13 @@ ######################################################################### # Get residuals -# ------------------------------------- +# ------------- resid = masker.fit_transform(fmri_glm.residuals) ######################################################################### # Plot distribution of residuals -# ------------------------------------- +# ------------------------------ # Note that residuals are not really distributed normally. @@ -112,7 +112,7 @@ ######################################################################### # Plot R-squared -# ------------------------------------- +# -------------- # Because we stored the residuals, we can plot the R-squared: the proportion # of explained variance of the GLM as a whole. Note that the R-squared is markedly # lower deep down the brain, where there is more physiological noise and we From fcf8320d7c8a788ff2e7d79e4763fe996a652cbf Mon Sep 17 00:00:00 2001 From: James Kent Date: Thu, 19 Dec 2019 11:40:24 -0600 Subject: [PATCH 15/32] fix grammar Co-Authored-By: bthirion --- examples/02_first_level_models/plot_predictions_residuals.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/02_first_level_models/plot_predictions_residuals.py b/examples/02_first_level_models/plot_predictions_residuals.py index ad55097b..0220fb22 100644 --- a/examples/02_first_level_models/plot_predictions_residuals.py +++ b/examples/02_first_level_models/plot_predictions_residuals.py @@ -4,7 +4,7 @@ Here we fit a First Level GLM with the `minimize_memory`-argument set to `False`. By doing so, the `FirstLevelModel`-object stores the residuals, which we can then inspect. -Also, the predicted time series can be extracted, which are useful to assess the quality of the model fit. +Also, the predicted time series can be extracted, which is useful to assess the quality of the model fit. """ From a86491907cacd229465b6f5aa569088ceef493d5 Mon Sep 17 00:00:00 2001 From: James Kent Date: Thu, 19 Dec 2019 11:51:39 -0600 Subject: [PATCH 16/32] fix code formatting and do not standardize --- .../plot_predictions_residuals.py | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/examples/02_first_level_models/plot_predictions_residuals.py b/examples/02_first_level_models/plot_predictions_residuals.py index ad55097b..db464d08 100644 --- a/examples/02_first_level_models/plot_predictions_residuals.py +++ b/examples/02_first_level_models/plot_predictions_residuals.py @@ -1,7 +1,7 @@ """ Predicted time series and residuals =================================== - + Here we fit a First Level GLM with the `minimize_memory`-argument set to `False`. By doing so, the `FirstLevelModel`-object stores the residuals, which we can then inspect. Also, the predicted time series can be extracted, which are useful to assess the quality of the model fit. @@ -12,7 +12,6 @@ # Import modules # -------------- from nistats.datasets import fetch_spm_auditory -from nilearn.plotting import plot_stat_map from nilearn import input_data, image import matplotlib.pyplot as plt from nilearn import plotting, masking @@ -28,8 +27,8 @@ mean_img = image.mean_img(fmri_img) mask = masking.compute_epi_mask(mean_img) -# Clean (demean) and smooth data -fmri_img = image.clean_img(fmri_img) +# Clean and smooth data +fmri_img = image.clean_img(fmri_img, standardize=False) fmri_img = image.smooth_img(fmri_img, 5.) # load events @@ -65,10 +64,10 @@ # ---------------------------- table = get_clusters_table(z_map, stat_threshold=3.1, - cluster_threshold=20).set_index('Cluster ID', drop=True) + cluster_threshold=20).set_index('Cluster ID', drop=True) table.head() -masker = input_data.NiftiSpheresMasker(table.loc[range(1,7), ['X', 'Y', 'Z']].values) +masker = input_data.NiftiSpheresMasker(table.loc[range(1, 7), ['X', 'Y', 'Z']].values) real_timeseries = masker.fit_transform(fmri_img) predicted_timeseries = masker.fit_transform(fmri_glm.predicted) @@ -78,13 +77,13 @@ # Plot predicted and actual time series for 6 most significant clusters # --------------------------------------------------------------------- for i in range(1, 7): - plt.subplot(2, 3, i) + plt.subplot(2, 3, i) plt.title('Cluster peak {}\n'.format(table.loc[i, ['X', 'Y', 'Z']].tolist())) plt.plot(real_timeseries[:, i-1], c='k', lw=2) plt.plot(predicted_timeseries[:, i-1], c='r', ls='--', lw=2) plt.xlabel('Time') plt.ylabel('Signal intensity') - + plt.gcf().set_size_inches(12, 7) plt.tight_layout() @@ -101,7 +100,7 @@ for i in range(1, 7): - plt.subplot(2, 3, i) + plt.subplot(2, 3, i) plt.title('Cluster peak {}\n'.format(table.loc[i, ['X', 'Y', 'Z']].tolist())) plt.hist(resid[:, i-1]) print('Mean residuals: {}'.format(resid[:, i-1].mean())) From 806f11d9e43e3c941c9049f180f72fb1ff071976 Mon Sep 17 00:00:00 2001 From: James Kent Date: Thu, 19 Dec 2019 11:54:42 -0600 Subject: [PATCH 17/32] change parameter timeseries to result_as_time_series --- nistats/first_level_model.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/nistats/first_level_model.py b/nistats/first_level_model.py index 38705354..4b6098ce 100644 --- a/nistats/first_level_model.py +++ b/nistats/first_level_model.py @@ -585,7 +585,7 @@ def compute_contrast(self, contrast_def, stat_type=None, return outputs if output_type == 'all' else output - def get_voxelwise_model_attribute_(self, attribute, timeseries=True): + def get_voxelwise_model_attribute_(self, attribute, result_as_time_series=True): """Transform RegressionResults instances within a dictionary (whose keys represent the autoregressive coefficient under the 'ar1' noise model or only 0.0 under 'ols' noise_model and values are the @@ -595,7 +595,7 @@ def get_voxelwise_model_attribute_(self, attribute, timeseries=True): ---------- attribute : str an attribute of a RegressionResults instance - timeseries : bool, optional + result_as_time_series : bool, optional whether the RegressionResult attribute has a value per timepoint of the input nifti image. @@ -618,7 +618,7 @@ def get_voxelwise_model_attribute_(self, attribute, timeseries=True): for design_matrix, labels, results in zip(self.design_matrices_, self.labels_, self.results_): - if timeseries: + if result_as_time_series: voxelwise_attribute = np.zeros((design_matrix.shape[0], len(labels))) else: voxelwise_attribute = np.zeros((1, len(labels))) @@ -671,7 +671,7 @@ def r_square(self): a list of Nifti1Images if FirstLevelModel is fit with a list of Nifti1Images, or a single Nifti1Image otherwise. """ - return self.get_voxelwise_model_attribute_('r_square', timeseries=False) + return self.get_voxelwise_model_attribute_('r_square', result_as_time_series=False) @replace_parameters({'mask': 'mask_img'}, end_version='next') From a47586fc1326c86578fcf780768395f64b42c9f1 Mon Sep 17 00:00:00 2001 From: James Kent Date: Wed, 8 Jan 2020 12:07:13 -0600 Subject: [PATCH 18/32] attempt to address @bthirion comments --- .../plot_predictions_residuals.py | 31 ++++++++++++++++++- 1 file changed, 30 insertions(+), 1 deletion(-) diff --git a/examples/02_first_level_models/plot_predictions_residuals.py b/examples/02_first_level_models/plot_predictions_residuals.py index 75a51778..f60d7cc1 100644 --- a/examples/02_first_level_models/plot_predictions_residuals.py +++ b/examples/02_first_level_models/plot_predictions_residuals.py @@ -18,6 +18,8 @@ from nistats.reporting import get_clusters_table from nistats.first_level_model import FirstLevelModel import pandas as pd +import numpy as np + # load fMRI data subject_data = fetch_spm_auditory() @@ -115,6 +117,33 @@ # Because we stored the residuals, we can plot the R-squared: the proportion # of explained variance of the GLM as a whole. Note that the R-squared is markedly # lower deep down the brain, where there is more physiological noise and we -# are further away from the receive coils. +# are further away from the receive coils. However, R-Squared should be interpreted +# with a grain of salt. The R-squared value will necessarily increase with +# the addition of more factors (such as rest, active, drift, motion) into the GLM. +# Additionally, we are looking at the overall fit of the model, so we are +# unable to say whether a voxel/region has a large R-squared value because +# the voxel/region is responsive to the experiment (such as active or rest) +# or because the voxel/region fits the noise factors (such as drift or motion) +# that could be present in the GLM. To isolate the influence of the experiment, +# we can use an F-test as shown in the next section. plotting.plot_stat_map(fmri_glm.r_square, bg_img=mean_img, threshold=.1, display_mode='z', cut_coords=7) + + +######################################################################### +# Calculate and Plot F-test +# ------------------------- +# The F-test tells you how well the GLM fits effects of interest such as +# the active and rest conditions together. This is different from R-squared, +# which tells you how well the overall GLM fits the data, including active, rest +# and all the other columns in the design matrix such as drift and motion. +design_matrix = fmri_glm.design_matrices_[0] +active = np.array([1 if c == 'active' else 0 for c in design_matrix.columns]) +rest = np.array([1 if c == 'rest' else 0 for c in design_matrix.columns]) +effects_of_interest = np.vstack((active, rest)) + +z_map = fmri_glm.compute_contrast(effects_of_interest, + output_type='z_score') + +plotting.plot_stat_map(z_map, + bg_img=mean_img, threshold=2.33, display_mode='z', cut_coords=7) From cdf21612f060fa9ce22003ecbbf900cf5d27bd4c Mon Sep 17 00:00:00 2001 From: James Kent Date: Wed, 8 Jan 2020 16:43:19 -0600 Subject: [PATCH 19/32] split imports statements --- .../plot_predictions_residuals.py | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/examples/02_first_level_models/plot_predictions_residuals.py b/examples/02_first_level_models/plot_predictions_residuals.py index f60d7cc1..a32a2749 100644 --- a/examples/02_first_level_models/plot_predictions_residuals.py +++ b/examples/02_first_level_models/plot_predictions_residuals.py @@ -12,13 +12,9 @@ # Import modules # -------------- from nistats.datasets import fetch_spm_auditory -from nilearn import input_data, image -import matplotlib.pyplot as plt -from nilearn import plotting, masking -from nistats.reporting import get_clusters_table -from nistats.first_level_model import FirstLevelModel +from nilearn import image +from nilearn import masking import pandas as pd -import numpy as np # load fMRI data @@ -44,6 +40,7 @@ # stores the residuals. # `signal_scaling` is set to False, so we keep the same scaling as the # original data in `fmri_img`. +from nistats.first_level_model import FirstLevelModel fmri_glm = FirstLevelModel(t_r=7, drift_model='cosine', @@ -57,6 +54,8 @@ ######################################################################### # Calculate and plot contrast # --------------------------- +from nilearn import plotting + z_map = fmri_glm.compute_contrast('active - rest') plotting.plot_stat_map(z_map, bg_img=mean_img, threshold=3.1) @@ -64,6 +63,8 @@ ######################################################################### # Extract the largest clusters # ---------------------------- +from nistats.reporting import get_clusters_table +from nilearn import input_data table = get_clusters_table(z_map, stat_threshold=3.1, cluster_threshold=20).set_index('Cluster ID', drop=True) @@ -78,6 +79,8 @@ ######################################################################### # Plot predicted and actual time series for 6 most significant clusters # --------------------------------------------------------------------- +import matplotlib.pyplot as plt + for i in range(1, 7): plt.subplot(2, 3, i) plt.title('Cluster peak {}\n'.format(table.loc[i, ['X', 'Y', 'Z']].tolist())) @@ -137,6 +140,8 @@ # the active and rest conditions together. This is different from R-squared, # which tells you how well the overall GLM fits the data, including active, rest # and all the other columns in the design matrix such as drift and motion. +import numpy as np + design_matrix = fmri_glm.design_matrices_[0] active = np.array([1 if c == 'active' else 0 for c in design_matrix.columns]) rest = np.array([1 if c == 'rest' else 0 for c in design_matrix.columns]) From b8c483dbc96361b6a100a4835b9eae19918674ea Mon Sep 17 00:00:00 2001 From: James Kent Date: Wed, 8 Jan 2020 16:44:37 -0600 Subject: [PATCH 20/32] always return list get_voxelwise_model_attribute_ --- nistats/first_level_model.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/nistats/first_level_model.py b/nistats/first_level_model.py index 4b6098ce..aab5195f 100644 --- a/nistats/first_level_model.py +++ b/nistats/first_level_model.py @@ -629,9 +629,6 @@ def get_voxelwise_model_attribute_(self, attribute, result_as_time_series=True): output.append(self.masker_.inverse_transform(voxelwise_attribute)) - if len(output) == 1: - return output[0] - else: return output @setattr_on_read From 72aec1dd2e5b8bbc0abcda3875f925f500bf2aa7 Mon Sep 17 00:00:00 2001 From: James Kent Date: Wed, 8 Jan 2020 19:47:53 -0600 Subject: [PATCH 21/32] change docstrings for output to always be a list --- nistats/first_level_model.py | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/nistats/first_level_model.py b/nistats/first_level_model.py index aab5195f..be25f5c6 100644 --- a/nistats/first_level_model.py +++ b/nistats/first_level_model.py @@ -601,9 +601,8 @@ def get_voxelwise_model_attribute_(self, attribute, result_as_time_series=True): Returns ------- - output : list or Nifti1Image - a list of Nifti1Images if FirstLevelModel is fit with - a list of Nifti1Images, or a single Nifti1Image otherwise. + output : list + a list of Nifti1Image(s) """ if self.minimize_memory: raise ValueError('To access voxelwise attributes like R-squared, residuals, ' @@ -638,9 +637,8 @@ def residuals(self): Returns ------- - output : list or Nifti1Image - a list of Nifti1Images if FirstLevelModel is fit with - a list of Nifti1Images, or a single Nifti1Image otherwise. + output : list + a list of Nifti1Image(s) """ return self.get_voxelwise_model_attribute_('resid') @@ -651,9 +649,8 @@ def predicted(self): Returns ------- - output : list or Nifti1Image - a list of Nifti1Images if FirstLevelModel is fit with - a list of Nifti1Images, or a single Nifti1Image otherwise. + output : list + a list of Nifti1Image(s) """ return self.get_voxelwise_model_attribute_('predicted') @@ -664,9 +661,8 @@ def r_square(self): Returns ------- - output : list or Nifti1Image - a list of Nifti1Images if FirstLevelModel is fit with - a list of Nifti1Images, or a single Nifti1Image otherwise. + output : list + a list of Nifti1Image(s) """ return self.get_voxelwise_model_attribute_('r_square', result_as_time_series=False) From a245635ca875fdc96607ff29d6062303c68bd2c6 Mon Sep 17 00:00:00 2001 From: James Kent Date: Wed, 8 Jan 2020 20:13:05 -0600 Subject: [PATCH 22/32] modify tests to treat output as list --- .../plot_predictions_residuals.py | 67 ++++++++++++------- nistats/tests/test_first_level_model.py | 6 +- 2 files changed, 44 insertions(+), 29 deletions(-) diff --git a/examples/02_first_level_models/plot_predictions_residuals.py b/examples/02_first_level_models/plot_predictions_residuals.py index a32a2749..3410d3a7 100644 --- a/examples/02_first_level_models/plot_predictions_residuals.py +++ b/examples/02_first_level_models/plot_predictions_residuals.py @@ -70,10 +70,13 @@ cluster_threshold=20).set_index('Cluster ID', drop=True) table.head() -masker = input_data.NiftiSpheresMasker(table.loc[range(1, 7), ['X', 'Y', 'Z']].values) +# get the 6 largest clusters' max x, y, and z coordinates +coords = table.loc[range(1, 7), ['X', 'Y', 'Z']].values +# extract time series from each coordinate +masker = input_data.NiftiSpheresMasker(coords) real_timeseries = masker.fit_transform(fmri_img) -predicted_timeseries = masker.fit_transform(fmri_glm.predicted) +predicted_timeseries = masker.fit_transform(fmri_glm.predicted[0]) ######################################################################### @@ -81,37 +84,44 @@ # --------------------------------------------------------------------- import matplotlib.pyplot as plt -for i in range(1, 7): - plt.subplot(2, 3, i) - plt.title('Cluster peak {}\n'.format(table.loc[i, ['X', 'Y', 'Z']].tolist())) - plt.plot(real_timeseries[:, i-1], c='k', lw=2) - plt.plot(predicted_timeseries[:, i-1], c='r', ls='--', lw=2) - plt.xlabel('Time') - plt.ylabel('Signal intensity') +# colors for each of the clusters +colors = ['blue', 'navy', 'purple', 'magenta', 'olive', 'teal'] +# plot the time series and corresponding locations +fig1, axs1 = plt.subplots(2, 6) +for i in range(0, 6): + # plotting time series + axs1[0, i].set_title('Cluster peak {}\n'.format(coords[i])) + axs1[0, i].plot(real_timeseries[:, i], c=colors[i], lw=2) + axs1[0, i].plot(predicted_timeseries[:, i], c='r', ls='--', lw=2) + axs1[0, i].set_xlabel('Time') + axs1[0, i].set_ylabel('Signal intensity', labelpad=0) + # plotting image below the time series + roi_img = plotting.plot_stat_map( + z_map, cut_coords=[coords[i][2]], threshold=3.1, figure=fig1, + axes=axs1[1, i], display_mode='z', colorbar=False, bg_img=mean_img) + roi_img.add_markers([coords[i]], colors[i], 300) +fig1.set_size_inches(24, 14) -plt.gcf().set_size_inches(12, 7) -plt.tight_layout() ######################################################################### # Get residuals # ------------- -resid = masker.fit_transform(fmri_glm.residuals) +resid = masker.fit_transform(fmri_glm.residuals[0]) ######################################################################### # Plot distribution of residuals # ------------------------------ # Note that residuals are not really distributed normally. +fig2, axs2 = plt.subplots(2, 3) +axs2 = axs2.flatten() +for i in range(0, 6): + axs2[i].set_title('Cluster peak {}\n'.format(coords[i])) + axs2[i].hist(resid[:, i], color=colors[i]) + print('Mean residuals: {}'.format(resid[:, i].mean())) - -for i in range(1, 7): - plt.subplot(2, 3, i) - plt.title('Cluster peak {}\n'.format(table.loc[i, ['X', 'Y', 'Z']].tolist())) - plt.hist(resid[:, i-1]) - print('Mean residuals: {}'.format(resid[:, i-1].mean())) - -plt.gcf().set_size_inches(12, 7) -plt.tight_layout() +fig2.set_size_inches(12, 7) +fig2.tight_layout() ######################################################################### @@ -129,7 +139,7 @@ # or because the voxel/region fits the noise factors (such as drift or motion) # that could be present in the GLM. To isolate the influence of the experiment, # we can use an F-test as shown in the next section. -plotting.plot_stat_map(fmri_glm.r_square, +plotting.plot_stat_map(fmri_glm.r_square[0], bg_img=mean_img, threshold=.1, display_mode='z', cut_coords=7) @@ -143,12 +153,17 @@ import numpy as np design_matrix = fmri_glm.design_matrices_[0] + +# contrast with a one for "active" and zero everywhere else active = np.array([1 if c == 'active' else 0 for c in design_matrix.columns]) + +# contrast with a one for "rest" and zero everywhere else rest = np.array([1 if c == 'rest' else 0 for c in design_matrix.columns]) -effects_of_interest = np.vstack((active, rest)) -z_map = fmri_glm.compute_contrast(effects_of_interest, +effects_of_interest = np.vstack((active, rest)) +# f-test for rest and activity +z_map_ftest = fmri_glm.compute_contrast(effects_of_interest, output_type='z_score') -plotting.plot_stat_map(z_map, - bg_img=mean_img, threshold=2.33, display_mode='z', cut_coords=7) +plotting.plot_stat_map(z_map_ftest, + bg_img=mean_img, threshold=3.1, display_mode='z', cut_coords=7) diff --git a/nistats/tests/test_first_level_model.py b/nistats/tests/test_first_level_model.py index 8cb420e2..09165151 100644 --- a/nistats/tests/test_first_level_model.py +++ b/nistats/tests/test_first_level_model.py @@ -606,7 +606,7 @@ def test_first_level_model_residuals(): noise_model='ols') model.fit(fmri_data, design_matrices=design_matrices) - resid = model.residuals + resid = model.residuals[0] mean_resids = model.masker_.transform(resid).mean(0) assert_array_almost_equal(mean_resids, 0) @@ -624,9 +624,9 @@ def test_first_level_model_predictions_r_square(): noise_model='ols') model.fit(fmri_data, design_matrices=design_matrices) - pred = model.predicted + pred = model.predicted[0] data = fmri_data[0] - r_square_3d = model.r_square + r_square_3d = model.r_square[0] y_predicted = model.masker_.transform(pred) y_measured = model.masker_.transform(data) From afef0a9857f6d5cc754aa2880516af406216536b Mon Sep 17 00:00:00 2001 From: James Kent Date: Wed, 8 Jan 2020 20:15:54 -0600 Subject: [PATCH 23/32] make _get_voxelwise_model_attribute private and improve documentation --- nistats/first_level_model.py | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/nistats/first_level_model.py b/nistats/first_level_model.py index be25f5c6..f9265f72 100644 --- a/nistats/first_level_model.py +++ b/nistats/first_level_model.py @@ -37,6 +37,7 @@ from .regression import (ARModel, OLSModel, SimpleRegressionResults, + RegressionResults ) from .utils import (_basestring, _check_run_tables, @@ -585,7 +586,7 @@ def compute_contrast(self, contrast_def, stat_type=None, return outputs if output_type == 'all' else output - def get_voxelwise_model_attribute_(self, attribute, result_as_time_series=True): + def _get_voxelwise_model_attribute(self, attribute, result_as_time_series=True): """Transform RegressionResults instances within a dictionary (whose keys represent the autoregressive coefficient under the 'ar1' noise model or only 0.0 under 'ols' noise_model and values are the @@ -594,7 +595,9 @@ def get_voxelwise_model_attribute_(self, attribute, result_as_time_series=True): Parameters ---------- attribute : str - an attribute of a RegressionResults instance + an attribute of a RegressionResults instance. + possible values include: resid, norm_resid, predicted, + SSE, r_square, MSE. result_as_time_series : bool, optional whether the RegressionResult attribute has a value per timepoint of the input nifti image. @@ -604,6 +607,13 @@ def get_voxelwise_model_attribute_(self, attribute, result_as_time_series=True): output : list a list of Nifti1Image(s) """ + # check if valid attribute is being accessed. + all_attributes = dict(vars(RegressionResults)).keys() + possible_attributes = [prop for prop in all_attributes if '__' not in prop] + if attribute not in possible_attributes: + msg = "attribute must be one of: {attr}".format(attr=possible_attributes) + raise ValueError(msg) + if self.minimize_memory: raise ValueError('To access voxelwise attributes like R-squared, residuals, ' 'and predictions, the `FirstLevelModel`-object needs to store ' @@ -640,7 +650,7 @@ def residuals(self): output : list a list of Nifti1Image(s) """ - return self.get_voxelwise_model_attribute_('resid') + return self._get_voxelwise_model_attribute('resid') @setattr_on_read def predicted(self): @@ -652,7 +662,7 @@ def predicted(self): output : list a list of Nifti1Image(s) """ - return self.get_voxelwise_model_attribute_('predicted') + return self._get_voxelwise_model_attribute('predicted') @setattr_on_read def r_square(self): @@ -664,7 +674,7 @@ def r_square(self): output : list a list of Nifti1Image(s) """ - return self.get_voxelwise_model_attribute_('r_square', result_as_time_series=False) + return self._get_voxelwise_model_attribute('r_square', result_as_time_series=False) @replace_parameters({'mask': 'mask_img'}, end_version='next') From 1cbd072cbacc37bcaa0139cadb96d4a10a43a273 Mon Sep 17 00:00:00 2001 From: James Kent Date: Wed, 8 Jan 2020 21:27:07 -0600 Subject: [PATCH 24/32] fix formatting of function call --- .../02_first_level_models/plot_predictions_residuals.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/examples/02_first_level_models/plot_predictions_residuals.py b/examples/02_first_level_models/plot_predictions_residuals.py index 3410d3a7..4645d705 100644 --- a/examples/02_first_level_models/plot_predictions_residuals.py +++ b/examples/02_first_level_models/plot_predictions_residuals.py @@ -162,8 +162,10 @@ effects_of_interest = np.vstack((active, rest)) # f-test for rest and activity -z_map_ftest = fmri_glm.compute_contrast(effects_of_interest, - output_type='z_score') +z_map_ftest = fmri_glm.compute_contrast( + effects_of_interest, + stat_type='F', + output_type='z_score') plotting.plot_stat_map(z_map_ftest, bg_img=mean_img, threshold=3.1, display_mode='z', cut_coords=7) From 1a8e0c599673cdd687b94488a6db9d585708020f Mon Sep 17 00:00:00 2001 From: James Kent Date: Wed, 8 Jan 2020 21:28:21 -0600 Subject: [PATCH 25/32] add empty line back in --- nistats/regression.py | 1 + 1 file changed, 1 insertion(+) diff --git a/nistats/regression.py b/nistats/regression.py index ed3e57d3..87456a3f 100644 --- a/nistats/regression.py +++ b/nistats/regression.py @@ -61,6 +61,7 @@ class OLSModel(object): df_resid : scalar Degrees of freedom of the residuals. Number of observations less the rank of the design. + df_model : scalar Degrees of freedome of the model. The rank of the design. """ From e801f8e955ed266c2c498fea06ee55334585781a Mon Sep 17 00:00:00 2001 From: James Kent Date: Thu, 9 Jan 2020 10:59:16 -0600 Subject: [PATCH 26/32] revert regression.py to master --- nistats/regression.py | 36 ++++++++++++++++++++++++++---------- 1 file changed, 26 insertions(+), 10 deletions(-) diff --git a/nistats/regression.py b/nistats/regression.py index 87456a3f..9607d1b6 100644 --- a/nistats/regression.py +++ b/nistats/regression.py @@ -333,15 +333,13 @@ def MSE(self): return self.SSE / self.df_resid -class SimpleRegressionResults(RegressionResults): +class SimpleRegressionResults(LikelihoodModelResults): """This class contains only information of the model fit necessary for contast computation. - Its intended to save memory when details of the model are unnecessary. """ def __init__(self, results): """See LikelihoodModelResults constructor. - The only difference is that the whitened Y and residual values are stored for a regression model. """ @@ -355,19 +353,37 @@ 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('minimize_memory should be set to False to make residuals or predictions.') + raise ValueError('can not use this method for simple results') - def resid(self): + def resid(self, Y): """ Residuals from the fit. """ - raise ValueError('minimize_memory should be set to False to make residuals or predictions.') + return Y - self.predicted - def norm_resid(self): - raise ValueError('minimize_memory should be set to False to make residuals or predictions.') + def norm_resid(self, Y): + """ + Residuals, normalized to have unit length. + Notes + ----- + Is this supposed to return "stanardized residuals," + residuals standardized + to have mean zero and approximately unit variance? + d_i = e_i / sqrt(MS_E) + Where MS_E = SSE / (n - k) + See: Montgomery and Peck 3.2.1 p. 68 + Davidson and MacKinnon 15.2 p 662 + """ + return self.resid(Y) * positive_reciprocal(np.sqrt(self.dispersion)) + + def predicted(self): + """ 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) \ No newline at end of file From bf670662cffab51c15899d097b2d4392e76eadf4 Mon Sep 17 00:00:00 2001 From: James Kent Date: Thu, 9 Jan 2020 11:00:30 -0600 Subject: [PATCH 27/32] make result_as_time_series mandatory --- nistats/first_level_model.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/nistats/first_level_model.py b/nistats/first_level_model.py index f9265f72..ef4e770c 100644 --- a/nistats/first_level_model.py +++ b/nistats/first_level_model.py @@ -586,7 +586,7 @@ def compute_contrast(self, contrast_def, stat_type=None, return outputs if output_type == 'all' else output - def _get_voxelwise_model_attribute(self, attribute, result_as_time_series=True): + def _get_voxelwise_model_attribute(self, attribute, result_as_time_series): """Transform RegressionResults instances within a dictionary (whose keys represent the autoregressive coefficient under the 'ar1' noise model or only 0.0 under 'ols' noise_model and values are the @@ -598,7 +598,7 @@ def _get_voxelwise_model_attribute(self, attribute, result_as_time_series=True): an attribute of a RegressionResults instance. possible values include: resid, norm_resid, predicted, SSE, r_square, MSE. - result_as_time_series : bool, optional + result_as_time_series : bool whether the RegressionResult attribute has a value per timepoint of the input nifti image. @@ -650,7 +650,7 @@ def residuals(self): output : list a list of Nifti1Image(s) """ - return self._get_voxelwise_model_attribute('resid') + return self._get_voxelwise_model_attribute('resid', result_as_time_series=True) @setattr_on_read def predicted(self): @@ -662,7 +662,7 @@ def predicted(self): output : list a list of Nifti1Image(s) """ - return self._get_voxelwise_model_attribute('predicted') + return self._get_voxelwise_model_attribute('predicted', result_as_time_series=True) @setattr_on_read def r_square(self): From b55c13b7dff4900cf7f9c56cea83ba5bbc875b18 Mon Sep 17 00:00:00 2001 From: James Kent Date: Thu, 9 Jan 2020 12:18:53 -0600 Subject: [PATCH 28/32] add newlines to docs --- nistats/regression.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/nistats/regression.py b/nistats/regression.py index 9607d1b6..036a34f0 100644 --- a/nistats/regression.py +++ b/nistats/regression.py @@ -336,10 +336,12 @@ def MSE(self): class SimpleRegressionResults(LikelihoodModelResults): """This class contains only information of the model fit necessary for contast computation. + Its intended to save memory when details of the model are unnecessary. """ def __init__(self, results): """See LikelihoodModelResults constructor. + The only difference is that the whitened Y and residual values are stored for a regression model. """ @@ -373,8 +375,11 @@ def norm_resid(self, Y): Is this supposed to return "stanardized residuals," residuals standardized to have mean zero and approximately unit variance? + d_i = e_i / sqrt(MS_E) + Where MS_E = SSE / (n - k) + See: Montgomery and Peck 3.2.1 p. 68 Davidson and MacKinnon 15.2 p 662 """ From 999caf99b25329fb6499c3cfc235a621415dd382 Mon Sep 17 00:00:00 2001 From: James Kent Date: Thu, 9 Jan 2020 12:19:17 -0600 Subject: [PATCH 29/32] add newline to end of file --- nistats/regression.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nistats/regression.py b/nistats/regression.py index 036a34f0..36a1f928 100644 --- a/nistats/regression.py +++ b/nistats/regression.py @@ -391,4 +391,4 @@ def predicted(self): beta = self.theta # the LikelihoodModelResults has parameters named 'theta' X = self.model.design - return np.dot(X, beta) \ No newline at end of file + return np.dot(X, beta) From 4b5d485296fab2b01780e884f33c6129db320f48 Mon Sep 17 00:00:00 2001 From: James Kent Date: Thu, 9 Jan 2020 12:21:49 -0600 Subject: [PATCH 30/32] fix missing newline --- nistats/regression.py | 1 + 1 file changed, 1 insertion(+) diff --git a/nistats/regression.py b/nistats/regression.py index 36a1f928..37e64657 100644 --- a/nistats/regression.py +++ b/nistats/regression.py @@ -370,6 +370,7 @@ def resid(self, Y): def norm_resid(self, Y): """ Residuals, normalized to have unit length. + Notes ----- Is this supposed to return "stanardized residuals," From 78995d02577a75448053f19f2f70b57e9ce8c206 Mon Sep 17 00:00:00 2001 From: James Kent Date: Fri, 10 Jan 2020 20:40:54 -0600 Subject: [PATCH 31/32] add James Kent to .mailmap --- .mailmap | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.mailmap b/.mailmap index 5c229607..85ac94db 100644 --- a/.mailmap +++ b/.mailmap @@ -21,6 +21,9 @@ Fabian Pedregosa Franz Liem Gael Varoquaux Greg Kiar +James D. Kent +James D. Kent +James D. Kent Fred Mertz Jan Margeta Jaques Grobler Jason Gors From 44b0e85ffa0b9abc49c9fa76edc9d96281c5cbf7 Mon Sep 17 00:00:00 2001 From: James Kent Date: Fri, 10 Jan 2020 20:42:17 -0600 Subject: [PATCH 32/32] add entry for the new attributes to FirstLevelModel --- doc/whats_new.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/doc/whats_new.rst b/doc/whats_new.rst index ebd28b3f..bc17e417 100644 --- a/doc/whats_new.rst +++ b/doc/whats_new.rst @@ -9,6 +9,9 @@ New --- +* :func:`nistats.first_level_model.FirstLevelModel` now has the attributes: ``residuals``, ``predicted``, and ``r_square`` + which returns a Niimg-like object in the same shape as the input Niimg-like object. + Additionally, there is an example showcasing the use of the attributes. * Use :func:`nistats.reporting.make_glm_report` to easily generate HTML reports from fitted first and second level models and contrasts. * New dataset fetcher, :func:`nistats.datasets.fetch_language_localizer_demo_dataset` , BIDS 1.2 compatible. * New example showcasing the use of a GLM to get beta maps for decoding experiments (aka beta-regression).