IBMA aggressive masks are inherited across fits, causing opaque errors

[tsalo]

Summary

When you fit the same IBMA estimator multiple times, each time the aggressive mask is updated based on the new data and the existing aggressive mask using a logical or, meaning that good voxels in an older dataset will be retained when fitting to the new dataset. This causes problems when some of those older voxels are not good for the new dataset.

This problem wasn't already apparent because all of our Jackknife tests on IBMAs had good data in all of the voxels, I believe. Or at least the aggressive mask didn't change from one subset of studies to another.

I'm not sure what the best solution is. Here are some ideas:

• Automatically replace the aggressive mask when fitting an estimator to a new dataset. I think this would just be a step in _preprocess_inputs to delete that key from the inputs_ dictionary if it exists at the beginning of the preprocessing.
• Add an argument to amend or replace the mask to the IBMA estimator initialization, so users can decide.

Additional details

• NiMARE version: ~0.0.12rc1 (current main branch)

This issue stems from #649 (comment).

What were you trying to do?

Run a Jackknife diagnostic on a DerSimonianLaird Estimator.

What did you expect to happen?

The aggressive mask would be recalculated from scratch on each subset Dataset within the Jackknife process.

What actually happened?

The same Estimator was used for each subset Dataset, and the aggressive mask was expanded based on the good voxels in each subset, so when study 1 was removed, that added X voxels, but then when study 2 was removed and study 1 was included, those voxels, even though they're "bad" for study 1, were included, leading to an SVD convergence failure in the second Jackknife meta-analysis.

Reproducing the bug

import os
import warnings

from nimare.dataset import Dataset
from nimare.extract import download_nidm_pain
from nimare.meta.ibma import DerSimonianLaird
from nimare.transforms import ImageTransformer
from nimare.utils import get_resource_path

# Ignore resampling-related warnings from nilearn for the sake of a clean terminal
warnings.filterwarnings(action="ignore", module="nilearn", category=RuntimeWarning)

dset_dir = download_nidm_pain()

dset_file = os.path.join(get_resource_path(), "nidm_pain_dset.json")
dset = Dataset(dset_file)
dset.update_path(dset_dir)

# Calculate missing statistical images from the available stats.
xformer = ImageTransformer(target=["varcope"])
dset = xformer.transform(dset)

# Initialize Estimator and fit to full Dataset
meta_ibma = DerSimonianLaird(resample=True)
results = meta_ibma.fit(dset)

# Fit the same Estimator to the same Dataset, but without study 1
temp_dset2 = dset.slice(dset.ids[1:])
results2 = meta_ibma.fit(temp_dset2)

# Fit to the Dataset minus study 2
temp_ids = list(dset.ids)
temp_ids = temp_ids[:1] + temp_ids[2:]
temp_dset3 = dset.slice(temp_ids)
results3 = meta_ibma.fit(temp_dset3)  # fails here

Here's the traceback:

# After the first fit
WARNING:nimare.base:Masking out 39113 additional voxels. The updated masker is available in the Estimator.masker attribute.
# Fit on dataset without study 1
WARNING:nimare.base:Masking out 37645 additional voxels. The updated masker is available in the Estimator.masker attribute.
# Fit on dataset without study 2
WARNING:nimare.base:Masking out 36039 additional voxels. The updated masker is available in the Estimator.masker attribute.
/Users/taylor/Documents/tsalo/PyMARE/pymare/stats.py:23: RuntimeWarning: divide by zero encountered in true_divide
  w = 1.0 / (v + tau2)
/Users/taylor/Documents/tsalo/PyMARE/pymare/estimators/estimators.py:202: RuntimeWarning: divide by zero encountered in true_divide
  w = 1.0 / v
Traceback (most recent call last):
  File "test_jackknife.py", line 35, in <module>
    results3 = meta_ibma.fit(temp_dset3)
  File "/Users/taylor/Documents/tsalo/NiMARE/nimare/base.py", line 314, in fit
    maps = self._fit(dataset)
  File "/Users/taylor/Documents/tsalo/NiMARE/nimare/meta/ibma.py", line 292, in _fit
    est.fit_dataset(pymare_dset)
  File "/Users/taylor/Documents/tsalo/PyMARE/pymare/estimators/estimators.py", line 88, in fit_dataset
    self.fit(*args, **all_kwargs)
  File "/Users/taylor/Documents/tsalo/PyMARE/pymare/estimators/estimators.py", line 215, in fit
    beta_dl, inv_cov = weighted_least_squares(y, v, X, tau2=tau_dl, return_cov=True)
  File "/Users/taylor/Documents/tsalo/PyMARE/pymare/stats.py", line 31, in weighted_least_squares
    precision = np.linalg.pinv(cov).T
  File "<__array_function__ internals>", line 5, in pinv
  File "/opt/miniconda3/lib/python3.8/site-packages/numpy/linalg/linalg.py", line 2002, in pinv
    u, s, vt = svd(a, full_matrices=False, hermitian=hermitian)
  File "<__array_function__ internals>", line 5, in svd
  File "/opt/miniconda3/lib/python3.8/site-packages/numpy/linalg/linalg.py", line 1660, in svd
    u, s, vh = gufunc(a, signature=signature, extobj=extobj)
  File "/opt/miniconda3/lib/python3.8/site-packages/numpy/linalg/linalg.py", line 97, in _raise_linalgerror_svd_nonconvergence
    raise LinAlgError("SVD did not converge")
numpy.linalg.LinAlgError: SVD did not converge