Gå till index

Analys och forskning med R

0% färdig
0/0 Steps
  1. Analys och forskning med R och Posit (Rstudio)
  2. Grunderna i R och Rstudio
    7 Ämnen
  3. Importera, exportera, spara och ladda data
    5 Ämnen
  4. Strängar och regular expressions (regex)
    1 Ämne
  5. Bearbetning av data med dplyr
    12 Ämnen
  6. Visualisera och presentera
    14 Ämnen
  7. Explorerande och deskriptiva analyser
    6 Ämnen
  8. Prediktionsmodeller
    12 Ämnen
  9. Klassisk regressionsanalys
    8 Ämnen
  10. Machine learning (ML) och Artificiell Intelligens (AI)
    9 Ämnen
  11. Skapa prediktionsmodeller med Tidymodels
    6 Ämnen
  12. Hypotestester och epidemiologiska mått
    5 Ämnen
Avsnitt Progress
0% färdig

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.

R
# Installera och aktivera mlbench
install.packages("mlbench")
library(mlbench)

# Aktivera data
data(PimaIndiansDiabetes)

# Kopiera data
pima <- PimaIndiansDiabetes

# Inspektera data
head(pima)
Resultat
  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      neg

Den 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.

R
pima$diabetes = recode_factor(pima$diabetes,
                              'pos'='Yes',
                              'neg'='No')
                              
head(pima)
Resultat
  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       No

Skapa 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.

R
# 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:

R
logistic_reg()
Resultat
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.

R
log_mod <-
  logistic_reg() |> 
  set_engine("glm") |> 
  set_mode("classification")

log_mod
Resultat
Logistic 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 (.):

R
log_fit <-
  log_mod |> 
  fit(diabetes ~ ., data=pima_train)

# Skriv ut resultatet
log_fit
Resultat
parsnip 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.4

Extrahera 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:

R
tidy(log_fit, exponentiate=T, conf.int=T)
R
  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.998

Kolumnen 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:

R
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:

R
pima_test <- unnest(pima_test, predicted)

# Inspektera data
head(pima_test)
Resultat
   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.616

Kolumnerna .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:

R
pima_test %>%
  roc_curve(diabetes, .pred_Yes)
R
   .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:

R
pima_test %>%
  roc_auc(diabetes, .pred_Yes)
Resultat
.metric    .estimator     .estimate
roc_auc    binary         0.840

AUC ROC är 0.840.

Kurvan för ROC erhålles med följande kod:

R
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.

R
# 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()
R
.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:

R
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.

Logga in för att läsa

Som medlem har du tillgång till allt nuvarande och kommande utbildningsmaterial. Du kan följa din utveckling och få intyg på avklarade moment.