06-hypothesis-tests.qmd

# מבחני השערות {#sec-hypothesis-tests}

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

מטרת הפרק היא לא לספק את התיאוריה. לצורך השלמות תיאורתיות, אני ממליץ לקרוא את [@walpole1993probability].

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

-   מבחני השערות על תוחלת

-   מבחני השערות על שונות

-   מבחני השערות על אי-תלות

-   מבחני השערות על התפלגות

## מבחן השערות על תוחלת

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

```{r setup}
#| warning: false
#| message: false
#| include: false

ggplot2::theme_set(ggplot2::theme_bw())
```

נדגים זאת בעזרת מדגם הפינגויינים שראינו בפרק הקודם:

```{r penguin means}
#| warning: false
#| message: false
#| include: false
library(palmerpenguins)
library(tidyverse)
penguins %>% 
  filter(!is.na(bill_length_mm)) %>% 
  group_by(species) %>% 
  summarize(avg_length = mean(bill_length_mm),
            sample_size = n())
```

לכל אחד מזני הפינגויינים ממוצע שונה לאורך המקור. לדוגמה, פינגויינים מסוג Adelie הם בעלי מקור באורך ממוצע 38.8 מ"מ. בפועל מדובר בממוצע המחושב על מדגם של 151 תצפיות, וסביר להניח שהתוחלת בכלל אוכלוסית הפינגויים מזן Adelie היא קצת שונה. לדוגמה, תוחלת של 38 מ"מ, יכולה להביא לכך שממוצע על פני דגימה של 151 תצפיות יהיה 38.8 מ"מ.

המטרה של מבחן השערות על התוחלת היא לעזור לנו לקבוע ברמת שגיאה של $\alpha$ (אשר בדרך כלל נקבעת כ-5% שגיאה), *האם התוחלת היא בעלת ערך מסוים* $\mu_0$, *בהינתן תוצאות המדגם* שלנו.

צורת הסתכלות אחרת היא בניית רווח בר סמך לתוחלת, או במילים אחרות, *מה הטווח הסביר שבו נמצאת התוחלת בהינתן תוצאות המדגם*.

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

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

על ההיסטוגרמה הלבשנו שני קווים אנכיים: הממוצע בפועל בירוק (38.3913), והערך שאותו נבחן כתוחלת (38) באמצעות מבחן השערות.

```{r adelie bill length histogram}
adelie <- penguins %>% 
  filter(species == "Adelie") %>% 
  filter(!is.na(bill_length_mm))

ggplot(adelie, aes(x = bill_length_mm)) + 
  geom_histogram(bins = 15) + 
  geom_vline(xintercept = 38, color = "#E41A1C", linewidth = 1, linetype = 2) + 
  geom_vline(xintercept = mean(adelie$bill_length_mm), color = "#4DAF4A", linewidth = 1) + 
  ggtitle("Distribution of Adelie bill length")
```

נשתמש בפקודת `t.test` על מנת לייצר רווח בר סמך (או מבחן השערות) לתוחלת אורך המקור של פינגויינים מזן Adelie.

נבנה את רווח בר הסמך ברמת ביטחון של 95% (כלומר רמת הוודאות של הדיווח שנקבל היא 95%, או 5% שגיאה). כמו כן, נבחן האם ייתכן שהתוחלת הינה 38, למרות שהתצפיות במדגם בעלות ממוצע של 38.8.

```{r adelie ttest}
adeliettest <- t.test(adelie$bill_length_mm, 
       alternative = "two.sided",
       mu = 38, conf.level = 0.95)

adeliettest
```

**פרשנות תוצאות המבחן: אנו דוחים את השערת ה-0 שהתוחלת הינה 38.** להלן הסבר מפורט לפלט:

1.  בראש הפלט מופיע סוג המבחן שהופעל (One Sample t-test), כלומר מבחן t עבור מדגם יחיד (יש גם מבחנים עבור שני מדגמים מזווגים, ושני מדגמים לא-מזווגים, מיד נדגים).
2.  בשורה השניה מופיע הדאטה שבו השתמשנו.
3.  בשורה השלישית ההשערה האלטרנטיבית שנבחנה (המכונה $H_1$ . במקרה זה הגדרנו בארגומנט `alternative` מבחן השערה דו-צדדי, כלומר: $H_1: \mu\neq38$.
4.  בשורה הרביעית ערך $T$ של המבחן (המכונה גם סטטיסטי המבחן), במקרה זה ערכו הוא 3.6513, ובנוסף מספר דרגות החופש (df), וערך המובהקות (p-value). ערך p-value במקרה זה הינו מתחת ל-0.05 מה שמעיד על כך שעלינו לדחות את השערת ה-0 ($H_0: \mu=38$).
5.  בשורה הבאה (ובזו שאחריה) מופיע הרווח בר סמך לתוחלת. במקרה זה הרווח בר סמך ברמת בטחון של 95% מעיד על כך שתוחלת אורך המקור של הפינגויין מהזן הנבחן נמצאת בין 38.36 לבין 39.22. ניתן לשים לב שהתוחלת אותה בחנו בהשערת ה-0 (\$\\mu=38\$) נמצאת מחוץ לטווח זה, וזה עולה בקנה אחד עם הדחיה שתיארנו בסעיף קודם בהתייחס לערך ה-p-value.
6.  לבסוף מופיע הממוצע עצמו של המשתנה (38.7914).

::: callout-tip
# מעבר ממבחן דו-צדדי למבחן חד-צדדי

היעזרו בתיעוד הפקודה `t-test` ושנו את הארגומנט `alternative` על מנת לבצע את המבחן הסטטיסטי:

$$
H_0: \mu=38
$$

$$
H_1: \mu>38
$$

מה תוצאת המבחן? מה המסקנה?
:::

## שני מדגמים לא-מזווגים

כעת נפנה למבחן t-test שבו נשתמש כאשר ברצוננו להשוות בין שני מדגמים, לא מזווגים. לדוגמה על מנת לבחן האם ניתן לדחות את ההשערה שאורך המקור של פינגויינים מזן Chinstrap ואורך המקור של פינגויינים מזן Gentoo שווים (כפי שראינו ערכיהם דומים יחסית).

ראשית, נציג את ההתפלגות של שני הזנים בתרשים הבא:

```{r distribution of chinstrap and gentoo bill length}
chin_gentoo <- penguins %>% 
  filter(!is.na(bill_length_mm)) %>% 
  filter(species != "Adelie")

ggplot(chin_gentoo,
       aes(x = bill_length_mm, color = species)) + 
  geom_density() + 
  ggtitle("Distribution of Chinstrap and Gentoo bill length")
```

התרשים ממחיש לנו הבדלים בין התפלגות אורך המקור של שני הזנים (ובהמשך הפרק נראה איך בוחנים הבדלים בין שתי התפלגויות), אך כעת נרצה לבחון את ההבדלים בין ממוצעי אורך המקור של כל זן. נשתמש בפקודה `t-test` אך בצורת הפעלה מעט שונה מזו שבה השתמשנו במקרה של מדגם יחיד.

```{r two sample ttest}
chin_gentoo_ttest <- t.test(formula = bill_length_mm ~ species,
                            data = chin_gentoo)
chin_gentoo_ttest
```

צורת ההפעלה שבה השתמשנו מגדירה נוסחה ודאטה. הנוסחה היא `bill_length_mm ~ species`, כלומר בחן את אורך המקור לפי המשתנה species (במקרה זה, חייבים שלמשתנה המסביר `species` יהיו בדיוק שתי רמות, אחרת נקבל שגיאה).

**קיבלנו דחיה של השערת ה-0, כלומר ניתן לדחות את ההשערה שתוחלת אורך המקור בין שני הזנים שווה.** ברמת מובהקות של 0.00773.

הרווח בר סמך מתייחס להפרש בין התוחלות, קרי, ברמת ביטחון של 95% ההפרש בין התוחלות נע בין 0.357 לבין 2.301 (לטובת Chinstrap שבעלי אורך מקור גדול יותר).

הפקודה שהרצנו משתמשת במבחן המניח שונות שונה, ניתן לשנות זאת על ידי הוספת הארגומנט `var.equal=TRUE`.

```{r two sample ttest non equal variance}
chin_gentoo_nonequal_var_ttest <- t.test(formula = bill_length_mm ~ species,
                                         data = chin_gentoo,
                                         var.equal = TRUE)
chin_gentoo_nonequal_var_ttest
```

## שני מדגמים מזווגים

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

נמחיש זאת על ידי השוואה בין אורך המקור לעומק המקור (`bill_length_mm`, `bill_depth_mm`). מדובר בשני מדגמים, אך הם קשורים ביניהם (משום שכל תצפית של אורך קשורה לתצפית מסוימת של עומק).

ראשית נשתמש בתרשים על מנת להמחיש את ההתפלגויות של המשתנים:

```{r length and depth distribution}
penguins %>% 
  mutate(specimen = seq_along(species)) %>% 
  select(specimen, bill_length_mm, bill_depth_mm) %>% 
  pivot_longer(cols = -specimen, names_to = "measurement_type", values_to = "mm") %>% 
  ggplot(aes(x = mm, color = measurement_type)) + 
  geom_density()
```

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

נשתמש שוב בפקודת `t.test` עם שינוי נוסף.

```{r paired sample ttest}
paired_ttest <- t.test(x = penguins$bill_length_mm,
                       y = penguins$bill_depth_mm, 
                       paired = TRUE)
paired_ttest
```

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

ניתן לראות שההבדל בין אורך לעומק המקור עומד על 26.1 ועד 27.4 מ"מ, ברווח בר סמך של 95%.

::: callout-tip
היעזרו בפקודה `t.test` על מנת לבחון את ההשערה שההבדל בין אורך לעומק המקור עומד על 26.5 מ"מ.

לשם כך השתמשו במבחן מזווג עם הארגומנט `mu` (מה צריך להיות ערכו של הארגומנט?).
:::

## מבחן לשיוויון שונויות

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

המבחן שבו נשתמש בוחן את היחס בין השונויות, וכן הוא מניח ששתי השונויות הן של משתנים מקריים הלקוחים מהתפלגות נורמלית. המבחן נקרא מבחן F. נדגים את המבחן על השונות של אורך המקור של Gentoo לעומת אורך המקור של Chinstrap.

```{r var test gentoo chinstrap}
chin_gentoo_vartest <- var.test(formula = bill_length_mm ~ species,
                                data = chin_gentoo)
chin_gentoo_vartest
```

שימו לב שאופן השימוש בפקודה מאוד דומה לאופן השימוש בפקודה `t.test` (ראו בחלק של מדגמים לא מזווגים). להלן הסבר על פלט הפקודה:

-   השורה הראשונה של הפלט מציינת שמדובר במבחן F.

-   בשורה השלישית אנחנו רואים את נתוני המבחן (הסטטיסטי, דרגות חופש במונה ובמכנה, וכן את ה-p-value המחושב) כפי שניתן לראות p-value \> 0.05, ולכן לא ניתן לדחות את השערת האפס שהשונויות שוות בין שני המדגמים.

-   בשורה הבאה ניתן לראות מה בחן המבחן (ההשערה האלטרנטיבית היא שהיחס בין השונויות שונה מ-1).

-   לאחר מכן מופיע רווח בר סמך ליחס בין השונויות $(0.778, 1.819)$. ניתן לראות רווח בר הסמך כולל את הערך 1 (שקול לאי-דחיית השערת ה-0).

-   בשורה האחרונה הפקודה מספקת לנו אמד ליחס בין השונויות (1.17).

## מבחן חי-בריבוע (טיב התאמה, אי-תלות)

מבחן חי-בריבוע (\$\\Chi\^2\$) יכול לשמש אותנו לשתי מטרות מרכזיות (שתיהן בוחנות את השונות):

1.  מבחן טיב התאמה
2.  מבחן אי-תלות

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

### מבחן טיב התאמה

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

נחזור לדוגמה על זן ה-Adelie ונבחן את ההשערה שהתפלגות אורך המקור לקוחה מהתפלגות נורמלית.

נלביש את ההתפלגות הנצפית וההתפלגות התיאורתית (עם ממוצע וסטיית תקן לפי הממוצע וסטיית התקן הנאמדים מהנתונים).

```{r distribution adelie bill length}

length_mu <- mean(adelie$bill_length_mm)
length_sd <- sd(adelie$bill_length_mm)

theoretical_bill_length <- tibble(x = seq(32, 46, 0.1)) %>% 
  mutate(length_dens = dnorm(x = x,
                             mean = length_mu,
                             sd = length_sd))

ggplot(adelie, aes(x = bill_length_mm)) + 
  geom_density(color = "blue") +
  geom_line(data = theoretical_bill_length, 
            inherit.aes = FALSE,
            aes(x = x, y = length_dens), color = "red") +
  ggtitle("Density of Adelie bill length with theoretical density")
```

השלבים שבוצעו בקוד כללו ראשית בניית טבלה של צפיפות של ההתפלגות הנורמלית עם ממוצע `r round(mean(adelie$bill_length_mm), 2)`, וסטיית תקן `r round(sd(adelie$bill_length_mm), 2)`, הם האמדים מהנתונים האמפיריים של תוחלת וסטיית תקן אורך המקור של פינגויינים מזן Adelie. הפונקציה `dnorm` משמשת לחישוב הצפיפות הנורמלית.

לאחר מכן השתמשנו ב`ggplot`, שאליו הוספנו שכבה נוספת של `geom_line` אשר לא משתמשת בדאטה הראשי (`adelie`) אלא בטבלה המשנית שחישבנו בשלב קודם (`theoretical_bill_length`). השתמשנו בארגומנטים `data = theoretical_bill_length` וב-`inherit.aes = FALSE` , על מנת להנחות את הפקודה של `geom_line` להשתמש בדאטה שונה, ובהגדרות אסטטיקה שונות.

התוצאה נראית כמו שני גרפי צפיפויות דומים אחד לשני, עם חריגה מסוימת (באדום הצפיפות התיאורתית ובכחול האמפירית). דומים, אבל לא זהים. האם הם "שונים מספיק" על מנת שנדחה את השערת האפס (שהם אותה ההתפלגות)?

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

#### תרשים Q-Q-Plot

בתרשים מסוג Q-Q-Plot המחשב מציב אחוזונים של שתי התפלגויות כנקודות במרחב (התפלגות אחת בציר X והשניה בציר Y). אם הנקודות המוצגות מתלבשות על אלכסון הדבר מעיד על כך שהן מגיעות מאותן התפלגות. ככל שהנקודות מתרחקות מהאלכסון, הדבר מעיד על כך שההתפלגויות אינן זהות. התרשים מהווה כלי ויזואלי שיכול לתמוך את התיאור של מבחן טיב ההתאמה.

```{r qqplot demonstration}
adelie %>% 
  ggplot(aes(sample = bill_length_mm)) + 
  geom_qq(distribution = stats::qnorm) + 
  stat_qq_line()
```

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

#### ביצוע מבחן טיב ההתאמה

במבחן טיב ההתאמה עלינו לחלק את ההתפלגות האמפירית לתאים סטטיסטיים (שבדרך כלל מקובל לסמן ב-\$O_i\$ מלשון Observed), ולהשוות למול תאים בהתפלגות התיאורתית (שמקובל לסמן באות $E_i$ מלשון Expected). אחד התנאים לשימוש במבחן הוא שבכל תא סטטיסטי מספר התצפיות הצפוי יהיה לפחות 5 תצפיות.

```{r showing goodness of fit}

distribution_breaks <- c(0, 35, 37, 39, 41, 50)

observed_expected_counts <- adelie %>% 
  mutate(length_groups = cut(bill_length_mm, 
                             breaks = distribution_breaks)) %>% 
  count(length_groups, name = "observed") %>% 
  mutate(probability_upper = 
           pnorm(q = distribution_breaks[-1], mean = length_mu, sd = length_sd),
         probability_lower = 
           pnorm(q = distribution_breaks[-6], mean = length_mu, sd = length_sd)
  ) %>% 
  mutate(expected_prob = (probability_upper - probability_lower)) %>% 
  mutate(expected_counts = expected_prob * sum(observed))

observed_expected_counts

goodness_of_fit_chisqtest <- chisq.test(x = observed_expected_counts$observed,
                                        p = observed_expected_counts$expected_prob,
                                        rescale.p = TRUE)
goodness_of_fit_chisqtest
```

ראשית נקבעו הגבולות של הקבוצות בהתפלגות (הדיסקרטיזציה, יצירת התפלגות בדידה), כך שההתפלגות תכיל את כל תוצאות ההתפלגות האפשריות מ-0 ועד 1.

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

לאחר מכן אנחנו משתמשים בפקודה `chisq.test`. שימו לב שלמעשה לא השתמשנו במספר התצפיות הצפוי אלא בהתפלגות שלהם (הארגומנט `p`), שמעיד מה ההסתברות לקבל דגימות בתאים הסטטיסטיים הנבחנים.

במקרה זה קיבלנו `p.value=0.801` כלומר לא ניתן לדחות את השערת האפס שההתפלגות הנדונה היא התפלגות נורמלית (לכן נאמר שיש בסיס להתייחס להתפלגות כהתפלגות נורמלית).

::: callout-note
במרבית ספרי הסטטיסטיקה תתיקלו בנוסחה הבאה עבור מבחן חי-בריבוע.

ראשית מחשבים את הסטטיסטי:\
$$
X=\sum_{i=1}^l\frac{\left(O_i-E_i\right)^2}{E_i}
$$

כאשר $l$ הוא מספר התאים הסטטיסטיים (בדוגמה הקודמת 5), ומספר דרגות החופש בהתפלגות יהיה $df=l-1$.

לאחר מכן מחשבים את ה-p-value לפי ערכו של הסטטיסטי:

$$
\Pr\left(X>\chi^2_{df}\right)
$$

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

```{r chi sq direct computation}
X <- ((observed_expected_counts$observed - 
        observed_expected_counts$expected_counts)^2) / 
  observed_expected_counts$expected_counts

p.value <- pchisq(q = sum(X), df = 4, 
                  lower.tail = F)
p.value
```
:::

### מבחן אי-תלות

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

```{r year and island penguin sampling}
penguins %>% 
  count(year, species) %>%
  group_by(year) %>% 
  mutate(prop = n/sum(n)) %>% 
  ggplot(aes(x = factor(year), y = prop, fill = species)) + 
  geom_col(position = position_fill()) + 
  scale_y_continuous(labels = scales::percent) + 
  ylab("Sample distribution") + xlab("") +
  ggtitle("Sample distribution of species over years") +
  geom_text(aes(label = glue::glue("{round(prop*100)}%")),
            position = position_fill(vjust = 0.5), show.legend = FALSE)
```

כפי שניתן להבחין בתרשים, ישנם הבדלים מסוימים. לדוגמה, בשנת 2007 רק 31% מהדגימות הם של פינגויינים מזן Gentoo לעומת שנת 2008 שם 40% מהדגימות הן מזן Gentoo. השערת האפס שלנו היא שהמשתנה שנה והמשתנה זן הפינגויין הם בלתי תלויים, וההשערה האלטרנטיבית הינה שהם תלויים.

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

```{r independence test variation one}
independent_chisqtest <- chisq.test(x = penguins$year, y = penguins$species)
independent_chisqtest
```

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

```{r independence test variation two}
penguin_species_year_counts <- penguins %>% 
  mutate(year = as.character(year)) %>% 
  count(year, species) %>% 
  pivot_wider(names_from = species, values_from = n) %>% 
  select(-year) %>% 
  as.matrix()

chisq.test(x = penguin_species_year_counts)
```

::: callout-note
# למה התוצאות זהות?

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

מבחינת פרשנות- ניתן לראות שערך ה-`p.value` הינו 0.5224, ולכן לא ניתן לדחות את השערת האפס (כלומר ייתכן שהמשתנים בלתי תלויים).

## מבחן פרופורציות

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

## מבחן פרופורציות עם מדגם יחיד

ראשית נדגים מבחן פרופורציות עם מדגם יחיד: נבצע מבחן סטטיסטי על פרופורצית הפינגויינים ששוקלים מעל 4 ק"ג, ונבחן האם פרופורציה זו שונה באופן מובהק מ-50%.

```{r proportion test}
penguins %>% 
  filter(!is.na(body_mass_g)) %>% 
  mutate(is_above_4kg = body_mass_g >= 4000) %>% 
  summarize(tot_above_4kg = sum(is_above_4kg),
            tot_under_4kg = sum(!is_above_4kg),
            total_specimens = n())
```

כפי שניתן לראות יש 165 פרטים עם משקל מתחת ל-4 ק"ג, ו-177 פרטים עם משקל מעל 4 ק"ג. סה"כ 342 פרטים. נשתמש בפקודה `prop.test` על מנת לבחון את ההשערה של חלוקה שווה (פרופורציה 50%).

```{r prop test penguins}
penguin_prop <- prop.test(x = 177, n = 342, p = 0.5)
penguin_prop
```

כפי שניתן לראות בפלט הפקודה, ערך ה-p-value הינו 0.552, ולכן לא ניתן לדחות את השערת האפס של חלוקה חצי/חצי. רווח בר סמך לפרופורציה של הפינגויינים בעלי משקל מעל 4 ק"ג עומד על בין 0.463 ועד 0.571.

::: callout-tip
# נקודות למחשבה

נסו לבצע את המבחן עם פרופורצית הפרטים מתחת למשקל 4 ק"ג. האם מתקבלת אותה התוצאה? מה שונה ומה דומה בפלט? הסבירו.
:::

## מבחן פרופורציות לשני מדגמים

נבצע מבחן המשווה את שיעור הפרטים בעלי משקל מעל 4.5 ק"ג המשווה בין זן ה-Gentoo ליתר הזנים.

```{r prop test two samples}
gentoo_weight <- penguins %>% 
  mutate(is_gentoo = species == "Gentoo") %>% 
  filter(!is.na(body_mass_g)) %>% 
  group_by(is_gentoo) %>% 
  summarize(above_4.5kg = sum(body_mass_g >= 4500),
            under_4.5kg = sum(body_mass_g < 4500),
            tot = n())

gentoo_weight

twosample_prop <- prop.test(x = c(11, 107), n = c(219, 123))
twosample_prop
```

כפי שניתן לראות על פי תוצאות המבחן, יש שיעור של כ-5% מהפינגויינים שאינם מזן Gentoo שיש להם משקל מעל 4.5 ק"ג, לעומת 87% מהפינגויינים שהם כן מזן Gentoo ובעלי משקל מעל 4.5 ק"ג. הבדל זה מובהק סטטיסטית עם p-value \< 2.2e-16 (כלומר קטן מהמספר הקטן ביותר שהמחשב יכול להציג).

רווח בר סמך להפרש בין אחוזים אלו עומד על בין 0.747 ועד 0.892 (בשורה בפלט זה מופיע בסימן שלילי משום שהפקודה מחסרת את הפרופורציה הראשונה מהשניה ובמקרה זה הפרופורציה הראשונה נמוכה מהשניה - 11 מתוך 219 לעומת 107 מתוך 123).

::: callout-tip
# הקשר שבין מבחני פרופורציות ומבחן חי-בריבוע לאי-תלות

נסו להשתמש במבחן חי-בריבוע `chisq.test` על מנת לבצע את מבחן הפרופורציות מהסעיף הקודם

רמז: עליכם לבנות מטריצה 2\*2 כאשר השורות הן מעל/מתחת למשקל שנקבע והעמודות הם שני הסוגים של הפינגויינים (כן/לא מסוג Gentoo.

האם קיבלתם את אותה התוצאה כמו במבחן הפרופורציות? (מה שונה ומה דומה?)

הוסיפו ארגומנט למבחן הפרופורציות: `correct = FALSE`, ארגומנט זה מבטל תיקון (מניפולציה) שמבצעת הפונקציה. האם כעת התוצאות זהות?
:::

## מבחן קולמוגורוב-סמירנוב (KS test) להתפלגות

מבחן קולמוגורוב-סמירנוב משמש להשוואה בין שתי התפלגויות, בדומה למבחן חי-בריבוע לטיב התאמה, אך הוא אינו דורש התפלגויות בדידות (כפי שמבחן חי-בריבוע דורש).

אפשר להפעיל אותו על מנת להשוות בין שתי התפלגויות רציפות נתונות או בין התפלגות רציפה נתונה לבין התפלגות רציפה תיאורתית. לדוגמה, נשוואה את אורך המקור של פינגויינים מזן Adelie להתפלגות רציפה תיאורתית נורמלית עם תוחלת `r round(length_mu, 2)`, וסטיית תקן `r round(length_sd, 2)`:

```{r kolmogorov-smirnov example}
ks_test_output <- ks.test(adelie$bill_length_mm, "pnorm", mean = length_mu, sd = length_sd)
ks_test_output
```

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

::: callout-note
# הודעת האזהרה של פלט המבחן `ks.test`

במסגרת המבחן שהרצנו קיבלנו הודעת אזהרה:

ties should not be present for the Kolmogorob-Smirnov test.

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

## סידור הפלט של מבחני השערות

בכל המקרים שבחנו בפרק זה הזנו למבחנים נתונים בפורמט מסודר של טבלאות (אם הפעלנו את ארגומנט ה-`data` וה-`formula` או בפורמט של וקטורים). תוצאות המבחנים הוצגו כפלט טקסטואלי, אך בפועל היו רשימה (list) עם רכיבים שונים שניתן להפריד אותם. לדוגמה במבחן האחרון שהרצנו (`ks.test`) ניתן היה לקרוא רק לערכו של ה-p-value על ידי שימוש בקוד: `ks_test_output$p.value`. ישנה חבילה הנקראת `broom` ובה פקודה הנקראת `tidy` שמסייעת לנו להפוך את הפלט של הפקודות השונות לטבלה (פחות או יותר אחידה).

נדגים זאת על חלק מהפלטים שראינו בפרק (כולל סידור של פורמט העמודות בטבלה באמצעות חבילת `gt`).

```{r example usage of broom no run}
#| warning: false
#| message: false
#| eval: false
library(broom)
bind_rows(
  tidy(adeliettest),
  tidy(paired_ttest),
  tidy(independent_chisqtest),
  tidy(penguin_prop),
  tidy(ks_test_output)
) %>% 
  gt::gt() %>% 
  gt::fmt_number(columns = c(1:2, 4:6), decimals = 2)
```

::: {style="zoom:80%; direction:ltr;"}
```{r example usage of broom}
#| warning: false
#| message: false
#| eval: true
#| echo: false
library(broom)
bind_rows(
  tidy(adeliettest),
  tidy(paired_ttest),
  tidy(independent_chisqtest),
  tidy(penguin_prop),
  tidy(ks_test_output)
) %>% 
  gt::gt() %>% 
  gt::fmt_number(columns = c(1:2, 4:6), decimals = 2)
```
:::

גישה זו שימושית במיוחד כאשר מריצים הרבה מבחנים זהים על קבוצות שונות, כפי שנדגים ב [פרק @sec-purrr-explained].

## סיכום

בפרק זה בחנו מבחנים סטטיסטיים מסוגים שונים:

התחלנו במבחני השערות על תוחלת של מדגם יחיד, מבחני השערות על תוחלת עבור שני מדגמים לא-מזווגים ועבור שני מדגמים מזווגים.

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

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

לבסוף למדנו כיצד ניתן להעביר את הפלט של המבחנים לפלט טבלאי מסודר.

::: end-page