{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Introduction to mathematical statistics \n", "\n", "Welcome to the lecture 11 in 02403\n", "\n", "During the lectures we will present both slides and notebooks. " ] }, { "cell_type": "code", "execution_count": 89, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import pandas as pd\n", "import scipy.stats as stats\n", "import statsmodels.api as sm\n", "import statsmodels.formula.api as smf\n", "import statsmodels.stats.power as smp\n", "import statsmodels.stats.proportion as smprop" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Two-way ANOVA" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To introduce two-way ANOVA we start by considering the following data (and do some data visuallisation)" ] }, { "cell_type": "code", "execution_count": 90, "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", " \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", " \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", "
valuetreatmentblock
02.8A1
13.6A2
23.4A3
32.3A4
45.5B1
56.3B2
66.1B3
75.7B4
85.8C1
98.3C2
106.9C3
116.1C4
\n", "
" ], "text/plain": [ " value treatment block\n", "0 2.8 A 1\n", "1 3.6 A 2\n", "2 3.4 A 3\n", "3 2.3 A 4\n", "4 5.5 B 1\n", "5 6.3 B 2\n", "6 6.1 B 3\n", "7 5.7 B 4\n", "8 5.8 C 1\n", "9 8.3 C 2\n", "10 6.9 C 3\n", "11 6.1 C 4" ] }, "execution_count": 90, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Typing the data:\n", "data = pd.DataFrame({\n", " 'value': [2.8, 3.6, 3.4, 2.3, 5.5, 6.3, 6.1, 5.7, 5.8, 8.3, 6.9, 6.1], \n", " 'treatment': np.repeat([\"A\", \"B\", \"C\"], 4),\n", " 'block': np.tile([\"1\", \"2\", \"3\", \"4\"], 3)})\n", "data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Visualising the data with box plots:" ] }, { "cell_type": "code", "execution_count": 91, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0.5, 0.98, 'Boxplots')" ] }, "execution_count": 91, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0YAAAHeCAYAAAC/jhsfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABJy0lEQVR4nO3deXxU9b3/8fckkwxZSMIShGjAAMqETS5BJIICFRAEr/a64NJbUOtPKy5UtFV7q2BV3LDYqPzcClZLoVKXFsqScgUKAYQgyhaUEIEfi2FLQhIYZjLn9wc3cwkhMJM5k8nMeT0fjzzCnJz5ns85ZzjfeZ/vzDk2wzAMAQAAAICFxYS7AAAAAAAIN4IRAAAAAMsjGAEAAACwPIIRAAAAAMsjGAEAAACwPIIRAAAAAMsjGAEAAACwPIIRAAAAAMsjGAEAAACwPIIRACAq2Gw2TZ48OdxlAAAiFMEIAHBes2bNks1mq/PTrl07DR06VAsXLgx3eUHZunWrJk+erO+//z7cpQAAwsge7gIAAJHj2WefVVZWlgzD0A8//KBZs2bpuuuu09///neNGTMm3OU1ytatWzVlyhQNGTJEF198cbjLAQCECcEIAOC3UaNGqV+/fr7H99xzjy644AL9+c9/jthgBACAxEfpAABBSEtLU0JCguz2/z3PVlVVpUmTJikzM1MOh0PdunXTq6++KsMwJEnHjx+X0+mU0+nU8ePHfc87cuSIOnTooCuvvFI1NTWSpPHjxys5OVk7d+7Utddeq6SkJGVkZOjZZ5/1tXcuX331lUaNGqWUlBQlJyfrmmuu0Zo1a3x/nzVrlm655RZJ0tChQ30fE1y2bJkkaf369br22mvVtm1bJSQkKCsrS3fffXfQ2w0A0PwwYgQA8Ft5ebkOHTokwzBUWlqqvLw8VVZW6ic/+YkkyTAM/fu//7u++OIL3XPPPerTp48WL16sxx9/XHv37tXvfvc7JSQk6IMPPtDAgQP161//Wq+99pokacKECSovL9esWbMUGxvrW2ZNTY1GjhypAQMG6OWXX9aiRYv0zDPPyOPx6Nlnn22w1i1btuiqq65SSkqKfvnLXyouLk5vv/22hgwZouXLl+uKK67Q1VdfrYcffli///3v9dRTTyk7O1uSlJ2drdLSUo0YMULp6el64oknlJaWpu+//16ffPJJCLcwACBsDAAAzmPmzJmGpHo/DofDmDVrlm++zz77zJBkPPfcc3Wef/PNNxs2m83YsWOHb9qTTz5pxMTEGCtWrDA+/vhjQ5Ixffr0Os8bN26cIcl46KGHfNO8Xq8xevRoIz4+3jh48KBvuiTjmWee8T2+8cYbjfj4eKO4uNg3bd++fUbLli2Nq6++2jetdtlffPFFnWV/+umnhiRj3bp1gW0sAEBE4qN0AAC/vfnmm8rPz1d+fr4++ugjDR06VD/72c98oyj/+Mc/FBsbq4cffrjO8yZNmiTDMOpcwW7y5Mnq0aOHxo0bpwceeECDBw+u97xaDz74oO/fNptNDz74oE6ePKl//vOfZ52/pqZGS5Ys0Y033qjOnTv7pnfo0EF33HGHVq5cqYqKinOua1pamiRp/vz5crvd55wXABD5CEYAAL/1799fw4YN07Bhw3TnnXdqwYIF6t69uy+o7Nq1SxkZGWrZsmWd59V+RG3Xrl2+afHx8frDH/6gkpISHTt2TDNnzpTNZqu3zJiYmDrhRpIuvfRSSWrwEtsHDx5UdXW1unXrVu9v2dnZ8nq92rNnzznXdfDgwbrppps0ZcoUtW3bVjfccINmzpwpl8t1zucBACITwQgA0GgxMTEaOnSo9u/fr++++y7g5y9evFiSdOLEiUY9P5RsNpvmzZun1atX68EHH9TevXt19913KycnR5WVleEuDwBgMoIRACAoHo9HklRZWalOnTpp3759OnbsWJ15ioqKJEmdOnXyTfvmm2/07LPP6q677tK//du/6Wc/+5nKy8vrte/1erVz584607799ltJavC+Q+np6UpMTNT27dvr/a2oqEgxMTHKzMyUpLOOUp1uwIABev7557V+/Xr96U9/0pYtWzRnzpxzPgcAEHkIRgCARnO73VqyZIni4+OVnZ2t6667TjU1NXrjjTfqzPe73/1ONptNo0aN8j1v/PjxysjI0Ouvv65Zs2bphx9+0C9+8YuzLuf09gzD0BtvvKG4uDhdc801Z50/NjZWI0aM0Oeff17n43Y//PCDZs+erUGDBiklJUWSlJSUJEkqKyur08bRo0frXRK8T58+ksTH6QAgCnG5bgCA3xYuXOgb/SktLdXs2bP13Xff6YknnlBKSoquv/56DR06VL/+9a/1/fff67LLLtOSJUv0+eefa+LEierSpYsk6bnnntPGjRu1dOlStWzZUr1799bTTz+t//qv/9LNN9+s6667zrfMFi1aaNGiRRo3bpyuuOIKLVy4UAsWLNBTTz2l9PT0Bmt97rnnlJ+fr0GDBumBBx6Q3W7X22+/LZfLpZdfftk3X58+fRQbG6uXXnpJ5eXlcjgc+tGPfqTZs2frrbfe0o9//GN16dJFx44d07vvvquUlJQ69QEAokSYr4oHAIgAZ7tcd4sWLYw+ffoYM2bMMLxer2/eY8eOGb/4xS+MjIwMIy4uzrjkkkuMV155xTdPYWGhYbfb61yC2zAMw+PxGJdffrmRkZFhHD161DCMU5frTkpKMoqLi40RI0YYiYmJxgUXXGA888wzRk1NTZ3n64zLdRuGYWzYsMG49tprjeTkZCMxMdEYOnSoUVBQUG/93n33XaNz585GbGys79LdGzZsMG6//XajY8eOhsPhMNq1a2eMGTPGWL9+vQlbFADQ3NgMw49bhwMAEAbjx4/XvHnzuNgBACDk+I4RAAAAAMsjGAEAAACwPIIRAAAAAMsjGAFBmDVrlmw2W53LAQMwz6xZs/h+EdBE6NNgdQQjAAAAAJZHMAIAAABgeQQjAAAAAJZHMIKlzJs3TzabTcuXL6/3t7fffls2m02bN2/WN998o/Hjx6tz585q0aKF2rdvr7vvvluHDx8+7zJsNpsmT55cb/rFF1+s8ePH15lWVlamiRMnKjMzUw6HQ127dtVLL70kr9fb2FUEAFgEfRpgLnu4CwCa0ujRo5WcnKy//OUvGjx4cJ2/zZ07Vz169FDPnj01bdo07dy5U3fddZfat2+vLVu26J133tGWLVu0Zs0a2Wy2oGuprq7W4MGDtXfvXt13333q2LGjCgoK9OSTT2r//v2aPn160MsAAEQv+jTAXAQjWEpCQoKuv/56zZs3T7///e8VGxsrSTpw4ICWL1/uOyv2wAMPaNKkSXWeO2DAAN1+++1auXKlrrrqqqBree2111RcXKyvvvpKl1xyiSTpvvvuU0ZGhl555RVNmjRJmZmZQS8HABCd6NMAc/FROljO2LFjVVpaqmXLlvmmzZs3T16vV2PHjpV0qrOpdeLECR06dEgDBgyQJG3YsMGUOj7++GNdddVVatWqlQ4dOuT7GTZsmGpqarRixQpTlgMAiF70aYB5GDGC5YwcOVKpqamaO3eurrnmGkmnPnLQp08fXXrppZKkI0eOaMqUKZozZ45KS0vrPL+8vNyUOr777jt98803Sk9PP+vfz1wuAABnok8DzEMwguU4HA7deOON+vTTT/XWW2/phx9+0KpVq/TCCy/45rn11ltVUFCgxx9/XH369FFycrK8Xq9GjhzZ6C+R1tTU1Hns9Xo1fPhw/fKXvzzr/LUdGgAADaFPA8xDMIIljR07Vh988IGWLl2qbdu2yTAM30cOjh49qqVLl2rKlCl6+umnfc/57rvv/Gq7VatWKisrqzPt5MmT2r9/f51pXbp0UWVlpYYNGxbcygAALI0+DTAH3zGCJQ0bNkytW7fW3LlzNXfuXPXv319ZWVmS5PvyqmEYdZ7j7xV1unTpUu+z1O+88069s2u33nqrVq9ercWLF9dro6ysTB6Px9/VAQBYGH0aYA5GjGBJcXFx+o//+A/NmTNHVVVVevXVV31/S0lJ0dVXX62XX35ZbrdbF154oZYsWaKSkhK/2v7Zz36m+++/XzfddJOGDx+ur7/+WosXL1bbtm3rzPf444/rb3/7m8aMGaPx48crJydHVVVV2rRpk+bNm6fvv/++3nMAADgTfRpgDoIRLGvs2LF67733ZLPZdOutt9b52+zZs/XQQw/pzTfflGEYGjFihBYuXKiMjIzztnvvvfeqpKRE77//vhYtWqSrrrpK+fn5vi/F1kpMTNTy5cv1wgsv6OOPP9Yf//hHpaSk6NJLL9WUKVOUmppq6voCAKIXfRoQPJtx5tgqAAAAAFgM3zECAAAAYHkEIwAAAACWRzACAAAAYHkEIwAAAACWRzACAAAAYHkEIwAAAACW1+T3MfJ6vdq3b59atmwpm83W1IsHAPjBMAwdO3ZMGRkZionhHNrZ0J8BQGTwt09r8mC0b98+ZWZmNvViAQCNsGfPHl100UXhLqNZoj8DgMhyvj6tyYNRy5YtJZ0qLCUlpakX3+TcbreWLFmiESNGKC4uLtzlIATYx9Zgtf1cUVGhzMxM3zEb9UVif2a113E4sa2bDtu66UTqtva3T2vyYFT7cYOUlJSI6UiC4Xa7lZiYqJSUlIh6AcF/7GNrsOp+5iNiDYvE/syqr+NwYFs3HbZ104n0bX2+Po0PjgMAAACwPIIRAAAAAMsjGAEAAACwPIIRAAAAAMsjGAEAAACwPIIRAAAAAMsjGAEAAACwPIIRAAAAAMsjGAEAAACwPIIRAAAAAMsjGAEAAACwPIIRAAAAAMsjGAEAAACwPHu4CwCAcKuurlZRUdE556k87lLBpmK1arteyQmOc87rdDqVmJhoZokAmgl/jhdSYMeMWhw7gPAiGAGwvKKiIuXk5Pg178t+zFNYWKi+ffsGVxSAZimQ44Xk3zGjFscOILwIRgAsz+l0qrCw8JzzbN9fpkc/3qTXbumlbh3SztsegOjkz/FCCuyYcXrbAMKHYATA8hITE897ljZm12E5/nVc2T0vU59ObZqoMgDNjT/HC4ljBhCJuPgCAAAAAMsjGAEAAACwPIIRAAAAAMsjGAEAAACwPIIRAAAAAMsjGAEAAACwPIIRAAAAAMsjGAEAAACwPIIRAAAAAMsjGAEAAACwvICCUU1NjX7zm98oKytLCQkJ6tKli37729/KMIxQ1QcAgOnozwAAZ7IHMvNLL72kGTNm6IMPPlCPHj20fv163XXXXUpNTdXDDz8cqhoBADAV/RkA4EwBBaOCggLdcMMNGj16tCTp4osv1p///Gd9+eWXISkOAIBQoD8DAJwpoGB05ZVX6p133tG3336rSy+9VF9//bVWrlyp1157rcHnuFwuuVwu3+OKigpJktvtltvtbmTZkaN2Ha2wrlbFPrYGj8fj+22FfR3t62jV/ozjVdOx2jEjnHhdN51I3db+1htQMHriiSdUUVEhp9Op2NhY1dTU6Pnnn9edd97Z4HOmTp2qKVOm1Ju+ZMkSJSYmBrL4iJafnx/uEhBi7OPotqdSkuxas2aN9m4OdzWhV11dHe4SQsrq/RnHq9Cz2jGjOeB13XQibVv726fZjAC+aTpnzhw9/vjjeuWVV9SjRw9t3LhREydO1GuvvaZx48ad9TlnO8OWmZmpQ4cOKSUlxd9FRyy32638/HwNHz5ccXFx4S4HIcA+toavdx/Rze+u17x7++myjq3DXU7IVVRUqG3btiovL4/KY7VV+zOOV03HaseMcOJ13XQidVv726cFNGL0+OOP64knntBtt90mSerVq5d27dqlqVOnNtiROBwOORyOetPj4uIiaoMGy2rra0Xs4+hmt9t9v62wn6N9Ha3en0VizZHGaseM5oDXddOJtG3tb60BXa67urpaMTF1nxIbGyuv1xtIMwAAhBX9GQDgTAGNGF1//fV6/vnn1bFjR/Xo0UNfffWVXnvtNd19992hqg8AANPRnwEAzhRQMMrLy9NvfvMbPfDAAyotLVVGRobuu+8+Pf3006GqDwAA09GfAQDOFFAwatmypaZPn67p06eHqBwAAEKP/gwAcKaAvmMEAAAAANGIYAQAAADA8ghGAAAAACyPYAQAAADA8ghGAAAAACyPYAQAAADA8ghGAAAAACyPYAQAAADA8ghGAAAAACyPYAQAAADA8ghGAAAAACyPYAQAAADA8ghGAAAAACyPYAQAAADA8ghGAAAAACyPYAQAAADA8ghGAAAAACyPYAQAAADA8ghGAAAAACyPYAQAAADA8ghGAAAAACyPYAQAAADA8ghGAAAAACyPYAQAAADA8ghGAAAAACyPYAQAAADA8ghGAAAAACyPYAQAAADA8ghGAAAAACyPYAQAAADA8ghGAAAAACyPYAQAAADA8uzhLgAAAKA5KDlUpSqXx5S2ig9W+X7b7ea93Upy2JXVNsm09gD8L4IRAACwvJJDVRr66jLT2500b5PpbX7x2BDCERACBCMAUc2sM8Bmnv3ljC/Q/NQeJ6aP7aOu7ZKDb++4S/OXrdaYIblKSnAE3Z4k7Sit1MS5G00b1QJQF8EIQNQKxRlgs87+csYXaJ66tktWzwtTg27H7XbrQLrUt1MrxcXFmVAZgFAjGAGIWmaeATbr7C9nfAEAaJ4IRgCinhlngDn7CwBAdONy3QAAAAAsj2AEAAAAwPIIRgAAAAAsj2AEAAAAwPIIRgAAAAAsj2AEAAAAwPIIRgAAAAAsL6BgdPHFF8tms9X7mTBhQqjqAwAgJOjTAACnC+gGr+vWrVNNTY3v8ebNmzV8+HDdcsstphcGAEAo0acBAE4XUDBKT0+v8/jFF19Uly5dNHjwYFOLAgAg1OjTAACna/R3jE6ePKmPPvpId999t2w2m5k1AQDQpOjTAAABjRid7rPPPlNZWZnGjx9/zvlcLpdcLpfvcUVFhSTJ7XbL7XY3dvERo3YdrbCuVsU+br48Ho/vd7D7x6z9bGZNodScawsFf/q0aOjPOF41zOz/m6HY1pFy/GhqvK6bTqRua3/rtRmGYTRmAddee63i4+P197///ZzzTZ48WVOmTKk3ffbs2UpMTGzMogHAL3sqpVc32fVYL48yk8NdzSnNsaazqa6u1h133KHy8nKlpKSEu5yQ86dPoz+LbpHwfzMSagSaI3/7tEYFo127dqlz58765JNPdMMNN5xz3rOdYcvMzNShQ4cs0dm63W7l5+dr+PDhiouLC3c5CAH2cfO1ZV+FbpyxRp/9fIB6ZAR3vDFrP5tZUyhVVFSobdu2lghG/vZp0dCfcbxqmNn/N0OxrSPl+NHUeF03nUjd1v72aY36KN3MmTPVrl07jR49+rzzOhwOORyOetPj4uIiaoMGy2rra0Xs4+bHbrf7fpu1b4Ldz6GoKRSac21m87dPi6b+LBJrDrVQ/d80c1tHyvEjXHhdN51I29b+1hrwxRe8Xq9mzpypcePG+f6DAgAQiejTAAC1Au4F/vnPf2r37t26++67Q1EPAABNhj4NaL6qq6tVVFR03vkqj7tUsKlYrdquV3JC/VHdMzmdTr4XiLMKOBiNGDFCjbxeAwAAzQp9GtB8FRUVKScnx+/5X/ZzvsLCQvXt27dxRSGq8bkBAAAANDtOp1OFhYXnnW/7/jI9+vEmvXZLL3XrkOZXu8DZEIwAAADQ7CQmJvo1shOz67Ac/zqu7J6XqU+nNk1QGaJVwBdfAAArWntgrV6veF1rD6wNdykAACAECEYAcB6GYShvY54Oeg8qb2Me30kBACAKEYwA4DwK9hVo65GtkqStR7aqYF9BmCsCAABmIxgBwDkYhqG8r/IUYzt1uIyxxSjvK0aNAACINgQjADiHgn0F2nJ4i7yGV5LkNbzacngLo0YAAEQZghEANODM0aJajBoBABB9CEYA0IAzR4tqMWoEAED0IRgBwFnUjhbZZDvr322yMWoEAEAUIRgBwFm4vW4dqDogQ2cPPoYMHag6ILfX3cSVAQCAULCHuwAAaI7iY+M1Z8wcHTlxRJLk8Xi0auUqDRw0UHb7qUNn6xatFR8bH84yAQCASQhGANCA9knt1T6pvSTJ7XarxF6i7NbZiouLC3NlAADAbHyUDgAAAIDlEYwAAAAAWB7BCAAAAIDlEYwAAAAAWB7BCAAAAIDlEYwAAAAAWB7BCAAAAIDlEYwAAAAAWB43eAUQ1Wz2CpVUbFdMi+Sg2vF4PNrn2adtR7bJbm/8obOkolI2e0VQtQBApCs5VKUql8eUtooPVvl+B3N8PlOSw66stkmmtYfmj2AEIKrFpa3VU1++YFp7by16K+g24tKukXRd8MUAQAQqOVSloa8uM73dSfM2md7mF48NIRxZCMEIQFRzl12haaPvUJd2wY8YrVq5SgMHDQzqjGRxaaUe/lNxULUAQCSrHSmaPraPugZ5bJakquMuzV+2WmOG5CopwRF0e5K0o7RSE+duNG1UC5GBYAQgqhmeFGWldFP3NqlBteN2u1ViL1F262zFxcU1uh3viXIZnoNB1QIA0aBru2T1vDC4Y7N06vh8IF3q26lVUMdngIsvAAAAALA8ghEAAAAAyyMYAQAAALA8ghEAAAAAyyMYAQAAALA8ghEAAAAAyyMYAQAAALA87mMEnEN1dbWKiorOOU/lcZcKNhWrVdv1SvbjxnJOp1OJiYlmlQgAMInNXqGSiu2KaRH8TUc9Ho/2efZp25FtQd0U+nQlFZWy2StMaQtAfQQj4ByKioqUk5Pj17wv+9lmYWGh+vbt2/iiAAAhEZe2Vk99+YKpbb616C1T24tLu0bSdaa2CeAUghFwDk6nU4WFheecZ/v+Mj368Sa9dksvdeuQ5lebAIDmx112haaNvkNd2pkzYrRq5SoNHDTQtBGj4tJKPfynYlPaAlAfwQg4h8TExPOO7sTsOizHv44ru+dl6tOpTRNVBgAwm+FJUVZKN3Vvkxp0W263WyX2EmW3zlZcXJwJ1UneE+UyPAdNaQtAfVx8AQAAAIDlEYwAAAAAWB7BCAAAAIDlEYwAAAAAWB7BCAAAAIDlcVU6WFrJoSpVuTxBtVF8sMr324xLsiY57MpqmxR0OwAAAPAfwQiWVXKoSkNfXWZae5PmbTKtrS8eG0I4AgAAaEIEI1hW7UjR9LF91DWIm/lVHXdp/rLVGjMkV0kJjqBq2lFaqYlzNwY9igUACK+1B9bq9YrX1eZAGw3KHBTucgD4gWAEy+vaLlk9L2z8zfzcbrcOpEt9O7Uy7SZ+AIDIZRiG8jbm6aD3oPI25mngRQNls9nCXVazYrNXqKRiu2JaNP7EZC2Px6N9nn3admSbKR9pl6SSikrZ7BWmtBUJqqurVVRUdN75Ko+7VLCpWK3arleyHyeDnU6nEhMTzSixSRCMAAAATFSwr0Bbj2yVJG09slUF+wo08MKBYa6qeYlLW6unvnzB1DbfWvSWqe3FpV0j6TpT22yuioqKlJOT4/f8L/s5X2Fhofr27du4osIg4GC0d+9e/epXv9LChQtVXV2trl27aubMmerXr18o6gMAICTozxAKhmEo76s8xdhi5DW8irHFKO+rPF2ZcSWjRqdxl12haaPvUJcgPspey+PxaNXKVRo4aKBpI0bFpZV6+E/FprQVCZxOpwoLC8873/b9ZXr040167ZZe6tYhza92I0lAr56jR49q4MCBGjp0qBYuXKj09HR99913atWqVajqAwDAdPRnCJWCfQXacniL77HX8GrL4S2MGp3B8KQoK6Wburdp/EfZa7ndbpXYS5TdOtu0j7R7T5TL8Bw0pa1IkJiY6NfITsyuw3L867iye16mPp3aNEFlTSugYPTSSy8pMzNTM2fO9E3LysoyvSgAMMNxd40kafPe8qDbqjru0vqDUvtdR4O6yMaO0sqga0Hw6M8QCmeOFtVi1AiIDAEFo7/97W+69tprdcstt2j58uW68MIL9cADD+jee+9t8Dkul0sul8v3uKLi1BfZ3G633G53I8uOHLXraIV1jTQej8f3O5j9Y+Y+NqsmnPLt/lOB6IlPzLqUul0f7lhnSkuOWKNZ7+PmXJsZrNqf0Sc1zIzj75mjRbVqR41W7F6hKzOuDGuNzYHZ6xGK13W0bGuzRep28bfWgILRzp07NWPGDD366KN66qmntG7dOj388MOKj4/XuHHjzvqcqVOnasqUKfWmL1myJKKuUhGs/Pz8cJeAM+yplCS7Vq5cqV3Bf8TZlH1sdk2W55Zu62xTuwRD8THBNfXDcenDHXb9Z1ePLkgIri1HrLR17XJtDa6ZkKqurg53CSFl9f6MPqm+YI+/hmHo/1b+X9lkkyGj3t9tsmnqv6bq/uT7Gz1qFC19RKjWw8zXdbRsa7PVbpc1a9Zo7+ZwV+M/f/s0m2EY9f/3NiA+Pl79+vVTQUGBb9rDDz+sdevWafXq1Wd9ztnOsGVmZurQoUNKSUnxd9ERy+12Kz8/X8OHD+dSzs3Mln0VunHGGn328wHqkdH416KZ+9ismmC+r3cf0c3vrte8e/vpso6tw11OyFVUVKht27YqLy+PymO1Vfsz+qSGBXv8PVlzUtd9fp2OnDjS4DxtWrTRghsWKD42Piw1Nhdmr0coXtfRsq3NFql9ob99WkAjRh06dFD37t3rTMvOztZf//rXBp/jcDjkcNT/PH5cXJylDspWW99IUHvlGrvdHtS+8d3E73DwN/EzqyaYz2r7JtrX0er9WSTWHGrB/h+Pi4vT3DFzfcHobFdKa92itZJaJIWtxuYiVOth5us6Wra12SJ1u/hba0DBaODAgdq+fXudad9++606deoUSDNA1OAmfkBkoj9DKLRPaq/2Se0lheZKaQBCK6BP3f/iF7/QmjVr9MILL2jHjh2aPXu23nnnHU2YMCFU9QHN2tlu4geg+aM/AwCcKaARo8svv1yffvqpnnzyST377LPKysrS9OnTdeedd4aqPiCkbPYKlVRsV0yLwL9ZaRiGXln3imIUI6+8ilGMXln3itIcaY0eNSqpqJTNXtGo5wLwH/0ZAKsoOVSlKpfHlLaKD1b5fpt1M11JSnLYldW28R8zNUvAazRmzBiNGTMmFLUATS4uba2e+vIFU9ryyqvi8mLdtuC2IGu6RtJ1ptQEoGH0ZwCiXcmhKg19dZnp7U6aZ9ZtMP7XF48NCXs4Mi/qARHIXXaFpo2+Q13aBTZiZBiGfr3y1yopL5FXp93ETzHKSs3S84Oeb9SoUXFppR7+U3HAzwMAADhT7UjR9LF91DXA9zpnbe+4S/OXrdaYIblB3ez8dDtKKzVx7kbTRrWCQTCCpRmeFGWldFP3NqkBPW/V3lUqLq8fYGpHjcpcZRp44cCA6/GeKJfhORjw8wAAABrStV2yel4Y2Huds3G73TqQLvXt1CoqLyoS5C0PAesxDEN5X+XJprOPCNlkU95XeQrgFmEAAAAIM4IRECC3160DVQfOemdzSTJk6EDVAbm97iauDAAAAI3FR+mAAMXHxmvOmDnnvYlfY+9sjqZXXV2toqKic86zfX+ZXAd2aNvmBHkPp51zXqfTqcTERBMrBM7Nn9ewJFUed6lgU7FatV2vZD+/H8DrGYBVEIyARuAmftGlqKhIOTk5fs17xwfnn6ewsFB9+/YNsirAf4G8hiXp5QDa5vUMwCoIRgAsz+l0qrCw8JzzVB53acEXqzV6aO55z7Q7nU4zywPOy5/XsHRq5PPRjzfptVt6qVuHNL/bBgArIBgBsLzExMTznhF3u906eqhUuf37MTKIZsef17Akxew6LMe/jiu752Xq06lNE1QGAJGDiy8AAAAAsDxGjAAAAIAoZbNXqKRiu2JaBH+DV4/Ho32efdp2ZJvvglPBKqmolM1eYUpbwSIYAQAAAFEqLm2tnvryBVPbfGvRW6a2F5d2jaTrTG2zMQhGAAAAQJRyl12haaPvUJd25owYnXmLkmAVl1bq4T8Vm9JWsAhGsKzj7hpJ0ua95UG1U3XcpfUHpfa7jirJz/uCNGRHaWVQzwcQfUoOVanK5TGlreKDVb7fZr2pkaQkh11ZbZNMaw+AeQxPirJSuql7m9Sg2wrFLUq8J8pleA6a0lawCEawrOL/CSFPfLLJhNbs+nDHOhPaOSXJwX9NAKdC0dBXl5ne7qR5Zhz36vrisSGEIwARjXdfsKwRPU7doLVLu2QlxMU2up3t+8s1ad4mTbu5l7p1CP5sDGdeAdSqHSmaPraPuprwMZiq4y7NX7ZaY4bkBj3CXWtHaaUmzt1o2qgWAIQLwQiW1TopXrf17xh0Ox7PqTcDXdKT1PPC4IMRAJypa7tkU44vbrdbB9Klvp1acT8uADgD9zECAAAAYHkEIwAAAACWRzACAAAAYHkEIwAAAACWRzACAABAxFp7YK1er3hdaw+sDXcpiHBclQ44h+rqahUVFZ1znu37y+Q6sEPbNifIezjtvG06nU4lJiaaVCGAaGezV6ikYrtiWphz1/p9nn3admSbaTd4LamolM1eYUpbQKAMw1Dexjwd9B5U3sY8DbxooGw2W7jLQoQiGAHnUFRUpJycHL/mveMD/9osLCxU3759g6gKgJXEpa3VU1++YGqbby16y9T24tKukXSdqW0C/ijYV6CtR7ZKkrYe2aqCfQUaeOHAMFeFSEUwAs7B6XSqsLDwnPNUHndpwRerNXporpL9uGGi0+k0qzwAFuAuu0LTRt+hLibc4NXj8WjVylUaOGigaSNGxaWVevhPxaa0BQTCMAzlfZWnGFuMvIZXMbYY5X2VpyszrmTUCI1CMALOITEx8byjO263W0cPlSq3fz9umAjAdIYnRVkp3dS9jTk3eC2xlyi7dbZpxyvviXIZnoOmtAUEomBfgbYc3uJ77DW82nJ4C6NGaDQuvgAAAICIcvpo0elqR40MwwhTZYhkBCMAAABElNrRIq/hrTP99FEjIFAEIwAAAESM2tEim87+PSKbbIwaoVEIRgAAAIgYbq9bB6oOyNDZg48hQweqDsjtdTdxZYh0XHwBAAAAESM+Nl5zxszRkRNHJJ39aoutW7RWfGx8OMtEBCIYAQAAIKK0T2qv9kntJYXmaouwJj5KBwAAAMDyCEYAAAAALI9gBAAAAMDyCEYAAAAALI9gBAAAAMDyCEYAAAAALI9gBAAAAMDyCEYAAAAALI9gBAAAAMDyCEYAAAAALI9gBAAAAMDyCEYAAAAALI9gBAAAAMDyCEYAAAAALI9gBAAAAMDyAgpGkydPls1mq/PjdDpDVRsAACFDnwYAOJ090Cf06NFD//znP/+3AXvATQAA0CzQpwEAagXcA9jtdrVv3z4UtQAA0KTo0wAAtQIORt99950yMjLUokUL5ebmaurUqerYsWOD87tcLrlcLt/jiooKSZLb7Zbb7W5EyZGldh2tsK5WxT62BqvtZ6usZyB9Wjj6M4/H4/ttxjJC8To2u8ZwYVs3HbZ102Fbn+JvuzbDMAx/G124cKEqKyvVrVs37d+/X1OmTNHevXu1efNmtWzZ8qzPmTx5sqZMmVJv+uzZs5WYmOjvogEATai6ulp33HGHysvLlZKSEu5yQiLQPi0c/dmeSunVTXY91sujzOSQLCJokVCjPyJhPSKhRn9EwnpEQo3+iIT1aIoa/e3TAgpGZyorK1OnTp302muv6Z577jnrPGc7w5aZmalDhw5FbWd7Orfbrfz8fA0fPlxxcXHhLgchwD62Bqvt54qKCrVt2zaqg9GZztenhaM/27KvQjfOWKPPfj5APTKCX0YoXsdm1xgubOumw7ZuOmzrU/zt04L6lmlaWpouvfRS7dixo8F5HA6HHA5HvelxcXGWeHNRy2rra0XsY2uwyn62wjqe6Xx9Wjj6s9qLQdjtdlOXYWbNoaqxqbkNmySp6IcqUy7CUXXcpfUHpfb7KpWUUP910xjfHzkhKfK3Na/rpsO2PsXfdoP6n19ZWani4mL953/+ZzDNAAAQdvRp1lZcWilJeuKTTSa2ateHO9aZ2N4pSQ6ungiEQkD/sx577DFdf/316tSpk/bt26dnnnlGsbGxuv3220NVHwAAIUGfhtON6HHq6oRd2iUrIS426Pa27y/XpHmbNO3mXurWITXo9molOezKaptkWnsA/ldAwej//b//p9tvv12HDx9Wenq6Bg0apDVr1ig9PT1U9QEAEBL0aThd66R43da/4avsBqr2Sltd0pPU80LzghGA0AkoGM2ZMydUdQAA0KTo0wAAp4sJdwEAAAAAEG4EIwAAAACWRzACAAAAYHkEIwAAAACWRzACAAAAYHkEIwAAAACWx62TAQCwiLUH1ur1itfV5kAbDcocFO5yYFHH3TWSpM17y01pr+q4S+sPSu13HVVSgsOUNneUVprSTrSJ9mMIwQgAAAswDEN5G/N00HtQeRvzNPCigbLZbOEuCxZU/D+h44lPNpnYql0f7lhnYnunJDl4q1zLCscQ9jYAABZQsK9AW49slSRtPbJVBfsKNPDCgWGuClY0okd7SVKXdslKiIsNur3t+8s1ad4mTbu5l7p1SA26vVpJDruy2iaZ1l6ks8IxhGAEAECUMwxDeV/lKcYWI6/hVYwtRnlf5enKjCuj7owvmr/WSfG6rX9H09rzeDySpC7pSep5oXnBCP/LKscQLr4AAECUK9hXoC2Ht8hreCVJXsOrLYe3qGBfQZgrAxAJrHIMIRgBABDFTj/Te7raM76GYYSpMgCRwErHEIIRAABR7MwzvbWi9YwvAHNZ6RhCMAIAIErVnum16ezfAbDJFnVnfAGYx2rHEIIRAABRyu1160DVARk6+5sWQ4YOVB2Q2+tu4soARAKrHUO4Kh0AAFEqPjZec8bM0ZETRySdunrXqpWrNHDQQNntp94CtG7RWvGx8eEsE0AzZbVjCMEIAIAo1j6pvdonnbpvjNvtVom9RNmtsxUXFxfmygBEAisdQ/goHQAAAADLIxgBAAAAsDyCEQAAAADLIxgBAAAAsDyCEQAAAADLIxgBAAAAsDyCEQAAAADLIxgBAAAAsDyCEQAAAADLIxgBAAAAsDyCEQAAAADLIxgBAAAAsDyCEQAAAADLIxgBAAAAsDyCEQAAAADLIxgBAAAAsDyCEQAAAADLIxgBAAAAsDyCEQAAAADLIxgBAAAAsDyCEQAAAADLIxgBAAAAsDyCEQAAAADLIxgBAAAAsDx7uAsAAABnd9xdI0navLfclPaqjru0/qDUftdRJSU4TGlzR2mlKe0AQLgRjAAAaKaK/yd0PPHJJhNbtevDHetMbO+UJAdvKQBENo5iAAA0UyN6tJckdWmXrIS42KDb276/XJPmbdK0m3upW4fUoNurleSwK6ttkmntAUA4EIwAAGimWifF67b+HU1rz+PxSJK6pCep54XmBSMAiAZcfAEAAACA5QUVjF588UXZbDZNnDjRpHIAAGh69GcAgEYHo3Xr1untt99W7969zawHAIAmRX8GAJAaGYwqKyt155136t1331WrVq3MrgkAgCZBfwYAqNWoiy9MmDBBo0eP1rBhw/Tcc8+dc16XyyWXy+V7XFFRIUlyu91yu92NWXxEqV1HK6yrVbGPrcFq+9kq62m1/qz24gsejydiao5UbOumw7Zu2LHjp45ZX+8+4ttOwag6cepeaG13HlRSC5PuhXawSlJo95+/7QYcjObMmaMNGzZo3Tr/7oEwdepUTZkypd70JUuWKDExMdDFR6z8/Pxwl4AQYx9bg1X2c3V1dbhLCDkr9md7KiXJrjVr1mjv5nBXE93Y1k2Hbd2w1T/YJMXq159vNbFVuz7c8ZWJ7Z2ybvVK7UowvVlJ/vdpAQWjPXv26JFHHlF+fr5atGjh13OefPJJPfroo77HFRUVyszM1IgRI5SSkhLI4iOS2+1Wfn6+hg8frri4uHCXgxBgH1uD1fZz7WhItLJqf/b17iPSpvUaMGCALuvYOtzlRDW2ddNhWzdsQNVJ9dpWqs7pSabcC+3bA+X65afb9PKPs3VpezPvhRari9uE7l5o/vZpAQWjwsJClZaWqm/fvr5pNTU1WrFihd544w25XC7Fxtbd6A6HQw5H/aG2uLg4S7y5qGW19bUi9rE1WGU/R/s6WrU/s9vtvt+RUnOkYls3HbZ1wy5Ii9OduVmmt3tp+1T16dTG9HZDxd/XRUDB6JprrtGmTZvqTLvrrrvkdDr1q1/9ql4nAgBAc0R/BgA4U0DBqGXLlurZs2edaUlJSWrTpk296QAANFf0ZwCAMwV1g1cAAAAAiAaNulz36ZYtW2ZCGQAAhBf9GQBYGyNGAAAAACyPYAQAAADA8ghGAAAAACyPYAQAAADA8ghGAAAAACyPYAQAAADA8ghGAAAAACyPYAQAAADA8ghGAAAAACyPYAQAAADA8ghGAAAAACyPYAQAAADA8ghGAAAAACyPYAQAAADA8ghGAAAAACzPHu4CIll1dbWKiorOOU/lcZcKNhWrVdv1Sk5wnHNep9OpxMREM0sEAAAm8qfvl6Tt+8vkOrBD2zYnyHs4za+2eR8AhBfBKAhFRUXKycnxa96X/ZinsLBQffv2Da4oAAAQMoH0/ZJ0xwf+t837ACC8CEZBcDqdKiwsPOc82/eX6dGPN+m1W3qpW4e087YHAACaL3/6funUJ0YWfLFao4fmnvcTI6e3DSB8CEZBSExMPO+ZnZhdh+X413Fl97xMfTq1aaLKAABAKPjT90uS2+3W0UOlyu3fT3FxcU1QGYBgcfEFAAAAAJZHMAIAAABgeQQjAAAAAJZHMAIAAABgeVx84RxKDlWpyuUJqo3ig1W+33Z7cJs7yWFXVtukoNoAAAAAUB/BqAElh6o09NVlprU3ad4mU9r54rEhhCMAAADAZASjBtSOFE0f20dd2yU3vp3jLs1ftlpjhuQqyc/7GJzNjtJKTZy7MegRLAAAAAD1EYzOo2u7ZPW8MLXRz3e73TqQLvXt1Ir7GAAAAPipurpaRUVF551v+/4yuQ7s0LbNCfIeTjvv/E6nU4mJiSZUiGhDMAIAAECzU1RUpJycHL/nv+MD/+YrLCz06ya9sB6CEQAAAJodp9OpwsLC885XedylBV+s1uihuUr242sLTqfTjPIQhQhGAAAAaHYSExP9Gtlxu906eqhUuf378bUFBIX7GAEAAACwPIIRAAAAAMsjGAEAAACwPIIRAAAAAMsjGAEAAACwPIIRAAAAAMsjGAEAAACwPIIRAAAAAMsjGAEAAACwPIIRAAAAAMsjGAEAAACwPIIRAAAAAMsjGAEAAACwPHu4C2jObPYKlVRsV0yL5Ea34fF4tM+zT9uObJPd3vjNXVJRKZu9otHPBwAAANAwgtE5xKWt1VNfvmBKW28teivoNuLSrpF0XfDFAAAAAKiDYHQO7rIrNG30HerSrnEjRoZh6Kl/PaWdFTvVOaWzXrjqBdlstka1VVxaqYf/VNyo5wIAAAA4t4CC0YwZMzRjxgx9//33kqQePXro6aef1qhRo0JRW9gZnhRlpXRT9zapjXr+qr2rtLNipyRpZ8VOlbnKNPDCgY1qy3uiXIbnYKOeCwCoz2p9GgDg3AK6+MJFF12kF198UYWFhVq/fr1+9KMf6YYbbtCWLVtCVV/EMgxDeV/lKcZ2ahPH2GKU91WeDMMIc2UAAIk+DQBQV0DB6Prrr9d1112nSy65RJdeeqmef/55JScna82aNaGqL2IV7CvQlsNb5DW8kiSv4dWWw1tUsK8gzJUBACT6NABAXY3+jlFNTY0+/vhjVVVVKTc3t8H5XC6XXC6X73FFxakrq7ndbrnd7sYuPuQ8Ho/vd6B1Goah32/4vWJsMb5gJJ0aNfr9ht/r8vTLA/6uUTD1ILRq9wf7JbpZbT9bZT1r+dOnRWp/djr6kqZjtWNGOLGtm06kHkP8rTXgYLRp0ybl5ubqxIkTSk5O1qeffqru3bs3OP/UqVM1ZcqUetOXLFmixMTEQBffZPZUSpJdK1eu1K4Ar73wnfs7ba3aWm+61/Bq65Gt+v3ffq9L4i5psnrQNPLz88NdApqAVfZzdXV1uEtoEoH0aZHan52uti9Zs2aN9m4OdzXWYJVjRnPAtg69SD2G+Nun2YwAv/Ry8uRJ7d69W+Xl5Zo3b57ee+89LV++vMGO5Gxn2DIzM3Xo0CGlpKQEsugmtWVfhW6csUaf/XyAemT4X6dhGPrPxf+pbUe2yVD9TWuTTdmts/XhtR8GNGrU2HoQem63W/n5+Ro+fLji4uLCXQ5CxGr7uaKiQm3btlV5eXmzPlYHK5A+LVL7s9N9vfuIbn53vebd20+XdWwd7nKimtWOGeHEtm46kXoM8bdPC3jEKD4+Xl27dpUk5eTkaN26dXr99df19ttvn3V+h8Mhh8NRb3pcXFyzfvHW3ozVbrcHVOfJmpP6ofqHs4YiSTJk6IfqH6RYKS7W/3YbWw+aTnN/TcMcVtnPVlhHKbA+LVL7s9PRlzS9SHp9RDq2dehF6jHE31qDvo+R1+utcwbN6uJj4zVnzBwdOXFE0qnPYK5auUoDBw30vZhat2it+Nj4cJYJADgL+jQAsK6AgtGTTz6pUaNGqWPHjjp27Jhmz56tZcuWafHixaGqLyK1T2qv9kntJZ0a3i2xlyi7dXZEJWsAiHb0aQCA0wUUjEpLS/XTn/5U+/fvV2pqqnr37q3Fixdr+PDhoaoPAICQoE8DAJwuoGD0/vvvh6oOAACaFH0aAOB0Ad3gFQAAAACiEcEIAAAAgOURjAAAAABYXtCX6wYAAOFVXV2toqKi8863fX+ZXAd2aNvmBHkPp/nVttPpVGJiYpAVAkDzRzACACDCFRUVKScnx+/57/jA/7YLCwvVt2/fRlQFAJGFYAQAQIRzOp0qLCw873yVx11a8MVqjR6aq+QEh99tA4AVEIwAAIhwiYmJfo3quN1uHT1Uqtz+/bjpOACcgYsvAAAAALA8ghEAAAAAyyMYAQAAALA8vmPUgOPuGknS5r3lQbVTddyl9Qel9ruOKsnPL7qezY7SyqDqAAAAANAwglEDiv8niDzxySYTWrPrwx3rTGhHSnKwywAAAACz8S67ASN6tJckdWmXrIS42Ea3s31/uSbN26RpN/dStw6pQdWU5LArq21SUG0AAAAAqI9g1IDWSfG6rX/HoNvxeDySpC7pSep5YXDBCAAAAEBocPEFAAAAAJZHMAIAAABgeQQjAAAAAJZHMAIAAABgeQQjAAAAAJZHMAIAAABgeQQjAAAAAJZHMAIAAABgeQQjAAAAAJZHMAIAAABgeQQjAAAAAJZHMAIAAABgeQQjAAAAAJZHMAIAAABgeQQjAAAAAJZHMAIAAABgeQQjAAAAAJZHMAIAAABgefZwFxDJqqurVVRUdM55tu8vk+vADm3bnCDv4bRzzut0OpWYmGhihQAAAMC5+fOeVgrsfa0Uee9tCUZBKCoqUk5Ojl/z3vHB+ecpLCxU3759g6wKAAAA8F8g72kl/97XSpH33pZgFASn06nCwsJzzlN53KUFX6zW6KG5Sk5wnLc9AAAAoCn5855WCux9bW27kYRgFITExMTzpmC3262jh0qV27+f4uLimqgyAAAAwD/+vKeVov99LRdfAAAAAGB5BCMAAAAAlkcwAgAAAGB5BCMAAAAAlkcwAgAAAGB5BCMAAAAAlkcwAgAAAGB5BCMAAAAAlkcwAgAAAGB5AQWjqVOn6vLLL1fLli3Vrl073Xjjjdq+fXuoagMAICTozwAAZwooGC1fvlwTJkzQmjVrlJ+fL7fbrREjRqiqqipU9QEAYDr6MwDAmeyBzLxo0aI6j2fNmqV27dqpsLBQV199tamFAQAQKvRnAIAzBfUdo/LycklS69atTSkGAIBwoD8DAAQ0YnQ6r9eriRMnauDAgerZs2eD87lcLrlcLt/j2s7nyJEjcrvdjV18xHC73aqurtbhw4cVFxcX7nIQAuxja7Dafj527JgkyTCMMFcSelbqz6z2Og4ntnXTYVs3nUjd1n73aUYj3X///UanTp2MPXv2nHO+Z555xpDEDz/88MNPBP6c7xgfDejP+OGHH36s8XO+47zNMAI/Hfjggw/q888/14oVK5SVlXXOec88w+b1enXkyBG1adNGNpst0EVHnIqKCmVmZmrPnj1KSUkJdzkIAfaxNVhtPxuGoWPHjikjI0MxMdF7Zwer9WdWex2HE9u66bCtm06kbmt/+7SAPkpnGIYeeughffrpp1q2bNl5OxFJcjgccjgcdaalpaUFstiokJKSElEvIASOfWwNVtrPqamp4S4hZKzen1npdRxubOumw7ZuOpG4rf3p0wIKRhMmTNDs2bP1+eefq2XLljpw4IBvQQkJCY2rEgCAJkZ/BgA4U0Cfj5gxY4bKy8s1ZMgQdejQwfczd+7cUNUHAIDp6M8AAGcK+KN0CIzD4dAzzzxT7+MXiB7sY2tgP0cXq/ZnvI6bDtu66bCtm060b+tGXXwBAAAAAKJJ9F5qCAAAAAD8RDACAAAAYHkEIwAAAACWRzACAAAAYHkEoxBbvXq1YmNjNXr06HCXApONHz9eNpvN99OmTRuNHDlS33zzTbhLg8kOHDighx56SJ07d5bD4VBmZqauv/56LV26NNylAX5bsWKFrr/+emVkZMhms+mzzz4Ld0lRaerUqbr88svVsmVLtWvXTjfeeKO2b98e7rKi0owZM9S7d2/fzUZzc3O1cOHCcJcV9V588UXZbDZNnDgx3KWYjmAUYu+//74eeughrVixQvv27Qt3OTDZyJEjtX//fu3fv19Lly6V3W7XmDFjwl0WTPT9998rJydH//3f/61XXnlFmzZt0qJFizR06FBNmDAh3OUBfquqqtJll12mN998M9ylRLXly5drwoQJWrNmjfLz8+V2uzVixAhVVVWFu7Soc9FFF+nFF19UYWGh1q9frx/96Ee64YYbtGXLlnCXFrXWrVunt99+W7179w53KSHB5bpDqLKyUh06dND69ev1zDPPqHfv3nrqqafCXRZMMn78eJWVldU567py5UpdddVVKi0tVXp6eviKg2muu+46ffPNN9q+fbuSkpLq/K2srExpaWnhKQwIgs1m06effqobb7wx3KVEvYMHD6pdu3Zavny5rr766nCXE/Vat26tV155Rffcc0+4S4k6lZWV6tu3r9566y0999xz6tOnj6ZPnx7uskzFiFEI/eUvf5HT6VS3bt30k5/8RH/4wx8se1NBK6isrNRHH32krl27qk2bNuEuByY4cuSIFi1apAkTJtQLRZIIRQDOq7y8XNKpN+wInZqaGs2ZM0dVVVXKzc0NdzlRacKECRo9erSGDRsW7lJCxh7uAqLZ+++/r5/85CeSTn3kqry8XMuXL9eQIUPCWxhMM3/+fCUnJ0s69TGVDh06aP78+YqJ4ZxDNNixY4cMw5DT6Qx3KQAikNfr1cSJEzVw4ED17Nkz3OVEpU2bNik3N1cnTpxQcnKyPv30U3Xv3j3cZUWdOXPmaMOGDVq3bl24Swkp3r2FyPbt2/Xll1/q9ttvlyTZ7XaNHTtW77//fpgrg5mGDh2qjRs3auPGjfryyy917bXXatSoUdq1a1e4S4MJGOEFEIwJEyZo8+bNmjNnTrhLiVrdunXTxo0btXbtWv385z/XuHHjtHXr1nCXFVX27NmjRx55RH/605/UokWLcJcTUowYhcj7778vj8ejjIwM3zTDMORwOPTGG28oNTU1jNXBLElJSeratavv8XvvvafU1FS9++67eu6558JYGcxwySWXyGazqaioKNylAIgwDz74oObPn68VK1booosuCnc5USs+Pt7XD+fk5GjdunV6/fXX9fbbb4e5suhRWFio0tJS9e3b1zetpqZGK1as0BtvvCGXy6XY2NgwVmgeRoxCwOPx6I9//KOmTZvmG03YuHGjvv76a2VkZOjPf/5zuEtEiNhsNsXExOj48ePhLgUmaN26ta699lq9+eabZ72iVFlZWdMXBaBZMwxDDz74oD799FP993//t7KyssJdkqV4vV65XK5wlxFVrrnmGm3atKnOe9p+/frpzjvv1MaNG6MmFEmMGIXE/PnzdfToUd1zzz31RoZuuukmvf/++7r//vvDVB3M5HK5dODAAUnS0aNH9cYbb6iyslLXX399mCuDWd58800NHDhQ/fv317PPPqvevXvL4/EoPz9fM2bM0LZt28JdIuCXyspK7dixw/e4pKREGzduVOvWrdWxY8cwVhZdJkyYoNmzZ+vzzz9Xy5YtfX1EamqqEhISwlxddHnyySc1atQodezYUceOHdPs2bO1bNkyLV68ONylRZWWLVvW+45cUlKS2rRpE3XfnSMYhcD777+vYcOGnfXjcjfddJNefvllffPNN1F7DXgrWbRokTp06CDp1IHD6XTq448/5gIbUaRz587asGGDnn/+eU2aNEn79+9Xenq6cnJyNGPGjHCXB/ht/fr1Gjp0qO/xo48+KkkaN26cZs2aFaaqok/tceHMfmDmzJkaP3580xcUxUpLS/XTn/5U+/fvV2pqqnr37q3Fixdr+PDh4S4NEYr7GAEAAACwPL5jBAAAAMDyCEYAAAAALI9gBAAAAMDyCEYAAAAALI9gBAAAAMDyCEYAAAAALI9gBAAAAMDyCEYAAABRYsiQIZo4cWKDf7/44os1ffp005ZndntAOBGMEPHO1wmYafz48brxxhubZFn+mDVrltLS0sJdBgAAQMQjGCHqGYYhj8cT7jIAAADQjBGMENHGjx+v5cuX6/XXX5fNZpPNZtOsWbNks9m0cOFC5eTkyOFwaOXKlfJ6vZo6daqysrKUkJCgyy67TPPmzfO1VVNTo3vuucf3927duun111/3/X3y5Mn64IMP9Pnnn/uWtWzZMn3//fey2Wz6y1/+oquuukoJCQm6/PLL9e2332rdunXq16+fkpOTNWrUKB08eLBO/e+9956ys7PVokULOZ1OvfXWW76/1bb7ySefaOjQoUpMTNRll12m1atXS5KWLVumu+66S+Xl5b56Jk+eHNoNDgBo9jwejx588EGlpqaqbdu2+s1vfiPDMM467+7du3XDDTcoOTlZKSkpuvXWW/XDDz/Umefvf/+7Lr/8crVo0UJt27bVj3/84waX/d577yktLU1Lly41dZ2AJmEAEaysrMzIzc017r33XmP//v3G/v37jX/+85+GJKN3797GkiVLjB07dhiHDx82nnvuOcPpdBqLFi0yiouLjZkzZxoOh8NYtmyZYRiGcfLkSePpp5821q1bZ+zcudP46KOPjMTERGPu3LmGYRjGsWPHjFtvvdUYOXKkb1kul8soKSkxJPna3rp1qzFgwAAjJyfHGDJkiLFy5Upjw4YNRteuXY3777/fV/tHH31kdOjQwfjrX/9q7Ny50/jrX/9qtG7d2pg1a5ZhGEaddufPn29s377duPnmm41OnToZbrfbcLlcxvTp042UlBRfPceOHWv6nQAAaDYGDx5sJCcnG4888ohRVFTk68veeecdwzAMo1OnTsbvfvc7wzAMo6amxujTp48xaNAgY/369caaNWuMnJwcY/Dgwb725s+fb8TGxhpPP/20sXXrVmPjxo3GCy+84Pv76e299NJLRps2bYy1a9c21eoCpiIYIeINHjzYeOSRR3yPv/jiC0OS8dlnn/mmnThxwkhMTDQKCgrqPPeee+4xbr/99gbbnjBhgnHTTTf5Ho8bN8644YYb6sxTG2Dee+8937Q///nPhiRj6dKlvmlTp041unXr5nvcpUsXY/bs2XXa+u1vf2vk5uY22O6WLVsMSca2bdsMwzCMmTNnGqmpqQ3WDwCwlsGDBxvZ2dmG1+v1TfvVr35lZGdnG4ZRN8gsWbLEiI2NNXbv3u2bt7af+fLLLw3DMIzc3FzjzjvvbHB5te398pe/NDp06GBs3rw5BGsFNA172IaqgBDr16+f7987duxQdXW1hg8fXmeekydP6t/+7d98j99880394Q9/0O7du3X8+HGdPHlSffr08Wt5vXv39v37ggsukCT16tWrzrTS0lJJUlVVlYqLi3XPPffo3nvv9c3j8XiUmpraYLsdOnSQJJWWlsrpdPpVFwDAWgYMGCCbzeZ7nJubq2nTpqmmpqbOfNu2bVNmZqYyMzN907p37660tDRt27ZNl19+uTZu3FinnzqbadOmqaqqSuvXr1fnzp3NXRmgCRGMELWSkpJ8/66srJQkLViwQBdeeGGd+RwOhyRpzpw5euyxxzRt2jTl5uaqZcuWeuWVV7R27Vq/lhcXF+f7d22HdOY0r9dbp553331XV1xxRZ12YmNjz9tubTsAAIRSQkLCeee56qqrtGDBAv3lL3/RE0880QRVAaFBMELEi4+Pr3cW7Ezdu3eXw+HQ7t27NXjw4LPOs2rVKl155ZV64IEHfNOKi4sDXpY/LrjgAmVkZGjnzp268847G92OWfUAAKLHmSf01qxZo0suuaTeibfs7Gzt2bNHe/bs8Y0abd26VWVlZerevbukU59aWLp0qe66664Gl9e/f389+OCDGjlypOx2ux577DGT1whoGgQjRLyLL75Ya9eu1ffff6/k5OSzjqa0bNlSjz32mH7xi1/I6/Vq0KBBKi8v16pVq5SSkqJx48bpkksu0R//+EctXrxYWVlZ+vDDD7Vu3TplZWXVWdbixYu1fft2tWnTpt7H3gIxZcoUPfzww0pNTdXIkSPlcrm0fv16HT16VI8++qjf615ZWamlS5fqsssuU2JiohITExtdEwAg8u3evVuPPvqo7rvvPm3YsEF5eXmaNm1avfmGDRumXr166c4779T06dPl8Xj0wAMPaPDgwb6Poz/zzDO65ppr1KVLF912223yeDz6xz/+oV/96ld12rryyiv1j3/8Q6NGjZLdbm+y+wsCZuJy3Yh4jz32mGJjY9W9e3elp6dr9+7dZ53vt7/9rX7zm99o6tSpys7O1siRI7VgwQJf8Lnvvvv0H//xHxo7dqyuuOIKHT58uM7okSTde++96tatm/r166f09HStWrWq0XX/7Gc/03vvvaeZM2eqV69eGjx4sGbNmlUniJ3PlVdeqfvvv19jx45Venq6Xn755UbXAwCIDj/96U91/Phx9e/fXxMmTNAjjzyi//N//k+9+Ww2mz7//HO1atVKV199tYYNG6bOnTtr7ty5vnmGDBmijz/+WH/729/Up08f/ehHP9KXX3551uUOGjRICxYs0H/9138pLy8vZOsHhIrNMBq4sD0AAAAAWAQjRgAAAAAsj2AEAAAAwPIIRgAAAAAsj2AEAAAAwPIIRgAAAAAsj2AEAAAAwPIIRgAAAAAsj2AEAAAAwPIIRgAAAAAsj2AEAAAAwPIIRgAAAAAsj2AEAAAAwPL+P1AMOsV2KbYwAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, axs = plt.subplots(1, 2, figsize=(10, 5))\n", "data.boxplot(\"value\", by=\"treatment\", ax=axs[0], showmeans=True)\n", "data.boxplot(\"value\", by=\"block\", ax=axs[1], showmeans=True)\n", "fig.suptitle(\"Boxplots\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Things to consider:\n", "\n", "Does it look like the value is different for different treatments?\n", "\n", "Does it look like the value is diiferent for different blocks?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Estimate parameters $\\mu$, $\\alpha_i$ and $\\beta_j$" ] }, { "cell_type": "code", "execution_count": 92, "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", " \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", " \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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
valuetreatmentblockoverall_mean
02.8A15.233333
13.6A25.233333
23.4A35.233333
32.3A45.233333
45.5B15.233333
56.3B25.233333
66.1B35.233333
75.7B45.233333
85.8C15.233333
98.3C25.233333
106.9C35.233333
116.1C45.233333
\n", "
" ], "text/plain": [ " value treatment block overall_mean\n", "0 2.8 A 1 5.233333\n", "1 3.6 A 2 5.233333\n", "2 3.4 A 3 5.233333\n", "3 2.3 A 4 5.233333\n", "4 5.5 B 1 5.233333\n", "5 6.3 B 2 5.233333\n", "6 6.1 B 3 5.233333\n", "7 5.7 B 4 5.233333\n", "8 5.8 C 1 5.233333\n", "9 8.3 C 2 5.233333\n", "10 6.9 C 3 5.233333\n", "11 6.1 C 4 5.233333" ] }, "execution_count": 92, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Compute the overall mean and add to dataframe:\n", "data['overall_mean'] = data[\"value\"].mean()\n", "data" ] }, { "cell_type": "code", "execution_count": 93, "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", " \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", " \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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
valuetreatmentblockoverall_meantreatment_mean
02.8A15.2333333.025
13.6A25.2333333.025
23.4A35.2333333.025
32.3A45.2333333.025
45.5B15.2333335.900
56.3B25.2333335.900
66.1B35.2333335.900
75.7B45.2333335.900
85.8C15.2333336.775
98.3C25.2333336.775
106.9C35.2333336.775
116.1C45.2333336.775
\n", "
" ], "text/plain": [ " value treatment block overall_mean treatment_mean\n", "0 2.8 A 1 5.233333 3.025\n", "1 3.6 A 2 5.233333 3.025\n", "2 3.4 A 3 5.233333 3.025\n", "3 2.3 A 4 5.233333 3.025\n", "4 5.5 B 1 5.233333 5.900\n", "5 6.3 B 2 5.233333 5.900\n", "6 6.1 B 3 5.233333 5.900\n", "7 5.7 B 4 5.233333 5.900\n", "8 5.8 C 1 5.233333 6.775\n", "9 8.3 C 2 5.233333 6.775\n", "10 6.9 C 3 5.233333 6.775\n", "11 6.1 C 4 5.233333 6.775" ] }, "execution_count": 93, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# compute the mean within each treatment-group and add to dataframe:\n", "data['treatment_mean'] = data.groupby(\"treatment\")['value'].transform('mean')\n", "data" ] }, { "cell_type": "code", "execution_count": 94, "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", " \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", " \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", " \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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
valuetreatmentblockoverall_meantreatment_meantreatment_alpha
02.8A15.2333333.025-2.208333
13.6A25.2333333.025-2.208333
23.4A35.2333333.025-2.208333
32.3A45.2333333.025-2.208333
45.5B15.2333335.9000.666667
56.3B25.2333335.9000.666667
66.1B35.2333335.9000.666667
75.7B45.2333335.9000.666667
85.8C15.2333336.7751.541667
98.3C25.2333336.7751.541667
106.9C35.2333336.7751.541667
116.1C45.2333336.7751.541667
\n", "
" ], "text/plain": [ " value treatment block overall_mean treatment_mean treatment_alpha\n", "0 2.8 A 1 5.233333 3.025 -2.208333\n", "1 3.6 A 2 5.233333 3.025 -2.208333\n", "2 3.4 A 3 5.233333 3.025 -2.208333\n", "3 2.3 A 4 5.233333 3.025 -2.208333\n", "4 5.5 B 1 5.233333 5.900 0.666667\n", "5 6.3 B 2 5.233333 5.900 0.666667\n", "6 6.1 B 3 5.233333 5.900 0.666667\n", "7 5.7 B 4 5.233333 5.900 0.666667\n", "8 5.8 C 1 5.233333 6.775 1.541667\n", "9 8.3 C 2 5.233333 6.775 1.541667\n", "10 6.9 C 3 5.233333 6.775 1.541667\n", "11 6.1 C 4 5.233333 6.775 1.541667" ] }, "execution_count": 94, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# also add the corresponding alpha for each treatment-group:\n", "data[\"treatment_alpha\"] = data[\"treatment_mean\"] - data[\"overall_mean\"]\n", "data" ] }, { "cell_type": "code", "execution_count": 95, "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", " \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", " \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", " \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", " \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", " \n", " \n", " \n", " \n", " \n", "
valuetreatmentblockoverall_meantreatment_meantreatment_alphablock_meanblock_beta
02.8A15.2333333.025-2.2083334.700000-0.533333
13.6A25.2333333.025-2.2083336.0666670.833333
23.4A35.2333333.025-2.2083335.4666670.233333
32.3A45.2333333.025-2.2083334.700000-0.533333
45.5B15.2333335.9000.6666674.700000-0.533333
56.3B25.2333335.9000.6666676.0666670.833333
66.1B35.2333335.9000.6666675.4666670.233333
75.7B45.2333335.9000.6666674.700000-0.533333
85.8C15.2333336.7751.5416674.700000-0.533333
98.3C25.2333336.7751.5416676.0666670.833333
106.9C35.2333336.7751.5416675.4666670.233333
116.1C45.2333336.7751.5416674.700000-0.533333
\n", "
" ], "text/plain": [ " value treatment block overall_mean treatment_mean treatment_alpha \\\n", "0 2.8 A 1 5.233333 3.025 -2.208333 \n", "1 3.6 A 2 5.233333 3.025 -2.208333 \n", "2 3.4 A 3 5.233333 3.025 -2.208333 \n", "3 2.3 A 4 5.233333 3.025 -2.208333 \n", "4 5.5 B 1 5.233333 5.900 0.666667 \n", "5 6.3 B 2 5.233333 5.900 0.666667 \n", "6 6.1 B 3 5.233333 5.900 0.666667 \n", "7 5.7 B 4 5.233333 5.900 0.666667 \n", "8 5.8 C 1 5.233333 6.775 1.541667 \n", "9 8.3 C 2 5.233333 6.775 1.541667 \n", "10 6.9 C 3 5.233333 6.775 1.541667 \n", "11 6.1 C 4 5.233333 6.775 1.541667 \n", "\n", " block_mean block_beta \n", "0 4.700000 -0.533333 \n", "1 6.066667 0.833333 \n", "2 5.466667 0.233333 \n", "3 4.700000 -0.533333 \n", "4 4.700000 -0.533333 \n", "5 6.066667 0.833333 \n", "6 5.466667 0.233333 \n", "7 4.700000 -0.533333 \n", "8 4.700000 -0.533333 \n", "9 6.066667 0.833333 \n", "10 5.466667 0.233333 \n", "11 4.700000 -0.533333 " ] }, "execution_count": 95, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# compute the mean within each person (block) and add to dataframe:\n", "data['block_mean'] = data.groupby(\"block\")['value'].transform('mean')\n", "\n", "# also add the corresponding beta for each person (block):\n", "data[\"block_beta\"] = data[\"block_mean\"] - data[\"overall_mean\"]\n", "\n", "data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: SST, SS(Tr), SS(Bl) and SSE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$Y_{ij} = \\mu + \\alpha_i + \\beta_j + \\epsilon_{ij}$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We start with residuals (and SSE):" ] }, { "cell_type": "code", "execution_count": 96, "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", " \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", " \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", " \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", " \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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
valuetreatmentblockoverall_meantreatment_meantreatment_alphablock_meanblock_betaresidual
02.8A15.2333333.025-2.2083334.700000-0.5333330.308333
13.6A25.2333333.025-2.2083336.0666670.833333-0.258333
23.4A35.2333333.025-2.2083335.4666670.2333330.141667
32.3A45.2333333.025-2.2083334.700000-0.533333-0.191667
45.5B15.2333335.9000.6666674.700000-0.5333330.133333
56.3B25.2333335.9000.6666676.0666670.833333-0.433333
66.1B35.2333335.9000.6666675.4666670.233333-0.033333
75.7B45.2333335.9000.6666674.700000-0.5333330.333333
85.8C15.2333336.7751.5416674.700000-0.533333-0.441667
98.3C25.2333336.7751.5416676.0666670.8333330.691667
106.9C35.2333336.7751.5416675.4666670.233333-0.108333
116.1C45.2333336.7751.5416674.700000-0.533333-0.141667
\n", "
" ], "text/plain": [ " value treatment block overall_mean treatment_mean treatment_alpha \\\n", "0 2.8 A 1 5.233333 3.025 -2.208333 \n", "1 3.6 A 2 5.233333 3.025 -2.208333 \n", "2 3.4 A 3 5.233333 3.025 -2.208333 \n", "3 2.3 A 4 5.233333 3.025 -2.208333 \n", "4 5.5 B 1 5.233333 5.900 0.666667 \n", "5 6.3 B 2 5.233333 5.900 0.666667 \n", "6 6.1 B 3 5.233333 5.900 0.666667 \n", "7 5.7 B 4 5.233333 5.900 0.666667 \n", "8 5.8 C 1 5.233333 6.775 1.541667 \n", "9 8.3 C 2 5.233333 6.775 1.541667 \n", "10 6.9 C 3 5.233333 6.775 1.541667 \n", "11 6.1 C 4 5.233333 6.775 1.541667 \n", "\n", " block_mean block_beta residual \n", "0 4.700000 -0.533333 0.308333 \n", "1 6.066667 0.833333 -0.258333 \n", "2 5.466667 0.233333 0.141667 \n", "3 4.700000 -0.533333 -0.191667 \n", "4 4.700000 -0.533333 0.133333 \n", "5 6.066667 0.833333 -0.433333 \n", "6 5.466667 0.233333 -0.033333 \n", "7 4.700000 -0.533333 0.333333 \n", "8 4.700000 -0.533333 -0.441667 \n", "9 6.066667 0.833333 0.691667 \n", "10 5.466667 0.233333 -0.108333 \n", "11 4.700000 -0.533333 -0.141667 " ] }, "execution_count": 96, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# calculate the residuals according to the model:\n", "data['residual'] = data['value'] - (data['overall_mean']+data['treatment_alpha']+data['block_beta'])\n", "data" ] }, { "cell_type": "code", "execution_count": 97, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.241666666666668\n" ] } ], "source": [ "data['sse_contribution'] = data['residual']**2\n", "\n", "SSE = data['sse_contribution'].sum()\n", "\n", "print(SSE)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now calculate SST:" ] }, { "cell_type": "code", "execution_count": 98, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "35.98666666666667\n" ] } ], "source": [ "data['sst_contribution'] = (data['value'] - data['overall_mean'])**2\n", "\n", "SST = data['sst_contribution'].sum()\n", "\n", "print(SST)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now calculate SS(tr) and SS(Bl):" ] }, { "cell_type": "code", "execution_count": 99, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "30.791666666666668\n" ] } ], "source": [ "data['sstr_contribution'] = (data['treatment_mean'] - data['overall_mean'])**2\n", "\n", "SSTr = data['sstr_contribution'].sum()\n", "\n", "print(SSTr)" ] }, { "cell_type": "code", "execution_count": 100, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.953333333333335\n" ] } ], "source": [ "data['ssbl_contribution'] = (data['block_mean'] - data['overall_mean'])**2\n", "\n", "SSBl = data['ssbl_contribution'].sum()\n", "\n", "print(SSBl)" ] }, { "cell_type": "code", "execution_count": 101, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "35.98666666666667\n", "35.98666666666667\n" ] } ], "source": [ "# sanity check:\n", "print(SST)\n", "print(SSTr + SSBl + SSE)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Projection set up" ] }, { "cell_type": "code", "execution_count": 102, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[1., 0., 0., 0., 0., 0.],\n", " [1., 0., 0., 1., 0., 0.],\n", " [1., 0., 0., 0., 1., 0.],\n", " [1., 0., 0., 0., 0., 1.],\n", " [1., 1., 0., 0., 0., 0.],\n", " [1., 1., 0., 1., 0., 0.],\n", " [1., 1., 0., 0., 1., 0.],\n", " [1., 1., 0., 0., 0., 1.],\n", " [1., 0., 1., 0., 0., 0.],\n", " [1., 0., 1., 1., 0., 0.],\n", " [1., 0., 1., 0., 1., 0.],\n", " [1., 0., 1., 0., 0., 1.]])" ] }, "execution_count": 102, "metadata": {}, "output_type": "execute_result" } ], "source": [ "one = np.repeat(1,4)\n", "zero = np.repeat(0,4)\n", "X0 = np.ones((12,1))\n", "col2 = np.append(np.append(zero,one),zero)\n", "col3 = np.append(np.append(zero,zero),one)\n", "\n", "X_tr = np.array([col2,col3]).T\n", "\n", "X_bl0 = np.c_[X0,X_tr]\n", "\n", "X_bl = np.array([[0,1,0,0,0,1,0,0,0,1,0,0],\n", " [0,0,1,0,0,0,1,0,0,0,1,0],\n", " [0,0,0,1,0,0,0,1,0,0,0,1]]).T\n", "\n", "X = np.c_[X0,X_tr,X_bl]\n", "X" ] }, { "cell_type": "code", "execution_count": 103, "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", " \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", " \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", " \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", " \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", " \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", " \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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
valuetreatmentblockoverall_meantreatment_meantreatment_alphablock_meanblock_betaresidualsse_contributionsst_contributionsstr_contributionssbl_contribution
02.8A15.2333333.025-2.2083334.700000-0.5333330.3083330.0950695.9211114.8767360.284444
13.6A25.2333333.025-2.2083336.0666670.833333-0.2583330.0667362.6677784.8767360.694444
23.4A35.2333333.025-2.2083335.4666670.2333330.1416670.0200693.3611114.8767360.054444
32.3A45.2333333.025-2.2083334.700000-0.533333-0.1916670.0367368.6044444.8767360.284444
45.5B15.2333335.9000.6666674.700000-0.5333330.1333330.0177780.0711110.4444440.284444
56.3B25.2333335.9000.6666676.0666670.833333-0.4333330.1877781.1377780.4444440.694444
66.1B35.2333335.9000.6666675.4666670.233333-0.0333330.0011110.7511110.4444440.054444
75.7B45.2333335.9000.6666674.700000-0.5333330.3333330.1111110.2177780.4444440.284444
85.8C15.2333336.7751.5416674.700000-0.533333-0.4416670.1950690.3211112.3767360.284444
98.3C25.2333336.7751.5416676.0666670.8333330.6916670.4784039.4044442.3767360.694444
106.9C35.2333336.7751.5416675.4666670.233333-0.1083330.0117362.7777782.3767360.054444
116.1C45.2333336.7751.5416674.700000-0.533333-0.1416670.0200690.7511112.3767360.284444
\n", "
" ], "text/plain": [ " value treatment block overall_mean treatment_mean treatment_alpha \\\n", "0 2.8 A 1 5.233333 3.025 -2.208333 \n", "1 3.6 A 2 5.233333 3.025 -2.208333 \n", "2 3.4 A 3 5.233333 3.025 -2.208333 \n", "3 2.3 A 4 5.233333 3.025 -2.208333 \n", "4 5.5 B 1 5.233333 5.900 0.666667 \n", "5 6.3 B 2 5.233333 5.900 0.666667 \n", "6 6.1 B 3 5.233333 5.900 0.666667 \n", "7 5.7 B 4 5.233333 5.900 0.666667 \n", "8 5.8 C 1 5.233333 6.775 1.541667 \n", "9 8.3 C 2 5.233333 6.775 1.541667 \n", "10 6.9 C 3 5.233333 6.775 1.541667 \n", "11 6.1 C 4 5.233333 6.775 1.541667 \n", "\n", " block_mean block_beta residual sse_contribution sst_contribution \\\n", "0 4.700000 -0.533333 0.308333 0.095069 5.921111 \n", "1 6.066667 0.833333 -0.258333 0.066736 2.667778 \n", "2 5.466667 0.233333 0.141667 0.020069 3.361111 \n", "3 4.700000 -0.533333 -0.191667 0.036736 8.604444 \n", "4 4.700000 -0.533333 0.133333 0.017778 0.071111 \n", "5 6.066667 0.833333 -0.433333 0.187778 1.137778 \n", "6 5.466667 0.233333 -0.033333 0.001111 0.751111 \n", "7 4.700000 -0.533333 0.333333 0.111111 0.217778 \n", "8 4.700000 -0.533333 -0.441667 0.195069 0.321111 \n", "9 6.066667 0.833333 0.691667 0.478403 9.404444 \n", "10 5.466667 0.233333 -0.108333 0.011736 2.777778 \n", "11 4.700000 -0.533333 -0.141667 0.020069 0.751111 \n", "\n", " sstr_contribution ssbl_contribution \n", "0 4.876736 0.284444 \n", "1 4.876736 0.694444 \n", "2 4.876736 0.054444 \n", "3 4.876736 0.284444 \n", "4 0.444444 0.284444 \n", "5 0.444444 0.694444 \n", "6 0.444444 0.054444 \n", "7 0.444444 0.284444 \n", "8 2.376736 0.284444 \n", "9 2.376736 0.694444 \n", "10 2.376736 0.054444 \n", "11 2.376736 0.284444 " ] }, "execution_count": 103, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data" ] }, { "cell_type": "code", "execution_count": 106, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([30.79166667, 3.95333333, 1.24166667, 35.98666667])" ] }, "execution_count": 106, "metadata": {}, "output_type": "execute_result" } ], "source": [ "H = X @ np.linalg.inv(X.T @ X) @ X.T\n", "H_bl0 = X_bl0 @ np.linalg.inv(X_bl0.T @ X_bl0) @ X_bl0.T\n", "H0 = X0@ np.linalg.inv(X0.T @ X0) @ X0.T\n", "I = np.identity(12)\n", "y = data[\"value\"]\n", "SST = y.T @ (I - H0) @ y\n", "SSTr = y.T @ (H_bl0 - H0) @ y\n", "SSBl = y.T @ (H - H_bl0) @ y\n", "SSE = y.T @ (I - H) @ y\n", "np.array([SSTr, SSBl, SSE, SST])\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: F-test and ANOVA table" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "TReatment groups:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$F_{Tr}=\\frac{SS(Tr)/(k-1)}{SSE/((k-1)(l-1))}$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Compute F-statistic (for treatment groups)\n", "FTr = (SSTr/(3-1)) / (SSE/((3-1)*(4-1)))\n", "print(FTr)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# compute critical value:\n", "stats.f.ppf(0.95, dfn = (3-1), dfd = (3-1)*(4-1) )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# p-value\n", "p_Tr = 1 - stats.f.cdf(FTr, dfn = (3-1), dfd = (3-1)*(4-1) )\n", "print(p_Tr)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Blocks:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$F_{Bl}=\\frac{SS(Bl)/(l-1)}{SSE/((k-1)(l-1))}$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Compute F-statistic (for block)\n", "FBl = (SSBl/(4-1)) / (SSE/((3-1)*(4-1)))\n", "print(FBl)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# compute critical value:\n", "stats.f.ppf(0.95, dfn = (4-1), dfd = (3-1)*(4-1) )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# p-value\n", "p_Bl = 1 - stats.f.cdf(FBl, dfn = (4-1), dfd = (3-1)*(4-1) )\n", "print(p_Bl)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Twoway ANOVA table" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Make the ANOVA table:\n", "fit = smf.ols(\"value ~ treatment + block\", data=data).fit()\n", "anova_table = sm.stats.anova_lm(fit)\n", "print(anova_table)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Model diagnostics" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$Y_{ij} = \\mu + \\alpha_i + \\beta_j + \\epsilon_{ij}$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# QQplot:\n", "sm.qqplot(data[\"residual\"], line='q',a=1/2)\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plot residuals grouped by either treatment or block:\n", "fig, axs = plt.subplots(1, 2, figsize=(10, 5))\n", "data.boxplot(\"residual\", by=\"treatment\", showmeans=True, ax=axs[0])\n", "data.boxplot(\"residual\", by=\"block\", showmeans=True, ax=axs[1])\n", "fig.suptitle(\"Boxplots of residuals\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Skive Fjord" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "SkiveAvg = pd.read_csv(\"../week1/skiveAvg.csv\", sep=';')\n", "SkiveAvg[\"vmp\"] = pd.Categorical(SkiveAvg[\"vmp\"])\n", "SkiveAvg[\"month\"] = pd.Categorical(SkiveAvg[\"month\"])\n", "SkiveAvg[\"year_f\"] = pd.Categorical(SkiveAvg[\"year\"])\n", "SkiveAvg" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fit = smf.ols(\"Nload ~ month + vmp\", data=SkiveAvg).fit()\n", "sm.stats.anova_lm(fit,typ=3)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sm.stats.anova_lm(fit,typ=3)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fit.summary(slim=True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#fig, ax =plt.subplots(1,2)\n", "plt.plot(np.arange(2,13),fit.params[1:12])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fit_year = smf.ols(\"Nload ~ month + year_f\", data=SkiveAvg).fit()\n", "sm.stats.anova_lm(fit_year,typ=1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sm.stats.anova_lm(fit,fit_year)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax =plt.subplots(1,2)\n", "ax[0].plot(np.arange(2,13),fit_year.params[1:12])\n", "ax[1].plot(np.arange(1983,2007),np.array(fit_year.params[12:36]))" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.16" } }, "nbformat": 4, "nbformat_minor": 2 }