Приложение-Множественные_сравнения.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0e967bde",
   "metadata": {},
   "source": [
    "### Приложение: множественные сравнения"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32f36f66",
   "metadata": {},
   "source": [
    "Байесовский подход "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9221b915",
   "metadata": {},
   "source": [
    "$$\n",
    "P(data | p) = Binom(p, s, N)\n",
    "\\\\\n",
    "P(p) = Uniform(0, 1)\n",
    "$$\n",
    "\n",
    "$$\n",
    "P(p | data) = \\frac{ P(data | p) P(p) }{P(data)}\n",
    "= Beta(p; \\alpha, \\beta)\n",
    "\\\\\n",
    "\\alpha = s + 1, \\qquad \\beta = N - s + 1\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d35e4e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "np.random.seed(7)\n",
    "\n",
    "from collections import namedtuple\n",
    "\n",
    "import scipy.stats as stats\n",
    "import plotly.express as px\n",
    "import plotly.graph_objects as go\n",
    "from plotly.subplots import make_subplots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7042388",
   "metadata": {},
   "outputs": [],
   "source": [
    "def posterior_for_binom_and_uniform_prior(n_success, n_trials):\n",
    "    alpha_prior = 1\n",
    "    beta_prior = 1\n",
    "    alpha_post = alpha_prior + n_success\n",
    "    beta_post = beta_prior + (n_trials - n_success)\n",
    "    return stats.beta(alpha_post, beta_post)\n",
    "\n",
    "pa_exact = 0.100\n",
    "pb_exact = 0.105\n",
    "pc_exact = 0.110\n",
    "\n",
    "N = 1000\n",
    "\n",
    "sa = stats.binom.rvs(n=N, p=pa_exact)\n",
    "sb = stats.binom.rvs(n=N, p=pb_exact)\n",
    "sc = stats.binom.rvs(n=N, p=pc_exact)\n",
    "print(sa, sb, sc)\n",
    "\n",
    "a_post_dist = posterior_for_binom_and_uniform_prior(sa, N)\n",
    "b_post_dist = posterior_for_binom_and_uniform_prior(sb, N)\n",
    "c_post_dist = posterior_for_binom_and_uniform_prior(sc, N)\n",
    "\n",
    "x = np.linspace(0, 1, 1001)\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Scatter(x=x, y=a_post_dist.pdf(x), mode='lines', name=\"A\"))\n",
    "fig.add_trace(go.Scatter(x=x, y=b_post_dist.pdf(x), mode='lines', name=\"B\"))\n",
    "fig.add_trace(go.Scatter(x=x, y=c_post_dist.pdf(x), mode='lines', name=\"C\"))\n",
    "fig.update_layout(title='Posterior Distributions',\n",
    "                  xaxis_title='p',\n",
    "                  yaxis_title='Prob Density',\n",
    "                  xaxis_range=[0, 0.3],\n",
    "                  hovermode=\"x\",\n",
    "                  height=450)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c417bd56",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_sample = 50000\n",
    "a_post_sample = a_post_dist.rvs(size=n_sample)\n",
    "b_post_sample = b_post_dist.rvs(size=n_sample)\n",
    "c_post_sample = c_post_dist.rvs(size=n_sample)\n",
    "\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Histogram(x=a_post_sample, histnorm='probability density', \n",
    "                           name='A',\n",
    "                           opacity=0.6))\n",
    "fig.add_trace(go.Histogram(x=b_post_sample, histnorm='probability density', \n",
    "                           name='B',\n",
    "                           opacity=0.6))\n",
    "fig.add_trace(go.Histogram(x=c_post_sample, histnorm='probability density', \n",
    "                           name='C',\n",
    "                           opacity=0.6))\n",
    "fig.update_layout(barmode='overlay')\n",
    "fig.update_layout(title='Posterior Samples',\n",
    "                  xaxis_title='n',\n",
    "                  yaxis_title='Prob Density',\n",
    "                  hovermode=\"x\",\n",
    "                  height=450)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b11be11",
   "metadata": {},
   "outputs": [],
   "source": [
    "diff_a_b = a_post_sample - b_post_sample\n",
    "diff_a_c = a_post_sample - c_post_sample\n",
    "\n",
    "diff_b_a = b_post_sample - a_post_sample\n",
    "diff_b_c = b_post_sample - c_post_sample\n",
    "\n",
    "diff_c_a = c_post_sample - a_post_sample\n",
    "diff_c_b = c_post_sample - b_post_sample\n",
    "\n",
    "prob_a_gt_b_and_c = len(a_post_sample[(diff_a_b > 0) & (diff_a_c > 0)]) / len(a_post_sample)\n",
    "prob_b_gt_a_and_c = len(b_post_sample[(diff_b_a > 0) & (diff_b_c > 0)]) / len(b_post_sample)\n",
    "prob_c_gt_b_and_a = len(c_post_sample[(diff_c_a > 0) & (diff_c_b > 0)]) / len(c_post_sample)\n",
    "\n",
    "print(f'Pr(Group A is the best) = Pr(p_a > p_b & p_a > p_c): {prob_a_gt_b_and_c}')\n",
    "print(f'Pr(Group B is the best) = Pr(p_b > p_a & p_b > p_c): {prob_b_gt_a_and_c}')\n",
    "print(f'Pr(Group C is the best) = Pr(p_c > p_b & p_c > p_a): {prob_c_gt_b_and_a}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4830940",
   "metadata": {},
   "source": [
    "Видно, что вероятность того, что C лучше A и B одновременно ниже, чем то, что C лучшая при попарном сравнении.\n",
    "\n",
    "$$\n",
    "P(p_c > p_a, p_b) < \\min(P(p_c > p_a), P(p_c > p_b))\n",
    "$$\n",
    "\n",
    "Это и есть проблема множественного сравнения."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2bbc6a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "diff_c_a = c_post_sample - a_post_sample\n",
    "prob_c_gt_a = len(diff_c_a[diff_c_a > 0]) / len(c_post_sample)\n",
    "print(f'Pr(p_c > p_a): {prob_c_gt_a}')\n",
    "\n",
    "diff_c_b = c_post_sample - b_post_sample\n",
    "prob_c_gt_b = len(diff_c_b[diff_c_b > 0]) / len(c_post_sample)\n",
    "print(f'Pr(p_c > p_b): {prob_c_gt_b}')\n",
    "\n",
    "prob_c_gt_b_and_a = len(diff_c_b[(diff_c_b > 0) & (diff_c_a > 0)]) / len(c_post_sample)\n",
    "print(f'Pr(p_c > p_b & p_c > p_a): {prob_c_gt_b_and_a}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6b84e8c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50591f8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def p_best_binomial(df_gr_n_s):\n",
    "    t = pd.DataFrame()\n",
    "    n_sample = 50000 * len(df_gr_n_s.index.unique())\n",
    "    for g in df_gr_n_s.index.unique():\n",
    "        d_post = posterior_for_binom_and_uniform_prior(df_gr_n_s['S'][g], df_gr_n_s['N'][g])\n",
    "        post_samp = d_post.rvs(size=n_sample)\n",
    "        t[g] = post_samp\n",
    "    return t.idxmax(axis=1).value_counts(normalize=True)\n",
    "\n",
    "\n",
    "df = pd.DataFrame({'gr': ['A', 'B', 'C'], 'N':N, 'S':[sa, sb, sc]})\n",
    "df['p_best'] = p_best_binomial(df.set_index('gr')[['N', 'S']])\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8be83c11",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0bae355a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "import plotly.graph_objects as go\n",
    "\n",
    "prior_conv = 0.1\n",
    "\n",
    "a_total = 1000\n",
    "a_conv = a_total * prior_conv\n",
    "\n",
    "b_total = 1000\n",
    "b_conv = b_total * (1.1 * prior_conv)\n",
    "\n",
    "c_total = 1000\n",
    "c_conv = c_total * (1.2 * prior_conv)\n",
    "\n",
    "x = np.linspace(0, 1, 1001)\n",
    "fig = go.Figure()\n",
    "\n",
    "\n",
    "# Uniform prior\n",
    "beta_prior = 1\n",
    "alpha_prior = 1\n",
    "\n",
    "fig.add_trace(go.Scatter(x=x, y=stats.beta.pdf(x, alpha_prior, beta_prior), mode='lines',\n",
    "                         name=f\"alpha_prior={alpha_prior}, beta_prior={beta_prior}\"))\n",
    "\n",
    "alpha_post = alpha_prior + a_conv\n",
    "beta_post = beta_prior + (a_total - a_conv)\n",
    "fig.add_trace(go.Scatter(x=x, y=stats.beta.pdf(x, alpha_post, beta_post), mode='lines',\n",
    "                         name=f\"<br>A: total = {a_total}, conv = {a_conv}<br>alpha_prior={alpha_prior}, beta_prior={beta_prior},<br>alpha_post={alpha_post}, beta_post={beta_post}\"))\n",
    "\n",
    "alpha_post = alpha_prior + b_conv\n",
    "beta_post = beta_prior + (b_total - b_conv)\n",
    "fig.add_trace(go.Scatter(x=x, y=stats.beta.pdf(x, alpha_post, beta_post), mode='lines',\n",
    "                         name=f\"<br>B: total = {b_total}, conv = {b_conv}<br>alpha_prior={alpha_prior}, beta_prior={beta_prior},<br>alpha_post={alpha_post}, beta_post={beta_post}\"))\n",
    "\n",
    "alpha_post = alpha_prior + c_conv\n",
    "beta_post = beta_prior + (c_total - c_conv)\n",
    "fig.add_trace(go.Scatter(x=x, y=stats.beta.pdf(x, alpha_post, beta_post), mode='lines',\n",
    "                         name=f\"<br>C: total = {c_total}, conv = {c_conv}<br>alpha_prior={alpha_prior}, beta_prior={beta_prior},<br>alpha_post={alpha_post}, beta_post={beta_post}\"))\n",
    "\n",
    "\n",
    "fig.update_layout(title='Posterior and Prior Distributions',\n",
    "                  xaxis_title='p',\n",
    "                  yaxis_title='Prob Density',\n",
    "                  hovermode=\"x\",\n",
    "                  height=550)\n",
    "\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c8988a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_sample = 50000\n",
    "alpha_prior = 1\n",
    "beta_prior = 1\n",
    "a = alpha_prior + a_conv\n",
    "b = beta_prior + (a_total - a_conv)\n",
    "post_sample_a = np.random.beta(a, b, n_sample)\n",
    "a = alpha_prior + b_conv\n",
    "b = beta_prior + (b_total - b_conv)\n",
    "post_sample_b = np.random.beta(a, b, n_sample)\n",
    "a = alpha_prior + c_conv\n",
    "b = beta_prior + (c_total - c_conv)\n",
    "post_sample_c = np.random.beta(a, b, n_sample)\n",
    "\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Histogram(x=post_sample_a, histnorm='probability density', \n",
    "                           name='A',\n",
    "                           opacity=0.6))\n",
    "fig.add_trace(go.Histogram(x=post_sample_b, histnorm='probability density', \n",
    "                           name='B',\n",
    "                           opacity=0.6))\n",
    "fig.add_trace(go.Histogram(x=post_sample_c, histnorm='probability density', \n",
    "                           name='C',\n",
    "                           opacity=0.6))\n",
    "fig.update_layout(barmode='overlay')\n",
    "fig.update_layout(title='Posterior Sample',\n",
    "                  xaxis_title='n',\n",
    "                  yaxis_title='Prob Density',\n",
    "                  hovermode=\"x\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96ada8ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "diff_c_a = post_sample_c - post_sample_a\n",
    "prob_c_ge_a = len(diff_c_a[diff_c_a > 0]) / len(post_sample_c)\n",
    "print(f'Pr(p_c > p_a): {prob_c_ge_a}')\n",
    "\n",
    "diff_c_b = post_sample_c - post_sample_b\n",
    "prob_c_ge_b = len(diff_c_b[diff_c_b > 0]) / len(post_sample_c)\n",
    "print(f'Pr(p_c > p_b): {prob_c_ge_b}')\n",
    "\n",
    "prob_c_ge_b_and_a = len(diff_c_b[(diff_c_b > 0) & (diff_c_a > 0)]) / len(post_sample_c)\n",
    "print(f'Pr(p_c > p_b & p_c > p_a): {prob_c_ge_b_and_a}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0affb96e",
   "metadata": {},
   "source": [
    "Видно, что вероятность того, что C лучше A и B одновременно ниже, чем то, что C лучшая при попарном сравнении.\n",
    "$$\n",
    "P(p_c > p_a, p_b) < \\min(P(p_c > p_a), P(p_c > p_b))\n",
    "$$\n",
    "Это и есть проблема множественного сравнения."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3cdf6461",
   "metadata": {},
   "source": [
    "С вероятностью $P(p_c > p_a, p_b)$ можно делать то же, что и с $P(p_c > p_a), P(p_c > p_b)$.  \n",
    "Моделировать, сколько нужно данных для достижения определенного уровня уверенности или записать функцию потерь и моделировать ее минимизацию. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1523e6f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d1c2435",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "e3d06a82",
   "metadata": {},
   "source": [
    "### Средние"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfdaf490",
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import namedtuple\n",
    "\n",
    "ConjugateNormalParams = namedtuple('ConjugateNormalParams', 'mu sigma k a b')\n",
    "\n",
    "def initial_parameters(mu, sigma, k=1/25, a=2, b=1):\n",
    "    return ConjugateNormalParams(mu=mu, sigma=sigma, k=k, a=a, b=b)\n",
    "\n",
    "def update_conj_parameters(x, conj_norm_pars):\n",
    "    mu_p = (x + conj_norm_pars.k * conj_norm_pars.mu) / (conj_norm_pars.k + 1)\n",
    "    sigma_p = conj_norm_pars.sigma\n",
    "    k_p = conj_norm_pars.k + 1\n",
    "    a_p = conj_norm_pars.a + 1/2\n",
    "    b_p = conj_norm_pars.b + conj_norm_pars.k / (conj_norm_pars.k + 1) * (x - conj_norm_pars.mu)**2 / (2 * conj_norm_pars.sigma**2)\n",
    "    return ConjugateNormalParams(mu=mu_p, sigma=sigma_p, k=k_p, a=a_p, b=b_p)\n",
    "\n",
    "def compute_posterior_parameters(sample, n_split):\n",
    "    means = reshape_and_compute_means(sample, n_split)\n",
    "    mean_1 = means[0]\n",
    "    sigma_1 = sample[0:n_split].std() / np.sqrt(n_split)\n",
    "    pars = []\n",
    "    pars.append(initial_parameters(mu=mean_1, sigma=sigma_1)) \n",
    "    for x in means[1:]:\n",
    "        new_pars = update_conj_parameters(x, pars[-1])\n",
    "        pars.append(new_pars)\n",
    "    return pars\n",
    "\n",
    "def reshape_and_compute_means(sample, n_split):\n",
    "    n_means = len(sample) // n_split\n",
    "    samp_reshaped = np.reshape(sample[0 : n_means * n_split], (n_means, n_split))\n",
    "    means = np.array([x.mean() for x in samp_reshaped])\n",
    "    return means\n",
    "\n",
    "def mu_marginal_distrib(conj_norm_pars):\n",
    "    df = 2 * conj_norm_pars.a\n",
    "    loc = conj_norm_pars.mu\n",
    "    scale = np.sqrt(conj_norm_pars.sigma**2 / conj_norm_pars.k * conj_norm_pars.b / conj_norm_pars.a)\n",
    "    return stats.t(df=df, loc=loc, scale=scale)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "170e1cd9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def p_best_means(df_gr_samp):\n",
    "    t = pd.DataFrame()\n",
    "    n_split = 25\n",
    "    n_sample = 30000 * len(df_gr_samp.index.unique())\n",
    "    for g in df_gr_samp.index.unique():\n",
    "        pars = compute_posterior_parameters(np.array(df_gr_samp['samp'][g]), n_split)\n",
    "        last_pars = pars[-1]\n",
    "        mu_dist = mu_marginal_distrib(last_pars)        \n",
    "        mu_samp = mu_dist.rvs(size=n_sample)\n",
    "        t[g] = mu_samp\n",
    "    return t.idxmax(axis=1).value_counts(normalize=True)\n",
    "\n",
    "\n",
    "n = 1000\n",
    "c_a_exact = 1.5\n",
    "exact_dist_a = stats.lomax(c=c_a_exact)\n",
    "samp_a = exact_dist_a.rvs(size=n)\n",
    "#c_b_exact = 2.8 \n",
    "c_b_exact = c_a_exact * 0.95\n",
    "exact_dist_b = stats.lomax(c=c_b_exact)\n",
    "samp_b = exact_dist_b.rvs(size=n)\n",
    "\n",
    "\n",
    "df = pd.concat([\n",
    "     pd.DataFrame({'gr':['A']*len(samp_a), 'samp':samp_a}),\n",
    "     pd.DataFrame({'gr':['B']*len(samp_b), 'samp':samp_b})])\n",
    "df_res = df.groupby('gr')['samp'].mean().rename('means').to_frame()\n",
    "df_res['p_best'] = p_best_means(df.set_index('gr'))\n",
    "df_res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab796c71",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4222b51d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "3892e7fe",
   "metadata": {},
   "source": [
    "# Ссылки"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eed3a57e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "350be6fd",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "9052afb6",
   "metadata": {},
   "source": [
    "В проверке статистических гипотез комбинирование результатов нескольких статистических тестов обычно делается разного рода поправками [[MultipleComp](https://en.wikipedia.org/wiki/Multiple_comparisons_problem), [FWER](https://en.wikipedia.org/wiki/Family-wise_error_rate), [Bonf](https://en.wikipedia.org/wiki/Bonferroni_correction)]. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "755bb480",
   "metadata": {},
   "source": [
    "[[MultipleComp](https://en.wikipedia.org/wiki/Multiple_comparisons_problem)] -   \n",
    "[[FWER](https://en.wikipedia.org/wiki/Family-wise_error_rate)] -  \n",
    "[[Bonf](https://en.wikipedia.org/wiki/Bonferroni_correction)] -   "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}