Combine data to investigate a hypothesis

This stuy case is based on the research question whether there is relation between town temperature and health risks. We would like to compare towns that have the most highest temperature with other towns in the Unitied States. The hypothesis is that the health risk value is higher for towns with higher temperature compared to ordinary towns. In particular the risk related to 'Adult Binch Drinking'. For this we need two datasets. The dataset of hottest towns and the dataset of USHealth_data.csv. The end goal is to run statistical comparion test. Before we can do such we first need to clean, organize and merge the data. The column 'Value' for instance contains not all float type data.

from bokeh.io import output_notebook
output_notebook()

Loading BokehJS ...

Load first dataset: hotest towns

First we load the text file hottest_towns.txt. This data we put in a dataframe. If we look with an editor in the file we see that the first two lines are not of interest, so we skip these. Furthermore there are no column labels so we add these.

import pandas as pd
df=pd.read_table('data/hottest_towns.txt', skiprows = 2, names=['place', 'days'])
df.head()

	place	days
0	Phoenix, Arizona	107
1	Las Vegas, Nevada	70
2	Riverside, California	24
3	Dallas, Texas	17
4	Austin, Texas	16

As we can see we need to split the City from the State. Let's do this

df[['City', 'State']] = df['place'].str.split(',', n=1, expand=True)
ht = df.drop('place', axis = 1) #drop column place
ht.head()

	days	City	State
0	107	Phoenix	Arizona
1	70	Las Vegas	Nevada
2	24	Riverside	California
3	17	Dallas	Texas
4	16	Austin	Texas

Load the health data

First we download data from http://www.bigcitieshealth.org/obesity-physical-activity and inspect it with .head()

df = pd.read_csv('data/BCHI-dataset_2019-03-04.csv', low_memory=False)
df.head(2)

	Indicator Category	Indicator	Year	Sex	Race/Ethnicity	Value	Place	BCHC Requested Methodology	Source	Methods	Notes	90% Confidence Level - Low	90% Confidence Level - High	95% Confidence Level - Low	95% Confidence Level - High
0	Behavioral Health/Substance Abuse	Opioid-Related Unintentional Drug Overdose Mor...	2010	Both	All	1.7	Washington, DC	Age-Adjusted rate of opioid-related mortality ...	D.C. Department of Health, Center for Policy, ...	NaN	This indicator is not exclusive of other drugs...	NaN	NaN	NaN	NaN
1	Behavioral Health/Substance Abuse	Opioid-Related Unintentional Drug Overdose Mor...	2010	Both	All	2.2	Fort Worth (Tarrant County), TX	Age-adjusted rate of opioid-related mortality ...	National Center for Health Statistics	NaN	This indicator is not exclusive of other drugs...	NaN	NaN	1.5	3.0

It seems that we need to do some cleaning as well in the place column. States are abbreviated. Let us split the column like we did in the hottest town data set and than change the state with the map function.

df[['City', 'State']] = df['Place'].str.split(',', n=1, expand=True)

df.head(2)

	Indicator Category	Indicator	Year	Sex	Race/Ethnicity	Value	Place	BCHC Requested Methodology	Source	Methods	Notes	90% Confidence Level - Low	90% Confidence Level - High	95% Confidence Level - Low	95% Confidence Level - High	City	State
0	Behavioral Health/Substance Abuse	Opioid-Related Unintentional Drug Overdose Mor...	2010	Both	All	1.7	Washington, DC	Age-Adjusted rate of opioid-related mortality ...	D.C. Department of Health, Center for Policy, ...	NaN	This indicator is not exclusive of other drugs...	NaN	NaN	NaN	NaN	Washington	DC
1	Behavioral Health/Substance Abuse	Opioid-Related Unintentional Drug Overdose Mor...	2010	Both	All	2.2	Fort Worth (Tarrant County), TX	Age-adjusted rate of opioid-related mortality ...	National Center for Health Statistics	NaN	This indicator is not exclusive of other drugs...	NaN	NaN	1.5	3.0	Fort Worth (Tarrant County)	TX

states = {'OH': 'Ohio', 'KY': 'Kentucky', 'AS': 'American Samoa', 'NV': 'Nevada', 
        'WY': 'Wyoming', 'NA': 'National', 'AL': 'Alabama', 'MD': 'Maryland', 
        'AK': 'Alaska', 'UT': 'Utah', 'OR': 'Oregon', 'MT': 'Montana', 'IL': 'Illinois', 
        'TN': 'Tennessee', 'DC': 'District of Columbia', 'VT': 'Vermont', 'ID': 'Idaho', 
        'AR': 'Arkansas', 'ME': 'Maine', 'WA': 'Washington', 'HI': 'Hawaii', 'WI': 'Wisconsin', 
        'MI': 'Michigan', 'IN': 'Indiana', 'NJ': 'New Jersey', 'AZ': 'Arizona', 'GU': 'Guam', 
        'MS': 'Mississippi', 'PR': 'Puerto Rico', 'NC': 'North Carolina', 'TX': 'Texas', 
        'SD': 'South Dakota', 'MP': 'Northern Mariana Islands', 'IA': 'Iowa', 'MO': 'Missouri', 
        'CT': 'Connecticut', 'WV': 'West Virginia', 'SC': 'South Carolina', 'LA': 'Louisiana', 
        'KS': 'Kansas', 'NY': 'New York', 'NE': 'Nebraska', 'OK': 'Oklahoma', 'FL': 'Florida', 
        'CA': 'California', 'CO': 'Colorado', 'PA': 'Pennsylvania', 'DE': 'Delaware', 
        'NM': 'New Mexico', 'RI': 'Rhode Island', 'MN': 'Minnesota', 'VI': 'Virgin Islands',
        'NH': 'New Hampshire', 'MA': 'Massachusetts', 'GA': 'Georgia', 'ND': 'North Dakota', 'VA': 'Virginia'}

df['State'] = df['State'].str.strip().map(states)

df.head(2)

	Indicator Category	Indicator	Year	Sex	Race/Ethnicity	Value	Place	BCHC Requested Methodology	Source	Methods	Notes	90% Confidence Level - Low	90% Confidence Level - High	95% Confidence Level - Low	95% Confidence Level - High	City	State
0	Behavioral Health/Substance Abuse	Opioid-Related Unintentional Drug Overdose Mor...	2010	Both	All	1.7	Washington, DC	Age-Adjusted rate of opioid-related mortality ...	D.C. Department of Health, Center for Policy, ...	NaN	This indicator is not exclusive of other drugs...	NaN	NaN	NaN	NaN	Washington	District of Columbia
1	Behavioral Health/Substance Abuse	Opioid-Related Unintentional Drug Overdose Mor...	2010	Both	All	2.2	Fort Worth (Tarrant County), TX	Age-adjusted rate of opioid-related mortality ...	National Center for Health Statistics	NaN	This indicator is not exclusive of other drugs...	NaN	NaN	1.5	3.0	Fort Worth (Tarrant County)	Texas

Inspect data

There are quit some missing values. We might want to drop the NaN's of the Values since this is our column of interest

df.isnull().sum()

Indicator Category                 0
Indicator                          0
Year                               0
Sex                                0
Race/Ethnicity                     0
Value                           3444
Place                              0
BCHC Requested Methodology         0
Source                          4964
Methods                        26968
Notes                          19851
90% Confidence Level - Low     31655
90% Confidence Level - High    31655
95% Confidence Level - Low     25657
95% Confidence Level - High    25616
City                               0
State                           2154
dtype: int64

df = df.dropna(subset = ['Value'])
df.isnull().sum()

Indicator Category                 0
Indicator                          0
Year                               0
Sex                                0
Race/Ethnicity                     0
Value                              0
Place                              0
BCHC Requested Methodology         0
Source                          4133
Methods                        24226
Notes                          19850
90% Confidence Level - Low     28214
90% Confidence Level - High    28214
95% Confidence Level - Low     22304
95% Confidence Level - High    22263
City                               0
State                           2154
dtype: int64

Let us inspect the indicator values. We were interested in adult binch drinking

set(df['Indicator'])

{'AIDS Diagnoses Rate (Per 100,000 people)',
 'All Types of Cancer Mortality Rate (Age-Adjusted; Per 100,000 people)',
 'All-Cause Mortality Rate (Age-Adjusted; Per 100,000 people)',
 'Asthma Emergency Department Visit Rate (Age-Adjusted; Per 10,000)',
 'Bike Score',
 'Chlamydia Rate (Per 100,000 People)',
 'Congenital Syphilis Rate (Per 100,000 Live Births)',
 'Diabetes Mortality Rate (Age-Adjusted; Per 100,000 people)',
 'Female Breast Cancer Mortality Rate (Age-Adjusted; Per 100,000 people)',
 'Firearm-Related Emergency Department Visit Rate (Age-Adjusted; Per 10,000 people)',
 'Firearm-Related Mortality Rate (Age-Adjusted; Per 100,000 people)',
 'Gonorrhea Rate (Per 100,000 People)',
 'HIV Diagnoses Rate (Per 100,000 people)',
 'HIV-Related Mortality Rate (Age-Adjusted; Per 100,000 people)',
 'Heart Disease Mortality Rate (Age-Adjusted; Per 100,000 people)',
 'Homicide Rate (Age-Adjusted; Per 100,000 people)',
 'Infant Mortality Rate (Per 1,000 live births)',
 'Life Expectancy at Birth (Years)',
 'Lung Cancer Mortality Rate (Age-Adjusted; Per 100,000 people)',
 'Median Household Income (Dollars)',
 'Motor Vehicle Mortality Rate (Age-Adjusted; Per 100,000 people)',
 'Opioid-Related Unintentional Drug Overdose Mortality Rate (Age-Adjusted; Per 100,000 people)',
 'Percent Foreign Born',
 'Percent Living Below 200% Poverty Level',
 'Percent Unemployed',
 'Percent Who Only Speak English at Home',
 'Percent Who Speak Spanish at Home',
 'Percent of 3 and 4 Year Olds Currently Enrolled in Preschool',
 'Percent of Adults 65 and Over Who Received Pneumonia Vaccine',
 'Percent of Adults Who Are Obese',
 'Percent of Adults Who Binge Drank',
 'Percent of Adults Who Currently Smoke',
 'Percent of Adults Who Meet CDC-Recommended Physical Activity Levels',
 'Percent of Adults Who Received Seasonal Flu Shot',
 'Percent of Children (Tested) Under Age 6 with Elevated Blood Lead Levels',
 'Percent of Children Living in Poverty',
 'Percent of Children Who Received Seasonal Flu Shot',
 'Percent of High School Graduates (Over Age 18)',
 'Percent of High School Students Who Are Obese',
 'Percent of High School Students Who Binge Drank',
 'Percent of High School Students Who Currently Smoke',
 'Percent of High School Students Who Meet CDC-Recommended Physical Activity Levels',
 'Percent of Households Whose Housing Costs Exceed 35% of Income',
 'Percent of Low Birth Weight Babies Born',
 'Percent of Mothers Under Age 20',
 'Percent of Population 65 and Over',
 'Percent of Population Under 18',
 'Percent of Population Uninsured',
 'Percent of Population with a Disability',
 'Persons Living with HIV/AIDS Rate (Per 100,000 people)',
 'Pneumonia and Influenza Mortality Rate (Age-Adjusted; Per 100,000 people)',
 'Primary and Secondary Syphilis Rate (Per 100,000 People)',
 'Race/Ethnicity (Percent)',
 'Rate of Laboratory Confirmed Infections Caused by Salmonella (Per 100,000 people)',
 'Rate of Laboratory Confirmed Infections Caused by Shiga Toxin-Producing E-Coli (Per 100,000 people)',
 'Sex (Percent)',
 'Suicide Rate (Age-Adjusted; Per 100,000 people)',
 'Total Population (People)',
 'Transit Score',
 'Tuberculosis Incidence Rate (Per 100,000 people)',
 'Walkability'}

df = df[df['Indicator'] == 'Percent of Adults Who Binge Drank']
len(df)

535

The question is wether the names of the cities and the states are comparable

print(sorted(set(df['City'])))
print("\n")
print(sorted(set(ht['City'])))

['Baltimore', 'Boston', 'Charlotte', 'Chicago', 'Columbus', 'Denver', 'Detroit', 'Fort Worth (Tarrant County)', 'Houston', 'Indianapolis (Marion County)', 'Kansas City', 'Las Vegas (Clark County)', 'Long Beach', 'Los Angeles', 'Miami (Miami-Dade County)', 'Minneapolis', 'New York City', 'Oakland (Alameda County)', 'Philadelphia', 'Phoenix', 'Portland (Multnomah County)', 'San Antonio', 'San Diego County', 'San Jose', 'Seattle', 'U.S. Total', 'Washington']


['Austin', 'Dallas', 'Houston', 'Kansas City', 'Las Vegas', 'Oklahoma City', 'Phoenix', 'Riverside', 'Sacramento', 'Salt Lake City', 'San Antonio']

df = df.dropna(subset = ['State'])
print((sorted(set(df['State']))))
print("\n")
print((sorted(set(ht['State']))))

['Arizona', 'California', 'Colorado', 'District of Columbia', 'Florida', 'Indiana', 'Maryland', 'Massachusetts', 'Michigan', 'Minnesota', 'Missouri', 'Nevada', 'New York', 'North Carolina', 'Ohio', 'Oregon', 'Pennsylvania', 'Texas', 'Washington']


[' Arizona ', ' California ', ' Missouri ', ' Nevada ', ' Oklahoma ', ' Texas ', ' Utah ']

As we can see we have to get rid of the (additional text) at some names. We can do such using regular expressions. The regular expression \s\([a-zA-Z\s\-]+\) matches any name starting with opening bracket (, followed by a name including spaces or '-' followed by a closing bracket )

df['City']= df['City'].str.replace(r"\s\([a-zA-Z\s\-]+\)", "")
print((sorted(set(df['City']))))

['Baltimore', 'Boston', 'Charlotte', 'Columbus', 'Denver', 'Detroit', 'Fort Worth', 'Houston', 'Indianapolis', 'Kansas City', 'Las Vegas', 'Long Beach', 'Los Angeles', 'Miami', 'Minneapolis', 'New York City', 'Oakland', 'Philadelphia', 'Phoenix', 'Portland', 'San Antonio', 'San Diego County', 'San Jose', 'Seattle', 'Washington']

Lastly we clean the States with strip() to make sure we have matching names

df['State'] = df['State'].str.strip()
ht['State'] = ht['State'].str.strip()

Now we cleaned and matched the city and states we can use this to merge the data

Merge both datasets to create a subset hottest towns health data

We join the datasets on City and State. Drop all the rows from that have NaN in the 'Value'

hh = pd.merge(df, ht, how='inner', left_on=['City', 'State'], right_on=['City', 'State'])
hh.head(2)

	Indicator Category	Indicator	Year	Sex	Race/Ethnicity	Value	Place	BCHC Requested Methodology	Source	Methods	Notes	90% Confidence Level - Low	90% Confidence Level - High	95% Confidence Level - Low	95% Confidence Level - High	City	State	days
0	Behavioral Health/Substance Abuse	Percent of Adults Who Binge Drank	2010	Both	Black	24.3	San Antonio, TX	BRFSS (or similar) How many times during the p...	2013 BRFSS Bexar County	NaN	Data includes Bexar County, TX, not just San A...	NaN	NaN	NaN	NaN	San Antonio	Texas	8
1	Behavioral Health/Substance Abuse	Percent of Adults Who Binge Drank	2010	Both	Hispanic	25.9	San Antonio, TX	BRFSS (or similar) How many times during the p...	2013 BRFSS Bexar County	NaN	Data includes Bexar County, TX, not just San A...	NaN	NaN	NaN	NaN	San Antonio	Texas	8

Create a dataset without the hottest towns

Now we need to create the second dataset without the hottest towns. We can use a trick to do such. We can create a list of towns and select all that is not in that list.

x = ht['City'].tolist()
nh =  df.where(~df['City'].isin(x))
nh.head(2)

	Indicator Category	Indicator	Year	Sex	Race/Ethnicity	Value	Place	BCHC Requested Methodology	Source	Methods	Notes	90% Confidence Level - Low	90% Confidence Level - High	95% Confidence Level - Low	95% Confidence Level - High	City	State
568	Behavioral Health/Substance Abuse	Percent of Adults Who Binge Drank	2010.0	Both	All	10.9	Miami (Miami-Dade County), FL	BRFSS (or similar) How many times during the p...	BRFSS	NaN	NaN	NaN	NaN	NaN	NaN	Miami	Florida
569	Behavioral Health/Substance Abuse	Percent of Adults Who Binge Drank	2010.0	Both	All	13.3	Charlotte, NC	BRFSS (or similar) How many times during the p...	2010 NC BRFSS (Mecklenburg Sample)	NaN	NaN	NaN	NaN	8.7	17.9	Charlotte	North Carolina

y_hot = nh['Value'].dropna()
y_nhot = hh['Value'].dropna()

y_hot.describe()

count    392.000000
mean      19.822194
std        7.395218
min        1.700000
25%       14.600000
50%       18.700000
75%       24.000000
max       54.000000
Name: Value, dtype: float64

y_nhot.describe()

count    103.000000
mean      18.218447
std       11.350605
min        6.300000
25%       12.800000
50%       15.300000
75%       21.150000
max       87.500000
Name: Value, dtype: float64

Visual inspection

A histogram can give insight in the distribution and whether both categories are comparable. From the following plot we can see that both distributions are comparable, but hottest town is a bit to the right if you ignore the outlier of not hottest towns

import numpy as np

from bokeh.layouts import gridplot
from bokeh.plotting import figure, show

def make_plot(data, color, label):
    measured = data
    hist, edges = np.histogram(measured, density=True, bins=50)
    p.quad(top=hist, bottom=0, left=edges[:-1], right=edges[1:],
           fill_color=color, line_color="white", alpha=0.5, legend_label=label)

    p.y_range.start = 0
    p.legend.location = "center_right"
    p.legend.background_fill_color = "#fefefe"
    p.xaxis.axis_label = 'Percent of Adults Who Binge Drank'
    p.yaxis.axis_label = 'Probability'
    p.grid.grid_line_color="white"
    return p

# Normal Distribution
p = figure(title="Distribution", tools='', background_fill_color="#fafafa")
p = make_plot(y_hot, 'red', 'hottest towns')
p = make_plot(y_nhot, 'green', 'normal towns')

show(p)

We still see a few outliers. We might want to delete them. Of course you first need to investigate wether you can disgard them. For demonstration purpose we drop the outliers.

y_nhot = y_nhot.where(y_nhot < 80).dropna()

Statistical tests

If we assume normality we can run the t-test to compare significant differences. If we do not we should use a non parametic test. From the plot we can see that our data is not normal.

#normal
from scipy.stats import ttest_ind
ttest_ind(np.array(y_hot), np.array(y_nhot), equal_var = False)
#not normal
from scipy.stats import kruskal
kruskal(np.array(y_hot), np.array(y_nhot))

Krukruskal(np.array(y_hot), np.array(y_nhot))

KruskalResult(statistic=14.72621798981739, pvalue=0.0001243056464991998)

Conclusion

The hotest town Value of Binch Drinking is significant higher than the regular towns.