fldm2_book.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FractureProof\n",
    "An Artificial Intelligence process for feature selection for multiple geographic layers wiothout the need of human selection for the purpose of informative policy analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The 'Mr FractureProofs Contemplative Woodcarvings' Process\n",
    "This FP option uses the following sequence:\n",
    "\n",
    "1. The raw data set is cleaned of observation labels, the outcome is defined, and a missing value methodology is applied to create a cohort for feature selection.  \n",
    "2. The absolue value of the eigenevctors for each feature in the cohort from components that explain significant variance in a PCA model  is collected. \n",
    "3. Gini impurity measures for each feature in the cohort from a RF classification model are collected. \n",
    "4. The features with below average eigenvectors and gini impurity are dropped from the cohort. \n",
    "5. RFE with cross-validation is used to identify the final list of features in the cohort. \n",
    "6. Variables are placed into multiple regression model with selected confounding variables\n",
    "7. GWR identifies weighted coefficients for the selected features for each 1st layer observation\n",
    "8. 1st layer observations are averaged by 2nd layer location boundaires. \n",
    "9. Each 2nd layer location is labeled with a categorical target based on the 1st layer feature with the highest coefficient.\n",
    "10. SVM are used to identify the 2nd layer feature with the highest average cofficients for each category.\n",
    "11. 2nd layer features are selected and placed into a second multiple regression model along side 1st layer features. \n",
    "12. Raw 1st and 2nd layer data is joined, processed, and test-train 50-50 split\n",
    "13. MLPs are run with 50-500 epochs based on loss and accuracy measures during training for each of the feature sets. \n",
    "14. C-statistics are calculated from the ROCs for comparison of accuracy in identifying true positives."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FractureProof Case Study: Diabetes in Florida\n",
    "The following notebook utilizes the Mr Fracture proof's Contemplative Woodcarvings option to conduct research on social and infastructural factors related to Diabetes mortality in Florida."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### *Finding Equity: Utilizing Artificial Intelligence to Identify Social and Infrastructural Predictors of Diabetes Mortality in Florida*\n",
    "Andrew S. Cistola, MPH\n",
    "<br>Department of Health Services Research, Management, and Policy\n",
    "<br>University of Florida\n",
    "![](_fig/fldm2_map.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Introduction**: DM outcomes represent one of the largest avoidable cost burdens with opportunity for improvement in the U.S. health care system. Improving health equity in the context of DM will require targeted community improvements, infrastructure investments, and policy interventions that are designed to maximize the impact of resource allocation through the use of available data and computational resources. \n",
    "\n",
    "**Methods**: By using an Artificial Intelligence approach to evaluate over 2000 socio-economic and infrastructural predictors of DM mortality, this study used a specific series of modeling techniques to identify significant predictors without human selection and compare their predictive ability with all possible factors when passed through artificial neural networks. \n",
    "\n",
    "**Results**: The final regression model using zip code and county level predictors had an R2 of 0.863. Significant predictors included: Population % White, Population % Householders, Population % Spanish spoken at home, Population % Divorced males, Population % With public health insurance coverage, Population % Employed with private health insurance coverage, Manufacturing-Dependent Designation, Low Education Designation, Population % Medicare Part A & B Female Beneficiaries, Number of Short Term General Hospitals with 50-99 Beds. Using a multi-layered perceptron to predict zip codes at risk the C-statistic for all 2000+ predictors was 0.7938 while the 13 selected predictors was 0.8232. \n",
    "\n",
    "**Discussion**: This indicates that these factors are highly relevant for DM mortality in Florida. This process was completed without the need of human variable selection and indicates how AI can be used for informative precision public health analyses for targeted population health management efforts. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prelimninary Step: Setup Workspace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Import python libraries\n",
    "import os # Operating system navigation\n",
    "from datetime import datetime\n",
    "from datetime import date\n",
    "\n",
    "### Import data science libraries\n",
    "import pandas as pd # Widely used data manipulation library with R/Excel like tables named 'data frames'\n",
    "import numpy as np # Widely used matrix library for numerical processes\n",
    "\n",
    "### Import statistics libraries\n",
    "import scipy.stats as st # Statistics package best for t-test, ChiSq, correlation\n",
    "import statsmodels.api as sm # Statistics package best for regression models\n",
    "\n",
    "### Import Visualization Libraries\n",
    "import matplotlib.pyplot as plt # Comprehensive graphing package in python\n",
    "import geopandas as gp # Simple mapping library for csv shape files with pandas like syntax for creating plots using matplotlib \n",
    "\n",
    "### Import scikit-learn libraries\n",
    "from sklearn.preprocessing import StandardScaler # Standard scaling for easier use of machine learning algorithms\n",
    "from sklearn.impute import SimpleImputer # Univariate imputation for missing data\n",
    "from sklearn.decomposition import PCA # Principal compnents analysis from sklearn\n",
    "from sklearn.ensemble import RandomForestRegressor # Random Forest regression component\n",
    "from sklearn.feature_selection import RFECV # Recursive Feature elimination with cross validation\n",
    "from sklearn.svm import LinearSVC # Linear Support Vector Classification from sklearn\n",
    "from sklearn.linear_model import LinearRegression # Used for machine learning with quantitative outcome\n",
    "from sklearn.linear_model import LogisticRegression # Used for machine learning with quantitative outcome\n",
    "from sklearn.model_selection import train_test_split # train test split function for validation\n",
    "from sklearn.metrics import roc_curve # Reciever operator curve\n",
    "from sklearn.metrics import auc # Area under the curve \n",
    "\n",
    "### Import PySAL Libraries\n",
    "import libpysal as ps # Spatial data science modeling tools in python\n",
    "from mgwr.gwr import GWR, MGWR # Geographic weighted regression modeling tools\n",
    "from mgwr.sel_bw import Sel_BW # Bandwidth selection for GWR\n",
    "\n",
    "### Import keras libraries\n",
    "from keras.models import Sequential # Uses a simple method for building layers in MLPs\n",
    "from keras.models import Model # Uses a more complex method for building layers in deeper networks\n",
    "from keras.layers import Dense # Used for creating dense fully connected layers\n",
    "from keras.layers import Input # Used for designating input layers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Set Directory\n",
    "os.chdir(directory) # Set wd to project repository"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 1: Data Processing of Predictors and Outcomes\n",
    "Dataset 1: Florida Deaprtment of Health Vital Statistics 113 Leading Mortality Causes 2014-2018 Zip Code 5-year Average\n",
    "Dataset 2: US Census American Community Survey 2014-2018 Zip Code 5-year Average"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Feature Processing\n",
    "Before features are input into the above models, features are selected if they have over 75% non-NA values. The remaining NA values are imputed with the median values for each feature. Users can also eliminate observations that do not fit a certain critera (ex. population under 100). The resulting values are then standard scaled. All observations missing a target value are dropped and descriptive statistics are calculated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Preprocess First Dataset\n",
    "df_d1 = pd.read_csv('fracture-proof/version_2/_data/FDOH_5Y2018_ZCTA.csv') # Import first dataset saved as csv in _data folder\n",
    "df_d1 = df_d1[df_d1['POPULATION'] > 500] # Susbet numeric column by condition\n",
    "df_d1 = df_d1.filter(['K00_K99_R1000', 'ZCTA']) # Drop or filter columns to keep only feature values and idenitifer\n",
    "df_d1 = df_d1.rename(columns = {'ZCTA': 'ID', 'K00_K99_R1000': 'quant'}) # Apply standard name to identifier and quantitative outcome\n",
    "df_d1.info() # Get class, memory, and column info: names, data types, obs\n",
    "\n",
    "### Preprocess Second Dataset\n",
    "df_d2 = pd.read_csv('fracture-proof/version_2/_data/ACS_5Y2018_ZCTA.csv') # Import second dataset saved as csv in _data folder\n",
    "df_d2 = df_d2.drop(columns = ['ST', 'FIPS']) # Drop or filter columns to keep only feature values and idenitifer\n",
    "df_d2 = df_d2.select_dtypes(exclude = ['int64']) # Drop all unwanted data types\n",
    "df_d2 = df_d2.rename(columns = {'ZCTA': 'ID'}) # Apply standard name to identifier used for joining datasets\n",
    "df_d2.info() # Get class, memory, and column info: names, data types, obs.\n",
    "\n",
    "### Join Datasets by ID and define targets\n",
    "df_XY = pd.merge(df_d1, df_d2, on = 'ID', how = 'inner') # Join datasets to create table with predictors and outcome\n",
    "df_XY = df_XY.dropna(subset = ['quant']) # Drop all outcome rows with NA values\n",
    "df_XY.info() # Get class, memory, and column info: names, data types, obs.\n",
    "\n",
    "### Create outcome table\n",
    "df_Y = df_XY.filter(['quant', 'ID']) # Create Outcome table\n",
    "df_Y = df_Y.set_index('ID') # Set identifier as index\n",
    "df_Y.info() # Get class, memory, and column info: names, data types, obs.\n",
    "\n",
    "### Create standard scaled predictor table\n",
    "df_X = df_XY.drop(columns = ['quant', 'ID']) # Drop Unwanted Columns\n",
    "df_X = df_X.replace([np.inf, -np.inf], np.nan) # Replace infitite values with NA\n",
    "df_X = df_X.dropna(axis = 1, thresh = 0.75*len(df_X)) # Drop features less than 75% non-NA count for all columns\n",
    "df_X = pd.DataFrame(SimpleImputer(strategy = 'median').fit_transform(df_X), columns = df_X.columns) # Impute missing data\n",
    "df_X = pd.DataFrame(StandardScaler().fit_transform(df_X.values), columns = df_X.columns) # Standard scale values by converting the normalized features into a tabular format with the help of DataFrame.\n",
    "df_X['ID'] = df_XY['ID'] # Save ID as column in predictor table\n",
    "df_X = df_X.set_index('ID') # Set identifier as index\n",
    "df_X.info() # Get class, memory, and column info: names, data types, obs.\n",
    "\n",
    "### Add feature labels\n",
    "df_l1 = pd.read_csv('fracture-proof/version_2/_data/ACS_5Y2018_labels.csv') # Import feature lables for first dataset saved as csv in _data folder\n",
    "df_l2 = pd.read_csv('fracture-proof/version_2/_data/FDOH_5Y2018_labels.csv') # Import feature lables for second dataset saved as csv in _data folder\n",
    "df_l1_l2 = pd.concat([df_l1, df_l2]) # Combine rows with same columns\n",
    "df_l1_l2 = df_l1_l2.filter(['Feature', 'Label']) # Keep only selected columns\n",
    "df_l1_l2 = df_l1_l2.set_index('Feature') # Set column as index\n",
    "df_l1_l2 = df_l1_l2.transpose() # Switch rows and columns\n",
    "df_l1_l2.info # Get class, memory, and column info: names, data types, obs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2: Identify Predictors with Open Box Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Principal Component Analysis (PCA)\n",
    "Using linear transformations on a normalized covariace matrix, PCA creates combinations of features in a regression model into principal components that explain a proportion of the variance observed in the dataset. The coefficients of the principal component models (eigenvectors) that explain a significant amount of variance can be compared and variables that do not explain significant variance can be dropped. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Principal Component Analysis\n",
    "degree = len(df_X.columns) - 1  # Save number of features -1 to get degrees of freedom\n",
    "pca = PCA(n_components = degree) # Pass the number of components to make PCA model based on degrees of freedom\n",
    "pca.fit(df_X) # Fit initial PCA model\n",
    "df_comp = pd.DataFrame(pca.explained_variance_) # Print explained variance of components\n",
    "df_comp = df_comp[(df_comp[0] > 1)] # Save eigenvalues above 1 to identify components\n",
    "components = len(df_comp.index) - 1 # Save count of components for Variable reduction\n",
    "pca = PCA(n_components = components) # you will pass the number of components to make PCA model\n",
    "pca.fit_transform(df_X) # finally call fit_transform on the aggregate data to create PCA results object\n",
    "df_pc = pd.DataFrame(pca.components_, columns = df_X.columns) # Export eigenvectors to data frame with column names from original data\n",
    "df_pc[\"Variance\"] = pca.explained_variance_ratio_ # Save eigenvalues as their own column\n",
    "df_pc = df_pc[df_pc[\"Variance\"] > df_pc[\"Variance\"].mean()] # Susbet by eigenvalues with above average exlained variance ratio\n",
    "df_pc = df_pc.abs() # Get absolute value of eigenvalues\n",
    "df_pc = df_pc.drop(columns = [\"Variance\"]) # Drop outcomes and targets\n",
    "df_p = pd.DataFrame(df_pc.max(), columns = [\"MaxEV\"]) # select maximum eigenvector for each feature\n",
    "df_p = df_p[df_p.MaxEV > df_p.MaxEV.mean()] # Susbet by above average max eigenvalues \n",
    "df_p = df_p.reset_index() # Add a new index of ascending values, existing index consisting of feature labels becomes column named \"index\"\n",
    "df_pca = df_p.rename(columns = {\"index\": \"Feature\"}) # Rename former index as features\n",
    "df_pca = df_pca.sort_values(by = [\"MaxEV\"], ascending = False) # Sort Columns by Value\n",
    "df_pca.info() # Get class, memory, and column info: names, data types, obs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Random Forests (RF)\n",
    "By aggregating decision trees from a bootstrapped sample of features, random forests can measure the importance of a given feature in predicting the outcome of interest. By calculating the change in prediction capability when the feature is removed, the importance value (Gini Impurity) of a feature can be compared to others. Features that are not important compared to the others can be dropped. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Random Forest Regressor\n",
    "forest = RandomForestRegressor(n_estimators = 1000, max_depth = 10) #Use default values except for number of trees. For a further explanation see readme included in repository. \n",
    "forest.fit(df_X, df_Y['quant']) # Fit Forest model, This will take time\n",
    "rf = forest.feature_importances_ # Output importances of features\n",
    "l_rf = list(zip(df_X, rf)) # Create list of variables alongside importance scores \n",
    "df_rf = pd.DataFrame(l_rf, columns = ['Feature', 'Gini']) # Create data frame of importances with variables and gini column names\n",
    "df_rf = df_rf[(df_rf['Gini'] > df_rf['Gini'].mean())] # Subset by Gini values higher than mean\n",
    "df_rf = df_rf.sort_values(by = ['Gini'], ascending = False) # Sort Columns by Value\n",
    "df_rf.info() # Get class, memory, and column info: names, data types, obs.\n",
    "\n",
    "### Fracture: Join RF and PCA \n",
    "df_fr = pd.merge(df_pca, df_rf, on = 'Feature', how = 'inner') # Join by column while keeping only items that exist in both, select outer or left for other options\n",
    "fracture = df_fr['Feature'].tolist() # Save features from data frame\n",
    "df_fr.info() # Get class, memory, and column info: names, data types, obs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Recursive Feature Elimination (RFE)\n",
    "Using regression methods, RFE creates a predictive model and removes weakest features by comparing prediction scores. RFE can cross-validate to remove cosnistently weaker features until a ideal feature set is defined. A minimum or predefined set of features can be used. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Recursive Feature Elimination\n",
    "recursive = RFECV(estimator = LinearRegression(), min_features_to_select = 5) # define selection parameters, in this case all features are selected. See Readme for more ifo\n",
    "recursive.fit(df_X[fracture], df_Y['quant']) # This will take time\n",
    "rfe = recursive.support_ # Save Boolean values as numpy array\n",
    "l_rfe = list(zip(df_X[fracture], rfe)) # Create list of variables alongside RFE value \n",
    "df_rfe = pd.DataFrame(l_rfe, columns = ['Feature', 'RFE']) # Create data frame of importances with variables and gini column names\n",
    "df_rfe = df_rfe.sort_values(by = ['RFE'], ascending = True) # Sort Columns by Value\n",
    "df_rfe = df_rfe[df_rfe['RFE'] == True] # Select Variables that were True\n",
    "df_rfe.info() # Get class, memory, and column info: names, data types, obs.\n",
    "\n",
    "### FractureProof: Join RFE with Fracture\n",
    "df_fp = pd.merge(df_fr, df_rfe, on = 'Feature', how = 'inner') # Join by column while keeping only items that exist in both, select outer or left for other options\n",
    "fractureproof = df_fp['Feature'].tolist() # Save chosen featres as list\n",
    "df_fp.info() # Get class, memory, and column info: names, data types, obs.\n",
    "\n",
    "### Get FractureProof feature labels\n",
    "df_lfp = df_l1_l2[fractureproof] # Save chosen featres as list\n",
    "df_lfp = df_lfp.transpose() # Switch rows and columns\n",
    "df_lfp = df_lfp.reset_index() # Reset index\n",
    "l_lfp = list(zip(df_lfp['Feature'], df_lfp['Label'])) # Create list of variables alongside RFE value \n",
    "df_lfp.info() # Get class, memory, and column info: names, data types, obs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 3: Create Informative Prediction Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Multiple regression modeling (MR)\n",
    "Using a multiple linear regression model, confounders are added and the selected features can be evaluated using R-squared, F-statistic and significance values. Features can also be compared for magnitude and direction. This model is created using the raw data before imputation and standard scaling. All observations with missing data are dropped. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Add confounders to multiple regression model\n",
    "mrfractureproof = df_X[fractureproof].columns.to_list() # Save list of selected variables for multiple regression model\n",
    "mrfractureproof.append('quant') # Add outcome to list of selected variables for multiple regression model\n",
    "mrfractureproof.append('DP05_0024PE') # Add confounder (Over 65) to list of selected variables for multiple regression model\n",
    "\n",
    "### Create Multiple Regression Model\n",
    "df_mrfp = df_XY[mrfractureproof] # Subset original nonscaled data for regression\n",
    "df_mrfp = df_mrfp.dropna() # Drop all rows with NA values\n",
    "X = df_mrfp.drop(columns = ['quant']) # Susbet predictors for regression\n",
    "Y = df_mrfp['quant'] # Subset quantitative outcome for regression\n",
    "mod = sm.OLS(Y, X) # Create linear regression model\n",
    "res = mod.fit() # Fit model to create result\n",
    "res.summary() # Print results of regression model\n",
    "\n",
    "### Add feature labels\n",
    "mrfractureproof.remove('quant') # Remove outcome to list of features used for collecting lables\n",
    "df_lmrfp = df_l1_l2[mrfractureproof] # Save selected features as list for collecting labels\n",
    "mrfractureproof.append('quant') # Add outcome to to list of selected variables for multiple regression model\n",
    "df_lmrfp = df_lmrfp.transpose() # Switch rows and columns\n",
    "df_lmrfp = df_lmrfp.reset_index() # Reset index\n",
    "l_lmrfp = list(zip(df_lmrfp['Feature'], df_lmrfp['Label'])) # Create list of variables alongside RFE value \n",
    "df_lmrfp.info() # Get class, memory, and column info: names, data types, obs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 4: Geographic Weighted Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Geographic Weighted Regression (GWR)\n",
    "Using GIS data for for observations from the first layer, regression models are caluclated from existing features with coefficients weighted based on location. Locations where the given features are highest are labeled for the purpose of identifying 2nd layer locations where sleected features have higher predictive weight. This process utilizes the PySal library (https://pysal.org)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Geojoin predictor and outcome table with polygons, Get centroid from coordinates\n",
    "gdf_XY = gp.read_file('fracture-proof/version_2/_data/cb_2018_us_zcta510_500k/cb_2018_us_zcta510_500k.shp') # Import shape files from folder with all other files downloaded\n",
    "gdf_XY['ID'] = gdf_XY['ZCTA5CE10'].astype('str') # Change data type of column in data frame\n",
    "gdf_XY['ID'] = gdf_XY['ID'].str.rjust(5, '0') # add leading zeros of character column using rjust() function\n",
    "gdf_XY['ID'] = 'ZCTA' + gdf_XY['ID'] # Combine string with column\n",
    "gdf_XY = gdf_XY.filter(['ID', 'geometry']) # Keep only selected columns\n",
    "gdf_XY = pd.merge(gdf_XY, df_XY, on = 'ID', how = 'inner') # Geojoins can use pandas merge as long as geo data is first passed in function\n",
    "gdf_XY['x'] = gdf_XY['geometry'].centroid.x # Save centroid coordinates as separate column\n",
    "gdf_XY['y'] = gdf_XY['geometry'].centroid.y # Save centroid coordinates as separate column\n",
    "gdf_XY['coordinates'] = list(zip(gdf_XY['x'], gdf_XY['y'])) # Save individual coordinates as column of paired list\n",
    "gdf_XY = gdf_XY.drop(columns = ['x', 'y', 'geometry']) # Drop Unwanted Columns\n",
    "gdf_XY.info() # Get class, memory, and column info: names, data types, obs.\n",
    "\n",
    "### Setup GWR table\n",
    "gdf_gwr = gdf_XY.set_index('ID') # Set ID column as index\n",
    "wood = gdf_gwr[fractureproof].columns.to_list() # Save fractureproof variables as list for GWR\n",
    "wood.append('quant') # Add outcome to list of gwr variables\n",
    "wood.append('coordinates') # Add coordinates to list of gwr variables\n",
    "gdf_gwr = gdf_gwr[wood] # Subset dataframe by sleetced variables\n",
    "gdf_gwr = gdf_gwr.dropna() # Drop all rows with NA values\n",
    "c = list(gdf_gwr[\"coordinates\"]) # save coordinates column as list\n",
    "x = gdf_gwr.drop(columns = ['quant', 'coordinates']).values # save selected features as numpy array\n",
    "y = gdf_gwr['quant'].values # save target as numpy array\n",
    "y = np.transpose([y]) # Transpose numpy array to fit GWR input\n",
    "gdf_gwr.info() # Get class, memory, and column info: names, data types, obs.\n",
    "\n",
    "### Create GWR model\n",
    "mgwr_selector = Sel_BW(c, y, x, multi = True) # create model to calibrate selector\n",
    "mgwr_bw = mgwr_selector.search(multi_bw_min = [2]) # search for selectors with minimum of 2 bandwidths, this may take a while\n",
    "mgwr_results = MGWR(c, y, x, mgwr_selector).fit() # fit MGWR model, this may take a while\n",
    "mgwr_results.summary() # Show MGWR summary\n",
    "\n",
    "### Export GWR results to new table\n",
    "wood.remove('quant') # Remove outcome to list of gwr variables\n",
    "wood.remove('coordinates') # Remove coordinates to list of gwr variables\n",
    "wood = ['Intercept'] + wood # Insert intercept label at front of gwr variable list\n",
    "df_gwr = pd.DataFrame(mgwr_results.params, columns = [wood]) # Create data frame of importances with variables and gini column names\n",
    "gdf_ID = gdf_gwr.reset_index() # Reset index on GWR inputs\n",
    "df_gwr['ID'] = gdf_ID['ID'] # Ad ID column from GWR inputs table\n",
    "df_gwr.info()  # Get class, memory, and column info: names, data types, obs.\n",
    "\n",
    "### Join first and second geographic layer lables\n",
    "df_layer = pd.read_csv('fracture-proof/version_2/_data/FIPS_ZCTA_key.csv') # Import layer key dataset saved as csv in _data folder\n",
    "df_layer = df_layer.filter(['FIPS', 'ZCTA']) # Keep only selected columns\n",
    "df_layer = df_layer.rename(columns = {'ZCTA': 'ID', 'FIPS': 'ID_2'}) # Rename geographic identifiers as standard features\n",
    "gdf_ID_2 = gdf_gwr.reset_index() # Reset Index as second geographic layer ID adn save as gdf for later\n",
    "df_gwr = pd.merge(gdf_gwr, df_layer, on = 'ID', how = 'left') # Join zip code geo weighted coefficients to county labels\n",
    "df_gwr = df_gwr.dropna() # Drop all rows with NA values\n",
    "df_gwr = df_gwr.set_index('ID') # Set first layer ID column as index\n",
    "df_gwr = df_gwr.drop(columns = ['coordinates', 'quant']) # Drop Unwanted Columns\n",
    "df_gwr = df_gwr.groupby(['ID_2'], as_index = False).mean() # Group 1st layer GWR coefficents by 2nd layer identifiers and calculate average\n",
    "df_gwr.info() # Get class, memory, and column info: names, data types, obs.\n",
    "\n",
    "### Create Multi-level categories based bandwidths\n",
    "df_bw = df_gwr.drop(columns = ['ID_2']) # Drop Unwanted Columns\n",
    "df_bw = df_bw.apply(st.zscore).abs() # Calculate absolute value of z-score for mean GWR coefficients in second layer\n",
    "df_bw['ID_2'] = df_gwr['ID_2'] # Save second layer identifiers from GWR dataset\n",
    "df_bw = df_bw.set_index('ID_2') # Set second layer identifiers as index\n",
    "bw = df_bw.idxmax(axis = 1) # Get second layer identifiers that have highest absolute value of z score\n",
    "l_bw = list(zip(df_bw.index, bw)) # Create list of variables alongside RFE value \n",
    "df_bw = pd.DataFrame(l_bw, columns = ['ID_2', 'multi']) # Create data frame of 1st layer features and 2nd layer identifiers\n",
    "df_bw['multi'] = df_bw['multi'].astype('category') # Save features as multi-level categoriacl variable with standard name\n",
    "df_bw['multi'] = df_bw['multi'].cat.codes # Convert string lable into numeric codes\n",
    "df_bw.info() # Get class, memory, and column info: names, data types, obs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 5: Data Processing of 2nd Geographic Layer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Feature Processing\n",
    "Before features are input into the above models, features are selected if they have over 75% non-NA values. The remaining NA values are imputed with the median values for each feature. Users can also eliminate observations that do not fit a certain critera (ex. population under 100). The resulting values are then standard scaled. All observations missing a target value are dropped and descriptive statistics are calculated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Preprocess Third Dataset\n",
    "df_d3 = pd.read_csv('fracture-proof/version_2/_data/AHRF_5Y2018_FIPS.csv') # Import third dataset saved as csv in _data folder\n",
    "df_d3 = df_d3.rename(columns = {'FIPS': 'ID_2'}) # Apply standard name to identifier used for joining 2nd layer datasets\n",
    "df_d3.info() # Get class, memory, and column info: names, data types, obs\n",
    "\n",
    "### Join Datasets by second layer identifier and define targets\n",
    "df_XY_2 = pd.merge(df_d3, df_bw, on = 'ID_2', how = 'inner') # Join datasets to create table with predictors and outcome\n",
    "df_XY_2.info() # Get class, memory, and column info: names, data types, obs.\n",
    "\n",
    "### Create second layer outcome table\n",
    "df_Y_2 = df_XY_2.filter(['multi', 'ID_2']) # Create Outcome table for second layer\n",
    "df_Y_2 = df_Y_2.set_index('ID_2') # Set second layer identifier as index\n",
    "df_Y_2.info() # Get class, memory, and column info: names, data types, obs.\n",
    "\n",
    "### Create second layer standard scaled predictor table\n",
    "df_X_2 = df_XY_2.drop(columns = ['multi', 'ID_2']) # Drop Unwanted Columns\n",
    "df_X_2 = df_X_2.replace([np.inf, -np.inf], np.nan) # Replace infitite values with NA\n",
    "df_X_2 = df_X_2.dropna(axis = 1, thresh = 0.75*len(df_X_2)) # Drop features less than 75% non-NA count for all columns\n",
    "df_X_2 = pd.DataFrame(SimpleImputer(strategy = 'median').fit_transform(df_X_2), columns = df_X_2.columns) # Impute missing data\n",
    "df_X_2 = pd.DataFrame(StandardScaler().fit_transform(df_X_2.values), columns = df_X_2.columns) # Standard scale values by converting the normalized features into a tabular format with the help of DataFrame.\n",
    "df_X_2['ID_2'] = df_XY_2['ID_2'] # Save ID as column in predictor table\n",
    "df_X_2 = df_X_2.set_index('ID_2') # Set identifier as index\n",
    "df_X_2.info() # Get class, memory, and column info: names, data types, obs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 6: Identify 2nd Layer Predictors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Support Vector Machines (SVM)\n",
    "Using the labeles creted by the GWR, support vector mahcines are used to identify 2nd layer features with the highest weights for the given location. Since the 2nd layer will often have signficiantly fewer obserrvations and predicts a multi-level categorical target, SVMs were slected for their ability to handle these constraints better than other available models. The result is a set of 2nd layer features that independently contribute or interact with the 1st layer features. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Support Vector Machines\n",
    "vector = LinearSVC() # Support vector machines with a linear kernel for multi-level categorical outrcomes\n",
    "vector.fit(df_X_2, df_Y_2['multi']) # fit model\n",
    "svm = vector.coef_ # Save coefficients for each category by feature\n",
    "df_svm = pd.DataFrame(svm, columns = df_X_2.columns, index = [fractureproof]) # Create data frame of coefficients by 2nd layer features and 1st layer features\n",
    "df_svm = df_svm.abs() # Get absolute value of all coefficients\n",
    "svm_max = df_svm.idxmax(axis = 1) # Get 2nd layer features that have highest values for each 1st layer feature\n",
    "l_svm_max = list(zip(df_svm.index, svm_max)) # Create list of 2nd layer features along 1st layer features\n",
    "df_svm_max = pd.DataFrame(l_svm_max, columns = ['GWR', 'Feature']) # Create data frame of 2nd layer features along 1st layer features\n",
    "carving = df_svm_max['Feature'].unique() # Print unique values in column to remove duplicate 2nd layer features and save as list\n",
    "df_svm_max.info() # Get class, memory, and column info: names, data types, obs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Principal Component Analysis (PCA)\n",
    "Using linear transformations on a normalized covariace matrix, PCA creates combinations of features in a regression model into principal components that explain a proportion of the variance observed in the dataset. PCA model is created and component loadings are used to eliminbate features accounting for less than 95% of variation in the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Principal Component Analysis\n",
    "degree = len(df_X_2[carving].columns) - 1  # Save number of features -1 to get degrees of freedom\n",
    "pca = PCA(n_components = degree) # Pass the number of components to make PCA model based on degrees of freedom\n",
    "pca.fit(df_X_2[carving]) # Fit initial PCA model\n",
    "\n",
    "### Variance ratios and component Loadings\n",
    "cvr = pca.explained_variance_ratio_.cumsum() # Save cumulative variance ratio\n",
    "comps = np.count_nonzero(cvr) - np.count_nonzero(cvr > 0.95) + 1 # Save number of components above threshold value\n",
    "load = pca.components_.T * np.sqrt(pca.explained_variance_) # Export component loadings\n",
    "df_load = pd.DataFrame(load, index = df_X_2[carving].columns) # Create data frame of component loading\n",
    "df_load = df_load.iloc[:, 0:comps] # Save columns by components above threshold\n",
    "df_load = df_load.abs() # get absolute value for column or data frame\n",
    "df_load = df_load[df_load > 0.5] # Subset by character\n",
    "df_load = df_load.dropna(thresh = 1) # Drop all rows without 1 non-NA value\n",
    "df_load = df_load.dropna(axis = 'columns', thresh = 1) # Drop all rows without 1 non-NA value\n",
    "woodcarving = df_load.index.to_list() # Save final set of 2nd layer features to list\n",
    "df_load.info() # Get class, memory, and column info: names, data types, obs.\n",
    "\n",
    "### Add 2nd layer feature labels\n",
    "df_l3 = pd.read_csv('fracture-proof/version_2/_data/AHRF_5Y2018_labels.csv') # Import dataset saved as csv in _data folder\n",
    "df_l_3 = df_l3.filter(['Feature', 'Label']) # Keep only selected columns\n",
    "df_l_3 = df_l_3.set_index('Feature') # Set column as index\n",
    "df_l_3 = df_l_3.transpose() # Switch rows and columns\n",
    "df_lwc = df_l_3[woodcarving] # Subset by 2nd layer selected featres\n",
    "df_lwc = df_lwc.transpose() # Switch rows and columns\n",
    "df_lwc = df_lwc.reset_index() # Reset index\n",
    "l_lwc = list(zip(df_lwc['Feature'], df_lwc['Label'])) # Create list of variables alongside RFE value \n",
    "df_lwc.info() # Get class, memory, and column info: names, data types, obs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 7: Create Informative Prediction Model with both geographic layers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Multiple regression modeling (MR)\n",
    "Using a multiple linear regression model, confounders are added and the selected features can be evaluated using R-squared, F-statistic and significance values. Features can also be compared for magnitude and direction. This model is created using the raw data before imputation and standard scaling. All observations with missing data are dropped. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Join Datasets by ID and define targets\n",
    "df_XY_f = pd.merge(df_XY_2, df_layer, on = 'ID_2', how = 'left') # Join datasets to create table with predictors and outcome\n",
    "df_XY_f = pd.merge(df_XY, df_XY_f, on = 'ID', how = 'inner') # Join datasets to create table with predictors and outcome\n",
    "df_XY_f = df_XY_f.drop(columns = ['ID_2', 'multi']) # Drop Unwanted Columns\n",
    "df_XY_f = df_XY_f.dropna(subset = ['quant']) # Drop all outcome rows with NA values\n",
    "df_XY_f.info() # Get class, memory, and column info: names, data types, obs.\n",
    "\n",
    "### Create Multiple Regression Model\n",
    "mrfractureproofswoodcarvings = mrfractureproof + woodcarving # Combine 2nd layer slected features with first layer regression model features\n",
    "mrfractureproofswoodcarvings.append('quant') # Remove outcome from 1st and 2nd layer regression feature list\n",
    "df_mrfpwc = df_XY_f[mrfractureproofswoodcarvings] # Subset full dataframe with 1st and 2nd layer regression model features\n",
    "df_mrfpwc = df_mrfpwc.dropna() # Drop all NA values from subset dataframe\n",
    "X = df_mrfpwc.drop(columns = ['quant']) # Create dataframe of predictors\n",
    "Y = df_mrfpwc['quant'] # Create dataframe of outcomes\n",
    "mod_f = sm.OLS(Y, X) # Create linear model\n",
    "res_f = mod_f.fit() # Fit model to create result\n",
    "res_f.summary() # Print results of regression model\n",
    "\n",
    "### Add feature labels\n",
    "df_lf = pd.concat([df_l1, df_l2, df_l3]) # Combine rows with same columns\n",
    "df_lf = df_lf.filter(['Feature', 'Label']) # Keep only selected columns\n",
    "df_lf = df_lf.set_index('Feature') # Set column as index\n",
    "df_lf = df_lf.transpose() # Switch rows and columns\n",
    "mrfractureproofswoodcarvings.remove('quant') # Remove outcome from 1st and 2nd layer regression feature list\n",
    "df_lmfpwc = df_lf[mrfractureproofswoodcarvings] # Save chosen featres as list\n",
    "df_lmfpwc = df_lmfpwc.transpose() # Switch rows and columns\n",
    "df_lmfpwc = df_lmfpwc.reset_index() # Reset index\n",
    "l_lmfpwc = list(zip(df_lmfpwc['Feature'], df_lmfpwc['Label'])) # Create list of variables alongside RFE value \n",
    "df_lmfpwc.info() # Get class, memory, and column info: names, data types, obs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 8: Predict Binary Outcome with Artificial Neural Networks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Feature Processing\n",
    "The raw data for the 1st and 2nd layers are processed using the steps in FractureProcess. The table is then randomly split 50-50 into test and train tables for evalatuing target prediction from the MLPs. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Create outcome table and define targets\n",
    "df_Y_f = df_XY_f.filter(['quant', 'ID']) # Create Outcome table\n",
    "df_Y_f['binary'] = np.where(df_Y_f['quant'] > df_Y_f['quant'].quantile(0.5), 1, 0) # Create binary outcome based on conditions\n",
    "df_Y_f = df_Y_f.set_index('ID') # Set identifier as index\n",
    "df_Y_f.info() # Get class, memory, and column info: names, data types, obs.\n",
    "\n",
    "### Create standard scaled predictor table\n",
    "df_X_f = df_XY_f.drop(columns = ['quant', 'ID']) # Drop Unwanted Columns\n",
    "df_X_f = df_X_f.replace([np.inf, -np.inf], np.nan) # Replace infitite values with NA\n",
    "df_X_f = df_X_f.dropna(axis = 1, thresh = 0.75*len(df_X_f)) # Drop features less than 75% non-NA count for all columns\n",
    "df_X_f = pd.DataFrame(SimpleImputer(strategy = 'median').fit_transform(df_X_f), columns = df_X_f.columns) # Impute missing data\n",
    "df_X_f = pd.DataFrame(StandardScaler().fit_transform(df_X_f.values), columns = df_X_f.columns) # Standard scale values by converting the normalized features into a tabular format with the help of DataFrame.\n",
    "df_X_f['ID'] = df_XY_f['ID'] # Save ID as column in predictor table\n",
    "df_X_f = df_X_f.set_index('ID') # Set identifier as index\n",
    "df_X_f.info() # Get class, memory, and column info: names, data types, obs.\n",
    "\n",
    "### Save FractureProof and Woodcarving feature list\n",
    "mrfractureproofscontemplativewoodcarvings = mrfractureproof + woodcarving # Combine 2nd layer slected features with first layer regression model features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Multi-layered Perceptrons (MLP)\n",
    "An artificical neural network consisting of 2 dense layers, and a binary activatrion layer is used to predict a binary outcome calculated based on the original quantitative target. Predictions are made with all possible features, 1st layer features only, 2nd layer features only, and the final list of 1st and 2nd layer features. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Receiver Operator Curve (ROC)\n",
    "In order to indetify whether the selected vfeatures provide real world practicality in improving prediction, ROCs are created and C-statistics are calcuklated to determine the amount opf true positives to false positives. This allows for easy comparison of whether the selected features are relevant to decision making. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Multi-Layered Perceptron with all predictors from all layers\n",
    "Y = df_Y_f.filter(['binary']) # Save binary outcome as MLP Input\n",
    "X = df_X_f # Save all predictors as MLP input\n",
    "X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size = 0.50) # Random 50/50 Train/Test Split\n",
    "input = X.shape[1] # Save number of columns as input dimension\n",
    "nodes = round(input / 2) # Number of input dimensions divided by two for nodes in each layer\n",
    "epochs = 50\n",
    "network = Sequential() # Build Network with keras Sequential API\n",
    "network.add(Dense(nodes, activation = 'relu', kernel_initializer = 'random_normal', input_dim = input)) # First dense layer\n",
    "network.add(Dense(nodes, activation = 'relu', kernel_initializer = 'random_normal')) # Second dense layer\n",
    "network.add(Dense(1, activation = 'sigmoid', kernel_initializer = 'random_normal')) # Output layer with binary activation\n",
    "network.compile(optimizer = 'adam', loss = 'binary_crossentropy', metrics = ['accuracy']) # Compile network with Adaptive moment estimation, and follow loss and accuracy\n",
    "final = network.fit(X_train, Y_train, batch_size = 10, epochs = epochs) # Fitting the data to the train outcome, with batch size and number of epochs\n",
    "Y_pred = network.predict(X_test) # Predict values from test data\n",
    "Y_pred = (Y_pred > 0.5) # Save predicted values close to 1 as boolean\n",
    "Y_test = (Y_test > 0.5) # Save test values close to 1 as boolean\n",
    "fpr, tpr, threshold = roc_curve(Y_test, Y_pred) # Create ROC outputs, true positive rate and false positive rate\n",
    "auc_a = auc(fpr, tpr) # Plot ROC and get AUC score\n",
    "e_a = epochs # Save epochs used for mlp\n",
    "print(auc_a) # Display object\n",
    "\n",
    "### Multi-Layered Perceptron with Mr. Fracture Proof predictors\n",
    "Y = df_Y_f.filter(['binary']) # Save binary outcome as MLP Input\n",
    "X = df_X_f[mrfractureproof] # Save selected predictors from all layers predictors as MLP input\n",
    "X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size = 0.50) # Random 50/50 Train/Test Split\n",
    "input = X.shape[1] # Save number of columns as input dimension\n",
    "nodes = round(input / 2) # Number of input dimensions divided by two for nodes in each layer\n",
    "epochs = 500\n",
    "network = Sequential() # Build Network with keras Sequential API\n",
    "network.add(Dense(nodes, activation = 'relu', kernel_initializer = 'random_normal', input_dim = input)) # First dense layer\n",
    "network.add(Dense(nodes, activation = 'relu', kernel_initializer = 'random_normal')) # Second dense layer\n",
    "network.add(Dense(1, activation = 'sigmoid', kernel_initializer = 'random_normal')) # Output layer with binary activation\n",
    "network.compile(optimizer = 'adam', loss = 'binary_crossentropy', metrics = ['accuracy']) # Compile network with Adaptive moment estimation, and follow loss and accuracy\n",
    "final = network.fit(X_train, Y_train, batch_size = 10, epochs = epochs) # Fitting the data to the train outcome, with batch size and number of epochs\n",
    "Y_pred = network.predict(X_test) # Predict values from test data\n",
    "Y_pred = (Y_pred > 0.5) # Save predicted values close to 1 as boolean\n",
    "Y_test = (Y_test > 0.5) # Save test values close to 1 as boolean\n",
    "fpr, tpr, threshold = roc_curve(Y_test, Y_pred) # Create ROC outputs, true positive rate and false positive rate\n",
    "auc_mrfp = auc(fpr, tpr) # Plot ROC and get AUC score\n",
    "e_mrfp = epochs # Save epochs used for mlp\n",
    "print(auc_mrfp) # Display object\n",
    "\n",
    "### Multi-Layered Perceptron with Woodcarving predictors\n",
    "Y = df_Y_f.filter(['binary']) # Save binary outcome as MLP Input\n",
    "X = df_X_f[woodcarving] # Save selected predictors from all layers predictors as MLP input\n",
    "X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size = 0.50) # Random 50/50 Train/Test Split\n",
    "input = X.shape[1] # Save number of columns as input dimension\n",
    "nodes = round(input / 2) # Number of input dimensions divided by two for nodes in each layer\n",
    "epochs = 500\n",
    "network = Sequential() # Build Network with keras Sequential API\n",
    "network.add(Dense(nodes, activation = 'relu', kernel_initializer = 'random_normal', input_dim = input)) # First dense layer\n",
    "network.add(Dense(nodes, activation = 'relu', kernel_initializer = 'random_normal')) # Second dense layer\n",
    "network.add(Dense(1, activation = 'sigmoid', kernel_initializer = 'random_normal')) # Output layer with binary activation\n",
    "network.compile(optimizer = 'adam', loss = 'binary_crossentropy', metrics = ['accuracy']) # Compile network with Adaptive moment estimation, and follow loss and accuracy\n",
    "final = network.fit(X_train, Y_train, batch_size = 10, epochs = epochs) # Fitting the data to the train outcome, with batch size and number of epochs\n",
    "Y_pred = network.predict(X_test) # Predict values from test data\n",
    "Y_pred = (Y_pred > 0.5) # Save predicted values close to 1 as boolean\n",
    "Y_test = (Y_test > 0.5) # Save test values close to 1 as boolean\n",
    "fpr, tpr, threshold = roc_curve(Y_test, Y_pred) # Create ROC outputs, true positive rate and false positive rate\n",
    "auc_wc = auc(fpr, tpr) # Plot ROC and get AUC score\n",
    "e_wc = epochs # Save epochs used for mlp\n",
    "print(auc_wc) # Display object\n",
    "\n",
    "### Multi-Layered Perceptron with Mr. Fracture Proof's Contemplative Woodcarving predictors\n",
    "Y = df_Y_f.filter(['binary']) # Save binary outcome as MLP Input\n",
    "X = df_X_f[mrfractureproofscontemplativewoodcarvings] # Save selected predictors from all layers predictors as MLP input\n",
    "X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size = 0.50) # Random 50/50 Train/Test Split\n",
    "input = X.shape[1] # Save number of columns as input dimension\n",
    "nodes = round(input / 2) # Number of input dimensions divided by two for nodes in each layer\n",
    "epochs = 500\n",
    "network = Sequential() # Build Network with keras Sequential API\n",
    "network.add(Dense(nodes, activation = 'relu', kernel_initializer = 'random_normal', input_dim = input)) # First dense layer\n",
    "network.add(Dense(nodes, activation = 'relu', kernel_initializer = 'random_normal')) # Second dense layer\n",
    "network.add(Dense(1, activation = 'sigmoid', kernel_initializer = 'random_normal')) # Output layer with binary activation\n",
    "network.compile(optimizer = 'adam', loss = 'binary_crossentropy', metrics = ['accuracy']) # Compile network with Adaptive moment estimation, and follow loss and accuracy\n",
    "final = network.fit(X_train, Y_train, batch_size = 10, epochs = epochs) # Fitting the data to the train outcome, with batch size and number of epochs\n",
    "Y_pred = network.predict(X_test) # Predict values from test data\n",
    "Y_pred = (Y_pred > 0.5) # Save predicted values close to 1 as boolean\n",
    "Y_test = (Y_test > 0.5) # Save test values close to 1 as boolean\n",
    "fpr, tpr, threshold = roc_curve(Y_test, Y_pred) # Create ROC outputs, true positive rate and false positive rate\n",
    "auc_mrfpctwc = auc(fpr, tpr) # Plot ROC and get AUC score\n",
    "e_mrfpctwc = epochs # Save epochs used for mlp\n",
    "print(auc_mrfpctwc) # Display object"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Disclaimer\n",
    "While the author (Andrew Cistola) is a Florida DOH employee and a University of Florida PhD student, these are NOT official publications by the Florida DOH, the University of Florida, or any other agency. \n",
    "No information is included in this repository that is not available to any member of the public. \n",
    "All information in this repository is available for public review and dissemination but is not to be used for making medical decisions. \n",
    "All code and data inside this repository is available for open source use per the terms of the included license. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### allocativ\n",
    "This repository is part of the larger allocativ project dedicated to prodiving analytical tools that are 'open source for public health.' Learn more at https://allocativ.com. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Last updated 13 November 2020 by DrewC!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}