import pandas as pd
import geopandas as gpd
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.api as sm
from shapely.geometry import Point
from shapely import wkt
import geopandas
try:
from pymc3 import *
import pymc3 as pm
except:
! pip install pymc3
from pymc3 import *
import pymc3 as pm
import arviz as az
from pymc3 import glm
from scipy.stats import poisson, norm, gamma
%matplotlib inline
## Reading into UNFILTERED Chronic Disease Indicator dataset.
us_cdi = pd.read_csv('us_cdi.csv')
us_cdi.info()
us_cdi.head()
us_cdi.describe()
#all values in dropped col are null -> drop
us_cdi[us_cdi['Response'].notnull()]
us_cdi = us_cdi.drop(['Response','StratificationCategory2','Stratification2','StratificationCategory3','Stratification3','StratificationCategoryID2','StratificationID2','StratificationCategoryID3','StratificationID3','ResponseID'], axis=1)
us_cdi.head()
us_cdi['Question'].value_counts()
us_cdi['DataSource'].value_counts()
us_cdi['Stratification1'].value_counts()
us_cdi['Topic'].value_counts()
us_cdi['DataValueType'].value_counts()
Figure A
ax = us_cdi['Topic'].value_counts().plot(kind='barh')
ax.set_ylabel('Disease Factors')
ax.set_xlabel('Count')
ax.set_title('Number of US Individuals By Disease-Related Factors / Topic')
df = us_cdi.pivot_table(values='DataValueAlt',index=['Topic','QuestionID','Question','DataValueUnit','DataValueType'],
columns='LocationAbbr',aggfunc='mean',dropna=True).round(2)
df
Figure B
The figure below shows the rates of the diseases grouped by gender over the course of 2008-2019. We want to explore the prominence of specific diseases such as asthma and diabetes, and the difference in prevalence among the two genders.
gen = us_cdi[us_cdi['StratificationCategory1'] == 'Gender']
gen1 = gen.groupby(['Stratification1','YearStart'])
gen1.mean().drop('LocationID',axis=1).round()
plt.figure(figsize=(16, 6))
sns.barplot(x='YearStart',y='DataValueAlt',data=gen,hue='Stratification1',ci=None,saturation=0.7)
gen = us_cdi[(us_cdi['StratificationCategory1'] == 'Gender')]# & (us_cdi['Topic'] =='Diabetes')]
gen2 = gen.groupby(['Stratification1','Topic'])
gen2.mean().drop('LocationID',axis=1).round()
plt.figure(figsize=(35, 10))
ax = sns.barplot(x='Topic',y='DataValueAlt',data=gen,hue='Stratification1',ci=None,saturation=0.7)
Figure C
The figure shows the prevalence of diseases according to gender. This allows us to see which diseases we can potentially do further research on and narrow our scope by not only the diseases, but as well as the questions associated with this disease.
gen = us_cdi[(us_cdi['StratificationCategory1'] == 'Gender')]# & (us_cdi['Topic'] =='Diabetes')]
gen2 = gen.groupby(['Stratification1','Topic'])
gen2.mean().drop('LocationID',axis=1).round()
plt.figure(figsize=(35, 10))
ax = sns.barplot(x='Topic',y='DataValueAlt',data=gen,hue='Stratification1',ci=None,saturation=0.7)
ax.set_xticklabels(ax.get_xticklabels(), rotation = 45, ha = 'right')
plt.tight_layout()
plt.show()
gen = us_cdi[(us_cdi['StratificationCategory1'] == 'Overall')]
gen2 = gen.groupby(['Stratification1','Topic'])
gen2.mean().drop('LocationID',axis=1).round()
plt.figure(figsize=(35, 10))
ax = sns.barplot(x='Topic',y='DataValueAlt',data=gen,hue='Stratification1',ci=None,saturation=0.7)
ax.set_xticklabels(ax.get_xticklabels(), rotation = 45, ha = 'right')
plt.tight_layout()
plt.show()
Figure D
The figure shows the prevalence of diseases according to race/ethnicity. This allows us to see which diseases we can potentially do further research on and narrow our scope by not only the diseases, but as well as the questions associated with this disease. We can further determine which diseases are more prevalent within a race/ethnicity group we may focus on.
race = us_cdi[us_cdi['StratificationCategory1'] == 'Race/Ethnicity']
race1 = gen.groupby(['Stratification1','YearStart'])
race1.mean().drop('LocationID',axis=1).round(0)
plt.figure(figsize=(16, 6))
sns.barplot(x='YearStart',y='DataValueAlt',data=race,hue='Stratification1',ci=None,saturation=0.7)
a = us_cdi[us_cdi[('StratificationCategory1'] == 'Race/Ethnicity') & us_cdi[('StratificationCategory1'] == 'Race/Ethnicity')]
a = a.loc[(a['Question'] == 'Mortality from total cardiovascular diseases') & (a['Topic'] == 'Cardiovascular Disease')]
a = a.loc[(a['DataValueType'] == 'Crude Rate')]
a1 = a.groupby(['Question','YearStart'])
#a.mean().drop('LocationID',axis=1).round(0)
plt.figure(figsize=(16, 6))
sns.barplot(x='YearStart',y='DataValue',data=a,hue='Question',ci=None,saturation=0.7)
temp = us_cdi['GeoLocation']
temp = temp.str.slice(start=8).str.strip(')')
temp = temp.str.split("//s+")
#temp = temp.astype(str)
#temp = temp.astype(int)
#temp = temp.str.strip(')').str.strip(' (').str.strip('POINT')
# temp['lon'] = temp.GeoLocation.x
# temp['lat'] = temp.point_object.y
temp
temp = us_cdi['GeoLocation']
temp = temp.str.slice(start=8).str.strip(')')
temp = temp.str.split("//s+")
temp
##
#import matplotlib.pyplot as plt
#us_cdi.plot(kind="scatter", x="longitude", y="latitude", alpha=0.4)
#plt.show()
#basemap
fig, ax = plt.subplots(figsize=(50, 10))
states = geopandas.read_file('cb_2018_us_state_500k.shp')
type(states)
states.boundary.plot(ax = ax, color="Black")
geo_df = gpd.GeoDataFrame(us_cdi, geometry='GeoLocation')
geo_df.head()
#us_cdi.plot(ax = ax, markersize=20, color='red')
Research Question Results
– Summarize and interpret your models. – Estimate any uncertainty in your GLM predictions, providing clear quantitative state- ments of the uncertainty in plain English.
Checkpoint
df = us_cdi.loc[us_cdi['Topic'] == 'Cardiovascular Disease']
df.head()
df['Question'].value_counts()
df = df.loc[(df['Question'] == 'Mortality from total cardiovascular diseases') & (df['Topic'] == 'Cardiovascular Disease')]
df = df.loc[(df['DataValueType'] == 'Age-adjusted Rate')]
df = df.loc[(df['StratificationCategory1'] == 'Race/Ethnicity')]
df
#drop response, datavalue symbol,
df = df[df['DataValue'].notna()]
df = df.drop(df.columns[[19,20,21,22,23]], axis=1)
df = df.drop(['DataValueFootnoteSymbol', 'DatavalueFootnote'], axis=1)
#df = df.drop(['Response', 'DataValueFootnoteSymbol', 'DatavalueFootnote', 'Stratification2', 'StratificationCategory2', 'StratificationCategory3', 'Stratification3', 'ResponseID'], axis=1)
#df = df.drop(df.columns[[17, 18,19,20,21,22,23,24,25]], axis=1)
df
d.columns.values
np.asarray(d.columns)
Comparing GLM and non-parametric methods: DT/RF/OSL
d2 = pd.get_dummies(df['LocationAbbr'], prefix='LocationAbbr')
df3 = pd.concat([df2, d2], axis=1)
df3
#d2.columns.values
import numpy as np
np.mean(df3['LowConfidenceLimit'] - df3['HighConfidenceLimit'])
from sklearn.model_selection import train_test_split
train, test = train_test_split(df3,test_size=.3, random_state = 4)
input_1
np.asarray(df3.YearEnd)
#Frequentist GLM - Poisson (race + year)
y = (df3.DataValue).astype(None)
poisson_model_raceyear = sm.GLM(
y, sm.add_constant(df3[['Stratification_American Indian or Alaska Native',
'Stratification_Asian or Pacific Islander',
'Stratification_Black, non-Hispanic', 'Stratification_Hispanic',
'Stratification_White, non-Hispanic', 'YearStart']]),
family=sm.families.Poisson()
)
poisson_results_raceyear = poisson_model_raceyear.fit()
print(poisson_results_raceyear.summary())
#Explanation: standard errors is too low for YearStart, which means that we have a problem with overconfidence. For our race variables,
#this issue is resolved and standard erros (while still low) aren't as unreaosnable. Overall, poisson does not seem to be GREAT
# when used in combination fo race and year. We want to try this with negative binomial and also with other features to differentiate.
#Frequentist GLM - Poisson (race ONLY)
y = (df3.DataValue).astype(None)
poisson_model_race = sm.GLM(
y, sm.add_constant(df3[['Stratification_American Indian or Alaska Native',
'Stratification_Asian or Pacific Islander',
'Stratification_Black, non-Hispanic', 'Stratification_Hispanic',
'Stratification_White, non-Hispanic']]),
family=sm.families.Poisson()
)
poisson_results_race = poisson_model_race.fit()
print(poisson_results_race.summary())
#Again: same problem, low standard of error, model is too confident. Can conclude poisson distribution is not a good fit for this data
# We will try this again with location once we use medicaid expenditure, etc. as features
#Frequent GLM - Negative Binomial (race + year)
y = (df3.DataValue).astype(None)
negbin_model_raceyear = sm.GLM(
y, sm.add_constant(df3[['Stratification_American Indian or Alaska Native',
'Stratification_Asian or Pacific Islander',
'Stratification_Black, non-Hispanic', 'Stratification_Hispanic',
'Stratification_White, non-Hispanic', 'YearStart']]),
family=sm.families.NegativeBinomial()
)
negbin_results_raceyear = negbin_model_raceyear.fit()
print(negbin_results_raceyear.summary())
#Frequentist Model - Negative Binomial (race only)
y = (df3.DataValue).astype(None)
negbin_model_race = sm.GLM(
y, sm.add_constant(df3[['Stratification_American Indian or Alaska Native',
'Stratification_Asian or Pacific Islander',
'Stratification_Black, non-Hispanic', 'Stratification_Hispanic',
'Stratification_White, non-Hispanic']]),
family=sm.families.NegativeBinomial()
)
negbin_results_race = negbin_model_race.fit()
print(negbin_results_race.summary())
#All of these worked terribly: log-likelihood was very negative (good model has larger log-likelihood). Chi-square statistic should
# roughly be n - p = 1885 - 5 = 1880 but instead is 123 LOL).
#Conclusion: our data does not fit a poisson or negative binomial distribution very well. We can condlude GLMs are not good for this data
# and nonparametric works a lot better
df4 = df3.copy()
df4 = df4.rename({'Stratification_American Indian or Alaska Native': "AI_AN"}, axis=1)
df4 = df4.rename({'Stratification_Asian or Pacific Islander': "Asian"}, axis=1)
df4 = df4.rename({'Stratification_Black, non-Hispanic': "Black"}, axis=1)
df4 = df4.rename({'Stratification_Hispanic': "Hispanic"}, axis=1)
df4 = df4.rename({'Stratification_White, non-Hispanic': "White"}, axis=1)
df4.head()
#Bayesian Model - poisson distribuition using Race + Year
formula = "DataValue ~ AI_AN + Asian + Black + Hispanic + White"
with pm.Model() as poisson_model:
glm.GLM.from_formula(formula, df4, family=glm.families.Poisson())
# draw posterior samples using NUTS sampling
poisson_trace = pm.sample(1000, cores=2, target_accept=0.97) #, return_inferencedata=True)
with poisson_model:
poisson_ppc = pm.sample_posterior_predictive(poisson_trace)
poisson_ppc['y'] = poisson_ppc['y'] + 0.0
ppc_poisson = az.from_pymc3(poisson_trace, posterior_predictive=poisson_ppc)
az.plot_ppc(ppc_poisson)
plt.xlabel('y = Age-Adjusted Mortality Rates (per 100,000 cases) from Total Cardiovascular Diseases')
plt.title('Bayesian Poisson Regression with SST covariate')
plt.show()
Decision Trees
#predicting mortality rate based on race
from sklearn.tree import DecisionTreeRegressor
race_tree_model = DecisionTreeRegressor()
X_train = train[['Stratification_American Indian or Alaska Native',
'Stratification_Asian or Pacific Islander',
'Stratification_Black, non-Hispanic', 'Stratification_Hispanic',
'Stratification_White, non-Hispanic', 'YearStart']]
X_test = test[['Stratification_American Indian or Alaska Native',
'Stratification_Asian or Pacific Islander',
'Stratification_Black, non-Hispanic', 'Stratification_Hispanic',
'Stratification_White, non-Hispanic', 'YearStart']]
y_train = train['DataValue']
y_test = test['DataValue']
race_tree_model.fit(X_train,y_train)
train["race_tree_pred"] = race_tree_model.predict(X_train)
test["race_tree_pred"] = race_tree_model.predict(X_test)
#RMSE for race predictor
train_rmse = np.mean((train["race_tree_pred"] - train["DataValue"]) ** 2) ** 0.5
test_rmse = np.mean((test["race_tree_pred"] - test["DataValue"]) ** 2) ** 0.5
print("Training set error for decision tree:", train_rmse)
print("Test set error for decision tree: ", test_rmse)
from sklearn.tree import plot_tree
plt.figure(figsize=(16, 4))
plot_tree(race_tree_model);
#predicting mortality rate based on location
from sklearn.tree import DecisionTreeRegressor
loc_tree_model = DecisionTreeRegressor()
y_train = train['DataValue']
y_test = train['DataValue']
X_train = train[['LocationAbbr_AK', 'LocationAbbr_AL', 'LocationAbbr_AR',
'LocationAbbr_AZ', 'LocationAbbr_CA', 'LocationAbbr_CO',
'LocationAbbr_CT', 'LocationAbbr_DC', 'LocationAbbr_DE',
'LocationAbbr_FL', 'LocationAbbr_GA', 'LocationAbbr_HI',
'LocationAbbr_IA', 'LocationAbbr_ID', 'LocationAbbr_IL',
'LocationAbbr_IN', 'LocationAbbr_KS', 'LocationAbbr_KY',
'LocationAbbr_LA', 'LocationAbbr_MA', 'LocationAbbr_MD',
'LocationAbbr_ME', 'LocationAbbr_MI', 'LocationAbbr_MN',
'LocationAbbr_MO', 'LocationAbbr_MS', 'LocationAbbr_MT',
'LocationAbbr_NC', 'LocationAbbr_ND', 'LocationAbbr_NE',
'LocationAbbr_NH', 'LocationAbbr_NJ', 'LocationAbbr_NM',
'LocationAbbr_NV', 'LocationAbbr_NY', 'LocationAbbr_OH',
'LocationAbbr_OK', 'LocationAbbr_OR', 'LocationAbbr_PA',
'LocationAbbr_RI', 'LocationAbbr_SC', 'LocationAbbr_SD',
'LocationAbbr_TN', 'LocationAbbr_TX', 'LocationAbbr_UT',
'LocationAbbr_VA', 'LocationAbbr_VT', 'LocationAbbr_WA',
'LocationAbbr_WI', 'LocationAbbr_WV', 'LocationAbbr_WY', 'YearStart']]
X_test = test[['LocationAbbr_AK', 'LocationAbbr_AL', 'LocationAbbr_AR',
'LocationAbbr_AZ', 'LocationAbbr_CA', 'LocationAbbr_CO',
'LocationAbbr_CT', 'LocationAbbr_DC', 'LocationAbbr_DE',
'LocationAbbr_FL', 'LocationAbbr_GA', 'LocationAbbr_HI',
'LocationAbbr_IA', 'LocationAbbr_ID', 'LocationAbbr_IL',
'LocationAbbr_IN', 'LocationAbbr_KS', 'LocationAbbr_KY',
'LocationAbbr_LA', 'LocationAbbr_MA', 'LocationAbbr_MD',
'LocationAbbr_ME', 'LocationAbbr_MI', 'LocationAbbr_MN',
'LocationAbbr_MO', 'LocationAbbr_MS', 'LocationAbbr_MT',
'LocationAbbr_NC', 'LocationAbbr_ND', 'LocationAbbr_NE',
'LocationAbbr_NH', 'LocationAbbr_NJ', 'LocationAbbr_NM',
'LocationAbbr_NV', 'LocationAbbr_NY', 'LocationAbbr_OH',
'LocationAbbr_OK', 'LocationAbbr_OR', 'LocationAbbr_PA',
'LocationAbbr_RI', 'LocationAbbr_SC', 'LocationAbbr_SD',
'LocationAbbr_TN', 'LocationAbbr_TX', 'LocationAbbr_UT',
'LocationAbbr_VA', 'LocationAbbr_VT', 'LocationAbbr_WA',
'LocationAbbr_WI', 'LocationAbbr_WV', 'LocationAbbr_WY', 'YearStart']]
loc_tree_model.fit(X_train,y_train)
train["loc_tree_pred"] = loc_tree_model.predict(X_train)
test["loc_tree_pred"] = loc_tree_model.predict(X_test)
#RMSE for race predictor
train_rmse = np.mean((train["loc_tree_pred"] - train["DataValue"]) ** 2) ** 0.5
test_rmse = np.mean((test["loc_tree_pred"] - test["DataValue"]) ** 2) ** 0.5
print("Training set error for decision tree:", train_rmse)
print("Test set error for decision tree: ", test_rmse)
from sklearn.tree import plot_tree
plt.figure(figsize=(30, 15))
plot_tree(loc_tree_model);
#closer look of nodes at the top
X_cols = ['LocationAbbr_AK', 'LocationAbbr_AL', 'LocationAbbr_AR',
'LocationAbbr_AZ', 'LocationAbbr_CA', 'LocationAbbr_CO',
'LocationAbbr_CT', 'LocationAbbr_DC', 'LocationAbbr_DE',
'LocationAbbr_FL', 'LocationAbbr_GA', 'LocationAbbr_HI',
'LocationAbbr_IA', 'LocationAbbr_ID', 'LocationAbbr_IL',
'LocationAbbr_IN', 'LocationAbbr_KS', 'LocationAbbr_KY',
'LocationAbbr_LA', 'LocationAbbr_MA', 'LocationAbbr_MD',
'LocationAbbr_ME', 'LocationAbbr_MI', 'LocationAbbr_MN',
'LocationAbbr_MO', 'LocationAbbr_MS', 'LocationAbbr_MT',
'LocationAbbr_NC', 'LocationAbbr_ND', 'LocationAbbr_NE',
'LocationAbbr_NH', 'LocationAbbr_NJ', 'LocationAbbr_NM',
'LocationAbbr_NV', 'LocationAbbr_NY', 'LocationAbbr_OH',
'LocationAbbr_OK', 'LocationAbbr_OR', 'LocationAbbr_PA',
'LocationAbbr_RI', 'LocationAbbr_SC', 'LocationAbbr_SD',
'LocationAbbr_TN', 'LocationAbbr_TX', 'LocationAbbr_UT',
'LocationAbbr_VA', 'LocationAbbr_VT', 'LocationAbbr_WA',
'LocationAbbr_WI', 'LocationAbbr_WV', 'LocationAbbr_WY']
plt.figure(figsize=(12, 7))
plot_tree(loc_tree_model, max_depth=2, fontsize=14, feature_names=X_cols);
#predicting mortality and location rate based on location
from sklearn.tree import DecisionTreeRegressor
loc_race_tree_model = DecisionTreeRegressor()
y_train = train['DataValue']
y_test = train['DataValue']
X_train = train[['Stratification_American Indian or Alaska Native',
'Stratification_Asian or Pacific Islander',
'Stratification_Black, non-Hispanic', 'Stratification_Hispanic',
'Stratification_White, non-Hispanic','LocationAbbr_AK', 'LocationAbbr_AL', 'LocationAbbr_AR',
'LocationAbbr_AZ', 'LocationAbbr_CA', 'LocationAbbr_CO',
'LocationAbbr_CT', 'LocationAbbr_DC', 'LocationAbbr_DE',
'LocationAbbr_FL', 'LocationAbbr_GA', 'LocationAbbr_HI',
'LocationAbbr_IA', 'LocationAbbr_ID', 'LocationAbbr_IL',
'LocationAbbr_IN', 'LocationAbbr_KS', 'LocationAbbr_KY',
'LocationAbbr_LA', 'LocationAbbr_MA', 'LocationAbbr_MD',
'LocationAbbr_ME', 'LocationAbbr_MI', 'LocationAbbr_MN',
'LocationAbbr_MO', 'LocationAbbr_MS', 'LocationAbbr_MT',
'LocationAbbr_NC', 'LocationAbbr_ND', 'LocationAbbr_NE',
'LocationAbbr_NH', 'LocationAbbr_NJ', 'LocationAbbr_NM',
'LocationAbbr_NV', 'LocationAbbr_NY', 'LocationAbbr_OH',
'LocationAbbr_OK', 'LocationAbbr_OR', 'LocationAbbr_PA',
'LocationAbbr_RI', 'LocationAbbr_SC', 'LocationAbbr_SD',
'LocationAbbr_TN', 'LocationAbbr_TX', 'LocationAbbr_UT',
'LocationAbbr_VA', 'LocationAbbr_VT', 'LocationAbbr_WA',
'LocationAbbr_WI', 'LocationAbbr_WV', 'LocationAbbr_WY', 'YearStart']]
X_test = test[['Stratification_American Indian or Alaska Native',
'Stratification_Asian or Pacific Islander',
'Stratification_Black, non-Hispanic', 'Stratification_Hispanic',
'Stratification_White, non-Hispanic','LocationAbbr_AK', 'LocationAbbr_AL', 'LocationAbbr_AR',
'LocationAbbr_AZ', 'LocationAbbr_CA', 'LocationAbbr_CO',
'LocationAbbr_CT', 'LocationAbbr_DC', 'LocationAbbr_DE',
'LocationAbbr_FL', 'LocationAbbr_GA', 'LocationAbbr_HI',
'LocationAbbr_IA', 'LocationAbbr_ID', 'LocationAbbr_IL',
'LocationAbbr_IN', 'LocationAbbr_KS', 'LocationAbbr_KY',
'LocationAbbr_LA', 'LocationAbbr_MA', 'LocationAbbr_MD',
'LocationAbbr_ME', 'LocationAbbr_MI', 'LocationAbbr_MN',
'LocationAbbr_MO', 'LocationAbbr_MS', 'LocationAbbr_MT',
'LocationAbbr_NC', 'LocationAbbr_ND', 'LocationAbbr_NE',
'LocationAbbr_NH', 'LocationAbbr_NJ', 'LocationAbbr_NM',
'LocationAbbr_NV', 'LocationAbbr_NY', 'LocationAbbr_OH',
'LocationAbbr_OK', 'LocationAbbr_OR', 'LocationAbbr_PA',
'LocationAbbr_RI', 'LocationAbbr_SC', 'LocationAbbr_SD',
'LocationAbbr_TN', 'LocationAbbr_TX', 'LocationAbbr_UT',
'LocationAbbr_VA', 'LocationAbbr_VT', 'LocationAbbr_WA',
'LocationAbbr_WI', 'LocationAbbr_WV', 'LocationAbbr_WY', 'YearStart']]
loc_race_tree_model.fit(X_train,y_train)
train["loc_race_tree_pred"] = loc_race_tree_model.predict(X_train)
test["loc_race_tree_pred"] = loc_race_tree_model.predict(X_test)
#RMSE for race predictor
train_rmse = np.mean((train["loc_race_tree_pred"] - train["DataValue"]) ** 2) ** 0.5
test_rmse = np.mean((test["loc_race_tree_pred"] - test["DataValue"]) ** 2) ** 0.5
print("Training set error for decision tree:", train_rmse)
print("Test set error for decision tree: ", test_rmse)
from sklearn.tree import plot_tree
plt.figure(figsize=(30, 15))
plot_tree(loc_tree_model);
#closer look at nodes
X_cols = ['Stratification_American Indian or Alaska Native',
'Stratification_Asian or Pacific Islander',
'Stratification_Black, non-Hispanic', 'Stratification_Hispanic',
'Stratification_White, non-Hispanic','LocationAbbr_AK', 'LocationAbbr_AL','LocationAbbr_AR',
'LocationAbbr_AZ', 'LocationAbbr_CA', 'LocationAbbr_CO',
'LocationAbbr_CT', 'LocationAbbr_DC', 'LocationAbbr_DE',
'LocationAbbr_FL', 'LocationAbbr_GA', 'LocationAbbr_HI',
'LocationAbbr_IA', 'LocationAbbr_ID', 'LocationAbbr_IL',
'LocationAbbr_IN', 'LocationAbbr_KS', 'LocationAbbr_KY',
'LocationAbbr_LA', 'LocationAbbr_MA', 'LocationAbbr_MD',
'LocationAbbr_ME', 'LocationAbbr_MI', 'LocationAbbr_MN',
'LocationAbbr_MO', 'LocationAbbr_MS', 'LocationAbbr_MT',
'LocationAbbr_NC', 'LocationAbbr_ND', 'LocationAbbr_NE',
'LocationAbbr_NH', 'LocationAbbr_NJ', 'LocationAbbr_NM',
'LocationAbbr_NV', 'LocationAbbr_NY', 'LocationAbbr_OH',
'LocationAbbr_OK', 'LocationAbbr_OR', 'LocationAbbr_PA',
'LocationAbbr_RI', 'LocationAbbr_SC', 'LocationAbbr_SD',
'LocationAbbr_TN', 'LocationAbbr_TX', 'LocationAbbr_UT',
'LocationAbbr_VA', 'LocationAbbr_VT', 'LocationAbbr_WA',
'LocationAbbr_WI', 'LocationAbbr_WV', 'LocationAbbr_WY']
plt.figure(figsize=(20, 15))
plot_tree(loc_tree_model, max_depth=7, fontsize=14, feature_names=X_cols);
df2
## OLS model from lab rather than sm library
def fit_OLS_model(df, target_variable, explanatory_variables, intercept = True):
"""
fits an OLS model from input data
Inputs:
df: pandas DataFrame
target_variable: string, name of the target variable
explanatory_variables: list of strings, names of the explanatory variables
intercept: bool, if True add intercept term
Outputs:
fitted_model: model containing OLS regression results
"""
target = df[target_variable]
inputs = df[explanatory_variables]
if intercept:
inputs = sm.add_constant(inputs)
fitted_model = sm.OLS(target, inputs).fit()
return(fitted_model)
OLS_X = df2.drop('DataValue', axis = 1)
OLS_Y_series = pd.to_numeric(df['DataValue'])
OLS_Y = OLS_Y_series.to_frame()
OLS_Y.head()
OLS_X.head()
#STRATIFICATION / INDEPENDENT COLUMNS
X_AIAN = OLS_X['Stratification_American Indian or Alaska Native']
X_API = OLS_X['Stratification_Asian or Pacific Islander']
X_BNH = OLS_X['Stratification_Black, non-Hispanic']
X_HIS = OLS_X['Stratification_Hispanic']
X_WNH = OLS_X['Stratification_White, non-Hispanic']
# TRYING TO COMBINE TO 1 SUMMARY WITH LAB METHOD
# OLS_X = df2.drop('DataValue', axis = 1)
# dv_numeric = pd.to_numeric(df['DataValue'])
# OLS_Y = OLS_Y_series.to_frame()
# OLS_Y.head()
df69 = df2
df69['DataValue'] = OLS_Y['DataValue']
AIAN = 'Stratification_American Indian or Alaska Native'
API = 'Stratification_Asian or Pacific Islander'
BNH = 'Stratification_Black, non-Hispanic'
HIS = 'Stratification_Hispanic'
WNH = 'Stratification_White, non-Hispanic'
df69.head()
gammas_1 = fit_OLS_model(df69, 'DataValue', ['Stratification_American Indian or Alaska Native', 'Stratification_Asian or Pacific Islander', 'Stratification_Black, non-Hispanic', 'Stratification_Hispanic', 'Stratification_White, non-Hispanic'])
print(gammas_1.summary())
print("The estimated causal effect on 'Mortality from total cardiovascular diseases' of 'Stratification_American Indian or Alaska Native' is {:.2f}".format(gammas_1.params[0]))
print("The estimated causal effect on 'Mortality from total cardiovascular diseases' of 'Stratification_Asian or Pacific Islander' is {:.2f}".format(gammas_1.params[1]))
# print("The estimated causal effect on 'Mortality from total cardiovascular diseases' of 'Stratification_Black, non-Hispanic' is {:.2f}".format(OLS_model3.params[0]))
# print("The estimated causal effect on 'Mortality from total cardiovascular diseases' of 'Stratification_Hispanic' is {:.2f}".format(OLS_model4.params[0]))
# print("The estimated causal effect on 'Mortality from total cardiovascular diseases' of 'Stratification_White, non-Hispanic' is {:.2f}".format(OLS_model5.params[0]))
#Notes/Old Work
# gammas_cols = df3['Stratification_American Indian or Alaska Native']
# cvd_rate_val = pd.to_numeric(df3['DataValue'])
# gammas_model = fit_OLS_model(df3, cvd_rate_val, gammas_cols)
# print(gammas_model.summary())
# ols_mod = sm.OLS(df2['Stratification_American Indian or Alaska Native'], df2['DataValue'])
# OLS Linear Model 1
# Dependent Var / Target: Mortality from CVDs (DataValue)
# Independent Var / Explanatory: Stratification_American Indian or Alaska Native
OLS_model1 = sm.regression.linear_model.OLS(OLS_Y, X_AIAN).fit()
print(OLS_model1.summary())
# OLS_model2 = sm.regression.linear_model.OLS(OLS_Y, X_API).fit()
# print(OLS_model2.summary())
# OLS_model3 = sm.regression.linear_model.OLS(OLS_Y, X_BNH).fit()
# print(OLS_model3.summary())
# OLS_model4 = sm.regression.linear_model.OLS(OLS_Y, X_HIS).fit()
# print(OLS_model4.summary())
# OLS_model5 = sm.regression.linear_model.OLS(OLS_Y, X_WNH).fit()
# print(OLS_model5.summary())
Random Forest
from sklearn.ensemble import RandomForestRegressor
race_forest_model = RandomForestRegressor(max_features = 1)
X_train = train[['Stratification_American Indian or Alaska Native',
'Stratification_Asian or Pacific Islander',
'Stratification_Black, non-Hispanic', 'Stratification_Hispanic',
'Stratification_White, non-Hispanic']]
X_test = test[['Stratification_American Indian or Alaska Native',
'Stratification_Asian or Pacific Islander',
'Stratification_Black, non-Hispanic', 'Stratification_Hispanic',
'Stratification_White, non-Hispanic']]
y_train = train['DataValue']
y_test = test['DataValue']
race_forest_model.fit(X_train,y_train)
train_pred = race_forest_model.predict(X_train)
test_pred = race_forest_model.predict(X_test)
train["race_forest_pred"] = train_pred
test["race_forest_pred"] = test_pred
train_rmse = np.mean((train["race_forest_pred"] - train["DataValue"]) ** 2) ** 0.5
test_rmse = np.mean((test["race_forest_pred"] - test["DataValue"]) ** 2) ** 0.5
print("Training set error for random forest:", train_rmse)
print("Test set error for random forest: ", test_rmse)
worst_pred_forest = test
worst_pred_forest["race absolute error"] = np.abs((worst_pred_forest["DataValue"] - worst_pred_forest['race_forest_pred']))
worst_forest_predicted_cars_df = worst_pred_forest.sort_values("race absolute error", ascending=False).iloc[0:2]
worst_forest_predicted_cars_df
from sklearn.ensemble import RandomForestRegressor
loc_forest_model = RandomForestRegressor(max_features = 1)
X_train = train[['LocationAbbr_AK', 'LocationAbbr_AL', 'LocationAbbr_AR',
'LocationAbbr_AZ', 'LocationAbbr_CA', 'LocationAbbr_CO',
'LocationAbbr_CT', 'LocationAbbr_DC', 'LocationAbbr_DE',
'LocationAbbr_FL', 'LocationAbbr_GA', 'LocationAbbr_HI',
'LocationAbbr_IA', 'LocationAbbr_ID', 'LocationAbbr_IL',
'LocationAbbr_IN', 'LocationAbbr_KS', 'LocationAbbr_KY',
'LocationAbbr_LA', 'LocationAbbr_MA', 'LocationAbbr_MD',
'LocationAbbr_ME', 'LocationAbbr_MI', 'LocationAbbr_MN',
'LocationAbbr_MO', 'LocationAbbr_MS', 'LocationAbbr_MT',
'LocationAbbr_NC', 'LocationAbbr_ND', 'LocationAbbr_NE',
'LocationAbbr_NH', 'LocationAbbr_NJ', 'LocationAbbr_NM',
'LocationAbbr_NV', 'LocationAbbr_NY', 'LocationAbbr_OH',
'LocationAbbr_OK', 'LocationAbbr_OR', 'LocationAbbr_PA',
'LocationAbbr_RI', 'LocationAbbr_SC', 'LocationAbbr_SD',
'LocationAbbr_TN', 'LocationAbbr_TX', 'LocationAbbr_UT',
'LocationAbbr_VA', 'LocationAbbr_VT', 'LocationAbbr_WA',
'LocationAbbr_WI', 'LocationAbbr_WV', 'LocationAbbr_WY']]
X_test = test[['LocationAbbr_AK', 'LocationAbbr_AL', 'LocationAbbr_AR',
'LocationAbbr_AZ', 'LocationAbbr_CA', 'LocationAbbr_CO',
'LocationAbbr_CT', 'LocationAbbr_DC', 'LocationAbbr_DE',
'LocationAbbr_FL', 'LocationAbbr_GA', 'LocationAbbr_HI',
'LocationAbbr_IA', 'LocationAbbr_ID', 'LocationAbbr_IL',
'LocationAbbr_IN', 'LocationAbbr_KS', 'LocationAbbr_KY',
'LocationAbbr_LA', 'LocationAbbr_MA', 'LocationAbbr_MD',
'LocationAbbr_ME', 'LocationAbbr_MI', 'LocationAbbr_MN',
'LocationAbbr_MO', 'LocationAbbr_MS', 'LocationAbbr_MT',
'LocationAbbr_NC', 'LocationAbbr_ND', 'LocationAbbr_NE',
'LocationAbbr_NH', 'LocationAbbr_NJ', 'LocationAbbr_NM',
'LocationAbbr_NV', 'LocationAbbr_NY', 'LocationAbbr_OH',
'LocationAbbr_OK', 'LocationAbbr_OR', 'LocationAbbr_PA',
'LocationAbbr_RI', 'LocationAbbr_SC', 'LocationAbbr_SD',
'LocationAbbr_TN', 'LocationAbbr_TX', 'LocationAbbr_UT',
'LocationAbbr_VA', 'LocationAbbr_VT', 'LocationAbbr_WA',
'LocationAbbr_WI', 'LocationAbbr_WV', 'LocationAbbr_WY']]
y_train = train['DataValue']
y_test = test['DataValue']
loc_forest_model.fit(X_train,y_train)
train_pred = loc_forest_model.predict(X_train)
test_pred = loc_forest_model.predict(X_test)
train["loc_forest_pred"] = train_pred
test["loc_forest_pred"] = test_pred
train_rmse = np.mean((train["loc_forest_pred"] - train["DataValue"]) ** 2) ** 0.5
test_rmse = np.mean((test["loc_forest_pred"] - test["DataValue"]) ** 2) ** 0.5
print("Training set error for random forest:", train_rmse)
print("Test set error for random forest: ", test_rmse)
worst_pred_forest = test
worst_pred_forest['loc absolute error'] = np.abs((worst_pred_forest["DataValue"] - worst_pred_forest['loc_forest_pred']))
worst_forest_predicted_cars_df = worst_pred_forest.sort_values("loc absolute error", ascending=False).iloc[0:2]
worst_forest_predicted_cars_df
from sklearn.ensemble import RandomForestRegressor
loc_race_forest_model = RandomForestRegressor(max_features = 1)
X_train = train[['LocationAbbr_AK', 'LocationAbbr_AL', 'LocationAbbr_AR',
'LocationAbbr_AZ', 'LocationAbbr_CA', 'LocationAbbr_CO',
'LocationAbbr_CT', 'LocationAbbr_DC', 'LocationAbbr_DE',
'LocationAbbr_FL', 'LocationAbbr_GA', 'LocationAbbr_HI',
'LocationAbbr_IA', 'LocationAbbr_ID', 'LocationAbbr_IL',
'LocationAbbr_IN', 'LocationAbbr_KS', 'LocationAbbr_KY',
'LocationAbbr_LA', 'LocationAbbr_MA', 'LocationAbbr_MD',
'LocationAbbr_ME', 'LocationAbbr_MI', 'LocationAbbr_MN',
'LocationAbbr_MO', 'LocationAbbr_MS', 'LocationAbbr_MT',
'LocationAbbr_NC', 'LocationAbbr_ND', 'LocationAbbr_NE',
'LocationAbbr_NH', 'LocationAbbr_NJ', 'LocationAbbr_NM',
'LocationAbbr_NV', 'LocationAbbr_NY', 'LocationAbbr_OH',
'LocationAbbr_OK', 'LocationAbbr_OR', 'LocationAbbr_PA',
'LocationAbbr_RI', 'LocationAbbr_SC', 'LocationAbbr_SD',
'LocationAbbr_TN', 'LocationAbbr_TX', 'LocationAbbr_UT',
'LocationAbbr_VA', 'LocationAbbr_VT', 'LocationAbbr_WA',
'LocationAbbr_WI', 'LocationAbbr_WV', 'LocationAbbr_WY',
'Stratification_American Indian or Alaska Native',
'Stratification_Asian or Pacific Islander',
'Stratification_Black, non-Hispanic', 'Stratification_Hispanic',
'Stratification_White, non-Hispanic']]
X_test = test[['LocationAbbr_AK', 'LocationAbbr_AL', 'LocationAbbr_AR',
'LocationAbbr_AZ', 'LocationAbbr_CA', 'LocationAbbr_CO',
'LocationAbbr_CT', 'LocationAbbr_DC', 'LocationAbbr_DE',
'LocationAbbr_FL', 'LocationAbbr_GA', 'LocationAbbr_HI',
'LocationAbbr_IA', 'LocationAbbr_ID', 'LocationAbbr_IL',
'LocationAbbr_IN', 'LocationAbbr_KS', 'LocationAbbr_KY',
'LocationAbbr_LA', 'LocationAbbr_MA', 'LocationAbbr_MD',
'LocationAbbr_ME', 'LocationAbbr_MI', 'LocationAbbr_MN',
'LocationAbbr_MO', 'LocationAbbr_MS', 'LocationAbbr_MT',
'LocationAbbr_NC', 'LocationAbbr_ND', 'LocationAbbr_NE',
'LocationAbbr_NH', 'LocationAbbr_NJ', 'LocationAbbr_NM',
'LocationAbbr_NV', 'LocationAbbr_NY', 'LocationAbbr_OH',
'LocationAbbr_OK', 'LocationAbbr_OR', 'LocationAbbr_PA',
'LocationAbbr_RI', 'LocationAbbr_SC', 'LocationAbbr_SD',
'LocationAbbr_TN', 'LocationAbbr_TX', 'LocationAbbr_UT',
'LocationAbbr_VA', 'LocationAbbr_VT', 'LocationAbbr_WA',
'LocationAbbr_WI', 'LocationAbbr_WV', 'LocationAbbr_WY',
'Stratification_American Indian or Alaska Native',
'Stratification_Asian or Pacific Islander',
'Stratification_Black, non-Hispanic', 'Stratification_Hispanic',
'Stratification_White, non-Hispanic']]
y_train = train['DataValue']
y_test = test['DataValue']
loc_race_forest_model.fit(X_train,y_train)
train_pred = loc_race_forest_model.predict(X_train)
test_pred = loc_race_forest_model.predict(X_test)
train["loc_race_forest_pred"] = train_pred
test["loc_race_forest_pred"] = test_pred
train_rmse = np.mean((train["loc_race_forest_pred"] - train["DataValue"]) ** 2) ** 0.5
test_rmse = np.mean((test["loc_race_forest_pred"] - test["DataValue"]) ** 2) ** 0.5
print("Training set error for random forest:", train_rmse)
print("Test set error for random forest: ", test_rmse)
worst_pred_forest = test
worst_pred_forest['loc and race absolute error'] = np.abs((worst_pred_forest["DataValue"] - worst_pred_forest['loc_race_forest_pred']))
worst_forest_predicted_cars_df = worst_pred_forest.sort_values("loc and race absolute error", ascending=False).iloc[0:2]
worst_forest_predicted_cars_df
from sklearn.ensemble import RandomForestRegressor
year_loc_race_forest_model = RandomForestRegressor(max_features = 1)
X_train = train[['LocationAbbr_AK', 'LocationAbbr_AL', 'LocationAbbr_AR',
'LocationAbbr_AZ', 'LocationAbbr_CA', 'LocationAbbr_CO',
'LocationAbbr_CT', 'LocationAbbr_DC', 'LocationAbbr_DE',
'LocationAbbr_FL', 'LocationAbbr_GA', 'LocationAbbr_HI',
'LocationAbbr_IA', 'LocationAbbr_ID', 'LocationAbbr_IL',
'LocationAbbr_IN', 'LocationAbbr_KS', 'LocationAbbr_KY',
'LocationAbbr_LA', 'LocationAbbr_MA', 'LocationAbbr_MD',
'LocationAbbr_ME', 'LocationAbbr_MI', 'LocationAbbr_MN',
'LocationAbbr_MO', 'LocationAbbr_MS', 'LocationAbbr_MT',
'LocationAbbr_NC', 'LocationAbbr_ND', 'LocationAbbr_NE',
'LocationAbbr_NH', 'LocationAbbr_NJ', 'LocationAbbr_NM',
'LocationAbbr_NV', 'LocationAbbr_NY', 'LocationAbbr_OH',
'LocationAbbr_OK', 'LocationAbbr_OR', 'LocationAbbr_PA',
'LocationAbbr_RI', 'LocationAbbr_SC', 'LocationAbbr_SD',
'LocationAbbr_TN', 'LocationAbbr_TX', 'LocationAbbr_UT',
'LocationAbbr_VA', 'LocationAbbr_VT', 'LocationAbbr_WA',
'LocationAbbr_WI', 'LocationAbbr_WV', 'LocationAbbr_WY',
'Stratification_American Indian or Alaska Native',
'Stratification_Asian or Pacific Islander',
'Stratification_Black, non-Hispanic', 'Stratification_Hispanic',
'Stratification_White, non-Hispanic', 'YearStart']]
X_test = test[['LocationAbbr_AK', 'LocationAbbr_AL', 'LocationAbbr_AR',
'LocationAbbr_AZ', 'LocationAbbr_CA', 'LocationAbbr_CO',
'LocationAbbr_CT', 'LocationAbbr_DC', 'LocationAbbr_DE',
'LocationAbbr_FL', 'LocationAbbr_GA', 'LocationAbbr_HI',
'LocationAbbr_IA', 'LocationAbbr_ID', 'LocationAbbr_IL',
'LocationAbbr_IN', 'LocationAbbr_KS', 'LocationAbbr_KY',
'LocationAbbr_LA', 'LocationAbbr_MA', 'LocationAbbr_MD',
'LocationAbbr_ME', 'LocationAbbr_MI', 'LocationAbbr_MN',
'LocationAbbr_MO', 'LocationAbbr_MS', 'LocationAbbr_MT',
'LocationAbbr_NC', 'LocationAbbr_ND', 'LocationAbbr_NE',
'LocationAbbr_NH', 'LocationAbbr_NJ', 'LocationAbbr_NM',
'LocationAbbr_NV', 'LocationAbbr_NY', 'LocationAbbr_OH',
'LocationAbbr_OK', 'LocationAbbr_OR', 'LocationAbbr_PA',
'LocationAbbr_RI', 'LocationAbbr_SC', 'LocationAbbr_SD',
'LocationAbbr_TN', 'LocationAbbr_TX', 'LocationAbbr_UT',
'LocationAbbr_VA', 'LocationAbbr_VT', 'LocationAbbr_WA',
'LocationAbbr_WI', 'LocationAbbr_WV', 'LocationAbbr_WY',
'Stratification_American Indian or Alaska Native',
'Stratification_Asian or Pacific Islander',
'Stratification_Black, non-Hispanic', 'Stratification_Hispanic',
'Stratification_White, non-Hispanic', 'YearStart']]
y_train = train['DataValue']
y_test = test['DataValue']
year_loc_race_forest_model.fit(X_train,y_train)
train_pred = year_loc_race_forest_model.predict(X_train)
test_pred = year_loc_race_forest_model.predict(X_test)
train["year_loc_race_forest_pred"] = train_pred
test["year_loc_race_forest_pred"] = test_pred
train_rmse = np.mean((train["year_loc_race_forest_pred"] - train["DataValue"]) ** 2) ** 0.5
test_rmse = np.mean((test["year_loc_race_forest_pred"] - test["DataValue"]) ** 2) ** 0.5
print("Training set error for random forest:", train_rmse)
print("Test set error for random forest:", test_rmse)
worst_pred_forest = test
worst_pred_forest['yr loc and race absolute error'] = np.abs((worst_pred_forest["DataValue"] - worst_pred_forest['year_loc_race_forest_pred']))
worst_forest_predicted_cars_df = worst_pred_forest.sort_values("yr loc and race absolute error", ascending=False).iloc[0:2]
worst_forest_predicted_cars_df