Logistic Regression ב-R: glm() עם family = binomial

כשהתוצאה היא כן/לא

רגרסיה לינארית חוזה מספר. אבל רבות מהשאלות ששווה למדל הן בינאריות: האם הלקוח נוטש, האם המטופל מחלים, האם לוחצים על האימייל. התאמת קו ישר לתוצאה של 0/1 נשברת מיד: הקו חוזה בשמחה הסתברויות של −0.3 או 1.4, שהן חסרות משמעות.

רגרסיה לוגיסטית פותרת את זה על ידי מידול ההסתברות של התוצאה דרך טרנספורמציית ה-log-odds (logit): log(p / (1 − p)) = intercept + slope × x. סקאלת ה-log-odds משתרעת על כל ציר המספרים, כך שמשוואה לינארית מתאימה לה באופן טבעי, והמיפוי חזרה דוחס כל חיזוי לתוך (0, 1) לאורך עקומת ה-S המוכרת. המחיר של הטריק: המקדמים חיים בסקאלת ה-log-odds, וכל המשחק של קריאת רגרסיה לוגיסטית הוא לתרגם אותם חזרה למשהו שבני אדם מבינים.

התאמה: glm() עם family = binomial

glm() (generalized linear model, מודל לינארי מוכלל) הוא האח הגדול של lm(); family = binomial בוחר ברגרסיה לוגיסטית. ב-mtcars, am מתעד את סוג תיבת ההילוכים (1 = ידנית, 0 = אוטומטית). האם מכוניות חסכוניות בדלק נוטות להיות ידניות?

fit <- glm(am ~ mpg, data = mtcars, family = binomial)
summary(fit)

שני דברים לפני שקוראים את הפלט. ראשית, family = binomial אינו אופציונלי: השמיטו אותו ו-glm() מתאים בשקט ריבועים פחותים רגילים. שנית, התוצאה חייבת להיות בינארית: 0/1, לוגית, או factor עם שתי רמות (R ממדל את ההסתברות של הרמה השנייה).

עכשיו ה-summary, בלוק אחר בלוק:

• Coefficients: ה-Estimate של mpg הוא בערך 0.31, והוא שיפוע של log-odds: כל mpg נוסף מוסיף 0.31 ל-log-odds של תיבה ידנית. חיובי פירושו "מעלה את ההסתברות", שלילי פירושו "מוריד". מעבר לסימן, האינטואיציה של אף אחד לא עובדת בסקאלה הזו, ובגלל זה קיים הסעיף הבא.
• z value ו-Pr(>|z|): אותו היגיון כמו מבחני ה-t של הרגרסיה (Estimate ÷ Std. Error, ואז ערך p לשאלה "האם זה יכול להיות אפס?"), רק עם קירוב נורמלי, ולכן z במקום t. כאן p ≈ 0.011: לא סביר שהקשר בין mpg לסוג תיבת ההילוכים הוא רעש.
• Null deviance מול Residual deviance: ה-deviance הוא מדד חוסר ההתאמה בעולם של glm (קטן יותר = טוב יותר). ה-null deviance (43.2 על 31 df) הוא המודל עם חותך בלבד; ה-residual deviance (29.7 על 30 df) הוא המודל שלכם. הירידה של כ-13.6 בתמורה לדרגת חופש של משתנה מסביר אחד היא המקבילה ב-glm ל"R-squared עלה".
• AIC: ציון להשוואת מודלים שמאזן בין התאמה למורכבות; הנמוך מנצח. חסר משמעות לבדו, שימושי בין מודלים מועמדים על אותם נתונים.

מ-log-odds ליחסי סיכויים: exp(coef())

העלאה באקספוננט מעבירה מקדמים מסקאלת ה-log-odds החיבורית לסקאלת הסיכויים הכפלית:

fit <- glm(am ~ mpg, data = mtcars, family = binomial)

print(exp(coef(fit)))
print(exp(confint(fit)))

exp(0.307) ≈ 1.36, והנה תבנית המשפט הכנה שכדאי לשנן: "כל mpg נוסף מכפיל את הסיכויים לתיבת הילוכים ידנית בכ-1.36." יחס סיכויים מעל 1 מעלה את הסיכויים, מתחת ל-1 מוריד אותם, ובדיוק 1 פירושו שאין אפקט. זו הסיבה שהשאלה המכריעה של רווח הסמך ליחסי סיכויים היא "האם הרווח לא כולל את 1?" (לא אפס; אפס היה הגבול בסקאלת ה-log-odds).

שימו לב לשפה: סיכויים אינם הסתברויות. סיכויים = p / (1 − p), כך שהסתברות של 0.75 היא סיכויים של 3. הכפלת הסיכויים ב-1.36 אינה זהה להכפלת ההסתברות ב-1.36, וכשהתוצאה נפוצה, הפער גדול. יחס סיכויים של 2 לתוצאה נדירה מתנהג כמו "בערך פי שניים סיכון"; לתוצאה שקורית ב-50% מהמקרים, הוא בהחלט לא. לעולם אל תדווחו על יחס סיכויים בניסוח של יחס סיכונים ("סביר פי 1.36") אלא אם התוצאה נדירה.

הסתברויות חזויות: המלכודת של type = "response"

הבאג הנפוץ ביותר ברגרסיה לוגיסטית בעולם האמיתי:

fit <- glm(am ~ mpg, data = mtcars, family = binomial)
new_cars <- data.frame(mpg = c(15, 20, 25, 30))

print(predict(fit, newdata = new_cars))                      # log-odds!
print(predict(fit, newdata = new_cars, type = "response"))   # probabilities

הקריאה הראשונה מחזירה את ברירת המחדל type = "link": חיזויים בסקאלת ה-log-odds, כולל ערכים שליליים. השנייה מחזירה הסתברויות אמיתיות. אם ה"הסתברויות" שלכם יוצאות אי פעם שליליות או גדולות מ-1, זו הסיבה. הריצו את הבלוק: למכונית של 15 mpg אין כמעט שום סיכוי להיות ידנית, מכונית של 30 mpg ככל הנראה ידנית, ועקומת ה-S מתכופפת באמצע.

סיווג: סף וטבלת הבלבול

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

fit <- glm(am ~ mpg, data = mtcars, family = binomial)

probs <- predict(fit, type = "response")   # fitted probabilities, one per car
predicted <- probs > 0.5

print(table(actual = mtcars$am, predicted = predicted))

התאים שעל האלכסון הם החלטות נכונות; שני התאים שמחוץ לאלכסון הם שתי הטעויות השונות (חיזוי של ידנית למכונית אוטומטית, והפוך). דיוק כולל לבדו יכול להחמיא למודל מאוד: אם 95% מהלקוחות לא נוטשים, "לחזות שאף אחד לא נוטש" משיג 95% ותופס אפס נוטשים. לכן תמיד הסתכלו על שני סוגי השגיאות. ו-0.5 הוא מוסכמה, לא חוק: כשלשתי הטעויות יש מחירים שונים, הזיזו את הסף בהתאם.

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

כמה משתנים מסבירים

בדיוק כמו lm(): הוסיפו איברים עם +, וכל פירוש מקבל את ההסתייגות "כשהאחרים נשארים קבועים":

fit2 <- glm(am ~ mpg + qsec, data = mtcars, family = binomial)
summary(fit2)
print(exp(coef(fit2)))

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

אזהרות

• הפרדה מלאה (complete separation). אם משתנה מסביר מפצל את התוצאה באופן מושלם (כל מכונית מעל mpg מסוים ידנית, כל אחת מתחת אוטומטית), מקדם הנראות המרבית רוצה להיות אינסופי. R מזהיר, glm.fit: fitted probabilities numerically 0 or 1 occurred, ומדווח על מקדמים עצומים עם שגיאות תקן אבסורדיות. אל תפרסמו את המספרים האלה; פשטו את המודל, השיגו עוד נתונים, או השתמשו בשיטה עם קנס (החבילות brglm2 או logistf).
• מספיק אירועים. המגבלה המכרעת היא מספר המקרים של התוצאה הנדירה יותר, לא מספר השורות הכולל. כלל אצבע ישן דורש סדר גודל של 10 עד 15 אירועים לכל משתנה מסביר; הדוגמאות עם 32 המכוניות כאן נועדו ללמד את המכניקה, לא לשמש תבנית לגודל מדגם שאפשר לפרסם.
• יחסי סיכויים אינם יחסי סיכונים כשהתוצאה נפוצה. כוסה למעלה, וחוזר כאן כי סוקרים יתפסו את זה גם אם אתם לא.

מה לקחת מכאן

• תוצאה בינארית → glm(y ~ x, data = df, family = binomial); לעולם אל תשכחו את ה-family.
• המקדמים הגולמיים הם log-odds; exp(coef(fit)) נותן יחסי סיכויים, וערך האפס של הרווחים שלהם הוא 1.
• תבנית המשפט: "כל עלייה ביחידה אחת ב-x מכפילה את הסיכויים לתוצאה ב-exp(b)."
• predict(..., type = "response") להסתברויות: ברירת המחדל מחזירה log-odds, הבלבול מספר אחת.
• סווגו עם סף ושפטו עם טבלת בלבול; דיוק לבדו יכול לשקר.
• שימו לב לאזהרות הפרדה, ספרו את האירועים שלכם, ואל תלבישו יחסי סיכויים כיחסי סיכונים.

בהמשך: המנגנון שמאחורי כל רווח שראיתם עד עכשיו, רווחי סמך עם t.test(), confint() ו-prop.test().