Diskriminantinė analizė
Tautvydas Merkevičius
2022-03-22
Užduotis:
Labai trumpai charakterizuoti duomenų rinkiniai, sprendžiamas klasifikavimo uždavinys abiems pasirinktiems duomenų rinkiniams. Pasirinktiems duomenims sudarykite LDA klasifikatorių pagal pavyzdį 128-131 psl. iš 5 skyriaus “Machine learning with R, tidyverse and mlr”. Padarykite vieną “hold-out” testą ir išveskite klasifikavimo lenteles “confusion matrix” mokymo ir testavimo duomenims.
Užduočiai pasiekti pasirenkami du skirtingi duomenų rinkiniai, patalpinti “UCI Machine Learning Repository” duomenų bazėje. Vienas iš jų – klasikinis rinkodaros banko duomenų rinkinys (“BANK”), kuriame yra turima informacija apie 45 000 asmenis. Kitame duomenų rinkinyje (“ADULT”) rasime informacijos apie maždaug 30 000 asmenų pajamas, kurių metines pajamas lemia įvairūs veiksniai: iššilavinimo lygis, amžius, lytis ir kita.
Reikalingos bibliotekos:
rm(list = ls())
library(caret)
library(mlr)
library(dplyr)
library(tidyr)
library(readxl)Pirmas duomenų rinkinys “bank”.
Nuskaitykime pirmą duomenų rinkinį ir jį pavadinkime “bankclient”. Duomenų rinkinį sudaro daugiau nei 45 000 eilučių ir 17 atributų. Ištrinamos nežinomos reikšmės, sudaroma duomenų lentelė ir “dplyr” paketo pagalba parenkamos klasifikavimo uždaviniui reikalingos reikšmės: amžius (age), paskutinė susisiekimo diena (day), sąskaitos likutis (balance), paskutinio skambučio trukmė (duration), kompanijos kontaktų skaičius (campaign), dienų skaičius kurios praėjo nuo paskutinio susisiekimo su klientu iš ankstesnės kompanijos (pdays), ankstesnis kliento kontaktų skaičius (previous), ar klientas sutartį pasirašo ar ne (y). Būtent šis uždavinys ir yra skirtas prognozuoti, ar klientas banke pasirašys terminuotą sutartį ar ne.
bankclient <- read.csv(file = 'C:/Users/Tautvydas/Desktop/DUOMENŲ GAVYBA/bank/bank-full.csv',
stringsAsFactors = FALSE,
strip.white = TRUE,
sep = ';')
data <- bankclient[complete.cases(bankclient), ] #Deleting NA values
str(data)## 'data.frame': 45211 obs. of 17 variables:
## $ age : int 58 44 33 47 33 35 28 42 58 43 ...
## $ job : chr "management" "technician" "entrepreneur" "blue-collar" ...
## $ marital : chr "married" "single" "married" "married" ...
## $ education: chr "tertiary" "secondary" "secondary" "unknown" ...
## $ default : chr "no" "no" "no" "no" ...
## $ balance : int 2143 29 2 1506 1 231 447 2 121 593 ...
## $ housing : chr "yes" "yes" "yes" "yes" ...
## $ loan : chr "no" "no" "yes" "no" ...
## $ contact : chr "unknown" "unknown" "unknown" "unknown" ...
## $ day : int 5 5 5 5 5 5 5 5 5 5 ...
## $ month : chr "may" "may" "may" "may" ...
## $ duration : int 261 151 76 92 198 139 217 380 50 55 ...
## $ campaign : int 1 1 1 1 1 1 1 1 1 1 ...
## $ pdays : int -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
## $ previous : int 0 0 0 0 0 0 0 0 0 0 ...
## $ poutcome : chr "unknown" "unknown" "unknown" "unknown" ...
## $ y : chr "no" "no" "no" "no" ...bankdata <- bankclient %>% dplyr::select(age, day, balance, duration, campaign, pdays, previous, y)
bankdata<- as_tibble(bankdata)
bankdata$y<- as.factor(bankdata$y)
str(bankdata)## tibble [45,211 x 8] (S3: tbl_df/tbl/data.frame)
## $ age : int [1:45211] 58 44 33 47 33 35 28 42 58 43 ...
## $ day : int [1:45211] 5 5 5 5 5 5 5 5 5 5 ...
## $ balance : int [1:45211] 2143 29 2 1506 1 231 447 2 121 593 ...
## $ duration: int [1:45211] 261 151 76 92 198 139 217 380 50 55 ...
## $ campaign: int [1:45211] 1 1 1 1 1 1 1 1 1 1 ...
## $ pdays : int [1:45211] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
## $ previous: int [1:45211] 0 0 0 0 0 0 0 0 0 0 ...
## $ y : Factor w/ 2 levels "yes","no": 2 2 2 2 2 2 2 2 2 2 ...Šiek tiek panagrinėkime turimą duomenų rinkinį. Matome, kiekvieno atributo minimalias, makslimalias reikšmes, kvartilius bei mediana. Kaip faktorius klasifikavimo uždaviniui pasirenkamas sutarties pasirašymo kintamasis, kuris gali įgyti taip (yes) – 1, ne (no) – 2.
summary(bankdata)## age day balance duration
## Min. :18.00 Min. : 1.00 Min. : -8019 Min. : 0.0
## 1st Qu.:33.00 1st Qu.: 8.00 1st Qu.: 72 1st Qu.: 103.0
## Median :39.00 Median :16.00 Median : 448 Median : 180.0
## Mean :40.94 Mean :15.81 Mean : 1362 Mean : 258.2
## 3rd Qu.:48.00 3rd Qu.:21.00 3rd Qu.: 1428 3rd Qu.: 319.0
## Max. :95.00 Max. :31.00 Max. :102127 Max. :4918.0
## campaign pdays previous y
## Min. : 1.000 Min. : -1.0 Min. : 0.0000 yes: 5289
## 1st Qu.: 1.000 1st Qu.: -1.0 1st Qu.: 0.0000 no :39922
## Median : 2.000 Median : -1.0 Median : 0.0000
## Mean : 2.764 Mean : 40.2 Mean : 0.5803
## 3rd Qu.: 3.000 3rd Qu.: -1.0 3rd Qu.: 0.0000
## Max. :63.000 Max. :871.0 Max. :275.0000bankUntidy <- gather(bankdata, "Variable", "Value", -y)
ggplot(bankUntidy, aes(y, Value)) +
facet_wrap(~ Variable, scales = "free_y") +
geom_boxplot() +
theme_bw()
Vadinamieji box plot grafikai, rodo grafinį skaitinių imties charakteristikų išsidėstymą, kurias galima palyginti. Stačiakampis brėžiamas nuo pirmojo iki trečiojo kvartilio. Stačiakampį į dvi dalis dalija brūkšnys, žymintis medianą. Nuo stačiakampio į viršų ir į apačią brėžiamos dvi tiesių atkarpos: apatinė – iki mažiausios imties reikšmės, o viršutinė – iki didžiausios imties reikšmės. Taip pat galima pastebėti išskirtis.
Sukurkime reikalingą modelį apsibrėždami klasifikatorių ir nurodykime mokymosi imtį. Šį kartą padarykime tai naudodami LDA (Linear Discriminant Analysis)
bankTask <- makeClassifTask(data = bankdata, target = "y")
lda <- makeLearner("classif.lda")
ldaModel <- train(lda, bankTask)Iš modelio išskaičiuokime reikalingą informaciją pasinaudodami funkcija getLearnerModel() ir gaukime duomenų lentelės reikšmes kiekvienu atveju naudojant prognozavimo funkciją. Matome, kad modelis išmoko vieną diskriminantinę funkciją (DF) – LD1. Prognozės funkcija grąžina šios funkcijos reikšmes kiekvienu atveju mūsų turimo “bankdata” duomenų rinkinio.
ldaModelData <- getLearnerModel(ldaModel)
ldaPreds <- predict(ldaModelData)$x
head(ldaPreds)## LD1
## 1 -0.0756109
## 2 0.5275920
## 3 0.8983154
## 4 0.6943838
## 5 0.4078993
## 6 0.6246096Sudarykime QDA (Quadratic Discriminant Analysis) modelį.
qda <- makeLearner("classif.qda")
qdaModel <- train(qda, bankTask)Patikrinkime, kaip sudaryti modeliai veikia k-fold metodu.
kFold <- makeResampleDesc(method = "RepCV", folds = 10, reps = 50,
stratify = TRUE)
ldaCV <- resample(learner = lda, task = bankTask, resampling = kFold,
measures = list(mmce, acc))
qdaCV <- resample(learner = qda, task = bankTask, resampling = kFold,
measures = list(mmce, acc))
ldaCV$aggr## mmce.test.mean acc.test.mean
## 0.1117954 0.8882046qdaCV$aggr## mmce.test.mean acc.test.mean
## 0.1273663 0.8726337LDA modelis teisingai klasifikavo vidutiniškai 89% sutarties pasirašymo prognozės, o QDA modeliui, pavyko teisingai klasifikuoti 87%. Išvada, kad LDA klasifikuoja geriau. Patikrinkime, kaip atrodys klasifikavimo lentelės (confusion matrix) ir įvertinkime kaip jie veiks turint naujus duomenis.
calculateConfusionMatrix(ldaCV$pred, relative = TRUE)## Relative confusion matrix (normalized by row/column):
## predicted
## true yes no -err.-
## yes 0.25/0.55 0.75/0.09 0.75
## no 0.03/0.45 0.97/0.91 0.03
## -err.- 0.45 0.09 0.11
##
##
## Absolute confusion matrix:
## predicted
## true yes no -err.-
## yes 66814 197636 197636
## no 55083 1941017 55083
## -err.- 55083 197636 252719newclient <- tibble(age = 40, day = 2, balance = 1000, duration = 10, campaign = 3, pdays = -1, previous = 0)
predict(ldaModel, newdata = newclient)## Prediction: 1 observations
## predict.type: response
## threshold:
## time: 0.00
## response
## 1 noModelis prognozuoja, kad naujas klientas terminuotos sutarties nepasirašys!
Papildomai atliekamas hold-out testas.
holdout <- makeResampleDesc(method = "Holdout", split = 2/3,
stratify = TRUE)
lda_holdoutCV <- resample(learner = lda, task = bankTask,
resampling = holdout, measures = list(mmce, acc))
qda_holdoutCV <- resample(learner = qda, task = bankTask,
resampling = holdout, measures = list(mmce, acc))
lda_holdoutCV$aggr## mmce.test.mean acc.test.mean
## 0.1092827 0.8907173qda_holdoutCV$aggr## mmce.test.mean acc.test.mean
## 0.1246765 0.8753235calculateConfusionMatrix(lda_holdoutCV$pred, relative = TRUE)## Relative confusion matrix (normalized by row/column):
## predicted
## true yes no -err.-
## yes 0.26/0.57 0.74/0.09 0.74
## no 0.03/0.43 0.97/0.91 0.03
## -err.- 0.43 0.09 0.11
##
##
## Absolute confusion matrix:
## predicted
## true yes no -err.-
## yes 467 1296 1296
## no 351 12957 351
## -err.- 351 1296 1647Matome, kad atlikus holdout testą LDA klasifikuoja geriau (89%).
Antras duomenų rinkinys “ADULT”
Šis duomenų rinkinys yra skirtas prognozuoti, ar asmens pajamos per metus yra didesnės už 50 000 dolerių ar mažesnės. Tai klasifikuoti padeda įvairūs faktoriai apie žmogų. Rinkinyje yra viso 14 atributų ir daugiau nei 32 000 duomenų eilučių (vartotojų). Šiam klasifikavimui naudojami duomenys nuskaitomi R kodu, ištrinamos nežinomos NA reikšmės, sudaroma duomenų lentelė.
adult <- read_excel("DUOMENŲ GAVYBA/adult.xlsx")
str(adult)## tibble [32,561 x 15] (S3: tbl_df/tbl/data.frame)
## $ age : num [1:32561] 39 50 38 53 28 37 49 52 31 42 ...
## $ workclass : chr [1:32561] "State-gov" "Self-emp-not-inc" "Private" "Private" ...
## $ fnlwgt : num [1:32561] 77516 83311 215646 234721 338409 ...
## $ education : chr [1:32561] "Bachelors" "Bachelors" "HS-grad" "11th" ...
## $ educationnum : num [1:32561] 13 13 9 7 13 14 5 9 14 13 ...
## $ maritalstatus: chr [1:32561] "Never-married" "Married-civ-spouse" "Divorced" "Married-civ-spouse" ...
## $ occupation : chr [1:32561] "Adm-clerical" "Exec-managerial" "Handlers-cleaners" "Handlers-cleaners" ...
## $ relationship : chr [1:32561] "Not-in-family" "Husband" "Not-in-family" "Husband" ...
## $ race : chr [1:32561] "White" "White" "White" "Black" ...
## $ sex : chr [1:32561] "Male" "Male" "Male" "Male" ...
## $ capitalgain : num [1:32561] 2174 0 0 0 0 ...
## $ capitalloss : num [1:32561] 0 0 0 0 0 0 0 0 0 0 ...
## $ hoursperweek : num [1:32561] 40 13 40 40 40 40 16 45 50 40 ...
## $ nativecountry: chr [1:32561] "United-States" "United-States" "United-States" "United-States" ...
## $ income : chr [1:32561] "<=50K" "<=50K" "<=50K" "<=50K" ...data1 <- adult[complete.cases(adult), ] #Deleting NA values
adultdata <- data1 %>% dplyr::select(age, fnlwgt, educationnum, capitalgain, capitalloss, hoursperweek, income)
adultdata<- as_tibble(adultdata)
adultdata$income<- as.factor(adultdata$income)
str(adultdata)## tibble [32,561 x 7] (S3: tbl_df/tbl/data.frame)
## $ age : num [1:32561] 39 50 38 53 28 37 49 52 31 42 ...
## $ fnlwgt : num [1:32561] 77516 83311 215646 234721 338409 ...
## $ educationnum: num [1:32561] 13 13 9 7 13 14 5 9 14 13 ...
## $ capitalgain : num [1:32561] 2174 0 0 0 0 ...
## $ capitalloss : num [1:32561] 0 0 0 0 0 0 0 0 0 0 ...
## $ hoursperweek: num [1:32561] 40 13 40 40 40 40 16 45 50 40 ...
## $ income : Factor w/ 2 levels "<=50K",">50K": 1 1 1 1 1 1 1 2 2 2 ...summary(adultdata)## age fnlwgt educationnum capitalgain
## Min. :17.00 Min. : 12285 Min. : 1.00 Min. : 0
## 1st Qu.:28.00 1st Qu.: 117827 1st Qu.: 9.00 1st Qu.: 0
## Median :37.00 Median : 178356 Median :10.00 Median : 0
## Mean :38.58 Mean : 189778 Mean :10.08 Mean : 1078
## 3rd Qu.:48.00 3rd Qu.: 237051 3rd Qu.:12.00 3rd Qu.: 0
## Max. :90.00 Max. :1484705 Max. :16.00 Max. :99999
## capitalloss hoursperweek income
## Min. : 0.0 Min. : 1.00 <=50K:24720
## 1st Qu.: 0.0 1st Qu.:40.00 >50K : 7841
## Median : 0.0 Median :40.00
## Mean : 87.3 Mean :40.44
## 3rd Qu.: 0.0 3rd Qu.:45.00
## Max. :4356.0 Max. :99.00Naudojami atributai: vartotojo amžius (age), pajamų svoris (fnlwgt), mokymosi metai (educationnum), kapitalo prieaugis (capitalgain), kapitalo nuostolis (capitalloss), darbo valandos per savaitę (hoursperweek), pajamos (income). Visi kintamieji kiekybiniai, išskyrus pajamas, kurios gali įgauti dvi reikšmes: arba >50K arba <50K.
Atvaizduojami boxplot grafikai, kad grafiškai būtų matomos įvairios duomenų matematinės charakteristikos.
adultUntidy <- gather(adultdata, "Variable", "Value", -income)
ggplot(adultUntidy, aes(income, Value)) +
facet_wrap(~ Variable, scales = "free_y") +
geom_boxplot() +
theme_bw()
Apibrėžiame pajamų klasifikavimo uždavinį ir sukuriami LDA ir QDA klasifikatoriai. Sukurti klasifikatoriai yra apmokomi.
adultTask <- makeClassifTask(data = adultdata, target = "income")## Warning in makeTask(type = type, data = data, weights = weights, blocking =
## blocking, : Provided data is not a pure data.frame but from class tbl_df, hence
## it will be converted.lda <- makeLearner("classif.lda")
ldaModel <- train(lda, adultTask)
ldaModelData <- getLearnerModel(ldaModel)
ldaPreds <- predict(ldaModelData)$x
head(ldaPreds)## LD1
## 1 0.6933577
## 2 0.1397711
## 3 -0.4135499
## 4 -0.4042299
## 5 0.3283618
## 6 0.8652966qda <- makeLearner("classif.qda")
qdaModel <- train(qda, adultTask)Sudaryto modelio validavimas pasinaudojant k-fold metodu.
kFold <- makeResampleDesc(method = "RepCV", folds = 10, reps = 50,
stratify = TRUE)
ldaCV <- resample(learner = lda, task = adultTask, resampling = kFold,
measures = list(mmce, acc))
qdaCV <- resample(learner = qda, task = adultTask, resampling = kFold,
measures = list(mmce, acc))
ldaCV$aggr## mmce.test.mean acc.test.mean
## 0.1975633 0.8024367qdaCV$aggr## mmce.test.mean acc.test.mean
## 0.2036669 0.7963331LDA modelis teisingai klasifikavo vidutiniškai 80.2% pajamų prognozės, o QDA modeliui, pavyko teisingai klasifikuoti 79.6%. Išvada, kad LDA klasifikuoja geriau.
Patikrinkime, kaip atrodys klasifikavimo lentelės (confusion matrix) ir įvertinkime kaip jie veiks turint naujus duomenis.
calculateConfusionMatrix(ldaCV$pred, relative = TRUE)## Relative confusion matrix (normalized by row/column):
## predicted
## true <=50K >50K -err.-
## <=50K 0.95/0.82 0.05/0.31 0.05
## >50K 0.67/0.18 0.33/0.69 0.67
## -err.- 0.18 0.31 0.20
##
##
## Absolute confusion matrix:
## predicted
## true <=50K >50K -err.-
## <=50K 1177649 58351 58351
## >50K 263292 128758 263292
## -err.- 263292 58351 321643newadult <- tibble(age = 40, fnlwgt = 0, educationnum = 1, capitalgain = 0, capitalloss = 1000, hoursperweek = 10)
predict(qdaModel, newdata = newadult)## Prediction: 1 observations
## predict.type: response
## threshold:
## time: 0.00
## response
## 1 <=50KGauname, kad naujo kliento pajamos bus mažesnės arba lygios 50 000 dolerių.
Papildomai atlikime hold-out testą:
holdout <- makeResampleDesc(method = "Holdout", split = 2/3,
stratify = TRUE)
lda_holdoutCV <- resample(learner = lda, task = adultTask,
resampling = holdout, measures = list(mmce, acc))
qda_holdoutCV <- resample(learner = qda, task = adultTask,
resampling = holdout, measures = list(mmce, acc))
lda_holdoutCV$aggr## mmce.test.mean acc.test.mean
## 0.2013083 0.7986917qda_holdoutCV$aggr## mmce.test.mean acc.test.mean
## 0.2057306 0.7942694calculateConfusionMatrix(lda_holdoutCV$pred, relative = TRUE)## Relative confusion matrix (normalized by row/column):
## predicted
## true <=50K >50K -err.-
## <=50K 0.95/0.81 0.05/0.33 0.05
## >50K 0.68/0.19 0.32/0.67 0.68
## -err.- 0.19 0.33 0.20
##
##
## Absolute confusion matrix:
## predicted
## true <=50K >50K -err.-
## <=50K 7836 404 404
## >50K 1781 833 1781
## -err.- 1781 404 2185Matome, kad atlikus holdout testą geriau klasifikuoja LDA (80.4%).