<a href="https://colab.research.google.com/github/dlsun/pods/blob/master/02-Categorical-Data/2_3_Multi_Way_Tables_and_Simpson's_Paradox.ipynb" target="_parent"></a>

## 2.3 Multi-Way Tables and Simpson's Paradox

In the previous lesson, we summarized two categorical variables by cross-tabulating their frequencies.


```python
import pandas as pd
data_dir = "http://dlsun.github.io/pods/data/"
df_titanic = pd.read_csv(data_dir + "titanic.csv")

def class_to_type(c):
  if c in ["1st", "2nd", "3rd"]:
    return "passenger"
  else:
    return "crew"
df_titanic["type"] = df_titanic["class"].map(class_to_type)

joint_type_survived = pd.crosstab(
    df_titanic["type"],
    df_titanic["survived"],
    normalize=True
)
joint_type_survived
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th>survived</th>
      <th>0</th>
      <th>1</th>
    </tr>
    <tr>
      <th>type</th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>crew</th>
      <td>0.307657</td>
      <td>0.095605</td>
    </tr>
    <tr>
      <th>passenger</th>
      <td>0.370186</td>
      <td>0.226552</td>
    </tr>
  </tbody>
</table>
</div>



Each number in this table represents a joint probability. For example:
$$ P(\text{crew}, \text{survived}) = 0.095605. $$

We might want to know whether crew members or passengers survived at higher rates. To do this, we have to compare the conditional probabilities
$$\begin{align}
P(\text{survived} | \text{crew}) & & \text{vs.} & & P(\text{survived} | \text{passenger}).
\end{align}$$

In the last lesson, we learned to calculate conditional distributions using broadcasting.


```python
survived_given_type = joint_type_survived.divide(
    joint_type_survived.sum(axis=1),
    axis=0
)
survived_given_type
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th>survived</th>
      <th>0</th>
      <th>1</th>
    </tr>
    <tr>
      <th>type</th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>crew</td>
      <td>0.762921</td>
      <td>0.237079</td>
    </tr>
    <tr>
      <td>passenger</td>
      <td>0.620349</td>
      <td>0.379651</td>
    </tr>
  </tbody>
</table>
</div>



From the table, it is apparent that passengers survived at much higher rates than crew members:
$$ P(\text{survived}|\text{crew}) = 0.237079 < 0.379651 = P(\text{survived}|\text{passenger}). $$



### 2.3.1 Communication Corner: Reporting Differences in Probabilities

How do we report the difference between the two probabilities $23.71\%$ and $37.97\%$ above? There are a number of ways:

1. As an _additive change_: "Passengers were 
$$ 37.97\% - 23.71\% = 14.26 \text{ percentage points} $$
more likely to survive than crew members."
2. As a _relative change_ (or _relative risk_): "Passengers were 
$$ 37.97\% \big/ 23.71\% = 1.60 \text{ times} $$
as likely to survive as crew members." 
3. We can translate relative changes to _percent changes_ by subtracting $1$ and multiplying by $100\%$. So we can rephrase the above as: "Passengers were 
$$ 100\% \times (1.60 - 1.00) = 60\% $$
more likely to survive than crew members."
4. As an _odds ratio_: "The odds of a passenger surviving was 
$$ \frac{37.97\% \big/ (100\% - 37.97\%)}{23.71\% \big/ (100\% - 23.71\%)} = 1.97 \text{ times} $$
greater than the odds of a crew member surviving."

Note that additive changes and percent changes should be compared to a baseline of 0.0 (and can be negative), while relative changes and odds ratios should be compared to a baseline of 1.0 (and cannot be negative).

Watch out: it is incorrect to say that passengers are $14.26\%$ more likely to survive than crew members, since the percent change is $60\%$. An additive change should always be reported in units of "percentage points".

### 2.3.2 Controlling for a Variable

But is this the whole story? We know that survival rates for males and females were very different. Will the trend between the survival rates for crew and passengers still hold after we _control_ for **gender**? 

To do this, let's determine the joint distribution of these two variables and a third variable, **gender**. In principle, the frequencies could be represented using a three-dimensional table, but it is difficult to visualize more than two dimensions on paper or on a screen. So we put two of the variables along one dimension and one variable along the other, creating a _three-way table_.


```python
joint_gender_type_survived = pd.crosstab(
    [df_titanic["gender"], df_titanic["type"]],
    df_titanic["survived"],
    normalize=True
)
joint_gender_type_survived
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>survived</th>
      <th>0</th>
      <th>1</th>
    </tr>
    <tr>
      <th>gender</th>
      <th>type</th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td rowspan="2" valign="top">female</td>
      <td>crew</td>
      <td>0.001359</td>
      <td>0.009062</td>
    </tr>
    <tr>
      <td>passenger</td>
      <td>0.057544</td>
      <td>0.153602</td>
    </tr>
    <tr>
      <td rowspan="2" valign="top">male</td>
      <td>crew</td>
      <td>0.306298</td>
      <td>0.086543</td>
    </tr>
    <tr>
      <td>passenger</td>
      <td>0.312642</td>
      <td>0.072950</td>
    </tr>
  </tbody>
</table>
</div>



Of course, we would have chosen any two of the variables to place along the rows, or had the two variables along the columns instead of the rows. The particular representation above was chosen because it makes it easy to survival rates for each gender and type, i.e.,
$$ P(\text{survived} | \textbf{gender}, \textbf{type}), $$
where **gender** is either "male" or "female" and **type** is either "crew" or "passenger". Recall that the conditional probability is calculated as 
$$ P(\text{survived} | \textbf{gender}, \textbf{type}) = \frac{P(\text{survived}, \textbf{gender}, \textbf{type})}{P(\textbf{gender}, \textbf{type})}. $$
The numerator is the joint distribution above. The denominator can be calculated by summing over the possible values of **survived**---in other words, across each row, over the columns. 


```python
joint_gender_type = joint_gender_type_survived.sum(axis=1)
joint_gender_type
```




    gender  type     
    female  crew         0.010421
            passenger    0.211146
    male    crew         0.392841
            passenger    0.385591
    dtype: float64



To obtain the conditional probabilities, we simply divide the joint distribution by the marginal.


```python
survived_given_gender_type = joint_gender_type_survived.divide(
    joint_gender_type,
    axis=0
)
survived_given_gender_type
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>survived</th>
      <th>0</th>
      <th>1</th>
    </tr>
    <tr>
      <th>gender</th>
      <th>type</th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td rowspan="2" valign="top">female</td>
      <td>crew</td>
      <td>0.130435</td>
      <td>0.869565</td>
    </tr>
    <tr>
      <td>passenger</td>
      <td>0.272532</td>
      <td>0.727468</td>
    </tr>
    <tr>
      <td rowspan="2" valign="top">male</td>
      <td>crew</td>
      <td>0.779700</td>
      <td>0.220300</td>
    </tr>
    <tr>
      <td>passenger</td>
      <td>0.810811</td>
      <td>0.189189</td>
    </tr>
  </tbody>
</table>
</div>



Now, let's compare the survival rates of passengers and crew members for females and males separately. 

- For females, crew members survived at a higher rate:
$$ P(\text{survived} | \text{female}, \text{crew}) = 0.869565 > 0.727468 = P(\text{survived} | \text{female}, \text{passenger}) $$
- For males, crew members survived at a higher rate:
$$ P(\text{survived} | \text{male}, \text{crew}) = 0.220300 > 0.189189 = P(\text{survived} | \text{male}, \text{passenger}) $$

But remember, we found earlier that passengers survived at a higher rate overall:
$$ P(\text{survived} | \text{crew}) < P(\text{survived} | \text{passenger}). $$

How is it possible that both male and female crew members survived at a higher rate, yet crew members survived at a lower rate overall? This surprising phenomenon is known as _Simpson's paradox_.

### 2.3.3 Simpson's Paradox

Simpson's paradox is a phenomenon where a trend disappears or reverses when the data is aggregated. In the Titanic data set, both male and female crew members survived at higher rates, but when we aggregated over gender, the trend reversed.

In order to investigate Simpson's paradox, we first reorganize the probabilities. First, we keep only the survival rate, dropping the death rate (since it is redundant; it is just one minus the survival rate).


```python
survived_given_gender_type[1]
```




    gender  type     
    female  crew         0.869565
            passenger    0.727468
    male    crew         0.220300
            passenger    0.189189
    Name: 1, dtype: float64



Next, we rearrange these probabilities into a two-way table, with gender along one dimension and type along the other. This can be achieved in code by "unstacking" a level of the index. (There are two "levels": **gender** and **type**.)


```python
survival_rates_by_gender_type = (survived_given_gender_type[1].
                                 unstack(level="type"))
survival_rates_by_gender_type
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th>type</th>
      <th>crew</th>
      <th>passenger</th>
    </tr>
    <tr>
      <th>gender</th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>female</td>
      <td>0.869565</td>
      <td>0.727468</td>
    </tr>
    <tr>
      <td>male</td>
      <td>0.220300</td>
      <td>0.189189</td>
    </tr>
  </tbody>
</table>
</div>



Caution: the probabilities in this table do not represent a distribution. They do not add up to 1.0. These probabilities originally came from the conditional distribution of **survived** given **gender** and **type**, but we dropped the death rates from the data.

Finally, we append the overall survival rates (which we calculated at the beginning of this lesson) to the last row of this `DataFrame`.


```python
survival_rates_by_gender_type.append(survived_given_type[1])
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th>type</th>
      <th>crew</th>
      <th>passenger</th>
    </tr>
    <tr>
      <th>gender</th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>female</td>
      <td>0.869565</td>
      <td>0.727468</td>
    </tr>
    <tr>
      <td>male</td>
      <td>0.220300</td>
      <td>0.189189</td>
    </tr>
    <tr>
      <td>1</td>
      <td>0.237079</td>
      <td>0.379651</td>
    </tr>
  </tbody>
</table>
</div>



The overall survival rates are weighted averages of the survival rates for each gender. If we look at the survival rates for crew members:

- The survival rate for female crew is 87.0%.
- The survival rate for male crew is 22.0%.
- The overall survival rate for all crew is 23.7%, which is between the gender-specific survival rates, but much closer to the survival rate for male crew.

Likewise, if we look at the survival rates for passengers:

- The survival rate for female crew is 72.7%.
- The survival rate for male crew is 18.9%.
- The overall survival rate for all crew is 38.0%, which is closer to the middle of the gender-specific survival rates.

Why would the survival rate for crew members be so close to the survival rate for male crew? To answer this question, let's examine the weights that go into this weighted average.

In mathematical notation, the overall survival rate can be decomposed as:
$$ \underbrace{P(\text{survived} | \textbf{type})}_{\text{overall survival rate}} = \sum_{\textbf{gender}} \underbrace{P(\textbf{gender} | \textbf{type})}_{\text{weight}} \underbrace{P(\text{survived} | \textbf{gender}, \textbf{type})}_{\text{gender-specific survival rate}}. $$
So we see that the weights are $P(\textbf{gender} | \textbf{type})$. 

First, we calculate this conditional distribution from the joint distribution of **gender** and **type**.


```python
joint_gender_type = pd.crosstab(
    df_titanic["gender"],
    df_titanic["type"],
    normalize=True
)

gender_given_type = joint_gender_type.divide(
    joint_gender_type.sum(axis=0),
    axis=1
)
gender_given_type
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th>type</th>
      <th>crew</th>
      <th>passenger</th>
    </tr>
    <tr>
      <th>gender</th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>female</td>
      <td>0.025843</td>
      <td>0.353834</td>
    </tr>
    <tr>
      <td>male</td>
      <td>0.974157</td>
      <td>0.646166</td>
    </tr>
  </tbody>
</table>
</div>



Notice that 97.4% of crew members were male! So the lower male survival rate is going to dominate the weighted average when we calculate the overall survival rate for crew members. On the other hand, the gender ratio for passengers was more balanced, so their overall survival rate will end up being closer to the middle of the male and female survival rates.

Now, we calculate the weighted average, using the conditional distribution of gender as "weights" that we multiply by the survival rates. Then, we sum over the genders to get the weighted averages---i.e., the overall survival rates.


```python
(gender_given_type * survival_rates_by_gender_type).sum(axis=0)
```




    type
    crew         0.237079
    passenger    0.379651
    dtype: float64



Check that these match the overall survival rates that we calculated above.

So the secret of Simpson's Paradox lies in two facts:

1. Survival rates were generally much lower for men than for women.
2. Because crew members were predominantly male, their survival rate was weighted towards the lower male survival rate, that their overall survival rate ended up being lower than the survival rate for passengers.

Simpson's Paradox means that we have to be careful when comparing proportions from a two-way table, such as survival rates for crew and passengers. When we control for a third variable, such as **gender**, the direction of the effect could change!

### 2.3.4 Exercises

1\. Calculate the _percent change_ in survival rates between passengers and crew members, controlling for where they embarked. Does there appear to be a Simpson's paradox effect here?

Exercise 2 asks you to work with the Florida Death Penalty data set, which is available at  `https://dlsun.github.io/pods/data/death_penalty.csv`. This data set contains information about the races of the defendant and the victim, as well as whether a death penalty verdict was rendered, in 674 homicide trials in Florida between 1976-1987.

2\. Calculate the _relative risk_ of a death penalty verdict for black defendants (relative to white defendants), adjusting for the race of the victim. How does this compare to what you found at the end of the last chapter? Can you explain intuitively what is going on?

ab2b31b86ca21938051d93f8c20b294fcec5434d
ipynb
Jupyter Notebook
TesterBook/02_Categorical_Data/.ipynb_checkpoints/2.3 Multi Way Tables and Simpson's Paradox-checkpoint.ipynb
bfkwong/Latex-Book-Parser
07ed107b785c7297eb4a41bbf337bc006a14f00c
MIT
TesterBook/02_Categorical_Data/.ipynb_checkpoints/2.3 Multi Way Tables and Simpson's Paradox-checkpoint.ipynb
bfkwong/Latex-Book-Parser
07ed107b785c7297eb4a41bbf337bc006a14f00c
MIT
2020-07-06T00:12:10.000Z
2020-07-06T01:53:58.000Z
TesterBook/02_Categorical_Data/2.3 Multi Way Tables and Simpson's Paradox.ipynb
bfkwong/Latex-Book-Parser
07ed107b785c7297eb4a41bbf337bc006a14f00c
MIT
Qwen/Qwen-72B
1. YES
2. YES