{ "cells": [ { "cell_type": "markdown", "id": "c8868457", "metadata": {}, "source": [ "# Permutation test for correlation\n", "\n", "In the previous examples we used permutation testing to assess the significance of a difference between groups (difference of means or mean difference).\n", "\n", "Permutation testing can also be used to assess the statistical significance of a correlation.\n", "\n", "As a reminder, a correlation can occur only in paired designs, as when two variables are correlated, it means that an individual's score on one variable is related to their score on the other variable. \n", "\n", "Correlations can be interesting in themselves (do students who score highly on English tests also score highly on maths tests?; do people who eat more broccli have greater bone density?). \n", "\n", "They can also reflect the fact that experimental measures often depend on factors other than the one we are manipulating (sometimes called confounding factors), which are what we try to control for by using a paired design. For example if we are interested in whether men earn more than women, we might use a paired design comparing brothers and sisters to take into account the very important effects of parental occupation and education on earnings which mean that high-earning brothers often have high-earning sisters. The fact that brothers' and sisters' earnings are correlated actually reflects the confounds that we want to 'cancel out' by using a paired design to test gender differences.\n", "\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": "5fb0416d", "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": "f0806e8c", "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", "We are interested in whether people who eat more broccoli have higher IQs.\n", "\n", "#### Question & design\n", "\n", "We hypothesise that those wo eat more broccoli ave higher IQs \n", "\n", "This is tecnically a *repeated measures design* as we have two measurements (broccoli consumption and IQ) for each individual. \n", "\n", "#### Hypotheses\n", "\n", "We can state our hypotheses as follows:\n", "\n", "$\\mathcal{H_o}:$ There is no relationsip between broccoli consumption and IQ\n", "* the correlation, Pearson's $r=0$\n", "\n", "$\\mathcal{H_a}:$ Those with higher broccoli consumption have higher IQ\n", "* the correlation, Pearson's $r>0$\n", "\n", "This is a one-tailed (directional) alternative hypothesis\n", "\n", "#### Data\n", "\n", "The following made-up data give weekly broccoli consumption in grams and IQ for 25 individuals:" ] }, { "cell_type": "code", "execution_count": 2, "id": "537072c8", "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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
broccoli_gIQ
0087
12891
20101
32092
4096
59295
68892
712894
8096
92299
1011499
11096
1214699
13255108
14131100
15255107
16390114
17402107
18216108
19719104
20395107
21485114
22553116
23682116
24815111
\n", "
" ], "text/plain": [ " broccoli_g IQ\n", "0 0 87\n", "1 28 91\n", "2 0 101\n", "3 20 92\n", "4 0 96\n", "5 92 95\n", "6 88 92\n", "7 128 94\n", "8 0 96\n", "9 22 99\n", "10 114 99\n", "11 0 96\n", "12 146 99\n", "13 255 108\n", "14 131 100\n", "15 255 107\n", "16 390 114\n", "17 402 107\n", "18 216 108\n", "19 719 104\n", "20 395 107\n", "21 485 114\n", "22 553 116\n", "23 682 116\n", "24 815 111" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "broccoli = pd.read_csv('https://raw.githubusercontent.com/jillxoreilly/StatsCourseBook_2024/main/data/broccoli.csv')\n", "broccoli" ] }, { "cell_type": "markdown", "id": "325c73df", "metadata": {}, "source": [ "Let's plot the data:" ] }, { "cell_type": "code", "execution_count": 3, "id": "356b33e8", "metadata": { "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj8AAAG1CAYAAAAWb5UUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAx6ElEQVR4nO3dfXRU9Z3H8U9gMsmMMTFgSNACRiA8E1IJBCkPhod6BNmDqBwQDoKiForCqqDFLVDEdkVQQJEFolgtK1ZYiq4Pq9hWWzAEFJUIQSBBQPJgCGCYZCYJd//IydQhaUogmbk39/06J+eY328y+X4v6Hy8v3vvL8wwDEMAAAA20SLUBQAAAAQT4QcAANgK4QcAANgK4QcAANgK4QcAANgK4QcAANgK4QcAANiKI9QFmE3fvn3l8/kUFxcX6lIAAMBFKioqktPp1O7du//lawk/F/B6vaqqqgp1GQAAoAEqKyt1sc9tJvxcoE2bNpKk7du3h7gSAABwsYYNG3bRr+WaHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCvs7QUAuCie8goVlZSpzFspV4RDcbEuuSPDQ12WKXBsrIXwAwD4lwpOeZSZna9Sj88/FuV2qn+PBMW3coewstDj2FgPy14AgHp5yitqfbhLUqnHp8zsfHnKK0JUWehxbKyJ8AMAqFdRSVmtD/capR6fikrKglyReXBsrInwAwCoV5m3sv55X/3zzRnHxpoIPwCAerki6r881OW07+WjHBtrIvwAAOoVF+tSlNtZ51yU26m4WFeQKzIPjo01EX4AAPVyR4arf4+EWh/yUW6n0nom2PqWbo6NNXE+DgDwL8W3cmt4arvqZ9n4KuVy8iybGhwb6yH8AAAuijsyXB3a8oFeF46NtbDsBQAAbMV04Wf16tWaPHlyrfHc3Fz16dNHx48fDxjftWuXunTpUutrx44dwSoZAABYiKmWvTZs2KCVK1cqNTU1YDwnJ0f333+/yspqPywqJydH7du318aNGwPGY2JimrRWAABgTaYIPwUFBZo/f7727NmjxMTEgLkXX3xRa9asUceOHXXy5MlaP3vw4EF17txZcXFxwSoXAABYmCmWvbKzsxUTE6Nt27YpOTk5YO6TTz7R0qVLNW/evDp/NicnR506dQpGmQAAoBkwxZmf9PR0paen1zlXs5yVmZlZa84wDH3zzTeKi4vTbbfdpoKCAiUlJWnOnDnq3bt3k9YMAACsyRRnfi7V8ePH5fF45PP59Otf/1qrV69Wq1atNGnSJB06dCjU5QEAABMyxZmfS9WuXTvt3r1bbrdbLVu2lCQtXbpUo0eP1quvvqpFixaFuEIAAGA2lg4/knTllVcGfN+iRQt16tRJBQUFIaoIAACYmaWXvf7yl7+oT58+AXeBVVZW6sCBA1wEDQAA6mTp8NO3b1+1bt1ac+fOVXZ2tnJycjRv3jydPn1ad999d6jLAwAAJmTp8BMVFaUNGzYoNjZW06ZN0/jx43X69Gm99tpruvrqq0NdHgAAMKEwwzCMUBdhJsOGDZMkbd++PcSVAACAi9WQz2/LX/AMAJ7yChWVlKnMWylXhENxsS65I627w3Zz6wcwG8IPAEsrOOVRZna+Sj0+/1iU26n+PRIU38odwsouTXPrBzAjS1/zA8DePOUVtYKCJJV6fMrMzpenvCJElV2a5tYPYFaEHwCWVVRSViso1Cj1+FRUUhbkii5Pc+sHMCvCDwDLKvNW1j/vq3/ebJpbP4BZEX4AWJYrov7LFl1Oa13W2Nz6AcyK8APAsuJiXYpyO+uci3I7FRfrCnJFl6e59QOYFeEHgGW5I8PVv0dCrcAQ5XYqrWeC5W4Pb279AGbFOVQAlhbfyq3hqe2qn4vjq5TLae3n4jS3fgAzIvwAsDx3ZLg6tG0+4aC59QOYDcteAADAVgg/AADAVgg/AADAVgg/AADAVgg/AADAVgg/AADAVgg/AADAVgg/AADAVgg/AADAVgg/AADAVgg/AADAVgg/AADAVgg/AADAVgg/AADAVgg/AADAVgg/AADAVgg/AADAVgg/AADAVgg/AADAVgg/AADAVgg/AADAVgg/AADAVgg/AADAVgg/AADAVgg/AADAVhyhLgAAANiDp7xCRSVlKvNWyhXhUFysS+7I8KDXQfgBAABNruCUR5nZ+Sr1+PxjUW6n+vdIUHwrd1BrYdkLAAA0KU95Ra3gI0mlHp8ys/PlKa8Iaj2EHwAA0KSKSspqBZ8apR6fikrKgloP4QcAADSpMm9l/fO++ucbG+EHAAA0KVdE/ZcYu5zBvQSZ8AMAAJpUXKxLUW5nnXNRbqfiYl1BrYfwAwAAmpQ7Mlz9eyTUCkBRbqfSeiYE/XZ3bnUHAABNLr6VW8NT21U/58dXKZeT5/wAAIBmzh0Zrg5tgx92LsSyFwAAsBXThZ/Vq1dr8uTJtcZzc3PVp08fHT9+PGDc6/Vq0aJFGjBggFJSUvTggw+quLg4WOUCAACLMVX42bBhg1auXFlrPCcnR1OnTlVZWe2HIC1cuFB///vftWrVKr3yyis6duyYHnrooWCUCwAALMgU4aegoED33nuvVqxYocTExIC5F198UXfeeadatWpV589t3bpVTzzxhPr27avevXtr+fLlysrK0t69e4NUPQAAsBJThJ/s7GzFxMRo27ZtSk5ODpj75JNPtHTpUs2bN6/Wz+3Zs0eS1L9/f/9YYmKi4uPjlZWV1bRFAwAASzLF3V7p6elKT0+vc27jxo2SpMzMzFpzBQUFio2NVURERMB4mzZtdPLkycYvFAAAWJ4pzvxcqrKyMjmdtZ8YGRERIa/XG4KKAACA2Vk6/ERGRsrnq71LrNfrlcsV3EdlAwAAa7B0+ElISNDp06drBaDCwkIlJCSEqCoAAGBmlg4/N9xwg86fP++/8FmSjhw5ooKCAvXt2zeElQEAALOydPiJj4/XqFGj9MQTTygzM1NffvmlHn74YfXr1099+vQJdXkAAMCELB1+JGnx4sUaMGCAfvnLX+qee+7R9ddfX+eDEgEAACQpzDAMI9RFmMmwYcMkSdu3bw9xJQAA4GI15PPb8md+AAAAGoLwAwAAbIXwAwAAbIXwAwAAbIXwAwAAbIXwAwAAbIXwAwAAbIXwAwAAbIXwAwAAbIXwAwAAbIXwAwAAbIXwAwAAbMUR6gIAq/OUV6iopExl3kq5IhyKi3XJHRlu2vcFALsj/ACXoeCUR5nZ+Sr1+PxjUW6n+vdIUHwrt+neFwDAshdwyTzlFbUCiiSVenzKzM6Xp7zCVO8LAKhG+AEuUVFJWa2AUqPU41NRSZmp3hcAUI3wA1yiMm9l/fO++ueD/b4AgGqEH+ASuSLqv2TO5by0S+qa6n0BANUIP8Aliot1KcrtrHMuyu1UXKzLVO8LAKhG+AEukTsyXP17JNQKKlFup9J6JlzybelN9b4AgGqcPwcuQ3wrt4antqt+Ho+vUi5n4zyPp6neFwBA+AEumzsyXB3aNn4oaar3BQC7Y9kLAADYCuEHAADYCstegEmxt1fDXXjMoqOcOlvq4xgCCED4AUyIvb0a7sJjVlpWobLyCnW7rpVOnS3XeYNjCKAay16AybC3V8NdeMwqKqt0orBU+cUefXGoSDFREZI4hgCqEX4Ak2Fvr4a78JidK6uUr6JKknTqjFeGYfjnOIYACD+AybC3V8NdeMwqq84HfF9RGfg9xxCwN8IPYDLs7dVwFx4zR8vA/7SFOwK/5xgC9kb4AUyGvb0a7sJjdoXLIWd4S0lSq5gIhYWF+ec4hgAIP4DJsLdXw114zMIdLXVtmygltHYruXOczpR6JXEMAVTj3C9gQuzt1XB1HTP/c344hgB+hPADmBR7ezVcXccs9srIEFUDwKxY9gIAALZC+AEAALZC+AEAALZC+AEAALZC+AEAALZC+AEAALZC+AEAALZC+AEAALZC+AEAALZC+AEAALZC+AEAALZC+AEAALZiqvCzevVqTZ48OWBs//79mjRpkvr06aOhQ4cqIyMjYP7EiRPq0qVLra8//vGPwSwdAABYhGl2dd+wYYNWrlyp1NRU/1hJSYmmTp2q4cOHa9GiRdq7d68WLVqkq666SuPGjZMk5eTkKCIiQh9++KHCwsL8P3vllVcGvQcAAGB+IQ8/BQUFmj9/vvbs2aPExMSAuTfeeENOp1MLFy6Uw+FQx44ddfToUa1bt84ffg4ePKjExES1adMmFOUDAACLCfmyV3Z2tmJiYrRt2zYlJycHzO3evVupqalyOP6R0dLS0pSbm6vi4mJJ1Wd+OnXqFNSaAQCAdYX8zE96errS09PrnMvPz1dSUlLAWM0Znu+++06tW7fWwYMHFRcXp4kTJyovL08dOnTQjBkzNGjQoCavHQAAWE/Iz/zUp7y8XE6nM2AsIiJCkuT1euXz+ZSXl6fS0lLNnj1ba9euVa9evTR9+nTt3LkzFCUDAACTC/mZn/pERkbK5/MFjHm9XkmS2+2W0+lUVlaWHA6HPyT17NlThw8fVkZGhgYMGBD0mgEAgLmZ+sxPQkKCCgsLA8Zqvo+Pj5f0jxD0Y0lJSSooKAhOkQAAwFJMHX5SU1O1Z88eVVVV+cd27typxMREtW7dWgcOHFBKSop2794d8HP79u3jImgAAFAnU4efcePGqbS0VPPnz9ehQ4e0ZcsWvfLKK7r//vslVZ/h6dy5sxYtWqTdu3fr8OHD+u1vf6u9e/fqgQceCHH1AADAjEwdflq3bq3169crNzdXY8eO1fPPP6+5c+dq7NixkqQWLVpozZo16tWrl2bPnq2xY8fqiy++0Msvv6wuXbqEuHoAAGBGYYZhGKEuwkyGDRsmSdq+fXuIKwEAABerIZ/fpr7bCwAAu/KUV6iopExl3kq5IhyKi3XJHRke6rKaBcIPAAAmU3DKo8zsfJV6/vG4lyi3U/17JCi+lTuElTUPpr7mBwAAu/GUV9QKPpJU6vEpMztfnvKKEFXWfBB+AAAwkaKSslrBp0apx6eikrIgV9T8EH4AADCRMm9l/fO++ufxrxF+AAAwEVdE/Zfjupxcrnu5CD8AAJhIXKxLUW5nnXNRbqfiYl1Brqj5IfwAAGAi7shw9e+RUCsARbmdSuuZwO3ujYBzZwAAmEx8K7eGp7arfs6Pr1IuJ8/5aUyEHwAATMgdGa4ObQk7TYFlLwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuOS/3BY8eOafv27Tpx4oQMw1Dbtm01ePBgde7cWZL06quvKjo6Wv/2b//WaMUCAABcrgaHn9LSUi1evFjbtm2TYRgBc88884xGjhypRx55RM8//7wWLVrUaIUCAAA0hgaFn6qqKs2YMUN79+7Vvffeq1GjRql9+/Zq2bKlvv32W73//vvKyMjQn//8ZyUnJ+vmm29uqroBAAAuSYPCz+bNm/XFF1/o1VdfVXJycsBc586d1blzZ8XExGjJkiUaOnRoY9YJAADQKBp0wfOWLVs0YcKEWsGnRmVlpdavX69rr71W//d//9coBQIAADSmBoWf3NxcDRgw4J/OHz58WFdccYWeeOIJ5ebmXnZxAAAAja1B4aeqqkoRERH/dL5Lly5655135HK5VFVVddnFNSee8godPXlWB/JO6ejJs/KUV4S6JAAAbKlB1/y0b99ee/fuVVpaWr2v+/zzz9WhQ4fLKqw5KTjlUWZ2vko9Pv9YlNup/j0SFN/KHcLKAACwnwad+Rk1apQ2bNiggoKCf/qa7777Tr///e81atSoyy6uOfCUV9QKPpJU6vEpMzufM0AAAARZg8LPxIkTdfXVV2v8+PF65513VF5e7p8rKyvTW2+9pQkTJiguLk533XVXoxdrRUUlZbWCT41Sj09FJWVBrggAAHtr0LKXy+XSunXr9PDDD+vf//3f5XA4dNVVVyksLEwlJSWqrKxU3759tWzZMkVGRjZVzZZS5q2sf95X/zwAAGhcDX7Cc9u2bfWHP/xBO3fu1EcffeTf3uInP/mJ0tPTNWDAAIWFhTVFrZbkiqj/ELucl7zDCAAAuASX9MkbFhamG2+8UTfeeGNj19PsxMW6FOV21rn0FeV2Ki7WFYKqAACwrwaFn65du170WZ2wsDB9/fXXl1RUc+KODFf/Hgl13u2V1jNB7sjwEFYHAID9NCj8zJw5kyWtSxDfyq3hqe1UVFKmMl+lXE6H4mJdBB8AAEKgQeFn1qxZTVVHs+eODFeHtoQdAABCrUG3ugMAAFidqcLP6tWrNXny5ICx/fv3a9KkSerTp4+GDh2qjIyMgPnz589r5cqVGjRokJKTkzVt2jQdPXo0mGUDAAALMU342bBhg1auXBkwVlJSoqlTp+q6667T5s2bNWvWLK1YsUKbN2/2v2b16tV6/fXX9eSTT2rTpk0KCwvT9OnT5fPV/WBBmBN7nwEAgiXkD5kpKCjQ/PnztWfPHiUmJgbMvfHGG3I6nVq4cKEcDoc6duyoo0ePat26dRo3bpx8Pp9eeuklPfrooxoyZIgk6dlnn9WgQYP0wQcfsMWGRbD3GQAgmEJ+5ic7O1sxMTHatm2bkpOTA+Z2796t1NRUORz/yGhpaWnKzc1VcXGxDhw4oHPnzgVstBodHa3u3bsrKysraD3g0rH3GQAg2EJ+5ic9PV3p6el1zuXn5yspKSlgrE2bNpKqN1DNz8+XVP3U6Qtfc/LkySaoFo3tYvY+4y45AEBjCvmZn/qUl5fL6XQGjEVEREiSvF6vysqqNwWt6zVerzc4ReKysPcZACDYTB1+IiMja124XBNq3G63f/PUul7jcrFthBWw9xkAINhMHX4SEhJUWFgYMFbzfXx8vH+5q67XJCQkBKdIXJaavc/qwt5nAICmYOrwk5qaqj179qiqqso/tnPnTiUmJqp169bq2rWroqKilJmZ6Z8/e/asvv76a/Xt2zcUJaOBavY+uzAAsfcZAKCpmHpNYdy4cVq/fr3mz5+ve++9V19++aVeeeUVLVq0SFL1tT6TJk3SM888o1atWunaa6/V0qVLlZCQoBEjRoS4elws9j4DAASTqcNP69attX79ei1ZskRjx45VXFyc5s6dq7Fjx/pf8+CDD6qyslJPPPGEysvLlZqaqoyMjFoXQcPc2PsMABAsYYZhGKEuwkyGDRsmSdq+fXuIKwEAABerIZ/fpj7zY2ee8orqZSBvpVwR1l4GaqxemtMxAQCEDuHHhJrTdg+N1UtzOiYAgNAy9d1edtSctntorF6a0zEBAIQe4cdkLma7B6torF6a0zEBAIQe4cdkmtN2D43VS3M6JgCA0CP8mExz2u6hsXppTscEABB6hB+TaU7bPTRWL83pmAAAQo/wYzLNabuHxuqlOR0TAEDosV5gQs1pu4fG6qU5HRMAQGgRfkyqOW330Fi9NKdjAgAIHZa9AACArRB+AACArRB+AACArRB+AACArRB+AACArRB+AACArRB+AACArRB+AACArRB+AACArRB+AACArRB+AACArRB+AACArRB+AACArRB+AACArRB+AACArRB+AACArRB+AACArRB+AACArRB+AACArRB+AACArRB+AACArRB+AACArRB+AACArRB+AACArRB+AACArThCXYBdFJZ4dPS7syot8ynK7VSHttFqE+sOdVkAANgO4ScI9ucW64/bDyrv5A/+sevaXqk7hiWpW2LrEFYGAID9sOzVxApLPLWCjyTlnfxBf9x+UIUlnhBVBgCAPRF+mtjR787WCj418k7+oKPfnQ1yRQAA2Bvhp4mVlvkuax4AADQuwk8Ti3I5L2seAAA0LsJPE+twTbSua3tlnXPXtb1SHa6JDnJFAADYG+GnibWJdeuOYUm1AtB1ba/UncO7cLs7AABBxq3uQdAtsbUeGJf8j+f8uJzqcA3P+QEAIBQIP0HSJtZN2AEAwARY9gIAALZiiTM/586d0/Lly/Xhhx+qtLRU/fv312OPPab27dtLkh5//HFt2bIl4Gfi4+P18ccfh6JcAABgYpYIP7Nnz1ZOTo4WLVqk9u3ba8OGDZowYYLefvttxcbGKicnRw888IAmTZrk/5mWLVuGsOLmzVNeoaKSMpV5K+WKcCgu1iV3ZHioywIA4KKYPvwcOHBAH3/8sdauXashQ4ZIkhYtWqRdu3Zp48aNeuCBB3To0CHNmDFDcXFxIa62+Ss45VFmdr5KPf94OGOU26n+PRIU34prmgAA5mf6a35yc3MlSX379vWPtWjRQl27dlVWVpby8vLk9XrVsWPHUJVoG57yilrBR5JKPT5lZufLU14RosoAALh4pg8/NWdz8vPzA8ZPnDih4uJiHTx4UGFhYXrllVeUnp6u4cOHa/Hixfrhh7r308KlKyopqxV8apR6fCoqKQtyRQAANJzpw09ycrI6duyoBQsW6OTJk/L5fNqwYYP2798vn8+nb775Ri1atNC1116rNWvWaN68efrrX/+qGTNm6Pz586Euv1kp81bWP++rfx4AADMw/TU/4eHheuGFF/TYY49p6NChcjgcGjp0qG6//Xbt27dPs2bN0t13363o6OptIpKSkhQXF6fx48frq6++UnJycog7aD5cEfX/dXE5Tf/XCQAA84cfSUpMTNSmTZt05swZhYWFKTo6Wg899JCuu+46//c/lpSUJKl6qYzw03jiYl2KcjvrXPqKcjsVF+sKQVUAADSM6Ze9SktLNWnSJO3bt08xMTGKjo7WDz/8oB07dmjQoEF6+OGHdc899wT8zFdffSVJ6tSpUyhKbrbckeHq3yNBUe7Aneij3E6l9UzgdncAgCWY/sxPVFSUwsLC9NRTT2nBggUyDEOLFy/WNddco9GjRysmJka/+MUv9OKLL2rUqFHKzc3Vb37zG40ePZo7wJpAfCu3hqe2q37Oj69SLifP+QEAWEuYYRhGqIv4VwoLC/Xkk09q586datGihW666SbNnTtXrVq1kiS9//77WrNmjY4cOaIrr7xSt956q2bPnq2IiIgG/65hw4ZJkrZv396oPQAAgKbTkM9vS4SfYCL8AABgPQ35/Db9sheaP7bLAAAEE+EHIcV2GQCAYDP93V5ovtguAwAQCoQfhAzbZQAAQoHwg5BhuwwAQCgQfhAybJcBAAgFwg9Cpma7jLqwXQYAoKkQfhAybJcBAAgF1hUQUmyXAQAINsIPQs4dGa4ObQk7AIDgYNkLAADYCuEHAADYCuEHAADYCuEHAADYCuEHAADYCuEHAADYCuEHAADYCuEHAADYCuEHAADYCuEHAADYCuEHAADYCuEHAADYCuEHAADYCru6m5CnvEJFJWUq81bKFeFQXKxL7sjQ7npuxpoAALgUhB+TKTjlUWZ2vko9Pv9YlNup/j0SFN/KTU0AAFwmlr1MxFNeUStkSFKpx6fM7Hx5yiuoCQCAy0T4MZGikrJaIaNGqcenopKyIFdkzpoAALgchB8TKfNW1j/vq3++KZixJgAALgfhx0RcEfVfguVyBv8SLTPWBADA5SD8mEhcrEtRbmedc1Fup+JiXUGuyJw1AQBwOQg/JuKODFf/Hgm1wkaU26m0ngkhubXcjDUBAHA5WLMwmfhWbg1PbVf9TB1fpVzO0D9Tx4w1AQBwqQg/JuSODFeHtuYKFmasCQCAS8GyFwAAsBXCDwAAsBWWvSyAfbUAAGg8hB+TY18tAAAaF8teJsa+WgAAND7Cj4mxrxYAAI2P8GNi7KsFAEDjI/yYGPtqAQDQ+Ag/Jsa+WgAAND7Cj4mxrxYAAI2PdROTY18tAAAaF+HHAthXCwCAxsOyFwAAsBVLhJ9z585p8eLFGjJkiG644QbNmDFD3377rX9+//79mjRpkvr06aOhQ4cqIyMjhNUCAAAzs0T4mT17tj744AMtXLhQb7zxhlq3bq0JEyaopKREJSUlmjp1qq677jpt3rxZs2bN0ooVK7R58+ZQlw0AAEzI9Nf8HDhwQB9//LHWrl2rIUOGSJIWLVqkXbt2aePGjXI4HHI6nVq4cKEcDoc6duyoo0ePat26dRo3blyIqwcAAGZj+jM/ubm5kqS+ffv6x1q0aKGuXbsqKytLu3fvVmpqqhyOf+S4tLQ05ebmqri4OOj1AgAAczN9+ImLi5Mk5efnB4yfOHFCxcXFys/PV0JCQsBcmzZtJEnfffddcIoEAACWYfrwk5ycrI4dO2rBggU6efKkfD6fNmzYoP3798vn86m8vFxOZ+BDACMiIiRJXq83FCUDAAATM334CQ8P1wsvvKCKigoNHTpUKSkpysrK0u23366oqChFRkbK5wvc+bwm9Ljd7lCUDAAATMz0FzxLUmJiojZt2qQzZ84oLCxM0dHReuihh3Tdddfp7NmzKiwsDHh9zffx8fGhKBcAAJiY6c/8lJaWatKkSdq3b59iYmIUHR2tH374QTt27NCgQYOUmpqqPXv2qKqqyv8zO3fuVGJiolq3bh3CygEAgBmZPvxERUUpLCxMTz31lHJycnTgwAE98MADuuaaazR69GiNGzdOpaWlmj9/vg4dOqQtW7bolVde0f333x/q0gEAgAmZPvxI0rJly3T11Vdr0qRJmjJlitq1a6eXX35ZDodDrVu31vr165Wbm6uxY8fq+eef19y5czV27NhQlw0AAEzIEtf8tGnTRitXrvyn871799amTZuCWBEAALAqS5z5AQAAaCyEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuEHwAAYCuOUBdgF57yChWVlKnMWylXhENxsS65I8NDXRYAALZD+AmCglMeZWbnq9Tj849FuZ3q3yNB8a3cIawMAAD7YdmriXnKK2oFH0kq9fiUmZ0vT3lFiCoDAMCeCD9NrKikrFbwqVHq8amopCzIFQEAYG+EnyZW5q2sf95X/zwAAGhchJ8m5oqo/7Iql5PLrgAACCbCTxOLi3Upyu2scy7K7VRcrCvIFQEAYG+EnybmjgxX/x4JtQJQlNuptJ4J3O4OAECQseYSBPGt3Bqe2q76OT++SrmcPOcHAIBQIfwEiTsyXB3aEnYAAAg1lr0AAICtEH4AAICtEH4AAICtEH4AAICtEH4AAICtEH4AAICtEH4AAICtEH4AAICtEH4AAICtEH4AAICtsL3FBQoLC1VVVaVhw4aFuhQAAHCRTp48qZYtW17Uaznzc4GIiAg5HGRCAACsxOFwKCIi4qJeG2YYhtHE9QAAAJgGZ34AAICtEH4AAICtEH4AAICtEH4AAICtEH4AAICtEH4AAICtEH4AAICtEH4AAICtEH4AAICtEH4AAICtEH4AAICtEH4AAICtEH6C4Pz581q5cqUGDRqk5ORkTZs2TUePHg11WQ22evVqTZ48OWBs//79mjRpkvr06aOhQ4cqIyMjYN7MvZ8+fVq//vWvNXjwYP30pz/VhAkTtHv3bv+8lXuTpOLiYj366KNKS0tTSkqK7rvvPh06dMg/b/X+auTm5iolJUVbtmzxj1m5txMnTqhLly61vv74xz9KsnZvNbZu3apbbrlFvXr10qhRo/Tuu+/656zaX2ZmZp1/bl26dNGwYcMkWbc3SaqoqNCzzz6roUOHKiUlRRMnTtRnn33mn7dcbwaa3KpVq4wBAwYYf/nLX4z9+/cb06ZNM0aMGGF4vd5Ql3bRXn75ZaNLly7GpEmT/GOnTp0y+vfvb8yfP984dOiQ8eabbxq9evUy3nzzTf9rzNz71KlTjTFjxhhZWVnG4cOHjcWLFxu9e/c2Dh06ZPneDMMw7rjjDmP8+PHGl19+aRw6dMiYNWuWMXDgQMPj8TSL/gzDMHw+n3HbbbcZSUlJxubNmw3DsP7fy+3btxu9evUyCgoKjMLCQv9XWVmZ5XszDMPYunWr0a1bN2PDhg1GXl6e8fzzzxtdu3Y1PvvsM0v35/V6A/68CgsLjb/97W9G9+7djTfeeMPSvRmGYaxYscIYOHCg8cknnxh5eXnG/PnzjZ/+9KdGfn6+JXsj/DQxr9drpKSkGBs3bvSPnTlzxujdu7fx9ttvh7Cyi5Ofn2/cc889Rp8+fYybb745IPysWbPGGDRokFFRUeEfW7ZsmfHzn//cMAxz956Xl2ckJSUZe/bs8Y+dP3/eGDFihPHcc89ZujfDqA4Ac+bMMQ4ePOgf279/v5GUlGR88cUXlu+vxrJly4zJkycHhB+r9/biiy8aY8aMqXPO6r2dP3/euOmmm4zf/e53AePTpk0z1qxZY/n+fszn8xmjRo0yZs+ebRiG9f/sxowZY/z2t7/1f//DDz8YSUlJxnvvvWfJ3lj2amIHDhzQuXPnlJaW5h+Ljo5W9+7dlZWVFcLKLk52drZiYmK0bds2JScnB8zt3r1bqampcjgc/rG0tDTl5uaquLjY1L3HxsZq7dq16tmzp38sLCxMhmHozJkzlu5Nqu5v+fLl6ty5syTp+++/V0ZGhhISEtSpUyfL9ydJWVlZ2rRpk/7zP/8zYNzqveXk5KhTp051zlm9tyNHjujEiRO69dZbA8YzMjJ0//33W76/H/vDH/6gkydP6vHHH5dk/T+7q666Sn/+8591/PhxVVVVadOmTXI6nerWrZsleyP8NLH8/HxJUtu2bQPG27Rpo5MnT4aipAZJT0/XsmXL1K5du1pz+fn5SkhICBhr06aNJOm7774zde/R0dEaMmSInE6nf+zdd9/Vt99+q5/97GeW7u1C//Ef/6GBAwfqvffe05IlS+R2uy3f39mzZzV37lw98cQTtWq0em8HDx5UcXGxJk6cqBtvvFETJkzQJ598Isn6veXl5UmSPB6P7rnnHg0YMEB33HGHPvroI0nW76+G1+vVmjVrNGXKFH/9Vu9t/vz5cjgcGjZsmHr16qVnn31Wzz33nNq3b2/J3gg/TaysrEySAj5kJSkiIkJerzcUJTWa8vLyOvuSqv/lt1Lve/bs0a9+9SsNGzZM6enpzaq3KVOmaPPmzRozZoxmzpyp7Oxsy/e3cOFC9enTp9YZBMnafy99Pp/y8vJUWlqq2bNna+3aterVq5emT5+unTt3Wro3SSotLZUkzZs3T6NHj9ZLL72kgQMHasaMGc2ivxp/+tOf5PV6A24QsXpvhw8fVnR0tF544QVt2rRJt912m+bNm6cDBw5YsjfHv34JLkdkZKSk6v+o1fyzVP0XwuVyhaqsRhEZGSmfzxcwVvMX2e12W6b3Dz/8UI888oiSk5O1fPlySc2nN0n+JZTFixdr7969eu211yzd39atW7V792699dZbdc5buTen06msrCw5HA7/B0XPnj11+PBhZWRkWLo3SQoPD5ck3XPPPRo7dqwkqVu3bvr666/18ssvW76/Glu3btXIkSMVGxvrH7NybydOnNCjjz6qDRs2qG/fvpKkXr166dChQ1q1apUle+PMTxOrOc1XWFgYMF5YWFjrNKHVJCQk1NmXJMXHx1ui99dee02zZs3S4MGDtW7dOv+/mFbvrbi4WG+//baqqqr8Yy1atFDHjh39NVq1v82bN6u4uNh/y21KSookacGCBRo1apSle5OqPywu/D/kpKQkFRQUWL63mhqSkpICxjt16qTjx49bvj9JOnXqlD7//HPdcsstAeNW7u3LL79URUWFevXqFTCenJysvLw8S/ZG+GliXbt2VVRUlDIzM/1jZ8+e1ddff+1P0FaVmpqqPXv2BHzA7ty5U4mJiWrdurXpe9+4caMWL16su+66S88991zAB47VeyssLNTDDz+sXbt2+ccqKir09ddfq2PHjpbu75lnntE777yjrVu3+r8k6cEHH9TatWst3duBAweUkpIS8LwpSdq3b586depk6d4kqXv37rriiiv0xRdfBIwfPHhQ7du3t3x/kvTZZ58pLCxM/fr1Cxi3cm814SUnJydg/ODBg+rQoYM1ewvJPWY2s3z5cqNfv37Ghx9+6H++wciRI03x7IaGmDdvXsCt7t9//72RmppqzJs3z/jmm2+MzZs3G7169TK2bNnif41Zez9y5IjRo0cPY+bMmbWezXH27FlL92YY1bcUT5s2zfj5z39uZGVlGTk5OcacOXOM1NRU48SJE5bv70I/vtXdyr1VVVUZd9xxhzF69GgjKyvLOHTokPHUU08ZPXv2NA4cOGDp3mq88MILRkpKivHWW28ZR48eNVavXm107drV+PTTT5tFf6tWrTJGjhxZa9zKvVVVVRkTJ040br75ZmPnzp1Gbm6u8eyzzxrdunUzPv/8c0v2RvgJgsrKSuPpp5820tLSjD59+hjTp083jh07FuqyGuzC8GMYhvHFF18Yd955p9GzZ0/jpptuMl599dWAebP2/uKLLxpJSUl1fs2bN88wDOv2VuPs2bPGggULjIEDBxq9e/c2pk2bFvDcH6v392M/Dj+GYe3eiouLjccff9wYOHCg0atXL2P8+PFGVlaWf97KvdV46aWXjPT0dKNHjx7GmDFjjA8++MA/Z/X+FixYYNx55511zlm5t9OnTxsLFy40hg4daqSkpBjjx483MjMz/fNW6y3MMAwjNOecAAAAgo9rfgAAgK0QfgAAgK0QfgAAgK0QfgAAgK0QfgAAgK0QfgAAgK0QfgAAgK0QfgAAgK0QfgA0mfT0dD322GOhLuOydOnSRatWrZIkZWZmqkuXLgF7FAGwHkeoCwAAq+jRo4c2bdqkTp06hboUAJeB8AMAFykqKkp9+vQJdRkALhPLXgCaVEVFhZ588kmlpqYqNTVV8+bN06lTpyRJjz32mKZMmaIFCxaob9++Gjt2rCorK+X1evXCCy/o5ptvVq9evTRy5EitXbtW58+fD3jv//3f/9Vtt92m5ORkDR06VEuXLpXP5/PP79u3T/fee69uuOEGpaWlac6cOTp58qR/vrCwUI8//riGDBmi3r176/bbb9f27dv/aS+XuuxVWFioOXPmqF+/fkpNTdWvf/1rPfvss0pPT2/Q+wBoHJz5AdCk3n33XfXu3Vu/+93vdOrUKT3zzDM6evSoXn/9dUnS7t27FRYWplWrVuncuXNq2bKlpk+frr1792rmzJnq1q2bMjMz9dxzz+nYsWNavHixJOn111/XggULdPvtt2vOnDk6fvy4nn76aZWUlOipp57SgQMHNGHCBP/vNgxDy5Yt07Rp07Rt2zadOXNGt99+u8LDwzVnzhzFxsZqy5Ytmjlzpp5++mmNGTOmUfr3+XyaMmWKPB6PfvWrXykqKkpr167V/v37FRcX1yi/A0DDEH4ANKno6GitX79eUVFRkqTY2FjNnDlTf/vb3yRJlZWVWrRokTp06CBJ+utf/6odO3Zo6dKl/gAycOBARUZGasWKFZoyZYquv/56rVq1SiNGjNCSJUv8v8vr9ep//ud/5PP5tHr1asXExOill15SRESEJCkhIUGzZ89WTk6O3n33XZ06dUrvvvuu2rVrJ0kaMmSI7r77bj399NMaPXq0WrS4/JPj27Zt05EjR7R582b17NlTkpSWlqbhw4df9nsDuDQsewFoUkOGDPEHH6n6DrDw8HDt2LFDkhQZGan27dv753ft2qWWLVvqlltuCXifmiCUmZmp3Nxcff/997UCxN13360//elPcjqd2rNnjwYPHuwPPpLUu3dvffTRR+rZs6d27dqllJQUf/D58e8pKirSkSNHGqX/Tz/9VO3atfMHH6n62qGbbrqpUd4fQMMRfgA0qauvvjrg+xYtWuiqq67S2bNnJUmtW7dWWFiYf/7MmTOKjY2VwxF4YrpmieiHH37Q6dOn/T/7z5w+fbre+TNnztSq7cf11tR3uUpKSuqso67fDSA4CD8AmtSFIaKqquqfBgJJiomJUUlJiSorKwPGCwsLJVUvm0VHR0uS/8LpGqdPn9bf//53nTt3TldeeWWteal6Wa2goEAxMTH6/vvva80XFRX5f09jiI+PV3Fxca3xusYABAfhB0CT2rFjR0CQef/991VZWan+/fvX+fp+/fqpqqpK77zzTsD4tm3bJEk33HCDrr/+esXGxta6M+utt97S9OnT5fV61bdvX33yyScBd3/l5OTovvvu01dffaXU1FR9/vnnOnbsWK3fExcX578G6XL169dPx44d0/79+/1jXq9XH3/8caO8P4CG44JnAE3q+++/16xZszR58mTl5eVp+fLlGjhwoAYMGOAPND82ePBg9e/fXwsWLFBhYaG6d++uXbt2ad26dRo7dqz/AYOzZs3Sb37zGy1cuFAjRoxQXl6ennvuOU2YMEGtWrXSjBkzNH78eE2fPl1TpkyRz+fTihUr1KNHDw0ePFjJycnatm2bpk6dql/+8peKjY3V1q1b9emnn+qpp55qlIudJWn06NFau3atZs6cqYceekjR0dF66aWXVFxcrGuuuaZRfgeAhiH8AGhSd955p8rLyzVz5kw5nU7deuutevTRRwOu8/mxsLAw/dd//ZdWrlyp3//+9zp16pR+8pOfaM6cOZo6dar/dXfddZfcbrcyMjL05ptvKj4+XtOmTdN9990nSerevbteffVVLVu2THPmzNEVV1yhIUOG6JFHHpHT6VRcXJz++7//W8uWLdOSJUtUUVGhrl27avXq1Ro2bFij9e9wOJSRkaElS5Zo4cKFcjgcGjNmjGJjY5Wbm9tovwfAxQszDMMIdREA0Fx98803OnLkiEaOHBkQ+MaNG6e2bdvq+eefD2F1gD1x5gcALkFVVZX+1f87hoWFyePx6KGHHtLEiRM1YsQIVVVV6e2331Z2drYeffTRIFUL4Mc48wMAlyA9PV0nTpyo9zXXXnutPvroI7333nvKyMjQ4cOHZRiGunfvrl/84hf62c9+FqRqAfwY4QcALkFOTk7AnWR1cTqd6tKlS5AqAnCxCD8AAMBWeM4PAACwFcIPAACwFcIPAACwFcIPAACwFcIPAACwFcIPAACwFcIPAACwlf8HMUS0R0CWl70AAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sns.scatterplot(data=broccoli, x='broccoli_g', y='IQ', alpha=0.5)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "e65b5443", "metadata": {}, "source": [ "We can see that there does seem to be a correlation. \n", "\n", "### Test statistic\n", "\n", "Here our test statistic is Pearson's $r$" ] }, { "cell_type": "code", "execution_count": 4, "id": "d93aa2ed", "metadata": { "tags": [] }, "outputs": [ { "data": { "text/plain": [ "0.8011532312132079" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "broccoli.broccoli_g.corr(broccoli.IQ)" ] }, { "cell_type": "markdown", "id": "9b7aafdf", "metadata": {}, "source": [ "The correlation is 0.80, which is actually very strong (remember they are made-up data!)." ] }, { "cell_type": "markdown", "id": "ed8de2f6", "metadata": {}, "source": [ "### Is this result statistically significant? \n", "\n", "That is, would the result (a correlation of 0.80) be very unlikely to occur due to chance, if the null hypothesis were true? \n", "\n", "Even if there was no relationship between broccoli consumption and IQ, it would sometimes happen that for 25 random people the ones with the highest IQ tend to also be the biggest broccoli-eaters, resulting in a positive correlation. The question is how often a positive correlation as large as $r$=0.80 would occur just due to chance if broccoli consuption and IQ were paired randomly.\n", "\n", "\n", "### Permutation\n", "\n", "The sample tells us several interesting things about broccoli consumption and IQ, regardless of whether broccoli consummption and IQ are related:\n", "\n", "**Independent features of Broccoli distribution (to be retained):**\n", "\n", "* Most people eat between 0 and 500g of broccoli a week\n", "* Quite a few people eat 0g of broccoli (they never eat broccoli, basically)\n", "\n", "It also tells us about some potential effects to do with the relationship between IQ and broccoli:\n", "\n", "**Possible correlation effects (to be shuffled out)\n", "* the quantity of broccoli eaten seems to be higher for individuals with higher IQ\n", "* none of the highest IQ people are the 0g of broccoli people\n", "\n", "\n", "### Visualizing random shuffles in correlation \n", "\n", "What we are going to do is shuffle the data around to create many new (re)samples preserving distribution within each variable (broccoli consumption and IQ - so for example there will always be 5 non-broccoli-eaters) but **pairing the datapoints** at random so one person's IQ is matched with another person's broccoli consumption. \n", "\n", "Using these simulated (shuffled) datasets we will work out how often we get a correlation of 0.80 or more. This is equivalent to determining how likely our correlation is to have occurred due to chance.\n", "\n" ] }, { "cell_type": "markdown", "id": "86435bf3", "metadata": {}, "source": [ "Below you see jointplots for the original data and three random permutations (joint plots are scatter plots with histograms of the distributions of x and y at the sides, created using the function `sns.jointplot()`)\n", "\n", "\n", "Notice that in the shuffled data, when each person's broccoli consupmtion is randomly paired with another person's IQ, the correlation goes away. However, the distribution of each variable by itself (the histograms in the margins) remain unchanged.\n", "\n", "" ] }, { "cell_type": "markdown", "id": "e9a292f4", "metadata": {}, "source": [ "You can see that the relationship between broccoli consumption and IQ looks less tight in the shuffled datasets.\n", "\n" ] }, { "cell_type": "markdown", "id": "23046d85", "metadata": {}, "source": [ "## Run 10,000 shuffles\n", "\n", "To get a fine-grained estimate of how likely it is to get a test statistic of r=0.8 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.\n" ] }, { "cell_type": "code", "execution_count": 5, "id": "216aabae-a11d-440d-af33-d4d35e100f61", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PermutationTestResult(statistic=0.8011532312132079, pvalue=9.999000099990002e-05, null_distribution=array([-0.08879756, 0.08141763, 0.1368379 , ..., 0.28560897,\n", " -0.49493745, -0.04183689]))\n" ] } ], "source": [ "def correlate(x, y):\n", " tmp = np.corrcoef(x,y)\n", " c = tmp[0][1] \n", " return c\n", "\n", "results = stats.permutation_test((broccoli.broccoli_g, broccoli.IQ), correlate, \n", " permutation_type='pairings', alternative='greater', n_resamples=10000)\n", "print(results)" ] }, { "cell_type": "markdown", "id": "3621e34c-11a3-4450-8778-59692c21c55e", "metadata": { "tags": [] }, "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 (correlation) from our actual data, ie r=0.8\n", "* **pValue** - the proportion of shuffles in which a result greater than 0.8 was obtained (due to chance) - note this is very small, $9.9 \\times 10^{-5}$ or about 1/10000 \n", "* **null distribution** - this contains all 10,000 correlation 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 correlation, 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:" ] }, { "cell_type": "code", "execution_count": 6, "id": "63f2fab4-9b0c-49ad-b6cc-515a5d4e42fc", "metadata": { "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkgAAAGgCAYAAABR4ZjdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAuL0lEQVR4nO3dfVjVZZ7H8c8RBGTQMlKwJ3V0kHEbxeQkbpJmg7vlTBPLtTNN4mOZTgaZBVnaijtamqbpNOakWLaNIzmYlbU7PuQ07WYG9qxgakq2ihihJMI5Avf+0XImfnCUczicw8P7dV1el9z37xy+Xw7Kh999n9/PZowxAgAAgEunQBcAAADQ2hCQAAAALAhIAAAAFgQkAAAACwISAACABQEJAADAgoAEAABgERzoAlqb+Ph4OZ1O9ejRI9ClAACAJjp16pRCQkKUn5/vk+cjIFk4HA7V1NQEugwAAOCB6upq+fLa1wQki549e0qSdu7cGeBKAABAU918880+fT72IAEAAFgQkAAAACwISAAAABYEJAAAAAsCEgAAgAUBCQAAwIKABAAAYEFAAgAAsCAgAQAAWBCQAAAALAhIAAAAFgQkAAAACwISAACABQEJAADAgoAEAABgQUACAACwCA50AQDgT+eqzqu6xridDw6yKTyssx8rAtAaEZAAdCjVNUaL1ue5nZ890e7HagC0VgQkAPAAZ6CAjoGABAAe4AwU0DGwSRsAAMCCgAQAAGBBQAIAALAgIAEAAFgQkAAAACwISAAAABa8zR8Avsdmk8ornG7njXF/DSQA7QcBCQC+p7bWaPGL+W7nH54Q78dqAAQKS2wAAAAWBCQAAACLgAek//3f/9WAAQMa/Nm0aZMkqaCgQKmpqYqLi9OoUaOUnZ1d7/G1tbVauXKlEhMTNXjwYE2ZMkVFRUWBaAUAALQTAd+DdODAAYWGhmrHjh2y2Wyu8a5du6qsrEyTJ0/WT3/6U82fP18fffSR5s+fr0svvVQpKSmSpFWrVmnjxo164oknFBUVpSVLlmjq1KnaunWrQkJCAtUWAABowwIekD7//HP17dtXPXv2bDC3fv16hYSEKCsrS8HBwerXr5+Kioq0Zs0apaSkyOl0at26dcrIyNDIkSMlScuXL1diYqK2b9+usWPH+rsdAADQDgR8ie3AgQPq379/o3P5+fmy2+0KDv57jktISNCRI0dUWlqqwsJCVVRUKCEhwTXfrVs3DRw4UHl57u+2DQAAcCEBD0iff/65SktLdeedd+of//Ef9etf/1rvvPOOJKm4uFjR0dH1jq8703T8+HEVFxdLknr16tXgmBMnTvihegAA0B4FdInN6XTq6NGj6tKlizIzMxUeHq7XXntNU6dO1fPPP6+qqqoG+4hCQ0MlSQ6HQ5WVlZLU6DFnzpzxTxMAAKDdCWhACgkJUV5enoKDg10h59prr9Xhw4eVnZ2tsLAwOZ31r2jrcDgkSeHh4QoLC5P0XdCq+3vdMV26dPFTFwAAoL0J+BJbeHh4gzNAMTExOnnypKKjo1VSUlJvru7jqKgo19JaY8dYl+YAAACaKqABqbCwUEOGDFF+fv3L+n/22Wfq37+/7Ha79u7dq5qaGtfc7t271bdvX0VGRio2NlYRERHas2ePa768vFz79+9XfDy3AwAAAN4JaECKiYnRj370I82fP1/5+fk6fPiwnnjiCX300UeaPn26UlJSdPbsWc2ZM0eHDh3S5s2btX79ek2bNk3Sd0t0qampWrp0qXbu3KnCwkI98MADio6OVlJSUiBbAwAAbVhA9yB16tRJq1ev1tKlSzVz5kyVl5dr4MCBev755zVgwABJ0tq1a7Vw4UIlJyerR48eyszMVHJysus50tPTVV1drblz56qqqkp2u13Z2dlcJBIAAHgt4BeKvOyyy/T444+7nR80aJBycnLczgcFBSkjI0MZGRktUR4AAOiAAr5JGwAAoLUhIAEAAFgQkAAAACwISAAAABYB36QNAL50ruq8qmuM23lj3M8BQB0CEoB2pbrGaNH6PLfzD0/gIrIALo4lNgAAAAsCEgAAgAVLbADaFPYYAfAHAhKANoU9RgD8gSU2AAAACwISAACABQEJAADAgoAEAABgQUACAACwICABAABYEJAAAAAsCEgAAAAWBCQAAAALrqQNwK8udquQ4CCbwsM6+7EiAGiIgATAry52q5DZE+1+rAYAGscSGwAAgAUBCQAAwIKABAAAYEFAAgAAsCAgAQAAWBCQAAAALAhIAAAAFgQkAAAACwISAACABQEJAADAgoAEAABgQUACAACwICABAABYEJAAAAAsCEgAAAAWBCQAAACL4EAXAADfZ7NJ5RVOt/PGGD9WA6CjIiABaFVqa40Wv5jvdv7hCfF+rAZAR8USGwAAgAUBCQAAwIIlNgDwo3NV51Vdc+F9VMFBNoWHdfZTRQAaQ0ACAD+qrjFatD7vgsfMnmj3UzUA3GGJDQAAwIKABAAAYEFAAgAAsCAgAQAAWBCQAAAALAhIAAAAFgQkAAAAC66DBAA+xM12gfaBgAQAPsTNdoH2gYAEAK3Mxc5CcSsSoOURkACglbnYWShuRQK0PDZpAwAAWLSqgHTkyBENGTJEmzdvdo0VFBQoNTVVcXFxGjVqlLKzs+s9pra2VitXrlRiYqIGDx6sKVOmqKioyN+lAwCAdqTVBKTz58/roYce0rlz51xjZWVlmjx5svr06aPc3FylpaVpxYoVys3NdR2zatUqbdy4UQsWLFBOTo5sNpumTp0qp9P9+j0AAMCFtJqA9Lvf/U4/+MEP6o29/PLLCgkJUVZWlvr166eUlBRNmjRJa9askSQ5nU6tW7dOaWlpGjlypGJjY7V8+XKdPHlS27dvD0QbAACgHWgVASkvL085OTlavHhxvfH8/HzZ7XYFB/99L3lCQoKOHDmi0tJSFRYWqqKiQgkJCa75bt26aeDAgcrLy/Nb/QAAoH0JeEAqLy9XZmam5s6dq169etWbKy4uVnR0dL2xnj17SpKOHz+u4uJiSWrwuJ49e+rEiRMtWDUAAGjPAv42/6ysLMXFxennP/95g7mqqiqFhITUGwsNDZUkORwOVVZWSlKjx5w5c6aFKgZwIeeqzqu6xv3VormSNIC2IKABacuWLcrPz9frr7/e6HxYWFiDzdYOh0OSFB4errCwMEnf7UWq+3vdMV26dGmhqgFcSHWN0aL17pe4uZI0gLYgoAEpNzdXpaWlGjVqVL3xefPmKTs7W1dccYVKSkrqzdV9HBUVperqatfYNddcU++Y2NjYli0eAAC0WwENSEuXLlVVVVW9sTFjxig9PV233nqr3njjDW3cuFE1NTUKCgqSJO3evVt9+/ZVZGSkunbtqoiICO3Zs8cVkMrLy7V//36lpqb6vR8AANA+BDQgRUVFNToeGRmpK6+8UikpKVq7dq3mzJmju+++W5988onWr1+v+fPnS/pu71FqaqqWLl2qyy67TFdeeaWWLFmi6OhoJSUl+bMVAADQjgR8k/aFREZGau3atVq4cKGSk5PVo0cPZWZmKjk52XVMenq6qqurNXfuXFVVVclutys7O7vBxm0AAICmanUB6cCBA/U+HjRokHJyctweHxQUpIyMDGVkZLR0aQAAoIMI+HWQAAAAWhsCEgAAgAUBCQAAwIKABAAAYEFAAgAAsGh172ID0LpxrzUAHQEBCYBHuNcagI6AJTYAAAALAhIAAIAFAQkAAMCCgAQAAGDBJm0AaGNsNqm8wul2PjjIpvCwzn6sCGh/CEgA0MbU1hotfjHf7fzsiXY/VgO0TyyxAQAAWBCQAAAALFhiA4B2hj1KQPMRkACgnWGPEtB8LLEBAABYEJAAAAAsCEgAAAAWBCQAAAALAhIAAIAFAQkAAMCCgAQAAGBBQAIAALAgIAEAAFgQkAAAACwISAAAABYEJAAAAAsCEgAAgAUBCQAAwIKABAAAYEFAAgAAsCAgAQAAWBCQAAAALAhIAAAAFgQkAAAAi+BAFwCgdTlXdV7VNcbtvDHu5wCgvSAgAainusZo0fo8t/MPT4j3YzUAEBgssQEAAFhwBgnoYFhCA4CLIyABHQxLaABwcSyxAQAAWBCQAAAALAhIAAAAFgQkAAAACzZpA0AHY7NJ5RVOt/PBQTaFh3X2Y0VA60NAAoAOprbWaPGL+W7nZ0+0+7EaoHViiQ0AAMCCgAQAAGBBQAIAALAgIAEAAFgQkAAAACwISAAAABZeBaS8vDxVVFQ0OldeXq433nijyc9VWlqqjIwMJSQkaMiQIbrnnnt06NAh13xBQYFSU1MVFxenUaNGKTs7u97ja2trtXLlSiUmJmrw4MGaMmWKioqKvGkLAABAkpcBacKECTp8+HCjc/v379cjjzzS5Of6zW9+o2PHjmnNmjX685//rLCwME2aNEmVlZUqKyvT5MmT1adPH+Xm5iotLU0rVqxQbm6u6/GrVq3Sxo0btWDBAuXk5Mhms2nq1KlyOt1fBA0AAOBCmnyhyIcfflgnTpyQJBljlJWVpYiIiAbHHT16VJdffnmTnrOsrExXXXWVfvOb3+hHP/qRJOnee+/VL37xCx08eFC7d+9WSEiIsrKyFBwcrH79+qmoqEhr1qxRSkqKnE6n1q1bp4yMDI0cOVKStHz5ciUmJmr79u0aO3ZsU9sDAABwafIZpH/6p3+SMUbGGNdY3cd1fzp16qS4uDg98cQTTXrO7t27a9myZa5w9PXXXys7O1vR0dHq37+/8vPzZbfbFRz89xyXkJCgI0eOqLS0VIWFhaqoqFBCQoJrvlu3bho4cKDy8vKa2hoAAEA9TT6DNHr0aI0ePVqSNH78eGVlZalfv34+K+Sxxx7Tyy+/rJCQED377LMKDw9XcXGxYmJi6h3Xs2dPSdLx48dVXFwsSerVq1eDY+rOdgEAAHjKqz1I//Ef/+HTcCRJEydOVG5urm677TbNmDFD+/btU1VVlUJCQuodFxoaKklyOByqrKyUpEaPcTgcPq0PAAB0HF7drLayslKrV6/Wrl27VFlZqdra2nrzNptNO3bs8Og5+/fvL0n67W9/q48++kgvvfSSwsLCGmy2rgs+4eHhCgsLkyQ5nU7X3+uO6dKli8d9AQAASF4GpIULFyo3N1fXX3+9fvzjH6tTJ+8up1RaWqrdu3frlltuUVBQkCSpU6dO6tevn0pKShQdHa2SkpJ6j6n7OCoqStXV1a6xa665pt4xsbGxXtUEAADgVUDatm2bHnjgAd1zzz3N+uQlJSV68MEHFRkZqeHDh0uSzp8/r/3792v06NG6/PLLtXHjRtXU1LgC1O7du9W3b19FRkaqa9euioiI0J49e1wBqby8XPv371dqamqzagMAAB2XV6d+qqurNWjQoGZ/8tjYWI0YMULz589Xfn6+Pv/8cz388MMqLy/XpEmTlJKSorNnz2rOnDk6dOiQNm/erPXr12vatGmSvtt7lJqaqqVLl2rnzp0qLCzUAw88oOjoaCUlJTW7PgAA0DF5dQZpxIgR+tvf/lbv7fXesNlsevrpp/XUU09p5syZ+vbbbxUfH68//vGPuuKKKyRJa9eu1cKFC5WcnKwePXooMzNTycnJrudIT09XdXW15s6dq6qqKtntdmVnZzfYuA0AANBUXgWkW2+9VfPmzdM333yjwYMHN7oh+vbbb2/Sc3Xt2lVZWVnKyspqdH7QoEHKyclx+/igoCBlZGQoIyOjSZ8PAADgYrwKSDNnzpQkbdmyRVu2bGkwb7PZmhyQAAAAWhuvAtLOnTt9XQcAAECr4VVAuvLKK31dBwAAQKvhVUB65plnLnrMfffd581TAwAABJzPA1JERIR69uxJQAIAAG2WVwGpsLCwwdi5c+e0d+9eZWVl6bHHHmt2YQAAAIHi3T1CGhEeHq7ExETNmDFDTz75pK+eFgAAwO98FpDq9OrVS4cPH/b10wIAAPiNV0tsjTHG6MSJE1qzZg3vcgMAAG2aVwEpNjZWNput0TljDEtsAACgTfMqIM2YMaPRgBQREaFRo0apT58+za0LAAAgYLwKSGlpab6uAwAAoNXweg+S0+nU5s2btWfPHpWXl6t79+6Kj49XcnKyQkNDfVkjAACAX3kVkMrLyzVhwgQVFhbqiiuuUI8ePXTkyBFt3bpVf/zjH7VhwwZ17drV17UCAAD4hVdv83/qqadUXFysl156SW+99ZZycnL01ltv6aWXXlJpaalWrFjh6zoBAAD8xquAtHPnTs2cOVPx8fH1xuPj45Wenq5t27b5pDgAAIBA8GqJraKiQldffXWjc1dffbVOnz7dnJoAXMC5qvOqrjFu54ODbAoP6+zHitDe2GxSeYXT7TzfY+gIvApIP/zhD7Vr1y7dcMMNDeZ27typ3r17N7swAI2rrjFatD7P7fzsiXY/VoP2qLbWaPGL+W7n+R5DR+BVQLrrrrs0a9YsOZ1O/fznP9fll1+ur7/+Wq+//ro2bdqkrKwsH5cJoKku9tu/Me7PPgEAvuNVQLr11lt19OhRrV69Wps2bXKNd+7cWTNmzNCvfvUrnxUIwDMX++3/4QnxbucAAN/xKiCdO3dO9957r1JTU/XRRx/pzJkzOnHihH71q1/pkksu8XWNAAAAfuXRu9gKCgp0++2364UXXpAkdevWTTfeeKNuvPFGPf3007rzzjt1+PDhlqgTAADAb5ockI4dO6ZJkybpzJkz6t+/f725kJAQPfroo6qoqNCdd96p4uJinxcKAADgL00OSM8995y6d++uV155RWPGjKk316VLF6Wmpio3N1fh4eFavXq1zwsFAADwlyYHpN27d+vuu+/WpZde6vaYyMhITZ48Wbt37/ZFbQAAAAHR5IB06tSpJl3fKCYmhiU2AADQpjU5IF122WUqKSm56HHffPPNBc8yAQAAtHZNDkh2u12bN2++6HFbtmzRj3/842YVBQAAEEhNDkjjx4/Xnj17tGjRIjkcjgbzTqdTixcv1jvvvKNx48b5tEgAAAB/avKFIn/yk5/okUce0eOPP65XX31Vw4cP11VXXaWamhodP35ce/bsUVlZme6//34lJia2ZM0AAAAtyqMraY8bN06xsbHKzs7Wzp07XWeSfvCDH2jEiBGaMmWKBg8e3CKFAgAA+IvHtxoZOnSohg4dKkkqKytTp06duL0IAABoV7y6F1ud7t27+6oOAACAVsOje7EBAAB0BAQkAAAACwISAACABQEJAADAgoAEAABgQUACAACwICABAABYEJAAAAAsCEgAAAAWBCQAAAALAhIAAIAFAQkAAMCCgAQAAGBBQAIAALAgIAEAAFgQkAAAACwISAAAABYEJAAAAAsCEgAAgEVwoAsAALQtNptUXuF0Ox8cZFN4WGc/VgT4HgEJAOCR2lqjxS/mu52fPdHux2qAlkFAAgD4FGeY0B4EPCCdPn1ay5Yt01//+ledPXtWAwYM0IMPPqj4+HhJUkFBgRYuXKjPPvtMl156qcaPH6+77rrL9fja2lo988wz2rRpk8rLyzV06FDNmzdPvXv3DlRLANChcYYJ7UHAN2nPmjVLH3/8sZYtW6Y///nP+od/+AfdddddOnz4sMrKyjR58mT16dNHubm5SktL04oVK5Sbm+t6/KpVq7Rx40YtWLBAOTk5stlsmjp1qpxO97+9AAAAXEhAzyAVFRXpf/7nf/SnP/1J1113nSRpzpw5+tvf/qatW7cqLCxMISEhysrKUnBwsPr166eioiKtWbNGKSkpcjqdWrdunTIyMjRy5EhJ0vLly5WYmKjt27dr7NixgWwPAAC0UQE9g9S9e3c999xzuvbaa11jNptNxhidOXNG+fn5stvtCg7+e45LSEjQkSNHVFpaqsLCQlVUVCghIcE1361bNw0cOFB5eXl+7QUAALQfAQ1I3bp108iRIxUSEuIa+8///E99+eWXGjFihIqLixUdHV3vMT179pQkHT9+XMXFxZKkXr16NTjmxIkTLVw9AABorwK+B+n79u7dq0cffVQ333yzRo8eraqqqnrhSZJCQ0MlSQ6HQ5WVlZLU6DEOh8M/RQMAgHan1QSkHTt26K677tKgQYO0bNkySVJYWFiDzdZ1wSc8PFxhYWGS1OgxXbp08UPVAACgPWoVAemll15SWlqabrzxRq1Zs8YVfKKjo1VSUlLv2LqPo6KiXEtrjR1jXZoDAABoqoAHpA0bNui3v/2txo0bp6effrrecpndbtfevXtVU1PjGtu9e7f69u2ryMhIxcbGKiIiQnv27HHNl5eXa//+/a7rKAEAAHgqoAHpyJEjevzxx5WUlKRp06aptLRUp06d0qlTp/Ttt98qJSVFZ8+e1Zw5c3To0CFt3rxZ69ev17Rp0yR9t/coNTVVS5cu1c6dO1VYWKgHHnhA0dHRSkpKCmRrAACgDQvodZD+8pe/6Pz589q+fbu2b99eby45OVmLFi3S2rVrtXDhQiUnJ6tHjx7KzMxUcnKy67j09HRVV1dr7ty5qqqqkt1uV3Z2doON2wAAAE0V0IA0ffp0TZ8+/YLHDBo0SDk5OW7ng4KClJGRoYyMDF+XBwAAOqiA34sNQH3nqs6rusa4nTfG/RwAwDcISEArU11jtGi9+yvBPzyBNyAAQEsL+LvYAAAAWhsCEgAAgAUBCQAAwIKABAAAYEFAAgAAsOBdbICf8TZ+AGj9CEiAn/E2fgBo/VhiAwAAsCAgAQAAWLDEBgDwK5tNKq9wup0PDrIpPKyzHysCGiIgAQD8qrbWaPGL+W7nZ0+0+7EaoHEssQEAAFgQkAAAACwISAAAABYEJAAAAAsCEgAAgAUBCQAAwIKABAAAYEFAAgAAsCAgAQAAWBCQAAAALAhIAAAAFgQkAAAACwISAACABQEJAADAgoAEAABgQUACAACwCA50AQAAfJ/NJpVXON3OBwfZFB7W2Y8VoSMiIAEAWpXaWqPFL+a7nZ890e7HatBRscQGAABgQUACAACwICABAABYEJAAAAAsCEgAAAAWBCQAAAALAhIAAIAFAQkAAMCCgAQAAGBBQAIAALAgIAEAAFgQkAAAACy4WS3gY+eqzqu6xridN8b9HACgdSAgAT5WXWO0aH2e2/mHJ8T7sRoAgDdYYgMAALAgIAEAAFgQkAAAACwISAAAABYEJAAAAAsCEgAAgAVv8wcAtCk2m1Re4XQ7HxxkU3hYZz9WhPaIgAQAaFNqa40Wv5jvdn72RLsfq0F7xRIbAACABQEJAADAolUFpFWrVmn8+PH1xgoKCpSamqq4uDiNGjVK2dnZ9eZra2u1cuVKJSYmavDgwZoyZYqKior8WTYAAGhnWk1AeuGFF7Ry5cp6Y2VlZZo8ebL69Omj3NxcpaWlacWKFcrNzXUds2rVKm3cuFELFixQTk6ObDabpk6dKqfT/QY+AACACwn4Ju2TJ09qzpw52rt3r/r27Vtv7uWXX1ZISIiysrIUHBysfv36qaioSGvWrFFKSoqcTqfWrVunjIwMjRw5UpK0fPlyJSYmavv27Ro7dmwgWgIAAG1cwM8g7du3T5dccolee+01DR48uN5cfn6+7Ha7goP/nuMSEhJ05MgRlZaWqrCwUBUVFUpISHDNd+vWTQMHDlRenvu7qQMAAFxIwM8gjR49WqNHj250rri4WDExMfXGevbsKUk6fvy4iouLJUm9evVqcMyJEydaoFoAANARBPwM0oVUVVUpJCSk3lhoaKgkyeFwqLKyUpIaPcbhcPinSAAA0O606oAUFhbWYLN1XfAJDw9XWFiYJDV6TJcuXfxTJAAAaHdadUCKjo5WSUlJvbG6j6OiolxLa40dEx0d7Z8iAQBAu9OqA5LdbtfevXtVU1PjGtu9e7f69u2ryMhIxcbGKiIiQnv27HHNl5eXa//+/YqPjw9EyQCAAKu7V5u7P+eqzge6RLQBAd+kfSEpKSlau3at5syZo7vvvluffPKJ1q9fr/nz50v6bu9Ramqqli5dqssuu0xXXnmllixZoujoaCUlJQW4egBAIHCvNvhCqw5IkZGRWrt2rRYuXKjk5GT16NFDmZmZSk5Odh2Tnp6u6upqzZ07V1VVVbLb7crOzm6wcRsAAKCpWlVAWrRoUYOxQYMGKScnx+1jgoKClJGRoYyMjJYsDXA5V3Ve1TXG7bwx7ucAAG1DqwpIQFtQXWO0aL37C5E+PIH9bwDQ1rXqTdoAAACBQEACAACwICABAABYEJAAAAAsCEgAAAAWBCQAAAAL3uYPfM/FrnEkcZ0jAOgICEjA91zsGkcS1zkCgI6AJTYAAAALziChQ+E2IQCApiAgoUPhNiEAgKZgiQ0AAMCCgAQAAGBBQAIAALAgIAEAAFgQkAAAACwISAAAABYEJAAAAAuugwQA6FBsNqm8wul2PjjIpvCwzn6sCK0RAQkA0KHU1hotfjHf7fzsiXY/VoPWiiU2AAAACwISAACABQEJAADAgoAEAABgQUACAACwICABAABYEJAAAAAsuA4SAADfw4UkIRGQAACohwtJQmKJDQAAoAHOIKFdOVd1XtU1xu28Me7nAACoQ0BCm9KUAHShU+MPT4hvibIAAO0MAQltSnWN0aL1eW7nCUAAAF9gDxIAAIAFAQkAAMCCgAQAAGBBQAIAALBgkzYAAB7gStsdAwEJAAAPcKXtjoElNgAAAAsCEgAAgAUBCQAAwIKABAAAYMEmbQAA/Ohi95TkXXCtAwEJAAAfuthlAC52U23eBdc6EJAAAPChi10GgJtqtw3sQQIAALDgDBJalYutzRvjfg4AAF8hIKFVqa4xWrQ+z+08p6YBtHfcyqR1ICABANCKcCuT1oE9SAAAABacQYJfsccIANAWEJDgV+wxAoDmudgepU42qfYCv2uyh6lp2kVAqq2t1TPPPKNNmzapvLxcQ4cO1bx589S7d+9Al9bhcIYIAFpWU66zxB6m5msXAWnVqlXauHGjnnjiCUVFRWnJkiWaOnWqtm7dqpCQkECX1640JQBxgTQAQFvX5gOS0+nUunXrlJGRoZEjR0qSli9frsTERG3fvl1jx44NcIXtC0tkAICOoM0HpMLCQlVUVCghIcE11q1bNw0cOFB5eXkEJIuLnQG62No1S2QA0LZdbA+TdPGfBR1hn5PNtPGfeNu2bVNaWpo+/vhjhYWFucbvv/9+VVVV6Q9/+INHz/eTn/xENTU16tWrl69LlSQZIxm5/5LbZJPN1iKfWpJUa4xOf+twO39p11DmLzDfGmpgvn3Pt4YamO/Y8754ju5dQy/wk65lftadOHFCQUFB+vTTT33yfG0+IL366qvKzMxUQUGBOnX6+2WdMjMzVVJSohdeeMGj54uPj5fT6VSPHj18XCkAAGgpp06dUkhIiPLz3e+D9USbX2KrO2vkdDrrnUFyOBzq0qWLx8/nqy8sAABou9r8lbTrlsJKSkrqjZeUlCg6OjoQJQEAgDauzQek2NhYRUREaM+ePa6x8vJy7d+/X/HxvKMKAAB4rs0vsYWEhCg1NVVLly7VZZddpiuvvFJLlixRdHS0kpKSAl0eAABog9p8QJKk9PR0VVdXa+7cuaqqqpLdbld2djYXiQQAAF5p8+9iAwAA8LU2vwcJAADA1whIAAAAFgQkAAAACwISAACABQEJAADAgoAEAABgQUACAACwICB5wOFwaP78+Ro+fLiGDBmi9PR0lZaWXvAxR48e1T333KP4+HiNGDFCc+fOVXl5uZ8qbj5vej579qzmzZunhIQEDR06VNOnT9exY8f8VHHzeNPv9z377LMaMGBAC1boW970+8EHH2j8+PEaOnSoEhMTNWfOHJ0+fdo/BXuhtrZWK1euVGJiogYPHqwpU6aoqKjI7fFlZWV68MEHZbfbZbfb9dhjj+ncuXN+rLh5PO334MGDuueeezRs2DANHz5c6enpOn78uB8rbj5Pe/6+119/XQMGDNBXX33VwlX6jqf9nj9/Xk899ZQSExMVFxen1NRUFRQU+LHi5vO051OnTmnWrFkaNmyYhg0bpvvvv1/FxcWefVKDJps9e7ZJSkoyeXl55uOPPza33367GTdunNvjnU6nSUpKMvfee685dOiQycvLM2PGjDFpaWl+rLp5PO3ZGGMmTZpkkpKSzPvvv28KCgrMHXfcYcaOHWtqamr8VLX3vOm3zscff2wGDhxoYmJiWrhK3/G03y+++MLExcWZBQsWmMOHD5u8vDzzs5/9zIwfP96PVXvmd7/7nRk+fLj561//agoKCsyUKVNMUlKScTgcjR6fmppq/vVf/9V89tln5t133zU33XSTyczM9HPV3vOk32+++cbccMMNZubMmebzzz83n376qUlNTTW33HKLqaqqCkD13vH0Na7z1VdfmaFDh5qYmBhz7NgxP1XbfJ72++ijj5qEhASza9cuc+jQITNjxgxzww03mPLycj9X7j1Pex43bpy54447zL59+8y+ffvML3/5S5OcnOzR5yQgNVFxcbGJjY01b7/9tmvsiy++MDExMebDDz9s9DH79+83MTExprCw0DW2fv16ExcX19Ll+oQ3Pb/33ntmwIAB9Xo+ePCgGTVqlDl8+HBLl9ws3vRbp6KiwowZM8ZMmDChzQQkb/pdtmyZGTNmjKmtrXWN5eXlmZiYGPPll1+2dMkeczgcZsiQIWbDhg2usTNnzphBgwaZrVu3Njj+gw8+MDExMebQoUOusXfeeccMGDDAFBcX+6Xm5vC035dfftlcd9119cLQiRMnTExMjHn33Xf9UnNzedpznZqaGvPrX//a9W+2rQQkT/v98ssvTUxMjNm1a1e942+66aZ2+xqfOXPGxMTEmJ07d7rGduzYYWJiYsw333zT5M/LElsT7d27V5I0bNgw11jfvn0VFRWlvLy8Rh9zySWXyGazadOmTXI6nfrmm2/0l7/8RXFxcf4oudm86fmdd95RTExMvWWm/v37a9euXfrhD3/YsgU3kzf91lm4cKFiYmL0i1/8okVr9CVv+r3tttu0ePFi2Wy2BnOtcZmtsLBQFRUVSkhIcI1169ZNAwcObLTH/Px89ejRQ/369XONXX/99bLZbK6vV2vmab/Dhw/X73//e4WGhjaYO3PmTIvW6iue9lxn9erVOn/+vKZNm+aPMn3G037/+7//W926ddONN95Y7/i33npLw4cP90vNzeVpz6GhoQoPD9eWLVt09uxZnT17Vq+++qr69OmjSy65pMmfl4DURCdPnlT37t0b/EfSs2dPnThxotHHXHHFFXr00UeVm5urwYMHa/jw4Tpz5oyWLVvmj5KbzZuejx49qt69e2vDhg0aO3asEhMTNXPmTJ08edIfJTeLN/1K0vbt2/X222/r3//931u6RJ/ypt9+/fo1CPhr1qxRjx49FBsb21Kleq1uz0GvXr3qjbvr8eTJkw2ODQkJ0aWXXnrB74HWwtN+r7rqqno/dCTpD3/4g0JDQ2W321uuUB/ytGdJ+uSTT7Ru3TotWbJEQUFBLV6jL3na79GjR3X11Vdr27Zt+pd/+RfdcMMNmjp1qg4fPuyXen3B055DQ0O1cOFCvf/++4qPj5fdbtdHH32kNWvWqFOnpsee4OaV3X589dVXuvnmm93O33///QoJCWkwHhoaKofD0ehjnE6nDh48qDFjxmjcuHEqKyvTk08+qQceeEDr1q3z6IVqCS3R89mzZ7Vv3z6VlZVp/vz5kqSlS5dqwoQJeu211xr9TdVfWqLfkydP6rHHHtOTTz6p7t27+6xWX2iJfq0WLVqkt99+WytXrlTnzp29rrWlVFZWSlKDPkNDQxs9Q1JZWdnsr0kgedqv1YsvvqgNGzbokUceUWRkZIvU6Gue9nzu3Dk99NBDeuihh9SnT5828cvb93na79mzZ/Xll19q1apVyszMVLdu3fTss8/qzjvv1JtvvtkmXmdPezbG6MCBAxoyZIjuvvtu1dTUaPny5ZoxY4b+9Kc/KSIiokmfl4D0/6KiovTmm2+6nX/77bfldDobjDscDnXp0qXRxzz//PN6//339eabb7p+S+nTp4/GjBmjt956Sz/96U99U7yXWqLnzp07y+Fw6Pe//73rVOYzzzyjxMREvfXWW7rlllt8U7wXfN2vMUazZ8/WLbfcUu/0dWvREq9vnfPnz+vf/u3f9Morr2jevHkaM2ZMs+ttCWFhYZK++2Wl7u+S+x7DwsLcfk3Cw8NbrlAf8bTfOsYYrVixQs8++6ymTZumSZMmtXSpPuNpzwsWLFCfPn10xx13+K1GX/K0386dO+vbb7/V8uXLXUvHy5cv18iRI/XKK6/o7rvv9k/hzeBpz2+88YY2bNigXbt2ucLQ6tWrddNNNyk3N1cTJ05s0uclIP2/zp0719t3YHXgwAGdPn1aTqezXootKSlRdHR0o4/Zu3evBg4cWO8Ubu/evdW9e3cdPXrUZ7V7qyV6jo6OVlRUVL113ssvv1yXXnppwN9G6+t+jx8/rnfffVcffPCBtmzZIkmqrq6WJA0ZMkTTpk3T9OnTfduEB1ri9ZW++430vvvuU35+vp566imNHTvWp3X7Ut0p+ZKSEl1zzTWu8ZKSkkaXBKOjo7Vjx456Y06nU6dPn1ZUVFTLFusDnvYrfRd2H3nkEW3dulWZmZm66667/FKrr3jac25urkJCQjRkyBBJUk1NjSTpZz/7mW677bZWv1Tuzfd0cHBwvf8LwsLCdPXVVwf8/+Sm8rTnvXv3qm/fvvXOFF1yySXq27evRz972YPUREOHDlVtbW29jZpffPGFTp48qfj4+EYfEx0drYMHD8oY4xo7efKkTp8+rT59+rR0yc3mTc/x8fE6fvy4SkpKXGMlJSUqKytT7969W7zm5vC036ioKG3btk2vvfaatmzZoi1btig9PV2StGXLllb/G6o3r6/T6dS0adP06aefau3ata06HElSbGysIiIitGfPHtdYeXm59u/f32iPdrtdxcXF9a6vUvfY6667ruULbiZP+5WkzMxM/dd//ZeeeuqpNheOJM973rZtm7Zu3er6N7tgwQJJ0nPPPaf777/fb3V7y9N+4+PjVV1drU8//dQ1VlVVpWPHjrX6/5PreNpzr169VFRUVG9ZvLKyUl999ZVnPXv/xruOZ9asWWb06NHmvffec10zJjU11TXvcDhMSUmJ67oMBw4cMNdee62ZM2eOOXTokPnwww/NL3/5S3PbbbcZp9MZqDY84mnPDofD3HrrreaOO+4wn3zyidm3b58ZN26c+ed//ueLXpOkNfC0X6vc3Nw28zZ/Yzzvd+XKlWbAgAFm69atpqSkpN6f1vr6Llu2zFx//fVmx44druunjBkzxjgcDlNdXW1KSkpMZWWlMcaY2tpac8cdd5jk5GTz8ccfm927d5ubbrrJzJ49O8BdNJ0n/dZ9v65du7bB61l3TFvgSc9W7733Xpt6m78xnvc7adIkc8stt5i8vDxz8OBBk5aWZoYPH25KS0sD2IVnPOn55MmT5vrrrzfTp083BQUFpqCgwEybNs2MGDHCnDlzpsmfk4DkgYqKCjNnzhwTHx9v4uPjzaxZs+pdU6HuH9p7773nGvvwww/N+PHjzdChQ82IESPM7Nmz29Q3pTc9l5SUmFmzZpnrrrvODBkyxNx3331t4hoyxnjX7/e1tYDkab9jxowxMTExjf5x9zUJtOrqavPkk0+ahIQEExcXZ6ZOner6YXjs2DETExNjcnNzXcd//fXXJi0tzcTFxZlhw4aZefPmtamLJnrS7+TJk92+nt//mrR2nr7G39cWA5Kn/X777bdm3rx5ZtiwYWbw4MFm8uTJ5uDBg4Eq3yue9nzo0CEzbdo0c/3115uEhARz3333efwa24z53voPAAAA2IMEAABgRUACAACwICABAABYEJAAAAAsCEgAAAAWBCQAAAALAhIAAIAFAQkAAMCCgAQAAGBBQAIAALAgIAEAAFj8H+/2NxtPMg3kAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sns.histplot(results.null_distribution)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "4b6633d8-fe62-449b-9a26-1863bd3b1ea8", "metadata": {}, "source": [ "We can see that the null distribution is nice and symmetrical (the average correlation in shuffled data is zero, whilst positive and negative correlations are equally likely)\n", "\n", "Note that our observed correlation, 0.8, is so rare that it was *never* obtained in the 10,000 shuffles" ] }, { "cell_type": "code", "execution_count": 8, "id": "ddddc13f-7a62-44a7-975e-36d09b55fd39", "metadata": { "tags": [] }, "outputs": [ { "data": { "text/plain": [ "10000" ] }, "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 >= 0.801)" ] }, { "cell_type": "markdown", "id": "0f6dc507", "metadata": {}, "source": [ "### The $𝑝$-valueΒΆ\n", "\n", "The probability that the test statistic (in this case, the correlation between broccoli consumption and IQ) 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 basically zero- we never get a correlation of 0.80 in our 10,000 random shuffles.\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, so the result is significant if $p=0.05$ or $p=0.01$. Our result would be considered highly statistically significant." ] }, { "cell_type": "markdown", "id": "604bf628-6d76-49e6-8a3e-eeb5abea174c", "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", "\n", "Firstly, we had to give the function `stats.permutation_test()` our two samples. In this case, the data are in 'wide form', ie the two measurements for each individual (broccoli consumption and IW) are contained in two different columns (broccoli.broccoli_g, broccoli.IQ)\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, `correlation()`, using a function definition, before running the permutation test\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 correlation so we shuffle the *pairings*, ie in shuffled data we pair each person's IQ with a random person's broccoli consumption" ] }, { "cell_type": "markdown", "id": "5fca8a34", "metadata": {}, "source": [ "## Recap\n", "\n", "To run a permutation test on a correlation, we shuffled up all the pairings so each person's IQ was paired with someone else's broccoli consumption. We did not switch any datapoints from the broccoli column into the the IQ column (!).\n", "\n", "For each shuffle we calculated the correlation between broccoli consumption and IQ\n", "\n", "Permutation testing in this way gives us a null distribution for the correlation. Values of the correlation coefficient that occur rarely in the null distriubtion are considered statistically significant.\n", " \n", "To run the permutation test with scipy.stats we need the option `permutation_type='pairings'`" ] }, { "cell_type": "code", "execution_count": null, "id": "267cc3b2", "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 }