Analysis of Simpson's paradox example dataset.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4dd74f9d",
   "metadata": {},
   "source": [
    "# Analysis of Simpson's paradox example dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dba34439",
   "metadata": {},
   "source": [
    "[Simpson's Paradox](https://plato.stanford.edu/entries/paradox-simpson/) is a statistical phenomenon where an association between two variables in a population emerges, disappears or reverses when the population is divided into subpopulations. It can have serious consequences for decision making because people easily confuse association with causal relationship and infer a wrong conclusion. Classical example of Simpson's paradox is the [University of California as Berkeley 1973 admission dataset](https://discovery.cs.illinois.edu/dataset/berkeley/). We shall briefly analyze this dataset here in order to resolve the paradox."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7d02160",
   "metadata": {},
   "source": [
    "## Was there a sex bias in Berkeley graduate admissions?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30431018",
   "metadata": {},
   "source": [
    "In the Fall of 1973, the University of California at Berkeley released data about their graduate admissions. At the first glance, the data showed a significant difference in the acceptance rate for male and female applicants. They had admitted 44.3% of men and only 34.6% of women.\n",
    "\n",
    "Is this difference really significant? Is there a bias in the application procedure? Is there something that applicants should complain about? Did we have sex-based discrimination at Berkley? We shall try to answer these questions in this article using some data analysis and visualization."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "502cc3b1",
   "metadata": {},
   "source": [
    "Firstly, let's see the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "fb247a42",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Year</th>\n",
       "      <th>Major</th>\n",
       "      <th>Gender</th>\n",
       "      <th>Admission</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1973</td>\n",
       "      <td>C</td>\n",
       "      <td>F</td>\n",
       "      <td>Rejected</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1973</td>\n",
       "      <td>B</td>\n",
       "      <td>M</td>\n",
       "      <td>Accepted</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1973</td>\n",
       "      <td>Other</td>\n",
       "      <td>F</td>\n",
       "      <td>Accepted</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1973</td>\n",
       "      <td>Other</td>\n",
       "      <td>M</td>\n",
       "      <td>Accepted</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1973</td>\n",
       "      <td>Other</td>\n",
       "      <td>M</td>\n",
       "      <td>Rejected</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12758</th>\n",
       "      <td>1973</td>\n",
       "      <td>Other</td>\n",
       "      <td>M</td>\n",
       "      <td>Accepted</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12759</th>\n",
       "      <td>1973</td>\n",
       "      <td>D</td>\n",
       "      <td>M</td>\n",
       "      <td>Accepted</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12760</th>\n",
       "      <td>1973</td>\n",
       "      <td>Other</td>\n",
       "      <td>F</td>\n",
       "      <td>Rejected</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12761</th>\n",
       "      <td>1973</td>\n",
       "      <td>Other</td>\n",
       "      <td>M</td>\n",
       "      <td>Rejected</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12762</th>\n",
       "      <td>1973</td>\n",
       "      <td>Other</td>\n",
       "      <td>M</td>\n",
       "      <td>Accepted</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>12763 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       Year  Major Gender Admission\n",
       "0      1973      C      F  Rejected\n",
       "1      1973      B      M  Accepted\n",
       "2      1973  Other      F  Accepted\n",
       "3      1973  Other      M  Accepted\n",
       "4      1973  Other      M  Rejected\n",
       "...     ...    ...    ...       ...\n",
       "12758  1973  Other      M  Accepted\n",
       "12759  1973      D      M  Accepted\n",
       "12760  1973  Other      F  Rejected\n",
       "12761  1973  Other      M  Rejected\n",
       "12762  1973  Other      M  Accepted\n",
       "\n",
       "[12763 rows x 4 columns]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "df=pd.read_csv(\"https://waf.cs.illinois.edu/discovery/berkeley.csv\")\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb64c146",
   "metadata": {},
   "source": [
    "Here we have application outcomes for 12,763 applicants. They were either accepted or rejected. The __Major__ column indicates which department the candidates applied for. Actually, we see just the labels, but this should be enough for the purpose of statistical analysis. We see, also, applicants' gender as well as the year of admission. The year data is always the same: 1973. We could ignore it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ac331ea5",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df.drop('Year',axis=1, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39e1116a",
   "metadata": {},
   "source": [
    "Let's see how many male and female applicants were accepted and rejected. For the purpose, we will use pandas function for cross-tabulation: `crosstab()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "230c9333",
   "metadata": {},
   "outputs": [],
   "source": [
    "ct=pd.crosstab(df.Gender, df.Admission)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "eb53f6cc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>Admission</th>\n",
       "      <th>Accepted</th>\n",
       "      <th>Rejected</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gender</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>F</th>\n",
       "      <td>1494</td>\n",
       "      <td>2827</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M</th>\n",
       "      <td>3738</td>\n",
       "      <td>4704</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Admission  Accepted  Rejected\n",
       "Gender                       \n",
       "F              1494      2827\n",
       "M              3738      4704"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# cross-tabulation\n",
    "ct"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dae6c5f2",
   "metadata": {},
   "source": [
    "We can see that much fewer women were accepted but also that much fewer women applied. It is hard to say whether the admission difference between male and female applicants is significant. We can test this doubt in two ways. One is to apply a statistical test. That is how most scientists would do it. The other approach is how we could do it without statistical tests - just with a simulation of possible outcomes. We will start with the latter approach."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c42bd5e",
   "metadata": {},
   "source": [
    "### Simulation of an unbiased admission"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6866b79c",
   "metadata": {},
   "source": [
    "We shall assume that the application procedure is unbiased and simulate possible admission results. For every applicant in the array, we should simulate admission in the same way regardless of the candidate's gender. In order to create such a simulation, we need an array of random numbers between 0 and 1 as long as the array of candidates is. If a number in the array is greater than the threshold, we would say the candidate is accepted. Otherwise, we would say he/she is rejected. For the threshold, we will use the ratio of the number of candidates who were accepted and the number of applicants. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "068b1333",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "63192667",
   "metadata": {},
   "outputs": [],
   "source": [
    "# total number of candidates and number of those who were accepted\n",
    "n=len(df)\n",
    "a=len(df[df.Admission=='Accepted'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66d1752c",
   "metadata": {},
   "source": [
    "The array **sim_adm** contains values 0 and 1 for rejected and accepted candidates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "176c701a",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim_adm=np.array(np.random.rand(n)+a/n,dtype=\"int\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fb49fb4",
   "metadata": {},
   "source": [
    "Let's see what the crosstabulation would look like in this particular simulated case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7d85e075",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>col_0</th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gender</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>F</th>\n",
       "      <td>2532</td>\n",
       "      <td>1789</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M</th>\n",
       "      <td>4940</td>\n",
       "      <td>3502</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "col_0      0     1\n",
       "Gender            \n",
       "F       2532  1789\n",
       "M       4940  3502"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.crosstab(df.Gender,sim_adm)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48134241",
   "metadata": {},
   "source": [
    "Again, we have much fewer female applicants and much fewer those who were admitted. In this simulated case 1806 women were accepted while in the observed dataset there are 1494 of them. Is this difference significant? We can't say that yet. One simulated outcome is not enough to infer anything. We will repeat the simulation 1000 times, record the number of accepted female applicants for all of them and display this data graphically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "8c65242e",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim_fa=np.zeros(1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "4dc7d64b",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(1000):\n",
    "    sim_adm=np.array(np.random.rand(n)+a/n,dtype=\"int\")\n",
    "    fem_acc=pd.crosstab(df.Gender,sim_adm).iloc[0,1]\n",
    "    sim_fa[i]=fem_acc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ee4e2b39",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "337e9aa1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAPlklEQVR4nO3dXYxcZ33H8e+vMeSiQJU0Tuo6VjdCBpFcYGDrUgHlrSVpItVJqyDnglpqJFOUtATRqg6tSqTKkqG8iKqFyjQRpqUkrgiNq9BCsKAREo2zidK8mYAhLlls2UtBIr0JjfPvxRwrU2dfZnd2dmcffz/SaM485zk7P49Xvz37+Mw4VYUkqS0/s9oBJEnLz3KXpAZZ7pLUIMtdkhpkuUtSg9atdgCACy64oCYmJlY7hiStKQ888MAPq2r9bPvGotwnJiaYmppa7RiStKYk+a+59rksI0kNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDRqLd6hKozCx6+4F5xzdc9UKJJFWnmfuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUoAXLPcmmJF9LcjjJY0ne243fkuQHSR7qblf2HXNzkiNJnkhy+Sj/AJKkFxrk4weeBd5fVQ8meSnwQJJ7un0fr6qP9E9OcimwHbgM+EXgq0leUVWnljO4JGluC565V9Xxqnqw234aOAxsnOeQbcDtVfVMVT0JHAG2LkdYSdJgFrXmnmQCeA1wXzd0Y5KHk9yW5LxubCPwVN9h08zywyDJziRTSaZmZmYWn1ySNKeByz3JS4AvADdV1U+ATwEvB7YAx4GPnp46y+H1goGqvVU1WVWT69evX2xuSdI8Bir3JC+iV+yfq6o7AarqRFWdqqrngE/z/NLLNLCp7/CLgWPLF1mStJBBrpYJcCtwuKo+1je+oW/aNcCj3fYBYHuSc5NcAmwGDi1fZEnSQga5WuYNwLuAR5I81I19ALguyRZ6Sy5HgXcDVNVjSfYDj9O70uYGr5SRpJW1YLlX1TeYfR39S/McsxvYPUQuSdIQfIeqJDXIcpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KDLHdJatC6hSYk2QR8FvgF4Dlgb1V9Isn5wB3ABHAUeGdV/bg75mbgeuAU8IdV9eWRpJeGNLHr7gXnHN1z1QokkZbXIGfuzwLvr6pXAa8HbkhyKbALOFhVm4GD3WO6fduBy4ArgE8mOWcU4SVJs1uw3KvqeFU92G0/DRwGNgLbgH3dtH3A1d32NuD2qnqmqp4EjgBblzm3JGkei1pzTzIBvAa4D7ioqo5D7wcAcGE3bSPwVN9h093YmV9rZ5KpJFMzMzNLiC5JmsvA5Z7kJcAXgJuq6ifzTZ1lrF4wULW3qiaranL9+vWDxpAkDWCgck/yInrF/rmqurMbPpFkQ7d/A3CyG58GNvUdfjFwbHniSpIGsWC5JwlwK3C4qj7Wt+sAsKPb3gHc1Te+Pcm5SS4BNgOHli+yJGkhC14KCbwBeBfwSJKHurEPAHuA/UmuB74PXAtQVY8l2Q88Tu9Kmxuq6tRyB5ckzW3Bcq+qbzD7OjrA2+c4Zjewe4hckqQh+A5VSWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1KBBrnOXxs4gH9Urnc08c5ekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ1asNyT3JbkZJJH+8ZuSfKDJA91tyv79t2c5EiSJ5JcPqrgkqS5DXLm/hngilnGP15VW7rblwCSXApsBy7rjvlkknOWK6wkaTALlntV3Qv8aMCvtw24vaqeqaongSPA1iHySZKWYJg19xuTPNwt25zXjW0EnuqbM92NvUCSnUmmkkzNzMwMEUOSdKallvungJcDW4DjwEe78cwyt2b7AlW1t6omq2py/fr1S4whSZrNksq9qk5U1amqeg74NM8vvUwDm/qmXgwcGy6iJGmxllTuSTb0PbwGOH0lzQFge5Jzk1wCbAYODRdRkrRY6xaakOTzwFuAC5JMAx8E3pJkC70ll6PAuwGq6rEk+4HHgWeBG6rq1EiSS5LmtGC5V9V1swzfOs/83cDuYUJJkobjO1QlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KDLHdJatC61Q4gjbuJXXcvOOfonqtWIIk0OM/cJalBlrskNchlGY2dQZZBJM3PM3dJapDlLkkNstwlqUELlnuS25KcTPJo39j5Se5J8p3u/ry+fTcnOZLkiSSXjyq4JGlug5y5fwa44oyxXcDBqtoMHOwek+RSYDtwWXfMJ5Ocs2xpJUkDWbDcq+pe4EdnDG8D9nXb+4Cr+8Zvr6pnqupJ4AiwdXmiSpIGtdQ194uq6jhAd39hN74ReKpv3nQ39gJJdiaZSjI1MzOzxBiSpNks9z+oZpaxmm1iVe2tqsmqmly/fv0yx5Cks9tSy/1Ekg0A3f3Jbnwa2NQ372Lg2NLjSZKWYqnlfgDY0W3vAO7qG9+e5NwklwCbgUPDRZQkLdaCHz+Q5PPAW4ALkkwDHwT2APuTXA98H7gWoKoeS7IfeBx4Frihqk6NKLskaQ4LlntVXTfHrrfPMX83sHuYUJKk4fgOVUlqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktSgdcMcnOQo8DRwCni2qiaTnA/cAUwAR4F3VtWPh4spSVqMocq989aq+mHf413Awarak2RX9/hPluF5pLE1sevuBecc3XPVCiSRekaxLLMN2Ndt7wOuHsFzSJLmMWy5F/CVJA8k2dmNXVRVxwG6+wtnOzDJziRTSaZmZmaGjCFJ6jfssswbqupYkguBe5J8a9ADq2ovsBdgcnKyhswhSeoz1Jl7VR3r7k8CXwS2AieSbADo7k8OG1KStDhLLvckP5vkpae3gXcAjwIHgB3dtB3AXcOGlCQtzjDLMhcBX0xy+uv8Y1X9W5L7gf1Jrge+D1w7fExJ0mIsudyr6nvAq2cZ/2/g7cOEkiQNx3eoSlKDLHdJapDlLkkNWo6PH5AGNsjb9CUNzzN3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KD/M86pBUyyH9UcnTPVSuQRGcDz9wlqUGWuyQ1yHKXpAa55q5l439+LY0Pz9wlqUGWuyQ1yHKXpAZZ7pLUIMtdkho0sqtlklwBfAI4B/i7qtozqufS6HkljLS2jKTck5wD/A3wG8A0cH+SA1X1+CieT2qFH1Gg5TKqZZmtwJGq+l5V/RS4Hdg2oueSJJ1hVMsyG4Gn+h5PA7/SPyHJTmBn9/B/kjwxoixLdQHww9UOsQhrKe9aygpjljcfmnf3WGUdwFrKO45Zf2muHaMq98wyVv/vQdVeYO+Inn9oSaaqanK1cwxqLeVdS1lhbeVdS1lhbeVdS1lhdMsy08CmvscXA8dG9FySpDOMqtzvBzYnuSTJi4HtwIERPZck6QwjWZapqmeT3Ah8md6lkLdV1WOjeK4RGtslozmspbxrKSusrbxrKSusrbxrKSupqoVnSZLWFN+hKkkNstwlqUFnVbknuS3JySSP9o3dkuQHSR7qblf27bs5yZEkTyS5vG/8dUke6fb9VZLZLv1c9qzd+B90eR5L8uFxyDpX3iR39L2uR5M8NA5558i6Jcl/dFmnkmwdh6zz5H11km92z/8vSV42DnmTbErytSSHu+/R93bj5ye5J8l3uvvzVjvvPFmv7R4/l2TyjGNW9XthUarqrLkBvwa8Fni0b+wW4I9mmXsp8J/AucAlwHeBc7p9h4BfpXc9/78Cv7lCWd8KfBU4t3t84ThknSvvGfs/Cvz5OOSd47X9yunnAq4Evj4OWefJez/w5m7794C/GIe8wAbgtd32S4Fvd5k+DOzqxncBH1rtvPNkfRXwSuDrwGTf/FX/XljM7aw6c6+qe4EfDTh9G3B7VT1TVU8CR4CtSTYAL6uqb1bvb/WzwNUrlPU9wJ6qeqabc3Icss6TF4DuLOadwOfHIe8cWQs4ffb7czz/voxxfW1fCdzbbd8D/M445K2q41X1YLf9NHCY3jvWtwH7umn7+p571fLOlbWqDlfVbO+YX/XvhcU4q8p9Hjcmebj79ff0r4uzfYTCxu42Pcv4SngF8KYk9yX59yS/3I2PY9Z+bwJOVNV3usfjmPcm4C+TPAV8BLi5Gx/HrACPAr/VbV/L828aHJu8SSaA1wD3ARdV1XHolSpw4TjlPSPrXMYi66Asd/gU8HJgC3Cc3vIBzP0RCgt+tMIIrQPOA14P/DGwvzsrHses/a7j+bN2GM+87wHeV1WbgPcBt3bj45gVeksxNyR5gN6Swk+78bHIm+QlwBeAm6rqJ/NNnWVsRfOupayLcdaXe1WdqKpTVfUc8Gl6n2gJc3+EwnS3feb4SpgG7qyeQ8Bz9D7MaByzApBkHfDbwB19w+OYdwdwZ7f9T4z39wFV9a2qekdVvY7eD87vdrtWPW+SF9Ery89V1enX9ES3fEF3f3pJcVXzzpF1Lqv+2i7GWV/up7/hOtfQ+3UXeh+XsD3JuUkuATYDh7pfKZ9O8vrurPl3gbtWKO4/A2/rcr8CeDG9T6kbx6yn/Trwrarq/7V1HPMeA97cbb8NOL2ENI5ZSXJhd/8zwJ8BfzsOebuvfStwuKo+1rfrAL0foHT3d/WNr0reebLOZSy/F+a02v+iu5I3emc4x4H/pffT9nrg74FHgIfp/eVt6Jv/p/TOiJ6g71+/gUl6PwS+C/w13Tt9VyDri4F/6J77QeBt45B1rrzd+GeA359l/ri9tm8EHqB3NcR9wOvGIes8ed9L7+qObwN7+p97lV/bN9JbkngYeKi7XQn8PHCQ3g/Ng8D5q513nqzXdK/zM8AJ4MurnXUpNz9+QJIadNYvy0hSiyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1KD/A/1748qw4+XaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(sim_fa,bins=np.linspace(1450,2150,35));"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e54af21f",
   "metadata": {},
   "source": [
    "What we see here is the distribution of the number of accepted female candidates in 1000 simulated unbiased admissions. This is what we should expect for the number of accepted women if the number of female applicants stays the same. The distribution of simulated outcomes is bell-shaped where the great majority of outcomes lie between 1700 and 1850. The real-data outcome (1494) is far smaller than the expected values. In other words, the probability that Berkeley accepts only 1494 women is very low, much lower than 1 in 1000."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c80ccd0c-16eb-4002-9091-81fb668898a2",
   "metadata": {},
   "source": [
    "It is clear that the acceptance rate for male and female candidates for the University of California at Berkeley was not the same. For women, it was significantly more difficult to enroll in 1973. Is this evidence hard enough to accuse Berkeley of discrimination? No! There are at least three reasons why we shouldn't run into the conclusion that the bias was intentional Firstly, events with the probability less than 0.1% happen all the time. It could be just a coincidence. Secondly, statistically significant differences without the cause can't be taken as serious evidence. We need some theory behind this data. Thirdly, there might be some confounding variable, which has nothing to do with Berkeley, influencing the acceptance rate."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c65b50dd",
   "metadata": {},
   "source": [
    "### Statistical test $\\chi^2$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6355db11-6d97-4696-9a26-174917c15bab",
   "metadata": {},
   "source": [
    "There is a commonly used statistical test for contingency tables called $\\chi^2$ or chi-square or just chi2 (if you want to google it). Pearson's $\\chi^2$ test determines whether there is a statistically significant difference between the expected frequencies and the observed frequencies in one or more categories of a contingency table. A function for this test can be found in scipy.stats library.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "8712805c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import chi2_contingency "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "de0b023b",
   "metadata": {},
   "outputs": [],
   "source": [
    "stat, p, dof, arr=chi2_contingency(ct)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58b8e0fe",
   "metadata": {},
   "source": [
    "Here we have three outputs: statistic $\\chi^2$ (__stat__), _p_-value (__p__) and the number of the degrees of freedom (__dof__). Right now, we are interested in _p_-value only. It stands for estimated probability that observed results can be this far from results expected under the assumption that acceptance rate does not depend on the gender. Using simulation, we have already estimated that this probability must be much lower than 0.001. Let's see what statistical test says."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "6f285cf7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.385628329926473e-26"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "882dcfe8",
   "metadata": {},
   "source": [
    "If this exponential form of the number is not convincing enough, we can print it with all the leading zeroes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "3a2e9268",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p=0.00000000000000000000000006\n"
     ]
    }
   ],
   "source": [
    "print(\"p=%0.26f\" % p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c762b125",
   "metadata": {},
   "source": [
    "Definitely, the probability of the observed outcome is extremely low. If we had repeated admission at Berkeley every day since the Big Bang, it still would be extremely unlikely that we ever obtain such an outcome. Can we now accuse Berkeley of discrimination? Such results can't be a coincidence. The answer is the same as before. No, we need at least \"a motive for murder\" like in criminalistic series. Instead, we can look for the confounding variables."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c36edab4",
   "metadata": {},
   "source": [
    "### Looking for the confounding variable"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15b13ce2",
   "metadata": {},
   "source": [
    "The only recorded variable that could be confounding is __Major__. Results that look biased for the whole university might be different if we stratify the sample and look at subgroups, i.e. departments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "c201c754",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Other    7924\n",
       "A        1246\n",
       "C         918\n",
       "D         792\n",
       "F         714\n",
       "B         585\n",
       "E         584\n",
       "Name: Major, dtype: int64"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df['Major'].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c25e822",
   "metadata": {},
   "source": [
    "Unfortunately, we have labels for the six biggest departments only (A-E). All the others are put in a vast category \"Other\". If the major is really the confounding variable, mixing outcomes into the \"Other\" category would just blur the results. We should exclude this category. Anyway, Simpson's paradox appears to be even more obvious if we take the candidates who applied for the labeled departments only. Therefore, we will exclude those who applied for \"Other\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "74def5f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "df2=df[df.Major!='Other']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "b40707cf",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Major</th>\n",
       "      <th>Gender</th>\n",
       "      <th>Admission</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>C</td>\n",
       "      <td>F</td>\n",
       "      <td>Rejected</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>B</td>\n",
       "      <td>M</td>\n",
       "      <td>Accepted</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>F</td>\n",
       "      <td>F</td>\n",
       "      <td>Accepted</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>A</td>\n",
       "      <td>M</td>\n",
       "      <td>Accepted</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>B</td>\n",
       "      <td>M</td>\n",
       "      <td>Accepted</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12747</th>\n",
       "      <td>A</td>\n",
       "      <td>M</td>\n",
       "      <td>Accepted</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12748</th>\n",
       "      <td>F</td>\n",
       "      <td>M</td>\n",
       "      <td>Rejected</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12752</th>\n",
       "      <td>A</td>\n",
       "      <td>M</td>\n",
       "      <td>Accepted</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12754</th>\n",
       "      <td>A</td>\n",
       "      <td>M</td>\n",
       "      <td>Accepted</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12759</th>\n",
       "      <td>D</td>\n",
       "      <td>M</td>\n",
       "      <td>Accepted</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>4839 rows × 3 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      Major Gender Admission\n",
       "0         C      F  Rejected\n",
       "1         B      M  Accepted\n",
       "6         F      F  Accepted\n",
       "9         A      M  Accepted\n",
       "11        B      M  Accepted\n",
       "...     ...    ...       ...\n",
       "12747     A      M  Accepted\n",
       "12748     F      M  Rejected\n",
       "12752     A      M  Accepted\n",
       "12754     A      M  Accepted\n",
       "12759     D      M  Accepted\n",
       "\n",
       "[4839 rows x 3 columns]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09b91d35",
   "metadata": {},
   "source": [
    "> Suggestion: Repeat simulation and statistical test for this truncated DataFrame. Do we still have an extremely low _p_-value again?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "1cb9e273",
   "metadata": {},
   "outputs": [],
   "source": [
    "dfA=df[df.Major=='A']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "e56ad906",
   "metadata": {},
   "outputs": [],
   "source": [
    "ctA=pd.crosstab(dfA.Gender, dfA.Admission)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "4dbfa68f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>Admission</th>\n",
       "      <th>Accepted</th>\n",
       "      <th>Rejected</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gender</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>F</th>\n",
       "      <td>89</td>\n",
       "      <td>19</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M</th>\n",
       "      <td>825</td>\n",
       "      <td>313</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Admission  Accepted  Rejected\n",
       "Gender                       \n",
       "F                89        19\n",
       "M               825       313"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ctA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "3adc9ea6",
   "metadata": {},
   "outputs": [],
   "source": [
    "stat, p, dof, arr=chi2_contingency(ctA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "9b4d2a47",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.034619660876707294"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbc94a02",
   "metadata": {},
   "source": [
    "The difference is still statistically significant, barely, but _p_<0.05. Although the significance is not so convincing, we see the same bias. But wait! There is a difference in acceptance rate but female candidates for the major A have a GREATER acceptance rate than the male! For women, acceptance rate is 87% while for men it is 72%. It seems that we have found a hole in that biased admission logic. Major could be the confounding variable."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d4c8b7f",
   "metadata": {},
   "source": [
    "You can do the same statistical analysis for all majors, but it could be better to visualize these contingency tables. There are a lot of ways to visualize a Simpson's paradox example dataset, which don't contribute much to the clarity. We want something better, a visualization that reveals what really matters."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01892ce1",
   "metadata": {},
   "source": [
    "### Visualization of the Berkeley dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f94492b8",
   "metadata": {},
   "source": [
    "This visualization has several steps since we have to create it from scratch. Already available charts won't be good enough. Be patient. We will do it step by step.\n",
    "\n",
    "Firstly, we would like to compare the acceptance rate for both sexes and all available majors regardless of the number of applicants. That is why on the y-axis we should put the acceptance rate given as a percentage. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "8679508a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# maximal height (percentage)\n",
    "h=100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "19b2115d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of F and M applicants per major\n",
    "noa=pd.crosstab(df2.Gender,df2.Major)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "cf86aba7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# total number of applicants per major\n",
    "toa=noa.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "c81cc3ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of majors\n",
    "nom=len(df2.Major.unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "b15e641e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# list of outcomes\n",
    "outcomes=df2.Admission.unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "37ae07ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of accepted candidates per major\n",
    "nos=pd.crosstab(df2['Gender'][df2.Admission==outcomes[1]],df2['Major'][df2.Admission==outcomes[1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9c29d81",
   "metadata": {},
   "source": [
    "We will display the percentage of accepted candidates as the light part of the bar. The rest up to 100% is dark. They are rejected. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "4fffd820",
   "metadata": {},
   "outputs": [],
   "source": [
    "# acceptance rate DataFrame\n",
    "ar=nos/noa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "353b7a3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# bar height (percentage of acceptance)\n",
    "ha=ar.transpose().to_numpy().flatten()*100"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c337b69d",
   "metadata": {},
   "source": [
    "The number of candidates will be displayed through the bars' width. The wider bar, the more candidates applied. Of course, the width should be scaled. Width 1.00 will be assigned to the major with the most applicants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "cd2ef097",
   "metadata": {},
   "outputs": [],
   "source": [
    "# relative weight (number of applicants)\n",
    "w=noa/max(toa)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "026bed7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set x-axis\n",
    "x=np.linspace(1,nom,nom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "2ba5793c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# distance between middles of adjacent bars\n",
    "interspace=1.05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "a68cd3a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# shifted bars' positions for male and female \n",
    "shift=pd.DataFrame([x-(toa/max(toa)-w.iloc[0,:])/2/interspace,x-(toa/max(toa)-2*w.iloc[0,:]-w.iloc[1,:])/2/interspace])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "559080c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# rearanged x-axis\n",
    "x2=shift.transpose().to_numpy().flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "0474939b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# bars' width\n",
    "wl=w.transpose().to_numpy().flatten()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4922192",
   "metadata": {},
   "source": [
    "Finally, we plot six groups of bars for each major. The left bar is for female candidates, and the right is for males."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "a1d918d8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6UAAAD4CAYAAAAU2UDyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAPaklEQVR4nO3da6zkd13H8c/XPdwrsbVn69oWF5JVrE0EsiIX460S8BJbEpq0BrMxjfsEFNQtWXhSnxCbSAgmXsgG0EWR2gBJG2PQZpF4R7elRsqKrVTLyto9gDeMAq1fH+wsOeye7Z6eObO/mdnXK2lm5j8zZ77Z+WVO3+f/n5nq7gAAAMAIXzd6AAAAAC5eohQAAIBhRCkAAADDiFIAAACGEaUAAAAMszJ6gCS5/PLLe/fu3aPHAAAAYAbuvffez3X36kbXzUWU7t69O0ePHh09BgAAADNQVf98ruscvgsAAMAwohQAAIBhRCkAAADDiFIAAACGEaUAAAAMI0oBAAAY5rxRWlXvqaqTVfWJddsuq6p7qurByeml6657c1U9VFWfqqpXzmpwAAAAFt9m9pT+VpJXnbHtYJIj3b0nyZHJ5VTVNUluSvIdk/v8elXt2LZpAQAAWCrnjdLu/pMkXzhj8/VJDk/OH05yw7rtd3T3l7r74SQPJXnx9owKAADAslnZ4v2u6O4TSdLdJ6pq52T7lUn+at3tjk+2naWq9ifZnyTPec5ztjjGhfO2qk3f9tb88gwnYbTuA6NH2JT1a/ZiWpOL8vxM68m8Ji2beVnPi77W5n0NzcvzvFmLvh62w7yvqdO2e2157jc2b+thnl9TrKHt/6CjjVZfb3TD7j7U3Xu7e+/q6uo2jwEAAMAi2GqUPlpVu5Jkcnpysv14kqvX3e6qJJ/d+ngAAAAss61G6d1J9k3O70ty17rtN1XV06rquUn2JPnr6UYEAABgWZ33PaVV9f4k35/k8qo6nuS2JLcnubOqbknySJIbk6S7H6iqO5N8MsljSV7X3Y/PaHYAAAAW3HmjtLtvPsdV153j9m9N8tZphgIAAODisN0fdAQAAACbJkoBAAAYRpQCAAAwjCgFAABgGFEKAADAMKIUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGEaUAgAAMIwoBQAAYBhRCgAAwDCiFAAAgGFEKQAAAMOIUgAAAIYRpQAAAAwjSgEAABhGlAIAADCMKAUAAGAYUQoAAMAwohQAAIBhRCkAAADDiFIAAACGEaUAAAAMI0oBAAAYRpQCAAAwjCgFAABgGFEKAADAMKIUAACAYVZGD7AoDjyJ2946sykAAACWiz2lAAAADDNVlFbVz1XVA1X1iap6f1U9vaouq6p7qurByeml2zUsAAAAy2XLUVpVVyb52SR7u/vaJDuS3JTkYJIj3b0nyZHJZQAAADjLtIfvriR5RlWtJHlmks8muT7J4cn1h5PcMOVjAAAAsKS2HKXd/S9J3pbkkSQnkvxHd/9Rkiu6+8TkNieS7Nzo/lW1v6qOVtXRtbW1rY4BAADAApvm8N1Lc2qv6HOTfHOSZ1XVazd7/+4+1N17u3vv6urqVscAAABggU1z+O4PJXm4u9e6+ytJPpTkZUkerapdSTI5PTn9mAAAACyjaaL0kSQvqapnVlUluS7JsSR3J9k3uc2+JHdNNyIAAADLamWrd+zuj1XVB5Lcl+SxJB9PcijJJUnurKpbcipcb9yOQQEAAFg+W47SJOnu25LcdsbmL+XUXlMAAAB4QtN+JQwAAABsmSgFAABgGFEKAADAMKIUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGEaUAgAAMIwoBQAAYJiV0QMspwOjBwAAAFgI9pQCAAAwjCgFAABgGFEKAADAMKIUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGGZl9ADAbBxYd/7WYVMAAMATs6cUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGEaUAgAAMIwoBQAAYBhRCgAAwDCiFAAAgGFWRg8AXAgHRg8AAAAbsqcUAACAYUQpAAAAw0wVpVX1DVX1gar6+6o6VlUvrarLquqeqnpwcnrpdg0LAADAcpl2T+mvJPlwdz8/yXcmOZbkYJIj3b0nyZHJZQAAADjLlqO0qp6d5HuTvDtJuvvL3f3vSa5Pcnhys8NJbphuRAAAAJbVNHtKn5dkLclvVtXHq+pdVfWsJFd094kkmZzu3IY5AQAAWELTROlKkhcl+Y3ufmGS/86TOFS3qvZX1dGqOrq2tjbFGAAAACyqaaL0eJLj3f2xyeUP5FSkPlpVu5Jkcnpyozt396Hu3tvde1dXV6cYAwAAgEW15Sjt7n9N8pmq+rbJpuuSfDLJ3Un2TbbtS3LXVBMCAACwtFamvP/PJHlfVT01yaeT/FROhe6dVXVLkkeS3DjlYwAAALCkporS7r4/yd4Nrrpump8LAADAxWHa7ykFAACALROlAAAADCNKAQAAGEaUAgAAMIwoBQAAYBhRCgAAwDCiFAAAgGFEKQAAAMOIUgAAAIYRpQAAAAwjSgEAABhGlAIAADCMKAUAAGAYUQoAAMAwohQAAIBhRCkAAADDiFIAAACGEaUAAAAMI0oBAAAYRpQCAAAwjCgFAABgGFEKAADAMKIUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGEaUAgAAMIwoBQAAYBhRCgAAwDCiFAAAgGFEKQAAAMOIUgAAAIaZOkqrakdVfbyqfn9y+bKquqeqHpycXjr9mAAAACyj7dhT+oYkx9ZdPpjkSHfvSXJkchkAAADOMlWUVtVVSX40ybvWbb4+yeHJ+cNJbpjmMQAAAFhe0+4pfUeSNyX5v3XbrujuE0kyOd250R2ran9VHa2qo2tra1OOAQAAwCLacpRW1Y8lOdnd927l/t19qLv3dvfe1dXVrY4BAADAAluZ4r4vT/LjVfUjSZ6e5NlV9TtJHq2qXd19oqp2JTm5HYMCAACwfLa8p7S739zdV3X37iQ3JflId782yd1J9k1uti/JXVNPCQAAwFKaxfeU3p7kFVX1YJJXTC4DAADAWaY5fPeruvujST46Of/5JNdtx88FAABguW1LlALAhXNg9AAAwDaaxeG7AAAAsCmiFAAAgGFEKQAAAMN4TykAXITm/Z25t44eAIALxp5SAAAAhhGlAAAADOPwXYAFNO+HXs6SwzoBYLmIUgBgDl3Mf3oBuLg4fBcAAIBhRCkAAADDiFIAAACGEaUAAAAMI0oBAAAYRpQCAAAwjCgFAABgGFEKAADAMKIUAACAYUQpAAAAw6yMHgAAAC6cA6MHAM5gTykAAADDiFIAAACGEaUAAAAMI0oBAAAYRpQCAAAwjCgFAABgGFEKAADAML6nFACAqS3Kt3/eOnoA4Cz2lAIAADCMKAUAAGAYUQoAAMAwohQAAIBhthylVXV1Vf1xVR2rqgeq6g2T7ZdV1T1V9eDk9NLtGxcAAIBlMs2e0seS/EJ3f3uSlyR5XVVdk+RgkiPdvSfJkcllAAAAOMuWo7S7T3T3fZPz/5XkWJIrk1yf5PDkZoeT3DDljAAAACypbXlPaVXtTvLCJB9LckV3n0hOhWuSnee4z/6qOlpVR9fW1rZjDAAAABbM1FFaVZck+WCSN3b3f272ft19qLv3dvfe1dXVaccAAABgAU0VpVX1lJwK0vd194cmmx+tql2T63clOTndiAAAACyraT59t5K8O8mx7n77uqvuTrJvcn5fkru2Ph4AAADLbGWK+748yU8m+buqun+y7S1Jbk9yZ1XdkuSRJDdONSEAAABLa8tR2t1/lqTOcfV1W/25AAAAXDy25dN3AQAAYCtEKQAAAMOIUgAAAIYRpQAAAAwjSgEAABhGlAIAADCMKAUAAGAYUQoAAMAwohQAAIBhRCkAAADDiFIAAACGEaUAAAAMI0oBAAAYRpQCAAAwjCgFAABgGFEKAADAMKIUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGGZl9AAAAMByOTB6gDPcOnoAnpA9pQAAAAwjSgEAABjG4bsAAMCSm7cDilnPnlIAAACGEaUAAAAMI0oBAAAYRpQCAAAwjCgFAABgGFEKAADAMKIUAACAYUQpAAAAw8wsSqvqVVX1qap6qKoOzupxAAAAWFwzidKq2pHk15L8cJJrktxcVdfM4rEAAABYXLPaU/riJA9196e7+8tJ7khy/YweCwAAgAW1MqOfe2WSz6y7fDzJd6+/QVXtT7J/cvGLVfWpGc0yQI0eYJYuT/K50UOMVAv59F7woYetk8V8fi5aW1wn8/EkW2uz9tV/4IX4vWM9DPck1sn2Plme+0VRyZy+nlxEa+hbznXFrKJ0o3/a/poL3YeSHJrR4zMjVXW0u/eOnoP5Zp2wGdYJm2GdsBnWCZthncyvWR2+ezzJ1esuX5XkszN6LAAAABbUrKL0b5LsqarnVtVTk9yU5O4ZPRYAAAALaiaH73b3Y1X1+iR/mGRHkvd09wOzeCwuOIdcsxnWCZthnbAZ1gmbYZ2wGdbJnKruPv+tAAAAYAZmdfguAAAAnJcoBQAAYBhRyqZU1aurqqvq+aNnYT5V1eNVdX9V/W1V3VdVLxs9E/Opqr6pqu6oqn+sqk9W1R9U1beOnov5se715IHJa8rPV5X/Z+Es69bK6f8Ojp6J+bPBOtk9eia+lveUsilVdWeSXUmOdPcvDh6HOVRVX+zuSybnX5nkLd39fYPHYs5UVSX5iySHu/udk20vSPL13f2nI2djfpzxerIzye8m+fPuvm3sZMyb9WsFzsU6mX/+6sh5VdUlSV6e5Jac+nofOJ9nJ/m30UMwl34gyVdOB2mSdPf9gpRz6e6TSfYnef3kjxoALJmZfCUMS+eGJB/u7n+oqi9U1Yu6+77RQzF3nlFV9yd5ek7tVf/BseMwp65Ncu/oIVgs3f3pyeG7O5M8Onoe5srp3z2n/VJ3/96oYZhb69fJw9396pHDcDZRymbcnOQdk/N3TC6LUs70P939giSpqpcmeW9VXdveIwBsD3tJ2chXf/fAE7BO5pwo5QlV1Tfm1B6va6uqk+xI0lX1JrHBuXT3X1bV5UlWk5wcPQ9z5YEkrxk9BIulqp6X5PF4PQFYSt5Tyvm8Jsl7u/tbunt3d1+d5OEk3zN4LubY5FOadyT5/OhZmDsfSfK0qvrp0xuq6ruqyodisaGqWk3yziS/6o+hAMvJnlLO5+Ykt5+x7YNJfiKJDyZhvfXv16gk+7r78YHzMIe6u6vq1UneMfnqhv9N8k9J3jhyLubO6deTpyR5LMlvJ3n70ImYV2e+p/TD3e1rYWDB+EoYAAAAhnH4LgAAAMOIUgAAAIYRpQAAAAwjSgEAABhGlAIAADCMKAUAAGAYUQoAAMAw/w8rLdnRmKmaIgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# barplot\n",
    "plt.figure(figsize=(16, 4))\n",
    "plt.bar(x2,height=h,width=wl/interspace,color=['darkred','darkblue'])\n",
    "plt.bar(x2,height=ha,width=wl/interspace,color=['red','blue'])\n",
    "plt.bar(x, 0, tick_label=noa.columns);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98d1c4c9",
   "metadata": {},
   "source": [
    "Take a look. What can we see here? For instance, the most of candidates applied for major A. Only a small fraction of the applicants were women. The fraction of women in major B was even smaller. What is more important, we can see that differences in acceptance rates are much greater between majors than between sexes. If we look acceptance rate for men and women for particular majors, the difference is not so big."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d313b27",
   "metadata": {},
   "source": [
    "Why do we notice such a difference in the acceptance rate for the whole university if it doesn't exist at the department level? It is simple. Women were reluctant to apply to majors with a high acceptance rate (i.e. A and B) so there were a greater number of men were accepted. They tended to apply more for the majors with low acceptance rates so the successful outcomes were not so numerous. Although we don't know what is behind the major labels, we can imagine that more women in 1973 applied for humanities than for engineering, while men tended to apply in higher numbers to engineering and science departments."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58336a0c",
   "metadata": {},
   "source": [
    "It looks like the major was the confounding variable, which had a decisive influence on the acceptance rate. Typically, men and women prefer different majors so their choice of major actually creates this picture of biased admission."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7b3b736",
   "metadata": {},
   "source": [
    "### Observational study"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebe4c5af-f1a7-4f8b-bd9d-f710ec76e5a0",
   "metadata": {},
   "source": [
    "Data obtained from the Berkeley admission office is not a result of an experiment. This data is just an observation of something that happened. All the performed analysis is an observational study. If our intention was to test if the University of California at Berkeley had discriminated against women during the admission process, we should have designed an experiment for this. In the experiment, we should prepare almost identical treatment and the control group. For instance, we should compare the acceptance rate of 100 randomly selected male candidates with 100 randomly selected female candidates. It is difficult to make experiments with the admission process. We should do what we can with observations. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e1952f1-5e0f-4e4f-9657-da809487c9b5",
   "metadata": {},
   "source": [
    "### Brief conclusion"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3537e375-b94b-42ea-9f6a-edda771acce9",
   "metadata": {},
   "source": [
    "Simpson’s paradox can have serious consequences because people easily confuse association with causal relationship and infer the wrong conclusion. We should be always cautious while interpreting observational data. Often, people don't recognize the confounders. In this example, it was easy to identify a confounder variable because it was the only available in the dataset, but in the general case, it could be challenging. Insights gained by splitting data according to factors may help us to explain observed associations."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a8d9945-3665-4a4b-aeb4-9356efd95ac2",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "### Bonus: Function for the visualization made above\n",
    "If you have a similar dataset that needs to be visualized in the same manner, here is a function that combines all previous cells. Put your DataFrame as an argument. Try, for example: `plt.figure(figsize=(16, 4)); lrp_bars(df2)` Just pay attention to the order of columns. The categorical variable that turned out to be the confounding variable is in the first column, the binary variable that separates bars into left and right parts is in the second, while the outcome was in the third.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "afa7b444",
   "metadata": {},
   "outputs": [],
   "source": [
    "def lrp_bars(f):\n",
    "    h=100\n",
    "    noa=pd.crosstab(f.iloc[:,1],f.iloc[:,0])\n",
    "    toa=noa.sum()\n",
    "    nom=len(noa.columns)\n",
    "    outcomes=f.iloc[:,2].unique()\n",
    "    nos=pd.crosstab(f.iloc[:,1][f.iloc[:,2]==outcomes[1]],f.iloc[:,0][f.iloc[:,2]==outcomes[1]])\n",
    "    ar=nos/noa\n",
    "    ha=ar.transpose().to_numpy().flatten()*100\n",
    "    w=noa/max(toa)\n",
    "    x=np.linspace(1,nom,nom)\n",
    "    interspace=1.05\n",
    "    shift=pd.DataFrame([x-(toa/max(toa)-w.iloc[0,:])/2/interspace,x-(toa/max(toa)-2*w.iloc[0,:]-w.iloc[1,:])/2/interspace])\n",
    "    x2=shift.transpose().to_numpy().flatten()\n",
    "    wl=w.transpose().to_numpy().flatten()\n",
    "    plt.bar(x2,height=h,width=wl/interspace,color=['darkred','darkblue'])\n",
    "    plt.bar(x2,height=ha,width=wl/interspace,color=['red','blue'])\n",
    "    plt.bar(x,0,tick_label=noa.columns)"
   ]
  }
 ],
 "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}