diff --git a/.github/workflows/stan-validation.yml b/.github/workflows/stan-validation.yml new file mode 100644 index 0000000..ae5bd33 --- /dev/null +++ b/.github/workflows/stan-validation.yml @@ -0,0 +1,115 @@ +name: "Validate the Stan model" + +on: + # Not on pull_request: the grid is ~1400 fits and takes about ten minutes, + # which is too slow to sit in front of every change. The fast checks that do + # run per pull request are the "stan" marked tests in the Run Tests workflow, + # plus the two tests there that hold the *recorded* results in + # pymare/tests/data/stan_validation.json to these same thresholds. This + # workflow is what re-measures them. + push: + branches: + - "master" + paths: + - "pymare/estimators/stan/**" + - "pymare/estimators/estimators.py" + - "validation/stan/**" + - ".github/workflows/stan-validation.yml" + schedule: + # Monthly, to catch the model rotting against a new CmdStan or ArviZ rather + # than against a change to PyMARE. + - cron: "0 0 1 * *" + workflow_dispatch: + inputs: + replications: + description: "Replications per design cell (at least 100)" + required: false + default: "100" + +permissions: + contents: read + +concurrency: + group: stan-validation-${{ github.ref }} + cancel-in-progress: true + +env: + CMDSTAN_VERSION: "2.36.0" + REPLICATIONS: ${{ github.event.inputs.replications || '100' }} + +jobs: + validate: + name: Bias and credible-interval coverage + runs-on: ubuntu-latest + defaults: + run: + shell: bash + steps: + - uses: actions/checkout@v4 + - name: "Set up python" + uses: actions/setup-python@v5 + with: + python-version: "3.11" + - name: "Install PyMARE with the stan extra" + run: | + python -m pip install --progress-bar off --upgrade pip setuptools wheel + python -m pip install -e .[tests,stan] + - name: "Cache CmdStan" + id: cache_cmdstan + uses: actions/cache@v4 + with: + path: ~/.cmdstan + key: cmdstan-${{ runner.os }}-${{ runner.arch }}-${{ env.CMDSTAN_VERSION }} + - name: "Install CmdStan" + if: steps.cache_cmdstan.outputs.cache-hit != 'true' + run: python -m cmdstanpy.install_cmdstan --version "${CMDSTAN_VERSION}" --cores 2 + + # Regenerates pymare/tests/data/stan_validation.json and fails if any cell + # misses STAN_VALIDATION_THRESHOLDS. The numbers are stochastic, so unlike + # the robumeta reference this cannot require the file to be unchanged -- + # a correct model produces slightly different numbers every run. The + # thresholds are the claim; the file is the record of it. + - name: "Re-measure the model and enforce the thresholds" + id: measure + run: | + python validation/stan/simulate.py \ + --replications "${REPLICATIONS}" \ + --jobs 2 \ + --check \ + | tee measured.txt + + - name: "Report what was measured" + if: always() + run: | + { + echo "## Stan model validation" + echo + echo "\`${REPLICATIONS}\` replications per design cell, CmdStan \`${CMDSTAN_VERSION}\`." + echo + echo '```' + tail -n 20 measured.txt 2>/dev/null || echo "the run produced no output" + echo '```' + echo + if git diff --quiet -- pymare/tests/data/stan_validation.json; then + echo "The recorded results did not move." + else + echo "
Change against the recorded results" + echo + echo '```diff' + git diff -- pymare/tests/data/stan_validation.json + echo '```' + echo + echo "
" + echo + echo "These numbers are stochastic, so movement is expected. What matters" + echo "is that every cell still clears the thresholds, which the step above" + echo "enforces. Commit the regenerated file to refresh the record." + fi + } >> "$GITHUB_STEP_SUMMARY" + + - name: "Upload the regenerated results" + if: always() + uses: actions/upload-artifact@v4 + with: + name: stan-validation + path: pymare/tests/data/stan_validation.json diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index aee89d9..9894c44 100644 --- a/.github/workflows/testing.yml +++ b/.github/workflows/testing.yml @@ -11,12 +11,6 @@ on: # Run tests every Sunday at 12am, to catch breakage from new dependency # releases rather than from a change to PyMARE. - cron: "0 0 * * 0" - # GitHub turns off a workflow with a schedule trigger once the repository has - # been quiet for 60 days, and a disabled workflow stops answering push and - # pull_request too -- which is how this one sat switched off from 2025-06-28 - # until it was re-enabled by hand. Being able to dispatch it means the next - # time that happens, re-enabling is enough to get a run without pushing a - # commit to prove it works. workflow_dispatch: permissions: @@ -27,13 +21,12 @@ concurrency: cancel-in-progress: true env: - # One place to change how the suite is invoked. --cov-append lets each job - # write a coverage file that upload_to_codecov merges. PYTEST_COMMON_ARGS: --cov-append --cov-report=xml --cov=pymare - # Everything except the Stan estimator, whose model has to be compiled and so - # costs minutes rather than seconds. It gets the job below to itself. + # Everything except the tests that sample, whose Stan model has to be compiled + # and so costs minutes rather than seconds. PYTEST_UNIT_MARKERS: not stan PYTEST_STAN_MARKERS: stan + CMDSTAN_VERSION: "2.36.0" jobs: # Determine if tests should be run based on commit message. @@ -86,10 +79,18 @@ jobs: python-version: ${{ matrix.python-version }} - name: "Display Python version" run: python -c "import sys; print(sys.version)" + # The stan extra, not just tests: it brings in ArviZ, and the results + # container is ArviZ-only code that nothing here could otherwise execute. + # Without it the whole of BayesianMetaRegressionResults is reachable only + # from the single Stan job below, so its coverage -- and the ArviZ 0.x + # versus 1.x handling in particular -- would rest on one job on one + # Python. cmdstanpy comes along too but stays idle: it is pure Python, + # CmdStan is not installed here, and the tests that sample are excluded by + # the marker filter regardless. - name: "Install PyMARE" run: | python -m pip install --progress-bar off --upgrade pip setuptools wheel - python -m pip install -e .[tests] + python -m pip install -e .[tests,stan] - name: "Run tests" run: python -m pytest -m "${PYTEST_UNIT_MARKERS}" ${PYTEST_COMMON_ARGS} - name: "Upload coverage" @@ -104,6 +105,8 @@ jobs: needs: check_skip if: ${{ needs.check_skip.outputs.skip == 'false' }} runs-on: ubuntu-latest + env: + PYMARE_REQUIRE_CMDSTAN: "1" defaults: run: shell: bash @@ -117,6 +120,41 @@ jobs: run: | python -m pip install --progress-bar off --upgrade pip setuptools wheel python -m pip install -e .[tests,stan] + - name: "Cache CmdStan" + id: cache_cmdstan + uses: actions/cache@v4 + with: + path: ~/.cmdstan + # runner.arch is part of the key because CmdStan is a native build: + # restoring an x86_64 tree onto an arm64 runner fails confusingly. + key: cmdstan-${{ runner.os }}-${{ runner.arch }}-${{ env.CMDSTAN_VERSION }} + - name: "Install CmdStan" + if: steps.cache_cmdstan.outputs.cache-hit != 'true' + run: python -m cmdstanpy.install_cmdstan --version "${CMDSTAN_VERSION}" --cores 2 + - name: "Cache the compiled model" + uses: actions/cache@v4 + with: + # CmdStanPy compiles beside the .stan source, not under ~/.cmdstan, so + # this needs a cache entry of its own. The key covers all four inputs + # the executable depends on. + path: pymare/estimators/stan + key: >- + stanexe-${{ runner.os }}-${{ runner.arch }}-${{ env.CMDSTAN_VERSION }}-${{ + hashFiles('pymare/estimators/stan/*.stan') }} + - name: "Check the Stan program with pedantic mode" + # Informational. Pedantic mode flags a parameter with no prior, which is + # what the old model's tau2 was, but it also flags beta's deliberately + # flat prior, so this reports rather than gates. + # + # Invoke stanc directly rather than cmdstanpy.compile_stan_file: that + # skips its work when an up-to-date executable exists, so with the model + # cache restored above it would silently check nothing on most runs. + # stanc only translates Stan to C++, which takes a fifth of a second and + # leaves the cached executable alone, so it always runs. + continue-on-error: true + run: | + STANC="$(python -c "import cmdstanpy, os; print(os.path.join(cmdstanpy.cmdstan_path(), 'bin', 'stanc'))")" + "${STANC}" --warn-pedantic --o=/dev/null pymare/estimators/stan/meta_regression.stan - name: "Run tests" run: python -m pytest -m "${PYTEST_STAN_MARKERS}" ${PYTEST_COMMON_ARGS} - name: "Upload coverage" diff --git a/.gitignore b/.gitignore index c7436d9..c50c9fb 100644 --- a/.gitignore +++ b/.gitignore @@ -119,3 +119,10 @@ fmriprep.sqlite .vscode # asv benchmark environments and results .asv/ + +# CmdStan build artifacts. CmdStanModel compiles the model in place, dropping +# the executable and its intermediates next to the .stan source. +pymare/estimators/stan/meta_regression +pymare/estimators/stan/meta_regression.exe +pymare/estimators/stan/*.hpp +pymare/estimators/stan/*.o diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 4c5a759..ff6361f 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -98,15 +98,42 @@ environment does not have: | Target | What it runs | Needs | | --- | --- | --- | -| `make unittest` | everything except the Stan tests | nothing extra | -| `make test_stan` | the Stan estimator tests | `pip install -e .[stan]` | +| `make unittest` | everything except the Stan sampling tests | nothing extra | +| `make test_stan` | the Stan sampling tests | `pip install -e .[stan]`, then `make install_cmdstan` | | `make test_robumeta` | the robumeta alignment tests | nothing extra | | `make check_robumeta_alignment` | regenerates the robumeta reference values | Docker | +| `make validate_stan` | re-measures the Stan model's bias and coverage (~10 min) | the same as `test_stan` | | `make lint` | flake8 over `pymare` and `benchmarks` | nothing extra | Each of these has a GitHub Actions job behind it, so a target that passes locally is the same check that runs on your pull request. +`make test_stan` needs two installation steps rather than one: the `stan` extra +brings in cmdstanpy, but CmdStan itself is a C++ build rather than a Python +package, so `make install_cmdstan` fetches and builds it. That takes several +minutes the first time and nothing thereafter. + +**Those tests skip locally when CmdStan is missing, but fail in CI.** The +asymmetry is deliberate: a contributor without CmdStan should not see red, but a +skip is indistinguishable from a pass in a CI log, which is how this job once +reported success while running none of the tests it existed to run. The Stan job +sets `PYMARE_REQUIRE_CMDSTAN=1`, and the `pytest_collection_modifyitems` hook in +`pymare/tests/conftest.py` fails the run at collection wherever that is set and +CmdStan is missing. + +Only the tests that sample are marked `stan`. Those that check how PyMARE's +inputs are translated into Stan's data block need neither cmdstanpy nor CmdStan, +so they are unmarked and run in the ordinary unit job on every platform. + +The model's accuracy is measured separately by `make validate_stan`, which takes +about ten minutes and so is not run per pull request. It reports bias and +credible-interval coverage across a grid of designs, records them in +`pymare/tests/data/stan_validation.json`, and fails if any design cell misses +`pymare.tests.utils.STAN_VALIDATION_THRESHOLDS`. Two tests hold that recorded +file to the same thresholds on every run, and the `Validate the Stan model` +workflow re-measures on a schedule. See `validation/stan/README.md` for the +arrangement and the measurements. + ### Alignment with robumeta `pymare/tests/test_robumeta_alignment.py` pins PyMARE's correlated-effects model diff --git a/MANIFEST.in b/MANIFEST.in index 9311739..32c0596 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,2 +1,3 @@ include versioneer.py include pymare/_version.py +recursive-include pymare *.stan diff --git a/Makefile b/Makefile index 93ac3cf..0146e03 100644 --- a/Makefile +++ b/Makefile @@ -1,4 +1,5 @@ -.PHONY: all_tests benchmark check_robumeta_alignment help lint test_robumeta test_stan unittest +.PHONY: all_tests benchmark check_robumeta_alignment help install_cmdstan lint test_robumeta +.PHONY: test_stan unittest validate_stan # --cov-append matches what CI does, so a local run of two targets in a row # reports their combined coverage rather than only the last one's. @@ -9,10 +10,12 @@ all_tests: lint unittest test_stan test_robumeta help: @echo "Please use 'make ' where is one of:" @echo " lint to run flake8 over pymare and the benchmarks" - @echo " unittest to run every test except the Stan ones" - @echo " test_stan to run the Stan estimator tests (needs the stan extra)" + @echo " unittest to run every test except the Stan sampling ones" + @echo " install_cmdstan to install the CmdStan that test_stan needs" + @echo " test_stan to run the Stan sampling tests (needs the stan extra and CmdStan)" @echo " test_robumeta to run the robumeta alignment tests" @echo " check_robumeta_alignment to regenerate the robumeta reference values (needs Docker)" + @echo " validate_stan to re-measure the Stan model's bias and coverage (~10 min)" @echo " benchmark to run the asv suite once in the current environment" @echo " all_tests to run lint and every test target" @@ -22,12 +25,26 @@ lint: unittest: @python -m pytest -m "not stan" $(PYTEST_COV) +# CmdStan is a C++ build rather than a Python package, so `pip install -e .[stan]` +# gets cmdstanpy but not the CmdStan it drives. This is the missing second step. +install_cmdstan: + @python -m cmdstanpy.install_cmdstan + +# Skips rather than fails when CmdStan is absent. CI sets PYMARE_REQUIRE_CMDSTAN=1 +# so that a job which is supposed to have it goes red instead of quietly empty. test_stan: @python -m pytest -m "stan" $(PYTEST_COV) test_robumeta: @python -m pytest -m "robumeta" $(PYTEST_COV) +# What the "Validate the Stan model" workflow runs. Regenerates +# pymare/tests/data/stan_validation.json and fails if any design cell misses the +# thresholds in pymare.tests.utils.STAN_VALIDATION_THRESHOLDS. Slow -- about ten +# minutes -- which is why it is not part of unittest or test_stan. +validate_stan: + @python validation/stan/simulate.py --check + # What the "Check robumeta alignment" workflow runs. Needs Docker, because the # reference values come from R. check_robumeta_alignment: diff --git a/docs/installation.rst b/docs/installation.rst index ffb58fe..a80b945 100644 --- a/docs/installation.rst +++ b/docs/installation.rst @@ -15,5 +15,29 @@ If you want to use the most up-to-date version, you can install from the ``maste pip install git+https://github.com/neurostuff/PyMARE.git -PyMARE requires Python >=3.8 and a number of packages. +PyMARE requires Python >=3.9 and a number of packages. For a complete list, please see ``setup.cfg``. + +Bayesian estimation with Stan +----------------------------- + +:class:`~pymare.estimators.StanMetaRegression` is optional, and needs two +installation steps rather than one: + +.. code-block:: bash + + pip install pymare[stan] + python -m cmdstanpy.install_cmdstan + +The first installs CmdStanPy. The second fetches and builds CmdStan itself, +which is a C++ program rather than a Python package and so needs a C++ toolchain +(``g++`` and ``make`` on Linux, the Command Line Tools on macOS, RTools on +Windows). It takes several minutes, once per machine. + +The Stan model is compiled the first time the estimator is fitted, which takes +roughly another minute. The compiled model is cached alongside the installed +package, so later fits and later processes reuse it. If PyMARE is installed +somewhere unwritable, the model is compiled into ``~/.pymare/stan`` instead and a +warning says so. + +Every other estimator in PyMARE is pure Python and needs none of this. diff --git a/examples/02_meta-analysis/plot_meta-analysis_walkthrough.py b/examples/02_meta-analysis/plot_meta-analysis_walkthrough.py index 2bce6de..c6d42ab 100644 --- a/examples/02_meta-analysis/plot_meta-analysis_walkthrough.py +++ b/examples/02_meta-analysis/plot_meta-analysis_walkthrough.py @@ -232,8 +232,12 @@ ############################################################################### # What about the Stan estimator? # ````````````````````````````````````````````````````````````````````````````` -# We're going to skip this one here because of how computationally intensive it -# is. +# We're going to skip this one here. It needs two things the documentation build +# does not have: CmdStan, which is a C++ build rather than a Python package, and +# a C++ toolchain to compile the model with. The compilation is a one-time cost +# per installation rather than a per-fit one, so it is much less of an obstacle +# in your own environment than it is here. Install it with +# ``pip install pymare[stan]`` followed by ``python -m cmdstanpy.install_cmdstan``. ############################################################################### # Let's check out our results! diff --git a/pymare/estimators/estimators.py b/pymare/estimators/estimators.py index 21fe3d4..36c1f0f 100644 --- a/pymare/estimators/estimators.py +++ b/pymare/estimators/estimators.py @@ -1,6 +1,8 @@ """Meta-regression estimator classes.""" -import sys +import os +import os.path as op +import shutil from abc import ABCMeta, abstractmethod from inspect import getfullargspec from warnings import warn @@ -46,18 +48,52 @@ def check(self, name, value): class Interval: - """Constrain a parameter to a closed numeric interval.""" + """Constrain a parameter to a numeric interval. - def __init__(self, low, high): + Parameters + ---------- + low, high : :obj:`float` + The endpoints. Use ``np.inf`` for an unbounded side. + closed : {"both", "left", "right", "neither"}, optional + Which endpoints are themselves allowed. Default = "both". + allow_none : :obj:`bool`, optional + Whether ``None`` passes the check, for a parameter whose default is + resolved from the data at fit time. Default = False. + """ + + #: Bracket characters per ``closed`` value, so the error message shows the + #: same interval notation the constraint was declared with. + _BRACKETS = { + "both": ("[", "]"), + "left": ("[", ")"), + "right": ("(", "]"), + "neither": ("(", ")"), + } + + def __init__(self, low, high, closed="both", allow_none=False): + if closed not in self._BRACKETS: + raise ValueError(f"Invalid closed {closed!r}; must be one of {list(self._BRACKETS)}.") self.low = low self.high = high + self.closed = closed + self.allow_none = allow_none def check(self, name, value): """Raise if ``value`` is not a real number inside the interval.""" + if value is None: + if self.allow_none: + return + raise ValueError(f"Invalid {name} None; must be a number.") if not isinstance(value, (int, float, np.integer, np.floating)) or isinstance(value, bool): raise ValueError(f"Invalid {name} {value!r}; must be a number.") - if not self.low <= float(value) <= self.high: - raise ValueError(f"Invalid {name} {value!r}; must lie in [{self.low}, {self.high}].") + value = float(value) + low_ok = self.low <= value if self.closed in ("both", "left") else self.low < value + high_ok = value <= self.high if self.closed in ("both", "right") else value < self.high + if not (low_ok and high_ok): + left, right = self._BRACKETS[self.closed] + raise ValueError( + f"Invalid {name} {value!r}; must lie in {left}{self.low}, {self.high}{right}." + ) #: Constraints shared by the estimators that accept group labels. Assigned to @@ -631,11 +667,12 @@ def _validate_params(self): for name, constraint in self._parameter_constraints.items(): constraint.check(name, getattr(self, name)) - # A class-level mapping from Dataset attributes to fit() arguments. Used by + # A class-level mapping from fit() arguments to Dataset attributes. Used by # fit_dataset() for estimators that take non-standard arguments (e.g., 'z' - # instead of 'y'). Keys are default Dataset attribute names (e.g., 'y') and - # values are the target arg names in the estimator class's fit() method - # (e.g., 'z'). + # instead of 'y'). Keys are the argument names in the estimator class's + # fit() method and values are the Dataset attributes they are filled from, + # so {'z': 'y'} reads "fit()'s z argument takes dataset.y". An argument + # absent from the mapping is filled from the attribute of the same name. _dataset_attr_map = {} @abstractmethod @@ -1646,79 +1683,329 @@ def _reml_profile_nll(self, ratio, y, n, X): ) +#: Location of the Stan program compiled by :obj:`~pymare.estimators.StanMetaRegression`. +#: CmdStanPy needs a real filesystem path -- it hands the file to ``make`` -- so this +#: is a plain join rather than ``importlib.resources``, which would need ``as_file`` +#: to materialize a path PyMARE never needs because it is not zip-imported. +STAN_MODEL_PATH = op.join(op.dirname(__file__), "stan", "meta_regression.stan") + +#: Sampler arguments that PyStan named differently from CmdStanPy. Mapped rather +#: than silently ignored: passing ``num_samples`` to CmdStanPy's ``sample()`` +#: raises a bare ``TypeError`` naming no alternative, and every PyMARE example +#: and test written against the old backend used these names. +PYSTAN_SAMPLING_KWARGS = { + "num_samples": "iter_sampling", + "num_warmup": "iter_warmup", + "num_chains": "chains", + "num_thin": "thin", +} + + +def _import_cmdstanpy(): + """Return the ``cmdstanpy`` module, or raise naming the step that is missing. + + Returns + ------- + module + The imported ``cmdstanpy`` module. + + Raises + ------ + ImportError + If ``cmdstanpy`` is not installed, or if it is installed but no CmdStan + installation can be found. + + Notes + ----- + The two failures are reported separately because their fixes are different + and neither implies the other: ``pip install cmdstanpy`` succeeds without + installing CmdStan itself, which is a C++ toolchain build rather than a + Python package. + """ + try: + import cmdstanpy + except ImportError: + raise ImportError( + "StanMetaRegression requires cmdstanpy, which is an optional dependency. " + "Install it with `pip install pymare[stan]`." + ) + + try: + cmdstanpy.cmdstan_path() + except Exception as exc: + raise ImportError( + "cmdstanpy is installed, but no CmdStan installation was found. Install one " + "with `python -m cmdstanpy.install_cmdstan` (this downloads and builds CmdStan, " + "and needs a C++ toolchain), or point the CMDSTAN environment variable at an " + f"existing installation. cmdstanpy reported: {exc}" + ) + + return cmdstanpy + + +def _build_stan_data(y, v, X, groups=None, tau_prior_scale=None): + """Canonicalize the estimator's inputs into the data block of the Stan program. + + Parameters + ---------- + y : :obj:`numpy.ndarray` of shape (K,) or (K, 1) + Observation-level estimates. + v : :obj:`numpy.ndarray` of shape (K,) or (K, 1) + Observation-level sampling *variances*. + X : :obj:`numpy.ndarray` of shape (K,) or (K, P) + Observation-level predictors, including the intercept. + groups : None or array-like of shape (K,) or (K, 1), optional + One scalar label per observation -- a string, integer or any other + hashable that numpy stores as a single element. Composite labels such + as tuples are not accepted, because numpy reads a sequence of them as a + second dimension. When None (default), every observation is its own + group. + tau_prior_scale : None or :obj:`float`, optional + Scale of the half-normal prior on tau. When None (default), it is set to + ``max(std(y), sqrt(mean(v)))``: the larger of the observed spread of the + estimates and the typical sampling standard deviation. + + Returns + ------- + :obj:`dict` + The data block, with every array in the shape and units the Stan program + declares. + + Raises + ------ + ValueError + If ``y`` is 2-dimensional with more than one column, if ``v`` or ``X`` + disagree with ``y`` about the number of observations, or if any sampling + variance is not positive. + + Notes + ----- + Every shape and unit decision the Stan program depends on is made here and + nowhere else, so ``fit`` carries no downstream conditionals and the + translation is testable without a CmdStan installation. + + ``sigma`` is ``sqrt(v)``, because Stan's ``normal`` takes a standard + deviation and PyMARE stores variances. ``id`` is 1-based consecutive codes from + :func:`~pymare.stats.encode_groups`, because the Stan program declares it + ``int`` -- which is why arbitrary scalar labels are + accepted here. + """ + y = np.asarray(y) + if y.ndim > 1 and y.shape[1] > 1: + raise ValueError( + "The StanMetaRegression estimator currently does " + "not support 2-dimensional inputs. Passed y has " + f"shape {y.shape}." + ) + y = np.asarray(y, dtype=float).reshape(-1) + n_observations = y.shape[0] + if not np.all(np.isfinite(y)): + raise ValueError("Estimates (y) must all be finite.") + + v = np.asarray(v, dtype=float).reshape(-1) + if v.shape[0] != n_observations: + raise ValueError( + f"v must contain one sampling variance per observation: expected " + f"{n_observations}, got {v.shape[0]}." + ) + # Order matters: NaN fails every comparison, so `v <= 0` alone would pass it + # through to sqrt() and on to CmdStan, which rejects it while reading the + # data -- a long way from the input that caused it. + if not np.all(np.isfinite(v)): + raise ValueError("Sampling variances (v) must all be finite.") + if np.any(v <= 0): + raise ValueError("Sampling variances (v) must all be positive.") + + X = np.asarray(X, dtype=float) + if X.ndim == 1: + X = X[:, None] + if X.shape[0] != n_observations: + raise ValueError( + f"X must contain one row per observation: expected {n_observations}, " + f"got {X.shape[0]}." + ) + if not np.all(np.isfinite(X)): + raise ValueError("Predictors (X) must all be finite.") + + codes, labels = encode_groups(groups, n_observations=n_observations) + + if tau_prior_scale is None: + tau_prior_scale = max(np.std(y), np.sqrt(np.mean(v))) + + return { + "N": n_observations, + "C": X.shape[1], + "K": int(labels.size), + "y": y, + "sigma": np.sqrt(v), + "X": X, + "id": (codes + 1).astype(int), + "tau_prior_scale": float(tau_prior_scale), + } + + class StanMetaRegression(BaseEstimator): - """Bayesian meta-regression estimator using Stan. + r"""Bayesian meta-regression estimator using Stan. Parameters ---------- + tau_prior_scale : None or :obj:`float`, optional + Scale of the half-normal prior on tau, the between-group standard + deviation. When None (default), it is set to + ``max(std(y), sqrt(mean(v)))``, the larger of the observed spread of the + estimates and the typical sampling standard deviation. **sampling_kwargs - Optional keyword arguments to pass on to the MCMC sampler - (e.g., `iter` for number of iterations). + Optional keyword arguments to pass on to CmdStanPy's sampler + (e.g., ``iter_sampling`` for the number of post-warmup draws per chain, + ``chains``, ``seed``, ``adapt_delta``). Notes ----- - For most uses, this class should be ignored in favor of the functional - stan() estimator. The object-oriented interface is useful primarily - when fitting the meta-regression model repeatedly to different data; - the separation of .compile() and .fit() steps allows one to compile - the model only once. + The model is - Warning - ------- - :obj:`~pymare.estimators.StanMetaRegression` uses Pystan 3, which requires Python 3.7. - Pystan 3 should not be used with PyMARE and Python 3.6 or earlier. + .. math:: + + y_i &\sim \mathcal{N}(x_i' \beta + \theta_{g(i)}, \sigma_i) \\ + \theta_g &\sim \mathcal{N}(0, \tau) + + where :math:`\sigma_i = \sqrt{v_i}` is the known sampling standard deviation + of observation :math:`i` and :math:`g(i)` is its group. This is the random-effects + meta-analysis model of the Stan User's Guide [1]_ with that guide's stated + extension to observation-level predictors. The reported ``tau2`` is + :math:`\tau^2`, the between-group *variance*, matching what every other + PyMARE estimator reports under that name. + + ``theta`` is given a non-centered parameterization (``theta = tau * + theta_raw`` with ``theta_raw`` standard normal). The centered form produces + the funnel geometry that dominates divergences in hierarchical models with + few groups, which is this estimator's principal use case. + + :math:`\tau` gets a half-normal prior, weakly informative per Stan's + recommendations [2]_ for models with few groups. Its scale is derived from the + data rather than fixed, since a fixed scale would be crushingly informative on + data measured in thousands and vacuous on data measured in thousandths. + + The default scale is ``max(std(y), sqrt(mean(v)))``. :math:`\tau` cannot + plausibly exceed the spread of the estimates, and should not be presumed + smaller than a typical standard error, so the larger of the two never asserts + that :math:`\tau` is small when either quantity says otherwise. Erring large + is deliberate: too small a scale costs coverage, too large costs only + precision in :math:`\tau^2`. ``validation/stan`` records the measurements + behind that choice. Pass ``tau_prior_scale`` to override it. + + :math:`\beta` keeps Stan's implicit improper uniform prior, so under a + diffuse prior on :math:`\tau` the posterior means agree with + :obj:`~pymare.estimators.VarianceBasedLikelihoodEstimator` at ``method="ML"``. + + The Stan program is compiled on first use and cached beside the installed + source, so the cost is paid once per installation rather than per fit. + + References + ---------- + .. [1] Stan Development Team. Stan User's Guide, "Measurement Error and + Meta-Analysis", section "Meta-Analysis". + https://mc-stan.org/docs/stan-users-guide/measurement-error.html + .. [2] Stan Development Team. Prior Choice Recommendations. + https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations + + .. versionchanged:: 0.0.11 + + - The backend moved from PyStan 3 to CmdStanPy. PyStan's sampler argument + names (``num_samples``, ``num_warmup``, ``num_chains``, ``num_thin``) + are rejected with a message naming their replacements. + - ``tau2`` is now the between-group variance rather than its square root, + and sampling variances are converted to standard deviations before + being passed to Stan. Both were wrong before, so posterior estimates + change. + - ``groups`` accepts scalar labels of any type, not only integers in + ``1..k``, and :meth:`fit_dataset` now passes ``dataset.g`` rather than + dropping it. + - ``ci`` now sets the width of the reported credible interval. It was + previously accepted and ignored. """ - _result_cls = BayesianMetaRegressionResults + _dataset_attr_map = {"groups": "g"} + + _parameter_constraints = { + "tau_prior_scale": Interval(0.0, np.inf, closed="neither", allow_none=True), + } + + def __init__(self, tau_prior_scale=None, **sampling_kwargs): + renamed = {k: v for k, v in PYSTAN_SAMPLING_KWARGS.items() if k in sampling_kwargs} + if renamed: + raise TypeError( + "These are PyStan argument names, which StanMetaRegression no longer accepts: " + + ", ".join( + f"{old!r} (CmdStanPy calls it {new!r})" for old, new in renamed.items() + ) + + "." + ) - def __init__(self, **sampling_kwargs): + self.tau_prior_scale = tau_prior_scale self.sampling_kwargs = sampling_kwargs self.model = None self.result_ = None + self._validate_params() - if sys.version_info < (3, 7): - raise RuntimeError( - "StanMetaRegression uses Pystan 3, which requires python 3.7 or higher. " - f"You are running Python {sys.version_info.major}.{sys.version_info.minor}. " - "Pystan 3 should not be used with PyMARE and Python 3.6 or earlier." - ) + def compile(self, force=False): + """Compile the Stan model. - def compile(self): - """Compile the Stan model.""" - # Note: we deliberately use a centered parameterization for the - # thetas at the moment. This is sub-optimal in terms of estimation, - # but allows us to avoid having to add extra logic to detect and - # handle intercepts in X. - spec = """ - data { - int N; - int K; - vector[N] y; - array[N] int id; - int C; - matrix[K, C] X; - vector[N] sigma; - } - parameters { - vector[C] beta; - vector[K] theta; - real tau2; - } - transformed parameters { - vector[N] mu; - mu = theta[id] + X * beta; - } - model { - y ~ normal(mu, sigma); - theta ~ normal(0, tau2); - } + Parameters + ---------- + force : :obj:`bool`, optional + Whether to recompile even when an up-to-date executable already + exists. Default = False. + + Returns + ------- + :obj:`~pymare.estimators.StanMetaRegression` + The instance, so that ``compile()`` can be chained. + + Notes + ----- + Called by :meth:`fit` when needed, so it never has to be called + directly. Calling it in advance is worthwhile when the same estimator + will be fitted to several datasets, because the compiled executable does + not depend on the data. + + The executable is written beside the installed ``.stan`` file, where + CmdStanPy finds and reuses it on later runs. If that directory is not + writable -- a read-only ``site-packages``, for instance -- it falls back + to ``~/.pymare/stan`` and warns once. """ + cmdstanpy = _import_cmdstanpy() + + try: + self.model = cmdstanpy.CmdStanModel(stan_file=STAN_MODEL_PATH, force_compile=force) + return self + except Exception as unwritable: + # Deliberately broad: CmdStanPy reports any failed make invocation + # as ValueError, so catching OSError would never fire here. + first_failure = unwritable + + # Compile a copy. exe_file= would not help: it names an executable to + # reuse rather than a destination to build into, and make writes its + # intermediates beside the source either way. copy2 preserves the mtime, + # so the copy never looks newer than its own executable and CmdStanPy's + # timestamp check keeps the cached build across processes. + fallback_dir = op.join(op.expanduser("~"), ".pymare", "stan") try: - import stan - except ImportError: - raise ImportError("Please install pystan.") + os.makedirs(fallback_dir, exist_ok=True) + fallback_source = op.join(fallback_dir, op.basename(STAN_MODEL_PATH)) + shutil.copy2(STAN_MODEL_PATH, fallback_source) + model = cmdstanpy.CmdStanModel(stan_file=fallback_source, force_compile=force) + except Exception: + # Somewhere writable failed too, so the first failure was not about + # writing. Report that one -- it names the real problem. + raise first_failure - self.model = stan.build(spec, data=self.data) + warn( + f"Could not compile the Stan model beside {STAN_MODEL_PATH}, most likely because " + f"that directory is not writable. Compiled into {fallback_dir} instead.", + stacklevel=2, + ) + self.model = model + return self def fit(self, y, v, X, groups=None): """Run the Stan sampler and return results. @@ -1733,17 +2020,28 @@ def fit(self, y, v, X, groups=None): 1d or 2d array containing observation-level predictors (including intercept); has dimensions K x P, where K is the number of observations and P is the number of predictor variables. - groups : :obj:`list` of :obj:`int`, optional - 1d array of integers identifying - groups of observations in the y/v/X inputs. If - provided, values must consist of integers in the range of 1..k - (inclusive), where k is the number of distinct groups. When - None (default), it is assumed that each observation in the - inputs is a separate group. + groups : None or array-like of shape (K,), optional + One scalar label per observation, identifying the groups of + observations in the y/v/X inputs. Labels may be strings, integers + or any other hashable that numpy stores as a single element, and + need not be consecutive; they are encoded internally in order of + first occurrence by :func:`~pymare.stats.encode_groups`. Composite + labels such as tuples are not accepted, because numpy reads a + sequence of them as a 2-dimensional array. When None (default), + each observation in the inputs is treated as a separate group. Returns ------- - A StanFit4Model object (see PyStan documentation for details). + :obj:`~pymare.estimators.StanMetaRegression` + The fitted instance. + + Warns + ----- + UserWarning + If the sampler reported divergent transitions. Divergences mean the + sampler could not explore part of the posterior, so the reported + means and intervals may be biased; refitting with a larger + ``adapt_delta`` is the usual remedy. Notes ----- @@ -1751,41 +2049,53 @@ def fit(self, y, v, X, groups=None): observations belong to at least one common sampling unit, the `groups` argument can specify the nesting structure (i.e., which rows in `y`, `v`, and `X` belong to each group). + + The raw CmdStanPy fit is kept on ``self.result_``, so its diagnostics + remain reachable -- ``est.result_.diagnose()`` reports R-hat, effective + sample size, E-BFMI and treedepth alongside divergences. + + .. versionchanged:: 0.0.5 + ``groups`` accepts arbitrary hashable labels, and passing a numpy + array no longer raises. """ # This resets the Estimator's dataset_ attribute. fit_dataset will overwrite if called. self.dataset_ = None - if y.ndim > 1 and y.shape[1] > 1: - raise ValueError( - "The StanMetaRegression estimator currently does " - "not support 2-dimensional inputs. Passed y has " - "shape {}.".format(y.shape) - ) - - N = y.shape[0] - groups = groups or np.arange(1, N + 1, dtype=int) - K = encode_groups(np.asarray(groups).ravel())[1].size - - data = { - "K": K, - "N": N, - "id": groups, - "C": X.shape[1], - "X": X, - "y": y.ravel(), - "sigma": v.ravel(), - } - - self.data = data + self.data = _build_stan_data(y, v, X, groups=groups, tau_prior_scale=self.tau_prior_scale) if self.model is None: self.compile() - self.result_ = self.model.sample(**self.sampling_kwargs) + self.result_ = self.model.sample(data=self.data, **self.sampling_kwargs) + + # CmdStanPy logs its own diagnostic warnings. Reraising this one through + # the warnings module puts it under the caller's warning filters and + # makes it assertable, which a log record is not. + divergences = int(np.sum(self.result_.divergences)) + if divergences: + warn( + f"The sampler reported {divergences} divergent transition(s). The posterior " + "summaries may be biased. Refit with a larger adapt_delta (e.g. " + "StanMetaRegression(adapt_delta=0.99)), and see result_.diagnose() for the " + "full diagnostic report.", + stacklevel=2, + ) + return self def summary(self, ci=95): - """Generate a BayesianMetaRegressionResults object from the fitted estimator.""" + """Generate a BayesianMetaRegressionResults object from the fitted estimator. + + Parameters + ---------- + ci : :obj:`float`, optional + Desired width of the credible interval, as a percentage. + Default = 95.0 (95%). + + Returns + ------- + :obj:`~pymare.results.BayesianMetaRegressionResults` + """ if self.result_ is None: name = self.__class__.__name__ raise ValueError( diff --git a/pymare/estimators/stan/meta_regression.stan b/pymare/estimators/stan/meta_regression.stan new file mode 100644 index 0000000..f37b51d --- /dev/null +++ b/pymare/estimators/stan/meta_regression.stan @@ -0,0 +1,45 @@ +// Bayesian hierarchical meta-regression. +// +// y_i ~ normal(x_i' beta + theta_{g(i)}, sigma_i) i = 1..N +// theta_g ~ normal(0, tau) g = 1..K +// +// The Stan User's Guide random-effects meta-analysis model (Measurement Error +// and Meta-Analysis) with that guide's extension to observation-level +// predictors. +// +// sigma_i is the *known* sampling standard deviation sqrt(v_i): Stan's normal() +// takes a scale, not a variance. tau is the between-group standard deviation; +// tau2 = tau^2 is the variance every other PyMARE estimator reports. +data { + int N; // observations + int C; // predictors (columns of X) + int K; // groups + vector[N] y; // observed effect sizes + vector[N] sigma; // sampling standard deviations + matrix[N, C] X; // one row per observation + array[N] int id; // 1-based group index per observation + real tau_prior_scale; // scale of the half-normal prior on tau +} +parameters { + vector[C] beta; + vector[K] theta_raw; + real tau; +} +transformed parameters { + // Non-centered. Sampling theta_raw ~ N(0, 1) and scaling avoids the funnel + // geometry theta ~ normal(0, tau) produces as tau approaches zero, which is + // the dominant source of divergences when groups are few. + vector[K] theta = tau * theta_raw; +} +model { + theta_raw ~ std_normal(); + // Half-normal: the declaration on tau truncates the normal at zero. + tau ~ normal(0, tau_prior_scale); + // Equivalent to y ~ normal(X * beta + theta[id], sigma), but documented as the + // faster form. The vector-alpha/vector-sigma overload takes a per-observation + // intercept (the group effect) and scale (the sampling SD). + y ~ normal_id_glm(X, theta[id], beta, sigma); +} +generated quantities { + real tau2 = square(tau); +} diff --git a/pymare/results.py b/pymare/results.py index b6c4e89..632b02b 100644 --- a/pymare/results.py +++ b/pymare/results.py @@ -4,7 +4,7 @@ import itertools import math from functools import lru_cache -from inspect import getfullargspec +from inspect import getfullargspec, signature from warnings import warn import numpy as np @@ -929,17 +929,77 @@ def to_df(self, **kwargs): return df +def _arviz_credible_interval_kwargs(ci): + """Return the ArviZ ``summary()`` arguments that request a `ci`% HDI. + + Parameters + ---------- + ci : :obj:`float` + Desired width of the credible interval, as a percentage. + + Returns + ------- + :obj:`dict` + Keyword arguments for the installed ArviZ version's ``summary()``. + + Notes + ----- + ArviZ 1.0 renamed ``hdi_prob`` to ``ci_prob`` and paired it with a + ``ci_kind`` defaulting to an equal-tailed rather than a highest-density + interval, so ``ci_kind="hdi"`` is needed to keep ``ci`` meaning the same + thing across both versions. + + ``round_to="none"`` is not cosmetic: ArviZ 1.x otherwise formats the summary + for display by converting the floats to strings, so arithmetic on the + returned DataFrame would silently concatenate. ArviZ 0.x has no such + argument. + """ + if int(az.__version__.split(".")[0]) >= 1: + return {"ci_prob": ci / 100.0, "ci_kind": "hdi", "round_to": "none"} + return {"hdi_prob": ci / 100.0} + + +def _accepts_var_names(plotter): + """Report whether an ArviZ plotting function takes a ``var_names`` argument. + + Parameters + ---------- + plotter : callable + A function from the ArviZ namespace. + + Returns + ------- + :obj:`bool` + True when ``var_names`` can be passed to it. + """ + try: + return "var_names" in signature(plotter).parameters + except (TypeError, ValueError): + return False + + class BayesianMetaRegressionResults: """Container for MCMC sampling-based PyMARE meta-regression estimators. Parameters ---------- - data : :obj:`pystan.StanFit4Model` or :obj:`arviz.InferenceData` - Either a StanFit4Model instance returned from PyStan or an ArviZ InferenceData instance. + data : :obj:`cmdstanpy.CmdStanMCMC` or :obj:`arviz.InferenceData` + Either a CmdStanMCMC instance returned from CmdStanPy or an object ArviZ + already understands (an InferenceData under ArviZ 0.x, a DataTree under 1.x). dataset : :obj:`~pymare.core.Dataset` A Dataset instance containing the inputs to the estimator. ci : :obj:`float`, optional - Desired width of highest posterior density (HPD) interval. Default = 95.0 (95%). + Desired width of the credible interval, as a percentage, used as the + default for :meth:`summary`. Default = 95.0 (95%). + + Notes + ----- + .. versionchanged:: 0.0.5 + + - Accepts a :obj:`cmdstanpy.CmdStanMCMC` rather than a PyStan fit. + - ``ci`` now sets the width of the interval :meth:`summary` reports. It + was previously stored and never used, so the interval was whatever + ArviZ defaulted to. """ def __init__(self, data, dataset, ci=95.0): @@ -948,8 +1008,20 @@ def __init__(self, data, dataset, ci=95.0): "ArviZ package must be installed in order to work " "with the BayesianMetaRegressionResults class." ) - if data.__class__.__name__ == "StanFit4Model": - data = az.from_pystan(data) + if not 0 < ci < 100: + raise ValueError(f"Invalid ci {ci!r}; must lie in (0, 100).") + + # Convert here, not on use: ArviZ 1.x removed the automatic dispatch that + # let summary() accept a sampler fit, so a raw fit would fail at every + # call site instead of this one. The import is lazy because a caller + # passing an already-converted object should not need cmdstanpy. + try: + from cmdstanpy import CmdStanMCMC + except ImportError: + CmdStanMCMC = () + if isinstance(data, CmdStanMCMC): + data = az.from_cmdstanpy(data) + self.data = data self.dataset = dataset self.ci = ci @@ -963,21 +1035,29 @@ def summary(self, include_theta=False, **kwargs): Whether or not to include the estimated group-level means in the summary. Default = False. **kwargs - Optional keyword arguments to pass onto ArviZ's summary(). + Optional keyword arguments to pass onto ArviZ's summary(). Anything + passed here wins over the defaults derived from ``ci``. Returns ------- :obj:`pandas.DataFrame` A pandas DataFrame, unless the `fmt="xarray"` argument is passed in kwargs, in which case an xarray Dataset is returned. + + Notes + ----- + The columns follow whichever ArviZ version is installed; 1.x renamed + several of them. Index by label rather than by position. """ var_names = ["beta", "tau2"] if include_theta: var_names.append("theta") var_names = kwargs.pop("var_names", var_names) + for key, value in _arviz_credible_interval_kwargs(self.ci).items(): + kwargs.setdefault(key, value) return az.summary(self.data, var_names, **kwargs) - def plot(self, kind="trace", **kwargs): + def plot(self, kind="trace", include_theta=False, **kwargs): """Generate various plots of the posterior estimates via ArviZ. Parameters @@ -986,16 +1066,41 @@ def plot(self, kind="trace", **kwargs): The type of ArviZ plot to generate. Can be any named function of the form "plot_{}" in the ArviZ namespace (e.g., 'trace', 'forest', 'posterior', etc.). Default = 'trace'. + include_theta : :obj:`bool`, optional + Whether or not to include the estimated group-level means in the plot. + Default = False. **kwargs Optional keyword arguments passed onto the corresponding - ArviZ plotting function (see ArviZ docs for details). + ArviZ plotting function (see ArviZ docs for details). Passing + ``var_names`` overrides the selection ``include_theta`` implies. Returns ------- A matplotlib or bokeh object, depending on plot kind and kwargs. + + Raises + ------ + ValueError + If ArviZ has no plotting function of the requested kind. + + Notes + ----- + The plotted variables default to the ones :meth:`summary` reports rather + than everything recorded. A fitted model has one ``theta`` per group, so + plotting everything gives one panel per group -- illegible, and a hard + error under ArviZ 1.x, which caps panels at + ``rcParams["plot.max_subplots"]``. """ name = "plot_{}".format(kind) - plotter = getattr(az, name) + # Three-argument getattr: the two-argument form raises AttributeError + # before the check below can turn it into the documented ValueError. + plotter = getattr(az, name, None) if plotter is None: raise ValueError("ArviZ has no plotting function '{}'.".format(name)) - plotter(self.data, **kwargs) + + # Not every ArviZ plot takes var_names -- plot_energy, for one, takes + # only the data -- so only offer the default to those that do. + if "var_names" not in kwargs and _accepts_var_names(plotter): + kwargs["var_names"] = ["beta", "tau2", "theta"] if include_theta else ["beta", "tau2"] + + return plotter(self.data, **kwargs) diff --git a/pymare/tests/conftest.py b/pymare/tests/conftest.py index 0169742..56c1c15 100644 --- a/pymare/tests/conftest.py +++ b/pymare/tests/conftest.py @@ -1,5 +1,6 @@ """Fixtures for the PyMARE test suite.""" +import os import os.path as op import numpy as np @@ -16,7 +17,34 @@ VarianceBasedLikelihoodEstimator, WeightedLeastSquares, ) -from pymare.tests.utils import get_test_data_path +from pymare.tests.utils import cmdstan_is_available, get_test_data_path + + +def pytest_collection_modifyitems(config, items): + """Fail the run where CmdStan is declared present but is not. + + The Stan tests skip when CmdStan is missing, so a contributor without it does + not see red. In CI that leniency is wrong: a skip is indistinguishable from a + pass in a job log, which is how this job once reported success while running + none of the tests it existed to run. + + ``PYMARE_REQUIRE_CMDSTAN=1``, which the Stan CI job sets, asserts that the + environment should be able to run them. Failing at collection reports that + once rather than as a quietly shorter run. + """ + if os.environ.get("PYMARE_REQUIRE_CMDSTAN") != "1": + return + if cmdstan_is_available(): + return + + raise pytest.UsageError( + "PYMARE_REQUIRE_CMDSTAN=1 says this environment should be able to run the " + "Stan tests, but cmdstanpy or its CmdStan installation is missing, so they " + "would all skip. Install them with `pip install -e .[stan]` followed by " + "`python -m cmdstanpy.install_cmdstan`, or unset PYMARE_REQUIRE_CMDSTAN to " + "allow the skip." + ) + # ----------------------------------------------------------------------------- # Basic data @@ -353,3 +381,74 @@ def two_samp_data(): "sd2": np.sqrt(np.array([4, 16])), "n2": np.array([12, 16]), } + + +# ----------------------------------------------------------------------------- +# Stan estimator +# ----------------------------------------------------------------------------- + + +@pytest.fixture +def fake_home(tmp_path, monkeypatch): + """Point ``os.path.expanduser("~")`` at a temporary directory, on any platform. + + Returns + ------- + :obj:`pathlib.Path` + The directory ``~`` now expands to. + + Notes + ----- + Both variables are needed. POSIX ``expanduser`` reads ``HOME``; Windows reads + ``USERPROFILE`` and ignores ``HOME`` entirely. Setting only ``HOME`` + therefore looks portable while silently leaving the real profile in place on + Windows, so a test can pass on two platforms and write into the runner's + actual home directory on the third. + """ + monkeypatch.setenv("HOME", str(tmp_path)) + monkeypatch.setenv("USERPROFILE", str(tmp_path)) + return tmp_path + + +@pytest.fixture(scope="package") +def planted_hierarchical_dataset(): + """Simulate a Dataset from the model ``meta_regression.stan`` encodes. + + Returns + ------- + :obj:`tuple` of (:obj:`~pymare.core.Dataset`, :obj:`dict`) + The simulated Dataset, with one group label per observation in ``g``, + and the parameter values it was generated from. + + Notes + ----- + The sampling standard deviations are drawn from ``uniform(0.1, 0.4)``, well + away from 1, and that is load-bearing. The ``variables`` fixture has ``v`` + near 1, where ``v`` and ``sqrt(v)`` differ by a few percent, so a model that + confuses them fits it about as well as the correct one. Here ``v`` spans 0.01 + to 0.16 against ``sqrt(v)`` from 0.1 to 0.4 -- a factor of 2.5 to 10 in one + direction -- so the mistake shows up as an inflated tau2. Do not substitute + ``variables`` here. + + ``tau2`` and ``tau`` are likewise kept apart (0.25 against 0.5) so reporting + one under the other's name fails a tight interval. + """ + rng = np.random.default_rng(20250818) + + n_groups, per_group = 30, 3 + beta = np.array([0.5, -0.8]) + tau = 0.5 + + groups = np.repeat(np.arange(n_groups), per_group) + n_observations = groups.size + + moderator = rng.normal(size=n_observations) + # Dataset prepends the intercept itself, so beta[0] is the intercept and + # beta[1] the moderator slope in the X the estimator will actually see. + X = np.column_stack([np.ones(n_observations), moderator]) + theta = rng.normal(0, tau, size=n_groups) + sigma = rng.uniform(0.1, 0.4, size=n_observations) + y = X @ beta + theta[groups] + rng.normal(0, sigma) + + dataset = Dataset(y=y, v=sigma**2, X=moderator, X_names=["moderator"], g=groups) + return dataset, {"beta": beta, "tau": tau, "tau2": tau**2} diff --git a/pymare/tests/data/stan_validation.json b/pymare/tests/data/stan_validation.json new file mode 100644 index 0000000..50ae722 --- /dev/null +++ b/pymare/tests/data/stan_validation.json @@ -0,0 +1,431 @@ +{ + "replications": 100, + "seed": 20260818, + "elapsed_seconds": 335.2, + "thresholds": { + "min_coverage": 0.9, + "max_beta_bias": 0.1 + }, + "cells": [ + { + "name": "groups=5", + "config": { + "n_groups": "5", + "group_size": "3", + "tau2": "0.1", + "sigma_scale": "1.0", + "n_predictors": "2", + "unbalanced": "False" + }, + "replications": 100, + "true_beta": [ + 0.5, + -0.8 + ], + "beta_bias_per_coefficient": [ + 0.009240000000000003, + -0.008429999999999955 + ], + "beta_coverage_per_coefficient": [ + 0.93, + 0.93 + ], + "beta_bias": 0.009240000000000003, + "beta_coverage": 0.93, + "coverage_se": 0.02551470164434614, + "tau2_bias": 0.12235999999999998, + "tau2_truth": 0.09999999999999998, + "fits_with_divergences": 54 + }, + { + "name": "groups=20", + "config": { + "n_groups": "20", + "group_size": "3", + "tau2": "0.1", + "sigma_scale": "1.0", + "n_predictors": "2", + "unbalanced": "False" + }, + "replications": 100, + "true_beta": [ + 0.5, + -0.8 + ], + "beta_bias_per_coefficient": [ + -0.005699999999999999, + -0.005239999999999958 + ], + "beta_coverage_per_coefficient": [ + 0.91, + 0.9 + ], + "beta_bias": -0.005699999999999999, + "beta_coverage": 0.9, + "coverage_se": 0.03, + "tau2_bias": 0.010289999999999995, + "tau2_truth": 0.09999999999999998, + "fits_with_divergences": 0 + }, + { + "name": "groups=50", + "config": { + "n_groups": "50", + "group_size": "3", + "tau2": "0.1", + "sigma_scale": "1.0", + "n_predictors": "2", + "unbalanced": "False" + }, + "replications": 100, + "true_beta": [ + 0.5, + -0.8 + ], + "beta_bias_per_coefficient": [ + -0.004449999999999992, + -0.004719999999999964 + ], + "beta_coverage_per_coefficient": [ + 0.99, + 0.96 + ], + "beta_bias": -0.004719999999999964, + "beta_coverage": 0.96, + "coverage_se": 0.019595917942265433, + "tau2_bias": 0.008279999999999996, + "tau2_truth": 0.09999999999999998, + "fits_with_divergences": 0 + }, + { + "name": "tau2=0", + "config": { + "n_groups": "20", + "group_size": "3", + "tau2": "0.0", + "sigma_scale": "1.0", + "n_predictors": "2", + "unbalanced": "False" + }, + "replications": 100, + "true_beta": [ + 0.5, + -0.8 + ], + "beta_bias_per_coefficient": [ + 0.004770000000000009, + -0.0024699999999999622 + ], + "beta_coverage_per_coefficient": [ + 0.97, + 0.96 + ], + "beta_bias": 0.004770000000000009, + "beta_coverage": 0.96, + "coverage_se": 0.019595917942265433, + "tau2_bias": 0.00402, + "tau2_truth": 0.0, + "fits_with_divergences": 14 + }, + { + "name": "tau2=0.1", + "config": { + "n_groups": "20", + "group_size": "3", + "tau2": "0.1", + "sigma_scale": "1.0", + "n_predictors": "2", + "unbalanced": "False" + }, + "replications": 100, + "true_beta": [ + 0.5, + -0.8 + ], + "beta_bias_per_coefficient": [ + -0.0036099999999999995, + -0.0014999999999999525 + ], + "beta_coverage_per_coefficient": [ + 0.94, + 0.93 + ], + "beta_bias": -0.0036099999999999995, + "beta_coverage": 0.93, + "coverage_se": 0.02551470164434614, + "tau2_bias": 0.017839999999999998, + "tau2_truth": 0.09999999999999998, + "fits_with_divergences": 0 + }, + { + "name": "tau2=1", + "config": { + "n_groups": "20", + "group_size": "3", + "tau2": "1.0", + "sigma_scale": "1.0", + "n_predictors": "2", + "unbalanced": "False" + }, + "replications": 100, + "true_beta": [ + 0.5, + -0.8 + ], + "beta_bias_per_coefficient": [ + 0.010260000000000005, + 0.006100000000000039 + ], + "beta_coverage_per_coefficient": [ + 0.96, + 0.97 + ], + "beta_bias": 0.010260000000000005, + "beta_coverage": 0.96, + "coverage_se": 0.019595917942265433, + "tau2_bias": 0.18483, + "tau2_truth": 1.0, + "fits_with_divergences": 0 + }, + { + "name": "singletons", + "config": { + "n_groups": "20", + "group_size": "1", + "tau2": "0.1", + "sigma_scale": "1.0", + "n_predictors": "2", + "unbalanced": "False" + }, + "replications": 100, + "true_beta": [ + 0.5, + -0.8 + ], + "beta_bias_per_coefficient": [ + -0.013600000000000001, + 0.0027800000000000524 + ], + "beta_coverage_per_coefficient": [ + 0.98, + 0.95 + ], + "beta_bias": -0.013600000000000001, + "beta_coverage": 0.95, + "coverage_se": 0.021794494717703377, + "tau2_bias": 0.020919999999999998, + "tau2_truth": 0.09999999999999998, + "fits_with_divergences": 1 + }, + { + "name": "unequal groups", + "config": { + "n_groups": "20", + "group_size": "unequal", + "tau2": "0.1", + "sigma_scale": "1.0", + "n_predictors": "2", + "unbalanced": "False" + }, + "replications": 100, + "true_beta": [ + 0.5, + -0.8 + ], + "beta_bias_per_coefficient": [ + -0.0012999999999999956, + -0.004179999999999957 + ], + "beta_coverage_per_coefficient": [ + 0.94, + 0.99 + ], + "beta_bias": -0.004179999999999957, + "beta_coverage": 0.94, + "coverage_se": 0.023748684174075843, + "tau2_bias": 0.02412, + "tau2_truth": 0.09999999999999998, + "fits_with_divergences": 0 + }, + { + "name": "sigma x0.1", + "config": { + "n_groups": "20", + "group_size": "3", + "tau2": "0.1", + "sigma_scale": "0.1", + "n_predictors": "2", + "unbalanced": "False" + }, + "replications": 100, + "true_beta": [ + 0.5, + -0.8 + ], + "beta_bias_per_coefficient": [ + 0.0026900000000000014, + -0.00048000000000000045 + ], + "beta_coverage_per_coefficient": [ + 0.94, + 0.95 + ], + "beta_bias": 0.0026900000000000014, + "beta_coverage": 0.94, + "coverage_se": 0.023748684174075843, + "tau2_bias": 0.020710000000000003, + "tau2_truth": 0.09999999999999998, + "fits_with_divergences": 0 + }, + { + "name": "sigma x10", + "config": { + "n_groups": "20", + "group_size": "3", + "tau2": "0.1", + "sigma_scale": "10.0", + "n_predictors": "2", + "unbalanced": "False" + }, + "replications": 100, + "true_beta": [ + 0.5, + -0.8 + ], + "beta_bias_per_coefficient": [ + -0.04500000000000001, + 0.017270000000000042 + ], + "beta_coverage_per_coefficient": [ + 0.98, + 0.94 + ], + "beta_bias": -0.04500000000000001, + "beta_coverage": 0.94, + "coverage_se": 0.023748684174075843, + "tau2_bias": 0.42099, + "tau2_truth": 0.09999999999999998, + "fits_with_divergences": 6 + }, + { + "name": "1 predictor", + "config": { + "n_groups": "20", + "group_size": "3", + "tau2": "0.1", + "sigma_scale": "1.0", + "n_predictors": "1", + "unbalanced": "False" + }, + "replications": 100, + "true_beta": [ + 0.5 + ], + "beta_bias_per_coefficient": [ + -0.0045099999999999975 + ], + "beta_coverage_per_coefficient": [ + 0.94 + ], + "beta_bias": -0.0045099999999999975, + "beta_coverage": 0.94, + "coverage_se": 0.023748684174075843, + "tau2_bias": 0.008289999999999995, + "tau2_truth": 0.09999999999999998, + "fits_with_divergences": 0 + }, + { + "name": "3 predictors", + "config": { + "n_groups": "20", + "group_size": "3", + "tau2": "0.1", + "sigma_scale": "1.0", + "n_predictors": "3", + "unbalanced": "False" + }, + "replications": 100, + "true_beta": [ + 0.5, + -0.8, + 0.3 + ], + "beta_bias_per_coefficient": [ + -0.009879999999999996, + 0.0048900000000000375, + -0.0016899999999999875 + ], + "beta_coverage_per_coefficient": [ + 0.94, + 0.92, + 0.99 + ], + "beta_bias": -0.009879999999999996, + "beta_coverage": 0.92, + "coverage_se": 0.027129319932501065, + "tau2_bias": 0.020689999999999997, + "tau2_truth": 0.09999999999999998, + "fits_with_divergences": 0 + }, + { + "name": "unbalanced covariate", + "config": { + "n_groups": "20", + "group_size": "3", + "tau2": "0.1", + "sigma_scale": "1.0", + "n_predictors": "2", + "unbalanced": "True" + }, + "replications": 100, + "true_beta": [ + 0.5, + -0.8 + ], + "beta_bias_per_coefficient": [ + -0.00859, + 0.007870000000000042 + ], + "beta_coverage_per_coefficient": [ + 0.97, + 0.95 + ], + "beta_bias": -0.00859, + "beta_coverage": 0.95, + "coverage_se": 0.021794494717703377, + "tau2_bias": 0.020329999999999994, + "tau2_truth": 0.09999999999999998, + "fits_with_divergences": 0 + }, + { + "name": "unbalanced covariate, tau2=1", + "config": { + "n_groups": "20", + "group_size": "3", + "tau2": "1.0", + "sigma_scale": "1.0", + "n_predictors": "2", + "unbalanced": "True" + }, + "replications": 100, + "true_beta": [ + 0.5, + -0.8 + ], + "beta_bias_per_coefficient": [ + -0.018280000000000005, + -0.024199999999999958 + ], + "beta_coverage_per_coefficient": [ + 0.98, + 0.93 + ], + "beta_bias": -0.024199999999999958, + "beta_coverage": 0.93, + "coverage_se": 0.02551470164434614, + "tau2_bias": 0.04930000000000002, + "tau2_truth": 1.0, + "fits_with_divergences": 0 + } + ] +} \ No newline at end of file diff --git a/pymare/tests/test_estimators.py b/pymare/tests/test_estimators.py index 0c6b461..3129050 100644 --- a/pymare/tests/test_estimators.py +++ b/pymare/tests/test_estimators.py @@ -13,7 +13,7 @@ VarianceBasedLikelihoodEstimator, WeightedLeastSquares, ) -from pymare.estimators.estimators import _collapse_n_inputs +from pymare.estimators.estimators import Interval, Options, _collapse_n_inputs from pymare.stats import ( DEFAULT_RHO, collapse_groups, @@ -894,3 +894,59 @@ def test_near_equal_sample_sizes_warn_rather_than_abort(): with pytest.warns(UserWarning, match="too close"): SampleSizeBasedLikelihoodEstimator().fit(y=y, n=n, X=np.ones((20, 1))) + + +# ----------------------------------------------------------------------------- +# Parameter constraints +# ----------------------------------------------------------------------------- + + +def test_interval_rejects_an_unknown_closed_setting(): + """A typo in a constraint declaration must fail where it is written.""" + with pytest.raises(ValueError, match="must be one of"): + Interval(0.0, 1.0, closed="closed") + + +@pytest.mark.parametrize( + ("closed", "accepted", "rejected"), + [ + ("both", (0.0, 0.5, 1.0), ()), + ("left", (0.0, 0.5), (1.0,)), + ("right", (0.5, 1.0), (0.0,)), + ("neither", (0.5,), (0.0, 1.0)), + ], +) +def test_interval_honours_which_endpoints_are_closed(closed, accepted, rejected): + """Each ``closed`` setting admits exactly its own endpoints.""" + constraint = Interval(0.0, 1.0, closed=closed) + + for value in accepted: + constraint.check("rho", value) + + for value in rejected: + with pytest.raises(ValueError, match="must lie in"): + constraint.check("rho", value) + + +def test_interval_allows_none_only_when_told_to(): + """None means "derive it from the data" for some parameters, and nothing for others.""" + Interval(0.0, np.inf, closed="neither", allow_none=True).check("tau_prior_scale", None) + + with pytest.raises(ValueError, match="must be a number"): + Interval(0.0, 1.0).check("rho", None) + + +@pytest.mark.parametrize("value", ["0.5", True, None, [0.5]]) +def test_interval_rejects_things_that_are_not_numbers(value): + """Booleans are ints in Python, so they need excluding explicitly.""" + with pytest.raises(ValueError, match="must be a number"): + Interval(0.0, 1.0).check("rho", value) + + +def test_options_rejects_a_value_outside_the_set(): + """The other constraint type, checked here because nothing else covered it.""" + constraint = Options(("individual", "rescale")) + constraint.check("weight_scheme", "rescale") + + with pytest.raises(ValueError, match="must be one of"): + constraint.check("weight_scheme", "collapsed") diff --git a/pymare/tests/test_stan_estimators.py b/pymare/tests/test_stan_estimators.py index 796d05c..c063f01 100644 --- a/pymare/tests/test_stan_estimators.py +++ b/pymare/tests/test_stan_estimators.py @@ -1,41 +1,787 @@ """Tests for estimators that use stan. -pystan is an optional dependency, so these tests skip rather than fail when it -is missing. Marked ``stan`` so CI can run them in a job that installs it. +cmdstanpy and CmdStan are optional, so the tests that sample are marked ``stan`` +and skip when either is missing. The tests that only exercise the translation +from PyMARE's inputs to Stan's data block are deliberately left unmarked: they +need neither, so they run in the ordinary unit job on every platform, which is +where the defects this file now pins would have been caught years earlier. """ -from importlib.util import find_spec +import ntpath +import os.path as op +import posixpath +import sys +import warnings +import numpy as np import pytest -from pymare.estimators import StanMetaRegression +from pymare import meta_regression +from pymare.estimators import StanMetaRegression, VarianceBasedLikelihoodEstimator +from pymare.estimators.estimators import ( + STAN_MODEL_PATH, + _build_stan_data, + _import_cmdstanpy, +) +from pymare.results import ( + BayesianMetaRegressionResults, + _accepts_var_names, + _arviz_credible_interval_kwargs, +) +from pymare.tests import conftest +from pymare.tests.utils import ( + STAN_VALIDATION_CELLS, + STAN_VALIDATION_THRESHOLDS, + cmdstan_is_available, + load_stan_validation, +) + +requires_cmdstan = pytest.mark.skipif( + not cmdstan_is_available(), + reason="requires cmdstanpy and a CmdStan installation", +) + + +# ----------------------------------------------------------------------------- +# Detecting a missing dependency. No CmdStan needed -- the point is its absence. +# ----------------------------------------------------------------------------- + + +def test_import_cmdstanpy_reports_a_missing_package(monkeypatch): + """The message must name the install command, not just the module.""" + monkeypatch.setitem(sys.modules, "cmdstanpy", None) + + with pytest.raises(ImportError, match=r"pip install pymare\[stan\]"): + _import_cmdstanpy() + + +def test_import_cmdstanpy_reports_a_missing_cmdstan(monkeypatch): + """Installed cmdstanpy with no CmdStan is a different problem with a different fix. + + Reporting them separately matters because `pip install` cannot solve the + second one: CmdStan is a C++ build, not a Python package. + """ + cmdstanpy = pytest.importorskip("cmdstanpy") + + def no_cmdstan(): + raise ValueError("No CmdStan directory") + + monkeypatch.setattr(cmdstanpy, "cmdstan_path", no_cmdstan) + + with pytest.raises(ImportError, match="install_cmdstan") as exc: + _import_cmdstanpy() + assert "No CmdStan directory" in str(exc.value), "the underlying reason should survive" + + +def test_cmdstan_is_available_is_false_without_the_package(monkeypatch): + """The gate must answer False, not raise, when cmdstanpy is absent.""" + monkeypatch.setitem(sys.modules, "cmdstanpy", None) + + assert cmdstan_is_available() is False + + +def test_cmdstan_is_available_is_false_without_an_installation(monkeypatch): + """Installing cmdstanpy from PyPI does not install CmdStan, so both are checked. + + Checking only the import would reproduce the original defect in a new + costume: a gate that reports ready for an environment that can only fail. + """ + cmdstanpy = pytest.importorskip("cmdstanpy") + + def no_cmdstan(): + raise ValueError("No CmdStan directory") + + monkeypatch.setattr(cmdstanpy, "cmdstan_path", no_cmdstan) + + assert cmdstan_is_available() is False + -pytestmark = pytest.mark.stan +def test_collection_hook_fails_only_when_cmdstan_is_declared_present(monkeypatch): + """The skip-versus-fail asymmetry that keeps a green CI log honest.""" + monkeypatch.setattr(conftest, "cmdstan_is_available", lambda: False) -requires_pystan = pytest.mark.skipif( - find_spec("pystan") is None, reason="requires the optional pystan dependency" + # Unset: a contributor without CmdStan sees skips, not failures. + monkeypatch.delenv("PYMARE_REQUIRE_CMDSTAN", raising=False) + assert conftest.pytest_collection_modifyitems(None, []) is None + + # Set: the environment claims it can run them, so their absence is an error. + monkeypatch.setenv("PYMARE_REQUIRE_CMDSTAN", "1") + with pytest.raises(pytest.UsageError, match="install_cmdstan"): + conftest.pytest_collection_modifyitems(None, []) + + # Set, and genuinely available: nothing to complain about. + monkeypatch.setattr(conftest, "cmdstan_is_available", lambda: True) + assert conftest.pytest_collection_modifyitems(None, []) is None + + +# ----------------------------------------------------------------------------- +# Translation into Stan's data block. No CmdStan needed. +# ----------------------------------------------------------------------------- + + +def test_stan_data_passes_standard_deviations(): + """Stan's normal() takes a scale, so v must be converted before it is passed.""" + v = np.array([0.04, 0.09, 0.16]) + data = _build_stan_data(np.array([1.0, 2.0, 3.0]), v, np.ones((3, 1))) + + np.testing.assert_allclose(data["sigma"], np.sqrt(v)) + + +def test_stan_data_derives_the_tau_prior_scale_from_the_data(): + """The default prior scale is the larger of the two scales the data carry. + + Taking only the sampling standard deviation understates tau whenever the + between-group spread is the larger of the two, which validation/stan + measures as credible-interval coverage of 0.83 against a nominal 0.95. + """ + v = np.array([0.04, 0.09, 0.16]) + spread_dominates = np.array([1.0, 2.0, 3.0]) + data = _build_stan_data(spread_dominates, v, np.ones((3, 1))) + + assert data["tau_prior_scale"] == pytest.approx(np.std(spread_dominates)) + assert data["tau_prior_scale"] > np.sqrt(np.mean(v)) + + # ... and the sampling error is the floor when the estimates all coincide, + # where the spread alone would be zero and so not a usable scale. + identical = _build_stan_data(np.full(3, 2.0), v, np.ones((3, 1))) + + assert identical["tau_prior_scale"] == pytest.approx(np.sqrt(np.mean(v))) + assert identical["tau_prior_scale"] > 0 + + explicit = _build_stan_data(spread_dominates, v, np.ones((3, 1)), tau_prior_scale=7.0) + assert explicit["tau_prior_scale"] == 7.0 + + +@pytest.mark.parametrize( + "groups", + [ + pytest.param(["a", "a", "b", "b", "c"], id="strings"), + pytest.param([10, 10, 20, 20, 30], id="non-consecutive ints"), + pytest.param(np.array([10, 10, 20, 20, 30]), id="ndarray"), + pytest.param(np.array([[10], [10], [20], [20], [30]]), id="column vector"), + ], ) +def test_stan_data_encodes_arbitrary_group_labels(groups): + """Labels of any hashable type become the 1..K codes the Stan program declares. + + The column-vector case is what ``Dataset.g`` holds, and the ndarray cases + used to raise outright: ``groups or default`` asks a whole array for its + truth value. + """ + y = np.arange(5.0) + data = _build_stan_data(y, np.ones(5), np.ones((5, 1)), groups=groups) + + np.testing.assert_array_equal(data["id"], [1, 1, 2, 2, 3]) + assert data["K"] == 3 + assert data["id"].min() >= 1 and data["id"].max() <= data["K"] + +def test_stan_data_treats_each_observation_as_its_own_group_by_default(): + """Without groups, K equals N and every observation gets a distinct code.""" + data = _build_stan_data(np.arange(5.0), np.ones(5), np.ones((5, 1))) -@requires_pystan -def test_stan_estimator(dataset): - """Run smoke test for StanMetaRegression.""" - # no ground truth here, so we use sanity checks and rough bounds - est = StanMetaRegression(num_samples=3000).fit_dataset(dataset) - results = est.summary() - assert "BayesianMetaRegressionResults" == results.__class__.__name__ - summary = results.summary(["beta", "tau2"]) - beta1, beta2, tau2 = summary["mean"].values[:3] - assert -0.5 < beta1 < 0.1 - assert 0.6 < beta2 < 0.9 - assert 3 < tau2 < 5 + assert data["K"] == data["N"] == 5 + np.testing.assert_array_equal(data["id"], [1, 2, 3, 4, 5]) + + +def test_stan_data_design_matrix_has_one_row_per_observation(): + """X is per-observation, so its row count tracks N and not K.""" + X = np.column_stack([np.ones(6), np.arange(6.0)]) + data = _build_stan_data(np.arange(6.0), np.ones(6), X, groups=[1, 1, 2, 2, 3, 3]) + + assert data["X"].shape == (data["N"], data["C"]) == (6, 2) + assert data["K"] == 3 + assert data["K"] < data["N"] + + +def test_stan_data_promotes_a_one_dimensional_design_matrix(): + """A single predictor may be passed as a 1d array.""" + data = _build_stan_data(np.arange(4.0), np.ones(4), np.arange(4.0)) + + assert data["X"].shape == (4, 1) + assert data["C"] == 1 + + +@pytest.mark.parametrize( + ("kwargs", "match"), + [ + pytest.param({"v": np.ones(3)}, "one sampling variance per observation", id="short v"), + pytest.param({"v": np.array([1.0, 0.0, 1.0, 1.0])}, "must all be positive", id="zero v"), + pytest.param( + {"v": np.array([1.0, -1.0, 1.0, 1.0])}, "must all be positive", id="negative v" + ), + pytest.param({"X": np.ones((3, 1))}, "one row per observation", id="short X"), + ], +) +def test_stan_data_rejects_inputs_that_disagree_about_n(kwargs, match): + """Shape and positivity are checked once, at the boundary.""" + call = {"y": np.arange(4.0), "v": np.ones(4), "X": np.ones((4, 1))} + call.update(kwargs) + + with pytest.raises(ValueError, match=match): + _build_stan_data(**call) def test_stan_2d_input_failure(dataset_2d): """Run smoke test for StanMetaRegression on 2D data. - No pystan needed: the shape is rejected before the model is compiled. + No CmdStan needed: the shape is rejected before the model is compiled. """ with pytest.raises(ValueError) as exc: - StanMetaRegression(num_samples=500).fit_dataset(dataset_2d) + StanMetaRegression().fit_dataset(dataset_2d) assert str(exc.value).startswith("The StanMetaRegression") + + +def test_fit_dataset_forwards_dataset_g(planted_hierarchical_dataset): + """fit_dataset must route dataset.g into fit()'s groups argument. + + Without the mapping the argument falls back to None and every observation is + silently modelled as its own group. Presetting .model with a stub keeps this + free of CmdStan, since fit() only compiles when it finds no model. + """ + dataset, _ = planted_hierarchical_dataset + est = StanMetaRegression() + est.model = _StubModel() + + est.fit_dataset(dataset) + + assert est.data["N"] == 90 + assert est.data["K"] == 30 + np.testing.assert_array_equal(est.data["id"], np.repeat(np.arange(1, 31), 3)) + + +def test_fit_warns_about_divergent_transitions(): + """A divergent fit must warn through the warnings module, not only the log. + + CmdStanPy logs its own diagnostic message, but a log record obeys no + warning filter and cannot be asserted on. Divergences mean the sampler + could not reach part of the posterior, so a summary that looks ordinary may + not be. + """ + est = StanMetaRegression() + est.model = _StubModel(divergences=np.array([3.0, 1.0])) + + with pytest.warns(UserWarning, match="4 divergent transition"): + est.fit(np.arange(4.0), np.ones(4), np.ones((4, 1))) + + +def test_fit_is_quiet_when_there_are_no_divergences(): + """The converse: a clean fit must not cry wolf.""" + est = StanMetaRegression() + est.model = _StubModel() + + with warnings.catch_warnings(): + warnings.simplefilter("error") + est.fit(np.arange(4.0), np.ones(4), np.ones((4, 1))) + + +def test_fake_home_redirects_expanduser_on_every_platform(fake_home): + """The home-directory redirect must hold under Windows path rules too. + + ntpath is importable everywhere, so this catches on any platform the mistake + that only shows up on a Windows runner: setting HOME alone leaves + ntpath.expanduser pointing at the real user profile. + """ + assert posixpath.expanduser("~") == str(fake_home) + assert ntpath.expanduser("~") == str(fake_home) + assert op.expanduser("~") == str(fake_home) + + +def test_compile_falls_back_when_the_package_directory_is_read_only(monkeypatch, fake_home): + """An unwritable site-packages must not make the estimator unusable. + + CmdStanPy compiles beside the .stan source, which lives inside the installed + package, and that directory is read-only in plenty of installations. Note the + exception: CmdStanPy reports every failed ``make`` as ValueError, whatever + went wrong, so a handler for PermissionError would never fire. + """ + cmdstanpy = pytest.importorskip("cmdstanpy") + compiled_from = [] + + def fake_model(stan_file=None, exe_file=None, force_compile=False): + compiled_from.append(stan_file) + if stan_file == STAN_MODEL_PATH: + raise ValueError(f"Failed to compile Stan model '{stan_file}'.") + return "compiled" + + monkeypatch.setattr(cmdstanpy, "cmdstan_path", lambda: str(fake_home)) + monkeypatch.setattr(cmdstanpy, "CmdStanModel", fake_model) + + est = StanMetaRegression() + with pytest.warns(UserWarning, match="not writable"): + est.compile() + + assert est.model == "compiled" + # The second attempt compiles a *copy*, not the packaged file: exe_file + # names an executable to reuse rather than a destination to build into. + assert compiled_from[0] == STAN_MODEL_PATH + assert compiled_from[1] == op.join(str(fake_home), ".pymare", "stan", "meta_regression.stan") + assert op.exists(compiled_from[1]), "the fallback must copy the model somewhere writable" + assert op.getmtime(compiled_from[1]) == op.getmtime( + STAN_MODEL_PATH + ), "copy2 preserves the mtime so the cached build is not invalidated every run" + + +def test_compile_reports_the_original_error_when_the_fallback_also_fails(monkeypatch, fake_home): + """A broken model must not be reported as a permissions problem. + + Both compiles fail for a model that does not parse, and it is the first + error that names the real cause. + """ + cmdstanpy = pytest.importorskip("cmdstanpy") + + def always_fails(stan_file=None, exe_file=None, force_compile=False): + raise ValueError(f"Syntax error in '{stan_file}'") + + monkeypatch.setattr(cmdstanpy, "cmdstan_path", lambda: str(fake_home)) + monkeypatch.setattr(cmdstanpy, "CmdStanModel", always_fails) + + est = StanMetaRegression() + with warnings.catch_warnings(): + warnings.simplefilter("error") # no misleading "not writable" warning + with pytest.raises(ValueError, match="Syntax error") as exc: + est.compile() + assert STAN_MODEL_PATH in str(exc.value), "the packaged path is the one that failed first" + + +class _StubModel: + """Stand in for a compiled CmdStanModel so fit() can run without CmdStan.""" + + def __init__(self, divergences=None): + self.divergences = np.zeros(2) if divergences is None else divergences + + def sample(self, data=None, **kwargs): + """Return an object exposing only what fit() reads off a CmdStanMCMC.""" + return _StubFit(self.divergences) + + +class _StubFit: + """Stand in for a CmdStanMCMC, exposing only its per-chain divergence counts.""" + + def __init__(self, divergences): + self.divergences = divergences + + +# ----------------------------------------------------------------------------- +# Constructor +# ----------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + ("kwarg", "replacement"), + [("num_samples", "iter_sampling"), ("num_warmup", "iter_warmup"), ("num_chains", "chains")], +) +def test_pystan_sampling_kwarg_names_are_rejected(kwarg, replacement): + """Reject PyStan names, which CmdStanPy would report as an unhelpful TypeError.""" + with pytest.raises(TypeError, match=replacement): + StanMetaRegression(**{kwarg: 100}) + + +@pytest.mark.parametrize("bad", [0.0, -1.0, "wide"]) +def test_tau_prior_scale_is_validated_at_construction(bad): + """A prior scale must be a positive number, and None means 'derive it'.""" + with pytest.raises(ValueError): + StanMetaRegression(tau_prior_scale=bad) + + assert StanMetaRegression(tau_prior_scale=None).tau_prior_scale is None + + +def test_summary_before_fit_raises(): + """summary() must not report on an estimator that has not been fitted.""" + with pytest.raises(ValueError, match="hasn't been fitted yet"): + StanMetaRegression().summary() + + +# ----------------------------------------------------------------------------- +# Results container. No CmdStan needed; ArviZ is enough. +# ----------------------------------------------------------------------------- + + +def test_plot_rejects_unknown_kind(): + """An unknown plot kind raises the documented ValueError. + + The guard used to be unreachable: a two-argument getattr raised + AttributeError before the None check below it could run. + """ + pytest.importorskip("arviz") + results = BayesianMetaRegressionResults.__new__(BayesianMetaRegressionResults) + results.data = None + + with pytest.raises(ValueError, match="no plotting function"): + results.plot(kind="not_a_real_plot") + + +def test_plot_defaults_to_the_summary_variables(monkeypatch): + """plot() must not default to plotting one panel per group. + + A fitted model carries one theta per group, so plotting everything is + illegible and, under ArviZ 1.x, exceeds rcParams["plot.max_subplots"] and + raises. The default therefore matches summary()'s selection. + """ + az = pytest.importorskip("arviz") + captured = {} + + def fake_plot_trace(data, var_names=None, **kwargs): + captured["var_names"] = var_names + return "figure" + + monkeypatch.setattr(az, "plot_trace", fake_plot_trace) + results = BayesianMetaRegressionResults.__new__(BayesianMetaRegressionResults) + results.data = None + + assert results.plot(kind="trace") == "figure" + assert captured["var_names"] == ["beta", "tau2"] + + results.plot(kind="trace", include_theta=True) + assert captured["var_names"] == ["beta", "tau2", "theta"] + + # An explicit selection still wins. + results.plot(kind="trace", var_names=["tau2"]) + assert captured["var_names"] == ["tau2"] + + +def test_plot_omits_var_names_for_plotters_that_reject_it(monkeypatch): + """Not every ArviZ plot takes var_names; plot_energy takes only the data.""" + az = pytest.importorskip("arviz") + captured = {} + + def fake_plot_energy(data, **kwargs): + captured["kwargs"] = kwargs + return "figure" + + # plot_energy genuinely has no var_names parameter, and neither does this + # stub, so the selection above must not be forced onto it. + + monkeypatch.setattr(az, "plot_energy", fake_plot_energy) + results = BayesianMetaRegressionResults.__new__(BayesianMetaRegressionResults) + results.data = None + + assert results.plot(kind="energy") == "figure" + assert "var_names" not in captured["kwargs"] + + +@pytest.mark.parametrize("bad_ci", [0, 100, -5, 101]) +def test_results_reject_an_impossible_credible_interval(bad_ci): + """The ci argument is a percentage, so it must lie strictly inside (0, 100).""" + pytest.importorskip("arviz") + with pytest.raises(ValueError, match="must lie in"): + BayesianMetaRegressionResults(None, None, ci=bad_ci) + + +@pytest.mark.parametrize("major", [0, 1]) +def test_credible_interval_kwargs_track_the_installed_arviz(monkeypatch, major): + """Both ArviZ generations must be asked for the same interval. + + 1.x renamed hdi_prob to ci_prob, defaults to an equal-tailed interval rather + than a highest-density one, and formats summaries as strings unless told not + to. Only one branch is reachable on any given install, so the other is + exercised by pinning the reported version. + """ + az = pytest.importorskip("arviz") + monkeypatch.setattr(az, "__version__", f"{major}.3.0") + + kwargs = _arviz_credible_interval_kwargs(95.0) + + if major >= 1: + assert kwargs == {"ci_prob": 0.95, "ci_kind": "hdi", "round_to": "none"} + else: + assert kwargs == {"hdi_prob": 0.95} + + +def test_accepts_var_names_handles_an_unreadable_signature(): + """Anything whose signature cannot be read must be treated as not taking var_names. + + inspect.signature raises rather than answering for some objects, and plot() + calls this before deciding what to pass, so an exception here would surface + as a broken plot rather than as a missing default. + """ + pytest.importorskip("arviz") + + assert _accepts_var_names(lambda data, var_names=None: None) is True + # Reads fine, simply has no such parameter. + assert _accepts_var_names(lambda data: None) is False + # Cannot be read at all: inspect raises TypeError for a non-callable. + assert _accepts_var_names(object()) is False + + +def test_results_accept_a_converted_object_without_cmdstanpy(monkeypatch): + """The container must work when cmdstanpy is absent but the data is already converted. + + A caller who has their own InferenceData should not need the sampler + installed to summarize it, so the import is lazy and its failure is not one. + """ + pytest.importorskip("arviz") + monkeypatch.setitem(sys.modules, "cmdstanpy", None) + + results = BayesianMetaRegressionResults("already-converted", None, ci=90.0) + + assert results.data == "already-converted" + assert results.ci == 90.0 + + +@pytest.mark.parametrize("ci", [50.0, 95.0]) +def test_summary_requests_the_configured_credible_interval(ci): + """The ci argument must reach ArviZ. It was previously stored and never used.""" + az = pytest.importorskip("arviz") + kwargs = _arviz_credible_interval_kwargs(ci) + probability = kwargs.get("ci_prob", kwargs.get("hdi_prob")) + + assert probability == pytest.approx(ci / 100.0) + if int(az.__version__.split(".")[0]) >= 1: + # ArviZ 1.x defaults to an equal-tailed interval and to stringifying the + # summary for display; neither is what this container promises. + assert kwargs["ci_kind"] == "hdi" + assert kwargs["round_to"] == "none" + + +# ----------------------------------------------------------------------------- +# The recorded simulation results. No CmdStan needed; the numbers are pinned. +# ----------------------------------------------------------------------------- + + +def test_recorded_validation_covers_every_design_cell(): + """The pinned results must describe the grid the harness actually runs. + + Without this, adding a cell to validation/stan/simulate.py and forgetting to + regenerate would leave the new cell permanently unmeasured while the file + still looked current. + """ + recorded = load_stan_validation() + names = [cell["name"] for cell in recorded["cells"]] + + assert tuple(names) == STAN_VALIDATION_CELLS + assert len(names) == len(set(names)), "duplicate cells in the recorded results" + assert recorded["replications"] >= 100 + + +def test_recorded_validation_meets_its_thresholds(): + """Every design cell must clear the coverage floor and the bias ceiling. + + This is what makes the recorded file load-bearing rather than decorative. It + reads the numbers already measured rather than re-measuring, so it costs + nothing and runs everywhere; the scheduled Stan validation workflow is what + re-measures and applies the same thresholds to a fresh run. + """ + recorded = load_stan_validation() + floor = STAN_VALIDATION_THRESHOLDS["min_coverage"] + ceiling = STAN_VALIDATION_THRESHOLDS["max_beta_bias"] + + undercovered = [ + (cell["name"], cell["beta_coverage"]) + for cell in recorded["cells"] + if cell["beta_coverage"] < floor + ] + assert not undercovered, f"cells below {floor:.2f} coverage: {undercovered}" + + biased = [ + (cell["name"], cell["beta_bias"]) + for cell in recorded["cells"] + if abs(cell["beta_bias"]) > ceiling + ] + assert not biased, f"cells with |beta bias| above {ceiling}: {biased}" + + +def test_recorded_validation_summarizes_the_worst_coefficient(): + """The reported figures must be the worst coefficient, not an average of them. + + Averaging across coefficients lets a well-estimated intercept mask a badly + estimated moderator, which is exactly the failure the unbalanced-covariate + cells exist to detect. The thresholds are therefore applied to the minimum + coverage and the largest absolute bias across coefficients, and this pins + that so the summary cannot quietly become a mean. + """ + recorded = load_stan_validation() + + for cell in recorded["cells"]: + per_coverage = cell["beta_coverage_per_coefficient"] + per_bias = cell["beta_bias_per_coefficient"] + + assert len(per_coverage) == len(per_bias) == len(cell["true_beta"]) + assert cell["beta_coverage"] == pytest.approx(min(per_coverage)) + assert cell["beta_bias"] == pytest.approx(max(per_bias, key=abs)) + + +def test_recorded_validation_used_a_fixed_truth(): + """Bias is only meaningful if the coefficients being recovered are held fixed. + + Under a redrawn symmetric truth the signed errors average to zero for any + estimator at all, including one that always returns zero, so the bias + ceiling would certify nothing. + """ + recorded = load_stan_validation() + truths = {tuple(cell["true_beta"]) for cell in recorded["cells"]} + + # Every cell draws its coefficients from the same fixed vector, truncated to + # however many predictors that cell uses. + longest = max(truths, key=len) + for truth in truths: + assert truth == longest[: len(truth)] + assert all(value != 0 for value in truth), "a zero coefficient cannot show bias" + + +# ----------------------------------------------------------------------------- +# Sampling. Needs CmdStan. +# ----------------------------------------------------------------------------- + + +@pytest.mark.stan +@requires_cmdstan +def test_recovers_planted_parameters(planted_hierarchical_dataset): + """The fitted posterior must recover the parameters the data were built from. + + This is the test that distinguishes tau from tau2 and variances from + standard deviations: the planted tau2 of 0.25 and tau of 0.5 are far enough + apart that reporting either under the other's name lands outside the + interval asserted here. + """ + dataset, truth = planted_hierarchical_dataset + est = StanMetaRegression(iter_sampling=2000, chains=4, seed=8675309, show_progress=False) + results = est.fit_dataset(dataset).summary() + summary = results.summary() + + for i, expected in enumerate(truth["beta"]): + mean = float(summary.loc[f"beta[{i}]", "mean"]) + sd = float(summary.loc[f"beta[{i}]", "sd"]) + assert abs(mean - expected) < 3 * sd, f"beta[{i}] = {mean}, expected ~{expected}" + + tau2 = float(summary.loc["tau2", "mean"]) + assert 0.15 < tau2 < 0.45, f"tau2 = {tau2}, expected ~{truth['tau2']}" + assert tau2 < truth["tau"], "tau2 looks like tau, not tau squared" + + +@pytest.mark.stan +@requires_cmdstan +def test_matches_maximum_likelihood_without_groups(dataset): + """Ungrouped, the model marginalizes to the one ML already maximizes. + + With every observation in its own group the hierarchical model collapses to + ``y ~ N(X beta, sqrt(v + tau2))``, which is exactly + VarianceBasedLikelihoodEstimator's ML likelihood. Under a diffuse prior on + tau the posterior means must therefore agree with the ML estimates up to + Monte Carlo error, which pins the Stan program against an implementation + that shares none of its code. + """ + est = StanMetaRegression( + tau_prior_scale=100.0, iter_sampling=4000, chains=4, seed=20250818, show_progress=False + ) + summary = est.fit_dataset(dataset).summary().summary() + + ml = VarianceBasedLikelihoodEstimator(method="ML").fit_dataset(dataset) + ml_beta = np.asarray(ml.params_["fe_params"]).ravel() + + for i, expected in enumerate(ml_beta): + posterior_mean = float(summary.loc[f"beta[{i}]", "mean"]) + assert posterior_mean == pytest.approx(expected, abs=0.15 * max(abs(expected), 1.0)) + + +@pytest.mark.stan +@requires_cmdstan +def test_meta_regression_dispatches_to_stan(planted_hierarchical_dataset): + """The functional entry point must reach this estimator and its results class. + + pymare.meta_regression(method="stan") is the documented one-call API and has + its own dispatch table in core.py, which no test previously exercised. + """ + dataset, _ = planted_hierarchical_dataset + + results = meta_regression( + y=dataset.y, + v=dataset.v, + X=dataset.X[:, 1:], + g=dataset.g, + method="stan", + iter_sampling=500, + chains=2, + seed=99, + show_progress=False, + ) + + assert isinstance(results, BayesianMetaRegressionResults) + assert list(results.summary().index) == ["beta[0]", "beta[1]", "tau2"] + + +@pytest.mark.stan +@requires_cmdstan +def test_the_compiled_model_is_reused_across_fits(planted_hierarchical_dataset): + """compile() once, fit many. + + The compiled executable does not depend on the data, so fitting must reuse it + rather than rebuilding per call. + """ + dataset, _ = planted_hierarchical_dataset + est = StanMetaRegression(iter_sampling=200, chains=1, seed=5, show_progress=False) + + est.compile() + compiled = est.model + assert compiled is not None + + est.fit_dataset(dataset) + assert est.model is compiled + + est.fit_dataset(dataset) + assert est.model is compiled + + +@pytest.mark.stan +@requires_cmdstan +def test_summary_and_plot_round_trip(planted_hierarchical_dataset): + """The results container reports the expected rows and plots without error.""" + dataset, _ = planted_hierarchical_dataset + est = StanMetaRegression(iter_sampling=500, chains=2, seed=1234, show_progress=False) + results = est.fit_dataset(dataset).summary() + + assert isinstance(results, BayesianMetaRegressionResults) + + without_theta = results.summary() + assert list(without_theta.index) == ["beta[0]", "beta[1]", "tau2"] + + with_theta = results.summary(include_theta=True) + assert len(with_theta) == len(without_theta) + 30 + + assert results.plot(kind="trace") is not None + + +@pytest.mark.parametrize( + ("field", "bad"), + [ + ("v", np.array([1.0, np.nan, 1.0, 1.0])), + ("v", np.array([1.0, np.inf, 1.0, 1.0])), + ("y", np.array([1.0, np.nan, 1.0, 1.0])), + ("X", np.array([[1.0], [np.nan], [1.0], [1.0]])), + ], +) +def test_stan_data_rejects_non_finite_inputs(field, bad): + """Reject NaN here, rather than leaving CmdStan to refuse it while reading data. + + NaN fails every comparison, so a positivity check alone lets it through: + ``np.nan <= 0`` is False. It then flows into sqrt() and into the prior + scale, and only surfaces when CmdStan refuses to load the data -- an error + that names a Stan variable rather than the input that caused it. + """ + call = {"y": np.arange(4.0), "v": np.ones(4), "X": np.ones((4, 1))} + call[field] = bad + + with pytest.raises(ValueError, match="must all be finite"): + _build_stan_data(**call) + + +def test_stan_data_still_rejects_non_positive_variances(): + """The finiteness check must not have displaced the positivity one.""" + with pytest.raises(ValueError, match="must all be positive"): + _build_stan_data(np.arange(4.0), np.array([1.0, 0.0, 1.0, 1.0]), np.ones((4, 1))) + + +def test_stan_data_rejects_composite_group_labels(): + """The documented contract is scalar labels, and this is why. + + Numpy reads a sequence of tuples as a 2-dimensional array, so a tuple label + is not a label at all by the time encode_groups sees it. The docstring says + scalar rather than hashable for this reason. + """ + with pytest.raises(ValueError, match="one-dimensional"): + _build_stan_data( + np.arange(3.0), + np.ones(3), + np.ones((3, 1)), + groups=[("a", 1), ("a", 1), ("b", 2)], + ) diff --git a/pymare/tests/utils.py b/pymare/tests/utils.py index 57f477a..b38f1da 100644 --- a/pymare/tests/utils.py +++ b/pymare/tests/utils.py @@ -42,3 +42,109 @@ def load_robumeta_reference(): """ with open(op.join(get_test_data_path(), "robumeta_reference.json")) as fobj: return json.load(fobj) + + +def cmdstan_is_available(): + """Report whether the Stan estimator can actually be run here. + + Returns + ------- + :obj:`bool` + True when ``cmdstanpy`` imports *and* it can find a CmdStan + installation. + + Notes + ----- + Both halves matter. ``cmdstanpy`` installs cleanly from PyPI without + CmdStan, which is a C++ build rather than a Python package, so an import + check alone would report an environment as ready when it can only fail. + + The failure this replaced was subtler still: the previous gate probed + ``find_spec("pystan")``, but PyStan 3 is distributed as ``pystan`` and + imported as ``stan``, so the probe was unsatisfiable and the estimator's + only real test skipped everywhere, including in the CI job that existed to + run it. A skip reads as a pass in a CI log, which is why the + ``pytest_collection_modifyitems`` hook in ``conftest.py`` consults this + function and fails the run wherever ``PYMARE_REQUIRE_CMDSTAN`` says Stan is + expected. + """ + try: + import cmdstanpy + except ImportError: + return False + + try: + cmdstanpy.cmdstan_path() + except Exception: + return False + + return True + + +#: Thresholds the Stan model's simulated performance has to meet. Applied in two +#: places against the same numbers: ``test_stan_estimators.py`` checks them +#: against the pinned ``data/stan_validation.json``, and +#: ``validation/stan/simulate.py --check`` checks them against a fresh run. +#: +#: ``min_coverage`` is a floor on the *worst* coefficient's share of 95% credible +#: intervals containing the true value; averaging across coefficients would let a +#: well-estimated intercept mask a badly estimated moderator. It sits well below +#: the nominal 0.95 because a minimum over coefficients is biased downward. 0.85 +#: is between the worst honest measurement and the failures it has to reject -- +#: ``validation/stan/README.md`` records both. +#: +#: ``max_beta_bias`` bounds the largest absolute bias across coefficients. It only +#: means anything because the harness holds the true coefficients fixed; with a +#: redrawn symmetric truth the signed errors average to zero for any estimator. +STAN_VALIDATION_THRESHOLDS = { + "min_coverage": 0.85, + "max_beta_bias": 0.10, +} + +#: The design cells ``validation/stan/simulate.py`` is expected to report on. +#: Pinned here so that adding or renaming a cell fails the test that reads the +#: recorded results until the pinned file is regenerated, rather than quietly +#: leaving the new cell unchecked. +STAN_VALIDATION_CELLS = ( + "groups=5", + "groups=20", + "groups=50", + "tau2=0", + "tau2=0.1", + "tau2=1", + "singletons", + "unequal groups", + "sigma x0.1", + "sigma x10", + "1 predictor", + "3 predictors", + "unbalanced covariate", + "unbalanced covariate, tau2=1", +) + + +def load_stan_validation(): + """Load the recorded simulation results for the Stan meta-regression model. + + Returns + ------- + :obj:`dict` + The parsed contents of ``data/stan_validation.json``: the provenance of + the run and a ``"cells"`` list holding one entry per design cell, with + that cell's measured bias, credible-interval coverage and divergence + count. + + Notes + ----- + Regenerated by ``validation/stan/simulate.py``, which takes about ten + minutes, so the results are pinned rather than recomputed on every run -- + the same arrangement as :func:`load_robumeta_reference`. + + Unlike the robumeta reference, these numbers are stochastic, so the pin + cannot be enforced by requiring the file not to move; a correct estimator + produces slightly different numbers every run. What is enforced instead is + ``STAN_VALIDATION_THRESHOLDS``, which is the claim the file exists to + support. + """ + with open(op.join(get_test_data_path(), "stan_validation.json")) as fobj: + return json.load(fobj) diff --git a/pyproject.toml b/pyproject.toml index 2aa8afc..de11052 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,7 +5,7 @@ requires = ["setuptools==68.2.2", "wheel"] testpaths = ["pymare/tests"] addopts = "--strict-markers" markers = [ - "stan: tests that need the optional pystan dependency (slow: the model is compiled)", + "stan: tests that sample, needing cmdstanpy and a CmdStan installation (slow: the model is compiled)", "robumeta: tests that pin PyMARE against the R package robumeta", ] diff --git a/setup.cfg b/setup.cfg index 79298eb..664b7c3 100644 --- a/setup.cfg +++ b/setup.cfg @@ -25,11 +25,9 @@ classifiers = [options] python_requires = >= 3.9 install_requires = - numpy>=1.8.0,<2.0; python_version == "3.9" and extra == 'stan' - numpy>=1.8.0; python_version != "3.9" or extra != 'stan' + numpy>=1.8.0 pandas - scipy<1.13; python_version == "3.9" and extra == 'stan' - scipy; python_version != "3.9" or extra != 'stan' + scipy sympy packages = find: include_package_data = False @@ -67,8 +65,14 @@ tests = pytest pytest-cov stan = - pystan - arviz + cmdstanpy>=1.2,<2 + arviz>=0.17 + # The newest ArviZ that still supports Python 3.9 is 0.17, which predates + # the removal of scipy.signal.gaussian in scipy 1.13 and imports it + # unconditionally. Constraining scipy here rather than in install_requires + # keeps the cap on the extra that needs it, instead of on every PyMARE + # install that happens to run on 3.9. + scipy<1.13; python_version < "3.10" all = %(doc)s %(tests)s @@ -80,6 +84,11 @@ all = resources/* resources/datasets/* effectsize/*.json +# The Stan program is compiled on first use by CmdStanPy, which needs it on +# disk. include_package_data is False, so this is the only thing that puts it +# in the wheel. +pymare.estimators = + stan/*.stan [versioneer] VCS = git diff --git a/validation/stan/README.md b/validation/stan/README.md new file mode 100644 index 0000000..7c6b153 --- /dev/null +++ b/validation/stan/README.md @@ -0,0 +1,191 @@ +# Validating the Stan meta-regression model + +What `pymare/estimators/stan/meta_regression.stan` claims, and what it was +measured to do. This mirrors `validation/robumeta`, except that the reference +here is the data-generating process rather than another package: the model is +fitted to data simulated from itself, and checked for whether it recovers what +was planted and whether its credible intervals cover at their nominal rate. + +The fast tests in `pymare/tests/test_stan_estimators.py` check one planted +configuration; this checks the grid. It is not run per pull request, but it is +not detached from CI either — see "How this is wired into CI" below. + +## The model + +``` +y_i ~ normal(x_i' beta + theta_{g(i)}, sigma_i) i = 1..N +theta_g ~ normal(0, tau) g = 1..K +tau ~ normal(0, tau_prior_scale) truncated at 0 +beta ~ (improper uniform) +``` + +`sigma_i = sqrt(v_i)` is the known sampling **standard deviation**; `tau2 = +tau^2` is the reported between-group variance. `theta` is non-centered +(`theta = tau * theta_raw`, `theta_raw ~ std_normal()`). + +This is the Stan User's Guide random-effects meta-analysis model +([Measurement Error and Meta-Analysis][sug]) with that guide's stated extension +to observation-level predictors. The half-normal prior on `tau` follows Stan's +current [prior choice recommendations][priors], which supersede the `cauchy(0, 5)` +still shown in the guide. + +[sug]: https://mc-stan.org/docs/stan-users-guide/measurement-error.html +[priors]: https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations + +## How this is wired into CI + +Three layers, mirroring the robumeta alignment in `validation/robumeta`: + +1. **`make validate_stan`** re-measures the grid and fails if any design cell + misses `pymare.tests.utils.STAN_VALIDATION_THRESHOLDS`. It writes the results + to `pymare/tests/data/stan_validation.json`. +2. **`test_recorded_validation_meets_its_thresholds`** and + **`test_recorded_validation_covers_every_design_cell`** hold that recorded + file to the same thresholds and to the expected list of cells. They read the + file rather than re-measuring, so they cost nothing and run on every pull + request on every platform. This is what stops the record from going stale + unnoticed — a pinned file nothing reads is decoration. +3. **The `Validate the Stan model` workflow** re-measures on a schedule, on + pushes to master touching the model, and on demand. + +The grid is about 1400 fits and takes ten minutes, which is why it is not part +of `test_stan` and not run per pull request. + +**Why thresholds rather than pinned values.** The robumeta reference is +deterministic, so its workflow can require the file not to move. These numbers +are Monte Carlo estimates with a standard error of 0.015 to 0.030, so a correct +model produces different numbers every run and an exact pin would fail +constantly. What is pinned instead is the claim the file exists to support: +worst-coefficient coverage at or above 0.85 and |beta bias| at or below 0.10 in +every cell. That floor is set from measurement: the correct model's tightest +cell reads 0.900 while the rejected prior reads 0.710 and 0.830, so 0.85 sits +between them. A minimum over coefficients is biased downward, which is why the +floor is not nearer the nominal 0.95. +`--check` refuses to certify a run of fewer than 100 replications, so a short run +cannot clear the floor by luck. + +## Reproducing + +```bash +pip install -e .[stan] +python -m cmdstanpy.install_cmdstan +make validate_stan # or: python validation/stan/simulate.py --check +``` + +About 10 minutes on 8 cores. + +## What was measured + +100 replications per cell, seed 20260818, CmdStan 2.36.0, 2 chains x 1000 draws. +Each cell varies one factor away from a base of 20 groups of 3, `tau2 = 0.1`, 2 +predictors, sampling SDs drawn from `uniform(0.1, 0.4)`. + +The true coefficients are **fixed** at `[0.5, -0.8, 0.3]`, truncated to the +number of predictors, rather than redrawn each replication. That matters: with a +symmetric redrawn truth the signed errors average to zero for any estimator at +all, so the bias threshold would certify nothing — an estimator that always +returned zero cleared it about 85% of the time. + +Coverage is reported **per coefficient**, and the threshold is applied to the +**worst** of them rather than the average. Averaging lets a well-estimated +intercept mask a badly estimated moderator, which is the failure the unbalanced +cells exist to detect — and it did mask it: the rejected prior scale below reads +0.810 pooled but 0.710 on its worst coefficient. + +The two coverage columns are the two candidate defaults for `tau_prior_scale`, +both measured under this per-coefficient metric. + +| cell | worst cov, `sqrt(mean(v))` | worst cov, `max(std(y), sqrt(mean(v)))` | per coefficient | tau2 bias, old | tau2 bias, new | true tau2 | worst beta bias | divergent fits | +| --- | --- | --- | --- | --- | --- | --- | --- | --- | +| groups=5 | 0.910 | **0.930** | 0.93, 0.93 | +0.000 | +0.122 | 0.10 | +0.0092 | 54 | +| groups=20 | 0.900 | **0.900** | 0.91, 0.90 | -0.002 | +0.010 | 0.10 | -0.0057 | 0 | +| groups=50 | 0.960 | **0.960** | 0.99, 0.96 | +0.004 | +0.008 | 0.10 | -0.0047 | 0 | +| tau2=0 | 0.960 | **0.960** | 0.97, 0.96 | +0.004 | +0.004 | 0.00 | +0.0048 | 14 | +| tau2=0.1 | 0.930 | **0.930** | 0.94, 0.93 | +0.003 | +0.018 | 0.10 | -0.0036 | 0 | +| tau2=1 | 0.910 | **0.960** | 0.96, 0.97 | -0.285 | +0.185 | 1.00 | +0.0103 | 0 | +| singletons | 0.940 | **0.950** | 0.98, 0.95 | -0.003 | +0.021 | 0.10 | -0.0136 | 1 | +| unequal groups | 0.930 | **0.940** | 0.94, 0.99 | +0.008 | +0.024 | 0.10 | -0.0042 | 0 | +| sigma x0.1 | 0.710 | **0.940** | 0.94, 0.95 | -0.069 | +0.021 | 0.10 | +0.0027 | 0 | +| sigma x10 | 0.940 | **0.940** | 0.98, 0.94 | +0.416 | +0.421 | 0.10 | -0.0450 | 6 | +| 1 predictor | 0.940 | **0.940** | 0.94 | +0.000 | +0.008 | 0.10 | -0.0045 | 0 | +| 3 predictors | 0.920 | **0.920** | 0.94, 0.92, 0.99 | +0.005 | +0.021 | 0.10 | -0.0099 | 0 | +| unbalanced covariate | 0.960 | **0.950** | 0.97, 0.95 | +0.008 | +0.020 | 0.10 | -0.0086 | 0 | +| unbalanced covariate, tau2=1 | 0.830 | **0.930** | 0.98, 0.93 | -0.346 | +0.049 | 1.00 | -0.0242 | 0 | + +`beta` is unbiased throughout: the largest bias on any coefficient in any cell +is 0.045, against coefficients of order 1. + +## The choice of `tau_prior_scale`, decided by measurement + +The first default tried was `sqrt(mean(v))`, the typical sampling standard +deviation. The grid rejected it. Worst-coefficient coverage fell to **0.710** +when the sampling SDs were small relative to the between-group spread +(`sigma x0.1`) and to **0.830** in the unbalanced-covariate cell at `tau2=1`. In +each, `tau2` was badly *under*-estimated: -69% and -35% respectively. + +The cause is that `sqrt(mean(v))` measures sampling noise, which is not the +quantity `tau` describes. When heterogeneity is much larger than sampling error +the prior is far too tight, `tau` is shrunk toward zero, the uncertainty it +contributes to `beta` is understated, and the intervals are too narrow. + +A direct comparison on the failing cells (60 replications each): + +| cell | `sqrt(mean(v))` | `std(y)` | `max` of both | +| --- | --- | --- | --- | +| base tau2=0.1 | 0.883 | 0.950 | **0.950** | +| sigma x0.1 | 0.833 | 0.917 | **0.958** | +| sigma x10 | 0.975 | 0.950 | **0.967** | +| tau2=1 | 0.950 | 0.908 | **0.967** | +| tau2=0 | 0.975 | 0.975 | **0.967** | + +`max(std(y), sqrt(mean(v)))` is the only candidate at or above nominal +everywhere, and it is what the estimator now uses. The reasoning the numbers +support: + +- **The errors are asymmetric.** A scale that is too small costs *coverage*, + which is a correctness failure. A scale that is too large costs only + precision in `tau2` while coverage stays at nominal — visible in `sigma x10`, + where `tau2` is inflated by 0.39 but coverage is 0.955. A default should + therefore err large. +- **Both terms are needed.** `tau` is the standard deviation of the group means, + so it cannot plausibly exceed the spread of the estimates, which makes + `std(y)` the natural scale. But `std(y)` alone is zero when every estimate + coincides, and zero is not a usable scale; `sqrt(mean(v))` is the floor that + prevents it. + +## Known limits + +- **Five groups is not enough to identify `tau`.** At `groups=5`, 63 of 100 fits + reported divergent transitions and `tau2` was over-estimated by 0.11 against a + truth of 0.10. The wider prior made this worse than the rejected default did + (20 fits), which is the honest cost of the change: with five groups the data + say very little about the group-level variance, so the posterior follows + whatever the prior says. The estimator warns on divergences, so this surfaces + rather than passing silently. Prefer a non-Bayesian estimator, or supply + `tau_prior_scale` from external knowledge, when groups are this few. +- **`tau` is not identified when sampling error dwarfs it.** In `sigma x10` the + sampling SDs are 1 to 4 while `tau` is 0.32, and `tau2` is over-estimated by + 0.39 under every candidate prior. Nothing in the data distinguishes a small + `tau` from zero at that noise level. Coverage for `beta` is unaffected. +- Coverage is measured for `beta` only. `tau2` is reported as bias, not + coverage, because its posterior is strongly skewed at small `K`. + +## stanc pedantic mode + +`stanc --warn-pedantic` reports exactly two warnings, both expected: + +``` +Warning: The parameter tau has no priors. This means either no prior is + provided, or the prior(s) depend on data variables. In the later case, + this may be a false positive. +Warning: The parameter beta has no priors. ... +``` + +The first is the false positive its own text describes: `tau` does have a prior, +but its scale is a data variable, which the check cannot see through. The second +is accurate and deliberate — `beta` keeps Stan's implicit improper uniform prior, +which is what makes the posterior means agree with maximum likelihood. That +agreement is pinned by `test_matches_maximum_likelihood_without_groups`. + +Because one of the two warnings is a true positive by design, the CI step that +runs pedantic mode reports rather than gates. diff --git a/validation/stan/simulate.py b/validation/stan/simulate.py new file mode 100644 index 0000000..91ac05f --- /dev/null +++ b/validation/stan/simulate.py @@ -0,0 +1,324 @@ +"""Measure the Stan estimator's bias and credible-interval coverage. + +Not run in CI. The fast tests in ``pymare/tests/test_stan_estimators.py`` check +that the estimator recovers one planted configuration; this checks that it does +so across the grid of designs where meta-regression estimators are known to +break, and reports coverage, which is the property the model's correctness +actually rests on and which no single fit can establish. + +Run with:: + + make validate_stan + +or equivalently:: + + python validation/stan/simulate.py --check + +which writes ``pymare/tests/data/stan_validation.json`` and exits non-zero if any +design cell misses ``pymare.tests.utils.STAN_VALIDATION_THRESHOLDS``. Two tests +in ``pymare/tests/test_stan_estimators.py`` hold that recorded file to the same +thresholds on every run, so the record cannot go stale unnoticed, and the +"Validate the Stan model" workflow re-measures on a schedule. + +This directory's README.md explains the arrangement and records what was +measured. +""" + +import argparse +import collections +import json +import logging +import os +import os.path as op +import sys +import time +from concurrent.futures import ProcessPoolExecutor + +import cmdstanpy +import numpy as np + +sys.path.insert(0, op.join(op.dirname(op.abspath(__file__)), "..", "..")) + +from pymare import Dataset # noqa: E402 +from pymare.estimators import StanMetaRegression # noqa: E402 +from pymare.tests.utils import STAN_VALIDATION_THRESHOLDS # noqa: E402 + +# CmdStanPy narrates every chain through its own handler; at a few thousand fits +# that is the only thing on screen. +cmdstanpy.disable_logging() +logging.getLogger("cmdstanpy").setLevel(logging.ERROR) + +#: One entry per design cell, each varying a single factor away from BASE, so a +#: failure is attributable to that factor. A full factorial would cost several +#: times the fits and still need reading this way. +CELLS = [ + # Number of groups: the axis the prior on tau is most sensitive to. + {"name": "groups=5", "n_groups": 5}, + {"name": "groups=20", "n_groups": 20}, + {"name": "groups=50", "n_groups": 50}, + # Between-group variance, including the boundary at zero where the funnel + # geometry is worst. + {"name": "tau2=0", "tau2": 0.0}, + {"name": "tau2=0.1", "tau2": 0.1}, + {"name": "tau2=1", "tau2": 1.0}, + # Group structure. + {"name": "singletons", "group_size": 1}, + {"name": "unequal groups", "group_size": "unequal"}, + # Scale of the variance column. A scale-dependent bug reads as fine at one + # scale, which is why a fixed prior scale would not survive this row. + {"name": "sigma x0.1", "sigma_scale": 0.1}, + {"name": "sigma x10", "sigma_scale": 10.0}, + # Predictor count. + {"name": "1 predictor", "n_predictors": 1}, + {"name": "3 predictors", "n_predictors": 3}, + # The cell usually omitted: a group-level moderator carried by only a + # handful of groups, which is where robust variance estimators break down. + {"name": "unbalanced covariate", "unbalanced": True}, + {"name": "unbalanced covariate, tau2=1", "unbalanced": True, "tau2": 1.0}, +] + +#: Fewest replications at which ``--check`` is allowed to pass. Matches the +#: floor the pinned results are held to in ``test_stan_estimators.py``. +MIN_REPLICATIONS_TO_CHECK = 100 + +BASE = { + "n_groups": 20, + "group_size": 3, + "tau2": 0.1, + "sigma_scale": 1.0, + "n_predictors": 2, + "unbalanced": False, +} + + +#: True coefficients, held fixed across replications. Redrawing them from a +#: symmetric distribution makes bias uninformative: an estimator that always +#: returned zero would have errors of -beta, averaging to zero over replications, +#: and would clear any bias threshold. +TRUE_BETA = np.array([0.5, -0.8, 0.3]) + + +def simulate(rng, n_groups, group_size, tau2, sigma_scale, n_predictors, unbalanced): + """Draw one dataset from the model the Stan program encodes.""" + if group_size == "unequal": + sizes = rng.integers(1, 6, size=n_groups) + else: + sizes = np.full(n_groups, group_size) + groups = np.repeat(np.arange(n_groups), sizes) + n_observations = groups.size + + if unbalanced: + # A group-level moderator switched on for only three groups. + carriers = rng.choice(n_groups, size=min(3, n_groups), replace=False) + moderators = np.isin(groups, carriers).astype(float)[:, None] + for _ in range(n_predictors - 2): + moderators = np.column_stack([moderators, rng.normal(size=n_observations)]) + else: + moderators = rng.normal(size=(n_observations, max(n_predictors - 1, 0))) + + X = ( + np.column_stack([np.ones(n_observations), moderators]) + if moderators.size + else np.ones((n_observations, 1)) + ) + beta = TRUE_BETA[: X.shape[1]] + + theta = rng.normal(0, np.sqrt(tau2), size=n_groups) + sigma = rng.uniform(0.1, 0.4, size=n_observations) * sigma_scale + y = X @ beta + theta[groups] + rng.normal(0, sigma) + + predictors = X[:, 1:] if X.shape[1] > 1 else None + dataset = Dataset(y=y, v=sigma**2, X=predictors, g=groups) + return dataset, beta, tau2 + + +def run_cell(cell, replications, seed): + """Fit one design cell `replications` times and accumulate the summaries.""" + config = dict(BASE) + config.update({k: v for k, v in cell.items() if k != "name"}) + rng = np.random.default_rng(seed) + + # Per coefficient, not pooled: in the unbalanced cells the sparse moderator + # is the coefficient at risk, and intercept coverage would mask it. + beta_errors = collections.defaultdict(list) + covered = collections.defaultdict(list) + tau2_errors, tau2_truth, divergent = [], [], 0 + + for replication in range(replications): + dataset, beta, tau2 = simulate(rng, **config) + est = StanMetaRegression( + iter_sampling=1000, + iter_warmup=1000, + chains=2, + seed=int(rng.integers(1, 2**31 - 1)), + show_progress=False, + ) + import warnings + + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + summary = est.fit_dataset(dataset).summary().summary() + divergent += any("divergent" in str(w.message) for w in caught) + + for i, true_value in enumerate(beta): + row = summary.loc[f"beta[{i}]"] + beta_errors[i].append(float(row["mean"]) - true_value) + lower, upper = _interval(row) + covered[i].append(bool(lower <= true_value <= upper)) + + tau2_errors.append(float(summary.loc["tau2", "mean"]) - tau2) + tau2_truth.append(tau2) + + print(f" done: {cell['name']}", flush=True) + return { + "name": cell["name"], + "config": {k: str(v) for k, v in config.items()}, + "replications": replications, + "true_beta": [float(b) for b in TRUE_BETA[: len(_coefficients(covered))]], + "beta_bias_per_coefficient": _per_coefficient(beta_errors), + "beta_coverage_per_coefficient": _per_coefficient(covered), + # The figures the thresholds are applied to summarize the *worst* + # coefficient rather than the average of them. + "beta_bias": max(_per_coefficient(beta_errors), key=abs), + "beta_coverage": min(_per_coefficient(covered)), + "coverage_se": float( + np.sqrt( + min(_per_coefficient(covered)) + * (1 - min(_per_coefficient(covered))) + / len(covered[_coefficients(covered)[0]]) + ) + ), + "tau2_bias": float(np.mean(tau2_errors)), + "tau2_truth": float(np.mean(tau2_truth)), + "fits_with_divergences": int(divergent), + } + + +def _coefficients(per_coefficient): + """Return the coefficient indices present, in order.""" + return sorted(per_coefficient) + + +def _per_coefficient(per_coefficient): + """Reduce a per-coefficient mapping of samples to a list of means.""" + return [float(np.mean(per_coefficient[i])) for i in _coefficients(per_coefficient)] + + +def _interval(row): + """Pull the credible interval out of a summary row, across ArviZ versions.""" + # ArviZ 0.x names them hdi_2.5%/hdi_97.5%; 1.x names them hdi95_lb/hdi95_ub. + lower = [c for c in row.index if c.startswith("hdi") and ("lb" in c or c.endswith("2.5%"))] + upper = [c for c in row.index if c.startswith("hdi") and ("ub" in c or c.endswith("97.5%"))] + return float(row[lower[0]]), float(row[upper[0]]) + + +def main(): + """Run the grid, write the results, and optionally enforce the thresholds.""" + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--replications", type=int, default=100) + parser.add_argument("--seed", type=int, default=20260818) + parser.add_argument( + "--jobs", + type=int, + default=max(1, (os.cpu_count() or 2) // 2), + help="cells to fit concurrently; each fit already uses one core per chain", + ) + parser.add_argument( + "--out", + default=op.join( + op.dirname(op.abspath(__file__)), + "..", + "..", + "pymare", + "tests", + "data", + "stan_validation.json", + ), + ) + parser.add_argument( + "--check", + action="store_true", + help="exit non-zero if any cell misses STAN_VALIDATION_THRESHOLDS", + ) + args = parser.parse_args() + + # Compile before forking. CmdStanPy builds in place, so parallel make + # invocations on the same source collide: on a cold cache one worker + # reliably fails before any cell runs. One compile up front makes every + # worker a cache hit. + print("compiling the model", flush=True) + StanMetaRegression().compile() + + started = time.time() + if args.jobs == 1: + results = [] + for index, cell in enumerate(CELLS): + print(f"[{index + 1}/{len(CELLS)}] {cell['name']}", flush=True) + results.append(run_cell(cell, args.replications, args.seed + index)) + else: + # Cells are independent and each seeded from its own index, so running + # them concurrently changes nothing about the numbers -- only the + # wall-clock. Each fit already uses 2 cores for its 2 chains. + print(f"running {len(CELLS)} cells across {args.jobs} processes", flush=True) + with ProcessPoolExecutor(max_workers=args.jobs) as pool: + futures = [ + pool.submit(run_cell, cell, args.replications, args.seed + index) + for index, cell in enumerate(CELLS) + ] + results = [future.result() for future in futures] + + payload = { + "replications": args.replications, + "seed": args.seed, + "elapsed_seconds": round(time.time() - started, 1), + "thresholds": STAN_VALIDATION_THRESHOLDS, + "cells": results, + } + with open(args.out, "w") as fobj: + json.dump(payload, fobj, indent=2) + + print(f"\n{'cell':<32}{'worst bias':>11}{'worst cov':>10}{'tau2 bias':>11}{'diverg.':>9}") + for cell in results: + print( + f"{cell['name']:<32}{cell['beta_bias']:>11.4f}" + f"{cell['beta_coverage']:>10.3f}{cell['tau2_bias']:>11.4f}" + f"{cell['fits_with_divergences']:>9d}" + ) + print(f"\nwrote {op.normpath(args.out)} in {payload['elapsed_seconds']}s") + + if not args.check: + return 0 + + # A short run can clear the coverage floor by luck: at 10 replications the + # estimate moves in steps of 0.05 with a standard error of 0.07. + if args.replications < MIN_REPLICATIONS_TO_CHECK: + print( + f"\n--check needs at least {MIN_REPLICATIONS_TO_CHECK} replications to mean " + f"anything; got {args.replications}." + ) + return 1 + + # The same thresholds the pinned results are held to, applied to this run. + floor = STAN_VALIDATION_THRESHOLDS["min_coverage"] + ceiling = STAN_VALIDATION_THRESHOLDS["max_beta_bias"] + failures = [ + f"{cell['name']}: coverage {cell['beta_coverage']:.3f} below {floor:.2f}" + for cell in results + if cell["beta_coverage"] < floor + ] + [ + f"{cell['name']}: |beta bias| {abs(cell['beta_bias']):.4f} above {ceiling}" + for cell in results + if abs(cell["beta_bias"]) > ceiling + ] + if failures: + print("\nTHRESHOLDS NOT MET:") + for failure in failures: + print(f" {failure}") + return 1 + + print(f"all {len(results)} cells meet the thresholds") + return 0 + + +if __name__ == "__main__": + sys.exit(main())