Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[ENH] Do not apply IBMA methods to voxels with zeros or NaNs #544

Merged
merged 4 commits into from
Jul 26, 2021
Merged

[ENH] Do not apply IBMA methods to voxels with zeros or NaNs #544

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
29b258e
Start implementing a masker-update step.
tsalo Jul 22, 2021
47bb838
Replace masker update with boolean masking.
tsalo Jul 22, 2021
8a5ec99
Squeeze.
tsalo Jul 22, 2021
3c7a9ab
Remove tau2 from WLS results.
tsalo Jul 23, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
40 changes: 35 additions & 5 deletions nimare/base.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -348,6 +348,12 @@ def _preprocess_input(self, dataset):
if isinstance(mask_img, str):
mask_img = nb.load(mask_img)

# Ensure that protected values are not included among _required_inputs
assert "aggressive_mask" not in self._required_inputs.keys(), "This is a protected name."

# A dictionary to collect masked image data, to be further reduced by the aggressive mask.
temp_image_inputs = {}

for name, (type_, _) in self._required_inputs.items():
if type_ == "image":
# If no resampling is requested, check if resampling is required
Expand All @@ -371,12 +377,22 @@ def _preprocess_input(self, dataset):

# An intermediate step to mask out bad voxels.
# Can be dropped once PyMARE is able to handle masked arrays or missing data.
bad_voxel_idx = np.where(temp_arr == 0)[1]
bad_voxel_idx = np.unique(bad_voxel_idx)
LGR.debug(f"Masking out {len(bad_voxel_idx)} 'bad' voxels")
temp_arr[:, bad_voxel_idx] = 0
nonzero_voxels_bool = np.all(temp_arr != 0, axis=0)
nonnan_voxels_bool = np.all(~np.isnan(temp_arr), axis=0)
good_voxels_bool = np.logical_and(nonzero_voxels_bool, nonnan_voxels_bool)

data = masker.transform(img4d)

temp_image_inputs[name] = data
if "aggressive_mask" not in self.inputs_.keys():
self.inputs_["aggressive_mask"] = good_voxels_bool
else:
# Remove any voxels that are bad in any image-based inputs
self.inputs_["aggressive_mask"] = np.logical_or(
self.inputs_["aggressive_mask"],
good_voxels_bool,
)

self.inputs_[name] = temp_arr
elif type_ == "coordinates":
# Try to load existing MA maps
if hasattr(self, "kernel_transformer"):
Expand All @@ -389,6 +405,20 @@ def _preprocess_input(self, dataset):
if all(f is not None for f in files):
self.inputs_["ma_maps"] = files

# Further reduce image-based inputs to remove "bad" voxels
# (voxels with zeros or NaNs in any studies)
if "aggressive_mask" in self.inputs_.keys():
n_bad_voxels = (
self.inputs_["aggressive_mask"].size - self.inputs_["aggressive_mask"].sum()
)
if n_bad_voxels:
LGR.warning(
f"Masking out {n_bad_voxels} additional voxels. "
"The updated masker is available in the Estimator.masker attribute."
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This was a mistake. The masker is not updated directly.

)
for name, raw_masked_data in temp_image_inputs.items():
self.inputs_[name] = raw_masked_data[:, self.inputs_["aggressive_mask"]]


class Transformer(NiMAREBase):
"""Transformers take in Datasets and return Datasets.
Expand Down
91 changes: 65 additions & 26 deletions nimare/meta/ibma.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@

from ..base import MetaEstimator
from ..transforms import p_to_z, t_to_z
from ..utils import boolean_unmask

LGR = logging.getLogger(__name__)

Expand Down Expand Up @@ -65,8 +66,8 @@ def _fit(self, dataset):
est.fit_dataset(pymare_dset)
est_summary = est.summary()
results = {
"z": est_summary.z,
"p": est_summary.p,
"z": boolean_unmask(est_summary.z.squeeze(), self.inputs_["aggressive_mask"]),
"p": boolean_unmask(est_summary.p.squeeze(), self.inputs_["aggressive_mask"]),
}
return results

Expand Down Expand Up @@ -144,8 +145,8 @@ def _fit(self, dataset):
est_summary = est.summary()

results = {
"z": est_summary.z,
"p": est_summary.p,
"z": boolean_unmask(est_summary.z.squeeze(), self.inputs_["aggressive_mask"]),
"p": boolean_unmask(est_summary.p.squeeze(), self.inputs_["aggressive_mask"]),
}
return results

Expand Down Expand Up @@ -214,11 +215,17 @@ def _fit(self, dataset):
est = pymare.estimators.WeightedLeastSquares(tau2=self.tau2)
est.fit_dataset(pymare_dset)
est_summary = est.summary()
# tau2 is an float, not a map, so it can't go in the results dictionary
results = {
"tau2": est_summary.tau2,
"z": est_summary.get_fe_stats()["z"].squeeze(),
"p": est_summary.get_fe_stats()["p"].squeeze(),
"est": est_summary.get_fe_stats()["est"].squeeze(),
"z": boolean_unmask(
est_summary.get_fe_stats()["z"].squeeze(), self.inputs_["aggressive_mask"]
),
"p": boolean_unmask(
est_summary.get_fe_stats()["p"].squeeze(), self.inputs_["aggressive_mask"]
),
"est": boolean_unmask(
est_summary.get_fe_stats()["est"].squeeze(), self.inputs_["aggressive_mask"]
),
}
return results

Expand Down Expand Up @@ -282,10 +289,16 @@ def _fit(self, dataset):
est.fit_dataset(pymare_dset)
est_summary = est.summary()
results = {
"tau2": est_summary.tau2,
"z": est_summary.get_fe_stats()["z"].squeeze(),
"p": est_summary.get_fe_stats()["p"].squeeze(),
"est": est_summary.get_fe_stats()["est"].squeeze(),
"tau2": boolean_unmask(est_summary.tau2.squeeze(), self.inputs_["aggressive_mask"]),
"z": boolean_unmask(
est_summary.get_fe_stats()["z"].squeeze(), self.inputs_["aggressive_mask"]
),
"p": boolean_unmask(
est_summary.get_fe_stats()["p"].squeeze(), self.inputs_["aggressive_mask"]
),
"est": boolean_unmask(
est_summary.get_fe_stats()["est"].squeeze(), self.inputs_["aggressive_mask"]
),
}
return results

Expand Down Expand Up @@ -345,10 +358,16 @@ def _fit(self, dataset):
est.fit_dataset(pymare_dset)
est_summary = est.summary()
results = {
"tau2": est_summary.tau2,
"z": est_summary.get_fe_stats()["z"].squeeze(),
"p": est_summary.get_fe_stats()["p"].squeeze(),
"est": est_summary.get_fe_stats()["est"].squeeze(),
"tau2": boolean_unmask(est_summary.tau2.squeeze(), self.inputs_["aggressive_mask"]),
"z": boolean_unmask(
est_summary.get_fe_stats()["z"].squeeze(), self.inputs_["aggressive_mask"]
),
"p": boolean_unmask(
est_summary.get_fe_stats()["p"].squeeze(), self.inputs_["aggressive_mask"]
),
"est": boolean_unmask(
est_summary.get_fe_stats()["est"].squeeze(), self.inputs_["aggressive_mask"]
),
}
return results

Expand Down Expand Up @@ -414,10 +433,16 @@ def _fit(self, dataset):
est.fit_dataset(pymare_dset)
est_summary = est.summary()
results = {
"tau2": est_summary.tau2,
"z": est_summary.get_fe_stats()["z"].squeeze(),
"p": est_summary.get_fe_stats()["p"].squeeze(),
"est": est_summary.get_fe_stats()["est"].squeeze(),
"tau2": boolean_unmask(est_summary.tau2.squeeze(), self.inputs_["aggressive_mask"]),
"z": boolean_unmask(
est_summary.get_fe_stats()["z"].squeeze(), self.inputs_["aggressive_mask"]
),
"p": boolean_unmask(
est_summary.get_fe_stats()["p"].squeeze(), self.inputs_["aggressive_mask"]
),
"est": boolean_unmask(
est_summary.get_fe_stats()["est"].squeeze(), self.inputs_["aggressive_mask"]
),
}
return results

Expand Down Expand Up @@ -498,10 +523,16 @@ def _fit(self, dataset):
est.fit_dataset(pymare_dset)
est_summary = est.summary()
results = {
"tau2": est_summary.tau2,
"z": est_summary.get_fe_stats()["z"].squeeze(),
"p": est_summary.get_fe_stats()["p"].squeeze(),
"est": est_summary.get_fe_stats()["est"].squeeze(),
"tau2": boolean_unmask(est_summary.tau2.squeeze(), self.inputs_["aggressive_mask"]),
"z": boolean_unmask(
est_summary.get_fe_stats()["z"].squeeze(), self.inputs_["aggressive_mask"]
),
"p": boolean_unmask(
est_summary.get_fe_stats()["p"].squeeze(), self.inputs_["aggressive_mask"]
),
"est": boolean_unmask(
est_summary.get_fe_stats()["est"].squeeze(), self.inputs_["aggressive_mask"]
),
}
return results

Expand Down Expand Up @@ -572,7 +603,10 @@ def _fit(self, dataset):
# Convert t to z, preserving signs
dof = self.parameters_["tested_vars"].shape[0] - self.parameters_["tested_vars"].shape[1]
z_map = t_to_z(t_map, dof)
images = {"t": t_map.squeeze(), "z": z_map.squeeze()}
images = {
"t": boolean_unmask(t_map.squeeze(), self.inputs_["aggressive_mask"]),
"z": boolean_unmask(z_map.squeeze(), self.inputs_["aggressive_mask"]),
}
return images

def correct_fwe_montecarlo(self, result, n_iters=10000, n_cores=-1):
Expand Down Expand Up @@ -639,5 +673,10 @@ def correct_fwe_montecarlo(self, result, n_iters=10000, n_cores=-1):
sign = np.sign(t_map)
sign[sign == 0] = 1
z_map = p_to_z(p_map, tail="two") * sign
images = {"logp_level-voxel": log_p_map.squeeze(), "z_level-voxel": z_map.squeeze()}
images = {
"logp_level-voxel": boolean_unmask(
log_p_map.squeeze(), self.inputs_["aggressive_mask"]
),
"z_level-voxel": boolean_unmask(z_map.squeeze(), self.inputs_["aggressive_mask"]),
}
return images
4 changes: 4 additions & 0 deletions nimare/tests/test_meta_ibma.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
import os.path as op
from contextlib import ExitStack as does_not_raise

import numpy as np
import pytest
from nilearn.input_data import NiftiLabelsMasker

Expand Down Expand Up @@ -156,7 +157,10 @@ def test_ibma_with_custom_masker(testdata_ibma, caplog, estimator, expectation,
# Only fit the estimator if it doesn't raise a ValueError
if expectation != "error":
assert isinstance(meta.results, nimare.results.MetaResult)
# There are five "labels", but one of them has no good data,
# so the outputs should be 4 long.
assert meta.results.maps["z"].shape == (5,)
assert np.isnan(meta.results.maps["z"][0])
assert meta.results.get_map("z").shape == (10, 10, 10)


Expand Down
35 changes: 35 additions & 0 deletions nimare/utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -897,3 +897,38 @@ def mni2tal(coords):
if use_dim == 1:
out_coords = out_coords.transpose()
return out_coords


def boolean_unmask(data_array, bool_array):
"""Unmask data based on a boolean array, with NaNs in empty voxels.

Parameters
----------
data_array : 1D or 2D :obj:`numpy.ndarray`
Masked data array.
bool_array : 1D :obj:`numpy.ndarray`
Boolean mask array. Must have the same number of ``True`` entries as elements in the
second dimension of ``data_array``.

Returns
-------
unmasked_data : 1D or 2D :obj:`numpy.ndarray`
Unmasked data array.
If 1D, first dimension is the same size as the first (and only) dimension of
``boolean_array``.
If 2D, first dimension is the same size as the first dimension of ``data_array``, while
second dimension is the same size as the first (and only) dimension of ``boolean_array``.
All elements corresponding to ``False`` values in ``boolean_array`` will have NaNs.
"""
assert data_array.ndim in (1, 2)
assert bool_array.ndim == 1
assert bool_array.sum() == data_array.shape[-1]

unmasked_data = np.full(
shape=bool_array.shape + data_array.T.shape[1:],
fill_value=np.nan,
dtype=data_array.dtype,
)
unmasked_data[bool_array] = data_array
unmasked_data = unmasked_data.T
return unmasked_data