{ "cells": [ { "cell_type": "markdown", "id": "77d7a001", "metadata": {}, "source": [ "# Chapter 15: Comparing means" ] }, { "cell_type": "code", "execution_count": 104, "id": "2673e835-2b2a-412d-bcef-6b862fa2f630", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The rpy2.ipython extension is already loaded. To reload it, use:\n", " %reload_ext rpy2.ipython\n" ] } ], "source": [ "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "from scipy.stats import norm, t, binom, scoreatpercentile, binomtest, f\n", "import pingouin as pg\n", "import matplotlib\n", "import statsmodels.api as sm\n", "import statsmodels.formula.api as smf\n", "from sklearn.metrics import r2_score, mean_squared_error\n", "from sklearn.model_selection import KFold\n", "\n", "import rpy2.robjects as ro\n", "from rpy2.robjects.packages import importr\n", "from rpy2.robjects import pandas2ri\n", "pandas2ri.activate()\n", "from rpy2.robjects.conversion import localconverter\n", "# import NHANES package\n", "base = importr('NHANES')\n", "%load_ext rpy2.ipython\n", "\n", "with localconverter(ro.default_converter + pandas2ri.converter):\n", " NHANES = ro.conversion.rpy2py(ro.r['NHANES'])\n", "\n", " \n", "NHANES = NHANES.drop_duplicates(subset='ID').replace(-2147483648, np.nan)\n", "NHANES_adult = NHANES.dropna(subset=['BPDiaAve']).query('Age > 17 and BPDiaAve > 0')\n", "\n", "rng = np.random.default_rng(1234)" ] }, { "cell_type": "markdown", "id": "c44a6659", "metadata": {}, "source": [ "## Binomial test for a single proportion" ] }, { "cell_type": "code", "execution_count": 48, "id": "c5cc9cb6-7b22-4a86-ac88-0565536e5649", "metadata": { "Rmd_chunk_options": "echo=FALSE", "jupyter": { "output_hidden": false }, "kernel": "R", "tags": [ "report_output" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "BinomTestResult(k=22, n=200, alternative='greater', proportion_estimate=0.11, pvalue=1.0)\n" ] } ], "source": [ "NHANES_adult['Hypertensive'] = NHANES_adult.BPDiaAve > 80\n", "\n", "NHANES_sample = NHANES_adult.sample(200, random_state=rng)\n", "\n", "# compute sign test for differences between first and second measurement\n", "npos = NHANES_sample.Hypertensive.sum()\n", "\n", "print(binomtest(npos, NHANES_sample.shape[0], alternative='greater'))" ] }, { "cell_type": "markdown", "id": "1fbad11c", "metadata": {}, "source": [ "## T-test output" ] }, { "cell_type": "code", "execution_count": 49, "id": "83608a00-d63d-46f7-9ab9-b7bae815c3a9", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Tdofalternativep-valCI95%cohen-dBF10power
T-test-14.293279199greater1.0[67.29, inf]1.0106873.983e-300.0
\n", "
" ], "text/plain": [ " T dof alternative p-val CI95% cohen-d BF10 \\\n", "T-test -14.293279 199 greater 1.0 [67.29, inf] 1.010687 3.983e-30 \n", "\n", " power \n", "T-test 0.0 " ] }, "execution_count": 49, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pg.ttest(NHANES_sample.BPDiaAve, 80, alternative='greater')" ] }, { "cell_type": "markdown", "id": "024dbd24-67e1-4ad6-8aea-240cfe9d07cb", "metadata": { "kernel": "Markdown", "tags": [ "report_output" ] }, "source": [ "## Bayes factor output" ] }, { "cell_type": "code", "execution_count": 50, "id": "6484b24c", "metadata": { "Rmd_chunk_options": "echo=FALSE, message=FALSE, warning=FALSE", "jupyter": { "output_hidden": false }, "kernel": "R", "tags": [ "report_output" ] }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "R[write to console]: t is large; approximation invoked.\n", "\n", "R[write to console]: t is large; approximation invoked.\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Bayes factor analysis\n", "--------------\n", "[1] Alt., r=0.707 -Inf" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "NHANES_sample = NHANES_adult.dropna(subset=['TVHrsDay', 'RegularMarij']).sample(200, random_state=rng)\n", "\n", "tv_recoder = {\n", " \"More_4_hr\": 5,\n", " \"4_hr\": 4,\n", " \"2_hr\": 2,\n", " \"1_hr\": 1,\n", " \"3_hr\": 3,\n", " \"0_to_1_hr\": 0.5,\n", " \"0_hrs\": 0\n", "}\n", "\n", "NHANES_sample['TVHrsNum'] = [tv_recoder[i] for i in NHANES_sample.TVHrsDay]\n", "NHANES_sample['RegularMarijNum'] = [1 if i == \"Yes\" else 0 for i in NHANES_sample['RegularMarij']]\n", "\n", "fig, ax = plt.subplots(1, 2, figsize=(12,6))\n", "\n", "sns.violinplot(data=NHANES_sample, x='RegularMarij', y='TVHrsNum', ax=ax[0])\n", "\n", "\n", "sns.violinplot(data=NHANES_sample, x='RegularMarij', y='TVHrsNum', ax=ax[1])\n", "sns.regplot(data=NHANES_sample, x='RegularMarijNum', y='TVHrsNum', ax=ax[1], \n", " scatter=False, ci=None, color='black')\n", "\n", "for x in range(2):\n", " ax[x].set_xlabel(\"Regular marijuana user\")\n", " ax[x].set_ylabel(\"TV hours per day\")\n" ] }, { "cell_type": "markdown", "id": "fab81cbb-b65c-4b5b-8a95-7596aa46e87b", "metadata": {}, "source": [ "## T-test result" ] }, { "cell_type": "code", "execution_count": 52, "id": "c8eb8c48-b621-48d1-9dcd-5841d477de43", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Tdofalternativep-valCI95%cohen-dBF10power
T-test1.470448198greater0.071514[-0.04, inf]0.2434570.9560.4288
\n", "
" ], "text/plain": [ " T dof alternative p-val CI95% cohen-d BF10 \\\n", "T-test 1.470448 198 greater 0.071514 [-0.04, inf] 0.243457 0.956 \n", "\n", " power \n", "T-test 0.4288 " ] }, "execution_count": 52, "metadata": {}, "output_type": "execute_result" } ], "source": [ "tt = pg.ttest(x=NHANES_sample.query('RegularMarij == \"Yes\"').TVHrsNum, \n", " y=NHANES_sample.query('RegularMarij == \"No\"').TVHrsNum,\n", " alternative='greater', correction=False)\n", "tt" ] }, { "cell_type": "markdown", "id": "f63207df", "metadata": {}, "source": [ "## Linear model summary" ] }, { "cell_type": "code", "execution_count": 53, "id": "66ffd6f5-14b8-4c24-9be6-3671c263e761", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: TVHrsNum R-squared: 0.011\n", "Model: OLS Adj. R-squared: 0.006\n", "Method: Least Squares F-statistic: 2.162\n", "Date: Thu, 01 Sep 2022 Prob (F-statistic): 0.143\n", "Time: 11:59:45 Log-Likelihood: -353.68\n", "No. Observations: 200 AIC: 711.4\n", "Df Residuals: 198 BIC: 718.0\n", "Df Model: 1 \n", "Covariance Type: nonrobust \n", "=======================================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "---------------------------------------------------------------------------------------\n", "Intercept 2.2467 0.116 19.432 0.000 2.019 2.475\n", "RegularMarij[T.Yes] 0.3470 0.236 1.470 0.143 -0.118 0.812\n", "==============================================================================\n", "Omnibus: 15.476 Durbin-Watson: 1.966\n", "Prob(Omnibus): 0.000 Jarque-Bera (JB): 10.929\n", "Skew: 0.452 Prob(JB): 0.00423\n", "Kurtosis: 2.297 Cond. No. 2.50\n", "==============================================================================\n", "\n", "Notes:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n" ] } ], "source": [ "# print summary of linear regression to perform t-test\n", "\n", "lmresult = smf.ols(formula='TVHrsNum ~ RegularMarij', data=NHANES_sample).fit()\n", "\n", "print(lmresult.summary())" ] }, { "cell_type": "markdown", "id": "d4d78775", "metadata": {}, "source": [ "## Bayes factor for mean differences" ] }, { "cell_type": "code", "execution_count": 54, "id": "7dd94a73", "metadata": { "Rmd_chunk_options": "echo=FALSE", "jupyter": { "output_hidden": false }, "kernel": "R", "tags": [ "report_output" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Bayes factor analysis\n", "--------------\n", "[1] Alt., r=0.707 0" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "NHANES_sample = NHANES_adult.dropna(\n", " subset=['BPSys1', 'BPSys2']).sample(200, random_state=rng)[['ID', 'BPSys1', 'BPSys2']]\n", "NHANES_sample_long = pd.melt(NHANES_sample, id_vars='ID', var_name='timepoint', value_name='BPSys')\n", "NHANES_sample_long.timepoint = NHANES_sample_long.timepoint.replace({'BPSys1': 'Time1', 'BPSys2': 'Time2'})\n", "fig, ax = plt.subplots(1, 2, figsize=(12,6))\n", "\n", "sns.violinplot(data=NHANES_sample_long, x='timepoint', y='BPSys', ax=ax[0])\n", "\n", "sns.violinplot(data=NHANES_sample_long, x='timepoint', y='BPSys', ax=ax[1])\n", "sns.lineplot(data=NHANES_sample_long, x='timepoint', y='BPSys', hue='ID',\n", " ax=ax[1], color='black', legend=False)\n", "ax[1].set_xlim((-1, 2))\n", "ax[1].set_ylim(ax[0].get_ylim())\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "id": "054242b1", "metadata": {}, "source": [ "## T-test output" ] }, { "cell_type": "code", "execution_count": 74, "id": "c732702f-0552-432a-89d1-df1ae398d0eb", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Tdofalternativep-valCI95%cohen-dBF10power
T-test0.466004398greater0.320734[-2.13, inf]0.04660.2460.119072
\n", "
" ], "text/plain": [ " T dof alternative p-val CI95% cohen-d BF10 \\\n", "T-test 0.466004 398 greater 0.320734 [-2.13, inf] 0.0466 0.246 \n", "\n", " power \n", "T-test 0.119072 " ] }, "execution_count": 74, "metadata": {}, "output_type": "execute_result" } ], "source": [ "tt = pg.ttest(x=NHANES_sample_long.query('timepoint == \"Time1\"').BPSys, \n", " y=NHANES_sample_long.query('timepoint == \"Time2\"').BPSys,\n", " alternative='greater', correction=False)\n", "tt" ] }, { "cell_type": "markdown", "id": "0854d27a", "metadata": {}, "source": [ "## Figure 15.3" ] }, { "cell_type": "code", "execution_count": 79, "id": "d7a4c431-6151-4bdf-8de9-702eac16751b", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 79, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAGdCAYAAAAIbpn/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAlsklEQVR4nO3df1BV54H/8Q9GuGrgXsUfXKigGFONMZhZasjND0eFSEjqasN00h/bmqxrqoumSrdN2CYY2M1g7a6adAmmW4vpJJasOzHZZBrdiBFnR6BKtIRkw0THRKJe6CaFqyRcKJzvH473m1sFAa+c8+D7NXNmPOc59/DJM0Y/nvNwiLIsyxIAAICBRtgdAAAAYLAoMgAAwFgUGQAAYCyKDAAAMBZFBgAAGIsiAwAAjEWRAQAAxqLIAAAAY420O8DV1tPTo9OnTysuLk5RUVF2xwEAAP1gWZbOnj2rpKQkjRjR+32XYV9kTp8+reTkZLtjAACAQWhqatLkyZN7HR/2RSYuLk7S+Ylwu902pwEAAP0RCASUnJwc+nu8N8O+yFx4nOR2uykyAAAY5nLLQljsCwAAjEWRAQAAxqLIAAAAY1FkAACAsSgyAADAWBQZAABgLIoMAAAwFkUGAAAYiyIDAACMRZEBAADGckyR2bBhg6KiorR27drQsY6ODuXl5Wn8+PGKjY1Vbm6umpub7QsJAAAcxRFF5tChQ3r++eeVlpYWdnzdunV6/fXXtXPnTlVVVen06dN64IEHbEoJAACcxvYic+7cOX33u9/Vv//7v2vcuHGh421tbdq2bZs2bdqkhQsXKj09XeXl5Tp48KBqampsTAzgarIsqb39/GZZdqcB4HS2F5m8vDzdf//9ysrKCjteV1enrq6usOMzZ85USkqKqqure71eMBhUIBAI2wCY4/PPpdjY89vnn9udBoDTjbTzi1dUVOidd97RoUOHLhrz+/2KiYnR2LFjw44nJCTI7/f3es2SkhIVFRVFOioAAHAg2+7INDU16Yc//KFeeukljRo1KmLXLSgoUFtbW2hramqK2LUBAICz2FZk6urq1NLSor/6q7/SyJEjNXLkSFVVVenZZ5/VyJEjlZCQoM7OTrW2toZ9rrm5WV6vt9frulwuud3usA0AAAxPtj1ayszM1Lvvvht27OGHH9bMmTP12GOPKTk5WdHR0aqsrFRubq4kqbGxUSdPnpTP57MjMgAAcBjbikxcXJxmz54dduz666/X+PHjQ8eXL1+u/Px8xcfHy+12a82aNfL5fLr99tvtiAwAABzG1sW+l7N582aNGDFCubm5CgaDys7O1nPPPWd3LAAA4BBRljW839QQCATk8XjU1tbGehnAAO3t57/1WpLOnZOuv97ePADs0d+/v21/jwwAAMBgUWQAAICxKDIAAMBYFBkAAGAsigwAADAWRQYAABiLIgMAAIxFkQEAAMaiyAAAAGNRZAAAgLEoMgAAwFgUGQAAYCyKDAAAMBZFBgAAGIsiAwAAjEWRAQAAxqLIAAAAY1FkAACAsSgyAADAWBQZAABgLIoMAAAwFkUGAAAYiyIDAACMRZEBAADGosgAAABjUWQAAICxKDIAAMBYFBkAAGAsigwAADAWRQYAABiLIgMAAIxFkQEAAMaiyAAAAGPZWmTKysqUlpYmt9stt9stn8+nN998MzQ+f/58RUVFhW0rV660MTEAAHCSkXZ+8cmTJ2vDhg268cYbZVmWXnjhBS1ZskRHjhzRzTffLElasWKFiouLQ58ZM2aMXXEBAIDD2FpkFi9eHLb/9NNPq6ysTDU1NaEiM2bMGHm9XjviAQAAh3PMGpnu7m5VVFSovb1dPp8vdPyll17ShAkTNHv2bBUUFOjzzz/v8zrBYFCBQCBsAwAAw5Otd2Qk6d1335XP51NHR4diY2O1a9cuzZo1S5L0ne98R1OmTFFSUpLq6+v12GOPqbGxUa+88kqv1yspKVFRUdFQxQcwABkZGfL7/b2Oe71e7dtXO4SJAJguyrIsy84AnZ2dOnnypNra2vSf//mf+tWvfqWqqqpQmfmyffv2KTMzU8eOHdMNN9xwyesFg0EFg8HQfiAQUHJystra2uR2u6/afweAy5syZUrYmre/VFhYqPff/1ixsef3z52Trr9+iMIBcJRAICCPx3PZv79tvyMTExOj6dOnS5LS09N16NAhPfPMM3r++ecvOjcjI0OS+iwyLpdLLpfr6gUGAACO4Zg1Mhf09PSE3VH5sqNHj0qSEhMThzARAABwKlvvyBQUFCgnJ0cpKSk6e/asduzYof3792vPnj06fvy4duzYofvuu0/jx49XfX291q1bp3nz5iktLc3O2AAAwCFsLTItLS36/ve/rzNnzsjj8SgtLU179uzRPffco6amJu3du1dbtmxRe3u7kpOTlZubqyeeeMLOyAAAwEFsLTLbtm3rdSw5OVlVVVVDmAYAAJjGcWtkAAAA+osiAwAAjEWRAQAAxqLIAAAAY1FkAACAsSgyAADAWBQZAABgLIoMAAAwFkUGAAAYiyIDAACMRZEBAADGosgAAABjUWQAAICxKDIAAMBYFBkAAGAsigwAADAWRQYAABiLIgMAAIxFkQEAAMaiyAAAAGNRZAAAgLEoMgAAwFgUGQAAYCyKDAAAMBZFBgAAGIsiAwAAjEWRAQAAxqLIAAAAY1FkAACAsSgyAADAWBQZAABgLIoMAAAwFkUGAAAYy9YiU1ZWprS0NLndbrndbvl8Pr355puh8Y6ODuXl5Wn8+PGKjY1Vbm6umpubbUwMAACcxNYiM3nyZG3YsEF1dXU6fPiwFi5cqCVLlui9996TJK1bt06vv/66du7cqaqqKp0+fVoPPPCAnZEBAICDjLTziy9evDhs/+mnn1ZZWZlqamo0efJkbdu2TTt27NDChQslSeXl5brppptUU1Oj22+/3Y7IAADAQRyzRqa7u1sVFRVqb2+Xz+dTXV2durq6lJWVFTpn5syZSklJUXV1da/XCQaDCgQCYRsAABiebC8y7777rmJjY+VyubRy5Urt2rVLs2bNkt/vV0xMjMaOHRt2fkJCgvx+f6/XKykpkcfjCW3JyclX+b8AAADYxfYiM2PGDB09elS1tbVatWqVli1bpvfff3/Q1ysoKFBbW1toa2pqimBaAADgJLaukZGkmJgYTZ8+XZKUnp6uQ4cO6ZlnntGDDz6ozs5Otba2ht2VaW5ultfr7fV6LpdLLpfrascGAAAOYPsdmb/U09OjYDCo9PR0RUdHq7KyMjTW2NiokydPyufz2ZgQAAA4ha13ZAoKCpSTk6OUlBSdPXtWO3bs0P79+7Vnzx55PB4tX75c+fn5io+Pl9vt1po1a+Tz+fiOJQAAIMnmItPS0qLvf//7OnPmjDwej9LS0rRnzx7dc889kqTNmzdrxIgRys3NVTAYVHZ2tp577jk7IwMAAAextchs27atz/FRo0aptLRUpaWlQ5QIAACYxHFrZAAAAPqLIgMAAIxFkQEAAMaiyAAAAGNRZAAAgLEoMgAAwFgUGQAAYCyKDAAAMBZFBgAAGIsiAwAAjEWRAQAAxrL1Zy0BMEdGRob8fn+v416vV7W1tUOYCAAoMgD6ye/3q7i4uNfxwsLCIUwDAOfxaAkAABiLIgMAAIxFkQEAAMaiyAAAAGNRZAAAgLEoMgAAwFgUGQAAYCzeIwPAKJF4MR8v9wOGD4oMAKNE4sV8vNwPGD54tAQAAIxFkQEAAMaiyAAAAGNRZAAAgLEoMgAAwFgUGQAAYCyKDAAAMBZFBgAAGIsiAwAAjEWRAQAAxqLIAAAAY9laZEpKSjR37lzFxcVp0qRJWrp0qRobG8POmT9/vqKiosK2lStX2pQYAAA4ia1FpqqqSnl5eaqpqdFbb72lrq4uLVq0SO3t7WHnrVixQmfOnAltGzdutCkxAABwElt/+vXu3bvD9rdv365Jkyaprq5O8+bNCx0fM2aMvF7vUMcDAAAO56g1Mm1tbZKk+Pj4sOMvvfSSJkyYoNmzZ6ugoECff/55r9cIBoMKBAJhGwAAGJ5svSPzZT09PVq7dq3uvPNOzZ49O3T8O9/5jqZMmaKkpCTV19frscceU2Njo1555ZVLXqekpERFRUVDFRsAANjIMUUmLy9PDQ0N+p//+Z+w44888kjo17fccosSExOVmZmp48eP64YbbrjoOgUFBcrPzw/tBwIBJScnX73gAADANo4oMqtXr9Ybb7yhAwcOaPLkyX2em5GRIUk6duzYJYuMy+WSy+W6KjkBAICz2FpkLMvSmjVrtGvXLu3fv1+pqamX/czRo0clSYmJiVc5HQAAcDpbi0xeXp527Nih1157TXFxcfL7/ZIkj8ej0aNH6/jx49qxY4fuu+8+jR8/XvX19Vq3bp3mzZuntLQ0O6MDAAAHsLXIlJWVSTr/0rsvKy8v10MPPaSYmBjt3btXW7ZsUXt7u5KTk5Wbm6snnnjChrQAAMBpbH+01Jfk5GRVVVUNURoAAGAaR71HBgAAYCAoMgAAwFgUGQAAYCyKDAAAMBZFBgAAGIsiAwAAjEWRAQAAxqLIAAAAY1FkAACAsSgyAADAWIMqMtOmTdOnn3560fHW1lZNmzbtikMBAAD0x6CKzEcffaTu7u6LjgeDQZ06deqKQwEAAPTHgH5o5H/913+Ffr1nzx55PJ7Qfnd3tyorKzV16tSIhQMAAOjLgIrM0qVLJUlRUVFatmxZ2Fh0dLSmTp2qf/3Xf41YOAAAgL4MqMj09PRIklJTU3Xo0CFNmDDhqoQCAADojwEVmQtOnDgR6RwAAAADNqgiI0mVlZWqrKxUS0tL6E7NBb/+9a+vOBgAAMDlDKrIFBUVqbi4WF/72teUmJioqKioSOcCAAC4rEEVma1bt2r79u363ve+F+k8AAAA/Tao98h0dnbqjjvuiHQWAACAARlUkfm7v/s77dixI9JZAAAABmRQj5Y6Ojr0y1/+Unv37lVaWpqio6PDxjdt2hSRcAAAAH0ZVJGpr6/XrbfeKklqaGgIG2PhLwAAGCqDKjJvv/12pHMAAAAM2KDWyAAAADjBoO7ILFiwoM9HSPv27Rt0IAAAgP4aVJG5sD7mgq6uLh09elQNDQ0X/TBJAACAq2VQRWbz5s2XPP7UU0/p3LlzVxQIAACgvyK6RuZv/uZv+DlLAABgyES0yFRXV2vUqFGRvCQAAECvBvVo6YEHHgjbtyxLZ86c0eHDh/Xkk09GJBgAAMDlDKrIeDyesP0RI0ZoxowZKi4u1qJFiyISDAAA4HIGVWTKy8sjnQMAAGDArmiNTF1dnV588UW9+OKLOnLkyIA/X1JSorlz5youLk6TJk3S0qVL1djYGHZOR0eH8vLyNH78eMXGxio3N1fNzc1XEhsAAAwTgyoyLS0tWrhwoebOnatHH31Ujz76qNLT05WZmak//vGP/b5OVVWV8vLyVFNTo7feektdXV1atGiR2tvbQ+esW7dOr7/+unbu3KmqqiqdPn36ojU6AADg2jSoIrNmzRqdPXtW7733nj777DN99tlnamhoUCAQ0KOPPtrv6+zevVsPPfSQbr75Zs2ZM0fbt2/XyZMnVVdXJ0lqa2vTtm3btGnTJi1cuFDp6ekqLy/XwYMHVVNTM5joAABgGBnUGpndu3dr7969uummm0LHZs2apdLS0ita7NvW1iZJio+Pl3T+0VVXV5eysrJC58ycOVMpKSmqrq7W7bffftE1gsGggsFgaD8QCAw6DwAAcLZB3ZHp6elRdHT0Rcejo6PV09MzqCA9PT1au3at7rzzTs2ePVuS5Pf7FRMTo7Fjx4adm5CQIL/ff8nrlJSUyOPxhLbk5ORB5QEAAM43qCKzcOFC/fCHP9Tp06dDx06dOqV169YpMzNzUEHy8vLU0NCgioqKQX3+goKCArW1tYW2pqamK7oeAABwrkE9Wvq3f/s3/fVf/7WmTp0auuPR1NSk2bNn68UXXxzw9VavXq033nhDBw4c0OTJk0PHvV6vOjs71draGnZXprm5WV6v95LXcrlccrlcA84AAADMM6gik5ycrHfeeUd79+7VBx98IEm66aabwtay9IdlWVqzZo127dql/fv3KzU1NWw8PT1d0dHRqqysVG5uriSpsbFRJ0+elM/nG0x0AAAwjAyoyOzbt0+rV69WTU2N3G637rnnHt1zzz2Szi/Uvfnmm7V161bdfffd/bpeXl6eduzYoddee01xcXGhdS8ej0ejR4+Wx+PR8uXLlZ+fr/j4eLndbq1Zs0Y+n++SC30BAMC1ZUBrZLZs2aIVK1bI7XZfNObxePSDH/xAmzZt6vf1ysrK1NbWpvnz5ysxMTG0vfzyy6FzNm/erK9//evKzc3VvHnz5PV69corrwwkNgAAGKYGdEfmD3/4g372s5/1Or5o0SL9y7/8S7+vZ1nWZc8ZNWqUSktLVVpa2u/rAgCAa8OA7sg0Nzdf8tuuLxg5cuSA3uwLAABwJQZUZL7yla+ooaGh1/H6+nolJiZecSgAAID+GFCRue+++/Tkk0+qo6PjorEvvvhC69ev19e//vWIhQMAAOjLgNbIPPHEE3rllVf01a9+VatXr9aMGTMkSR988IFKS0vV3d2tn/70p1clKAAAwF8aUJFJSEjQwYMHtWrVKhUUFIQW60ZFRSk7O1ulpaVKSEi4KkEBAAD+0oBfiDdlyhT97ne/05/+9CcdO3ZMlmXpxhtv1Lhx465GPgAAgF4N6s2+kjRu3DjNnTs3klkAAAAGZFA/NBIAAMAJKDIAAMBYFBkAAGAsigwAADDWoBf7AsC1LCMjQ36/v9dxr9er2traIUwEXJsoMgAwCH6/X8XFxb2OFxYWDmEa4NrFoyUAAGAsigwAADAWRQYAABiLIgMAAIxFkQEAAMaiyAAAAGNRZAAAgLEoMgAAwFgUGQAAYCyKDAAAMBZFBgAAGIsiAwAAjEWRAQAAxqLIAAAAY1FkAACAsSgyAADAWBQZAABgLIoMAAAwFkUGAAAYiyIDAACMZWuROXDggBYvXqykpCRFRUXp1VdfDRt/6KGHFBUVFbbde++99oQFAACOY2uRaW9v15w5c1RaWtrrOffee6/OnDkT2n77298OYUIAAOBkI+384jk5OcrJyenzHJfLJa/XO0SJAACASRy/Rmb//v2aNGmSZsyYoVWrVunTTz/t8/xgMKhAIBC2AQCA4cnRRebee+/Vb37zG1VWVupnP/uZqqqqlJOTo+7u7l4/U1JSIo/HE9qSk5OHMDEAABhKtj5aupxvfetboV/fcsstSktL0w033KD9+/crMzPzkp8pKChQfn5+aD8QCFBmAAAYphx9R+YvTZs2TRMmTNCxY8d6PcflcsntdodtAABgeDKqyHzyySf69NNPlZiYaHcUAADgALY+Wjp37lzY3ZUTJ07o6NGjio+PV3x8vIqKipSbmyuv16vjx4/rJz/5iaZPn67s7GwbUwMAAKewtcgcPnxYCxYsCO1fWNuybNkylZWVqb6+Xi+88IJaW1uVlJSkRYsW6Z/+6Z/kcrnsigwAABzE1iIzf/58WZbV6/iePXuGMA0AADCNUWtkAAAAvowiAwAAjEWRAQAAxqLIAAAAY1FkAACAsSgyAADAWBQZAABgLIoMAAAwFkUGAAAYiyIDAACMRZEBAADGosgAAABjUWQAAICxKDIAAMBYFBkAAGAsigwAADAWRQYAABiLIgMAAIxFkQEAAMaiyAAAAGONtDsAgL5lZGTI7/f3Ou71elVbWzuEiQDAOSgygMP5/X4VFxf3Ol5YWDiEaQDAWXi0BAAAjEWRAQAAxqLIAAAAY1FkAACAsSgyAADAWBQZAABgLIoMAAAwFkUGAAAYiyIDAACMRZEBAADGosgAAABj2VpkDhw4oMWLFyspKUlRUVF69dVXw8Yty1JhYaESExM1evRoZWVl6cMPP7QnLAAAcBxbi0x7e7vmzJmj0tLSS45v3LhRzz77rLZu3ara2lpdf/31ys7OVkdHxxAnBQAATmTrT7/OyclRTk7OJccsy9KWLVv0xBNPaMmSJZKk3/zmN0pISNCrr76qb33rW0MZFQAAOJBj18icOHFCfr9fWVlZoWMej0cZGRmqrq7u9XPBYFCBQCBsAwAAw5Otd2T64vf7JUkJCQlhxxMSEkJjl1JSUqKioqKrmg3or4yMjD5/v3q9XtXW1g5hIgAYXhxbZAaroKBA+fn5of1AIKDk5GQbE+Fa5vf7VVxc3Ot4YWHhEKYBgOHHsY+WvF6vJKm5uTnseHNzc2jsUlwul9xud9gGAACGJ8cWmdTUVHm9XlVWVoaOBQIB1dbWyufz2ZgMAAA4ha2Pls6dO6djx46F9k+cOKGjR48qPj5eKSkpWrt2rf75n/9ZN954o1JTU/Xkk08qKSlJS5cutS80AABwDFuLzOHDh7VgwYLQ/oW1LcuWLdP27dv1k5/8RO3t7XrkkUfU2tqqu+66S7t379aoUaPsigwAABzE1iIzf/58WZbV63hUVJSKi4v7XCwJAACuXY5dIwMAAHA5FBkAAGCsYfceGQAX48V8AIYrigxwDeDFfACGKx4tAQAAY1FkAACAsSgyAADAWBQZAABgLIoMAAAwFkUGAAAYiyIDAACMxXtkAMAmvKgQuHIUGQCwCS8qBK4cj5YAAICxKDIAAMBYFBkAAGAsigwAADAWRQYAABiLIgMAAIxFkQEAAMaiyAAAAGNRZAAAgLEoMgAAwFgUGQAAYCyKDAAAMBZFBgAAGIsiAwAAjEWRAQAAxhppdwAAwOBlZGTI7/f3Ou71elVbWzuEiYChRZEBAIP5/X4VFxf3Ol5YWDiEaYChx6MlAABgLIoMAAAwFkUGAAAYy9FF5qmnnlJUVFTYNnPmTLtjAQAAh3D8Yt+bb75Ze/fuDe2PHOn4yAAAYIg4vhWMHDlSXq/X7hgAAMCBHP1oSZI+/PBDJSUladq0afrud7+rkydP9nl+MBhUIBAI2wAAwPDk6DsyGRkZ2r59u2bMmKEzZ86oqKhId999txoaGhQXF3fJz5SUlKioqGiIk2I44kVjAOB8ji4yOTk5oV+npaUpIyNDU6ZM0X/8x39o+fLll/xMQUGB8vPzQ/uBQEDJyclXPSuGH140BgDO5+gi85fGjh2rr371qzp27Fiv57hcLrlcriFMBQAA7OL4NTJfdu7cOR0/flyJiYl2RwEAAA7g6CLzD//wD6qqqtJHH32kgwcP6hvf+Iauu+46ffvb37Y7GgAAcABHP1r65JNP9O1vf1uffvqpJk6cqLvuuks1NTWaOHGi3dEAAIADOLrIVFRU2B0BAAA4mKMfLQEAAPSFIgMAAIxFkQEAAMaiyAAAAGNRZAAAgLEoMgAAwFgUGQAAYCyKDAAAMBZFBgAAGIsiAwAAjEWRAQAAxnL0z1oCAFx9GRkZ8vv9vY57vV7V1tZe9WsAg0GRAYBrnN/vV3Fxca/jhYWFQ3INYDB4tAQAAIxFkQEAAMaiyAAAAGNRZAAAgLEoMgAAwFgUGQAAYCyKDAAAMBbvkRkGLvciKunaexkVL+cCzMP/txgMiswwcLkXUUnX3suoeDkXYB7+v8Vg8GgJAAAYiyIDAACMRZEBAADGosgAAABjUWQAAICxKDIAAMBYfPv1FeD9LQDgLLyLJvKcPqcUmSvA+1sAwFl4F03kOX1OebQEAACMRZEBAADGosgAAABjGVFkSktLNXXqVI0aNUoZGRn6/e9/b3ckAADgAI4vMi+//LLy8/O1fv16vfPOO5ozZ46ys7PV0tJidzQAAGAzxxeZTZs2acWKFXr44Yc1a9Ysbd26VWPGjNGvf/1ru6MBAACbOfrbrzs7O1VXV6eCgoLQsREjRigrK0vV1dWX/EwwGFQwGAztt7W1SZICgUDE8/X09OiLL7647DlX42s7MYeTXG5O+jMfXMP+awQCUne3/Tm4xrV1DYSza04vXNOyrL5PtBzs1KlTliTr4MGDYcd//OMfW7fddtslP7N+/XpLEhsbGxsbG9sw2JqamvrsCo6+IzMYBQUFys/PD+339PTos88+0/jx4xUVFWVjsqETCASUnJyspqYmud1uu+MMC8xp5DGnkcecRh5zGnn9nVPLsnT27FklJSX1eT1HF5kJEybouuuuU3Nzc9jx5uZmeb3eS37G5XLJ5XKFHRs7duzViuhobreb//EijDmNPOY08pjTyGNOI68/c+rxeC57HUcv9o2JiVF6eroqKytDx3p6elRZWSmfz2djMgAA4ASOviMjSfn5+Vq2bJm+9rWv6bbbbtOWLVvU3t6uhx9+2O5oAADAZo4vMg8++KD++Mc/qrCwUH6/X7feeqt2796thIQEu6M5lsvl0vr16y96xIbBY04jjzmNPOY08pjTyIv0nEZZ1uW+rwkAAMCZHL1GBgAAoC8UGQAAYCyKDAAAMBZFBgAAGIsiM4x89NFHWr58uVJTUzV69GjdcMMNWr9+vTo7O8POq6+v1913361Ro0YpOTlZGzdutCmxGZ5++mndcccdGjNmTK8vVzx58qTuv/9+jRkzRpMmTdKPf/xj/fnPfx7aoAYpLS3V1KlTNWrUKGVkZOj3v/+93ZGMcuDAAS1evFhJSUmKiorSq6++GjZuWZYKCwuVmJio0aNHKysrSx9++KE9YQ1QUlKiuXPnKi4uTpMmTdLSpUvV2NgYdk5HR4fy8vI0fvx4xcbGKjc396KXteL/KysrU1paWuildz6fT2+++WZoPJLzSZEZRj744AP19PTo+eef13vvvafNmzdr69at+sd//MfQOYFAQIsWLdKUKVNUV1enn//853rqqaf0y1/+0sbkztbZ2alvfvObWrVq1SXHu7u7df/996uzs1MHDx7UCy+8oO3bt6uwsHCIk5rh5ZdfVn5+vtavX6933nlHc+bMUXZ2tlpaWuyOZoz29nbNmTNHpaWllxzfuHGjnn32WW3dulW1tbW6/vrrlZ2drY6OjiFOaoaqqirl5eWppqZGb731lrq6urRo0SK1t7eHzlm3bp1ef/117dy5U1VVVTp9+rQeeOABG1M72+TJk7VhwwbV1dXp8OHDWrhwoZYsWaL33ntPUoTnMyI/3RGOtXHjRis1NTW0/9xzz1njxo2zgsFg6Nhjjz1mzZgxw454RikvL7c8Hs9Fx3/3u99ZI0aMsPx+f+hYWVmZ5Xa7w+YZ5912221WXl5eaL+7u9tKSkqySkpKbExlLknWrl27Qvs9PT2W1+u1fv7zn4eOtba2Wi6Xy/rtb39rQ0LztLS0WJKsqqoqy7LOz190dLS1c+fO0Dn/+7//a0myqqur7YppnHHjxlm/+tWvIj6f3JEZ5tra2hQfHx/ar66u1rx58xQTExM6lp2drcbGRv3pT3+yI6Lxqqurdcstt4S9pDE7O1uBQCD0rw+c19nZqbq6OmVlZYWOjRgxQllZWaqurrYx2fBx4sQJ+f3+sDn2eDzKyMhgjvupra1NkkJ/dtbV1amrqytsTmfOnKmUlBTmtB+6u7tVUVGh9vZ2+Xy+iM8nRWYYO3bsmH7xi1/oBz/4QeiY3++/6K3IF/b9fv+Q5hsumNP++7//+z91d3dfcr6Yq8i4MI/M8eD09PRo7dq1uvPOOzV79mxJ5+c0JibmojVyzGnf3n33XcXGxsrlcmnlypXatWuXZs2aFfH5pMgY4PHHH1dUVFSf2wcffBD2mVOnTunee+/VN7/5Ta1YscKm5M41mDkFMPzl5eWpoaFBFRUVdkcx3owZM3T06FHV1tZq1apVWrZsmd5///2Ifx3H/6wlSD/60Y/00EMP9XnOtGnTQr8+ffq0FixYoDvuuOOiRbxer/eileEX9r1eb2QCG2Cgc9oXr9d70XfdXItz2h8TJkzQddddd8nfg8xVZFyYx+bmZiUmJoaONzc369Zbb7UplRlWr16tN954QwcOHNDkyZNDx71erzo7O9Xa2hp2F4Hft32LiYnR9OnTJUnp6ek6dOiQnnnmGT344IMRnU+KjAEmTpyoiRMn9uvcU6dOacGCBUpPT1d5eblGjAi/6ebz+fTTn/5UXV1dio6OliS99dZbmjFjhsaNGxfx7E41kDm9HJ/Pp6efflotLS2aNGmSpPNz6na7NWvWrIh8jeEiJiZG6enpqqys1NKlSyWdv5VfWVmp1atX2xtumEhNTZXX61VlZWWouAQCgdC/inExy7K0Zs0a7dq1S/v371dqamrYeHp6uqKjo1VZWanc3FxJUmNjo06ePCmfz2dHZCP19PQoGAxGfj4juCAZNvvkk0+s6dOnW5mZmdYnn3xinTlzJrRd0NraaiUkJFjf+973rIaGBquiosIaM2aM9fzzz9uY3Nk+/vhj68iRI1ZRUZEVGxtrHTlyxDpy5Ih19uxZy7Is689//rM1e/Zsa9GiRdbRo0et3bt3WxMnTrQKCgpsTu5MFRUVlsvlsrZv3269//771iOPPGKNHTs27Lu+0LezZ8+Gfh9KsjZt2mQdOXLE+vjjjy3LsqwNGzZYY8eOtV577TWrvr7eWrJkiZWammp98cUXNid3plWrVlkej8fav39/2J+bn3/+eeiclStXWikpKda+ffusw4cPWz6fz/L5fDamdrbHH3/cqqqqsk6cOGHV19dbjz/+uBUVFWX993//t2VZkZ1PiswwUl5ebkm65PZlf/jDH6y77rrLcrlc1le+8hVrw4YNNiU2w7Jlyy45p2+//XbonI8++sjKycmxRo8ebU2YMMH60Y9+ZHV1ddkX2uF+8YtfWCkpKVZMTIx12223WTU1NXZHMsrbb799yd+Ty5Ytsyzr/LdgP/nkk1ZCQoLlcrmszMxMq7Gx0d7QDtbbn5vl5eWhc7744gvr7//+761x48ZZY8aMsb7xjW+E/SMR4f72b//WmjJlihUTE2NNnDjRyszMDJUYy4rsfEZZlmUN9jYRAACAnfiuJQAAYCyKDAAAMBZFBgAAGIsiAwAAjEWRAQAAxqLIAAAAY1FkAACAsSgyAADAWBQZAABgLIoMAAAwFkUGAAAYiyIDAACM9f8ARAh0L4xuaiUAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "NHANES_sample.diff = NHANES_sample.BPSys1 - NHANES_sample.BPSys2\n", "\n", "sns.histplot(NHANES_sample.diff, bins=30, discrete=True, color='gray')\n", "plt.plot([NHANES_sample.diff.mean()] * 2, [0, 40], color='blue')" ] }, { "cell_type": "markdown", "id": "300aa26a", "metadata": {}, "source": [ "## Sign test\n" ] }, { "cell_type": "code", "execution_count": 91, "id": "49596473", "metadata": { "Rmd_chunk_options": "echo=FALSE", "jupyter": { "output_hidden": false }, "kernel": "R", "tags": [ "report_output" ] }, "outputs": [ { "data": { "text/plain": [ "BinomTestResult(k=96, n=200, alternative='two-sided', proportion_estimate=0.48, pvalue=0.6207289280968307)" ] }, "execution_count": 91, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# compute sign test for differences between first and second measurement\n", "npos = sum(NHANES_sample.diff > 0)\n", "binomtest(npos, NHANES_sample.shape[0], .5)\n" ] }, { "cell_type": "markdown", "id": "b8a28564", "metadata": {}, "source": [ "## Paired t-test\n" ] }, { "cell_type": "code", "execution_count": 82, "id": "a0b31206-8842-4998-9970-82884f5ebe01", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Tdofalternativep-valCI95%cohen-dBF10power
T-test1.840435199greater0.033597[0.09, inf]0.04660.8280.161561
\n", "
" ], "text/plain": [ " T dof alternative p-val CI95% cohen-d BF10 \\\n", "T-test 1.840435 199 greater 0.033597 [0.09, inf] 0.0466 0.828 \n", "\n", " power \n", "T-test 0.161561 " ] }, "execution_count": 82, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ttpaired = pg.ttest(x=NHANES_sample.BPSys1, \n", " y=NHANES_sample.BPSys2, paired=True,\n", " alternative='greater', correction=False)\n", "ttpaired" ] }, { "cell_type": "markdown", "id": "69fd2886", "metadata": {}, "source": [ "## Bayes factor result" ] }, { "cell_type": "code", "execution_count": 94, "id": "fc72bcd8", "metadata": { "Rmd_chunk_options": "echo=FALSE", "jupyter": { "output_hidden": false }, "kernel": "R", "tags": [ "report_output" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Bayes factor analysis\n", "--------------\n", "[1] Alt., r=0.707 : 0.4138236 ±0.05%\n", "\n", "Against denominator:\n", " Null, mu = 0 \n", "---\n", "Bayes factor type: BFoneSample, JZS\n", "\n" ] } ], "source": [ "%%R -i NHANES_sample\n", "\n", "# compute Bayes factor for paired t-test\n", "ttestBF(x = NHANES_sample$BPSys1, y = NHANES_sample$BPSys2, paired = TRUE)" ] }, { "cell_type": "markdown", "id": "14a5a076", "metadata": {}, "source": [ "## Figure 15.4" ] }, { "cell_type": "code", "execution_count": 117, "id": "f854fdc3-1656-44c2-a67d-f6f8eaad8188", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'systolic blood pressure')" ] }, "execution_count": 117, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGwCAYAAABPSaTdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA4sUlEQVR4nO3deViVdf7/8ddBdmPxuIQYZGOp6Cj51VTSErcADReYyYVJS9I2m5LSLscxbRszcU2LpnFcujLNNC2bkdwKMzMFaZp+jtu4pqiFCGggyP37w/HUGaThwIFzuHk+rutc17k/930+531z3XZefe7Pfd8WwzAMAQAAmJSHqwsAAACoSYQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgap6uLsAdlJWV6dSpUwoICJDFYnF1OQAAoBIMw1BBQYFCQ0Pl4VHx+A1hR9KpU6cUFhbm6jIAAEAVnDhxQjfddFOF6wk7kgICAiRd/WMFBga6uBoAAFAZ+fn5CgsLs/2OV4SwI9lOXQUGBhJ2AACoY/7XFBQmKAMAAFMj7AAAAFMj7AAAAFMj7AAAAFMj7AAAAFMj7AAAAFNzadjJyMhQfHy8QkNDZbFYtG7dOrv1Fovluq9Zs2bZtsnNzVVSUpICAwMVHBys5ORkFRYW1vKeAAAAd+XSsHPx4kVFRkZq0aJF111/+vRpu9df//pXWSwWJSYm2rZJSkrSt99+q02bNmnDhg3KyMjQuHHjamsXAACAm7MYhmG4ugjp6ijOBx98oCFDhlS4zZAhQ1RQUKAtW7ZIkvbt26d27dpp9+7d6tKliyRp48aNGjBggE6ePKnQ0NDr9lNcXKzi4mLb8rU7MF64cIGbCgIAUEfk5+crKCjof/5+15k5O2fOnNHHH3+s5ORkW9vOnTsVHBxsCzqS1K9fP3l4eGjXrl0V9jVjxgwFBQXZXjwXCwAA86ozYWfZsmUKCAhQQkKCrS0nJ0fNmjWz287T01NWq1U5OTkV9jV58mRduHDB9jpx4kSN1Q0AAFyrzjwb669//auSkpLk6+tb7b58fHzk4+PjhKoAAIC7qxNhZ/v27dq/f79WrVpl1x4SEqKzZ8/atZWWlio3N1chISG1WSIAwM0ZhqGioiKX13BtzqiPj8//fIBlTfP19XV5DbWhToSdxYsXq3PnzoqMjLRrj4qKUl5enjIzM9W5c2dJ0tatW1VWVqZu3bq5olQA/8EPS3n15YfFXRUVFSkmJsbVZbiV9PR0+fn5ubqMGufSsFNYWKhDhw7Zlo8cOaLs7GxZrVaFh4dLujrTevXq1Zo9e3a5z0dERCg2NlZjx45VWlqaSkpKNH78eA0fPrzCK7EA1A5+WMqrLz8sgLtxadjZs2ePevfubVtOSUmRJI0ePVpLly6VJK1cuVKGYWjEiBHX7eOdd97R+PHj1bdvX3l4eCgxMVELFiyo8doBAHWLr6+v0tPTXVpDUVGRBg8eLElav369U+ahVoerv7+2uM19dlypstfpA6g8dziN5Y4/LJzGqt9+/PFH24gnI33VV9nf7zoxZwdA3WOxWNzqP+S+vr5uVQ+A2lNn7rMDAABQFYQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgai4NOxkZGYqPj1doaKgsFovWrVtXbpt9+/Zp0KBBCgoKUsOGDXXHHXfo+PHjtvVFRUV6/PHH1bhxY91www1KTEzUmTNnanEvAACAO/N05ZdfvHhRkZGRGjNmjBISEsqtP3z4sHr27Knk5GQ9//zzCgwM1LfffitfX1/bNhMmTNDHH3+s1atXKygoSOPHj1dCQoJ27NhRm7viUoZhqKioyOU1FBcXS5J8fHxksVhcWo8k+fr6ukUdAADXcmnYiYuLU1xcXIXrp0yZogEDBujVV1+1tbVq1cr2/sKFC1q8eLFWrFihPn36SJKWLFmiiIgIffnll+revft1+y0uLrb9MEtSfn5+dXfFpYqKihQTE+PqMtxOenq6/Pz8XF0GAMDF3HbOTllZmT7++GO1bt1aMTExatasmbp162Z3qiszM1MlJSXq16+fra1t27YKDw/Xzp07K+x7xowZCgoKsr3CwsJqclcAAIALuXRk55ecPXtWhYWFeuWVV/TSSy9p5syZ2rhxoxISErRt2zb16tVLOTk58vb2VnBwsN1nb7zxRuXk5FTY9+TJk5WSkmJbzs/Pr9OBx9fXV+np6S6toaioSIMHD5YkrV+/3u5Uo6u4Qw0AANdz27BTVlYmSRo8eLAmTJggSbr99tv1xRdfKC0tTb169apy3z4+PvLx8XFKne7AYrG41ekaX19ft6oHAFC/ue1prCZNmsjT01Pt2rWza4+IiLBdjRUSEqLLly8rLy/PbpszZ84oJCSktkoFAABuzG3Djre3t+644w7t37/frv3AgQO6+eabJUmdO3eWl5eXtmzZYlu/f/9+HT9+XFFRUbVaLwAAcE8uPY1VWFioQ4cO2ZaPHDmi7OxsWa1WhYeHa+LEiRo2bJjuvvtu9e7dWxs3btRHH32kTz/9VJIUFBSk5ORkpaSkyGq1KjAwUE888YSioqIqvBILAADULy4NO3v27FHv3r1ty9cmDY8ePVpLly7V0KFDlZaWphkzZuj3v/+92rRpozVr1qhnz562z8ydO1ceHh5KTExUcXGxYmJi9Prrr9f6vgAAAPfk0rATHR0twzB+cZsxY8ZozJgxFa739fXVokWLtGjRImeXBwAATMBt5+wAAAA4A2EHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYWpXCTmlpqTZv3qw333xTBQUFkqRTp06psLDQqcUBAABUl6ejHzh27JhiY2N1/PhxFRcXq3///goICNDMmTNVXFystLS0mqgTAACgShwe2XnyySfVpUsXnT9/Xn5+frb2oUOHasuWLU4tDgAAoLocHtnZvn27vvjiC3l7e9u1t2zZUt99953TCgMAAHAGh0d2ysrKdOXKlXLtJ0+eVEBAgFOKAgAAcBaHw84999yjefPm2ZYtFosKCws1bdo0DRgwwJm1AQAAVJvDp7FSU1MVGxurdu3aqaioSCNHjtTBgwfVpEkTvfvuuzVRIwAAQJU5HHbCwsL09ddfa9WqVfr6669VWFio5ORkJSUl2U1YBgAAcAcOhZ2SkhK1bdtWGzZsUFJSkpKSkmqqLgAAAKdwaM6Ol5eXioqKaqoWAAAAp3N4gvLjjz+umTNnqrS0tNpfnpGRofj4eIWGhspisWjdunV26x944AFZLBa7V2xsrN02ubm5SkpKUmBgoIKDg5WcnMydnAEAgI3Dc3Z2796tLVu26JNPPlGHDh3UsGFDu/Vr166tdF8XL15UZGSkxowZo4SEhOtuExsbqyVLltiWfXx87NYnJSXp9OnT2rRpk0pKSvTggw9q3LhxWrFihQN7BQAAzMrhsBMcHKzExESnfHlcXJzi4uJ+cRsfHx+FhIRcd92+ffu0ceNG7d69W126dJEkvfbaaxowYIBSU1MVGhrqlDoBAEDd5XDY+fkoS2349NNP1axZMzVq1Eh9+vTRSy+9pMaNG0uSdu7cqeDgYFvQkaR+/frJw8NDu3bt0tChQ6/bZ3FxsYqLi23L+fn5NbsTQC0zDIP5dZLd34C/x1W+vr6yWCyuLgOoVQ6HndoUGxurhIQE3XLLLTp8+LD+8Ic/KC4uTjt37lSDBg2Uk5OjZs2a2X3G09NTVqtVOTk5FfY7Y8YMPf/88zVdPuAyRUVFiomJcXUZbmXw4MGuLsEtpKenc5sQ1DsOh51bbrnlF/+v4N///ne1Cvq54cOH29536NBBHTt2VKtWrfTpp5+qb9++Ve538uTJSklJsS3n5+crLCysWrUCACrGaONVjDaWVxujjQ6HnaeeespuuaSkRHv37tXGjRs1ceJEZ9V1Xb/61a/UpEkTHTp0SH379lVISIjOnj1rt01paalyc3MrnOcjXZ0H9N8TnQGzWnR3nnwaGK4uwyUMQ7pcdvW9t4dUX8/eFF+x6PGMYJfWwGhjeYw2XlUbo40Oh50nn3zyuu2LFi3Snj17ql3QLzl58qR++OEHNW/eXJIUFRWlvLw8ZWZmqnPnzpKkrVu3qqysTN26davRWoC6wqeBId8Grq7CdThhI0n1M+wC1zhtzk5cXJwmT57s0ATmwsJCHTp0yLZ85MgRZWdny2q1ymq16vnnn1diYqJCQkJ0+PBhTZo0Sbfeeqvt/w4iIiIUGxursWPHKi0tTSUlJRo/fryGDx/OlVgA4KauxF9x8xmjNciQdOU/7xtIqqejjSqVGnxUe/8X5rTD7f3335fVanXoM3v27FHv3r1ty9fm0YwePVpvvPGG/vGPf2jZsmXKy8tTaGio7rnnHr344ot2p6DeeecdjR8/Xn379pWHh4cSExO1YMEC5+wUAMD5PFV/w44kebm6gPrH4cOtU6dOdhOJDMNQTk6Ozp07p9dff92hvqKjo2UYFQ+vpqen/88+rFYrNxAEAAAVcjjsDBkyxG7Zw8NDTZs2VXR0tNq2beusugAAAJzC4bAzbdq0mqgDAACgRjj8INCsrCx98803tuX169dryJAh+sMf/qDLly87tTgAAIDqcjjsPPzwwzpw4ICkqzcQHDZsmPz9/bV69WpNmjTJ6QUCAABUh8Nh58CBA7r99tslSatXr1avXr20YsUKLV26VGvWrHF2fQAAANXicNgxDENlZVdvSbp582YNGDBAkhQWFqbvv//eudUBAABUk8Nhp0uXLnrppZf09ttv67PPPtPAgQMlXb0h4I033uj0AgEAAKrD4bAzb948ZWVlafz48ZoyZYpuvfVWSVdvKnjnnXc6vUAAAIDqcPjS844dO9pdjXXNrFmz1KBBPX4ADwAAcEsOj+ycOHFCJ0+etC1/9dVXeuqpp7R8+XJ5eXEPbAAA4F4cDjsjR47Utm3bJEk5OTnq37+/vvrqK02ZMkUvvPCC0wsEAACoDofDzj//+U917dpVkvTee+/p17/+tb744gu98847Wrp0qbPrAwAAqBaHw05JSYntqeObN2/WoEGDJElt27bV6dOnnVsdAABANTkcdtq3b6+0tDRt375dmzZtUmxsrCTp1KlTaty4sdMLBAAAqA6Hw87MmTP15ptvKjo6WiNGjFBkZKQk6cMPP7Sd3gIAAHAXDl96Hh0dre+//175+flq1KiRrX3cuHHy9/d3anEAAADV5fDIjnT1kRGZmZl68803VVBQIEny9vYm7AAAALfj8MjOsWPHFBsbq+PHj6u4uFj9+/dXQECAZs6cqeLiYqWlpdVEnQAAAFXi8MjOk08+qS5duuj8+fPy8/OztQ8dOlRbtmxxanEAAADV5fDIzvbt2/XFF1/I29vbrr1ly5b67rvvnFYYAACAMzg8slNWVqYrV66Uaz958qQCAgKcUhQAAICzOBx27rnnHs2bN8+2bLFYVFhYqGnTpmnAgAHOrA0AAKDaHD6NlZqaqtjYWLVr105FRUUaOXKkDh48qCZNmujdd9+tiRoBAACqzOGwExYWpq+//lqrVq3S119/rcLCQiUnJyspKcluwjIA1zEMw/a+uPxZZ9QzPz8Gfn5sAPWFQ2GnpKREbdu21YYNG5SUlKSkpKSaqgtANRQXF9veP57R6Be2RH1TXFzMPdFQ7zg0Z8fLy0tFRUU1VQsAAIDTOXwa6/HHH9fMmTP1l7/8RZ6eDn8cQC3w8fGxvV9093n5NHBhMXC54is/jfD9/NgA6guH08ru3bu1ZcsWffLJJ+rQoYMaNmxot37t2rVOKw5A1VgsFtt7nwaSL2EH//HzYwOoLxwOO8HBwUpMTKyJWgAAAJzO4bCzZMmSmqgDAACgRlR50s3Zs2e1f/9+SVKbNm3UrFkzpxUFAADgLA7fQTk/P1/333+/WrRooV69eqlXr15q0aKFfve73+nChQs1USMAAECVORx2xo4dq127dmnDhg3Ky8tTXl6eNmzYoD179ujhhx+uiRoBAACqzOHTWBs2bFB6erp69uxpa4uJidFbb72l2NhYpxYHAABQXQ6P7DRu3FhBQUHl2oOCgtSoEXdqBQAA7sXhkZ0//vGPSklJ0dtvv62QkBBJUk5OjiZOnKipU6c6vUAAQN1n90yuUtfVATfxs2OgNp7X5nDYeeONN3To0CGFh4crPDxcknT8+HH5+Pjo3LlzevPNN23bZmVlOa9SAECd9fPntTX4iLtc4ie18bw2h8POkCFDaqAMAACAmuFw2Jk2bVpN1AEAMLGfP5PrSvyVatzlDaZQ+tMIX208r43DDQBQ4+yeyeUpfn1gUxvPa3P4aiwAAIC6hLADAABMjbADAABMzaVhJyMjQ/Hx8QoNDZXFYtG6desq3PaRRx6RxWLRvHnz7Npzc3OVlJSkwMBABQcHKzk5WYWFhTVbOAAAqDMqNUUsJSWl0h3OmTOn0ttevHhRkZGRGjNmjBISEirc7oMPPtCXX36p0NDQcuuSkpJ0+vRpbdq0SSUlJXrwwQc1btw4rVixotJ1AAAA86pU2Nm7d6/dclZWlkpLS9WmTRtJ0oEDB9SgQQN17tzZoS+Pi4tTXFzcL27z3Xff6YknnlB6eroGDhxot27fvn3auHGjdu/erS5dukiSXnvtNQ0YMECpqanXDUcAAKB+qVTY2bZtm+39nDlzFBAQoGXLltmehXX+/Hk9+OCDuuuuu5xaXFlZme6//35NnDhR7du3L7d+586dCg4OtgUdSerXr588PDy0a9cuDR069Lr9FhcX293NMz8/36l1AwAA9+HwnJ3Zs2drxowZdg/9bNSokV566SXNnj3bqcXNnDlTnp6e+v3vf3/d9Tk5OWrWrJldm6enp6xWq3Jycirsd8aMGQoKCrK9wsLCnFo3AABwHw6Hnfz8fJ07d65c+7lz51RQUOCUoiQpMzNT8+fP19KlS51+w6HJkyfrwoULtteJEyec2j8AAHAfDoedoUOH6sEHH9TatWt18uRJnTx5UmvWrFFycvIvTjJ21Pbt23X27FmFh4fL09NTnp6eOnbsmJ5++mm1bNlSkhQSEqKzZ8/afa60tFS5ubm2J7Jfj4+PjwIDA+1eAADAnBy+YXdaWpqeeeYZjRw5UiUlJVc78fRUcnKyZs2a5bTC7r//fvXr18+uLSYmRvfff78efPBBSVJUVJTy8vKUmZlpmxy9detWlZWVqVu3bk6r5ZcYhqGioqJa+S539vO/AX+Pn/j6+tbKrdABABVzOOz4+/vr9ddf16xZs3T48GFJUqtWrdSwYUOHv7ywsFCHDh2yLR85ckTZ2dmyWq0KDw9X48aN7bb38vJSSEiI7SqwiIgIxcbGauzYsUpLS1NJSYnGjx+v4cOH19qVWEVFRYqJiamV76orBg8e7OoS3EZ6err8/PxcXQYA1GtVfhRbw4YNZbVabe+rYs+ePerdu7dt+dr9fEaPHq2lS5dWqo933nlH48ePV9++feXh4aHExEQtWLCgSvUAAADzcTjslJWV2a68unan4oCAAD399NOaMmWKPDwqPw0oOjpahmFUevujR4+Wa7NarW5zA8GL/5ckedTTR/kahlRWevW9h6dUn0/dlJWqYdY7rq4CAPAfDv8yT5kyRYsXL9Yrr7yiHj16SJI+//xzTZ8+XUVFRXr55ZedXmSd4eEpNfBydRUu5O3qAgAAKMfhsLNs2TL95S9/0aBBg2xtHTt2VIsWLfTYY4/V77ADAADcjsOXnufm5qpt27bl2tu2bavc3FynFAUAAOAsDoedyMhILVy4sFz7woULFRkZ6ZSiAAAAnMXh01ivvvqqBg4cqM2bNysqKkrS1WdUnThxQn/729+cXiAAAEB1ODyy06tXLx04cEBDhw5VXl6e8vLylJCQoP379zv9QaAAAADVVaXrpENDQ5mIDAAA6oQqhZ28vDwtXrxY+/btkyS1b99eY8aMUVBQkFOLAwAAqC6HT2Pt2bNHrVq10ty5c5Wbm6vc3FzNmTNHrVq1UlZWVk3UCAAAUGUOj+xMmDBBgwYN0ltvvSVPz6sfLy0t1UMPPaSnnnpKGRkZTi8SAACgqhwOO3v27LELOtLVp55PmjRJXbp0cWpxAAAA1eXwaazAwEAdP368XPuJEycUEBDglKIAAACcxeGwM2zYMCUnJ2vVqlU6ceKETpw4oZUrV+qhhx7SiBEjaqJGAACAKnP4NFZqaqosFotGjRql0tKrT7n28vLSo48+qldeecXpBQIAAFSHw2HH29tb8+fP14wZM3T48GFJUqtWreTv7+/04gAAAKqrSvfZkSR/f3916NDBmbUAAAA4XaXCTkJCQqU7XLt2bZWLAQAAcLZKhR3ujAwAAOqqSoWdJUuW1HQdAAAANaLKc3bOnj2r/fv3S5LatGmjZs2aOa0oAAAAZ3E47OTn5+vxxx/XypUrdeXKFUlSgwYNNGzYMC1atIhTXgCAX1bq6gJcyJB05T/vG0iyuLAWV6rlY8DhsDN27Fjt3btXGzZsUFRUlCRp586devLJJ/Xwww9r5cqVTi8SAGAeDT5q4OoSUM84HHY2bNig9PR09ezZ09YWExOjt956S7GxsU4tDgAAoLocDjuNGze+7qmqoKAgNWrUyClFAQDMxdfXV+np6a4uw+WKioo0ePBgSdL69evl6+vr4opcrzb+Bg6HnT/+8Y9KSUnR22+/rZCQEElSTk6OJk6cqKlTpzq9QADVU3zFoqsTBeofw5Aul1197+0hWerp/Iirx4BrWSwW+fn5uboMt+Lr68vfpJZUKux06tRJlp/9V+LgwYMKDw9XeHi4JOn48ePy8fHRuXPn9PDDD9dMpQCq5PGMYFeXAAAuVamwM2TIkBouAwAAoGZUKuxMmzatpusA4ETMj7iK+RHl8TdAfVTlmwoCcF/MjyiP+RFA/eXh6gIAAABqEiM71WQYP7vK5UqJ6wqB+/jZcWB3fAAAXIKwU03FxcW29w33rnBhJXBHxcXF8vf3d3UZAFCvcRoLAACYmsMjO4mJierataueffZZu/ZXX31Vu3fv1urVq51WXF3g4+Nje3+x00ipgZcLq4FbuFJiG+X7+fEBAHANh8NORkaGpk+fXq49Li5Os2fPdkZNdcrPb7aoBl6EHdix1Ndb9gKAG3H4NFZhYaG8vb3LtXt5eSk/P98pRQEAADiLw2GnQ4cOWrVqVbn2lStXql27dk4pCgAAwFkcPo01depUJSQk6PDhw+rTp48kacuWLXr33Xfr3XwdAADg/hwOO/Hx8Vq3bp3+9Kc/6f3335efn586duyozZs3q1evXjVRIwAAQJVV6T47AwcO1MCBA51dCwAAgNNxnx0AAGBqlRrZsVqtOnDggJo0aaJGjRr94uW0ubm5TisOAACguioVdubOnauAgABJ0rx582qyHgAAAKeqVNgZPXr0dd9XV0ZGhmbNmqXMzEydPn1aH3zwgYYMGWJbP336dK1cuVInTpyQt7e3OnfurJdfflndunWzbZObm6snnnhCH330kTw8PJSYmKj58+frhhtucFqdAACg7qrUnJ38/PxKvxxx8eJFRUZGatGiRddd37p1ay1cuFDffPONPv/8c7Vs2VL33HOPzp07Z9smKSlJ3377rTZt2qQNGzYoIyND48aNc6gOAABgXpUa2QkODv6ft703DEMWi0VXrlyp9JfHxcUpLi6uwvUjR460W54zZ44WL16sf/zjH+rbt6/27dunjRs3avfu3erSpYsk6bXXXtOAAQOUmpqq0NDQStcCAADMqVJhZ9u2bTVdx/90+fJl/fnPf1ZQUJAiIyMlSTt37lRwcLAt6EhSv3795OHhoV27dmno0KHX7au4uFjFxcW2ZR5zAQCAeVUq7LjyZoEbNmzQ8OHDdenSJTVv3lybNm1SkyZNJEk5OTlq1qyZ3faenp6yWq3KycmpsM8ZM2bo+eefr9G6AQCAe6jSTQXz8vK0ePFi7du3T5LUvn17jRkzRkFBQU4tTpJ69+6t7Oxsff/993rrrbd03333adeuXeVCjiMmT56slJQU23J+fr7CwsKcUS4AAHAzDt9UcM+ePWrVqpXmzp2r3Nxc5ebmas6cOWrVqpWysrKcXmDDhg116623qnv37lq8eLE8PT21ePFiSVJISIjOnj1rt31paalyc3MVEhJSYZ8+Pj4KDAy0ewEAAHNyOOxMmDBBgwYN0tGjR7V27VqtXbtWR44c0b333qunnnqqBkq0V1ZWZptvExUVpby8PGVmZtrWb926VWVlZXaXpwMAgPrL4dNYe/bs0VtvvSVPz58+6unpqUmTJtlNFK6MwsJCHTp0yLZ85MgRZWdny2q1qnHjxnr55Zc1aNAgNW/eXN9//70WLVqk7777Tr/97W8lSREREYqNjdXYsWOVlpamkpISjR8/XsOHD+dKLAAAIKkKIzuBgYE6fvx4ufYTJ07Y7rJcWXv27FGnTp3UqVMnSVJKSoo6deqk5557Tg0aNNC//vUvJSYmqnXr1oqPj9cPP/yg7du3q3379rY+3nnnHbVt21Z9+/bVgAED1LNnT/35z392dLcAAIBJOTyyM2zYMCUnJys1NVV33nmnJGnHjh2aOHGiRowY4VBf0dHRMgyjwvVr1679n31YrVatWLHCoe8FAAD1h8NhJzU1VRaLRaNGjVJpaakkycvLS48++qheeeUVpxcIAABQHQ6HHW9vb82fP18zZszQ4cOHJUmtWrWSv7+/04sDAACoLofn7IwZM0YFBQXy9/dXhw4d1KFDB/n7++vixYsaM2ZMTdQIAABQZQ6HnWXLlunHH38s1/7jjz9q+fLlTikKAADAWSp9Gis/P1+GYcgwDBUUFMjX19e27sqVK/rb3/5WrbsaAwAA1IRKh51rTz63WCxq3bp1ufUWi4XnTQEAALdT6bCzbds2GYahPn36aM2aNbJarbZ13t7euvnmm7mRHwAAcDuVDjvXnnx+5MgRhYeHy2Kx1FhRAAAAzuLwBOV9+/Zpx44dtuVFixbp9ttv18iRI3X+/HmnFgcAAFBdDoediRMnKj8/X5L0zTffKCUlRQMGDNCRI0eUkpLi9AIBAACqw+GbCh45ckTt2rWTJK1Zs0bx8fH605/+pKysLA0YMMDpBQIAAFSHwyM73t7eunTpkiRp8+bNuueeeyRdfUbVtREfAAAAd+HwyE7Pnj2VkpKiHj166KuvvtKqVaskSQcOHNBNN93k9AIBAACqw+GRnYULF8rT01Pvv/++3njjDbVo0UKS9Pe//12xsbFOLxAAAKA6HB7ZCQ8P14YNG8q1z5071ykFAQAAOJPDIzu9evXS8uXLr/t8LAAAAHfjcNjp1KmTnnnmGYWEhGjs2LH68ssva6IuAAAAp3A47MybN0+nTp3SkiVLdPbsWd19991q166dUlNTdebMmZqoEQAAoMocDjuS5OnpqYSEBK1fv14nT57UyJEjNXXqVIWFhWnIkCHaunWrs+sEAACokiqFnWu++uorTZs2TbNnz1azZs00efJkNWnSRPfee6+eeeYZZ9UIAABQZQ5fjXX27Fm9/fbbWrJkiQ4ePKj4+Hi9++67iomJsT0c9IEHHlBsbKxSU1OdXjAAAIAjHA47N910k1q1aqUxY8bogQceUNOmTctt07FjR91xxx1OKRAAAKA6HA47W7Zs0V133fWL2wQGBmrbtm1VLgoAAMBZHA47Xbp00aVLl+Tv7y9JOnbsmD744AO1a9fO9pysequs1NUVuI5h/LT/Hp7Sf05p1kv1+TgAADfkcNgZPHiwEhIS9MgjjygvL0/dunWTl5eXvv/+e82ZM0ePPvpoTdRZJzTMesfVJQAAgP/i8NVYWVlZttNY77//vm688UYdO3ZMy5cv14IFC5xeIAAAQHU4PLJz6dIlBQQESJI++eQTJSQkyMPDQ927d9exY8ecXqC78/X1VXp6uqvLcLmioiINHjxYkrR+/Xr5+vq6uCL3wN8BAFzP4bBz6623at26dRo6dKjS09M1YcIESVcvSQ8MDHR6ge7OYrHIz8/P1WW4FV9fX/4mAAC34fBprOeee07PPPOMWrZsqW7duikqKkrS1VGeTp06Ob1AAACA6nB4ZOc3v/mNevbsqdOnTysyMtLW3rdvXw0dOtSpxQEA4CyGYaioqMilNfz8+11di3R1JN5SD66edTjsSFJISIhCQkLs2rp27eqUggAAqAlFRUWKiYlxdRk21+Y5ulJ6enq9mHZQrWdjAQAAuLsqjewAAFDXuMPVs4ZhqLi4WJLk4+Pj8lNI9eWKUcIOgBrB/Ijy6sv8CHflLlfPXnsCAWoPYQdAjWB+RHn1ZX4E4G6YswMAAEyNkR0ANYL5EeXVl/kRgLsh7ACoEcyPAOAuOI0FAABMjbADAABMjbADAABMjbADAABMjbADAABMzaVhJyMjQ/Hx8QoNDZXFYtG6dets60pKSvTss8+qQ4cOatiwoUJDQzVq1CidOnXKro/c3FwlJSUpMDBQwcHBSk5OVmFhYS3vCQAAcFcuDTsXL15UZGSkFi1aVG7dpUuXlJWVpalTpyorK0tr167V/v37NWjQILvtkpKS9O2332rTpk3asGGDMjIyNG7cuNraBQAA4OZcep+duLg4xcXFXXddUFCQNm3aZNe2cOFCde3aVcePH1d4eLj27dunjRs3avfu3erSpYsk6bXXXtOAAQOUmpqq0NDQGt8HAADg3urUnJ0LFy7IYrEoODhYkrRz504FBwfbgo4k9evXTx4eHtq1a1eF/RQXFys/P9/uBQAAzKnOhJ2ioiI9++yzGjFihAIDAyVJOTk5atasmd12np6eslqtysnJqbCvGTNmKCgoyPYKCwur0doBAIDr1ImwU1JSovvuu0+GYeiNN96odn+TJ0/WhQsXbK8TJ044oUoAAOCO3P7ZWNeCzrFjx7R161bbqI4khYSE6OzZs3bbl5aWKjc3VyEhIRX26ePjIx8fnxqrGQAAuA+3Htm5FnQOHjyozZs3q3Hjxnbro6KilJeXp8zMTFvb1q1bVVZWpm7dutV2uQAAwA25dGSnsLBQhw4dsi0fOXJE2dnZslqtat68uX7zm98oKytLGzZs0JUrV2zzcKxWq7y9vRUREaHY2FiNHTtWaWlpKikp0fjx4zV8+HCuxAIAAJJcHHb27Nmj3r1725ZTUlIkSaNHj9b06dP14YcfSpJuv/12u89t27ZN0dHRkqR33nlH48ePV9++feXh4aHExEQtWLCgVuoHAADuz6VhJzo6WoZhVLj+l9ZdY7VatWLFCmeWBQAATMSt5+wAAABUF2EHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYGmEHAACYmkvDTkZGhuLj4xUaGiqLxaJ169bZrV+7dq3uueceNW7cWBaLRdnZ2eX6KCoq0uOPP67GjRvrhhtuUGJios6cOVM7OwAAANyeS8POxYsXFRkZqUWLFlW4vmfPnpo5c2aFfUyYMEEfffSRVq9erc8++0ynTp1SQkJCTZUMAADqGE9XfnlcXJzi4uIqXH///fdLko4ePXrd9RcuXNDixYu1YsUK9enTR5K0ZMkSRURE6Msvv1T37t2dXjMAAKhb6vScnczMTJWUlKhfv362trZt2yo8PFw7d+6s8HPFxcXKz8+3ewEAAHOq02EnJydH3t7eCg4Otmu/8cYblZOTU+HnZsyYoaCgINsrLCyshisFAACuUqfDTlVNnjxZFy5csL1OnDjh6pIAAEANcemcneoKCQnR5cuXlZeXZze6c+bMGYWEhFT4OR8fH/n4+NRChQAAwNXq9MhO586d5eXlpS1bttja9u/fr+PHjysqKsqFlQEAAHfh0pGdwsJCHTp0yLZ85MgRZWdny2q1Kjw8XLm5uTp+/LhOnTol6WqQka6O6ISEhCgoKEjJyclKSUmR1WpVYGCgnnjiCUVFRXElFgAAkOTikZ09e/aoU6dO6tSpkyQpJSVFnTp10nPPPSdJ+vDDD9WpUycNHDhQkjR8+HB16tRJaWlptj7mzp2re++9V4mJibr77rsVEhKitWvX1v7OAAAAt2QxDMNwdRGulp+fr6CgIF24cEGBgYGuLqdO+vHHHxUTEyNJSk9Pl5+fn4srAgCYXWV/v+v0nB0AAID/hbADAABMjbADAABMjbADAABMjbADAABMjbADAABMjbADAABMrU4/GwtXGYahoqIil9bw8+93dS3X+Pr6ymKxuLoMAICLEXZMoKioyHZDP3cwePBgV5cgiZsbAgCu4jQWAAAwNUZ2TMDX11fp6ekurcEwDBUXF0uSfHx83OL0ka+vr6tLAAC4AcKOCVgsFrc4XePv7+/qEgAAKIfTWAAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNR46rkkwzAkSfn5+S6uBAAAVNa13+1rv+MVIexIKigokCSFhYW5uBIAAOCogoICBQUFVbjeYvyvOFQPlJWV6dSpUwoICJDFYnF1OXVWfn6+wsLCdOLECQUGBrq6HEASxyXcD8ek8xiGoYKCAoWGhsrDo+KZOYzsSPLw8NBNN93k6jJMIzAwkH/AcDscl3A3HJPO8UsjOtcwQRkAAJgaYQcAAJgaYQdO4+Pjo2nTpsnHx8fVpQA2HJdwNxyTtY8JygAAwNQY2QEAAKZG2AEAAKZG2AEAAKZG2EGFoqOj9dRTT7m6DMCGYxKu1rJlS82bN6/Wvu/TTz+VxWJRXl5erX2nGRF24HKnT5/WyJEj1bp1a3l4ePBjBpdbu3at+vfvr6ZNmyowMFBRUVFKT093dVkAqoiwgyq5fPmy0/oqLi5W06ZN9cc//lGRkZFO6xf1izOPyYyMDPXv319/+9vflJmZqd69eys+Pl579+512ncAqD2EHUiSLl68qFGjRumGG25Q8+bNNXv2bLv1LVu21IsvvqhRo0YpMDBQ48aNu+7wanZ2tiwWi44ePWpre+uttxQWFiZ/f38NHTpUc+bMUXBwsF3f8+fP16hRoyp122/UD648JufNm6dJkybpjjvu0G233aY//elPuu222/TRRx/V8F7D1aKjozV+/HiNHz9eQUFBatKkiaZOnVrhU7XnzJmjDh06qGHDhgoLC9Njjz2mwsJCu2127Nih6Oho+fv7q1GjRoqJidH58+clXX0244wZM3TLLbfIz89PkZGRev/998t9z44dO9SxY0f5+vqqe/fu+uc//2m3fs2aNWrfvr18fHzUsmXLcv9e6jvCDiRJEydO1Geffab169frk08+0aeffqqsrCy7bVJTUxUZGam9e/dq6tSplep3x44deuSRR/Tkk08qOztb/fv318svv1wTuwCTcadjsqysTAUFBbJarVXeH9Qdy5Ytk6enp7766ivNnz9fc+bM0V/+8pfrbuvh4aEFCxbo22+/1bJly7R161ZNmjTJtj47O1t9+/ZVu3bttHPnTn3++eeKj4/XlStXJEkzZszQ8uXLlZaWpm+//VYTJkzQ7373O3322Wd23zNx4kTNnj1bu3fvVtOmTRUfH6+SkhJJUmZmpu677z4NHz5c33zzjaZPn66pU6dq6dKlNfMHqosM1HsFBQWGt7e38d5779nafvjhB8PPz8948sknDcMwjJtvvtkYMmSI3ee2bdtmSDLOnz9va9u7d68hyThy5IhhGIYxbNgwY+DAgXafS0pKMoKCgq5bS69evWzfifrLnY5JwzCMmTNnGo0aNTLOnDlTrf2C++vVq5cRERFhlJWV2dqeffZZIyIiwjCMq8fd3LlzK/z86tWrjcaNG9uWR4wYYfTo0eO62xYVFRn+/v7GF198YdeenJxsjBgxwjCMn47plStX2tZf+7ewatUqwzAMY+TIkUb//v3t+pg4caLRrl27Suxx/cDIDnT48GFdvnxZ3bp1s7VZrVa1adPGbrsuXbo43Pf+/fvVtWtXu7b/Xgb+mzsdkytWrNDzzz+v9957T82aNXP4+1D3dO/eXRaLxbYcFRWlgwcP2kZjfm7z5s3q27evWrRooYCAAN1///364YcfdOnSJUk/jexcz6FDh3Tp0iX1799fN9xwg+21fPlyHT582G7bqKgo2/tr/xb27dsnSdq3b5969Ohht32PHj0qrLk+8nR1Aag7GjZsaLfs4XE1Kxs/O5d9bVgVqA01fUyuXLlSDz30kFavXq1+/fpVuR+Y09GjR3Xvvffq0Ucf1csvvyyr1arPP/9cycnJunz5svz9/eXn51fh56/N7fn444/VokULu3U8N8u5GNmBWrVqJS8vL+3atcvWdv78eR04cOAXP9e0aVNJVy8dvyY7O9tumzZt2mj37t12bf+9DPw3dzgm3333XT344IN69913NXDgQEd3AXXYz487Sfryyy912223qUGDBnbtmZmZKisr0+zZs9W9e3e1bt1ap06dstumY8eO2rJly3W/p127dvLx8dHx48d166232r3CwsLK1XDNtX8LERERkqSIiAjt2LHDbvsdO3aodevW5WqurxjZgW644QYlJydr4sSJaty4sZo1a6YpU6bY/i+5Itf+QU6fPl0vv/yyDhw4UO4KgCeeeEJ333235syZo/j4eG3dulV///vf7YaIpZ9+kAoLC3Xu3DllZ2fL29tb7dq1c+q+om5w9TG5YsUKjR49WvPnz1e3bt2Uk5MjSfLz8+OKwXrg+PHjSklJ0cMPP6ysrCy99tpr17266dZbb1VJSYlee+01xcfHa8eOHUpLS7PbZvLkyerQoYMee+wxPfLII/L29ta2bdv029/+Vk2aNNEzzzyjCRMmqKysTD179tSFCxe0Y8cOBQYGavTo0bZ+XnjhBTVu3Fg33nijpkyZoiZNmmjIkCGSpKefflp33HGHXnzxRQ0bNkw7d+7UwoUL9frrr9fo36lOcfWkIbiHgoIC43e/+53h7+9v3Hjjjcarr75qN1m4okl5n3/+udGhQwfD19fXuOuuu4zVq1fbTQY1DMP485//bLRo0cLw8/MzhgwZYrz00ktGSEiIXT+Syr1uvvnmmtthuD1XHpO9evW67jE5evTomt1puFyvXr2Mxx57zHjkkUeMwMBAo1GjRsYf/vAH24Tl/z7u5syZYzRv3tzw8/MzYmJijOXLl5ebJP/pp58ad955p+Hj42MEBwcbMTExtvVlZWXGvHnzjDZt2hheXl5G06ZNjZiYGOOzzz4zDOOnCcofffSR0b59e8Pb29vo2rWr8fXXX9vV/f777xvt2rUzvLy8jPDwcGPWrFk1+neqayyGUcHNA4AaMnbsWP3rX//S9u3bXV0KIIljEj+Jjo7W7bffXquPhEDN4zQWalxqaqr69++vhg0b6u9//7uWLVvG8CpcimMSqF8IO6hxX331lV599VUVFBToV7/6lRYsWKCHHnrI1WWhHuOYBOoXTmMBAABT49JzAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdAKZw+fJlV5cAwE0RdgC4pYKCAiUlJalhw4Zq3ry55s6dq+joaD311FOSpJYtW+rFF1/UqFGjFBgYqHHjxkmS1qxZo/bt28vHx0ctW7Ys97Rqi8WidevW2bUFBwdr6dKlkqSjR4/KYrFo5cqVuvPOO+Xr66tf//rX+uyzz2p6lwHUEMIOALeUkpKiHTt26MMPP9SmTZu0fft2ZWVl2W2TmpqqyMhI7d27V1OnTlVmZqbuu+8+DR8+XN98842mT5+uqVOn2oKMIyZOnKinn35ae/fuVVRUlOLj4/XDDz84ae8A1CaejQXA7RQUFGjZsmVasWKF+vbtK0lasmSJQkND7bbr06ePnn76adtyUlKS+vbtq6lTp0qSWrdurf/3//6fZs2apQceeMChGsaPH6/ExERJ0htvvKGNGzdq8eLFmjRpUjX2DIArMLIDwO38+9//VklJibp27WprCwoKUps2bey269Kli93yvn371KNHD7u2Hj166ODBg7py5YpDNURFRdnee3p6qkuXLtq3b59DfQBwD4QdAHVWw4YNHf6MxWLRfz//uKSkxFklAXBDhB0AbudXv/qVvLy8tHv3blvbhQsXdODAgV/8XEREhHbs2GHXtmPHDrVu3VoNGjSQJDVt2lSnT5+2rT948KAuXbpUrq8vv/zS9r60tFSZmZmKiIio0v4AcC3m7ABwOwEBARo9erQmTpwoq9WqZs2aadq0afLw8JDFYqnwc08//bTuuOMOvfjiixo2bJh27typhQsX6vXXX7dt06dPHy1cuFBRUVG6cuWKnn32WXl5eZXra9GiRbrtttsUERGhuXPn6vz58xozZkyN7C+AmsXIDgC3NGfOHEVFRenee+9Vv3791KNHD0VERMjX17fCz/zf//2f3nvvPa1cuVK//vWv9dxzz+mFF16wm5w8e/ZshYWF6a677tLIkSP1zDPPyN/fv1xfr7zyil555RVFRkbq888/14cffqgmTZrUxK4CqGEW479PXgOAG7p48aJatGih2bNnKzk5uca+5+jRo7rlllu0d+9e3X777TX2PQBqD6exALilvXv36l//+pe6du2qCxcu6IUXXpAkDR482MWVAahrCDsA3FZqaqr2798vb29vde7cWdu3b+dUEgCHcRoLAACYGhOUAQCAqRF2AACAqRF2AACAqRF2AACAqRF2AACAqRF2AACAqRF2AACAqRF2AACAqf1/I1xyAwvQIAkAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "nPerGroup = 36\n", "noiseSD = 10\n", "meanSysBP = 140\n", "effectSize = 0.8\n", "DrugDf = pd.DataFrame({'group': ['drug1', 'drug2', 'placebo'] * nPerGroup})\n", "\n", "DrugDf['SysBP'] = rng.normal(loc=meanSysBP, scale=noiseSD, size=DrugDf.shape[0]) - (DrugDf.group == 'drug1') * effectSize * noiseSD\n", "\n", "sns.boxplot(data=DrugDf, x='group', y='SysBP')\n", "plt.ylabel('systolic blood pressure')" ] }, { "cell_type": "markdown", "id": "add73229", "metadata": {}, "source": [ "## Figure 15.5" ] }, { "cell_type": "code", "execution_count": 112, "id": "8ea003ea-ada4-42a5-b35c-4b092bd15bb7", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'density')" ] }, "execution_count": 112, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABmMUlEQVR4nO3deXxU9b3/8deZNftGyAYBIjuyg1JQqyiKS12virtSq/daaFV+XaQVqNaKWrVai1K8tdpeFa211qrFIoorIquiIggSw5KEJWRfZ+b8/pjMJIEkJJNZksn7+bjnMTNnzvJJ7DVvv9/v+X4N0zRNRERERKKEJdIFiIiIiASTwo2IiIhEFYUbERERiSoKNyIiIhJVFG5EREQkqijciIiISFRRuBEREZGoYot0AeHm8XjYt28fiYmJGIYR6XJERESkA0zTpKKigpycHCyW9ttmel242bdvH7m5uZEuQ0RERAKwe/du+vfv3+4xvS7cJCYmAt5fTlJSUoSrERERkY4oLy8nNzfX/3e8Pb0u3Pi6opKSkhRuREREepiODCnRgGIRERGJKgo3IiIiElUUbkRERCSq9LoxNyIiIgBut5uGhoZIlyHNOByOYz7m3REKNyIi0quYpklRURGlpaWRLkWOYLFYyMvLw+FwdOk6CjciItKr+IJNRkYGcXFxmtC1m/BNsltYWMiAAQO69M9F4UZERHoNt9vtDzZ9+vSJdDlyhL59+7Jv3z5cLhd2uz3g62hAsYiI9Bq+MTZxcXERrkRa4+uOcrvdXbpORMPNe++9x/nnn09OTg6GYfDKK6+0e/zLL7/MmWeeSd++fUlKSmLq1Km8+eab4SlWRESihrqiuqdg/XOJaLipqqpi3LhxLFmypEPHv/fee5x55pm88cYbbNiwgenTp3P++eezadOmEFcqIiIiPUVEx9ycc845nHPOOR0+/pFHHmnx+d577+Wf//wn//rXv5gwYUKQqxMREZGeqEePufF4PFRUVJCWltbmMXV1dZSXl7fYREREehrTNLn55ptJS0vDMAw2b94c6ZK6rR4dbh588EEqKyu5/PLL2zxm8eLFJCcn+7fc3NzQFOOqh7K9UFoQmuuLiEivtmLFCp5++mlee+01CgsLGT16dJvHdnZM65FefvllzjrrLPr06dMjg1SPDTfPPfccd911Fy+++CIZGRltHjd//nzKysr82+7du0NT0J518LtR8NeLQ3N9ERHp1Xbu3El2djbTpk0jKysLm63tkSWdHdPa2vknn3wy999/f6DlRlSPnOdm+fLl/OAHP+Bvf/sbM2bMaPdYp9OJ0+kMfVGOxscK66tDfy8REQka0zSpaejao8eBirVbO/SE0A033MAzzzwDeJ8oGjhwIPn5+W0e39kxrUe69tprAdq9R3fW48LN888/z/e//32WL1/OeeedF+lymtjjva8NCjciIj1JTYObUQsjM63Il3fPJM5x7D/Fjz76KIMHD2bZsmWsW7cOq9Uahup6roiGm8rKSnbs2OH/vGvXLjZv3kxaWhoDBgxg/vz57N27l7/85S+Atyvq+uuv59FHH2XKlCkUFRUBEBsbS3JyckR+Bj97rPdV4UZERIIsOTmZxMRErFYrWVlZkS6n24touFm/fj3Tp0/3f543bx4A119/PU8//TSFhYUUFDQN0F22bBkul4s5c+YwZ84c/37f8RHlaGy5cdeD2wXWHtcoJiLSK8XarXx598yI3VuCL6J/gU877TRM02zz+yMDy+rVq0NbUFfYm03l3VAN1qTI1SIiIh1mGEaHuoak5+ixT0t1OzYn0DgoTF1TIiIiEaOoGiyG4e2aqq+E+qpIVyMiIr3Ysca0HktJSQkFBQXs27cPgG3btgGQlZXVI8b8qOUmmHxdUw01ka1DRER6tfXr1zNhwgT/0kTz5s1jwoQJLFy4sEPnv/rqq0yYMMH/VPIVV1zBhAkTWLp0achqDibDbG/QSxQqLy8nOTmZsrIykpKCPC7mkbFQ+i3cuBJyTwzutUVEpMtqa2vZtWsXeXl5xMTERLocOUJ7/3w68/dbLTfB5HtiSt1SIiIiEaNwE0z+bikNKBYRkdAoKCggISGhza35FCqtef/999s9PxpoQHEwOTTmRkREQisnJ6fdhSxzcnLaPX/y5Mk9biHMzlK4CSZfy426pUREJERsNhtDhgwJ+PzY2Ngund8TqFsqmNQtJSIiEnEKN8HkULgRERGJNIWbYPJ3SynciIiIRIrCTTCpW0pERCTiFG6CyTfPjcKNiIhIxCjcBJM91vuqbikREQky0zS5+eabSUtLwzCMqH+cuysUboJJ3VIiIhIiK1as4Omnn+a1116jsLCQ0aNHt3nse++9x/nnn09OTg6GYfDKK6906l433HADhmG02M4+++wWx5SUlHD11VeTlJRESkoKN954I5WVlYH8aEGncBNM6pYSEZEQ2blzJ9nZ2UybNo2srCxstranqquqqmLcuHEsWbIk4PudffbZFBYW+rfnn3++xfdXX301X3zxBStXruS1117jvffe4+abbw74fsGkSfyCSd1SIiI9j2lG7j9K7XFgGMc87IYbbuCZZ54BwDAMBg4cSH5+fpvHn3POOZxzzjldKs3pdJKVldXqd1u3bmXFihWsW7eOyZMnA/DYY49x7rnn8uCDDx5zluRQU7gJJruv5UYzFIuI9BgN1XBvhP4Y/2JfU6t/Ox599FEGDx7MsmXLWLduHVarNeSlrV69moyMDFJTUzn99NO555576NOnDwBr1qwhJSXFH2wAZsyYgcViYe3atVx88cUhr689CjfB5NA8NyIiEnzJyckkJiZitVrbbE0JprPPPptLLrmEvLw8du7cyS9+8QvOOecc1qxZg9VqpaioiIyMjBbn2Gw20tLSKCoqCnl9x6JwE0y+biktnCki0nPY47wtKJG6dzd0xRVX+N+PGTOGsWPHMnjwYFavXs0ZZ5wRwco6RuEmmNQtJSLS8xhGh7qGerPjjjuO9PR0duzYwRlnnEFWVhb79+9vcYzL5aKkpCQsLUvHoqelgkndUiIiEoX27NnDoUOHyM7OBmDq1KmUlpayYcMG/zFvv/02Ho+HKVOmRKpMP7XcBJOvedHTAO4GsNojW4+IiPRKlZWV7Nixw/95165dbN68mbS0NAYMGHDMc++66y7+67/+i6ysLHbu3MnPfvYzhgwZwsyZMwEYOXIkZ599NjfddBNLly6loaGBuXPncsUVV0T8SSlQy01wNe871Vw3IiISIevXr2fChAlMmDABgHnz5jFhwgQWLlx4zHOtViufffYZF1xwAcOGDePGG29k0qRJvP/++zidTv9xzz77LCNGjOCMM87g3HPP5eSTT2bZsmUh+5k6wzBN04x0EeFUXl5OcnIyZWVlJCUlBffipgl3p4HpgXlfQVJ2cK8vIiJdUltby65du8jLyyMmJibS5cgR2vvn05m/32q5CSbDaDaoWC03IiIikaBwE2z+x8EVbkREJPgKCgpISEhocysoKGj3/Pfff7/d86OBBhQHmyMOqtATUyIiEhI5OTntrgh+rAG9kydPjvoVxRVugk3dUiIiEkI2m40hQ4YEfH5sbGyXzu8J1C0VbOqWEhERiSiFm2DTRH4iIiIRpXATbFqCQUREJKIUboJNi2eKiIhElMJNsPm7pdRyIyIiEgkKN8Gmp6VERCQETNPk5ptvJi0tDcMwov5x7q5QuAk2dUuJiEgIrFixgqeffprXXnuNwsJCRo8e3eax7733Hueffz45OTkYhsErr7xy1DGmabJw4UKys7OJjY1lxowZfP311x2uZ9CgQRiG0WK77777Whzz2WefccoppxATE0Nubi4PPPBAh6/fFQo3weZobLlRt5SIiATRzp07yc7OZtq0aWRlZWGztT1VXVVVFePGjWPJkiVtHvPAAw/w+9//nqVLl7J27Vri4+OZOXMmtbW1Ha7p7rvvprCw0L/96Ec/8n9XXl7OWWedxcCBA9mwYQO//e1v+dWvfhWWxTU1iV+w+VYGV7eUiEiPYJomNa7ItLbH2mIxDOOYx91www0888wzABiGwcCBA8nPz2/z+HPOOYdzzjmnze9N0+SRRx7hzjvv5MILLwTgL3/5C5mZmbzyyitcccUVHao/MTGRrKysVr979tlnqa+v56mnnsLhcHD88cezefNmHn74YW6++eYOXT9QCjfB5htQrG4pEZEeocZVw5TnpkTk3muvWkuc7z+K2/Hoo48yePBgli1bxrp167BarV26765duygqKmLGjBn+fcnJyUyZMoU1a9Z0ONzcd999/PrXv2bAgAFcddVV3H777f4WpTVr1vDd734Xh8PhP37mzJncf//9HD58mNTU1C79DO1RuAk2u56WEhGR4EpOTiYxMRGr1dpmS0lnFBUVAZCZmdlif2Zmpv+7Y/nxj3/MxIkTSUtL46OPPmL+/PkUFhby8MMP+++Rl5d31PV93ync9CTqlhIR6VFibbGsvWptxO7dU82bN8//fuzYsTgcDv77v/+bxYsX43Q6I1iZwk3wORRuRER6EsMwOtQ1FE18rT/FxcVkZ2f79xcXFzN+/PiArjllyhRcLhf5+fkMHz6crKwsiouLWxzj+xyM1qf26GmpYLNrbSkREene8vLyyMrKYtWqVf595eXlrF27lqlTpwZ0zc2bN2OxWMjIyABg6tSpvPfeezQ0NPiPWblyJcOHDw9plxSo5Sb41C0lIiIRVllZyY4dO/yfd+3axebNm0lLS2PAgAEYhsFtt93GPffcw9ChQ8nLy2PBggXk5ORw0UUXHfP6a9asYe3atUyfPp3ExETWrFnD7bffzjXXXOMPLldddRV33XUXN954Iz//+c/5/PPPefTRR/nd734Xqh/bT+Em2Pzz3CjciIhIZKxfv57p06f7P/vGx1x//fU8/fTTAPzsZz+jqqqKm2++mdLSUk4++WRWrFhBTEzMMa/vdDpZvnw5v/rVr6irqyMvL4/bb7+9xTic5ORk/vOf/zBnzhwmTZpEeno6CxcuDPlj4ACGaZpmyO/SjZSXl5OcnExZWRlJSUnBv0HZHvjd8WCxw8KDwb++iIgErLa2ll27dpGXl9ehP+ISXu398+nM3++IjrnpyPTQR1q9ejUTJ07E6XQyZMgQfwLtNnzdUp4GcDe0f6yIiIgEXUTDTUemh25u165dnHfeeUyfPp3Nmzdz22238YMf/IA333wzxJV2gq9bCjTXjYiIBF1BQQEJCQltbgUFBV2+x7333tvm9dub+bi7iOiYm2NND32kpUuXkpeXx0MPPQTAyJEj+eCDD/jd737HzJkzWz2nrq6Ouro6/+fy8vKuFX0sVgcYFjA93lmKY1NCez8REelVcnJy2l0RPCcnp8v3+J//+R8uv/zyVr+Lje3+c/P0qAHFa9asaTFVNHincr7tttvaPGfx4sXcddddIa6sGcMAezzUV+iJKRERCTqbzcaQIUNCeo+0tDTS0tJCeo9Q6lHz3BQVFbU6VXR5eTk1Na2v5TR//nzKysr82+7du0NfqENLMIiIdGcejyfSJUgrgvWMU49quQmE0+kM/zTQ9sYmOy2eKSLSrTgcDiwWC/v27aNv3744HI4OrcotoWeaJgcOHMAwDOx2e5eu1aPCTVtTOSclJXWvPkB746DiBrXciIh0JxaLhby8PAoLC9m3b1+ky5EjGIZB//79u7zqeY8KN1OnTuWNN95osW/lypUBTxUdMg4twSAi0l05HA4GDBiAy+XC7XZHuhxpxm63dznYQITDzbGmh54/fz579+7lL3/5C+Advf2HP/yBn/3sZ3z/+9/n7bff5sUXX+T111+P1I/QOv8SDOqWEhHpjnxdH13t/pDuKaIDitevX8+ECROYMGEC4J0eesKECSxcuBCAwsLCFs/r5+Xl8frrr7Ny5UrGjRvHQw89xP/+7/+2+Rh4xPjDjbqlREREwi2iLTennXZauyOjW5t9+LTTTmPTpk0hrCoI1C0lIiISMT3qUfAeQy03IiIiEaNwEwoacyMiIhIxCjehoG4pERGRiFG4CQXNcyMiIhIxCjehoBmKRUREIkbhJhTULSUiIhIxCjehoG4pERGRiFG4CQV1S4mIiESMwk0oqFtKREQkYhRuQkHdUiIiIhGjcBMK6pYSERGJGIWbUHA0ttyoW0pERCTsFG5CQWtLiYiIRIzCTSj4uqU8LnDVR7YWERGRXkbhJhR83VIADeqaEhERCSeFmyD5dHcp0x9czaw/rgGrAwyr9wuFGxERkbCyRbqAaLLrYBV1DW4wDG/rTV25BhWLiIiEmVpugiTe6c2JlXUu7w7/4+DhCzcbijewef/msN1PRESkO1LLTZAkNIabqno3pmli+J+YCn24Kasr475P7uO1b17Dalj550X/ZGDSwJDfV0REpDtSy02QxDu9Y2zcHpM6l6fZXDehfRz8w70fcsmrl/DaN69572+6eXbrsyG9p4iISHemcBMk8Y6mRrDKOlfIZymucdVw95q7+Z+3/of91fsZlDSI2yfdDsArO16hvL48JPcVERHp7hRugsRiMYhzeFtvqupczSbyC0231BOfPsHftv8NgKtHXs2L57/I7ONnMzR1KDWuGv6+/e8hua+IiEh3p3ATRC0GFYe4W+qjvR8BsOA7C7jjxDuItcViGAbXjrwWgOe+eg6XxxWSe4uIiHRnCjdB5B9UXOcOabdUVUMVX5d+DcBpuae1+O7c484lLSaNoqoi3vr2raDfW0REpLtTuAki36Dilt1SwW+5+fzg53hMD9nx2WTEZbT4zml1Mmv4LAD++uVfg35vERGR7k7hJoh8g4pbdksFf8zNpwc+BWBc33Gtfn/58MuxW+x8dvAzzXsjIiK9jsJNEDV1S7lCOonfscJNemw65x13HqDWGxER6X0UboKoxYBie2PLTZDDjWmafHbgM6DtcANw7SjvwOK3Ct5iX+W+oNYgIiLSnSncBFF88wHFjsYxN0Hulvq2/FtK60pxWp2MSBvR5nHDUocxJXsKHtPD8m3Lg1qDiIhId6ZwE0QJvgHF9aHrlvJ1SY3qMwq71d7usZcOuxTwzmIsIiLSWyjcBFE4uqWONd6mucmZkwH4+vDXmrFYRER6DYWbIGoxoDhE3VKdCTfpsekMTBqIiamnpkREpNdQuAmi+BZPSwV/+YWqhip2lO4AOhZuACZkTABgY/HGoNUhIiLSnSncBFHLbqnghxvf5H058Tn0jevboXMmZkwEYNP+TUGrQ0REpDtTuAki/4DiED0t1ZkuKZ+Jmd5ws+XgFurcdUGrRUREpLtSuAki3wzFVSEaUOwLN2P7ju3wOQMSB9Anpg8Nnga+OPhF0GoRERHprhRugqhlt1RwHwXv6OR9RzIMw996s3G/xt2IiEj0U7gJolaflvK4wFXf5Wt3dPK+1vjG3Wwo3tDlOkRERLo7hZsg8j8tVe/GY4tr+iIIK4N3ZvK+I/labj7d/yluj7vLtYiIiHRnCjdB5Gu5Aah2W8DwDjAOxqDiQAYT+wxLHUacLY6Khgr/o+QiIiLRSuEmiGLsFiyG931VvRscvkHFNV2+dlfCjc1iY3zGeEDjbkREJPop3ASRYRhtzHXTtW6pQCbvO5Im8xMRkd5C4SbIWgwq9j0x1cVuqe2Ht+MxPWTGZXZ48r4jTcqcBHjDjWmaXapHRESkO1O4CbIWLTeO4Mx1U1hZCEBuYm7A1xidPhqbxcb+mv3srdzbpXpERES6M4WbIGtaX8odtCUYCqu84SY7Pjvga8TaYhnVZxSgcTciIhLdIh5ulixZwqBBg4iJiWHKlCl88skn7R7/yCOPMHz4cGJjY8nNzeX222+ntrY2TNUeW9MSDMHrliqqKgIgKz6rS9fxzXejcTciIhLNIhpuXnjhBebNm8eiRYvYuHEj48aNY+bMmezfv7/V45977jnuuOMOFi1axNatW/nTn/7ECy+8wC9+8YswV9423xIMweyWCnq4UcuNiIhEsYiGm4cffpibbrqJ2bNnM2rUKJYuXUpcXBxPPfVUq8d/9NFHnHTSSVx11VUMGjSIs846iyuvvPKYrT3h1HJAcXC7pboabnxPTO0q20VJbUmXriUiItJdRSzc1NfXs2HDBmbMmNFUjMXCjBkzWLNmTavnTJs2jQ0bNvjDzDfffMMbb7zBueee2+Z96urqKC8vb7GFUnzzcBOT5N1ZU9qlaxZVe1tuujLmBiAlJoXByYMB/OtUiYiIRJuIhZuDBw/idrvJzMxssT8zM5OioqJWz7nqqqu4++67Ofnkk7Hb7QwePJjTTjut3W6pxYsXk5yc7N9ycwN/4qgjmp6WckN8hndnVevdbB1R3VBNWV0Z0PVwA/gHFX9V8lWXryUiItIdRXxAcWesXr2ae++9l8cff5yNGzfy8ssv8/rrr/PrX/+6zXPmz59PWVmZf9u9e3dIa2wxoDihcU6aygMBX8/XapNgTyDBkdDl+oanDQdgW8m2Ll9LRESkO7Id+5DQSE9Px2q1Ulxc3GJ/cXExWVmtjy1ZsGAB1157LT/4wQ8AGDNmDFVVVdx888388pe/xGI5Oqs5nU6cTmfwf4A2+Ftu6l1BabkpqgzOYGIf34riarkREZFoFbGWG4fDwaRJk1i1apV/n8fjYdWqVUydOrXVc6qrq48KMFart6Wku8y622LMTUJjl1tlcTtntM/XchPscLOncg8V9RVBuaaIiEh3EtFuqXnz5vHkk0/yzDPPsHXrVm655RaqqqqYPXs2ANdddx3z58/3H3/++efzxBNPsHz5cnbt2sXKlStZsGAB559/vj/kRFqLp6Wad0sFGL6CMYFfc8nOZP+11DUlIiLRKGLdUgCzZs3iwIEDLFy4kKKiIsaPH8+KFSv8g4wLCgpatNTceeedGIbBnXfeyd69e+nbty/nn38+v/nNbyL1Ixyl1QHFrhqorwRnYqev51t6IVgtN+Add1NYVci2w9uYnDU5aNcVERHpDiIabgDmzp3L3LlzW/1u9erVLT7bbDYWLVrEokWLwlBZYFoMKHYmgD3euyp45f6Awk2wHgNvbkTaCFbvXq1xNyIiEpV61NNSPUGLMTfQrGsqsEHFwZqduLkRqd5xN+qWEhGRaKRwE2Qtll+ALj0xZZpmaMJNH2+4+br0axrcDUG7roiISHegcBNkvgHFdS4PLrcHEhrDTQAtN4frDlPnrsPAIDMu89gndFBOfA6J9kRcHhfflH0TtOuKiIh0Bwo3QebrlgKoqnN3Kdz4Wm36xPbBYXUEpT4AwzD8k/lp3I2IiEQbhZsgc9gsOKzeX2tXJ/IL9mPgzWkyPxERiVYKNyEQ32IJBl/LTeeXYAjFeBsftdyIiEi0UrgJgaa5bpqFmwBabkIZbkamjQS8T0x1l9mdRUREgkHhJgRazFLs65YKYAkGX7dUVlzww81xycdhs9ioaKhgX9W+oF9fREQkUhRuQqDl+lLNuqU62ULia7nJTgj+mBu71c6QlCGAuqZERCS6KNyEQIslGBKOWIKhE0I5oBhgeKp33I0m8xMRkWiicBMCLZZgcMR7l2CATj0O3uBp4EC1dxByKMbcAIzs4x13s7Vka0iuLyIiEgkKNyFw1CzFAcx1c6D6ACYmdoudtJi0YJcIqOVGRESik8JNCBy9vlTnn5jydUllxmViMULzj8n3OHhhVSFldWUhuYeIiEi4KdyEQMKR4Sa+84tnhnIwsU+iI5F+Cf0Atd6IiEj0ULgJgRYDiiGgbqlQPgbenG+mYo27ERGRaKFwEwItBhQDJDQuetmJbqlQTuDXnC/cqOVGRESihcJNCPjH3NQf2S3V8SUYwtEtBWq5ERGR6KNwEwItll+AZt1SHZ+lOFzdUsNShwGQX5ZPg7shpPcSEREJB4WbEDhqQHEXuqVCNYGfT3Z8Ngn2BFymi13lu0J6LxERkXBQuAmBpkfBGwcUN++W6sASDNUN1ZTXlwOhH3NjGIZ/GYavD38d0nuJiIiEg8JNCPgGFB/VLdXBJRh8rTaJ9kQSHAkhqbG5oalDAYUbERGJDgo3IdB8Ej/TNL1LMPhCSgceB/dP4BefGbIam/OHm1KFGxER6fkUbkIgrnH5BZfHpM7l8e7sxER+4Rpv4zM0xRtudhzeEZb7iYiIhJLCTQjEO6z+94EsweB/UirE4218fC03+6r2UdnJlctFRES6G4WbELBZLcTYvb/aqgBmKfaFm3C13CQ7k8mI89a3o1StNyIi0rMp3IRIwpFz3cR3PNwUV3nnwwlXyw00td5sP7w9bPcUEREJBYWbEDlqluJOdEsVV3vDTWZceAYUAwxL8U7mpyemRESkpwso3LzzzjvBriPqxDvamqX42EswHK47DEBaTFpIamuNnpgSEZFoEVC4Ofvssxk8eDD33HMPu3fvDnZNUeGoWYrjO7YEg8vjoqyuDIDUmNSQ1Xek5nPdmB2YaFBERKS7Cijc7N27l7lz5/LSSy9x3HHHMXPmTF588UXq6+uDXV+PFX/UyuAd65YqrSsFwMAgxZkSouqOlpech9WwUl5fzv7qji8TISIi0t0EFG7S09O5/fbb2bx5M2vXrmXYsGH88Ic/JCcnhx//+Md8+umnwa6zx2laPLNzSzCU1JYA3ieYrBZrm8cFm9PqZEDSAEBdUyIi0rN1eUDxxIkTmT9/PnPnzqWyspKnnnqKSZMmccopp/DFF18Eo8Ye6ejFM5stwVBX0eZ5h2vDP97GxzeZnwYVi4hITxZwuGloaOCll17i3HPPZeDAgbz55pv84Q9/oLi4mB07djBw4EAuu+yyYNbao8QfGW6aL8FQ1fagYl+4Ced4Gx+tMSUiItHAFshJP/rRj3j++ecxTZNrr72WBx54gNGjR/u/j4+P58EHHyQnJydohfY08UfOcwPerqn6Su9cN30Gt3qer1sqIi03jeFGE/mJiEhPFlC4+fLLL3nssce45JJLcDqdrR6Tnp7eqx8ZTzhyQDFAQiYc3tXuoGLfY+CpzvC33PjmutlZuhOXx4XNEtD/PERERCIqoG6pRYsWcdlllx0VbFwuF++99x4ANpuNU089tesV9lBHDSgGSDj24pklNY0tN7Hhb7npl9iPWFss9Z56CioKwn5/ERGRYAgo3EyfPp2SkpKj9peVlTF9+vQuFxUNjhpQDB1agiGSLTcWw8KQlCGAxt2IiEjPFVC4MU0TwzCO2n/o0CHi4+O7XFQ08M1Q7F9+AbzdUtBut1Qkx9yABhWLiEjP16lBFZdccgkAhmFwww03tOiWcrvdfPbZZ0ybNi24FfZQrQ4o7kC3VCSflgI9Di4iIj1fp8JNcnIy4G25SUxMJDY21v+dw+HgO9/5DjfddFNwK+yhAu6WinS40RpTIiLSw3Uq3Pz5z38GYNCgQfzkJz9RF1Q7mpZfaD6guP1uKbfH7V9+IdLdUnsq9lDdUE2cPS4idYiIiAQq4KelFGza52+5qXc1LUSZ0P4SDKV1pZh494dzXanm0mLS6BPTBxOTnaU7I1KDiIhIV3S45WbixImsWrWK1NRUJkyY0OqAYp+NGzcGpbiezDfmxjShut7t/RzfbAmG2jKITWlxjq9LKtmZHNE5ZoamDuVQ4SG+Lv2aMX3HRKwOERGRQHT4L+iFF17oH0B80UUXhaqeqBHnsGIY3nBTVefyhhtHnDfgVO33TuYXO6HFOZF8DLy5oalD+bjwYw0qFhGRHqnD4WbRokWtvpfWGYZBvMNGZZ2LyjoXGb4v+gzxhptDOyGnZbiJ9GPgPsNSvTMVK9yIiEhPFNCYm927d7Nnzx7/508++YTbbruNZcuWdfpaS5YsYdCgQcTExDBlyhQ++eSTdo8vLS1lzpw5ZGdn43Q6GTZsGG+88Uan7xsOvkHF1fXNBhX71pQ6dPT6Td0l3PgGFW87vK1pvJCIiEgPEVC4ueqqq/zrRhUVFTFjxgw++eQTfvnLX3L33Xd3+DovvPAC8+bNY9GiRWzcuJFx48Yxc+ZM9u9v/Wmi+vp6zjzzTPLz83nppZfYtm0bTz75JP369Qvkxwi5Vue66eOdAbi1cBPpx8B9BicPxmJYKK0r5WDNwYjWIiIi0lkBhZvPP/+cE088EYAXX3yRMWPG8NFHH/Hss8/y9NNPd/g6Dz/8MDfddBOzZ89m1KhRLF26lLi4OJ566qlWj3/qqacoKSnhlVde4aSTTmLQoEGceuqpjBs3LpAfI+RaneumnXDja7mJdLiJscUwIHEAoK4pERHpeQIKNw0NDf7BxW+99RYXXHABACNGjKCwsLBD16ivr2fDhg3MmDGjqRiLhRkzZrBmzZpWz3n11VeZOnUqc+bMITMzk9GjR3PvvffidrtbPR6grq6O8vLyFlu4+JZgaL3lZudRj4P7Wm4i3S0FTeNuth/eHuFKREREOiegcHP88cezdOlS3n//fVauXMnZZ58NwL59++jTp0+HrnHw4EHcbjeZmZkt9mdmZlJUVNTqOd988w0vvfQSbrebN954gwULFvDQQw9xzz33tHmfxYsXk5yc7N9yc3M7+FN2Xby/5aZZ+ErLAwyoK4eqAy2O9z0tpXAjIiISuIDCzf33388f//hHTjvtNK688kp/t9Crr77q764KBY/HQ0ZGBsuWLWPSpEnMmjWLX/7ylyxdurTNc+bPn09ZWZl/2717d8jqO1JC44DiyrqGpp02J6R4u3yO7Joqqeke3VKgZRhERKTnCmimuNNOO42DBw9SXl5OamrTH+Kbb76ZuLiOTdefnp6O1WqluLi4xf7i4mKysrJaPSc7Oxu73Y7VavXvGzlyJEVFRdTX1+NwOI46x+l0tljgM5z6JHjve6Ci7ogvhkDpt95wM7BpodHuMs8NNLXc7CzdSYOnAbvFHuGKREREOiaglhsAq9XaItiAd82pjIyMNs5oyeFwMGnSJFatWuXf5/F4WLVqFVOnTm31nJNOOokdO3bg8Xj8+7Zv3052dnarwSbSspNjACgsq235RSuDij2mJ+LrSjWXk5BDnC2OBk8DBeUFkS5HRESkwwIKN8XFxVx77bXk5ORgs9mwWq0tto6aN28eTz75JM888wxbt27llltuoaqqitmzZwNw3XXXMX/+fP/xt9xyCyUlJdx6661s376d119/nXvvvZc5c+YE8mOEXHayd9X0tsNN09pNZXVleExvaEuJSQlHee2yGBZ/15TG3YiISE8SULfUDTfcQEFBAQsWLCA7O7vddabaM2vWLA4cOMDChQspKipi/PjxrFixwj/IuKCgAIulKX/l5uby5ptvcvvttzN27Fj69evHrbfeys9//vOA7h9q2Snelpuio8LN0RP5+Z6USnIkdZsuoGGpw/j0wKdsP7ydc/LOiXQ5IiIiHRJQuPnggw94//33GT9+fJcLmDt3LnPnzm31u9WrVx+1b+rUqXz88cddvm845DS23BSV1+L2mFgtjSHQ13JT8g143GCxcqj2ENA9uqR8/IOKNdeNiIj0IAF1S+Xm5mpa/g7om+jEajFwe8yWg4qTc8HqBHc9lHmf3uousxM3p8fBRUSkJwoo3DzyyCPccccd5OfnB7mc6GK1GGQmep+YKiyrafrCYjmqa8ofbrrBk1I+Q1K8LUyFVYWU14dv8kMREZGuCCjczJo1i9WrVzN48GASExNJS0trsUmT7JS2BhX7wo13UHFJXeOimbHd5/eX7EwmK977WP6Ow0cvFyEiItIdBTTm5pFHHglyGdErq4OPg3fHlhvwdk0VVRWx/fB2JmZOjHQ5IiIixxRQuLn++uuDXUfUyvGFm9Kall8cEW58i2Z2pwHFAENThvLenvc07kZERHqMgCfx27lzJ3feeSdXXnkl+/fvB+Df//43X3zxRdCKiwZZx5rr5uARLTfdaEAxNA0q1hNTIiLSUwQUbt59913GjBnD2rVrefnll6msrATg008/ZdGiRUEtsKfzt9yUtdFyU7YbGmr8LTfdLdw0X2NKT8iJiEhPEFC4ueOOO7jnnntYuXJli2UPTj/99B4zB024tDmgOK4PxCQDJpTs8rfc9Inp2Krq4TIoeRA2i42qhir2Ve2LdDkiIiLHFFC42bJlCxdffPFR+zMyMjh48GCXi4omvvWl9lfU4XI3rYmFYfhbbzwHt/vXlepuLTd2i53Byd4nu7aXaNyNiIh0fwGFm5SUFAoLC4/av2nTJvr169floqJJeoITm28iv8pWVgcHyg98idt0A93vaSlAa0yJiEiPElC4ueKKK/j5z39OUVERhmHg8Xj48MMP+clPfsJ1110X7Bp7NKvFIDPJ23qzr7T1QcUlh7yDdRPtidit3WNdqeb8g4pLNahYRES6v4DCzb333suIESPIzc2lsrKSUaNGccoppzBt2jTuvPPOYNfY4/m6ptpaQPNw6S6ge03g15xabkREpCcJaJ4bh8PBk08+ycKFC9myZQuVlZVMmDCBoUOHBru+qJCdEgvfHm7zianDlfsgxdktu6SgqeXm2/JvqXXVEmOLiXBFIiIibetwuJk3b1673zd/Surhhx8OvKIolN3WLMVp3pabElcV4Ox2g4l9+sb2JdWZyuG6w+wo3cHo9NGRLklERKRNHQ43mzZtavF548aNuFwuhg8fDsD27duxWq1MmjQpuBVGgey25rpxJkBiNiXWKqD7zU7sYxgGI9JGsKZwDVtLtirciIhIt9bhcPPOO+/43z/88MMkJibyzDPPkJrqbW04fPgws2fP5pRTTgl+lT1cduMsxUcNKAboM4TD5Z8D3e8x8OZG9PGGm68OfRXpUkRERNoV0IDihx56iMWLF/uDDUBqair33HMPDz30UNCKixZtDigG6DOYw1bvP4bu2nIDMDJtJABflSjciIhI9xZQuCkvL+fAgQNH7T9w4AAVFRVdLiraZKf4JvKrbTmRH0CfIZRYrUA3b7lJGwF4n5hye9wRrkZERKRtAYWbiy++mNmzZ/Pyyy+zZ88e9uzZw9///nduvPFGLrnkkmDX2OOlxzuxWw08pnem4hayxlLia7lxdt+Wm4FJA4m1xVLrriW/PD/S5YiIiLQpoHCzdOlSzjnnHK666ioGDhzIwIEDueqqqzj77LN5/PHHg11jj2dpNpHfUYOKcyZw2NLYctON16W0GBaGp3oHj28t2RrhakRERNoWULiJi4vj8ccf59ChQ2zatIlNmzZRUlLC448/Tnx8fLBrjAo5bQwqNp2JlDZ2S6UdLgh7XZ0xPM0bbjSoWEREurOAJvHziY+PZ+zYscGqJapltTGouLy+HJfhfZ+6f1u4y+oUDSoWEZGeIKCWG+k836DifUd0Sx2uPQxAgseDY++mo87rTkb08Q4q3lqyFdPsxn1oIiLSqynchImvW+rIlpuS2hIAUt1u2LcRPJ6jzu0uhqYMxWbYKK8vp7Dq6FXhRUREugOFmzDxdUvtOyLc+FpuUj1AbRmU7Ax3aR3msDo4LuU4QIOKRUSk+1K4CRNfy01hactuqZI6b8tNmj3Bu2PvhrDW1Vm++W407kZERLorhZsw8Y25OVBZR0OzifxKa0sBSI3L8O7Ysz7cpXWKf1CxnpgSEZFuSuEmTNLiHDisFkwTisubuqZK60oBSEnK9e7oIS036pYSEZHuSuEmTCwWo9XHwf3hps8w746iLdDQyhpU3YQv3BRXF/vHC4mIiHQnCjdhlN3KoGJ/uEkZBHHp4GmA4s8jUF3HJDgSyE30tjKp9UZERLojhZsw8oWb5oOKfeEmOSYF+k3y7uzm4240qFhERLozhZswyk5pfGKqecuNb0CxMxX6T/bu7ObjbjSoWEREujOFmzDyt9yUHd1yk+JMgX4TvTv39oyWG3VLiYhId6RwE0bZR8xS7PK4qKivACAlJgVyGsNNyTdQXRKJEjtkZB9vy8235d9S3VAd4WpERERaUrgJoyMHFJfXl2PiXaMpyZEEcWmQNth78N6NEamxI9Jj00mPTcfEZPvh7ZEuR0REpAWFmzDyhZuDlXXUuzz+LqkkRxI2S+MC7f5xN+qaEhERCYTCTRilxTtIdNowTdh1sMo/mDjFmdJ0kO+JqZ4yqFhPTImISDejcBNGhmEwMjsJgC8Ly1oOJvbp19hys2c9mGZ4C+yEUX1GAfD5we47J4+IiPROCjdhNiqnMdzsK28KNzEpTQdkjQarA2pK4PCu8BfYQWPSxwCwo3SHBhWLiEi3onATZqP8LTflrbfc2JyQPc77Pv/D8BbXCZnxmWTEZeAxPXx56MtIlyMiIuKncBNmvpabrYUVrY+5ARh8hvd1x8rwFRYAX+vNloNbIlyJiIhIE4WbMBuSkYDVYlBSVU9hxSGglXAzZIb3dedqcLvCWl9nKNyIiEh3pHATZjF2K0P6JgCwt+IgAMnO5JYH9ZsIsWlQVwZ71oW7xA4b23csoHAjIiLdi8JNBPi6pg5UeWchTo1JbXmAxQqDT/e+78ZdU6P6jMLAoKiqiAPVByJdjoiICKBwExEjsxMBKKsvA1rplgIYeqb39evuG27i7fEMTvHOqKzWGxER6S4UbiJgVLa3G6rW07iuVGvhxtdyU/QZVBSHqbLOU9eUiIh0N90i3CxZsoRBgwYRExPDlClT+OSTTzp03vLlyzEMg4suuii0BQaZt+XGg8eoAtoINwkZkD3e+37nqnCV1mmj00cDsOWAwo2IiHQPEQ83L7zwAvPmzWPRokVs3LiRcePGMXPmTPbv39/uefn5+fzkJz/hlFNOCVOlwdMnwUlmMhiGdwbiVsMN9IiuqbHp3pabzw99jsf0RLgaERGRbhBuHn74YW666SZmz57NqFGjWLp0KXFxcTz11FNtnuN2u7n66qu56667OO6449q9fl1dHeXl5S227uC4LO+rwxKL3Wpv/SD/I+Fvd9tHwgenDCbWFktVQxW7yrrvjMoiItJ7RDTc1NfXs2HDBmbMmOHfZ7FYmDFjBmvWrGnzvLvvvpuMjAxuvPHGY95j8eLFJCcn+7fc3Nyg1N5VueneVhurmdD2Qf0mQ0wy1JbCvo3hKayTbBabf52pzw58FuFqREREIhxuDh48iNvtJjMzs8X+zMxMioqKWj3ngw8+4E9/+hNPPvlkh+4xf/58ysrK/Nvu3bu7XHcwZKS4AfC44to+yGprGljcjbumNJmfiIh0JxHvluqMiooKrr32Wp588knS09M7dI7T6SQpKanF1h2kJDYAUFPrxOVuZ6zKkMZxN914vhtfuNEK4SIi0h3YInnz9PR0rFYrxcUtH3UuLi4mKyvrqON37txJfn4+559/vn+fx+MNBjabjW3btjF48ODQFh0kNlsNAG5XHPmHqhiSkdj6gUMa15natwkqD0BC3zBV2HG+x8G3H95OjauGWFtshCsSEZHeLKItNw6Hg0mTJrFqVdOjzh6Ph1WrVjF16tSjjh8xYgRbtmxh8+bN/u2CCy5g+vTpbN68uduMp+mIsvpSAEx3HF8WVrR9YGIWZHlbRtj5dugLC0BmXCbpsem4TTdbD22NdDkiItLLRbxbat68eTz55JM888wzbN26lVtuuYWqqipmz54NwHXXXcf8+fMBiImJYfTo0S22lJQUEhMTGT16NA6HI5I/SqeU1pUCjeFm3zGe4OrmXVOGYWjcjYiIdBsR7ZYCmDVrFgcOHGDhwoUUFRUxfvx4VqxY4R9kXFBQgMUS8QwWdKW1pQCY7ni+LDxGuBl6JnzwMHz9H3DVgc0Z+gI7aWzfsbyz+x2FGxERibiIhxuAuXPnMnfu3Fa/W716dbvnPv3008EvKAz8LTeuDrTc5E6BxByo2Afb34RRF4S+wE7STMUiItJdRF+TSA/hCzd44jlYWcf+itq2D7ZYYcyl3vefvRDy2gIxus9oDAz2Ve3jYM3BSJcjIiK9mMJNhPjCTU5iHwC2tjeoGGDcFd7X7W9CdUkIKwtMgiOB45K9s0V/euDTCFcjIiK9mcJNBJim6Q83w/p6xxYds2sq83jIHA2eBvjyldAWGKAJmRMA2FC8IcKViIhIb6ZwEwFVDVW4PN61osZk5wCwZW/psU8cO8v7+mn37JqanDkZgPVF6yNciYiI9GYKNxHga7WJscZw0nHZAHy08xBuj9n+iWMuBQzY/TGUdL9FKn3hZtvhbVTUH6ObTUREJEQUbiLAF25SYlIYn5tCotNGaXUDn+8ta//EpBw47lTv+y1/C22RAciMzyQ3MReP6WHT/k2RLkdERHophZsI8IcbZwo2q4VpQ7yDit/bfuDYJ/u6pj57AcxjtPREgL9rqlhdUyIiEhkKNxHgCzfJzmQAvjvMu17U+1934BHqkeeDLRYO7YC9G0NVYsAmZ3nDzYYiDSoWEZHIULiJAN/sxKnOVAC+O9QbbjYWHKaitqH9k52JMOI87/tuOOeNr+Xmi0NfUN1QHeFqRESkN1K4iYAjW25y0+LIS4/H5TFZs/PQsS/g65r6/O/gPkYYCrOchBxy4nNwm242798c6XJERKQXUriJAF+4SY1J9e87ZWg6AO993YFxN4NPh7h0qD4IO1Yd+/gwm5Q5CdC4GxERiQyFmwhoPqDY55ShnRh3Y7XB2Mu979c9GeTqus437kbhRkREIkHhJgJaCzdTB/fBZjH49lA13x6qOvZFTrwJMGDHW7D/q5DUGSjfuJstB7dQ46qJcDUiItLbKNxEgG9AcfNwk+C0MXGgt5vqvY603qQd1zSw+OPHg1xh1+Qm5pIRm4HL4+KzA59FuhwREellFG4ioLWWG4BTfY+Ed2S+G4Cpc7yvny6Hqu6zErdhGEzK0rgbERGJDIWbCGg+Q3FzvkHFH+08RIPbc+wLDZgKORPAXQfr/hTkKrvG1zWlRTRFRCTcFG7CrMZVQ527Dji65eb4nGRS4+xU1rnYvLv02BczDJg61/t+3ZPQUBvcYrvAN6j4swOfUe+uj3A1IiLSmyjchJlvvI3dYifOFtfiO6vF4OTGp6Y6tBQDwKgLIakfVB2Az18KZqldkpeUR1pMGnXuOrYc3BLpckREpBdRuAmz5uNtDMM46vum+W46OIbGaocTb/a+X7Ok26w3ZRhG03w3RRp3IyIi4aNwE2ZHzk58JN9SDJ/tKeVwVQe7cyZdD/Z42P8lfPNOMMoMCi2iKSIikaBwE2atzU7cXFZyDMMzEzFNWPllcccuGpsKE67xvl+zJAhVBscJWScAsGn/Jmpd3Wc8kIiIRDeFmzBr6zHw5i6ckAPASxv2dPzC3/kf/JP67e0eTygNSRlCZlwmde461hWti3Q5IiLSSyjchFlrE/gd6ZIJ/bEY8El+CfkHOzBbMXgn9fMtybDq7q4VGSSGYXBK/1MAeH/v+xGuRkREeguFmzDrSMtNVnKMf62pv2/sROvN9F+AxQ7frIad3WPszSn9GsPNnvcxu8lgZxERiW4KN2F2rAHFPpdN7g/A3zfswe3pYChIHQQn3Oh9/9avwNOBiQBD7DvZ38FmsbGncg/55fmRLkdERHoBhZswO9aAYp8ZIzNJirGxr6yWj3Z2YmmFU34CjgQo3Axb/xl4oUESZ4/zPzX1wd4PIlyNiIj0Bgo3YdaRbimAGLuVC8f3Azo5sDihb9Osxat+De6GAKoMruZdUyIiIqGmcBNmHRlQ7HPpJG/X1IrPiyir6URImTYX4tKhZCds+msAVQaXb1Dx+uL1VDdUR7gaERGJdgo3YdbRlhuAsf2TGZaZQJ3Lw2uf7ev4TZyJ8N2fet+vvh/qIxsoBiUNon9Cfxo8DawtXBvRWkREJPop3IRRnbuOapc3aBy5InhrDMPgskm5QCe7pgAmz4aUAVBZBB9HdmI/PRIuIiLhpHATRodqDgHgsDhItCd26JwLJ+RgtRhsKihlx/6Kjt/M5oTTF3jfv/cQHM7vZLXBdXK/kwHvoGI9Ei4iIqGkcBNGB2u8Tz2lx6a3umhmazISY5g+3Dvnzd8623oz5jIYdAq4auC1eRFdVPOErBNwWp0UVhWys3RnxOoQEZHop3ATRs3DTWdc6uuaWr+Hmnp3x080DPjeI2B1ws5V8PnfO3XfYIq1xfrXmlLXlIiIhJLCTRj5wk2f2D6dOu+MkRkMSIvjUFU9z679tnM3TR8C3/2J9/2KO6C6pHPnB5H/kXCFGxERCSGFmzDyjbnpbMuN3WphzvTBACx99xtqGzrRegNw0m2QPhyqDsBbizp3bhD5BhVvKt5ERX0nxg+JiIh0gsJNGAXaLQVw8YT+9EuJ5WBlHc+tLejcyTYHnP+o9/3Gv0D+h52+fzDkJuYyKGkQLtPFmn1rIlKDiIhEP4WbMOpKuHHYLMyZPgSApe/u7HzrzcCpMOkG7/t/3Qquuk7XEAzf7f9dAFYVrIrI/UVEJPop3ITRwdrAxtz4XDqpPznJMeyvqOOFdbs7f4EZv4L4DDj0NayMTPfU2YPOBuCd3e9Q46qJSA0iIhLdFG7CKNAxNz4Om4VbGltvnli9kzpXJ1tvYlPhwj943699Arb9O6A6umJ0+mj6JfSjxlXDu3veDfv9RUQk+inchIlpml3qlvK5fHJ/spJiKCqv5cX1nZz3BmDYTPjOHO/7V34IZXsDriUQhmFwTt45AKzYtSKs9xYRkd5B4SZMKhsqqXN7x7n0iQmsWwrAabNyy2neJ6eeeGdH51tvAGYsguzxUFMCL98EngCu0QW+cPP+nvf11JSIiASdwk2Y+FptEu2JxNhiunStWSfkkpHoZF9ZLc939skp8C7NcOlT4EiAbz+E937bpXo6a2jKUAYnD6beU8/bBW+H9d4iIhL9FG7CJNAJ/FoTY7fyozOGAvDQf7ZTXF7b+Yv0GeydvRjg3fsh/4Mu19VRhmFwdp53YPG/88M/7kdERKKbwk2YdHUw8ZGuOnEA43JTqKhzcde/vgjsImMvg/FXg+mBl74PZQGM4QmQ76mpj/d9TElt5GZNFhGR6KNwEybBGEzcnNVisPjiMVgtBm9sKWLV1uLALnTOA5BxPFQWw/NXQF1lUOo7lkHJgxiZNhK36eatb98Kyz1FRKR36BbhZsmSJQwaNIiYmBimTJnCJ5980uaxTz75JKeccgqpqamkpqYyY8aMdo/vLoIdbgBG5STxg5PzAFj4zy+oqnN1/iLOBLhqOcT3haItYR1gfG7euQD8e5e6pkREJHgiHm5eeOEF5s2bx6JFi9i4cSPjxo1j5syZ7N+/v9XjV69ezZVXXsk777zDmjVryM3N5ayzzmLv3vA+0txZwRxz09ytM4bSLyWWvaU1PPLW9sAukjIArnjeu3r4tjfCtv7UzEEzAdhQvIHiqgBbnkRERI4Q8XDz8MMPc9NNNzF79mxGjRrF0qVLiYuL46mnnmr1+GeffZYf/vCHjB8/nhEjRvC///u/eDweVq3q3tP5+2YnDmbLDUCcw8Y9F40G4KkP8/l8b1lgF8o9AS563Pv+o8dgwzNBqrBt2QnZTMiYgInJm/lvhvx+IiLSO0Q03NTX17NhwwZmzJjh32exWJgxYwZr1nRsYcXq6moaGhpIS0tr9fu6ujrKy8tbbJEQ7AHFzU0fkcF5Y7Nxe0x+8Y8tNLg9gV1ozKVw6h3e96/Pg52hf0zbN7B4Rb4m9BMRkeCIaLg5ePAgbrebzMzMFvszMzMpKirq0DV+/vOfk5OT0yIgNbd48WKSk5P9W25ubpfrDkQoxtw0t+h7o0iMsfHZnjLu//dXgV/otDtg9H+BxwXLr4ZvPwpeka04a9BZWAwLWw5u4ZvSb0J6LxER6R0i3i3VFffddx/Lly/nH//4BzExrU+MN3/+fMrKyvzb7t0BLDjZRW6P2/+4c6jCTUZSDA9eNg6A//1gF29sKQzsQoYBFz0BQ2ZAQzU8ezns2RDESltKj03n1P6nArB82/KQ3UdERHqPiIab9PR0rFYrxcUtB5MWFxeTlZXV7rkPPvgg9913H//5z38YO3Zsm8c5nU6SkpJabOF2uO4wHtODxbCQ6kwN2X1mHp/Ff596HAA//dun7DwQ4GPdNifM+j8YdArUV8D/XQyFnwWx0pauGHEFAK/ufJWqhqqQ3UdERHqHiIYbh8PBpEmTWgwG9g0Onjp1apvnPfDAA/z6179mxYoVTJ48ORyldolvvE2qMxWrxRrSe/30rOFMyUujqt7N//x1Q2CPhwPYY+HK5dD/RKgtg79eBPu70N3Vju9kf4dBSYOoaqjitZ2vheQeIiLSe0S8W2revHk8+eSTPPPMM2zdupVbbrmFqqoqZs+eDcB1113H/Pnz/cfff//9LFiwgKeeeopBgwZRVFREUVERlZXhmXwuEKEeb9OczWrhsasmkJHo5Ov9lcx/eQumaQZ2MWcCXPOSd5HN6kPwlwug6POg1gtgMSzMGj4L8HZNBVyviIgI3SDczJo1iwcffJCFCxcyfvx4Nm/ezIoVK/yDjAsKCigsbBo/8sQTT1BfX8+ll15Kdna2f3vwwQcj9SMcUzjDDUBGYgxLrp6I1WLw6qf7+NMHuwK/WEwyXPsPyBztncX46XOh4OPgFdvogiEXEGuLZUfpDtYXrw/69UVEpPeIeLgBmDt3Lt9++y11dXWsXbuWKVOm+L9bvXo1Tz/9tP9zfn4+pmketf3qV78Kf+EdFKoJ/NpzwqA05p8zAoB7Xt/Kyxu7sG5UXBrc8BrkTvF2Uf3lItj+n+AU2ijJkcR5x50HwPKvNLBYREQC1y3CTbQLd8uNz40n5zH7pEEA/PSlz/jPFx17vL5Vsalw7Ssw5Exw1cDyK+GzF4NSp88Vw70Di98ueJv91a3PUC0iInIsCjdhEMoJ/NpjGAYLzhvFpZP64/aYzH1uEx/uOBj4BR1xcOXzMOZy7zw4L98EH/4egjRGZnjacCZmTMRlunhp+0tBuaaIiPQ+CjdhEKqlFzrCYjG475IxzDw+k3q3h5v+sp5NBYcDv6DVDhf/Eabc4v28cgG8cgs01AalXt9j4S9tf4kGT0NQrikiIr2Lwk0YRKpbysdmtfD7KydwytB0quvd3PDndV0LOBYLnL0YznkADCt8+jw8fR6UBzhxYDMzBsygT0wfDtQcYFVB914vTEREuieFmzCIxIDiIzltVv547SQmDUylrKaBq55cy9tfdWElbsOAKf8N1/wdYlJg73p4cjrs7dpsxnarnUuHXQrAM58/o8fCRUSk0xRuQqzOXUdFfQUQuZYbnziHjb98/0S+O6wvNQ1ubvrLBl5c18XlKAZPh5vehr4joKIQnjoH1v1vl8bhXDniSmJtsXx+6HPe3fNu1+oTEZFeR+EmxHyDiR0WB4n2xAhXA/FOG3+6fjKXTOyH22Pys79/xmOrvu5aC0mfwXDjShh+Lrjr4PX/By9eCzWBdX31ie3DlSOuBGDJ5iV4zABXORcRkV5J4SbEmo+3MQwjwtV42a0WHrpsHD88bTAAD63czvyXt1Db4A78ojFJMOtZmHkvWOyw9V+w9BQoWBvQ5WYfP5t4ezxflXylsTciItIpCjchFunBxG0xDIOfnT2Cuy44HsOA5et2c9nSNRQcqg78ohYLTJ0DN/4HUvOgbDf8+Rx49wFwd+7Jp5SYFK4ZeQ0Aj29+HLenC8FLRER6FYWbEOsOg4nbc/20QTw9+0RS4+xs2VvGeY+937XJ/gD6TYT/fg/GXAamG975DTx5OhRt6dRlrjv+OhIdiewo3cGb+W92rSYREek1FG5CLFIT+HXGqcP68vqPT2HigBQqal3c/NcN3PvGVhrcXRjrEpMElzzp3WJToegzWHYavHMvuOo7dIkkRxLXj7oegCc+fQKXJ8AVzkVEpFdRuAmx7totdaSclFiW3zyVG0/OA2DZe99w4R8+5PO9ZYFf1DBg7OXww7Uw8nzvrMbv3u8NObs/6dAlrhl1DSnOFPLL83n9m9cDr0VERHoNhZsQ6ynhBsBhs7Dge6NYes1EUuLsfFlYzoVLPuS3b37VtcHGiZlw+V/h0j9DXB/Y/wX86Ux45YdQ2f4aUvH2eGaPng3A0k+XatZiERE5JoWbEPMtvdBdx9y05uzR2ay8/VTOG5ON22Oy5J2dfO+xD9jwbRdmNTYMGH0JzPkEJngHCrP5WXhsEnz8BLjb7nK6YvgVpMWksadyD89++WzgNYiISK+gcBNiPWHMTWv6JjpZcvVEnrh6IukJTnbsr+S/nviIeS9sprCsJvALx6fDhUvgxrcgexzUlcOKO2DpybBtRauT/8XZ47h14q0APP7p4+yp2BP4/UVEJOop3ISQaZocqD4A9Lxw43POmGzemvddLp3UH4CXN+1l+oOreeSt7VTXd2GAb+4JcNM78L3feQccH9gKz8/yrlG1e91Rh1885GImZ06mxlXDPWvv0bIMIiLSJoWbEKpoqKDe430yqE9Mz+mWOlJKnIMHLxvHP+ecxOSBqdQ2eHjkra85/cF3eWFdQeBPVVmsMPn78ONNcNKtYHXCtx/Cn2bAC9fA/q/8hxqGwcKpC7Fb7Hy490NW5K8I0k8nIiLRRuEmhHyDiRPticTYYiJcTdeNy03hb/8zlT9cNYF+KbEUldfy879vYfqDq1n+SQH1rgBDTmwqnHk3/HijdzyOYfHOcPz4d+DF66HocwDykvO4eezNANz3yX2U1XXhSS4REYlaCjch5Btv05MGEx+LYRh8b2wOq/7fqdx53kjSE5zsOVzDHS97Q85zawsCf7Iqub93PM4tH3kfHceEL1+BpSfB81fBvk3cOPpGBicPpqS2hIc3PBzMH01ERKKEwk0I9aTHwDsrxm7lB6ccx/s/m+4POXtLa/jFP7Zw0n1v8/B/trG/ojawi2eMhFn/5w05x18CGLDtdVh2Gva/XszC/mcD8PLXL7Ou6OjxOSIi0rsp3IRQNIcbn1hHU8hZ8L1R9EuJ5VBVPb9/ewcn3fc2817czGd7SgMbAJx5PFz2Z5izFsbOAsMK+e8z8bWfc1m9FYCFHy6gvL48yD+ViIj0ZAo3IdQbwo1PrMPKjSfn8e5PT2PJVROZOCCFBrfJyxv3csEfPuTc33/AMx/lU1YdwCR8fYfDJcvg1k9h2o/AmcRthfnkNLjYU7mXX7x8CZ6DXwf/hxIRkR5J4SaEiqq8C1BG05ibY7FZLZw3NpuXf3gS//jhNC4cn4PDZmFrYTmLXv2CE+59i1uXb2L1tv2df8oqJRfOugfmfUnSmffyuxobTo+Hd+uK+eP/nQ7PXABf/rPDa1eJiEh0MsxeNmFIeXk5ycnJlJWVkZSUFNJ7XfTKRews28mSM5bw3f7fDem9urPS6npe2bSX5et281VRhX9/WryDc8dkccG4fkwemIrFYnTuwh43r3y0mAU7X8AwTZYUH+CUmlrvEg9jLofxV0H22CD/NCIiEgmd+futcBMi1Q3VTH1+Kh7TwzuXv9MruqaOxTRNPttTxt837uH1zwo5VNXUwpKVFMOZozI56/hMpuT1wWHreKPir9f8mhe3v0iiYeeFg5Xklhc3fZk5GsZd6V36ISknmD+OiIiEkcJNO8IVbjbv38y1/76WvrF9efvyt0N2n57K5fbw4c5DvLp5H//5ooiKuqbZjhNjbJw+IoPTR2RwytC+pMU72r1Wvbue2Stm89nBzxieOoy/Dvs+sVv+Bl+9Dm5fgDJg4DRvyBl5IST0DeFPJyIiwaZw045whZvntj7H4k8W893+32XJGUtCdp9oUNvgZs3OQ/znyyJWflnMwcqmFh3DgLH9Uzh1WF9OHZbO2P4p2K1Ht+oUVRUx67VZlNSWcHK/k3l0+qM46irhi5dhy0tQsKbpYMMKg06CEefDiPMguV84fkwREekChZt2hCvcLPhwAa/seIX/HvvfzJ0wN2T3iTZuj8nm3Yf5z5fFvLvtQIsxOgDxDisn5KUx9bg+TB3ch+NzkrE2jtXZtH8TN//nZmrdtZza/1R+d9rvsFvt3hPL9sAX//AGncLNLW+aM9Ebcoad7X383Ojk2B8REQk5hZt2hCvcXPrqpWw7vI1Hpz/K6QNOD9l9ol1RWS3vfX2Ad7cd4MOdByk94lHyBKeNCQNSmDwwjRMGpeJybGfeez+mzl3H9NzpPHTaQ9gt9pYXLfnG22W19TXYvRZo9v8CSf1g6JkwdCbkfRecCaH/IUVE5JgUbtoRjnBT567jO89+B5fpYuWlK8mKzwrJfXobj8dka1E5a3Ye4uNvDrH2m5IWY3UArBaDAf32UJKwFA8NTMk8jT/MeIgYWxvjdiqKYdsbsO3fsOs9cNU0fWexQ+6JcNx0GDwdciZ4F/sUEZGwU7hpRzjCzecHP+fK168k1ZnKu7PexVA3R0i4PSZfFZWz4dvDrM8/zPr8EvaVeZd8sMZvI7b/XzAsbjyVYxlq3MSYnHRG90vi+JxkhmYm4LQdEVQaaiD/A9j+Jnz9JpQWtPw+JhkGngSDTvZumaMVdkREwkThph3hCDcvbnuRX3/8a6blTOOPZ/4xJPeQ1hWW1fDp7lI27y7j/b3v8q3tCQzDjbsml5o912C6kgFvC8/gvvEMz0piRFYiI7ISGZaZSL+UWO98O6bp7b765h3Y+Q7seh+OXIU8JhkGTIXcKTDgO96WHXtsBH5qEZHop3DTjnCEm7vW3MVL21/ixtE3ctuk20JyD+mYj/au4Sfv/oSKhnJiLCnk1P03u/b0pbzW1erxMXYLQzISGJqRyOC+8RzXN4Hj+sYzKNVJzIHP4dsPvK07366B+paDnbHYIWc89D8B+k2C/pMhZaAGKIuIBIHCTTvCEW6ueO0Kvjj0BQ+d+hBnDTorJPeQjttdsZtb37mVrw9/jc1i45dTfslJmefxVWEFW4vK2VZUwbaiCr45UEV9G0tCGAbkJMeSlx7PwD5x5KU6GW3N57jqLaSVbMS25xOo2n/0iXHp3qCTMx6yx0P2OO9kggo8IiKdonDTjlCHmwZ3A1Oem0KDp4E3LnmD3MTcoN9DOq+6oZo7P7yTld+uBOCCwRfwsxN+RrIz2X+My+2hoKSar/dXsmN/JTsPVPLNgSq+OVDZZkuPT3q8nUnJ5Zzs2MkoczsDa78irWIbFk8rC4XG94WssZA1GjLHeF/7DAWrLag/s4hINFG4aUeow81XJV9x2b8uI9GeyIdXfqjBxN2IaZos+2wZSzYvwcQkLSaN+SfOZ+agme3+czJNk0NV9ew6WEX+wSq+PVRN/qEq8g9VsbukhrKa1lc6d1LPKONbxlq+YaI9n7GWbxngKcDK0a1DptUJfYdhZIyCjJHQdyRkjIDkAWDR+rYiIgo37Qh1uPnH1/9g4UcLOTHrRP40809Bv7503eb9m1n00SK+KfsGgFP7n8qd37kz4Ef2y2oa2F1S7d0OV7P3cA17S2vY0/ha0azVJ4Y6RhoFjLQUMMrIZ6SlgBFGAfFGXavXbrA4qUzIoz51CKQPx5k1nPicEdj7DgFHfED1ioj0RAo37Qh1uPnNx79h+bbl3HD8Dfy/yf8v6NeX4Kh31/OnLX9i2ZZluDwuYm2xXDPyGq4//voWXVXBUF7bQGFpLfvKaigsraWwrIaislqKymspLq+luKya5LpChhl7GGbsZrjF+3qcUYjTaLs77IDRh2J7P0qd/amMz6UhaRBmWh7WPseRlJJGapyD1HgHqXF2Yu1WtSKKSI+mcNOOUIeba964hk8PfMr9p9zPucedG/TrS3DtLN3Jrz76FZsPbAYg3h7PtaOu5dpR15LkCN0M1keqqXezv6KW/RV1FJfXsr+8joMVVXgO5RNTuoOkql1k1OaT7d7HIKOQPkZFu9c7bCaw2+zbuGVQbPSlzJlFTWw2NXH9iElIJTnWTkqcneQ4O8mxTVtSTNP7xBgbtlbW8hIRCTeFm3aEMty4PC6mPjeVWnctr170KnnJeUG9voSGaZq8XfA2j3/6ONsPbwcg0Z7IFSOu4L+G/Rf9ErrPwpoej0lpTQOHDxZRVbgd94GvMQ7n46z4lviqAlLr9pLoLj3mdcrNWArNPhSafdhnplFo9qGINIpM71ZsplJOHGAQ57CSFGMnKdZGYow38DS92kh02khw2kho3Jfg/+z9Lt5pI86hliMR6RqFm3aEMtzsOLyDi1+9mDhbHGuuWoPF0H/x9iQe08Nb377F45sfZ2fZTgAMDKb1m8ZlQy/ju7nfPXqdqu6orhJKv4XD32IezqfhUD6ukm8xyvZgq9yLve5why5TYzooNlPZTwr7zRQOmCnsN1M5QDIHzeTGz8kcJhEX7T/pZRgQ77AR77QS77Q1vXfYiHPaiHdYiXN4Q1Bc4/5Yh9X72WEl1m5ret94bKzdSozdotAk0kt05u+3nj0Noq0lWwEYkTZCwaYHshgWzhp0FmcMOIO3d7/Ni9te5OPCj/lw74d8uPdD0mPTOXPgmZw+4HQmZU7qvkHHmeBd3TzzeAzA0bj51VVC+V7vSunle6FsL5TvgfJCqCiE8n1QW0qsUc8go5hBFB/zlpWWJMotyRw2kjlkJnHQk8h+TwJFDQkcMhMpIZHS+kQO1yWwj0RqcQblRzUMiLF5Q0+M3Rt8Yu3eLcZhJdZuIcZuJcbm/c5ptxBj8x4b4/uucZ+z2avT5t3vPOKzw6owJdITKNwE0ZeHvgRgVJ9REa5EusJqsXLmwDM5c+CZ7C7fzd+//jv/2PEPDtYc5Pmvnuf5r54n0ZHIqf1P5dTcU5mcOZn02PRIl91xzgToO9y7taW+2ht0KouhoqjZ637v+6r93vdVB8D0kOApJ8FTTg67W16njfzntjqptydTZ0+mxppEtTWJSksilUY8ZSRQbsZz2IyjxBNHiSuWg+4YDjbEsr8hhrIGC/Uu7+P0pgk1DW5qGtxB+uUcm9NmwWmz4LBZve8bQ4/Tbm36zmrBYfNuzsZXh9XaYp/dajQe591vtxqN+33HW7D7Xq3e74/8zm41sFss3iVDRMRP3VJBdMOKG9hQvIHfnPwbLhh8QVCvLZHV4G5gTeEa3i54m3d2v0NJbUmL7/OS85icOZnJmZMZ03cM/RP6947/wve4oabUG3aqDngDT/UhqDrofa0+CFWHoKak8XMJtDaxYWdYnZgxSZiOJNyORNyORBps8dTbEqm3xlNriaPOEkeNEUu1EUcVMVQTSyUxlHtiqPA4KXPHUOGxU+MyqHO5qW1wU9vgaXzf8rXO5aG7/1vSajH8QccXemwWiz802Xz7Ld6AZGsMSjaL4d9vszYd2/x7m9X7vbXx+jard5/NYnjPb7ye97PFf5yvJqul6Zo2i9G439K4v+X3Vv9no3f8/490isbctCNU4cZjepj2/DSqGqp4+YKXGZo6NGjXlu7F7XHz6YFPWVWwio8LP/YPQm4u0ZHIiLQRjEwbyYi0EeQl5zEoaRAJjoQIVNyNmCbUV3qDTk2pN/TUHPaGntpS7z7fa02pd7HSmjKoLTt64dJgsMWCI847Z5A93vvqiGt6b4/FtMfhscXissbSYImhwRpDg+Gk3nBS5391UGs6qMVJDXZqTQfVpp0aj506N9S7PNS7Pd5X3+Y+el+D27vVNX7vcpv+/b5jXR4Ttyf6/7VtMfCGJV8IsrYMQ74g1DwkNf9sMQxsVu9ra/tsFgOLxcBqtDzX9521ne+tFvzXbdrX8ljDwL+v+XWO2m8xsBhN17MYza5nAcNoPNYwsFhocYyl8VoWXw1Gy+sbBlEVEnvcmJslS5bw29/+lqKiIsaNG8djjz3GiSee2Obxf/vb31iwYAH5+fkMHTqU+++/n3PPjexj1wXlBVQ1VBFjjdFTUlHOarEyMXMiEzMnAlBWV8aG4g2sL17PxuKNbD+8nYr6CtYVrWNd0boW56bHpjMoaRC5iblkx2eTnZBNTnwO2fHZ9IntQ5w9LhI/UvgYBjgTvVtqJ8/1uKGuAurKoba82WuFN/jUVXi32nJvgPJ9rquA+qrGfZXeBU/NxlmiXTXerfpQ2yUD1sYtoJFCVifYYsAeA7bG9/7N6d3iGt9bnWBzeL+zOlrua/bqsdhxWey4sOEy7LgMOw3YaMD32rTVm1bqsdFgWr3vPRYaPCYut0mDx0NDY2BqcJu43J7G73z7vAHL5fG0+N7duM/VeB2Xx3fcEe/dzY8xcXuawpnLbeI22w5qHhPvWm/h63GMSoZBi3DUPBhZLM3eN9tvGM1C1xHHGM3CndHiXFp8HpaZyK8uOD5iP3fEw80LL7zAvHnzWLp0KVOmTOGRRx5h5syZbNu2jYyMjKOO/+ijj7jyyitZvHgx3/ve93juuee46KKL2LhxI6NHj47AT+C1r3IfsbZYhqYOxWaJ+K9VwijZmczpA07n9AGnA94urJ1lO9l6aCtbS7ay/fB28svyOVR7iIM1BzlYc5D1xetbvVacLY702HTSY9NJi0kj2ZlMsjOZFGcKKc4UEh2JxNvjSXQkkmBPIMGRQJwtjhhbTPQPYrdYITbFu3WFaYKr1juuqL7SG3wafO+rG99XtXxtqGn5vqGm8RpV3tfm+xpqWna9ueu8WxBbniy0MlC8M6wO7yr2VlvL9xY7WO2tfLZ631ts3s+25p9t3ldL4/EWa7PPtsbPzfYZzT9b8RiNm2nBbbF5XzHwYMGFFQ8W3I2by2x8b1pwGwZu07vP0/idC7z7PAYuLLhNw7vfNLzvseAy8e7zGN79GI1hzcTTGLbcpjfAuU0Tj++zB//7pn3NzvGAx2y8TrP9TddsPP+I7zwmTe+b3cs0m+5hNh7ju7fvvI4wTXCZJmCGNSjWuVpfhDhcIt4tNWXKFE444QT+8Ic/AODxeMjNzeVHP/oRd9xxx1HHz5o1i6qqKl577TX/vu985zuMHz+epUuXHvN+oRxz4/a4Ka0rpU9sn6BeV6JDRX0F+WX55Jfns7dyL0VVRRRWFbKvch9FVUXUumu7dP1YW6x/c1qdOK1OYmwxOK1OHFYHDosDu9WOw+LAYXVgs9iwW+zYLXZsFlvTZtiwWqxYjcbN0vRqMSxYDWtjU3nTZ4thwcDwvhreVwve9waGf5+Bt4ncv7/xHAyavms83ve5+fG+7337mn9udnibxzbX/Ppt7mulRb+18446xuMBV2OocdU2e1/vfW2oA09d4/76Zq/14G7wH2t4Grz7PI3fNf/sbmja527cPC7vdTwN4HZ595ldHOPUK1jAaG0z2tjv244898jjDbz/427+XeNri/3Njmtrv/9z43/EGBZMw8DEu4GBx/+/TQMTMA2L9xUDE+97/3ctXg1Ms/F94z1NDDym75imDdNoXB3P+2qaTff13cMDOBL7MfW824P6T6nHdEvV19ezYcMG5s+f799nsViYMWMGa9asafWcNWvWMG/evBb7Zs6cySuvvNLq8XV1ddTVNa3bU15e3vXC22C1WBVspE2JjkTG9B3DmL5jjvrONE2qXdX+lp0DNQc4XHuYsroy/1ZaV0plQyUV9RVUNlRSVV9FZUMlZuO/smpcNdS4asL9Y0mk+PrKWrAAMY2bhIencQsC84jXHmxcjZ2pBDfcdEZEw83Bgwdxu91kZma22J+ZmclXX33V6jlFRUWtHl9UVNTq8YsXL+auu+4KTsEiIWIYBvH2eOLt8QxMGtjh80zTpNZdS42rhuqGaqpd1dS6aqlz1/lfa1w1NHgaqHfXezdPPQ3uBho8DbhMFw3uBlweFy7Thdvj9r93eVx4TA9ujxu36d08pse7z3Q3Npt7Xz2mBw+epn2YTfsbx7eYeD+bpukPZL73zY/x/t/RxzT/mX3HNj/G//2RfxlMjvruqGOOuEabv+8g/dXpjs9xBOtn6/7MFi+dOifQ+3X5VxuifzZm6K5tj0sLyXU7KuoHh8yfP79FS095eTm5ubkRrEgkeAzD8HdFpcVE9l8mIiLdRUTDTXp6OlarleLiljOgFhcXk5WV1eo5WVlZnTre6XTidAZnNlQRERHp/iL6eIXD4WDSpEmsWrXKv8/j8bBq1SqmTp3a6jlTp05tcTzAypUr2zxeREREepeId0vNmzeP66+/nsmTJ3PiiSfyyCOPUFVVxezZswG47rrr6NevH4sXLwbg1ltv5dRTT+Whhx7ivPPOY/ny5axfv55ly5ZF8scQERGRbiLi4WbWrFkcOHCAhQsXUlRUxPjx41mxYoV/0HBBQQEWS1MD07Rp03juuee48847+cUvfsHQoUN55ZVXIjrHjYiIiHQfEZ/nJtxCOc+NiIiIhEZn/n5H+ZSmIiIi0tso3IiIiEhUUbgRERGRqKJwIyIiIlFF4UZERESiisKNiIiIRBWFGxEREYkqCjciIiISVRRuREREJKpEfPmFcPNNyFxeXh7hSkRERKSjfH+3O7KwQq8LNxUVFQDk5uZGuBIRERHprIqKCpKTk9s9ptetLeXxeNi3bx+JiYkYhhHwdcrLy8nNzWX37t1aoyoM9PsOL/2+w0u/7/DS7zu8gvX7Nk2TiooKcnJyWiyo3Zpe13JjsVjo379/0K6XlJSk/+cII/2+w0u/7/DS7zu89PsOr2D8vo/VYuOjAcUiIiISVRRuREREJKoo3ATI6XSyaNEinE5npEvpFfT7Di/9vsNLv+/w0u87vCLx++51A4pFREQkuqnlRkRERKKKwo2IiIhEFYUbERERiSoKNyIiIhJVFG4CtGTJEgYNGkRMTAxTpkzhk08+iXRJUWnx4sWccMIJJCYmkpGRwUUXXcS2bdsiXVavcd9992EYBrfddlukS4lae/fu5ZprrqFPnz7ExsYyZswY1q9fH+myopLb7WbBggXk5eURGxvL4MGD+fWvf92htYrk2N577z3OP/98cnJyMAyDV155pcX3pmmycOFCsrOziY2NZcaMGXz99dchqUXhJgAvvPAC8+bNY9GiRWzcuJFx48Yxc+ZM9u/fH+nSos67777LnDlz+Pjjj1m5ciUNDQ2cddZZVFVVRbq0qLdu3Tr++Mc/Mnbs2EiXErUOHz7MSSedhN1u59///jdffvklDz30EKmpqZEuLSrdf//9PPHEE/zhD39g69at3H///TzwwAM89thjkS4tKlRVVTFu3DiWLFnS6vcPPPAAv//971m6dClr164lPj6emTNnUltbG/xiTOm0E0880ZwzZ47/s9vtNnNycszFixdHsKreYf/+/SZgvvvuu5EuJapVVFSYQ4cONVeuXGmeeuqp5q233hrpkqLSz3/+c/Pkk0+OdBm9xnnnnWd+//vfb7HvkksuMa+++uoIVRS9APMf//iH/7PH4zGzsrLM3/72t/59paWlptPpNJ9//vmg318tN51UX1/Phg0bmDFjhn+fxWJhxowZrFmzJoKV9Q5lZWUApKWlRbiS6DZnzhzOO++8Fv87l+B79dVXmTx5MpdddhkZGRlMmDCBJ598MtJlRa1p06axatUqtm/fDsCnn37KBx98wDnnnBPhyqLfrl27KCoqavHvlOTkZKZMmRKSv529buHMrjp48CBut5vMzMwW+zMzM/nqq68iVFXv4PF4uO222zjppJMYPXp0pMuJWsuXL2fjxo2sW7cu0qVEvW+++YYnnniCefPm8Ytf/IJ169bx4x//GIfDwfXXXx/p8qLOHXfcQXl5OSNGjMBqteJ2u/nNb37D1VdfHenSol5RURFAq387fd8Fk8KN9Bhz5szh888/54MPPoh0KVFr9+7d3HrrraxcuZKYmJhIlxP1PB4PkydP5t577wVgwoQJfP755yxdulThJgRefPFFnn32WZ577jmOP/54Nm/ezG233UZOTo5+31FG3VKdlJ6ejtVqpbi4uMX+4uJisrKyIlRV9Js7dy6vvfYa77zzDv379490OVFrw4YN7N+/n4kTJ2Kz2bDZbLz77rv8/ve/x2az4Xa7I11iVMnOzmbUqFEt9o0cOZKCgoIIVRTdfvrTn3LHHXdwxRVXMGbMGK699lpuv/12Fi9eHOnSop7v72O4/nYq3HSSw+Fg0qRJrFq1yr/P4/GwatUqpk6dGsHKopNpmsydO5d//OMfvP322+Tl5UW6pKh2xhlnsGXLFjZv3uzfJk+ezNVXX83mzZuxWq2RLjGqnHTSSUdNbbB9+3YGDhwYoYqiW3V1NRZLyz97VqsVj8cToYp6j7y8PLKyslr87SwvL2ft2rUh+dupbqkAzJs3j+uvv57Jkydz4okn8sgjj1BVVcXs2bMjXVrUmTNnDs899xz//Oc/SUxM9PfNJicnExsbG+Hqok9iYuJR45ni4+Pp06ePxjmFwO233860adO49957ufzyy/nkk09YtmwZy5Yti3RpUen888/nN7/5DQMGDOD4449n06ZNPPzww3z/+9+PdGlRobKykh07dvg/79q1i82bN5OWlsaAAQO47bbbuOeeexg6dCh5eXksWLCAnJwcLrroouAXE/Tnr3qJxx57zBwwYIDpcDjME0880fz4448jXVJUAlrd/vznP0e6tF5Dj4KH1r/+9S9z9OjRptPpNEeMGGEuW7Ys0iVFrfLycvPWW281BwwYYMbExJjHHXec+ctf/tKsq6uLdGlR4Z133mn139fXX3+9aZrex8EXLFhgZmZmmk6n0zzjjDPMbdu2haQWwzQ1NaOIiIhED425ERERkaiicCMiIiJRReFGREREoorCjYiIiEQVhRsRERGJKgo3IiIiElUUbkRERCSqKNyIiIhIVFG4ERERkaiicCMiIiJRReFGREREoorCjYj0eAcOHCArK4t7773Xv++jjz7C4XCwatWqCFYmIpGghTNFJCq88cYbXHTRRXz00UcMHz6c8ePHc+GFF/Lwww9HujQRCTOFGxGJGnPmzOGtt95i8uTJbNmyhXXr1uF0OiNdloiEmcKNiESNmpoaRo8eze7du9mwYQNjxoyJdEkiEgEacyMiUWPnzp3s27cPj8dDfn5+pMsRkQhRy42IRIX6+npOPPFExo8fz/Dhw3nkkUfYsmULGRkZkS5NRMJM4UZEosJPf/pTXnrpJT799FMSEhI49dRTSU5O5rXXXot0aSISZuqWEpEeb/Xq1TzyyCP89a9/JSkpCYvFwl//+lfef/99nnjiiUiXJyJhppYbERERiSpquREREZGoonAjIiIiUUXhRkRERKKKwo2IiIhEFYUbERERiSoKNyIiIhJVFG5EREQkqijciIiISFRRuBEREZGoonAjIiIiUUXhRkRERKLK/wczjTJ+CHfmlQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fdata = pd.DataFrame({'x': np.arange(0.1, 10, .1)})\n", "df = [[1, 1], [1, 50], [10, 50]]\n", "labels = []\n", "for df1, df2 in df:\n", " label = f'f_{df1}_{df2}'\n", " fdata[label] = f.pdf(fdata.x, df1, df2)\n", " sns.lineplot(data=fdata, x='x', y=label, palette='colorblind')\n", " labels.append(label)\n", "\n", "plt.legend(labels)\n", "plt.ylabel('density')" ] }, { "cell_type": "markdown", "id": "6039300f", "metadata": {}, "source": [ "## ANOVA result\n", "\n", "Statsmodels automatically generates dummy variables when it sees that an independent variable is categorical." ] }, { "cell_type": "code", "execution_count": 119, "id": "687e7cd1-ba2b-439d-8cb3-1ca99709e0e2", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: SysBP R-squared: 0.161\n", "Model: OLS Adj. R-squared: 0.145\n", "Method: Least Squares F-statistic: 10.07\n", "Date: Thu, 01 Sep 2022 Prob (F-statistic): 0.000100\n", "Time: 12:52:54 Log-Likelihood: -413.59\n", "No. Observations: 108 AIC: 833.2\n", "Df Residuals: 105 BIC: 841.2\n", "Df Model: 2 \n", "Covariance Type: nonrobust \n", "====================================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "------------------------------------------------------------------------------------\n", "Intercept 131.9393 1.883 70.066 0.000 128.205 135.673\n", "group[T.drug2] 10.4696 2.663 3.931 0.000 5.189 15.750\n", "group[T.placebo] 10.2214 2.663 3.838 0.000 4.941 15.502\n", "==============================================================================\n", "Omnibus: 0.439 Durbin-Watson: 2.146\n", "Prob(Omnibus): 0.803 Jarque-Bera (JB): 0.587\n", "Skew: 0.021 Prob(JB): 0.746\n", "Kurtosis: 2.641 Cond. No. 3.73\n", "==============================================================================\n", "\n", "Notes:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n" ] } ], "source": [ "lmresultAnovaBasic = smf.ols(formula='SysBP ~ group', data=DrugDf).fit()\n", "\n", "print(lmresultAnovaBasic.summary())" ] } ], "metadata": { "Rmd_chunk_options": { "output": { "bookdown::gitbook": { "includes": { "in_header": "google_analytics.html" }, "lib_dir": "book_assets" }, "html_document": "default", "pdf_document": "default" } }, "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.9.15" }, "sos": { "kernels": [ [ "SoS", "sos", "", "" ], [ "R", "ir", "", "" ], [ "Markdown", "markdown", "", "" ] ] } }, "nbformat": 4, "nbformat_minor": 5 }