diff --git a/fitlins/interfaces/nistats.py b/fitlins/interfaces/nistats.py index 773d7f5d..93ec590c 100644 --- a/fitlins/interfaces/nistats.py +++ b/fitlins/interfaces/nistats.py @@ -3,6 +3,8 @@ import pandas as pd from functools import partial + + from nipype.interfaces.base import LibraryBaseInterface, SimpleInterface, isdefined from .abstract import ( @@ -322,7 +324,7 @@ def _run_interface(self, runtime): compute_fixed_effects, _compute_fixed_effects_params, ) - + from ..interfaces.pymare_extension import pymare_model spec = self.inputs.spec smoothing_fwhm = self.inputs.smoothing_fwhm smoothing_type = self.inputs.smoothing_type @@ -332,7 +334,6 @@ def _run_interface(self, runtime): raise NotImplementedError( "Only the iso smoothing type is available for the nistats estimator." ) - effect_maps = [] variance_maps = [] stat_maps = [] @@ -342,12 +343,11 @@ def _run_interface(self, runtime): spec_metadata = spec['metadata'].to_dict('records') out_ents = spec['entities'].copy() # Same for all out_ents.setdefault("contrast", spec['contrasts'][0]['name']) - + # Only keep files which match all entities for contrast stat_metadata = _flatten(self.inputs.stat_metadata) input_effects = _flatten(self.inputs.effect_maps) input_variances = _flatten(self.inputs.variance_maps) - filtered_effects = [] filtered_variances = [] names = [] @@ -361,7 +361,6 @@ def _run_interface(self, runtime): names.append(m['contrast']) contrasts = prepare_contrasts(spec['contrasts'], spec['X'].columns) - is_cifti = filtered_effects[0].endswith('dscalar.nii') if is_cifti: fname_fmt = os.path.join(runtime.cwd, '{}_{}.dscalar.nii').format @@ -369,9 +368,8 @@ def _run_interface(self, runtime): fname_fmt = os.path.join(runtime.cwd, '{}_{}.nii.gz').format model_type = spec["model"].get("type", "") - - # Do not fit model for meta-analyses - if model_type != 'Meta': + + if model_type.lower() != "meta": if len(filtered_effects) < 2: raise RuntimeError( "At least two inputs are required for a 't' for 'F' " "second level contrast" @@ -386,6 +384,13 @@ def _run_interface(self, runtime): else: model = level2.SecondLevelModel(smoothing_fwhm=smoothing_fwhm) model.fit(filtered_effects, design_matrix=spec['X']) + else: + if is_cifti: + model = pymare_model(is_cifti=True) + model.fit(filtered_effects, filtered_variances, spec['X']) + else: + model = pymare_model(is_cifti=False) + model.fit(filtered_effects, filtered_variances, spec['X']) for name, weights, cont_ents, contrast_test in contrasts: contrast_metadata.append( @@ -397,40 +402,7 @@ def _run_interface(self, runtime): } ) - # Pass-through happens automatically as it can handle 1 input - if model_type == 'Meta': - # Index design identity matrix on non-zero contrasts weights - con_ix = weights[0].astype(bool) - # Index of all input files "involved" with that contrast - dm_ix = spec['X'].iloc[:, con_ix].any(axis=1) - - contrast_imgs = np.array(filtered_effects)[dm_ix] - variance_imgs = np.array(filtered_variances)[dm_ix] - if is_cifti: - ffx_cont, ffx_var, ffx_t = _compute_fixed_effects_params( - np.squeeze( - [nb.load(fname).get_fdata(dtype='f4') for fname in contrast_imgs] - ), - np.squeeze( - [nb.load(fname).get_fdata(dtype='f4') for fname in variance_imgs] - ), - precision_weighted=False, - ) - img = nb.load(filtered_effects[0]) - maps = { - 'effect_size': dscalar_from_cifti(img, ffx_cont, "effect_size"), - 'effect_variance': dscalar_from_cifti(img, ffx_var, "effect_variance"), - 'stat': dscalar_from_cifti(img, ffx_t, "stat"), - } - - else: - ffx_res = compute_fixed_effects(contrast_imgs, variance_imgs) - maps = { - 'effect_size': ffx_res[0], - 'effect_variance': ffx_res[1], - 'stat': ffx_res[2], - } - else: + if model_type.lower() != "meta": if is_cifti: contrast = compute_contrast( labels, estimates, weights, contrast_type=contrast_test @@ -452,6 +424,16 @@ def _run_interface(self, runtime): second_level_stat_type=contrast_test, output_type='all', ) + else: + if is_cifti: + dict_maps = model.compute_contrast(weights) + img = nb.load(filtered_effects[0]) + maps = { + map_type: dscalar_from_cifti(img, dict_maps[map_type], map_type) + for map_type in dict_maps.keys() + } + else: + maps = model.compute_contrast(weights) for map_type, map_list in ( ('effect_size', effect_maps), diff --git a/fitlins/interfaces/pymare_extension.py b/fitlins/interfaces/pymare_extension.py new file mode 100644 index 00000000..0a2677b2 --- /dev/null +++ b/fitlins/interfaces/pymare_extension.py @@ -0,0 +1,80 @@ +import nibabel as nb +import numpy as np +from nilearn.input_data import NiftiMasker +from pymare import estimators +from nilearn.image import mean_img + + +def remove_zero_var_voxels_(effect_data, variance_data): + nonzero_var_mask = np.sum(variance_data == 0, 0) == 0 + variance_data_masked = variance_data[:, nonzero_var_mask] + effect_data_masked = effect_data[:, nonzero_var_mask] + return effect_data_masked, variance_data_masked, nonzero_var_mask + + + +class pymare_model: + def __init__( + self, + smoothing_fwhm=None, + is_cifti=False + ): + self.smoothing_fwhm = smoothing_fwhm + self.is_cifti = is_cifti + + def fit( + self, + filtered_effects=None, + filtered_variances=None, + design_matrix=None + ): + self.design_matrix_ = design_matrix.to_numpy() + if self.is_cifti: + effect_data = np.squeeze( + [nb.load(effect).get_fdata(dtype='f4') for effect in filtered_effects] + ) + variance_data = np.squeeze( + [nb.load(variance).get_fdata(dtype='f4') for variance in filtered_variances] + ) + self.effect_data_, self.variance_data_, self.nonzero_var_mask_ = \ + remove_zero_var_voxels_(effect_data, variance_data) + + else: + self.masker_ = NiftiMasker( + smoothing_fwhm=self.smoothing_fwhm, + ) + sample_map = mean_img(filtered_effects) + self.masker_.fit(sample_map) + effect_data = self.masker_.transform(filtered_effects) + variance_data = self.masker_.transform(filtered_variances) + self.effect_data_, self.variance_data_, self.nonzero_var_mask_ = \ + remove_zero_var_voxels_(effect_data, variance_data) + + self.wls_ = estimators.WeightedLeastSquares() + self.wls_.fit(y=self.effect_data_, + v=self.variance_data_, + X=self.design_matrix_ + ) + + def compute_contrast( + self, + con_val=None, + ): + outputs = {} + tmp = np.einsum('ij, jkl->kl', con_val, self.wls_.params_['inv_cov']) + outputs['effect_variance'] = np.einsum('ij, jk->k', con_val, tmp) + outputs['effect_size'] = np.einsum('ij, jk->k', con_val, self.wls_.params_['fe_params']) + outputs['stat'] = outputs['effect_size'] / outputs['effect_variance']**.5 + + outputs_array = {} + for image_type, image in outputs.items(): + outputs_array[image_type] = np.zeros(self.nonzero_var_mask_.shape) + outputs_array[image_type][self.nonzero_var_mask_] = image + + if self.is_cifti is False: + outputs_nifti = {} + for image_type, image in outputs_array.items(): + outputs_nifti[image_type] = self.masker_.inverse_transform(image) + return outputs_nifti + else: + return outputs_array diff --git a/setup.cfg b/setup.cfg index 9c16e749..9ddf82cd 100644 --- a/setup.cfg +++ b/setup.cfg @@ -31,6 +31,7 @@ install_requires = tables>=3.2.1 pybids~=0.15.4 jinja2 + pymare>=0.0.3 [options.extras_require] duecredit = duecredit