רווחי סמך ב-R: t.test(), confint() ו-prop.test()

מה באמת אומר "ביטחון של 95%"

רווח סמך הופך אומדן נקודתי ("ממוצע המדגם הוא 5.61") לטווח כן ("הממוצע האמיתי נמצא ככל הנראה בין 5.29 ל-5.93"). אבל הביטוי ביטחון של 95% הוא מהמפורשים ביותר באופן שגוי בסטטיסטיקה, אז בואו נבין אותו נכון כבר בהתחלה.

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

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

רווח סמך לממוצע, בדרך הקלה

כל t.test() נושא איתו רווח סמך; אפשר להריץ את המבחן רק כדי לקצור אותו:

times <- c(5.1, 4.9, 6.2, 5.8, 5.5, 5.9, 6.1, 5.2, 5.4, 6.0)

result <- t.test(times)
print(result$conf.int)
print(result$estimate)

ממוצע 5.61, רווח סמך של 95% בערך מ-5.29 עד 5.93. הרווח עושה את מה שממוצע חשוף לא יכול: הוא מראה כמה דיוק עשר תצפיות באמת קונות לכם.

אותו רווח סמך ידנית

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

times <- c(5.1, 4.9, 6.2, 5.8, 5.5, 5.9, 6.1, 5.2, 5.4, 6.0)

n <- length(times)
m <- mean(times)
se <- sd(times) / sqrt(n)          # standard error of the mean
t_crit <- qt(0.975, df = n - 1)    # t critical value, 2.5% in each tail

margin <- t_crit * se
print(t_crit)
print(margin)
print(c(m - margin, m + margin))

השוו לפלט של t.test(): אותו רווח עד הספרה האחרונה. שלושה חלקים נעים:

• שגיאת התקן sd(x)/sqrt(n) מודדת כמה ממוצע המדגם מתנדנד (ראו סטטיסטיקה תיאורית להבחנה בין סטיית תקן לשגיאת תקן).
• qt(0.975, df) הוא הערך הקריטי של t: לכיסוי של 95% משאירים 2.5% בכל זנב, ומכאן 0.975. עם df = 9 הוא בערך 2.26, רחב יותר מה-1.96 של ההתפלגות הנורמלית: המחיר של מדגם קטן על אמידת סטיית התקן מאותם נתונים.
• מרווח הטעות הוא המכפלה שלהם: המספר עם ה-"±" שמופיע בכותרות.

שינוי הרמה: 90%, 99% והפשרה

conf.level שולט בכיסוי, ואיתו ברוחב:

times <- c(5.1, 4.9, 6.2, 5.8, 5.5, 5.9, 6.1, 5.2, 5.4, 6.0)

print(t.test(times, conf.level = 0.90)$conf.int)
print(t.test(times, conf.level = 0.95)$conf.int)
print(t.test(times, conf.level = 0.99)$conf.int)

יותר ביטחון עולה יותר רוחב: רווח של 99% חייב להיות רחב מספיק כדי לצדוק 99 פעמים מתוך 100, ולכן הוא נמתח; רווח של 90% צר יותר אבל מפספס פי שניים יותר מ-95%. אין ארוחות חינם, יש רק כפתור כיוון. 95% הוא מוסכמה טהורה: ברירת מחדל ששווה לשמור עליה אלא אם יש לכם סיבה, לא חוק טבע.

רווח סמך לפרופורציה: prop.test()

נניח ש-47 מתוך 120 משתמשים שנסקרו מאמצים פיצ'ר חדש. מה הטווח הסביר לשיעור האימוץ האמיתי?

result <- prop.test(47, 120)
print(result$estimate)
print(result$conf.int)

פרופורציה במדגם 0.39, רווח סמך של 95% בערך מ-0.30 עד 0.49, כך ש"אימוץ של כ-40%" הוא הוגן רק עם ההילה של ±9 נקודות שמוצמדת אליו. prop.test() משתמשת בקירוב טוב יותר מנוסחת ספרי הלימוד p ± 1.96 × sqrt(p(1−p)/n) (זה רווח מסוג Wilson עם תיקון רציפות), וזה חשוב במיוחד ליד 0 או 1: רווח ספרי הלימוד יכול לחרוג מחוץ ל-[0, 1]; זה לא. לספירות קטנות מאוד, binom.test(47, 120)$conf.int נותנת את הגרסה המדויקת.

רווחי סמך למקדמי מודל: confint()

מודלים מותאמים מקבלים רווחים דרך פונקציה גנרית אחת:

fit <- lm(mpg ~ wt, data = mtcars)
print(confint(fit))
print(confint(fit, level = 0.99))

כל שורה תוחמת מקדם אחד של הרגרסיה הלינארית: רווח ה-95% של השיפוע של wt משתרע בערך בין −6.5 ובין −4.2 mpg לכל 1000 ליברות. זה אומר יותר מערך ה-p בסיכום: לא רק שההשפעה שונה מאפס, אלא שהיא לפחות כ-4 mpg ואולי כ-6.5.

ב-glm() לוגיסטי, confint(fit) עובדת באותה צורה אבל מחזירה גבולות ב-log-odds; העלו באקספוננט כדי לקבל רווחים ל-odds ratio, exp(confint(fit)), וזכרו שערך הייחוס של "אין השפעה" הופך ל-1 במקום 0.

גודל מדגם: חוק השורש הריבועי

הרוחב קטן לפי sqrt(n), ולזה יש תוצאה שקל לזכור: נתונים פי ארבעה, רווח בחצי. ראו את זה קורה עם נתונים מדומים: אותה התפלגות, מדגם אחד גדול פי 4 (ה-seed הופך את זה לשחזור; ראו מספרים אקראיים):

set.seed(42)
small <- rnorm(50, mean = 100, sd = 15)
large <- rnorm(200, mean = 100, sd = 15)   # 4 times the n

ci_small <- t.test(small)$conf.int
ci_large <- t.test(large)$conf.int

print(diff(ci_small))   # width with n = 50
print(diff(ci_large))   # width with n = 200 - about half

הרוחבים נוחתים קרוב ליחס של 2:1 (רעש הדגימה מונע ממנו להיות מדויק). השורש הריבועי הוא הסיבה שדיוק נעשה יקר: 100 התצפיות הראשונות קונות לכם יותר צמצום מ-300 הבאות יחד, וחציית רווח שכבר צר תמיד עולה פי 4 ממה ששילמתם עד עכשיו.

רווח סמך וערך p: שתי השקפות על אותו מבחן

רווח סמך ומבחן השערות הם אותו מידע בבגדים שונים. רווח הסמך של 95% מכיל בדיוק את ערכי הפרמטר שמבחן דו-צדדי ברמת 0.05 לא היה דוחה. לכן:

• רווח סמך להפרש ממוצעים לא כולל 0 ⇔ מבחן t אומר p < 0.05.
• רווח סמך ל-odds ratio לא כולל 1 ⇔ ה-p של המקדם קטן מ-0.05.

כשיש לכם את הרווח, בדרך כלל יש לכם את הסיכום הטוב יותר: הוא נותן את אותה הכרעה לגבי מובהקות, ובנוסף את גודל ההשפעה ביחידות אמיתיות. "p = 0.03" אומר שקיים הבדל; "רווח סמך של 95%: 0.2 עד 7.6" אומר שהוא קיים ושהוא עשוי להיות זניח או עצום, ולעיתים קרובות זה הממצא שחשוב.

מה לקחת מכאן

• 95% מתאר את שיעור הפגיעה של ההליך לאורך זמן, לא את ההסתברות שהרווח הספציפי הזה תפס את האמת.
• ממוצע: t.test(x)$conf.int, או ידנית כממוצע ± qt(0.975, n−1) × שגיאת תקן.
• פרופורציה: prop.test(x, n)$conf.int; מקדמי מודל: confint(fit) (העלו באקספוננט בשביל odds ratios של glm).
• ביטחון גבוה יותר = רווח רחב יותר; 95% הוא מוסכמה, לא חוק.
• הרוחב קטן לפי sqrt(n): נתונים פי 4 חוצים את הרווח.
• רווח הסמך מכיל כל ערך שמבחן ברמת 0.05 לא היה דוחה, ובניגוד לערך p, הוא מראה את גודל ההשפעה.

הבא בתור: ארגז הכלים לסימולציה שעומד מאחורי כל זה: מספרים אקראיים עם rnorm, runif, sample ו-set.seed.