Regularized Linear Regression with Model Stacking

# import packages

# basic packages
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import random

# general packages
import warnings
warnings.filterwarnings('ignore')
%matplotlib inline

# model packages
from sklearn import linear_model
from sklearn.linear_model import ElasticNet, Lasso,  Ridge, LassoLarsIC, LinearRegression, LogisticRegression
from sklearn.ensemble import RandomForestRegressor,  GradientBoostingRegressor
from sklearn.metrics import mean_squared_error
import xgboost
from sklearn.model_selection import train_test_split

# Model Tune packages
from sklearn.model_selection import GridSearchCV

# CV Folder
from sklearn.model_selection import cross_validate
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import ShuffleSplit

# load train/test data
train = pd.read_csv('../input/train.csv')
test = pd.read_csv('../input/test.csv')
df = pd.concat((train, test)).reset_index(drop=True)

# basic info
print('Shape of Train: ',train.shape)
print('-'*40)
print('Shape of Test: ',test.shape)
print('-'*40)
print(train.columns)
print('-'*40)
print(test.columns)

# saleprice
g = sns.distplot(train.SalePrice,color = 'g')
g.set_title('Sale Price Distribution')
plt.show()

# scale sale price
df['SalePrice'] = np.log1p(df['SalePrice'])

# plot and cjeck
g1 = sns.distplot(df[:1460].SalePrice,color = 'g')
g1.set_title('Scaled Sale Price Distribution')

plt.show()

# skewness and kurtosis
print('Skewness is: ', train.SalePrice.skew())
print('-'*40)
print('Kurtosis is: ',train.SalePrice.kurt())

# Using Corraltion matrix to analyze numeric data first
_,ax = plt.subplots(figsize=(14,12))
sns.heatmap(train.drop('Id',axis = 1).corr(),cmap = 'Greens')
ax.set_title('SalePrice vs All Numeric Variables \n The darker color, the higher relationship.')
plt.show()

# list relationship of all numeric variables with SalePrice
_,ax = plt.subplots(figsize = (14,12))
g = train.drop(['Id'],axis = 1).corr()['SalePrice'].sort_values().plot.barh(color = 'g')
ax.set_title('Relationship of Survived vs All Numeric Variables')

# add annotation
for p in g.patches:
    g.annotate(str(round(p.get_width(),2)), (p.get_width() * 1.01,p.get_y()))

plt.show()

# plot all missing values
_,ax = plt.subplots(1,2,figsize = (14,8))
g1 = df.isnull().sum().sort_values(ascending = False).head(30).plot.barh(ax = ax[0],color = 'g')
g2 = train.isnull().sum().sort_values(ascending = False).head(20).plot.barh(ax = ax[1],color = 'g')
g1.set_title('Missing Values for All Data')
g2.set_title('Missing Values fro Train')

# add annotation
for p in g1.patches:
    g1.annotate(str(round(p.get_width(),2)), (p.get_width() * 1.01,p.get_y()))

plt.show()

# drop columns with more than 600 missings values, eg ('PoolQC','MiscFeature','Alley','Fence','FireplaceQu')
df = df.drop(['PoolQC','MiscFeature','Alley','Fence','FireplaceQu'],axis = 1)
print(df.shape)

# LotFrontage : immutate with neighborhood variables
df['LotFrontage'] = df.groupby(['Neighborhood'])['LotFrontage'].transform(lambda x: x.fillna(x.median()))

for col in ('GarageType', 'GarageFinish', 'GarageQual', 'GarageCond'):
    df[col] = df[col].fillna('None')
for col in ('GarageYrBlt', 'GarageArea', 'GarageCars'):
    df[col] = df[col].fillna(0)
for col in ('BsmtFinSF1', 'BsmtFinSF2', 'BsmtUnfSF','TotalBsmtSF', 'BsmtFullBath', 'BsmtHalfBath'):
    df[col] = df[col].fillna(0)
for col in ('BsmtQual', 'BsmtCond', 'BsmtExposure', 'BsmtFinType1', 'BsmtFinType2'):
    df[col] = df[col].fillna('None')
df["MasVnrType"] = df["MasVnrType"].fillna("None")
df["MasVnrArea"] = df["MasVnrArea"].fillna(0)

# MSZoning
print(df['MSZoning'].describe()) # top is RL
df['MSZoning'] = df['MSZoning'].fillna('RL')

# Utilities
df['Utilities'].value_counts() # there are 2916 allpub and 1 Nosewa and 2 NAs, drop this variable
df = df.drop('Utilities',axis = 1)

# Functional
print(df['Functional'].value_counts()) # top is Typ
df['Functional'] = df['Functional'].fillna('Typ')

# Exterior1st / Exterior2nd
print(df['Exterior1st'].value_counts()) # top is VinylSd and ImStucc only has 1 row, can be deleted when get_dummies
print(df['Exterior2nd'].value_counts()) # top is VinylSd and Other only has 1 row, can be deleted when get_dummies
df['Exterior1st'] = df['Exterior1st'].fillna('VinylSd')
df['Exterior2nd'] = df['Exterior2nd'].fillna('VinylSd')

# KitchenQual
print(df['KitchenQual'].value_counts())
df['KitchenQual'] = df['KitchenQual'].fillna('TA')

# SaleType
print(df['SaleType'].value_counts())
df['SaleType'] = df['SaleType'].fillna('WD')

# Electrical
print(df['Electrical'].value_counts()) # top is SBrkr and Mix has only 1 row, deletes when get_dummies
df['Electrical'] = df['Electrical'].fillna('SBrkr')

_,ax = plt.subplots(figsize = (8,6))
g1 = df.isnull().sum().sort_values(ascending = False).head(30).plot.barh(color = 'g')
g1.set_title('Missing Values for All Data')
for p in g1.patches:
    g1.annotate(str(round(p.get_width(),2)), (p.get_width() * 1.01,p.get_y()))   
plt.show()

_,ax = plt.subplots(3,2,figsize = (14,14))
sns.regplot(x = df.GrLivArea,y = df.SalePrice,ax = ax[0,0])
sns.regplot(x = df.TotalBsmtSF,y = df.SalePrice,ax = ax[0,1])
sns.boxplot(x = df.OverallCond,y = df.SalePrice,ax = ax[1,0])
sns.boxplot(x = df.OverallQual,y = df.SalePrice,ax = ax[1,1])
sns.regplot(x = df['1stFlrSF'],y = df.SalePrice,ax = ax[2,0])
sns.regplot(x = df['2ndFlrSF'],y = df.SalePrice, ax = ax[2,1])

_,ax = plt.subplots(2,2,figsize = (14,10))
sns.boxplot(y = df.GarageArea,x = df.GarageCars,ax = ax[1,0])
sns.regplot(x = df['GarageArea'],y = df.SalePrice, ax = ax[0,1])
sns.boxplot(y = df.SalePrice,x = df.GarageCars,ax = ax[0,0])
sns.regplot(x = df['GarageArea']/df.GarageCars,y = df.SalePrice, ax = ax[1,1])

ax[0,0].set_title('GarageCars vs SalePrice')
ax[0,1].set_title('GarageAarea vs SalePrice')
ax[1,0].set_title('GarageCars vs GarageArea')
ax[1,1].set_title('GarageCars vs GarageArea Per Cars')

# drop outliers from GrLivArea and TotalBsmtSF plot
df = df.drop(df[df['Id'] == 1299].index)
df = df.drop(df[df['Id'] == 524].index)
# watch out there is one extremely point in test data too.

scaled = ["1stFlrSF","2ndFlrSF","3SsnPorch","BsmtFinSF1","BsmtFinSF2","BsmtUnfSF","EnclosedPorch",
          "GarageArea","GrLivArea","LotArea","LotFrontage","LowQualFinSF","MasVnrArea","MiscVal",
          "OpenPorchSF","PoolArea","ScreenPorch","TotalBsmtSF", "WoodDeckSF"]

# ckech kurtosis and skewness
def skew_kur(df) :
    # select numeric variables
    num = []
    for c in df.columns :
        if c in scaled:
            num.append(c)
    # skewness
    skewness = []
    for c in num :
        skewness.append(df[c].skew())
    # kurtosis  
    kurtosis = []
    for c in num :
        kurtosis.append(df[c].kurtosis())
    # norm dataframe
    norm = pd.DataFrame({'Variable' : num,
                        'skewness' : skewness,
                        'kurtosis' : kurtosis})
    return norm

norm = skew_kur(df)
skew_kur(df)

# consider -0.8 to 0.8 for skewness and -3.0 to 3.0 for kurtosis as the acceptable ranege to ensure 95% CI.
temp = norm[(abs(norm['kurtosis']) > 3) | (abs(norm['skewness']) > 0.8)]['Variable']

df[temp] = df[temp].apply(lambda x: np.log1p(x))

# double check
skew_kur(df)
# there are still some extremely value but much better

# HasPool : if poolarea != 0 then haspool = 1
df['HasPool'] = 0
df[df.PoolArea != 0]['HasPool'] = 1

# HasGarage : if garageyrblt =0 then hasgarage = 0
df['HasGarage'] = 0
df[df.GarageYrBlt != 0]['HasGarage'] = 1

# Month
_,ax = plt.subplots(1,2,figsize = (14,6))
sns.countplot(df.MoSold,ax =ax[0])
sns.boxplot(x = df.MoSold,y = df.SalePrice,ax = ax[1])

# Built year
_,ax = plt.subplots(3,1,figsize = (14,10))
sns.regplot(x = df.YearBuilt.astype(int),y = df.SalePrice,ax = ax[0])
sns.boxplot(x = df.YearBuilt.astype(int),y = df.SalePrice,ax =ax[1])
sns.countplot(df.YearBuilt.astype(int),ax =ax[2])

_,ax = plt.subplots(1,2,figsize = (14,6))
sns.countplot(df.YrSold,ax = ax[0])
sns.regplot(x = df.YrSold.astype(int),y = df.SalePrice,ax = ax[1] )

# Hot Mon
df['HotMon'] = 0
df[df.MoSold == 5|6|7]['HotMon'] =1

# HouseAge : YearSold - YearBuilt
df['HouseAge'] = df.YrSold.astype(int) - df.YearBuilt.astype(int)

# Neighborhood
_,ax = plt.subplots(2,1,figsize = (14,8))
sns.countplot(df.Neighborhood,ax = ax[0])
sns.boxplot(x = df.Neighborhood,y = df.SalePrice,ax = ax[1])

# consider average of each Neighborhood > 12.5 as 'ClassI' and others as 'Others'
print(df.groupby('Neighborhood')['SalePrice'].mean().sort_values(ascending = False).head())

# classI
df['ClassI'] = 0
df[(df.Neighborhood == 'NoRidge')|(df.Neighborhood =='NridgHt')|(df.Neighborhood =='StoneBr')]['ClassI'] = 1

# TotalArea = basement + 1stfloor + 2ndfloor
df['TotalArea'] = df['1stFlrSF'] +df['2ndFlrSF'] +df['TotalBsmtSF']

# street tp 0/1
df.Street.replace(['Grvl', 'Pave'],[0,1],inplace = True)

df_dummy = pd.get_dummies(df)
print(df_dummy.shape)

# drop some columns
# Exterior1st :ImStucc
# Exterior2nd :Other
# Electrical :Mix
df_dummy.drop(['Exterior1st_ImStucc','Exterior2nd_Other','Electrical_Mix'],axis = 1,inplace = True)
df_dummy.shape

# seperate train and test
train_dummy = df_dummy[:train.shape[0]-2] # 2 outliers
test_dummy = df_dummy[train.shape[0]-2:]
print(train_dummy.shape,test_dummy.shape)

# set up train/ test /traget / data dataset
train_Price = train_dummy['SalePrice'] # target
train_data = train_dummy.drop(['SalePrice','Id'],axis = 1) # data
print(train_Price.shape,train_data.shape)

# train test split for scoring
x_train,x_test,y_train,y_test = train_test_split(train_data,train_Price,
                                               test_size = 0.25,random_state = 13)
print(x_train.info())

# set up dataframe to see rmse
models = pd.DataFrame({
    'Model': ['Linear', 'Lasso', 'Ridge',
              'Elastic Net', 'Gradient Boosting', 'Random Forest',
              'Xgboost'],
    'Result': [ np.sqrt(mean_squared_error(LinearRegression().fit(x_train,y_train).predict(x_test),y_test)),
              np.sqrt(mean_squared_error(Lasso().fit(x_train,y_train).predict(x_test),y_test)),
              np.sqrt(mean_squared_error(Ridge().fit(x_train,y_train).predict(x_test),y_test)),
              np.sqrt(mean_squared_error(ElasticNet().fit(x_train,y_train).predict(x_test),y_test)),
              np.sqrt(mean_squared_error(GradientBoostingRegressor().fit(x_train,y_train).predict(x_test),y_test)),
              np.sqrt(mean_squared_error(RandomForestRegressor().fit(x_train,y_train).predict(x_test),y_test)),
              np.sqrt(mean_squared_error(xgboost.XGBRegressor().fit(x_train,y_train).predict(x_test),y_test))
              ]})
print(models.sort_values(by='Result'))
g = models.plot(kind = 'bar',title = 'Models RMSE Score \n The Lower the Better')
g.set_xlabel(list(models['Model']))

# split cv folder
cv_split = ShuffleSplit(n_splits = 10, test_size = .3,
                                        train_size = .6, random_state = 13)

# models will test
mod = [LinearRegression(),Lasso(),Ridge(),ElasticNet(),GradientBoostingRegressor(),RandomForestRegressor(),xgboost.XGBRegressor()]

# Cross validation
cv_score = list()
for model in mod :
    cv_result = np.sqrt(-cross_val_score(model, train_data, train_Price, cv  = cv_split,scoring='mean_squared_error'))
    cv_score.append(cv_result.mean())

cv_model = pd.DataFrame({
        'Model': ['Linear', 'Lasso', 'Ridge',
              'Elastic Net', 'Gradient Boosting', 'Random Forest',
              'Xgboost'],
        'CV Result' :cv_score})

print(cv_model.sort_values(by='CV Result'))

# Lasso
Lasso().get_params()
param_grid = {'alpha' : [0.0005,0.005,0.05,0.2,0.4],
              'tol': [0.1,0.01,0.001, 0.0001],
              'random_state' : [13]}
tuned_lasso = GridSearchCV(Lasso(),param_grid=param_grid, scoring = 'neg_mean_squared_error', cv = cv_split)
tuned_lasso.fit(train_data,train_Price)
print("Best Hyper Parameters:\n",tuned_lasso.best_params_)

# Ridge
Ridge().get_params()
param_grid = {'alpha' : [0.005,0.05,0.2,0.4,0.6,1],
              'tol': [0.1,0.01,0.001, 0.0001],
              'random_state' : [13]}
tuned_ridge = GridSearchCV(Ridge(),param_grid=param_grid, scoring = 'neg_mean_squared_error', cv = cv_split)
tuned_ridge.fit(train_data,train_Price)
print("Best Hyper Parameters:\n",tuned_ridge.best_params_)

# ElasticNet
ElasticNet().get_params()
param_grid = {'alpha' : [0.1,1],
              'l1_ratio' :[0,0.4,0.7,1],
              'tol': [0.1,0.01],
              'random_state' : [13]}
tuned_elnet = GridSearchCV(ElasticNet(),param_grid=param_grid, scoring = 'neg_mean_squared_error', cv = cv_split)
tuned_elnet.fit(train_data,train_Price)
print("Best Hyper Parameters:\n",tuned_ridge.best_params_)

# GradientBoostingRegressor
GradientBoostingRegressor().get_params()
param_grid = {'learning_rate': [0.01],
              'min_samples_split':[5],
              'min_samples_leaf':[5],
              'max_depth':[3],
              'max_features':['sqrt'],
              'subsample':[.5,],
              'n_estimators':[3000,4000],
              'random_state':[13]}
tuned_gbm = GridSearchCV(GradientBoostingRegressor(),param_grid=param_grid, scoring = 'neg_mean_squared_error', cv = cv_split)
tuned_gbm.fit(train_data,train_Price)
print("Best Hyper Parameters:\n",tuned_gbm.best_params_)

# RandomForestRegressor
RandomForestRegressor().get_params()
param_grid = {'n_estimators': [2000,3000],
              'max_features': ['sqrt'],
              'criterion' : ['mse'],
              'n_jobs' : [-1],
              'random_state' : [13]}
tuned_rf = GridSearchCV(RandomForestRegressor(),param_grid=param_grid, scoring = 'neg_mean_squared_error', cv = cv_split)
tuned_rf.fit(train_data,train_Price)
print("Best Hyper Parameters:\n",tuned_rf.best_params_)

# xgboost.XGBRegressor
xgboost.XGBRegressor().get_params()
param_grid = {"colsample_bytree" : [0.4],
              "gamma" : [0.4],
              "learning_rate" : [0.05],
              "max_depth" : [5],              
              "min_child_weight":[2],
              "n_estimators" : [2500],
              "subsample" : [0.6],
              "random_state" :[13],
              "nthread" : [-1]}
tuned_xgb = GridSearchCV(xgboost.XGBRegressor(),param_grid=param_grid, scoring = 'neg_mean_squared_error', cv = cv_split)
tuned_xgb.fit(train_data,train_Price)
print("Best Hyper Parameters:\n",tuned_xgb.best_params_)

# model will tested
ln_mod = LinearRegression()
las_mod = Lasso(alpha = 0.0005,tol = 0.01,random_state = 13)
rid_mod = Ridge(alpha = 0.1,tol = 0.1,random_state = 13)
enet_mod = ElasticNet(alpha = 0.1,tol = 0.1,random_state = 13)
gbm_mod = GradientBoostingRegressor(max_features = 'sqrt',n_estimators = 4000,random_state = 13,
                                learning_rate = 0.01,max_depth = 3,min_samples_leaf = 5,min_samples_split = 5,subsample = 0.5)
rf_mod = RandomForestRegressor(criterion = 'mse',max_features = 'sqrt',n_estimators = 3000,
                                    min_samples_split = 8,random_state = 13)
xgb_mod = xgboost.XGBRegressor(colsample_bytree = .4,max_depth = 5,learning_rate = 0.05,
                                     n_estimators = 2500,gamma = 0.4,
                                     random_state = 13,subsample = .6, min_child_weight = 2)

mod_tuned = [ln_mod,las_mod,rid_mod,enet_mod,gbm_mod,rf_mod,xgb_mod]

# Cross validation
cv_score = list()
for model in mod_tuned :
    cv_result = np.sqrt(-cross_val_score(model, train_data, train_Price, cv  = cv_split,scoring='mean_squared_error'))
    cv_score.append(cv_result.mean())

tuned_cv_model = pd.DataFrame({
                                'Model': ['Linear','Lasso', 'Ridge',
                                      'Elastic Net', 'Gradient Boosting', 'Random Forest',
                                      'Xgboost'],
                                'Tuned CV Result' :cv_score})
print(tuned_cv_model)

_,ax = plt.subplots(figsize=(8,6))
g = tuned_cv_model.plot(ax = ax)
cv_model.plot(ax = g)
models.plot(ax = g)
g.set_title('Tuned Result vs CV result vs Original Result')
g.set_xlabel(list(models['Model']))

# Lasso Ridge GBM
# m1 = 0.6*lasso + 0.35*gbm + 0.05*ridge
m1 = 0.6*las_mod.fit(x_train,y_train).predict(x_test) + 0.35*gbm_mod.fit(x_train,y_train).predict(x_test) + 0.05*rid_mod.fit(x_train,y_train).predict(x_test)
r = np.sqrt(mean_squared_error(m1,y_test))
print('Combined Model RMSE is : ',r)

# combined results from 6 single models together
train_index = random.sample(range(1,len(train_data)),round(0.6*len(train_data)))
test_index = random.sample(range(1,len(test)),round(0.6*len(test)))
t_train_stacking = np.vstack(
    (las_mod.fit(train_data.iloc[train_index],train_Price.iloc[train_index]).predict(train_data),
     rid_mod.fit(train_data.iloc[train_index],train_Price.iloc[train_index]).predict(train_data),
     enet_mod.fit(train_data.iloc[train_index],train_Price.iloc[train_index]).predict(train_data),
     gbm_mod.fit(train_data.iloc[train_index],train_Price.iloc[train_index]).predict(train_data),
     rf_mod.fit(train_data.iloc[train_index],train_Price.iloc[train_index]).predict(train_data),
     xgb_mod.fit(train_data.iloc[train_index],train_Price.iloc[train_index]).predict(train_data)))
train_stacking = pd.DataFrame(t_train_stacking.T)
train_stacking.shape

# split into train and test part
x_train_stacking,x_test_stacking,y_train_stacking,y_test_stacking = train_test_split(train_stacking,train_Price,
                                                                                    test_size = .3,random_state = 13)

# Using xgb as Level 2 model
l2_mod.fit(x_train_stacking,y_train_stacking)
predictions = l2_mod.predict(x_test_stacking)
print('Level 2 Modeling RMSE Score:',np.sqrt(mean_squared_error(predictions,y_test_stacking)))