<jupyter_start><jupyter_text>Power Outages* **See the main project notebook for instructions to be sure you satisfy the rubric!*** See Project 03 for information on the dataset.* A few example prediction questions to pursue are listed below. However, don't limit yourself to them!    * Predict the severity (number of customers, duration, or demand loss) of a major power outage.    * Predict the cause of a major power outage.    * Predict the number and/or severity of major power outages in the year 2020.    * Predict the electricity consumption of an area.Be careful to justify what information you would know at the "time of prediction" and train your model using only those features.  Summary of Findings IntroductionIn this project, we explore the major power outage data in the continental U.S. from January 2000 to July 2016 and use the classification prediction model to predict if an outage will be serious or not. Setting the median of customer affected as threshold, we used a Binarizer by labeling the outage with the customer affected larger than the median as 1 and 0 vice versa. We assigned the results as "Outcomes" in our dataset which also stands as our targeting feature. We decided to run the baseline model using one hot ecoding on categorical variable and dropped the features which has large amount (>300) of nan values since missingness does not help us to predict our outcome. Randomly splitting the data into 70% training set and 30% testing set, the baseline model gives us a train accuracy around 81% and test accuracy around 70% which can definitely be improved. In final model, we added the status of the population in our model and log transform IND.CUST.PCT in order to improve our model. We obtain a test accuracy around 76% which is definitely an improvement. Lastly, we did a fairness analysis to to check if our model is fair when evaluating the severity in coastal and inland states using permutation test on  Baseline ModelWe dropped three features which had large amount of nan values thus leaving us 50 features to predict if the outage is serious or not. After dropping these features, we also dropped the rows that contained nan values to make it easier for us to predict our model. Later, we converted the features that are supposed to be quantitative from objects to float. We then apply one hot encoding to the categorical features and put the adjusted data into Column Transformer and Pipeline using the Decision Tree Classifier which gives us the test accuracy around 84% and test accuracy around 70% which definitely can be improved.Below are the name of the features and whether it is nominal(N), ordinal(O), or quantitative(Q): 'YEAR': O, 'MONTH': N, NERC.REGION': N 'CLIMATE.REGION': N, 'ANOMALY.LEVEL': Q, 'CAUSE.CATEGORY': N, 'CAUSE.CATEGORY': N, 'OUTAGE.DURATION': Q, 'CUSTOMERS.AFFECTED': Q, 'RES.PRICE': Q, 'COM.PRICE: Q, 'IND.PRICE': Q, 'TOTAL.PRICE': Q, 'RES.SALES': Q, 'COM.SALES': Q, 'IND.SALES': Q, 'TOTAL.SALES': Q, 'RES.PERCEN': Q, 'COM.PERCEN': Q, 'IND.PERCEN': Q, 'RES.CUSTOMERS': Q, 'COM.CUSTOMERS': Q, 'IND.CUSTOMERS': Q, 'TOTAL.CUSTOMERS': Q ,'RES.CUST.PCT': Q, 'COM.CUST.PCT': Q, 'IND.CUST.PCT': Q, 'PC.REALGSP.REL':Q, 'PC.REALGSP.CHANGE': Q, 'UTIL.REALGSP': Q, 'TOTAL.REALGSP': Q, 'UTIL.CONTRI': Q, 'PI.UTIL.OFUSA': Q,'POPULATION': Q, 'POPPCT_URBAN': Q, 'POPPCT_UC': Q, 'POPDEN_URBAN': Q, 'POPDEN_UC':Q, 'POPDEN_RURAL': Q, 'AREAPCT_URBAN': Q, 'AREAPCT_UC': Q, 'PCT_LAND': Q, 'PCT_WATER_TOT': Q, 'PCT_WATER_INLAND': Q, 'OUTAGE.START': O,'OUTAGE.RESTORATION': O, 'Outcome': O Final ModelIn the final model, we added the feature POPULATION_CATEGORICAl which used the quatiles of population as thresholds to determine if the population is Very_Little_People, Little_People, Many_People or Very_Many_People. The reason for us to add this feature was because . After running the model, we obtain the train accuracy around 83% and test accuracy around 72% which is only a slightly improvement. We also find out one quantitative feature, IND.CUST.PCT, with large standard deviation and highly skewed which may cause a problem to our prediction; therefore, we log transform this feature using FunctionTransformer and repeated the steps as baseline model. After using GridSearchCV, we obtain the best max_depth to be 3 and obtain the test accuray around 76% which is an improvement comparing to our test accuracy in the baseline model. Fairness EvaluationWe evaluate our model for "fairness" on the 'U.S.State' subset of the data using a permutation test. We want to explore whether our model is fairer for inland states or costal state. Thus, we binarize our subset by setting the States that are inland as 'inland' and the States that are costal as 'costal'. Now my dataset is split into subsets: a costal subset and an inland subset. These are our X and Y.Our null and alternative tests are states as below with the significance level of 0.01:- Null Hypothesis: our model is fair, as the precision for my two subsets are roughly the same- Alternative Hypothesis: our model is unfair, as the precision for the 'coastal' subset is different from the "inland" subsetWe first created a column 'is_coastal', using classification from relavent [Wikipedia page] (https://en.wikipedia.org/wiki/List_of_U.S._states_and_territories_by_coastline). Here, our methodology is to shuffle the 'is_coastal' column and verify if our observative accuracy is similar enough to other permutated result by checking it poistion in the simulated distribution. Once we got the subsets of the data with no Null values, we ran our final model on the subset of the data and found the difference in accuracy score for the coastal and inland group in the model. Then we store this difference of score in an array. We ran this proccess of getting the test statistic 1000 times to get a distribution of test statistics. We found our p-value based on the proportion of test statistics above a pre-determined threshold. We chose our threshold to be 0.01. Comparing to our siginificance level 0.49 (<0.01), we fail to reject the null hypothesis and conclude that our model is fair when predicting the seriousness of the outage for coastal and inland states.  Code<jupyter_code>import matplotlib.pyplot as plt
import numpy as np
import os
import pandas as pd
import seaborn as sns
import sklearn.metrics as metrics

from sklearn.preprocessing import Binarizer 
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import FunctionTransformer
from sklearn.preprocessing import OneHotEncoder
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.tree import DecisionTreeClassifier
from sklearn.model_selection import train_test_split
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import accuracy_score
from sklearn.neighbors import KNeighborsClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn import metrics 

%matplotlib inline
%config InlineBackend.figure_format = 'retina'  # Higher resolution figures<jupyter_output><empty_output><jupyter_text>Cleaning and EDAIt is vital for us to do data cleaning before performing EDA and further testing since the original data contains invalid column names and some invalid rows.<jupyter_code>#Read in Excel file
df = pd.read_excel('outage.xlsx')
df.head()
# Get all the variables
variables = df.loc[4]
# Set the column names to variables
df.columns = variables
# Extract the units for the variables
units = df.loc[5].values[1:]
# Drop the useless columns 
df = df.drop([0, 1, 2, 3, 4, 5])
df = df.drop(columns = 'variables')
# Reset the index
df = df.set_index('OBS')

# Combine OUTAGE START DATE and OUTAGE START TIME to one column
df['OUTAGE.START'] = pd.to_datetime(df['OUTAGE.START.DATE'].apply(lambda x: str(x)) + " " + df['OUTAGE.START.TIME'].apply(lambda x: str(x)), errors='coerce')
# Combine OUTAGE RESTORATION DATE and OUTAGE RESTORATION TIME to one column
df['OUTAGE.RESTORATION'] = pd.to_datetime(df['OUTAGE.RESTORATION.DATE'].apply(lambda x: str(x)) + " " + df['OUTAGE.RESTORATION.TIME'].apply(lambda x: str(x)), errors='coerce')

# Drop the orginial columns 
df = df.drop(columns = ['OUTAGE.START.DATE', 'OUTAGE.START.TIME', 'OUTAGE.RESTORATION.DATE','OUTAGE.RESTORATION.TIME' ])
# Convert all numerical features to float
for col in df.columns:
    if col not in ['YEAR', 'MONTH', 'U.S._STATE', 'POSTAL.CODE', 'NERC.REGION','CLIMATE.REGION', 'CLIMATE.CATEGORY', 'CAUSE.CATEGORY',
                   'CAUSE.CATEGORY.DETAIL', 'HURRICANE.NAMES','OUTAGE.START','OUTAGE.RESTORATION']:
        df[col] = df[col].astype(float)
df.head(10)<jupyter_output><empty_output><jupyter_text>Baseline Model<jupyter_code># Drop the rows with nan values in CUSTOMERS.AFFECTED
df = df.dropna(subset = ['CUSTOMERS.AFFECTED'])
# Get the median of CUSTOMERS.AFFECTED
cusaff_median = df['CUSTOMERS.AFFECTED'].median()
cusaff_median
# Create a binarizer
transformer = Binarizer(threshold = cusaff_median).fit(df[['CUSTOMERS.AFFECTED']])
# Create outcomes based on data
outcomes = transformer.transform(df[['CUSTOMERS.AFFECTED']])
# Adding outcomes to our dataframe
df['Outcome'] = outcomes
df[['CUSTOMERS.AFFECTED','Outcome']].head()<jupyter_output><empty_output><jupyter_text>Determine which features to drop<jupyter_code>df.isna().sum()<jupyter_output><empty_output><jupyter_text>We observed that features like HURRICANE.NAMES, CAUSE.CATEGORY.DETAIL, and DEMAND.LOSS.MW have more than 300 nan values which will not be sufficient if we included them into our prediction model. Thus we decided to drop these three features.<jupyter_code>df = df.drop(columns = ['HURRICANE.NAMES', 'CAUSE.CATEGORY.DETAIL', 'DEMAND.LOSS.MW'])<jupyter_output><empty_output><jupyter_text>We first drop the nan values in the cleaned dataset. Then, keep the numerical features as original and apply one hot encoding to the categorical features. Later, we put the adjusted data into Column Transformer and Pipeline using the Decision Tree Classifier which gives us the test accuracy around 81% and test accuracy around 70%.<jupyter_code>states = df[['U.S._STATE']] # store it for fairness test
df = df.drop(columns = ['U.S._STATE','POSTAL.CODE'])<jupyter_output><empty_output><jupyter_text>Then, we dropped the 'U.S._STATE','POSTAL.CODE' because they can cause overfitting.<jupyter_code># Drop all the remain nans in the dataset
df = df.dropna()
# do a training test split
X = df.drop('Outcome', axis=1)
y = df['Outcome']

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3)
steps = []
for col in df.columns:
    data_type = type(df[col].iloc[0])
    if data_type == str:
        steps.append((col,OneHotEncoder(handle_unknown='ignore'),[col]))

col_t = ('col_trans',ColumnTransformer(steps))
pl = Pipeline([col_t,('tree', DecisionTreeClassifier())])

pl.fit(X_train, y_train)
train_accuracy = pl.score(X_train, y_train)
test_accuracy = pl.score(X_test, y_test)
train_accuracy, test_accuracy
# look at the importance of each feature in our tree model
# _: attributes created when model is fit. Contains info that it learned from the data
tree = pl['tree']
dict(zip(df.drop('Outcome', axis=1).columns, tree.feature_importances_))<jupyter_output><empty_output><jupyter_text>The baseline model has max tree depth around 13 - 17 and 150- 170 nodes.<jupyter_code># maximum depth
print("Max tree depth:")
print(tree.tree_.max_depth)
# node count
print("Number of nodes:")
print(tree.tree_.node_count)<jupyter_output>Max tree depth:
14
Number of nodes:
163<jupyter_text>Final Model POPULATION is impactful when predicting our model. As we can see from the histogram, population that is in between quantile 1 and quantile 3 usually have a higher chance to encounter serious outrage. This is reasonable since we define the seriousness of the outage using the customers that were being affected.Therefore, we can generalize our model to make better predicitions if we categorize POUPLATION to different category instead of being extra specific on listing out the population. We used the quatiles of population as thresholds to determine if the population is Very_Little_People, Little_People, Many_People or Very_Many_People. After running the model with this new feature again, we obtain the test accuracy around 83% and test accuracy around 72%<jupyter_code>sns.histplot(data=df, x="POPULATION", hue="Outcome")
df.columns
Q1 = df['POPULATION'].quantile(q = 0.25)
Q2 = df['POPULATION'].quantile(q = 0.5)
Q3 = df['POPULATION'].quantile(q = 0.75)
def convert(x):
    if x <= Q1:
        return "Very_Little_People"
    elif x <= Q2:
        return "Little_People"
    elif x <= Q3:
        return "Many_People"
    else:
        return "Very_Many_People"
df['POPULATION'] = df['POPULATION'].apply(convert)
# do a training test split
X = df.drop('Outcome', axis=1)
y = df['Outcome']

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3)
steps = []
for col in df.columns:
    data_type = type(df[col].iloc[0])
    if data_type == str:
        steps.append((col,OneHotEncoder(handle_unknown='ignore'),[col]))

col_t = ('col_trans',ColumnTransformer(steps))
pl = Pipeline([col_t,('tree', DecisionTreeClassifier())])

pl.fit(X_train, y_train)
train_accuracy = pl.score(X_train, y_train)
test_accuracy = pl.score(X_test, y_test)
train_accuracy, test_accuracy<jupyter_output><empty_output><jupyter_text>From baseline model, we identify the feature 'IND.CUST.PCTA' as one of the most impactful feature when predicting our outcomes. As we can see from the plot below, there is a skewness in the data with large values, thus, applying log transformation to this feature will make this feature more usable for our model.<jupyter_code>df['IND.CUST.PCT'].plot()<jupyter_output><empty_output><jupyter_text>We can see a higher accuracy after applying log transformation to IND.CUST.PCT. The testing accuracy jumps from around 70% to 76% after we apply grid search to find the best parameters!<jupyter_code># do a training test split
X = df.drop('Outcome', axis=1)
y = df['Outcome']

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3)
steps = [('log_trans',FunctionTransformer(np.log),['IND.CUST.PCT'])]
for col in df.columns:
    data_type = type(df[col].iloc[0])
    if data_type == str:
        steps.append((col,OneHotEncoder(handle_unknown='ignore'),[col]))

col_t = ('col_trans',ColumnTransformer(steps))
pl = Pipeline([col_t,('tree', DecisionTreeClassifier())])

pl.fit(X_train, y_train)
train_accuracy = pl.score(X_train, y_train)
test_accuracy = pl.score(X_test, y_test)
train_accuracy, test_accuracy<jupyter_output><empty_output><jupyter_text>We then used gridsearch to avoid overfitting.<jupyter_code>parameter = {'tree__max_depth': [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,None]}
clf = GridSearchCV(estimator=pl, param_grid=parameter, cv = 5)
clf.fit(X_train, y_train)
clf.best_params_
test_accuracy = clf.score(X_test, y_test)
clf.best_score_, test_accuracy<jupyter_output><empty_output><jupyter_text>Fairness Evaluation In this section we want to check if our model is fair when evaluating the severity in coastal and inland states. (based on this [https://en.wikipedia.org/wiki/List_of_U.S._states_and_territories_by_coastline])- Null Hypothesis: our model is fair, as the precision for my two subsets are roughly the same- Alternative Hypothesis: our model is unfair, as the precision for the 'coastal' subset is different from the "inland" subset<jupyter_code># accuracy
clf.score(X_train,y_train)
# confusion matrix
preds = clf.predict(X_test)
metrics.confusion_matrix(y_test, preds) / len(preds)
# Recall: proportion of the number of true severe outage events predicted 
# divided by all occurence of severe outages.

metrics.recall_score(y_test, preds)
# Specificity: proportion of non-servere outages that were correctly identified.

metrics.recall_score(y_test, preds, pos_label=0)
# Precision: proportion of predicted severe outages that were actually severe

metrics.precision_score(y_test, preds)
# FDR: proportion of predicted severe outages that are not severe

1 - metrics.precision_score(y_test, preds)<jupyter_output><empty_output><jupyter_text>Calculate Parity Measures* A = { coastal city = 'coastal' & inland city = 'inland'}* Y = severe outage (1.0) or non-severe outage (0.0)* C = did we predict it as severe (1.0) or not (0.0)<jupyter_code># get the U.s. State we have in the df
results = pd.concat([df, states], axis=1) # concat df with states
np.unique(results['U.S._STATE'])
results['is_coastal'] = (results['U.S._STATE'].isin(['Alabama', 'Arkansas', 'California', 
       'Connecticut', 'Delaware', 'Florida', 'Georgia', 
       'Indiana', 'Louisiana', 'Maine', 'Maryland', 'Massachusetts', 
       'Mississippi', 'New Hampshire', 'New Jersey', 'New York', 'North Carolina',
       'South Carolina', 'Texas', 'Virginia', 'Washington'])).replace({True:'coastal', False:'inland'})
# drop states
results = results.drop(columns = ['U.S._STATE'])
# Drop all the remain nans in the dataset
results = results.dropna()

# do a training test split
X = results.drop('Outcome', axis=1)
y = results['Outcome']

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3)
steps = [('log_trans',FunctionTransformer(np.log),['IND.CUST.PCT'])]
for col in results.columns:
    data_type = type(results[col].iloc[0])
    if data_type == str:
        steps.append((col,OneHotEncoder(handle_unknown='ignore'),[col]))

col_t = ('col_trans',ColumnTransformer(steps))
pl = Pipeline([col_t,('tree', DecisionTreeClassifier())])

parameter = {'tree__max_depth': [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,None]}
clf = GridSearchCV(estimator=pl, param_grid=parameter, cv = 5)
clf.fit(X_train, y_train)
clf.best_params_
results
# Demographic Parity 
# Inland States have a better chance to encounter serious outage

preds = clf.predict(X_test)
results = X_test
results['prediction'] = preds
results['tag'] = y_test
results.groupby('is_coastal').prediction.mean().to_frame()
# Accuracy Parity

(
    results
    .groupby('is_coastal')
    .apply(lambda x: metrics.accuracy_score(x.tag, x.prediction))
    .rename('accuracy')
    .to_frame()
)<jupyter_output><empty_output><jupyter_text>Is this difference in accuracy significant?* Use a permutation test:    - are the distributions of severe outages "the same" for coastal/ inland groups?    - test-statistic: "the same" is measured in terms of different notions of parity.         - (using 'metrics.accuracy_score()' to obtain the difference in two prediction scores)* Set a significance level of 0.01<jupyter_code>obs = results.groupby('is_coastal').apply(lambda x: metrics.accuracy_score(x.tag, x.prediction)).diff().iloc[-1]

metrs = []
# permutate 1000 times
for _ in range(100):
    s = (
        results[['is_coastal', 'prediction', 'tag']]
        .assign(is_coastal=results.is_coastal.sample(frac=1.0, replace=False).reset_index(drop=True))
        .groupby('is_coastal')
        .apply(lambda x: metrics.accuracy_score(x.tag, x.prediction))
        .diff()
        .iloc[-1]
    )
    
    metrs.append(s)
print(pd.Series(metrs <= obs).mean())
pd.Series(metrs).plot(kind='hist', bins=20, title='Permutation Test for severity prediction between coastal/ inland groups')
plt.scatter(obs, 0, c='r');<jupyter_output>0.49
{'fraction_non_alphanumeric': 0.10126975900492355, 'fraction_numerical': 0.01135009069707178, 'mean_word_length': 5.353638459005597, 'pattern_counts': {'":': 0, '<': 64, '<?xml version=': 0, '>': 60, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '219995', 'license': 'MIT', 'n_tokens_mistral': 6259, 'n_tokens_neox': 5942, 'n_words': 2177, 'path': 'outages_proj5.ipynb', 'repo_name': 'DylanTao/DSC80', 'size': '19295'}