import pandas as pd
import numpy as np
import warnings
warnings.filterwarnings("ignore")
import warnings
warnings.filterwarnings("ignore", category=DeprecationWarning)

import matplotlib.pyplot as plt
import seaborn as sns
sns.set()

covid_vacc_zipcode = pd.read_csv('people/coverage-by-modzcta-allages.csv')
covid_vacc_zipcode.head()

import pgeocode

nomi = pgeocode.Nominatim('US')
# get latitude and longitude from zip code
covid_vacc_zipcode[['LATITUDE', 'LONGITUDE']] = covid_vacc_zipcode['MODZCTA'].apply(lambda x: nomi.query_postal_code(str(x))[-3:-1])

covid_vacc_zipcode.head()

# get education by zipcodes
educational_attainment = pd.read_csv('people/ACSST5Y2019.S1501_data_with_overlays_2021-06-22T105606.csv', header=1)

# get cols we need for analysis
cols = ['id',
        'Geographic Area Name',
        'Estimate!!Total!!AGE BY EDUCATIONAL ATTAINMENT!!Population 18 to 24 years!!Less than high school graduate',
       'Estimate!!Total!!AGE BY EDUCATIONAL ATTAINMENT!!Population 18 to 24 years!!High school graduate (includes equivalency)',
       'Estimate!!Total!!AGE BY EDUCATIONAL ATTAINMENT!!Population 18 to 24 years!!Some college or associate\'s degree',
       'Estimate!!Total!!AGE BY EDUCATIONAL ATTAINMENT!!Population 18 to 24 years!!Bachelor\'s degree or higher',
       'Estimate!!Total!!AGE BY EDUCATIONAL ATTAINMENT!!Population 25 years and over!!High school graduate (includes equivalency)',
       'Estimate!!Total!!AGE BY EDUCATIONAL ATTAINMENT!!Population 25 years and over!!Some college, no degree',
       'Estimate!!Total!!AGE BY EDUCATIONAL ATTAINMENT!!Population 25 years and over!!Associate\'s degree',
       'Estimate!!Total!!AGE BY EDUCATIONAL ATTAINMENT!!Population 25 years and over!!Bachelor\'s degree',
       'Estimate!!Total!!AGE BY EDUCATIONAL ATTAINMENT!!Population 25 years and over!!Graduate or professional degree']

# get filtered dataset
educational_attainment = educational_attainment[cols]

# column filter
def filter_column(df, expressions):
    '''
    df - dataframe
    expression - str, boolean logic operators
                e.g 'Dog|dog|cat|Cat'
    '''
    return df.filter(regex=expressions)

# educational attainment
less_than_high_school = filter_column(educational_attainment, 'Less than').sum(axis=1)
high_school = filter_column(educational_attainment,'High').sum(axis=1)
some_college_associates = filter_column(educational_attainment, 'Some|Associate').sum(axis=1)
bachelors = filter_column(educational_attainment, 'Bachelor').sum(axis=1)
graduate_professional_deg = filter_column(educational_attainment, 'Graduate').sum(axis=1)

# create dataframe for educational attainment
degrees = pd.DataFrame({'raw_zip_code': educational_attainment['Geographic Area Name'],'less_than_high_school': less_than_high_school, 'high_school':high_school, 'some_college_associates':some_college_associates, 'bachelors':bachelors, 'graduate_professional_deg':graduate_professional_deg})


# clean zip code column
degrees.loc[:, 'raw_zip_code'] = degrees.apply(lambda row: row.raw_zip_code[-5:].strip(), axis=1).astype('int')

# merge covid dataset and education attainment on zip code
covid_vacc_zipcode_education = covid_vacc_zipcode.merge(degrees, how='inner', left_on='MODZCTA', right_on='raw_zip_code')
covid_vacc_zipcode_education.head()


# calculate educational attainment percentage by zipcode
# get percentage of educational attainment in neighborhood
covid_vacc_zipcode_education[['less_than_high_school %', 'high_school %', 'some_college_associates %', 'bachelors %', 'graduate_professional_deg %']] = covid_vacc_zipcode_education[['less_than_high_school', 'high_school', 'some_college_associates', 'bachelors', 'graduate_professional_deg']].div(covid_vacc_zipcode_education['POP_DENOMINATOR'], axis=0)
covid_vacc_zipcode_education.head()

# median income
median_income_nyc = pd.read_csv('people/Median Incomes.csv')

# rename data column to income
median_income_nyc.rename(columns={'Data':'Income'}, inplace=True)


median_income_nyc[median_income_nyc['Location'].str[0]=='Z'].index[0]
# find zip code by filtering by Z
# clean zip code
median_income_nyc.loc[3900:, 'Location'] = median_income_nyc[median_income_nyc['Location'].str[0]=='Z'].apply(lambda zc: zc.Location[-5:].strip(), axis=1)

# get zip code only
median_income_nyc = median_income_nyc.loc[3900:, :].copy()

# convert income column to float
median_income_nyc['Income'] = median_income_nyc['Income'].astype('float')

# convert object to int
median_income_nyc.loc[:, 'Location'] = median_income_nyc.loc[:, 'Location'].astype('int')

# filter to get "All Household" and most recent median household income which is 2019
# then merge dataset to covid_vacc_zipcode_education
covid_vacc_zipcode_education_income = covid_vacc_zipcode_education.merge(median_income_nyc.loc[(median_income_nyc['Household Type']=='All Households') & (median_income_nyc['TimeFrame']==2019), :], how='inner', left_on='MODZCTA', right_on='Location')

covid_vacc_zipcode_education_income.head()

covid_death = pd.read_csv('coronavirus-data/totals/data-by-modzcta.csv')

covid_death.head()

# convert to int
covid_death['MODIFIED_ZCTA'] = covid_death['MODIFIED_ZCTA'].astype('int')

# merge datasets
covid_vacc_zipcode_education_income_deaths = covid_vacc_zipcode_education_income.merge(covid_death, how='inner', left_on='MODZCTA', right_on='MODIFIED_ZCTA')

# this dataset contains lat and longitude so we could've just used that for tableau
covid_vacc_zipcode_education_income_deaths.head()

# fortunately, we have no missing datapoints
covid_vacc_zipcode_education_income_deaths.isna().sum()

DATE                           0
NEIGHBORHOOD_NAME_x            0
BOROUGH                        0
MODZCTA                        0
Label                          0
AGE_GROUP                      0
POP_DENOMINATOR_x              0
COUNT_PARTIALLY_CUMULATIVE     0
COUNT_FULLY_CUMULATIVE         0
COUNT_1PLUS_CUMULATIVE         0
PERC_PARTIALLY                 0
PERC_FULLY                     0
PERC_1PLUS                     0
LATITUDE                       0
LONGITUDE                      0
raw_zip_code                   0
less_than_high_school          0
high_school                    0
some_college_associates        0
bachelors                      0
graduate_professional_deg      0
less_than_high_school %        0
high_school %                  0
some_college_associates %      0
bachelors %                    0
graduate_professional_deg %    0
Location                       0
Household Type                 0
TimeFrame                      0
DataFormat                     0
Income                         0
Fips                           0
MODIFIED_ZCTA                  0
NEIGHBORHOOD_NAME_y            0
BOROUGH_GROUP                  0
label                          0
lat                            0
lon                            0
COVID_CASE_COUNT               0
COVID_CASE_RATE                0
POP_DENOMINATOR_y              0
COVID_DEATH_COUNT              0
COVID_DEATH_RATE               0
PERCENT_POSITIVE               0
TOTAL_COVID_TESTS              0
dtype: int64

plt.figure(figsize=(20,15))
scatter = plt.scatter(covid_vacc_zipcode_education_income_deaths['PERC_FULLY'], covid_vacc_zipcode_education_income_deaths['Income'], c=covid_vacc_zipcode_education_income_deaths['Income'], s=covid_vacc_zipcode_education_income_deaths['COVID_DEATH_COUNT'])
plt.colorbar()
plt.xlabel('Fully Vacc %')
plt.ylabel('Income')
plt.title('Income vs. COVID 19 Full Vaccinated')

handles, labels = scatter.legend_elements(prop="sizes", alpha=0.6)
legend2 = plt.legend(handles, labels, loc="upper left", title="COVID-19 Death Count")
plt.show()

plt.figure(figsize=(20,15))
scatter = plt.scatter(covid_vacc_zipcode_education_income_deaths['PERC_1PLUS'], covid_vacc_zipcode_education_income_deaths['Income'], c=covid_vacc_zipcode_education_income_deaths['Income'], s=covid_vacc_zipcode_education_income_deaths['COVID_DEATH_COUNT'])
plt.colorbar()
plt.xlabel('1 Plus Vacc %')
plt.ylabel('Income')
plt.title('Income vs. COVID 19 1 Plus')

handles, labels = scatter.legend_elements(prop="sizes", alpha=0.6)

legend2 = plt.legend(handles, labels, loc="upper left", title="COVID-19 Death Count")
plt.show()

plt.figure(figsize=(10,10))
scatter = plt.scatter(covid_vacc_zipcode_education_income_deaths['COVID_DEATH_RATE'], covid_vacc_zipcode_education_income_deaths['Income'], c=covid_vacc_zipcode_education_income_deaths['Income'], s=covid_vacc_zipcode_education_income_deaths['PERC_FULLY'])
plt.colorbar()
plt.xlabel('Death Rate')
plt.ylabel('Income')
plt.title('Income vs. COVID 19 Death Rate')

handles, labels = scatter.legend_elements(prop="sizes", alpha=0.6)
legend2 = plt.legend(handles, labels, loc="upper right", title="COVID-19 FULLY VACC %")
plt.show()

import geopandas as gpd
from geopandas import GeoDataFrame
from shapely.geometry import Point

nyc_neighborhood = gpd.read_file('Neighborhood Tabulation Areas (NTA)/geo_export_ab6fc730-beb8-4922-8a7a-24de1a3fdc50.shp')

geometry = [Point(xy) for xy in zip(covid_vacc_zipcode_education_income_deaths['LONGITUDE'], covid_vacc_zipcode_education_income_deaths['LATITUDE'])]
gdf = GeoDataFrame(covid_vacc_zipcode_education_income_deaths, geometry=geometry)

import matplotlib as mpl
cmap = mpl.colors.ListedColormap(['red', 'orange',
                                  'green'])

import contextily as ctx


ax = nyc_neighborhood.boundary.plot(figsize=(20,15), color='grey')

plt.axis('off')
plt.title('COVID-19 VACCINATION NYC', fontsize=20)

# # plt.legend(labels=covid_vacc_income_education_deaths['COVID_DEATH_COUNT'].sort_values(ascending=True), scatterpoints=1)

# plot scatter again on shapefile to get size legend
scatter = plt.scatter('LONGITUDE', 'LATITUDE', data=covid_vacc_zipcode_education_income_deaths, s='COVID_DEATH_COUNT', c='PERC_FULLY',cmap=cmap)
plt.legend(*scatter.legend_elements("sizes", num=6, alpha=0.6), title='COVID-19 DEATH COUNT')

plt.colorbar(cmap=cmap, boundaries=covid_vacc_zipcode_education_income_deaths['PERC_FULLY'].sort_values(ascending=True).to_list())

ctx.add_basemap(ax=ax, crs='epsg:4326')

from sklearn.model_selection import train_test_split
import re
# create functions for training regression models
# split into test and train datasets
def modelBuilder(X, y, seed=7, test_size=0.1):
    '''returns X_train, X_test, y_train, y_test '''
    seed = seed
    test_size = test_size
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=test_size, random_state=seed)

    return X_train, y_train, X_test, y_test

def modelsTrainer(X_train, y_train, X_test, y_test, models):
    '''Returns results'''
    # create an empty list for results
    results = []
    # loop over models
    for model in models:
        model_regression = model
        # if value error, make sure function is taking correct parameters in the correct order
        model_regression.fit(X_train, y_train)
        y_pred = model_regression.predict(X_test)
        # bias
        training_score = model_regression.score(X_train, y_train)

        # variance
        testing_score = model_regression.score(X_test, y_test)

        # percent difference between training and testing score
        delta = training_score - testing_score
        # get raw model name
        model_name = re.search('[a-zA-Z]+', str(model))[0]
        # append scores to results
        results.append({model_name:[training_score, testing_score, delta]})


    return results

from sklearn.ensemble import RandomForestRegressor
# import all of the models we want to test
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import ElasticNet
from sklearn.linear_model import Lasso, Ridge
from xgboost import XGBRegressor

X = covid_vacc_zipcode_education_income_deaths[['POP_DENOMINATOR_x', 'Income','bachelors', 'COVID_DEATH_COUNT']]
y = covid_vacc_zipcode_education_income_deaths['COUNT_FULLY_CUMULATIVE']


X_train, y_train, X_test, y_test = modelBuilder(X, y)

# models to train
models = [LinearRegression(), RandomForestRegressor(), ElasticNet(), GradientBoostingRegressor(), XGBRegressor(verbosity=0), Lasso(), Ridge()]

# fit models
# returns training score, testing score and difference between the two
modelsTrainer(X_train, y_train, X_test, y_test, models)

[{'LinearRegression': [0.8349011295818145,
   0.805264936984036,
   0.0296361925977785]},
 {'RandomForestRegressor': [0.9393786736093015,
   0.7816819810432987,
   0.15769669256600283]},
 {'ElasticNet': [0.8349011295434671,
   0.8052659443078068,
   0.029635185235660222]},
 {'GradientBoostingRegressor': [0.9854864931948837,
   0.8559366915021283,
   0.12954980169275543]},
 {'XGBRegressor': [0.9781925378371675,
   0.8717985612917483,
   0.10639397654541927]},
 {'Lasso': [0.8349011295798674, 0.8052651639636368,
0.02963596561623061]},
 {'Ridge': [0.8349011295818097, 0.805264948229599,
0.02963618135221069]}]

from sklearn.model_selection import cross_val_score, KFold

cv = KFold(n_splits=10)

# we can use cross-validation to confirm our linear models have low variance and does well generalizing data
for model in models:
  cross_val2 = cross_val_score(model, X_train, y_train, cv=cv)
  # linear models have low standard deviation meaning their variance is low
  # we can conclude our linear model(s) will do good generalizing data
  print(re.search('[a-zA-Z]+', str(model))[0], '-', 'mean:', cross_val2.mean(), 'std:', cross_val2.std())

LinearRegression - mean: 0.7988492776333851 std: 0.09774188277510705
RandomForestRegressor - mean: 0.7321956559173095 std:
0.16819931769007415
ElasticNet - mean: 0.7988515079441996 std: 0.09774302215729834
GradientBoostingRegressor - mean: 0.6875606157113585 std:
0.23012268956944482
XGBRegressor - mean: 0.7207995095977946 std: 0.19388589251272895
Lasso - mean: 0.7988493286692742 std: 0.09774221635911819
Ridge - mean: 0.7988493084839189 std: 0.09774189636351902

lasso = Lasso()

lasso.fit(X_train, y_train)

# points are randomly dispersed around the horizontal axis
# a linear model is our best route
from yellowbrick.regressor import ResidualsPlot
visualizer = ResidualsPlot(lasso, hist=False, qqplot=True)
visualizer.fit(X_train, y_train)  # Fit the training data
visualizer.score(X_test, y_test)  # Evaluate the model for variance
visualizer.show()

<matplotlib.axes._subplots.AxesSubplot at 0x1449bbeb8>

from sklearn.metrics import mean_squared_error


mse = []
for alpha in [0, 0.0001, 0.001, 0.01, 0.1, 1]:
  model = Lasso(alpha=alpha)
  model.fit(X_train, y_train)
  mse.append(mean_squared_error(y_test, model.predict(X_test)))
# lowest mean squared error is when alpha is 1
plt.plot([0, 0.0001, 0.001, 0.01, 0.1, 1], mse)

[<matplotlib.lines.Line2D at 0x1365faba8>]

from sklearn.model_selection import GridSearchCV
# use grid search to tune our alpha
param_grids = {'alpha': [0, 0.0001, 0.001, 0.01, 0.1, 1]}

# setup grid search
grid_search_lasso = GridSearchCV(estimator=lasso, param_grid=param_grids, cv=cv)

grid_search_lasso.fit(X_train, y_train)

grid_search_lasso.best_score_

/Users/pythagoras/anaconda3/lib/python3.7/site-
packages/sklearn/model_selection/_search.py:813: DeprecationWarning:
The default of the `iid` parameter will change from True to False in
version 0.22 and will be removed in 0.24. This will change numeric
results when test-set sizes are unequal.
  DeprecationWarning)

0.8001171769697994

# best alpha is indeed 1
grid_search_lasso.best_params_

{'alpha': 1}

from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline

# make pipeline
pipe = make_pipeline(StandardScaler(), GridSearchCV(estimator=lasso, param_grid=param_grids, cv=cv))
pipe.fit(X_train, y_train)

# As you can see our model barely improved
pipe.score(X_test, y_test)

/Users/pythagoras/anaconda3/lib/python3.7/site-
packages/sklearn/model_selection/_search.py:813: DeprecationWarning:
The default of the `iid` parameter will change from True to False in
version 0.22 and will be removed in 0.24. This will change numeric
results when test-set sizes are unequal.
  DeprecationWarning)

0.8053618445370625

from sklearn.externals import joblib

# save model for later use
joblib.dump(lasso, 'lasso_grid_search.pkl')

/Users/pythagoras/anaconda3/lib/python3.7/site-
packages/sklearn/externals/joblib/__init__.py:15: DeprecationWarning:
sklearn.externals.joblib is deprecated in 0.21 and will be removed in
0.23. Please import this functionality directly from joblib, which can
be installed with: pip install joblib. If this warning is raised when
loading pickled models, you may need to re-serialize those models with
scikit-learn 0.21+.
  warnings.warn(msg, category=DeprecationWarning)

['lasso_grid_search.pkl']