How to ML-fy Your Small Business

ML is everywhere and those who do not know it, it is Machine Learning. The Machine Learning is out of the labs and being democratized by technologies. So, it can be used by all and sundry and it is not costly too. It is available in the micro-fraction of the cost it was used to be. Objective of this article is to understand how a small business can reap the benefits of machine learning and AI without hiccups. 

Before delving into machine learning and AI for small business, we must also know why a small business should even care about it? Can't it live without it. Yes, definitely a business can survive and be profitable without AI. It is a matter of competitive advantages if a company wants to have. AI will make small businesses what online eCommerce has done to the marketing function of small businesses. The eCommerce has uplifted the marketing and sales performance of hundreds and thousands of small businesses through pan-country reach and delivery systems. Now, any small business can take rest as far as marketing is concerned because eCommerce companies can take your products to consumers, they do marketing, they convince customers, they run sales promotion schemes....... almost everything, isn't it. Let us take an instance, what can you do if your competitor are also on same platform, selling same stuff, at same or lower price than you? Can you lower the price by capital financing the discounts, or delay the break-even or move out of the market. All these options are very scary. 

Here, ML-fy your business and this is what you can do as a logical step to sustain, survive or grow your business in fierce competitions. ML-fying your business means you delve into data generated in your small businesses and explore the opportunities of maximization and minimization. With the information of maximization and minimization opportunities, small business can foresee future business environment and internal performance. If that is known, you can have better planning, better approach to deal with competition. Over and above, your customer will feel the difference which means more customer loyalty and more referrals.

Now, we have got feeling of how ML-fication of small business can benefits in terms of minization and maximization opportunities.Following are the broad steps in implementing machine learning into business--

Let us understand this process of machine learning implementation with an example. Suppose there a sweet manufacturing company-- Sheetal Sweets. This company manufactures sweet "RASGULLA". It is a family owned business. This company is facing a problem that she has issues of either stockout situation or over-inventory of raw materials especially sugar and it is very frequent. Stockout situation of sugar leads to halt in manufacturing process and delay in sales order fulfillment. Oversupply leads to blocking of working capital which means interest has to be paid to bank on overdraft facility.

Let us dive deeper into the problem. Suppose working capital employed in business is Rs.30 lac. Bank charges interest rate @ 12% p.a. If working capital is paid in 6 months on an average, the interest burden on company is Rs.180000 for six month. On yearly basis, it is Rs.360000. Even if it is not borrowed from bank, we can assume the opportunity cost of employing the working capital. If Rs.30 lac amount is available with Sheetal Sweets, it should earn at least 12% return on business operations because bank might be offering same interest rates on bank deposits. Let us assume that working capital is completely absorbed in RASGULLA manufacturing.

The stockouts and over inventory adversely impact the working capital management. In case of stockouts, Sheetal Sweets may loose customer or sales and if there is over-inventory, the working capital rotation is delayed. So ultimately, either way it is harmful. Thus, Sheetal Sweets needs the predictability capability so that it can plan its inventory in optimal manner that minimizes the stockouts and over-inventory situations.

PREDICTABILITY CAPABILITY DEVELOPMENT

For predictability capability can be built through various measures but our focus in this article is to use machine learning for predicting the next day sales, next week sales, next fortnight sales and next month sales. There are range of models in machine learning that can predict the sales in coming periods but we will focus only on Time Series method. In time series method of forecasting, the variable data like sales is available for past periods and with suitable model, we try to predict the sales for coming periods.

                                   Picture 2 Process of Forecasting Model Development

The picture 2 shows the development process of forecasting model. Sheetal Sweets needs to collect its sales data for past periods. It is better to have daily sales data of at least 3 years or more. After collecting data, it should be stored into nice format in MS Excel. In Excel file, it sales data should be filled in cells corresponding to the periods they belong to.

Now, a small programming code in R can be used to build a forecasting model. The forecasting model programming in R can be asked separately via dropping the requirement in comment box. This model should be reviewed time to time to see how close the predicted values of model are to the actual sales realization. The lesser the gap, the better the model is.

Self Learning Technologies for eCommerce

The eCommerce business cannot survive without incorporation of Self Learning Technologies (SLTs) due to changes in scope of decision making from years and months to hours and minutes. Self Learning Technologies can be adopted with budget as minimum as $10000 and budget as big as $1000000. 

There are various ways to implement SLTs from as simple as hiring a freelancer on temporary engagement basis, crowdsourcing, renting the Machine Learning-as-a-Service (MLaaS) and as complex as buying SLT solution, ML cloud services and in-house SLT development. The allocation of human resources depends on approach of SLT implementation chosen. For simple SLT implementation, a data scientist and digital marketer may suffice but for complex SLT implementation, a dedicated team of data engineers, data scientists, business analytics professional and digital marketers should be deployed. 

The SLT deployment may take three months to one year or more which again depends on what kind of SLTs are being implemented. Simple and intermediate SLTs implementation may take three months to nine months and complex SLTs implementation may take twelve months to eighteen months. The decision of SLT implementation is very critical. 

For simple to intermediate SLT implementation, crowdsourcing and MLaaS are better options and for complex SLT implementation, renting MLaaS to buying SLT Solutions from Microsoft, Google, IBM etc. or build SLT solution in-house are suggested. If company wants to keep the SLT intellectual property confidential, limited to itself and to do full blown business transformation; it should go build SLT solution in-house or if company wants to deploy SLT implementation faster and battle tested SLT solution without actually owning it; the company can choose buy SLT solutions from reputed vendors.    

For more detail on how to implement self learning technologies (SLTs) in eCommerce, please click link below---

https://imojo.in/e9s5oc

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.