paired_stratification .ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "22da6d5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import namedtuple\n",
    "import scipy.stats as sps\n",
    "import statsmodels.stats.api as sms\n",
    "from tqdm.notebook import tqdm as tqdm_notebook\n",
    "from collections import defaultdict\n",
    "from statsmodels.stats.proportion import proportion_confint\n",
    "import numpy as np\n",
    "import itertools\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set(font_scale=1.5, palette='Set2')\n",
    "ExperimentComparisonResults = namedtuple('ExperimentComparisonResults', \n",
    "                                        ['pvalue', 'effect', 'ci_length', 'left_bound', 'right_bound'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48534773",
   "metadata": {},
   "source": [
    "# Примеры спаренных критериев "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8894020",
   "metadata": {},
   "source": [
    "**Absolute t-test**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e2eff740",
   "metadata": {},
   "outputs": [],
   "source": [
    "def absolute_ttest(control, test):\n",
    "    mean_control = np.mean(control)\n",
    "    mean_test = np.mean(test)\n",
    "    var_mean_control  = np.var(control) / len(control)\n",
    "    var_mean_test  = np.var(test) / len(test)\n",
    "    \n",
    "    difference_mean = mean_test - mean_control\n",
    "    difference_mean_var = var_mean_control + var_mean_test\n",
    "    difference_distribution = sps.norm(loc=difference_mean, scale=np.sqrt(difference_mean_var))\n",
    "\n",
    "    left_bound, right_bound = difference_distribution.ppf([0.025, 0.975])\n",
    "    ci_length = (right_bound - left_bound)\n",
    "    pvalue = 2 * min(difference_distribution.cdf(0), difference_distribution.sf(0))\n",
    "    effect = difference_mean\n",
    "    return ExperimentComparisonResults(pvalue, effect, ci_length, left_bound, right_bound)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "65ebef3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def paired_ttest(control, test):\n",
    "    mean_control = np.mean(control)\n",
    "    mean_test = np.mean(test)\n",
    "    \n",
    "    difference_mean = mean_test - mean_control\n",
    "    difference_mean_var = np.var(test - control) / len(control)\n",
    "    difference_distribution = sps.norm(loc=difference_mean, scale=np.sqrt(difference_mean_var))\n",
    "\n",
    "    left_bound, right_bound = difference_distribution.ppf([0.025, 0.975])\n",
    "    ci_length = (right_bound - left_bound)\n",
    "    pvalue = 2 * min(difference_distribution.cdf(0), difference_distribution.sf(0))\n",
    "    effect = difference_mean\n",
    "    return ExperimentComparisonResults(pvalue, effect, ci_length, left_bound, right_bound)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "318de56f",
   "metadata": {},
   "source": [
    "**Relative t-test**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b7dfd9db",
   "metadata": {},
   "outputs": [],
   "source": [
    "def paired_relative_ttest(control, test):\n",
    "    mean_control = np.mean(control)\n",
    "    mean_test = np.mean(test)\n",
    "    \n",
    "    var_mean_control = np.var(control) / len(control)\n",
    "    difference_mean = mean_test - mean_control\n",
    "    difference_mean_var  = np.var(test - control) / len(test)\n",
    "    \n",
    "    covariance = -np.cov(test - control, control)[0, 1] / len(control)\n",
    "\n",
    "    relative_mu = difference_mean / mean_control\n",
    "    relative_var = difference_mean_var / (mean_control ** 2) \\\n",
    "                    + var_mean_control * ((difference_mean ** 2) / (mean_control ** 4))\\\n",
    "                    - 2 * (difference_mean / (mean_control ** 3)) * covariance\n",
    "    relative_distribution = sps.norm(loc=relative_mu, scale=np.sqrt(relative_var))\n",
    "    left_bound, right_bound = relative_distribution.ppf([0.025, 0.975])\n",
    "    \n",
    "    ci_length = (right_bound - left_bound)\n",
    "    pvalue = 2 * min(relative_distribution.cdf(0), relative_distribution.sf(0))\n",
    "    effect = relative_mu\n",
    "    return ExperimentComparisonResults(pvalue, effect, ci_length, left_bound, right_bound)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bfe47cd",
   "metadata": {},
   "source": [
    "**Relative cuped**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7c689ada",
   "metadata": {},
   "outputs": [],
   "source": [
    "def paired_relative_cuped(control, test, control_before, test_before):\n",
    "    theta = np.cov(test - control, test_before - control_before)[0, 1] /\\\n",
    "                np.var(test_before - control_before)\n",
    "\n",
    "    control_cup = control - theta * control_before\n",
    "    test_cup = test - theta * test_before\n",
    "\n",
    "    mean_den = np.mean(control)\n",
    "    mean_num = np.mean(test_cup - control_cup)\n",
    "    var_mean_den  = np.var(control) / len(control)\n",
    "    var_mean_num  = np.var(test_cup - control_cup) / len(control_cup)\n",
    "\n",
    "    cov = np.cov(test_cup - control_cup, control)[0, 1] / len(control)\n",
    "\n",
    "    relative_mu = mean_num / mean_den\n",
    "    relative_var = var_mean_num / (mean_den ** 2)  + var_mean_den * ((mean_num ** 2) / (mean_den ** 4))\\\n",
    "                - 2 * (mean_num / (mean_den ** 3)) * cov\n",
    "    \n",
    "    relative_distribution = sps.norm(loc=relative_mu, scale=np.sqrt(relative_var))\n",
    "    left_bound, right_bound = relative_distribution.ppf([0.025, 0.975])\n",
    "    \n",
    "    ci_length = (right_bound - left_bound)\n",
    "    pvalue = 2 * min(relative_distribution.cdf(0), relative_distribution.sf(0))\n",
    "    effect = relative_mu\n",
    "    return ExperimentComparisonResults(pvalue, effect, ci_length, left_bound, right_bound)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e469be2c",
   "metadata": {},
   "source": [
    "**Парная пост-нормировка**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d0c993ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "def paired_post_normed_bootstrap(control, test, control_before, test_before, test_type='absolute'):\n",
    "    # Функция от средних, которую надо посчитать на каждой выборке\n",
    "    absolute_func = lambda C, T, C_b, T_b: T - (T_b / C_b) * C\n",
    "    relative_func = lambda C, T, C_b, T_b: (T / C) / (T_b / C_b) - 1\n",
    "    \n",
    "    boot_func = absolute_func if test_type == 'absolute' else relative_func\n",
    "    stat_sample = []\n",
    "    \n",
    "    batch_sz = 100\n",
    "    \n",
    "    #В теории boot_samples_size стоить брать не меньше размера выборки. Но на практике можно и меньше.\n",
    "    boot_samples_size = len(control)\n",
    "    for i in range(0, boot_samples_size, batch_sz):\n",
    "        N = len(control)\n",
    "\n",
    "        # Надо помнить, что мы семплируем именно юзеров\n",
    "        # Поэтому, если мы взяли n раз i элемент в выборке control\n",
    "        # То надо столько же раз взять i элемент в выборке control_before\n",
    "        # Поэтому будем семплировать индексы\n",
    "        indices = np.arange(N)\n",
    "\n",
    "        indices_sample = np.random.choice(control_indices, size=(len(control), batch_sz), replace=True)\n",
    "\n",
    "        \n",
    "        C   = np.mean(control[indices_sample], axis=0)\n",
    "        T   = np.mean(test[indices_sample], axis=0)\n",
    "        C_b = np.mean(control_before[indices_sample], axis=0)\n",
    "        T_b = np.mean(test_before[indices_sample], axis=0)\n",
    "        assert len(T) == batch_sz\n",
    "        stat_sample += list(boot_func(C, T, C_b, T_b))\n",
    "\n",
    "    stat_sample = np.array(stat_sample)\n",
    "    # считаем истинный эффект\n",
    "    effect = boot_func(np.mean(control), np.mean(test), np.mean(control_before), np.mean(test_before))\n",
    "    left_bound, right_bound = np.quantile(stat_sample, [0.025, 0.975])\n",
    "    \n",
    "    ci_length = (right_bound - left_bound)\n",
    "    # pvalue - процент статистик, которые лежат левее или правее 0.\n",
    "    pvalue = 2 * min(np.mean(stat_sample > 0), np.mean(stat_sample < 0))\n",
    "    return ExperimentComparisonResults(pvalue, effect, ci_length, left_bound, right_bound)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "755fdd75",
   "metadata": {},
   "source": [
    "# Алгоритм сплитования в парной стратификации"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b66e0a0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def splitter(before_metrics):\n",
    "    size = len(before_metrics)\n",
    "    \n",
    "    # отсортируем массив \n",
    "    sorted_array = np.sort(before_metrics)[::-1]\n",
    "    control = []\n",
    "    test = []\n",
    "    for i in range(0, size, 2):\n",
    "        if np.random.rand() < 0.5:\n",
    "            control.append(sorted_array[i])\n",
    "            test.append(sorted_array[i + 1])\n",
    "        else:\n",
    "            control.append(sorted_array[i + 1])\n",
    "            test.append(sorted_array[i])\n",
    "\n",
    "    return np.array(control), np.array(test)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb2c7162",
   "metadata": {},
   "source": [
    "# сравнение мощностей\n",
    "\n",
    "### парная стратификация применена"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "05e4c2a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def relative_ttest(control, test):\n",
    "    mean_control = np.mean(control)\n",
    "    mean_test = np.mean(test)\n",
    "    var_mean_control  = np.var(control) / len(control)\n",
    "    var_mean_test  = np.var(test) / len(test)\n",
    "\n",
    "    difference_mean = mean_test - mean_control\n",
    "    difference_mean_var = var_mean_control + var_mean_test\n",
    "    difference_distribution = sps.norm(loc=difference_mean, scale=np.sqrt(difference_mean_var))\n",
    "\n",
    "    left_bound, right_bound = difference_distribution.ppf([0.025, 0.975])\n",
    "    left_bound = left_bound / np.mean(control)   # Деление на среднее\n",
    "    right_bound = right_bound / np.mean(control) # Деление на среднее\n",
    "\n",
    "    ci_length = (right_bound - left_bound)\n",
    "    pvalue = 2 * min(difference_distribution.cdf(0), difference_distribution.sf(0))\n",
    "    effect = difference_mean\n",
    "    return ExperimentComparisonResults(pvalue, effect, ci_length, left_bound, right_bound)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b2b22eaa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "08b23f900838449aaa7e69c2b10e7a4d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/30000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Реальный уровень значимости для обычного критрия: 0.0; доверительный интервал: [0.0, 0.00013]\n",
      "Реальный уровень значимости для спаренного критерия: 0.0513; доверительный интервал: [0.04886, 0.05385]\n",
      "Мощность спаренного критрерия VS мощность обычного критерия: 0.8623 VS. 0.0\n"
     ]
    }
   ],
   "source": [
    "bad_cnt_paired = 0\n",
    "paired_power = 0\n",
    "bad_cnt = 0\n",
    "power = 0\n",
    "N = 30000\n",
    "for i in tqdm_notebook(range(N)):\n",
    "    before = sps.expon(scale=1000).rvs(4000)\n",
    "    C_b, T_b = splitter(before)\n",
    "\n",
    "    C = C_b + sps.norm(loc=0, scale=100).rvs(2000)\n",
    "    T = T_b + sps.norm(loc=0, scale=100).rvs(2000)\n",
    "    T *= 1.01\n",
    "\n",
    "    _, _, _, left_bound, right_bound = relative_ttest(C, T)\n",
    "    _, _, _, left_bound_paired, right_bound_paired = paired_relative_ttest(C, T)\n",
    "\n",
    "    if left_bound > 0.01 or right_bound < 0.01:\n",
    "        bad_cnt += 1\n",
    "        \n",
    "    if left_bound_paired > 0.01 or right_bound_paired < 0.01:\n",
    "        bad_cnt_paired += 1\n",
    "\n",
    "    if left_bound > 0:\n",
    "        power += 1\n",
    "\n",
    "    if left_bound_paired > 0:\n",
    "        paired_power += 1\n",
    "    \n",
    "left_real_level, right_real_level = proportion_confint(count = bad_cnt, nobs = N, \n",
    "                                                         alpha=0.05, method='wilson')\n",
    "left_real_level_paired, right_real_level_paired = proportion_confint(count = bad_cnt_paired, nobs = N, \n",
    "                                                         alpha=0.05, method='wilson')\n",
    "print(f\"Реальный уровень значимости для обычного критрия: {round(bad_cnt / N, 4)};\"\\\n",
    "      f\" доверительный интервал: [{round(left_real_level, 5)}, {round(right_real_level, 5)}]\")\n",
    "print(f\"Реальный уровень значимости для спаренного критерия: {round(bad_cnt_paired / N, 4)};\"\\\n",
    "      f\" доверительный интервал: [{round(left_real_level_paired, 5)}, {round(right_real_level_paired, 5)}]\")\n",
    "print(f\"Мощность спаренного критрерия VS мощность обычного критерия: {paired_power / N} VS. {power / N}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9dc2431",
   "metadata": {},
   "source": [
    "### стратификации не было"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "22c67489",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9b2d16768d8c4ad4acdec8bc3a39e6c7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/30000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Реальный уровень значимости для обычного критрия: 0.0505; доверительный интервал: [0.04808, 0.05304]\n",
      "Реальный уровень значимости для спаренного критерия: 0.0507; доверительный интервал: [0.04824, 0.05321]\n",
      "Мощность спаренного критрерия VS мощность обычного критерия: 0.05503333333333333 VS. 0.0482\n"
     ]
    }
   ],
   "source": [
    "bad_cnt_paired = 0\n",
    "paired_power = 0\n",
    "bad_cnt = 0\n",
    "power = 0\n",
    "N = 30000\n",
    "for i in tqdm_notebook(range(N)):\n",
    "    before = sps.expon(scale=1000).rvs(4000)\n",
    "    C_b, T_b = before[:2000], before[2000:]\n",
    "\n",
    "    C = C_b + sps.norm(loc=0, scale=100).rvs(2000)\n",
    "    T = T_b + sps.norm(loc=0, scale=100).rvs(2000)\n",
    "    T *= 1.01\n",
    "\n",
    "    _, _, _, left_bound, right_bound = relative_ttest(C, T)\n",
    "    _, _, _, left_bound_paired, right_bound_paired = paired_relative_ttest(C, T)\n",
    "\n",
    "    if left_bound > 0.01 or right_bound < 0.01:\n",
    "        bad_cnt += 1\n",
    "        \n",
    "    if left_bound_paired > 0.01 or right_bound_paired < 0.01:\n",
    "        bad_cnt_paired += 1\n",
    "\n",
    "    if left_bound > 0:\n",
    "        power += 1\n",
    "\n",
    "    if left_bound_paired > 0:\n",
    "        paired_power += 1\n",
    "    \n",
    "left_real_level, right_real_level = proportion_confint(count = bad_cnt, nobs = N, \n",
    "                                                         alpha=0.05, method='wilson')\n",
    "left_real_level_paired, right_real_level_paired = proportion_confint(count = bad_cnt_paired, nobs = N, \n",
    "                                                         alpha=0.05, method='wilson')\n",
    "print(f\"Реальный уровень значимости для обычного критрия: {round(bad_cnt / N, 4)};\"\\\n",
    "      f\" доверительный интервал: [{round(left_real_level, 5)}, {round(right_real_level, 5)}]\")\n",
    "print(f\"Реальный уровень значимости для спаренного критерия: {round(bad_cnt_paired / N, 4)};\"\\\n",
    "      f\" доверительный интервал: [{round(left_real_level_paired, 5)}, {round(right_real_level_paired, 5)}]\")\n",
    "print(f\"Мощность спаренного критрерия VS мощность обычного критерия: {paired_power / N} VS. {power /N}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9304636b",
   "metadata": {},
   "source": [
    "## Пример, когда от сплитования может стать хуже"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "0f50e641",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2d75a5bbdc834cc6a828adc04d77602d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "925ab55786c84ba4800d2271da1ff4f3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ширина доверительного интервала в случае неспаренного критерия VS спаренного критерия: 123.8 VS. 174.96\n",
      "pvalue сравнения: 0.0\n"
     ]
    }
   ],
   "source": [
    "ci_NOT_splitted_length = []\n",
    "ci_splitted_length = []\n",
    "N = 1000\n",
    "# сплитование\n",
    "for i in tqdm_notebook(range(N)):\n",
    "    before = sps.expon(scale=1000).rvs(4000)\n",
    "    # уже отсортированы\n",
    "    C_b, T_b = splitter(before)\n",
    "\n",
    "    C = C_b\n",
    "    T = -1 * T_b\n",
    "\n",
    "    _, _, ci_splitted_sample, _, _ = paired_ttest(C, T)\n",
    "    ci_splitted_length.append(ci_splitted_sample)\n",
    "###################################\n",
    "    \n",
    "# не было сплитования\n",
    "for i in tqdm_notebook(range(N)):\n",
    "    before = sps.expon(scale=1000).rvs(4000)\n",
    "    C_b, T_b = before[:2000], before[2000:]\n",
    "\n",
    "    C = C_b\n",
    "    T = -1 * T_b\n",
    "\n",
    "    _, _, ci_NOT_splitted_sample, _, _ = absolute_ttest(C, T)\n",
    "\n",
    "    ci_NOT_splitted_length.append(ci_NOT_splitted_sample)\n",
    "print(\"ширина доверительного интервала в случае неспаренного критерия VS спаренного критерия:\"\\\n",
    "      f\" {round(np.mean(ci_NOT_splitted_length), 2)} VS. {round(np.mean(ci_splitted_length), 2)}\")\n",
    "print(f\"pvalue сравнения: {absolute_ttest(np.array(ci_NOT_splitted_length), np.array(ci_splitted_length)).pvalue}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c470986f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}