From 87ac3fdf8b29642a061b7a41ce9ff2c380e4d519 Mon Sep 17 00:00:00 2001 From: James Kent Date: Wed, 19 Aug 2026 00:20:23 -0500 Subject: [PATCH 1/7] wip: refactor/update STAN --- .github/workflows/testing.yml | 53 +- .gitignore | 7 + CONTRIBUTING.md | 27 +- MANIFEST.in | 1 + Makefile | 14 +- docs/installation.rst | 26 +- .../plot_meta-analysis_walkthrough.py | 8 +- pymare/estimators/estimators.py | 494 ++++++++++++++---- pymare/estimators/stan/meta_regression.stan | 47 ++ pymare/results.py | 132 ++++- pymare/tests/conftest.py | 85 ++- pymare/tests/test_stan_estimators.py | 484 ++++++++++++++++- pymare/tests/utils.py | 37 ++ pyproject.toml | 2 +- setup.cfg | 21 +- validation/stan/README.md | 149 ++++++ validation/stan/results.json | 273 ++++++++++ validation/stan/simulate.py | 220 ++++++++ 18 files changed, 1936 insertions(+), 144 deletions(-) create mode 100644 pymare/estimators/stan/meta_regression.stan create mode 100644 validation/stan/README.md create mode 100644 validation/stan/results.json create mode 100644 validation/stan/simulate.py diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index aee89d9..996cd53 100644 --- a/.github/workflows/testing.yml +++ b/.github/workflows/testing.yml @@ -30,10 +30,17 @@ 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. They get the job below to + # themselves. The tests that only check how PyMARE's inputs are translated + # into Stan's data block are deliberately unmarked, so they run here on every + # platform without needing CmdStan. PYTEST_UNIT_MARKERS: not stan PYTEST_STAN_MARKERS: stan + # Pinned rather than "latest": an unpinned version makes the cache key move on + # its own and makes a red run ambiguous between a PyMARE change and a Stan + # release. Same reasoning as the pinned R image in validation/robumeta. + CMDSTAN_VERSION: "2.36.0" jobs: # Determine if tests should be run based on commit message. @@ -104,6 +111,13 @@ jobs: needs: check_skip if: ${{ needs.check_skip.outputs.skip == 'false' }} runs-on: ubuntu-latest + env: + # Read by the pytest_collection_modifyitems hook in + # pymare/tests/conftest.py, which turns a missing or broken CmdStan into a + # failure rather than a skip. The previous gate probed for a module name + # that never existed, so this job passed for years while running none of + # the tests it exists to run, and nothing in a green log said so. + PYMARE_REQUIRE_CMDSTAN: "1" defaults: run: shell: bash @@ -117,6 +131,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..904fd31 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -98,8 +98,8 @@ 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 lint` | flake8 over `pymare` and `benchmarks` | nothing extra | @@ -107,6 +107,29 @@ environment does not have: 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, and that is exactly how the +Stan job passed for years while running none of the tests it existed to run -- +its gate probed for a module name that PyStan 3 never provided. The Stan job now +sets `PYMARE_REQUIRE_CMDSTAN=1`, and the `pytest_collection_modifyitems` hook in +`pymare/tests/conftest.py` fails the run outright, at collection, wherever that +is set and CmdStan is missing. + +Only the tests that actually sample are marked `stan`. The ones 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, in `validation/stan/`, which +reports bias and credible-interval coverage across a grid of designs. That is +not run in CI; see its README for what it measured and how to rerun it. + ### 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..6b0f83f 100644 --- a/Makefile +++ b/Makefile @@ -1,4 +1,4 @@ -.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 test_stan unittest # --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,8 +9,9 @@ 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 " benchmark to run the asv suite once in the current environment" @@ -22,6 +23,13 @@ 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) 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 698ad76..992d57c 100644 --- a/pymare/estimators/estimators.py +++ b/pymare/estimators/estimators.py @@ -1,6 +1,7 @@ """Meta-regression estimator classes.""" -import sys +import os +import os.path as op from abc import ABCMeta, abstractmethod from inspect import getfullargspec, signature from warnings import warn @@ -43,18 +44,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 @@ -639,11 +674,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 @@ -1435,79 +1471,328 @@ def _reml_nll(self, theta, y, n, X): return ll_ + 0.5 * np.log(np.linalg.det(F)) +#: 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 hashable label per observation. 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 that ``fit`` carries no downstream conditionals and the + translation can be tested without a CmdStan installation -- which is what + the estimator's own tests could not do while the translation lived inside + ``fit`` next to a call to the sampler. + + Two of those decisions are corrections rather than conveniences. ``sigma`` + is ``sqrt(v)``, because Stan's ``normal`` is parameterized by 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``; arbitrary labels, including strings and + non-consecutive integers, are therefore 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 " + "shape {}.".format(y.shape) + ) + y = y.reshape(-1) + n_observations = y.shape[0] + + 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]}." + ) + 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]}." + ) + + 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` is given a half-normal prior whose scale is taken from the data. + Stan's prior choice recommendations [2]_ suggest a half-normal(0, 1) or + half-t(4, 0, 1) when the number of groups is small enough that the data say + little about the group-level variance, on data scaled to unit variance. + PyMARE cannot rescale a caller's data, so the scale is derived from it + instead, which makes the prior equivariant: a fixed scale would be crushingly + informative on data measured in thousands and vacuous on data measured in + thousandths. + + The default is ``max(std(y), sqrt(mean(v)))`` rather than either term alone. + :math:`\tau` is the standard deviation of the group means, so it cannot + plausibly exceed the spread of the estimates themselves; and it should not be + presumed smaller than a typical standard error. Taking the larger of the two + means the prior never asserts that :math:`\tau` is small when either quantity + says otherwise. That asymmetry is what matters: ``validation/stan`` measures + credible-interval coverage falling to 0.83 when the scale is too small, while + a scale that is too large costs only precision in :math:`\tau^2` and leaves + coverage at nominal. Using ``sqrt(mean(v))`` alone, which was the first + default tried, undercovers whenever the between-group spread is much larger + than the sampling error. Using ``std(y)`` alone is zero when every estimate + coincides, which is not a usable scale. Pass ``tau_prior_scale`` explicitly to + override it, including to make it diffuse. :math:`\beta` keeps Stan's + implicit improper uniform prior, so with a diffuse prior on :math:`\tau` the + posterior means agree with + :obj:`~pymare.estimators.VarianceBasedLikelihoodEstimator` at ``method="ML"``. + + A QR reparameterization of ``X`` was considered and not adopted. It improves + the geometry when predictors are strongly correlated and, under a flat prior + on :math:`\beta`, leaves the posterior unchanged, but it costs a matrix + inverse and a back-transform and is incompatible with the ``normal_id_glm`` + form the model uses. PyMARE designs typically carry one to three predictors, + where the conditioning it addresses is rare. + + The Stan program ships as a source file and is compiled on first use, with + the executable cached beside it so that later processes reuse it. Shipping + a precompiled binary instead would require building CmdStan at wheel-build + time and publishing one wheel per platform, which is not a reasonable trade + for one optional estimator in an otherwise pure-Python package. + + 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.5 + + - 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. It was previously the + between-group standard deviation, because the parameter was passed to + Stan's ``normal`` where a scale is expected. + - Sampling variances are now converted to standard deviations before + being passed to Stan. They previously were not, so the model treated + ``v`` as ``sqrt(v)``. + - ``groups`` accepts any hashable labels and no longer has to be + integers in ``1..k``. + - :meth:`fit_dataset` now passes ``dataset.g`` as ``groups``. It + previously dropped it silently. + - ``ci`` now sets the width of the reported credible interval. It was + previously accepted and ignored. """ - _result_cls = BayesianMetaRegressionResults + _dataset_attr_map = {"groups": "g"} - def __init__(self, **sampling_kwargs): + _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() + ) + + "." + ) + + 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: - import stan - except ImportError: - raise ImportError("Please install pystan.") + self.model = cmdstanpy.CmdStanModel(stan_file=STAN_MODEL_PATH, force_compile=force) + except (PermissionError, OSError): + fallback_dir = op.join(op.expanduser("~"), ".pymare", "stan") + os.makedirs(fallback_dir, exist_ok=True) + warn( + f"Could not compile the Stan model beside {STAN_MODEL_PATH}, most likely " + f"because it is not writable. Compiling into {fallback_dir} instead.", + stacklevel=2, + ) + self.model = cmdstanpy.CmdStanModel( + stan_file=STAN_MODEL_PATH, + exe_file=op.join(fallback_dir, "meta_regression"), + force_compile=force, + ) - self.model = stan.build(spec, data=self.data) + return self def fit(self, y, v, X, groups=None): """Run the Stan sampler and return results. @@ -1522,17 +1807,26 @@ 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 hashable label per observation, identifying the groups of + observations in the y/v/X inputs. Labels may be of any hashable + type and need not be consecutive; they are encoded internally in + order of first occurrence by + :func:`~pymare.stats.encode_groups`. 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 ----- @@ -1540,41 +1834,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..a0c49de --- /dev/null +++ b/pymare/estimators/stan/meta_regression.stan @@ -0,0 +1,47 @@ +// 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 +// +// This is the Stan User's Guide random-effects meta-analysis model (Measurement +// Error and Meta-Analysis, section "Meta-Analysis") with the guide's stated +// extension to trial-specific predictors: the per-observation effects are given +// a regression on X. sigma_i is the *known* sampling standard deviation of +// observation i, i.e. sqrt(v_i) -- Stan's normal() is parameterized by a scale, +// not a variance. tau is the between-group standard deviation, and tau2 = tau^2 +// the between-group variance that 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 by tau avoids the + // funnel geometry that theta ~ normal(0, tau) produces when tau is near zero. + // That geometry is the dominant source of divergences in small-K hierarchical + // models, which is exactly this estimator's use case. + 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 the GLM form has + // hand-derived gradients and is documented as the faster of the two. The + // vector-alpha/vector-sigma overload takes a per-observation intercept + // (the group effect) and a per-observation 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..29dd0c7 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,79 @@ 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 split the library into ``arviz_base``/``arviz_stats``/``arviz_plots`` + and renamed the interval arguments: ``hdi_prob`` became ``ci_prob``, paired + with a ``ci_kind`` that defaults to an equal-tailed rather than a + highest-density interval. Requesting ``ci_kind="hdi"`` keeps the reported + interval the same kind across both versions, which is what the ``ci`` + argument has always meant here. + + ``round_to="none"`` is not cosmetic. ArviZ 1.x defaults to ``"auto"``, which + formats the summary for display by converting the floats to strings; a + caller doing arithmetic on the returned DataFrame would silently get + concatenation instead. 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 +1010,21 @@ 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 explicitly. ArviZ 1.x removed the automatic dispatch that used + # to let summary() accept a sampler fit directly, so a fit stored raw + # here would fail at every call site rather than at this one. Import + # lazily: cmdstanpy is an optional dependency, and a caller who passes + # an already-converted object should not need it installed. + 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 +1038,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 +1069,43 @@ 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 same ones :meth:`summary` reports, + rather than to everything the sampler recorded. A fitted model has one + ``theta`` per group, so plotting everything means one panel per group: + illegible at any size, and a hard error under ArviZ 1.x, which caps a + figure at ``rcParams["plot.max_subplots"]`` panels. Pass + ``include_theta=True`` for the group-level means, and raise that + rcParam if there are many groups. """ 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..8069e24 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,37 @@ 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 that a contributor without + it does not see red. In CI that leniency is the wrong default: a skip is + indistinguishable from a pass in a job log, and that is precisely how the + Stan job reported success for years while running none of the tests it + existed to run -- its gate probed ``find_spec("pystan")``, but PyStan 3 + installs a module named ``stan``, so the condition was unsatisfiable. + + Setting ``PYMARE_REQUIRE_CMDSTAN=1``, as the Stan CI job does, asserts that + the environment is supposed to be able to run them. Failing here, at + collection, reports that once and unmissably 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 +384,55 @@ def two_samp_data(): "sd2": np.sqrt(np.array([4, 16])), "n2": np.array([12, 16]), } + + +# ----------------------------------------------------------------------------- +# Stan estimator +# ----------------------------------------------------------------------------- + + +@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 rather than arbitrary. The ``variables`` + fixture has ``v`` near 1 throughout, where ``sqrt(v)`` and ``v`` are within a + few percent of each other -- so a model that passes variances where standard + deviations belong fits it about as well as the correct one. That is how the + original defect survived. Here ``v`` spans 0.01 to 0.16 while ``sqrt(v)`` + spans 0.1 to 0.4, a factor of 2.5 to 10 in a consistent direction, so the + mistake shows up as a badly inflated tau2. Do not reuse ``variables`` for + this. + + ``tau2`` and ``tau`` are likewise kept well apart (0.25 against 0.5) so that + reporting the standard deviation under the name of the variance fails a + tight interval rather than landing inside it. + """ + 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/test_stan_estimators.py b/pymare/tests/test_stan_estimators.py index 796d05c..6b6a97c 100644 --- a/pymare/tests/test_stan_estimators.py +++ b/pymare/tests/test_stan_estimators.py @@ -1,41 +1,483 @@ """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 os.path as op +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 _build_stan_data +from pymare.results import BayesianMetaRegressionResults +from pymare.tests.utils import cmdstan_is_available -pytestmark = pytest.mark.stan +requires_cmdstan = pytest.mark.skipif( + not cmdstan_is_available(), + reason="requires cmdstanpy and a CmdStan installation", +) + + +# ----------------------------------------------------------------------------- +# 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))) -requires_pystan = pytest.mark.skipif( - find_spec("pystan") is None, reason="requires the optional pystan dependency" + 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"] -@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 + +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))) + + 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. + + _dataset_attr_map was empty, so the argument fell back to its None default + and every fit_dataset() call silently modelled each observation as its own + group. Presetting .model with a stub keeps this test free of CmdStan: 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_compile_falls_back_when_the_package_directory_is_read_only(monkeypatch, tmp_path): + """An unwritable site-packages must not make the estimator unusable. + + CmdStanPy compiles beside the .stan source, which is inside the installed + package. That directory is read-only in plenty of ordinary installations, + and the resulting error would otherwise surface from the middle of fit(). + """ + cmdstanpy = pytest.importorskip("cmdstanpy") + monkeypatch.setenv("HOME", str(tmp_path)) + attempts = [] + + def fake_model(stan_file=None, exe_file=None, force_compile=False): + attempts.append(exe_file) + if exe_file is None: + raise PermissionError("read-only file system") + return "compiled" + + # Stub the CmdStan lookup as well as the compiler, so this exercises the + # fallback itself rather than requiring a real CmdStan to get that far. + monkeypatch.setattr(cmdstanpy, "cmdstan_path", lambda: str(tmp_path)) + monkeypatch.setattr(cmdstanpy, "CmdStanModel", fake_model) + + est = StanMetaRegression() + with pytest.warns(UserWarning, match="not writable"): + est.compile() + + assert est.model == "compiled" + assert attempts[0] is None + assert attempts[1] == op.join(str(tmp_path), ".pymare", "stan", "meta_regression") + + +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("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") + from pymare.results import _arviz_credible_interval_kwargs + + 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" + + +# ----------------------------------------------------------------------------- +# 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 class docstring has always promised this, but it could not be done: the + old compile() read self.data, which only fit() assigned, so calling it + directly raised AttributeError and every fit recompiled. Under CmdStanPy the + executable does not depend on the data, so the promise is now keepable. + """ + 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 diff --git a/pymare/tests/utils.py b/pymare/tests/utils.py index 57f477a..9455417 100644 --- a/pymare/tests/utils.py +++ b/pymare/tests/utils.py @@ -42,3 +42,40 @@ 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 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 22cc7fe..d6679d7 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 wrapt packages = find: @@ -68,8 +66,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 @@ -81,6 +85,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..5f373a5 --- /dev/null +++ b/validation/stan/README.md @@ -0,0 +1,149 @@ +# 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. + +Not run in CI. The fast tests in `pymare/tests/test_stan_estimators.py` check one +planted configuration; this checks the grid. + +## 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 + +## Reproducing + +```bash +pip install -e .[stan] +python -m cmdstanpy.install_cmdstan +python validation/stan/simulate.py --replications 100 +``` + +About 10 minutes on 8 cores. Results are written to `results.json`, which is +committed, so a change to the model can be diffed against it. + +## 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)`. + +Coverage is the fraction of 95% credible intervals for `beta` containing the +planted value; nominal is 0.950 and the Monte Carlo standard error is about +0.015 to 0.030. The two coverage columns are the two candidate defaults for +`tau_prior_scale` (see below). + +| cell | coverage, `sqrt(mean(v))` | coverage, `max(std(y), sqrt(mean(v)))` | tau2 bias, old | tau2 bias, new | true tau2 | beta bias | fits with divergences | +| --- | --- | --- | --- | --- | --- | --- | --- | +| groups=5 | 0.910 | **0.925** | -0.002 | +0.114 | 0.10 | -0.0021 | 63 | +| groups=20 | 0.940 | **0.950** | +0.004 | +0.017 | 0.10 | -0.0029 | 0 | +| groups=50 | 0.960 | **0.955** | +0.006 | +0.010 | 0.10 | -0.0022 | 0 | +| tau2=0 | 0.970 | **0.965** | +0.004 | +0.004 | 0.00 | +0.0002 | 11 | +| tau2=0.1 | 0.930 | **0.940** | +0.007 | +0.020 | 0.10 | -0.0065 | 0 | +| tau2=1 | 0.880 | **0.935** | -0.334 | +0.085 | 1.00 | +0.0151 | 0 | +| singletons | 0.955 | **0.960** | -0.006 | +0.016 | 0.10 | +0.0068 | 3 | +| unequal groups | 0.930 | **0.935** | +0.009 | +0.023 | 0.10 | -0.0001 | 0 | +| sigma x0.1 | 0.810 | **0.925** | -0.070 | +0.016 | 0.10 | +0.0014 | 0 | +| sigma x10 | 0.935 | **0.955** | +0.379 | +0.386 | 0.10 | -0.0188 | 9 | +| 1 predictor | 0.950 | **0.950** | +0.002 | +0.010 | 0.10 | +0.0089 | 0 | +| 3 predictors | 0.943 | **0.940** | +0.004 | +0.018 | 0.10 | +0.0000 | 0 | +| unbalanced covariate | 0.965 | **0.965** | +0.008 | +0.020 | 0.10 | +0.0080 | 0 | +| unbalanced covariate, tau2=1 | 0.870 | **0.945** | -0.300 | +0.153 | 1.00 | -0.0081 | 0 | + +`beta` is unbiased throughout: the largest bias in any cell is 0.019, 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. Coverage fell to **0.810** when the sampling +SDs were small relative to the between-group spread (`sigma x0.1`), to 0.880 at +`tau2=1`, and to 0.870 in the unbalanced-covariate cell at `tau2=1`. In each, +`tau2` was badly *under*-estimated: -70%, -33% and -30% 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/results.json b/validation/stan/results.json new file mode 100644 index 0000000..b537e23 --- /dev/null +++ b/validation/stan/results.json @@ -0,0 +1,273 @@ +{ + "replications": 100, + "seed": 20260818, + "elapsed_seconds": 583.2, + "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, + "beta_bias": -0.0021201448821405065, + "beta_rmse": 0.12639884294523676, + "beta_coverage": 0.925, + "coverage_se": 0.01862458053218917, + "tau2_bias": 0.11401999999999997, + "tau2_truth": 0.09999999999999998, + "fits_with_divergences": 63 + }, + { + "name": "groups=20", + "config": { + "n_groups": "20", + "group_size": "3", + "tau2": "0.1", + "sigma_scale": "1.0", + "n_predictors": "2", + "unbalanced": "False" + }, + "replications": 100, + "beta_bias": -0.002885654421913514, + "beta_rmse": 0.05071385238690838, + "beta_coverage": 0.95, + "coverage_se": 0.015411035007422448, + "tau2_bias": 0.017379999999999993, + "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, + "beta_bias": -0.002236553636109281, + "beta_rmse": 0.03725240459012178, + "beta_coverage": 0.955, + "coverage_se": 0.014658615214269054, + "tau2_bias": 0.010019999999999996, + "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, + "beta_bias": 0.000234083355221673, + "beta_rmse": 0.02752582425327561, + "beta_coverage": 0.965, + "coverage_se": 0.012995191418367032, + "tau2_bias": 0.004270000000000001, + "tau2_truth": 0.0, + "fits_with_divergences": 11 + }, + { + "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, + "beta_bias": -0.00645483896524104, + "beta_rmse": 0.0624116445444163, + "beta_coverage": 0.94, + "coverage_se": 0.01679285562374667, + "tau2_bias": 0.020149999999999998, + "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, + "beta_bias": 0.015128182286879798, + "beta_rmse": 0.1552137746779895, + "beta_coverage": 0.935, + "coverage_se": 0.01743201078476031, + "tau2_bias": 0.08488, + "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, + "beta_bias": 0.006795324568292994, + "beta_rmse": 0.08982016398564083, + "beta_coverage": 0.96, + "coverage_se": 0.013856406460551024, + "tau2_bias": 0.015779999999999995, + "tau2_truth": 0.09999999999999998, + "fits_with_divergences": 3 + }, + { + "name": "unequal groups", + "config": { + "n_groups": "20", + "group_size": "unequal", + "tau2": "0.1", + "sigma_scale": "1.0", + "n_predictors": "2", + "unbalanced": "False" + }, + "replications": 100, + "beta_bias": -8.534368402361576e-05, + "beta_rmse": 0.06195725104514044, + "beta_coverage": 0.935, + "coverage_se": 0.01743201078476031, + "tau2_bias": 0.022759999999999992, + "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, + "beta_bias": 0.0013507767768359108, + "beta_rmse": 0.047671429829435605, + "beta_coverage": 0.925, + "coverage_se": 0.01862458053218917, + "tau2_bias": 0.015859999999999996, + "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, + "beta_bias": -0.018845679375840186, + "beta_rmse": 0.29437805323971644, + "beta_coverage": 0.955, + "coverage_se": 0.014658615214269054, + "tau2_bias": 0.38577000000000006, + "tau2_truth": 0.09999999999999998, + "fits_with_divergences": 9 + }, + { + "name": "1 predictor", + "config": { + "n_groups": "20", + "group_size": "3", + "tau2": "0.1", + "sigma_scale": "1.0", + "n_predictors": "1", + "unbalanced": "False" + }, + "replications": 100, + "beta_bias": 0.008909029793187874, + "beta_rmse": 0.07746874870350999, + "beta_coverage": 0.95, + "coverage_se": 0.021794494717703377, + "tau2_bias": 0.009759999999999993, + "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, + "beta_bias": 6.009360521037991e-06, + "beta_rmse": 0.052885373073410256, + "beta_coverage": 0.94, + "coverage_se": 0.013711309200802093, + "tau2_bias": 0.018249999999999995, + "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, + "beta_bias": 0.008010131782226839, + "beta_rmse": 0.15405818433913557, + "beta_coverage": 0.965, + "coverage_se": 0.012995191418367032, + "tau2_bias": 0.02037999999999999, + "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, + "beta_bias": -0.00812323831412477, + "beta_rmse": 0.4654110785668132, + "beta_coverage": 0.945, + "coverage_se": 0.016120638945153514, + "tau2_bias": 0.15338000000000002, + "tau2_truth": 1.0, + "fits_with_divergences": 0 + } + ] +} \ No newline at end of file diff --git a/validation/stan/simulate.py b/validation/stan/simulate.py new file mode 100644 index 0000000..3223db3 --- /dev/null +++ b/validation/stan/simulate.py @@ -0,0 +1,220 @@ +"""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:: + + python validation/stan/simulate.py --replications 100 --out results.json + +The measured output is recorded in this directory's README.md. +""" + +import argparse +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 + +# 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 varies a single factor away from the base +#: configuration, which is the arrangement that attributes a failure to a +#: factor; a full factorial over five factors would cost 5x the fits and still +#: need this reading to interpret. +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}, +] + +BASE = { + "n_groups": 20, + "group_size": 3, + "tau2": 0.1, + "sigma_scale": 1.0, + "n_predictors": 2, + "unbalanced": False, +} + + +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 = rng.normal(size=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) + + beta_errors, covered, 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.append(float(row["mean"]) - true_value) + lower, upper = _interval(row) + covered.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, + "beta_bias": float(np.mean(beta_errors)), + "beta_rmse": float(np.sqrt(np.mean(np.square(beta_errors)))), + "beta_coverage": float(np.mean(covered)), + "coverage_se": float(np.sqrt(np.mean(covered) * (1 - np.mean(covered)) / len(covered))), + "tau2_bias": float(np.mean(tau2_errors)), + "tau2_truth": float(np.mean(tau2_truth)), + "fits_with_divergences": int(divergent), + } + + +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 and write the results.""" + 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__)), "results.json")) + args = parser.parse_args() + + 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), + "cells": results, + } + with open(args.out, "w") as fobj: + json.dump(payload, fobj, indent=2) + + print(f"\n{'cell':<32}{'beta bias':>11}{'coverage':>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 {args.out} in {payload['elapsed_seconds']}s") + + +if __name__ == "__main__": + main() From 2be61efab256e977b0505abda2e6ed1bca4f33ff Mon Sep 17 00:00:00 2001 From: James Kent Date: Wed, 19 Aug 2026 00:54:30 -0500 Subject: [PATCH 2/7] validate the results of stan --- .github/workflows/stan-validation.yml | 115 ++++++++++++++++++ CONTRIBUTING.md | 22 +++- Makefile | 11 +- .../tests/data/stan_validation.json | 6 +- pymare/tests/test_stan_estimators.py | 58 ++++++++- pymare/tests/utils.py | 70 +++++++++++ validation/stan/README.md | 40 +++++- validation/stan/simulate.py | 79 +++++++++++- 8 files changed, 384 insertions(+), 17 deletions(-) create mode 100644 .github/workflows/stan-validation.yml rename validation/stan/results.json => pymare/tests/data/stan_validation.json (98%) 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/CONTRIBUTING.md b/CONTRIBUTING.md index 904fd31..93a550b 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -102,6 +102,7 @@ environment does not have: | `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 @@ -126,9 +127,24 @@ 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, in `validation/stan/`, which -reports bias and credible-interval coverage across a grid of designs. That is -not run in CI; see its README for what it measured and how to rerun it. +The model's accuracy is measured separately, by `validation/stan/simulate.py`, +which reports bias and credible-interval coverage across a grid of designs and +records them in `pymare/tests/data/stan_validation.json`. It follows the same +three-layer arrangement as the robumeta alignment: + +1. `make validate_stan` regenerates that file and fails if any design cell + misses `pymare.tests.utils.STAN_VALIDATION_THRESHOLDS`. +2. Two tests in `test_stan_estimators.py` hold the *recorded* file to those same + thresholds and to the expected list of design cells. They read the file + rather than re-measuring, so they cost nothing and run everywhere — which is + what stops the pin from quietly going stale. +3. The `Validate the Stan model` workflow re-measures on a schedule, on pushes + to master that touch the model, and on demand. + +The grid takes about ten minutes, which is why it is not part of `test_stan` and +not run per pull request. Unlike the robumeta reference, these numbers are +stochastic, so the pin cannot be enforced by requiring the file not to move; the +thresholds are the claim, and the file is the record of it. ### Alignment with robumeta diff --git a/Makefile b/Makefile index 6b0f83f..0146e03 100644 --- a/Makefile +++ b/Makefile @@ -1,4 +1,5 @@ -.PHONY: all_tests benchmark check_robumeta_alignment help install_cmdstan 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. @@ -14,6 +15,7 @@ help: @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" @@ -36,6 +38,13 @@ test_stan: 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/validation/stan/results.json b/pymare/tests/data/stan_validation.json similarity index 98% rename from validation/stan/results.json rename to pymare/tests/data/stan_validation.json index b537e23..86915b0 100644 --- a/validation/stan/results.json +++ b/pymare/tests/data/stan_validation.json @@ -1,7 +1,11 @@ { "replications": 100, "seed": 20260818, - "elapsed_seconds": 583.2, + "elapsed_seconds": 402.9, + "thresholds": { + "min_coverage": 0.9, + "max_beta_bias": 0.1 + }, "cells": [ { "name": "groups=5", diff --git a/pymare/tests/test_stan_estimators.py b/pymare/tests/test_stan_estimators.py index 6b6a97c..ab7acad 100644 --- a/pymare/tests/test_stan_estimators.py +++ b/pymare/tests/test_stan_estimators.py @@ -17,7 +17,12 @@ from pymare.estimators import StanMetaRegression, VarianceBasedLikelihoodEstimator from pymare.estimators.estimators import _build_stan_data from pymare.results import BayesianMetaRegressionResults -from pymare.tests.utils import cmdstan_is_available +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(), @@ -359,6 +364,57 @@ def test_summary_requests_the_configured_credible_interval(ci): 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 + checks 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 enforces the same thresholds against a fresh run. + + The thresholds are not decoration either: the first prior scale tried here + produced coverage of 0.810 in the ``sigma x0.1`` cell, which this floor + rejects. + """ + 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}" + + # ----------------------------------------------------------------------------- # Sampling. Needs CmdStan. # ----------------------------------------------------------------------------- diff --git a/pymare/tests/utils.py b/pymare/tests/utils.py index 9455417..5e7cd2f 100644 --- a/pymare/tests/utils.py +++ b/pymare/tests/utils.py @@ -79,3 +79,73 @@ def cmdstan_is_available(): return False return True + + +#: Thresholds the Stan model's simulated performance has to meet, checked in two +#: places against the same numbers: ``test_stan_estimators.py`` asserts them +#: against the pinned ``data/stan_validation.json``, and +#: ``validation/stan/simulate.py --check`` asserts them against a fresh run. +#: +#: ``min_coverage`` is the floor for the fraction of 95% credible intervals for +#: beta that contain the planted value. It sits below the nominal 0.95 because +#: the estimate is a Monte Carlo one: at 100 replications its standard error is +#: about 0.02, so a correct estimator lands below nominal routinely. 0.90 is +#: roughly nominal minus two standard errors, which the first prior tried here +#: violated outright at 0.810 -- so the floor is loose enough not to fire on +#: noise and tight enough to have caught the defect it was written for. +#: +#: ``max_beta_bias`` is generous for the same reason: beta is unbiased in theory, +#: and the largest bias measured across the grid is 0.019 against coefficients of +#: order 1. +STAN_VALIDATION_THRESHOLDS = { + "min_coverage": 0.90, + "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/validation/stan/README.md b/validation/stan/README.md index 5f373a5..bbf6d74 100644 --- a/validation/stan/README.md +++ b/validation/stan/README.md @@ -6,8 +6,9 @@ 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. -Not run in CI. The fast tests in `pymare/tests/test_stan_estimators.py` check one -planted configuration; this checks the grid. +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 @@ -31,16 +32,45 @@ 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: +coverage at or above 0.90 and |beta bias| at or below 0.10 in every cell. That +floor is roughly nominal minus two standard errors — loose enough not to fire on +noise, and tight enough that it rejects the 0.810 the first prior scale produced. +`--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 -python validation/stan/simulate.py --replications 100 +make validate_stan # or: python validation/stan/simulate.py --check ``` -About 10 minutes on 8 cores. Results are written to `results.json`, which is -committed, so a change to the model can be diffed against it. +About 10 minutes on 8 cores. ## What was measured diff --git a/validation/stan/simulate.py b/validation/stan/simulate.py index 3223db3..26d8783 100644 --- a/validation/stan/simulate.py +++ b/validation/stan/simulate.py @@ -8,9 +8,20 @@ Run with:: - python validation/stan/simulate.py --replications 100 --out results.json + make validate_stan -The measured output is recorded in this directory's README.md. +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 @@ -29,6 +40,7 @@ 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. @@ -65,6 +77,10 @@ {"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, @@ -166,7 +182,7 @@ def _interval(row): def main(): - """Run the grid and write the results.""" + """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) @@ -176,7 +192,23 @@ def main(): 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__)), "results.json")) + 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() started = time.time() @@ -201,6 +233,7 @@ def main(): "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: @@ -213,8 +246,42 @@ def main(): f"{cell['beta_coverage']:>10.3f}{cell['tau2_bias']:>11.4f}" f"{cell['fits_with_divergences']:>9d}" ) - print(f"\nwrote {args.out} in {payload['elapsed_seconds']}s") + print(f"\nwrote {op.normpath(args.out)} in {payload['elapsed_seconds']}s") + + if not args.check: + return 0 + + # A handful of replications can clear the coverage floor by luck -- at 10 + # replications the estimate moves in steps of 0.05 and its standard error is + # 0.07 -- so a short run must not be able to certify the model. + 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__": - main() + sys.exit(main()) From 5ff27ce9afd21a144e3298dccccfed7db32d1c73 Mon Sep 17 00:00:00 2001 From: James Kent Date: Wed, 19 Aug 2026 02:05:44 -0500 Subject: [PATCH 3/7] [TST] cover the ArviZ-only and constraint code paths Codecov reported 82.47% patch coverage against a 92.18% target. The cause was structural rather than a few missed lines: .codecov.yml ignores pymare/tests/, so only source counts, and BayesianMetaRegressionResults is ArviZ-only code that no unit job could execute because only the Stan job installed ArviZ. The whole results container was reachable from one job, on one Python, on one platform. The unit job now installs the stan extra as well. cmdstanpy comes with it but stays idle -- it is pure Python, CmdStan is not installed there, and the tests that sample are excluded by the marker filter regardless. The effect is that the ArviZ 0.x versus 1.x handling is now exercised across the whole matrix rather than resting on a single job. Also adds the tests Interval and Options never had. Interval had grown `closed` and `allow_none` for tau_prior_scale with nothing checking either. Co-Authored-By: Claude Opus 5 (1M context) --- .github/workflows/testing.yml | 10 +++++- pymare/tests/test_estimators.py | 58 ++++++++++++++++++++++++++++++++- 2 files changed, 66 insertions(+), 2 deletions(-) diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index 996cd53..0510026 100644 --- a/.github/workflows/testing.yml +++ b/.github/workflows/testing.yml @@ -93,10 +93,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" diff --git a/pymare/tests/test_estimators.py b/pymare/tests/test_estimators.py index 0218190..7ab9857 100644 --- a/pymare/tests/test_estimators.py +++ b/pymare/tests/test_estimators.py @@ -11,7 +11,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, collapse_groups_by_n @@ -750,3 +750,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") From 040c22ae87b6530605c95c596360809118d3c9f0 Mon Sep 17 00:00:00 2001 From: James Kent Date: Wed, 19 Aug 2026 02:06:06 -0500 Subject: [PATCH 4/7] [FIX] address review findings on the Stan migration Six review comments, each verified against the code before fixing. Two were worse than reported. The read-only fallback in compile() never worked. CmdStanPy reports every failed make invocation as ValueError, so catching PermissionError/OSError could not fire in the one situation the fallback exists for. Independently, exe_file= names an executable to reuse rather than a destination to build into, so the fallback would have failed even had it been reached: make writes its intermediates beside the source. It now compiles a copy in ~/.pymare/stan, using copy2 so the preserved mtime keeps the cached build across processes, and re-raises the original error when that also fails -- a model that does not parse should not be reported as a permissions problem. Verified end to end against a real chmod 555 directory. The test that covered it had asserted PermissionError because that is what a read-only filesystem sounds like. It passed while the code under it could not run, which is the same shape of defect as the skip gate this branch started from, so the test now pins the exception CmdStanPy actually raises. The validation harness could not detect what it was written to detect. It redrew the true coefficients from a symmetric normal on every replication, so the signed errors averaged to zero for any estimator at all: one that always returned zero cleared the bias ceiling 84.6% of the time. It also pooled coverage across coefficients, which let a well-estimated intercept mask a badly estimated moderator -- exactly the failure the unbalanced-covariate cells exist to probe. The truth is now fixed, coverage is reported per coefficient, and the thresholds apply to the worst one. Under the sharper metric the prior scale this branch already rejected reads 0.710 rather than 0.810, so the pooling was understating it. The coverage floor moved to 0.85, chosen by measuring the rejected prior under the new metric rather than by judgement: it reads 0.710 and 0.830 in two cells while the current model's tightest honest cell reads 0.900. A minimum over coefficients is biased downward, so a floor nearer nominal would flake. Parallel workers raced to compile the same model. With a cold cache and four workers one of them reliably failed with "Failed to compile Stan model" before any cell ran, which is what a fresh validation runner would have hit. The model is now compiled once in the parent before the pool starts. NaN sampling variances passed the positivity check, since NaN fails every comparison, and surfaced later as a CmdStan data-loading error naming a Stan variable rather than the input responsible. y, v and X are now all checked for finiteness at the boundary. The groups docstring promised any hashable label, but numpy reads a sequence of tuples as a second dimension, so composite labels are rejected by encode_groups. The contract is narrowed to scalar labels rather than widening shared code that other estimators depend on. The sixth comment, that no CI job ran the unmarked results tests, was already resolved by the preceding commit. Co-Authored-By: Claude Opus 5 (1M context) --- pymare/estimators/estimators.py | 74 ++++-- pymare/tests/data/stan_validation.json | 304 +++++++++++++++++++------ pymare/tests/test_stan_estimators.py | 279 +++++++++++++++++++++-- pymare/tests/utils.py | 34 ++- validation/stan/README.md | 75 +++--- validation/stan/simulate.py | 60 ++++- 6 files changed, 664 insertions(+), 162 deletions(-) diff --git a/pymare/estimators/estimators.py b/pymare/estimators/estimators.py index 992d57c..8ebf99e 100644 --- a/pymare/estimators/estimators.py +++ b/pymare/estimators/estimators.py @@ -2,6 +2,7 @@ import os import os.path as op +import shutil from abc import ABCMeta, abstractmethod from inspect import getfullargspec, signature from warnings import warn @@ -1543,8 +1544,11 @@ def _build_stan_data(y, v, X, groups=None, tau_prior_scale=None): 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 hashable label per observation. When None (default), every - observation is its own group. + 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 @@ -1585,8 +1589,10 @@ def _build_stan_data(y, v, X, groups=None, tau_prior_scale=None): "not support 2-dimensional inputs. Passed y has " "shape {}.".format(y.shape) ) - y = y.reshape(-1) + 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: @@ -1594,6 +1600,11 @@ def _build_stan_data(y, v, X, groups=None, tau_prior_scale=None): 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.") @@ -1605,6 +1616,8 @@ def _build_stan_data(y, v, X, groups=None, tau_prior_scale=None): 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) @@ -1778,20 +1791,37 @@ def compile(self, force=False): try: self.model = cmdstanpy.CmdStanModel(stan_file=STAN_MODEL_PATH, force_compile=force) - except (PermissionError, OSError): - fallback_dir = op.join(op.expanduser("~"), ".pymare", "stan") + return self + except Exception as unwritable: + # Deliberately broad. CmdStanPy reports *any* failed make invocation + # as ValueError, including for the read-only package directory this + # fallback exists for, so catching OSError would never fire. + first_failure = unwritable + + # Compile a copy instead. Passing exe_file= would not work: that names an + # executable to reuse, not a destination to build into, so a read-only + # source directory fails there too -- make writes its intermediates + # beside the source. copy2 preserves the modification time, so the copy + # is not perpetually 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: os.makedirs(fallback_dir, exist_ok=True) - warn( - f"Could not compile the Stan model beside {STAN_MODEL_PATH}, most likely " - f"because it is not writable. Compiling into {fallback_dir} instead.", - stacklevel=2, - ) - self.model = cmdstanpy.CmdStanModel( - stan_file=STAN_MODEL_PATH, - exe_file=op.join(fallback_dir, "meta_regression"), - force_compile=force, - ) + 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: + # Compiling somewhere writable failed too, so the first failure was + # not about writing. Report that one: it names the real problem, + # usually an error in the model or the C++ toolchain. + raise first_failure + 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): @@ -1808,12 +1838,14 @@ def fit(self, y, v, X, groups=None): (including intercept); has dimensions K x P, where K is the number of observations and P is the number of predictor variables. groups : None or array-like of shape (K,), optional - One hashable label per observation, identifying the groups of - observations in the y/v/X inputs. Labels may be of any hashable - type and need not be consecutive; they are encoded internally in - order of first occurrence by - :func:`~pymare.stats.encode_groups`. When None (default), each - observation in the inputs is treated as a separate group. + 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 ------- diff --git a/pymare/tests/data/stan_validation.json b/pymare/tests/data/stan_validation.json index 86915b0..50ae722 100644 --- a/pymare/tests/data/stan_validation.json +++ b/pymare/tests/data/stan_validation.json @@ -1,7 +1,7 @@ { "replications": 100, "seed": 20260818, - "elapsed_seconds": 402.9, + "elapsed_seconds": 335.2, "thresholds": { "min_coverage": 0.9, "max_beta_bias": 0.1 @@ -18,13 +18,24 @@ "unbalanced": "False" }, "replications": 100, - "beta_bias": -0.0021201448821405065, - "beta_rmse": 0.12639884294523676, - "beta_coverage": 0.925, - "coverage_se": 0.01862458053218917, - "tau2_bias": 0.11401999999999997, + "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": 63 + "fits_with_divergences": 54 }, { "name": "groups=20", @@ -37,11 +48,22 @@ "unbalanced": "False" }, "replications": 100, - "beta_bias": -0.002885654421913514, - "beta_rmse": 0.05071385238690838, - "beta_coverage": 0.95, - "coverage_se": 0.015411035007422448, - "tau2_bias": 0.017379999999999993, + "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 }, @@ -56,11 +78,22 @@ "unbalanced": "False" }, "replications": 100, - "beta_bias": -0.002236553636109281, - "beta_rmse": 0.03725240459012178, - "beta_coverage": 0.955, - "coverage_se": 0.014658615214269054, - "tau2_bias": 0.010019999999999996, + "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 }, @@ -75,13 +108,24 @@ "unbalanced": "False" }, "replications": 100, - "beta_bias": 0.000234083355221673, - "beta_rmse": 0.02752582425327561, - "beta_coverage": 0.965, - "coverage_se": 0.012995191418367032, - "tau2_bias": 0.004270000000000001, + "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": 11 + "fits_with_divergences": 14 }, { "name": "tau2=0.1", @@ -94,11 +138,22 @@ "unbalanced": "False" }, "replications": 100, - "beta_bias": -0.00645483896524104, - "beta_rmse": 0.0624116445444163, - "beta_coverage": 0.94, - "coverage_se": 0.01679285562374667, - "tau2_bias": 0.020149999999999998, + "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 }, @@ -113,11 +168,22 @@ "unbalanced": "False" }, "replications": 100, - "beta_bias": 0.015128182286879798, - "beta_rmse": 0.1552137746779895, - "beta_coverage": 0.935, - "coverage_se": 0.01743201078476031, - "tau2_bias": 0.08488, + "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 }, @@ -132,13 +198,24 @@ "unbalanced": "False" }, "replications": 100, - "beta_bias": 0.006795324568292994, - "beta_rmse": 0.08982016398564083, - "beta_coverage": 0.96, - "coverage_se": 0.013856406460551024, - "tau2_bias": 0.015779999999999995, + "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": 3 + "fits_with_divergences": 1 }, { "name": "unequal groups", @@ -151,11 +228,22 @@ "unbalanced": "False" }, "replications": 100, - "beta_bias": -8.534368402361576e-05, - "beta_rmse": 0.06195725104514044, - "beta_coverage": 0.935, - "coverage_se": 0.01743201078476031, - "tau2_bias": 0.022759999999999992, + "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 }, @@ -170,11 +258,22 @@ "unbalanced": "False" }, "replications": 100, - "beta_bias": 0.0013507767768359108, - "beta_rmse": 0.047671429829435605, - "beta_coverage": 0.925, - "coverage_se": 0.01862458053218917, - "tau2_bias": 0.015859999999999996, + "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 }, @@ -189,13 +288,24 @@ "unbalanced": "False" }, "replications": 100, - "beta_bias": -0.018845679375840186, - "beta_rmse": 0.29437805323971644, - "beta_coverage": 0.955, - "coverage_se": 0.014658615214269054, - "tau2_bias": 0.38577000000000006, + "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": 9 + "fits_with_divergences": 6 }, { "name": "1 predictor", @@ -208,11 +318,19 @@ "unbalanced": "False" }, "replications": 100, - "beta_bias": 0.008909029793187874, - "beta_rmse": 0.07746874870350999, - "beta_coverage": 0.95, - "coverage_se": 0.021794494717703377, - "tau2_bias": 0.009759999999999993, + "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 }, @@ -227,11 +345,25 @@ "unbalanced": "False" }, "replications": 100, - "beta_bias": 6.009360521037991e-06, - "beta_rmse": 0.052885373073410256, - "beta_coverage": 0.94, - "coverage_se": 0.013711309200802093, - "tau2_bias": 0.018249999999999995, + "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 }, @@ -246,11 +378,22 @@ "unbalanced": "True" }, "replications": 100, - "beta_bias": 0.008010131782226839, - "beta_rmse": 0.15405818433913557, - "beta_coverage": 0.965, - "coverage_se": 0.012995191418367032, - "tau2_bias": 0.02037999999999999, + "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 }, @@ -265,11 +408,22 @@ "unbalanced": "True" }, "replications": 100, - "beta_bias": -0.00812323831412477, - "beta_rmse": 0.4654110785668132, - "beta_coverage": 0.945, - "coverage_se": 0.016120638945153514, - "tau2_bias": 0.15338000000000002, + "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 } diff --git a/pymare/tests/test_stan_estimators.py b/pymare/tests/test_stan_estimators.py index ab7acad..07b5c42 100644 --- a/pymare/tests/test_stan_estimators.py +++ b/pymare/tests/test_stan_estimators.py @@ -8,6 +8,7 @@ """ import os.path as op +import sys import warnings import numpy as np @@ -15,8 +16,17 @@ from pymare import meta_regression from pymare.estimators import StanMetaRegression, VarianceBasedLikelihoodEstimator -from pymare.estimators.estimators import _build_stan_data -from pymare.results import BayesianMetaRegressionResults +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, @@ -30,6 +40,78 @@ ) +# ----------------------------------------------------------------------------- +# 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 + + +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) + + # 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. # ----------------------------------------------------------------------------- @@ -195,22 +277,26 @@ def test_fit_is_quiet_when_there_are_no_divergences(): def test_compile_falls_back_when_the_package_directory_is_read_only(monkeypatch, tmp_path): """An unwritable site-packages must not make the estimator unusable. - CmdStanPy compiles beside the .stan source, which is inside the installed - package. That directory is read-only in plenty of ordinary installations, - and the resulting error would otherwise surface from the middle of fit(). + CmdStanPy compiles beside the .stan source, which lives inside the installed + package, and that directory is read-only in plenty of ordinary + installations. + + The failure is raised as ValueError, not PermissionError: CmdStanPy reports + every failed ``make`` the same way, whatever went wrong. An earlier version + of this test asserted PermissionError because that is what a read-only + filesystem sounds like, and it passed while the fallback it was meant to + cover could never fire. """ cmdstanpy = pytest.importorskip("cmdstanpy") monkeypatch.setenv("HOME", str(tmp_path)) - attempts = [] + compiled_from = [] def fake_model(stan_file=None, exe_file=None, force_compile=False): - attempts.append(exe_file) - if exe_file is None: - raise PermissionError("read-only file system") + compiled_from.append(stan_file) + if stan_file == STAN_MODEL_PATH: + raise ValueError(f"Failed to compile Stan model '{stan_file}'.") return "compiled" - # Stub the CmdStan lookup as well as the compiler, so this exercises the - # fallback itself rather than requiring a real CmdStan to get that far. monkeypatch.setattr(cmdstanpy, "cmdstan_path", lambda: str(tmp_path)) monkeypatch.setattr(cmdstanpy, "CmdStanModel", fake_model) @@ -219,8 +305,37 @@ def fake_model(stan_file=None, exe_file=None, force_compile=False): est.compile() assert est.model == "compiled" - assert attempts[0] is None - assert attempts[1] == op.join(str(tmp_path), ".pymare", "stan", "meta_regression") + # 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(tmp_path), ".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, tmp_path): + """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") + monkeypatch.setenv("HOME", str(tmp_path)) + + 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(tmp_path)) + 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: @@ -347,12 +462,61 @@ def test_results_reject_an_impossible_credible_interval(bad_ci): 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") - from pymare.results import _arviz_credible_interval_kwargs - kwargs = _arviz_credible_interval_kwargs(ci) probability = kwargs.get("ci_prob", kwargs.get("hdi_prob")) @@ -415,6 +579,45 @@ def test_recorded_validation_meets_its_thresholds(): 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. + + An earlier version of the harness redrew beta from a symmetric normal on + every replication. The signed errors then averaged to zero for *any* + estimator -- one that always returned zero cleared the bias ceiling about + 85% of the time -- so the threshold certified 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. # ----------------------------------------------------------------------------- @@ -537,3 +740,49 @@ def test_summary_and_plot_round_trip(planted_hierarchical_dataset): 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 5e7cd2f..83fd8be 100644 --- a/pymare/tests/utils.py +++ b/pymare/tests/utils.py @@ -86,19 +86,31 @@ def cmdstan_is_available(): #: against the pinned ``data/stan_validation.json``, and #: ``validation/stan/simulate.py --check`` asserts them against a fresh run. #: -#: ``min_coverage`` is the floor for the fraction of 95% credible intervals for -#: beta that contain the planted value. It sits below the nominal 0.95 because -#: the estimate is a Monte Carlo one: at 100 replications its standard error is -#: about 0.02, so a correct estimator lands below nominal routinely. 0.90 is -#: roughly nominal minus two standard errors, which the first prior tried here -#: violated outright at 0.810 -- so the floor is loose enough not to fire on -#: noise and tight enough to have caught the defect it was written for. +#: ``min_coverage`` is the floor for the *worst* coefficient's share of 95% +#: credible intervals containing the true value. Both halves of that sentence +#: were chosen by measurement rather than taste. #: -#: ``max_beta_bias`` is generous for the same reason: beta is unbiased in theory, -#: and the largest bias measured across the grid is 0.019 against coefficients of -#: order 1. +#: Reporting the worst coefficient rather than the average matters because +#: averaging lets a well-estimated intercept mask a badly estimated moderator: +#: the prior scale this model originally shipped with measures 0.810 pooled but +#: 0.710 on its worst coefficient, so pooling understated the defect. +#: +#: The floor is 0.85 rather than something nearer the nominal 0.95 because a +#: minimum over coefficients is biased downward -- each coefficient's coverage is +#: a Monte Carlo estimate with a standard error near 0.022 at 100 replications, +#: and the smallest of two or three such estimates sits below nominal routinely. +#: The correct model's tightest cell measures 0.900, and the rejected prior +#: measures 0.710 and 0.830 in two cells, so 0.85 sits between them: about two +#: standard errors below the worst honest result, and far enough above the +#: failures to catch a regression of that kind. +#: +#: ``max_beta_bias`` applies to the largest absolute bias across coefficients. +#: The largest measured anywhere in the grid is 0.045, against coefficients of +#: order 1. It is only meaningful because the harness holds the true +#: coefficients fixed; when it redrew them from a symmetric normal, an estimator +#: that always returned zero cleared this ceiling about 85% of the time. STAN_VALIDATION_THRESHOLDS = { - "min_coverage": 0.90, + "min_coverage": 0.85, "max_beta_bias": 0.10, } diff --git a/validation/stan/README.md b/validation/stan/README.md index bbf6d74..3eeed53 100644 --- a/validation/stan/README.md +++ b/validation/stan/README.md @@ -56,9 +56,12 @@ 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: -coverage at or above 0.90 and |beta bias| at or below 0.10 in every cell. That -floor is roughly nominal minus two standard errors — loose enough not to fire on -noise, and tight enough that it rejects the 0.810 the first prior scale produced. +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 and +the rejected prior reads 0.710 and 0.830, so 0.85 sits between them — about two +standard errors below the worst honest result and clear of the failures. 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. @@ -78,38 +81,48 @@ About 10 minutes on 8 cores. 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)`. -Coverage is the fraction of 95% credible intervals for `beta` containing the -planted value; nominal is 0.950 and the Monte Carlo standard error is about -0.015 to 0.030. The two coverage columns are the two candidate defaults for -`tau_prior_scale` (see below). - -| cell | coverage, `sqrt(mean(v))` | coverage, `max(std(y), sqrt(mean(v)))` | tau2 bias, old | tau2 bias, new | true tau2 | beta bias | fits with divergences | -| --- | --- | --- | --- | --- | --- | --- | --- | -| groups=5 | 0.910 | **0.925** | -0.002 | +0.114 | 0.10 | -0.0021 | 63 | -| groups=20 | 0.940 | **0.950** | +0.004 | +0.017 | 0.10 | -0.0029 | 0 | -| groups=50 | 0.960 | **0.955** | +0.006 | +0.010 | 0.10 | -0.0022 | 0 | -| tau2=0 | 0.970 | **0.965** | +0.004 | +0.004 | 0.00 | +0.0002 | 11 | -| tau2=0.1 | 0.930 | **0.940** | +0.007 | +0.020 | 0.10 | -0.0065 | 0 | -| tau2=1 | 0.880 | **0.935** | -0.334 | +0.085 | 1.00 | +0.0151 | 0 | -| singletons | 0.955 | **0.960** | -0.006 | +0.016 | 0.10 | +0.0068 | 3 | -| unequal groups | 0.930 | **0.935** | +0.009 | +0.023 | 0.10 | -0.0001 | 0 | -| sigma x0.1 | 0.810 | **0.925** | -0.070 | +0.016 | 0.10 | +0.0014 | 0 | -| sigma x10 | 0.935 | **0.955** | +0.379 | +0.386 | 0.10 | -0.0188 | 9 | -| 1 predictor | 0.950 | **0.950** | +0.002 | +0.010 | 0.10 | +0.0089 | 0 | -| 3 predictors | 0.943 | **0.940** | +0.004 | +0.018 | 0.10 | +0.0000 | 0 | -| unbalanced covariate | 0.965 | **0.965** | +0.008 | +0.020 | 0.10 | +0.0080 | 0 | -| unbalanced covariate, tau2=1 | 0.870 | **0.945** | -0.300 | +0.153 | 1.00 | -0.0081 | 0 | - -`beta` is unbiased throughout: the largest bias in any cell is 0.019, against -coefficients of order 1. +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. Coverage fell to **0.810** when the sampling -SDs were small relative to the between-group spread (`sigma x0.1`), to 0.880 at -`tau2=1`, and to 0.870 in the unbalanced-covariate cell at `tau2=1`. In each, -`tau2` was badly *under*-estimated: -70%, -33% and -30% respectively. +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 diff --git a/validation/stan/simulate.py b/validation/stan/simulate.py index 26d8783..394b3f2 100644 --- a/validation/stan/simulate.py +++ b/validation/stan/simulate.py @@ -25,6 +25,7 @@ """ import argparse +import collections import json import logging import os @@ -91,6 +92,14 @@ } +#: True coefficients, held fixed across replications rather than redrawn. +#: Redrawing them from a symmetric distribution makes the pooled bias +#: uninformative: the error of an estimator that always returned zero would be +#: -beta, whose mean over replications is zero, so it would clear any bias +#: threshold. Fixing the truth means a bias estimate measures the estimator. +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": @@ -114,7 +123,7 @@ def simulate(rng, n_groups, group_size, tau2, sigma_scale, n_predictors, unbalan if moderators.size else np.ones((n_observations, 1)) ) - beta = rng.normal(size=X.shape[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 @@ -131,7 +140,12 @@ def run_cell(cell, replications, seed): config.update({k: v for k, v in cell.items() if k != "name"}) rng = np.random.default_rng(seed) - beta_errors, covered, tau2_errors, tau2_truth, divergent = [], [], [], [], 0 + # Per coefficient, not pooled. Pooling hides the case this grid exists to + # probe: in the unbalanced cells the sparse moderator is the coefficient at + # risk, and good intercept coverage would mask bad coverage for 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) @@ -151,9 +165,9 @@ def run_cell(cell, replications, seed): for i, true_value in enumerate(beta): row = summary.loc[f"beta[{i}]"] - beta_errors.append(float(row["mean"]) - true_value) + beta_errors[i].append(float(row["mean"]) - true_value) lower, upper = _interval(row) - covered.append(bool(lower <= true_value <= upper)) + covered[i].append(bool(lower <= true_value <= upper)) tau2_errors.append(float(summary.loc["tau2", "mean"]) - tau2) tau2_truth.append(tau2) @@ -163,16 +177,36 @@ def run_cell(cell, replications, seed): "name": cell["name"], "config": {k: str(v) for k, v in config.items()}, "replications": replications, - "beta_bias": float(np.mean(beta_errors)), - "beta_rmse": float(np.sqrt(np.mean(np.square(beta_errors)))), - "beta_coverage": float(np.mean(covered)), - "coverage_se": float(np.sqrt(np.mean(covered) * (1 - np.mean(covered)) / len(covered))), + "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. @@ -211,6 +245,14 @@ def main(): ) args = parser.parse_args() + # Compile before any worker starts. CmdStanPy builds in place, and parallel + # make invocations on the same source collide: with a cold cache and four + # workers, one of them reliably fails with "Failed to compile Stan model" + # before a single 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 = [] @@ -239,7 +281,7 @@ def main(): with open(args.out, "w") as fobj: json.dump(payload, fobj, indent=2) - print(f"\n{'cell':<32}{'beta bias':>11}{'coverage':>10}{'tau2 bias':>11}{'diverg.':>9}") + 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}" From 8447ed574fb372162b6426c0783c0158bfe060de Mon Sep 17 00:00:00 2001 From: James Kent Date: Wed, 19 Aug 2026 06:40:34 -0500 Subject: [PATCH 5/7] [FIX] make the compile-fallback tests work on Windows The Windows unit job failed on test_compile_falls_back_when_the_package_directory_is_read_only. The estimator is fine: expanduser("~") resolves the user profile on Windows, which is what the fallback wants. The test was wrong. It redirected the home directory by setting HOME. POSIX expanduser reads HOME, but Windows reads USERPROFILE and ignores HOME entirely, falling back to HOMEDRIVE/HOMEPATH and then to leaving "~" unexpanded. So the redirect silently did nothing there: the assertion compared a temporary path against the runner's real profile. Worse, the sibling test that drives both compiles to failure was creating .pymare/stan and copying the model into the runner's actual home directory, because the same redirect was equally ineffective. Both now use a fake_home fixture that sets HOME and USERPROFILE together. A test asserts the redirect holds under ntpath as well as posixpath, so this Windows-only failure mode is caught on every platform -- reverting the fixture to the HOME-only version fails on Linux. Co-Authored-By: Claude Opus 5 (1M context) --- pymare/tests/conftest.py | 24 ++++++++++++++++++++++++ pymare/tests/test_stan_estimators.py | 28 +++++++++++++++++++++------- 2 files changed, 45 insertions(+), 7 deletions(-) diff --git a/pymare/tests/conftest.py b/pymare/tests/conftest.py index 8069e24..73a6c05 100644 --- a/pymare/tests/conftest.py +++ b/pymare/tests/conftest.py @@ -391,6 +391,30 @@ def two_samp_data(): # ----------------------------------------------------------------------------- +@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, falling back to + ``HOMEDRIVE``/``HOMEPATH`` and then to leaving ``~`` unexpanded. Setting + only ``HOME`` therefore looks like it works everywhere while silently + leaving the real profile in place on Windows -- which is what let a test of + the Stan compile fallback assert against a temporary path on Linux and macOS + while writing into the CI runner's actual home directory on Windows. + """ + 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. diff --git a/pymare/tests/test_stan_estimators.py b/pymare/tests/test_stan_estimators.py index 07b5c42..e76af92 100644 --- a/pymare/tests/test_stan_estimators.py +++ b/pymare/tests/test_stan_estimators.py @@ -7,7 +7,9 @@ where the defects this file now pins would have been caught years earlier. """ +import ntpath import os.path as op +import posixpath import sys import warnings @@ -274,7 +276,21 @@ def test_fit_is_quiet_when_there_are_no_divergences(): est.fit(np.arange(4.0), np.ones(4), np.ones((4, 1))) -def test_compile_falls_back_when_the_package_directory_is_read_only(monkeypatch, tmp_path): +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 Linux and macOS the + mistake that only shows up on a Windows runner: setting HOME alone leaves + ntpath.expanduser pointing at the real user profile, so a test asserting a + temporary path passes on two platforms and fails on the third -- after + having written into the runner's actual home directory. + """ + 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 @@ -288,7 +304,6 @@ def test_compile_falls_back_when_the_package_directory_is_read_only(monkeypatch, cover could never fire. """ cmdstanpy = pytest.importorskip("cmdstanpy") - monkeypatch.setenv("HOME", str(tmp_path)) compiled_from = [] def fake_model(stan_file=None, exe_file=None, force_compile=False): @@ -297,7 +312,7 @@ def fake_model(stan_file=None, exe_file=None, force_compile=False): raise ValueError(f"Failed to compile Stan model '{stan_file}'.") return "compiled" - monkeypatch.setattr(cmdstanpy, "cmdstan_path", lambda: str(tmp_path)) + monkeypatch.setattr(cmdstanpy, "cmdstan_path", lambda: str(fake_home)) monkeypatch.setattr(cmdstanpy, "CmdStanModel", fake_model) est = StanMetaRegression() @@ -308,26 +323,25 @@ def fake_model(stan_file=None, exe_file=None, force_compile=False): # 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(tmp_path), ".pymare", "stan", "meta_regression.stan") + 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, tmp_path): +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") - monkeypatch.setenv("HOME", str(tmp_path)) 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(tmp_path)) + monkeypatch.setattr(cmdstanpy, "cmdstan_path", lambda: str(fake_home)) monkeypatch.setattr(cmdstanpy, "CmdStanModel", always_fails) est = StanMetaRegression() From 41c0b96b950db2342a56002a8aa228b0e1800180 Mon Sep 17 00:00:00 2001 From: James Kent Date: Wed, 19 Aug 2026 09:14:48 -0500 Subject: [PATCH 6/7] [DOC] tighten the docstrings and comments added by this branch A focus pass over the prose this branch introduced. Removes 88 lines without losing a fact, and fixes two that were wrong. Two factual errors. The class docstring said credible-interval coverage fell to 0.83 under the rejected prior scale; that was the pooled figure, and the per-coefficient measurement the branch now uses reads 0.710. The versionchanged note still advertised "any hashable labels" for groups after the contract had been narrowed to scalar labels, which is the opposite of what the code does. One explanation, one place. The three-layer validation arrangement was written out in full in both CONTRIBUTING.md and validation/stan/README.md, and the 0.810-versus-0.710 measurement appeared in three files. The README is the canonical record; the others now state the invariant and point at it. Prose duplicated across files goes stale in the copies nobody edits. Cut process narration. Several docstrings described what an earlier version of the same code or test had done -- "an earlier version of this test asserted PermissionError", "which was the first default tried". That belongs in the commit history, not in a docstring a reader meets years later. The invariant is what survives. Cut two rejected alternatives that were never implemented: a QR reparameterization and precompiled platform wheels. Recording what was removed and why is worth the lines; speculating about roads not taken is not, and the packaging rationale was already stated where the model path is resolved. Verified: 387 tests pass, the Stan model recompiles from the edited source and still samples, the docs build clean for both touched modules, and the rendered page keeps its math, references and versionchanged note. Co-Authored-By: Claude Opus 5 (1M context) --- CONTRIBUTING.md | 48 +++----- pymare/estimators/estimators.py | 115 ++++++++------------ pymare/estimators/stan/meta_regression.stan | 28 +++-- pymare/results.py | 41 +++---- pymare/tests/conftest.py | 52 ++++----- pymare/tests/test_stan_estimators.py | 47 +++----- pymare/tests/utils.py | 37 ++----- validation/stan/README.md | 11 +- validation/stan/simulate.py | 35 +++--- 9 files changed, 163 insertions(+), 251 deletions(-) diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 93a550b..ff6361f 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -114,37 +114,25 @@ 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, and that is exactly how the -Stan job passed for years while running none of the tests it existed to run -- -its gate probed for a module name that PyStan 3 never provided. The Stan job now +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 outright, at collection, wherever that -is set and CmdStan is missing. - -Only the tests that actually sample are marked `stan`. The ones 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 `validation/stan/simulate.py`, -which reports bias and credible-interval coverage across a grid of designs and -records them in `pymare/tests/data/stan_validation.json`. It follows the same -three-layer arrangement as the robumeta alignment: - -1. `make validate_stan` regenerates that file and fails if any design cell - misses `pymare.tests.utils.STAN_VALIDATION_THRESHOLDS`. -2. Two tests in `test_stan_estimators.py` hold the *recorded* file to those same - thresholds and to the expected list of design cells. They read the file - rather than re-measuring, so they cost nothing and run everywhere — which is - what stops the pin from quietly going stale. -3. The `Validate the Stan model` workflow re-measures on a schedule, on pushes - to master that touch the model, and on demand. - -The grid takes about ten minutes, which is why it is not part of `test_stan` and -not run per pull request. Unlike the robumeta reference, these numbers are -stochastic, so the pin cannot be enforced by requiring the file not to move; the -thresholds are the claim, and the file is the record of it. +`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 diff --git a/pymare/estimators/estimators.py b/pymare/estimators/estimators.py index 8ebf99e..9789570 100644 --- a/pymare/estimators/estimators.py +++ b/pymare/estimators/estimators.py @@ -1570,17 +1570,15 @@ def _build_stan_data(y, v, X, groups=None, tau_prior_scale=None): Notes ----- Every shape and unit decision the Stan program depends on is made here and - nowhere else, so that ``fit`` carries no downstream conditionals and the - translation can be tested without a CmdStan installation -- which is what - the estimator's own tests could not do while the translation lived inside - ``fit`` next to a call to the sampler. - - Two of those decisions are corrections rather than conveniences. ``sigma`` - is ``sqrt(v)``, because Stan's ``normal`` is parameterized by 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``; arbitrary labels, including strings and - non-consecutive integers, are therefore accepted here. + nowhere else, so ``fit`` carries no downstream conditionals and the + translation is testable without a CmdStan installation. + + Two of those decisions are conversions, not conveniences. ``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: @@ -1672,44 +1670,25 @@ class StanMetaRegression(BaseEstimator): the funnel geometry that dominates divergences in hierarchical models with few groups, which is this estimator's principal use case. - :math:`\tau` is given a half-normal prior whose scale is taken from the data. - Stan's prior choice recommendations [2]_ suggest a half-normal(0, 1) or - half-t(4, 0, 1) when the number of groups is small enough that the data say - little about the group-level variance, on data scaled to unit variance. - PyMARE cannot rescale a caller's data, so the scale is derived from it - instead, which makes the prior equivariant: a fixed scale would be crushingly - informative on data measured in thousands and vacuous on data measured in - thousandths. - - The default is ``max(std(y), sqrt(mean(v)))`` rather than either term alone. - :math:`\tau` is the standard deviation of the group means, so it cannot - plausibly exceed the spread of the estimates themselves; and it should not be - presumed smaller than a typical standard error. Taking the larger of the two - means the prior never asserts that :math:`\tau` is small when either quantity - says otherwise. That asymmetry is what matters: ``validation/stan`` measures - credible-interval coverage falling to 0.83 when the scale is too small, while - a scale that is too large costs only precision in :math:`\tau^2` and leaves - coverage at nominal. Using ``sqrt(mean(v))`` alone, which was the first - default tried, undercovers whenever the between-group spread is much larger - than the sampling error. Using ``std(y)`` alone is zero when every estimate - coincides, which is not a usable scale. Pass ``tau_prior_scale`` explicitly to - override it, including to make it diffuse. :math:`\beta` keeps Stan's - implicit improper uniform prior, so with a diffuse prior on :math:`\tau` the - posterior means agree with + :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"``. - A QR reparameterization of ``X`` was considered and not adopted. It improves - the geometry when predictors are strongly correlated and, under a flat prior - on :math:`\beta`, leaves the posterior unchanged, but it costs a matrix - inverse and a back-transform and is incompatible with the ``normal_id_glm`` - form the model uses. PyMARE designs typically carry one to three predictors, - where the conditioning it addresses is rare. - - The Stan program ships as a source file and is compiled on first use, with - the executable cached beside it so that later processes reuse it. Shipping - a precompiled binary instead would require building CmdStan at wheel-build - time and publishing one wheel per platform, which is not a reasonable trade - for one optional estimator in an otherwise pure-Python package. + 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 ---------- @@ -1721,19 +1700,16 @@ class StanMetaRegression(BaseEstimator): .. versionchanged:: 0.0.5 - - 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. It was previously the - between-group standard deviation, because the parameter was passed to - Stan's ``normal`` where a scale is expected. - - Sampling variances are now converted to standard deviations before - being passed to Stan. They previously were not, so the model treated - ``v`` as ``sqrt(v)``. - - ``groups`` accepts any hashable labels and no longer has to be - integers in ``1..k``. - - :meth:`fit_dataset` now passes ``dataset.g`` as ``groups``. It - previously dropped it silently. + - 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. """ @@ -1793,16 +1769,14 @@ def compile(self, force=False): 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, including for the read-only package directory this - # fallback exists for, so catching OSError would never fire. + # Deliberately broad: CmdStanPy reports any failed make invocation + # as ValueError, so catching OSError would never fire here. first_failure = unwritable - # Compile a copy instead. Passing exe_file= would not work: that names an - # executable to reuse, not a destination to build into, so a read-only - # source directory fails there too -- make writes its intermediates - # beside the source. copy2 preserves the modification time, so the copy - # is not perpetually newer than its own executable and CmdStanPy's + # 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: @@ -1811,9 +1785,8 @@ def compile(self, force=False): shutil.copy2(STAN_MODEL_PATH, fallback_source) model = cmdstanpy.CmdStanModel(stan_file=fallback_source, force_compile=force) except Exception: - # Compiling somewhere writable failed too, so the first failure was - # not about writing. Report that one: it names the real problem, - # usually an error in the model or the C++ toolchain. + # Somewhere writable failed too, so the first failure was not about + # writing. Report that one -- it names the real problem. raise first_failure warn( diff --git a/pymare/estimators/stan/meta_regression.stan b/pymare/estimators/stan/meta_regression.stan index a0c49de..f37b51d 100644 --- a/pymare/estimators/stan/meta_regression.stan +++ b/pymare/estimators/stan/meta_regression.stan @@ -3,13 +3,13 @@ // y_i ~ normal(x_i' beta + theta_{g(i)}, sigma_i) i = 1..N // theta_g ~ normal(0, tau) g = 1..K // -// This is the Stan User's Guide random-effects meta-analysis model (Measurement -// Error and Meta-Analysis, section "Meta-Analysis") with the guide's stated -// extension to trial-specific predictors: the per-observation effects are given -// a regression on X. sigma_i is the *known* sampling standard deviation of -// observation i, i.e. sqrt(v_i) -- Stan's normal() is parameterized by a scale, -// not a variance. tau is the between-group standard deviation, and tau2 = tau^2 -// the between-group variance that every other PyMARE estimator reports. +// 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) @@ -26,20 +26,18 @@ parameters { real tau; } transformed parameters { - // Non-centered: sampling theta_raw ~ N(0, 1) and scaling by tau avoids the - // funnel geometry that theta ~ normal(0, tau) produces when tau is near zero. - // That geometry is the dominant source of divergences in small-K hierarchical - // models, which is exactly this estimator's use case. + // 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 the GLM form has - // hand-derived gradients and is documented as the faster of the two. The - // vector-alpha/vector-sigma overload takes a per-observation intercept - // (the group effect) and a per-observation scale (the sampling SD). + // 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 { diff --git a/pymare/results.py b/pymare/results.py index 29dd0c7..632b02b 100644 --- a/pymare/results.py +++ b/pymare/results.py @@ -944,17 +944,15 @@ def _arviz_credible_interval_kwargs(ci): Notes ----- - ArviZ 1.0 split the library into ``arviz_base``/``arviz_stats``/``arviz_plots`` - and renamed the interval arguments: ``hdi_prob`` became ``ci_prob``, paired - with a ``ci_kind`` that defaults to an equal-tailed rather than a - highest-density interval. Requesting ``ci_kind="hdi"`` keeps the reported - interval the same kind across both versions, which is what the ``ci`` - argument has always meant here. - - ``round_to="none"`` is not cosmetic. ArviZ 1.x defaults to ``"auto"``, which - formats the summary for display by converting the floats to strings; a - caller doing arithmetic on the returned DataFrame would silently get - concatenation instead. ArviZ 0.x has no such argument. + 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"} @@ -1013,11 +1011,10 @@ def __init__(self, data, dataset, ci=95.0): if not 0 < ci < 100: raise ValueError(f"Invalid ci {ci!r}; must lie in (0, 100).") - # Convert explicitly. ArviZ 1.x removed the automatic dispatch that used - # to let summary() accept a sampler fit directly, so a fit stored raw - # here would fail at every call site rather than at this one. Import - # lazily: cmdstanpy is an optional dependency, and a caller who passes - # an already-converted object should not need it installed. + # 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: @@ -1088,13 +1085,11 @@ def plot(self, kind="trace", include_theta=False, **kwargs): Notes ----- - The plotted variables default to the same ones :meth:`summary` reports, - rather than to everything the sampler recorded. A fitted model has one - ``theta`` per group, so plotting everything means one panel per group: - illegible at any size, and a hard error under ArviZ 1.x, which caps a - figure at ``rcParams["plot.max_subplots"]`` panels. Pass - ``include_theta=True`` for the group-level means, and raise that - rcParam if there are many groups. + 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) # Three-argument getattr: the two-argument form raises AttributeError diff --git a/pymare/tests/conftest.py b/pymare/tests/conftest.py index 73a6c05..56c1c15 100644 --- a/pymare/tests/conftest.py +++ b/pymare/tests/conftest.py @@ -23,17 +23,14 @@ 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 that a contributor without - it does not see red. In CI that leniency is the wrong default: a skip is - indistinguishable from a pass in a job log, and that is precisely how the - Stan job reported success for years while running none of the tests it - existed to run -- its gate probed ``find_spec("pystan")``, but PyStan 3 - installs a module named ``stan``, so the condition was unsatisfiable. - - Setting ``PYMARE_REQUIRE_CMDSTAN=1``, as the Stan CI job does, asserts that - the environment is supposed to be able to run them. Failing here, at - collection, reports that once and unmissably rather than as a quietly - shorter run. + 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 @@ -402,13 +399,11 @@ def fake_home(tmp_path, monkeypatch): Notes ----- - Both variables are needed. POSIX ``expanduser`` reads ``HOME``; Windows - reads ``USERPROFILE`` and ignores ``HOME`` entirely, falling back to - ``HOMEDRIVE``/``HOMEPATH`` and then to leaving ``~`` unexpanded. Setting - only ``HOME`` therefore looks like it works everywhere while silently - leaving the real profile in place on Windows -- which is what let a test of - the Stan compile fallback assert against a temporary path on Linux and macOS - while writing into the CI runner's actual home directory on Windows. + 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)) @@ -428,18 +423,15 @@ def planted_hierarchical_dataset(): Notes ----- The sampling standard deviations are drawn from ``uniform(0.1, 0.4)``, well - away from 1, and that is load-bearing rather than arbitrary. The ``variables`` - fixture has ``v`` near 1 throughout, where ``sqrt(v)`` and ``v`` are within a - few percent of each other -- so a model that passes variances where standard - deviations belong fits it about as well as the correct one. That is how the - original defect survived. Here ``v`` spans 0.01 to 0.16 while ``sqrt(v)`` - spans 0.1 to 0.4, a factor of 2.5 to 10 in a consistent direction, so the - mistake shows up as a badly inflated tau2. Do not reuse ``variables`` for - this. - - ``tau2`` and ``tau`` are likewise kept well apart (0.25 against 0.5) so that - reporting the standard deviation under the name of the variance fails a - tight interval rather than landing inside it. + 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) diff --git a/pymare/tests/test_stan_estimators.py b/pymare/tests/test_stan_estimators.py index e76af92..c063f01 100644 --- a/pymare/tests/test_stan_estimators.py +++ b/pymare/tests/test_stan_estimators.py @@ -235,10 +235,9 @@ def test_stan_2d_input_failure(dataset_2d): def test_fit_dataset_forwards_dataset_g(planted_hierarchical_dataset): """fit_dataset must route dataset.g into fit()'s groups argument. - _dataset_attr_map was empty, so the argument fell back to its None default - and every fit_dataset() call silently modelled each observation as its own - group. Presetting .model with a stub keeps this test free of CmdStan: fit() - only compiles when it finds no model. + 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() @@ -279,11 +278,9 @@ def test_fit_is_quiet_when_there_are_no_divergences(): 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 Linux and macOS the - mistake that only shows up on a Windows runner: setting HOME alone leaves - ntpath.expanduser pointing at the real user profile, so a test asserting a - temporary path passes on two platforms and fails on the third -- after - having written into the runner's actual home directory. + 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) @@ -294,14 +291,9 @@ def test_compile_falls_back_when_the_package_directory_is_read_only(monkeypatch, """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 ordinary - installations. - - The failure is raised as ValueError, not PermissionError: CmdStanPy reports - every failed ``make`` the same way, whatever went wrong. An earlier version - of this test asserted PermissionError because that is what a read-only - filesystem sounds like, and it passed while the fallback it was meant to - cover could never fire. + 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 = [] @@ -566,13 +558,9 @@ 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 - checks the numbers already measured rather than re-measuring, so it costs + 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 enforces the same thresholds against a fresh run. - - The thresholds are not decoration either: the first prior scale tried here - produced coverage of 0.810 in the ``sigma x0.1`` cell, which this floor - rejects. + re-measures and applies the same thresholds to a fresh run. """ recorded = load_stan_validation() floor = STAN_VALIDATION_THRESHOLDS["min_coverage"] @@ -616,10 +604,9 @@ def test_recorded_validation_summarizes_the_worst_coefficient(): def test_recorded_validation_used_a_fixed_truth(): """Bias is only meaningful if the coefficients being recovered are held fixed. - An earlier version of the harness redrew beta from a symmetric normal on - every replication. The signed errors then averaged to zero for *any* - estimator -- one that always returned zero cleared the bias ceiling about - 85% of the time -- so the threshold certified nothing. + 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"]} @@ -718,10 +705,8 @@ def test_meta_regression_dispatches_to_stan(planted_hierarchical_dataset): def test_the_compiled_model_is_reused_across_fits(planted_hierarchical_dataset): """compile() once, fit many. - The class docstring has always promised this, but it could not be done: the - old compile() read self.data, which only fit() assigned, so calling it - directly raised AttributeError and every fit recompiled. Under CmdStanPy the - executable does not depend on the data, so the promise is now keepable. + 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) diff --git a/pymare/tests/utils.py b/pymare/tests/utils.py index 83fd8be..b38f1da 100644 --- a/pymare/tests/utils.py +++ b/pymare/tests/utils.py @@ -81,34 +81,21 @@ def cmdstan_is_available(): return True -#: Thresholds the Stan model's simulated performance has to meet, checked in two -#: places against the same numbers: ``test_stan_estimators.py`` asserts them +#: 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`` asserts them against a fresh run. +#: ``validation/stan/simulate.py --check`` checks them against a fresh run. #: -#: ``min_coverage`` is the floor for the *worst* coefficient's share of 95% -#: credible intervals containing the true value. Both halves of that sentence -#: were chosen by measurement rather than taste. +#: ``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. #: -#: Reporting the worst coefficient rather than the average matters because -#: averaging lets a well-estimated intercept mask a badly estimated moderator: -#: the prior scale this model originally shipped with measures 0.810 pooled but -#: 0.710 on its worst coefficient, so pooling understated the defect. -#: -#: The floor is 0.85 rather than something nearer the nominal 0.95 because a -#: minimum over coefficients is biased downward -- each coefficient's coverage is -#: a Monte Carlo estimate with a standard error near 0.022 at 100 replications, -#: and the smallest of two or three such estimates sits below nominal routinely. -#: The correct model's tightest cell measures 0.900, and the rejected prior -#: measures 0.710 and 0.830 in two cells, so 0.85 sits between them: about two -#: standard errors below the worst honest result, and far enough above the -#: failures to catch a regression of that kind. -#: -#: ``max_beta_bias`` applies to the largest absolute bias across coefficients. -#: The largest measured anywhere in the grid is 0.045, against coefficients of -#: order 1. It is only meaningful because the harness holds the true -#: coefficients fixed; when it redrew them from a symmetric normal, an estimator -#: that always returned zero cleared this ceiling about 85% of the time. +#: ``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, diff --git a/validation/stan/README.md b/validation/stan/README.md index 3eeed53..7c6b153 100644 --- a/validation/stan/README.md +++ b/validation/stan/README.md @@ -56,12 +56,11 @@ 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 and -the rejected prior reads 0.710 and 0.830, so 0.85 sits between them — about two -standard errors below the worst honest result and clear of the failures. A -minimum over coefficients is biased downward, which is why the floor is not -nearer the nominal 0.95. +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. diff --git a/validation/stan/simulate.py b/validation/stan/simulate.py index 394b3f2..91ac05f 100644 --- a/validation/stan/simulate.py +++ b/validation/stan/simulate.py @@ -48,10 +48,9 @@ cmdstanpy.disable_logging() logging.getLogger("cmdstanpy").setLevel(logging.ERROR) -#: One entry per design cell. Each varies a single factor away from the base -#: configuration, which is the arrangement that attributes a failure to a -#: factor; a full factorial over five factors would cost 5x the fits and still -#: need this reading to interpret. +#: 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}, @@ -92,11 +91,10 @@ } -#: True coefficients, held fixed across replications rather than redrawn. -#: Redrawing them from a symmetric distribution makes the pooled bias -#: uninformative: the error of an estimator that always returned zero would be -#: -beta, whose mean over replications is zero, so it would clear any bias -#: threshold. Fixing the truth means a bias estimate measures the estimator. +#: 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]) @@ -140,9 +138,8 @@ def run_cell(cell, replications, seed): config.update({k: v for k, v in cell.items() if k != "name"}) rng = np.random.default_rng(seed) - # Per coefficient, not pooled. Pooling hides the case this grid exists to - # probe: in the unbalanced cells the sparse moderator is the coefficient at - # risk, and good intercept coverage would mask bad coverage for it. + # 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 @@ -245,11 +242,10 @@ def main(): ) args = parser.parse_args() - # Compile before any worker starts. CmdStanPy builds in place, and parallel - # make invocations on the same source collide: with a cold cache and four - # workers, one of them reliably fails with "Failed to compile Stan model" - # before a single cell runs. One compile up front makes every worker a - # cache hit. + # 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() @@ -293,9 +289,8 @@ def main(): if not args.check: return 0 - # A handful of replications can clear the coverage floor by luck -- at 10 - # replications the estimate moves in steps of 0.05 and its standard error is - # 0.07 -- so a short run must not be able to certify the model. + # 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 " From 62cbccaa4226341e2a59f8ee8b08e360dd8275c4 Mon Sep 17 00:00:00 2001 From: James Kent Date: Wed, 19 Aug 2026 09:58:01 -0500 Subject: [PATCH 7/7] [DOC] trim workflow comments and correct the versionchanged target Working-tree edits from the IDE, plus the wrapping they needed to pass lint. Thins out the comments in the testing workflow, and corrects the versionchanged note on StanMetaRegression from 0.0.5 to 0.0.11: the latest tag is 0.0.10, so 0.0.11 is the next release rather than a version long past. The docstring reflow arrived under-indented by one space, which numpydoc reads as a nested block (D207); rewrapped at the surrounding indentation. Co-Authored-By: Claude Opus 5 (1M context) --- .github/workflows/testing.yml | 21 +-------------------- pymare/estimators/estimators.py | 9 ++++----- 2 files changed, 5 insertions(+), 25 deletions(-) diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index 0510026..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,19 +21,11 @@ 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 tests that sample, whose Stan model has to be compiled - # and so costs minutes rather than seconds. They get the job below to - # themselves. The tests that only check how PyMARE's inputs are translated - # into Stan's data block are deliberately unmarked, so they run here on every - # platform without needing CmdStan. + # and so costs minutes rather than seconds. PYTEST_UNIT_MARKERS: not stan PYTEST_STAN_MARKERS: stan - # Pinned rather than "latest": an unpinned version makes the cache key move on - # its own and makes a red run ambiguous between a PyMARE change and a Stan - # release. Same reasoning as the pinned R image in validation/robumeta. CMDSTAN_VERSION: "2.36.0" jobs: @@ -120,11 +106,6 @@ jobs: if: ${{ needs.check_skip.outputs.skip == 'false' }} runs-on: ubuntu-latest env: - # Read by the pytest_collection_modifyitems hook in - # pymare/tests/conftest.py, which turns a missing or broken CmdStan into a - # failure rather than a skip. The previous gate probed for a module name - # that never existed, so this job passed for years while running none of - # the tests it exists to run, and nothing in a green log said so. PYMARE_REQUIRE_CMDSTAN: "1" defaults: run: diff --git a/pymare/estimators/estimators.py b/pymare/estimators/estimators.py index 9789570..86de8dc 100644 --- a/pymare/estimators/estimators.py +++ b/pymare/estimators/estimators.py @@ -1573,9 +1573,8 @@ def _build_stan_data(y, v, X, groups=None, tau_prior_scale=None): nowhere else, so ``fit`` carries no downstream conditionals and the translation is testable without a CmdStan installation. - Two of those decisions are conversions, not conveniences. ``sigma`` is - ``sqrt(v)``, because Stan's ``normal`` takes a standard deviation and PyMARE - stores variances. ``id`` is 1-based consecutive codes from + ``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. @@ -1585,7 +1584,7 @@ def _build_stan_data(y, v, X, groups=None, tau_prior_scale=None): raise ValueError( "The StanMetaRegression estimator currently does " "not support 2-dimensional inputs. Passed y has " - "shape {}.".format(y.shape) + f"shape {y.shape}." ) y = np.asarray(y, dtype=float).reshape(-1) n_observations = y.shape[0] @@ -1698,7 +1697,7 @@ class StanMetaRegression(BaseEstimator): .. [2] Stan Development Team. Prior Choice Recommendations. https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations - .. versionchanged:: 0.0.5 + .. 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``)