{ "cells": [ { "cell_type": "markdown", "id": "e5a4bea8", "metadata": {}, "source": [ "# Permutation test for unpaired or independent samples data\n", "\n", "We now look at the case of unpaired data - data in which we wish to compare two groups but there are no particular links between individuals in one group and individuals in another group\n", "\n", "Examples of experimental designs using independent samples would be (say) a comparison of science and humanities students, or first- and final-year students, where the individuals in each group are distinct.\n", "\n", "### Set up Python libraries\n", "\n", "As usual, run the code cell below to import the relevant Python libraries" ] }, { "cell_type": "code", "execution_count": 1, "id": "83b773a6", "metadata": { "tags": [] }, "outputs": [], "source": [ "# Set-up Python libraries - you need to run this but you don't need to change it\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import scipy.stats as stats\n", "import pandas as pd\n", "import seaborn as sns\n", "sns.set_theme(style='white')\n", "import statsmodels.api as sm\n", "import statsmodels.formula.api as smf\n", "import warnings \n", "warnings.simplefilter('ignore', category=FutureWarning)" ] }, { "cell_type": "markdown", "id": "8d4519d3", "metadata": {}, "source": [ "## Toy example\n", "\n", "[A toy example is an example with a very small dataset, just to show how it works]\n", "\n", "\n", "#### Question & design\n", "\n", "We hypothesise that cat owners eat more yoghurt than dog owners \n", "\n", "We decide on an **independent samples** design in which a group of cat owners are compared to an unrelated group of dog owners. \n", "\n", "#### Hypotheses\n", "\n", "We can state our hypotheses as follows:\n", "\n", "$\\mathcal{H_o}:$ The difference in the mean number of yoghurts eaten by cat owners and dog owners is zero\n", "* mean(cat owners' yoghurt consumption) - mean(dog owners' yoghurt consumption) = 0\n", "\n", "$\\mathcal{H_a}:$ Cat owners eat more yoghurt\n", "* mean(cat owners' yoghurt consumption) > mean(dog owners' yoghurt consumption) = 0\n", "\n", "This is a one-tailed (directional) alternative hypothesis\n", "\n", "#### Data\n", "\n", "We survey five cat owners and four dog owners and ask how many times in the last week they ate yogurt. We obtain the following data:" ] }, { "cell_type": "code", "execution_count": 2, "id": "7a66358a", "metadata": { "tags": [] }, "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", "
PetYoghurt
0cat0
1cat7
2cat5
3cat2
4cat13
5cat1
6cat8
7cat6
8cat3
9cat5
10dog0
11dog0
12dog6
13dog4
14dog1
15dog4
16dog3
17dog2
\n", "
" ], "text/plain": [ " Pet Yoghurt\n", "0 cat 0\n", "1 cat 7\n", "2 cat 5\n", "3 cat 2\n", "4 cat 13\n", "5 cat 1\n", "6 cat 8\n", "7 cat 6\n", "8 cat 3\n", "9 cat 5\n", "10 dog 0\n", "11 dog 0\n", "12 dog 6\n", "13 dog 4\n", "14 dog 1\n", "15 dog 4\n", "16 dog 3\n", "17 dog 2" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "yog = pd.read_csv('https://raw.githubusercontent.com/jillxoreilly/StatsCourseBook_2024/main/data/yoghurt.csv')\n", "yog" ] }, { "cell_type": "markdown", "id": "d0c68564", "metadata": {}, "source": [ "Let's plot the data:" ] }, { "cell_type": "code", "execution_count": 3, "id": "a254b247", "metadata": { "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAN4AAAF3CAYAAADKA/UlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAeoUlEQVR4nO3de1TUdf7H8dcAoYmJaSZumQIJKoLgclFTQsw2RVdA3W5eMjSLsE23vGCJqOtlNUUU6qiY27Zu1IKTZbSZXdxjpiArGt5WIkRNcPOKKNfP7w9/TI2gMTDD+zszr8c5nRPf78C8+R6ffOf6GZ1SSoGIWpSD9ABE9ojhEQlgeEQCGB6RAIZHJIDhEQlgeEQCnKQHaGmBgYGorKxEp06dpEchG3Pu3Dk4OzsjJyfnVy9rd+FVVFSgpqZGegyyQdXV1Wjs61HsLrx7770XALBz507hScjWDB06tNGX5X08IgEMj0gAwyMSwPCIBDA8IgEMj0gAwyMSwPCIBDA8IgEMj0gAw9MYvV6P4OBguLi4IDg4GHq9XnoksgCGpyF6vR5RUVHIzs5GeXk5srOzER0dzfhsEMPTkCVLltTbppTC0qVLBaYhS2J4GpKfn2/SdrJeDE9DfHx8TNpO1ovhaUh8fDx0Op3RNp1Oh/j4eKGJyFIYnoZERkYiMzPT6FHNrVu3YvTo0dKjkZnZ3TvQtS4yMhKRkZHSY5CF8YxHJIDhEQlgeEQCGB6RAIZHJIDhEQlgeEQCGB6RAIZHJIDhEQlgeEQCGB6RAIZHJEBT4aWmpmLChAlG27744guMGTMGAQEBCA8Px/Lly3H9+nWhCYnMQzPhbd68GcnJyUbbcnJyEBcXh9/97nfQ6/VYsGABsrKykJiYKDQlkXmIh1dSUoIpU6ZgzZo1cHd3N9r33nvvoX///njuuefQrVs3hIaGYsaMGdi2bRsqKyuFJiZqPvE3wubn58PV1RXbtm1DSkoKTp8+bdj37LPPwsGh/t+G6upqlJWVoUOHDi05KpHZiIcXHh6O8PDwBvf17t3b6OvKykq8/fbb8PHxYXRk1cTDa6zq6mrMmjULJ06cwN///nfpcYiaxSrCKysrw8svv4y9e/ciOTkZffv2lR6JqFk0H15paSmmTp2KU6dOYcOGDejfv7/0SETNpunwLl26hEmTJqGsrAxbtmyBt7e39EhEZqHp8JYuXYri4mJs3LgRHTp0wLlz5wz7OnToAEdHR8HpiJpOs+HV1tbik08+QVVVFSZNmlRv/86dO3H//fcLTEbUfJoKb9myZYb/d3BwwMGDBwWnIbIc8VeuENkjhkckgOERCWB4RAIYHpEAhkckgOERCWB4RAIYHpEAhkckgOERCWB4RAIYHpEAhkckgOERCWB4RAIYHpEAhkckgOERCWB4RAIYHpEAhkckgOERCWB4RAIYHpEAhkckgOERCWB4RAIYHpEAhkckgOERCWB4RAIYHpEATYWXmpqKCRMmGG07cuQIxo8fD39/f4SFhSEtLU1oOiLz0Ux4mzdvRnJystG2CxcuYPLkyejevTsyMjIwffp0rFmzBhkZGUJTEpmH+Gegl5SUYN68edi/fz/c3d2N9r3//vtwdnbGggUL4OTkBE9PTxQVFWHDhg0YM2aM0MREzSd+xsvPz4erqyu2bduGvn37Gu3LyclBUFAQnJx+/vvQv39/FBYW4qeffmrpUYnMRvyMFx4ejvDw8Ab3nT17Fl5eXkbb7r33XgDAmTNn0LFjR4vPR2QJ4me827l+/TqcnZ2NtrVq1QoAUFFRITESkVloOrzWrVujsrLSaFtdcG3atJEYicgsNB2em5sbSktLjbbVfd25c2eJkYjMQtPhBQUFYf/+/aipqTFs27NnD9zd3Xn/jqyapsMbM2YMysrKMG/ePJw4cQKZmZn461//imnTpkmPRtQsmg6vY8eO2LhxIwoLCxEVFYV169Zh1qxZiIqKkh6NqFnEn074pWXLltXb5ufnh/T0dIFpiCxH02c8IlvF8IgEMDwiAQyPSADDIxLA8IgEMDwiAQyPSADDIxLA8IgEMDwiAQyPSADDIxLA8IgEMDwiAQyPSADDIxLA8IgEMDwiAQyPSADDIxLA8IgEMDwiAQyPSADDIxLA8IgEMDwiAQyPSADDIxLA8IgEMDwiAQyPSADDIxJgFeFVVVVh9erVCAsLQ0BAAJ566ink5uZKj0XUZFYR3ptvvomMjAwsXrwYer0eHh4emDp1KkpKSqRHI2oSqwhv586dGDlyJAYNGoRu3bphzpw5KCsrw4EDB6RHM4ler0dwcDBcXFwQHBwMvV4vPRIJsYrw2rdvjy+//BKnTp1CTU0N0tPT4ezsjF69ekmP1mh6vR5RUVHIzs5GeXk5srOzER0dzfjslFWEN2/ePDg5OWHo0KHw9fXF6tWrkZSUhAceeEB6tEZbsmRJvW1KKSxdulRgGpJmcnjr1q275X2rU6dOYeHChc0e6mYFBQVo164dUlJSkJ6ejujoaMyePRtHjx41+3VZSn5+vknbybaZHF5KSsotw8vLy8MHH3zQ7KF+6fTp03j11Vfxpz/9CY888gh8fX2xcOFCeHl5Ye3atWa9Lkvy8fExaTvZNqfGXOiJJ55AXl4egBs3jx5//PFbXtbX19c8k/2/gwcPoqqqqt7P7du3L3bt2mXW67Kk+Ph4REdHQyll2KbT6RAfHy84FUlpVHh//vOfkZWVBaUUUlJSMGbMGLi5uRldxsHBAe3atcOjjz5q1gG7dOkCADh27Bj8/PwM248fP45u3bqZ9bosKTIyEpmZmVi6dCny8/Ph4+OD+Ph4jB49Wno0EqBTv/wT3AhPPPEE1qxZg86dO1tqJiO1tbWYMGECzp8/j4SEBLi5uUGv12P9+vXYsmUL/P39Tfp5Q4cOBXDjKQoiczLl35bJ9/GKi4uxb98+06dqIgcHB6SmpqJ///6YO3cuoqOj8e2332Lz5s0mR0ekFY26qflLrVq1QqtWrSwxyy25uroiISEBCQkJLXq9RJZicnjTpk3D/PnzcfToUfTo0QP33HNPvcsEBQWZZTgiW2VyeHVnndTUVAA3Hpmro5SCTqfDkSNHzDQekW0yObx33nnHEnMQ2RWTwwsODrbEHER2xeTwGvOi3sjIyCaMQmQ/TA5vzpw5DW7X6XRwdHSEo6MjwyP6FSaH19CTg+Xl5di/fz/Wr1+PlJQUswxGZMtMDu++++5rcHuPHj1QVVWFRYsWYcuWLc0ejMiWmfX9eF5eXnybC1EjmC28yspKvP/+++jYsaO5fiSRzTL5pmZ4eLjRk+bAjRcyX7hwARUVFZg9e7bZhiOyVU16Hu/m8ACgbdu2GDJkCAYOHGiWwYhsmcnhLVu2zBJz2BW9Xo8lS5YYvS+PT8HYF5PDq3P16lVcuXIFtbW19fb95je/adZQtqxutbE6dauNZWZmMj47YnJ4J0+exMyZM2/76CVfJH1rt1ttjOHZD5PDS0xMxMmTJzFt2jR07doVDg5WsUKgZnC1MQKaEF5ubi4SEhL417mJfHx8kJ2d3eB2sh8mn65cXFzQqVMnS8xiF+Lj4+s9KszVxuyPyeGNHj0a77zzDmpqaiwxj82rW23sl5+hsHXrVq42ZmcadVNz7ty5hv+vrq7Gv//9bwwbNgx+fn648847jS6r0+kafACBfhYZGcmb6nauUeHt3bvX6Ou6NTUPHjxY77INPblORMYaFd4XX3xh6TmI7AqfCyASYJYXSddxcHBAmzZt0K1bN0yYMIHL/BHdgslnvFGjRuHcuXMoLy9HcHAwIiIi0L9/f1RWVuLMmTPo3r07zpw5g0mTJmHPnj2WmJnI6pl8xrt48SJ69+6NtLQ0uLi4GLZXVFRg6tSp6NSpE9asWYM5c+YgNTUVAwYMMOvARLbA5DPep59+iueee84oOuDG0u6TJ0/GRx99BAAYOXIkDh8+bJ4piWxMkx5cuXr1aoPbr1y5gurqagCAk5MTn1ogugWTwxs4cCBWrVpV7x0IR48eRVJSEh566CEAwI4dO+Dp6WmeKYlsjMn38eLj4zFx4kRER0eja9eu6NChA86fP4/i4mJ4eHhg3rx5+Oyzz7BlyxasWbPGEjMTWT2Tw+vUqRM+/PBDbNu2DXv37sX58+fh6emJ2NhYjBo1Co6OjvDw8EB6errRJ7gS0c+a9A50Z2dnjB07FhERESgrK0P79u1xxx13GPY/+OCDZhuQyBY1KbycnBysWLEChw4dQt0nOfv5+WHGjBno37+/WQckskUmP7iSm5uLZ555BpcvX0ZsbCwSEhLwwgsv4OLFi5gyZQr+85//WGJO6PV6jBgxAr6+voiIiEBWVpZFroeoJZh8xktKSkJgYCDS0tLg6Oho2B4XF4eYmBisXbsWmzZtMuuQH374IeLj4zF79myEhYXh448/xsyZM+Hm5oaAgACzXhdRSzD5jHfo0CFMnDjRKDrgxus0x48f3+BbhZpDKYU1a9Zg0qRJmDRpErp164YXX3wRAwcOxL59+8x6XdL0er3RG2Qb85FoZMxajqHJZzwXFxfDk+Q3q6qqMtznM5fvv/8ep0+fxqhRo4y2p6WlmfV6pHHZv+azpmNo8hmvX79+eOutt+q9eqWsrAzr169HYGCg2YYDgB9++AHAjY8Ci4mJwYABAzBu3Dibe4/g7Zb9o8axpmPYqPBWrFiBqqoqAMDMmTNRVFSERx55BHPnzsWqVaswd+5cDBs2DD/88ANmzJhh1gHLysoAALNnz8bIkSOxadMmPPTQQ4iNjbWpdz9w2b/ms6Zj2Kjw0tLSMGbMGBw7dgzdu3dHeno6QkJCsGvXLmzatAm7du1CSEgI3n//ffTs2dOsA9Y9PxgTE4OoqCj06tULL7/8MkJDQ/H222+b9bok3Wp5Py7713jWdAwbFd6GDRtQVlaGsWPHYv369fD09ERSUhJ2796N7777Drt370ZSUpJFnjivW9/Fy8vLaPuDDz6IU6dOmf36pHDZv+azpmPYqPAGDx6Mjz/+GOPGjcPq1avx1FNPobi42NKzAQB69+4NFxcX5OXlGW0/fvw4HnjggRaZoSVw2b/ms6ZjqFMmPgyZm5uLBQsWoLi4GLNnz0ZoaGi9y5j7Q0tSU1OxceNGLFy4EH5+fti+fTuSk5OxefNmhISEmPSzhg4dCqDhz3Inag5T/m2Z/HRCv379kJmZicmTJyMxMbHBy5j7Q0tiY2Nx5513YvXq1SgpKYGnpyfWrl1rcnREWmFyePn5+Vi4cCHy8vIwYsQIDB482BJz1TN58mRMnjy5Ra6LyNIaHV5FRQWSkpLwt7/9De3bt8e6devwyCOPWHI2IpvVqPC++eYbJCQkoLi4GBEREXj99dfRvn17C49GZLsaFd6zzz6Le+65BykpKYY7kETUdI0Kb9SoUXjttdfg6upq6XmI7EKjwluxYoWl5yCyK/zsBCIBDI9IAMMjEsDwiAQwPCIBDI9IAMMjEsDwNEav16NHjx5wdHSEo6MjevToodmVsqjpGJ6G1K2SdeLECdTW1qK2thYnTpxAVFQU47MxDE9DGlolq44WV8qipmN4NzH3uqCmuN1qWJZcKUvyd7ZXTfrQElum0+lw4cKFWy7aa0ne3t63/OwJb29vnDt3zuzX6eTkhLvvvtvsP5duj+E1oLq6WiS8uLg4TJkypcEzUFxcnMhMZBm8qakhw4cPx8aNG+Hu7g4HBwc4ODjA3d0dmzZtwmOPPSY9HpkRz3gaM3z4cAwfPlx6DLIwnvGIBDA8IgEMj0gAwyMSwPCIBDA8IgEMj0gAwyMSwPCIBDA8IgEMj0gAwyMSwPCIBFhVeIWFhQgICEBmZqb0KETNYjXhVVVV4ZVXXkF5ebn0KETNZjXhrV27Fi4uLtJjkMZY63KIVhFednY20tPTsXz5culRSEOseTlEzYd3+fJlzJo1C6+99hq6dOkiPQ5piDUvh6j58BYsWAB/f3+MGjVKehTSGKnlEM1B0+Hp9Xrk5ORgwYIF0qOQBvn4+DRpnxZoOryMjAz89NNPCAsLQ0BAAAICAgAACQkJiIiIEJ6OpMXHx0On091yn5ZpepWxlStX4vr160bbHn30Ubz00ksYMWKE0FSkFZGRkcjMzMSsWbNQUFAAAPD09MSKFSswevRo4eluT9Phde7cucHtHTt2xH333dfC05AWRUZGIjIyUnoMk2n6piaRrdL0Ga8hx44dkx6BqNl4xiMSwPCIBDA8IgEMj0gAwyMSwPCIBDA8IgEMj0gAwyMSwPCIBDA8IgEMj0gAwyMSwPA0LCsrCyNGjICnpydGjBiBrKws6ZFalF6vR3BwMFxcXBAcHKz5lcNMwfA0KisrCzExMThw4ACuXbuGAwcOYMqUKXYTX93SfdnZ2SgvL0d2djaio6NtJj6Gp1Fr166tt00phXXr1glM0/IaWrpPKaX5Zfsai+Fp1K3e8GsvbwS+1fJ8Wl+2r7EYnkZ5e3ubtN3W3Gp5Pq0v29dYDE+jpk+fXm/pOp1Oh+nTpwtN1LIaWrpPp9Npftm+xmJ4GjV8+HBs3LgRAQEBaNOmDQICApCWlobHHntMerQWUbd03y8f1dy6davml+1rLKtb7MieDB8+HMOHD5ceQ4y1Lt3XGDzjEQlgeEQCGB6RAIZHJIDhEQlgeEQCGB41i1JKeoQWZ47fmc/jUbPodDpcuHAB1dXV0qO0CCcnJ9x9993N/zlmmIXsXHV1td2EZy68qUkkgOERCWB4RAKsIryLFy9i/vz5CA0NRb9+/fDkk08iJydHeiyiJrOK8GbOnIm8vDysWrUK//znP+Hj44OYmBgUFBRIj0bUJJoPr6ioCLt370ZCQgICAwPh4eGBefPmoXPnzvj444+lx7Moe19l7Ga2dDw0H97dd9+N9evXo0+fPoZtOp0OSilcunRJcDLLsvdVxm5ma8dD8+G1a9cODz/8MJydnQ3bsrKycPLkSQwaNEhwMsuy91XGbmZrx0Pz4d1s//79iI+Px9ChQxEeHi49jsXY+ypjN7O142FV4X3++eeIiYmBn58fVq1aJT2ORdn7KmM3s7XjYTXhvfvuu5g+fTpCQ0OxYcMGtG7dWnoki7L3VcZuZmvHwyrC27JlCxYtWoSnn34aSUlJRvf3bJW9rzJ2M1s7Hpp/kXRhYSGWLFmCYcOGYdq0afjpp58M+1q3bo277rpLcDrLsvdVxm5mS8dD8+H961//QlVVFXbs2IEdO3YY7YuKisKyZcuEJiNqOs2H9/zzz+P555+XHoPIrKziPh6RrWF4RAIYHpEAhkckgOERCWB4RAIYHpEAhkckgOERCWB4RAIYHpEAhkckgOERCWB4LciWlqezJHs4Tgyvhdja8nSWYi/HieG1EFtbns5S7OU4MbwWYmvL01mKvRwnhtdCbG15Okuxl+PE8FqIrS1PZyn2cpwYXguxteXpLMVejpPmFzuyJba0PJ0l2cNx4hmPSADDIxLA8IgEMDwiAQyPSADDIxLA8IgEMDwiAQyPSADDIxLA8IgEMDwiAQyPSIBVhFdbW4vk5GQMHjwYffv2xbPPPouioiLpsYiazCrCS01NxXvvvYfFixcjPT0dOp0OU6dORWVlpfRoTWYPK2mZg6nHyVqOq+bDq6ysxKZNmzB9+nQ8/PDD6NmzJ1avXo2SkhLs2LFDerwmsZeVtJrL1ONkTcdVp5RS0kPczsGDBzFu3Dh8+umncHd3N2x/8skn4e3tjQULFpj083x9fVFTU4MuXbrc8jK1tbWw5GHJzc3FlStX6m2/66670K9fP4tdb0N0Oh0cHJr399dSx8vU49QSx/V2x+vHH3+Eo6MjDh069Ks/R/PvQD979iwA1Avl3nvvxY8//mjyz2vVqtWv3kRt7j/EX3P16tUGt5eXl8PR0dGi120Jljpeph4n6ePq5OQEZ2fnxl3WwrM027Vr1wCg3i/UqlUrXLp0yeSfl5OTY5a5iJpD8/fxWrduDQD1zlIVFRW48847JUYiajbNh1d3E7O0tNRoe2lpKdzc3CRGImo2zYfXs2dPtG3bFnv37jVsu3z5Mg4fPozAwEDByYiaTvP38ZydnTF+/HisXLkSHTp0wH333YcVK1bAzc0Nw4YNkx6PqEk0Hx4AvPTSS6iursZrr72G69evIygoCGlpaY1+BIlIazT/PB6RLdL8fTwiW8TwiAQwPCIBDI9IAMMjEsDwiAQwPCtw4cIFfPDBB9JjiPD29kZmZqb0GGbH8KzAX/7yF2zbtk16DDIjhmcF+BoH28PwWkh5eTkWL16MQYMGISAgAE8//TQOHjwIAMjIyEBkZCT8/Pzg7++PCRMmID8/HwAwZ84cbN26Ffv27YO3t7fkr2BxZ8+exQsvvICAgACEhYVh+/btRvu/+uor/OEPf0BAQAAGDRqEZcuWoaKiwrD//PnzmDFjBgIDAxESEoIVK1Zg4sSJWLt2bUv/Kr9OUYt47rnnVHh4uPr666/VDz/8oObNm6cCAwPVJ598onx8fNTWrVvVqVOn1IEDB9TYsWPV6NGjlVJKXb58Wf3xj39Ujz/+uCotLZX9JSyoqqpKRUREqMcff1x99913Kjc3V40ePVp5eXmpjIwMtWPHDtWzZ0+1bt06VVBQoL744gsVGhqq4uLilFJK1dTUqLFjx6qoqCiVm5urvvvuOzV+/Hjl7e2tkpOThX+7+hheC/j++++Vl5eX2rVrl2FbRUWFWrJkidq7d6/aunWr0eXT09NVz549DV/Pnj1bjR8/vqXGFbFr1y7l5eWlioqKDNsOHz5sCG/s2LFq+vTpRt+zc+dO5eXlpU6cOKH27NmjvLy8VEFBgWH///73P+Xr66vJ8Kzi3QnW7tixYwAAf39/wzZnZ2fMnTsXAFBQUIDU1FQUFRWhsLAQR44cQW1trcSoYo4fPw5XV1c88MADhm29evUyrDJw/PhxREREGH1PUFAQgBvH9+zZs3B1dYWHh4dhf8eOHY0WyNIS3sdrAU5ON/6+6XS6evu2b9+O3//+9ygqKoKfnx9effVVzJkzp6VH1ATVwINIdcdOKVXv+NXU1Bgu4+joaFV/rBheC/D09AQAo2XfqqurERYWhuTkZIwdOxbLly/H008/jaCgIBQXFwP4+R9iQ8Hamt69e+Py5cv473//a9hWWFhoWK7Py8sL+/fvN/qeuoWrPD090bNnT1y5cgUFBQWG/RcvXtTsiuMMrwW4u7vj0UcfRWJiIvbs2YPCwkLMnz8flZWV6NatG3Jzc5Gfn4+TJ09i8+bNePfddwH8vMBTmzZtUFpaagjSFoWEhKBv376YNWsWDhw4gEOHDmHOnDmGpQNjYmLw2WefISUlBYWFhfjyyy+xaNEiDBkyBJ6enggJCYG/v7/h+48ePYpXXnkF165d0+YfLuH7mHbjypUr6vXXX1chISHK399fTZo0SR05ckSdPHlSPfPMM8rf318FBweriRMnqu3btysvLy+1d+9epZRSBw8eVIMHD1Z+fn6qpKRE+DexnPPnz6uZM2eqgIAANWDAALV582Y1cOBAlZGRoZRS6qOPPlIjR45UPj4+avDgwWr58uXq2rVrhu8/e/asio2NVX379lXBwcEqKSlJhYWFqbfeekvqV7olvgOdbML58+eRl5eHQYMG4Y477gBw4xZDSEgIEhISEBkZKTvgTfioJtkEJycnzJgxA0888QSefPJJVFVVGdblCQ0NlR6vHp7xyGZ8++23SEpKwrFjx6DT6fDb3/4Wr7zyiiZf8cPwiATwUU0iAQyPSADDIxLA8IgEMDwyCz5GZxqGZ4cmTJgAb29vo//69OmDsLAwJCYmmvSBn2fPnsW0adNw+vRpC05se/gEup3q3bs3EhISDF9XVVUhPz8fq1atwpEjR/CPf/yjUa9x/Oabb/DVV1/h9ddft+S4Nofh2am2bdsavT8QuPH+tqtXryI5ORl5eXn19pP58KYmGenTpw8A4MyZMwCAzz//HNHR0fD19cVDDz2ExYsXo7y8HACQmZlpeDPv0KFD7fZ9hE3B8MhIYWEhAKBr16746KOP8OKLL8LDwwMpKSmIi4vDtm3bEBsbC6UUwsLC8MILLwAA1q1bh9jYWMnRrQpvatoppRSqq6sNX1+6dAn79u3Dm2++CX9/f/Tp0wdxcXEYPHgwVq5cabhc9+7d8cwzz+Drr79GWFiYYamGXr164f7772/x38NaMTw7lZ2dDR8fH6NtDg4OGDBgABYtWoTvv//e8IjlLwMNCgpC27ZtsXv3boSFhbXw1LaD4dkpHx8fJCYmArixtESrVq3QpUsXtG3bFgAMyywkJiYaLvdLpaWlLTesDWJ4dsrFxQW+vr633N+uXTsAwKxZsxAcHFxvv6urq8Vmswd8cIUa5OHhgY4dO+LUqVPw9fU1/Ofm5oY33ngDhw8fBgDDmihkGp7xqEGOjo6YMWMG5s+fD0dHRwwZMgSXL19GamoqSkpKDPcP686MO3bsQGhoqGFFNbo9hke3NG7cOLi4uGDjxo1IT09HmzZt0K9fP6xcuRJdu3YFcGN1sIEDB+KNN97Anj17sH79euGprQPfgU4kgDfQiQQwPCIBDI9IAMMjEsDwiAQwPCIBDI9IAMMjEsDwiAQwPCIBDI9IwP8BEw0rcxqo+3sAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=(9,4))\n", "plt.subplot(1,4,1)\n", "sns.barplot(data=yog, x='Pet', y='Yoghurt', color=[0.9,0.9,0.9], errorbar=None)\n", "sns.stripplot(data=yog, x='Pet', y='Yoghurt', color=[0,0,0]) # this bit plots the dots\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "f792710a", "metadata": {}, "source": [ "We note that the mean number of yoghurts eaten is higher for cat owners, but this may be skewed by one cat-owning yogurt enthusiast.\n", " \n", "### Test statistic\n", "\n", "We need to summarize the difference in yoghurt consumption as a number. \n", "\n", "This mean difference, mean(Cat owners) - mean(Dog owners)), is calculated as follows:" ] }, { "cell_type": "code", "execution_count": 4, "id": "8591939b-3d91-4862-bfd7-e9f681a313cf", "metadata": { "tags": [] }, "outputs": [ { "data": { "text/plain": [ "2.5" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "yog.query('Pet == \"cat\"').Yoghurt.mean() - yog.query('Pet == \"dog\"').Yoghurt.mean()" ] }, { "cell_type": "markdown", "id": "cc794ca8", "metadata": {}, "source": [ "On average, cat owners eat yoghurt 2.5 more times per week than dog owners. Possibly because cats stealing yoghurt are cute.\n", "\n", "### Is the result statistically significant? \n", "\n", "Even if there were no true pet-based differences in yoghurt consumption, if we picked a sample of random cat- and dog-owners and measured their yoghurt consumption we wouldn't expect the difference in average yoghurt consupmtion to be *exactly zero* - just by chance the cat owners would sometimes eat more yogurt, and vice versa.\n", "\n", "Is it plausible that the pet-based difference we observed (cat owners eat yoghurt 2.5 more times per week than dog owners) arose due to chance from a situation in which there is no true pet-based difference in yoghurt consumption (ie, the *null hypothesis is true*)?\n", "\n", "In that case, it would be just chance whether a person is labelled as a cat- or dog-owner.\n", "\n", "\n", "### Permutation\n", "\n", "To work out how likely it would be to get a test statistic (difference in mean yogurt consumtion bbetwee the cat- and dog groups) as large as 2.5, we are going to do is shuffle the data around to create many new (re)samples in which, **people are randomly moved between groups**\n", "\n", "Within each resample, some individuals will retain their true pet identity, and some will be relabelled; which individuals are which will differ between resamples.\n", "\n", "For each resample, we then calculate the mean difference in yogurt consumption, mean(Cat-Dog).\n", "\n", "After repeating the process thousands of times, we can ask on what proportion of random resamples we get a values of mean(Cat-Dog) at least as large as our observed difference, 2.5.\n", "\n", "In this case as data are unpaired, we don't need to preserve pared realationships - ie we can swap the pet-labels for any two people in the sample of 18. However, we do still want to preserve the non-pet-based features of the dataset, namely the sample sizes (10 cat owners and 8 dog owners).\n", "\n", "Shuffling within our sample of 18 will also preserve the overall features of the dataset:\n", "\n", "**Non pet-based effects (to be preserved):**\n", "* Most people eat yoghurt less than once a day\n", "* Some people eat loads of yoghurt\n", " * ie the distriution of differences has positive skew \n", "\n", "**Potential pet-based effects (to be \"shuffled out\")**:\n", "* The number of times yoghurt is eaten tends to be higher for cat owners\n", "* The outlier individual who ate loads of yoghurt was a cat owner\n", "\n" ] }, { "cell_type": "markdown", "id": "f6f76266", "metadata": {}, "source": [ "### Run the simulation\n", "\n", "To generate new simulated datasets, we will shuffle around the datapoints in our original dataset. \n", "\n", "Which ones can we shuffle?\n", "\n", "* We assume there are no pet-based differences, so we can swap the pet-labels for any individuals in the sample\n", "* We stick to the original sample sizes, ie 10 people should always be labelled as cat-owners and 8 as dog-owners\n", "\n", "Conceptually, we would achieve this by randomly re-ordering the 18 values in the 'Yogurt' column so that a different 10 fall into rows laelled 'cat'" ] }, { "cell_type": "markdown", "id": "811525bd", "metadata": {}, "source": [ "### Visualizing randoms shuffles in the unpaired data\n", "\n", "Let's see how the difference of means changes over a few random shuffles.\n", "\n", "Below I generate 4 random shuffles of our yoghurt data (in which individuals are randomly relabelled as cat- or dog-owners, maintaining the correct number of cat- and dog-owners), and plot the outcomes:\n", "\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "id": "2d7a4d09", "metadata": {}, "source": [ "Note that **each individual can move between the 'cat' and 'dog' group at random**, but there are always 10 people in the 'cat' group and 8 people in the 'dog' group.\n", "\n", "You can see that the difference of means (difference in heights of the bars) depends on whether more of the high-yoghurt individuals are assigned as cat owners, and also depends quite a lot on which group the 13-yoghurt individual ends up in.\n", "\n", "## Run 10,000 shuffles\n", "\n", "To get a fine-grained estimate of how likely it is to get a test statistic of 2.5 due to chance, we run a *lot* of shuffles - say 10,000.\n", "\n", "We can do this using a function from the library `scipy.stats`, called `scipy.stats.permutation_test()`\n", "\n", "First of all we will run this (using the code block below) and learn about the output. Afterwards we will have a closer look at the Python syntax including defining a function." ] }, { "cell_type": "code", "execution_count": 5, "id": "4d598dbe-4c5c-41f5-9b34-17615be43c7e", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PermutationTestResult(statistic=2.5, pvalue=0.063993600639936, null_distribution=array([-1.1 , 1.6 , -1.1 , ..., -1.775, -0.875, 1.375]))\n" ] } ], "source": [ "def dMeans(x, y):\n", " return np.mean(x)-np.mean(y)\n", "\n", "results = stats.permutation_test((yog.query('Pet == \"cat\"').Yoghurt, yog.query('Pet == \"dog\"').Yoghurt), dMeans, \n", " permutation_type='independent', alternative='greater', n_resamples=10000)\n", "print(results)" ] }, { "cell_type": "markdown", "id": "18b16556", "metadata": {}, "source": [ "\n", "The results are displayed below the code box; there are number of outputs from the test (contained in the parentheses)\n", "\n", "* **statistic** - the test statistic (mean difference) from our actual data, ie 2.5\n", "* **pValue** - the proportion of shuffles in which a result greater than 2.5 was obtained (due to chance)\n", "* **null distribution** - this contains all 10,000 difference of mean values, for the 10,000 shuffles\n", "\n", "### The Null Distribution\n", "\n", "The **null distribution** is the set of 10,000 values we get for the difference of mean yoghurt consumption for (cat owners - dog owners), from our 10,000 shuffles (the equivalent of the red numbers on the 4-permutation example above)\n", "\n", "Let's plot the null distribution in a histogram:\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "480c0e15", "metadata": { "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj8AAAGgCAYAAABFdAY8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAlP0lEQVR4nO3df1SW9f3H8RdCNz9ClhBwm8syzZFrYsW9LzsT4dhoZ3ra5HB2OhVsqTlbDVYt2EnYxJOkJdN0O8wF2FzOIIdjy7azyNZp2+FwwM3ahrTpIdwmdzgwMH7ct/z4/uHxzntoIt5w3ff9eT7O6Rzu67pueF/O2LP7+tz3FTI6OjoqAAAAQ0yzegAAAICpRPwAAACjED8AAMAoxA8AADAK8QMAAIxC/AAAAKMQPwAAwChhVg/gb1JSUuR2uxUfH2/1KAAAYJxOnjwpm82m5ubmSx5L/PwPl8ul4eFhq8cAAACXYWhoSOP93Gbi538kJCRIkg4ePGjxJAAAYLzuvPPOcR/Lmh8AAGAU4gcAABiF+AEAAEYhfgAAgFGIHwAAYBTiBwAAGIX4AQAARiF+AACAUYgfAABgFOIHAAAYhfgBAABGIX4AAIBRiB8AAGAU4gcAABiF+AEAAEYhfgAAgFHCrB4AAKxS9+ZRObv7JUn22CityJhn8UQApgLxA8BYzu5+HXeetnoMAFOMy14AAMAoxA8AADAK8QMAAIxC/AAAAKOw4BkAPgbvCAOCD/EDAB+Dd4QBwYfLXgAAwCjEDwAAMArxAwAAjEL8AAAAoxA/AADAKMQPAAAwCvEDAACMQvwAAACjED8AAMAoxA8AADAK8QMAAIxC/AAAAKMQPwAAwCjEDwAAMArxAwAAjEL8AAAAoxA/AADAKJbHz3/+8x996lOfGvPPvn37JElHjhxRTk6OFi1apIyMDFVVVXk9f2RkRDt27FBaWpqSk5O1atUqtbe3W3EqAAAgAIRZPcC7776r8PBwvf766woJCfFsnz59uk6dOqWVK1fqC1/4gjZs2KDDhw9rw4YNuuaaa5SdnS1JKi8vV3V1tTZt2qTExERt2bJFa9as0YEDB2Sz2aw6LQAGqXvzqJzd/ZIke2yUVmTMs3giAB/H8vj5xz/+oTlz5ighIWHMvt27d8tms6mkpERhYWGaO3eu2tvbVVFRoezsbLndbu3atUsFBQVKT0+XJG3btk1paWmqr6/X8uXLp/p0ABjI2d2v487TVo8BYJwsv+z17rvvat68C/9XUnNzsxwOh8LCPmq01NRUtbW1qaurS62trerr61Nqaqpnf0xMjBYsWKCmpqZJnx0AAAQev3jlJz4+Xvfdd5/ee+893XDDDXr44YeVlpYmp9Op+fPnex1/7hWiEydOyOl0SpJmzpw55piOjo6pOQEAfuf8y1ASl6IAeLM0ftxut9577z1FRkaqsLBQUVFR+vWvf601a9bohRde0ODg4Jh1O+Hh4ZIkl8ulgYEBSbrgMT09PVNzEgD8DpehAHwcS+PHZrOpqalJYWFhnoC59dZbdezYMVVVVSkiIkJut9vrOS6XS5IUFRWliIgISWcj6tzX546JjIycorMAMFEsFL40/owA37P8sldUVNSYbfPnz9cf//hH2e12dXZ2eu079zgxMVFDQ0OebbNnz/Y6JikpaRKnBuALvEJzafwZAb5n6YLn1tZW3XbbbWpubvba/re//U3z5s2Tw+HQoUOHNDw87NnX0NCgOXPmKC4uTklJSYqOjlZjY6Nnf29vr1paWpSSkjJl5wEAAAKHpfEzf/583XzzzdqwYYOam5t17Ngxbdq0SYcPH9ZDDz2k7OxsffjhhyoqKtLRo0e1f/9+7d69W2vXrpV09rJZTk6OysrKdPDgQbW2tuqxxx6T3W5XZmamlacGAAD8lKWXvaZNm6adO3eqrKxMjz76qHp7e7VgwQK98MIL+tSnPiVJqqysVGlpqbKyshQfH6/CwkJlZWV5vkd+fr6GhoZUXFyswcFBORwOVVVV8QGHAPwa70gDrGP5mp/Y2Fg9/fTTF92/cOFC1dTUXHR/aGioCgoKVFBQMBnjAcCkYC0PYB3LP+QQAABgKhE/AADAKMQPAAAwiuVrfgDg4/AhfwB8jfgB4NdYGAzA17jsBQAAjEL8AAAAoxA/AADAKMQPAAAwCvEDAACMQvwAAACjED8AAMAoxA8AADAK8QMAAIxC/AAAAKMQPwAAwCjEDwAAMArxAwAAjMJd3QFMiro3j8rZ3S9JssdGaUXGPIsnAoCziB8Ak8LZ3a/jztNWjwEAY3DZCwAAGIX4AQAARiF+AACAUVjzA+CynL+QWWIxM4DAQ/wAuCwsZAYQ6LjsBQAAjEL8AAAAoxA/AADAKMQPAAAwCvEDAACMwru9ACDAcR814PIQPwAQ4Mbz8QMEEvAR4gcADMDnMwEfYc0PAAAwCvEDAACMQvwAAACjED8AAMAoxA8AADAK8QMAAIxC/AAAAKMQPwAAwCjEDwAAMArxAwAAjEL8AAAAoxA/AADAKMQPAAAwCvEDAACMQvwAAACjED8AAMAoxA8AADAK8QMAAIxC/AAAAKMQPwAAwCh+FT9tbW267bbbtH//fs+2I0eOKCcnR4sWLVJGRoaqqqq8njMyMqIdO3YoLS1NycnJWrVqldrb26d6dAAAECD8Jn7OnDmjJ554Qv39/Z5tp06d0sqVK3XjjTeqtrZWeXl52r59u2praz3HlJeXq7q6Whs3blRNTY1CQkK0Zs0aud1uK04DAAD4Ob+Jnx/+8Ie6+uqrvba9/PLLstlsKikp0dy5c5Wdna0HHnhAFRUVkiS3261du3YpLy9P6enpSkpK0rZt2/T++++rvr7eitMAAAB+zi/ip6mpSTU1NXrmmWe8tjc3N8vhcCgsLMyzLTU1VW1tberq6lJra6v6+vqUmprq2R8TE6MFCxaoqalpyuYHAACBw/L46e3tVWFhoYqLizVz5kyvfU6nU3a73WtbQkKCJOnEiRNyOp2SNOZ5CQkJ6ujomMSpAQBAoLI8fkpKSrRo0SLdfffdY/YNDg7KZrN5bQsPD5ckuVwuDQwMSNIFj3G5XJM0MQAACGRhlz5k8tTV1am5uVmvvPLKBfdHRESMWbh8LmqioqIUEREh6ezan3NfnzsmMjJykqYGAACBzNL4qa2tVVdXlzIyMry2r1+/XlVVVbruuuvU2dnpte/c48TERA0NDXm2zZ492+uYpKSkyR0eAAAEJEvjp6ysTIODg17b7rrrLuXn52vZsmV69dVXVV1dreHhYYWGhkqSGhoaNGfOHMXFxWn69OmKjo5WY2OjJ356e3vV0tKinJycKT8fAADg/yyNn8TExAtuj4uL06xZs5Sdna3KykoVFRXpwQcf1DvvvKPdu3drw4YNks6u9cnJyVFZWZliY2M1a9YsbdmyRXa7XZmZmVN5KgAAIEBYGj+XEhcXp8rKSpWWliorK0vx8fEqLCxUVlaW55j8/HwNDQ2puLhYg4ODcjgcqqqqGrMIGgAAQPLD+Hn33Xe9Hi9cuFA1NTUXPT40NFQFBQUqKCiY7NEAAEAQ8Lv4ATB16t48Kmf32VvK2GOjtCJjnsUTAcDkI34Agzm7+3XcedrqMeCnzo9jiUBG8CB+AAAXRBwjWFn+Cc8AAABTifgBAABGIX4AAIBRiB8AAGAU4gcAABiF+AEAAEYhfgAAgFGIHwAAYBTiBwAAGIX4AQAARiF+AACAUYgfAABgFOIHAAAYhfgBAABGIX4AAIBRiB8AAGAU4gcAABiF+AEAAEYhfgAAgFGIHwAAYJQwqwcAMDnq3jwqZ3e/JMkeG6UVGfMsnggA/APxAwQpZ3e/jjtPWz0GAPgdLnsBAACjED8AAMAoxA8AADAK8QMAAIxC/AAAAKMQPwAAwCjEDwAAMArxAwAAjEL8AAAAoxA/AADAKMQPAAAwCvf2AgIQNy0FgIkjfoAAxE1LAWDiuOwFAACMQvwAAACjcNkLAHBFWIOGQEP8AACuCGvQEGi47AUAAIxC/AAAAKMQPwAAwCgTip+mpib19fVdcF9vb69effXVKxoKAABgskwofr72ta/p2LFjF9zX0tKiJ5988oqGAgAAmCzjfrfXd7/7XXV0dEiSRkdHVVJSoujo6DHHvffee7r22mt9NyEAAIAPjfuVny9+8YsaHR3V6OioZ9u5x+f+mTZtmhYtWqRNmzZNyrAAAABXatyv/CxdulRLly6VJOXm5qqkpERz586dtMEAAAAmw4Q+5PDFF1/09RwAAABTYkLxMzAwoJ07d+r3v/+9BgYGNDIy4rU/JCREr7/+uk8GBAAA8KUJxU9paalqa2v12c9+VrfccoumTZv4xwV1dXVp8+bN+sMf/iCXyyWHw6HCwkLNm3f23jBHjhxRaWmp/va3v+maa65Rbm6uVq9e7Xn+yMiIfvSjH2nfvn3q7e3VHXfcofXr1+uGG26Y8EwAACB4TSh+XnvtNT322GP6xje+ccUDfPOb39S0adNUUVGhqKgobd++XQ888IDq6+s1ODiolStX6gtf+II2bNigw4cPa8OGDbrmmmuUnZ0tSSovL1d1dbU2bdqkxMREbdmyRWvWrNGBAwdks9mueD4AABBcJhQ/Q0NDWrhw4RX/8FOnTumTn/ykvvnNb+rmm2+WJD388MP6yle+on/+859qaGiQzWZTSUmJwsLCNHfuXLW3t6uiokLZ2dlyu93atWuXCgoKlJ6eLknatm2b0tLSVF9fr+XLl1/xjAAAILhM6HrV4sWL9dZbb13xD58xY4a2bt3qCZ///ve/qqqqkt1u17x589Tc3CyHw6GwsI8aLTU1VW1tberq6lJra6v6+vqUmprq2R8TE6MFCxaoqanpiucDAADBZ0Kv/Cxbtkzr169Xd3e3kpOTFRkZOeaYFStWXNb3/N73vqeXX35ZNptNP/7xjxUVFSWn06n58+d7HZeQkCBJOnHihJxOpyRp5syZY44594GMAAAA55tQ/Dz66KOSpLq6OtXV1Y3ZHxISctnx8/Wvf1333HOPXnrpJT3yyCPau3evBgcHx6zbCQ8PlyS5XC4NDAxI0gWP6enpuayfDwAAzDCh+Dl48KCv5/C8u+upp57S4cOHtWfPHkVERMjtdnsd53K5JElRUVGKiIiQJLndbs/X54650KtRAAAAE4qfWbNm+eSHd3V1qaGhQV/60pcUGhoqSZo2bZrmzp2rzs5O2e12dXZ2ej3n3OPExEQNDQ15ts2ePdvrmKSkJJ/MCAAAgsuE4udHP/rRJY/51re+dcljOjs79Z3vfEdxcXH63Oc+J0k6c+aMWlpatHTpUl177bWqrq7W8PCwJ44aGho0Z84cxcXFafr06YqOjlZjY6Mnfnp7e9XS0qKcnJyJnBoAAAhyPo+f6OhoJSQkjCt+kpKStHjxYm3YsEEbN25UTEyMdu7cqd7eXj3wwAMKDw9XZWWlioqK9OCDD+qdd97R7t27tWHDBkln1/rk5OSorKxMsbGxmjVrlrZs2SK73a7MzMyJnBoAAAhyE4qf1tbWMdv6+/t16NAhlZSU6Hvf+964vk9ISIiee+45/eAHP9Cjjz6q06dPKyUlRT//+c913XXXSZIqKytVWlqqrKwsxcfHq7CwUFlZWZ7vkZ+fr6GhIRUXF2twcFAOh0NVVVV8wCEAALigCcXPhURFRSktLU2PPPKInn32Wf3yl78c1/OmT5+ukpISlZSUXHD/woULVVNTc9Hnh4aGqqCgQAUFBRMZGwAAGGbiN+W6iJkzZ+rYsWO+/rYAAAA+4bNXfkZHR9XR0aGKigqfvRsMAADA1yYUP0lJSQoJCbngvtHRUT377LNXNBQAAMBkmVD8PPLIIxeMn+joaGVkZOjGG2+80rkAAAAmxYTiJy8vz9dzAAAATIkJr/lxu93av3+/Ghsb1dvbqxkzZiglJUVZWVme+28BAAD4mwnFT29vr772ta+ptbVV1113neLj49XW1qYDBw7o5z//ufbu3avp06f7elYAAIArNqG3uv/gBz+Q0+nUnj179MYbb6impkZvvPGG9uzZo66uLm3fvt3XcwIAAPjEhOLn4MGDevTRR5WSkuK1PSUlRfn5+Xrttdd8MhwAAICvTeiyV19fn66//voL7rv++uv1wQcfXMlMAIAgU/fmUTm7+yVJ9tgorciYZ/FEMNmEXvm56aab9Pvf//6C+w4ePKgbbrjhioYCAAQXZ3e/jjtP67jztCeCAKtM6JWf1atX6/HHH5fb7dbdd9+ta6+9Vv/973/1yiuvaN++fRe9TxcAAIDVJhQ/y5Yt03vvvaedO3dq3759nu1XXXWVHnnkEd1zzz0+GxAAAMCXJhQ//f39evjhh5WTk6PDhw+rp6dHHR0duueee/SJT3zC1zMCAAD4zGWt+Tly5IhWrFihn/70p5KkmJgYLVmyREuWLNFzzz2n++67jzu6AwAAvzbu+PnXv/6lBx54QD09PZo3z3uVvs1m07p169TX16f77rtPTqfT54MCAAD4wrjj5/nnn9eMGTP0y1/+UnfddZfXvsjISOXk5Ki2tlZRUVHauXOnzwcFAADwhXHHT0NDgx588EFdc801Fz0mLi5OK1euVENDgy9mAwAA8Llxx8/JkyfH9fk98+fP57IXAADwW+OOn9jYWHV2dl7yuO7u7o99dQgAAMBK444fh8Oh/fv3X/K4uro63XLLLVc0FAAAwGQZd/zk5uaqsbFRmzdvlsvlGrPf7XbrmWee0R/+8Afdf//9Ph0SAADAV8b9IYef+cxn9OSTT+rpp5/Wr371K33uc5/TJz/5SQ0PD+vEiRNqbGzUqVOn9O1vf1tpaWmTOTMAIAhx81NMlcv6hOf7779fSUlJqqqq0sGDBz2vAF199dVavHixVq1apeTk5EkZFAAQ3M7d/BSYbJd9e4s77rhDd9xxhyTp1KlTmjZtGre0AAAAAWNC9/Y6Z8aMGb6aAwAAYEpc1r29AAAAAh3xAwAAjEL8AAAAoxA/AADAKMQPAAAwCvEDAACMQvwAAACjXNHn/ADwrfM/3l/iI/4BYDIQP4Af4eP9AWDycdkLAAAYhfgBAABGIX4AAIBRWPMDTJLzFy+zcBkA/AfxA0wSFi8DgH/ishcAADAK8QMAAIxC/AAAAKMQPwAAwCjEDwAAMArxAwAAjEL8AAAAoxA/AADAKMQPAAAwCvEDAACMQvwAAACjED8AAMAolsfPBx98oO9///tasmSJbr/9dt17771qbm727D9y5IhycnK0aNEiZWRkqKqqyuv5IyMj2rFjh9LS0pScnKxVq1apvb19qk8DAAAECMvj5/HHH9fbb7+trVu36he/+IU+/elPa/Xq1Tp27JhOnTqllStX6sYbb1Rtba3y8vK0fft21dbWep5fXl6u6upqbdy4UTU1NQoJCdGaNWvkdrstPCsEu7o3j2rn/ne0c/87qnvzqNXjAAAuQ5iVP7y9vV1/+tOf9NJLL+n222+XJBUVFemtt97SgQMHFBERIZvNppKSEoWFhWnu3Llqb29XRUWFsrOz5Xa7tWvXLhUUFCg9PV2StG3bNqWlpam+vl7Lly+38vQQxJzd/TruPG31GACACbD0lZ8ZM2bo+eef16233urZFhISotHRUfX09Ki5uVkOh0NhYR81Wmpqqtra2tTV1aXW1lb19fUpNTXVsz8mJkYLFixQU1PTlJ4LAAAIDJbGT0xMjNLT02Wz2Tzbfvvb3+r48eNavHixnE6n7Ha713MSEhIkSSdOnJDT6ZQkzZw5c8wxHR0dkzw9AAAIRJZe9vpfhw4d0rp163TnnXdq6dKl2rRpk1cYSVJ4eLgkyeVyaWBgQJIueExPT8/UDA0AsEzdm0fl7O73PLbHRmlFxjwLJ0Ig8Jv4ef311/XEE08oOTlZW7dulSRFRESMWbjscrkkSVFRUYqIiJAkud1uz9fnjomMjJyiyQEAVmH9HSbC8nd7SdKePXuUl5enJUuWqKKiwhMydrtdnZ2dXseee5yYmOi53HWhY/73chkAAIDkB/Gzd+9ePfXUU7r//vv13HPPeV3CcjgcOnTokIaHhz3bGhoaNGfOHMXFxSkpKUnR0dFqbGz07O/t7VVLS4tSUlKm9DwAAEBgsDR+2tra9PTTTyszM1Nr165VV1eXTp48qZMnT+r06dPKzs7Whx9+qKKiIh09elT79+/X7t27tXbtWkln1/rk5OSorKxMBw8eVGtrqx577DHZ7XZlZmZaeWoAAMBPWbrm53e/+53OnDmj+vp61dfXe+3LysrS5s2bVVlZqdLSUmVlZSk+Pl6FhYXKysryHJefn6+hoSEVFxdrcHBQDodDVVVVYxZBAwAASBbHz0MPPaSHHnroY49ZuHChampqLro/NDRUBQUFKigo8PV4AAAgCFm+5gcAAGAqET8AAMAoxA8AADAK8QMAAIxC/AAAAKMQPwAAwCjEDwAAMArxAwAAjEL8AAAAoxA/AADAKMQPAAAwCvEDAACMQvwAAACjED8AAMAoxA8AADAK8QMAAIxC/AAAAKMQPwAAwCjEDwAAMArxAwAAjEL8AAAAoxA/AADAKMQPAAAwCvEDAACMQvwAAACjED8AAMAoYVYPAPibujePytndL0myx0ZpRcY8iycCAPgS8QP8D2d3v447T1s9BgBgknDZCwAAGIX4AQAARiF+AACAUYgfAABgFOIHAAAYhfgBAABGIX4AAIBRiB8AAGAU4gcAABiF+AEAAEYhfgAAgFG4txeMcf4NSyVuWgqYhBsW43zED4zBDUsBc/HvP87HZS8AAGAU4gcAABiF+AEAAEYhfgAAgFGIHwAAYBTiBwAAGIX4AQAARiF+AACAUYgfAABgFOIHAAAYhfgBAABGIX4AAIBR/Cp+ysvLlZub67XtyJEjysnJ0aJFi5SRkaGqqiqv/SMjI9qxY4fS0tKUnJysVatWqb29fSrHBgAAAcRv4uenP/2pduzY4bXt1KlTWrlypW688UbV1tYqLy9P27dvV21treeY8vJyVVdXa+PGjaqpqVFISIjWrFkjt9s91acAAAACQJjVA7z//vsqKirSoUOHNGfOHK99L7/8smw2m0pKShQWFqa5c+eqvb1dFRUVys7Oltvt1q5du1RQUKD09HRJ0rZt25SWlqb6+notX77cilMCAAB+zPJXfv7+97/rE5/4hH79618rOTnZa19zc7McDofCwj5qtNTUVLW1tamrq0utra3q6+tTamqqZ39MTIwWLFigpqamKTsHAAAQOCx/5Wfp0qVaunTpBfc5nU7Nnz/fa1tCQoIk6cSJE3I6nZKkmTNnjjmmo6NjEqYFAASrujePytndL0myx0ZpRcY8iyfCZLE8fj7O4OCgbDab17bw8HBJksvl0sDAgCRd8Jienp6pGRIAEBSc3f067jxt9RiYApZf9vo4ERERYxYuu1wuSVJUVJQiIiIk6YLHREZGTs2QAAAgoPh1/NjtdnV2dnptO/c4MTHRc7nrQsfY7fapGRIAAAQUv44fh8OhQ4cOaXh42LOtoaFBc+bMUVxcnJKSkhQdHa3GxkbP/t7eXrW0tCglJcWKkQEAgJ/z6/jJzs7Whx9+qKKiIh09elT79+/X7t27tXbtWkln1/rk5OSorKxMBw8eVGtrqx577DHZ7XZlZmZaPD0AAPBHfr3gOS4uTpWVlSotLVVWVpbi4+NVWFiorKwszzH5+fkaGhpScXGxBgcH5XA4VFVVNWYRNAAAgORn8bN58+Yx2xYuXKiampqLPic0NFQFBQUqKCiYzNEAAECQ8OvLXgAAAL5G/AAAAKMQPwAAwCh+teYHGC8+hh6AFfjdExyIHwQkPoYegBX43RMcuOwFAACMQvwAAACjED8AAMAoxA8AADAK8QMAAIxC/AAAAKMQPwAAwCjEDwAAMArxAwAAjEL8AAAAoxA/AADAKMQPAAAwCvEDAACMQvwAAACjED8AAMAoxA8AADAK8QMAAIwSZvUAAAAEs7o3j8rZ3e95bI+N0oqMeRZOBOIHAIBJ5Ozu13HnaavHwHm47AUAAIxC/AAAAKMQPwAAwCjEDwAAMAoLnjGlzn/Xw8Xe8TCeYwAAmCjiB1NqPO964J0RAIDJxGUvAABgFOIHAAAYhcteAABYjE+BnlrEDwAAFmOt49TishcAADAK8QMAAIxC/AAAAKMQPwAAwCjEDwAAMArxAwAAjEL8AAAAoxA/AADAKMQPAAAwCvEDAACMQvwAAACjcG8v+AQ35QMABAriB2OcHzLjjRhuygcACBTED8YgZAAAwYw1PwAAwCjEDwAAMAqXvQwzkfU8AAD/MFm/w0170wrxYxjW8wBA4BrP73DetHJpQXHZa2RkRDt27FBaWpqSk5O1atUqtbe3Wz0WAABT7lzIHHee9no1Bx8JivgpLy9XdXW1Nm7cqJqaGoWEhGjNmjVyu91WjwYAAPxMwF/2crvd2rVrlwoKCpSeni5J2rZtm9LS0lRfX6/ly5dbPOGVM+1aLAAAkyng46e1tVV9fX1KTU31bIuJidGCBQvU1NQUFPEz3muxLGYGAFgtEP6/KGR0dHTU6iGuxGuvvaa8vDy9/fbbioiI8Gz/9re/rcHBQf3kJz+5rO/3mc98RsPDw5o5c6avR52w/sEzGhr+6H+msNAQRUVc9bHHTfYxzMiMzMiMzOi7Gaf6PCYyo8s9rOHzkiE0JEThtlCf/Hxf6OjoUGhoqP76179e8tiAf+VnYGBAkmSz2by2h4eHq6en57K/X3h4uN+tFRrvX5zxHOerYyb6HGb0zXOY0TfPYUbfPIcZp+77WjnjhULHVz/fF8LCwsa0wEWPneRZJt25V3vcbrfXKz8ul0uRkZGX/f2am5t9NhsAAPA/Af9ur3OXpzo7O722d3Z2ym63WzESAADwYwEfP0lJSYqOjlZjY6NnW29vr1paWpSSkmLhZAAAwB8F/GUvm82mnJwclZWVKTY2VrNmzdKWLVtkt9uVmZlp9XgAAMDPBHz8SFJ+fr6GhoZUXFyswcFBORwOVVVVjXvhEwAAMEfAv9UdAADgcgT8mh8AAIDLQfwAAACjED8AAMAoxA8AADAK8QMAAIxC/AAAAKMQPwAAwCjEDyakublZt9xyi9dtReA7HR0devzxx/X5z39eDodDq1ev1j//+U+rxwoKIyMj2rFjh9LS0pScnKxVq1apvb3d6rGC0gcffKDvf//7WrJkiW6//Xbde++93Dx6CrS1tem2227T/v37rR7FbxE/uGynT59WYWGhRkZGrB4lKLndbn3jG99QV1eXfvKTn2jv3r2aPn26vv71r6u7u9vq8QJeeXm5qqurtXHjRtXU1CgkJERr1qyR2+22erSg8/jjj+vtt9/W1q1b9Ytf/EKf/vSntXr1ah07dszq0YLWmTNn9MQTT6i/v9/qUfwa8YPLVlJSouuvv97qMYJWc3Oz/vGPf+jZZ5/VrbfeqptvvlnPPvus+vv79cYbb1g9XkBzu93atWuX8vLylJ6erqSkJG3btk3vv/++6uvrrR4vqLS3t+tPf/qT1q9fr5SUFN10000qKipSYmKiDhw4YPV4QeuHP/yhrr76aqvH8HvEDy7Lr371K/3lL3/RunXrrB4laN188816/vnnlZiY6LV9dHRUPT09Fk0VHFpbW9XX16fU1FTPtpiYGC1YsEBNTU0WThZ8ZsyYoeeff1633nqrZ1tISAh/jydRU1OTampq9Mwzz1g9it8LihubYmr8+9//VmlpqcrLy/kvi0kUHx+v9PR0r20/+9nP5HK59PnPf96iqYKD0+mUJM2cOdNre0JCgjo6OqwYKWjFxMSM+Xv829/+VsePH9fixYstmip49fb2qrCwUMXFxWP+fmMs4geSzobNnXfeedH9b731lgoLC3XPPfcoJSVF//73v6dwuuByqT/rP/7xj4qPj/c8fu2117Rt2zbl5uYqKSlpKkYMWgMDA5Ikm83mtT08PJxXIybZoUOHtG7dOt15551aunSp1eMEnZKSEi1atEh333231aMEBOIHkqTExET95je/uej+ffv2qb+/X3l5eVM4VXC61J91bGys5+uXXnpJTz31lJYtW6Ynn3xyKsYLahEREZLOrv0597UkuVwuRUZGWjVW0Hv99df1xBNPKDk5WVu3brV6nKBTV1en5uZmvfLKK1aPEjBCRkdHR60eAv5v6dKl6uzs1FVXXSXp7PqTgYEBhYeH67Of/awqKystnjD4lJWVqaKiQrm5uSoqKlJISIjVIwW8d955R1/96ldVX1+v2bNne7bfe++9SkpK0vr16y2cLjjt2bNHpaWlyszMVFlZ2ZhX3XDlcnNz9ec//9nrz7a/v182m02zZ8/Wq6++auF0/olXfjAuL774ooaGhjyP33//feXm5mrjxo36v//7PwsnC05btmxRZWWlCgsLtXr1aqvHCRpJSUmKjo5WY2OjJ356e3vV0tKinJwci6cLPnv37tVTTz2l3NxcrVu3TtOm8R6byVBWVqbBwUGvbXfddZfy8/O1bNkyi6byb8QPxmXWrFlej0NDQyWdvYTzv+9KwpVpbGxUZWWlcnNz9eUvf1knT5707IuKimKx+RWw2WzKyclRWVmZYmNjNWvWLG3ZskV2u12ZmZlWjxdU2tra9PTTTyszM1Nr165VV1eXZ19ERISmT59u4XTB5WK/g+Pi4sb87sZZxA/gZ859BsqLL76oF1980Wvft771LdZdXaH8/HwNDQ2puLhYg4ODcjgcqqqq4nKMj/3ud7/TmTNnVF9fP+YzlLKysrR582aLJgNY8wMAAAzDBVgAAGAU4gcAABiF+AEAAEYhfgAAgFGIHwAAYBTiBwAAGIX4AQAARiF+AACAUYgfAABgFOIHAAAYhfgBAABG+X9aYcQs3M5ylwAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sns.histplot(results.null_distribution, bins=np.arange(-5,5,0.1))\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "8ede169d-d08c-4614-bc72-de564aabbdc6", "metadata": {}, "source": [ "The histogram looks a bit spikey because for this dataset (which is small and has several repeated values) there are a limited number of possible values for the means of a groups of 8 and a group of 10.\n", "\n", "In any case we can see that the observed difference between cat and dog owners' mean consumption (+2.5) is quite rare but not completely extreme." ] }, { "cell_type": "markdown", "id": "4e7f3153", "metadata": {}, "source": [ "### The $𝑝$-valueΒΆ\n", "\n", "We can calculate the exact proportion of cases in which the difference of means for cat- minus dog-owners exceeds the value we observed in our original sample, 2.5. " ] }, { "cell_type": "code", "execution_count": 8, "id": "d6344afc-b707-4f4b-9350-675a5118151d", "metadata": { "tags": [] }, "outputs": [ { "data": { "text/plain": [ "0.0639" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# check which values is the null distribution are greater than 4.6 and code as 1/0\n", "# then use the numpy function np.mean() to get the avergae of those ones and zeros (ie he proportion of 1s)\n", "np.mean(results.null_distribution >= 2.5)" ] }, { "cell_type": "markdown", "id": "5689e6df-b5ee-4c13-9767-7b70e0b02ee6", "metadata": {}, "source": [ "This proportion is 0.064 or 6.4%.\n", "\n", "This tells us that if we simulate a situation in which pet ownership does not determine the number of times yoghurt is eaten (but preserving some other important features of the dataset like the high skew), there is still almost a 5% chance that we would get an apparent pet-based difference as large as the one we observed in our 'real' data.\n", "\n", "The probability that the test statistic (in this case, the difference in mean times yoghurt is eaten) would be observed if the null hypothesis were true, is sometimes called the $𝑝$-value.\n", "\n", "Our permutation test shows that the $𝑝$-value associated with the observed difference of means is 0.047 or 4.7%.\n", "\n", "The result is considered statistically significant if $𝑝$ is smaller than some predetermined level, known as $\\alpha$. Usually $\\alpha=0.05$ or $\\alpha=0.01$ is used; this result is not quite significant if $\\alpha=0.05$ but not significant if $\\alpha=0.01$. " ] }, { "cell_type": "markdown", "id": "24efb854-8ae2-4879-ae68-a3eba5b3d722", "metadata": {}, "source": [ "## Syntax of `stats.permutation_test()`\n", "\n", "Let's have a look at the Python code to run the permutation test.\n", "\n", "Firstly, we had to give the function `stats.permutation_test()` our two samples. In this case, the data are in 'long form', ie the pet category is defined in one column of the dataset, so cat owners and dog owners are in different rows. We therefore identify our categories useing `df.query()` with the column we want (Yoghurt) tacked on the end:\n", "\n", "`(yog.query('Pet == \"cat\"').Yoghurt, yog.query('Pet == \"dog\"').Yoghurt)`\n", "\n", "Secondly, to tell `stats.permutation_test()` the test statistic we want to get the null distribution of, we had to pass it a *function* that calculates that test statistic. We created that function, `dMeans()`, using a function definition, before running the permutation test\n", "\n", "\n", "### Arguments\n", "\n", "When running the permutation test, we have to specify some arguments\n", "\n", "* **n_resamples** the number of shuffles to do.\n", "* **alternative** - whether to do a one-tailed test (and in whcih direction: 'greater' counts how many shuffles result in test statistic values greater than the one obtained from unshuffled data) or a two tailed test\n", "* **permutation_type** - pairings, samples, or independent. These are for test of correlation, paired data and unpaired data respectively.\n", "\n", "Here we are testing for a difference of means in *independent* (unpaired) samples, so we choose `permutation_type='independent'` - this randomly shuffles all the datapoints and assignes them randomly to the groups cat/dog." ] }, { "cell_type": "markdown", "id": "a25bc28d-0a21-4428-8044-173683688fde", "metadata": {}, "source": [ "## Recap\n", "\n", "To run a permutation test on unpaired data, we mixed up all the individuals in the sample regardless of category (pet ownership), and randomly split them into two groups. We kept the original group sizes (10 and 8 in this case), as we want to retain the characteristics of the original dataset that one group is larger\n", "\n", "For each shuffle we calculated the mean yoghurt consuption for each 'pet' group. \n", "\n", "Permutation testing in this way gives us a null distribution for the difference of means. Values of difference of means that occur rarely in the null distribution are considered statistically significant.\n", " \n", "To run the permutation test with scipy.stats we need the option `permutation_type='independent'`" ] }, { "cell_type": "code", "execution_count": null, "id": "65731b03-9c05-41a2-9386-64aa63dde21c", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.7" } }, "nbformat": 4, "nbformat_minor": 5 }