I detta kapitel används data Pimaindianer. Diabetes är vanligare bland Pimaindianer och i dessa data har 768 individer undersökts avseende diabetes, blodsocker, etc. Data finns i paketet mlbench. Du behöver installera mlbench först och därefter aktivera data. Vi kommer dessutom skapa en kopia av data som objektet pima.
# Installera och aktivera mlbench
install.packages("mlbench")
library(mlbench)
# Aktivera data
data(PimaIndiansDiabetes)
# Kopiera data
pima <- PimaIndiansDiabetes
# Inspektera data
head(pima) pregnant glucose pressure triceps insulin mass pedigree age diabetes
1 6 148 72 35 0 33.6 0.627 50 pos
2 1 85 66 29 0 26.6 0.351 31 neg
3 8 183 64 0 0 23.3 0.672 32 pos
4 1 89 66 23 94 28.1 0.167 21 neg
5 0 137 40 35 168 43.1 2.288 33 pos
6 5 116 74 0 0 25.6 0.201 30 negDen sista kolumnen (diabetes) anger om personen har diabetes (pos) eller inte (neg). I nästa steg kodar vi om kolumnen diabetes till Yes eller Nej istället för pos respektive neg.
pima$diabetes = recode_factor(pima$diabetes,
'pos'='Yes',
'neg'='No')
head(pima) pregnant glucose pressure triceps insulin mass pedigree age diabetes
1 6 148 72 35 0 33.6 0.627 50 Yes
2 1 85 66 29 0 26.6 0.351 31 No
3 8 183 64 0 0 23.3 0.672 32 Yes
4 1 89 66 23 94 28.1 0.167 21 No
5 0 137 40 35 168 43.1 2.288 33 Yes
6 5 116 74 0 0 25.6 0.201 30 NoSkapa träningsdata och testdata
Vi använder sedan funktionen initial_split() för att skapa träningsdata och testdata. Vi anger prop=0.7 för att 70% av data skall användas till träning. Med funktionen training() och testing() extraherar vi sedan träningsdata och testdata till separata filer.
# Skapar träning och test
pima_split <- initial_split(pima, strata=diabetes, prop=0.7)
pima_train <- training(pima_split)
pima_test <- testing(pima_split)Notera i koden ovan att vi angav strata=diabetes, vilket innebär att vår split stratifieras på kolumnen diabetes, vilket leder till att andelen med Yes och No blir lika stor i träningsdata och testdata. Detta är önskvärt för att andelen som drabbas av diabetes skall vara lika stor i träningsdata och testdata.
Specificera modell och metod
Vi börjar med att specificera vilken typ av modell vi vill använda. Detta görs med hjälp av funktioner i parsnip-paketet. Vi börjar med att specificera att vi vill använda logistisk regression, vilket vi gör med funktionen logistic_reg() i parsnip. Om vi endast skriver ut funktionens namn och kör koden erhåller vi följande resultat:
logistic_reg()Logistic Regression Model Specification (classification)
Computational engine: glm Av ovanstånde kan utläsas att en logistisk modell har specificerats för att göra klassificering (classification), vilket innebär att utfallet är en kategorisk variabel. Detta kommer göras med glm (generalized linear model), vilket är den klassiska formen av logistisk regression. Metoden ("engine") GLM är alltså standardmetoden för logistisk regression i Tidymodels. Vi kan dock specificera andra metoder för logistisk regression. Vi gör detta genom att addera funktionen set_engine() och specificera en metod. Alla alternativen för logistic_reg() finns i dokumentationen här.
Vi har ännu inte sparat vår specifikation. Nedan specificerar vi en logistisk regression med hjälp av GLM och den skall användas för klassifikation. Det sistnämnda specificeras med funktionen set_mode(). Denna funktion är särskilt viktig om modellen kan användas både för klassificering (predicera kategoriska utfall) och regression (predicera kontinuerliga utfall). Vi sparar specifikationen i objektet log_mod.
log_mod <-
logistic_reg() |>
set_engine("glm") |>
set_mode("classification")
log_modLogistic Regression Model Specification (classification)
Computational engine: glm Specificera en formel
Vi går vidare med att specificera själva formeln med funktionen fit() och adderar detta till objektet log_mod. Notera att formeln diabetes ~ . betyder att vi skapar en modell där diabetes är utfallet och alla andra variabler är prediktorer, vilket specificeras med punkten (.):
log_fit <-
log_mod |>
fit(diabetes ~ ., data=pima_train)
# Skriv ut resultatet
log_fitparsnip model object
Call: stats::glm(formula = diabetes ~ ., family = stats::binomial,
data = data)
Coefficients:
(Intercept) pregnant glucose pressure triceps insulin mass pedigree age
8.591e+00 -1.499e-01 -3.117e-02 1.470e-02 5.995e-06 1.321e-03 -1.056e-01 -7.798e-01 -2.357e-02
Degrees of Freedom: 536 Total (i.e. Null); 528 Residual
Null Deviance: 694.2
Residual Deviance: 500.4 AIC: 518.4Extrahera resultat
Resultatet blir prydligare med funktionen tidy() från paketet broom. Vi anger att argumenten exponentiate och conf.int är TRUE för att erhålla exponentierade koefficienter samt konfidensintervaller:
tidy(log_fit, exponentiate=T, conf.int=T) term estimate std.error statistic p.value conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 5382. 0.873 9.84 7.85e-23 1040. 32053.
2 pregnant 0.861 0.0386 -3.88 1.05e- 4 0.797 0.928
3 glucose 0.969 0.00431 -7.24 4.61e-13 0.961 0.977
4 pressure 1.01 0.00648 2.27 2.32e- 2 1.00 1.03
5 triceps 1.00 0.00843 0.000711 9.99e- 1 0.984 1.02
6 insulin 1.00 0.00102 1.29 1.97e- 1 0.999 1.00
7 mass 0.900 0.0192 -5.51 3.52e- 8 0.866 0.933
8 pedigree 0.458 0.365 -2.14 3.27e- 2 0.222 0.930
9 age 0.977 0.0109 -2.16 3.04e- 2 0.956 0.998Kolumnen estimate är odds ratio (oddskvot). Kolumnerna conf.low och conf.high är nedre respektive övre gräns för 95% konfidensintervall.
Predicera
I nästa steg används modellen för att predicera på testdata (pima_test). Detta görs med funktionen predict(). Resultatet sparas i en ny kolumn i pima_test. Den nya kolumnen kallas predicted:
pima_test$predicted <- predict(log_fit,
new_data = pima_test,
type="prob")Om vi granskar pima_test ser vi att den nya kolumnen faktiskt består av två kolumner, nämligen .pred_Yes och .pred_No. Det innebär att detta är en "nestad" kolumn eftersom den består av flera kolumner.
Nestade kolumner är förvisso effektiva lagringssätt men i detta fall något opraktiskt. Därför bryter vi ut kolumnerna till separata kolumner med hjälp av funktionen unnest(), enligt följande:
pima_test <- unnest(pima_test, predicted)
# Inspektera data
head(pima_test) pregnant glucose pressure triceps insulin mass pedigree age diabetes .pred_Yes .pred_No
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <dbl> <dbl>
1 6 148 72 35 0 33.6 0.627 50 Yes 0.719 0.281
2 0 137 40 35 168 43.1 2.29 33 Yes 0.894 0.106
3 5 116 74 0 0 25.6 0.201 30 No 0.137 0.863
4 3 78 50 32 88 31 0.248 26 Yes 0.0620 0.938
5 10 115 0 0 0 35.3 0.134 29 No 0.597 0.403
6 10 168 74 0 0 38 0.537 34 Yes 0.885 0.115
7 1 189 60 23 846 30.1 0.398 59 Yes 0.760 0.240
8 5 166 72 19 175 25.8 0.587 51 Yes 0.658 0.342
9 7 100 0 0 0 30 0.484 32 Yes 0.361 0.639
10 0 118 84 47 230 45.8 0.551 31 Yes 0.384 0.616Kolumnerna .pred_Yes och .pred_No är de predicerade sannolikheterna för utfallen Yes respektive No enligt den logistiska regressionsmodellen. Summan av sannolikheten för Yes och No är alltid 1.0.
Utvärdera modellen
Om vi endast betraktar sannolikheten för Yes (.pred_Yes) så kan vi förslagsvis välja att predicera att individen har diabetes om sannolikheten är >0.5 (större än 0.5). Gränsen 0.5 används ofta som standard i många paket och funktioner eftersom värden >0.5 innebär att händelsens inträffande (i detta fall diabetes) är mer sannolik än att händelsen uteblir. Denna gränsen kallas ofta decision threshold. Man kan självklart sänka decision threshold till 0.4, 0.35, 0.2, eller ännu lägre. Ju lägre tröskeln är desto fler individer med diabetes kommer att upptäckas men desto fler individer får falskpositiva diagnoser. Det finns enkla funktioner för att utvärdera sensitivitet och specifictet för ett stort antal decision thresholds. Med funktionen roc_curve() beräknar vi sensitivitet och specifictet för över 200 decision thresholds:
pima_test %>%
roc_curve(diabetes, .pred_Yes) .threshold specificity sensitivity
<dbl> <dbl> <dbl>
1 -Inf 0 1
2 0.00171 0 1
3 0.0133 0.00667 1
4 0.0176 0.0133 1
5 0.0190 0.0200 1
6 0.0210 0.0267 1
7 0.0260 0.0333 1
8 0.0279 0.0400 1
9 0.0282 0.0467 1
10 0.0297 0.0533 1
....
Vi kan automatiskt beräkna andrean under ROC-kurvan med följande kod:
pima_test %>%
roc_auc(diabetes, .pred_Yes).metric .estimator .estimate
roc_auc binary 0.840AUC ROC är 0.840.
Kurvan för ROC erhålles med följande kod:
pima_test %>%
roc_curve(diabetes, .pred_Yes) %>%
ggplot(aes(x = 1 - specificity, y = sensitivity)) +
geom_path() +
geom_abline(lty = 3) +
coord_equal()
Vi kan för intresses skull göra om modellen fast denna gången med extreme gradient boosting (XGBoost), som är en mycket potent modell. XGBoost finns i funktionen boost_tree(). Den har hyperparametrar (exempelvis tree_depth) som egentligen behöver optimeras men det görs inte i detta exemplet. Hyperparametrarna får ett enda standardvärde för enkelhetens skull. På rad 13 skapar vi en ny data frame kallad resultat, vilken endast innehåller kolumnen diabetes från pima_test. På rad 16 använder vi får XGBoost modell för att predicera sannolikheten för diabetes och resultatet skickas till kolumnen predicted2 i resultat. På Rad 21 unnestar vi kolumnen predicted2. Därefter beräknas AUC-ROC och ROC-kurvan.
# Specificera modell och metod
xgb_mod <-
boost_tree(tree_depth = 4,
trees = 600,
learn_rate = 0.005,
min_n = 5,
sample_size = 0.5) |>
set_engine('xgboost') |>
set_mode('classification') |>
fit(diabetes ~ ., data=pima_train)
# Skapa en data frame med faktiska diagnoser
resultat <- data.frame(diabetes = pima_test$diabetes)
# Predicera diagnos
resultat$predicted2 <- predict(xgb_mod,
new_data = pima_test,
type="prob")
# Unnesta kolumnen predicted2
resultat <- unnest(resultat, predicted2)
# ROC AUC
resultat %>%
roc_auc(diabetes, .pred_Yes)
# ROC-kurvan
resultat |>
roc_curve(diabetes, .pred_Yes) |>
ggplot(aes(x = 1 - specificity, y = sensitivity)) +
geom_path() +
geom_abline(lty = 3) +
coord_equal().metric .estimator .estimate
roc_auc binary 0.844
ROC för XGBoost utan hyperparametertuning var 0.844 jämfört med 0.840 för logistisk regression.
Variable importance
Machine learning av denna typen har inbyggda metoder för att jämföra prediktorernas relativa betydelse. Detta kallas relative variable importance. Paketet vip (variable importance) kan användas för att extrahera denna informationen från Tidymodels objekt, enligt följande:
library(vip)
xgb_mod %>% vip(geom = "point")
Frågor och svar
Var finns listan över tillgängliga modeller och metoder? Klicka här för att komma till parsnips dokumentation.