Mendel's model¶

In [42]:
import numpy as np
from matplotlib.pyplot import subplots

# simulating the test statistic
num_repetitions = 10000
L = []
for i in range(num_repetitions):
    n = np.random.choice([1, 0], size = 929, replace=True, p = [0.75, 0.25])
    L.append(n.sum())

L = np.array(L)
L = np.abs(L/929 - 0.75)

fig, ax = subplots(figsize = (4, 4))
ax.hist(L, density = True);
ax.scatter(0.0088, 0, color='red');
No description has been provided for this image
In [43]:
# Calculating the empirical p-value for Mendel's model
empirical_pvalue = np.sum(L >= 0.0088)/num_repetitions
empirical_pvalue
Out[43]:
np.float64(0.5439)

Smoking and birth weight¶

In [44]:
import pandas as pd

baby = pd.read_csv('baby.csv')
births = baby[['Maternal Smoker', 'Birth Weight']]
births
Out[44]:
Maternal Smoker Birth Weight
0 False 120
1 False 113
2 True 128
3 True 108
4 False 136
... ... ...
1169 False 113
1170 False 128
1171 True 130
1172 False 125
1173 False 117

1174 rows × 2 columns

In [45]:
births.groupby('Maternal Smoker')['Birth Weight'].count().reset_index()
Out[45]:
Maternal Smoker Birth Weight
0 False 715
1 True 459
In [46]:
# Checking the birth weights of non smokers
births[births['Maternal Smoker'] == False]['Birth Weight']
Out[46]:
0       120
1       113
4       136
5       138
6       132
       ... 
1168    132
1169    113
1170    128
1172    125
1173    117
Name: Birth Weight, Length: 715, dtype: int64
In [27]:
fig, ax = subplots(figsize = (4, 4))
ax.hist(births[births['Maternal Smoker'] == False]['Birth Weight'], alpha=0.5, label = 'Non Smoker')
ax.hist(births[births['Maternal Smoker'] == True]['Birth Weight'], alpha=0.5, label = 'Smoker')
ax.legend();
No description has been provided for this image
In [47]:
# Finding the average weight of the groups
means = births.groupby('Maternal Smoker')['Birth Weight'].mean().reset_index()
means
Out[47]:
Maternal Smoker Birth Weight
0 False 123.085315
1 True 113.819172
In [49]:
# Finding the observed value
observed_value = means.iloc[1, 1] - means.iloc[0, 1]
observed_value
Out[49]:
np.float64(-9.266142572024918)
In [50]:
# Finding a random permutation of labels
shuffle = births.sample(frac=1.0, replace = False)['Maternal Smoker']
shuffle
Out[50]:
121     True
320    False
420    False
196    False
70     False
       ...  
865     True
35     False
185    False
132    False
969    False
Name: Maternal Smoker, Length: 1174, dtype: bool
In [39]:
ons = births.assign(Shuffled_label=shuffle.reset_index().iloc[:, 1])
ons
Out[39]:
Maternal Smoker Birth Weight Shuffled_label
0 False 120 True
1 False 113 False
2 True 128 True
3 True 108 True
4 False 136 False
... ... ... ...
1169 False 113 True
1170 False 128 True
1171 True 130 False
1172 False 125 True
1173 False 117 False

1174 rows × 3 columns

In [41]:
# Checking the test statistic for one shuffle
means = ons.groupby('Shuffled_label')['Birth Weight'].mean()
means.reset_index().iloc[0, 1] - means.reset_index().iloc[1, 1]
Out[41]:
np.float64(0.9561893444246437)
In [52]:
# Repeated shuffling to find empirical p-value
num_repetitions = 10000

L=[]
for i in range(num_repetitions):
    shuffle = births.sample(frac=1.0, replace=False)['Maternal Smoker']
    ons = births.assign(Shuffled_label = shuffle.reset_index().iloc[:, 1])
    means = ons.groupby('Shuffled_label')['Birth Weight'].mean()
    L.append(means.reset_index().iloc[0, 1] - means.reset_index().iloc[1, 1])

L = np.array(L)
L
Out[52]:
array([ 0.34090224, -0.47828816,  0.24073922, ...,  0.21569846,
        1.42838643,  0.02610418])
In [53]:
fig, ax = subplots(figsize = (4, 4))
ax.hist(L, density=True);
ax.set_xlabel('Difference between group means')
ax.set_title('Prediction under null hypothesis')
Out[53]:
Text(0.5, 1.0, 'Prediction under null hypothesis')
No description has been provided for this image
In [55]:
# Calculating empirical p-value
empirical_pvalue = np.sum(L <= observed_value)
empirical_pvalue
Out[55]:
np.int64(0)

RCT for lower back pain¶

In [56]:
bta = pd.read_csv('bta.csv')
bta
Out[56]:
Group Result
0 Control 1.0
1 Control 1.0
2 Control 0.0
3 Control 0.0
4 Control 0.0
5 Control 0.0
6 Control 0.0
7 Control 0.0
8 Control 0.0
9 Control 0.0
10 Control 0.0
11 Control 0.0
12 Control 0.0
13 Control 0.0
14 Control 0.0
15 Control 0.0
16 Treatment 1.0
17 Treatment 1.0
18 Treatment 1.0
19 Treatment 1.0
20 Treatment 1.0
21 Treatment 1.0
22 Treatment 1.0
23 Treatment 1.0
24 Treatment 1.0
25 Treatment 0.0
26 Treatment 0.0
27 Treatment 0.0
28 Treatment 0.0
29 Treatment 0.0
30 Treatment 0.0
In [57]:
means = bta.groupby('Group')['Result'].mean()
means
Out[57]:
Group
Control      0.125
Treatment    0.600
Name: Result, dtype: float64
In [58]:
obout = pd.read_csv('observed_outcomes.csv')
obout
Out[58]:
Group Outcome if assigned treatment Outcome if assigned control
0 Control Unknown 1
1 Control Unknown 1
2 Control Unknown 0
3 Control Unknown 0
4 Control Unknown 0
5 Control Unknown 0
6 Control Unknown 0
7 Control Unknown 0
8 Control Unknown 0
9 Control Unknown 0
10 Control Unknown 0
11 Control Unknown 0
12 Control Unknown 0
13 Control Unknown 0
14 Control Unknown 0
15 Control Unknown 0
16 Treatment 1 Unknown
17 Treatment 1 Unknown
18 Treatment 1 Unknown
19 Treatment 1 Unknown
20 Treatment 1 Unknown
21 Treatment 1 Unknown
22 Treatment 1 Unknown
23 Treatment 1 Unknown
24 Treatment 1 Unknown
25 Treatment 0 Unknown
26 Treatment 0 Unknown
27 Treatment 0 Unknown
28 Treatment 0 Unknown
29 Treatment 0 Unknown
30 Treatment 0 Unknown
In [59]:
bta
Out[59]:
Group Result
0 Control 1.0
1 Control 1.0
2 Control 0.0
3 Control 0.0
4 Control 0.0
5 Control 0.0
6 Control 0.0
7 Control 0.0
8 Control 0.0
9 Control 0.0
10 Control 0.0
11 Control 0.0
12 Control 0.0
13 Control 0.0
14 Control 0.0
15 Control 0.0
16 Treatment 1.0
17 Treatment 1.0
18 Treatment 1.0
19 Treatment 1.0
20 Treatment 1.0
21 Treatment 1.0
22 Treatment 1.0
23 Treatment 1.0
24 Treatment 1.0
25 Treatment 0.0
26 Treatment 0.0
27 Treatment 0.0
28 Treatment 0.0
29 Treatment 0.0
30 Treatment 0.0
In [101]:
num_repetitions = 20000
L = []
for i in range(num_repetitions):
    shuffle = bta.sample(frac=1.0, replace=False)['Group']
    shuffled_bta = bta.assign(Shuffled_groups=shuffle.reset_index().iloc[:, 1])
    means = shuffled_bta.groupby('Shuffled_groups')['Result'].mean().reset_index()
    L.append(means.iloc[0, 1] - means.iloc[1, 1])
L = np.array(L)
L = np.abs(L)
In [102]:
fig, ax = subplots(figsize = (4, 4))
ax.hist(L, bins = np.arange(0, 0.7, 0.1));
ax.scatter(0.475, 0, color='red');
No description has been provided for this image
In [103]:
empirical_pvalue = np.sum(L >= 0.475)/num_repetitions
empirical_pvalue
Out[103]:
np.float64(0.00905)

Deflategate¶

In [104]:
balls = pd.read_csv('deflategate.csv')
balls
Out[104]:
Team Blakeman Prioleau
0 Patriots 11.50 11.80
1 Patriots 10.85 11.20
2 Patriots 11.15 11.50
3 Patriots 10.70 11.00
4 Patriots 11.10 11.45
5 Patriots 11.60 11.95
6 Patriots 11.85 12.30
7 Patriots 11.10 11.55
8 Patriots 10.95 11.35
9 Patriots 10.50 10.90
10 Patriots 10.90 11.35
11 Colts 12.70 12.35
12 Colts 12.75 12.30
13 Colts 12.50 12.95
14 Colts 12.55 12.15
In [105]:
average_balls = balls.assign(Combined = (balls['Blakeman']+balls['Prioleau'])/2)
average_balls
Out[105]:
Team Blakeman Prioleau Combined
0 Patriots 11.50 11.80 11.650
1 Patriots 10.85 11.20 11.025
2 Patriots 11.15 11.50 11.325
3 Patriots 10.70 11.00 10.850
4 Patriots 11.10 11.45 11.275
5 Patriots 11.60 11.95 11.775
6 Patriots 11.85 12.30 12.075
7 Patriots 11.10 11.55 11.325
8 Patriots 10.95 11.35 11.150
9 Patriots 10.50 10.90 10.700
10 Patriots 10.90 11.35 11.125
11 Colts 12.70 12.35 12.525
12 Colts 12.75 12.30 12.525
13 Colts 12.50 12.95 12.725
14 Colts 12.55 12.15 12.350
In [107]:
L = np.ones(11)*12.5
L = np.append(L, np.ones(4)*13)
L
Out[107]:
array([12.5, 12.5, 12.5, 12.5, 12.5, 12.5, 12.5, 12.5, 12.5, 12.5, 12.5,
       13. , 13. , 13. , 13. ])
In [109]:
balls2 = average_balls.assign(Drop = L - average_balls['Combined'])
balls2
Out[109]:
Team Blakeman Prioleau Combined Drop
0 Patriots 11.50 11.80 11.650 0.850
1 Patriots 10.85 11.20 11.025 1.475
2 Patriots 11.15 11.50 11.325 1.175
3 Patriots 10.70 11.00 10.850 1.650
4 Patriots 11.10 11.45 11.275 1.225
5 Patriots 11.60 11.95 11.775 0.725
6 Patriots 11.85 12.30 12.075 0.425
7 Patriots 11.10 11.55 11.325 1.175
8 Patriots 10.95 11.35 11.150 1.350
9 Patriots 10.50 10.90 10.700 1.800
10 Patriots 10.90 11.35 11.125 1.375
11 Colts 12.70 12.35 12.525 0.475
12 Colts 12.75 12.30 12.525 0.475
13 Colts 12.50 12.95 12.725 0.275
14 Colts 12.55 12.15 12.350 0.650
In [110]:
means = balls2.groupby('Team')['Drop'].mean().reset_index()
means
Out[110]:
Team Drop
0 Colts 0.468750
1 Patriots 1.202273
In [122]:
num_repetitions = 30000
L = []
for i in range(num_repetitions):
    shuffle = balls2.sample(frac=1.0, replace=False)['Team']
    shuffled_balls = balls2.assign(Shuffled_teams=shuffle.reset_index().iloc[:, 1])
    means = shuffled_balls.groupby('Shuffled_teams')['Drop'].mean().reset_index()
    L.append(means.iloc[1, 1] - means.iloc[0, 1])

L = np.array(L)
In [123]:
fig, ax = subplots(figsize = (4, 4))
ax.hist(L);
ax.scatter(0.73352, 0, color='red');
No description has been provided for this image
In [124]:
empirical_pvalue = np.sum(L >= 0.73352)/num_repetitions
empirical_pvalue
Out[124]:
np.float64(0.0036666666666666666)
In [ ]: