From 624fba0b04d5bcd46d6141d6fdf639bd376f3a17 Mon Sep 17 00:00:00 2001 From: JulioAPeraza Date: Tue, 2 Apr 2024 11:39:32 -0400 Subject: [PATCH 01/11] Add correction term for multiple contrasts in Stouffer's combination test --- pymare/estimators/combination.py | 74 +++++++++++++++++++++++--- pymare/tests/test_combination_tests.py | 38 +++++++++++++ 2 files changed, 106 insertions(+), 6 deletions(-) diff --git a/pymare/estimators/combination.py b/pymare/estimators/combination.py index f2ffabc..3d198e7 100644 --- a/pymare/estimators/combination.py +++ b/pymare/estimators/combination.py @@ -1,4 +1,5 @@ """Estimators for combination (p/z) tests.""" + import warnings from abc import abstractmethod @@ -110,17 +111,78 @@ class StoufferCombinationTest(CombinationTest): """ # Maps Dataset attributes onto fit() args; see BaseEstimator for details. - _dataset_attr_map = {"z": "y", "w": "v"} - - def fit(self, z, w=None): + _dataset_attr_map = {"z": "y", "w": "n", "g": "v"} + + def _inflation_term(self, z, w, g): + """Calculate the variance inflation term for each group. + + This term is used to adjust the variance of the combined z-scorem when + multiple sample from the same group are present. + + Parameters + ---------- + z : array-like + Array of z-values. + w : array-like + Array of weights. + g : array-like + Array of group labels. + + Returns + ------- + sigma : float + The variance inflation term. + """ + # Centering + all_samples_same = np.all(np.equal(z, z[0]), axis=0).all() + if not all_samples_same: + # Only center if the samples are not all the same, to prevent division by zero + # when calculating the correlation matrix. + mean_sample = z.mean(0) + z = z - mean_sample # This centering is problematic for N=2 + + # Use the value for one voxel, as all voxels have the same group and weight + groups = g[:, 0] + weights = w[:, 0] + + # Loop over groups + unique_groups = np.unique(groups) + + sigma = 0 + for group in unique_groups: + group_indices = np.where(groups == group)[0] + group_z = z[group_indices] + + # For groups with only one sample the contribution to the summand is 0 + n_samples = len(group_indices) + if n_samples < 2: + continue + + # Calculate the within group correlation matrix and select the non-diagonal elements + corr = np.corrcoef(group_z, rowvar=True) + upper_indices = np.triu_indices(n_samples, k=1) + non_diag_corr = 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): """Fit the estimator to z-values, optionally with weights.""" - return super().fit(z, w=w) + return super().fit(z, w=w, g=g) - def p_value(self, z, w=None): + def p_value(self, z, w=None, g=None): """Calculate p-values.""" if w is None: w = np.ones_like(z) - cz = (z * w).sum(0) / np.sqrt((w**2).sum(0)) + + # Calculate the variance inflation term + sigma = self._inflation_term(z, w, g) if g is not None else 0 + + variance = (w**2).sum(0) + sigma + + cz = (z * w).sum(0) / np.sqrt(variance) return ss.norm.sf(cz) diff --git a/pymare/tests/test_combination_tests.py b/pymare/tests/test_combination_tests.py index feacbc9..1b123ff 100644 --- a/pymare/tests/test_combination_tests.py +++ b/pymare/tests/test_combination_tests.py @@ -1,4 +1,5 @@ """Tests for pymare.estimators.combination.""" + import numpy as np import pytest import scipy.stats as ss @@ -33,6 +34,43 @@ def test_combination_test(Cls, data, mode, expected): assert np.allclose(z, expected, atol=1e-5) +def test_stouffer_adjusted(): + """Test StoufferCombinationTest with weights and groups.""" + # Test with weights and groups + data = np.array( + [ + [2.1, 0.7, -0.2, 4.1, 3.8], + [1.1, 0.2, 0.4, 1.3, 1.5], + [-0.6, -1.6, -2.3, -0.8, -4.0], + [2.5, 1.7, 2.1, 2.3, 2.5], + [3.1, 2.7, 3.1, 3.3, 3.5], + [3.6, 3.2, 3.6, 3.8, 4.0], + ] + ) + weights = np.tile(np.array([4, 3, 4, 10, 15, 10]), (data.shape[1], 1)).T + groups = np.tile(np.array([0, 0, 1, 2, 2, 2]), (data.shape[1], 1)).T + + results = StoufferCombinationTest("directed").fit(z=data, w=weights, g=groups).params_ + z = ss.norm.isf(results["p"]) + + z_expected = np.array([5.00088912, 3.70356943, 4.05465924, 5.4633001, 5.18927878]) + assert np.allclose(z, z_expected, atol=1e-5) + + # Test with weights and no groups. Limiting cases. + # Limiting case 1: all correlations are one. + n_maps_l1 = 5 + common_sample = np.array([2.1, 0.7, -0.2]) + data_l1 = np.tile(common_sample, (n_maps_l1, 1)) + groups_l1 = np.tile(np.array([0, 0, 0, 0, 0]), (data_l1.shape[1], 1)).T + + results_l1 = StoufferCombinationTest("directed").fit(z=data_l1, g=groups_l1).params_ + z_l1 = ss.norm.isf(results_l1["p"]) + + 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) + + @pytest.mark.parametrize("Cls,data,mode,expected", _params) def test_combination_test_from_dataset(Cls, data, mode, expected): """Test CombinationTest Estimators with PyMARE Datasets.""" From aa389e77eed611787159543a45ab60fc15469cfd Mon Sep 17 00:00:00 2001 From: JulioAPeraza Date: Tue, 2 Apr 2024 12:19:24 -0400 Subject: [PATCH 02/11] Update RTD yml --- .readthedocs.yml | 5 +++++ pymare/estimators/combination.py | 4 ++-- pymare/tests/test_combination_tests.py | 20 ++++++++++---------- 3 files changed, 17 insertions(+), 12 deletions(-) diff --git a/.readthedocs.yml b/.readthedocs.yml index f994a7e..4e696ba 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -5,6 +5,11 @@ # Required version: 2 +build: + os: "ubuntu-22.04" + tools: + python: "3.9" + # Build documentation in the docs/ directory with Sphinx sphinx: configuration: docs/conf.py diff --git a/pymare/estimators/combination.py b/pymare/estimators/combination.py index 3d198e7..61f18e9 100644 --- a/pymare/estimators/combination.py +++ b/pymare/estimators/combination.py @@ -116,8 +116,8 @@ class StoufferCombinationTest(CombinationTest): def _inflation_term(self, z, w, g): """Calculate the variance inflation term for each group. - This term is used to adjust the variance of the combined z-scorem when - multiple sample from the same group are present. + This term is used to adjust the variance of the combined z-score when + multiple sample come from the same group are present. Parameters ---------- diff --git a/pymare/tests/test_combination_tests.py b/pymare/tests/test_combination_tests.py index 1b123ff..9f5f7e1 100644 --- a/pymare/tests/test_combination_tests.py +++ b/pymare/tests/test_combination_tests.py @@ -34,6 +34,16 @@ def test_combination_test(Cls, data, mode, expected): assert np.allclose(z, expected, atol=1e-5) +@pytest.mark.parametrize("Cls,data,mode,expected", _params) +def test_combination_test_from_dataset(Cls, data, mode, expected): + """Test CombinationTest Estimators with PyMARE Datasets.""" + dset = Dataset(y=data) + est = Cls(mode).fit_dataset(dset) + results = est.summary() + z = ss.norm.isf(results.p) + assert np.allclose(z, expected, atol=1e-5) + + def test_stouffer_adjusted(): """Test StoufferCombinationTest with weights and groups.""" # Test with weights and groups @@ -69,13 +79,3 @@ 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) - - -@pytest.mark.parametrize("Cls,data,mode,expected", _params) -def test_combination_test_from_dataset(Cls, data, mode, expected): - """Test CombinationTest Estimators with PyMARE Datasets.""" - dset = Dataset(y=data) - est = Cls(mode).fit_dataset(dset) - results = est.summary() - z = ss.norm.isf(results.p) - assert np.allclose(z, expected, atol=1e-5) From 3c8e58aa397645f220f48798a68443045261f939 Mon Sep 17 00:00:00 2001 From: JulioAPeraza Date: Tue, 2 Apr 2024 12:22:42 -0400 Subject: [PATCH 03/11] Update .readthedocs.yml --- .readthedocs.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.readthedocs.yml b/.readthedocs.yml index 4e696ba..616b148 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -15,7 +15,6 @@ sphinx: configuration: docs/conf.py python: - version: 3.7 install: - method: pip path: . From 1973650dfaf91b77887a8d5ecc29051d22f47215 Mon Sep 17 00:00:00 2001 From: JulioAPeraza Date: Tue, 2 Apr 2024 12:22:44 -0400 Subject: [PATCH 04/11] Update combination.py --- pymare/estimators/combination.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/pymare/estimators/combination.py b/pymare/estimators/combination.py index 61f18e9..8bcaf5b 100644 --- a/pymare/estimators/combination.py +++ b/pymare/estimators/combination.py @@ -134,12 +134,11 @@ def _inflation_term(self, z, w, g): The variance inflation term. """ # Centering + # Only center if the samples are not all the same, to prevent division by zero + # when calculating the correlation matrix. + # This centering is problematic for N=2 all_samples_same = np.all(np.equal(z, z[0]), axis=0).all() - if not all_samples_same: - # Only center if the samples are not all the same, to prevent division by zero - # when calculating the correlation matrix. - mean_sample = z.mean(0) - z = z - mean_sample # This centering is problematic for N=2 + z = z if all_samples_same else z - z.mean(0) # Use the value for one voxel, as all voxels have the same group and weight groups = g[:, 0] From 919068490c582549bda4ad2e0a3d5d0c70a9ba33 Mon Sep 17 00:00:00 2001 From: JulioAPeraza Date: Tue, 2 Apr 2024 12:23:53 -0400 Subject: [PATCH 05/11] Update .readthedocs.yml --- .readthedocs.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.readthedocs.yml b/.readthedocs.yml index 616b148..1d35214 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -20,4 +20,3 @@ python: path: . extra_requirements: - doc - system_packages: true From c76a18878ab3465497a8d4a55beb5221ca4fe579 Mon Sep 17 00:00:00 2001 From: JulioAPeraza Date: Tue, 2 Apr 2024 12:30:45 -0400 Subject: [PATCH 06/11] Update setup.cfg --- setup.cfg | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.cfg b/setup.cfg index 17959f6..a74a38a 100644 --- a/setup.cfg +++ b/setup.cfg @@ -57,7 +57,7 @@ doc = sphinx>=3.5 sphinx-argparse sphinx-copybutton - sphinx_gallery==0.10.1 + sphinx_gallery sphinx_rtd_theme sphinxcontrib-bibtex tests = From 23b1e2e758714baf8906f41511d0227391319360 Mon Sep 17 00:00:00 2001 From: JulioAPeraza Date: Tue, 2 Apr 2024 12:49:38 -0400 Subject: [PATCH 07/11] Update testing.yml --- .github/workflows/testing.yml | 45 ++++++++++++++++------------------- 1 file changed, 21 insertions(+), 24 deletions(-) diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index 050df84..db78a6c 100644 --- a/.github/workflows/testing.yml +++ b/.github/workflows/testing.yml @@ -6,10 +6,7 @@ on: - "master" pull_request: branches: - - '*' - schedule: - # Run tests every Sunday at 12am - - cron: '0 0 * * 0' + - "*" concurrency: group: environment-${{ github.ref }} @@ -28,46 +25,46 @@ jobs: - id: result_step uses: mstachniuk/ci-skip@master with: - commit-filter: '[skip ci];[ci skip];[skip github]' - commit-filter-separator: ';' + commit-filter: "[skip ci];[ci skip];[skip github]" + commit-filter-separator: ";" run_tests: needs: check_skip if: ${{ needs.check_skip.outputs.skip == 'false' }} runs-on: ${{ matrix.os }} strategy: - fail-fast: false - matrix: - os: ["ubuntu-latest", "macos-latest"] - python-version: ["3.6", "3.7", "3.8", "3.9"] - include: - # ubuntu-20.04 is used only to test python 3.6 - - os: "ubuntu-20.04" - python-version: "3.6" - exclude: - # ubuntu-latest does not support python 3.6 - - os: "ubuntu-latest" - python-version: "3.6" + fail-fast: false + matrix: + os: ["ubuntu-latest", "macos-latest"] + python-version: ["3.6", "3.7", "3.8", "3.9"] + include: + # ubuntu-20.04 is used only to test python 3.6 + - os: "ubuntu-20.04" + python-version: "3.6" + exclude: + # ubuntu-latest does not support python 3.6 + - os: "ubuntu-latest" + python-version: "3.6" name: ${{ matrix.os }} with Python ${{ matrix.python-version }} defaults: run: shell: bash steps: - uses: actions/checkout@v2 - - name: 'Set up python' + - name: "Set up python" uses: actions/setup-python@v2 with: - python-version: ${{ matrix.python-version }} - - name: 'Display Python version' + python-version: ${{ matrix.python-version }} + - name: "Display Python version" shell: bash {0} run: python -c "import sys; print(sys.version)" - - name: 'Install PyMARE' + - name: "Install PyMARE" shell: bash {0} run: pip install -e .[tests,stan] - - name: 'Run tests' + - name: "Run tests" shell: bash {0} run: python -m pytest --pyargs pymare --cov=pymare - - name: 'Upload coverage to CodeCov' + - name: "Upload coverage to CodeCov" uses: codecov/codecov-action@v1 if: success() From cd291f1b2d8564c5979621babbf9851845e4f480 Mon Sep 17 00:00:00 2001 From: JulioAPeraza Date: Tue, 2 Apr 2024 12:52:21 -0400 Subject: [PATCH 08/11] Update testing.yml --- .github/workflows/testing.yml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index db78a6c..53e6630 100644 --- a/.github/workflows/testing.yml +++ b/.github/workflows/testing.yml @@ -7,6 +7,9 @@ on: pull_request: branches: - "*" + schedule: + # Run tests every Sunday at 12am + - cron: "0 0 * * 0" concurrency: group: environment-${{ github.ref }} From 91880ee179d7c465e6f16d4246c4d123dec7770e Mon Sep 17 00:00:00 2001 From: JulioAPeraza Date: Tue, 2 Apr 2024 12:56:38 -0400 Subject: [PATCH 09/11] Run black --- pymare/__init__.py | 1 + pymare/datasets/__init__.py | 1 + pymare/datasets/metadat.py | 1 + pymare/effectsize/__init__.py | 1 + pymare/effectsize/expressions.py | 1 + pymare/estimators/__init__.py | 1 + pymare/results.py | 1 + pymare/tests/conftest.py | 1 + pymare/tests/test_core.py | 1 + pymare/tests/test_datasets.py | 1 + pymare/tests/test_effectsize_base.py | 1 + pymare/tests/test_effectsize_expressions.py | 1 + pymare/tests/test_estimators.py | 1 + pymare/tests/test_results.py | 1 + pymare/tests/test_stan_estimators.py | 1 + pymare/tests/test_stats.py | 1 + pymare/tests/test_utils.py | 1 + pymare/utils.py | 1 + 18 files changed, 18 insertions(+) diff --git a/pymare/__init__.py b/pymare/__init__.py index 6f112e9..8127a73 100644 --- a/pymare/__init__.py +++ b/pymare/__init__.py @@ -1,4 +1,5 @@ """PyMARE: Python Meta-Analysis & Regression Engine.""" + import sys import warnings diff --git a/pymare/datasets/__init__.py b/pymare/datasets/__init__.py index d4c0a6e..01d3f4d 100644 --- a/pymare/datasets/__init__.py +++ b/pymare/datasets/__init__.py @@ -1,4 +1,5 @@ """Open meta-analytic datasets.""" + from .metadat import michael2013 __all__ = [ diff --git a/pymare/datasets/metadat.py b/pymare/datasets/metadat.py index f3fcfd5..4592520 100644 --- a/pymare/datasets/metadat.py +++ b/pymare/datasets/metadat.py @@ -1,4 +1,5 @@ """Datasets from metadat.""" + import json import os.path as op diff --git a/pymare/effectsize/__init__.py b/pymare/effectsize/__init__.py index d0a5424..7dbfd8e 100644 --- a/pymare/effectsize/__init__.py +++ b/pymare/effectsize/__init__.py @@ -1,4 +1,5 @@ """Tools for converting between effect-size measures.""" + from .base import ( OneSampleEffectSizeConverter, TwoSampleEffectSizeConverter, diff --git a/pymare/effectsize/expressions.py b/pymare/effectsize/expressions.py index 2a56a61..6cc1eeb 100644 --- a/pymare/effectsize/expressions.py +++ b/pymare/effectsize/expressions.py @@ -1,4 +1,5 @@ """Statistical expressions.""" + import json from collections import defaultdict from itertools import chain diff --git a/pymare/estimators/__init__.py b/pymare/estimators/__init__.py index c964654..5b48193 100644 --- a/pymare/estimators/__init__.py +++ b/pymare/estimators/__init__.py @@ -1,4 +1,5 @@ """Estimators for meta-analyses and meta-regressions.""" + from .combination import FisherCombinationTest, StoufferCombinationTest from .estimators import ( DerSimonianLaird, diff --git a/pymare/results.py b/pymare/results.py index b84b1f6..31c7d3e 100644 --- a/pymare/results.py +++ b/pymare/results.py @@ -1,4 +1,5 @@ """Tools for representing and manipulating meta-regression results.""" + import itertools from functools import lru_cache from inspect import getfullargspec diff --git a/pymare/tests/conftest.py b/pymare/tests/conftest.py index 6a3e988..cce93c8 100644 --- a/pymare/tests/conftest.py +++ b/pymare/tests/conftest.py @@ -1,4 +1,5 @@ """Data for tests.""" + import numpy as np import pytest diff --git a/pymare/tests/test_core.py b/pymare/tests/test_core.py index 8ce58b1..12fd1b4 100644 --- a/pymare/tests/test_core.py +++ b/pymare/tests/test_core.py @@ -1,4 +1,5 @@ """Tests for pymare.core.""" + import numpy as np import pandas as pd import pytest diff --git a/pymare/tests/test_datasets.py b/pymare/tests/test_datasets.py index fcf312f..26859e1 100644 --- a/pymare/tests/test_datasets.py +++ b/pymare/tests/test_datasets.py @@ -1,4 +1,5 @@ """Tests for the pymare.datasets module.""" + import pandas as pd from pymare import datasets diff --git a/pymare/tests/test_effectsize_base.py b/pymare/tests/test_effectsize_base.py index 0616fa8..515b8cb 100644 --- a/pymare/tests/test_effectsize_base.py +++ b/pymare/tests/test_effectsize_base.py @@ -1,4 +1,5 @@ """Tests for pymare.effectsize.base.""" + import numpy as np import pandas as pd import pytest diff --git a/pymare/tests/test_effectsize_expressions.py b/pymare/tests/test_effectsize_expressions.py index b28a675..9c33c9d 100644 --- a/pymare/tests/test_effectsize_expressions.py +++ b/pymare/tests/test_effectsize_expressions.py @@ -1,4 +1,5 @@ """Tests for pymare.effectsize.expressions.""" + import pytest from sympy import Symbol from sympy.core.sympify import SympifyError diff --git a/pymare/tests/test_estimators.py b/pymare/tests/test_estimators.py index 339e163..8b6fece 100644 --- a/pymare/tests/test_estimators.py +++ b/pymare/tests/test_estimators.py @@ -1,4 +1,5 @@ """Tests for pymare.estimators.estimators.""" + import numpy as np import pytest diff --git a/pymare/tests/test_results.py b/pymare/tests/test_results.py index 2baf3bd..5b00742 100644 --- a/pymare/tests/test_results.py +++ b/pymare/tests/test_results.py @@ -1,4 +1,5 @@ """Tests for pymare.results.""" + import numpy as np import pytest diff --git a/pymare/tests/test_stan_estimators.py b/pymare/tests/test_stan_estimators.py index c685ecd..71efc28 100644 --- a/pymare/tests/test_stan_estimators.py +++ b/pymare/tests/test_stan_estimators.py @@ -1,4 +1,5 @@ """Tests for estimators that use stan.""" + import pytest from pymare.estimators import StanMetaRegression diff --git a/pymare/tests/test_stats.py b/pymare/tests/test_stats.py index c350f3f..f01d434 100644 --- a/pymare/tests/test_stats.py +++ b/pymare/tests/test_stats.py @@ -1,4 +1,5 @@ """Tests for pymare.stats.""" + from pymare import stats diff --git a/pymare/tests/test_utils.py b/pymare/tests/test_utils.py index 474b8f8..7ff794c 100644 --- a/pymare/tests/test_utils.py +++ b/pymare/tests/test_utils.py @@ -1,4 +1,5 @@ """Tests for pymare.utils.""" + import os.path as op import numpy as np diff --git a/pymare/utils.py b/pymare/utils.py index b159a48..d03ea9a 100644 --- a/pymare/utils.py +++ b/pymare/utils.py @@ -1,4 +1,5 @@ """Miscellaneous utility functions.""" + import os.path as op import numpy as np From a718ac4c485fced7c67f414a2e394c7509edacc5 Mon Sep 17 00:00:00 2001 From: JulioAPeraza Date: Tue, 2 Apr 2024 15:03:34 -0400 Subject: [PATCH 10/11] Make sure solutions and symbols match --- pymare/effectsize/base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pymare/effectsize/base.py b/pymare/effectsize/base.py index ac15b9a..68d8642 100644 --- a/pymare/effectsize/base.py +++ b/pymare/effectsize/base.py @@ -68,7 +68,7 @@ def solve_system(system, known_vars=None): # solver will return a dict if there's only one non-dummy expression if isinstance(solutions, dict): - solutions = [list(solutions.values())] + solutions = [[solutions[s] for s in symbols]] # Prepare the dummy list and data args in a fixed order dummy_list = list(dummies) From 539974edea25fe4249b3584eef0b4a74c3bdd391 Mon Sep 17 00:00:00 2001 From: JulioAPeraza Date: Fri, 5 Apr 2024 11:59:45 -0400 Subject: [PATCH 11/11] Update combination.py --- pymare/estimators/combination.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/pymare/estimators/combination.py b/pymare/estimators/combination.py index 8bcaf5b..6f9f41a 100644 --- a/pymare/estimators/combination.py +++ b/pymare/estimators/combination.py @@ -117,15 +117,15 @@ def _inflation_term(self, z, w, g): """Calculate the variance inflation term for each group. This term is used to adjust the variance of the combined z-score when - multiple sample come from the same group are present. + multiple sample come from the same study. Parameters ---------- - z : array-like + z : :obj:`numpy.ndarray` of shape (n, d) Array of z-values. - w : array-like + w : :obj:`numpy.ndarray` of shape (n, d) Array of weights. - g : array-like + g : :obj:`numpy.ndarray` of shape (n, d) Array of group labels. Returns @@ -133,14 +133,13 @@ def _inflation_term(self, z, w, g): sigma : float The variance inflation term. """ - # Centering # Only center if the samples are not all the same, to prevent division by zero # when calculating the correlation matrix. # This centering is problematic for N=2 all_samples_same = np.all(np.equal(z, z[0]), axis=0).all() z = z if all_samples_same else z - z.mean(0) - # Use the value for one voxel, as all voxels have the same group and weight + # Use the value from one feature, as all features have the same groups and weights groups = g[:, 0] weights = w[:, 0] @@ -157,7 +156,7 @@ def _inflation_term(self, z, w, g): if n_samples < 2: continue - # Calculate the within group correlation matrix and select the non-diagonal elements + # Calculate the within group correlation matrix and sum the non-diagonal elements corr = np.corrcoef(group_z, rowvar=True) upper_indices = np.triu_indices(n_samples, k=1) non_diag_corr = corr[upper_indices] @@ -168,7 +167,7 @@ def _inflation_term(self, z, w, g): return sigma def fit(self, z, w=None, g=None): - """Fit the estimator to z-values, optionally with weights.""" + """Fit the estimator to z-values, optionally with weights and groups.""" return super().fit(z, w=w, g=g) def p_value(self, z, w=None, g=None): @@ -176,9 +175,10 @@ def p_value(self, z, w=None, g=None): if w is None: w = np.ones_like(z) - # Calculate the variance inflation term + # 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 + # The sum of diagonal elements of sigma is given by (w**2).sum(0). variance = (w**2).sum(0) + sigma cz = (z * w).sum(0) / np.sqrt(variance)