ANOVA ב-R: aov(), טבלת ANOVA ו-Tukey HSD

למה ANOVA ולא ערימה של מבחני t

מבחן t משווה שני ממוצעים. כשיש שלוש קבוצות או יותר, הפיתוי הוא להריץ מבחן t על כל זוג, וזו בדיוק המלכודת. כל מבחן ברמת 0.05 נושא סיכון של 5% לתוצאה חיובית שגויה, והסיכונים מצטברים: עם 4 קבוצות מדובר ב-6 מבחנים בין זוגות ובסיכוי של כ-26% לפחות לתוצאה "מובהקת" מדומה אחת; עם 5 קבוצות (10 מבחנים), כ-40%. אתם פשוט מייצרים תגליות מתוך רעש.

ANOVA (ANalysis Of VAriance, ניתוח שונות) פותר את זה בכך שהוא שואל שאלה אחת עם ערך p אחד: האם כל ממוצעי הקבוצות זהים, או שלפחות אחד מהם שונה? למרות השם, הוא משווה ממוצעים, רק שהוא עושה זאת בעזרת ניתוח של השונות: אם ממוצעי הקבוצות מפוזרים יותר ממה שהרעש בתוך הקבוצות יכול להסביר, קורה שם משהו אמיתי.

ANOVA חד-כיווני עם aov()

PlantGrowth בנוי בדיוק לזה: משקלי צמחים בקבוצת ביקורת ובשני תנאי טיפול. מתאימים עם aov(), ואז, וזה חשוב, מדפיסים את הטבלה עם summary():

fit <- aov(weight ~ group, data = PlantGrowth)
summary(fit)

קראו את הנוסחה כ"האם weight תלוי ב-group?". משתנה הקיבוץ חייב להיות factor. PlantGrowth$group כבר כזה, אבל אם הקבוצות שלכם מקודדות כמספרים (מינונים, מזהי אצווה), עטפו אותן: aov(y ~ factor(dose), ...). אחרת aov() מתאים בשקט קו רגרסיה דרך קודי הקבוצות במקום להשוות את ממוצעי הקבוצות: מודל שגוי, בלי שום הודעת שגיאה.

קריאת טבלת ה-ANOVA

לטבלה שתי שורות, ה-factor והשאריות, וחמש עמודות. תא אחר תא:

            Df Sum Sq Mean Sq F value Pr(>F)
group        2  3.766  1.8832   4.846 0.0159
Residuals   27 10.492  0.3886

• Df: דרגות חופש. שורת ה-factor מקבלת קבוצות − 1 (3 קבוצות → 2); שורת השאריות מקבלת תצפיות − קבוצות (30 − 3 = 27). בדיקה מהירה ש-R ראה את מבנה הניסוי שהתכוונתם אליו.
• Sum Sq: השונות, מחולקת לשתי ערימות. שורת group היא הערימה שבין הקבוצות: כמה רחוק ממוצעי הקבוצות יושבים מהממוצע הכללי. Residuals היא הערימה שבתוך הקבוצות: כמה הצמחים משתנים סביב הממוצע של הקבוצה שלהם. יחד הן מסתכמות לשונות הכוללת בנתונים.
• Mean Sq: כל Sum Sq מחולק ב-Df שלו, וכך ערימות השונות הופכות לשיעורים ברי השוואה לכל דרגת חופש. ה-Mean Sq של השאריות (0.389) הוא רמת הרעש.
• F value: היחס: Mean Sq של group חלקי Mean Sq של Residuals (1.8832 / 0.3886 ≈ 4.85). אם כל ממוצעי הקבוצות היו באמת שווים, היחס הזה היה מסתובב סביב 1. ככל ש-F גדל, ההבדלים בין הקבוצות עולים על מה שהרעש יכול להסביר.
• Pr(>F): ערך ה-p: ההסתברות לקבל F גדול כל כך אם שלושת הממוצעים האמיתיים היו זהים. כאן 0.016, קטן מספיק ברמה המקובלת של 0.05 כדי לדחות את "כולם שווים". כמו תמיד, זו לא ההסתברות שהשערת האפס נכונה, והוא לא אומר דבר על אילו קבוצות שונות או בכמה.

הנקודה האחרונה היא המגבלה המכרעת: F מובהק אומר "יש הבדל כלשהו, איפשהו". לא יותר.

אילו קבוצות שונות? TukeyHSD()

שאלת ההמשך דורשת מבחן בדיעבד (post-hoc). מבחן ה-Honest Significant Difference של Tukey בודק כל זוג ושומר על שיעור הטעות ברמת המשפחה על 5% בכל ההשוואות יחד:

fit <- aov(weight ~ group, data = PlantGrowth)
TukeyHSD(fit)

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

• diff: ההפרש המשוער בין הממוצעים (הקבוצה השנייה בשם פחות הראשונה).
• lwr, upr: רווח הסמך של 95% ברמת המשפחה עבור ההפרש הזה.
• p adj: ערך ה-p, כבר מתוקנן לכך שבוצעו שלוש השוואות.

עבור PlantGrowth: השורה trt2-trt1 מראה הפרש של כ-0.87 עם p adj ≈ 0.012 ורווח שלא נוגע באפס: טיפול 2 מצמיח יותר מטיפול 1. לשתי השורות של טיפול מול ביקורת (trt1-ctrl, trt2-ctrl) יש רווחים שחוצים את האפס ו-p adj הרבה מעל 0.05: אי אפשר להבחין בין אף אחד מהטיפולים לבין הביקורת במדגם הזה. כלל האצבע זהה לרווחי סמך בכל מקום: הרווח לא כולל את האפס ⇔ p adj מתחת ל-0.05.

שימו לב שה-F הכולל אמר "קיים הבדל", ואילו Tukey מאתר אותו בזוג אחד בדיוק. המבנה הדו-שלבי הוא כל הרעיון: מבחן כולל הוגן אחד, ואחריו עבודת בילוש בין זוגות עם תיקון נכון.

ANOVA דו-כיווני: שני factors והאינטראקציה ביניהם

עם שני משתני קיבוץ, קריאה אחת בודקת את שניהם, וגם אם יש ביניהם אינטראקציה. ToothGrowth מצליב את סוג התוסף (supp) עם המינון (0.5, 1, 2 מ"ג, מספרי, ולכן צריך factor()):

fit2 <- aov(len ~ supp * factor(dose), data = ToothGrowth)
summary(fit2)

supp * factor(dose) מתרחב לשלושה אפקטים, שורה אחת בטבלה לכל אחד:

• supp: בממוצע על פני המינונים, האם לסוג התוסף יש השפעה? (כן: p ≈ 0.0002.)
• factor(dose): בממוצע על פני התוספים, האם למינון יש השפעה? (בהחלט: ה-p זעיר.)
• supp:factor(dose): האינטראקציה: האם ההשפעה של התוסף משתנה בהתאם למינון? כאן p ≈ 0.022, כלומר כן. מיץ תפוזים עדיף על חומצה אסקורבית במינונים נמוכים, אבל הפער נסגר במינון של 2 מ"ג.

אינטראקציה מובהקת היא תווית אזהרה על האפקטים הראשיים: "לסוג התוסף יש השפעה" נכון רק בממוצע, והממוצע מסתיר סיפור שתלוי במינון. כשהאינטראקציה מובהקת, תארו את הצירופים (בעזרת TukeyHSD(fit2, "supp:factor(dose)") או טבלה של ממוצעי קבוצות דרך סיכומים לפי קבוצה) במקום לדווח רק על האפקטים הראשיים. השתמשו ב-+ במקום * רק כשאתם רוצים במכוון מודל בלי אינטראקציה.

ההנחות, וחלופה מבוססת דירוגים

המתמטיקה של ANOVA נשענת על שלוש הנחות, מהקשיחה לגמישה:

• אי-תלות: התצפיות לא משפיעות זו על זו. שום דבר לא מתקן הפרה בדיעבד; זו תכונה של מבנה המחקר.
• שונויות שוות בין הקבוצות: ה-Mean Sq של השאריות הוא אומדן רעש משותף אחד, ולכן הקבוצות צריכות להיות רועשות במידה דומה. בדיקה בשורה אחת: bartlett.test(weight ~ group, data = PlantGrowth) (ערך p גדול פירושו שאין ראיה לשונויות לא שוות).
• שאריות בערך נורמליות: חשוב הכי פחות כשהמבנה מאוזן והקבוצות בגודל סביר; הסתכלו על היסטוגרמה או boxplot של השאריות דרך residuals(fit).

כשהנתונים מוטים מאוד, סודרים או מלאים בחריגים, החלופה מבוססת הדירוגים ל-ANOVA חד-כיווני היא kruskal.test(weight ~ group, data = PlantGrowth), האח הרב-קבוצתי של wilcox.test(). זו שורה אחת, והיא מוותרת על מעט עוצמה תמורת עמידות.

מה לקחת מכאן

• ANOVA שואל שאלה אחת על ממוצעים של 3 קבוצות או יותר עם ערך p אחד: הפתרון לבעיית ריבוי מבחני t.
• fit <- aov(y ~ group, data = df); summary(fit), ומשתנה הקיבוץ חייב להיות factor.
• ערך F הוא האות שבין הקבוצות חלקי הרעש שבתוכן; Pr(>F) קטן פירושו "יש הבדל כלשהו, איפשהו", ולא יותר.
• TukeyHSD(fit) מוצא אילו זוגות שונים, עם תיקון להשוואות מרובות מובנה.
• a * b מתאים שני factors וגם את האינטראקציה ביניהם; אינטראקציה מובהקת פירושה שהאפקטים הראשיים לבדם לא מספרים את כל הסיפור.
• בדקו שוויון שונויות (bartlett.test) ונורמליות של השאריות; kruskal.test() היא החלופה מבוססת הדירוגים.

הבא בתור: מהשוואת ממוצעי קבוצות למידול של קשר, רגרסיה לינארית עם lm().