Maskinl?ring i R med tidymodels
Luigi Maglanoc
1 Intro
Bygge maskinl?ringsmodeller med tidymodels.
Tidymodels er en rekke pakker basert p? tidyverse prinsipper.
Tidyverse er i seg selv en samling av pakker for ?
for ? utf?re data analyse i R (klikk her).
F?rst m? vi laste inn data. Pakken “datasets” inneholder en rekke datasett
library(datasets)
data(iris)
datasettet best?r av 4 ulike (blomster) egenskaper for 150 blomster fra iris planteslekten.
Vi kan kikke p? de f?rste 6 radene med head kommandoen
egenskaper = variabler i en maskinl?ringskontekst
head(iris)
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
#head(penguins)
Her har vi lengden og bredden til “Sepal” (bergerblad), samt lengden og bredden til “petal” (kronblad)
Vi kan se p? fordelingen av blomsterartene med table funksjonen
table(iris$Species)
##
## setosa versicolor virginica
## 50 50 50
Det er alts? tre arter av irisblomsten: setosa, virginica og versicolor
Vi kan bruke f?lgende kommando for ? se p? datatypene til variablene
str(iris)
## 'data.frame': 150 obs. of 5 variables:
## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
M?let v?rt er ? bygge en maskinl?ringsmodell, eller klassifiseringsmodell
som kan skille mellom de ulike irisartene basert p? de 4 egenskapene
2 Exploratory data Analysis (EDA)
Det man b?r gj?re f?r man bygger selve maskinl?ringsmodellene er ? utf?re det noen kaller for
Eksplorerende dataanalyse, som er rett ? slett ? unders?ke ulike aspekter ved datasettet,
enten med beskrivende statistikk og/eller datavisualisering.
La oss plotte ulike aspekter ved dataen for ? se om det er tegn p? forskjeller i artene basert p? egenskapene
Vi kan starte med “Sepal.length”
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ? ggplot2 3.3.2 ? purrr 0.3.4
## ? tibble 3.1.2 ? dplyr 1.0.4
## ? tidyr 1.1.2 ? stringr 1.4.0
## ? readr 1.3.1 ? forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
ggplot(data=iris, aes(x=Species, y=Sepal.Length))+
geom_boxplot()
Det er enkelt ? fargelegge boksene basert p? artene
ggplot(data=iris, aes(x=Species, y=Sepal.Length, fill=Species))+
geom_boxplot()
Det er ganske tydelig forskjeller i “Sepal.length” blant de tre artene. Vi ser st?rst lengde p? “virginica” og minst lengde p? “setosa”.
La oss bytte ut “Sepal.length” med “Sepal.Width”
ggplot(data=iris, aes(x=Species, y=Sepal.Width, fill=Species))+
geom_boxplot()
Vi kan gjenta samme kode for ? individuelt plotte de ?vrige to variablene, men det er ogs? mulig
? plotte alt i samme figur. Dette krever at vi omstrukturer data fra “long” til “wide” format.
Det er flere m?ter ? gj?re det i R, men siden vi forholder oss til tidyverse kommer vi til ?
ta i bruk pivot_longer funksjonen fra tidyr
I “wide” format har vi en kolonne per variabel og en rad per observasjon.
Dvs. jo flere variabler vi har, jo lengre datarammen blir.
I “long” format samler vi dataen, slik at vi har en kolonne som angir hvilken variabel det gjelder,
og en annen “kolonne” som angir verdien til variabelen.
Antall rader er alts? antall observasjon ganger antall variabler vi
velger ? omstrukturere
facet_plot <- iris %>%
tidyr::pivot_longer(!Species, names_to = "variable", values_to = "value")
facet_plot
## # A tibble: 600 x 3
## Species variable value
## <fct> <chr> <dbl>
## 1 setosa Sepal.Length 5.1
## 2 setosa Sepal.Width 3.5
## 3 setosa Petal.Length 1.4
## 4 setosa Petal.Width 0.2
## 5 setosa Sepal.Length 4.9
## 6 setosa Sepal.Width 3
## 7 setosa Petal.Length 1.4
## 8 setosa Petal.Width 0.2
## 9 setosa Sepal.Length 4.7
## 10 setosa Sepal.Width 3.2
## # … with 590 more rows
N? kan vi dytte iris datarammen inn i ggplot.
ggplot(data=facet_plot, aes(x=variable, y=value, fill=Species))+
geom_boxplot()
Denne figuren hadde v?rt litt uoversiktlig hvis det var mange variabler, og det var store forskjeller
i verdiene.
En l?sning er ? dele figure opp i subplots med facet_wrap
ggplot(data=facet_plot, aes(x=variable, y=value, fill=Species))+
geom_boxplot()+
facet_wrap(~variable, scales="free")
Dette gir oss en ganske grei oversikt p? forskjellene i blomsterartene" basert p? “egenskapene”, og en anelse p? hvilke egenskaper kan v?re nyttige i klassifiseringen.
Her ser det ut til at alle kan v?re gode egenskaper,
Vi kan ogs? unders?ke? forholdet mellom egenskapene med “scatterplots”
Dett er en “base-R” funksjon (dvs. ikke en del av ggplot2 pakken).
pairs(iris[,1:4], pch = 19, lower.panel = NULL)
H?yt korrelasjon blant egenskapene (kolinearitet) er ikke ideelt for maskinl?ring.
En m?te ? h?ndtere dette p? er ? fjerne en av egenskapene som er h?yt korrelert.
3 dataforberedelse (klassifisere to kategorier)
For ? gj?re ting litt enklere s? kan vi starte med ? klassifisere to av blomsterartene, s?kalt bin?r klassifisering
iris_bin <- iris %>%
filter(Species != "setosa")
table(iris_bin$Species)
##
## setosa versicolor virginica
## 0 50 50
For ? forhindre problemer med tidymodels b?r vi fjerne faktorniv?et fra selve kolonnen
levels(iris_bin$Species)
## [1] "setosa" "versicolor" "virginica"
iris_bin$Species <- factor(iris_bin$Species)
table(iris_bin$Species)
##
## versicolor virginica
## 50 50
4 Trening- og testsett
Det man vanligvis gj?r f?r selve maskinl?ringssteget er ? dele opp datasettet i en trening og testsett.
Maskinl?ringsmodellene bygges p? treningssettet, og testsettet brukes til prediksjon/modellevaluering
Oppdelingen er gjort slik at treningssettet er s? stor som mulig, og er spesielt viktig med sm? datasett.
Her velger jeg litt vilk?rlig at 60% av datasettet skal v?re treningssettet
library(tidymodels)
## ── Attaching packages ────────────────────────────────────── tidymodels 0.1.2 ──
## ? broom 0.7.5 ? recipes 0.1.15
## ? dials 0.0.9 ? rsample 0.0.9
## ? infer 0.5.4 ? tune 0.1.3
## ? modeldata 0.1.0 ? workflows 0.2.3
## ? parsnip 0.1.5 ? yardstick 0.0.8
## ── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
## x scales::discard() masks purrr::discard()
## x workflows::extract_recipe() masks tune::extract_recipe()
## x dplyr::filter() masks stats::filter()
## x recipes::fixed() masks stringr::fixed()
## x dplyr::lag() masks stats::lag()
## x yardstick::spec() masks readr::spec()
## x recipes::step() masks stats::step()
set.seed(123) # for reproduserbarhet
train_test_split <- rsample::initial_split(iris_bin, prop=0.6, strata=Species)
train <- training(train_test_split)
test <- testing(train_test_split)
blomsterartfordelingen i treningsettet og testsettet
table(train$Species)
##
## versicolor virginica
## 30 30
table(test$Species)
##
## versicolor virginica
## 20 20
5 Bin?r klassifisering (med XGBoost)
En popul?r algoritme i Kaggle ML konkurranser de siste 5-7 ?rene er XGBoost, som st?r for (e)Xtreme gradient boosting
Dette er en salgs “ensemble” metode, som betyr at flere modeller aggregeres sammen til ? bygge den endelige modellen. XGBoost l?rer ved bruk av s?kalte beslutningstr?r
I motsetning til LDA har XGBoost flere hyperparametre. Kort sagt er disse konfigurasjoner som ikke kan direkte estimeres fra dataen.
For ? illustrere hvordan ? endre verdiene, gj?r vi det for “nrounds” (antall treningsrunder) med trees parameteren.
library(xgboost)
##
## Attaching package: 'xgboost'
## The following object is masked from 'package:dplyr':
##
## slice
xgb.fit <-
boost_tree(trees=50) %>%
set_engine("xgboost") %>%
set_mode("classification") %>%
fit(Species ~ ., data=train)
Det er flere m?ter ? bygge modeller med tidymodels
Med boost_tree() angir vi at vi vil bruke em trebasertalgoritme,
og med set_engine velger vi en spesifikk implementering.
Med set_mode angir vi hva slags modell vi bygger
xgboost og enkelte andre maskinl?ringsalgoritmer
kan bygges opp til ? skille mellom kategorier, eller ? predikere kontinuerlige
verdier (av og til kjent som “regresjonsmodeller”)
fit funksjonen “bygger” selve modellen p? datasettet vi mater inn,
og utvalgt “target” (Species)
xgb.pred = predict(xgb.fit, test)
xgb.pred
## # A tibble: 40 x 1
## .pred_class
## <fct>
## 1 versicolor
## 2 versicolor
## 3 versicolor
## 4 versicolor
## 5 versicolor
## 6 versicolor
## 7 versicolor
## 8 versicolor
## 9 versicolor
## 10 versicolor
## # … with 30 more rows
Det som egentlig skjer n?r man anvender en klassifiseringsalgoritme er at den predikerer
sannsynligheten for ? tilh?re de ulike blomsterartene. F?lgende kode lar oss inspisere
nettop dette.Her ser vi f.eks at
Legg merke til at summen av hver rad er 1.
Det man m? gj?re n? er ? spesifisere en “terskel” (cutoff) for ? f? til selve klassifiseringen.
Defaulten i de aller fleste bin?rklassifiseringssettinger er 0.5, eller 50%.
Dvs. hvis det er st?rre enn 50% sjanse for at en blomsterobservasjon er av “versicolr” art, s? blir den “tagget” som en
“veriscolor”. Under grensn blir den da tagget som “setosa”.
Siden summen av hver rad er 1, tar vi bare utgangspunkt i en av kolonnene.
Defaulten er f?rste kolonnen
Det kan v?re lurt ? velge andre terskler hvis datasettet er veldig ujevn, og kan dermed “vekte” klassifiseringen mot en kategori med f?rre observasjoner
xgb.pred.prob = predict(xgb.fit, test, type="prob")
xgb.pred.prob
## # A tibble: 40 x 2
## .pred_versicolor .pred_virginica
## <dbl> <dbl>
## 1 0.571 0.429
## 2 0.974 0.0262
## 3 0.974 0.0262
## 4 0.990 0.0101
## 5 0.974 0.0262
## 6 0.990 0.0101
## 7 0.974 0.0262
## 8 0.987 0.0128
## 9 0.978 0.0215
## 10 0.978 0.0215
## # … with 30 more rows
Hva gj?r vi med de predikerte og faktiske blomsterartene?
Vi lager enn krysstabell (“confusion matrix”) med conf_mat som teller opp de predikerte kategoriene opp mot de faktiske kategoriene.
Kolonnene angir den faktiske kategorien, men sradene angir den predikerte kategorien (dvs. Det maskinl?ringsmodellen tror).
Diagonalen viser antall observasjoner som har blitt korrekt klassifisert.
# Samle relevante data i en dataramme
xgb.res <- cbind(test$Species, xgb.pred, xgb.pred.prob)
# Gj?re om kolonnenavnene
names(xgb.res) <- c("obs", "pred", "pred_versicolor", "pred_virginica")
# Lage krysstabell
xgb.conf <- xgb.res %>%
conf_mat(obs, pred)
xgb.conf
## Truth
## Prediction versicolor virginica
## versicolor 19 3
## virginica 1 17
Men hvordan kan vi tolke hvorvidt vi har en bra maskinl?ringsmodell, og i tillegg, hvordan kan vi sammenligne maskinl?ringsmodeller, eller andre datasett?
6 Modellevalueringsm?l
6.1 N?yaktighet (accuracy)
N?yaktighet er antall korrekte prediksjoner delt p? totalt antall observasjone
sum(diag(xgb.conf$table))/sum(xgb.conf$table)
## [1] 0.9
Med to kategorier s? er i prinsippet det laveste vi kan f? basert p? sjanse 50%.
6.2 Presisjon
Presisjon er antall korrekte prediksoner i en kategori (f.eks. versicolor) delt p? total antall prediksjoner for samme kategori
Dvs. Hvor mye er riktig predikert for en spesifikk kategori
diag(xgb.conf$table) / rowSums(xgb.conf$table)
## versicolor virginica
## 0.8636364 0.9444444
6.3 Tilbakekalling (recall)
Det er antall korrekte prediksjoner for en spesifikk kategori delt p? total antall reele observasjoner i samme kategori
diag(xgb.conf$table) / colSums(xgb.conf$table) #tidymodels
## versicolor virginica
## 0.95 0.85
Kan vekte enten preisjon eller tilbakekalling for en av kategoriene
Balansegang mellom presisjon og tilbakekalling
Ogs? balansengang innad hver modellevalueringsm?l p? tvers av kategoriene
Tidymodels har innebygde funksjoner for disse evalueringsm?lene
xgb.res %>%
accuracy(obs, pred)
## # A tibble: 1 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.9
xgb.res %>%
recall(obs, pred)
## # A tibble: 1 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 recall binary 0.95
xgb.res %>%
precision(obs, pred)
## # A tibble: 1 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 precision binary 0.864
Legg merke til at vi her bare f?r ett presisjons- og tilbakekallingsestimat.
fordi den tar referanseniv?et. Du kan velge referansen med
event_level parameteren:
levels(xgb.res$obs)
## [1] "versicolor" "virginica"
xgb.res %>%
recall(obs, pred, event_level = "second")
## # A tibble: 1 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 recall binary 0.85
xgb.res %>%
precision(obs, pred, event_level = "second")
## # A tibble: 1 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 precision binary 0.944
xgb.res %>%
metrics(obs, pred)
## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.9
## 2 kap binary 0.8
Man kan ogs? f? ut andre evalueringsm?l hvis man i tillegg
legger inn predikert sannsynlighet for prediksjonene til
referanseniv?et (dvs. f?rste niv? i faktorvariabelen)
#metrics(lda.res, truth=obs, estimate=pred, pred_AV)
xgb.res %>%
metrics(obs, pred, pred_versicolor)
## # A tibble: 4 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.9
## 2 kap binary 0.8
## 3 mn_log_loss binary 0.340
## 4 roc_auc binary 0.968
Det er egentlig samme prosess for ? klassifisere “setosa” mot “virginca” eller mot “versicolor”.
Det er bare ? justere p? linjen hvor man dropper en av de tre tekstsjangerne
7 Feature importance
Feature importance er en slags form for bidragsestimat.
N?r man har mange egenskaper kan dette brukes til variabelutvelgelse:: fjerne egenskaper som ikke bidrar noe,
eller faktisk svekker modellen.
I det minste for man innsikt i hvilke egenskaper er viktige i modellene.
library(vip)
##
## Attaching package: 'vip'
## The following object is masked from 'package:utils':
##
## vi
vip(xgb.fit)
Vi kan hente ut disse “feature importance” sk?rene og plotte dem i
ggplot hvis vi vil justere p? figuren:
xgb_features <- vip::vi(xgb.fit)
xgb_features
## # A tibble: 4 x 2
## Variable Importance
## <chr> <dbl>
## 1 Petal.Length 0.905
## 2 Petal.Width 0.0749
## 3 Sepal.Width 0.0128
## 4 Sepal.Length 0.00734
ggplot(data=xgb_features, aes(x=reorder(Variable, Importance),
y=Importance))+
geom_col(aes(fill=Variable))+
coord_flip()+
theme(axis.text.y=element_blank(),
axis.ticks.y=element_blank())+
labs(x="Bidragsestimat",
y="",
title="Binary classification",
subtitle="XGboost",
caption="first try")
Det var i korte trekk en m?te ? utf?re maskinl?ring med tidymodels i R.
Men, det kan hende at vi utf?re noen preprosseringssteg
f?r vi bygger maskinl?ringsmodellene.
En typisk preprosesseringssteg er ? z-sk?re normalisere egenskapene,
dvs. at hver variabel f?r en gjennomsnitt p? 0 og en standardavvik p? 1.
For ? gj?re dette m? vi bruke recipe funksjonen.
iris_recipe <- recipe(Species ~ ., data = train) %>%
step_normalize(all_predictors())
Deretter m? vi initialisere en maskinl?ringsalgoritme
Her skal jeg fort vise dere hvordan dere velger en annen maskinl?ringsalgoritme
library("discrim")
##
## Attaching package: 'discrim'
## The following object is masked from 'package:dials':
##
## smoothness
lda_mod <-
discrim_linear() %>%
set_engine("MASS")
Og dytte v?r “oppskrift” og algoritmen i en s?kalt “workflow”
lda_wfl <- workflow() %>%
add_recipe(iris_recipe) %>%
add_model(lda_mod)
Ogs? kan vi enkelt mate inn “arbeidsflyten” i fit funksjonen
for ? bygge maskinl?ringsmodellen, og predikere testsettet med predict
lda_fit_norm <- fit(lda_wfl, data = train)
Det som egentlig skjer med workflow() metoden er f?lgende:
train_norm <- iris_recipe %>%
prep(training = train) %>%
bake(train)
test_norm <- iris_recipe %>%
prep(training = train) %>%
bake(test)
train_norm
## # A tibble: 60 x 5
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## <dbl> <dbl> <dbl> <dbl> <fct>
## 1 1.35 0.949 -0.200 -0.583 versicolor
## 2 0.342 0.949 -0.446 -0.359 versicolor
## 3 -1.17 -1.58 -1.06 -0.807 versicolor
## 4 -0.835 -0.173 -0.446 -0.807 versicolor
## 5 -2.01 -2.42 -1.67 -1.48 versicolor
## 6 -0.499 0.388 -0.814 -0.359 versicolor
## 7 -0.331 -1.86 -1.06 -1.48 versicolor
## 8 0.846 0.668 -0.569 -0.583 versicolor
## 9 -1.00 0.388 -0.446 -0.359 versicolor
## 10 -0.667 -0.453 -0.937 -1.48 versicolor
## # … with 50 more rows
test_norm
## # A tibble: 40 x 5
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## <dbl> <dbl> <dbl> <dbl> <fct>
## 1 1.18 0.668 0.0450 -0.359 versicolor
## 2 0.510 -0.173 -0.323 -0.359 versicolor
## 3 0.174 1.23 -0.200 -0.134 versicolor
## 4 -2.18 -1.29 -1.92 -1.48 versicolor
## 5 0.678 0.108 -0.323 -0.807 versicolor
## 6 -1.68 -0.453 -1.18 -0.583 versicolor
## 7 -0.162 0.108 -0.200 -0.583 versicolor
## 8 -1.00 0.108 -1.55 -0.807 versicolor
## 9 0.342 0.108 -0.691 -0.807 versicolor
## 10 0.678 0.388 -0.569 -0.583 versicolor
## # … with 30 more rows
Det fine med workflow() er at man slipper ? manuelt
passe p? at slike transformeringer p? testsettet er basert p? treningssettet
dvs. ikke basert p? hele datasettet (f?r oppdelingen), som er viktig for ? unng? “data leakage”
“Data leakage” kan f?re til modeller som virker bedre enn det de egentlig er.
Denne arbeidsflyten (med workflow() ) kan enkelt anvendes p? nye datasett, og lar deg i
prinsippet skrive mer konsist kode n?r det er flere preprosesseringssteg
Den kan i tillegg integreres med videre steg, som f.eks. kryss-validering, hyperparameteroptimalisering, m.fl
I tillegg er det mulig ? mate inn en liste med ulike kombinasjoner av
preprosesseringer, samt modeller, som er fint hvis man vil finne hvilke kombinasjoner
f?rer til best modellersultat.