diff --git a/app/Console/Commands/SendExperimentFunnelReportCommand.php b/app/Console/Commands/SendExperimentFunnelReportCommand.php index 5b3cd9a2..cd2caaa7 100644 --- a/app/Console/Commands/SendExperimentFunnelReportCommand.php +++ b/app/Console/Commands/SendExperimentFunnelReportCommand.php @@ -4,7 +4,9 @@ namespace App\Console\Commands; use App\Features\SubscriptionExperiment; use App\Services\Discord\DiscordWebhook; +use App\Services\Stats\BinomialProportion; use App\Services\Stats\ExperimentFunnelCollector; +use App\Services\Stats\ProportionSignificance; use Carbon\CarbonImmutable; use Illuminate\Console\Command; @@ -22,8 +24,10 @@ class SendExperimentFunnelReportCommand extends Command SubscriptionExperiment::PAY_NOW => 'pay_now', ]; - public function __construct(private ExperimentFunnelCollector $collector) - { + public function __construct( + private ExperimentFunnelCollector $collector, + private ProportionSignificance $significance, + ) { parent::__construct(); } @@ -78,7 +82,7 @@ class SendExperimentFunnelReportCommand extends Command foreach (self::LABELS as $key => $label) { $row = $report['variants'][$key]; $mature = $row['assignedMature'] > 0; - $money = $revenue && $mature; + $showMoney = $revenue && $mature; $lines[] = sprintf( '%-8s %5d %5d %5d %5d %5d %6s %7s %7s %7s %7s %7s', @@ -89,11 +93,11 @@ class SendExperimentFunnelReportCommand extends Command $row['assignedMature'], $row['convertedMature'], $mature ? ((int) round($row['conversionRate'] * 100)).'%' : 'pend', - $money && $row['arpuCents'] !== null ? $this->money($row['arpuCents'], $currency) : '—', - $money ? $this->money($row['mrrCents'], $currency) : '—', + $showMoney && $row['arpuCents'] !== null ? $this->money($row['arpuCents'], $currency) : '—', + $showMoney ? $this->money($row['mrrCents'], $currency) : '—', $mature ? $this->money($row['costCents'], $currency) : '—', $mature ? $this->money($row['wastedCostCents'], $currency) : '—', - $money ? $this->money($row['contributionMarginCents'], $currency) : '—', + $showMoney ? $this->money($row['contributionMarginCents'], $currency) : '—', ); } @@ -101,19 +105,18 @@ class SendExperimentFunnelReportCommand extends Command } /** - * Conversion-rate uncertainty per variant (95% Wilson interval) plus a - * two-proportion z-test between the two leaders, Bonferroni-corrected for - * the three pairwise comparisons, so "check significance before calling a - * winner" has the numbers behind it instead of just the warning. + * Per-variant conversion-rate uncertainty (95% Wilson interval) plus the + * leader-vs-runner-up verdict from {@see ProportionSignificance} — a Fisher + * exact test and a Newcombe difference interval, Bonferroni-corrected — so + * "check significance before calling a winner" has the numbers behind it. * * @param array{startedAt: ?CarbonImmutable, currency: string, revenueAvailable: bool, costPerConnectionCents: int, variants: array>} $report * @return list */ private function significanceLines(array $report): array { - $z = 1.96; // two-sided 95% $lines = ['', 'Significance (95% Wilson CI on Conv%, n = MatU):']; - $scored = []; + $arms = []; foreach (self::LABELS as $key => $label) { $row = $report['variants'][$key]; @@ -126,154 +129,43 @@ class SendExperimentFunnelReportCommand extends Command continue; } - [$low, $high] = $this->wilsonInterval($k, $n, $z); + [$low, $high] = $this->significance->wilsonInterval($k, $n); $lines[] = sprintf(' %-8s %6s [%6s – %6s] (n=%d)', $label, $this->percent($k / $n), $this->percent($low), $this->percent($high), $n); - $scored[] = ['label' => $label, 'k' => $k, 'n' => $n, 'rate' => $k / $n]; + $arms[] = new BinomialProportion($label, $k, $n); } - if (count($scored) < 2) { + if (count($arms) < 2) { $lines[] = 'Not enough matured variants to compare yet.'; return $lines; } - usort($scored, fn (array $a, array $b): int => $b['rate'] <=> $a['rate']); - [$leader, $runnerUp] = [$scored[0], $scored[1]]; - - // Decide with Fisher's exact test, not a normal-approximation z: at the - // small conversion counts this report has, the pooled z overstates the - // evidence (expected cell counts fall below the np>=5 the z-test needs). - // Fisher is exact at any n. Bonferroni-correct for the 3 pairwise arms. - $alpha = 0.05 / 3; - $delta = ($leader['rate'] - $runnerUp['rate']) * 100; - [$diffLow, $diffHigh] = $this->newcombeDiffInterval($leader['k'], $leader['n'], $runnerUp['k'], $runnerUp['n'], $z); - $fisherP = $this->fisherExactTwoSided( - $leader['k'], $leader['n'] - $leader['k'], - $runnerUp['k'], $runnerUp['n'] - $runnerUp['k'], - ); - $significant = $fisherP < $alpha; + usort($arms, fn (BinomialProportion $a, BinomialProportion $b): int => $b->rate() <=> $a->rate()); + [$leader, $runnerUp] = [$arms[0], $arms[1]]; + $result = $this->significance->compare($leader, $runnerUp); $lines[] = sprintf( 'Leader %s vs %s: Δ %+.1f pts (95%% CI %+.1f … %+.1f pts, Newcombe).', - $leader['label'], $runnerUp['label'], $delta, $diffLow * 100, $diffHigh * 100, + $leader->label, $runnerUp->label, + ($leader->rate() - $runnerUp->rate()) * 100, $result['diffLow'] * 100, $result['diffHigh'] * 100, ); $lines[] = sprintf( 'Fisher exact p=%.3f %s α=%.3f (Bonferroni×3) -> %s.%s', - $fisherP, $significant ? '<' : '≥', $alpha, - $significant ? 'significant' : 'not significant', - $significant ? '' : ' Keep running.', + $result['fisherP'], $result['significant'] ? '<' : '≥', $result['alpha'], + $result['significant'] ? 'significant' : 'not significant', + $result['significant'] ? '' : ' Keep running.', ); - $pooled = ($leader['k'] + $runnerUp['k']) / ($leader['n'] + $runnerUp['n']); - $minExpected = min($leader['n'], $runnerUp['n']) * min($pooled, 1 - $pooled); - if ($minExpected < 5.0) { + if ($result['minExpectedCount'] < 5.0) { $lines[] = sprintf( '(Small sample: min expected conversions %.1f < 5, so the normal-approx z=%.2f overstates — exact test used.)', - $minExpected, $this->twoProportionZ($leader['k'], $leader['n'], $runnerUp['k'], $runnerUp['n']), + $result['minExpectedCount'], $result['z'], ); } return $lines; } - /** - * Wilson score interval for a binomial proportion — accurate for small n - * and near 0/1, where the normal approximation misbehaves. - * - * @return array{0: float, 1: float} lower and upper bound, clamped to [0, 1] - */ - private function wilsonInterval(int $k, int $n, float $z): array - { - $p = $k / $n; - $z2 = $z * $z; - $denom = 1 + $z2 / $n; - $center = ($p + $z2 / (2 * $n)) / $denom; - $margin = ($z / $denom) * sqrt($p * (1 - $p) / $n + $z2 / (4 * $n * $n)); - - return [max(0.0, $center - $margin), min(1.0, $center + $margin)]; - } - - private function twoProportionZ(int $kA, int $nA, int $kB, int $nB): float - { - $pooled = ($kA + $kB) / ($nA + $nB); - $se = sqrt($pooled * (1 - $pooled) * (1 / $nA + 1 / $nB)); - - return $se > 0.0 ? (($kA / $nA) - ($kB / $nB)) / $se : 0.0; - } - - /** - * Newcombe (Wilson-based) 95% interval for the difference of two proportions - * pA − pB. The correct object when asking "is A better than B" — overlapping - * marginal intervals do NOT imply the difference includes 0. - * - * @return array{0: float, 1: float} - */ - private function newcombeDiffInterval(int $kA, int $nA, int $kB, int $nB, float $z): array - { - $pA = $kA / $nA; - $pB = $kB / $nB; - [$lA, $uA] = $this->wilsonInterval($kA, $nA, $z); - [$lB, $uB] = $this->wilsonInterval($kB, $nB, $z); - - $lower = ($pA - $pB) - sqrt(($pA - $lA) ** 2 + ($uB - $pB) ** 2); - $upper = ($pA - $pB) + sqrt(($uA - $pA) ** 2 + ($pB - $lB) ** 2); - - return [$lower, $upper]; - } - - /** - * Two-sided Fisher exact test p-value for the 2x2 table [[a, b], [c, d]] - * (a/c = conversions, b/d = non-conversions). Sums the hypergeometric - * probabilities of every table, with the same margins, no more likely than - * the observed one. Exact at any sample size — no normal approximation. - */ - private function fisherExactTwoSided(int $a, int $b, int $c, int $d): float - { - $rowA = $a + $b; - $rowB = $c + $d; - $col = $a + $c; - $total = $rowA + $rowB; - - if ($rowA === 0 || $rowB === 0 || $col === 0 || $col === $total) { - return 1.0; - } - - $logProbObserved = $this->hypergeometricLogProb($a, $rowA, $rowB, $col); - $p = 0.0; - for ($x = max(0, $col - $rowB); $x <= min($col, $rowA); $x++) { - $logProb = $this->hypergeometricLogProb($x, $rowA, $rowB, $col); - if ($logProb <= $logProbObserved + 1e-7) { - $p += exp($logProb); - } - } - - return min(1.0, $p); - } - - private function hypergeometricLogProb(int $x, int $rowA, int $rowB, int $col): float - { - return $this->logChoose($rowA, $x) + $this->logChoose($rowB, $col - $x) - $this->logChoose($rowA + $rowB, $col); - } - - private function logChoose(int $n, int $k): float - { - if ($k < 0 || $k > $n) { - return -INF; - } - - return $this->logFactorial($n) - $this->logFactorial($k) - $this->logFactorial($n - $k); - } - - private function logFactorial(int $n): float - { - $sum = 0.0; - for ($i = 2; $i <= $n; $i++) { - $sum += log($i); - } - - return $sum; - } - private function percent(float $rate): string { return number_format($rate * 100, 1).'%'; diff --git a/app/Services/Stats/BinomialProportion.php b/app/Services/Stats/BinomialProportion.php new file mode 100644 index 00000000..d6cf923d --- /dev/null +++ b/app/Services/Stats/BinomialProportion.php @@ -0,0 +1,22 @@ +trials > 0 ? (float) $this->successes / $this->trials : 0.0; + } +} diff --git a/app/Services/Stats/ProportionSignificance.php b/app/Services/Stats/ProportionSignificance.php new file mode 100644 index 00000000..344fe98f --- /dev/null +++ b/app/Services/Stats/ProportionSignificance.php @@ -0,0 +1,149 @@ +=5 rule. + */ +final class ProportionSignificance +{ + /** Two-sided 95% standard-normal quantile. */ + private const Z_95 = 1.96; + + /** @var list memoised log-factorials, index i = log(i!) */ + private array $logFactorials = [0.0, 0.0]; + + /** + * Compare the two leading arms: Newcombe difference interval, Fisher exact + * p-value, and a family-wise (Bonferroni) corrected significance flag. `z` + * and `minExpectedCount` are returned for the small-sample caveat only. + * + * @return array{alpha: float, diffLow: float, diffHigh: float, fisherP: float, significant: bool, minExpectedCount: float, z: float} + */ + public function compare(BinomialProportion $leader, BinomialProportion $runnerUp, float $familyAlpha = 0.05, int $comparisons = 3): array + { + $alpha = $familyAlpha / $comparisons; + [$diffLow, $diffHigh] = $this->newcombeDiffInterval($leader, $runnerUp); + $fisherP = $this->fisherExactTwoSided( + $leader->successes, $leader->trials - $leader->successes, + $runnerUp->successes, $runnerUp->trials - $runnerUp->successes, + ); + + $pooled = ($leader->successes + $runnerUp->successes) / ($leader->trials + $runnerUp->trials); + $minExpectedCount = min($leader->trials, $runnerUp->trials) * min($pooled, 1 - $pooled); + + return [ + 'alpha' => $alpha, + 'diffLow' => $diffLow, + 'diffHigh' => $diffHigh, + 'fisherP' => $fisherP, + 'significant' => $fisherP < $alpha, + 'minExpectedCount' => $minExpectedCount, + 'z' => $this->twoProportionZ($leader, $runnerUp), + ]; + } + + /** + * Wilson score interval for a binomial proportion — accurate for small n + * and near 0/1, where the normal approximation misbehaves. + * + * @return array{0: float, 1: float} lower and upper bound, clamped to [0, 1] + */ + public function wilsonInterval(int $successes, int $trials, float $z = self::Z_95): array + { + $p = $successes / $trials; + $z2 = $z * $z; + $denom = 1 + $z2 / $trials; + $center = ($p + $z2 / (2 * $trials)) / $denom; + $margin = ($z / $denom) * sqrt($p * (1 - $p) / $trials + $z2 / (4 * $trials * $trials)); + + return [max(0.0, $center - $margin), min(1.0, $center + $margin)]; + } + + /** + * Newcombe (Wilson-based) 95% interval for the difference pA − pB. The + * correct object for "is A better than B": overlapping marginal intervals do + * NOT imply the difference includes 0. + * + * @return array{0: float, 1: float} + */ + public function newcombeDiffInterval(BinomialProportion $a, BinomialProportion $b, float $z = self::Z_95): array + { + $pA = $a->rate(); + $pB = $b->rate(); + [$lA, $uA] = $this->wilsonInterval($a->successes, $a->trials, $z); + [$lB, $uB] = $this->wilsonInterval($b->successes, $b->trials, $z); + + $lower = ($pA - $pB) - sqrt(($pA - $lA) ** 2 + ($uB - $pB) ** 2); + $upper = ($pA - $pB) + sqrt(($uA - $pA) ** 2 + ($pB - $lB) ** 2); + + return [$lower, $upper]; + } + + /** Pooled two-proportion z statistic — descriptive only, not the decision test. */ + public function twoProportionZ(BinomialProportion $a, BinomialProportion $b): float + { + $pooled = ($a->successes + $b->successes) / ($a->trials + $b->trials); + $se = sqrt($pooled * (1 - $pooled) * (1 / $a->trials + 1 / $b->trials)); + + return $se > 0.0 ? ($a->rate() - $b->rate()) / $se : 0.0; + } + + /** + * Two-sided Fisher exact p-value for the 2x2 table [[a, b], [c, d]] + * (a/c = successes, b/d = failures). Sums the hypergeometric probabilities + * of every same-margin table no more likely than the observed one. Exact at + * any sample size — no normal approximation. + */ + public function fisherExactTwoSided(int $a, int $b, int $c, int $d): float + { + $rowA = $a + $b; + $rowB = $c + $d; + $col = $a + $c; + $total = $rowA + $rowB; + + if ($rowA === 0 || $rowB === 0 || $col === 0 || $col === $total) { + return 1.0; + } + + $logProbObserved = $this->hypergeometricLogProb($a, $rowA, $rowB, $col); + $p = 0.0; + for ($x = max(0, $col - $rowB); $x <= min($col, $rowA); $x++) { + $logProb = $this->hypergeometricLogProb($x, $rowA, $rowB, $col); + if ($logProb <= $logProbObserved + 1e-7) { + $p += exp($logProb); + } + } + + return min(1.0, $p); + } + + private function hypergeometricLogProb(int $x, int $rowA, int $rowB, int $col): float + { + return $this->logChoose($rowA, $x) + $this->logChoose($rowB, $col - $x) - $this->logChoose($rowA + $rowB, $col); + } + + private function logChoose(int $n, int $k): float + { + if ($k < 0 || $k > $n) { + return -INF; + } + + return $this->logFactorial($n) - $this->logFactorial($k) - $this->logFactorial($n - $k); + } + + private function logFactorial(int $n): float + { + for ($i = count($this->logFactorials); $i <= $n; $i++) { + $this->logFactorials[$i] = $this->logFactorials[$i - 1] + log($i); + } + + return $this->logFactorials[$n]; + } +} diff --git a/tests/Unit/ProportionSignificanceTest.php b/tests/Unit/ProportionSignificanceTest.php new file mode 100644 index 00000000..77e89701 --- /dev/null +++ b/tests/Unit/ProportionSignificanceTest.php @@ -0,0 +1,66 @@ +stats = new ProportionSignificance; +}); + +it('computes the Wilson score interval', function () { + [$low, $high] = $this->stats->wilsonInterval(6, 52); + expect($low)->toEqualWithDelta(0.0540, 0.0005) + ->and($high)->toEqualWithDelta(0.2297, 0.0005); + + [$low2, $high2] = $this->stats->wilsonInterval(5, 179); + expect($low2)->toEqualWithDelta(0.0120, 0.0005) + ->and($high2)->toEqualWithDelta(0.0640, 0.0005); +}); + +it('keeps the Wilson interval inside [0, 1] at the k=0 and k=n boundaries', function () { + expect($this->stats->wilsonInterval(0, 30)[0])->toBe(0.0) + ->and($this->stats->wilsonInterval(30, 30)[1])->toBe(1.0); +}); + +it('computes a two-sided Fisher exact p-value', function () { + // Real experiment table: reduced 6/52 vs pay_now 5/179. + expect($this->stats->fisherExactTwoSided(6, 46, 5, 174))->toEqualWithDelta(0.0182, 0.0005); + // Clean separation 5/5 vs 0/5 → 2 * C(5,5)/C(10,5) = 2/252. + expect($this->stats->fisherExactTwoSided(5, 0, 0, 5))->toEqualWithDelta(0.007936, 0.00001); + // Degenerate margin (no successes anywhere) → p = 1. + expect($this->stats->fisherExactTwoSided(0, 10, 0, 10))->toBe(1.0); +}); + +it('computes the Newcombe difference interval and the descriptive z', function () { + $reduced = new BinomialProportion('reduced', 6, 52); + $payNow = new BinomialProportion('pay_now', 5, 179); + + [$low, $high] = $this->stats->newcombeDiffInterval($reduced, $payNow); + expect($low)->toEqualWithDelta(0.0164, 0.0005) + ->and($high)->toEqualWithDelta(0.2029, 0.0005) + ->and($this->stats->twoProportionZ($reduced, $payNow))->toEqualWithDelta(2.607, 0.01); +}); + +it('calls the borderline real comparison NOT significant under Bonferroni', function () { + $result = $this->stats->compare( + new BinomialProportion('reduced', 6, 52), + new BinomialProportion('pay_now', 5, 179), + ); + + // Fisher p ≈ 0.018 exceeds the corrected bar 0.05/3 ≈ 0.0167 → not a winner yet. + expect($result['significant'])->toBeFalse() + ->and($result['fisherP'])->toEqualWithDelta(0.0182, 0.0005) + ->and($result['alpha'])->toEqualWithDelta(0.0167, 0.0005) + ->and($result['minExpectedCount'])->toEqualWithDelta(2.48, 0.05); +}); + +it('calls a clean separation significant', function () { + $result = $this->stats->compare( + new BinomialProportion('a', 5, 5), + new BinomialProportion('b', 0, 5), + ); + + // Fisher p ≈ 0.008 clears the corrected bar. + expect($result['significant'])->toBeTrue() + ->and($result['fisherP'])->toEqualWithDelta(0.007936, 0.00001); +});