Money and death

We return to the death penalty.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
# Make plots look a little bit more fancy
plt.style.use('fivethirtyeight')

In this case, we are going to analyze whether people with higher incomes are more likely to favor the death penalty.

To do this, we are going to analyze the results from a sample of the US General Social Survey from 2002.

If you are running on your laptop, download the data file GSS2002.csv.

# Read the data into a data frame
gss = pd.read_csv('GSS2002.csv')
gss

Each row corresponds to a single respondent.

Show the column names:

gss.columns

Index(['ID', 'Region', 'Gender', 'Race', 'Education', 'Marital', 'Religion',
       'Happy', 'Income', 'PolParty', 'Politics', 'Marijuana', 'DeathPenalty',
       'OwnGun', 'GunLaw', 'SpendMilitary', 'SpendEduc', 'SpendEnv',
       'SpendSci', 'Pres00', 'Postlife'],
      dtype='object')

We want to work with only two columns from this data frame. These are “Income”, and “DeathPenalty”.

“Income” gives the income bracket of the respondent. “DeathPenalty” is the answer to a question about whether they “Favor” or “Oppose” the death penalty.

First make a list with the names of the columns that we want.

cols = ['Income', 'DeathPenalty']
cols

['Income', 'DeathPenalty']

Next make a new data frame by indexing the data frame with this list.

The new data frame has only the columns we selected.

money_death = gss[cols]
money_death

There are many missing question responses, indicated by NaN. To make our life easier, we drop the respondents who didn’t specify an income bracket, and those who did not give an answer to the death penalty question. We use Pandas dropna method of the data frame, to drop all rows that have any missing values in the row.

money_death = money_death.dropna()
money_death

Get the income column.

income = money_death['Income']

Show the unique values:

income.value_counts()

40000-49999      88
30000-34999      78
50000-59999      72
25000-29999      60
35000-39999      54
60000-74999      51
20000-22499      44
12500-14999      44
130000-149999    43
22500-24999      40
110000-129999    38
17500-19999      37
15000-17499      36
10000-124999     36
1000-2999        32
8000-9999        32
75000-89999      26
3000-3999        19
under 1000       17
5000-5999        16
4000-4999        13
90000-109999     11
7000-7999         9
6000-6999         8
Name: Income, dtype: int64

These are strings. We want to get income as a number. We estimate this by recoding the “Income” column. We replace the string, giving the income bracket, with the average of the minimum and maximum in the range.

We can do this with a recoder function. We have not covered functions yet, so do not worry about the details of this function.

def recode_income(value):
    if value == 'under 1000':
        return 500
    low_str, high_str = value.split('-')
    low, high = int(low_str), int(high_str)
    return np.mean([low, high])

Here is what the recoder function gives with the lowest income bracket.

recode_income('under 1000')

500

Here is the return from a higher bracket:

recode_income('90000-109999')

99999.5

Use this function to recode the “Income” strings into numbers. Again, we have not covered the apply method yet, so don’t worry about the details.

income_ish = income.apply(recode_income)
income_ish

0        32499.5
1        82499.5
5        44999.5
8        67499.5
13       67499.5
17       44999.5
18       82499.5
21       18749.5
22       54999.5
31       32499.5
32       54999.5
33       82499.5
35         500.0
36        7499.5
37       67499.5
42       32499.5
45       37499.5
46         500.0
52       18749.5
55       37499.5
58        1999.5
62       54999.5
64       13749.5
74      119999.5
77       82499.5
78       37499.5
81       32499.5
92       21249.5
93       67499.5
95       67499.5
          ...   
2671     82499.5
2677      1999.5
2678     16249.5
2684       500.0
2689      3499.5
2690     23749.5
2692      8999.5
2696      3499.5
2697     32499.5
2699     27499.5
2702      8999.5
2706     67499.5
2709     13749.5
2714     13749.5
2715     44999.5
2716    139999.5
2717      3499.5
2723     23749.5
2725     44999.5
2726     16249.5
2727     13749.5
2729       500.0
2737     67499.5
2741     16249.5
2742     27499.5
2747      8999.5
2753     37499.5
2756     32499.5
2760     23749.5
2764     67499.5
Name: Income, Length: 904, dtype: float64

Now get the results of the answer to the death penalty question.

death = money_death['DeathPenalty']
death.value_counts()

Favor     622
Oppose    282
Name: DeathPenalty, dtype: int64

We will identify the rows for respondents who are in favor of the death penalty. To do this, we make a Boolean vector:

death == 'Favor'

0        True
1        True
5        True
8        True
13       True
17       True
18       True
21      False
22       True
31       True
32      False
33      False
35      False
36      False
37       True
42       True
45       True
46       True
52       True
55       True
58       True
62       True
64       True
74      False
77       True
78       True
81       True
92       True
93       True
95      False
        ...  
2671     True
2677    False
2678     True
2684     True
2689     True
2690    False
2692     True
2696    False
2697     True
2699     True
2702    False
2706    False
2709    False
2714     True
2715     True
2716     True
2717    False
2723     True
2725     True
2726    False
2727     True
2729     True
2737    False
2741     True
2742     True
2747     True
2753     True
2756     True
2760     True
2764    False
Name: DeathPenalty, Length: 904, dtype: bool

Use this vector to select the income values for the respondents in favor of the death penalty. Show the distribution of values.

favor_income = income_ish[death == 'Favor']
favor_income.hist();

Likewise select incomes for those opposed. Show the distribution.

oppose_income = income_ish[death == 'Oppose']
oppose_income.hist();

Calculate the difference in mean income between the groups. This is the difference we observe.

actual_diff = np.mean(favor_income) - np.mean(oppose_income)
actual_diff

4535.163012246019

We want to know whether this difference in income is compatible with random sampling. That is, we want to know whether a difference this large is plausible, if the incomes are in fact random samples from the same population.

To estimate how variable the mean differences can be, for such random sampling, we simulate this sampling by pooling the income values that we have, from the two groups, and the permuting them.

First, we get the number of respondents in favor of the death penalty.

n_favor = len(favor_income)
n_favor

622

Then we pool the in-favor and oppose groups.

pooled = np.append(favor_income, oppose_income)

To do the random sampling we permute the values, so the pooled vector is a random mixture of the two groups.

np.random.shuffle(pooled)

Treat the first n_favor observations from this shuffled vector as our simulated in-favor group. The rest are our simulated oppose group.

fake_favor = pooled[:n_favor]
fake_oppose = pooled[n_favor:]

Calculate the difference in means for this simulation.

fake_diff = np.mean(fake_favor) - np.mean(fake_oppose)
fake_diff

3143.6351793573704

Now it is your turn. Do this simulation 10000 times, to build up the distribution of differences compatible with random sampling.

Use the Brexit ages notebook for inspiration.

differences = np.zeros(10000)
for i in np.arange(10000):
    # Permute the pooled incomes
    np.random.shuffle(pooled)
    # Make a fake favor sample

    # Make a fake opposed sample

    # Calculate the mean difference for the fake samples

    # Put the mean difference into the differences array.

When you have that working, do a histogram of the differences.

# Your code here

You can get an idea of where the actual difference we saw sits on this histogram, and therefore how likely that difference is, assuming the incomes come from the same underlying population of incomes.

To be more specific, count how many of the differences you calculated were greater than or equal to the actual difference.

# Your code here

Now calculate the proportion of these differences, to give an estimate of the probability of seeing a difference this large, if the incomes all come from the same underlying population:

# Your code here