Linear Regression ב-R: lm(), summary() ו-predict()

הרעיון: קו הריבועים הפחותים

רגרסיה לינארית מתאימה קו ישר דרך ענן של נקודות: y = intercept + slope × x. מכל הקווים האפשריים, lm() בוחר את זה שממזער את סכום ריבועי השאריות, כשהשארית היא הפער האנכי בין נקודה לקו. ההעלאה בריבוע גורמת לפספוסים גדולים להיחשב באופן לא פרופורציונלי, ולכן חריג פראי אחד יכול להטות את כל ההתאמה.

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

fit <- lm(mpg ~ wt, data = mtcars)
print(coef(fit))

קראו את הנוסחה כ"מדלו את mpg כפונקציה של wt". שני המקדמים הם הקו המותאם: mpg ≈ 37.3 − 5.3 × weight. לשיפוע יש יחידות אמיתיות: כל 1000 ליברות נוספות של מכונית (wt נמדד ביחידות של 1000 ליברות) עולות כ-5.3 מיילים לגלון. החותך (37.3 mpg במשקל אפס) הוא רק המקום שבו הקו חוצה את האפס; אף מכונית לא שוקלת כלום, אז אל תפרשו אותו יותר מדי.

מעבר על summary()

summary(fit) הוא הפלט שכל קורס סטטיסטיקה מבקש לפרש. הריצו אותו, ואז עברו עליו בלוק אחר בלוק:

fit <- lm(mpg ~ wt, data = mtcars)
summary(fit)

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

Residuals: סיכום חמשת המספרים של השאריות (actual − predicted). רוצים שהחציון יהיה קרוב ל-0 ושתהיה סימטריה גסה בין Min/Max ובין 1Q/3Q; אסימטריה חזקה רומזת שמודל הקו הישר מפספס משהו.

Coefficients: לב הפלט, שורה אחת לכל איבר:

• Estimate: הערך המותאם. עבור wt, −5.34: כל 1000 ליברות נוספות קשורות לכ-5.3 mpg פחות. תמיד תרגמו את השיפוע למשפט עם יחידות; המשפט הזה הוא כל התוכן המעשי של המודל.
• Std. Error: כמה האומדן היה מתנדנד בין מדגמים חוזרים. אומדנים שנמצאים בטווח של כמה שגיאות תקן מאפס רעועים.
• t value: Estimate ÷ Std. Error: כמה שגיאות תקן המקדם רחוק מאפס (−9.56 כאן).
• Pr(>|t|): ערך ה-p לשאלה "האם ייתכן שהמקדם הזה הוא באמת אפס?". עבור wt הוא כ-1.3e-10: אם למשקל באמת לא היה קשר לינארי ל-mpg, שיפוע תלול כזה כמעט אף פעם לא היה מופיע במדגם של 32. ערך p קטן = ראיה שהקשר קיים, לא הוכחה שהמודל נכון, ולא מדד לחשיבות (אפקט זעיר שנאמד בדיוק מקבל גם הוא ערך p זעיר).
• Signif. codes / כוכביות: קיצור חזותי לעמודת ערכי ה-p. נוח; לא מוסיף מידע.

Residual standard error: 3.05 on 30 degrees of freedom: הגודל הטיפוסי של פספוס בחיזוי, ביחידות של המשתנה התלוי עצמו: החיזויים סוטים בדרך כלל בכ-3 mpg. שפטו את זה מול הסקאלה של mpg (שנעה בערך בין 10 ל-34).

Multiple R-squared: 0.75: המשקל מסביר כ-75% מהשונות ב-mpg. Adjusted R-squared (0.74) מחשב את זה מחדש עם קנס על כל משתנה מסביר, כי הגרסה הגולמית יכולה רק לעלות כשמוסיפים משתנים, אפילו רעש אקראי. כשמשווים מודלים עם מספר שונה של משתנים מסבירים, הגרסה המתוקננת היא הכנה. והתנגדו לרפלקס ש"מודל טוב = R² גבוה": אפקט שימושי באמת יכול לחיות במודל עם R² נמוך (תוצאה רועשת, גורם אחד מתוך רבים), בעוד ש-R² גבוה יכול לנבוע מהתאמת יתר או ממשתנה שדלף.

F-statistic: 91.4 ... p-value: 1.29e-10: המבחן של המודל כולו: האם המודל הזה מנצח את "פשוט לחזות את הממוצע לכולם"? עם משתנה מסביר אחד הוא משכפל את מבחן ה-t של השיפוע (שימו לב ש-9.56² ≈ 91.4); עם כמה משתנים מסבירים הוא הופך למבחן המשותף לכך שלפחות מקדם אחד שונה מאפס. המנגנון שלו הוא אותו פירוק שונות של ANOVA.

רגרסיה מרובה: החזקת האחרים קבועים

הוסיפו משתנים מסבירים עם +:

fit2 <- lm(mpg ~ wt + hp, data = mtcars)
summary(fit2)

הפירוש משתנה בדרך מכרעת אחת. כל Estimate הוא עכשיו האפקט של אותו משתנה מסביר כשהאחרים נשארים קבועים: המקדם של wt (בערך −3.9, ירידה מ-−5.3) הוא המחיר ב-mpg של משקל נוסף בהשוואה בין מכוניות עם אותו כוח סוס. ה-−5.3 של הרגרסיה הפשוטה כלל בשקט את העובדה שמכוניות כבדות יותר נוטות גם להיות חזקות יותר; הרגרסיה המרובה מפרידה את זה. זו גם הסיבה שמקדמים זזים כשמוסיפים משתנים: אם המשתנה המסביר החדש מתואם עם ישן, הגדרת התפקיד של הישן משתנה. R-squared מטפס לכ-0.83, וכאן הגרסה המתוקננת היא ההשוואה ההוגנת מול המודל עם משתנה מסביר יחיד.

חיזויים: predict()

המודל המותאם הוא פונקציה; predict() מעריכה אותה. בנו data frame של newdata ששמות העמודות שלו תואמים בדיוק את המשתנים המסבירים:

fit <- lm(mpg ~ wt, data = mtcars)

new_cars <- data.frame(wt = c(2.5, 3.5))
print(predict(fit, newdata = new_cars))

print(predict(fit, newdata = new_cars, interval = "confidence"))
print(predict(fit, newdata = new_cars, interval = "prediction"))

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

• interval = "confidence": אי ודאות לגבי הממוצע: "עבור כל המכוניות ששוקלות 2500 ליברות, היכן נמצא ה-mpg הממוצע?". צר, ומצטמצם ככל שיש יותר נתונים.
• interval = "prediction": הטווח שבו מכונית חדשה בודדת במשקל הזה צפויה לנחות. רחב בהרבה, כי מכונית בודדת נושאת פיזור משלה סביב הקו, פיזור שאף כמות נתונים לא ממצעת אותו.

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

אבחון ומלכודת האקסטרפולציה

summary() אומר לכם מה המודל אומד; גרפי השאריות אומרים לכם אם להאמין לו. בסשן אינטראקטיבי:

par(mfrow = c(2, 2))
plot(fit)   # four diagnostic plots

למה לשים לב: Residuals vs Fitted צריך להיות ענן חסר צורה. עקומה פירושה שהקשר לא ישר; משפך (פיזור שגדל עם הערכים המותאמים) פירושו שונות לא קבועה, ושגיאות התקן שלכם שגויות. נקודות ה-Q-Q plot צריכות להיצמד לקו: זנבות כבדים פירושם שחריגים מעוותים את ההתאמה. Scale-Location הוא שוב בדיקת המשפך. Residuals vs Leverage מסמן נקודות משפיעות, תצפיות שלבדן גוררות את המקדמים (ב-mtcars, מכוניות אקזוטיות כמו Chrysler Imperial נוטות להופיע כאן). גרף פיזור מהיר של הנתונים הגולמיים לפני ההתאמה תופס את רוב זה מוקדם.

ולבסוף, המלכודת שאף אבחון לא תופס: אקסטרפולציה. המודל למד ממכוניות ששוקלות בערך 1500 עד 5400 ליברות. תנו ל-predict() ערך wt של 8 והוא יחזיר בעליזות mpg שלילי: המתמטיקה מאריכה את הקו לנצח, אבל הראיות נעצרות בקצה הנתונים. חזו רק בתוך הטווח שעליו התאמתם (או קרוב אליו).

מה לקחת מכאן

• fit <- lm(y ~ x, data = df) מתאים את קו הריבועים הפחותים; coef(fit) היא המשוואה, summary(fit) הדוח המלא.
• קראו את ה-Estimates כמשפטים עם יחידות; Pr(>|t|) שואל "האם זה יכול להיות אפס?", לא "האם זה חשוב?".
• Residual standard error הוא הפספוס הטיפוסי ביחידות אמיתיות; adjusted R-squared הוא המספר ההוגן להשוואת מודלים.
• ברגרסיה מרובה כל מקדם פירושו "כשהאחרים נשארים קבועים", ומקדמים זזים כשמצטרפים משתנים מסבירים מתואמים.
• predict(fit, newdata, interval = ...): "confidence" לממוצע, "prediction" למקרה חדש אחד, הרחב מביניהם.
• בדקו את plot(fit) לעקומות, משפכים ונקודות משפיעות; לעולם אל תסמכו על חיזויים מחוץ לטווח הנתונים.

בהמשך: כשהתוצאה היא כן/לא ולא מספר, רגרסיה לוגיסטית עם glm().