ranklens.stats

ranklens.stats

Statistics of comparing two runs: is the difference real or is it noise?

Which hypothesis. Everything here is about one number: the mean over queries of the per-query difference of one metric at one cutoff (for example NDCG@10), B minus A.

• permutation_test tests H0: the per-query differences are symmetric around zero, that is, for every query it makes no difference which run is called A and which B. Rejecting H0 means such a mean difference is unlikely to come from chance alone.
• paired_bootstrap tests nothing: it says how precisely that mean is measured if the queries at hand are a sample of a larger population of queries.
• benjamini_hochberg keeps many comparisons honest: over 40 segments a couple of "significant" ones appear by chance, so reports show q-values.
• ppi_compare compares on cheap (proxy) judgements corrected by a gold sample: a systematic error of the cheap judgements would bias the difference, and no resampling of queries could see it.
• minimum_detectable_effect answers the question to ask first: is this many queries enough to see the effect we care about?
• The unit of observation is the query. The inference is about the population of queries, not about documents, sessions or users.

What these do not answer. Not "do the two rankings differ": two runs can put documents in completely different orders and still have the same NDCG — for the orders themselves use ranklens.metrics.rbo. And not "will users notice": that is an online experiment.

MIN_QUERIES module-attribute

MIN_QUERIES = 20

Below this the interval is so wide that it says little; a warning, not an error.

mde module-attribute

mde = minimum_detectable_effect

Short alias of minimum_detectable_effect.

MIN_GOLD module-attribute

MIN_GOLD = 30

Fewer gold queries give a correction too noisy to rely on; a warning is emitted.

paired_bootstrap

paired_bootstrap(
    deltas: Sequence[float] | Deltas,
    *,
    n_resamples: int = 10000,
    alpha: float = 0.05,
    seed: int = 0,
) -> BootstrapInterval

Percentile confidence interval of the mean difference.

Queries are resampled with replacement n_resamples times; the interval is the alpha/2 and 1 - alpha/2 quantiles of the resampled means. Resampling happens in memory-bounded chunks (see ranklens.stats.chunking).

Source code in src/ranklens/stats/bootstrap.py

def paired_bootstrap(
    deltas: Sequence[float] | Deltas,
    *,
    n_resamples: int = 10_000,
    alpha: float = 0.05,
    seed: int = 0,
) -> BootstrapInterval:
    """Percentile confidence interval of the mean difference.

    Queries are resampled with replacement ``n_resamples`` times; the interval is the
    ``alpha/2`` and ``1 - alpha/2`` quantiles of the resampled means. Resampling happens
    in memory-bounded chunks (see `ranklens.stats.chunking`).
    """
    values = np.asarray(deltas, dtype=np.float64)
    if values.ndim != 1:
        raise ValueError(f"deltas must be one-dimensional, got shape {values.shape}")
    if not 0.0 < alpha < 1.0:
        raise ValueError(f"alpha must be in (0, 1), got {alpha}")
    if n_resamples < 1:
        raise ValueError(f"n_resamples must be >= 1, got {n_resamples}")
    if len(values) == 0:
        raise ValueError("deltas must not be empty")

    means = resample_means(values, n_resamples=n_resamples, seed=seed)
    low, high = np.quantile(means, [alpha / 2, 1.0 - alpha / 2])
    return BootstrapInterval(
        delta=float(values.mean()),
        low=float(low),
        high=float(high),
        alpha=alpha,
        n_queries=len(values),
        n_resamples=n_resamples,
        seed=seed,
    )

paired_deltas

paired_deltas(
    per_query_a: Mapping[QueryId, float],
    per_query_b: Mapping[QueryId, float],
    *,
    min_queries: int = MIN_QUERIES,
) -> tuple[tuple[QueryId, ...], Deltas]

Queries measured in both runs, sorted, and their per-query differences b - a.

Queries present in only one of the runs are dropped: a paired test needs pairs. With no common queries at all there is nothing to compare — that is InsufficientSampleError. With fewer than min_queries the numbers are computable but unreliable, so they are produced with a SmallSampleWarning.

Source code in src/ranklens/stats/bootstrap.py

def paired_deltas(
    per_query_a: Mapping[QueryId, float],
    per_query_b: Mapping[QueryId, float],
    *,
    min_queries: int = MIN_QUERIES,
) -> tuple[tuple[QueryId, ...], Deltas]:
    """Queries measured in both runs, sorted, and their per-query differences ``b - a``.

    Queries present in only one of the runs are dropped: a paired test needs pairs.
    With no common queries at all there is nothing to compare — that is
    `InsufficientSampleError`. With fewer than ``min_queries`` the numbers are
    computable but unreliable, so they are produced with a `SmallSampleWarning`.
    """
    common = tuple(sorted(per_query_a.keys() & per_query_b.keys()))
    if not common:
        raise InsufficientSampleError(0, min_queries)
    if len(common) < min_queries:
        warnings.warn(
            f"comparison on {len(common)} paired queries: the interval will be wide and the "
            f"p-value unstable; {min_queries} or more are recommended (see the MDE calculator)",
            SmallSampleWarning,
            stacklevel=2,
        )
    deltas = np.fromiter(
        (per_query_b[query] - per_query_a[query] for query in common),
        dtype=np.float64,
        count=len(common),
    )
    return common, deltas

adjust

adjust(
    results: Sequence[Comparison],
) -> tuple[Comparison, ...]

Same comparisons with q_value filled in; significant then uses it.

Works for ComparisonResult and PPIComparison alike: BH needs only p-values.

Source code in src/ranklens/stats/multiple.py

def adjust(results: Sequence[Comparison]) -> tuple[Comparison, ...]:
    """Same comparisons with ``q_value`` filled in; ``significant`` then uses it.

    Works for `ComparisonResult` and `PPIComparison` alike: BH needs only p-values.
    """
    q_values = benjamini_hochberg([result.p_value for result in results])
    return tuple(
        dataclasses.replace(result, q_value=q_value)
        for result, q_value in zip(results, q_values, strict=True)
    )

benjamini_hochberg

benjamini_hochberg(
    p_values: Sequence[float] | NDArray[float64],
) -> tuple[float, ...]

q-values of the Benjamini-Hochberg procedure, in the order of the input.

For p-values sorted ascending, q(i) = min over j >= i of (m / j) * p(j), capped at 1. A comparison is a discovery at level alpha when its q-value is below alpha.

Source code in src/ranklens/stats/multiple.py

def benjamini_hochberg(p_values: Sequence[float] | npt.NDArray[np.float64]) -> tuple[float, ...]:
    """q-values of the Benjamini-Hochberg procedure, in the order of the input.

    For p-values sorted ascending, q(i) = min over j >= i of (m / j) * p(j), capped at 1.
    A comparison is a discovery at level alpha when its q-value is below alpha.
    """
    values = np.asarray(p_values, dtype=np.float64)
    if values.ndim != 1:
        raise ValueError(f"p_values must be one-dimensional, got shape {values.shape}")
    if values.size == 0:
        return ()
    if not np.all((values >= 0.0) & (values <= 1.0)):
        raise ValueError("p_values must all lie in [0, 1]")

    m = values.size
    ascending = np.argsort(values, kind="stable")
    ranks = np.arange(m, 0, -1)  # from m down to 1, matching the reversed order
    # walk from the largest p-value down, keeping the running minimum
    adjusted = np.minimum.accumulate(values[ascending][::-1] * m / ranks)[::-1]
    result = np.empty(m, dtype=np.float64)
    result[ascending] = np.minimum(adjusted, 1.0)
    return tuple(float(value) for value in result)

permutation_test

permutation_test(
    deltas: Sequence[float] | Deltas,
    *,
    n_permutations: int = 10000,
    seed: int = 0,
) -> float

Two-sided p-value against H0: the per-query differences are symmetric around zero.

In other words, under H0 it makes no difference for any query which run is called A and which B, so the sign of each difference is arbitrary. The test flips signs at random and counts how often the mean difference comes out at least as extreme as the observed one. The p-value is (b + 1) / (m + 1): the observed assignment belongs to the reference set, so it is never 0 — with 10k permutations the smallest reportable value is about 1e-4.

A small p-value says the mean difference of this metric is unlikely to be noise of the query sample. It says nothing about how differently the two runs order documents (see ranklens.metrics.rbo) or about the effect on users.

Source code in src/ranklens/stats/permutation.py

def permutation_test(
    deltas: Sequence[float] | Deltas,
    *,
    n_permutations: int = 10_000,
    seed: int = 0,
) -> float:
    """Two-sided p-value against H0: the per-query differences are symmetric around zero.

    In other words, under H0 it makes no difference for any query which run is called A
    and which B, so the sign of each difference is arbitrary. The test flips signs at
    random and counts how often the mean difference comes out at least as extreme as the
    observed one. The p-value is ``(b + 1) / (m + 1)``: the observed assignment belongs to
    the reference set, so it is never 0 — with 10k permutations the smallest reportable
    value is about 1e-4.

    A small p-value says the mean difference of **this metric** is unlikely to be noise
    of the query sample. It says nothing about how differently the two runs order
    documents (see `ranklens.metrics.rbo`) or about the effect on users.
    """
    values = np.asarray(deltas, dtype=np.float64)
    if values.ndim != 1:
        raise ValueError(f"deltas must be one-dimensional, got shape {values.shape}")
    if n_permutations < 1:
        raise ValueError(f"n_permutations must be >= 1, got {n_permutations}")
    if len(values) < 2:
        raise InsufficientSampleError(len(values), 2)

    rng = np.random.default_rng(seed)
    observed = abs(float(values.mean()))
    n = len(values)
    at_least_as_extreme = 0
    for size in chunk_sizes(n_permutations, n):
        flipped = np.where(rng.random(size=(size, n)) < 0.5, -values, values)
        means = np.abs(flipped.mean(axis=1))
        at_least_as_extreme += int((means >= observed - _TIES_TOLERANCE).sum())
    return (at_least_as_extreme + 1) / (n_permutations + 1)

minimum_detectable_effect

minimum_detectable_effect(
    sd: float,
    n_queries: int,
    *,
    alpha: float = 0.05,
    power: float = 0.8,
) -> float

Smallest true difference that would be detected with probability power.

MDE = (z(1 - alpha/2) + z(power)) * sd / sqrt(n), where sd is the standard deviation of the per-query differences. A smaller difference is not "absent": it is simply below the resolution of this many queries.

Source code in src/ranklens/stats/power.py

def minimum_detectable_effect(
    sd: float, n_queries: int, *, alpha: float = 0.05, power: float = 0.8
) -> float:
    """Smallest true difference that would be detected with probability ``power``.

    MDE = (z(1 - alpha/2) + z(power)) * sd / sqrt(n), where ``sd`` is the standard
    deviation of the per-query differences. A smaller difference is not "absent":
    it is simply below the resolution of this many queries.
    """
    _check(alpha, power)
    if sd < 0:
        raise ValueError(f"sd must be >= 0, got {sd}")
    if n_queries < 1:
        raise ValueError(f"n_queries must be >= 1, got {n_queries}")
    return float(_factor(alpha, power) * sd / np.sqrt(n_queries))

required_queries

required_queries(
    delta: float,
    sd: float,
    *,
    alpha: float = 0.05,
    power: float = 0.8,
) -> int

How many paired queries are needed to detect a difference of delta.

Never fewer than 2: below that no test can tell a difference from noise, even when sd is 0.

Source code in src/ranklens/stats/power.py

def required_queries(delta: float, sd: float, *, alpha: float = 0.05, power: float = 0.8) -> int:
    """How many paired queries are needed to detect a difference of ``delta``.

    Never fewer than 2: below that no test can tell a difference from noise, even when
    ``sd`` is 0.
    """
    _check(alpha, power)
    if delta <= 0:
        raise ValueError(f"delta must be > 0, got {delta}")
    if sd < 0:
        raise ValueError(f"sd must be >= 0, got {sd}")
    return max(2, int(np.ceil((_factor(alpha, power) * sd / delta) ** 2)))

standard_deviation

standard_deviation(deltas: Sequence[float]) -> float

Standard deviation of per-query differences, as the two functions above expect.

Source code in src/ranklens/stats/power.py

def standard_deviation(deltas: Sequence[float]) -> float:
    """Standard deviation of per-query differences, as the two functions above expect."""
    values = np.asarray(deltas, dtype=np.float64)
    if values.size < 2:
        raise ValueError("at least two differences are needed to estimate a spread")
    return float(values.std(ddof=1))

ppi_compare

ppi_compare(
    proxy_a: MetricResult,
    proxy_b: MetricResult,
    gold_a: MetricResult,
    gold_b: MetricResult,
    *,
    lam: float | None = None,
    alpha: float = 0.05,
) -> PPIComparison

Difference B - A of one metric, measured on proxy judgements and corrected on gold.

proxy_* are the metric of each run under proxy judgements (all queries), gold_* under gold judgements (a random sample of queries). Every gold query must also be measured under proxy judgements. lam fixes the weight of the proxy part: 0 is gold only, 1 is classic PPI; None estimates it (PPI++).

With fewer than two proxy-only queries the spread of the proxy part cannot be measured, so the weight is 0 whatever lam says: an interval that ignores that spread would be too narrow. Queries measured under proxy judgements in only one of the runs are not used.

Source code in src/ranklens/stats/ppi.py

def ppi_compare(
    proxy_a: MetricResult,
    proxy_b: MetricResult,
    gold_a: MetricResult,
    gold_b: MetricResult,
    *,
    lam: float | None = None,
    alpha: float = 0.05,
) -> PPIComparison:
    """Difference B - A of one metric, measured on proxy judgements and corrected on gold.

    ``proxy_*`` are the metric of each run under proxy judgements (all queries),
    ``gold_*`` under gold judgements (a random sample of queries). Every gold query
    must also be measured under proxy judgements. ``lam`` fixes the weight of the
    proxy part: 0 is gold only, 1 is classic PPI; None estimates it (PPI++).

    With fewer than two proxy-only queries the spread of the proxy part cannot be
    measured, so the weight is 0 whatever ``lam`` says: an interval that ignores that
    spread would be too narrow. Queries measured under proxy judgements in only one of
    the runs are not used.
    """
    if not 0.0 < alpha < 1.0:
        raise ValueError(f"alpha must be in (0, 1), got {alpha}")
    metrics = {proxy_a.metric, proxy_b.metric, gold_a.metric, gold_b.metric}
    if len(metrics) != 1:
        raise ValueError(f"all four results must be of one metric, got {sorted(metrics)}")

    proxy = _deltas(proxy_a.per_query, proxy_b.per_query)
    gold = _deltas(gold_a.per_query, gold_b.per_query)
    missing = gold.keys() - proxy.keys()
    if missing:
        raise ValueError(
            f"{len(missing)} gold queries have no proxy measurement, e.g. {min(missing)!r}: "
            "the correction needs both judgements on the same query"
        )
    if len(gold) < 2:
        raise InsufficientSampleError(len(gold), 2, unit="gold queries")
    if len(gold) < MIN_GOLD:
        warnings.warn(
            f"{len(gold)} gold queries: the correction is too noisy to rely on; "
            f"{MIN_GOLD} or more are recommended",
            SmallSampleWarning,
            stacklevel=2,
        )

    queries_g = sorted(gold)
    queries_u = sorted(proxy.keys() - gold.keys())
    g, p_g = _array(gold, queries_g), _array(proxy, queries_g)
    p_u = _array(proxy, queries_u)
    weight = _optimal_lambda(g, p_g, p_u) if lam is None else float(lam)
    if len(queries_u) < 2:
        weight = 0.0  # no proxy-only spread to weigh against the gold one

    delta, variance = _estimate(g, p_g, p_u, weight)
    half = _NORMAL.inv_cdf(1.0 - alpha / 2.0) * math.sqrt(variance)
    if variance > 0:
        p_value = 2.0 * (1.0 - _NORMAL.cdf(abs(delta) / math.sqrt(variance)))
    else:  # no spread at all: the difference is exactly what was measured
        p_value = 1.0 if delta == 0 else 0.0
    mean_a = _estimate(
        _array(gold_a.per_query, queries_g),
        _array(proxy_a.per_query, queries_g),
        _array(proxy_a.per_query, queries_u),
        weight,
    )[0]
    return PPIComparison(
        metric=proxy_a.metric,
        mean_a=mean_a,
        mean_b=mean_a + delta,
        delta=delta,
        ci_low=delta - half,
        ci_high=delta + half,
        alpha=alpha,
        p_value=min(1.0, p_value),
        lam=weight,
        n_proxy=len(queries_u),
        n_gold=len(queries_g),
    )