PythonCode_IncomePrediction.py

# -*- coding: utf-8 -*-
"""IncomePredict.ipynb

Automatically generated by Colaboratory.

#**Import Library**
"""

import pandas as pd
import numpy as np
import seaborn as sn
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.style as style

from sklearn import preprocessing
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.cluster import KMeans, Birch
from sklearn.preprocessing import MinMaxScaler, StandardScaler

from __future__ import print_function
from sklearn.datasets import make_blobs
from sklearn.metrics import silhouette_samples, silhouette_score

"""# **Data Preprocessing**
Training Set
"""

from google.colab import files
uploaded = files.upload()

"""###Data Cleaning"""

df_train = pd.read_csv("CensusCanada2016Training.csv",encoding='latin-1') #Reads csv as pandas dataframe
df_train.head()

df_train.describe()

df_train.shape

#df_train.info()
df = df_train.copy()

#renaming column names
df.columns = ['Total_population', 'Total_households', 'Median_income', 'THC', 
              'THC_before1961','THC1961_1980','THC1981_1990','THC1991_2000','THC2001_2005',
              'TH_houses','TH_apartment','TH_other_dewelling',
              'TH_tenure','TH_owner','TH_renter']
df.info()

#create several columns serve for later calculation
df['THC_before1991'] = df['THC_before1961'] + df['THC1961_1980'] + df['THC1981_1990']
df['THC_before2001'] = df['THC_before1991'] + df['THC1991_2000']
df['THC_before2006'] = df['THC_before2001'] + df['THC2001_2005']
df['THC2006_2016'] = df['Total_households'] - df['THC_before2006']

"""## Derive Input Variables"""

# Derive additional columns
# % increase in total household for period of construction
df['%increase_1991_2000'] = df.THC1991_2000 / df.THC_before1991
df['%increase_2001_2005'] = df.THC2001_2005 / df.THC_before2001
df['%increase_2006_2016'] = df.THC2006_2016 / df.THC_before2006

# % Households by Structure Type
df['%Houses'] = df.TH_houses / df.Total_households
df['%Apartment'] = df.TH_apartment / df.Total_households
df['%Other_dwelling'] = df.TH_other_dewelling / df.Total_households

# % Households by Tenure
df['%Owner'] = df.TH_owner / df.Total_households
df['%Renter'] = df.TH_renter / df.Total_households

#Drop unused columns
df2 = df.drop(columns = ['THC','THC_before1961','THC1961_1980','THC1981_1990','THC1991_2000','THC2001_2005',
                          'TH_houses','TH_apartment','TH_other_dewelling','TH_tenure','TH_owner','TH_renter',
                          'THC_before1991','THC_before2001','THC_before2006','THC2006_2016'])
df2.info()

#replace null value with 0
df2.replace([np.nan, np.inf, -np.inf], 0, inplace=True)
np.where(pd.isnull(df2))

df2.describe()

#outlier in %Owner & %Renter has value greater than 1
outlier = df2[(df2['%Owner'] > 1) | (df2['%Renter'] > 1)]
outlier

#drop outlier value
df2.drop(outlier.index, inplace=True)
df2 = df2.reset_index(drop=True)
df2

df2.describe()

"""## Correlation analysis"""

#correlation matrix
input = df2.drop(columns ='Median_income')
corr_matrix = input.corr()
corr_matrix

fig, ax = plt.subplots(figsize=(11,8))
sn.set(font_scale=1)
sn.heatmap(corr_matrix, annot=True, fmt=".2f", cmap="Blues") #Blues, flare, crest, vlag, 
plt.show()

sn.set_style("ticks")
sn.pairplot(input)

#Total population vs Total households
fig, ax = plt.subplots(figsize=(5,5))
sc = sn.scatterplot(x="Total_population", y="Total_households", data=df2)
sc.set_title("Total population vs. Total households")

#"%Renter vs. %Owner"
fig, ax = plt.subplots(figsize=(5,5))
sc = sn.scatterplot(x="%Renter", y="%Owner", data=df2)
sc.set_title("%Owner vs %Renter ")

#"%Houses vs. %Apartment"
fig, ax = plt.subplots(figsize=(5,5))
sc = sn.scatterplot(x="%Houses", y="%Apartment", data=df2)
sc.set_title("%Houses vs. %Apartment")

#drop highly correlated variables
df3 = df2.drop(columns = ['Total_population','%Apartment','%Other_dwelling', '%Renter'])
input = df3.drop(columns ='Median_income')

input.corr()

"""#**Clustering Model**

## Silhouette analysis
"""

#target variable
y = df3.loc[:,'Median_income']

#scale the dataset use Min-Max Scaler
scaler = preprocessing.MinMaxScaler()
X = input
X = scaler.fit_transform(X)

#Elbow Method:
ks = range(1,7)
inertias = []

for k in ks:
    model = KMeans(n_clusters=k)
    model.fit(X)
    inertias.append(model.inertia_)

plt.plot(ks, inertias, '-o')
plt.xlabel('number of clusters, k')
plt.ylabel('inertia')
plt.grid(False)
plt.xticks(ks)

#silhouette analysis
range_n_clusters = range(2,9)

for n_clusters in range_n_clusters:
    # Create a subplot with 1 row and 2 columns
    fig, (ax1) = plt.subplots(1, 1)
    fig.set_size_inches(7, 5)

    # The silhouette coefficient can range from -1, 1 but in this example all
    # lie within [-0.1, 1]
    ax1.set_xlim([-0.1, 1])
    # The (n_clusters+1)*10 is for inserting blank space between silhouette
    # plots of individual clusters, to demarcate them clearly.
    ax1.set_ylim([0, len(X) + (n_clusters + 1) * 10])

    # Initialize the clusterer with n_clusters value and a random generator
    # seed of 11 for reproducibility.
    clusterer = KMeans(n_clusters=n_clusters, random_state=11)
    cluster_labels = clusterer.fit_predict(X)

    # The silhouette_score gives the average value for all the samples.
    # This gives a perspective into the density and separation of the formed
    # clusters
    silhouette_avg = silhouette_score(X, cluster_labels)
    print("For n_clusters =", n_clusters,
          "The average silhouette_score is :", silhouette_avg)

    # Compute the silhouette scores for each sample
    sample_silhouette_values = silhouette_samples(X, cluster_labels)

    y_lower = 10
    for i in range(n_clusters):
        # Aggregate the silhouette scores for samples belonging to
        # cluster i, and sort them
        ith_cluster_silhouette_values = \
            sample_silhouette_values[cluster_labels == i]

        ith_cluster_silhouette_values.sort()

        size_cluster_i = ith_cluster_silhouette_values.shape[0]
        y_upper = y_lower + size_cluster_i

        color = cm.nipy_spectral(float(i) / n_clusters)
        ax1.fill_betweenx(np.arange(y_lower, y_upper),
                          0, ith_cluster_silhouette_values,
                          facecolor=color, edgecolor=color, alpha=0.7)

        # Label the silhouette plots with their cluster numbers at the middle
        ax1.text(-0.05, y_lower + 0.5 * size_cluster_i, str(i))

        # Compute the new y_lower for next plot
        y_lower = y_upper + 10  # 10 for the 0 samples

    ax1.set_title("The silhouette plot for the various clusters.")
    ax1.set_xlabel("The silhouette coefficient values")
    ax1.set_ylabel("Cluster label")

    # The vertical line for average silhouette score of all the values
    ax1.axvline(x=silhouette_avg, color="red", linestyle="--")

    ax1.set_yticks([])  # Clear the yaxis labels / ticks
    ax1.set_xticks([-0.1, 0, 0.2, 0.4, 0.6, 0.8, 1])

    plt.suptitle(("Silhouette analysis for KMeans clustering on sample data "
                  "with n_clusters = %d" % n_clusters),
                 fontsize=14, fontweight='bold')

plt.show()

"""## K-Means clustering model"""

#select k=3
kmeans = KMeans(n_clusters=3, random_state=5)
pipeline = make_pipeline(kmeans)
pipeline.fit(X)
labels = pipeline.predict(X)
kmeans_cluster = pd.DataFrame(data=labels)

#add columns for kmeans cluster
df4 = df3.copy()
df4['kmeans_cluster'] = kmeans_cluster[0]
#rename the cluster to 1, 2, 3
df4.replace({'kmeans_cluster':{0: 1, 1: 2, 2:3}}, inplace=True)
df4['kmeans_cluster'].value_counts()

pd.Series.sort_index(df4['kmeans_cluster'].value_counts())

"""## Clusters Profile"""

#Exploring Cluster
df4['kmeans_cluster'].value_counts().plot(kind='bar', title= "Type of clusters", color=(0.2, 0.4, 0.6, 0.6))

#aggregate statistics by cluster
df_cluster_summary = df4.groupby('kmeans_cluster').describe().T.reset_index()
df_cluster_summary = df_cluster_summary.rename(columns={'level_0':'column','level_1':'metric'})
# mean value for each cluster
df_cluster_summary = df_cluster_summary[df_cluster_summary['metric'] == "mean"]
df_cluster_summary = df_cluster_summary.set_index('column')
df_cluster_summary

# join into single summary dataset
df_profile_overall = df4.describe().T
df_profile = df_cluster_summary.join(df_profile_overall) # joins on Index
df_profile

# %House by cluster
sn.boxplot(x="kmeans_cluster", y="%Houses", data=df4, palette='Set2', hue="kmeans_cluster")
plt.legend(loc='upper left', title='Cluster')
plt.title("% of Houses of total households by Cluster")

# %Owner by cluster
sn.boxplot(x="kmeans_cluster", y="%Owner", data=df4, palette='Set2', hue="kmeans_cluster") #, color='skyblue', palette='hls', Set2, Blues, Paired 
plt.legend(loc='upper left', title='Cluster')
plt.title("% of Owner of total households by Cluster")

Cluster_1 = df4[df4['kmeans_cluster'] == 1]
Cluster_2 = df4[df4['kmeans_cluster'] == 2]
Cluster_3 = df4[df4['kmeans_cluster'] == 3]

plt.style.use('default')
sn.set_style("ticks")
#Median Income by cluster
plt.figure(figsize=(8,5))
plt.hist(Cluster_1['Median_income'], bins=100, alpha=0.5, label="Cluster_1")
plt.hist(Cluster_2['Median_income'], bins=100, alpha=0.5, label="Cluster_2")
plt.hist(Cluster_3['Median_income'], bins=100, alpha=0.5, label="Cluster_3")
plt.xlabel("Median Income", size=14)
plt.title("Median Income by Cluster")
plt.legend(loc='upper right')
plt.xlim([0, 200000])

#Total Households by cluster
plt.figure(figsize=(8,5))
plt.hist(Cluster_1['Total_households'], bins=100, alpha=0.5, label="Cluster_1")
plt.hist(Cluster_2['Total_households'], bins=100, alpha=0.5, label="Cluster_2")
plt.hist(Cluster_3['Total_households'], bins=100, alpha=0.5, label="Cluster_3")
plt.xlabel("Total householdse", size=14)
plt.title("Total households by Cluster")
plt.legend(loc='upper right')
plt.xlim([0, 7000])

#%Houses by cluster
plt.figure(figsize=(8,5))
plt.hist(Cluster_1['%Houses'], bins=50, alpha=0.5, label="Cluster_1")
plt.hist(Cluster_2['%Houses'], bins=50, alpha=0.5, label="Cluster_2")
plt.hist(Cluster_3['%Houses'], bins=50, alpha=0.5, label="Cluster_3")
plt.xlabel("%Houses", size=14)
plt.title("%Houses by Cluster")
plt.legend(loc='upper right')

#%Owner by cluster
plt.figure(figsize=(8,5))
plt.hist(Cluster_1['%Owner'], bins=50, alpha=0.5, label="Cluster_1")
plt.hist(Cluster_2['%Owner'], bins=50, alpha=0.5, label="Cluster_2")
plt.hist(Cluster_3['%Owner'], bins=50, alpha=0.5, label="Cluster_3")
plt.xlabel("%Owner", size=14)
plt.title("%Owner by Cluster")
plt.legend(loc='upper right')

for c in df4:
    grid= sn.FacetGrid(df4, col='kmeans_cluster')
    grid.map(plt.hist, c)    
    plt.show()

"""## BIRCH clustering model"""

#BIRCH clustering
#threshod=0.1
brc = Birch(n_clusters=3, threshold=0.1)
brc.fit(X)
BIRCH = brc.predict(X)
BIRCH_cluster = pd.DataFrame(data=BIRCH)

#add columns for BIRCH cluster
df6 = df4.copy()
df6['BIRCH'] = BIRCH_cluster[0]

#rename the cluster to 1, 2, 3
df6.replace({'BIRCH':{0: 2, 2:3}}, inplace=True)
df6['BIRCH'].value_counts()

pd.Series.sort_index(df6['BIRCH'].value_counts())

df6.head()

#aggregate statistics by cluster
df_cluster_summary2 = df6.groupby('BIRCH').describe().T.reset_index()
df_cluster_summary2 = df_cluster_summary2.rename(columns={'level_0':'column','level_1':'metric'})
# mean value for each cluster
df_cluster_summary2 = df_cluster_summary2[df_cluster_summary2['metric'] == "mean"]
df_cluster_summary2 = df_cluster_summary2.set_index('column')
df_cluster_summary2

df_cluster_summary2.T

#Compare to K-Means
df_cluster_summary

# join into single summary dataset
df_profile_overall2 = df6.describe().T
df_profile2 = df_cluster_summary2.join(df_profile_overall2) # joins on Index
df_profile2

# %House by cluster
sn.boxplot(x="BIRCH", y="%Houses", data=df6, palette='Paired', hue = "BIRCH") #, color='skyblue', palette='hls', Set2, Paired

# %Owner by cluster
sn.boxplot(x="BIRCH", y="%Owner", data=df6, palette='Paired', hue = "BIRCH") #, color='skyblue', palette='hls', Set2, Blues, Paired

Cluster_1_B = df6[df6['BIRCH'] == 1]
Cluster_2_B = df6[df6['BIRCH'] == 2]
Cluster_3_B = df6[df6['BIRCH'] == 3]

#Median Income by cluster
plt.figure(figsize=(8,5))
plt.hist(Cluster_1_B['Median_income'], bins=100, alpha=0.5, label="Cluster_1_BIRCH")
plt.hist(Cluster_2_B['Median_income'], bins=100, alpha=0.5, label="Cluster_2_BIRCH")
plt.hist(Cluster_3_B['Median_income'], bins=100, alpha=0.5, label="Cluster_3_BIRCH")
plt.xlabel("Median Income", size=14)
plt.title("Median Income by Cluster")
plt.legend(loc='upper right')
plt.xlim([0, 250000])

#Total Household by cluster
plt.figure(figsize=(8,5))
plt.hist(Cluster_1['Total_households'], bins=70, alpha=0.5, label="Cluster_1_BIRCH")
plt.hist(Cluster_2['Total_households'], bins=70, alpha=0.5, label="Cluster_2_BIRCH")
plt.hist(Cluster_3['Total_households'], bins=70, alpha=0.5, label="Cluster_3_BIRCH")
plt.xlabel("Total_households", size=14)
plt.title("BIRCH: Total_households by Cluster")
plt.legend(loc='upper right')
plt.xlim([0, 8000])

#%House by cluster
plt.figure(figsize=(8,5))
plt.hist(Cluster_1_B['%Houses'], bins=50, alpha=0.5, label="Cluster_1_BIRCH")
plt.hist(Cluster_2_B['%Houses'], bins=50, alpha=0.5, label="Cluster_2_BIRCH")
plt.hist(Cluster_3_B['%Houses'], bins=50, alpha=0.5, label="Cluster_3_BIRCH")
plt.xlabel("%Houses", size=14)
plt.title("BIRCH: %Houses by Cluster")
plt.legend(loc='upper right')

#%Owner by cluster
plt.figure(figsize=(8,5))
plt.hist(Cluster_1_B['%Owner'], bins=50, alpha=0.5, label="Cluster_1_BIRCH")
plt.hist(Cluster_2_B['%Owner'], bins=50, alpha=0.5, label="Cluster_2_BIRCH")
plt.hist(Cluster_3_B['%Owner'], bins=50, alpha=0.5, label="Cluster_3_BIRCH")
plt.xlabel("%Owner", size=14)
plt.title("BIRCH: %Owner by Cluster")
plt.legend(loc='upper right')

"""# **Segmentation Modeling: KNN + Regression Tree + Random Forest**

### Split training and validation sets
"""

# dataframe for part 2
df_p2 = df4.copy()
df_p2

X = df_p2.drop(columns=["Median_income"])
y = df_p2.drop(columns=["Total_households", "%increase_1991_2000", "%increase_2001_2005", "%increase_2006_2016", "%Houses", "%Owner"])

test_size = 0.25
seed = 88

# split data into training and test set
X_train, X_test, y_train, y_test = train_test_split(X,y,test_size=test_size, random_state=seed)

X_train.shape

X_test.shape

X_train

y_train

y_test

"""### Clustering training set and validation set 



"""

# identify segments, ie. clusters for training set 
X_cluster1_train = X_train[X_train['kmeans_cluster'] == 1.0]
X_cluster2_train = X_train[X_train['kmeans_cluster'] == 2.0]
X_cluster3_train = X_train[X_train['kmeans_cluster'] == 3.0]

X_cluster1_train= X_cluster1_train.drop(columns="kmeans_cluster")
X_cluster2_train= X_cluster2_train.drop(columns="kmeans_cluster")
X_cluster3_train= X_cluster3_train.drop(columns="kmeans_cluster")

y_cluster1_train = y_train[y_train['kmeans_cluster'] == 1.0]
y_cluster2_train = y_train[y_train['kmeans_cluster'] == 2.0]
y_cluster3_train = y_train[y_train['kmeans_cluster'] == 3.0]

y_cluster1_train= y_cluster1_train.drop(columns="kmeans_cluster")
y_cluster2_train= y_cluster2_train.drop(columns="kmeans_cluster")
y_cluster3_train= y_cluster3_train.drop(columns="kmeans_cluster")

# overview of the training set
train_summary = X_train.groupby('kmeans_cluster').describe().T.reset_index()
train_summary = train_summary.rename(columns={'level_0':'column','level_1':'metric'})

train_summary = train_summary[train_summary['metric'] == "mean"]
train_summary = train_summary.set_index('column')
train_summary

X_train = X_train.drop(columns="kmeans_cluster")
y_train = y_train.drop(columns="kmeans_cluster")

# identify segments, ie. clusters for validation set 
X_cluster1_test = X_test[X_test['kmeans_cluster'] == 1.0]
X_cluster2_test = X_test[X_test['kmeans_cluster'] == 2.0]
X_cluster3_test = X_test[X_test['kmeans_cluster'] == 3.0]

X_cluster1_test= X_cluster1_test.drop(columns="kmeans_cluster")
X_cluster2_test= X_cluster2_test.drop(columns="kmeans_cluster")
X_cluster3_test= X_cluster3_test.drop(columns="kmeans_cluster")

y_cluster1_test = y_test[y_test['kmeans_cluster'] == 1.0]
y_cluster2_test = y_test[y_test['kmeans_cluster'] == 2.0]
y_cluster3_test = y_test[y_test['kmeans_cluster'] == 3.0]

y_cluster1_test= y_cluster1_test.drop(columns="kmeans_cluster")
y_cluster2_test= y_cluster2_test.drop(columns="kmeans_cluster")
y_cluster3_test= y_cluster3_test.drop(columns="kmeans_cluster")

# overview of the validation set
validation_summary = X_test.groupby('kmeans_cluster').describe().T.reset_index()
validation_summary = validation_summary.rename(columns={'level_0':'column','level_1':'metric'})

validation_summary = validation_summary[validation_summary['metric'] == "mean"]
validation_summary = validation_summary.set_index('column')
validation_summary

X_test = X_test.drop(columns="kmeans_cluster")
y_test = y_test.drop(columns="kmeans_cluster")

"""### Training and Validation Set Size"""

from tabulate import tabulate
data = [['Total', X_train.shape[0], X_test.shape[0]],
['Cluster 1', X_cluster1_train.shape[0], X_cluster1_test.shape[0]],
['Cluster 2', X_cluster2_train.shape[0], X_cluster2_test.shape[0]],
['Cluster 3', X_cluster3_train.shape[0], X_cluster3_test.shape[0]]]
print (tabulate(data, headers=["Dataset", "Training Set Size", "Validation Set Size"]))

"""## KNN

### A Global KNN model on training set
"""

from numpy import arange
from pandas import read_table
from pandas import set_option
from pandas.plotting import scatter_matrix

from sklearn.preprocessing import MinMaxScaler
from sklearn.model_selection import train_test_split
from sklearn.model_selection import KFold
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import GridSearchCV
from sklearn.neighbors import KNeighborsRegressor
from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import r2_score
from sklearn.metrics import confusion_matrix

# Rescale Training X
scaler = MinMaxScaler().fit(X_train)
rescaled_X_train = scaler.transform(X_train)

# Rescale Test X
scaler = MinMaxScaler().fit(X_test)
rescaled_X_test = scaler.transform(X_test)

seed = 88
num_folds = 5

scoring = 'neg_mean_absolute_error'

k_values = np.arange(1,30) # try k values

param_grid = dict(n_neighbors=k_values)

model = KNeighborsRegressor()
kfold = KFold(n_splits = num_folds, shuffle = True, random_state = seed)

grid = GridSearchCV(estimator = model, param_grid = param_grid, scoring = scoring, cv = kfold)
grid_result = grid.fit(rescaled_X_train, y_train)

print("Best: %f using %s" %(grid_result.best_score_, grid_result.best_params_))

# all of the means, standard deviations, and parameters
means = grid_result.cv_results_['mean_test_score']
stds = grid_result.cv_results_['std_test_score']
params = grid_result.cv_results_['params']

plt.figure(figsize = (8,5))
plt.plot(k_values, means)
plt.xlabel("Number of Neighbors K")
plt.ylabel("Negative Mean Absolute Error")

best_neighbour= grid_result.best_estimator_.get_params()['n_neighbors']

model = KNeighborsRegressor(n_neighbors = best_neighbour)

# take all the training data and fit the model
model.fit(rescaled_X_train, y_train)

# calculate training MAE
y_pred_train = model.predict(rescaled_X_train)
train_MAE_global = round(mean_absolute_error(y_train, y_pred_train),2)
train_MAE_global

# calculate test MAE
y_pred_test = model.predict(rescaled_X_test)
test_MAE_global = round(mean_absolute_error(y_test, y_pred_test),2)
test_MAE_global

"""### KNN Segmentation model on cluster 1"""

# Rescale Training X for Cluster 1
scaler = MinMaxScaler().fit(X_cluster1_train)
rescaled_X_train1 = scaler.transform(X_cluster1_train)

# Rescale Test X for Cluster 1
scaler = MinMaxScaler().fit(X_cluster1_test)
rescaled_X_test1 = scaler.transform(X_cluster1_test)

seed = 88
num_folds = 5

scoring = 'neg_mean_absolute_error'

k_values = np.arange(1,30) # try k values

param_grid = dict(n_neighbors=k_values)

model = KNeighborsRegressor()
kfold = KFold(n_splits = num_folds, shuffle = True, random_state = seed)

grid = GridSearchCV(estimator = model, param_grid = param_grid, scoring = scoring, cv = kfold)
grid_result = grid.fit(rescaled_X_train1, y_cluster1_train)

print("Best: %f using %s" %(grid_result.best_score_, grid_result.best_params_))

# all of the means, standard deviations, and parameters
means = grid_result.cv_results_['mean_test_score']
stds = grid_result.cv_results_['std_test_score']
params = grid_result.cv_results_['params']

plt.figure(figsize = (8,5))
plt.plot(k_values, means)
plt.xlabel("Number of Neighbors K")
plt.ylabel("Negative Mean Absolute Error")

best_neighbour1= grid_result.best_estimator_.get_params()['n_neighbors']
model1 = KNeighborsRegressor(n_neighbors = best_neighbour1)

# take cluster 1 training data and fit the model
model1.fit(rescaled_X_train1, y_cluster1_train)

# calculate training MAE
y_pred_train1 = model1.predict(rescaled_X_train1)
train_MAE1 = round(mean_absolute_error(y_cluster1_train, y_pred_train1),2)
train_MAE1

# calculate test MAE
y_pred_test1 = model1.predict(rescaled_X_test1)
test_MAE1 = round(mean_absolute_error(y_cluster1_test, y_pred_test1),2)
test_MAE1

"""### KNN Segmentation model on cluster 2"""

# Rescale Training X for Cluster 2
scaler = MinMaxScaler().fit(X_cluster2_train)
rescaled_X_train2 = scaler.transform(X_cluster2_train)

# Rescale Test X for Cluster 2
scaler = MinMaxScaler().fit(X_cluster2_test)
rescaled_X_test2 = scaler.transform(X_cluster2_test)

seed = 88
num_folds = 5

scoring = 'neg_mean_absolute_error'

k_values = np.arange(1,30) # try k values

param_grid = dict(n_neighbors=k_values)

model = KNeighborsRegressor()
kfold = KFold(n_splits = num_folds, shuffle = True, random_state = seed)

grid = GridSearchCV(estimator = model, param_grid = param_grid, scoring = scoring, cv = kfold)
grid_result = grid.fit(rescaled_X_train2, y_cluster2_train)

print("Best: %f using %s" %(grid_result.best_score_, grid_result.best_params_))

# all of the means, standard deviations, and parameters
means = grid_result.cv_results_['mean_test_score']
stds = grid_result.cv_results_['std_test_score']
params = grid_result.cv_results_['params']

plt.figure(figsize = (8,5))
plt.plot(k_values, means)
plt.xlabel("Number of Neighbors K")
plt.ylabel("Negative Mean Absolute Error")

best_neighbour2= grid_result.best_estimator_.get_params()['n_neighbors']
model2 = KNeighborsRegressor(n_neighbors = best_neighbour2)

# take cluster 2 training data and fit the model
model2.fit(rescaled_X_train2, y_cluster2_train)

# calculate training MAE
y_pred_train2 = model2.predict(rescaled_X_train2)
train_MAE2 = round(mean_absolute_error(y_cluster2_train, y_pred_train2),2)
train_MAE2

# calculate test MAE
y_pred_test2 = model2.predict(rescaled_X_test2)
test_MAE2 = round(mean_absolute_error(y_cluster2_test, y_pred_test2),2)
test_MAE2

"""### KNN Segmentation model on cluster 3"""

# Rescale Training X for Cluster 3
scaler = MinMaxScaler().fit(X_cluster3_train)
rescaled_X_train3 = scaler.transform(X_cluster3_train)

# Rescale Test X for Cluster 3
scaler = MinMaxScaler().fit(X_cluster3_test)
rescaled_X_test3 = scaler.transform(X_cluster3_test)

seed = 88
num_folds = 5

scoring = 'neg_mean_absolute_error'

k_values = np.arange(1,30) # try k values

param_grid = dict(n_neighbors=k_values)

model = KNeighborsRegressor()
kfold = KFold(n_splits = num_folds, shuffle = True, random_state = seed)

grid = GridSearchCV(estimator = model, param_grid = param_grid, scoring = scoring, cv = kfold)
grid_result = grid.fit(rescaled_X_train3, y_cluster3_train)

print("Best: %f using %s" %(grid_result.best_score_, grid_result.best_params_))

# all of the means, standard deviations, and parameters
means = grid_result.cv_results_['mean_test_score']
stds = grid_result.cv_results_['std_test_score']
params = grid_result.cv_results_['params']

plt.figure(figsize = (8,5))
plt.plot(k_values, means)
plt.xlabel("Number of Neighbors K")
plt.ylabel("Negative Mean Absolute Error")

best_neighbour3= grid_result.best_estimator_.get_params()['n_neighbors']
model3 = KNeighborsRegressor(n_neighbors = best_neighbour3)

# take cluster 3 training data and fit the model
model3.fit(rescaled_X_train3, y_cluster3_train)

# calculate training MAE
y_pred_train3 = model3.predict(rescaled_X_train3)
train_MAE3 = round(mean_absolute_error(y_cluster3_train, y_pred_train3),2)
train_MAE3

# calculate test MAE
y_pred_test3 = model3.predict(rescaled_X_test3)
test_MAE3 = round(mean_absolute_error(y_cluster3_test, y_pred_test3),2)
test_MAE3

"""### Compare KNN Global Model and Segmentation Models"""

print("Global KNN n_neighbors =", best_neighbour, '\n'
      "Cluster 1 KNN n_neighbors =", best_neighbour1, '\n'
      "Cluster 2 KNN n_neighbors =", best_neighbour2, '\n'
      "Cluster 3 KNN n_neighbors =", best_neighbour3)

from tabulate import tabulate
data = [['Training Data', train_MAE_global, test_MAE_global],
['Cluster 1',train_MAE1, test_MAE1],
['Cluster 2', train_MAE2, test_MAE2],
['Cluster 3', train_MAE3, test_MAE3]]
print (tabulate(data, headers=["Dataset", "Training MAE", "Testing MAE"]))

"""## Regression Tree"""

!pip install dmba

import graphviz
import json
from sklearn import tree
from sklearn.tree import DecisionTreeClassifier, DecisionTreeRegressor,export_graphviz
from sklearn.model_selection import train_test_split, cross_val_score, GridSearchCV
from sklearn.metrics import r2_score, accuracy_score, mean_absolute_error
from dmba import classificationSummary
from sklearn.model_selection import train_test_split

"""### A Global Regression Tree model on training set

"""

#min_samples_leaf: the minimum number of samples required to be at a leaf node. Minimum 50 as per assignment requirements
RegressionTree = DecisionTreeRegressor(min_samples_leaf=50 ,random_state=88)

RegressionTree.fit(X_train, y_train)

y_predicted = RegressionTree.predict(X_test)

expected_y  = y_test

print("Max depth: " + str(RegressionTree.tree_.max_depth))
print("Node count: " + str(RegressionTree.tree_.node_count))

print("Tree description: " )
RegressionTree.get_params()

dot_data = tree.export_graphviz(RegressionTree, out_file=None, 
                                feature_names=X_train.columns,  
                                filled=True)
graphviz.Source(dot_data, format="png")

# Evaluating the Model
print("Mean absolute errors for:")
rt_global_training_mae = mean_absolute_error(y_train, RegressionTree.predict(X_train))
rt_global_testing_mae =mean_absolute_error(expected_y, y_predicted)
print("   Global training: "+ str(rt_global_training_mae))
print("   Global testing: " + str(rt_global_testing_mae))

rt_global_r2 = r2_score(expected_y, y_predicted)
rt_global_r2

"""### Regression Tree Segmentation model on cluster 1

"""

#min_samples_leaf: the minimum number of samples required to be at a leaf node. Minimum 50 as per assignment requirements
RegressionTree_cluster1 = DecisionTreeRegressor(min_samples_leaf=50 ,random_state=88)

RegressionTree_cluster1.fit(X_cluster1_train, y_cluster1_train)

y_predicted_cluster1 = RegressionTree_cluster1.predict(X_cluster1_test)

expected_y_cluster1  = y_cluster1_test.Median_income

rt_cluster1_err = expected_y_cluster1 - y_predicted_cluster1

print("Max depth: " + str(RegressionTree_cluster1.tree_.max_depth))
print("Node count: " + str(RegressionTree_cluster1.tree_.node_count))

print("Tree description: " )
RegressionTree_cluster1.get_params()

dot_data = tree.export_graphviz(RegressionTree_cluster1, out_file=None, 
                                feature_names=X_cluster1_train.columns,  
                                filled=True)
graphviz.Source(dot_data, format="png")

# Evaluating the Model
print("Mean absolute errors for:")
rt_cluster1_training_mae = mean_absolute_error(y_cluster1_train, RegressionTree_cluster1.predict(X_cluster1_train))
rt_cluster1_testing_mae =mean_absolute_error(expected_y_cluster1, y_predicted_cluster1)
print("   Cluster 1 training: "+ str(rt_cluster1_training_mae))
print("   Cluster 1 testing: " + str(rt_cluster1_testing_mae))

"""### Regression Tree Segmentation model on cluster 2"""

#min_samples_leaf: the minimum number of samples required to be at a leaf node. Minimum 50 as per assignment requirements
RegressionTree_cluster2 = DecisionTreeRegressor(min_samples_leaf=50 ,random_state=88)

RegressionTree_cluster2.fit(X_cluster2_train, y_cluster2_train)

y_predicted_cluster2 = RegressionTree_cluster2.predict(X_cluster2_test)
expected_y_cluster2  = y_cluster2_test

import statistics
statistics.median(expected_y_cluster2)

print("Max depth: " + str(RegressionTree_cluster2.tree_.max_depth))
print("Node count: " + str(RegressionTree_cluster2.tree_.node_count))

dot_data = tree.export_graphviz(RegressionTree_cluster2, out_file=None, 
                                feature_names=X_cluster2_train.columns,  
                                filled=True)
graphviz.Source(dot_data, format="png")

# Evaluating the Model
print("Mean absolute errors for:")
rt_cluster2_training_mae = mean_absolute_error(y_cluster2_train, RegressionTree_cluster2.predict(X_cluster2_train))
rt_cluster2_testing_mae =mean_absolute_error(expected_y_cluster2, y_predicted_cluster2)
print("   Cluster 2 training: "+ str(rt_cluster2_training_mae))
print("   Cluster 2 testing: " + str(rt_cluster2_testing_mae))

"""### Regression Tree Segmentation model on cluster 3

"""

#min_samples_leaf: the minimum number of samples required to be at a leaf node. Minimum 50 as per assignment requirements
RegressionTree_cluster3 = DecisionTreeRegressor(min_samples_leaf=50 ,random_state=88)

RegressionTree_cluster3.fit(X_cluster3_train, y_cluster3_train)

y_predicted_cluster3 = RegressionTree_cluster3.predict(X_cluster3_test)
expected_y_cluster3  = y_cluster3_test

statistics.median(expected_y_cluster3)

print("Max depth: " + str(RegressionTree_cluster3.tree_.max_depth))
print("Node count: " + str(RegressionTree_cluster3.tree_.node_count))

dot_data = tree.export_graphviz(RegressionTree_cluster3, out_file=None, 
                                feature_names=X_cluster3_train.columns,  
                                filled=True)
graphviz.Source(dot_data, format="png")

# Evaluating the Model
print("Mean absolute errors for:")
rt_cluster3_training_mae = mean_absolute_error(y_cluster3_train, RegressionTree_cluster3.predict(X_cluster3_train))
rt_cluster3_testing_mae =mean_absolute_error(expected_y_cluster3, y_predicted_cluster3)
print("   Cluster 3 training: "+ str(rt_cluster3_training_mae))
print("   Cluster 3 testing: " + str(rt_cluster3_testing_mae))

"""## Random Forest

"""

from sklearn.ensemble import RandomForestRegressor # for building the model

"""### A Global Random Forest Model on training set"""

# Initializing the Random Forest Regression model with 10 decision trees
rf = RandomForestRegressor(n_estimators = 128, random_state = 0)

# Fitting the Random Forest Regression model to the data
rf.fit(X_train, y_train)

y_predicted = rf.predict(X_test)
mean_y_global_rf = y_predicted.mean()
y_predicted

y_expected = y_test.Median_income
y_expected

expected_test = y_expected.mean()
expected_test

rf_global_err = y_expected - y_predicted

# Evaluating the Model
print("Mean absolute errors for:")
global_rf_train_mae = mean_absolute_error(y_train, rf.predict(X_train))
global_rf_test_mae =mean_absolute_error(y_expected.to_numpy(), y_predicted)
print("   training: "+ str(global_rf_train_mae))
print("   testing: " + str(global_rf_test_mae))

rf_global_r2 = r2_score(y_expected, y_predicted)

median_income_global_rf = statistics.mean(y_predicted)
median_income_global_rf

"""### Random Forest Segmentation Model on Cluster 1

"""

# Initializing the Random Forest Regression model with 10 decision trees
rf_cluster1 = RandomForestRegressor(n_estimators = 128, random_state = 0)

# Fitting the Random Forest Regression model to the data
rf_cluster1.fit(X_cluster1_train, y_cluster1_train)

y_predicted_cluster1 = rf_cluster1.predict(X_cluster1_test)
mean_y_cluster1_rf = y_predicted_cluster1.mean()

y_expected_cluster1 = y_cluster1_test.Median_income

expected_cluster1 = y_expected_cluster1.mean()
expected_cluster1

rf_cluster1_err = y_expected_cluster1 - y_predicted_cluster1

# Evaluating the Model
print("Mean absolute errors for:")
cluster1_rf_train_mae = mean_absolute_error(y_cluster1_train, rf.predict(X_cluster1_train))
cluster1_rf_test_mae =mean_absolute_error(y_expected_cluster1, y_predicted_cluster1)
print("   training: "+ str(cluster1_rf_train_mae))
print("   testing: " + str(cluster1_rf_test_mae))

rf_cluster1_r2 = r2_score(y_expected_cluster1, y_predicted_cluster1)
rf_cluster1_r2

rf_cluster1.feature_importances_

"""### Random Forest Segmentation Model on Cluster 2

"""

# Initializing the Random Forest Regression model with 10 decision trees
rf_cluster2 = RandomForestRegressor(n_estimators = 128, random_state = 0)

# Fitting the Random Forest Regression model to the data
rf_cluster2.fit(X_cluster2_train, y_cluster2_train)

y_predicted_cluster2 = rf_cluster2.predict(X_cluster2_test)
mean_y_cluster2_rf = y_predicted_cluster2.mean()

y_expected_cluster2 = y_cluster2_test.Median_income

expected_cluster2 = y_expected_cluster2.mean()
expected_cluster2

rf_cluster2_err = y_expected_cluster2 - y_predicted_cluster2

# Evaluating the Model
print("Mean absolute errors for:")
cluster2_rf_train_mae = mean_absolute_error(y_cluster2_train, rf.predict(X_cluster2_train))
cluster2_rf_test_mae =mean_absolute_error(y_expected_cluster2, y_predicted_cluster2)
print("   training: "+ str(cluster2_rf_train_mae))
print("   testing: " + str(cluster2_rf_test_mae))

rf_cluster2_r2 = r2_score(y_expected_cluster2, y_predicted_cluster2)
rf_cluster2_r2

"""### Random Forest Segmentation Model on Cluster 3"""

# Initializing the Random Forest Regression model with 10 decision trees
rf_cluster3 = RandomForestRegressor(n_estimators = 128, random_state = 0)

# Fitting the Random Forest Regression model to the data
rf_cluster3.fit(X_cluster3_train, y_cluster3_train)

y_predicted_cluster3 = rf_cluster3.predict(X_cluster3_test)
mean_y_cluster3_rf = y_predicted_cluster3.mean()

y_expected_cluster3 = y_cluster3_test.Median_income

expected_cluster3 = y_expected_cluster3.mean()
expected_cluster3

rf_cluster3_err = y_expected_cluster3 - y_predicted_cluster3

# Evaluating the Model
print("Mean absolute errors for:")
cluster3_rf_train_mae = mean_absolute_error(y_cluster3_train, rf.predict(X_cluster3_train))
cluster3_rf_test_mae =mean_absolute_error(y_expected_cluster3, y_predicted_cluster3)
print("   training: "+ str(cluster3_rf_train_mae))
print("   testing: " + str(cluster3_rf_test_mae))

rf_cluster3_r2 = r2_score(y_expected_cluster3, y_predicted_cluster3)
rf_cluster3_r2

"""## **Model Comparison**"""

from tabulate import tabulate
data = [
    ["KNN",'Training Data', train_MAE_global, test_MAE_global],
    ["",'Cluster 1', train_MAE1, test_MAE1],
    ["",'Cluster 2', train_MAE2, test_MAE2],
    ["",'Cluster 3', train_MAE3, test_MAE3],
    ["---------------", "-------------", "-----------", "--------"],
    ["Regression Tree",'Training Data', round(rt_global_training_mae,2), round(rt_global_testing_mae,2)],
    ["",'Cluster 1', round(rt_cluster1_training_mae,2), round(rt_cluster1_testing_mae,2)],
    ["",'Cluster 2', round(rt_cluster2_training_mae,2), round(rt_cluster2_testing_mae,2)],
    ["",'Cluster 3', round(rt_cluster3_training_mae,2), round(rt_cluster3_testing_mae,2)],
    ["---------------", "-------------", "-----------", "-------"],
    ["Random Forest",'Training Data', round(global_rf_train_mae,2), round(global_rf_test_mae,2)],
    ["",'Cluster 1', round(cluster1_rf_train_mae,2), round(cluster1_rf_test_mae,2)],
    ["",'Cluster 2', round(cluster2_rf_train_mae,2), round(cluster2_rf_test_mae,2)],
    ["",'Cluster 3', round(cluster3_rf_train_mae,2), round(cluster3_rf_test_mae,2)]
    ]
print (tabulate(data, headers=["Model", "Dataset", "Training MAE", "Testing MAE"]))

"""###Compare 3 model MAEs on 3 Clusters"""

plt.figure(figsize=(12, 8), dpi=80)
N = 3
ind = np.arange(N) 
width = 0.25
  
knn_mae = [test_MAE1, test_MAE2, test_MAE3]
bar1 = plt.bar(ind, knn_mae, width, color = 'lightskyblue')
  
tree_mae = [rt_cluster1_testing_mae, rt_cluster2_testing_mae, rt_cluster3_testing_mae]
bar2 = plt.bar(ind+width, tree_mae, width, color='pink')
  
rf_mae = [cluster1_rf_test_mae, cluster2_rf_test_mae, cluster3_rf_test_mae]
bar3 = plt.bar(ind+width*2, rf_mae, width, color = 'yellowgreen')


plt.xlabel("Clusters")
plt.ylabel('Test MAE')
plt.title("Comparing Test MAE")
plt.xticks(ind+width,['Cluster 1', 'Cluster 2', 'Cluster 3'])
plt.ylim(6000, 19000)
plt.legend( (bar1, bar2, bar3), ('KNN', 'Regression Tree', 'Random Forest') )
plt.show()

"""### Feature Importance"""

#Random Forest on Cluster 2 Feature Importance
rf_cluster2.feature_importances_
col_map = plt.get_cmap('Paired')
color_list = ["lightsteelblue","lightsteelblue", "lightsteelblue", "cornflowerblue","cornflowerblue", "royalblue"]
sorted_idx = rf_cluster2.feature_importances_.argsort()
plt.figure(figsize=(6, 8), dpi=80)
plt.barh(X_cluster2_train.columns[sorted_idx], rf_cluster2.feature_importances_[sorted_idx],color=color_list)
plt.xlim(0,0.7)
plt.xlabel("Cluster 2 - Random Forest Model Feature Importance")
plt.show()

#Regression Tree on Cluster 3 Feature Importance
RegressionTree_cluster3.feature_importances_
col_map = plt.get_cmap('Paired')
color_list = ["lightsteelblue","lightsteelblue", "lightsteelblue", "cornflowerblue","cornflowerblue", "royalblue"]
sorted_idx = RegressionTree_cluster3.feature_importances_.argsort()
plt.figure(figsize=(6, 8), dpi=80)
plt.barh(X_cluster3_train.columns[sorted_idx], RegressionTree_cluster3.feature_importances_[sorted_idx],color=color_list)
plt.xlim(0,0.7)
plt.xlabel("Cluster 3 - Regression Tree Feature Importance")
plt.show()

"""# Test Set"""

from google.colab import files
uploaded = files.upload()

"""## Data PreProcessing and Kmeans Clustering

"""

df_test = pd.read_csv("CensusCanada2016Test.csv",encoding='latin-1') #Reads csv as pandas dataframe
df_test.head()

df_test.shape

df_test.describe()

#renaming column names
df_test.columns = ['Total_population', 'Total_households', 'THC', 
              'THC_before1961','THC1961_1980','THC1981_1990','THC1991_2000','THC2001_2005',
              'TH_houses','TH_apartment','TH_other_dewelling',
              'TH_tenure','TH_owner','TH_renter']
df_test.info()

#create several columns serve for later calculation
df_test['THC_before1991'] = df_test['THC_before1961'] + df_test['THC1961_1980'] + df_test['THC1981_1990']
df_test['THC_before2001'] = df_test['THC_before1991'] + df_test['THC1991_2000']
df_test['THC_before2006'] = df_test['THC_before2001'] + df_test['THC2001_2005']
df_test['THC2006_2016'] = df_test['Total_households'] - df_test['THC_before2006']

# Derive additional columns
df_test['%increase_1991_2000'] = df_test.THC1991_2000 / df_test.THC_before1991
df_test['%increase_2001_2005'] = df_test.THC2001_2005 / df_test.THC_before2001
df_test['%increase_2006_2016'] = df_test.THC2006_2016 / df_test.THC_before2006
df_test['%Houses'] = df_test.TH_houses / df_test.Total_households
df_test['%Owner'] = df_test.TH_owner / df_test.Total_households

#replace null value with 0
df_test.replace([np.nan, np.inf, -np.inf], 0, inplace=True)
np.where(pd.isnull(df_test))

#Drop unused columns
df_test = df_test.drop(columns = ['THC','THC_before1961','THC1961_1980','THC1981_1990','THC1991_2000','THC2001_2005',
                          'TH_houses','TH_apartment','TH_other_dewelling','TH_tenure','TH_owner','TH_renter',
                          'THC_before1991','THC_before2001','THC_before2006','THC2006_2016','Total_population'])
df_test.info()

#K-Means Clustering
scaler = preprocessing.MinMaxScaler() 
X2 = scaler.fit_transform(df_test)
kmeans = KMeans(n_clusters=3, random_state=9)
pipeline = make_pipeline(kmeans)
pipeline.fit(X2)
labels = pipeline.predict(X2)
kmeans_cluster = pd.DataFrame(data=labels)
df_test2 = df_test.copy()
df_test2['kmeans_cluster'] = kmeans_cluster[0]
df_test2.replace({'kmeans_cluster':{0: 1, 1: 2, 2:3}}, inplace=True)
df_test2['kmeans_cluster'].value_counts()
pd.Series.sort_index(df_test2['kmeans_cluster'].value_counts())

df_test2

"""## Predictions via Segmentation Modelling """

# create df for each training cluster
X_test_c1 = df_test2[df_test2['kmeans_cluster'] == 1.0]
X_test_c2 = df_test2[df_test2['kmeans_cluster'] == 2.0]
X_test_c3 = df_test2[df_test2['kmeans_cluster'] == 3.0]

X_test_c1= X_test_c1.drop(columns="kmeans_cluster")
X_test_c2= X_test_c2.drop(columns="kmeans_cluster")
X_test_c3= X_test_c3.drop(columns="kmeans_cluster")

X_test_c1.index

#aggregate statistics by cluster
df_cluster_summary = df_test2.groupby('kmeans_cluster').describe().T.reset_index()
df_cluster_summary = df_cluster_summary.rename(columns={'level_0':'column','level_1':'metric'})
df_cluster_summary = df_cluster_summary[df_cluster_summary['metric'] == "mean"]
df_cluster_summary = df_cluster_summary.set_index('column')
df_cluster_summary

t1 = pd.concat([X_test_c1, X_test_c2, X_test_c3])
t1['index_column'] = t1.index
Index_col = t1.loc[:,'index_column']
Index_col

"""### KNN Segmentation model on cluster 1"""

# Rescale Test X for Cluster 1
scaler = MinMaxScaler().fit(X_test_c1)
knn_rescaled_X_test1 = scaler.transform(X_test_c1)

# calculate test MAE
predicted_cluster1 = model1.predict(knn_rescaled_X_test1)
income_pred = predicted_cluster1

mean_income_cluster1 = predicted_cluster1.mean()

print("Mean predicted median income for cluster 1 is: " + str(mean_income_cluster1))

"""### Random Forest Segmentation Model on Cluster 2

"""

predicted_cluster2 = rf_cluster2.predict(X_test_c2)
mean_income_cluster2 = predicted_cluster2.mean()

income_pred = np.append(income_pred, predicted_cluster2)

print("Mean predicted median income for cluster 2 is: " + str(mean_income_cluster2))

"""### Regression Tree Segmentation model on cluster 3

"""

RegressionTree_cluster3.fit(X_cluster3_train, y_cluster3_train)

predicted_cluster3 = RegressionTree_cluster3.predict(X_test_c3)

mean_income_cluster3 = predicted_cluster3.mean()

income_pred = np.append(income_pred, predicted_cluster3)

print("Mean predicted median income for cluster 3 is: " + str(mean_income_cluster3))

"""### Overview of All Predictions"""

predictions_data = [["cluster 1", mean_income_cluster1], ["cluster 2",mean_income_cluster2], ["cluster 3", mean_income_cluster3]]
print(tabulate(predictions_data, headers=["Clusters", "Predicted Median income"]))

income_pred.shape

df_output = pd.DataFrame()
df_output["index"] = Index_col
df_output["predict"] = income_pred
df_output = df_output.sort_values('index')
df_output

#write to text file
file = open("Team_13_predictions.txt", "w")
for row in df_output["predict"]:
    file.write(str(row) + "\n")
file.close()