From b7ee23942e11c6a12386283142cf784c00e28b03 Mon Sep 17 00:00:00 2001 From: PJ Date: Tue, 18 Aug 2026 20:13:05 +0530 Subject: [PATCH] feat(analyze): add the gehan generalized wilcoxon test The rank-sum carried over to right-censored samples: every pair of runs is scored by which one outlived the other, and a pair censoring cannot order counts as half rather than as a difference neither run supports. The effect size and the p-value are the same statistic, and with nothing censored both are exactly what the rank-sum reports. --- cmd/internal-tools/analyze/gehan.go | 97 +++++++++ cmd/internal-tools/analyze/gehan_test.go | 256 +++++++++++++++++++++++ 2 files changed, 353 insertions(+) create mode 100644 cmd/internal-tools/analyze/gehan.go create mode 100644 cmd/internal-tools/analyze/gehan_test.go diff --git a/cmd/internal-tools/analyze/gehan.go b/cmd/internal-tools/analyze/gehan.go new file mode 100644 index 0000000..5c496c2 --- /dev/null +++ b/cmd/internal-tools/analyze/gehan.go @@ -0,0 +1,97 @@ +package main + +import "math" + +// gehanResult is one pairwise comparison of two arms of right-censored runs: an +// effect size, the run pairs that have no order between them, and the test of +// the same statistic against the null of equal hazards. +type gehanResult struct { + FirstSize int + SecondSize int + Statistic float64 + A12 float64 + Unordered int + PValue float64 +} + +// outlives orders two runs the only way right-censoring allows. A run censored +// at step t violated at no step up to t and stopped for a reason of its own, so +// it outlives a violation at or before t and nothing orders it against a +// violation after t or against another censored run. Comparing the two step +// counts as plain numbers instead reads a run the wall clock stopped at step 12 +// as one that violated at step 12. +func outlives(left, right observation) int { + switch { + case left.Event && right.Event: + switch { + case left.Steps > right.Steps: + return 1 + case left.Steps < right.Steps: + return -1 + } + case left.Event: + if right.Steps >= left.Steps { + return -1 + } + case right.Event: + if left.Steps >= right.Steps { + return 1 + } + } + return 0 +} + +// atRiskWeight is Gehan's weight: an event counts for as many runs as were still +// at risk when it happened. It is what makes the weighted log-rank statistic the +// same quantity as the pairwise count below, so the effect size and the p-value +// are one statistic rather than two that can disagree. +func atRiskWeight(atRisk float64) float64 { return atRisk } + +// gehanTest is the Gehan-Breslow generalized Wilcoxon test: the rank-sum +// carried over to right-censored samples by scoring every pair of runs by which +// one outlived the other and leaving the pairs censoring cannot order out of the +// count. Gehan (1965), "A Generalized Wilcoxon Test for Comparing Arbitrarily +// Singly-Censored Samples", Biometrika 52(1-2), 203-223; Breslow (1970). +// +// Statistic is that count, U, and A12 is it over the number of pairs: the share +// of run pairs in which the first arm survived longer, an unordered pair +// counting as half. With nothing censored the two are exactly the Mann-Whitney U +// and the Vargha-Delaney A12 the uncensored rank-sum reports. Where censoring +// leaves a pair unordered, the half it contributes is the null value, so an +// unordered pair can only pull the effect size toward 0.5 and can never +// manufacture a direction. +// +// The p-value is the same statistic standardized: the weighted log-rank with +// Gehan's weight has this U for its statistic, and its variance is the +// conditional hypergeometric one summed over event times, which is what keeps +// the test honest when the arms censor on different schedules. The permutation +// variance Gehan originally paired with the statistic does not. +func gehanTest(first, second []observation) gehanResult { + result := gehanResult{ + FirstSize: len(first), + SecondSize: len(second), + Statistic: math.NaN(), + A12: math.NaN(), + PValue: math.NaN(), + } + if len(first) == 0 || len(second) == 0 { + return result + } + outlived := 0.0 + for _, left := range first { + for _, right := range second { + switch outlives(left, right) { + case 1: + outlived++ + case 0: + outlived += 0.5 + result.Unordered++ + } + } + } + result.Statistic = outlived + result.A12 = outlived / float64(len(first)*len(second)) + test := weightedLogRank([]string{"first", "second"}, [][]observation{first, second}, atRiskWeight) + result.PValue = test.PValue + return result +} diff --git a/cmd/internal-tools/analyze/gehan_test.go b/cmd/internal-tools/analyze/gehan_test.go new file mode 100644 index 0000000..fe6af21 --- /dev/null +++ b/cmd/internal-tools/analyze/gehan_test.go @@ -0,0 +1,256 @@ +package main + +import ( + "math" + "testing" +) + +func tiedPairs(first, second []float64) int { + tied := 0 + for _, left := range first { + for _, right := range second { + if left == right { + tied++ + } + } + } + return tied +} + +func events(values []float64) []observation { + items := make([]observation, 0, len(values)) + for _, value := range values { + items = append(items, observation{Steps: value, Event: true}) + } + return items +} + +// Every ordering censoring supports and every ordering it does not. +func TestOutlives_OrdersOnlyWhatTheCensoringSupports(t *testing.T) { + cases := []struct { + name string + left, right observation + want int + }{ + {"two violations", observation{30, true}, observation{10, true}, 1}, + {"two violations the other way", observation{10, true}, observation{30, true}, -1}, + {"two violations at the same step", observation{10, true}, observation{10, true}, 0}, + {"censored after the violation", observation{30, false}, observation{10, true}, 1}, + {"censored on the violation's own step", observation{10, false}, observation{10, true}, 1}, + {"censored before the violation", observation{10, false}, observation{30, true}, 0}, + {"violation before the censoring", observation{10, true}, observation{30, false}, -1}, + {"violation after the censoring", observation{30, true}, observation{10, false}, 0}, + {"both censored", observation{10, false}, observation{30, false}, 0}, + } + for _, test := range cases { + if got := outlives(test.left, test.right); got != test.want { + t.Errorf("%s: %v against %v ordered %+d, want %+d", test.name, test.left, test.right, got, test.want) + } + } +} + +// With nothing censored the test is the rank-sum, so its statistic and effect +// size have to be the ones the rank-sum reports on the same numbers, ties +// included. +func TestGehanTest_ReducesToTheRankSumWhenNothingIsCensored(t *testing.T) { + cases := [][2][]float64{ + {chorioamnionTerm, chorioamnionEarly}, + {{1, 2, 3, 4}, {3, 4, 5, 6}}, + {{40, 40, 40}, {40, 40, 40, 40}}, + {{5, 6, 7}, {1, 2}}, + {{3}, {9}}, + } + for _, test := range cases { + result := gehanTest(events(test[0]), events(test[1])) + reference := rankSum(test[0], test[1]) + if result.Statistic != reference.Statistic { + t.Errorf("u %v over %v and %v, want the rank-sum's %v", + result.Statistic, test[0], test[1], reference.Statistic) + } + if result.A12 != reference.A12 { + t.Errorf("a12 %v over %v and %v, want the rank-sum's %v", + result.A12, test[0], test[1], reference.A12) + } + if want := tiedPairs(test[0], test[1]); result.Unordered != want { + t.Errorf("%d unordered pair(s) over %v and %v, want the %d tied ones and no others", + result.Unordered, test[0], test[1], want) + } + } +} + +// The failure the flattening produced: an arm the wall clock stopped at step 12 +// says nothing about step 100, so there is no difference to find and no +// direction to report. +func TestGehanTest_RunsStoppedBeforeEveryViolationOrderNothing(t *testing.T) { + stopped := []observation{{12, false}, {12, false}, {12, false}, {12, false}} + violated := []observation{{100, true}, {100, true}, {100, true}} + result := gehanTest(stopped, violated) + + if result.Unordered != 12 || result.A12 != 0.5 { + t.Errorf("%d of 12 pairs unordered, a12 %v, want all of them and 0.5", result.Unordered, result.A12) + } + if result.PValue < 0.05 { + t.Errorf("p %v, want no difference between arms never observed over the same steps", result.PValue) + } +} + +// A censored run outliving a violation is evidence, and it is the only kind the +// wall-clock case leaves: four runs still clean at step 12 against three +// violations by step 5. +func TestGehanTest_CensoringLeavesTheEvidenceItDoesSupport(t *testing.T) { + stopped := []observation{{12, false}, {12, false}, {12, false}, {12, false}} + violated := []observation{{5, true}, {5, true}, {5, true}} + result := gehanTest(stopped, violated) + + if result.Unordered != 0 || result.Statistic != 12 || result.A12 != 1 { + t.Errorf("result %+v, want every pair ordered for the arm that had not violated", result) + } + if result.PValue > 0.05 { + t.Errorf("p %v, want the arms to separate", result.PValue) + } +} + +func riskAndDeaths(group []observation, steps float64) (float64, float64) { + atRisk, deaths := 0.0, 0.0 + for _, item := range group { + if item.Steps >= steps { + atRisk++ + } + if item.Steps == steps && item.Event { + deaths++ + } + } + return atRisk, deaths +} + +// gehanReference is Gehan's statistic and its conditional variance written +// straight from the definitions, +// +// S = sum over event times of (Y2*d1 - Y1*d2) +// V = sum over event times of d(Y-d)/(Y-1) * Y1*Y2 +// +// which is an independent calculation rather than a second call into the code +// under test. +func gehanReference(first, second []observation) (float64, float64) { + pooled := append(append([]observation{}, first...), second...) + statistic, variance := 0.0, 0.0 + for _, steps := range distinctSteps(pooled) { + firstAtRisk, firstDeaths := riskAndDeaths(first, steps) + secondAtRisk, secondDeaths := riskAndDeaths(second, steps) + deaths := firstDeaths + secondDeaths + if deaths == 0 { + continue + } + atRisk := firstAtRisk + secondAtRisk + statistic += secondAtRisk*firstDeaths - firstAtRisk*secondDeaths + if atRisk > 1 { + variance += deaths * (atRisk - deaths) / (atRisk - 1) * firstAtRisk * secondAtRisk + } + } + return statistic, variance +} + +// The effect size and the p-value have to be the same statistic seen twice, or +// the report can carry a direction its p-value does not support. Counting run +// pairs and accumulating over risk sets are two routes to Gehan's statistic, and +// they are tied by S = mn - 2U. +func TestGehanTest_PairCountAndRiskSetAgreeOnOneStatistic(t *testing.T) { + cases := []struct { + name string + first, second []observation + }{ + {"6-mp against placebo", gehanSixMercaptopurine, gehanPlacebo}, + {"maintained against nonmaintained", amlMaintained, amlNonmaintained}, + {"stopped short against violating late", []observation{{12, false}, {12, false}, {14, false}}, + []observation{{5, true}, {100, true}, {100, true}}}, + {"censoring tied with an event", []observation{{20, false}, {20, true}, {35, true}}, + []observation{{20, true}, {20, false}, {9, true}}}, + } + for _, test := range cases { + result := gehanTest(test.first, test.second) + statistic, variance := gehanReference(test.first, test.second) + pairs := float64(len(test.first) * len(test.second)) + + if got := pairs - 2*result.Statistic; math.Abs(got-statistic) > 1e-9 { + t.Errorf("%s: pair count gives a statistic of %v, the risk sets give %v", test.name, got, statistic) + } + expected := chiSquareUpperTail(statistic*statistic/variance, 1) + if math.Abs(result.PValue-expected) > 1e-12 { + t.Errorf("%s: p %v, want %v from statistic %v over variance %v", + test.name, result.PValue, expected, statistic, variance) + } + } +} + +// The 6-MP trial is the dataset the test is named for. Its log-rank result is +// checked elsewhere against the published one; here the generalized Wilcoxon +// has to reach the same conclusion, with the maintained arm outliving the +// placebo arm on both routes. +func TestGehanTest_SeparatesThePublishedLeukaemiaTrial(t *testing.T) { + result := gehanTest(gehanSixMercaptopurine, gehanPlacebo) + if result.A12 <= 0.5 { + t.Errorf("a12 %v, want the 6-mp arm to outlive the placebo arm", result.A12) + } + if result.PValue > 0.001 { + t.Errorf("p %v, want the arms to separate as the log-rank has them separate", result.PValue) + } + logRankResult := logRank([]string{"6-mp", "placebo"}, + [][]observation{gehanSixMercaptopurine, gehanPlacebo}) + if logRankResult.PValue > 0.001 { + t.Fatalf("log-rank p %v: the comparison being made is not the one this test assumes", logRankResult.PValue) + } +} + +func censoredAt(count int, steps float64) []observation { + items := make([]observation, 0, count) + for index := 0; index < count; index++ { + items = append(items, observation{Steps: steps}) + } + return items +} + +// A specification whose violations are all obligations reported when the run +// ends puts every event on one step, and the comparison collapses to a single +// two-by-two table of violated against clean. The statistic there is the +// Mantel-Haenszel chi-square of that table, +// +// (N-1)(ad-bc)^2 / ((a+b)(c+d)(a+c)(b+d)) +// +// and the (Y-d)/(Y-1) term in the variance is what carries the (N-1)/N that +// separates it from the Pearson chi-square. Nothing else in the pipeline +// exercises that term hard, because tied events are otherwise rare. +func TestGehanTest_EveryViolationOnOneStepIsTheMantelHaenszelTable(t *testing.T) { + cases := [][4]int{{30, 20, 10, 40}, {25, 25, 15, 35}, {20, 20, 12, 28}} + for _, test := range cases { + firstEvents, firstCensored, secondEvents, secondCensored := test[0], test[1], test[2], test[3] + first := append(events(repeated(firstEvents, 400)), censoredAt(firstCensored, 400)...) + second := append(events(repeated(secondEvents, 400)), censoredAt(secondCensored, 400)...) + + total := float64(firstEvents + firstCensored + secondEvents + secondCensored) + crossProduct := float64(firstEvents*secondCensored - firstCensored*secondEvents) + chiSquare := (total - 1) * crossProduct * crossProduct / + (float64(firstEvents+firstCensored) * float64(secondEvents+secondCensored) * + float64(firstEvents+secondEvents) * float64(firstCensored+secondCensored)) + + got := gehanTest(first, second).PValue + want := chiSquareUpperTail(chiSquare, 1) + if math.Abs(got-want) > 1e-12 { + t.Errorf("%v: p %v, want the table's %v from chi-square %v", test, got, want, chiSquare) + } + } +} + +func repeated(count int, value float64) []float64 { + values := make([]float64, 0, count) + for index := 0; index < count; index++ { + values = append(values, value) + } + return values +} + +func TestGehanTest_EmptyArmHasNoComparison(t *testing.T) { + result := gehanTest(nil, []observation{{5, true}}) + if !math.IsNaN(result.PValue) || !math.IsNaN(result.A12) || !math.IsNaN(result.Statistic) { + t.Errorf("result %+v, want everything undefined", result) + } +}