CASE STUDY 4: House Price Prediction

HOUSE PRICE PREDICTION

We have a dataset related to house prices and attributes of each house sale. It is a kaggle competition dataset.

In univariate analysis, we will try to understand the overall structure of data, each variable and some univariate hypotheses analysis.

First, let us understand the overall structure of the data set by running following chunk of commands—

setwd("E:/r/R/kaggle/houseprice")

getwd()
## [1] "E:/r/R/kaggle/houseprice"
train=read.csv("train.csv", header = T, stringsAsFactors = TRUE)

str(train)
## 'data.frame':    1460 obs. of  81 variables:
##  $ Id           : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ MSSubClass   : int  60 20 60 70 60 50 20 60 50 190 ...
##  $ MSZoning     : Factor w/ 5 levels "C (all)","FV",..: 4 4 4 4 4 4 4 4 5 4 ...
##  $ LotFrontage  : int  65 80 68 60 84 85 75 NA 51 50 ...
##  $ LotArea      : int  8450 9600 11250 9550 14260 14115 10084 10382 6120 7420 ...
##  $ Street       : Factor w/ 2 levels "Grvl","Pave": 2 2 2 2 2 2 2 2 2 2 ...
##  $ Alley        : Factor w/ 2 levels "Grvl","Pave": NA NA NA NA NA NA NA NA NA NA ...
##  $ LotShape     : Factor w/ 4 levels "IR1","IR2","IR3",..: 4 4 1 1 1 1 4 1 4 4 ...
##  $ LandContour  : Factor w/ 4 levels "Bnk","HLS","Low",..: 4 4 4 4 4 4 4 4 4 4 ...
##  $ Utilities    : Factor w/ 2 levels "AllPub","NoSeWa": 1 1 1 1 1 1 1 1 1 1 ...
##  $ LotConfig    : Factor w/ 5 levels "Corner","CulDSac",..: 5 3 5 1 3 5 5 1 5 1 ...
##  $ LandSlope    : Factor w/ 3 levels "Gtl","Mod","Sev": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Neighborhood : Factor w/ 25 levels "Blmngtn","Blueste",..: 6 25 6 7 14 12 21 17 18 4 ...
##  $ Condition1   : Factor w/ 9 levels "Artery","Feedr",..: 3 2 3 3 3 3 3 5 1 1 ...
##  $ Condition2   : Factor w/ 8 levels "Artery","Feedr",..: 3 3 3 3 3 3 3 3 3 1 ...
##  $ BldgType     : Factor w/ 5 levels "1Fam","2fmCon",..: 1 1 1 1 1 1 1 1 1 2 ...
##  $ HouseStyle   : Factor w/ 8 levels "1.5Fin","1.5Unf",..: 6 3 6 6 6 1 3 6 1 2 ...
##  $ OverallQual  : int  7 6 7 7 8 5 8 7 7 5 ...
##  $ OverallCond  : int  5 8 5 5 5 5 5 6 5 6 ...
##  $ YearBuilt    : int  2003 1976 2001 1915 2000 1993 2004 1973 1931 1939 ...
##  $ YearRemodAdd : int  2003 1976 2002 1970 2000 1995 2005 1973 1950 1950 ...
##  $ RoofStyle    : Factor w/ 6 levels "Flat","Gable",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ RoofMatl     : Factor w/ 8 levels "ClyTile","CompShg",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ Exterior1st  : Factor w/ 15 levels "AsbShng","AsphShn",..: 13 9 13 14 13 13 13 7 4 9 ...
##  $ Exterior2nd  : Factor w/ 16 levels "AsbShng","AsphShn",..: 14 9 14 16 14 14 14 7 16 9 ...
##  $ MasVnrType   : Factor w/ 4 levels "BrkCmn","BrkFace",..: 2 3 2 3 2 3 4 4 3 3 ...
##  $ MasVnrArea   : int  196 0 162 0 350 0 186 240 0 0 ...
##  $ ExterQual    : Factor w/ 4 levels "Ex","Fa","Gd",..: 3 4 3 4 3 4 3 4 4 4 ...
##  $ ExterCond    : Factor w/ 5 levels "Ex","Fa","Gd",..: 5 5 5 5 5 5 5 5 5 5 ...
##  $ Foundation   : Factor w/ 6 levels "BrkTil","CBlock",..: 3 2 3 1 3 6 3 2 1 1 ...
##  $ BsmtQual     : Factor w/ 4 levels "Ex","Fa","Gd",..: 3 3 3 4 3 3 1 3 4 4 ...
##  $ BsmtCond     : Factor w/ 4 levels "Fa","Gd","Po",..: 4 4 4 2 4 4 4 4 4 4 ...
##  $ BsmtExposure : Factor w/ 4 levels "Av","Gd","Mn",..: 4 2 3 4 1 4 1 3 4 4 ...
##  $ BsmtFinType1 : Factor w/ 6 levels "ALQ","BLQ","GLQ",..: 3 1 3 1 3 3 3 1 6 3 ...
##  $ BsmtFinSF1   : int  706 978 486 216 655 732 1369 859 0 851 ...
##  $ BsmtFinType2 : Factor w/ 6 levels "ALQ","BLQ","GLQ",..: 6 6 6 6 6 6 6 2 6 6 ...
##  $ BsmtFinSF2   : int  0 0 0 0 0 0 0 32 0 0 ...
##  $ BsmtUnfSF    : int  150 284 434 540 490 64 317 216 952 140 ...
##  $ TotalBsmtSF  : int  856 1262 920 756 1145 796 1686 1107 952 991 ...
##  $ Heating      : Factor w/ 6 levels "Floor","GasA",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ HeatingQC    : Factor w/ 5 levels "Ex","Fa","Gd",..: 1 1 1 3 1 1 1 1 3 1 ...
##  $ CentralAir   : Factor w/ 2 levels "N","Y": 2 2 2 2 2 2 2 2 2 2 ...
##  $ Electrical   : Factor w/ 5 levels "FuseA","FuseF",..: 5 5 5 5 5 5 5 5 2 5 ...
##  $ X1stFlrSF    : int  856 1262 920 961 1145 796 1694 1107 1022 1077 ...
##  $ X2ndFlrSF    : int  854 0 866 756 1053 566 0 983 752 0 ...
##  $ LowQualFinSF : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ GrLivArea    : int  1710 1262 1786 1717 2198 1362 1694 2090 1774 1077 ...
##  $ BsmtFullBath : int  1 0 1 1 1 1 1 1 0 1 ...
##  $ BsmtHalfBath : int  0 1 0 0 0 0 0 0 0 0 ...
##  $ FullBath     : int  2 2 2 1 2 1 2 2 2 1 ...
##  $ HalfBath     : int  1 0 1 0 1 1 0 1 0 0 ...
##  $ BedroomAbvGr : int  3 3 3 3 4 1 3 3 2 2 ...
##  $ KitchenAbvGr : int  1 1 1 1 1 1 1 1 2 2 ...
##  $ KitchenQual  : Factor w/ 4 levels "Ex","Fa","Gd",..: 3 4 3 3 3 4 3 4 4 4 ...
##  $ TotRmsAbvGrd : int  8 6 6 7 9 5 7 7 8 5 ...
##  $ Functional   : Factor w/ 7 levels "Maj1","Maj2",..: 7 7 7 7 7 7 7 7 3 7 ...
##  $ Fireplaces   : int  0 1 1 1 1 0 1 2 2 2 ...
##  $ FireplaceQu  : Factor w/ 5 levels "Ex","Fa","Gd",..: NA 5 5 3 5 NA 3 5 5 5 ...
##  $ GarageType   : Factor w/ 6 levels "2Types","Attchd",..: 2 2 2 6 2 2 2 2 6 2 ...
##  $ GarageYrBlt  : int  2003 1976 2001 1998 2000 1993 2004 1973 1931 1939 ...
##  $ GarageFinish : Factor w/ 3 levels "Fin","RFn","Unf": 2 2 2 3 2 3 2 2 3 2 ...
##  $ GarageCars   : int  2 2 2 3 3 2 2 2 2 1 ...
##  $ GarageArea   : int  548 460 608 642 836 480 636 484 468 205 ...
##  $ GarageQual   : Factor w/ 5 levels "Ex","Fa","Gd",..: 5 5 5 5 5 5 5 5 2 3 ...
##  $ GarageCond   : Factor w/ 5 levels "Ex","Fa","Gd",..: 5 5 5 5 5 5 5 5 5 5 ...
##  $ PavedDrive   : Factor w/ 3 levels "N","P","Y": 3 3 3 3 3 3 3 3 3 3 ...
##  $ WoodDeckSF   : int  0 298 0 0 192 40 255 235 90 0 ...
##  $ OpenPorchSF  : int  61 0 42 35 84 30 57 204 0 4 ...
##  $ EnclosedPorch: int  0 0 0 272 0 0 0 228 205 0 ...
##  $ X3SsnPorch   : int  0 0 0 0 0 320 0 0 0 0 ...
##  $ ScreenPorch  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ PoolArea     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ PoolQC       : Factor w/ 3 levels "Ex","Fa","Gd": NA NA NA NA NA NA NA NA NA NA ...
##  $ Fence        : Factor w/ 4 levels "GdPrv","GdWo",..: NA NA NA NA NA 3 NA NA NA NA ...
##  $ MiscFeature  : Factor w/ 4 levels "Gar2","Othr",..: NA NA NA NA NA 3 NA 3 NA NA ...
##  $ MiscVal      : int  0 0 0 0 0 700 0 350 0 0 ...
##  $ MoSold       : int  2 5 9 2 12 10 8 11 4 1 ...
##  $ YrSold       : int  2008 2007 2008 2006 2008 2009 2007 2009 2008 2008 ...
##  $ SaleType     : Factor w/ 9 levels "COD","Con","ConLD",..: 9 9 9 9 9 9 9 9 9 9 ...
##  $ SaleCondition: Factor w/ 6 levels "Abnorml","AdjLand",..: 5 5 5 1 5 5 5 5 1 5 ...
##  $ SalePrice    : int  208500 181500 223500 140000 250000 143000 307000 200000 129900 118000 ...
dim(train)
## [1] 1460   81
So, this dataset contains 1460 rows and 80 variables. The first column is just a simple index therefore it is ignored. The 80th variable is SalePrice variable which is a target variable. Other 79 variables are predictor variables. These variables reflect various aspects of a house. To understand each variable better, we can cacluate summary statistics by following command chunk—

summary(train)
##        Id           MSSubClass       MSZoning     LotFrontage 
##  Min.   :   1.0   Min.   : 20.0   C (all):  10   Min.   : 21.00
##  1st Qu.: 365.8   1st Qu.: 20.0   FV     :  65   1st Qu.: 59.00
##  Median : 730.5   Median : 50.0   RH     :  16   Median : 69.00
##  Mean   : 730.5   Mean   : 56.9   RL     :1151   Mean   : 70.05
##  3rd Qu.:1095.2   3rd Qu.: 70.0   RM     : 218   3rd Qu.: 80.00
##  Max.   :1460.0   Max.   :190.0                  Max.   :313.00
##                                                  NA's   :259   
##     LotArea        Street      Alley      LotShape  LandContour
##  Min.   :  1300   Grvl:   6   Grvl:  50   IR1:484   Bnk:  63 
##  1st Qu.:  7554   Pave:1454   Pave:  41   IR2: 41   HLS:  50 
##  Median :  9478               NA's:1369   IR3: 10   Low:  36 
##  Mean   : 10517                           Reg:925   Lvl:1311 
##  3rd Qu.: 11602                                             
##  Max.   :215245                                             
##                                                             
##   Utilities      LotConfig    LandSlope   Neighborhood   Condition1
##  AllPub:1459   Corner : 263   Gtl:1382   NAmes  :225   Norm   :1260
##  NoSeWa:   1   CulDSac:  94   Mod:  65   CollgCr:150   Feedr  :  81
##                FR2    :  47   Sev:  13   OldTown:113   Artery :  48
##                FR3    :   4              Edwards:100   RRAn   :  26
##                Inside :1052              Somerst: 86   PosN   :  19
##                                          Gilbert: 79   RRAe   :  11
##                                          (Other):707   (Other):  15
##    Condition2     BldgType      HouseStyle   OverallQual 
##  Norm   :1445   1Fam  :1220   1Story :726   Min.   : 1.000
##  Feedr  :   6   2fmCon:  31   2Story :445   1st Qu.: 5.000
##  Artery :   2   Duplex:  52   1.5Fin :154   Median : 6.000
##  PosN   :   2   Twnhs :  43   SLvl   : 65   Mean   : 6.099
##  RRNn   :   2   TwnhsE: 114   SFoyer : 37   3rd Qu.: 7.000
##  PosA   :   1                 1.5Unf : 14   Max.   :10.000
##  (Other):   2                 (Other): 19                 
##   OverallCond      YearBuilt     YearRemodAdd    RoofStyle 
##  Min.   :1.000   Min.   :1872   Min.   :1950   Flat   :  13
##  1st Qu.:5.000   1st Qu.:1954   1st Qu.:1967   Gable  :1141
##  Median :5.000   Median :1973   Median :1994   Gambrel:  11
##  Mean   :5.575   Mean   :1971   Mean   :1985   Hip    : 286
##  3rd Qu.:6.000   3rd Qu.:2000   3rd Qu.:2004   Mansard:   7
##  Max.   :9.000   Max.   :2010   Max.   :2010   Shed   :   2
##                                                           
##     RoofMatl     Exterior1st   Exterior2nd    MasVnrType    MasVnrArea 
##  CompShg:1434   VinylSd:515   VinylSd:504   BrkCmn : 15   Min.   :   0.0
##  Tar&Grv:  11   HdBoard:222   MetalSd:214   BrkFace:445   1st Qu.:   0.0
##  WdShngl:   6   MetalSd:220   HdBoard:207   None   :864   Median :   0.0
##  WdShake:   5   Wd Sdng:206   Wd Sdng:197   Stone  :128   Mean   : 103.7
##  ClyTile:   1   Plywood:108   Plywood:142   NA's   :  8   3rd Qu.: 166.0
##  Membran:   1   CemntBd: 61   CmentBd: 60                 Max.   :1600.0
##  (Other):   2   (Other):128   (Other):136                 NA's   :8     
##  ExterQual ExterCond  Foundation  BsmtQual   BsmtCond    BsmtExposure
##  Ex: 52    Ex:   3   BrkTil:146   Ex  :121   Fa  :  45   Av  :221 
##  Fa: 14    Fa:  28   CBlock:634   Fa  : 35   Gd  :  65   Gd  :134 
##  Gd:488    Gd: 146   PConc :647   Gd  :618   Po  :   2   Mn  :114 
##  TA:906    Po:   1   Slab  : 24   TA  :649   TA  :1311   No  :953 
##            TA:1282   Stone :  6   NA's: 37   NA's:  37   NA's: 38 
##                      Wood  :  3                                   
##                                                                   
##  BsmtFinType1   BsmtFinSF1     BsmtFinType2   BsmtFinSF2   
##  ALQ :220     Min.   :   0.0   ALQ :  19    Min.   :   0.00
##  BLQ :148     1st Qu.:   0.0   BLQ :  33    1st Qu.:   0.00
##  GLQ :418     Median : 383.5   GLQ :  14    Median :   0.00
##  LwQ : 74     Mean   : 443.6   LwQ :  46    Mean   :  46.55
##  Rec :133     3rd Qu.: 712.2   Rec :  54    3rd Qu.:   0.00
##  Unf :430     Max.   :5644.0   Unf :1256    Max.   :1474.00
##  NA's: 37                      NA's:  38                   
##    BsmtUnfSF       TotalBsmtSF      Heating     HeatingQC CentralAir
##  Min.   :   0.0   Min.   :   0.0   Floor:   1   Ex:741    N:  95 
##  1st Qu.: 223.0   1st Qu.: 795.8   GasA :1428   Fa: 49    Y:1365 
##  Median : 477.5   Median : 991.5   GasW :  18   Gd:241           
##  Mean   : 567.2   Mean   :1057.4   Grav :   7   Po:  1           
##  3rd Qu.: 808.0   3rd Qu.:1298.2   OthW :   2   TA:428           
##  Max.   :2336.0   Max.   :6110.0   Wall :   4                     
##                                                                   
##  Electrical     X1stFlrSF      X2ndFlrSF     LowQualFinSF 
##  FuseA:  94   Min.   : 334   Min.   :   0   Min.   :  0.000
##  FuseF:  27   1st Qu.: 882   1st Qu.:   0   1st Qu.:  0.000
##  FuseP:   3   Median :1087   Median :   0   Median :  0.000
##  Mix  :   1   Mean   :1163   Mean   : 347   Mean   :  5.845
##  SBrkr:1334   3rd Qu.:1391   3rd Qu.: 728   3rd Qu.:  0.000
##  NA's :   1   Max.   :4692   Max.   :2065   Max.   :572.000
##                                                           
##    GrLivArea     BsmtFullBath     BsmtHalfBath        FullBath 
##  Min.   : 334   Min.   :0.0000   Min.   :0.00000   Min.   :0.000
##  1st Qu.:1130   1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:1.000
##  Median :1464   Median :0.0000   Median :0.00000   Median :2.000
##  Mean   :1515   Mean   :0.4253   Mean   :0.05753   Mean   :1.565
##  3rd Qu.:1777   3rd Qu.:1.0000   3rd Qu.:0.00000   3rd Qu.:2.000
##  Max.   :5642   Max.   :3.0000   Max.   :2.00000   Max.   :3.000
##                                                                 
##     HalfBath       BedroomAbvGr    KitchenAbvGr   KitchenQual
##  Min.   :0.0000   Min.   :0.000   Min.   :0.000   Ex:100   
##  1st Qu.:0.0000   1st Qu.:2.000   1st Qu.:1.000   Fa: 39   
##  Median :0.0000   Median :3.000   Median :1.000   Gd:586   
##  Mean   :0.3829   Mean   :2.866   Mean   :1.047   TA:735   
##  3rd Qu.:1.0000   3rd Qu.:3.000   3rd Qu.:1.000           
##  Max.   :2.0000   Max.   :8.000   Max.   :3.000           
##                                                           
##   TotRmsAbvGrd    Functional    Fireplaces    FireplaceQu   GarageType
##  Min.   : 2.000   Maj1:  14   Min.   :0.000   Ex  : 24    2Types :  6
##  1st Qu.: 5.000   Maj2:   5   1st Qu.:0.000   Fa  : 33    Attchd :870
##  Median : 6.000   Min1:  31   Median :1.000   Gd  :380    Basment: 19
##  Mean   : 6.518   Min2:  34   Mean   :0.613   Po  : 20    BuiltIn: 88
##  3rd Qu.: 7.000   Mod :  15   3rd Qu.:1.000   TA  :313    CarPort:  9
##  Max.   :14.000   Sev :   1   Max.   :3.000   NA's:690    Detchd :387
##                   Typ :1360                               NA's   : 81
##   GarageYrBlt   GarageFinish   GarageCars      GarageArea     GarageQual
##  Min.   :1900   Fin :352     Min.   :0.000   Min.   :   0.0   Ex  :   3
##  1st Qu.:1961   RFn :422     1st Qu.:1.000   1st Qu.: 334.5   Fa  :  48
##  Median :1980   Unf :605     Median :2.000   Median : 480.0   Gd  :  14
##  Mean   :1979   NA's: 81     Mean   :1.767   Mean   : 473.0   Po  :   3
##  3rd Qu.:2002                3rd Qu.:2.000   3rd Qu.: 576.0   TA  :1311
##  Max.   :2010                Max.   :4.000   Max.   :1418.0   NA's:  81
##  NA's   :81                                                           
##  GarageCond  PavedDrive   WoodDeckSF      OpenPorchSF     EnclosedPorch 
##  Ex  :   2   N:  90     Min.   :  0.00   Min.   :  0.00   Min.   :  0.00
##  Fa  :  35   P:  30     1st Qu.:  0.00   1st Qu.:  0.00   1st Qu.:  0.00
##  Gd  :   9   Y:1340     Median :  0.00   Median : 25.00   Median :  0.00
##  Po  :   7              Mean   : 94.24   Mean   : 46.66   Mean   : 21.95
##  TA  :1326              3rd Qu.:168.00   3rd Qu.: 68.00   3rd Qu.:  0.00
##  NA's:  81              Max.   :857.00   Max.   :547.00   Max.   :552.00
##                                                                         
##    X3SsnPorch      ScreenPorch        PoolArea        PoolQC 
##  Min.   :  0.00   Min.   :  0.00   Min.   :  0.000   Ex  :   2
##  1st Qu.:  0.00   1st Qu.:  0.00   1st Qu.:  0.000   Fa  :   2
##  Median :  0.00   Median :  0.00   Median :  0.000   Gd  :   3
##  Mean   :  3.41   Mean   : 15.06   Mean   :  2.759   NA's:1453
##  3rd Qu.:  0.00   3rd Qu.:  0.00   3rd Qu.:  0.000           
##  Max.   :508.00   Max.   :480.00   Max.   :738.000           
##                                                               
##    Fence      MiscFeature    MiscVal             MoSold   
##  GdPrv:  59   Gar2:   2   Min.   :    0.00   Min.   : 1.000
##  GdWo :  54   Othr:   2   1st Qu.:    0.00   1st Qu.: 5.000
##  MnPrv: 157   Shed:  49   Median :    0.00   Median : 6.000
##  MnWw :  11   TenC:   1   Mean   :   43.49   Mean   : 6.322
##  NA's :1179   NA's:1406   3rd Qu.:    0.00   3rd Qu.: 8.000
##                           Max.   :15500.00   Max.   :12.000
##                                                           
##      YrSold        SaleType    SaleCondition    SalePrice   
##  Min.   :2006   WD     :1267   Abnorml: 101   Min.   : 34900
##  1st Qu.:2007   New    : 122   AdjLand:   4   1st Qu.:129975
##  Median :2008   COD    :  43   Alloca :  12   Median :163000
##  Mean   :2008   ConLD  :   9   Family :  20   Mean   :180921
##  3rd Qu.:2009   ConLI  :   5   Normal :1198   3rd Qu.:214000
##  Max.   :2010   ConLw  :   5   Partial: 125   Max.   :755000
##                 (Other):   9
Let us organize the dataframe in more formatted and comprehensive way—

train.numeric=train[,c(81,4,5,27,35,37,38,39,44:47,63,67:72,76)]

train.numericfactor=train[,c(2, 18,19, 48:53,55,57,62,77)]

train.numericfactor=as.data.frame(lapply(train.numericfactor, as.factor))

PeriodBuilt=2018-train$YearBuilt
PeriodRemod=2018-train$YearRemodAdd
PeriodGarage=2018-train$GarageYrBlt
PeriodSold=2018-train$YrSold

train.period=cbind.data.frame(PeriodBuilt,PeriodGarage,PeriodRemod, PeriodSold)

train.factor=train[,c(3,6:17, 22:26,28:34,36,40:43,54,56,58,59,61, 64:66,73:75, 79,80)]

train.combined=cbind(train.numeric,train.numericfactor,train.period, train.factor)

dim(train.combined)
## [1] 1460   80
Let us understand the distribution of each numeric variable by observing their histograms—

library(e1071)

for (i in colnames(train.combined[,c(1:20)])){

  skew.measure1=round(skewness(train.combined[[i]]),digits = 2)

  skew.measure2=skewness(log(train.combined[[i]]))

  xname=names(train.combined[i])

  p1=hist(scale(train.combined[[i]], center=T, scale=T), n=100, xlab = xname, main=paste("Histogram of", xname, "with skewness", skew.measure1))

  p2=hist(log(train.combined[[i]]), n=100, xlab = xname, main=paste("Histogram of","log", xname, "with skewness", skew.measure2))


}


As we can observe the graphs, few variables with log transformation are normalized. We retain the log tranformation of these numeric variables—

train.combined$SalePrice=log(train.combined$SalePrice)

train.combined$LotArea=log(train.combined$LotArea)

train.combined$X1stFlrSF=log(train.combined$X1stFlrSF)

train.combined$GrLivArea=log(train.combined$GrLivArea)
Now, columns 1 to 38 are nummeric columns and 39 to 81 are factor columns. The correlation is bivariate analysis. The independent variables should be screened out is there is no significant correlation with dependent variable “SalePrice”. Let us the correlation of SalePrice with each of the numeric columns by plotting scatter plots—

for (i in 2:length(train.combined[,c(2:20)])) {
  a <- cor.test(train.combined$SalePrice, train.combined[,i])
  print(paste(colnames(train.combined)[i], " est:", a$estimate, " p=value:", a$p.value))

}
## [1] "LotFrontage  est: 0.355878470385983  p=value: 3.55822388129835e-37"
## [1] "LotArea  est: 0.39991774112559  p=value: 3.47429417639357e-57"
## [1] "MasVnrArea  est: 0.43080852419511  p=value: 1.10796556420336e-66"
## [1] "BsmtFinSF1  est: 0.372023073567088  p=value: 3.84595104338152e-49"
## [1] "BsmtFinSF2  est: 0.00483241053204665  p=value: 0.853629797314542"
## [1] "BsmtUnfSF  est: 0.22198505352521  p=value: 9.31852717731222e-18"
## [1] "TotalBsmtSF  est: 0.612133975369787  p=value: 7.53455105447948e-151"
## [1] "X1stFlrSF  est: 0.608946654903631  p=value: 6.97929348105704e-149"
## [1] "X2ndFlrSF  est: 0.319299984347205  p=value: 5.86690165149455e-36"
## [1] "LowQualFinSF  est: -0.037962803137427  p=value: 0.147104011505013"
## [1] "GrLivArea  est: 0.730254851198229  p=value: 1.5984652370953e-243"
## [1] "GarageArea  est: 0.650887555902007  p=value: 1.10625537128625e-176"
## [1] "WoodDeckSF  est: 0.334135073957513  p=value: 2.05532422562011e-39"
## [1] "OpenPorchSF  est: 0.321052972019767  p=value: 2.34462590169273e-36"
## [1] "EnclosedPorch  est: -0.149050281427676  p=value: 1.05080829167697e-08"
## [1] "X3SsnPorch  est: 0.0549002264676204  p=value: 0.0359475707337999"
## [1] "ScreenPorch  est: 0.121207604896864  p=value: 3.40897065981049e-06"
## [1] "PoolArea  est: 0.0697978060096918  p=value: 0.00763178070309401"
Now, we can visualize the correlation between SalePrice and other numeric variables through scatter plots—

for (i in 2:length(train.combined[,c(2:20)])) {

  p <- plot(scale(train.combined[,i], center = T, scale = T), train.combined$SalePrice,
            xlab =colnames(train.combined)[i], ylab = "SalePrice", main=paste("Cor between Sale Price and", colnames(train.combined)[i]))
 
  p
}


We need to impute missing values in LotFrontage and MasVnrArea variables through regression—

mode.LotFrontage=lm(train.combined$LotFrontage~.,data=train.combined[,c(1:20)])

summary(mode.LotFrontage)
##
## Call:
## lm(formula = train.combined$LotFrontage ~ ., data = train.combined[,
##     c(1:20)])
##
## Residuals:
##     Min      1Q  Median      3Q     Max
## -73.485  -7.729  -0.465   6.828 210.494
##
## Coefficients: (1 not defined because of singularities)
##                 Estimate Std. Error t value Pr(>|t|) 
## (Intercept)   -1.876e+02  2.813e+01  -6.668 3.97e-11 ***
## SalePrice     -5.547e+00  2.366e+00  -2.344 0.019238 *
## LotArea        2.750e+01  1.229e+00  22.379  < 2e-16 ***
## MasVnrArea     7.095e-03  3.222e-03   2.202 0.027850 *
## BsmtFinSF1     5.110e-03  2.169e-03   2.356 0.018644 *
## BsmtFinSF2     5.484e-03  3.782e-03   1.450 0.147305 
## BsmtUnfSF      4.420e-03  2.125e-03   2.080 0.037743 *
## TotalBsmtSF           NA         NA      NA       NA 
## X1stFlrSF      9.269e+00  5.598e+00   1.656 0.098050 .
## X2ndFlrSF      4.036e-03  4.349e-03   0.928 0.353623 
## LowQualFinSF   1.170e-02  1.048e-02   1.116 0.264562 
## GrLivArea     -1.427e-01  6.700e+00  -0.021 0.983008 
## GarageArea     9.983e-03  3.144e-03   3.175 0.001537 **
## WoodDeckSF    -1.571e-02  4.615e-03  -3.405 0.000685 ***
## OpenPorchSF    3.904e-03  8.431e-03   0.463 0.643431 
## EnclosedPorch  4.457e-03  8.592e-03   0.519 0.604010 
## X3SsnPorch     2.988e-02  1.769e-02   1.689 0.091462 .
## ScreenPorch   -2.608e-02  9.432e-03  -2.765 0.005776 **
## PoolArea       7.337e-02  1.389e-02   5.281 1.53e-07 ***
## MiscVal       -6.408e-04  2.729e-03  -0.235 0.814387 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.49 on 1176 degrees of freedom
##   (265 observations deleted due to missingness)
## Multiple R-squared:  0.4893, Adjusted R-squared:  0.4815
## F-statistic:  62.6 on 18 and 1176 DF,  p-value: < 2.2e-16
pred=predict(mode.LotFrontage, newdata=train.combined[,c(1:20)])
## Warning in predict.lm(mode.LotFrontage, newdata = train.combined[,
## c(1:20)]): prediction from a rank-deficient fit may be misleading
train.combined$LotFrontage=ifelse(is.na(train.combined$LotFrontage),pred,
                              train.combined$LotFrontage)

train.combined$LotFrontage=round(train.combined$LotFrontage, digits = 0)
We can see there are some numerical variables which have majority value zero. These variables can be used to create binanry variables derived from their respective numeric variable—

train.combined$MasVnrArea_yes=ifelse(train.combined$MasVnrArea>0,1,0)

train.combined$BsmtFinSF1_yes=ifelse(train.combined$BsmtFinSF1>0,1,0)

train.combined$BsmtUnfSF_yes=ifelse(train.combined$BsmtUnfSF>0,1,0)

train.combined$X2ndFlrSF_yes=ifelse(train.combined$X2ndFlrSF>0,1,0)

train.combined$GarageArea_yes=ifelse(train.combined$GarageArea>0,1,0)

train.combined$WoodDeckSF_yes=ifelse(train.combined$WoodDeckSF>0,1,0)

train.combined$OpenPorchSF_yes=ifelse(train.combined$OpenPorchSF>0,1,0)

train.combined$EnclosedPorch_yes=ifelse(train.combined$EnclosedPorch>0,1,0)

train.combined$X3SsnPorch_yes=ifelse(train.combined$X3SsnPorch>0,1,0)

train.combined$ScreenPorch_yes=ifelse(train.combined$ScreenPorch>0,1,0)

train.combined$PoolArea_yes = ifelse(train.combined$PoolArea>0,1,0)

train.combined$MiscVal_yes=ifelse(train.combined$MiscVal>0,1,0)
Anova is a technique to compare multiple group means. As in our datframe, we have lot of categorical and dummy variables, we can check group means various categorical variables with respect to SalePrice are significantly different or nor—

for(i in colnames(train.combined[,c(21:33,38:92)])){

boxplots=boxplot(train.combined$SalePrice~train.combined[[i]],xlab=colnames(train.combined[i]), main=paste("Boxplot SalePrice Vs", colnames(train.combined[i])))

boxplots

}


Before model building we need to work out missing values in data frame. To impute missing values, we can use missForest package—

library(missForest)
## Loading required package: randomForest
## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.
## Loading required package: foreach
## Loading required package: itertools
## Loading required package: iterators
set.seed(1000)

train.withoutMV=randomForest::rfImpute(x=train.combined[,2:92],
                                       y=train.combined[,1],
                                       ntree=500,
                                       iter=5)
##      |      Out-of-bag   |
## Tree |      MSE  %Var(y) |
##  500 |  0.01983    12.44 |
##      |      Out-of-bag   |
## Tree |      MSE  %Var(y) |
##  500 |  0.01979    12.41 |
##      |      Out-of-bag   |
## Tree |      MSE  %Var(y) |
##  500 |   0.0195    12.23 |
##      |      Out-of-bag   |
## Tree |      MSE  %Var(y) |
##  500 |  0.01956    12.27 |
##      |      Out-of-bag   |
## Tree |      MSE  %Var(y) |
##  500 |  0.01985    12.45 |
train.withoutMV
train.combined[, 1]
<dbl>
LotFrontage
<dbl>
LotArea
<dbl>
MasVnrArea
<dbl>
BsmtFinSF1
<int>
BsmtFinSF2
<int>
BsmtUnfSF
<int>
TotalBsmtSF
<int>
12.24769 65.00000 9.041922 196.00000 706 0 150 856
12.10901 80.00000 9.169518 0.00000 978 0 284 1262
12.31717 68.00000 9.328123 162.00000 486 0 434 920
11.84940 60.00000 9.164296 0.00000 216 0 540 756
12.42922 84.00000 9.565214 350.00000 655 0 490 1145
11.87060 85.00000 9.554993 0.00000 732 0 64 796
12.63460 75.00000 9.218705 186.00000 1369 0 317 1686
12.20607 77.00000 9.247829 240.00000 859 32 216 1107
11.77452 51.00000 8.719317 0.00000 0 0 952 952
11.67844 50.00000 8.911934 0.00000 851 0 140 991
Next123456...146Previous
1-10 of 1,460 rows | 1-8 of 92 columns
Let us check how many missing values after imputation—

sum(is.na(train.withoutMV))
## [1] 0
The output shows that no missing values are left in the dataset train.withoutMV.Now let us do some model building with randomForest technique—

train.withoutMV$SalePrice=train.withoutMV$`train.combined[, 1]`

train.withoutMV$`train.combined[, 1]`=NULL
Now, it is time to build model. First, we will try a simple random forest model—

set.seed(1001)

model.rf=randomForest(SalePrice~.,
                      data=train.withoutMV,
                      ntree=500,
                      mtry=5,
                      importance=TRUE,
                      proximity=TRUE)

pred=predict(model.rf,
             newdata=train.withoutMV,
             type="response")

train.withoutMV$pred.rf=exp(pred)

train.withoutMV$SalePrice=exp(train.withoutMV$SalePrice)

library(Metrics)

rmse.rf=rmse(train.withoutMV$SalePrice, train.withoutMV$pred.rf)

rmse.rf
## [1] 16054.31
mape.rf=Metrics::mape(train.withoutMV$SalePrice, train.withoutMV$pred.rf)

mape.rf
## [1] 0.04494797
bias.rf=bias(train.withoutMV$SalePrice, train.withoutMV$pred.rf)

bias.rf
## [1] 2760.606
percent_bias.rf=percent_bias(train.withoutMV$SalePrice, train.withoutMV$pred.rf)

percent_bias.rf
## [1] -0.0028082
plot(model.rf)


varImpPlot(model.rf)

CASE STUDY 3: Loan Prediction with Random Forest

MACHINE LEARNING FOR LOAN PREDICTION

The loan approval depends on many parameters and based on those parameters, loan application is either rejected or accepted. Loan dataset is used in this case study. This dataset is a hackathon  dataset available on website : https//:www.analyticsvidhya.com. Following is the link to access the dataset---

https://datahack.analyticsvidhya.com/contest/practice-problem-loan-prediction-iii/

The dataset is divided into two sets-- train and test. The train dataset is used for model building and test dataset is used for submitting the predictions to the hackthon for evaluation. Purpose of model building is to automate the loan approval process through machine learning algorithm.

PREPARATION OF WORKING DIRECTORY, LIBRARIES AND DATASET

So, lets understand the model building for this case study.  It starts with loading the proper working directory---

# Loading the working directory

setwd(E:/r/R/files)
getwd()

To perform the analysis in R, required libraries or packages should be called in working directory---

# Loading the libraries

library(dplyr)
library(e1071)
library(Hmisc)
library(ggplot2)
library(randomForest)
library(caret)
library(pROC)

There are two datasets-- train and test in csv format. The data is loaded in working directory---

# Loading the datasets

train=read.csv("train.csv", header = T, stringsAsFactors = F)
test=read.csv("test.csv", header=T, stringsAsFactors = F)

In train dataset, there are 614 applicants and 13 variables and in test datset, there are 367 applicants an 12 variables.  Before combining, a "DataType" variable is added to both datasets with categorical values "Train" and "Test" respectively.  In test dataset, Loan_Status variable is absent therefore it is added in test dataset.

For model building, it is better to have more number of observations, hence we combine both datasets. To combine two datasets successfully , the datasets should have equal number of variables.  and than, both train and test datasets are combined  ---

train$DataType="Train"
test$DataType="Test"
test$Loan_Status=NA

df=rbind(train,test)

Now, the combined dataset "df" has following structure---

> str(df)
'data.frame': 981 obs. of  14 variables:
 $ Loan_ID                     : chr  "LP001002" "LP001003" "LP001005" "LP001006" ...
 $ Gender                        : chr  "Male" "Male" "Male" "Male" ...
 $ Married                       : chr  "No" "Yes" "Yes" "Yes" ...
 $ Dependents                 : chr  "0" "1" "0" "0" ...
 $ Education                    : chr  "Graduate" "Graduate" "Graduate" "Not Graduate" ...
 $ Self_Employed           : chr  "No" "No" "Yes" "No" ...
 $ ApplicantIncome        : int  5849 4583 3000 2583 6000 5417 2333 3036 4006 12841 ...
 $ CoapplicantIncome     : num  0 1508 0 2358 0 ...
 $ LoanAmount               : int  NA 128 66 120 141 267 95 158 168 349 ...
 $ Loan_Amount_Term   : int  360 360 360 360 360 360 360 360 360 360 ...
 $ Credit_History             : int  1 1 1 1 1 1 1 0 1 1 ...
 $ Property_Area              : chr  "Urban" "Rural" "Urban" "Urban" ...
 $ Loan_Status                  : chr  "Y" "N" "Y" "Y" ...
 $ DataType                       : chr  "Train" "Train" "Train" "Train" ...

MISSING VALUES TREATMENT AND VARIABLES MANIPULATION

The dateset contains few missing values which has to be treated carefully and data formats needs to be changed with appropriate transformations. The summary of missing values is as follows---

>MV_summary=sapply(df, function(x) sum(is.na(x)))
 
>MV_summary
          Loan_ID          Gender       Married      Dependents    Education     Self_Employed
                0                 0                 0                 0                    0                   0
  ApplicantIncome CoapplicantIncome        LoanAmount  Loan_Amount_Term    Credit_History   
                0               0                                   27                    20                                  79               

Property_Area     Loan_Status          DataType
               0                367                           0
                             
We can see that Loan Amount, Loan Amount Term and Credit History have 27, 20 and 79 missing values respectively. Loan Status also has 367 missing values which are deliberately kept missing for submission in Hackathon.

There are many character type variables in the dataset. These variables should be converted into numeric values so that statistical computation can be done. For example, The variable Gender has two options-- male and female. In numeric form male can be converted into "1" and female into "0". Similar transformations are carried out for Married, Dependents, Education, Credit History, Loan Status and Property Area.

In Dependent variable, there are 25 missing values but they were not visible into missing value summary (MV_summary) which may be due to some mistakes while doing data feeding---

> table(df$Dependents)

       0     1     2      3+
 25 545 160 160  91

As we can see that 0 dependent category is the mode of Dependents variable hence we can replace missing value with 0 dependent category by following command---

>df$Dependents=ifelse(df$Dependents=="", 0, df$Dependents)

We can also observe that there is 3+ Dependent category in Dependents variable. The value 3+ can not be used for statistical computation so it is converted from 3+ to 3 and categories are converted to numeric type by following commands---

>df$Dependents=substr(df$Dependents,1,1)

>df$Dependents=substr(df$Dependents,1,1)

To indicate, if applicant has dependent or not, a new variable "depend_yes" is introduced. If there is dependent(s) than depend_yes variable will take value 1 or else zero---

>df$depend_yes=ifelse(df$Dependents>0,1,0)

Variable "Education" is a dichotomous variable having value "Graduate" and "Non Graduate". By following commands, Education variable is converted into numeric form and "Education" variable is removed from dataframe and new variable "Graduate" is introduced---

>df=df%>%
  mutate(graduate=as.numeric(Education=="Graduate"))

>df$Education=NULL

Similarly, "Self_Employed" variable is converted from character type to numeric type by following commands---

>df=df%>%
  mutate(SE_yes=as.numeric(Self_Employed=="Yes"))

>df$Self_Employed=NULL

In ApplicantIncome variable, there are two applicants with zero income. We fill up these cells with median of ApplicantIncome which is more appropriate than mean of ApplicantIncome because mean is sensitive to extreme values while median is not---

>df$ApplicantIncome[df$ApplicantIncome==0]=median(df$ApplicantIncome,na.rm = TRUE)

A new variable TotalIncome is introduced by adding ApplicantIncome and CoapplicantIncome---

>df$TotalIncome=df$ApplicantIncome+df$CoapplicantIncome

To indicate whether there is CoapplicantIncome or not, a new variable CoappIncome_yes is introduced---

>df$CoappIncome_yes=ifelse(df$CoapplicantIncome>0,1,0)

In original dattaframe, LoanAmount variable in given in 1000s which should be converted into full figures which is done by following command---

>df$LoanAmount=1000*df$LoanAmount

There are missing values in LoanAmount variable. It can be assumed that LoanAmount should vary as per Loan_Amount_Term therefore we replace NAs in LoanAmount by groupwise median LoanAmount for each Loan_Amount_Term by following commands---

>df$LoanAmount[is.na(df$LoanAmount)] <- ave(df$LoanAmount,
                                 df$Loan_Amount_Term,
                                 FUN=function(x)median(x,
                                                     na.rm = T))[is.na(df$LoanAmount)]  

There are 20 missing values in  Loan_Amount_Term. The Loan_Amount_Term 360 is highest frequency term hence all missing values are replaced with 360 by following command---

>df$Loan_Amount_Term[is.na(df$Loan_Amount_Term)]=360

In Credit_History variable, there are 79 missing variables. It has to be filled carefully. Let us observe distribution of Credit_History with Loan_Status---

> table(df$Loan_Status, df$Credit_History)

   

        0    1

  N   82  97

  Y   7    378

We can see that Applicants who have Credit_History=1 , their Loan_Status is "Y". Hence we can fill Credit_History of those Applicant as 1 where Loan_Status is "Y"---

>df$Credit_History[is.na(df$Credit_History)]=ifelse(df$Loan_Status=="Y",1,df$Credit_History)

The table below shows the distribution of credit history as per Education and Loan_Status---

> table(df$Loan_Status, df$Credit_History, df$Education)

, ,  = Graduate

   

       0    1

  N  57  73

  Y   6    307

, ,  = Not Graduate

   

        0    1

  N  25  24

  Y   1   71

 We can see that Credit_History is 1 for most of the cases where Loan_status is "Y" and Education is "Graduate". So, we can replace those NAs of Credit_History variable where Loan_Status is "Y" and Education is "Graduate". by following command---

>df$Credit_History[is.na(df$Credit_History)]=ifelse(df$Loan_Status=="Y" &  

  df$Education=="Graduate",1,df$Credit_History)

Similarly, we observe the distribution of Credit_History with Loan_Status, Education and SE_yes---

> table(df$Loan_Status, df$Credit_History, df$Education, df$SE_yes)

, ,  = Graduate,  = 0

   

        0   1

  N  50  71

  Y   7   287

, ,  = Not Graduate,  = 0

   

       0    1

  N  23  22

  Y   2   70

, ,  = Graduate,  = 1

   

       0    1

  N  10   9

  Y   2    44

, ,  = Not Graduate,  = 1

   

        0   1

  N   2   5

  Y   0  10

We can observe that Applicant who is Graduate, either self employed or not, the Credit_History is most likely 1 if Loan_Status is also "Y". We replace NAs with following commands---

> df$Credit_History[is.na(df$Credit_History)]=ifelse(df$Loan_Status=="Y" & df$Education=="Graduate" & df$SE_yes==1,1,df$Credit_History)

> df$Credit_History[is.na(df$Credit_History)]=ifelse(df$Loan_Status=="Y" & df$Education=="Graduate" & df$SE_yes==0,1,df$Credit_History)

> which(is.na(df$Credit_History))

integer(0)

> table(df$Credit_History)

  0    1 

163  818 

Loan_Status variable should also be changed from character variable to numeric variable which is done by following command---

>df$Loan_Status=ifelse(df$Loan_Status=="Y",1,0)

For  approval of LoanAmount, the LoanAmount should not exceed a certain multiple of TotalIncome. According to literature survey, it is found that Loan to Income ratio should not exceed 40 times of Income in case of Self Employed  applicants and 60 times in case of non Self Employed applicants. Hereby, we have created a variable "LoantoIncomeRatio", LoanIncomeRatio_for_self and LoanIncomeRatio_for_nonself= ifelse(df$LoantoIncomeRatio<=60 & df$SE_yes==0,1,0)

---

>df$LoantoIncomeRatio=df$LoanAmount/df$TotalIncome

>df$LoanIncomeRatio_for_self= ifelse(df$LoantoIncomeRatio<=40 & df$SE_yes==1,1,0)

>df$LoanIncomeRatio_for_nonself= ifelse(df$LoantoIncomeRatio<=60 & df$SE_yes==0,1,0)

There are three types of Property_Area in the dataset--- Urban, Semiurban and Rural. This variable is mutated into two new variables-- Area_Urban and Area_Semiurban. 

Finally, the Loan_Status variable is mutated into numeric variable by following command---

>df$Loan_Status=ifelse(df$Loan_Status=="Y",1,0)

Now, let us check whether there is any missing value left by following commands---

> MV_summary=sapply(df, function(x) sum(is.na(x)))

> MV_summary

                    Loan_ID                   Dependents            ApplicantIncome 

                          0                          0                            0 

          CoapplicantIncome            LoanAmount          Loan_Amount_Term 

                          0                          0                            0 

             Credit_History                 Loan_Status          DataType 

                          0                         367                         0 

                       male                       married                  depend_yes 

                          0                          0                            0 

                     SE_yes                    TotalIncome           LoantoIncomeRatio 

                          0                          0                            0 

   LoanIncomeRatio_for_self       LoanIncomeRatio_for_nonself     area_urban 

                          0                          0                                                    0 

             area_semiurban               graduate 

                          0                          0 

Except for Loan_Status, all variables have no missing values.

The ApplicantIncome and CoapplicantIncome variables are redundant now so, we can remove them from dataset---

df=df[,-c(3,4)]

The TotalIncome and LoanAmount variables are continuous variables. Let us observe their distribution through histogram attached below---

> hist(df$LoanAmount, n=100)

> hist(df$TotalIncome, n=100)

We can see that TotalIncome and LoanAmount variables have skewed distribution towards higher side. For modeling purpose, the variables should be normally distributed like bell shaped curve. To remove skewness, we apply log transformation and observe the histograms again---

> df$LoanAmount=log(df$LoanAmount)

> df$TotalIncome=log(df$TotalIncome)

> hist(df$LoanAmount, n=100)

> hist(df$TotalIncome, n=100)

                               



Now, the dataset "df" is divided into two datasets--- "df_train" and "df_test". The first 614 rows are taken in "df_train" and last 367 rows are taken in "df_test".

>df_train=df[c(1:614),]
>df_test=df[615:981,]

MODEL BUILDING WITH TRAINING DATASET

For model building, we can start with logistic regression as the response variable i.e. Loan_Status is binary. The formula for model is as follows---

> formula=as.formula(Loan_Status~ LoanAmount
                                                            +Loan_Amount_Term
                                                            +Credit_History           
                                                            +male                     
                                                            +married                 
                                                            +Dependents               
                                                            +depend_yes               
                                                            +graduate                 
                                                            +SE_yes                   
                                                            +LoantoIncomeRatio
                                                            +area_urban               
                                                            +area_semiurban           
                                                            +TotalIncome
                                                            +LoanIncomeRatio_for_self
                                                            +LoanIncomeRatio_for_nonself)

Logistic regression is special case of linear regression where dependent variable is categorical variable. Following commands are used to generate logistic model---

>model_glm=glm(formula,data=df_train,family="binomial" )
>summary(model_glm)

Model fit index for logistic regression model is AIC which should be as minimum as possible. for model_glm, AIC is 597.9.  There are only three variables which are significant in the regression equation i.e. married, area_semiurban and Credit_History. To fine tune model, we will iteratively remove the variable with highest p-value one by one and observe the summary of the model until we are left with equation having only significant variables. To keep distinction between initial formula and final formula, there are named as "formula" and "formula.edited".

After number of iterations, we have got following results---

>formula.edited=as.formula(Loan_Status~
                  Credit_History           
                  +married                 
                  +area_semiurban           
                  +LoanIncomeRatio_for_self
                  +LoanIncomeRatio_for_nonself)


>model_glm=glm(formula.edited,data=df_train,family="binomial" )

>summary(model_glm)

Now, we can see that there are five predictor variables found significant i.e. LoanIncomeRatio_for_nonself, LoanIncomeRatio_for_self, area_urban, married and Credit_History.

PREDICTION OF LOAN_STATUS

After modeling, we need to predict the Loan_Status variable for all cases of df_train dataset. The predict function estimates the probability for each loan applicant. Following command adds "predict_glm" column to df_train dataset which contains probabilities---

>df_train$Predict_glm=predict(model_glm, type="response", newdata = df_train)

> str(df_train$Predict_glm)

    num [1:614] 0.692 0.79 0.799 0.79 0.692 ...

Random forest technique is advance level modeling technique in which number of models are build based on certain parameters and best model is chosen algorithmically. Following command applies RandomForest technique to fine tune the model and predict the probabilities for each applicant---

> model_rf=randomForest(formula.edited, data=df_train, importance=TRUE, proximity=TRUE,eig=TRUE, ntree=10000,type="classification")

>df_train$predict_rf=predict(model_rf, type="response", newdata=df_train)

>str(df_train$predict_rf)

 num [1:614] 0.694 0.733 0.728 0.733 0.694 ...

After modeling with Random Forest, the cutoff for prediction probabilities should ascertained to classify the Applicant into either Loan_Status=Y or Loan_Status=N. Iteratively, we will change the cutoff and see how predicted classification improves. The optimal cutoff found is 0.68 where maximum accuracy, sensitivity and specificity is achieved---

> cutoff=0.68

> df_train$cf=ifelse(df_train$predict_rf<cutoff,0,1)

> df_confusion_matrix=cbind(Survived=df_train$Loan_Status,Predicted=df_train$cf)

> df_confusion_matrix=as.data.frame(df_confusion_matrix)

> df_confusion_matrix$Survived=as.factor(df_confusion_matrix$Survived)

> df_confusion_matrix$Predicted=as.factor(df_confusion_matrix$Predicted)

> confusionMatrix(df_confusion_matrix$Survived, df_confusion_matrix$Predicted)

Confusion Matrix and Statistics

                       Reference

Prediction          0     1 

         0                93  99

         1                20 402

                                          

               Accuracy              : 0.8062          

                 95% CI               : (0.7727, 0.8367)

    No Information Rate       : 0.816           

    P-Value [Acc > NIR]      : 0.7526          

                                          

                                 Kappa : 0.4922          

     Mcnemar's Test P-Value : 8.662e-13       

                                          

            Sensitivity                : 0.8230          

            Specificity                : 0.8024          

         Pos Pred Value            : 0.4844          

         Neg Pred Value           : 0.9526          

             Prevalence               : 0.1840          

         Detection Rate             : 0.1515          

   Detection Prevalence         : 0.3127          

      Balanced Accuracy         : 0.8127          

                     'Positive' Class : 0    

CONCLUSION

The loan applicant with good credit history, married, having property in semiurban area, having Loan to Income Ratio upto 40 times for self employed and upto 60 times for non self employed have higher chances of getting loan approved. The loan applicant's loan status does not depend on gender, education, number of dependents, loan amount and total income. 

CASE STUDY 2: Who Did Survive in Titanic Disaster?

Titanic is the classic man made disaster. The brilliant craft of man kind, the great Titanic ship, sunk in the sea after striking an iceberg in its first voyage. The disaster had more grievous impact on passengers due to lack of life boats in the ship. 

There were many classes of passengers, males, females, variety of age groups from children to old age citizens in the ship. The dataset we are going to analyze in this case study contains the passengers details with information of whether passenger survived the disaster or not. Let us look at the data set--

The dataset contains details of 891 passengers. Each passenger's detail is explained through 12 columns or variables. The variables and their definitions are explained in following table---

The column "survival" tells that whether passenger survived the disaster or not (1= Survived, 0=Did not survive). Out of 891 passengers, 342 passengers survived and 549 passengers died. 

Questions which should be answered to find out the survival pattern among the passengers, corresponding hypotheses, relevant tests, results and interpretation to be tested are as following---

a) Whether "Ticket class" has any influence on chances of survival?

Ho : The survival of ith passenger is independent of Ticket class of ith passenger  

H1  : The survival of ith passenger is dependent of Ticket class of ith passenger

Statistical tool: Chi Squared Test of independence

Result and Interpretation:

> table(titanic$Survived, titanic$Pclass)
 
      1      2    3
  0  80    97  372
  1  136  87  119

We can observe from the the table that passengers of class 1 has got more survivors proportion while passengers of class 3 has more deaths.

> chisq.test(table(titanic$Survived,titanic$Pclass))

Pearson's Chi-squared test

data:  table(titanic$Survived, titanic$Pclass)
X-squared = 102.89, df = 2, p-value < 2.2e-16

The Chi-squared test result also shows that passengers' class has influence on survival chances of passengers. The class 1 passengers are safer and class 3 passengers are more vulnerable.

b) Whether gender of passenger has any influence on chances of survival?

Ho : The survival of ith passenger is independent of gender of ith passenger  

H1  : The survival of ith passenger is dependent of gender of ith passenger

Statistical tool: Chi Squared Test of independence

Result and Interpretations:

> table(titanic$Survived, titanic$Sex)
 
         female  male
  0     81         468
  1     233       109

We can observe that if passenger is female, chances of survival are more than the fellow male passengers.

> chisq.test(table(titanic$Survived, titanic$Sex))

Pearson's Chi-squared test with Yates' continuity correction

data:  table(titanic$Survived, titanic$Sex)
X-squared = 260.72, df = 1, p-value < 2.2e-16

The Chi-squared test result shows that the chances of survival depends on sex of the passenger.

c) Whether age of passenger has any influence on chances of survival?

Ho : There is no significant difference among the group means of survivors and non survivors.  

H1  : There is significant difference among the group means of survivors and non survivors.

Statistical tool: Two sample t-test

Before doing two sample t test, the variance of age of survived passengers and dead passengers has to be compared to see whether variance of both dataset is equal or not. The following is the result of F test---

Result and Interpretations---

> var.test(titanic_survived$Age, titanic_dead$Age)

F test to compare two variances

data:  titanic_survived$Age and titanic_dead$Age

F = 1.1129, num df = 289, denom df = 423, p-value = 0.317

alternative hypothesis: true ratio of variances is not equal to 1

95 percent confidence interval:

 0.9023453 1.3785538

sample estimates:

ratio of variances 

          1.112932  

The p-value is greater than 0.05 hence age variable of both datasets have homogeneous variance.

> t.test(titanic_survived$Age, titanic_dead$Age, var.equal = TRUE, paired = FALSE)

Two Sample t-test

data:  titanic_survived$Age and titanic_dead$Age
t = -2.0667, df = 712, p-value = 0.03912
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -4.450798 -0.114181
sample estimates:
mean of x mean of y
 28.34369  30.62618

We can observe that the difference in age of titanic survivor and titanic non-survivor is significant at 0.05 significance level. 

d) The influence of having siblings or spouse and number of siblings on survival chances of passengers need to be studied thoroughly. It will be done in two parts---

i. Whether having a sibling or spouse have influence on chances of survival?

Ho : The survival of ith passenger is independent of having siblings or spouse  

H1  : The survival of ith passenger is dependent of having siblings or spouse

Statistical tool: Chi Squared Test of independence

Result and Interpretations---

> table(titanic$Survived, titanic$SibSp_yes)
 
      0      1
  0  398  151
  1  210  132

> chisq.test(table(titanic$Survived, titanic$SibSp_yes))

Pearson's Chi-squared test with Yates' continuity correction

data:  table(titanic$Survived, titanic$SibSp_yes)
X-squared = 11.456, df = 1, p-value = 0.0007128

The p-vslue is much less than 0.05 hence null hypothesis is rejected and it can be interpreted that there is dependence of passengers' survival chances on having siblings or spouse.

ii. Whether number of siblings or spouse have influence on chances of survival?

To answer this question we can observe following tables---

> table(titanic$SibSp)

  0   1      2    3    4     5   8
608 209  28  16  18   5   7

> table(titanic$Survived, titanic$SibSp)
 
      0      1       2   3    4    5   8
  0  398  97    15  12  15  5   7
  1  210  112  13   4    3   0   0

We can understand very well that having no sibling or spouse clearly reduces chances of passenger survival and having only one sibling or spouse has greater chances of survival than having no sibling or more than one sibling.

e) The influence of having parents or children and number of children on survival chances of passengers need to be studied thoroughly. It will be done in two parts---

i. Whether having parents or children have influence on chances of survival?

Ho : The survival of ith passenger is independent of having parents or children  

H1  : The survival of ith passenger is dependent of having parents or children

Statistical tool: Chi Squared Test of independence

Result and Interpretations---

> table(titanic$Survived, titanic$parch_yes)
 
      0    1
  0 445 104
  1 233 109

> chisq.test(table(titanic$Survived, titanic$parch_yes))

Pearson's Chi-squared test with Yates' continuity correction

data:  table(titanic$Survived, titanic$parch_yes)
X-squared = 18.656, df = 1, p-value = 1.565e-05

It is evident from Chi-squared test that having parents or children influence chances of passenger's survival.

ii. Whether number of children have influence on chances of survival?

The tables below shows how number of children influence the chances of passenger survival---

> table(titanic$Parch)

  0   1       2    3   4   5   6
678 118  80   5   4   5   1

> table(titanic$Survived, titanic$Parch)
 
      0     1     2    3   4   5   6
  0 445  53  40   2   4   4   1
  1 233  65  40   3   0   1   0

We can observe that having no children or parents increases chances of non survival and having only one sibling or parent increases chances of survival but having children more than one decreases chances of survival.

g Whether fare of passenger has any influence on chances of survival?

Ho : There is no significant difference among the group means of survivors and non survivors.  

H1  : There is significant difference among the group means of survivors and non survivors.

Statistical tool: Two sample t-test

Result and Interpretations---

> var.test(titanic_survived$Fare, titanic_dead$Fare)

F test to compare two variances

data:  titanic_survived$Fare and titanic_dead$Fare
F = 4.5017, num df = 341, denom df = 548, p-value < 2.2e-16
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 3.725366 5.463382
sample estimates:
ratio of variances
          4.501697

The F test result shows that variance in fares of survived and non survived passengers is not homogeneous. Hence, two sample t test is applied with due consideration of non homogeneous variances.

> t.test(titanic_survived$Fare, titanic_dead$Fare, var.equal = FALSE, paired=FALSE)

Welch Two Sample t-test

data:  titanic_survived$Fare and titanic_dead$Fare
t = 6.8391, df = 436.7, p-value = 2.699e-11
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 18.72592 33.82912
sample estimates:
mean of x mean of y
 48.39541  22.11789 

The p-value is less than 0.05 hence it can be concluded that the mean fares of survived passengers and non survived passengers have significant difference. We can see that mean fare for  survived passengers is $48.40 and mean fare for non survived passengers is $22.12. There is clear evidence that passengers paying lower fares had discriminatory low opportunity to survive.

LOGISTIC REGRESSION MODELING

The passenger on Titanic either survived or died after striking the iceberg disaster. In terms of modeling, the outcome of the disaster is binary. Other variables in dataset are independent variables which are influencing the probability of passenger's survival.

The variables which were considered as independent variables are Pclass, Age, SibSp, Parch, SibSp_yes, parch_yes and female_yes. Pclass variable has three categories. Age variable is a continuous variable. SibSp variable is number of siblings or spouse. Parch variable is number of children or parents. SIbSp_yes variable is dichotomous variable has two categories i.e. passenger having sibling or spouse and passenger not having sibling or spouse. Similarly, Parch_yes variable represents whether passenger is having parent or children or not along with him. Female_yes variable represents whether passenger is male or female.

Logistic regression modeling results are as follows---

> summary(glm_model)

Call:
glm(formula = Survived ~ . - Embarked - Fare, family = "binomial",
    data = titanic_sub)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-2.7440  -0.6168  -0.4328   0.5865   2.6719

Coefficients:
                          Estimate    Std. Error      z value      Pr(>|z|) 
(Intercept)         2.090042   0.444176       4.705        2.53e-06 ***
Pclass              -1.093903    0.121616     -8.995          < 2e-16 ***
Age                 -0.037635    0.008076     -4.660         3.16e-06 ***
SibSp              -0.779355    0.202681     -3.845          0.00012 ***
Parch              -0.454346    0.210032     -2.163           0.03052 *
SibSp_yes       0.886692     0.331693       2.673           0.00751 **
parch_yes        0.903163    0.411734        2.194          0.02827 *
female_yes      2.718656    0.199072      13.657          < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1186.66  on 890  degrees of freedom
Residual deviance:  776.49  on 883  degrees of freedom
AIC: 792.49.

The model contains variables with significant regression coefficients. The following visualization shows how passengers are divided into survivor and non survivors---

The score axis tells  the probability of survival. Red points indicate the passengers we died in Titanic disaster. We can observe very well that red points are more concentrated towards 0 score which means chances of survival are very bleak or no chance of survival at all. 

The lift curve attached below also shows the accuracy of model classification. The line graph should be as distant as possible fro diagonal. It seems that classification model needs more fine tuning.

FINE TUNING THE MODEL

There are many ways to fine tune the model. Here, we are going to use random forest technique for fine tuning the model. In this technique, 500 decision tree models are created and than best model is selected which has minimum BIC and AIC values. Following result is obtained in random forest modeling---\

randomForest(formula = formula, data = df_for_model, importance = TRUE,      proximity = TRUE, eig = TRUE, ntree = 500, type = "classification") 

               Type of random forest: regression

                     Number of trees: 500

No. of variables tried at each split: 1

          Mean of squared residuals: 0.1411216

                    % Var explained: 40.33

When prediction is estimated for chances of survival of passenger, random forest model creates a vector of probabilities assigned with each passenger. The probabilities are not deterministic hence we need to find a cutoff score. This cutoff score will be a threshold. passengers with probabilities below the cutoff value will be assigned the value "Not Survived=0" and passengers with probabilities equal to or above cutoff value will be assigned the value "Survived=1". There are two ways to find optimal cutoff where maximum accuracy in classification can be obtained--- a) Iterative method and b) ROC method. In iterative method, the cutoff options are changes and model accuracy is calculated again and again. The cutoff is selected where maximum accuracy is obtained. In ROC method, algorithm automatically calculates the cutoff where accuracy is maximum. The ROC method is less time consuming hence we select ROC method to find out optimal cutoff and following results are obtained---

   

  Survived

         |        Predicted

        \/         0       1

         0        516   33

         1        107   235

Area under the curve: 0.8945 

Threshold       Accuracy 

0.3521408      0.8271605 

The area under the curve covers 89% of plot with threshold 0.35 and accuracy 83%. This model is better than logistic regression model. Following table shows the comparison of results----

The accuracy of random forest modeling is 84.29% while accuracy of logistic regression model is 77.89%. Random forest model is definitely chosen because of higher accuracy.

========================================================================

 R Script

=======================================================================

#####$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$#####
######################################################################
######################################################################
###              WHO DID SURVIVE THE TITANIC DISASTER?             ###
######################################################################
######################################################################
#####@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@#####


######################################################################
######################################################################
###                        GETTING STARTED                         ###
######################################################################
######################################################################


########################################################
## Loading the working directory

setwd()

getwd()

#########################################################
## Loading the libraries

library(dplyr)
library(e1071)
library(dplyr)
library(Hmisc)
library(ggplot2)
library(randomForest)
library(caret)
library(pROC)

#########################################################
## Loading the datasets

train=read.csv("train.csv", header = T, stringsAsFactors = F)
test=read.csv("test.csv", header=T, stringsAsFactors = F)

str(train)
str(test)

#########################################################
## Combining the datasets

train$DataType="Train"
test$DataType="Test"
test$Survived=NA

df=rbind(train,test)


######################################################################
######################################################################
###              MISSING VALUE ANALYSIS AND INPUTATIOn             ###
######################################################################
######################################################################

#########################################################
## Calculating the missing values for each column

MV_summary=sapply(df, function(x) sum(is.na(x)))
MV_summary

#########################################################
## Regression for estimating missing value in age variable

str(df)

model_age=lm(Age~Fare, data=df)

summary(model_age)

View(df)

View(test)

df$Sex=ifelse(df$Sex=="male",1,0)


model_age_2=lm(Age~+Pclass+Sex+SibSp+Parch, data=df)

summary(model_age_2)


df$Age[is.na(df$Age)] <- predict(model_age_2, na.rm = T, newdata=df)

df$Age=round(df$Age, digits = 0)

df$Fare=round(df$Fare, digits = 0)

MV_summary=sapply(df, function(x) sum(is.na(x)))
MV_summary

df$Fare[is.na(df$Fare)] <- mean(df$Fare, na.rm = T)

summary(df)

df[925,6]=2

which(df$Age<0)

which(df$Age==0)

which(df$Fare==0)

str(df)

######################################################################
######################################################################
###      Modeling of Response Variable with Predictors             ###
######################################################################
######################################################################



df_for_model=df[c(1:891),c(2,3,5,6,7,8,10)]


formula=as.formula(Survived~Age+Sex+SibSp+Pclass)

#########################################################
## Applying linear regression model


model_lm=lm(formula, data=df_for_model)

summary(model_lm)

#########################################################
## Applying logistic regression model

model_lg=glm(formula, data=df_for_model, family="binomial")

summary(model_lg)

df_for_model$Predicted=predict(model_lg, type="response", newdata = df_for_model)

View(df_for_model)

#########################################################
## Plotting the model predicted with actual variable

ggplot(df_for_model,aes(y=Survived,x=Predicted,color=factor(Survived)))+
  geom_point()+geom_jitter()

#########################################################
## Applying Random Forest technique for modeling

model_rf=randomForest(formula, data=df_for_model, importance=TRUE, proximity=TRUE,eig=TRUE, ntree=500, type="classification")

model_rf

summary(model_rf)

plot(model_rf)

df_for_model$predict_rf=predict(model_rf, type="response", newdata=df_for_model)

#########################################################
## Plotting the random forest model

ggplot(df_for_model,aes(y=Survived,x=predict_rf,color=factor(Survived)))+
  geom_point()+geom_jitter()

#########################################################
## Optimizing the cutoff value by interative procedure

cutoff=0.3521408

df_for_model$P0.10=ifelse(df_for_model$predict_rf<cutoff,0,1)

table(df_for_model$Survived, df_for_model$P0.10)

df_confusion_matrix=cbind(Survived=df_for_model$Survived,Predicted=df_for_model$P0.10)

df_confusion_matrix=as.data.frame(df_confusion_matrix)


df_confusion_matrix$Survived=as.factor(df_confusion_matrix$Survived)
df_confusion_matrix$Predicted=as.factor(df_confusion_matrix$Predicted)

confusionMatrix(df_confusion_matrix$Survived, df_confusion_matrix$Predicted)

model_roc=roc(df_for_model$Survived,df_for_model$predict_rf)

model_roc

plot(model_roc, print.thres="best", print.thres.best.method="closest.topleft", xlim=c(0,1))

result.coords <- coords(model_roc, "best", best.method="closest.topleft", ret=c("threshold", "accuracy"))
print(result.coords)

cutoff=0.3521408

df_for_model$Predicted=ifelse(df_for_model$Predicted<cutoff,0,1)

table(df_for_model$Survived, df_for_model$Predicted)