{
"cells": [
{
"cell_type": "markdown",
"id": "c2bf7028-42c0-4310-8bf8-e7e5a9bf27fd",
"metadata": {},
"source": [
"## Permutation Tests\n",
"### Exact Tests\n",
"Consider the following experiment from Efron and Tibshirani's [An Introduction to the Bootstrap](https://books.google.com/books?id=MWC1DwAAQBAJ&printsec=frontcoverhttps://books.google.com/books?id=MWC1DwAAQBAJ&printsec=frontcover). A new medical treatment is intended to prolong life after a form of surgery. Sixteen mice are randomly assigned to either a treatment group or control group under the constraint that only seven treatments are available. All mice receive the surgery, but only the treatment group will receive the treatment being studied. The survival time of each mouse after surgery is recorded below."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "823b24e2-d26a-4f79-94cb-cc1eb27c8e5b",
"metadata": {},
"outputs": [],
"source": [
"# survival times measured in days\n",
"import numpy as np\n",
"x = np.array([94, 197, 16, 38, 99, 141, 23]) # treatment group\n",
"y = np.array([52, 104, 146, 10, 51, 30, 40, 27, 46]) # control group"
]
},
{
"cell_type": "markdown",
"id": "e19f78a9-e1c0-4d8a-aeb6-3fa944752948",
"metadata": {},
"source": [
"The difference in mean lifetime after treatment between the two groups suggests that the treatment has a prolonging effect, as hypothesized."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "5b3b1ead-55f7-456f-b015-9e36d0380c9d",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"30.63492063492064\n"
]
}
],
"source": [
"def statistic(x, y):\n",
" return np.mean(x) - np.mean(y)\n",
"print(statistic(x, y))"
]
},
{
"cell_type": "markdown",
"id": "9d016372-3ebe-49e5-a70c-923d8e9a1650",
"metadata": {},
"source": [
"It is possible that the treatment has no effect in reality; perhaps the apparent prolonging effect of the treatment is due to the inherent variability in survival times and chance alone. This possibility is typically assessed using Student's t-test. A common formulation of the test begins with the null hypothesis that the survival times `x` and `y` are sampled at random from normal distributions $X$ and $Y$ with means $\\mu_x$ and $\\mu_y$ and a common standard deviation, $\\sigma$. To test the null hypothesis that $\\mu_x = \\mu_y$ against the alternative that $\\mu_x > \\mu_y$, we perform the independent sample t-test with `stats.ttest_ind`."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "b993796d-05c7-4b75-bc70-7cd974da76b0",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Ttest_indResult(statistic=1.1208453991208167, pvalue=0.14060629239765005)"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from scipy import stats\n",
"stats.ttest_ind(x, y, equal_var=True, alternative='greater')"
]
},
{
"cell_type": "markdown",
"id": "33e2913b-ac57-493a-a6c8-20f11ab9b933",
"metadata": {},
"source": [
"The probability of observing such an extreme test statistic under the null hypothesis (due to chance alone) is greater than 14%, so these data do not seem inconsistent with the null hypothesis. The *point estimate* of the statistic (~30 days) suggested a life-prolonging effect, but such a value of the statistic could quite easily have been observed due to chance alone.\n",
"\n",
"Although the t-test tends to be rather robust to violations of its underlying assumptions (e.g., $X$ and $Y$ do not need to be strictly normally distributed for the test to be reasonably accurate), it is possible to perform a hypothesis test which requires almost no such assumptions at all. \n",
"\n",
"Instead, let the null hypothesis be that the observations from the samples `x` and `y` were all drawn independently[**†**](#cite_note-2) from a single distribution ($X = Y = Z$), and test this against the alternative that the two samples were drawn from distinct distributions that would tend to produce a greater value of `statistic` (in this case such that $\\mu_x > \\mu_y$).\n",
"\n",
"[**†**](#cite_note-2) Actually, only exchangeability is required [[4]](https://en.wikipedia.org/wiki/Exchangeable_random_variables)."
]
},
{
"cell_type": "markdown",
"id": "96079705",
"metadata": {},
"source": [
"The complete population of mice survival times in the study is really:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "8888960d-bb51-4549-a01c-d1e70f880744",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ 94 197 16 38 99 141 23 52 104 146 10 51 30 40 27 46]\n"
]
}
],
"source": [
"z = np.concatenate([x, y])\n",
"print(z)"
]
},
{
"cell_type": "markdown",
"id": "6204af8f-b858-4d5a-8bed-b4ef1bc00c30",
"metadata": {},
"source": [
"Since the mice were randomly divided into the two groups under the constraint that there were only seven treatments available, any selection of seven mice from `z` to form the treatment group `x` was equally likely; the remaining mice would form the control group `y`. Furthermore, if the null hypothesis is true, the mice survival times would be *unaffected by the grouping*. Therefore, each value of the statistic obtained from the possible groupings is equally likely.\n",
"\n",
"We begin our hypothesis test by calculating the value of `statistic` for all possible *permutations*[**‡**](#cite_note-3) of mice into the the two groups, forming an exact null distribution.\n",
"\n",
"[**‡**](#cite_note-3) Here and below, we will refer to the the ways of rearranging samples as \"permutations\" even when the word is not stricly appropriate in the technical sense. "
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "c7fed236-d06e-478b-b580-cbdd468730eb",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0, 0.5, 'Observed Frequency')"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEGCAYAAABy53LJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAd4ElEQVR4nO3dfZxdVX3v8c+XYHiw0EgTaMxDJ9IATX0hhCGhCFWeNEEkYHmUXihY0yjprXq1BPFay+vevlBQbimYGErAeIUAlofUhgJigSvXQAaEkACRIYAMiRDsLU+hYPB3/9hrYOfkPOw9zJ45M/N9v17nNXuvvdY6v3OYzI+99t5rKSIwMzMrarvBDsDMzIYWJw4zMyvFicPMzEpx4jAzs1KcOMzMrJTtBzuAgTB27Njo6OgY7DDMzIaU+++//4WIGFdbPiISR0dHB11dXYMdhpnZkCLp6XrlHqoyM7NSnDjMzKwUJw4zMyvFicPMzEpx4jAzs1KcOMzMrBQnDjMzK8WJw8zMSnHiMDOzUkbEk+Nm/a1jwb/ULX/qgo8NcCRmA89nHGZmVooTh5mZlVJp4pA0S9I6Sd2SFtQ5LkmXpOOrJU3PHVsi6XlJaxr0/UVJIWlslZ/BzMy2VlnikDQKuAyYDUwDTpU0rababGBqes0FFuaOXQXMatD3JOAo4Bf9G7WZmbVS5cXxGUB3RKwHkLQMmAM8kqszB1gaEQGslDRG0viI2BgRd0vqaND3xcBfAzdXF75Zeb5obiNBlUNVE4Bncvs9qaxsna1IOhZ4NiIealFvrqQuSV2bNm0qHrWZmTVVZeJQnbLoQ523K0s7A+cBX2315hGxOCI6I6Jz3LhtFrAyM7M+qnKoqgeYlNufCGzoQ528PYEpwEOSeus/IGlGRPzyHUdsVqPR0JPZSFblGccqYKqkKZJGA6cAy2vqLAdOT3dXHQS8GBEbG3UYEQ9HxO4R0RERHWSJZ7qThpnZwKkscUTEFmA+cCvwKHBdRKyVNE/SvFRtBbAe6AYuBz7b217SNcBPgb0l9Uj6VFWxmplZcZVOORIRK8iSQ75sUW47gLMbtD21QP8d7zBEMzMryU+Om5lZKU4cZmZWihOHmZmV4sRhZmalOHGYmVkpXsjJbAB4DisbTnzGYWZmpThxmJlZKU4cZmZWiq9xmA0iX/uwochnHGZmVooTh5mZleLEYWZmpfgahxlesMmsDJ9xmJlZKU4cZmZWihOHmZmV4sRhZmalOHGYmVkpThxmZlZKpYlD0ixJ6yR1S1pQ57gkXZKOr5Y0PXdsiaTnJa2paXOhpMdS/RsljanyM5iZ2dYqe45D0ijgMuAooAdYJWl5RDySqzYbmJpeM4GF6SfAVcClwNKarm8Hzo2ILZK+DpwLnFPV57Dhxc9rmL1zVT4AOAPojoj1AJKWAXOAfOKYAyyNiABWShojaXxEbIyIuyV11HYaEbfldlcCJ1T2CcwGSbME5wkQbbBVOVQ1AXgmt9+TysrWaeYs4JZ6ByTNldQlqWvTpk0lujQzs2aqTByqUxZ9qFO/c+k8YAvw/XrHI2JxRHRGROe4ceOKdGlmZgVUOVTVA0zK7U8ENvShzjYknQEcAxyRhrnMzGyAVHnGsQqYKmmKpNHAKcDymjrLgdPT3VUHAS9GxMZmnUqaRXYx/NiI2FxF4GZm1lhlZxzprqf5wK3AKGBJRKyVNC8dXwSsAI4GuoHNwJm97SVdA3wYGCupB/ibiLiC7E6rHYDbJQGsjIh5VX0Os3bjVQNtsGkkjPR0dnZGV1fXYIdh74Bvo+07JxTrK0n3R0RnbbmfHDczs1KcOMzMrBSvAGg2zPmaiPU3n3GYmVkpThxmZlaKE4eZmZXixGFmZqW0TByS3j8QgZiZ2dBQ5IxjkaT7JH3WiyaZmVnLxBERhwCnkU1G2CXpaklHVR6ZmZm1pULXOCLiceArZJMLfgi4JC3f+okqgzMzs/ZT5BrHvpIuBh4FDgc+HhF/kLYvrjg+MzNrM0WeHL8UuBz4ckS81lsYERskfaWyyMzMrC0VSRxHA69FxJsAkrYDdoyIzRHxvUqjMzOztlMkcfwIOBJ4Je3vDNwGHFxVUGZWPc9hZX1V5OL4jhHRmzRI2ztXF5KZmbWzIonjVUnTe3ckHQC81qS+mZkNY0WGqj4HXC9pQ9ofD5xcWUQ2YnmVP7OhoWXiiIhVkvYB9gYEPBYRv648MjMza0tFJzk8ENgX2B84VdLpRRpJmiVpnaRuSQvqHJekS9Lx1TVDYkskPS9pTU2b3STdLunx9PM9BT+DmZn1gyIPAH4PuAg4hCyBHAhss3h5nXajgMuA2cA0soQzrababGBqes0FFuaOXQXMqtP1AuCOiJgK3JH2zcxsgBS5xtEJTIuIKNn3DKA7ItYDSFoGzAEeydWZAyxNfa+UNEbS+IjYGBF3S+qo0+8c4MNp+7vAnWRToZiZ2QAoMlS1BvjdPvQ9AXgmt9+TysrWqbVHRGwESD9370NsZmbWR0XOOMYCj0i6D3i9tzAijm3RTnXKas9aitTpE0lzyYa/mDx5cn90aWZmFEscX+tj3z1kU7H3mghs6EOdWs/1DmdJGg88X69SRCwGFgN0dnb2SzIyM7Ni63HcBTwFvCttrwIeKND3KmCqpCmSRgOnAMtr6iwHTk93Vx0EvNg7DNXEcuCMtH0GcHOBWMzMrJ8Uuavq08APgO+kognATa3aRcQWYD5wK9mU7NdFxFpJ8yTNS9VWAOuBbrIZeD+be99rgJ8Ce0vqkfSpdOgC4ChJjwNHpX0zMxsgRYaqzia7Q+peyBZ1klTognRErCBLDvmyRbntSP3Xa3tqg/JfAUcUeX8zM+t/Re6qej0i3ujdkbQ9/XQB28zMhp4iieMuSV8GdkprjV8P/HO1YZmZWbsqkjgWAJuAh4G/IBt68sp/ZmYjVJFJDn9DduH68urDMbPB5gWerJWWiUPSk9S5phER76skIjMza2tF56rqtSNwIrBbNeGYmVm7K/IA4K9yr2cj4n8Bh1cfmpmZtaMiQ1XTc7vbkZ2B7FJZRGZm1taKDFV9M7e9hWz6kZMqicbMzNpekbuqDhuIQMzMbGgoMlT1hWbHI+Jb/ReOmZm1u6J3VR3I2zPbfhy4m60XYDIzsxGi6EJO0yPiZQBJXwOuj4g/rzIwG74aPWBmZkNDkSlHJgNv5PbfADoqicbMzNpekTOO7wH3SbqR7Any44GllUZlZmZtq8hdVf9T0i3AoanozIj4WbVhmZlZuyoyVAWwM/BSRPw90CNpSoUxmZlZGyuydOzfAOcA56aidwH/u8qgzMysfRU54zgeOBZ4FSAiNuApR8zMRqwiieONtDZ4AEh6d7UhmZlZOyuSOK6T9B1gjKRPAz+i4KJOkmZJWiepW9KCOscl6ZJ0fHV+QsVGbSXtJ2mlpAcldUmaUSQWMzPrH03vqpIk4FpgH+AlYG/gqxFxe6uOJY0CLgOOAnqAVZKWR8QjuWqzganpNRNYCMxs0fYbwN9GxC2Sjk77Hy7+kc3M7J1omjgiIiTdFBEHAC2TRY0ZQHdErAeQtAyYA+QTxxxgaRoKWylpjKTxZA8YNmobwK6p/W8DG0rGZWZm70CRBwBXSjowIlaV7HsCW89n1UN2VtGqzoQWbT8H3CrpIrKhtoPrvbmkucBcgMmTJ5cM3czMGilyjeMwsuTxRLoO8bCk1QXaqU5Z7drljeo0a/sZ4PMRMQn4PHBFvTePiMUR0RkRnePGjSsQrpmZFdHwjEPS5Ij4Bdl1iL7oASbl9iey7bBSozqjm7Q9A/irtH098I99jM/MzPqg2RnHTQAR8TTwrYh4Ov8q0PcqYKqkKZJGA6fw9tTsvZYDp6e7qw4CXoyIjS3abgA+lLYPBx4vEIuZmfWTZtc48sNF7yvbcURskTQfuBUYBSyJiLWS5qXji4AVwNFAN7AZOLNZ29T1p4G/l7Q98J+k6xhmZjYwmiWOaLBdWESsIEsO+bJFue0Azi7aNpX/BDigL/GYWd81WkflqQs+NsCR2GBrljg+IOklsjOPndI2aT8iYtfGTc3MbLhqmDgiYtRABmLDj1f6Mxueik6rbmZmBjhxmJlZSU4cZmZWihOHmZmV0uzJ8Zdpchuu76oyMxuZmt1VtQuApPOBXwLfI7sV9zS8AqCZ2YhVZKjqoxHx7Yh4OSJeioiFwJ9UHZiZmbWnIonjTUmnSRolaTtJpwFvVh2YmZm1pyKJ45PAScBz6XViKjMzsxGo5UJOEfEU2ep7Zmbb8BxWI0/LMw5Je0m6Q9KatL+vpK9UH5qZmbWjIkNVlwPnAr8GiIjVZOtjmJnZCFQkcewcEffVlG2pIhgzM2t/RRLHC5L2JD0MKOkEYGOlUZmZWdtqeXGcbKGlxcA+kp4FniR7CNDMzEagIonj6Yg4UtK7ge0i4uWqg7KhxetumI0sRYaqnpS0GDgIeKXieMzMrM0VOePYG/g42ZDVFZJ+CCxLa3+bmdXl5zuGr5ZnHBHxWkRcFxGfAPYHdgXuKtK5pFmS1knqlrSgznFJuiQdXy1pepG2kv4yHVsr6RtFYjEzs/5R5IwDSR8CTgZmA6vIpiBp1WYUcBlwFNADrJK0PCIeyVWbDUxNr5nAQmBms7aSDiN7kn3fiHhd0u7FPqqZmfWHlolD0pPAg8B1wJci4tWCfc8AuiNifepnGdkf/HzimAMsjYgAVkoaI2k80NGk7WeACyLidYCIeL5gPGZm1g+aDlWl//O/MiKOj4hrSiQNgAnAM7n9nlRWpE6ztnsBh0q6V9Jdkg5sEPtcSV2SujZt2lQibDMza6Zp4oiIN4HD+ti36nVZsE6zttsD7yG7y+tLwHWStqkfEYsjojMiOseNG1c8ajMza6rINY7/K+lS4FrgrTOOiHigRbseYFJufyKwoWCd0U3a9gA3pOGt+yT9BhgL+LTCzGwAFEkcB6ef5+fKAji8RbtVwFRJU4BnySZGrF3HYzkwP13DmAm8GBEbJW1q0vam9N53StqLLMm8UOBzmJlZPyiyHkefhqoiYouk+cCtwChgSUSslTQvHV8ErACOBrqBzcCZzdqmrpcAS9I0728AZ6SzDzMzGwBF7qraA/g74L0RMVvSNOCPIuKKVm0jYgVZcsiXLcptB9mDhYXapvI3gD9t9d5mZlaNIkNVVwFXAuel/Z+TXe9omThsePGcVNYf/ET50FdkrqqxEXEd8BvIhpGANyuNyszM2laRxPGqpN/h7fU4DgJerDQqMzNrW0WGqr5AdvfTnpLuAcYBJ1QalZmZta0id1U9kOaq2pvswbx1EfHryiMzM7O21HKoStKJwE7pdtjjgGvzs9iamdnIUuQax3+PiJclHQJ8FPgu2Sy2ZmY2AhVJHL13UH0MWBgRN5M9rW1mZiNQkcTxrKTvkK3BsULSDgXbmZnZMFTkrqqTgFnARRHxH2m9jC9VG5YNJj/oZ2bNFLmrarOkp4DZkmYB90TEbZVHZmYjip8oHzqK3FX1VbIL4r9DNn35lZK+UnVgZmbWnooMVZ0K7B8R/wkg6QLgAeB/VBmYmZm1pyIXuZ8Cdszt7wA8UUk0ZmbW9hqecUj6B7L5qV4H1kq6Pe0fBfxkYMIzM7N202yoqiv9vB+4MVd+Z2XRmJlZ22uYOCLiuwCSdgR+n+xs44neax1mZjYyNRuq2p5s5b+zgKfJrodMlHQlcJ4nOhz6/LyGmfVFs4vjFwK7AVMi4oCI2B/YExgDXDQAsZmZWRtqljiOAT4dES/3FkTES8BngKOLdC5plqR1krolLahzXJIuScdX52fdLdD2i5JC0tgisZiZWf9oljgiIqJO4Zuk1QCbkTQKuAyYDUwDTpU0rababGBqes0lzbrbqq2kSWR3d/2iVRxmZta/mt1V9Yik0yNiab5Q0p8CjxXoewbQHRHrU7tlwBzgkVydOcDSlKBWShqT5sLqaNH2YuCvgZsLxGFmQ5inImk/zRLH2cANks4iuyU3gAOBnYDjC/Q9AXgmt98DzCxQZ0KztpKOBZ6NiIckNXxzSXPJzmKYPHlygXDNzKyIZrfjPgvMlHQ48Idky8beEhF3FOy73l/12iGuRnXqlkvaGTgP+EirN4+IxcBigM7OzpZDa2ZmVkyR2XF/DPy4D333AJNy+xOBDQXrjG5QvicwBeg925gIPCBpRkT8sg8xmplZSVUuyLQKmCppiqTRwCnA8po6y4HT091VBwEvRsTGRm0j4uGI2D0iOiKigyzxTHfSMDMbOEVmx+2TiNgiaT5wKzAKWBIRayXNS8cXASvIbu3tBjYDZzZrW1WsZmZWXGWJAyAiVpAlh3zZotx2kF2EL9S2Tp2Odx7l8OcnxM2sP3ntcDMzK8WJw8zMSnHiMDOzUpw4zMysFCcOMzMrxYnDzMxKceIwM7NSKn2Ow8ysKs2eT/LMudXyGYeZmZXixGFmZqV4qGoY8dQiZjYQfMZhZmalOHGYmVkpThxmZlaKE4eZmZXii+NmNuw0ulHEz3f0D59xmJlZKU4cZmZWihOHmZmV4sRhZmalVJo4JM2StE5St6QFdY5L0iXp+GpJ01u1lXShpMdS/RsljanyM5iZ2dYqSxySRgGXAbOBacCpkqbVVJsNTE2vucDCAm1vB94fEfsCPwfOreozmJnZtqo845gBdEfE+oh4A1gGzKmpMwdYGpmVwBhJ45u1jYjbImJLar8SmFjhZzAzsxpVJo4JwDO5/Z5UVqROkbYAZwG31HtzSXMldUnq2rRpU8nQzcyskSofAFSdsihYp2VbSecBW4Dv13vziFgMLAbo7Oysfd8hzbPgmtlgqjJx9ACTcvsTgQ0F64xu1lbSGcAxwBERMaySgplZu6sycawCpkqaAjwLnAJ8sqbOcmC+pGXATODFiNgoaVOjtpJmAecAH4qIzRXGb2bDjKci6R+VJY6I2CJpPnArMApYEhFrJc1LxxcBK4CjgW5gM3Bms7ap60uBHYDbJQGsjIh5VX0OMzPbWqWTHEbECrLkkC9blNsO4OyibVP57/dzmGZmVoKfHDczs1KcOMzMrBQnDjMzK8WJw8zMSnHiMDOzUrx0rJmNeH6+oxwnjjbmqUXMrB15qMrMzErxGYeZWUkjfWjLicPMrAEPF9fnxNEG/MtpZkOJr3GYmVkpPuMwM+snI+XahxPHAPFwlJkNF04c/cwJwsyGO1/jMDOzUnzGYWZWseF27cNnHGZmVorPOPrI1zLMbKTyGYeZmZVSaeKQNEvSOkndkhbUOS5Jl6TjqyVNb9VW0m6Sbpf0ePr5nio/g5mZba2yoSpJo4DLgKOAHmCVpOUR8Uiu2mxganrNBBYCM1u0XQDcEREXpISyADinqs/hISkzq0p/XTRv9neqigvwVV7jmAF0R8R6AEnLgDlAPnHMAZZGRAArJY2RNB7oaNJ2DvDh1P67wJ1UmDjMzAZau/8Pa5WJYwLwTG6/h+ysolWdCS3a7hERGwEiYqOk3eu9uaS5wNy0+4qkdX35EHWMBV7op74GkuMeWEM1bhi6sTvuOvT1d9T89+oVVpk4VKcsCtYp0rapiFgMLC7TpghJXRHR2d/9Vs1xD6yhGjcM3dgd98Cp8uJ4DzAptz8R2FCwTrO2z6XhLNLP5/sxZjMza6HKxLEKmCppiqTRwCnA8po6y4HT091VBwEvpmGoZm2XA2ek7TOAmyv8DGZmVqOyoaqI2CJpPnArMApYEhFrJc1LxxcBK4CjgW5gM3Bms7ap6wuA6yR9CvgFcGJVn6GBfh/+GiCOe2AN1bhh6MbuuAeIshuazMzMivGT42ZmVooTh5mZleLEUYCkayU9mF5PSXowlXdIei13bNEgh7oNSV+T9GwuxqNzx85NU7qsk/TRwYyzlqQLJT2WpqK5UdKYVD4UvvOmU+20C0mTJP2bpEclrZX0V6m84e9Mu0j/Dh9O8XWlsraejkjS3rnv9EFJL0n63FD4vmv5GkdJkr5JdvfX+ZI6gB9GxPsHOayGJH0NeCUiLqopnwZcQ/aE/3uBHwF7RcSbAx5kHZI+Avw43SjxdYCIOKfdv/M0Xc7PyU2XA5xaM9VOW0i3s4+PiAck7QLcDxwHnESd35l2IukpoDMiXsiVfQP499x0RO+JiLacVSL9njxL9mDzmbT5913LZxwlSBLZP6prBjuWfjAHWBYRr0fEk2R3ts0Y5JjeEhG3RcSWtLuS7FmeoeCtqXYi4g2gd7qcthMRGyPigbT9MvAo2awNQ9UcsmmISD+PG7xQWjoCeCIinh7sQPrCiaOcQ4HnIuLxXNkUST+TdJekQwcrsBbmpyGfJbnT90bTvbSjs4Bbcvvt/J0Ppe/1LelMbn/g3lRU73emnQRwm6T70/RCUDMdEVB3OqI2cQpb/w9ou3/fW3HiSCT9SNKaOq/8/y2eytb/sTcCkyNif+ALwNWSdh3IuKFl7AuBPYH9Urzf7G1Wp6sBHbcs8p1LOg/YAnw/FbXFd97EoH+vZUn6LeCfgM9FxEs0/p1pJx+MiOlkM2yfLemPBzugopQ91HwscH0qGgrf91a8AmASEUc2Oy5pe+ATwAG5Nq8Dr6ft+yU9AewFdFUY6jZaxd5L0uXAD9NukSlhKlXgOz8DOAY4Is2g3DbfeROD/r2WIeldZEnj+xFxA0BEPJc7nv+daRsRsSH9fF7SjWRDhM9JGp8mP23n6YhmAw/0fs9D4fuu5TOO4o4EHouInt4CSePSRS4kvY9sXZH1gxRfXekfUK/jgTVpezlwiqQdJE0hi/2+gY6vEUmzyKbLPzYiNufK2/07LzLVTltI1+yuAB6NiG/lyhv9zrQFSe9OF/OR9G7gI2QxDpXpiLYauWj377sen3EUVzsmCfDHwPmStgBvAvMi4t8HPLLmviFpP7LhkqeAvwBI079cR7bGyRbg7Ha5oyq5FNgBuD37+8bKiJhHm3/nLabLaTcfBP4L8LDSLebAl4FT6/3OtJE9gBvT78X2wNUR8a+SVjG40xG1JGlnsjvu8t9p3X+j7cy345qZWSkeqjIzs1KcOMzMrBQnDjMzK8WJw8zMSnHiMDOzUpw4bEiRdKdqZvJNM4x+u0WbzorjuiZNGfH5mvLj0oSSfelzvzIzpUoaI+mzZetJeq+kH/RXfRv+nDhsqLmG7JmavHrP2AwYSb8LHBwR+0bExTWHjwP6lDjIpqAoM8X2GKBl4qitFxEbIuKEfqxvw5wThw01PwCOkbQDvDU533uBn0haKKlL2doSf1uvsaRXctsnSLoqbY+T9E+SVqXXB+u03VHSlcrWgfiZpMPSoduA3ZWtpXBorv7BZHMSXZiO7Zle/5om5/s/kvZJdU9M83Q9JOnu9NT5+cDJqe3JNbH8oaT70rHVkqYCFwB7prILJf2WpDskPZBi7p0DrLZeh6Q1JfrN1x8l6aLU/2pJf1n4v6QNXRHhl19D6gX8CzAnbS8ALkzbu6Wfo4A7gX3T/p1kazdAtu5Bbz8nAFel7auBQ9L2ZLJpOGrf978BV6btfcieTt4R6ADWNIj1KuCE3P4dwNS0PZNszRGAh4EJaXtM+vlnwKUN+v0H4LS0PRrYqTYOsqeqd03bY8mmzledem/tF+w3X/8zZHNdbZ//b+DX8H55yhEbinqHq25OP89K5Scpm2J7e2A82RDR6oJ9HglMS9NYAOwqaZfI1qnodQjZH1Yi4jFJT5NNsPhSkTdQNgvtwcD1uffZIf28B7gqTQNzQ4HufgqcJ2kicENEPJ7r8623BP5O2cyxvyGb3n2Pfug370hgUaS1U6KNpn+x6jhx2FB0E/AtSdOBnSJbwW4K8EXgwIj4f2kIasc6bfNz7OSPbwf8UUS81uR9m/4FLWA74D8iYr9tgoqYJ2km8DHgwTR3UUMRcbWke1P9WyX9OdtO9ngaMA44ICJ+rWzVvHrfSdl+80SbTxtv/c/XOGzIiYhXyIaflvD2RfFdgVeBFyXtQTZ1dT3PSfoDSduRzUTa6zZgfu9Ogz/cd5P9MUbSXmRDWutahPsysEuK+yXgSUknpj4k6QNpe8+IuDcivgq8QDY1+1ttaymbGXh9RFxCNivsvnXq/zbwfEoahwG/VxtTH/vNuw2Yp2zZASTt1uL7sGHAicOGqmuAD5AtzUpEPAT8DFhLllDuadBuAdl6Bz8mWzSn138FOtMF3keAeXXafhsYJelh4FrgzyJbH6SZZcCX0sX0PckSz6ckPZRi7b1gfWG6wLyGLEE9BPwb2fDZNhfHgZOBNcpmtd0HWBoRvwLuSRfZLyRb/KpTUld638fSd1Vbr2y/ef9Idq1ndfpMn2zxfdgw4NlxzcysFJ9xmJlZKU4cZmZWihOHmZmV4sRhZmalOHGYmVkpThxmZlaKE4eZmZXy/wHyovxJGFZxogAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"from itertools import combinations\n",
"import matplotlib.pyplot as plt\n",
"\n",
"def null_distribution(z, nx):\n",
" # z is the population of mice survival times in the study\n",
" # nx is the number of mice in the treatment group\n",
" z = set(z)\n",
" null_distribution = []\n",
" for x in combinations(z, nx):\n",
" y = z - set(x)\n",
" stat = statistic(list(x), list(y))\n",
" null_distribution.append(stat)\n",
" return null_distribution\n",
"\n",
"null_dist = null_distribution(z, len(x))\n",
"plt.hist(null_dist, density=True, bins=50)\n",
"plt.xlabel(\"Value of test statistic\")\n",
"plt.ylabel(\"Observed Frequency\")"
]
},
{
"cell_type": "markdown",
"id": "13914661-4e02-47db-82d0-a2302cc1e695",
"metadata": {},
"source": [
"We complete the hypothesis test by comparing the observed value of the test statistic to the rest of the null distribution."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "d0b5ba29-952e-4652-8f35-58928e76e6ad",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.1409965034965035\n"
]
}
],
"source": [
"pvalue = np.sum(null_dist >= statistic(x, y) ) / len(null_dist)\n",
"print(pvalue)"
]
},
{
"cell_type": "markdown",
"id": "90fd3def-2996-4366-a661-76c3216cd634",
"metadata": {},
"source": [
"Approximately 14% of the values of the null distribution are greater than the observed value of the statistic, so there is a 14% probability of observing such an extreme value of the statistic even if the the treatment had no effect at all. \n",
"\n",
"Given data and a statistic function, `stats.permutation_test` performs the same test automatically."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "8ff08473-77e4-4bc9-864e-2f06f19f2948",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"PermutationTestResult(statistic=30.63492063492064, pvalue=0.1409965034965035, null_distribution=array([ 30.63492063, 38. , 51.20634921, ..., -11.01587302,\n",
" -45.55555556, -34.88888889]))\n"
]
}
],
"source": [
"# `alternative` is 'greater' because we are interested in the percentage of values in the\n",
"# null distribution that are greater than the observed value of the test statistic.\n",
"# `n_resamples` is `np.inf` to ensure that all possible permutations are used\n",
"# Note that `(x, y)`, a tuple, is a single argument.\n",
"res = stats.permutation_test((x, y), statistic, alternative='greater', n_resamples = np.inf)\n",
"assert res.pvalue == pvalue\n",
"print(res)"
]
},
{
"cell_type": "markdown",
"id": "8b5bb238-7db9-40eb-8221-d703c3d26b70",
"metadata": {},
"source": [
"It returns the observed value of the test statistic, the null distribution, and the $p$-value. They are related as above:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "f303201f-cf56-46c1-8ee3-428afca93eb3",
"metadata": {},
"outputs": [],
"source": [
"assert np.sum(res.null_distribution >= res.statistic ) / len(res.null_distribution) == res.pvalue"
]
},
{
"cell_type": "markdown",
"id": "ae116861-ded9-4728-a1b1-4c9c34c50fdd",
"metadata": {},
"source": [
"Note that the exact $p$-value from the permutation test matches the $p$-value from the t-test quite closely. (As we shall see, Ronald Fisher introduced permutation tests primarily to support the use of the t-test in applications where the underlying normality assumptions were not strictly true [[5]](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2458144/https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2458144/).)"
]
},
{
"cell_type": "markdown",
"id": "ae96fea0-1255-48ac-9d7f-e56f1db10ae2",
"metadata": {},
"source": [
"### Randomized Tests\n",
"The number of possible permutations grows rather quickly as the number of observations increases. Specifically, if $n_x$ and $n_y$ are the number of observations in `x` and `y`, respectively, than the number of possible permutations is $\\frac{(n_x + n_y)!}{n_x! n_y!}$. "
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "d864ec56-9438-45f2-a41a-cb30e3c1ac35",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"11440\n"
]
}
],
"source": [
"from math import factorial as f\n",
"n_x = len(x)\n",
"n_y = len(y)\n",
"assert len(res.null_distribution) == f(n_x + n_y) / (f(n_x) * f(n_y))\n",
"print(len(res.null_distribution))"
]
},
{
"cell_type": "markdown",
"id": "432cff2d-f08e-4be3-99f2-c97a8f751e83",
"metadata": {},
"source": [
"When the number of possible permutations is too large, it is common to use a randomly-sampled subset of the possible permutations instead. As with `monte_carlo_test` the maximum number of resamples used by `permutation_test` is controlled using the `n_resamples` parameter."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "e9625f6e-d527-4e78-89fc-e24a603a2811",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.1494"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# use only 4999 randomly-sampled permutations\n",
"res = stats.permutation_test((x, y), statistic, alternative='greater', n_resamples = 4999)\n",
"res.pvalue"
]
},
{
"cell_type": "markdown",
"id": "e1937a92-5ad8-4164-a6c0-5ee01e6702b0",
"metadata": {},
"source": [
"If the number of distinct permutations of the data is less than or equal to `n_resamples`, `permutation_test` performs an exact test, computing the value of the test statistic for each distinct permutation exactly once. If the number of distinct permutations exceeds `n_resamples`, `permutation_test` computes the value of the statistic for `n_resamples` random permutations, and the $p$-value is computed as:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "8e7e5525-5cbf-43df-86b2-7af081ad96e4",
"metadata": {},
"outputs": [],
"source": [
"pvalue = (np.sum(res.null_distribution >= res.statistic ) + 1) / (len(res.null_distribution) + 1)\n",
"assert pvalue == res.pvalue"
]
},
{
"cell_type": "markdown",
"id": "459d7ca3-f719-4b8e-a855-8e47846c3ec0",
"metadata": {},
"source": [
"Note that `1` is added to both the numerator and denominator when performing the randomized test [[6]](https://www.degruyter.com/document/doi/10.2202/1544-6115.1585/html). This can be thought of as including the observed value of the test statistic in the null distribution, and it ensures that the $p$-value of a randomized test is never zero."
]
},
{
"cell_type": "markdown",
"id": "cb30657f-b7db-4762-95c2-1eebb2a05a64",
"metadata": {},
"source": [
"A wide variety of common hypothesis tests can be performed as permutation tests. We continue with several other examples to explore the flexibility of `permutation_test`, beginning with [Independent-Sample Tests](https://nbviewer.org/github/scipy/scipy-cookbook/blob/main/ipython/ResamplingAndMonteCarloMethods/resampling_tutorial_2a.ipynb)."
]
}
],
"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.10.5"
}
},
"nbformat": 4,
"nbformat_minor": 5
}