08-advanced-models.qmd

# מודלים מונחים {#sec-advanced-models}

## הקדמה

ניתן לחלק את עולם המודלים הסטטיסטיים לשני סוגים, המובחנים ביניהם בדרך שבה בונים אותם:

-   מודלים מונחים (Supervised learning)

-   מודלים לא מונחים (Unsupervised learning)

מודלים מונחים הינם מודלים שנבנים תוך שימוש במשתנה תלוי (כפי שראינו בפרק הקודם של רגרסיה לינארית), ולעומת זאת, מודלים לא מונחים הינם מודלים שבהם אין משתנה תלוי (נניח לדוגמה כשרוצים לחלק את מרחב התצפיות לאשכולות שונים).

בפרק זה נתמקד במידול מונחה, מנקודת מבט מעשית, תוך שימוש בחבילות במשפחה של `tidymodels`.

::: callout-note
## קריאה נוספת לצורך העמקה בתיאוריה

ישנן משפחות רבות של מודלים מונחים, בין המוכרים כלולים גם מודלי שפה (כגון ChatGPT), אך גם מודלים לזיהוי נטישה, זיהוי הונאות, רגרסיות למיניהן, למידה עמוקה, ועוד.

הספר הנוכחי אינו מתמקד בתיאוריה העומדת בבסיס המודלים השונים, אלא מהווה עזר מעשי להפעלת המודלים. אתם מוזמנים לעיין בספרים הבאים לצורך העמקה תיאורתית:

-   The Elements of statistical learning [@hastie2009]

-   An Introduction to Statistical Learning [@james2021]
:::

## תיאור שלבי העבודה

כפי שציינו במבוא לספר זה, תהליך ניתוח הנתונים כולל טרנספורמציה על משתנים, וויז'ואליזציות של המשתנים והקשרים שלהם, ומידול. שלב המידול עצמו כולל מספר תתי-שלבים:

1.  חלוקה אקראית של המדגם לקבוצת אימון (training set), קבוצת ולידציה (validation set) וקבוצת מבחן (test set).
2.  בניית מתכון (recipe) להכנת הדאטה למודל. בניית המתכון משמשת להגדרת כל הטרנספורמציות הדרושות לביצוע על הדאטה, לפני שהוא נכנס למודל.
3.  הגדרת מודל/ים: בשלב זה מגדירים את המודלים בהם הולכים לעשות שימוש. גם מבחינת האלגוריתם של המודל (סוג המודל) וגם מבחינת מטרת המודל (רגרסיה לעומת סיווג).
4.  הגדרת תהליך אימון המודל (workflow) ובניית המודל הראשוני (fitting).
5.  כיוונון המודל (tuning).
6.  הערכת ביצועי המודל (evaluation).
7.  שימוש במודל לחיזוי תצפיות חדשות (prediction).

בשלב 1 אנחנו מפרקים את הנתונים לשלוש קבוצות, כאשר מטרת קבוצת האימון הינה לבנות את המודל, מטרת קבוצת הולידציה הינה לכוון את "היפר-פרמטרים" של המודל (hyper-paramters)ומטרת קבוצת המבחן היא לבחון את שגיאת המודל (על דאטה שהמודל לא תלוי בו משום שלא נלמד על בסיסו). קבוצת האימון תשתמש אותנו בשלבים 2-5, קבוצת הולידציה תשמש אותנו בשלב 6, וקבוצת המבחן תשמש אותנו בשלב 7 בלבד. שלבים 2-6 עשויים לחזור על עצמם מספר פעמים במהלך פיתוח מודלים (עם גרסאות שונות למודל או מודלים שונים לגמרי). השלב האחרון הינו השלב שבו אנחנו מפיקים ערך מהמודל, ומשתמשים בו לחיזוי עבור תצפיות חדשות.

כעת נפרק כל אחד מהשלבים ונסביר עליו בפירוט. על מנת להדגים את השלבים נשתמש בדאטה שמלווה אותנו מתחילת הספר (Palmer penguins), כאשר הבעיה שנפתור בפרק זה היא סיווג של תצפית של פינגויין לזכר או נקבה לפי כל המשתנים האחרים שעומדים לרשותנו.

## חלוקת המדגם

ראשית נחלק את המדגם לשלוש הקבוצות: אימון, ולידציה, ומבחן.

אנחנו קוראים את החבילות הנדרשות, מגדירים Seed, כך שנוכל להגריל חלוקה אקראית אך עקבית (תמיד כשנריץ את הקוד נקבל את אותה החלוקה האקראית), ונפצל את הדאטה לשלושת הקבוצות.

```{r initial splitting of the data}
#| warning: false
#| message: false
library(tidymodels) # reading the meta-package tidymodels
library(palmerpenguins) # reading the palmer penguins data set which include penguins

set.seed(42) # set an initial seed

# Modify the dependent var to a dummy
penguins_new <- penguins %>% 
  mutate(sex = factor(sex))

# Create the split
penguin_split <- initial_validation_split(penguins_new, 
                                          prop = c(0.6, 0.2), 
                                          strata = sex)
penguin_split
```

ראשית לפני התחלה, שינינו את המשתנה sex למשתנה פקטור. שינוי זה יקל עלינו בהמשך בהפעלת מודלים לסיווג. משתנים קטגוריאליים אחרים המהווים משתנים בלתי-תלויים (משתנים מסבירים) יטופלו בהמשך.

כפי שניתן לראות, המדגם חולק לקבוצת אימון בגודל של 205 תצפיות (המהווה כ-60% מהדאטה המקורי), קבוצת ולידציה בגודל 69 תצפיות (המהווה כ-20% מהדאטה המקורי), וקבוצת מבחן בגודל 70 תצפיות (המהווה כ-20%). גודל הקבוצות הוגדר באמצעות שימוש בארגומנט `prop`. מקובל לקבוע את קבוצת האימון כקבוצה הגדולה ביותר. ערכים מקובלים לקבוצת האימון והולידציה הם בדרך כלל 70-80% במצטבר, כאשר גודל קבוצת המבחן יהווה כ-20-30% מהדאטה בהתאמה. עם זאת, אין הגדרות חד-משמעיות בהיבט זה, וניתן לחרוג ממספרים אלו כתלות בצורך.

השתמשנו בארגומנט `strata` על מנת לבצע את הדגימה באופן ששומר על הפרופורציות של מין התצפיות בתוך כל תת-קבוצה. ארגומנט זה שימושי במיוחד כאשר הקבוצות אינן מאוזנות (קבוצה מסוימת קטנה משמעותית מאחרות).

על מנת לחלץ את הדאטאות עצמם של כל אחת מהקבוצות נשתמש בפונקציות הבאות:

```{r extract sets}
penguin_train <- training(penguin_split)
penguin_valid <- validation(penguin_split)
penguin_test <- testing(penguin_split)
```

את המשך הלמידה על הדאטה, כולל בניית המודל, נבצע על קבוצת האימון בלבד `penguin_train`.

::: callout-tip
## תרגיל: הבנת הדאטה לפני צלילה למודל

בדרך כלל לאחר פיצול הדאטה נבצע ויז'ואליזציות שונות כדי להבין איך משתנים שונים משפיעים על המשתנה התלוי (לפני שרואים את השפעתם במודל עצמו). בשלב זה העבודה שתוארה ב[פרק @sec-ggplot2] תהא משמעותית מאוד. לצורך תרגול, בחנו כיצד המשתנים השונים שבדאטה שלנו (penguin_train) משפיעים על מין התצפית.

1.  השתמשו בויז'ואליזציות שונות כפי שנלמדו והודגמו בפרק על ויז'ואליזציות. לאור הבדיקה, אילו משתנים הייתם מצפים שיהיו משמעותיים במודלים שנפתח?
2.  ישנן שלוש תצפיות בעלות ערך חסר (ללא מין) בקבוצת האימון. לא נוכל להשתמש בתצפיות אלו לצורך אימון המודל - הסבירו מדוע.
:::

## בניית מתכון

כפי שמיד תראו, חבילת tidymodels מאמצת ז'רגון של מטבח כמו בניית מתכון - recipe, סחיטת מיצים - juice, אפייה - bake, גזר לבן - parsnip, ועוד.\
בשלב זה נבנה "מתכון" שהמטרה שלו להגדיר את השינויים שעובר הדאטה לפני שהוא נכנס למודל. במובן מסויים, התהליך מאוד דומה לתהליך שתואר ב[פרק @sec-data-munging-tidyverse] על הכנת נתונים, אך יש הבדל מהותי: מתכון זה יופעל אוטומטית בכל פעם שנרצה להפעיל את המודל מחדש על תצפיות חדשות, וכן הוא יופעל גם במהלך שלב הכיוונים (וגם ניתן יהיה להגדיר היפר-פרמטרים כחלק מהמתכון, ולאפשר לשלב הכיוונון לבדוק גם אותם).

בבניית המתכון הבא נפעיל מספר צעדים, החל מהפקודה `recipe` לאתחול המקום ודרך צעדים שונים (`step`) שמיד נסביר אותם.

```{r building a recipe}
# Excuse the pun (for the name...)
penguin_recipe <- recipe(sex ~ ., data = penguin_train) %>% 
  step_naomit(bill_length_mm:sex) %>% 
  step_normalize(all_numeric_predictors()) %>% 
  step_dummy(species, island) %>% 
  step_select(-year) %>% 
  step_corr(bill_length_mm:body_mass_g)

penguin_recipe
```

במתכון שלעיל הגדרנו ארבעה צעדים (פונקציות שמתחילות ב-`step_*`).

בצעד הראשון אנחנו משמיטים ערכים חסרים מהדאטה באמצעות `step_na_omit`. בדאטה שלנו תצפיות בעלות ערך חסר מכילות ערך חסר בכל המשתנים הרלוונטיים, ולכן בחרנו להשמיט אותן לגמרי. באופן כללי, בבעיות אחרות כדאי לנסות לזקוף ערכים חסרים על ידי שימוש בצעדים כגון `step_impute_bag`, `step_impute_knn`, `step_impute_linear`, ועוד לפני שמסירים את התצפיות לגמרי. ביצוע צעד של זקיפת נתונים יאפשר להתמודד עם חיזוי של תצפיות חדשות גם אם הן מגיעות עם ערך חסר, בעוד שצעד הסרה ימנע מאיתנו לספק תחזית.

הצעד הבא, `step_normalize`, הופך את כל המשתנים המספריים (`all_nominal_predictors()`) למשתנים מנורמלים (כלומר בעלי תוחלת 0 וסטיית תקן 1). צעד זה גם ישמור את התוחלת וסטיית התקן המחושבת על קבוצת האימון, וישתמש בהן לנרמול של תצפיות חדשות בשלב החיזוי (חשבו למה).

הצעד הבא, `step_dummy`,הופך משתנים מסוג מחרוזת או פקטור למשתני דמי (בעלי ערכי 0 או 1). במקרה שלנו אלו משתני סוג הפינגויין, האי שבו נמדדה התצפית, ומין הפינגויין (שכבר הומר למשתנה פקטור בתחילת הקוד לפני הפיצול הראשוני לקבוצות). ניתן לשים לב שהמשתנה אי התצפית פוצל לשני משתני דמי (במקור ישנן שלוש רמות Dream, Torgersen, ו-Biscoe, כאשר השניים הראשונים הם 0 המשמעות היא שהתצפית הגיעה מהשלישי), וכך גם המשתנה של סוג הפינגויין (שפוצל לשני משתני דמי Chinstrap, Gentoo, והשלישי Adelie אינו בדאטה ומשתמע מהשניים האחרים).

הצעד `step_select`מוריד את משתנה השנה (שמתעד מתי נמדדה התצפית). הנחת העבודה בהסרת המשתנה היא שהשנה שבה נמדדה התצפית ומין הפינגויינים אינם קשורים סטטיסטית (בהנחה שהתפלגות מין הפינגויינים נשארת קבועה בין השנים). הפעלת הפוקנציה מאוד דומה לפונקציה `select` שבה דנו בפרק הכנת הנתונים.

הצעד האחרון שבחרנו לבצע הינו `step_corr`. צעד זה בוחן קורלציות בין המשתנים ובמידה והוא מזהה משתנה שבעל קורלציות גבוהות למשתנים אחרים, הוא מסיר אותו (במקרה זה לצעד זה אין השפעה, ואף משתנה אינו מוסר משום שאין מולטיקולינאריות, כפי שמיד נראה).

נציין שרשימת הצעדים האפשריים גדולה במיוחד ומבוססת על best practices של הנדסת נתונים (כפי שנקראים בעגה המקצועית - feature engineering). מומלץ להיעזר [ברשימת הצעדים](https://recipes.tidymodels.org/articles/recipes.html) המתעדכנת באתר הרשמי. כמו כן, תמיד ניתן להגדיר צעדים חדשים נוספים, לדוגמה על ידי שימוש ב-`step_mutate` להגדרת משתנה נוסף שהוא תוצאת חישוב של משתנים אחרים, או באופן מתקדם יותר על ידי [פיתוח של פונקציה חדשה](https://www.tidymodels.org/learn/develop/recipes/).

על מנת לבחון את תוצאת המתכון, ניתן להשתמש בפונקציות הכנה ואפייה, באופן הבא:

```{r prep and bake}
penguin_recipe %>% 
  prep() %>% 
  bake(new_data = NULL) %>% 
  glimpse()
```

באמצעות הפונקציה `prep` אנחנו מנחים את הפעלת הצעדים על בסיס קבוצת האימון, והפונקציה `bake` מריצה אותם בפועל. אולי זה נראה ככפילות מסוימת אבל בעצם יכלנו להריץ את `bake` על כל דאטה (כגון קבוצת הולידציה, קבוצת המבחן, או תצפיות חדשות לגמרי). אם מפעילים את `bake` עם `new_data = NULL` אז היא פשוט תפעיל הכל על קבוצת האימון.

## הגדרת המודל

הבעיה עמה אנחנו מתמודדים בפרק זה היא בעיית סיווג (classification), בהתאם לכך, ישנם מודלי סיווג שבאפשרותנו להשתמש בהם, כגון: רגרסיה לוגיסטית, עצי החלטה, יערות אקראיים, boosting, bagging, ועוד. כפי שציינו ניתן להתעמק במודלים השונים, בשני הספרים שצויינו בתחילת הפרק.\
נדגים הגדרה של שלושה מודלים שונים: רגרסיה לוגיסטית, יער אקראי, ו-boosting.

```{r definition of models}
# Logistic regression
penguin_logistic <- logistic_reg(
  mode = "classification",
  engine = "glm",
)

# Random forest
penguin_forest <- rand_forest(
  mode = "classification",
  engine = "ranger",
  mtry = 3,
  trees = 50,
  min_n = 15)

# Boosting
penguin_boost <- boost_tree(
  mode = "classification",
  engine = "xgboost")
```

שימו לב לכמה נקודות מעניינות:

-   עד כה כל הפעולות שעשינו היו זהות לכל המודלים (פיצול והכנת הדאטה), וזו הפעולה הראשונה שבה אנחנו מבחינים בין שלושה מודלים שונים. שלושת המודלים יקבלו את אותו הדאטה בשלב האימון. העוצמה הגדולה של `tidymodels` היא ביכולת שלה לאפשר לנו ממשק לחבילות מרובות (כפי שמיד נסביר, אנחנו מפעילים פה פונקציות מחבילות `stats`, `ranger`, ו-`xgboost`), תוך שימוש בשפה אחידה. במקור, לכל חבילה שבה אנחנו משתמשים בהגדרה לעיל יש ממשק שונה שצריך ללמוד אותו אם מפעילים ישירות את החבילות, אבל במקרה זה `tidymodels` מספקת לנו מעטפת אחידה.

-   כל מודל מוגדר עם ארגומנט `mode` שמבחין בין בעיות סיווג לבעיות רגרסיה. במקרה זה אנחנו עובדים עם בעיית סיווג, ולכן הארגומנט הוגדר כ-`mode = classification`.

-   בחבילת `parsnip` שהינה חבילה בתוך `tidymodels` ישנם ממשקים שונים להפעלת פונקציות מידול שונות, ובמקרה זה השתמשנו בשלושה ממשקים:

    -   הפונקציה `logistic_reg` שבאמצעות `engine = "glm"` הנחנו אותה לעבוד עם פונקציית `glm` (generalized linear models) שנמצאת בחבילת הבסיס `stats`. מעיון בתיעוד של `logistic_reg` תוכלו לראות שהיא מספקת ממשקים למעל עשר חבילות שונות של רגרסיה לוגיסטית (הכללות שונות של המודל הבסיסי או אלגוריתמים בעלי יעילות טובה יותר).

    -   הפונקציה `rand_forest` גם היא מספקת ממשקים לפוקנציות ואלגוריתמים שונים של יערות אקראיים, ובמקרה זה הנחנו את הפונקציה להוות ממשק לחבילת `ranger`.

    -   הפונקציה `boost_tree` שמייצרת מודלים המורכבים מהרבה תתי מודלים אחרים (ensemble), ובמקרה זה תעבוד כממשק לחבילת `xgboost`.

-   כל פונקציית ממשק כזו תפעיל פונקציות מחבילות אחרות (בדוגמה שלנו מחבילת stats, ranger, ו-xgboost), ולכן עלינו לוודא שהתקנו את החבילות הנדרשות.

-   פונקציות הממשק יודעות להעביר ארגומנטים שונים שרלוונטיים לפונקציות שהן מפעילות. לדוגמה לחבילת `ranger` אנחנו מעבירים את הארגומנט `mtry = 3` שמסביר מכמה משתנים יש לבנות כל עץ שמשתתף ביער שלנו. עיון בתיעוד הפונקציות בחבילות של האלגוריתמים יראה אילו ארגומנטים קיימים.

## הגדרת תהליך ובניית המודל

כעת, על מנת להפעיל את האלגוריתמים השונים לצורך בניית המודלים, נגדיר תהליך בניה ונבנה את המודל בפועל באופן הבא:

```{r define workflow}
penguin_workflow <- workflow() %>% 
  add_recipe(penguin_recipe)

penguin_logistic_fit <- penguin_workflow %>% 
  add_model(penguin_logistic) %>% 
  fit(data = penguin_train)

penguin_forest_fit <- penguin_workflow %>% 
  add_model(penguin_forest) %>% 
  fit(data = penguin_train)

penguin_boost_fit <- penguin_workflow %>% 
  add_model(penguin_boost) %>% 
  fit(data = penguin_train)
```

לפעמים מודל הוא "קופסה שחורה" ולא ניתן להתעמק בו (לדוגמה, מודלי למידה עמוקה, ובפרט מודלי שפה כגון ChatGPT הם קופסה שחורה, בהיבט שקשה לחקור את המודל אך קל לחזות בתוצאותיו).\
עם זאת, לעיתים ההתעמקות בתכונות של מודל יכולה לשפוך אור על הקשר שבין המשתנים (כפי שראינו ב[פרק @sec-linear-regression]. לשם כך עומדות לרשותנו פונקציות שונות המאפשרות לחלץ את המודל וללמוד את תכונותיו. נדגים חקירה של מודל הרגרסיה הלוגיסטית ושל מודל ה-Boosting.

### חקירה של מודל הרגרסיה הלוגיסטית

```{r extract fit and summary of logistic}
logistic_model_extract <- penguin_logistic_fit %>% 
  extract_fit_engine() 

logistic_model_extract %>% 
  summary()
```

הפלט דומה לפלט הפקודה `summary` כאשר מריצים אותה על מודל הרגרסיה הלינארית שראינו בפרק הקודם אך עם מספר הבדלים, לדוגמה:

-   במקום Residuals אנחנו רואים Deviance Residuals (חישוב שונה של שאריות המבוסס על פונקציית נראות במקום הפרש בין ערך התצפית בפועל לערך החזוי של התצפית).

-   במקום מבחן $F$ וערך של $R^2$ אנחנו רואים ערך של AIC (Akaike Information Criteria) - קריטריון אחר המאפשר להשוות בין מודלים.

-   משמעות המקדמים שונה (Estimate): ברגרסיה לוגיסטית, האקספוננט של המקדם הוא ה-Odds Ratio, כלומר פי כמה שינוי של יחידה במשתנה המסביר מעלה את הסבירות לשינוי הסיווג.

נציג את ערכי האקספוננט של המקדמים. נשתמש בפונקציה `tbl_regression` מחבילת `gtsummary`, עם הארגומנט `exponentiate = TRUE` על מנת להציג את האקספוננט של המקדמים במקום המקדמים עצמם.

```{r exp estimates}
#| warning: false
#| message: false
gtsummary::tbl_regression(logistic_model_extract, exponentiate = TRUE)
```

לדוגמה, עיון במקדם של אורך המקור (`bill_length_mm`) שערכו כ-3.62 (בטבלת הפלט הגולמית), מעיד על כך ששינוי של 1מ"מ באורך המקור, מעלה את ה-Odds Ratio (היחס בין הסבירות להשתייך לקבוצת הזכרים לסבירות להשתייך לקבוצת הנקבות) פי 37 (עמודת OR בטבלת הפלט המעובדת). ממצא זה מובהק סטטיסטית (p-value\<0.001).

המשתנה המשמעותי ביותר כפי שעולה מתוצאות מודל הרגרסיה הלוגיסטית הינו עומק המקור (`bill_depth_mm`) עם ערך Odds Ratio של 183.

::: callout-tip
## תרגיל: אינטראקציות ברגרסיה לוגיסטית

בפרק הקודם שדן ברגרסיה לינארית הדגמנו כיצד פועלות אינטראקציות ברגרסיה לינארית. כמו כן, ראינו שישנם קשרים בין הזנים השונים של הפינגויינים לבין תכונותיהם, וכעת ראינו גם שיש קשר בין מין הפינגויינים לתכונותיהם.

1.  העריכו כיצד הממצאים הללו עשויים להשפיע על ביצועי מודל הרגרסיה הלוגיסטית שהוצג?
2.  אילו אינטראקציות הייתם שוקלים להוסיף למודל על מנת לשפרו?
3.  הוסיפו את האינטראקציות הללו על ידי שימוש בפקודה `step_interact` במתכון שהוצג, והריצו מחדש את מודל הרגרסיה הלוגיסטית.
4.  נתחו את תוצאות המודל שקיבלתם לאחר האינטראקציה, מה אתם מסיקים מהמודל החדש?
:::

מומלץ להעמיק בתכונות הרגרסיה הלוגיסטית, לדוגמה בספר של @hastie2009.\
נדגים כעת כיצד אנו מחלצים את מודל ה-boosting ולומדים ממנו על השפעת המשתנים על מין הפינגויין.

### חקירה של מודל ה-Boosting

ראשית, נחלץ את מודל ה-boosting:

```{r extract fit and summary of boosting}
bst_model_extract <- penguin_boost_fit %>% 
  extract_fit_engine()

bst_model_extract
```

בשונה ממודל הרגרסיה הלוגיסטית, למודל boosting אין "מקדמים" במובן של מקדמים של מודל רגרסיה לוגיסטית. עם זאת, ניתן להשתמש בו כדי לחלץ את חשיבות המשתנים השונים, קרי, דירוג של המשתנים על פי מידת ההשפעה שיש להם על רמת הדיוק של העצים של המודל (תתי-המודלים מהם מורכב מודל ה-boosting). חילוץ נתונים אלו דורש היכרות והפעלה של פקודות מחבילת `xgboost`. במקרה זה אנחנו מציגים את ה-Gain של המשתנים השונים, כלומר התרומה היחסית שלהם בבעיית הסיווג (אחוזים גבוהים יותר משמעם תרומה גבוהה יותר למודל החיזוי).

```{r extract variable importance}
importance_matrix <- xgboost::xgb.importance(model = bst_model_extract)
importance_matrix %>% 
  mutate(Feature = forcats::fct_reorder(Feature, Gain)) %>% 
  ggplot(aes(y = Feature, x = Gain)) + 
  geom_col() + 
  ggtitle("Feature Gain in the Penguin Sex Boosting Model") +
  theme_bw()
```

כעת נעבור לשלב כיוונון המודלים, ונלמד מהם היפר-פרמטרים (Hyper parameters).

## כיוונון מודלים

ככל שעולה מידת המורכבות של מודלים, כך מתווספים להם פרמטרים הנקראים "היפר-פרמטרים" (Hyper parameters). יש הבדל מהותי בין המונח "פרמטר" (המתייחס לערך של האוכלוסיה שאנחנו אומדים, כגון תוחלת או סטיית תקן), לבין היפר-פרמטר. כאשר המונח היפר-פרמטר מתייחס לפרמטרים של האלגוריתם (ולא של האוכלוסיה). להלן מספר דוגמאות להיפר-פרמטרים:

-   מספר העצים שישתתפו במודל (נניח מסוג יערות אקראיים או מסוג boosting).

-   מספר האיטרציות המירבי שאלגוריתם צריך לעשות במהלך אימון המודל.

-   מספר הפיצ'רים המירבי שמורשים להשתתף בכל עץ במודל יערות אקראיים.

-   עומק העץ המירבי במודל יערות אקראיים.

לעיתים, ההיפר-פרמטר יהיה תוצר של צעדי המתכון שהכנו. לדוגמה אם השתמשנו בדיסקרטיזציה (`step_discretize` הפיכה של משתנה רציף למשתנה בדיד בעל מספר ערכים מוגבל), אז מספר הערכים שאנחנו מגדירים למשתנה הבדיד הינו היפר-פרמטר (הארגומנט `num_breaks` בדוגמה). ערכים שונים שלו עשויים להניב מודלים בעלי ביצועים שונים (מעט מדי ערכים יפגעו ביכולת של המשתנה להסביר את המשתנה התלוי, ויותר מדי ערכים יכניסו למודל הרבה פיצ'רים שיביאו לתופעה של התאמת-יתר over-fitting, גם היא תניב מודל בעל ביצועיים נחותים).

על כן, מהסיבות הללו ועוד, נשאלת השאלה- איך נכוון (tune) את המודל שלנו, על מנת שנצליח לבחור בסט הערכים המיטבי (שיביא את ביצועי המודל הטובים ביותר). מדובר בבעיה מורכבת שכן ברגע שיש לנו שילוב של היפר-פרמטרים, ישנם הרבה אפשרויות לבחינה (לדוגמה, חמישה פרמטרים שלכל אחד חמישה ערכים אפשריים מהווים $5^5$ צירופים שונים.

נדגים כעת כיצד אנחנו מכווננים את מודל היער האקראי בהתייחס להיפר-פרמטרים mtry, trees, n_min המתארים את מספר המשתנים המסבירים שמופיעים בכל עץ, מספר העצים הכללי במודל, וגודל המדגם המינימלי לכל קודקוד בעץ (שמתחתיו לא ניתן לפצל את הקוקוד).

כיוונון המודל מתבצע על ה-Validation set, ולבסוף כאשר נמצא את סט ערכי ההיפר פרמטרים המיטבי, נבצע אימון מחודש בהתאם לסט הערכים המיטבי, על גבי ה-Train set, ונאמוד את השגיאה באמצעות ה-Test set.

הפקודה הבאה מגדירה מודל, שבו רכיבים מסוימים הינם רכיבים אשר יעברו כוונון. הפונקציה `tune()` משמשת כמעין "ממלא מקום". הפונקציות בהן נשתמש בהמשך ידעו להחליף את `tune()` בערכים שונים במסגרת החיפוש המתבצע בתהליך הכיוונון.

```{r defining the tune placeholders}
# Random forest - tunning
# We'r tunning on the following hyper-parameters (mtry, trees, min_n)
penguin_forest_tune_spec <- rand_forest(
  mode = "classification",
  engine = "ranger",
  mtry = tune(),
  trees = tune(),
  min_n = tune())

penguin_forest_tune_spec
```

כעת נגדיר את רשת החיפוש (Search grid). המונח רשת מבטא את העובדה שישנם מספר היפר-פרמטרים וההצלבה ביניהם מייצרת מעין "רשת" שבה מתבצע החיפוש.

```{r define the search grid}
penguin_grid <- grid_regular(mtry(range = c(1, 5)), 
                             trees(), 
                             min_n(),
                             levels = 3)

penguin_grid
```

לחלק מהפונקציות בהם השתמשנו (`trees`, `min_n`) יש ערכי ברירת מחדל המגדירים את טווח החיפוש. \
עבור `mtry` הגדרנו את טווח החיפוש. הארגומנט `levels` מגדיר עבור כמה ערכים בטווח זה יבוצע החיפוש. כלומר, ברשת שהגדרנו יבוצע חיפוש ברשת של $3^3=27$ נקודות (כל ההצלבות האפשריות בין ההיפר-פרמטרים).

הפקודה הבאה תחלק את קבוצת ה-validation ל-`v=5` קבוצות באקראי. על כל אחת מקבוצות אלו יבנו 27 מודלים (כלומר בסך הכל 135 מודלים יבנו).

```{r define the vfold cv sets}
penguin_folds <- penguin_valid %>% 
  na.omit() %>% 
  vfold_cv(, v = 5, strata = sex)
```

השתמשנו ב-`na.omit` על מנת להשמיט ערכים חסרים מהמשתנה `sex`. כעת נבצע את ההרצה בפועל של בניית כלל המודלים:

```{r run the validation}
penguin_forest_tune_wf <- workflow() %>% 
  add_model(penguin_forest_tune_spec) %>% 
  add_recipe(penguin_recipe)

penguin_forest_tune_res <- penguin_forest_tune_wf %>% 
  tune_grid(
    resamples = penguin_folds,
    grid = penguin_grid
  )

penguin_forest_tune_res
```

כעת ניתן לנתח את רמת הדיוק כפונקציה של ערכי ההיפר-פרמטרים:

```{r hyper-parameter analysis}
penguin_forest_metrics <- penguin_forest_tune_res %>% 
  collect_metrics() %>% 
  mutate(mtry = glue::glue("mtry={mtry}")) %>% 
  mutate(across(c(mtry, min_n), as.factor))

penguin_forest_metrics

penguin_forest_metrics %>% 
  filter(.metric == "accuracy") %>% 
  ggplot(aes(x = trees, y = mean, color = min_n)) + 
  geom_line() + 
  geom_point() +
  facet_wrap(~ mtry) + 
  theme_bw() + 
  ylab("Accuracy") + 
  ggtitle("Model accuracy as a function of trees, min_n, and mtry")

penguin_forest_metrics %>% 
  filter(.metric == "roc_auc") %>% 
  ggplot(aes(x = trees, y = mean, color = min_n)) + 
  geom_line() + 
  geom_point() +
  facet_wrap(~ mtry) + 
  theme_bw() + 
  ylab("AUC") + 
  ggtitle("Model AUC as a function of trees, min_n, and mtry")
```

::: callout-tip
## מהם AUC ו-Accuracy?

TBD
:::

```{r select best}
penguin_forest_tune_res %>% 
  select_best(metric = "accuracy")

penguin_forest_tune_res %>% 
  select_best(metric = "roc_auc")
```

מניתוח שני התרשימים, נראה שהדיוק הגבוה ביותר במונחי accuracy מתקבל בעבור סט הערכים:

-   `min_n`= 21

-   `mtry`=5

-   `trees`=2000

במונחי AUC_ROC הדיוק הגבוה ביותר מתקבל כאשר מספר העצים הוא 1000 וערך min_n הינו 2, אך ביצועים אלו דומים גם לביצועים בעבור סט הערכים לעיל, ולכן יש העדפה קלה לבחירתו.

ייתכן שאם היינו מגדילים את מספר העצים במודל היתה מתקבלת רמת דיוק גבוהה אף יותר (משום שנראה שיש מגמת עליה קלה ברמת הדיוק כאשר מספר העצים גדל).

כעת נאמן את המודל עם ערכים אלו.

## הערכת ביצועי המודל

TBD

## שימוש במודל לחיזוי תצפיות חדשות

TBD

## קריאה נוספת

לצורך קריאה והעמקה נוספת בתהליך המידול, מומלץ לקרוא את [@kuhn2022tidy] המפרט על השימוש בחבילת `tidymodels` עם דוגמאות רבות. הספר [זמין באינטרנט](https://www.tmwr.org/).

## סיכום

::: end-page