# Zip Code Plotting

## Beginning

This is a quick test to see if I can plot some GeoJSON that defines a zipcode map for the state of Oregon. I got the file from GitHub but the README indicates that the original source (see below) was the the U.S. Census Bureau. The files are from 2010, so they're a little old, but I don't imagine that zip codes update that often anyway. Maybe as a future exercise I'll try to replicated what was done with the 2019 update.

### Imports

# python
from functools import partial
from pathlib import Path

import os

# pypi
from dotenv import load_dotenv

import geopandas
import hvplot.pandas
import matplotlib.pyplot as pyplot
import pandas

# my stuff
from graeae import EmbedHoloviews


### Set Up

SLUG = "zip-code-plotting"
Embed = partial(EmbedHoloviews, folder_path=f"files/posts/{SLUG}")


## Middle

This should be pretty straightforward. The documentation for GeoPandas indicates that it can handle GeoJSON directly.

load_dotenv("posts/.env")
path = Path(os.environ["OREGON_ZIP_CODES"]).expanduser()


Let's see what's there.

print(geoframe.head())

  STATEFP10 ZCTA5CE10  GEOID10 CLASSFP10 MTFCC10 FUNCSTAT10    ALAND10  \
0        41     97833  4197833        B5   G6350          S  228152974
1        41     97840  4197840        B5   G6350          S  295777905
2        41     97330  4197330        B5   G6350          S  199697439
3        41     97004  4197004        B5   G6350          S  113398767
4        41     97023  4197023        B5   G6350          S  330220870

AWATER10   INTPTLAT10    INTPTLON10 PARTFLG10  \
0         0  +44.9288886  -118.0148791         N
1  10777783  +44.8847111  -116.9184395         N
2    814864  +44.6424890  -123.2562655         N
3     71994  +45.2549625  -122.4493774         N
4   2345079  +45.2784758  -122.3231876         N

geometry
0  MULTIPOLYGON (((-118.15748 44.99903, -118.1575...
1  POLYGON ((-116.98971 44.88256, -116.98957 44.8...
2  POLYGON ((-123.18294 44.64477, -123.18293 44.6...
3  POLYGON ((-122.48691 45.22227, -122.48713 45.2...
4  POLYGON ((-122.07580 45.10889, -122.07594 45.1...


This next bit is a little slow.

plot = geoframe.hvplot(hover_cols=["ZCTA5CE10"], legend=False, tools=["hover", "wheel_zoom"],).opts(
title="Oregon Zip Codes",
width=800,
height=700,
fontscale=2,
)
outcome = Embed(plot=plot, file_name="oregon_zip_codes")()

print(outcome)


Well, on the one hand it kind of works, but it has a lot of weird holes and I can't get the tools to appear so you can zoom in. I also don't really want the zip-codes to be interpreted as a heat map. But I guess it's a start.

# The Geopandas DataFrame Coordinates Example

## Beginning

### Imports

# python
from functools import partial

# pypi
import geopandas
import hvplot.pandas
import matplotlib.pyplot as pyplot
import pandas

# my stuff
from graeae import EmbedHoloviews


### Set Up

Embed = partial(EmbedHoloviews, folder_path="files/posts/the-geopandas-dataframe-coordinates-example")


## Middle

### Longitudes and Latitudes

data = pandas.DataFrame(
{'City': ['Buenos Aires', 'Brasilia', 'Santiago', 'Bogota', 'Caracas'],
'Country': ['Argentina', 'Brazil', 'Chile', 'Colombia', 'Venezuela'],
'Latitude': [-34.58, -15.78, -33.45, 4.60, 10.48],
'Longitude': [-58.66, -47.91, -70.66, -74.08, -66.86]})

geoframe = geopandas.GeoDataFrame(data,
geometry=geopandas.points_from_xy(
data.Longitude, data.Latitude))

print(geoframe.head())

           City    Country  Latitude  Longitude                     geometry
0  Buenos Aires  Argentina    -34.58     -58.66  POINT (-58.66000 -34.58000)
1      Brasilia     Brazil    -15.78     -47.91  POINT (-47.91000 -15.78000)
2      Santiago      Chile    -33.45     -70.66  POINT (-70.66000 -33.45000)
3        Bogota   Colombia      4.60     -74.08    POINT (-74.08000 4.60000)
4       Caracas  Venezuela     10.48     -66.86   POINT (-66.86000 10.48000)


### Plot It

# read in the plotting data

# restrict it to South America.
ax = world[world.continent == 'South America'].plot(
color='white', edgecolor='black')

# We can now plot our GeoDataFrame.
geoframe.plot(ax=ax, color='red')

pyplot.savefig("files/posts/the-geopandas-dataframe-coordinates-example/south_america.png")


### With Holoviews

Since this is going to be interactive we can add more to the plot.

print(world.head())

     pop_est      continent                      name iso_a3  gdp_md_est  \
0     920938        Oceania                      Fiji    FJI      8374.0
1   53950935         Africa                  Tanzania    TZA    150600.0
2     603253         Africa                 W. Sahara    ESH       906.5
3   35623680  North America                    Canada    CAN   1674000.0
4  326625791  North America  United States of America    USA  18560000.0

geometry
0  MULTIPOLYGON (((180.00000 -16.06713, 180.00000...
1  POLYGON ((33.90371 -0.95000, 34.07262 -1.05982...
2  POLYGON ((-8.66559 27.65643, -8.66512 27.58948...
3  MULTIPOLYGON (((-122.84000 49.00000, -122.9742...
4  MULTIPOLYGON (((-122.84000 49.00000, -120.0000...


Let's add the population estimate to the countries and the city to the points.

world_plot = world[world.continent == "South America"].hvplot(
color="white",
x="Longitude",
y="Latitude",
legend=False,
hover_cols=["pop_est", "name"])

# NOTE: the overlay uses the pandas.DataFrame, not the GeoDataFrame
points = data.hvplot.points("Longitude", "Latitude", hover_cols=["City"])
plot = (world_plot * points).opts(
title="South American Cities",
width=700,
height=890,
)
outcome = Embed(plot=plot, file_name="south_america", height_in_pixels=1000)()

print(outcome)


This was mostly a check to see if I could get it working so I'll have to think about adding more explanations later.

### Notes

• To get matplotlib to work I had to install descartes (pip install descartes)
• To get hvplot to work I had to install geoviews (pip install geoviews)

# Secrets of Mental Math

## Citation

• Benjamin A, Shermer M, Benjamin A, 3M Cloud Library. Secrets of mental math: the mathemagician’s guide to lightning calculation and amazing math tricks [Internet]. Place of publication not identified: Crown Publishing Group; 2006 [cited 2020 Nov 9].

# Bibliography for Think Complexity

## Citation

1. Downey A. Think complexity: complexity science and computational modeling. Second edition. Beijing ; Boston: O’Reilly; 2018. 185 p.

This is a creative commons book that is for sale but is also available to read online or as a free PDF from Green Tea Press. There is also a github repository here.

# CPR and Blood Thinners

## Beginning

We want to test the hypothesis that blood thinners administered to patients that have been give CPR has some kind of effect (for good or for ill) on survival rates.

• $H_0$: Bood thinners don't have an effect on survival rate - the proportions will be the same for each group ($p_{treatment} = p_{control}$).
• $H_A$: Blood thinners impact survival, either positive or negative ($p_{treatment} - p_{control} \neq 0$).

We're going to use a significance level of 0.05.

ALPHA = 0.05


### Set Up

#### Imports

from argparse import Namespace
from functools import partial

# pypi
from numpy.random import default_rng
from tabulate import tabulate

import holoviews
import hvplot.pandas
import numpy
import pandas

# my stuff
from graeae import EmbedHoloviews, Timer


#### Random Generator

numpy_random = default_rng()


#### The Timer

TIMER = Timer()


#### Plotting

SLUG = "cpr-and-blood-thinners"
path = f"files/posts/{SLUG}"
Embed = partial(EmbedHoloviews, folder_path=path)

Plot = Namespace(
width=990,
height=780,
fontscale=2,
tan="#ddb377",
blue="#4687b7",
red="#ce7b6d",
color_cycle = holoviews.Cycle(["#4687b7", "#ce7b6d"])
)


#### The Table Printer

TABLE = partial(tabulate, headers="keys", tablefmt="orgtbl")


## Middle

### The Study Data

data = pandas.DataFrame.from_dict(
dict(
Survived=[11, 14],
Died = [39, 26]))
data.index = ["Control", "Treatment"]
row = data.sum()
row.name = "Total"
data = data.append(row)
data["Total"] = data.sum(axis="columns")
print(TABLE(data))

Survived Died Total
Control 11 39 50
Treatment 14 26 40
Total 25 65 90

The fifty people in the Control Group didn't receive blood thinners and the forty in the Treatment Group did.

Now I'll plot it to see what the differences look like.

# put the groups into a column
plotter = data.reset_index()

# get rid of the totals
del(plotter["Total"])
plotter = plotter.iloc[:-1]

plotter = plotter.melt(id_vars=["index"],
value_vars=["Survived", "Died"])

plotter = plotter.rename(columns=dict(
index="Group",
variable="Outcome",
value="Count"))
plot = plotter.hvplot.bar(x="Group", y="Count", by="Outcome").opts(
title="Survived or Died",
width=Plot.width,
height=Plot.height,
color=Plot.color_cycle,
fontscale=Plot.fontscale,
)

outcome = Embed(plot=plot, file_name="survived_or_died")()

print(outcome)


Hard to say, but it looks like more of the Treatment Group survived.

#### The Porportions

We are interested in the difference in survival rates so we need to calculate the proprtion that survived in each group and take the difference.

survival_rates = (data["Survived"] / data.Total).reset_index().rename(
columns={"index": "Group", 0: "Survival Rate"})
print(TABLE(survival_rates, showindex=False))

Group Survival Rate
Control 0.22
Treatment 0.35
Total 0.277778
survival_rates = survival_rates.set_index("Group")
TEST_STATISTIC = survival_rates.loc["Treatment"] - survival_rates.loc["Control"]
print(f"Difference: {TEST_STATISTIC}")

Difference: Survival Rate    0.13
dtype: float64


It appears that 13% more people survive when they are given blood thinners than when they aren't.

### Simulation

Now I'll run a simulation to build up the null distribution of the differences.

#### Set Up The Urn

The urn will be based on the total row on the contingency table. The that number survived will be represented as ones and the number that died will be represented as zeros.

total = data.loc["Total"]
survived = numpy.ones(total.Survived)
died = numpy.zeros(total.Died)
urn = numpy.append(survived, died)


Now we'll run the simulation to get the distribution.

TRIALS = 10**6
CONTROL_GROUP = data.Total.loc["Control"]
TREATMENT_GROUP = data.Total.loc["Treatment"]
survived = data["Survived"]

with TIMER:
controls = numpy.array([numpy_random.choice(urn,
CONTROL_GROUP,
replace=False).sum()
for trial in range(TRIALS)])
treatments = survived.loc["Total"] - controls

control_proportions = controls/CONTROL_GROUP
treatment_proportions = treatments/TREATMENT_GROUP

differences = treatment_proportions - control_proportions
differences = pandas.DataFrame.from_dict({"Point Estimate": differences})

2020-10-01 19:38:55,370 graeae.timers.timer start: Started: 2020-10-01 19:38:55.370038
2020-10-01 19:39:31,437 graeae.timers.timer end: Ended: 2020-10-01 19:39:31.437364
2020-10-01 19:39:31,438 graeae.timers.timer end: Elapsed: 0:00:36.067326


Plot the distribution.

plot = differences.hvplot.hist("Point Estimate").opts(
title="Null Distribution of Point Estimate",
width=Plot.width,
height=Plot.height,
color=Plot.tan,
fontscale=Plot.fontscale,
)
outcome = Embed(plot=plot, file_name="null_distribution")()

print(outcome)


Now we want to try a two-sided hypothesis test.

right_side = (differences >= TEST_STATISTIC.values).sum()
total = right_side * 2
p_value = total / len(differences)
print(f"One-Sided p-value: {(right_side/len(differences)).values[0]}")
print(f"Two-sided p-value: {p_value.values[0]}")

One-Sided p-value: 0.129026
Two-sided p-value: 0.258052

print(f"We reject the Null Hypothesis: {p_value.values[0] < ALPHA}")

We reject the Null Hypothesis: False


Even a single-sided test wouldn't have provided enough evidence to reject the null hypothesis, but with the two-sided version it looks really like the difference was from chance.

## End

Although it initially looked like the use of blood thinners after someone has been given CPR helped survival rates our data does not contain evidence to say that it does.

# Opportunity Cost Neglect

## Citation

• Frederick S, Novemsky N, Wang J, Dhar R, Nowlis S. Opportunity cost neglect. Journal of Consumer Research. 2009 Dec 1;36(4):553-61.

## Abstract

### And

To properly consider the opportunity costs of a purchase, consumers must actively generate the alternatives that it would displace.

### But

The current research suggests that consumers often fail to do so. Even under conditions promoting cognitive effort, various cues to consider opportunity costs reduce purchase rates and increase the choice share of more affordable options. Sensitivity to such cues varies with chronic dispositional differences in spending attitudes.

### Therefore

We discuss the implications of these results for the marketing strategies of economy and premium brands.

# Opportunity Cost

## Beginning

We want to answer the question of whether reminding college students that spending money now deprives them of money in the future will affect their spending behavior.

### Imports

# python
from argparse import Namespace
from functools import partial

# pypi
from expects import (
equal,
expect,
)
from numpy.random import default_rng
from tabulate import tabulate

import holoviews
import hvplot.pandas
import numpy
import pandas

# my stuff
from graeae import EmbedHoloviews, Timer


### Set Up

#### Randomness

numpy_random = default_rng()


#### Table

TABLE = partial(tabulate, tablefmt="orgtbl", headers="keys")


#### Timer and Plotting

TIMER = Timer()

SLUG = "opportunity-cost"
Embed = partial(EmbedHoloviews,
folder_path=f"files/posts/{SLUG}")

Plot = Namespace(
width=990,
height=780,
fontscale=2,
tan="#ddb377",
blue="#4687b7",
red="#ce7b6d",
color_cycle = holoviews.Cycle(["#4687b7", "#ce7b6d"])
)


## Middle

### Hypotheses

• $H_0:$ Null Hypothesis. Reminding students that they can save money for future purchases won't impact their spending decisions.
• $H_A:$ Alternative Hypothesis. Reminding students that they can save money for future purchases will reduce their spending.

I'll us a 95% significance level.

ALPHA = "0.05"


### The Study Data

150 Students were given the following statement:

Imagine that you have been saving some extra money on the side to make some purchases and on your most recent visit to the video store you come across a special sale on a new video. This video is one with your favorite actor or actress, and your favorite type of movie (such as comedy, drama, thriller, etc.). This particular video that you are considering is one you have been thinking about buying for a long time. It is available for a special sale price of $14.99. What would you do in this situation? Please circle one of the options below. Video store? Are we going to have to invent a time machine to run this study? No, we'll just swipe the data published by Frederick et al. in 2009. They conducted their study at Arizona State University, which perhaps still had some DVD stores around back then. The study was building on prior work that showed that people focus exclusively on explicitly presented details and ignore the facts that the explicit statements imply but don't state. (Also, I just looked it up and there's quite a few video stores around me (here in Portland) I guess I just haven't been to one in a while). The options that the subjects were given to circle depended on which group they belonged to. The control group (75 students) were given these options: A. Buy this entertaining video. B. Not buy this entertaining video. The treatment group (also 75 students) were given these options: A. Buy this entertaining video. B. Not buy this entertaining video. Keep the$14.99 for other purchases.

Here's the outcome of the study.

Outcome = Namespace(
control="Control Group",
treatment="Treatment Group",
total="Total"
)

data = pandas.DataFrame.from_dict({Outcome.buy:[56, 41],
data.index = [Outcome.control, Outcome.treatment]
row = data.sum()
row.name = Outcome.total
data = data.append(row)
data[Outcome.total] = data.sum(axis="columns")

print(TABLE(data))

Control Group 56 19 75
Treatment Group 41 34 75
Total 97 53 150
Total = Namespace(
group = data.loc[Outcome.control][Outcome.total],
students = data.loc[Outcome.total][Outcome.total],
)


Now a little plotting.

# put the groups into a column
plotter = data.reset_index()

# get rid of the totals
del(plotter["Total"])
plotter = plotter.iloc[:-1]

# move the outcome headers into a column
plotter = plotter.melt(id_vars=["index"],
plotter = plotter.rename(columns=dict(
index="Group",
variable="Outcome",
value="Count"))
plot = plotter.hvplot.bar(x="Group", y="Count", by="Outcome").opts(
width=Plot.width,
height=Plot.height,
color=Plot.color_cycle,
fontscale=Plot.fontscale,
)


print(outcome)


It looks like there was a significant difference since not only did the majority of the treatment group opt not to buy the DVD, but in the control a significant majority did.

Now as row-proportions.

proportions = data.divide(data[Outcome.total], axis="rows")
print(TABLE(proportions))

Control Group 0.746667 0.253333 1
Treatment Group 0.546667 0.453333 1
Total 0.646667 0.353333 1

Looking at the proportions it looks like quite a bit more didn't buy the DVD, but let's run the experiment and see.

### Point Estimate of Effect

The thing we are interested in is whether the wording of the second option to not buy the DVD made a difference, so our statistic of interest is the difference in proportions of the control and treatment group participants who didn't buy the DVD.

\begin{align} \hat{p}_{control} &= \frac{\textrm{Control Group that wouldn't buy}}{\textrm{size of Control Group}}\\ \hat{p}_{treatment} &= \frac{\textrm{Treatment Group that wouldn't buy}}{\textrm{size of Treament Group}}\\ \hat{p} &= \hat{p}_{treatment} - \hat{p}_{control} \end{align}
POINT_ESTIMATE = (proportions[Outcome.dont_buy].loc[Outcome.treatment]
print(f"{POINT_ESTIMATE * 100:.2f} %")

20.00 %


About twenty percent more of the treatment group said they would abstain from the purchase than control group.

### Simulating Random Chance

In a previous post looking at Gender Discrimination Inference I split the Urn into a 50-50 split of males and females to see if there was gender bias in choosing whether to promote them. In this case it's going to work a little differently.

We are asking here whether our treatment had an effect. If it didn't then we would expect that the distribution of "buy" and "don't buy" that we saw represents the distribution of the underlying population of ASU students, so we need to split our urn up to match the counts in the "Buy" and "Don't Buy" columns and then randomly split it into two equal groups. If the difference we saw was the result of random chance then we would expect that the difference between "buy" and "don't buy" group would be equal most of the time. This doesn't seem as intuitive a way to set it up as the previous method, but we'll see how it goes.

#### Set Up the Urn

buy = data.loc[Outcome.total][Outcome.buy]

didnt = numpy.ones(dont)

URN = numpy.append(bought, didnt)

expect(len(URN)).to(equal(STUDENTS))
expect(dont).to(equal(53))
expect(URN.sum()).to(equal(dont))


#### Simulate

The Process:

1. Shuffle the urn
2. Split the shuffled urn in half (control and treatment)
3. Count the proportion of each half that said they would not buy the DVD
4. Record the differences between the proportions of the control that wouldn't buy and the treatment that wouldn't buy

But, since I'm doing this with numpy it works a little different.

My process:

1. Randomly choose half the urn (without replacement) to get the control group
2. Sum the choices (this gives the count of the control that said they wouldn't buy)
3. Subtract the previous sum from the total that said they wouldn't buy (to get the treatment count)
4. Calculate the proportion of each group that said they wouldn't buy
5. Find the difference in proportion for each group
TRIALS = 1 * 10**6

with TIMER:
controls = numpy.array([numpy_random.choice(URN,
Total.group,
replace=False).sum()
for trial in range(TRIALS)])
treatments = Total.dont_buy - controls

control_proportions = controls/Total.group
treatment_proportions = treatments/Total.group

differences = control_proportions - treatment_proportions
simulation = pandas.DataFrame.from_dict({
"Point Estimate": differences,
"Control": control_proportions,
"Treatment": treatment_proportions,
})

2020-09-26 20:02:52,224 graeae.timers.timer start: Started: 2020-09-26 20:02:52.224515
2020-09-26 20:03:29,174 graeae.timers.timer end: Ended: 2020-09-26 20:03:29.174373
2020-09-26 20:03:29,175 graeae.timers.timer end: Elapsed: 0:00:36.949858


#### Plot the Distribution

plot = simulation.hvplot.hist("Point Estimate",
bins=25).opts(
title="Opportunity Cost Don't Buy Difference",
width=Plot.width,
height=Plot.height,
color=Plot.tan,
)
outcome = Embed(plot=plot, file_name="opportunity_cost_simulation")()

print(outcome)


The distribution appears to be reasonably normal, if a bit peaked in the middle.

matched = len(simulation[simulation["Point Estimate"] >= POINT_ESTIMATE])/len(simulation)
print(f"Percent of trials >= Point Estimate of Study: {100 * matched:0.3f} %")

Percent of trials >= Point Estimate of Study: 0.817 %


According to our simulation, less than one percent of the time we would see a difference like the study found by random chance alone.

#### Test Our Hypothesis

print(f"Null Hypothesis is {matched >= POINT_ESTIMATE}")

Null Hypothesis is False


So we'll conclude that the Null Hypothesis is false and that the Alternative Hypothesis that telling students about the opportunity cost of buying a DVD does have an effecte in getting them to not buy the DVD.

## End

So here we have a walk through of Inference using simulation and an experiment with Control and Treatment groups. Although the conclusion reached is that reminding students of the money they would have in the future if they didn't spend it is "causal", since this was an experiment, I'm not 100% convinced that asking studnents what the

### In the Abstract

1. Create an urn where the ones and zeros are equal to the totals for each of the outcome columns
2. Sample from the urn a simulated control group and treatment group (matching the sizes of the original groups)
3. Find the difference between the proportion of the control group and the treatment group that match the "success" outcome
4. Calculate the fraction of the differences that are equal to or greater than the differences in the study
5. Check if the fraction in the simulation meets your confidence level

# Gender Discrimination Inference

## Beginning

### Imports

# python
from argparse import Namespace
from functools import partial

# pypi
from expects import (
equal,
expect,
)
from numpy.random import default_rng
from tabulate import tabulate

import holoviews
import hvplot.pandas
import numpy
import pandas

# my stuff
from graeae import EmbedHoloviews, Timer


### Randomness

numpy_random = default_rng()


### Plotting, Tables, and a Timer

#### Plotting

Embed = partial(EmbedHoloviews,
folder_path="files/posts/gender-discrimination/")
Plot = Namespace(
width=990,
height=780,
fontscale=2,
tan="#ddb377",
blue="#4687b7",
red="#ce7b6d",
color_cycle = holoviews.Cycle(["#4687b7", "#ce7b6d"])
)


#### The Table

TABLE = partial(tabulate, tablefmt="orgtbl", headers="keys")


#### The Timer

TIMER = Timer()


## Middle

### The Original Study

This starts with a study where bank managers were each given a personnel file and asked to decide if they would promote the person represented by the file to branch manager. The files were identical but half of them were filled out as male and half as female. The researchers wanted to know if these manager were biased in favor of men. The table below shows their resuts.

study = pandas.DataFrame({"Promoted": [21, 14],
"Not Promoted":[3, 10]},
index=["Male", "Female"])
row = study.sum()
row.name = "Total"
study = study.append(row)
column = study.sum(axis="columns")
study["Total"] = column
print(TABLE(study))

Promoted Not Promoted Total
Male 21 3 24
Female 14 10 24
Total 35 13 48

Let's look at the values (without the totals) as a plot.

plotter = study.reset_index()
del(plotter["Total"])
plotter = plotter.iloc[:-1]
plotter = plotter.melt(id_vars=["index"],
value_vars=["Promoted", "Not Promoted"])
plotter = plotter.rename(columns=dict(index="Gender",
value="Count",
variable="Decision"))
plot = plotter.hvplot.bar(x="Gender", y="Count", by="Decision").opts(
title="Outcome by Gender",
height=Plot.height,
width=Plot.width,
fontscale=Plot.fontscale,
color=Plot.color_cycle,
)
outcome = Embed(plot=plot, file_name="promotions_by_gender")()

print(outcome)


It looks like a considerable majority of the males were promoted whereas for the females only a slight majority were promoted.

plot = plotter.hvplot.bar(x="Decision", y="Count", by="Gender").opts(
title="Male Vs Female by Decision",
height=Plot.height,
width=Plot.width,
fontscale=Plot.fontscale,
color=Plot.color_cycle,
)
outcome = Embed(plot=plot, file_name="decision_by_gender")()

print(outcome)


This plot doesn't make the disparity look quite as extreme as the previous plot did, I think.

#### About the Variables

I started using them without explicity stating it but we have two variable here - Promoted and Not Promoted are part of the decision variable and Male and Female are part of the gender variable.

Anyway… so more males were chosen for promotion than female, but by what proportion?

print(TABLE(study/study.loc["Total"]))

Promoted Not Promoted Total
Male 0.6 0.230769 0.5
Female 0.4 0.769231 0.5
Total 1 1 1

So 60% of the promoted were male and 40% of the promoted were female.

print(f"{.6/.4:.2f}")

1.50


If you were promoted you were one and a half times more likely to be male. Another way to look at it is to ask - What proportion of each gender was promoted?

fractions = study/study.Total.values
print(TABLE(fractions))

Promoted Not Promoted Total
Male 0.875 0.125 0.5
Female 0.583333 0.416667 0.5
Total 1.45833 0.541667 1

So about 88% of the males were promoted while only 58% of the females were promoted.

POINT_ESTIMATE = fractions.Promoted.loc['Male'] - fractions.Promoted.loc['Female']
print(f"{POINT_ESTIMATE:0.2f}")

0.29


There was a 30% difference between the rate of male promotion and the rate of female promotion. The question we have now is - could this difference have happened by random chance or is this evidence of bias?.

### The Experiment

We have a point estimate of the difference of 0.29 - is this evidence of bias?

#### Our Hypotheses

• $H_0$: Null Hypothesis. The variables gender and =decision are independent and the observed difference was due to chance.
• $H_A$: Alternative Hypothesis. The variables gender and decision are not independent and the difference was not due to chance, but rather women were less likely to be promoted than men.

I'm going to pick an arbitrary confidence interval of 0.95.

ALPHA = 0.05


#### The Simulation

The basic method here is we'll create an "urn" with an equal number of "balls" for men and women (24 of each in this case) and then randomly select 35 balls representing the number that were promoted and see the difference in the fraction of males promoted vs the number of females. To make the math simple I'll run it a number of times that is a power of 10.

• Some Setup

To make the counting easier I'll set males to be 1 and females to be 0 (so the number of males promoted is the sum and the females is the remainder).

NUMBER_OF_MALES = NUMBER_OF_FEMALES = 24
MALE = 1
FEMALE = 0
PROMOTED = 35
MALES = numpy.ones(NUMBER_OF_MALES)
FEMALES = numpy.zeros(NUMBER_OF_FEMALES)
URN = numpy.append(MALES, FEMALES)

expect(len(URN)).to(equal(NUMBER_OF_MALES + NUMBER_OF_FEMALES))
expect(URN.sum()).to(equal(NUMBER_OF_MALES))

• The Trials
TRIALS = 1 * 10**7
males_promoted = [numpy_random.choice(URN, PROMOTED, replace=False).sum()
for trial in range(TRIALS)]
females_promoted = [PROMOTED - males for males in males_promoted]
proportion_males_promoted = [males/NUMBER_OF_MALES for males in males_promoted]
proportion_females_promoted = [females/NUMBER_OF_FEMALES for females in females_promoted]
pairs = zip(proportion_males_promoted, proportion_females_promoted)
differences = [males - females for males, females in pairs]


### Looking at the Distribution

data = pandas.DataFrame.from_dict({
"Point Estimate": differences,
"Males Promoted": males_promoted,
"Females Promoted": females_promoted,
"Proportion Males Promoted": proportion_males_promoted,
"Proportion Females Promoted": proportion_females_promoted,
})

plot = data.hvplot.hist("Point Estimate").opts(
title="Distribution of Differences in Gender Promotion",
width=Plot.width,
height=Plot.height,
fontscale=Plot.fontscale,
)

outcome = Embed(plot=plot, file_name="difference_distribution")()

print(outcome)


The distribution looks mostly normal. As you might guess the distribution is centered around 0 (cases where exactly the same number of males and females were promoted - although 35 were promoted each time so the trials are never exactly 0). The cases where the difference in proportion is as great or greater than as it was in the original study are in the rightmost two bins.

I should note that because the sample size is so small, there's a weird distribution of the points - there's not enough variation in the differences to make a smooth curve (thus the gaps in the histogram).

But anyway, what proportion of our simulations had as much or more of a difference than that found in the study?

proportion = len(data[data["Point Estimate"] >= POINT_ESTIMATE])/TRIALS
print(f"Proportion: {100 * proportion:0.3f} %")

Proportion: 2.450 %

print(f"Probability of our study's difference in promotion between genders by chance alone: {proportion:0.2f}.")
print(f"Our tolerance was {ALPHA:0.2f}.")

Probability of our study's difference in promotion between genders by chance alone: 0.02.
Our tolerance was 0.05.


So we reject the null hypothesis and conclude that there is a statistically significant chance that the number of women promoted vs men in the original study was the result of bias.

## End

Well, that was the replication (sort of) of this problem from Introductory Statistics with Randomization and Simulation. The point of it was to both review Hypothesis Testing and see how it can be done using simulation rather than a p-test.

### In Abstract

1. Is your problem that you suspect some kind of bias in outcomes for two groups?
2. Get the Point Estimate for the value you want to test
3. State the Null Hypothesis that what happened could happen by chance and the Alternative Hypothesis that there was bias involved.
4. Decide on a tolerance level for the probability that it happened by chance.
5. Set up your urn
6. Simulate random selections for a large number of trials.
7. Calculate the proporiton of the trials that were greater than the original studies Point Estimate.
8. Make a conclusion whether the original outcome could have happened by random chance or not.

# Introductory Statistics with Randomization and Simulation

## Citation

• Diez DM, Barr CD, Çetinkaya-Rundel M. Introductory statistics with randomization and simulation. OpenIntro.org; 2014. 348 p.