Lilla Forskarskolan: Forskningsmetoder och Analys med R Skapa prediktionsmodeller med Tidymodels Jämför flera recept och flera modeller med olika grid search-metoder

Avsnitt 11, Ämne 6

Startad

← Previous Next→

Jämför flera recept och flera modeller med olika grid search-metoder

Avsnitt Progress

0% färdig

I detta kapitel används tidymodels för att skapa och jämföra olika modeller i ett arbetsflöde (workflow). Liksom i föregående kapitel kommer vi optimera modellerna med korsvalidering. Följande kommer göras:

Alla modeller utvärderas genom resampling (korsvalidering).
Två olika recept utvärderas.
Tre olika modeller utvärderas.
- Två av modellerna genomgår optimering av hyperparametrar. För en av modellerna görs det genom manuell definiering av en grid (kombination av hyperparametrar). För den andra modellen görs det genom att låta R skapa en grid.

Vi börjar med att aktivera nödvändiga paket och aktivera data från en studie på patienter med tjocktarmscancer (coloncancer).

library(tidyverse) # för hantering av data
library(tidymodels) # grundpaketet för prediktionsmodellering
library(dials) # för att optimera hyperparametrar
library(finetune) # ytterligare funktioner för att optimera hyperparametrar

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.

library(survival)
data(cancer) # aktivera all dataset med cancerstudier
head(colon, 10) # inspektera första 10 rader i studien om coloncancer

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, likt föregående kapitel, modifiera colon data med följande steg nedan:

Vi kodar om status 1 till Yes och 0 till No och gör variabeln till en faktor med funktionen recode_factor().
Vi kodar om sex så det står Kvinna och Man istället för 0 respektive 1.
Vi tar bort kolumnerna time, id och study eftersom vi inte skall använda de i vår modell.
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.

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

Kartlägg och imputera missing data

Vi fortsätter genom att undersöka förekomsten av missing data i colon2. Detta gör vi med funktionen vis_miss() från paketet visdat, enligt följande:

library(visdat)
vis_miss(colon2)

Som framgår av resultatet ovan erhåller vi en figur där kolumnerna och raderna colon2 visualiseras i en rektangel. Det framgår att det saknas totalt 0.3% av alla värden i colon2. Även om detta är en ovanligt låg siffra väljer vi ändå att imputera dessa värden med funktionen missRanger() som finns i paketet med samma namn. Resultatet skickar vi till colon3, som nu är en dataframe utan missing data (alla missing-värden har imputerats).

library(missRanger)
colon3 <- missRanger(data=colon2, maxiter=10, pmm.k = 5)

vis_miss(colon3)

Skapa träningsdata, testdata och resampling (korsvalidering)

Vi delar nu colon3 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 på utfallet (strata=status), så att andelen som dör är lika stor i träningsdata och testdata.

# Skapar träningsdata och testdata med set.seed() först
set.seed(1)

colon_split <- initial_split(colon3, 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 att vi angett set.seed(1), vilket eliminerar slumpens effekt när vi genererar test- och träningsdata. Det innebär att vi kan generera precis samma observationer i träningsdata och testdata varje gång vi använder initial_split(). Detta kan vara önskvärt för att få resultat som kan reproduceras exakt.

Skapa recept

Vi skapar två recept. Det första receptet kallar vi base_recept. Se kommentarerna i koden nedan för att se vilka pre-processing steg som görs. Det andra receptet kallas extended_recipe och bygger på base_recept men har tilläget step_pca(). Om du inte är bekant med metoden PCA så kan du för tillfället bara konstatera att PCA är en metod som reducerar antalet numeriska prediktorer, vilket ibland leder till bättre modeller.

# Grundläggande recept
base_recipe <- recipe(status ~ ., data=colon_train) |> 
  # Normalisering av numeriska prediktorer
  step_normalize(all_numeric()) |> 
  # Skapar dummy-variabler av kategoriska variabler
  step_dummy(all_nominal_predictors()) |> 
  # Skapar tolerans för nya variabelnivåer i testdata
  step_novel(all_nominal_predictors())

# Utökat recept baserat på base_recipe
extended_recipe <-
  base_recipe |> 
  step_pca(all_numeric_predictors())

Skapa modeller

Vi skapar tre modeller med följande inställningar:

Logistisk regression.
- Hyperparametrar: inga.
- Metod: GLM (vanlig logistisk regression)
Random forest.
- Hyperparametrar: 2 st. Vi skapar vår grid manuellt med funktionen crossing().
- Metod: Ranger-paketet.
Xtreme gradient boosting.
- Hyperparametrar: 8 st. Vi skapar vår grid med funktionen grid_latin_hypercube()., som utifrån antal tillåtna kombinationer (size=n) maximerar vilka kombinationer som undersökas.
- Metod: Ranger-paketet.

# 1. LOGISTISK REGRESSION =====================================
glm_spec <-
  logistic_reg() %>%
  set_engine('glm') %>%
  set_mode('classification')


# 2. RANDOM FOREST ============================================
rf_spec <- rand_forest(mtry = tune(),
                       trees = tune()) |> 
  set_engine("ranger") |> 
  set_mode("classification")

# Med funktionen crossing() skapas en grid manuellt.
# Ange manuellt vilka värden som ska testas för varje parameter
rf_grid <- crossing(mtry = c(3, 6, 9),
                    trees=c(100, 200))


# 3. XTREME GRADIENT BOOSTING ================================
# Definiera modellen
xgb_spec <-
  boost_tree(tree_depth = tune(),
             mtry = tune(),
             trees = tune(),
             learn_rate = tune(),
             min_n = tune(),
             loss_reduction = tune(),
             sample_size = tune(),
             stop_iter = tune()) %>%
  set_engine('xgboost') %>%
  set_mode('classification')

# Skapa en grid med grid_latin_hypercube()
xgb_grid <- grid_latin_hypercube(
  trees(),
  min_n(),
  learn_rate(), 
  loss_reduction(),
  tree_depth(),
  sample_prop(),
  stop_iter(),
  finalize(mtry(), colon_train),
  size = 10)

Skapar workflow

Vi skapar nu en uppsättning med workflows med funktionen workflow_set(). Här skapas två listor. Första listan är preproc och där anges recepten. Nästa argument är models och där anges modellerna. Argumentet cross sätts till TRUE eftersom alla modeller ska testas med alla recept.

Vi behöver addera våra grids till vårt workflow med funktionen option_add(). Här anges vilken grid som ska adderas till ett specifikt workflow. Som framgår nedan är namnet på varje workflow dess pre-processor och model separerade med ett understreck (exempelvis base_xgb).

my_wf <- 
   workflow_set(
      preproc = list(base = base_recipe,
                     extended = extended_recipe),
      models = list(glm = glm_spec,
                    rf = rf_spec,
                    xgb = xgb_spec),
      cross = TRUE
   ) %>%
   # Adderar specifik grid-search till specifikt arbetsflöde
   option_add(id = "base_xgb",     grid = xgb_grid) %>%
   option_add(id = "extended_xgb", grid = xgb_grid) %>%
   option_add(id = "base_rf",      grid = rf_grid) %>%
   option_add(id = "extended_rf",  grid = rf_grid)

Kör modellerna

För att köra modellerna använder vi parallellisering, vilket innebär att vi utnyttjar flera av datorns kärnor om det kan göras. Parallelliseringen kräver paketen doParallel och parallel. Paketet tictoc används för att mäta tiden det krävs för att köra alla modellerna.

Tänk på att inte använda alla kärnor i din dator. Låt alltid 2 kärnor vara lediga för andra processer i datorn. I nedanstående fall har datorn 8 kärnor och därför används 6 kärnor till beräkningarna nedan. Antal kärnor specificeras med funktionen makePSOCKcluster().

library(doParallel)
library(parallel)
library(tictoc)

# Starta parallellisering
mitt_kluster <- makePSOCKcluster(6)
registerDoParallel(mitt_kluster) 
tic()

# Kör modellerna
my_wf_results <- 
   workflow_map(object=my_wf,
                "tune_grid",
                seed = 1,
                resamples = my_cross_validation,
                metrics = yardstick::metric_set(roc_auc),
                verbose = TRUE,
                control   = control_grid(save_pred = TRUE,
                           parallel_over = "everything",
                           save_workflow = TRUE)
                )
                
# Stoppa parallellisering
stopCluster(mitt_kluster)
toc()

Resultat

i	No tuning parameters. `fit_resamples()` will be attempted
i 1 of 6 resampling: base_glm
✔ 1 of 6 resampling: base_glm (24.4s)
i 2 of 6 tuning:     base_rf
✔ 2 of 6 tuning:     base_rf (2m 52.9s)
i 3 of 6 tuning:     base_xgb
✔ 3 of 6 tuning:     base_xgb (6m 40.6s)
i	No tuning parameters. `fit_resamples()` will be attempted
i 4 of 6 resampling: extended_glm
✔ 4 of 6 resampling: extended_glm (11.6s)
i 5 of 6 tuning:     extended_rf
✔ 5 of 6 tuning:     extended_rf (3m 51.6s)
i 6 of 6 tuning:     extended_xgb
✔ 6 of 6 tuning:     extended_xgb (17m 25.1s)

1893.746 sec elapsed

Körningen tog 1893 sekunder.

Jämför modellerna

Med funktionen autoplot() får vi en snabb överblick över modellernas förmåga att predicera utfallet, baserat på det precisionsmått vi valt (AUC ROC).

autoplot(my_wf_results)

Vi kan extrahera de bästa modellerna med följande kod

rank_results(my_wf_results,
             rank_metric = "roc_auc",
             select_best = TRUE) %>% 
   select(rank, mean, model, wflow_id, .config)

Resultat

   rank  mean model        wflow_id     .config              
  <int> <dbl> <chr>        <chr>        <chr>                
1     1 0.798 rand_forest  base_rf      Preprocessor1_Model6 
2     2 0.735 boost_tree   base_xgb     Preprocessor1_Model01
3     3 0.700 logistic_reg base_glm     Preprocessor1_Model1 
4     4 0.694 boost_tree   extended_xgb Preprocessor1_Model02
5     5 0.681 logistic_reg extended_glm Preprocessor1_Model1 
6     6 0.681 rand_forest  extended_rf  Preprocessor1_Model1

Den bästa modellen är random forest med grundreceptet (base_rf). Den bästa konfigurationen heter Preprocessor1_Model6.

För att se vilken konfiguration detta är skriver vi:

best_model <- my_wf_results %>%
  extract_workflow_set_result("base_rf") %>%
  select_best(metric = "roc_auc")

best_model

   mtry   trees    .config             
      9    200     Preprocessor1_Model6

modellen har mtry = 9 och trees = 200.

Erinra att träningen gjordes med korsvalidering, vilket gör att modellen aldrig tränats på all träningsdata (en fold har använts för att testa modellen, se Korsvalidering). Den bästa modellen tränas därför en sista gång på all träningsdata:

my_finalized_wf <- 
  my_wf_results |> 
  extract_workflow("base_rf") |> 
  finalize_workflow(best_model) |> 
  last_fit(split = colon_split,
           metrics = yardstick::metric_set(roc_auc))

Objektet my_finalized_wf innehåller den slutgiltiga modellen. Med funktionen collect_metrics() kan vi extrahera information om hur bra modellen är på testdata. Observera att vi inte behöver specificera testdata eftersom tidymodels använder testadata automatiskt (objektet my_finalized_wf har tillgång till testdata eftersom vi specifierade split=colon_split ovan).

collect_metrics(my_finalized_wf)

  .metric  .estimator .estimate .config             
  roc_auc  binary         0.792 Preprocessor1_Model1

ROC-kurvan har en area på 0.793.

Vi kan nu studera själva prediktionerna genom att använda funktionen collect_predictions():

pred_data <- collect_predictions(my_finalized_wf

Resultat

# A tibble: 558 × 6
   id               .pred_Yes .pred_No  .row status .config             
   <chr>                <dbl>    <dbl> <int> <fct>  <chr>               
 1 train/test split    0.873     0.127     2 Yes    Preprocessor1_Model1
 2 train/test split    0.0959    0.904     4 No     Preprocessor1_Model1
 3 train/test split    0.559     0.441     8 Yes    Preprocessor1_Model1
 4 train/test split    0.578     0.422    13 Yes    Preprocessor1_Model1
 5 train/test split    0.560     0.440    14 Yes    Preprocessor1_Model1
 6 train/test split    0.327     0.673    20 No     Preprocessor1_Model1
 7 train/test split    0.292     0.708    21 No     Preprocessor1_Model1
 8 train/test split    0.333     0.667    22 No     Preprocessor1_Model1
 9 train/test split    0.454     0.546    30 No     Preprocessor1_Model1
10 train/test split    0.605     0.395    32 Yes    Preprocessor1_Model1

ROC-kurvan skapas med följande kod:

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

För att klargöra om modellen är välkalibrerad använder vi funktionen cal_plot_breaks() från paketet probably. Observera att denna funktionen kräver att testversionen av probably är installerad. Den installeras med följande kommando:

install.packages("remotes")
remotes::install_github("tidymodels/probably")
library(probably)

pred_data %>% 
  cal_plot_breaks(status, .pred_Yes)