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