<jupyter_start><jupyter_text>Smoke Alarm install model this notebook represents the current smoke alarm install model<jupyter_code>import pandas as pd
import os
import sys
from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt 
import missingno as msno
import geopandas as gpd
# import stored methods 
path = Path.cwd().parent.parent
LoadAndCleanACS = path /'src' /'data'/ 'LoadAndCleanACS.py'
LoadAndCleanARCP = path /'src' /'data'/ 'LoadAndCleanARCP.py'

# Run methods to obtain clean datasets 
%run $LoadAndCleanACS
%run $LoadAndCleanARCP


base_path = Path.cwd().parent.parent
shp_path =  base_path / 'Data' / 'Master Project Data' / 'Tiger_censusBlocks_2016'
shp_files = shp_path.glob('**/*.zip')
shapes = gpd.GeoDataFrame()
for z in shp_files:
    shapes = shapes.append( gpd.read_file( 'zip://' + str(z) ) )  

# to allow for all variables to be displayed in jupyter
pd.set_option('display.max_rows', 500)
pd.set_option('display.max_columns', 500)
def StandardizeColumnNames(df):
    """
    Standardizes column names
    """
    df.columns = map(str.lower, df.columns)
    df.columns = df.columns.str.replace(', ', '_')
    df.columns = df.columns.str.replace('-', '_')
    df.columns = df.columns.str.replace('/', '_')
    df.columns = df.columns.str.replace('(', '_')
    df.columns = df.columns.str.replace(')', '_')
    df.columns = df.columns.str.replace(' ', '_')
    #print(df.columns)
    return df<jupyter_output><empty_output><jupyter_text>Data<jupyter_code>input_loc =  path /'Data'/ 'Master Project Data'
output_loc = path /'Data'/ 'processed'
arc_path = input_loc / 'ARC Preparedness Data.csv'
arc = pd.read_csv(arc_path, 
                  dtype = {'GEOID': str, 'Zip': str})
arc = StandardizeColumnNames(arc)
arc.dropna(inplace = True)
# trim geoid leading saftey marks 
arc['geoid'] = arc['geoid'].str[2:]<jupyter_output><empty_output><jupyter_text>EDA  - remove all houses that don't have a previous smoke detector record - Determine the median number of house visist- Visualize visit distribution - use ACS data to determine % of blocks visited - determine % blocks visited with >15 visits<jupyter_code>arc
#block level
counts = arc['geoid'].value_counts()
counts_median = counts.median()
counts.describe()<jupyter_output><empty_output><jupyter_text>EDA- Geograpic Level repeat block Level analysis at various levels<jupyter_code># county
county_counts =  arc['geoid'].str[:7].value_counts()
print('County Level')
print(county_counts.describe())
# state 
state_counts =  arc['geoid'].str[:2].value_counts()
print('\n State Level')
print(state_counts.describe())
print(arc['pre_existing_alarms'].describe())
print(arc['pre_existing_alarms_tested_and_working'].describe())<jupyter_output>count    537984.000000
mean          1.507301
std           1.552540
min           0.000000
25%           0.000000
50%           1.000000
75%           2.000000
max           9.000000
Name: pre_existing_alarms, dtype: float64
count    537984.000000
mean          0.829638
std           1.336587
min           0.000000
25%           0.000000
50%           0.000000
75%           1.000000
max           9.000000
Name: pre_existing_alarms_tested_and_working, dtype: float64<jupyter_text>Confidence Interval Motivation A commonly used formula for a binomial confidence interval relies on approximating the distribution of error about a binomially-distributed observation, ${\displaystyle {\hat {p}}}$, with a normal distribution. This approximation is based on the central limit theorem and is unreliable when the sample size is small or the success probability is close to 0 or 1.Using the normal approximation, the success probability p is estimated as${\displaystyle {\hat {p}}\pm z{\sqrt {\frac {{\hat {p}}\left(1-{\hat {p}}\right)}{n}}},}$Source https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval<jupyter_code>def CreateConfidenceIntervals(num_surveys,percentage):
# this function takes the cleaned data and adds a confidence interval 

    z =	1.960 # corresponds to 95% confidence interval
    
    CI =  z * np.sqrt(
                     (percentage * (100 - percentage) ) / 
                      num_surveys  )

    return CI<jupyter_output><empty_output><jupyter_text>Feature Engineering Create Binary variables out of the smoke_alarms_present and smoke_alarms_tested_and_working variables   We'll then create a new dataset with the aggregated number and percantage of working smoke detectors in each census geography<jupyter_code>def CreateSingleLevelSmokeAlarmModel(df,geo_level, acs = ACS.copy() ):
# This function takes the arc data  into a dataset containing the percentage 
# and number of smoke detectors by census geography
#
# Inputs 
# arc-  the arc dataset
#
# geo_level- String var indcating what census geography to aggregate on. current levels are:
# State,County,Block,State
#
# The resultant dataset will have the following values:
#
#   num_surveys - total number of surveys conducted
#
#   detectors_found -   houses with at least one smoke detector in the home
#
#   detectors_workding - houses with at least one tested and working smoke detector in the home
#
#   Note: for variables the suffixes 
#       _total- indicates raw counts 
#        _prc  - indicates percentage: (_total / num_surveys * 100)
#
   
    # dict with relevant length of GEOID for tract geography
    geo_level_dict = {'State':2,'County':5,'Tract':11,'Block':12}
    
    df['geoid'] = df['geoid'].str[: geo_level_dict[geo_level]]
    acs.index =  acs.index.str[:geo_level_dict[geo_level]]
    acs.drop_duplicates(inplace = True)
    ## binarize pre_existing_alarms and _tested_and_working
    #  values will now be: 0 if no detectors present and 1 if any number were present
    df['pre_existing_alarms'].where(df['pre_existing_alarms'] < 1, other = 1, inplace = True) 
    df['pre_existing_alarms_tested_and_working'].where(
                                                        df['pre_existing_alarms_tested_and_working'] < 1,
                                                            other = 1, 
                                                            inplace = True)

    ## create detectors dataset
    # This happens by grouping data both on pre_existing alarms and then _tested_and working alarms 
    # and then merging the two into the final dataset

    detectors =  df.groupby('geoid')['pre_existing_alarms'].agg({np.size ,
                                                                  np.sum,
                                                                  lambda x: np.sum(x)/np.size(x)* 100 })

    detectors.rename({'size':'num_surveys','sum':'detectors_found_total','<lambda_0>':'detectors_found_prc'},
                     axis =1,
                     inplace = True)

    detectors['detectors_found_prc'] = detectors['detectors_found_prc'].round(2)
    
  
    
    d2 =  df.groupby('geoid')['pre_existing_alarms_tested_and_working'].agg({np.size,np.sum, 
                                                                              lambda x: np.sum(x)/np.size(x)* 100 })
    
    d2.rename({'size':'num_surveys2','sum':'detectors_working_total','<lambda_0>':'detectors_working_prc'},
                     axis =1,
                     inplace = True)

    
    d2['detectors_working_prc'] = d2['detectors_working_prc'].round(2)
    

    detectors = detectors.merge(d2,how = 'left', on ='geoid')

    detectors['detectors_found_CI'] = CreateConfidenceIntervals(detectors['num_surveys'].values,
                                                                detectors['detectors_found_prc'].values )
                                                                
    detectors['detectors_working_CI'] = CreateConfidenceIntervals(detectors['num_surveys'].values,
                                                                detectors['detectors_working_prc'].values )  
    
    
    
    
    
    # rearrange columns 
    column_order = ['num_surveys',	
                    'detectors_found_total',
                    'detectors_found_prc', 
                    'detectors_found_CI',
                    'detectors_working_total',
                    'detectors_working_prc',
                    'detectors_working_CI']
    
    detectors = detectors[column_order]
    
    detectors = detectors[~pd.isna(detectors.index)]
# fix block model to ensure blocks that weren't visited are added to the model 
    detectors = detectors.reindex(detectors.index.union(acs.index.unique()),fill_value = 0)
    detectors = detectors[~pd.isna(detectors.index)]
   

# test if there are missing values in resultant 

    return detectors
 arc_state = CreateSingleLevelSmokeAlarmModel(arc.copy(),'State',ACS.copy())
 arc_state.head()
arc_county = CreateSingleLevelSmokeAlarmModel(arc.copy(),'County',ACS.copy())
arc_county.describe()
arc_tract = CreateSingleLevelSmokeAlarmModel(arc.copy(),'Tract',ACS.copy())
arc_tract.describe()
arc_block = CreateSingleLevelSmokeAlarmModel(arc.copy(),'Block',ACS.copy())
arc_block.describe()
arc_block
arc_block
arc_county
plt.figure()
plt.hist(np.log(arc_tract['num_surveys']),bins=6)
plt.xlabel("X")
plt.ylabel("Y")
plt.xticks(np.arange(8))
plt.show()
plt.figure()
plt.hist(np.log(arc_block['num_surveys']+1),bins=6)
plt.xlabel("X")
plt.ylabel("Y")
plt.xticks(np.arange(8))
plt.show()
np.arange(8)
arc_block[arc_block['num_surveys'] > 30].shape[0]/ arc_block.shape[0]
ACS.index.isin(arc_block.index).all()
blockGEOID = []
EstimateSurveyNumber = []
EstimateMean = []
EstimateCI   = []
EstimateGeography = []



for index, row in arc_block.iterrows():
    blockGEOID.append(index) 
    # append block GEOID-Data  
    if  len(row) != 0 and  row['num_surveys'] > 30:
        EstimateGeography.append('Block')
        EstimateSurveyNumber.append(row['num_surveys'])
        EstimateMean.append(row['detectors_working_prc'])
        EstimateCI.append(row['detectors_working_CI'])
        continue
    
    if index[:11] in arc_tract.index:
        curr_tract = arc_tract[arc_tract.index  ==  index[:11]].values[0]
        if len(curr_tract) != 0 and curr_tract[0] > 30:
            EstimateGeography.append('Tract')
            EstimateSurveyNumber.append(curr_tract[0])
            EstimateMean.append(curr_tract[-2])
            EstimateCI.append(curr_tract[-1])
            continue
    if index[:5] in arc_county.index:
        curr_county = arc_county[arc_county.index  ==  index[:5]].values[0]
        if len(curr_county) != 0 and  curr_county[0]  > 30:
            EstimateGeography.append('County')
            EstimateSurveyNumber.append(curr_county[0])
            EstimateMean.append(curr_county[-2])
            EstimateCI.append(curr_county[-1])
            continue
    else:
        curr_state = arc_state[arc_state.index  ==  index[:2]].values[0]
        EstimateGeography.append('State')
        EstimateSurveyNumber.append(curr_state[0])
        EstimateMean.append(curr_state[-2])
        EstimateCI.append(curr_state[-1])
    
    # zip and create geoid
    MultiLevelModel = pd.DataFrame( 
     [{'geoid': geoid,'estimate_geography': geo, 'num_surveys': ns , 'detectors_prc_mean': dpm , 'detectors_prc_ci': dci} for 
     geoid,geo,ns,dpm,dci in zip(blockGEOID,EstimateGeography,EstimateSurveyNumber,EstimateMean,EstimateCI) ]
    )
   

import statsmodels.api as sm
p = Path.cwd()
data_path =  p.parent.parent / 'Data' / 'Model Outputs' / 'SmokeAlarmModelOutput.csv'

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

# Read in NFIRS dataframe
MultiLevelModel = pd.read_csv(data_path,
                    dtype = col_dtypes)


# Print out the statistics
MultiLevelModel=MultiLevelModel[['geoid', 'num_surveys','detectors_prc_mean']]
MultiLevelModel['geoid']=MultiLevelModel['geoid'].str[2:]
MultiLevelModel=MultiLevelModel.set_index('geoid')
Output_Var=MultiLevelModel[['detectors_prc_mean']]
#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) 
ACS_SmokeAlarm = ACS[ACS.index.isin(MultiLevelModel.index.tolist())]
Data_Matrix = ACS_SmokeAlarm.merge(Output_Var, how = 'left', left_index = True, right_index = True)
X=Data_Matrix.loc[:, Data_Matrix.columns != 'detectors_prc_mean']
y=Data_Matrix.loc[:, Data_Matrix.columns == 'detectors_prc_mean']
X = sm.add_constant(X)
results = sm.OLS(y,X).fit()
results.summary()

# Using Skicit-learn to split data into training and testing sets
from sklearn.model_selection import train_test_split
from sklearn import preprocessing
# Split the data into training and testing sets
train_features, test_features, train_labels, test_labels = train_test_split(X, 
                y, test_size = 0.30, random_state = 0)
scaler = preprocessing.StandardScaler().fit(train_features)
scaler.transform(train_features)
scaler.transform(test_features)
# Import the model we are using
from sklearn.ensemble import RandomForestRegressor
rf = RandomForestRegressor(n_estimators = 40, random_state = 0)
# Train the model on training data
model=rf.fit(train_features, train_labels['detectors_prc_mean'].ravel())
importances = model.feature_importances_
indices = np.argsort(importances)[::-1]
print("Feature ranking:")
for f in range(len(test_features.columns)):
    print("%d. %s (%f)" % (f + 1, test_features.columns[indices[f]], importances[indices[f]]))
        
# Use the forest's predict method on the test data
predictions = model.predict(test_features)
# Calculate the absolute errors
errors = abs(predictions - test_labels['detectors_prc_mean'].ravel())
# Print out the mean absolute error (mae)
print('Mean Absolute Error:', round(np.mean(errors), 2))
print('Median Absolute Error:', round(np.median(errors), 2))
print('Standard Deviation of Absolute Error:', round(np.std(errors), 2))


#out_path =  path / 'Data' /'Model Outputs'/ 'SmokeAlarmModelOutput.csv'

#MultiLevelModel.to_csv(out_path)
ACS.index.isin(MultiLevelModel['geoid'].str[2:]).all()
MultiLevelModelModel
MultiLevelModel
ACS.size
ACS[~ACS.index.isin(MultiLevelModel.index)]
x=ACS[~ACS.index.isin(MultiLevelModel['geoid'])]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11474030980584565, 'fraction_numerical': 0.017943505992850634, 'mean_word_length': 3.3915050784856877, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 31, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '537134', 'license': 'MIT', 'n_tokens_mistral': 4751, 'n_tokens_neox': 4390, 'n_words': 974, 'path': 'Code/Models/Smoke_Alarm_Model.ipynb', 'repo_name': 'dschon/rcp2', 'size': '14267'}