/ 
Stats-LikeR-0.313

Changes for version 0.313 - 2026-09-01

  • `merge`: a `left.on`/`right.on` join died when the right frame reused the key's name
    • The result of a join carries one key column, under the left name — R's convention, and not pandas', which keeps both. That left one collision unaccounted for: with `left.on`/`right.on`, a *non-key* column on the right can be named after the *left key*, and it then collides with the output key column rather than with a left data column. The suffix rule only looked at the left frame's data columns, so nothing was renamed and the collision guard fired:
      • merge($parents, $children, 'left.on' => 'name', 'right.on' => 'parent');
      • merge: output column 'name' collides; adjust 'suffixes'
    • That is `tests/reg-tests-1d.R`'s parents/children join, and both references perform it: R suffixes the right-hand copy under `no.dups = TRUE` (the default since R 3.5.0 — "if a `by.x` column name matches one of `y`, the y version gets suffixed as well"), and pandas keeps both key columns so the question never arises. `merge` now suffixes it too, giving R's names exactly:
    • | | columns | |---|---| | R 4.6.1 | `name`, `sex.x`, `age.x`, `name.y`, `sex.y`, `age.y` | | `merge` before | *croaks* | | `merge` now | `name`, `sex.x`, `age.x`, `name.y`, `sex.y`, `age.y` |
    • Every input this changes used to croak, so no join that worked before returns anything different. If the suffixes still leave two output columns sharing a name, `merge` still dies rather than hand back a frame with a column missing — which is R's behaviour too (`suffixes = c(".z", ".z")` is an error there).
  • `merge` is now cross-validated against R's and pandas' own merge suites
    • `t/merge.R.pandas.t` (329 tests) takes its cases from the references' test suites rather than inventing them, in the manner of the other `t/*.R.scipy.t` files:
    • **R 4.6.1** — 28 frames and 83 cases: the examples in `src/library/base/man/merge.Rd` (authors/books, and the `incomparables` example), plus R's own regression cases for `merge` from `reg-tests-1a.R` (PR#1510 `by.x`/`by.y` with multiple matches; the Cartesian product that did not make column names unique in 2.3.0; "merging when NA is a level"; the two character matrices that failed pre-2.0.0; merge on zero-row frames, not allowed ≤ 2.4.0), `reg-tests-1b.R` (the 2.15.0 and 2.15.1 suffixes regressions; the `women` zero-row merges that failed in 2.7.0), `reg-tests-1d.R` (the `by.y` naming case above) and `reg-tests-2.R` (the authors/books joins moved out of `merge.Rd`, and the 2002 case where every column is a join key).
    • **pandas 2.2.3** — 31 frames and 60 cases from `pandas/tests/reshape/merge/test_merge.py` (`test_intelligently_handle_join_key` GH#733, `test_merge_overlap`, `test_merge_different_column_key_names`, `test_merge_same_order_left_right` GH#35382, `test_left_merge_empty_dataframe`, all ten parametrisations of `test_merge_empty` GH#52777, `test_merge_on_ints_floats`, `test_merge_non_unique_index_many_to_many`, `test_merge_suffix`), `test_merge_cross.py` (all four cross-join tests) and `test_multi.py` (`test_merge_na_keys`, `test_merge_multiple_cols_with_mixed_cols_index` GH#29522).
    • Each frozen case runs through every input/output shape — AoH, HoA and HoH on the left crossed with the same on the right and with both output shapes, 18 joins per case — and its output column names are checked separately, because a join whose answer has no rows is exactly what several of the reference cases are about and a row-by-row comparison cannot see the names. Beyond the tables the file pins row order against pandas' `sort=False`, R's `by`/`by.x`/`by.y` and pandas' `left_on`/`right_on` spellings, every error path the references document, and double-valued keys with dyadic literals.
    • The two tables are generated by `t/merge.R.pandas.R` and `t/merge.R.pandas.py`, committed beside the test; the test itself never runs R or python, and needs neither installed. Both blocks reproduce byte-identically from the generators, and the pandas block is byte-identical under pandas 2.2.3 and 3.0.4, so nothing frozen there is specific to a release.
    • `t/merge.t` also gained the cases a data frame cannot express, and so no reference case can reach: an AoH row missing the key column entirely, a missing non-key cell, the `0` / `'0'` / `'0.0'` key identity, a NUL byte inside a key, a wide character, `Inf`/`-Inf`/`NaN` keys, and blessed frames.
  • An `undef` join key matches nothing — which is *not* "the pandas `NaN` rule"
    • The documentation for `merge` credited its treatment of a missing join key to pandas, and `LikeR.xs` credited it to R's default. Both attributions were wrong, in the same direction: both references match a missing key to a missing key. On `merge.Rd`'s own `incomparables` example the two agree with each other exactly and disagree with `merge`:
    • | join | R 4.6.1 | pandas 2.2.3 | `merge` | |---|---|---|---| | `by = c("k1","k2")` | 3 rows | 3 rows | 2 rows | | `by = "k1"` | 6 rows | 6 rows | 2 rows |
    • The behaviour is right and unchanged — it is SQL's rule for a `NULL` key, which is `merge(..., incomparables = NA)` in R, the line `merge.Rd` itself runs, and what `merge`'s documentation describes everywhere else. Only the credit was wrong; R's default is `incomparables = NULL`. The corrected wording says which rule it is and that it is not either reference's default, and `t/merge.R.pandas.t` now asserts the divergence — with R's own `incomparables = NA` answers frozen beside the generalised ones — so changing it later has to be deliberate.
  • 11 leak checks reported Devel::Cover's own counters as leaks
    • `t/vif_hoslem.t` failed under `cover` with 37 leaks attributed to a line of `hosmer_lemeshow` that allocates nothing (`next unless defined ... && looks_like_number ...`), and the leaked SVs were bare `IV`s holding values like `98113961175368` — pointers. They are Devel::Cover's per-line counters, which are allocated inside whichever block happens to be running and which `Test::LeakTrace` then counts; the same test reports 0 leaks on a plain perl.
    • Around ninety test files already guard their leak checks with `$INC{'Devel/Cover.pm'}` for this reason. Eleven did not: `t/age_standardize.t`, `t/dunn_test.t`, `t/effect_sizes.t`, `t/friedman_test.t`, `t/glm_families.t`, `t/ks_test.R.scipy.t`, `t/mcnemar_test.t`, `t/prop_test.t`, `t/tied.frames.t`, `t/vif_hoslem.t` and — for an import it never used — `t/_parse_csv.t`. Four of them failed; the rest passed only by luck, since whether the counters land inside the measured block depends on which lines get their first coverage there, which moves with file order. All of them are guarded now, in the two forms the suite already uses, and the whole suite passes under `Devel::Cover` (143 files, 35008 tests) as well as without it (143 files, 35577 tests — the difference is the leak checks, which still run and still report 0 leaks outside coverage mode).
  • `pt` near zero and `pf`'s upper tail lost up to eight digits to a cancellation
    • `incbeta()`, the regularized incomplete beta every t, F and binomial tail is built on, took only `x` and re-formed `1 - x` by subtraction. Its reflected branch — `I_x(a,b) = 1 - I_{1-x}(b,a)`, taken exactly when `x` is the side near `1` — then needs that complement, and forming it as `1.0 - x` is catastrophic cancellation there. Once `|1 - x|` fell below `2^-53` it collapsed to `0` and took the whole tail with it:
    • | call | returned | correct to 21 digits | |---|---|---| | `pf(1e-12, 1, 1e6, 'lower.tail' => 0)` | `1` exactly | `0.999999202115638638` | | `pt(1e-6, 1e6)` | `0.5` exactly | `0.500000398942180682` | | `pt(-1e-8, 1)` | `0.5` exactly | `0.499999996816901138` |
    • Every caller had the complement exactly, and was throwing it away: `pf` builds `x = df1·f/D` and `1-x = df2/D` over one denominator, `pt` has `x = df/(df + t²)` against `t²/(df + t²)`, and `pbinom`'s lower tail is `I_{1-p}(n-k, k+1)` with `p` itself to hand. So the fix is to pass both — `incbeta_xy(a, b, x, y)`, which is R's own split (`bratio()` takes `x` and `y` as separate arguments for this reason), with `incbeta()` kept as the one-argument wrapper for the callers that genuinely have only `x`, such as a bisection midpoint.
    • Worst error against `mpmath` at `mp.dps = 80`, over a 372-point grid:
    • | | before | after | |---|---|---| | `pt`, either tail | `3.99e-07` absolute | `1.11e-16` (1 ulp) | | `pf` upper tail | `3.19e-08` relative | `2.74e-11` | | `pf` lower tail | `2.01e-12` relative | `2.46e-13` |
    • The two tails also add up again: `pf(1e-12, 1, 1e6)`'s lower and upper summed to `1 + 7.98e-07` before, which is what first showed the bug.
    • This reached `t_test`, whose p-value is `1 -` that tail: a statistic small enough on enough degrees of freedom had its p-value pinned at exactly `1` instead of `1 - 4e-7`. `d_pf` now calls `pf()` rather than repeating its expression, so the two can no longer drift apart on the complement argument.
    • The one place `incbeta_xy` does not help is the Clopper-Pearson upper bound for a handful of successes in ~1e9 trials, which still carries ~2e-9 of relative error. The cancellation there is in the continued fraction's own argument during bisection, not in a complement a caller could have supplied, so `t/binom_test.R.scipy.t`'s existing note — that fixing it properly means porting `bratio()` — still stands.
  • `t_test` rejected `conf.level`, the spelling its own documentation lists
    • t_test(\@x, \@y, 'conf.level' => 0.99);
    • t_test: unknown argument 'conf.level'
    • `t_test` accepted only the underscored `conf_level` and `var_equal`, while its parameter table in this file has always documented the argument as `conf.level`, and while every sibling in the module — `var_test`, `wilcox_test`, `prop_test`, `cmh_test`, `glm` and the rest — already took both spellings. `var.equal`, which is what R calls it, was refused too. Both dotted forms now work, and the parameter table records the aliases.
  • `rank`, `wilcox_test` and `ks_test` are 1.6× to 2.3× faster
    • These were still sorting through `qsort()`, whose comparator the compiler cannot inline; the module's own `LIKER_DEFINE_SORT()` introsort and `nv_sort()` were already used by the three internal rankers but not by these. Ordering *n* records costs O(*n* log *n*) comparisons and an indirect call on each is most of what such a sort costs — the figure `LIKER_DEFINE_SORT`'s own comment records is 210µs against `qsort()`'s 355µs on 5000 NVs.
    • Measured on 20,000 doubles:
    • | | before | after | speedup | |---|---|---|---| | `wilcox_test` (two samples) | 6.11 ms | 2.64 ms | 2.3× | | `rank` | 2.87 ms | 1.43 ms | 2.0× | | `ks_test` (two samples) | 3.70 ms | 2.37 ms | 1.6× |
    • `rank_and_count_ties()`, which `wilcox_test`, `ks_test` and five other functions share, had been sorting `RankInfo` records through `cmp_nv3` — a comparator that reads a bare `NV`. That worked only because `val` is the struct's first member and a pointer to a struct is a pointer to its first member: true, but fragile as well as slow. It has a generated ordering of its own now.
    • What is left in `rank` is no longer the sort. Reading the same 20,000 values back out of an `AV` in plain perl costs 0.37 ms, so the gather loop and the `newSVnv` per result are now most of the call.
  • Five new cross-validation files, from R's and SciPy's own test suites
    • 4,444 tests, taking their cases from the references' suites and documented examples rather than inventing them, in the manner of the existing `t/*.R.scipy.t` files. Expected values are frozen literals; the generators are committed beside each test and are never run by it, so nothing here needs R, python or `mpmath` at install time.
    • | file | tests | sources | |---|---|---| | `t/t_test.R.scipy.t` | 2437 | `t.test.Rd`; `reg-tests-1a.R:4529` (one group of size one), `reg-tests-2.R:3199`, `reg-tests-1e.R:1985`; SciPy's `TestTTest_1samp`, `TestTTestIndMore`, `TestTTestRel`, `TestTTestCI` | | `t/tukey_aov_prcomp.R.t` | 940 | a 637-point `ptukey`/`qtukey` grid; PlantGrowth and chickwts `TukeyHSD`; mtcars `anova`/`vif`; USArrests `prcomp`; `scale` | | `t/p_adjust.R.t` | 767 | every method in R's own `p.adjust.methods`, on `p.adjust.Rd`'s own p-vector | | `t/friedman_mcnemar_prop_cmh.R.t` | 274 | Hollander & Wolfe (1973) p.140ff; Agresti (1990) p.350; Fleiss (1981) p.139; Agresti's Rabbits and `UCBAdmissions` | | `t/pf_pt_tails.R.mpmath.t` | 26 | `mpmath` at `mp.dps = 80`, plus R on the same grid |
    • `t_test` had no cross-validation at all before this, which is how the `conf.level` croak above survived. Each documented case is crossed over the whole argument space its function exposes — alternative × `var_equal` × `mu` × `conf.level` × `paired`, `correct`, `exact`, `p` — because a reference case exercised only at its defaults pins one code path out of dozens.
    • Two things fell out of writing them.
    • **`t/tukey.t`'s tolerances understate the `ptukey` port by ten orders of magnitude.** It checks `qtukey` to an absolute `1e-3` and the `TukeyHSD` columns to `1e-4`, where the Copenhaver & Holland port actually agrees with R to `3.0e-14` and `2.5e-12` over the new grid. A `1e-3` limit would not notice the port being replaced by a normal approximation, which is the regression it exists to catch. The new file checks it properly; the old one is left alone.
    • **A tail probability needs an absolute tolerance, not a relative one, across NV widths.** `ptukey`'s internals are deliberately plain `double`, exactly as R's `src/nmath/ptukey.c` has them, so the quadrature does not move with perl's `NV` — but its *argument* `q = |diff| / se` comes from the `aov` mean square, which is computed in `NV`. On the chickwts `horsebean-casein` comparison, where `p adj` is `3.07e-08`:
    • | NV width | `p adj` | relative | absolute | |---|---|---|---| | `double` | `3.0701967967949884e-08` | — | — | | x87 `long double` | `3.0701966635682254e-08` | `4.34e-08` | `1.3e-16` | | `long double` | `3.0701956643675032e-08` | `3.69e-07` | `1.1e-14` | | `__float128` | `3.0701956643675032e-08` | `3.69e-07` | `1.1e-14` |
    • Every width agrees to about `1e-14` absolute, which is all a probability in the `1e-8` tail can be asked for. Conversely `pf` and `pt` come out three to four orders *more* accurate on the wider widths (`2.74e-11` → `1.1e-15` for `pf`'s upper tail), which is the evidence that what is left there is the continued fraction's convergence and not another cancellation — a cancellation does not improve with the working precision.
  • `qf` is more accurate than R's own `qf` in the far lower tail
    • Writing the file above turned up a disagreement in the other direction, so it is recorded rather than quietly reconciled. Asked for the F quantile deep in the lower tail, R bottoms out at a fixed resolution:
    • | | `qf(1e-8, 1, 2)` | |---|---| | `mpmath`, `mp.dps = 80` | `2.0000000000000003e-16` | | `qf` | `2.0000000000000000e-16` | | R 4.6.1 `qf` | `4.4408920985006262e-16` — which is `2^-51` |
    • R is out by 122% there. Over the 108-point grid in `t/pf_pt_tails.R.mpmath.t`, `qf` is within `7.6e-15` relative of 80-digit truth at every point and R is out by as much as `1.22`, on 17 of them. The test asserts `qf` against `mpmath` and separately asserts that R really is the worse of the two, so that "fixing" `qf` towards R would fail loudly instead of passing quietly.
    • Also recorded there: `qtukey` is *not* an exact inverse of `ptukey`, in R or here. R's `qtukey.c` is a secant iteration that stops once successive iterates differ by less than a hardcoded `const static double eps = 0.0001` — an absolute `1e-4` in `q`, not a relative tolerance on `p` — so `ptukey(qtukey(p))` recovers `p` only to about `1e-7`. R's own worst round trip over that grid is `1.2744607e-07`, and this port's is the same to the digits printed, which is the strongest evidence in the file that the port really is faithful rather than merely close.

Documentation

CONTRIBUTING.md
README.md

Modules

Stats::LikeR
Get basic statistical functions, like in R, but with Perl using XS for performance

Other files