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)

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.

• 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)

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))

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.

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

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

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.

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

• \(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"

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(
    buy="Buy DVD",
    dont_buy="Not Buy DVD",
    control="Control Group",
    treatment="Treatment Group",
    total="Total"
)

data = pandas.DataFrame.from_dict({Outcome.buy:[56, 41],
                                   Outcome.dont_buy: [19, 34]})
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))

Total = Namespace(
    buy = data.loc[Outcome.total][Outcome.buy],
    dont_buy = data.loc[Outcome.total][Outcome.dont_buy],
    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"],
                       value_vars=[Outcome.buy, Outcome.dont_buy])
plotter = plotter.rename(columns=dict(
    index="Group",
    variable="Outcome",
    value="Count"))
plot = plotter.hvplot.bar(x="Group", y="Count", by="Outcome").opts(
    title="Buy Or Don't Buy DVD",
    width=Plot.width,
    height=Plot.height,
    color=Plot.color_cycle,
    fontscale=Plot.fontscale,
)

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

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))

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

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.

POINT_ESTIMATE = (proportions[Outcome.dont_buy].loc[Outcome.treatment]
                  - proportions[Outcome.dont_buy].loc[Outcome.control])
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.

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.

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

• Introductory Statistics with Randomization and Simulation
• Opportunity Cost Neglect

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. There was actually two types of jobs, one described by the authors as "simple" and one as "complex". This only looks at the "simple" job. The table below shows their results.

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))

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.

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

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).

output = stats.probplot(data["Point Estimate"])
theoretical, actual = output[0]
slope, intercept, r = output[1]

first = theoretical.min()
last = theoretical.max()

line = holoviews.Curve([(first, slope * first + intercept), (last, slope * last + intercept)],
                       color=Plot.red)

qq_frame = pandas.DataFrame.from_dict({
    "Theoretical Quantiles": theoretical,
    "Differences In Hiring Rates": actual,
})

scatter = qq_frame.hvplot.scatter(x="Theoretical Quantiles",
                                  y="Differences In Hiring Rates", datashade=True)
plot = (scatter * line).opts(
    title="Hiring Probability Plot",
    width=Plot.width,
    height=Plot.height,
    fontscale=Plot.fontscale,
)
outcome = Embed(plot=plot, file_name="differences_qq")()

print(outcome)

Hmm… I'm not sure I'm doing that right. But it looks reasonable, I suppose.

lowest = data["Point Estimate"].min()
data_sorted = data["Point Estimate"].sort_values()

max_loc = int(len(data_sorted) * 0.95)
highest = data_sorted.iloc[max_loc]

print(f"Point Estimate: {POINT_ESTIMATE: 0.3f}")
print(f"({lowest: 0.3f} , {highest: 0.3f})")

Point Estimate:  0.292
(-0.542 ,  0.208)

Looking at the confidence interval we can see that the original point estimate is higher than 95 % of the values, but the original book example used a hypothesis test so let's finish it up.

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.

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
• Rosen B, Jerdee TH. Influence of sex role stereotypes on personnel decisions. Journal of Applied Psychology. 1974 Feb;59(1):9 Link to PDF

BASE = "COVID-19-Weekly-Report-{date}-FINAL.pdf"

DATE = "2020-05-27"
load_dotenv(Path("~/.local/share/env").expanduser())
path = Path(os.environ["OREGON_COVID"]).expanduser()
filename = path/BASE.format(date=DATE)
assert filename.is_file()
tables = tabula.read_pdf(filename, pages="all")

print(type(tables))
print(len(tables))

<class 'list'>
15

Okay, so it found 15 tables, which should have been loaded as pandas DataFrames. Since a list doesn't really provide an easy way to get metadata, maybe I'll just traverse the list and see what's there.

for table in tables:
    print(table.head(1))
    print("*" * 20)

    Sex  Cases % of total cases Cases per 10,000a  Deaths Case fatality (%)  \
0  Male  1,874            47.5%               8.7      85              4.5%   

   Hospitalized  
0           402  
********************
  Age group Cases % of total cases  Cases per 10,000a  Deaths  \
0       0-9    53             1.3%                1.1       0   

  Case fatality (%)  Hospitalized  
0                0%             6  
********************
  Unnamed: 0  Race  Cases % of total cases Cases per 10,000b  Deaths  \
0      White   NaN  1,946            49.3%               6.2   109.0   

  Case fatality (%)  Hospitalized  
0              5.6%         452.0  
********************
  Ethnicity Case count % of total cases Cases per 10,000a  Deaths  \
0  Hispanic      1,317            33.4%              23.7      14   

  Case fatality (%)  Hospitalized  
0              1.1%           178  
********************
                Facility name     County First reported  Total casesa  \
0  Healthcare at Foster Creek  Multnomah      3/24/2020         119.0   

   Total deathsa  
0           30.0  
********************
   ZIP code Number of cases Cases per 10,000
0     97002             1–9              n/a
********************
   97053 1–9  n/a
0  97054   0  0.0
********************
   97148   0  0.0
0  97201  13  8.4
********************
   97324   0   0.0
0  97325  10  15.1
********************
   97411    0  0.0
0  97415  1–9  n/a
********************
   97501 1–9  n/a
0  97502  13  4.7
********************
   97813  1–9  n/a
0  97814  1–9  n/a
********************
  Hospital                    Max Patients    Max Patients.1
0      NaN  (Both suspected and confirmed)  (Confirmed only)
********************
                               Hospital
0  Asante Rogue Regional Medical Center
********************
  Providence Newberg Medical Center
0       Providence Seaside Hospital
********************

Most of the tables appear to be either fragments or case-counts based on different categorizations.

COLUMNS = Namespace(
    sex="Sex",
    cases="Cases",
    cases_per_10000 = "Cases per 10,000a",
    percent = "% of total cases",
)

tables_2 = dict(sex=tables[0])

print(tables_2["sex"])

             Sex  Cases % of total cases Cases per 10,000a  Deaths  \
0           Male  1,874            47.5%               8.7      85   
1         Female  2,071            52.4%               9.9      63   
2     Non-Binary      1               0%               n/a       0   
3  Not available      3             0.1%               n/a       0   
4          Total  3,949           100.0%               9.3     148   

  Case fatality (%)  Hospitalized  
0              4.5%           402  
1              3.0%           344  
2                0%             0  
3                0%             1  
4              3.7%           747  

One thing to notice is that there are "n/a", ",", and "%" strings in there that prevent some columns from being treated as numbers.

tables_2["sex"].info()
males_females = tables[0].iloc[:2]

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 5 entries, 0 to 4
Data columns (total 7 columns):
 #   Column             Non-Null Count  Dtype 
---  ------             --------------  ----- 
 0   Sex                5 non-null      object
 1   Cases              5 non-null      object
 2   % of total cases   5 non-null      object
 3   Cases per 10,000a  5 non-null      object
 4   Deaths             5 non-null      int64 
 5   Case fatality (%)  5 non-null      object
 6   Hospitalized       5 non-null      int64 
dtypes: int64(2), object(5)
memory usage: 408.0+ bytes

So to use most of the columns you have to do something to convert them back to numbers. I'll start by looking at male and female cases adjusted by population.

def bar_plot(data: pandas.DataFrame, x_column: str, y_column: str,
             title:str, filename: str) -> str:
    """Plot the data as a bar-plot

    Args:
     data: thing to plot
     x_column: name of column for the x-axis
     y_column: name of the column with the bar height
     title: Title for the plot
     filename: Name to save the plot

    Returns:
     html: the html to use to embed the plot into org
    """

    plot = data.hvplot.bar(x=x_column,
                           y=y_column,
                           color=Plot.color_cycle).opts(
                               width=Plot.width,
                               height=Plot.height,
                               fontscale=Plot.font_scale,
                               title=title,
)

    outcome = Embed(plot=plot, file_name=filename)()
    return outcome

gender_bar = partial(bar_plot, data=males_females, x_column=COLUMNS.sex)