Data Science עם R - רגרסיה לוגיסטית וסיווג
עדי שריד / adi@sarid-ins.co.il
ביחידה זו נדון ברגרסיה לוגיסטית.
רגרסיה לוגיסטית דומה בהרבה מובנים לרגרסיה לינארית: היא מתבססת על מבנה לינארי “מסוים”, וגם בהיבט הפרקטי, ב-R מריצים אותה בצורה מאוד דומה, ומפרשים אותה בצורה דומה. אבל יש לה הבדל משמעותי - במקום לחזות ערך מספרי “רחב” היא חוזה הסתברות, כלומר ערך מספרי בין 0 ל-1. לכן היא מתאימה לסיווג - נניח ההסתברות שלקוח ינטוש, או שליד פוטנציאלי יקנה את המוצר.
מה הבעיה עם רגרסיה לינארית לבעיות סיווג?
בבעיות סיווג, אנחנו מנסים לחזות את ההסתברות למאורע מסוים. \[p(Y=1|X)\] או בכתיב מקוצר, פשוט \[p(X)\] לכן אם ננסה להשתמש במודל רגרסיה לינארית נקבל את הנוסחה הבאה:
\[p(X) = \beta_0 + \beta_1X_1+\ldots\beta_p X_p\]
כאשר \(Y\) הוא המשתנה התלוי (לקוח ינטוש), ו-x הם המשתנים המסבירים. אבל אם נאמוד את פרמטרי ה-\(\beta\) ונבנה מודל רגרסיה לינארית, בהחלט ייתכן שנקבל עבור ערכים מסוימים הסתברויות שליליות או גבוהות מ-1!
זה כמובן לא ייתכן, אפשרות אחת היא להשתמש בפונקציית “עיגול” ולהפוך ערכים שליליים ל-0 וערכים מעל 1 ל-1. עם זאת, אם מראש היינו מתאימים מודל שתוצאותיו מוגבלות בין 0-1 ייתכן שביצועיו היו טובים יותר.
כאן נכנסת לביטוי פונקציית ה-logit.
\[p(X)=\frac{e^{\beta_0 + \beta_1 X_1 + \ldots + \beta_p x_p}}{1 + e^{\beta_0 + \beta_1 X_1 + \ldots + \beta_p x_p}}\]
פונקציה זו תמיד מקבלת ערכים בין 0 לבין 1. הרכיב שבחזקת \(e\) הוא לינארי, ודומה לרכיב שבו השתמשנו ברגרסיה לינארית.
בהימורי ספורט (מרוצי סוסים וכדומה) מקובל למדוד את היחס לזכיה (ולא את ההסתברות לזכיה). היחס לזכיה במונחים שהגדרנו יהיה:
\[\frac{p(X)}{1-p(X)} = e^{\beta_0 + \beta_1 X_1 + \ldots}\]
זה מלמד אותנו משהו חשוב בנוגע למשמעות של המקדמים \(\beta\). להזכירכם, ברגרסיה לינארית, מקדמים אלו שימשו אותנו כדי למדוד את השינוי הממוצע במשתנה התלוי \(y\) כאשר המשתנה הבלתי תלוי גדל ביחידה אחת. ברגרסיה לוגיסטית, הגודל \(e^{\beta_1}\) מלמד אותנו פי כמה ישתפר הסיכוי כאשר המשתנה הבלתי תלוי \(X_1\) גדל ביחידה אחת.
כדי למצוא את המקדמים ברגרסיה לוגיסטית, משתמשים בשיטה הנקראת Maximum Likelihood, אך לא נתמקד בתיאוריה מאחורי השיטה.
כדי להפעיל מודל של רגרסיה לוגיסטית ב-R כל מה שנדרש הוא לעשות כמה שינויים קלים בפקודות המוכרות:
library(tidyverse)## Warning: package 'tidyverse' was built under R version 3.4.4## -- Attaching packages ---------------------------------------------------------------- tidyverse 1.2.1 --## <U+221A> ggplot2 3.0.0 <U+221A> purrr 0.2.5
## <U+221A> tibble 1.4.2 <U+221A> dplyr 0.7.6
## <U+221A> tidyr 0.8.1 <U+221A> stringr 1.2.0
## <U+221A> readr 1.1.1 <U+221A> forcats 0.2.0## Warning: package 'ggplot2' was built under R version 3.4.4## Warning: package 'tibble' was built under R version 3.4.4## Warning: package 'tidyr' was built under R version 3.4.4## Warning: package 'purrr' was built under R version 3.4.4## Warning: package 'dplyr' was built under R version 3.4.4## -- Conflicts ------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()# read the no show data of patients
appointments <- read_csv("data-files/Medical_Appointments_No_Shows_KaggleV2-May-2016.csv") %>%
mutate(no_show = `No-show` == "Yes") # change the Yes/No into True/False (which is like having 1-0)## Parsed with column specification:
## cols(
## PatientId = col_double(),
## AppointmentID = col_integer(),
## Gender = col_character(),
## ScheduledDay = col_datetime(format = ""),
## AppointmentDay = col_datetime(format = ""),
## Age = col_integer(),
## Neighbourhood = col_character(),
## Scholarship = col_integer(),
## Hipertension = col_integer(),
## Diabetes = col_integer(),
## Alcoholism = col_integer(),
## Handcap = col_integer(),
## SMS_received = col_integer(),
## `No-show` = col_character()
## )## Warning: package 'bindrcpp' was built under R version 3.4.4# split to train and test set
appointments <- appointments %>%
mutate(is_train = runif(NROW(appointments)) <= 0.8)
# build the linear regression model
appointments_model <- glm(formula =
no_show ~ Gender + Age + Scholarship + Hipertension + Diabetes +
Alcoholism + Handcap + SMS_received,
family = binomial,
data = appointments %>% filter(is_train))
summary(appointments_model)##
## Call:
## glm(formula = no_show ~ Gender + Age + Scholarship + Hipertension +
## Diabetes + Alcoholism + Handcap + SMS_received, family = binomial,
## data = appointments %>% filter(is_train))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.9646 -0.6812 -0.6087 -0.5297 2.1198
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.3717824 0.0196477 -69.819 < 2e-16 ***
## GenderM -0.0236974 0.0181802 -1.303 0.192414
## Age -0.0065475 0.0004383 -14.939 < 2e-16 ***
## Scholarship 0.1862869 0.0273711 6.806 1e-11 ***
## Hipertension -0.1043517 0.0275337 -3.790 0.000151 ***
## Diabetes 0.0768353 0.0382977 2.006 0.044828 *
## Alcoholism 0.1371949 0.0500152 2.743 0.006087 **
## Handcap 0.0062107 0.0545697 0.114 0.909386
## SMS_received 0.6618315 0.0172679 38.327 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 89206 on 88333 degrees of freedom
## Residual deviance: 87336 on 88325 degrees of freedom
## AIC: 87354
##
## Number of Fisher Scoring iterations: 4מבנה התוצאות מאוד דומה. המקדמים כעת מייצגים את החזקות, כלומר כאשר:
Estimate of Scholarship = 0.17266 - כלומר בעלי מלגה יותר סביר שיבריזו מפגישה עם רופא, פי 1.2. או במילים אחרות, לגבי כל המקדמים, אפשר לגשת לחישוב המקדמים ברגרסיה באופן הבא:
exp(appointments_model$coefficients)## (Intercept) GenderM Age Scholarship Hipertension
## 0.2536545 0.9765812 0.9934739 1.2047679 0.9009084
## Diabetes Alcoholism Handcap SMS_received
## 1.0798642 1.1470517 1.0062300 1.9383392המדד R-Square התחלף במדד שנקרא AIC (קריטריון Akaike). כאשר המדד נמוך יותר המשמעות היא שהמודל טוב יותר. השאריות “residuals” התחלפו ב“deviance”. גם פה, ככל שהערך קטן יותר, המשמעות היא שהמודל יותר מתאים לנתונים. גם פה המקדמים השונים מוצגים עם מבחן השערה המדווח האם הם שונים מ-0 באופן מובהק, כלומר אכן משפיעים על משתנה המטרה.
תרגיל
נסו לשפר את מודל הרגרסיה הנ“ל על ידי הורדה ו/או הוספה של משתנים, כאשר אתם משתמשים בקריטריון deviance או AIC. באפשרותכם גם להיעזר בפונקציה
MASS::stepAICבאמצעות פונקציית
mutateהוסיפו ל-appointments משתנה המתאר את תוצאת החישוב של המודל המעודכן. הציגו בתרשים boxplot את ההתפלגות של ערכי תוצאת החישוב בהפרדה לפי המשתנה no_show. באפשרותכם להיעזר בקוד הבא:
# use the help of ?predict.glm to discover what do you need to put in the type argument
# (we want the probability of not showing to the appointment).
appointments <- appointments %>%
mutate(probability_no_show = predict(XXX,
newdata = appointments %>% filter(!is_train),
type = XXX))- ברגרסיה לוגיסטית מקבלים הסתברויות, אבל לרוב, אנחנו מעוניינים לסווג אותם באופן ברור יותר, כדי שניתן יהיה לקבל החלטה. סיווג אפשרי לדוגמה הוא להתייחס לכל מי שהסתברותו החזויה לאי-הגעה הינה מעל 0.2 כמישהו שהולך לא להגיע (0.2 הינו סף לדוגמה, אך באותה המידה זה יכל להיות סף אחר). הוסיפו משתנה ל-appointments הנקרא predicted_no_show0.2 שיקבל TRUE החל מערך 0.2 להסתברות החזויה.
appointments <- appointments %>%
mutate(XXX = XXX >= XXX)אחד מהכלים העומדים לרשותנו כדי להבין את טיב המודל, בנוסף ל-AIC ול-Deviance הוא “מטריצת בלבול”. מטריצת בלבול מציגה את כל הסיווגים הנכונים והלא נכונים מכל הסוגים. היא נראית כך:
|Val/Predict|FALSE|TRUE|
| FALSE | TN | FP |
| TRUE | FN | TP |כאשר:
- FP = False Positive = טעות מסוג ראשון
- FN = False Negative = טעות מסוג שני
או במילים אחרות: 
מקור: http://www.statisticssolutions.com/to-err-is-human-what-are-type-i-and-ii-errors/
- כעת נחשב מטריצת בלבול בהתבסס על וקטור הסיווג שחישבנו בסעיף הקודם. היעזרו בקוד הבא:
# This code will compute a confusion matrix using absolute numbers
confusion_abs <- appointments %>%
filter(!is_train) %>%
count(XXX, XXX) %>%
spread(predicted_no_show0.2, XXX)
# It is also useful to have it in percentage of row
confusion_prc <- appointments %>%
filter(!is_train) %>%
count(XXX, XXX) %>%
group_by(no_show) %>%
mutate(prop = n/sum(n)) %>%
select(-n) %>%
spread(predicted_no_show0.2, XXX)
מה שיעור הטעות מסוג ראשון? מה שיעור הטעות מסוג שני?
מה קורה כשמריצים את הקוד הנ“ל עם ערכים קיצוניים יותר (נניח 0.1 או 0.9)?
- חברת הביטוח הרפואי מעוניינת לבצע פעולות לעידוד הגעה. לשם כך היא מעוניינת לתמרץ את האנשים באמצעות שיחת טלפון אישית מאחות המרפאה יום לפני הביקור. החברה מעוניינת לקבוע סף תחזית שהחל ממנו תבצע את ההתקשרויות. כלכלני החברה חישבו שעלות התקשרות למטופל היא 5 שקלים, אך עלות אי-הגעה לפגישה עם רופא היא 25 שקלים. נתבקשתם להציע סף הסתברותי שהחל ממנו תבוצע השיחה למטופל. הניחו שאם מתבצעת שיחת הטלפון אז המטופל מגיע בסבירות של 100%. נבצע את התרגיל בשלבים.
# first, lets build a function that takes as an argument the threshold, and returns the total cost
compute_cost <- function(XXX){
appointments_cost <- appointments %>%
filter(!is_train) %>%
mutate(predicted_outcome = probability_no_show >= threshold) %>%
mutate(per_record_cost =
5*(predicted_outcome == 1 & no_show == 1) +
25*(XXX == 0 & XXX == XXX)
0*(XXX = XXX & XXX == XXX) +
5*(XXX = XXX & XXX = XXX))
return(sum(XXX))
}
# to get the function loaded into the environment just run the code
# (ctrl+L when standing on the first line of the function or last line } of the function.
# R will then run your code and "compile" the function.)
# Now for the second part, let's run a for loop in different probability threshold to look for the
# minimum cost.
# There are better ways than for loops, like functional programming, but that's a more advanced topic
# that we're not going to cover now.
# initialize a variable which will contain the cost
cost <- NULL
# now compute the cost as a function of varying thresholds
for (threshold in seq(from = XXX, to = XXX, by = XXX)){
cost <- c(cost,
compute_cost(XXX))
}
# now make a tibble containing the cost and the thresholds you used
cost_benefit <- tibble(threshold = seq(XXX),
cost = cost)
# plot the cost as a function of threshold. When should the nurse call?
ggplot(cost_benefit, aes(XXX)) +
geom_point() + geom_line()בקוד זה בדקנו את כל הdataset, למרות שעלינו להשתמש בפועל ב-test set. עדכנו את הקוד, האם המלצתכם משתנה?
חשבו, מה הייתם ממליצים אילו העלות של אי-הגעה לרופא היתה 50 ש“ח במקום 25 ש”ח? מה הייתם ממליצים אם היא היתה 5 ש“ח?
- מונח מאוד חשוב בבעיות סיווג הוא ROC. ראשי התיבות מייצגים Reciever Operating Characteristic וצורת הניתוח פותחה במלחמת העולם השנייה כדי לנתח את יכולות הצבא האמריקאי לזהות מטוסים יפנים. מאז משתמשים בצורת הניתוח הזו (למעשה בסוג הגרף הזה) כדי לנתח בעיות סיווג שמתבססות על סף סיווג, כמו ברגרסיה לוגיסטית. בסעיף זה נלמד על עקומות ROC ונראה איך ניתן לצייר אחת המתאימה למודל שלנו.
עקומות ROC מייצגות בציר y את הזיכוי לזהות אירוע (במקרה שלנו, לזהות אי-הגעה לפגישה), כנגד הסיכוי לטעות ולסווג אירוע בטעות (במקרה שלנו לסווג מישהו כ“לא מגיע” למרות שבעצם הוא כן יגיע).
בשפת “מטריצת הבלבול שלנו”, אנחנו נצייר בציר y את ה-True Positive כנגד False Positive. עבור כל סף שנקבע יהיו לנו שני מספרים כאלו, וכאשר נשנה את הסף נקבל זוג חדש.
כדי לצייר את ה-ROC, השתמשו בהנחיות הבאות:
# first arrange the appointments dataset by the predicted probability vector, in a descending order
# then, using mutate compute the cumulative sum of real no shows divided by the total number of no shows
# using another mutate, compute the cumulative sum of show-ups divided by the total number of show-ups
appointments_roc <- appointments %>%
arrange(desc(XXX)) %>%
mutate(tpr = cumsum(XXX)/sum(XXX)) %>%
mutate(fpr = cumsum(!XXX)/sum(!XXX))
# can you explain why this works? hint, think about what your are doing when you are sorting
# the probability column.
# plot the resulting line
ggplot(appointments_roc, aes(x = fpr, y = tpr)) +
geom_line() +
xlab("False positive rate (1 - Specificity)") +
ylab("True positive rate (Sensitivity)") +
ggtitle("An ROC for our medical appointment logistic regression model")
geom_abline(intercept = 0, slope = 1)
ggplot(appointments_roc, aes(x = fpr, y = tpr)) +
geom_line() +
xlab("False positive rate (1 - Specificity)") +
ylab("True positive rate (Sensitivity)") +
scale_x_continuous(labels = scales::percent) +
scale_y_continuous(labels = scales::percent) +
ggtitle("An ROC for our medical appointment logistic regression model") +
geom_abline(intercept = 0, slope = 1)
חשבו: מה המשמעות של הקו y=x? איזה “כלל הכרעה” מתלכד עם קו זה? מה היה קורא אילו היינו מקבלים עקום שנמצא כולו מתחת לקו y=x?
רגרסיה לוגיסטית לסיווג מספר קטגוריות
איך משתמשים ברגרסיה לוגיסטית כדי לתת סיווג ליותר מקטגוריה אחת? בעקרון מדובר בהכללה פשוטה. נניח ויש לנו שלושה סיווגים, setosa, versicolor, virginica. כל שעלינו לעשות הוא לבנות שני מודלים של רגרסיה לוגיסטית: \[P(Y=\text{setosa}|X)\] \[P(Y=\text{versicolor}|X)\] המודל השלישי קשור לשניים הקודמים: \[P(Y=\text{virginica}|X) = 1-P(Y=\text{setosa}|X)-P(Y=\text{versicolor}|X)\]
R יודע לעשות זאת בעצמו, לדוגמה:
iris_mult <- nnet::multinom(formula = Species ~ ., data = iris)## # weights: 18 (10 variable)
## initial value 164.791843
## iter 10 value 16.177348
## iter 20 value 7.111438
## iter 30 value 6.182999
## iter 40 value 5.984028
## iter 50 value 5.961278
## iter 60 value 5.954900
## iter 70 value 5.951851
## iter 80 value 5.950343
## iter 90 value 5.949904
## iter 100 value 5.949867
## final value 5.949867
## stopped after 100 iterationssummary(iris_mult)## Call:
## nnet::multinom(formula = Species ~ ., data = iris)
##
## Coefficients:
## (Intercept) Sepal.Length Sepal.Width Petal.Length Petal.Width
## versicolor 18.69037 -5.458424 -8.707401 14.24477 -3.097684
## virginica -23.83628 -7.923634 -15.370769 23.65978 15.135301
##
## Std. Errors:
## (Intercept) Sepal.Length Sepal.Width Petal.Length Petal.Width
## versicolor 34.97116 89.89215 157.0415 60.19170 45.48852
## virginica 35.76649 89.91153 157.1196 60.46753 45.93406
##
## Residual Deviance: 11.89973
## AIC: 31.89973אבל ישנן שיטות טובות יותר לצורך מציאת סיווגים של מספר קטגוריות, ובפרק הבא נתמקד בשתיים מהן: LDA ו-QDA.