MATH2349 Semester 1, 2018

Assignment 3

Required packages

Provide the packages required to reproduce the report. Make sure you fulfilled the minimum requirement #10.

library(readr)
library(tidyr)
library(dplyr)
library(lubridate)
library(MVN)
library(car)
library(forecast)
library(caret) 

Executive Summary

To achieve the high accuracy and precision of Machine learning models, data needs to be in clean format. To prepare such data there are number of tasks involved in it. The tasks includes data cleaning, data transformation, dealing with missing values and outliers and correct data inconsistencies. Another important factor is that the data has to be in tidy format.The data has around 19% of the missing values with some of the columns having more than 50% of the data is missing.Caaping method used to deal with the outliers present in the data.Furthermore, there is a huge variation between the ranges of the varibles which normalised using BoxCox transformation.

Data

Two datasets has been selected from kaggle website.[7] The data is about bank loan application. First file(application_train), contains the static data of all applications, one row represents one loan application. Second file(bureau) contains the applicants previous credits provided by other financial institutes & that were reported to credit bureau(for clients who have a loan application in first file). For every record in first file there are multiple records in second file.below are the steps to get the data. * Due to large number of volume in both the file, the data has been filtered on the basis of “SK_ID_CURR” column which is common in both the files and stored in separate CSV files to perform data preprocessing tasks.

##original source files
loan_app_source<-read.csv("application_train.csv")
bureau_source<-read.csv("bureau.csv")

##filter the data
loan_app_filtered<- loan_app_source %>% filter(SK_ID_CURR<=101156)
bureau_filtered<-bureau_source %>% filter(SK_ID_CURR<=101156)

#export filtered data
write.csv(loan_app_filtered, file = "loan_app_filt.csv",row.names=FALSE)
write.csv(bureau_filtered, file = "bureau_filt.csv",row.names=FALSE)

• In above chunk data has been filtered using “SK_ID_CURR<=101156” condition which will contain approximately 10% of the source data.

• Filtered data has been exported in separate CSV file and then imported to perform data preprocessing tasks.

• Using left join merged these two data frames

loan_app<-read.csv("loan_app_filt.csv")
bureau<-read.csv("bureau_filt.csv")
#combined two datasets
combined<- loan_app %>% left_join(bureau,by="SK_ID_CURR") 

Understand

To understand the data it is important to check the dimensions and data types of all the variables.

• Used dim function, which shows that the dataset has 4647 observations and 138 variables.
• Used sapply function to get the data type of all variables.

#check the dimention of the data frame
dim(combined)

[1] 4496  109

#check the data types
sapply(combined[1:20],class) # showing only first 20 variables as the compelete lsit long

                SK_ID_CURR                     TARGET         NAME_CONTRACT_TYPE 
                 "integer"                  "integer"                   "factor" 
               CODE_GENDER               FLAG_OWN_CAR            FLAG_OWN_REALTY 
                  "factor"                   "factor"                   "factor" 
              CNT_CHILDREN           AMT_INCOME_TOTAL                 AMT_CREDIT 
                 "integer"                  "numeric"                  "numeric" 
             AMT_ANNUITY.x            AMT_GOODS_PRICE            NAME_TYPE_SUITE 
                 "numeric"                  "numeric"                   "factor" 
          NAME_INCOME_TYPE        NAME_EDUCATION_TYPE         NAME_FAMILY_STATUS 
                  "factor"                   "factor"                   "factor" 
         NAME_HOUSING_TYPE REGION_POPULATION_RELATIVE                 DAYS_BIRTH 
                  "factor"                  "numeric"                  "integer" 
             DAYS_EMPLOYED          DAYS_REGISTRATION 
                 "integer"                  "integer" 

head(combined,5)

The results of the above chunk shows that the data fulfills the minimum requirements 2-4

• requirement 2- the data included multiple datatypes int, factor, numeric.
• requirement 3- contains variables on which data type conversion can be applied.numeric->factor and numeric-> date (code is in next chunk)
• requirement 4->contains factor variable that needs to be ordered(code is in next chunk)

#datatype conversion
#numeric->factor
combined$FLAG_MOBIL<-factor(combined$FLAG_MOBIL)
combined$FLAG_EMP_PHONE<-factor(combined$FLAG_EMP_PHONE)
combined$FLAG_WORK_PHONE<-factor(combined$FLAG_WORK_PHONE)
#numeric->date
a<-ddays(combined$DAYS_BIRTH)
##assume the date on whoch data is genereted is 18th May 2018
birth_date<-as.Date("2018-05-18")+a
class(birth_date)

[1] "Date"

#factor variable from unordered to ordered
##check the levels of factor variable
levels(combined$WEEKDAY_APPR_PROCESS_START)

[1] "FRIDAY"    "MONDAY"    "SATURDAY"  "SUNDAY"    "THURSDAY"  "TUESDAY"   "WEDNESDAY"

#order the levels
combined$WEEKDAY_APPR_PROCESS_START<- ordered(combined$WEEKDAY_APPR_PROCESS_START,levels=c("MONDAY","TUESDAY","WEDNESDAY","THURSDAY","FRIDAY","SATURDAY","SUNDAY"))
#re check whether the levels has been ordered or not
levels(combined$WEEKDAY_APPR_PROCESS_START)

[1] "MONDAY"    "TUESDAY"   "WEDNESDAY" "THURSDAY"  "FRIDAY"    "SATURDAY"  "SUNDAY"   

• numeric->factor: the flag columns had values either 0 or 1. After converting it into factor it created 2 levels
• numeric-> date: DAYS_BIRTH column has client’s age in number of days. using ddays function from lubridate package and subtracting the value from current date gives date of birth of the client which is further used to convert it into age.
• ordered factor: weekday_appr_process_start.x column has 7 levels contains the weekdays in unordered format. using ordered function defined the order fo the levels as Monday<Tuesday<Wednesday<Thursday<Friday<Saturday<Sunday.

Tidy & Manipulate Data I

The data is in tidy format as it shows the characteristics of tidy data as below

• Every value in a data belongs to an observation and a variable.
• Variable contains values that measure the same attribute across unit.(like income_total,amount_credit)
• Every variable is column as able to extract needed variables. Created new variables using existing variable was smooth and possible because of tidy data structure.

Tidy & Manipulate Data II

To arrive at a decision of whether to approve the credit or to disapprove two more important predictors can be created from existing data, AGE and AMOUNT remain in hand after paying any ongoing credit. Hence tried creating new variables on the basis of available information.

#create new column in data frame"AGE" wchich will store th client's age
combined$AGE<-round(interval(start = birth_date, end = "2018-05-18") / 
  duration(num = 1, units = "years"))
#create new variable which will show the remaining annual amount from income after decucting the loan amount
combined<-transform(combined,INCOME_MINUS_ANUITY= AMT_INCOME_TOTAL- AMT_ANNUITY.x)
cols<-c('SK_ID_CURR','DAYS_BIRTH','AGE','AMT_INCOME_TOTAL','AMT_ANNUITY.x','INCOME_MINUS_ANUITY')
tail(combined[cols],5)

• Created new column ‘AGE’ which is computed on the basis of birth date computed in above chunk. column shows the age of the client in years

• INCOME_MINUS_ANNUITY shows the amount remaining after paying the loan.

Scan I

Below chunk shows the amount of missing values in data and how to deal with them.

#missing values scan
# number of columns having missing values
missing_val_cols<-colnames(combined)[colSums(is.na(combined)) > 0]
length(missing_val_cols)

[1] 74

#total percentage of missing values
sum(is.na(combined)) / (nrow(combined) *ncol(combined))

[1] 0.1965631

#colSums(is.na(combined)) # as the result is long 
# Handeling missing values
#remove the columns where there are more than 50% values are missing
mis_cols_large<-colnames(combined)[colSums(is.na(combined)) > 2200] 
combined<- combined %>% select(-one_of(mis_cols_large))
dim(combined)

[1] 4647  109

#replace missing values with the Median
for(i in which(sapply(combined,is.numeric))){
  combined[is.na(combined[,i]), i] <- median(combined[,i], na.rm = TRUE)
}
missing_val_cols<-colnames(combined)[colSums(is.na(combined)) > 0]
length(missing_val_cols)

[1] 3

dim(combined)

[1] 4647  109

#still all the factor column contains 151 missing values
combined<-combined[complete.cases(combined),]
print('complete cases ')

[1] "complete cases "

dim(combined)

[1] 4496  109

• Out of 140 variables 74 variables have missing values.
• Total percentage of missing values is 19.6.
• colSums shows that some of the columns have more than 50% values are missing.
• removed columns(43) having missing value count more than 50%. remained with 109 variables.
• for remaining columns imputed with the median of the columns.
• still there were three factor columns have missing values. *Used complete.cases()to ignore 151 NA value rows[1]

#inconsistancies check
colnames(combined)[apply(combined, 2, function(x) any(is.infinite(x)))=='TRUE']

character(0)

colnames(combined)[apply(combined, 2, function(x) any(is.finite(x)))=='TRUE']

character(0)

colnames(combined)[apply(combined, 2, function(x) any(is.nan(x)))=='TRUE']

character(0)

##check for edit rules

Scan II

To perform outlier detection created subset with 11 numeric columns.

##get the numeric columns
num_cols <- unlist(lapply(combined, is.numeric)) 
combined_sub<-combined[ , num_cols] ##90 NUmeric variables
# creating subset with 12 features to perfrom outliers detection and data normalization
cols<-c("AMT_INCOME_TOTAL","AMT_CREDIT","AMT_ANNUITY.x","AMT_GOODS_PRICE","AGE",
        "DAYS_EMPLOYED","EXT_SOURCE_2","EXT_SOURCE_3","APARTMENTS_AVG",
        "TOTALAREA_MODE","DAYS_CREDIT_ENDDATE")
combined_sub<-combined_sub[cols]
combined_sub<-combined_sub[complete.cases(combined_sub),]
# detect ouliers
par(mar=c(10,6,4,1)+.2)
Boxplot(combined_sub,las=2) 

  [1] "1521" "1522" "1523" "1872" "225"  "226"  "227"  "228"  "229"  "267"  "872"  "1627"
 [13] "1628" "1629" "1630" "1631" "1632" "1633" "1634" "1635" "273"  "274"  "275"  "276" 
 [25] "4052" "4053" "4054" "872"  "4031" "4032" "872"  "1627" "1628" "1629" "1630" "1631"
 [37] "1632" "1633" "1634" "1635" "1255" "1256" "1257" "1258" "1259" "1260" "1261" "1262"
 [49] "1263" "1264" "40"   "41"   "42"   "43"   "53"   "54"   "55"   "56"   "98"   "99"  
 [61] "2895" "2896" "2897" "2898" "2899" "2900" "2901" "2902" "4607" "3994" "2592" "2593"
 [73] "2594" "2595" "2596" "4163" "4629" "4630" "4631" "1868" "2573" "2574" "2575" "2576"
 [85] "2577" "2578" "2579" "2580" "2581" "2582" "1223" "1224" "1225" "1226" "1227" "1228"
 [97] "1229" "1230" "1231" "4629" "2573" "2574" "2575" "2576" "2577" "2578" "2579" "2580"
[109] "2581" "2582" "3185" "3432" "4063" "3451" "2928" "4525" "700"  "933"  "657"  "4159"

• According to Boxplot there are 120 outliers in these 12 variables.
• Most of the outliers observed in three columns.
• To avoid simply removing these outliers records from data entirely used capping function to replace the values above 95th percentile and below 5th percentile.[2]

#Used capping method to replace outliers in the 
cap <- function(x){
  quantiles <- quantile( x, c(.05, 0.25, 0.75, .95 ) )
  x[ x < quantiles[2] - 1.5*IQR(x) ] <- quantiles[1]
  x[ x > quantiles[3] + 1.5*IQR(x) ] <- quantiles[4]
  x
}
comb_sub_wo_out<-apply(combined_sub, 2, FUN=cap)
#boxplot after outlier capping
par(mar=c(10,6,4,1)+.2)
Boxplot(comb_sub_wo_out,las=2)

 [1] "40"  "41"  "42"  "43"  "53"  "54"  "55"  "56"  "98"  "99"  "1"   "2"   "3"   "4"  
[15] "5"   "6"   "7"   "8"   "78"  "79"  "64"  "65"  "66"  "67"  "68"  "69"  "98"  "99" 
[29] "100" "101" "1"   "2"   "3"   "4"   "5"   "6"   "7"   "8"   "78"  "79"  "64"  "65" 
[43] "66"  "67"  "68"  "69"  "98"  "99"  "100" "101"

• The outliers visible in previous box plot has been capped.

Transform

First, selected five variables to see the distribution of them in the data.Plotted histograms to see the distribution.

#check the need for normalization
colnames <- colnames(comb_sub_wo_out)
par(mfrow = c(3, 2))
for (i in 1:5) {
  hist(unlist(comb_sub_wo_out[,i]), main = colnames[i])
}

• it is clear from the histograms the range of the data for each variable is not same. In fact, there is a huge difference. For example, age variable has range from 20-17 whereas AMT_GOOS_PRICE has range from 0-1400000.
• If we run the model on this data, results will be biased.
• Hence, need to normalize the data.
• Used BoxCox transformation to normalize the data.

boxcoxnorm <- apply(comb_sub_wo_out, 2, function(x) BoxCox(x, lambda = "auto"))

NA/Inf replaced by maximum positive value

#histograms after normalisation
par(mfrow = c(3, 2))
for (i in 1:5) {
  hist(unlist(boxcoxnorm[,i]), main = colnames[i])
}

• It is clear from the above histograms that the data has been normlised.
• The scale of the data has been changed.
• there is not much difference between the ranges for every variable.[3]

Feature Selection and Extraction

The dataset consists of large set of variables and most of the variables are correlated hence there si need to find few important variables. Principal component analysis is one of the method which extract important variables from large set.

• To perform PCA data should be numeric and standardized which we did in previous chunk.
• Since, the data which will be used by model consists of 79 numeric variables. i.e more than 1000 plots to analyze the relationship between variables.
• selected only numeric columns of the data
• identified and removed zero variance columns from data as prcomp cannot rescale a constant/zero column to unit variance.
• removed ID related columns from data as they no make sense when it comes to PCA.
• used prcomp to apply PCA on the data, which by default centers the variable to have mean equals to zero. scale.=T parameter normalize the variables to have SD equals to 1.
• This function results five useful measures(“sdev”,“rotation”,“center”,“scale”,“x”)[4,6]

#select only numeric columns
num_cols <- unlist(lapply(combined, is.numeric)) 
combined_mod<-combined[ , num_cols] ##4496
combined_mod <- combined_mod[complete.cases(combined_mod),]
#remove zero variance columns from data
combined_mod<-combined_mod[ , apply(combined_mod, 2, var) != 0]
#remove ID related columns form the data
combined_mod<-combined_mod %>% select(-c(SK_ID_CURR,SK_ID_BUREAU))
dim(combined_mod)

[1] 4496   76

#apply PCA
pca_combined <- prcomp(combined_mod, scale. = T)
#compute standard deviation of each principal component
std_dev <- pca_combined$sdev
#compute variance
pr_var <- std_dev^2
#check variance of first 10 components
pr_var[1:10]

 [1] 9.971499 4.350355 3.673241 3.291408 2.884132 2.759728 2.538433 2.322973 1.949930
[10] 1.786376

#aim is to find maximum variance.higher is the variance, higher will be the information.
#proportion of variance explained
prop_varex <- pr_var/sum(pr_var)
prop_varex[1:10]

 [1] 0.13120394 0.05724151 0.04833212 0.04330800 0.03794911 0.03631220 0.03340044 0.03056543
 [9] 0.02565698 0.02350495

• The proportion of variance shows that the 1st Component explains 12% variance, 2nd component explains 6.3 % variance and so on.
• To decide how many componenets should be selected used scatter plot.which shows the factors which explains the most of variability.[4,6]

par(mfrow = c(1, 2))
plot(prop_varex, xlab = "Principal Component",
     ylab = "Proportion of Variance Explained",
     type = "b")
#cumulative scree plot
plot(cumsum(prop_varex), xlab = "Principal Component",
       ylab = "Cumulative Proportion of Variance Explained",
       type = "b")

The plots above shows that ~58 components exlains around 98% variance.Same is confirmed with cumulative varaince plot.Therefore, the number of components as 58 can be used to perform modelling.[6]

References:

1.Week6 Demo slides,class worksheets

2.Week7 Demo slides,class worksheets

3.Week9 Demo slides,class worksheets

4.Week10 Demo slides,class worksheets

5.www.rdocumentation.org

6.https://www.analyticsvidhya.com/blog/2016/03/practical-guide-principal-component-analysis-python/

7.https://www.kaggle.com/c/home-credit-default-risk/data