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 ska vi använda tidymodels för att skapa ett komplett arbetsflöde som inbegriper val av prediktionsmodell, optimering av hyperparametrar med korsvalidering, samt selektion och utvärdering av bästa modellen. Vi börjar med att aktivera tre paket med följande kommandon:

R
# Aktivera tidymodels för prediktionsmodellering
library(tidymodels)

# Aktivera tidyverse för hantering av data
library(tidyverse)

# Aktivera vip för att extrahera variable importance
library(vip)

Data för denna övningen kommer från paketet survival, där data på individer med coloncancer finns att tillgå. Vi aktiverar survival och laddar alla dataset med cancerstudier. Därefter inspekterar vi colon, vilket är data vi skall arbeta med.

R
# Aktivera survival-paketet
library(survival)

# Aktivera all dataset med cancerstudier
data(cancer)

# Inspektera första 10 rader i studien om coloncancer
head(colon, 10)
Resultat
   id study      rx sex age obstruct perfor adhere nodes status differ extent surg node4 time etype
1   1     1 Lev+5FU   1  43        0      0      0     5      1      2      3    0     1 1521     2
2   1     1 Lev+5FU   1  43        0      0      0     5      1      2      3    0     1  968     1
3   2     1 Lev+5FU   1  63        0      0      0     1      0      2      3    0     0 3087     2
4   2     1 Lev+5FU   1  63        0      0      0     1      0      2      3    0     0 3087     1
5   3     1     Obs   0  71        0      0      1     7      1      2      2    0     1  963     2
6   3     1     Obs   0  71        0      0      1     7      1      2      2    0     1  542     1
7   4     1 Lev+5FU   0  66        1      0      0     6      1      2      3    1     1  293     2
8   4     1 Lev+5FU   0  66        1      0      0     6      1      2      3    1     1  245     1
9   5     1     Obs   1  69        0      0      0    22      1      2      3    1     1  659     2
10  5     1     Obs   1  69        0      0      0    22      1      2      3    1     1  523     1

I colon finns kolumnen status som är studiens utfallsmått. Värdet 1 indikerar att patienten dött och värdet 0 indikerar att patienten överlevt. Vi kommer modifiera colon data med följande steg nedan:

  1. Vi kodar om status 1 till Yes och 0 till No och gör variabeln till en faktor med funktionen recode_factor().
  2. Vi kodar om sex så det står Kvinna och Man istället för 0 respektive 1.
  3. Vi tar bort kolumnerna time, id och study eftersom vi inte skall använda de i vår modell.
  4. Variabler med 0/1 kodas om till No/Yes på ett finurligt sätt med mutate(across()). Notera att vi med across() kan välja flera kolumner som vi vill tillämpa funktionen recode_factor() på, och därefter representeras kolumnen av .x i funktionen. Se rad 5 nedan.
R
colon2 <- colon |> 
  mutate(status = recode_factor(status, "1"="Yes", "0"="No"),
         sex = recode_factor(sex, "0"="Kvinna", "1"="Man")) |> 
  select(-time, -id, -study) |> 
  mutate(across(c(obstruct, adhere, perfor, surg, node4), ~recode_factor(.x, "0"="No", "1"="Yes")))
  
# Inspektera data
head(colon2, 10)
Resultat
        rx    sex age obstruct perfor adhere nodes status differ extent surg node4 etype
1  Lev+5FU    Man  43       No     No     No     5    Yes      2      3   No   Yes     2
2  Lev+5FU    Man  43       No     No     No     5    Yes      2      3   No   Yes     1
3  Lev+5FU    Man  63       No     No     No     1     No      2      3   No    No     2
4  Lev+5FU    Man  63       No     No     No     1     No      2      3   No    No     1
5      Obs Kvinna  71       No     No    Yes     7    Yes      2      2   No   Yes     2
6      Obs Kvinna  71       No     No    Yes     7    Yes      2      2   No   Yes     1
7  Lev+5FU Kvinna  66      Yes     No     No     6    Yes      2      3  Yes   Yes     2
8  Lev+5FU Kvinna  66      Yes     No     No     6    Yes      2      3  Yes   Yes     1
9      Obs    Man  69       No     No     No    22    Yes      2      3  Yes   Yes     2
10     Obs    Man  69       No     No     No    22    Yes      2      3  Yes   Yes     1

Notera att vi hädanefter använder colon2.

Skapa träningsdata, testdata och resampling (korsvalidering)

Vi delar nu colon2 i träningsdata (70% av data) och testdata (30%). Eftersom studien syftar till att undersöka händelsen död väljer vi att stratifiera delningen, så att andelen som dör är lika stor i träningsdata och testdata. Detta är önskvärt för att modellen skall testas och tränas på realistiska förhållanden. Om många observationer finns blir det med stor sannolikhet lika stor andel som dör i träningsdata och testdata (eftersom träningsdata och testdata skapas helt slumpmässigt). Denna funktionen är därför särskilt användbar om antalet observationer är få, varvid funktionen säkerställer att det slumpmässiga urvalet har lika stora andelar dödsfall i träningsdata och testdata.

R
# Skapar träningsdata och testdata med set.seed() först
set.seed(1)
colon_split <- initial_split(colon2, strata=status, prop=0.7)
colon_train <- training(colon_split)
colon_test  <- testing(colon_split)

# Skapar korsvalidering av träningsdata
my_cross_validation <- vfold_cv(colon_train, 
                                v = 5,
                                repeats=10,
                                strata = status)

Notera även att vi angett set.seed(1), vilket på svenska kallas att så ett frö. När man sår ett frö i dessa sammanhang så tar man bort slumpen i den funktionen som körs efter set.seed(). I detta fall innebär det endast att vi kan erhålla precis samma observationer i träningsdata och testdata varje gång vi gör om delningen. Det är önskvärt för att få reproducerbara resultat. Du kan använda vilken siffra som helst i set.seed().

Resampling betyder att samma datamängd återanvänds för att träna och utvärdera en modell. I detta fall används en typ av resampling som kallas korsvalidering (se Korsvalidering). Korsvalideringen görs i träningsdata och med hjälp av funktionen vfold_cv(), som slumpmässigt delar in träningsdata in i 5 lika stora delar (folds) och varje modell kommer tränas i alla folds förutom en, vilken används för att testa modellen. Detta upprepas 10 gånger (repeats=10) för varje modell, vilket medför att alla modeller tränas och utvärderas i slumpmässigt valda delar av träningsdata. Den modellen som presterar bäst under korsvalideringen är med största sannolikhet även bäst på testdata (och därmed den bästa modellen).

Specificera en modell och metod

I nästa steg specificerar vi en modell, nämligen extreme gradient boosting (XGBoost). Vi gör detta genom att använda funktionen boost_tree(), där samtliga hyperparametrar ställs till tune(), vilket innebär att de ska optimeras via resampling. Vi specificerar att metoden med set_engine() och typ av prediktion med set_mode(). Vi skriver därefter ut modellobjektet (xgb_spec).

R
xgb_spec <-
  boost_tree(mtry = tune(),
             tree_depth = tune(),
             trees = tune(),
             learn_rate = tune(),
             min_n = tune(),
             loss_reduction = tune(),
             sample_size = tune(),
             stop_iter = tune()) |> 
  set_engine('xgboost') |> 
  set_mode('classification')

xgb_spec
Resultat
Boosted Tree Model Specification (classification)

Main Arguments:
  trees = tune()
  min_n = tune()
  tree_depth = tune()
  learn_rate = tune()
  loss_reduction = tune()
  sample_size = tune()
  stop_iter = tune()

Computational engine: xgboost 

Hyperparametrarna i xgboost är som följer

  • mtry: The number of predictors that will be randomly sampled at each split when creating the tree models.
  • trees: The number of trees contained in the ensemble. Default: 15.
  • min_n: The minimum number of data points in a node that are required for the node to be split further. Default: 1.
  • tree_depth: The maximum depth of the tree (i.e. number of splits). Default: 6.
  • learn_rate: The rate at which the boosting algorithm adapts from iteration-to-iteration. Default: 0.3.
  • loss_reduction: The reduction in the loss function required to split further. Type: double, default: 0.0.
  • sample_size: The amount of data exposed to the fitting routine. Type: double, default: 1.0.
  • stop_iter: The number of iterations without improvement before stopping. Type: integer, default: Inf.

Grid search

I nästa steg ska vi definiera en matris med kombinationer av hyperparametrar. Dessa matriser kallas för grid. Korsvalideringens syfte är att söka igenom matrisen efter bästa modellen, vilket kallas grid search. Vi kan skapa en grid med flera olika funktioner, varav en är grid_latin_hyperube(). Med denna funktionen skapas kombinationer som, baserat på antal kombinationer vi tillåter, täcker så "stort område som möjligt".

R
xgb_grid <- grid_latin_hypercube(
  trees(),
  min_n(),
  learn_rate(), 
  loss_reduction(),
  tree_depth(),
  sample_prop(),
  stop_iter(),
  finalize(mtry(), colon_train),
  size = 20)

Vår grid har totalt 50 olika kombinationer av hyperparametrar (size=50), vilket innebär att vi testar 50 olika modeller. Notera att parametern mtry (antal variabler) specificeras i funktionen finalize(), vilket förklaras av att mtry inte kan anta vilket värde som helst (antal kolumner kan vara för få).

R
xgb_grid
Resultat
   trees min_n learn_rate loss_reduction tree_depth sample_size stop_iter  mtry
   <int> <int>      <dbl>          <dbl>      <int>       <dbl>     <int> <int>
 1    55    24   4.15e- 2       2.28e- 2          6       0.358        14    10
 2   413    32   1.35e- 5       5.05e- 6         14       0.887        17     7
 3  1843    33   7.88e- 3       2.08e- 8          2       0.491        15     7
 4  1214    27   9.13e- 8       1.57e- 3          4       0.227        10     5
 5  1495    26   3.26e-10       5.35e- 7          4       0.596        15     6
 6   476    20   5.08e-10       1.33e- 8         13       0.828        11     2
 7  1117     8   5.41e- 3       2.92e- 7          8       0.731        17    12
 8   351     4   3.27e- 4       2.31e- 9         10       0.864        12    12
 9   310    35   9.31e- 9       1.98e+ 1         13       0.809        12     4
10  1450    40   1.34e- 6       2.03e-10          9       0.245         4    13
11  1269    35   1.13e- 8       4.29e- 1          4       0.615         4     7
12   990    22   8.24e- 2       2.13e- 5         11       0.532         8     2
13   128    23   5.87e-10       3.18e- 9         10       0.786        18     9
14  1347    13   5.02e- 9       8.81e- 2         12       0.722         7     9
15    94    38   6.35e- 7       9.03e- 5          7       0.168         5     4
16  1745    21   9.00e- 3       3.48e- 4          3       0.385        14     5
17  1791     6   1.93e- 4       1.23e- 6         10       0.413         3    12
18  1399    15   2.79e- 7       3.69e- 5          5       0.675         9    10
19   672    12   6.34e- 5       9.73e- 1          6       0.777        14     1
20  1531    12   2.25e-10       4.90e+ 0          3       0.266        15    11
21   893    19   1.54e- 6       1.15e+ 1          9       0.539         6    10
22   264    36   5.98e- 8       3.36e+ 0         14       0.208        10     8
23   233    38   1.25e- 9       2.41e- 3          9       0.998         6     3
24  1657    10   1.40e- 3       7.56e- 8          3       0.315        19     5
25   928     6   2.55e- 9       7.20e- 1          6       0.338        13     8
26  1418    33   2.85e- 4       1.40e- 7          1       0.970        11     6
27  1059     7   1.71e- 8       2.38e- 6          7       0.709         8     1
28   835    30   3.09e- 6       5.97e-10          8       0.131         3     8
29   162     9   2.75e- 8       1.66e-10         13       0.759        16     6
30   519     2   9.92e- 5       4.33e- 3          6       0.427        11    11
31  1153    23   3.08e- 5       1.07e- 4         13       0.922         7     2
32  1966     5   7.63e- 6       1.85e+ 0          2       0.398         8    10
33  1171    18   3.14e- 9       3.13e-10         12       0.652         9     1
34  1312    28   9.25e-10       2.25e- 4         12       0.576        13     3
35   590    20   5.51e- 5       1.34e- 9         15       0.189        18     6
36   699    25   3.85e- 8       1.04e- 2          8       0.289        13    11
37   389    14   5.52e- 2       7.44e- 6          7       0.474        16    11
38  1632     3   6.07e- 6       1.62e- 6          2       0.893        10    12
39   614    27   8.39e- 4       3.97e- 2         11       0.115        12    13
40  1011    16   4.35e- 7       1.07e+ 1         11       0.505        19     6
41   530    16   4.86e- 4       6.27e- 9          5       0.211         7     3
42   869    39   2.56e- 3       2.60e- 7         14       0.557        19     8
43  1821    16   2.37e- 7       2.43e- 1          9       0.678        20     4
44   799    29   4.12e- 6       1.63e- 2          5       0.935         6     3
45  1916    31   2.44e- 5       9.56e- 2          4       0.136        20     4
46    32    11   2.04e- 2       4.02e- 8         12       0.964         5     7
47  1713     8   1.74e- 7       6.99e- 9         15       0.838         9     3
48  1945    36   1.66e- 2       1.35e- 5          3       0.640         5     9
49  1567    18   1.50e-10       4.94e- 4          2       0.457        17     8
50   740    29   1.93e- 3       1.16e- 3          8       0.318        18    11

Skapa ett recept

I nästa steg skapar vi ett recept där vi definierar vad som är studiens utfall samt vilka prediktorer vi vill använda. Vi definierar också vilka pre-processing steg vi vill göra. Se kommentarerna med förklaringar.

R
xgb_recipe <- recipe(status ~ ., data=colon_train) |> 
  # step_normalize normaliserar  alla kontinuerliga prediktorer
  step_normalize(all_numeric()) |> 
  # step_dummy skapar dummy-variabler av alla kategoriska variabler. Detta krävs av xgboost
  step_dummy(all_nominal_predictors()) |> 
  # step_novel är en metod som behövs för att hantera situationen som kan uppstå om
  # testdata innehåller nivåer på kategoriska variabler som inte fanns i träningsdata
  step_novel(all_nominal_predictors())

Skapa ett arbetsflöde

Därefter skapar vi ett arbetsflöde med receptet och modellen:

R
xgb_wf <- workflow() |> 
  add_recipe(xgb_recipe) |> 
  add_model(xgb_spec)

xgb_wf
Resultat
══ Workflow ═══════════════════════════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: boost_tree()

── Preprocessor ───────────────────────────────────────────────────────────────────────────────────────
3 Recipe Steps

step_normalize()
step_dummy()
step_novel()

── Model ──────────────────────────────────────────────────────────────────────────────────────────────
Boosted Tree Model Specification (classification)

Main Arguments:
  mtry = tune()
  trees = tune()
  min_n = tune()
  tree_depth = tune()
  learn_rate = tune()
  loss_reduction = tune()
  sample_size = tune()
  stop_iter = tune()

Computational engine: xgboost 

Träna modellerna

Därefter genomförs träningen (dvs modellerna skapas). Observera argumentet metrics = yardstick::metric_set(roc_auc), som innebär att vi från paketet yardstick hämtar funktionen metric_set och sätter den till roc_auc (area under ROC-kurvan). Det blir således arean under ROC-kurvan som används för att jämföra modellerna.

R
xgb_res <- tune_grid(
  xgb_wf,
  resamples = my_cross_validation,
  grid      = xgb_grid,
  metrics   = yardstick::metric_set(roc_auc),
  control   = control_grid(save_pred = TRUE,
                           parallel_over = "everything",
                           save_workflow = TRUE))

Parallellisera träningen

Om du vill påskynda ovanstående steg kan du göra det genom att parallellisera körningen på flera av datorns kärnor. Detta görs med paketen parallell och doParallell. För att mäta tiden det tar att göra körningen använder vi paketet tictoc. Koden för en parallelliserad beräkning på 6 kärnor är som följer:

R

# Aktivera paketen (glöm inte att installera paketen)
library(tictoc)
library(parallel)
library(doParallel)

# Skapa ett kluster i R med 6 kärnor, så att körningen görs på 6 kärnor
mitt_kluster <- makePSOCKcluster(6)

# Aktivera ditt kluster
registerDoParallel(mitt_kluster) 

# Starta klockan
tic()

# Gör körningen
xgb_res <- tune_grid(
  xgb_wf,
  resamples = my_cross_validation,
  grid      = xgb_grid,
  metrics   = yardstick::metric_set(roc_auc),
  control   = control_grid(save_pred = TRUE,
                           parallel_over = "everything",
                           save_workflow = TRUE))

# Stoppa klustret
stopCluster(mitt_kluster)

# Stoppa klockan
toc()

Utvärdera modellerna

Resultaten inspekteras sedan med funktionen collect_metrics(). Den skapar den tabell med information om varje modells hyperparametrar samt modellens resultat under korsvalideringen. Kolumnen mean visar genomsnittlig ROC AUC under korsvalideringen. Vi sorterar alla modellerna på fallande ROC AUC (bästa modellen hamnar högst upp):

R
collect_metrics(xgb_res) %>% arrange(desc(mean))
Resultat
    mtry trees min_n tree_depth  learn_rate loss_reduct…¹ sampl…² stop_…³ .metric .esti…⁴  mean     n std_err 
 1    12   985     8          6 0.0719           9.45e- 7   0.695       4 roc_auc binary  0.743    50 0.00337 
 2     6  1713     2          2 0.0432           2.83e-10   0.517       7 roc_auc binary  0.733    50 0.00323 
 3     2   794     5          3 0.0567           4.73e- 7   0.904      14 roc_auc binary  0.732    50 0.00334 
 4     6   716     5          9 0.000956         1.62e- 2   0.880       9 roc_auc binary  0.727    50 0.00399 
 5     9   159     6          5 0.00270          2.68e- 6   0.796      16 roc_auc binary  0.712    50 0.00411 
 6     4  1394     9         13 0.000000417      1.68e- 9   0.862      13 roc_auc binary  0.707    50 0.00425 
 7     1   447     4         15 0.0000679        2.17e+ 0   0.986       7 roc_auc binary  0.705    50 0.00383 
 8     6   857    12          7 0.00000888       1.19e-10   0.745      12 roc_auc binary  0.704    50 0.00430
 9     9  1207    12          4 0.000126         3.87e- 2   0.946       6 roc_auc binary  0.698    50 0.00415
10    10  1280    14         12 0.00000152       1.17e- 8   0.666       6 roc_auc binary  0.695    50 0.00452 
# … with 40 more rows, and abbreviated variable names ¹​loss_reduction, ²​sample_size, ³​stop_iter, ⁴​.estimator
# ℹ Use `print(n = ...)` to see more rows

För att granska hur olika hyperparametrar påverkade ROC-AUC skriver vi följande, tämligen invecklade, kod:

R
xgb_res_plot <-  
  xgb_res |> 
  collect_metrics() |> 
  filter(.metric == "roc_auc") |> 
  select(mean, mtry:sample_size) |> 
  pivot_longer(mtry:sample_size,
               values_to = "value",
               names_to = "parameter") |>
  ggplot(aes(value, mean, color = parameter)) +
  geom_point(alpha = 0.8, show.legend = FALSE) +
  facet_wrap(~parameter, scales = "free_x") +
  labs(x = NULL, y = "AUC")

xgb_res_plot

För att se de 5 bästa modellerna i fallande ordning skriver vi show_best():

R
show_best(xgb_res, "roc_auc")
Resultat
   mtry trees min_n tree_depth learn_rate loss_reduction sample_size stop_…¹ .metric .esti…²  mean     n 
1    12   985     8          6   0.0719         9.45e- 7       0.695       4 roc_auc binary  0.743    50 
2     6  1713     2          2   0.0432         2.83e-10       0.517       7 roc_auc binary  0.733    50 
3     2   794     5          3   0.0567         4.73e- 7       0.904      14 roc_auc binary  0.732    50 
4     6   716     5          9   0.000956       1.62e- 2       0.880       9 roc_auc binary  0.727    50 
5     9   159     6          5   0.00270        2.68e- 6       0.796      16 roc_auc binary  0.712    50 
# … with abbreviated variable names ¹​stop_iter, ²​.estimator

För att extrahera den bästa modellen skriver vi select_best():

R
best_model <- select_best(xgb_res, "roc_auc")
best_model
Resultat
   mtry trees min_n tree_depth learn_rate loss_reduction sample_size stop_iter .config              
1    12   985     8          6     0.0719    0.000000945       0.695         4 Preprocessor1_Model09

Slutligen kör vi om den bästa modellen på hela träningsdatasettet (utan korsvalidering), och utvärderar modellen på testdata. Detta kan göras med funktionen last_fit(), där vi anger colon_split, som innehåller både träningsdata och testdata.

R
final_model <- 
  xgb_res %>% 
  extract_workflow("xgb_wf") %>%
  finalize_workflow(best_model) %>%
  last_fit(split = colon_split,
           metrics = yardstick::metric_set(roc_auc))

Denna modellens förmåga att predicera på testdata kan nu utvärderas:

R
collect_metrics(final_model)
Resultat
  .metric .estimator .estimate .config             
  roc_auc binary         0.707 Preprocessor1_Model1

ROC AUC var 0.707 för bästa modellen.

Vi kan extrahera prediktionerna på testdata med följande kod:

R
collect_predictions(final_model)
Resultat
   id               .pred_Yes .pred_No  .row status .config             
 1 train/test split     0.720   0.280      2 Yes    Preprocessor1_Model1
 2 train/test split     0.120   0.880      4 No     Preprocessor1_Model1
 3 train/test split     0.475   0.525      8 Yes    Preprocessor1_Model1
 4 train/test split     0.940   0.0603    13 Yes    Preprocessor1_Model1
 5 train/test split     0.905   0.0946    14 Yes    Preprocessor1_Model1
 6 train/test split     0.503   0.497     20 No     Preprocessor1_Model1
 7 train/test split     0.436   0.564     21 No     Preprocessor1_Model1
 8 train/test split     0.612   0.388     22 No     Preprocessor1_Model1
 9 train/test split     0.562   0.438     30 No     Preprocessor1_Model1
10 train/test split     0.831   0.169     32 Yes    Preprocessor1_Model1
# … with 548 more rows

En ROC-kurva kan enkelt skapas med följande kod:

R
collect_predictions(final_model) %>%
  roc_curve(status, .pred_Yes) %>%
  ggplot(aes(x = 1 - specificity, y = sensitivity)) +
  geom_path() +
  geom_abline(lty = 3) +
  coord_equal() +
  theme_bw()

Variable importance

Vi kan extrahera variablernas betydelse (relative variable importance) med paketet vip, enligt följande:

R
library(vip)

final_model %>%
  extract_fit_parsnip() %>%
  vip(geom = "point")

Referenser

Hvitfeldt E, Frick H (2022). censored: 'parsnip' Engines for Survival Models. R package version 0.1.1, https://github.com/tidymodels/censored.

Tidymodels

Kuhn M, Vaughan D (2022). parsnip: A Common API to Modeling and Analysis Functions. https://github.com/tidymodels/parsnip, https://parsnip.tidymodels.org/.

Vaughan D, Couch S (2022). workflows: Modeling Workflows. https://github.com/tidymodels/workflows, https://workflows.tidymodels.org.

XGBoost: A Scalable Tree Boosting System

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.