From bc28e6925e3a5a066139ca246393baf6763d97b2 Mon Sep 17 00:00:00 2001 From: JulioAPeraza Date: Wed, 24 Apr 2024 16:39:13 -0400 Subject: [PATCH 1/4] Support precomputed correlation matrix for calculating variance inflation term in Stouffers --- pymare/estimators/combination.py | 28 +++++++++++++++++++------- pymare/tests/test_combination_tests.py | 25 +++++++++++++++++++++++ 2 files changed, 46 insertions(+), 7 deletions(-) diff --git a/pymare/estimators/combination.py b/pymare/estimators/combination.py index 6f9f41a..5fd4fc4 100644 --- a/pymare/estimators/combination.py +++ b/pymare/estimators/combination.py @@ -113,7 +113,7 @@ class StoufferCombinationTest(CombinationTest): # Maps Dataset attributes onto fit() args; see BaseEstimator for details. _dataset_attr_map = {"z": "y", "w": "n", "g": "v"} - def _inflation_term(self, z, w, g): + def _inflation_term(self, z, w, g, corr=None): """Calculate the variance inflation term for each group. This term is used to adjust the variance of the combined z-score when @@ -127,6 +127,8 @@ def _inflation_term(self, z, w, g): Array of weights. g : :obj:`numpy.ndarray` of shape (n, d) Array of group labels. + corr : :obj:`numpy.ndarray` of shape (n, n), optional + The correlation matrix of the z-values. If None, it will be calculated. Returns ------- @@ -157,26 +159,38 @@ def _inflation_term(self, z, w, g): continue # Calculate the within group correlation matrix and sum the non-diagonal elements - corr = np.corrcoef(group_z, rowvar=True) + if corr is None: + if z.shape[1] < 2: + raise ValueError("The number of features must be greater than 1.") + group_corr = np.corrcoef(group_z, rowvar=True) + else: + group_corr = corr[group_indices][:, group_indices] + upper_indices = np.triu_indices(n_samples, k=1) - non_diag_corr = corr[upper_indices] + non_diag_corr = group_corr[upper_indices] w_i, w_j = weights[upper_indices[0]], weights[upper_indices[1]] sigma += (2 * w_i * w_j * non_diag_corr).sum() return sigma - def fit(self, z, w=None, g=None): + def fit(self, z, w=None, g=None, corr=None): """Fit the estimator to z-values, optionally with weights and groups.""" - return super().fit(z, w=w, g=g) + return super().fit(z, w=w, g=g, corr=corr) - def p_value(self, z, w=None, g=None): + def p_value(self, z, w=None, g=None, corr=None): """Calculate p-values.""" if w is None: w = np.ones_like(z) + if g is None and corr is not None: + warnings.warn("Correlation matrix provided without groups. Ignoring.") + + if g is not None and corr is not None and g.shape[0] != corr.shape[0]: + raise ValueError("Group labels must have the same length as the correlation matrix.") + # Calculate the variance inflation term, sum of non-diagonal elements of sigma. - sigma = self._inflation_term(z, w, g) if g is not None else 0 + sigma = self._inflation_term(z, w, g, corr=corr) if g is not None else 0 # The sum of diagonal elements of sigma is given by (w**2).sum(0). variance = (w**2).sum(0) + sigma diff --git a/pymare/tests/test_combination_tests.py b/pymare/tests/test_combination_tests.py index 9f5f7e1..ac6e956 100644 --- a/pymare/tests/test_combination_tests.py +++ b/pymare/tests/test_combination_tests.py @@ -79,3 +79,28 @@ def test_stouffer_adjusted(): sigma_l1 = n_maps_l1 * (n_maps_l1 - 1) # Expected inflation term z_expected_l1 = n_maps_l1 * common_sample / np.sqrt(n_maps_l1 + sigma_l1) assert np.allclose(z_l1, z_expected_l1, atol=1e-5) + + # Test with correlation matrix and groups. + data_corr = data - data.mean(0) + corr = np.corrcoef(data_corr, rowvar=True) + results_corr = ( + StoufferCombinationTest("directed").fit(z=data, w=weights, g=groups, corr=corr).params_ + ) + z_corr = ss.norm.isf(results_corr["p"]) + + z_corr_expected = np.array([5.00088912, 3.70356943, 4.05465924, 5.4633001, 5.18927878]) + assert np.allclose(z_corr, z_corr_expected, atol=1e-5) + + # Test with no correlation matrix and groups, but only one feature. + with pytest.raises(ValueError): + StoufferCombinationTest("directed").fit(z=data[:, :1], w=weights[:, :1], g=groups) + + # Test with correlation matrix and groups of different shapes. + with pytest.raises(ValueError): + StoufferCombinationTest("directed").fit(z=data, w=weights, g=groups, corr=corr[:-2, :-2]) + + # Test with correlation matrix and no groups. + results1 = StoufferCombinationTest("directed").fit(z=_z1, corr=corr).params_ + z1 = ss.norm.isf(results1["p"]) + + assert np.allclose(z1, [4.69574], atol=1e-5) From 71463d733eccb73c83dffe882c7511ca639bbdb9 Mon Sep 17 00:00:00 2001 From: JulioAPeraza Date: Wed, 24 Apr 2024 16:49:51 -0400 Subject: [PATCH 2/4] Update test_results.py --- pymare/tests/test_results.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pymare/tests/test_results.py b/pymare/tests/test_results.py index 5b00742..284e45e 100644 --- a/pymare/tests/test_results.py +++ b/pymare/tests/test_results.py @@ -87,7 +87,10 @@ def test_combination_test_results_from_arrays(dataset): # fit overwrites dataset_ attribute with None assert fitted_estimator.dataset_ is None + # fit_dataset overwrites it with the Dataset + dataset.n = dataset.v # StoufferCombinationTest now maps n to w + dataset.v = None fitted_estimator.fit_dataset(dataset) assert isinstance(fitted_estimator.dataset_, Dataset) # fit sets it back to None From a7fd12037e2bd311c6957fcad6c3ffbde8bba865 Mon Sep 17 00:00:00 2001 From: JulioAPeraza Date: Wed, 24 Apr 2024 16:50:44 -0400 Subject: [PATCH 3/4] Switch to macos-12 per https://github.com/actions/setup-python/issues/850 --- .github/workflows/testing.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index c792186..8d9804f 100644 --- a/.github/workflows/testing.yml +++ b/.github/workflows/testing.yml @@ -38,7 +38,7 @@ jobs: strategy: fail-fast: false matrix: - os: ["ubuntu-latest", "macos-latest"] + os: ["ubuntu-latest", "macos-12"] python-version: ["3.8", "3.9", "3.10", "3.11"] name: ${{ matrix.os }} with Python ${{ matrix.python-version }} From 306dd7b2a42f0d2564c766b1bf517782a2c36d8c Mon Sep 17 00:00:00 2001 From: JulioAPeraza Date: Wed, 24 Apr 2024 16:57:16 -0400 Subject: [PATCH 4/4] Update test_results.py --- pymare/tests/test_results.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/pymare/tests/test_results.py b/pymare/tests/test_results.py index 284e45e..edaa92e 100644 --- a/pymare/tests/test_results.py +++ b/pymare/tests/test_results.py @@ -89,9 +89,7 @@ def test_combination_test_results_from_arrays(dataset): assert fitted_estimator.dataset_ is None # fit_dataset overwrites it with the Dataset - dataset.n = dataset.v # StoufferCombinationTest now maps n to w - dataset.v = None - fitted_estimator.fit_dataset(dataset) + fitted_estimator.fit_dataset(Dataset(dataset.y)) assert isinstance(fitted_estimator.dataset_, Dataset) # fit sets it back to None fitted_estimator.fit(z=dataset.y)