<jupyter_start><jupyter_text>**NEXT STEP: add in fdid from nfirs and merge to ACS then run this.**<jupyter_code>from pathlib import Path
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from imblearn import under_sampling, over_sampling 
from sklearn import preprocessing
from sklearn.model_selection import train_test_split
from sklearn.neural_network import MLPClassifier
from sklearn.metrics import classification_report
from sklearn.metrics import log_loss
from sklearn.metrics import plot_confusion_matrix
from sklearn.metrics import confusion_matrix
from sklearn.metrics import roc_auc_score, average_precision_score
from sklearn.model_selection import cross_val_score


%matplotlib inline
plt.style.use('dark_background')
pd.set_option('display.max_columns',500)
sns.set()

import random

SEED = 27
random.seed(SEED)
np.random.seed(SEED)

p = Path.cwd()
data_path = p.parent.parent / 'data' / 'Master Project Data' 
nfirs_path =  data_path / 'NFIRS Fire Incident Data.csv'

# List the columns you want to download from the NFIRS csv
cols_to_use = ['state','fdid','inc_date','oth_inj','oth_death','prop_loss',
               'cont_loss','tot_loss','geoid']

# Specify particular data type for geoid column
col_dtypes = {'geoid':str}

# Read in NFIRS dataframe
nfirs = pd.read_csv(nfirs_path,
                    dtype = col_dtypes,
                    usecols = cols_to_use,
                    encoding='latin-1')

# Convert inc_date column values to python datetime type
nfirs['inc_date'] = pd.to_datetime(nfirs['inc_date'], infer_datetime_format=True)



#Read in ACS dataframe
ACS_path = data_path  / 'ACS 5YR Block Group Data.csv'
ACS = pd.read_csv(ACS_path,
                  dtype = {'GEOID':'object'},
                  index_col = 1)
            
#NFIRS Munging

# Ensure correct calculation of tot_loss column 
nfirs['tot_loss'] = nfirs['prop_loss'] + nfirs['cont_loss']

# Create mask for new death or injuries variable
death_inj_mask = (nfirs['oth_death'] > 0) | (nfirs['oth_inj'] > 0)

# By default assigns values of death or injuries column as not severe
nfirs['had_death_injury'] = 'no_death_injury'

# Applies filter to death or injuries column to label the death or injuries instances correctly
nfirs.loc[death_inj_mask,'had_death_injury'] = 'had_death_injury'

# Create new NFIRS variables based on specified thresholds of existing variables in dataframe
nfirs['had_inj'] = np.where(nfirs['oth_inj']>0,'had_inj','no_inj')
nfirs['had_death'] = np.where(nfirs['oth_death']>0,'had_death','no_death')

# Extract just the numeric portion of the geoid
nfirs['geoid'] =  nfirs['geoid'].str[2:-1] 
nfirs =  nfirs[~ nfirs['geoid'].isna()]

# Add a year column to be used to groupby in addition to geoid
nfirs['year'] = nfirs['inc_date'].dt.year.astype('str')
nfirs.set_index('geoid',inplace = True)

# ACS Munging

# Ensures GEOID variable is in the correct format and sets it as the dataframe index
ACS['GEOID'] = ACS['GEOID'].str[2:-1]   
ACS.set_index(['GEOID'],inplace = True)

# Removes extraneous columns in the dataframe
if 'Unnamed: 0' in ACS.columns:
    ACS.drop('Unnamed: 0','columns',inplace= True)

if 'NAME' in ACS.columns:
    ACS.drop('NAME','columns',inplace= True)

if 'inc_pcincome' in ACS.columns:
    ACS.drop('inc_pcincome','columns',inplace= True)

# Creates vector of total populations for each census tract to be used to normalize total fires per year variable
tot_pop = ACS[['tot_population']].groupby('GEOID').sum()

# Drop all total count columns in ACS and keeps all percentage columns
cols = ACS.columns.to_list()
for col in cols:
    if  col.find('tot') != -1: 
        ACS.drop(col,'columns', inplace = True)

# Integer indexing for all rows, but gets rid of county_name, state_name, and in_poverty       
ACS = ACS.iloc[:,3:]
ACS = ACS.groupby('GEOID').mean()

# Remove missing values from dataframe
ACS.replace([np.inf, -np.inf], np.nan,inplace = True)
ACS.dropna(inplace = True)

## Adjust total fires per year by the population counts

# Drop fires without fatalities
# nfirs_mort = nfirs.loc[nfirs['had_death'] == "had_death"]

# Creates dataframe that shows the number of deaths in each census block each year
fires = pd.crosstab(nfirs.index, nfirs['had_death_injury'])
fires.index.rename('GEOID',inplace = True)

# Grab total population values pulled from ACS dataframe and assign to each census block in NFIRS dataframe
fires = fires.merge(tot_pop, how = 'left', left_index = True, right_index = True)

# Remove resulting NaN/infinity values following merge
fires.replace([np.inf, -np.inf], np.nan,inplace = True)
fires.dropna(inplace = True)

# drop rows with no population count
fires = fires[fires['tot_population'] != 0 ]

# population adjustment
fires = fires.div(fires['tot_population'], axis = 'index')*1000
# fires = fires.loc[:,:'2016']

# view fires by year across geoids; displays additional information regarding # of fires in higher percentile categories
fires.describe(percentiles=[.75, .85, .9 ,.95, .99])
print(fires)
# drop unused variables
fires = fires.drop(columns=['no_death_injury', 'tot_population'])

# define variables to indicate census blocks in the top 10% percent of fire risk scores
top10 = fires > fires.quantile(.9)
fires.describe()
# function to create histogram of with # of fires on x-axis and # of census blocks on y-axis
def plotFires(df):

    figsize = (20, 16)
    cols = 4
    rows = 2
    f, axs = plt.subplots(cols,rows,figsize= figsize)


    cases = df.columns.to_list()
    for  case in enumerate(cases):
        ax = plt.subplot(cols,rows,case[0]+1)
        ax.set_title('All Fires {}'.format(str(case[1])) )
        plt.hist(df[case[1]],bins=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,15,20,40,80,100])

# Find correlated features in ACS dataset and identify the highly correlated relationships

# Based on feature importance experiments, select features with consistence importance across annual predictions

# Create ACS correlation matrix
corr = ACS.corr()

# Generate a mask for the upper triangle
mask = np.triu(np.ones_like(corr, dtype=np.bool))

# Set up the matplotlib figure
f, ax = plt.subplots(figsize=(11, 9))

# Generate a custom diverging colormap
cmap = sns.diverging_palette(220, 10, as_cmap=True)

#Draw the heatmap with the mask and correct aspect ratio
sns.heatmap(corr, mask=mask, cmap=cmap, vmax=.3, center=0,
            square=True, linewidths=.5, cbar_kws={"shrink": .5})
 
# Filtering out lower/upper triangular duplicates 
corr_high = corr[abs(corr) > 0.7].stack().reset_index()
corr_high = corr_high[corr_high['level_0'].astype(str)!=corr_high['level_1'].astype(str)]
corr_high['ordered-cols'] = corr_high.apply(lambda x: '-'.join(sorted([x['level_0'],x['level_1']])),axis=1)
corr_high = corr_high.drop_duplicates(['ordered-cols'])
corr_high.drop(['ordered-cols'], axis=1, inplace=True)
corr_high.columns = ['Pair Var 1', 'Pair Var 2', 'Corr Value']

# Display highly correlated pairs
print(corr_high.sort_values(by=['Corr Value'], ascending=False))

# From highly correlated pairs, remove one of the Pair Vars from the ACS dataset except for the 'mort' variables
ACS = ACS.drop(['house_pct_vacant', 'did_not_work_past_12_mo', 'house_pct_non_family', 'house_pct_rent_occupied',
                           'race_pct_nonwhite', 'race_pct_nonwhitenh', 'house_pct_incomplete_plumb',
                           'house_pct_incomplete_kitchen', 'race_pct_whitenh'], axis=1) 

# Based on feature importance experiments, select features with consistent importance across annual predictions
ACS = ACS[['house_yr_pct_earlier_1939', 'house_pct_occupied', 'house_pct_family_married', 'race_pct_black',
          'worked_past_12_mo', 'heat_pct_fueloil_kerosene', 'educ_bachelors', 'house_pct_live_alone', 
          'educ_some_col_no_grad', 'house_pct_ownd_occupied', 'house_w_home_equity_loan', 'house_val_175K_200K',
           'house_val_200K_250K']]

# Function to upsample or downsample our dataframe features if we have unbalanced classes

def resample_df(X,y,upsample=True,seed = SEED):
    from sklearn.utils import resample
    # check which of our two classes is overly represented 
    if np.mean(y) > .5:
        major,minor = 1,0
    else:
        major,minor = 0, 1
    
    # Add Class feature to dataframe equal to our existing dependent variable
    X['Class'] = y
    
    df_major = X[X.Class == major ]
    df_minor = X[X.Class == minor ]
    

    if upsample:      
    
        df_minor_resampled = resample(df_minor,
                                     replace = True,
                                     n_samples = df_major.shape[0], 
                                     random_state = seed)
    
    
   
        combined = pd.concat([df_major,df_minor_resampled])
        
        # Debug
        #print('minor class {}, major class {}'.format(df_minor_resampled.shape[0],
                                                       #df_major.shape[0]))
    
        
    else: # downsample
         
        df_major_resampled = resample(df_major,
                                     replace = False,
                                     n_samples = df_minor.shape[0],
                                     random_state = seed)
        
        
        combined = pd.concat([df_major_resampled,df_minor])
        
        #print('minor class {}, major class {}'.format(df_minor.shape[0],
                                                      #df_major_resampled.shape[0]))


    
    
    y_out = combined['Class']
    X_out = combined.drop('Class', axis =1)
    return X_out , y_out
top10.mean()
# Function to train model that predicts whether each census block is in the top 10% percent of fire risk scores
def train_model(top10, fires, ACS = pd.DataFrame(), nyears = 4, modeltype='LogisticRegression', resample = False, seed = SEED):
    from scipy.stats import zscore
    
    
    # Define model types & parameters 
    
    if modeltype =='LogisticRegression':
        from sklearn.linear_model import LogisticRegression
        model = LogisticRegression(warm_start=True,
                                   class_weight = 'balanced',
                                   max_iter = 1000)

        
    elif modeltype =='BalBagged':
        from imblearn.ensemble import BalancedBaggingClassifier
        from sklearn.tree import DecisionTreeClassifier
        model = BalancedBaggingClassifier(base_estimator=DecisionTreeClassifier(),
                                          n_estimators=80, sampling_strategy='auto',
                                          random_state=0)
        
    elif modeltype =='BalRF':
        from imblearn.ensemble import BalancedRandomForestClassifier
        model = BalancedRandomForestClassifier(n_estimators=80, sampling_strategy='auto',
                                               max_depth=10, random_state=0,
                                              max_features=None, min_samples_leaf=40)

    elif modeltype =='Bagged':
        from sklearn.ensemble import BaggingClassifier
        from sklearn.tree import DecisionTreeClassifier
        model = BaggingClassifier(base_estimator=DecisionTreeClassifier(),
                                    n_estimators=40,
                                    random_state=0)
    
    elif modeltype =='RF':
        from imblearn.ensemble import BalancedRandomForestClassifier
        model = BalancedRandomForestClassifier(n_estimators=60,
                                          warm_start = False,
                                          max_depth = 10,
                                            random_state = 0)
    elif modeltype =='RFtest':
        from sklearn.ensemble import RandomForestClassifier
        model = RandomForestClassifier(n_estimators=60, class_weight = "balanced", max_depth = 10 )
    
    
    # Create framework to predict whether a given census block has a fire risk score in the 90th percentile 
    # based on the specific number of previous years' data
        
#     years = top10.columns
       
#         X = fires.iloc[:,X_start:X_end].copy()
    X = fires.copy()

#     L = X.shape[1] 
#     X.columns = ['year-{}'.format(L - year) for year in range(L)]
#     sm = np.sum(X, axis = 1 )
#     mu = np.mean(X, axis = 1)
#     mx = np.max(X, axis =1)
#     X['Sum']  = sm
#     X['Mean'] = mu
#     X['Max']  = mx
    # set y to target variable
    y = top10.had_death_injury

    # merge in ACS Data into X unless NFIRS-Only model
    if not ACS.empty:
#         X = X[['Sum','Mean','Max']] # drop all other NFIRS columns that have low feature importance scores
        X = X.merge(ACS, how ='left',left_index = True, right_index = True)
        X = X.dropna()
        # drop target variable
        X = X.drop(columns="had_death_injury")
        y = y[y.index.isin(X.index)]
        print(X.shape, y.shape)



    # Create 80/20 training/testing set split
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = .2 )

    # Perform resampling if data classes are unbalanced
    if resample:
        X_train, y_train = resample_df(X_train,y_train)


    # Perform cross-validation 

    #scaler = preprocessing.StandardScaler().fit(X)
    #scaler.transform(X)
    #print ('Cross Val Score:')
    #print(cross_val_score(model, X, y))


    # Standardize features by removing the mean and scaling to unit variance

    scaler = preprocessing.StandardScaler().fit(X_train)
    scaler.transform(X_train)
    scaler.transform(X_test)


    # Fit model to training set

    print('Predicting')
    model = model.fit(X_train,y_train)

    # Calculate training set performance

    train_prediction_probs = model.predict_proba(X_train)
    train_predictions = model.predict(X_train)
    print (confusion_matrix(y_train, train_predictions))
    print (roc_auc_score(y_train, train_prediction_probs[:,1]))


    # Calculate test set performance

    test_prediction_probs = model.predict_proba(X_test)
    test_predictions = model.predict(X_test)
    print (confusion_matrix(y_test, test_predictions))
    print (roc_auc_score(y_test, test_prediction_probs[:,1]))
    print (classification_report(y_test,test_predictions))
    print (log_loss(y_test,test_predictions))


    #Calculate feature importance for each model

    if modeltype=="LogisticRegression":
        feature_importance = {}
        for coef, feat in zip(abs(model.coef_[0]),X_test.columns.tolist()):
            feature_importance[feat] = coef
        print("Feature ranking:")
        print (feature_importance)
    else:
        if modeltype=="RFtest" or modeltype=="BalRF":
            importances = model.feature_importances_
        elif modeltype=="Bagged":
            importances = np.mean([model.estimators_[i].feature_importances_ for i 
                           in range(len(model.estimators_))], axis=0)
        elif modeltype=="BalBagged":
            importances = np.mean([model.estimators_[i].steps[1][1].feature_importances_ for i 
                           in range(len(model.estimators_))], axis=0)

        indices = np.argsort(importances)[::-1]
        print("Feature ranking:")
        for f in range(len(X_test.columns)):
            print("%d. %s (%f)" % (f + 1, X_test.columns[indices[f]], importances[indices[f]]))

    # Increment sliding prediction window
    #start_pointer += 1
#     end_pointer += 1
#     y_pointer += 1

    
            
    return model, X_test, y_test
        

# Train NFIRS Only Model and output prediction performance metrics for each year
# mdl,X_test,y_test=train_model(top10.loc[:,'2009':'2016'],fires,nyears = 5, modeltype='RF', resample = True)

# Train NFIRS + ACS Model and output prediction performance metrics for each year
mdl, X_test, y_test = train_model(top10, fires, ACS = ACS, nyears=1, modeltype='LogisticRegression', resample = False)
# Train NFIRS Only Model and output prediction performance metrics for each year
# mdl,X_test,y_test=train_model(top10.loc[:,'2009':'2016'],fires,nyears = 5, modeltype='RF', resample = True)

# Train NFIRS + ACS Model and output prediction performance metrics for each year
mdl, X_test, y_test = train_model(top10, fires, ACS = ACS, nyears=1, modeltype='RFtest', resample = False)
# Train NFIRS Only Model and output prediction performance metrics for each year
# mdl,X_test,y_test=train_model(top10.loc[:,'2009':'2016'],fires,nyears = 5, modeltype='RF', resample = True)

# Train NFIRS + ACS Model and output prediction performance metrics for each year
mdl, X_test, y_test = train_model(top10, fires, ACS = ACS, nyears=1, modeltype='Bagged', resample = False)
# Train NFIRS Only Model and output prediction performance metrics for each year
# mdl,X_test,y_test=train_model(top10.loc[:,'2009':'2016'],fires,nyears = 5, modeltype='RF', resample = True)

# Train NFIRS + ACS Model and output prediction performance metrics for each year
mdl, X_test, y_test = train_model(top10, fires, ACS = ACS, nyears=1, modeltype='BalBagged', resample = False)<jupyter_output>(69945, 13) (69945,)
Predicting
{'fraction_non_alphanumeric': 0.10605614737213818, 'fraction_numerical': 0.015596492260608557, 'mean_word_length': 3.206981926219361, 'pattern_counts': {'":': 5, '<': 4, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 29, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '161055', 'license': 'MIT', 'n_tokens_mistral': 5640, 'n_tokens_neox': 5190, 'n_words': 1400, 'path': 'Code/Models/NFIRS_Fatalities_Block_level.ipynb', 'repo_name': 'dschon/rcp2', 'size': '16991'}