ממוצע, חציון וסטיית תקן ב-R

ממוצע וחציון

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

mpg <- mtcars$mpg

print(mean(mpg))
print(median(mpg))

mean() מחברת הכול ומחלקת במספר הערכים. median() ממיינת את הערכים ובוחרת את האמצעי (או ממצעת את שני האמצעיים כשמספר הערכים זוגי). עבור 32 המכוניות ב-mtcars, צריכת הדלק הממוצעת היא בערך 20.1 mpg והחציון 19.2: קרובים זה לזה, וזה אומר לכם שהנתונים לא מוטים במיוחד. כשהם חלוקים הרבה, גם זה מידע; נגיע לזה בסוף.

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

x <- c(4, 8, NA, 6)

print(mean(x))
print(mean(x, na.rm = TRUE))

זה מכוון: R מסרבת להעמיד פנים בשקט שהערך החסר לא קיים. na.rm = TRUE אומר "חשבי על הערכים שיש לך". כמעט כל פונקציית סיכום במאמר הזה מקבלת אותו. הסיפור המלא של איך NA מתפשט נמצא בערכים חסרים.

סטיית תקן ושונות

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

mpg <- mtcars$mpg

print(sd(mpg))
print(var(mpg))
print(sd(mpg)^2)

סטיית תקן של כ-6 mpg פירושה שמכוניות נמצאות בדרך כלל בטווח של כ-6 mpg מהממוצע של 20.1. מכיוון ש-sd() נמדדת ביחידות של הנתונים עצמם, היא זו שמדווחים עליה; var() מופיעה בעיקר בתוך נוסחאות אחרות.

הנה הפרט שחשוב בקורסים: sd() ו-var() מחשבות את הסטטיסטי של המדגם: הן מחלקות את סכום ריבועי הסטיות ב-n − 1, לא ב-n:

x <- c(2, 4, 4, 4, 5, 5, 7, 9)
n <- length(x)

sample_var <- sum((x - mean(x))^2) / (n - 1)
population_var <- sum((x - mean(x))^2) / n

print(sample_var)      # what var(x) gives you
print(var(x))
print(population_var)  # divide by n instead

למה n − 1? כי אמדתם את הממוצע מאותם נתונים, והסטיות סביב הממוצע המשוער הזה קטנות מעט מדי באופן שיטתי; חלוקה ב-n − 1 מתקנת את זה. מכיוון שהנתונים שלכם הם כמעט תמיד מדגם מתוך משהו גדול יותר, גרסת ה-n − 1 היא מה שאתם רוצים. אם באמת יש לכם את האוכלוסייה כולה (כל התלמידים בכיתה, כל המוצרים בקטלוג), הכפילו: var(x) * (n - 1) / n.

שגיאת התקן של הממוצע

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

mpg <- mtcars$mpg

se <- sd(mpg) / sqrt(length(mpg))
print(se)

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

summary(): סקירה בקריאה אחת

summary() נותנת לכם את סיכום חמשת המספרים ועוד הממוצע בבת אחת, והיא עובדת על data frames שלמים:

print(summary(mtcars$mpg))

print(summary(mtcars[, c("mpg", "wt", "hp")]))

כשקוראים לה על data frame, היא מסכמת כל עמודה: עמודות מספריות מקבלות מינימום, רבעונים, ממוצע ומקסימום, ו-factors מקבלים ספירה לכל רמה. זה הדבר הראשון שמריצים על כל מערך נתונים שזה עתה נטען: ערכים בלתי אפשריים (גיל שלילי, מקסימום של 9999) קופצים מיד לעין.

אחוזונים, טווח ו-IQR

quantile() מכלילה את החציון לכל נקודת חיתוך:

mpg <- mtcars$mpg

print(quantile(mpg))
print(quantile(mpg, probs = c(0.1, 0.5, 0.9)))

print(range(mpg))
print(IQR(mpg))
print(max(mpg) - min(mpg))

בלי ארגומנטים היא מחזירה את המינימום, הרבעונים והמקסימום. העבירו probs = לנקודות חיתוך ספציפיות: האחוזונים ה-10 וה-90 למעלה תוחמים את המקום שבו נמצא רוב הנתונים. IQR() (המרחק בין האחוזון ה-25 ל-75) הוא מדד פיזור שבניגוד ל-sd() לא נגרר על ידי חריגים. range() מחזירה את המינימום והמקסימום כזוג.

השכיח: mode() של R לא עושה את זה

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

x <- c(3, 7, 7, 2, 7, 5, 3)

print(mode(x))   # "numeric" - not what you wanted

counts <- table(x)
print(counts)
print(names(which.max(counts)))

הביטוי שכדאי לזכור: table(x) סופרת כמה פעמים כל ערך מופיע, which.max() מוצאת את הספירה הגדולה ביותר, ו-names() שולפת את הערך עצמו. שימו לב שהוא חוזר כמחרוזת תווים (שמות בטבלה תמיד כאלה); עטפו אותו ב-as.numeric() אם אתם צריכים לחשב איתו. אם שני ערכים שווים בספירה, which.max() מחזירה בשקט רק את הראשון: בדקו את הטבלה בעצמכם כשתיקו סביר.

ממוצע מול חציון: על מה לדווח

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

salaries <- c(42000, 45000, 47000, 48000, 51000, 54000, 58000)
print(mean(salaries))
print(median(salaries))

# The CEO joins the dataset:
salaries <- c(salaries, 900000)
print(mean(salaries))
print(median(salaries))

ערך אחד שנוסף דוחף את הממוצע מכ-49,300 ליותר מ-155,000, מספר שלא מתאר אף אחד בנתונים, ואילו החציון זז רק מ-48,000 ל-49,500. זו הסיבה שמדווחים על הכנסות, מחירי דירות ומשכי אשפוז כחציונים: נתונים מוטים עם זנב ארוך הופכים את הממוצע למטעה. בנתונים סימטריים בערך השניים מסכימים והממוצע בסדר (ויעיל יותר סטטיסטית). היסטוגרמה מהירה אומרת לכם באיזה מצב אתם נמצאים, והשוואת הממוצע לחציון היא בעצמה בדיקת הטיה בשורה אחת.

מה לקחת מכאן

• mean(x) ו-median(x): הוסיפו na.rm = TRUE כשיש ערכים חסרים.
• sd(x) ו-var(x) מחשבות את הסטטיסטי של המדגם (המכנה n − 1), וזה מה שאתם רוצים.
• שגיאת התקן היא sd(x) / sqrt(length(x)): היא מודדת כמה טוב אתם מכירים את הממוצע, לא כמה הנתונים מפוזרים.
• summary() על data frame טרי היא בדיקת השפיות המהירה ביותר ב-R.
• השכיח הוא names(which.max(table(x))): mode() של R עוסקת בסוגי אחסון.
• נתונים מוטים או חריגים: דווחו על החציון. נתונים סימטריים: הממוצע בסדר.

הבא בתור: מדידה של כמה שני משתנים נעים יחד: מתאם עם cor() ו-cor.test().