מתאם ב-R: cor(), cor.test() ומטריצות מתאם

מתאם בין שני משתנים: cor()

מתאם שואל: כשמשתנה אחד עולה, האם השני נוטה לעלות גם הוא (חיובי), לרדת (שלילי), או לעשות מה שבא לו (קרוב לאפס)? ב-R זו קריאה אחת:

print(cor(mtcars$wt, mtcars$mpg))

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

איך קוראים את r

המקדם תמיד נופל בין −1 ל-+1. הסימן נותן את הכיוון; הגודל נותן את העוצמה:

| |r| | פירוש מקובל | | --- | --- | | 0.0 עד 0.2 | זניח | | 0.2 עד 0.4 | חלש | | 0.4 עד 0.6 | בינוני | | 0.6 עד 0.8 | חזק | | 0.8 עד 1.0 | חזק מאוד |

התייחסו לטווחים האלה כנקודת פתיחה לשיחה, לא כחוק: בפיזיקה r של 0.6 הוא אכזבה, בפסיכולוגיה הוא שיא של קריירה. שתי תכונות שכדאי להפנים: ל-r אין יחידות (מתאם בין משקל בטונות ל-mpg נותן את אותו r כמו משקל בקילוגרמים מול mpg, כי r מחושב על ערכים מתוקננים), ו-r מודד רק קשר לינארי: קשר מושלם בצורת U יכול להיות עם r ≈ 0.

והמשפט שחייבים לומר: מתאם אינו סיבתיות. מכירות גלידה ומקרי טביעה נמצאים במתאם חזק לאורך חודשי השנה, לא כי גלידה מטביעה אנשים, אלא כי הקיץ מניע את שניהם. מתאם אומר לכם ששני משתנים נעים יחד; הוא שותק לגבי הלמה. אולי x מניע את y, אולי y מניע את x, אולי משהו שלישי (העונה, במקרה הזה) מניע את שניהם. כדי להכריע ביניהם צריך ניסויים או חשיבה סיבתית זהירה, לא r גדול יותר.

Spearman ו-Kendall: כש-Pearson הוא הכלי הלא נכון

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

x <- c(1, 2, 3, 4, 5, 6, 7, 8)
y <- c(1, 4, 9, 16, 25, 36, 49, 300)  # increasing, with a wild last value

print(cor(x, y))                        # Pearson - dragged by the outlier
print(cor(x, y, method = "spearman"))   # ranks: perfectly monotonic
print(cor(x, y, method = "kendall"))

Spearman מחליף כל ערך בדירוג שלו ואז מחשב Pearson על הדירוגים. מכיוון ש-y תמיד עולה כש-x עולה, כל הדירוגים מסתדרים ו-Spearman מדווח בדיוק 1: הגודל של החריג מפסיק להיות חשוב, רק המיקום שלו. פנו ל-Spearman כשהנתונים סודרים (סולמות בסקרים), מוטים מאוד, או כשהקשר מונוטוני אבל לא ישר. Kendall עונה על שאלה דומה בעזרת זוגות תואמים ולא תואמים; הוא עמיד יותר במדגמים קטנים אבל איטי יותר, ו-Spearman הוא ברירת המחדל הנפוצה.

מטריצת המתאם

כדי לסרוק קשרים בין משתנים רבים בבת אחת, העבירו ל-cor() כמה עמודות מספריות:

m <- cor(mtcars[, c("mpg", "wt", "hp", "disp")])
print(round(m, 2))

כל משתנה מול כל משתנה אחר, עם 1 באלכסון (כל דבר נמצא במתאם מושלם עם עצמו) ותמונת מראה משני צדיו. עיגול לשתי ספרות אחרי הנקודה חשוב יותר ממה שנשמע: המטריצה הלא מעוגלת היא קיר של ספרות, וכל הטעם במטריצה הוא לסרוק אותה. כאן הסריקה מראה ש-mpg נמצא במתאם שלילי עם שלושתם (מכוניות כבדות, חזקות ועם מנוע גדול שורפות יותר דלק), ואילו wt, hp ו-disp כולם במתאם חיובי חזק זה עם זה: אשכול של משתני "מכונית גדולה" שיהיה חשוב כשתגיעו לרגרסיה לינארית ולכאבי הראש של המולטיקולינאריות.

ערכים חסרים: הארגומנט use

כשיש נתונים חסרים, ברירת המחדל של cor() היא להחזיר NA במקום לנחש:

d <- mtcars[, c("mpg", "wt", "hp")]
d$wt[c(3, 7)] <- NA  # knock two values out

print(cor(d)[1, 2])
print(cor(d, use = "complete.obs")[1, 2])
print(cor(d, use = "pairwise.complete.obs")[1, 2])

• use = "complete.obs" משמיט כל שורה שמכילה NA כלשהו, ואז מחשב את כל המטריצה מהשורות שנשארו. עקבי, אבל בזבזני: wt חסר מסיר את השורה הזו גם מהזוג mpg ו-hp.
• use = "pairwise.complete.obs" מחשב כל תא מכל השורות שבהן הזוג הזה קיים. הוא שומר יותר נתונים, אבל תאים שונים נשענים על תת-קבוצות שונות, ולעיתים נדירות מאוד זה יכול לייצר מטריצה שאינה עקבית פנימית.

לכמה NA בודדים כל אחד מהם בסדר; רק ציינו באיזה השתמשתם.

האם זה מובהק? cor.test()

cor() נותנת מספר אבל שום תחושה אם הוא יכול להיות רעש. cor.test() מריצה את מבחן ההשערות:

print(cor.test(mtcars$wt, mtcars$mpg))

עברו על הפלט חלק אחר חלק:

• t = −9.56, df = 30: סטטיסטי המבחן. השערת האפס היא שהמתאם האמיתי הוא אפס; ה-r שנצפה מומר לסטטיסטי t עם n − 2 דרגות חופש (32 מכוניות − 2).
• p-value = 1.29e-10: אם המתאם האמיתי היה אפס, ההסתברות לראות r רחוק כל כך מאפס במדגם של 32 היא בערך 0.0000000001. זו ראיה מכרעת שהקשר אמיתי, אבל זכרו: ערך ה-p מדבר על השאלה אם r שונה מאפס, לא על השאלה אם הקשר גדול או סיבתי.
• 95 percent confidence interval: −0.93 to −0.74: הטווח הסביר של המתאם האמיתי. לעיתים קרובות שימושי יותר מערך ה-p: אפילו הקצה האופטימי של הרווח הזה הוא מתאם שלילי חזק.
• sample estimates: cor = −0.87: אותו מספר ש-cor() נתנה לכם.

למדגמים קטנים מגיע כאן כבוד מיוחד: עם n = 10, מתאמים של ±0.5 מופיעים במקרה לעיתים קרובות באופן מדאיג, ורווח הסמך הרחב יגיד לכם את זה. דווחו על הרווח, לא רק על ערך ה-p.

לראות את זה: תמיד לשרטט

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

plot(mtcars$wt, mtcars$mpg)              # one pair
pairs(mtcars[, c("mpg", "wt", "hp", "disp")])  # every pair in the matrix

pairs() מצייר רשת של תרשימי פיזור שתואמת את מטריצת המתאם שלכם: הדרך המהירה ביותר לבדוק שהמספרים אומרים את מה שאתם חושבים. לגרפיקה מלוטשת של מטריצה בסגנון מפת חום, החבילה corrplot (install.packages("corrplot"), ואז corrplot(cor(m))) היא הכלי הסטנדרטי.

מה לקחת מכאן

• cor(x, y) נותנת את r של Pearson: הסימן הוא הכיוון, הגודל הוא העוצמה, תמיד בטווח [−1, 1], בלי יחידות.
• מתאם מודד תנועה משותפת לינארית ולא אומר דבר על סיבתיות: משתנה שלישי נסתר הוא תמיד מועמד.
• method = "spearman" לדירוגים: נתונים סודרים, חריגים, קשרים מונוטוניים אבל מעוקלים.
• cor(df) על עמודות מספריות נותנת את המטריצה; עגלו אותה עם round(, 2), ושימו לב לארגומנט use = כשיש נתונים חסרים.
• cor.test(x, y) מוסיפה את ערך ה-p ורווח סמך: דווחו על הרווח.
• תמיד הסתכלו על תרשים הפיזור לפני שאתם מאמינים למספר.

הבא בתור: כשהשאלה מתחדדת מ"האם הם נעים יחד?" ל"האם הממוצע של הקבוצה הזו שונה מהממוצע של זו?", מבחן t.