# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: 'Python 3.9.4 64-bit (''.venv'': venv)'
#     name: python3
# ---

# + [markdown] colab_type="text" id="aBJaE--wKYOJ"
# # Female headed households in South Africa
#
# In  this notebook we will work on the data provided by Zindi during the  competition **Women in Big Data South Africa - Women-Headed Households in South Africa**. We will create our own model to find out which features have an impact on the income of a female-headed household. 

# + colab={} colab_type="code" id="gV8Fpz7vKYOP"
# import packages 
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from scipy.stats import boxcox
import statsmodels.api as sm
from sklearn.model_selection import train_test_split, KFold
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score

# %matplotlib inline
plt.rcParams['figure.figsize'] = [10, 6]

# + [markdown] colab_type="text" id="g6T5r-RjKYON"
# # Loading the data
# We will load the provided csv file from Zindi. In this notebook we will only work with the train data (Train.csv), as the test data (Test.csv) has no target variable and we can't evaluate our predictions. Additionally we will also load the file with the descriptions of the features (descriptions.csv).

# + colab={"base_uri": "https://localhost:8080/", "height": 309} colab_type="code" id="gfGrrqGlKYOV" outputId="c918b49a-67a7-4b1b-81bf-918684a2ca9b"
# import data
df = pd.read_csv('data/Train.csv')
df.head()
# -

# get some general information about the dataset
df.info()

# get a statistical overview of the dataframe
df.describe()

# check for null values
df.isnull().sum().sort_values(ascending=False)

# check for columns with single values
df.nunique().sort_values().head(10)

# print the columns with single values for a visual check
df[['dw_13', 'lan_13', 'dw_12', 'pw_08', 'pw_07']]

# The columns ['dw_13', 'lan_13', 'dw_12', 'pw_08', 'pw_07'] have only 0 in them --> we could drop these features, Also we can drop some other non-relevant features (ward, coordinates, features not given in percentage, Nightlight).

# +
# drop the non-numerical features
df = df.drop(['ward', 'ADM4_PCODE'], axis=1)

# drop non-percentage features: (total_households, total_individuals, lat, lon, NL and all-Zero values)
df = df.drop(['total_households', 'total_individuals', 'lat', 'lon', 'NL', 'dw_13', 'lan_13', 'dw_12', 'pw_08', 'pw_07'], axis=1)
# -

# check new dataframe after dropping the features with only 0 in them
df.shape

# + colab={"base_uri": "https://localhost:8080/", "height": 1000} colab_type="code" id="ytHU7Sh2KYOc" outputId="4233dd3a-0f6b-41e9-d494-2b3d98f883bf"
# load the description file and have a look at them
var_desc = pd.read_csv('data/variable_descriptions.csv')
pd.set_option('display.max_colwidth', 200) # So that we can see the full descriptions
var_desc

# + [markdown] colab_type="text" id="5chJlQBxKYOi"
# # Quick EDA
# Before we dive into modelling, let's take a quick look at the data and see what we have to work with. 

# + [markdown] colab_type="text" id="kO6UB7uQKYOj"
# ## Looking at the target variable
# The purpose of this project is to predict the percentage of households in a ward that are woman-headed and have an annual household income of <R19,600. To put that in context, that's less than 250USD per month. 

# + colab={"base_uri": "https://localhost:8080/", "height": 623} colab_type="code" id="fhCt8mEIKYOl" outputId="9fb84b2e-74d7-42fd-c87b-b8099e2798d8"
# Plot the distribution
ax = df['target'].plot(kind='hist')
plt.xlabel('Percentage of households that are w-headed and have income under R19.6k/month');

# + [markdown] colab_type="text" id="aBJaE--wKYOJ"
# # Baseline Model
#
# As a starter we will create a baseline model to get an idea how a more complex model could look like.
# -

# Get a very reduced dataframe for the baseline model
# we choose the highly correlated variable school attendance = yes (psa_00) for our baseline model 
df_base = df[['psa_00', 'target']]
df_base.head()

# plot the baseline model
df_base.plot(x='psa_00', y='target', kind='scatter', alpha=0.3);

# + [markdown] colab_type="text" id="88Ek8kY9KYPg"
# ## Modelling the baseline model
# Before we go deeper, we will build a quick model FIRST, to see a baseline to compare to. So, in this section we'll try a simple linear model based on only one highly correlated feature. 

# +
# Modelling by hand: regression line formula

# as we can see from the scatter plot, a line would run through (0 | -5), we only need a slope
# points would be: (0.2 | 10), (0.3 | 20) --> the slope is 100
# y = 100 * x -5

# put this line in the plot:
x = df_base['psa_00']
y = 100 * x - 5

df_base.plot(x='psa_00', y='target', kind='scatter', alpha=0.3)
plt.plot(x, y, '-r', label='y = 100 * x')
plt.show();

# + [markdown] colab_type="text" id="wL09aomvQ-eC"
# ### Evaluating Predictions
#
# We can evaluate our predictions with the "R-Squared Score", which is a measure of how well our model explains the variation in the target variable. The closer to 1, the better. 

# +
# try to calculate RMSE for this line:
from sklearn.metrics import mean_squared_error
import math

mse = mean_squared_error(df_base['target'], y)

rmse = math.sqrt(mse)

print("RMSE of baseline model:", round(rmse, 2))
print("Error in % of baseline model:", round((rmse/np.average(y))*100, 2))
print("Mean of the target:", round(np.average(y), 2))
print("Standard deviation of the target:", round(np.std(y), 2))
# -

# # Complex model

# In order to model the complex model we consider all the remaining 51 features,

# +
# Train test split with same random seed
# Defining X and y
features = df.columns.tolist()
features.remove('target')
X = df[features]
y = df.target

# Splitting the dataset
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=150, shuffle=True)

# Check the shape of the data sets
print("X_train:", X_train.shape)
print("y_train:", y_train.shape)
print("X_test:", X_test.shape)
print("y_test:", y_test.shape)
# -

# train model
lr = LinearRegression()
lr.fit(X_train,y_train)

# predict target values
y_pred = lr.predict(X_test)

# check error for predictions
from sklearn.metrics import r2_score
score = r2_score(y_test, y_pred)
print("r2 score is ", round(score, 2))
print("mean_sqrd_error is:", round(mean_squared_error(y_test, y_pred), 2))
print("root_mean_squared error is:", round(np.sqrt(mean_squared_error(y_test, y_pred)), 2))

# # Final Model

# After trying a first very basic Baseline model and then the complex model, we now make our model simpler by reducing the features. The features were chosen after checking the correlations.

# + [markdown] colab_type="text" id="Y5898TVfKYOs"
# # Looking at a few interesting correlations

# + colab={"base_uri": "https://localhost:8080/", "height": 204} colab_type="code" id="fVMaahPiKYOv" outputId="5bdd7b41-98f9-423e-ee7d-9c849032035d"
# Strong +ve correlations
df.corr()['target'].sort_values(ascending=False).head(10)

# + colab={"base_uri": "https://localhost:8080/", "height": 204} colab_type="code" id="d61CC8IlKYO4" outputId="86e36bad-b955-4045-b6a7-e160b674fdab"
# And negative correlations
df.corr()['target'].sort_values().head(10)
# -

# create a list of positive and negative correlated features above a treshold of |0.5|
# negatively correlated features (- 0.5 threshold)
list1 = list((df.corr()['target'].sort_values().head(8)).index)
# positively correlated features (+ 0.5 threshold)
list2 = list((df.corr()['target'].sort_values(ascending=False).head(6)).index)
# list2 = list((df.corr()['target'].sort_values(ascending=False).head(6)).index)
# combine both lists
corr_feature_list = list1 + list2
corr_feature_list

# plot the list of highly correlated features
corr = df[corr_feature_list].corr()
matrix = np.triu(corr)
sns.heatmap(corr, annot=True, mask=matrix);

# Car00 and Car01 / Stv00 and Stv01 / Lln00 and Lln01 perfectly correlated (- 1) --> we can remove each one of them

# select features with high correlations
#feature_list = ['pw_00', 'pw_01', 'psa_00', 'psa_01', 'car_00', 'lln_00', 'stv_00', 'pg_00', 'pg_03', 'lan_00', 'lan_01', 'target']
df_final = df[corr_feature_list]
df_final.drop(['car_01', 'stv_01', 'lln_01'], axis=1, inplace=True)
df_final.head()


# ## Transformation

# write function to transform features with boxcox
def convert_zeros(x):
    '''
    function to convert zeros to a postive number 
    so that it can be transformed with the boxcox'''
    if x == 0.0:
        return 0.0000001
    else :
        return x


# check for 0 before and convert them
print (df['pw_00'].isin([0]).sum())
print(df['pw_00'].apply(convert_zeros).isin([0]).sum())

# plot transformed features
for col in df_final.columns:
    if col != 'target':
        fig, axes = plt.subplots(1,2)
        df_final[col] = df_final[col].apply(convert_zeros)
        sns.histplot(df_final[col], ax=axes[0])
        plt.xlabel(col)
        sns.histplot(boxcox(df_final[col]), ax=axes[1])
        plt.xlabel('boxcox_transformed  ' + col)

# +
# apply the boxcox transformation on 

for col in df_final.columns:
    if col != 'target' and  col != 'psa_00' and  col != 'psa_01' and col != 'car_00' and col != 'pg_00':
        #df_final[col] = df_final[col].apply(convert_zeros)
        df_final[col] = boxcox(df_final[col])[0].reshape(-1,1);
# -

# check data after log transformation
df_final.head()

# +
# do the cross validation manually
from sklearn.model_selection import KFold

# Using this to test a model on 5 different splits
kf = KFold(n_splits=5, shuffle=False)

ycol = 'target'
in_cols = df_final.drop('target', axis=1).columns

scores = []
for train, test in kf.split(df_final):
  lr = LinearRegression()
  lr.fit(df_final[in_cols].iloc[train], df_final[ycol].iloc[train])
  rmse = np.sqrt(mean_squared_error(df_final[ycol].iloc[test], lr.predict(df_final[in_cols].iloc[test])))
  scores.append(rmse)
  print(rmse)

print("Average score in 5-fold CV:", np.mean(scores))

# +
# Train test split with same random seed
# Defining X and y
features = df_final.columns.tolist()
features.remove('target')
X = df_final[features]
y = df_final.target

# Splitting the dataset
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=150, shuffle=True)

# Check the shape of the data sets
print("X_train:", X_train.shape)
print("y_train:", y_train.shape)
print("X_test:", X_test.shape)
print("y_test:", y_test.shape)

# +
lr = LinearRegression()
lr.fit(X_train,y_train)

# predict target values
y_pred = lr.predict(X_test)
# -

# check error for predictions
from sklearn.metrics import r2_score
score = r2_score(y_test, y_pred)
print("r2 score is ", round(score, 2))
print("mean_sqrd_error is:", round(mean_squared_error(y_test, y_pred), 2))
print("root_mean_squared error is:", round(np.sqrt(mean_squared_error(y_test, y_pred)), 2))

# +
# Plot the errors
from yellowbrick.regressor import ResidualsPlot

#visualizer = ResidualsPlot(gs, hist=False, qqplot=True)
visualizer = ResidualsPlot(gs)

visualizer.fit(X_train, y_train)  # Fit the training data to the visualizer
visualizer.score(X_test, y_test)  # Evaluate the model on the test data
visualizer.show()
# -

#



StarcoderdataJupyterScriptsDedupFiltered