From 70bb14e662b22e89558a9d0489a7d697689e7208 Mon Sep 17 00:00:00 2001 From: Ondra Zahradnik Date: Fri, 12 Jun 2026 09:43:13 +0200 Subject: [PATCH] Keep sequential confidence intervals finite and warn on negative variance estimates The O'Brien-Fleming alpha spending function rounds the adjusted confidence level to 4 decimal places which produces exactly 1.0 early in the experiment. The t-quantile of confidence level 1.0 is infinite, making all confidence intervals inf or NaN. The adjusted confidence level is now capped at 0.9999 so confidence intervals stay finite (and appropriately wide) from day one. Metrics evaluated from pre-aggregated goals can produce negative variance estimates when sum_sqr_value is not a sum of squared per-unit goal totals, e.g. when a dimensional goal sums multiple dimension values per unit but squares were computed per dimension value. This produced silent NaN standard deviations, p-values and confidence intervals. The evaluation now emits an explicit warning naming the affected metrics and the likely cause. Also adds a regression test pinning that evaluate_by_unit computes sum_sqr_value from per-unit totals when a unit is split into multiple rows. Co-Authored-By: Claude Fable 5 --- src/epstats/toolkit/experiment.py | 162 ++++++++++++++---- src/epstats/toolkit/statistics.py | 54 ++++-- .../testing/resources/evaluations_metrics.csv | 4 +- tests/epstats/toolkit/test_experiment.py | 111 +++++++++++- tests/epstats/toolkit/test_statistics.py | 76 ++++++-- 5 files changed, 340 insertions(+), 67 deletions(-) diff --git a/src/epstats/toolkit/experiment.py b/src/epstats/toolkit/experiment.py index 7d0e70d..bb18e13 100644 --- a/src/epstats/toolkit/experiment.py +++ b/src/epstats/toolkit/experiment.py @@ -1,4 +1,5 @@ import logging +import warnings from collections import Counter from dataclasses import dataclass from datetime import datetime @@ -17,7 +18,9 @@ from .statistics import DEFAULT_CONFIDENCE_LEVEL, DEFAULT_POWER, Statistics from .utils import get_utc_timestamp, goals_wide_to_long -check_evaluation_errors_metric = get_prometheus_metric("check_evaluation_errors_total", PrometheusCounter) +check_evaluation_errors_metric = get_prometheus_metric( + "check_evaluation_errors_total", PrometheusCounter +) class Evaluation: @@ -25,7 +28,9 @@ class Evaluation: Results of an experiment evaluation. """ - def __init__(self, metrics: pd.DataFrame, checks: pd.DataFrame, exposures: pd.DataFrame): + def __init__( + self, metrics: pd.DataFrame, checks: pd.DataFrame, exposures: pd.DataFrame + ): self.metrics = metrics self.checks = checks self.exposures = exposures @@ -168,10 +173,20 @@ def __init__( self.metrics = metrics self._check_metric_ids_unique() self.checks = checks - self.date_from = datetime.strptime(date_from, "%Y-%m-%d").date() if date_from is not None else None - self.date_to = datetime.strptime(date_to, "%Y-%m-%d").date() if date_to is not None else None + self.date_from = ( + datetime.strptime(date_from, "%Y-%m-%d").date() + if date_from is not None + else None + ) + self.date_to = ( + datetime.strptime(date_to, "%Y-%m-%d").date() + if date_to is not None + else None + ) self.date_for = ( - datetime.strptime(date_for, "%Y-%m-%d").date() if date_for is not None else datetime.today().date() + datetime.strptime(date_for, "%Y-%m-%d").date() + if date_for is not None + else datetime.today().date() ) self.confidence_level = confidence_level self.variants = variants @@ -201,7 +216,9 @@ def _check_metric_ids_unique(self): id_counts = Counter(metric.id for metric in self.metrics) for id_, count in id_counts.items(): if count > 1: - raise ValueError(f"Metric ids must be unique. Id={id_} found more than once.") + raise ValueError( + f"Metric ids must be unique. Id={id_} found more than once." + ) def _update_dimension_to_value(self): """ @@ -531,7 +548,11 @@ def _exposures_fce_agg(goals: pd.DataFrame, exp_id: str, unit_type: str): Evaluates checks from already aggregated goals. """ checks_df = ( - goals[(goals["unit_type"] == unit_type) & (goals["agg_type"] == "global") & (goals["goal"] == "exposure")] + goals[ + (goals["unit_type"] == unit_type) + & (goals["agg_type"] == "global") + & (goals["goal"] == "exposure") + ] .groupby("exp_variant_id") .agg(exposures=("count", "sum")) .reset_index() @@ -544,16 +565,22 @@ def _exposures_fce_by_unit(goals: pd.DataFrame, exp_id: str, unit_type: str): """ Evaluates checks from already aggregated goals. """ - checks_df = goals[(goals["unit_type"] == unit_type) & (goals["agg_type"] == "unit")][ - [("exp_variant_id", ""), ("exposure", "count")] - ] + checks_df = goals[ + (goals["unit_type"] == unit_type) & (goals["agg_type"] == "unit") + ][[("exp_variant_id", ""), ("exposure", "count")]] checks_df = checks_df.droplevel(0, axis=1) checks_df.columns = ["exp_variant_id", "exposures"] - d = checks_df.groupby("exp_variant_id").agg(exposures=("exposures", "sum")).reset_index() + d = ( + checks_df.groupby("exp_variant_id") + .agg(exposures=("exposures", "sum")) + .reset_index() + ) d["exp_id"] = exp_id return d - def _evaluate(self, goals: pd.DataFrame, metrics_column_fce, checks_fce, exposures_fce): + def _evaluate( + self, goals: pd.DataFrame, metrics_column_fce, checks_fce, exposures_fce + ): metrics = self._evaluate_metrics(goals, metrics_column_fce) checks = self._evaluate_checks(goals, checks_fce) exposures = self._evaluate_exposures(goals, exposures_fce) @@ -570,10 +597,16 @@ def _evaluate_checks(self, goals: pd.DataFrame, check_fce) -> pd.DataFrame: r["exp_id"] = self.id res.append(r) except Exception as e: - self._logger.warning(f"Cannot evaluate check [{c.id} in experiment [{self.id}] because of {e}") + self._logger.warning( + f"Cannot evaluate check [{c.id} in experiment [{self.id}] because of {e}" + ) check_evaluation_errors_metric.inc() - c = pd.concat(res, axis=0) if res != [] else pd.DataFrame([], columns=Evaluation.check_columns()) + c = ( + pd.concat(res, axis=0) + if res != [] + else pd.DataFrame([], columns=Evaluation.check_columns()) + ) c["timestamp"] = round(get_utc_timestamp(datetime.now()).timestamp()) return c[Evaluation.check_columns()] @@ -624,7 +657,9 @@ def _fix_missing_agg(self, goals: pd.DataFrame) -> pd.DataFrame: ) for dimension in self.get_dimension_columns(): - empty_df[dimension] = np.repeat([g.dimension_to_value.get(dimension, "") for g in ngs], lnvs) + empty_df[dimension] = np.repeat( + [g.dimension_to_value.get(dimension, "") for g in ngs], lnvs + ) # join to existing data and use zeros for only missing variants and goals m = ( @@ -677,10 +712,14 @@ def _fix_missing_by_unit(self, goals: pd.DataFrame) -> pd.DataFrame: ) for dimension in self.get_dimension_columns(): - empty_df[dimension] = np.repeat([g.dimension_to_value.get(dimension, "") for g in ngs], lnvs) + empty_df[dimension] = np.repeat( + [g.dimension_to_value.get(dimension, "") for g in ngs], lnvs + ) # join to existing data and use zeros for only missing variants and goals - m = pd.concat([g, empty_df], axis=0).fillna({d: "" for d in self.get_dimension_columns()}) + m = pd.concat([g, empty_df], axis=0).fillna( + {d: "" for d in self.get_dimension_columns()} + ) return m[ [ "exp_id", @@ -711,9 +750,15 @@ def _get_required_sample_size( # to the sample size. In such case we do not evaluate the required sample size. # TODO: add suport for value() denominator metrics, # parser will return an additional column equal to count or count_unique. - sample_size = metric_row["count"] if metric_id not in metrics_with_value_denominator else np.nan + sample_size = ( + metric_row["count"] + if metric_id not in metrics_with_value_denominator + else np.nan + ) - if metric_row["exp_variant_id"] == self.control_variant or pd.isna(minimum_effect): + if metric_row["exp_variant_id"] == self.control_variant or pd.isna( + minimum_effect + ): return pd.Series([np.nan, sample_size, np.nan], index) metric_id = metric_row["metric_id"] @@ -734,7 +779,9 @@ def _get_required_sample_size( index, ) - def _get_required_sample_sizes(self, metrics: pd.DataFrame, n_variants: int) -> pd.DataFrame: + def _get_required_sample_sizes( + self, metrics: pd.DataFrame, n_variants: int + ) -> pd.DataFrame: controls = { r["metric_id"]: {"mean": r["mean"], "std": r["std"]} for _, r in metrics.iterrows() @@ -743,7 +790,9 @@ def _get_required_sample_sizes(self, metrics: pd.DataFrame, n_variants: int) -> minimum_effects = {m.id: m.minimum_effect for m in self.metrics} metrics_with_value_denominator = { - m.id for m in self.metrics if m.denominator.startswith("value(") and not isinstance(m, SimpleMetric) + m.id + for m in self.metrics + if m.denominator.startswith("value(") and not isinstance(m, SimpleMetric) } return metrics.apply( @@ -757,7 +806,9 @@ def _get_required_sample_sizes(self, metrics: pd.DataFrame, n_variants: int) -> axis=1, ) - def _get_power_from_required_sample_sizes(self, metrics: pd.DataFrame, n_variants: int) -> pd.Series: + def _get_power_from_required_sample_sizes( + self, metrics: pd.DataFrame, n_variants: int + ) -> pd.Series: return metrics.apply( lambda metric_row: Statistics.power_from_required_sample_size_per_variant( n_variants=n_variants, @@ -785,6 +836,28 @@ def _get_false_positive_risk(self, metric_row: pd.Series) -> float: def _get_false_positive_risks(self, metrics: pd.DataFrame) -> pd.Series: return metrics.apply(self._get_false_positive_risk, axis=1) + def _warn_on_negative_variance(self, var: np.array) -> None: + """ + Warns about metrics with a negative variance estimate of shape (metrics, variants). + + Negative variance estimate means that `sum_sqr_value` is not a sum of squared + per-unit goal totals. This usually happens when a metric matches multiple + pre-aggregated rows per unit, e.g. when a dimensional goal sums multiple + dimension values per unit and `sum_sqr_value` was squared per dimension value, + not per unit. + """ + affected = np.unique(np.nonzero(var < 0)[0]) + if affected.size: + metric_names = [self.metrics[i].name for i in affected] + warnings.warn( + f"Negative variance estimate of metrics {metric_names}. `sum_sqr_value` of these metrics" + " is not a sum of squared per-unit goal totals, this usually happens when a metric matches" + " multiple pre-aggregated rows per unit, e.g. when a dimensional goal sums multiple dimension" + " values per unit. Standard deviations, p-values and confidence intervals of affected metrics" + " will be NaN. Pre-aggregate `sum_sqr_value` from per-unit goal totals or evaluate using" + " `evaluate_by_unit`." + ) + def _evaluate_metrics(self, goals: pd.DataFrame, column_fce) -> pd.DataFrame: if not self.metrics: return pd.DataFrame([], columns=Evaluation.metric_columns()) @@ -805,26 +878,45 @@ def _evaluate_metrics(self, goals: pd.DataFrame, column_fce) -> pd.DataFrame: # There could be division by zero here which is expected as we return # nan or inf values to the caller. mean = sum_value / count - std = np.sqrt((sum_sqr_value - sum_value * sum_value / count) / (count - 1)) + var = (sum_sqr_value - sum_value * sum_value / count) / (count - 1) + std = np.sqrt(var) + + self._warn_on_negative_variance(var) # sequential testing correction if self.date_from is not None and self.date_to is not None: # Parameters - test_length = (self.date_to - self.date_from).days + 1 # test length in days - actual_day = (self.date_for - self.date_from).days + 1 # day(s) since beginning of the test - actual_day = min(actual_day, test_length) # actual day of evaluation must be in interval [1, test_length] + test_length = ( + self.date_to - self.date_from + ).days + 1 # test length in days + actual_day = ( + self.date_for - self.date_from + ).days + 1 # day(s) since beginning of the test + actual_day = min( + actual_day, test_length + ) # actual day of evaluation must be in interval [1, test_length] # confidence level adjustment - applied when actual_day < test_length (test is still running) - confidence_level = Statistics.obf_alpha_spending_function(self.confidence_level, test_length, actual_day) + confidence_level = Statistics.obf_alpha_spending_function( + self.confidence_level, test_length, actual_day + ) else: confidence_level = self.confidence_level # no change - stats = np.dstack((count, mean, std, sum_value, np.ones(count.shape) * confidence_level)) + stats = np.dstack( + (count, mean, std, sum_value, np.ones(count.shape) * confidence_level) + ) stats = np.dstack( ( - np.repeat([m.id for m in self.metrics], n_variants).reshape(metrics, n_variants, -1), - np.repeat([m.name for m in self.metrics], n_variants).reshape(metrics, n_variants, -1), - np.tile(goals["exp_variant_id"].unique(), metrics).reshape(metrics, n_variants, -1), + np.repeat([m.id for m in self.metrics], n_variants).reshape( + metrics, n_variants, -1 + ), + np.repeat([m.name for m in self.metrics], n_variants).reshape( + metrics, n_variants, -1 + ), + np.tile(goals["exp_variant_id"].unique(), metrics).reshape( + metrics, n_variants, -1 + ), stats, ) ) @@ -836,11 +928,15 @@ def _evaluate_metrics(self, goals: pd.DataFrame, column_fce) -> pd.DataFrame: # multiple variants (comparisons) correction - applied when we have multiple treatment variants if n_variants > 2: - c = Statistics.multiple_comparisons_correction(c, n_variants, metrics, confidence_level) + c = Statistics.multiple_comparisons_correction( + c, n_variants, metrics, confidence_level + ) c["exp_id"] = self.id c["timestamp"] = round(get_utc_timestamp(datetime.now()).timestamp()) - c[["minimum_effect", "sample_size", "required_sample_size"]] = self._get_required_sample_sizes(c, n_variants) + c[["minimum_effect", "sample_size", "required_sample_size"]] = ( + self._get_required_sample_sizes(c, n_variants) + ) c["power"] = self._get_power_from_required_sample_sizes(c, n_variants) c["false_positive_risk"] = self._get_false_positive_risks(c) return c[Evaluation.metric_columns()] diff --git a/src/epstats/toolkit/statistics.py b/src/epstats/toolkit/statistics.py index 0384d06..84e6fb2 100644 --- a/src/epstats/toolkit/statistics.py +++ b/src/epstats/toolkit/statistics.py @@ -9,6 +9,11 @@ DEFAULT_CONFIDENCE_LEVEL = 0.95 DEFAULT_POWER = 0.8 +# Maximum confidence level the sequential analysis adjustment can produce. +# It must stay below 1.0, otherwise t-quantiles and confidence intervals +# computed from the adjusted confidence level become infinite. +MAX_SEQUENTIAL_CONFIDENCE_LEVEL = 0.9999 + class Statistics: """ @@ -105,12 +110,16 @@ def ttest_evaluation(cls, stats: np.array, control_variant: str) -> pd.DataFrame # There could be division by zero here which is expected as we return # nan or inf values to the caller. # np.round() in case of roundoff errors, e.g. f = 9.999999998 => trunc(round(f, 5)) = 10 - f = np.trunc(np.round(num / den, 5)) # (rounded & truncated) degrees of freedom + f = np.trunc( + np.round(num / den, 5) + ) # (rounded & truncated) degrees of freedom # t-quantile with warnings.catch_warnings(): warnings.filterwarnings("ignore", category=RuntimeWarning) - t_quantile = st.t.ppf(conf_level + (1 - conf_level) / 2, f) # right quantile + t_quantile = st.t.ppf( + conf_level + (1 - conf_level) / 2, f + ) # right quantile # relative difference and test statistics with np.errstate(divide="ignore", invalid="ignore"): @@ -120,7 +129,10 @@ def ttest_evaluation(cls, stats: np.array, control_variant: str) -> pd.DataFrame rel_diff = (mean_treat - mean_cont) / np.abs(mean_cont) # standard error for relative difference rel_se = ( - np.sqrt((mean_treat * std_cont) ** 2 / (mean_cont**2 * count_cont) + (std_treat**2 / count_treat)) + np.sqrt( + (mean_treat * std_cont) ** 2 / (mean_cont**2 * count_cont) + + (std_treat**2 / count_treat) + ) / mean_cont ) test_stat = rel_diff / rel_se @@ -203,18 +215,28 @@ def multiple_comparisons_correction( index_to = (m + 1) * n_variants - 1 # p-value adjustment - pvals = df.loc[index_from:index_to, "p_value"].to_list() # select old p-values - adj_pvals = multipletests(pvals=pvals, alpha=alpha, method="holm")[1] # compute adjusted p-values + pvals = df.loc[ + index_from:index_to, "p_value" + ].to_list() # select old p-values + adj_pvals = multipletests(pvals=pvals, alpha=alpha, method="holm")[ + 1 + ] # compute adjusted p-values # confidence interval adjustment # we set ratio to 1 when test_stat is so big that pvals are zero, no reason to update ci adj_ratio = np.nan_to_num(pvals / adj_pvals, nan=1) # adjustment ratio adj_alpha = adj_ratio * alpha # adjusted level alpha - f = df.loc[index_from:index_to, "degrees_of_freedom"].to_list() # degrees of freedom - se = df.loc[index_from:index_to, "standard_error"].to_list() # standard error + f = df.loc[ + index_from:index_to, "degrees_of_freedom" + ].to_list() # degrees of freedom + se = df.loc[ + index_from:index_to, "standard_error" + ].to_list() # standard error - t_quantile = st.t.ppf(np.ones(n_variants - 1) - adj_alpha + adj_alpha / 2, f) # right t-quantile + t_quantile = st.t.ppf( + np.ones(n_variants - 1) - adj_alpha + adj_alpha / 2, f + ) # right t-quantile adj_conf_int = se * t_quantile # adjusted confidence interval # replace (unadjusted) p-values and confidence intervals with new adjusted ones @@ -223,7 +245,9 @@ def multiple_comparisons_correction( return df @classmethod - def obf_alpha_spending_function(cls, confidence_level: int, total_length: int, actual_day: int) -> int: + def obf_alpha_spending_function( + cls, confidence_level: int, total_length: int, actual_day: int + ) -> int: """ [O'Brien-Fleming alpha spending function](https://online.stat.psu.edu/stat509/lesson/9/9.6/). We adjust confidence level in time in experiment. Confidence level in this setting is @@ -236,13 +260,14 @@ def obf_alpha_spending_function(cls, confidence_level: int, total_length: int, a Returns: adjusted confidence level with respect to actual day of the experiment and total - length of the experiment. + length of the experiment, capped at `MAX_SEQUENTIAL_CONFIDENCE_LEVEL` so that + confidence intervals stay finite even early in the experiment. """ alpha = 1 - confidence_level t = actual_day / total_length # t in (0, 1] q = st.norm.ppf(1 - alpha / 2) # quantile of normal distribution alpha_adj = 2 - 2 * st.norm.cdf(q / np.sqrt(t)) - return np.round(1 - alpha_adj, decimals=4) + return min(np.round(1 - alpha_adj, decimals=4), MAX_SEQUENTIAL_CONFIDENCE_LEVEL) @staticmethod def required_sample_size_per_variant( @@ -425,11 +450,14 @@ def power_from_required_sample_size_per_variant( if n_variants < 2 or required_sample_size_per_variant == 0: return np.nan - required_sample_size_ratio = sample_size_per_variant / required_sample_size_per_variant + required_sample_size_ratio = ( + sample_size_per_variant / required_sample_size_per_variant + ) alpha = (1 - required_confidence_level) / (n_variants - 1) return st.norm.cdf( - np.sqrt(required_sample_size_ratio) * (st.norm.ppf(1 - alpha / 2) + st.norm.ppf(required_power)) + np.sqrt(required_sample_size_ratio) + * (st.norm.ppf(1 - alpha / 2) + st.norm.ppf(required_power)) - st.norm.ppf(1 - alpha / 2) ) diff --git a/src/epstats/toolkit/testing/resources/evaluations_metrics.csv b/src/epstats/toolkit/testing/resources/evaluations_metrics.csv index d29777e..f349203 100644 --- a/src/epstats/toolkit/testing/resources/evaluations_metrics.csv +++ b/src/epstats/toolkit/testing/resources/evaluations_metrics.csv @@ -39,8 +39,8 @@ test-srm-negative,b,1,Click-through Rate,5100,1100,1100,0.2156862745,0.015686274 test-srm-negative,c,1,Click-through Rate,5300,900,900,0.1698113208,-0.0301886792,-0.1509433962,0.3755023838,5.86106e-05,0.0352411269,-4.2831603193,10148.5415473186,1.85954e-05,3.71908e-05,0.5000005708,0.9749999715,N/A,0.0790012661,0.95,N/A,5300,N/A,N/A,N/A test-srm-one-variant,a,1,Click-through Rate,0,0,0,N/A,N/A,N/A,N/A,N/A,N/A,N/A,N/A,N/A,N/A,N/A,0.95,N/A,N/A,0.95,N/A,0,N/A,N/A,N/A test-srm-one-variant,b,1,Click-through Rate,5100,1100,1100,0.2156862745,-inf,N/A,0.4113379167,N/A,N/A,N/A,N/A,N/A,N/A,N/A,0.95,N/A,N/A,0.95,N/A,5100,N/A,N/A,N/A -test-sequential-v1,a,1,Average Bookings,2000,44400,2308800,22.2,0.0,0.0,25.727241311,0.6618909455,0.0366471534,0.0,3998.0,1.0,1.0,1.0,1.0,N/A,inf,1.0,N/A,2000,N/A,N/A,N/A -test-sequential-v1,b,1,Average Bookings,2020,35000,1700000,17.3267326733,-4.8732673267,-0.2195165462,23.2730881112,0.5990824184,0.0308725819,-7.1104045326,3970.1185027169,0.0,0.0,1.0,1.0,N/A,inf,1.0,N/A,2020,N/A,N/A,N/A +test-sequential-v1,a,1,Average Bookings,2000,44400,2308800,22.2,0.0,0.0,25.727241311,0.6618909455,0.0366471534,0.0,3998.0,1.0,1.0,1.0,0.9999,N/A,0.1427231163,0.9999,N/A,2000,N/A,N/A,N/A +test-sequential-v1,b,1,Average Bookings,2020,35000,1700000,17.3267326733,-4.8732673267,-0.2195165462,23.2730881112,0.5990824184,0.0308725819,-7.1104045326,3970.1185027169,0.0,0.0,1.0,0.9999,N/A,0.1202347778,0.9999,N/A,2020,N/A,N/A,N/A test-sequential-v2,a,1,Average Bookings,6000,133200,6926400,22.2,0.0,0.0,25.7229523651,0.2205567595,0.0211547166,0.0,11998.0,1.0,1.0,1.0,0.9796082526,N/A,0.0490655917,0.9796082526,N/A,6000,N/A,N/A,N/A test-sequential-v2,b,1,Average Bookings,6060,105000,5100000,17.3267326733,-4.8732673267,-0.2195165462,23.2692467167,0.1996278587,0.0178213388,-12.317623755,11914.3163356925,0.0,0.0,1.0,0.9796082526,N/A,0.0413342977,0.9796082526,N/A,6060,N/A,N/A,N/A test-sequential-v3,a,1,Average Bookings,8000,177600,9235200,22.2,0.0,0.0,25.7224163977,0.1654106763,0.0183201403,0.0,15998.0,1.0,1.0,1.0,0.95,N/A,0.0359095319,0.95,N/A,8000,N/A,N/A,N/A diff --git a/tests/epstats/toolkit/test_experiment.py b/tests/epstats/toolkit/test_experiment.py index 5a09ff7..e29fcd1 100644 --- a/tests/epstats/toolkit/test_experiment.py +++ b/tests/epstats/toolkit/test_experiment.py @@ -1,3 +1,5 @@ +import numpy as np +import pandas as pd import pytest from src.epstats.toolkit.check import SimpleSrmCheck, SrmCheck @@ -74,7 +76,9 @@ def test_update_dimension_to_value(unit_type): def test_binary_valued(dao, metrics, checks, unit_type): - experiment = Experiment("test-conversion", "a", metrics, checks, unit_type=unit_type) + experiment = Experiment( + "test-conversion", "a", metrics, checks, unit_type=unit_type + ) evaluate_experiment_agg(experiment, dao) @@ -151,6 +155,111 @@ def test_real_by_unit(dao, unit_type): evaluate_experiment_by_unit(experiment, dao) +def test_by_unit_sums_unit_rows_before_squaring(unit_type): + """ + A unit can be split into multiple rows of the same goal and dimension value, e.g. + when unit goals are pre-aggregated at a finer granularity than the experiment + dimensions. `evaluate_by_unit` must sum the values of the same unit before + squaring them, otherwise variance estimates can go negative resulting in NaN + standard deviations, p-values and confidence intervals. + """ + goals = pd.DataFrame( + { + "exp_id": "test-by-unit-sums", + "exp_variant_id": ["a", "a", "a", "a", "a", "b", "b", "b", "b", "b"], + "unit_type": unit_type, + "agg_type": "unit", + "unit_id": ["u1", "u1", "u1", "u2", "u2", "u3", "u3", "u3", "u4", "u4"], + "goal": [ + "exposure", + "conversion", + "conversion", + "exposure", + "conversion", + "exposure", + "conversion", + "conversion", + "exposure", + "conversion", + ], + "product": ["", "p_1", "p_1", "", "p_1", "", "p_1", "p_1", "", "p_1"], + "count": 1, + "sum_value": [1, 2, 3, 1, 1, 1, 2, 2, 1, 2], + } + ) + experiment = Experiment( + "test-by-unit-sums", + "a", + [ + Metric( + 1, + "Conversions per Exposure of Product p_1", + "value(test_unit_type.unit.conversion(product=p_1))", + "count(test_unit_type.unit.exposure)", + ), + ], + [], + unit_type=unit_type, + ) + metrics = experiment.evaluate_by_unit(goals).metrics.set_index( + ["metric_id", "exp_variant_id"] + ) + + # per-unit conversion totals are a: (5, 1), b: (4, 2); squaring the individual + # rows instead of the per-unit totals would give a negative variance in variant a + assert metrics.loc[(1, "a"), "mean"] == pytest.approx(3.0) + assert metrics.loc[(1, "a"), "std"] == pytest.approx(np.sqrt(8.0)) + assert metrics.loc[(1, "b"), "mean"] == pytest.approx(3.0) + assert metrics.loc[(1, "b"), "std"] == pytest.approx(np.sqrt(2.0)) + assert np.isfinite(metrics["p_value"].astype(float)).all() + + +def test_agg_warns_on_negative_variance(unit_type): + """ + Pre-aggregated goals can carry `sum_sqr_value` that is not a sum of squared + per-unit totals, e.g. when a dimensional goal is summed over multiple dimension + values per unit but squared per dimension value. The variance estimate is then + negative and the evaluation must warn about it instead of failing silently + with NaN results. + """ + goals = pd.DataFrame( + { + "exp_id": "test-agg-negative-variance", + "exp_variant_id": ["a", "a", "a", "b", "b", "b"], + "unit_type": unit_type, + "agg_type": "unit", + "goal": ["exposure", "conversion", "conversion"] * 2, + "product": ["", "p_1", "p_2"] * 2, + # `sum_sqr_value` of conversion is squared per (unit, product) row, the + # sum of squared per-unit totals would be much higher + "count": [1000, 1000, 1000] * 2, + "sum_sqr_count": [1000, 1000, 1000] * 2, + "sum_value": [1000, 1000, 1000] * 2, + "sum_sqr_value": [1000, 1.5, 1.5] * 2, + "count_unique": [1000, 1000, 1000] * 2, + } + ) + experiment = Experiment( + "test-agg-negative-variance", + "a", + [ + Metric( + 1, + "Conversions per Exposure", + "value(test_unit_type.unit.conversion)", + "count(test_unit_type.unit.exposure)", + ), + ], + [], + unit_type=unit_type, + ) + with pytest.warns(UserWarning, match="Negative variance estimate"): + metrics = experiment.evaluate_agg(goals).metrics.set_index( + ["metric_id", "exp_variant_id"] + ) + assert np.isnan(float(metrics.loc[(1, "a"), "std"])) + + def test_different_control_variant(dao, checks, unit_type): """ This test tests situation when control variant is not the first one as usual. diff --git a/tests/epstats/toolkit/test_statistics.py b/tests/epstats/toolkit/test_statistics.py index 028e6ba..6399b82 100644 --- a/tests/epstats/toolkit/test_statistics.py +++ b/tests/epstats/toolkit/test_statistics.py @@ -1,5 +1,6 @@ import numpy as np import pytest +import scipy.stats as st from statsmodels.stats.power import TTestIndPower from src.epstats.toolkit.statistics import Statistics @@ -24,13 +25,13 @@ def _assert_sample_sizes_equal_within_tolerance(x, y, n_variants): @pytest.mark.parametrize( "test_length, actual_day, expected", [ - (14, 0, 1.00), - (14, 1, 1.00), - (14, 2, 1.00), - (14, 3, 1.00), + (14, 0, 0.9999), + (14, 1, 0.9999), + (14, 2, 0.9999), + (14, 3, 0.9999), (14, 14, 0.95), - (7, 1, 1.00), - (28, 4, 1.00), + (7, 1, 0.9999), + (28, 4, 0.9999), (28, 8, 0.9998), (28, 28, 0.95), ], @@ -40,6 +41,21 @@ def test_obf_alpha_spending_function(test_length, actual_day, expected): assert alpha == expected +@pytest.mark.parametrize("test_length, actual_day", [(14, 1), (28, 4), (365, 1)]) +def test_obf_alpha_spending_function_keeps_confidence_intervals_finite( + test_length, actual_day +): + """ + Early in the experiment, the adjusted confidence level must stay below 1.0, + otherwise t-quantiles and confidence intervals become infinite. + """ + confidence_level = Statistics.obf_alpha_spending_function( + 0.95, test_length, actual_day + ) + assert confidence_level < 1.0 + assert np.isfinite(st.t.ppf(confidence_level + (1 - confidence_level) / 2, 100)) + + @pytest.mark.parametrize( # expected from https://bookingcom.github.io/powercalculator "n_variants, minimum_effect, mean, std, expected", @@ -57,7 +73,9 @@ def test_obf_alpha_spending_function(test_length, actual_day, expected): (4, 0.10, 0.2, 2.0, 208576), ], ) -def test_required_sample_size_per_variant_equal_variance(n_variants, minimum_effect, mean, std, expected): +def test_required_sample_size_per_variant_equal_variance( + n_variants, minimum_effect, mean, std, expected +): sample_size_per_variant = Statistics.required_sample_size_per_variant( n_variants=n_variants, minimum_effect=minimum_effect, @@ -75,8 +93,12 @@ def test_required_sample_size_per_variant_equal_variance(n_variants, minimum_eff nobs1=None, ) - _assert_sample_sizes_equal(sample_size_per_variant, round(expected_from_statsmodels)) - _assert_sample_sizes_equal_within_tolerance(sample_size_per_variant, expected, n_variants) + _assert_sample_sizes_equal( + sample_size_per_variant, round(expected_from_statsmodels) + ) + _assert_sample_sizes_equal_within_tolerance( + sample_size_per_variant, expected, n_variants + ) @pytest.mark.parametrize( @@ -86,7 +108,9 @@ def test_required_sample_size_per_variant_equal_variance(n_variants, minimum_eff (0.05, 0.3, 2.0, 2.5), ], ) -def test_required_sample_size_per_variant_unequal_variance(minimum_effect, mean, std, std_2): +def test_required_sample_size_per_variant_unequal_variance( + minimum_effect, mean, std, std_2 +): sample_size_per_variant = Statistics.required_sample_size_per_variant( n_variants=2, minimum_effect=minimum_effect, @@ -109,7 +133,9 @@ def test_required_sample_size_per_variant_unequal_variance(minimum_effect, mean, nobs1=None, ) - _assert_sample_sizes_equal(sample_size_per_variant, round(expected_from_statsmodels)) + _assert_sample_sizes_equal( + sample_size_per_variant, round(expected_from_statsmodels) + ) @pytest.mark.parametrize( @@ -122,7 +148,9 @@ def test_required_sample_size_per_variant_unequal_variance(minimum_effect, mean, (4, 0.10, 0.1, None, 19596), ], ) -def test_required_sample_size_per_variant_bernoulli(n_variants, minimum_effect, mean, std, expected): +def test_required_sample_size_per_variant_bernoulli( + n_variants, minimum_effect, mean, std, expected +): sample_size_per_variant = Statistics.required_sample_size_per_variant_bernoulli( n_variants=n_variants, minimum_effect=minimum_effect, @@ -145,8 +173,12 @@ def test_required_sample_size_per_variant_bernoulli(n_variants, minimum_effect, nobs1=None, ) - _assert_sample_sizes_equal(sample_size_per_variant, round(expected_from_statsmodels)) - _assert_sample_sizes_equal_within_tolerance(sample_size_per_variant, expected, n_variants) + _assert_sample_sizes_equal( + sample_size_per_variant, round(expected_from_statsmodels) + ) + _assert_sample_sizes_equal_within_tolerance( + sample_size_per_variant, expected, n_variants + ) @pytest.mark.parametrize( @@ -157,7 +189,9 @@ def test_required_sample_size_per_variant_bernoulli(n_variants, minimum_effect, (2, 0.1, 10.1, None, Statistics.required_sample_size_per_variant_bernoulli), ], ) -def test_required_sample_size_per_variant_raises_exception(n_variants, minimum_effect, mean, std, f): +def test_required_sample_size_per_variant_raises_exception( + n_variants, minimum_effect, mean, std, f +): args = {"minimum_effect": minimum_effect, "mean": mean, "n_variants": n_variants} if std is not None: @@ -177,7 +211,9 @@ def test_required_sample_size_per_variant_raises_exception(n_variants, minimum_e (np.nan, np.nan, np.nan, np.isnan), ], ) -def test_required_sample_size_per_variant_not_valid(minimum_effect, mean, std, expected): +def test_required_sample_size_per_variant_not_valid( + minimum_effect, mean, std, expected +): assert expected( Statistics.required_sample_size_per_variant( minimum_effect=minimum_effect, @@ -198,7 +234,9 @@ def test_required_sample_size_per_variant_not_valid(minimum_effect, mean, std, e (4, 300000), ], ) -def test_power_from_required_sample_size_per_variant(n_variants, sample_size_per_variant): +def test_power_from_required_sample_size_per_variant( + n_variants, sample_size_per_variant +): mean = 0.2 std = 2.0 minimum_effect = 0.05 @@ -256,6 +294,8 @@ def test_power_from_required_sample_size_per_variant_is_nan(args): def test_false_positive_risk(): - false_positive_risk = Statistics.false_positive_risk(null_hypothesis_rate=1 - 0.33, power=0.8, p_value=0.025) + false_positive_risk = Statistics.false_positive_risk( + null_hypothesis_rate=1 - 0.33, power=0.8, p_value=0.025 + ) assert false_positive_risk == 0.05966162065894922