align metafor and PyMARE

.github/workflows/alignment.yml

@@ -0,0 +1,106 @@
+name: "Check alignment with R packages"
+
+# metafor and clubSandwich, one matrix leg each. robumeta has its own workflow
+# rather than a third leg here: its check is named "Check robumeta alignment /
+# PyMARE vs robumeta" and may be a required check on the repository, which
+# folding it in would silently rename. The two share everything that matters --
+# validation/compare_reference.py and the regenerate.sh contract -- so the
+# duplication is a workflow header, not logic.
+
+on:
+  push:
+    branches:
+      - "master"
+  pull_request:
+    branches:
+      - "*"
+  schedule:
+    # Monthly, to catch the harness rotting rather than PyMARE moving: the R
+    # version and the package versions are pinned, so a failure here on an
+    # unchanged tree means a pinned image can no longer be built.
+    - cron: "0 0 1 * *"
+  # 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. Being able to dispatch it means re-enabling is enough to
+  # get a run without pushing a commit to prove it works.
+  workflow_dispatch:
+
+permissions:
+  contents: read
+
+concurrency:
+  group: alignment-${{ github.ref }}
+  cancel-in-progress: true
+
+jobs:
+  alignment:
+    name: PyMARE vs ${{ matrix.label }}
+    runs-on: ubuntu-latest
+    strategy:
+      # Independent references, so a failure in one should still report the
+      # other rather than being masked by it.
+      fail-fast: false
+      matrix:
+        include:
+          - package: metafor
+            label: metafor
+            marker: metafor
+            # rma.uni, escalc and permutest. validation/metafor/regenerate.sh
+            # rewrites all three; each is compared on its own so a failure
+            # names which reference moved.
+            references: >-
+              metafor_reference.json
+              metafor_escalc_reference.json
+              metafor_permutest_reference.json
+          - package: clubsandwich
+            label: clubSandwich
+            marker: clubsandwich
+            references: clubsandwich_reference.json
+    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"
+        run: |
+          python -m pip install --progress-bar off --upgrade pip setuptools wheel
+          python -m pip install -e .[tests]
+
+      # The pinned reference values are only worth anything if they still match
+      # what the R package prints, so regenerate them from the pinned image.
+      - name: "Regenerate the reference values from ${{ matrix.label }}"
+        run: validation/${{ matrix.package }}/regenerate.sh
+
+      # Compared numerically rather than with git diff: the values are written
+      # at full double precision, and which BLAS kernel R's image picks depends
+      # on the runner's CPU, so two runners disagree in the last bits with the R
+      # and package versions identical. See the tolerances in
+      # validation/compare_reference.py.
+      - name: "Require the reference values to still match ${{ matrix.label }}"
+        run: |
+          status=0
+          for reference in ${{ matrix.references }}; do
+              git show "HEAD:pymare/tests/data/${reference}" > "${RUNNER_TEMP}/${reference}"
+              python validation/compare_reference.py \
+                  "${RUNNER_TEMP}/${reference}" \
+                  "pymare/tests/data/${reference}" \
+                  | tee -a "$GITHUB_STEP_SUMMARY" || status=1
+          done
+          exit "$status"
+
+      # Run the alignment tests against the values just regenerated, so this job
+      # checks PyMARE against the R package itself and not only against the pin.
+      - name: "Check PyMARE against the regenerated values"
+        if: always()
+        run: python -m pytest -m ${{ matrix.marker }} -v
+
+      - name: "Upload the regenerated reference values"
+        if: failure()
+        uses: actions/upload-artifact@v4
+        with:
+          name: ${{ matrix.package }}-reference
+          path: pymare/tests/data/*_reference.json

.github/workflows/robumeta-alignment.yml

@@ -1,5 +1,10 @@
 name: "Check robumeta alignment"
 
+# Separate from the "Check alignment with R packages" workflow, which covers
+# metafor and clubSandwich, only so that this job's check name does not change:
+# it predates that workflow and may be a required check on the repository. The
+# steps are the same ones, against the same shared comparator.
+
 on:
   push:
     branches:
@@ -54,11 +59,13 @@ jobs:
       # Compared numerically rather than with git diff: the values are written at
       # full double precision, and which BLAS kernel R's image picks depends on
       # the runner's CPU, so two runners disagree in the last bits with the R and
-      # robumeta versions identical. See the tolerances in compare_reference.py.
+      # robumeta versions identical. See the tolerances in
+      # validation/compare_reference.py, which the metafor and clubSandwich legs
+      # of the "Check alignment with R packages" workflow share.
       - name: "Require the reference values to still match robumeta"
         run: |
           git show HEAD:pymare/tests/data/robumeta_reference.json > "${RUNNER_TEMP}/pinned.json"
-          python validation/robumeta/compare_reference.py \
+          python validation/compare_reference.py \
               "${RUNNER_TEMP}/pinned.json" \
               pymare/tests/data/robumeta_reference.json \
               | tee -a "$GITHUB_STEP_SUMMARY"

CONTRIBUTING.md

@@ -93,18 +93,27 @@ and so on. Fixtures live in `pymare/tests/conftest.py` and helpers that are
 neither fixtures nor tests live in `pymare/tests/utils.py`, so a test file holds
 only tests.
 
-Two groups of tests are marked, because they need something the default
-environment does not have:
+Some tests are marked, because they need something the default environment does
+not have, or because they cost more than a pull request should:
 
 | Target | What it runs | Needs |
 | --- | --- | --- |
 | `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 test_metafor` | the metafor alignment tests | nothing extra |
+| `make test_clubsandwich` | the clubSandwich alignment tests | nothing extra |
 | `make check_robumeta_alignment` | regenerates the robumeta reference values | Docker |
+| `make check_metafor_alignment` | regenerates the metafor reference values | Docker |
+| `make check_clubsandwich_alignment` | regenerates the clubSandwich 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 |
 
+The three `test_*` alignment targets need nothing extra because they read pinned
+numbers; only regenerating those numbers needs R, and that is what the
+`check_*_alignment` targets do. `make unittest` runs the alignment tests too, so
+you do not have to remember them.
+
 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.
 
@@ -134,20 +143,54 @@ 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
+### Alignment with R packages
 
-`pymare/tests/test_robumeta_alignment.py` pins PyMARE's correlated-effects model
-against the R package [robumeta][link_robumeta], over every combination of model,
-rho and variance column that both implementations can express. robumeta cannot be
-a test dependency, so its output is pinned in
-`pymare/tests/data/robumeta_reference.json`.
+Most of what PyMARE computes has a reference implementation in R, and six test
+modules pin PyMARE against one:
 
-`make check_robumeta_alignment` regenerates that file inside a Docker image with
-pinned R and robumeta versions, and fails if any number moved. The
-`Check robumeta alignment` workflow runs the same script on every pull request,
-so a change to the estimator that breaks agreement shows up as a failing check
-rather than as a stale pin. If you changed the estimator on purpose, rerun the
-script and commit the regenerated file.
+| Module | Pins against | Covers |
+| --- | --- | --- |
+| `test_robumeta_alignment.py` | [robumeta][link_robumeta] | the correlated-effects working model (`weight_scheme="rescale"`) |
+| `test_metafor_alignment.py` | [metafor][link_metafor] `rma.uni` | the Knapp-Hartung adjustment, over the whole inference path |
+| `test_metafor_random_effects.py` | metafor `rma.uni`, `confint` | tau^2, Cochran's Q, `I^2`, `H`, and the Q-profile interval |
+| `test_metafor_escalc.py` | metafor `escalc` | the effect-size converters |
+| `test_metafor_permutest.py` | metafor `permutest` | the exact permutation test |
+| `test_clubsandwich_alignment.py` | [clubSandwich][link_clubsandwich] | the CR2 covariance and its Satterthwaite degrees of freedom |
+
+None of those packages can be a test dependency, so their output is pinned under
+`pymare/tests/data/*_reference.json`. Each `validation/<package>/` directory
+holds the R script that produced its file, a Dockerfile pinning the R and package
+versions, a `regenerate.sh` that runs one against the other, and a README
+recording what agrees, to what tolerance, and what does not and why. **Read the
+README before changing an estimator**: it is where the known divergences are
+written down, and several of them are deliberate.
+
+`make check_<package>_alignment` regenerates the files and fails if any number
+moved. The `Check alignment with R packages` and `Check robumeta alignment`
+workflows run the same scripts on every pull request, so a change that breaks
+agreement shows up as a failing check rather than as a stale pin. If you changed
+an estimator on purpose, rerun the script and commit the regenerated file.
+
+The regenerated file is compared numerically, by
+`validation/compare_reference.py`, rather than with `git diff`. The numbers are
+written at full double precision and R reaches them through linear algebra whose
+last bits depend on which BLAS kernel its image picks for the CPU it runs on, so
+two machines with identical R and package versions produce files that differ in
+the last digits. The tolerances are in that script, set below the ones the
+alignment tests themselves rely on.
+
+#### Divergences recorded as strict xfails
+
+Some alignment tests are marked `xfail(strict=True)`, because measuring a
+quantity against its reference implementation turned up a defect rather than a
+divergence. The marker names the expression or code path and what it should be.
+`strict` is load-bearing: correcting the defect turns the test green, pytest
+reports XPASS as a *failure*, and the marker has to be removed in the same
+change. So if a fix of yours makes one of these pass, delete its marker -- do not
+work around it.
+
+That is not hypothetical: it is how the marker covering the two-sample Cohen's d
+variance was retired when [PR #144][link_pr144] landed.
 
 ### Benchmarks
 
@@ -203,4 +246,7 @@ You're awesome.
 [link_stemmrolemodels]: https://github.com/KirstieJane/STEMMRoleModels
 [link_zenodo]: https://github.com/neurostuff/PyMARE/blob/master/.zenodo.json
 [link_robumeta]: https://cran.r-project.org/package=robumeta
+[link_metafor]: https://cran.r-project.org/package=metafor
+[link_clubsandwich]: https://cran.r-project.org/package=clubSandwich
+[link_pr144]: https://github.com/neurostuff/PyMARE/pull/144
 [link_asv]: https://asv.readthedocs.io/en/stable/

Makefile

@@ -6,7 +6,7 @@
 # reports their combined coverage rather than only the last one's.
 PYTEST_COV := --cov-append --cov-report=xml --cov=pymare
 
-all_tests: lint unittest test_stan test_robumeta test_metafor
+all_tests: lint unittest test_stan test_robumeta test_metafor test_clubsandwich
 
 help:
 	@echo "Please use 'make <target>' where <target> is one of:"
@@ -16,8 +16,10 @@ 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 "  test_metafor               to run the metafor alignment tests"
+	@echo "  test_clubsandwich          to run the clubSandwich alignment tests"
 	@echo "  check_robumeta_alignment   to regenerate the robumeta reference values (needs Docker)"
 	@echo "  check_metafor_alignment    to regenerate the metafor reference values (needs Docker)"
+	@echo "  check_clubsandwich_alignment  to regenerate the clubSandwich reference values (needs Docker)"
 	@echo "  validate_stan              to re-measure the Stan model's bias and coverage (~10 min)"
 	@echo "  validate_knapp_hartung     to re-measure the small-sample corrections (~20 min)"
 	@echo "  benchmark                  to run the asv suite once in the current environment"
@@ -45,6 +47,9 @@ test_robumeta:
 test_metafor:
 	@python -m pytest -m "metafor" $(PYTEST_COV)
 
+test_clubsandwich:
+	@python -m pytest -m "clubsandwich" $(PYTEST_COV)
+
 # Re-measures the Type I error of the small-sample corrections and fails if any cell
 # misses the thresholds the default rests on. Not wired into CI and nothing is
 # pinned from it: these are Monte Carlo estimates, so re-measuring is the honest
@@ -59,19 +64,41 @@ validate_knapp_hartung:
 validate_stan:
 	@python validation/stan/simulate.py --check
 
-# What the "Check robumeta alignment" workflow runs. Needs Docker, because the
-# reference values come from R.
+# What the alignment workflows run. Needs Docker, because the reference values
+# come from R.
+#
+# Each target regenerates its package's reference files in place and then
+# compares them against the pinned copies from git. Numerically, through
+# validation/compare_reference.py, rather than with `git diff --exit-code`,
+# which the first two of these used to use: the numbers are written at full
+# double precision and R reaches them through linear algebra whose last bits
+# depend on which BLAS kernel its image picks for the CPU it runs on, so two
+# machines with identical R and package versions produce files that differ in
+# the 16th digit. A byte comparison reads that as drift.
+#
+# Note that these rewrite the working tree, so the comparison is against HEAD
+# rather than against the files on disk.
+COMPARE_REFERENCE = validation/compare_reference.py
+
+define compare_against_head
+	pinned="$$(mktemp)"; git show HEAD:pymare/tests/data/$(1) > "$$pinned"; \
+		python $(COMPARE_REFERENCE) "$$pinned" pymare/tests/data/$(1); \
+		status=$$?; rm -f "$$pinned"; exit $$status
+
+endef
+
 check_robumeta_alignment:
 	@validation/robumeta/regenerate.sh
-	@git diff --exit-code -- pymare/tests/data/robumeta_reference.json \
-		&& echo "The pinned robumeta reference values still match robumeta."
+	@$(call compare_against_head,robumeta_reference.json)
 
-# Regenerates the metafor reference values that pin the Knapp-Hartung adjustment.
-# Needs Docker, for the same reason as the robumeta target.
 check_metafor_alignment:
 	@validation/metafor/regenerate.sh
-	@git diff --exit-code -- pymare/tests/data/metafor_reference.json \
-		&& echo "The pinned metafor reference values still match metafor."
+	@$(foreach reference,metafor_reference.json metafor_escalc_reference.json \
+		metafor_permutest_reference.json,$(call compare_against_head,$(reference)))
+
+check_clubsandwich_alignment:
+	@validation/clubsandwich/regenerate.sh
+	@$(call compare_against_head,clubsandwich_reference.json)
 
 # A smoke test of the benchmark suite, not a measurement: --quick takes one
 # sample per benchmark. The Benchmark workflow is what measures, by timing a

pymare/effectsize/base.py

@@ -250,6 +250,19 @@ class OneSampleEffectSizeConverter(EffectSizeConverter):
     summaries, and are _not_ individual data points. E.g., do not pass in
     a vector of point estimates as `m` and a scalar for the SDs `sd`.
     The lengths of all inputs must match.
+
+    .. versionchanged:: 0.0.13
+
+        The sampling variance of a raw correlation (``'R'``) is now
+        ``(1 - r**2)**2 / (n - 1)``, matching ``metafor::escalc(measure="COR")``.
+        It was ``(1 - r**2) / (n - 2)``, which is the squared standard error of
+        ``r`` under the null hypothesis of *no* correlation rather than its
+        sampling variance at the observed value. The two agree near ``r = 0``
+        and diverge by a factor of ``(n - 1) / ((n - 2)(1 - r**2))`` -- 51x at
+        ``r = 0.99``. Since that factor depends on the data, the old expression
+        did not simply inflate variances: it reweighted studies against one
+        another, pulling a pooled estimate toward those with the weakest
+        correlations. ``'ZR'`` is unaffected and remains the measure to prefer.
     """
 
     _type = 1
@@ -272,7 +285,12 @@ def to_dataset(self, measure="RM", **kwargs):
                     a bias correction applied.
                 - 'D': Cohen's d. Note that no bias correction is applied
                     (use 'SM' instead).
-                - 'R': Raw correlation coefficient.
+                - 'R': Raw correlation coefficient. Prefer 'ZR' for
+                    meta-analysis: the sampling variance of a raw correlation
+                    depends strongly on the correlation itself, so studies are
+                    weighted very unequally by how large their correlations
+                    happen to be, which is the problem the Fisher transform
+                    exists to remove.
                 - 'ZR': Fisher z-transformed correlation coefficient.
         **kwargs
             Optional keyword arguments to pass onto the Dataset

pymare/effectsize/expressions.json

@@ -15,7 +15,7 @@
         "description": "Cohen's d (one-sample)"
     },
     {
-        "expression": "v_d - ((n - 1)/(n - 3)) * (1 / n + d**2) - d**2 / j**2 * n",
+        "expression": "v_d - (((n - 1)/(n - 3)) * (1 / n + d**2) - d**2 / j**2)",
         "type": 1,
         "description": "Variance of Cohen's d"
     },
@@ -30,7 +30,7 @@
         "description": "Standardized mean (Hedges's g)"
     },
     {
-        "expression": "v_sm - ((n - 1)/(n - 3)) * j**2 * (1 / n + d**2) - d**2",
+        "expression": "v_sm - (((n - 1)/(n - 3)) * j**2 * (1 / n + d**2) - d**2)",
         "type": 1,
         "description": "Variance of standardized mean"
     },
@@ -40,7 +40,7 @@
         "description": "Raw correlation coefficient"
     },
     {
-        "expression": "v_r - (1 - r**2) / (n - 2)",
+        "expression": "v_r - (1 - r**2)**2 / (n - 1)",
         "type": 1,
         "description": "Variance of raw correlation coefficient"
     },
@@ -60,7 +60,7 @@
         "description": "Raw mean difference"
     },
     {
-        "expression": "v_rmd - (sd1**2 / n1) + (sd2**2 / n2)",
+        "expression": "v_rmd - ((sd1**2 / n1) + (sd2**2 / n2))",
         "type": 2,
         "description": "Variance of raw mean difference"
     },