סטטיסטיקה בסיסית והסקה סטטיסטית ב- R

1. סטטיסטיקה בסיסית והסקה סטטיסטית ב-R

2. Everything differs!!! "צפויים להימצא הבדלים בין x ל-y" היא אמירה טריוויאלית הסטטיסטיקאי שבכם שואל "האם ההבדלים שנמצאו גדולים מהצפוי באקראי" הביולוג שבכם שואל "למה ההבדלים הם לכיוון ובדרגה שמצאתי"

3. מדדים לנטייה מרכזית* • ממוצע חשבוני, גיאומטרי, הרמוני Arithmetic mean: Σxi/n Geometric mean: (x1*x2*…*xn)1/n Harmonic mean: * = Moments of central tendency

4. מדדים לנטייה מרכזית ב-R Arithmetic mean: Σxi/n 1. ממוצע חשבוני הפונקציה “mean” data<-c(2,3,4,5,6,7,8) mean(data) • [1] 5 דוגמא: 2. ממוצע גיאומטרי Geometric mean: (x1*x2*…*xn)1/n data<-c(2,3,4,5,6,7,8) exp(mean(log(data))) [1] 4.549163 דוגמא:

5. מדדים לנטייה מרכזית* • א. ממוצע (mean) • ב. חציון (median) • ג. שכיח (mode) * = Moments of central tendency data<-c(2,3,4,5,6,7,8) median(data) • [1] 5 דוגמא:

6. מדדים לנטייה מרכזית http://www.statmethods.net/management/functions.html • ממוצע • שונות* הממוצע הוא מדד טוב יותר למה שקורה באוכלוסיה\מדגם כשהשונות קטנה או גדולה? *Variance = Σ(xi-μ)2 / n data<-c(2,3,4,5,6,7,8) var(data) • [1] 4.666667 דוגמא:

7. מדדים לנטייה מרכזית • ממוצע • שונות המומנט השני של נטייה מרכזית הוא מדד לפיזור הנתונים מסביב למומנט הראשון דוגמאות למומנט השני הן, למשל, השונות, סטיית התקן, שגיאת התקן, ה-coefficient of variation, ורווח הסמך (של 90, 95, 99% או מה שלא יהיה)

8. מדדים לנטייה מרכזית #for: data<-c(2,3,4,5,6,7,8) גודל מדגם: שונות: סטיית תקן: שגיאת תקן: coefficient of variation: length(data) var(data) sd(data) se<-(sd(data)/length(data)^0.5) se [1] 0.8164966 CV<-sd(data)/mean(data) CV [1] 0.4320494

9. מדדים לנטייה מרכזית • ממוצע • שונות • הטייה (Skew) התפלגות שכיחויות מוטה אינה סימטרית! האם בהתפלגות שכיחויות מוטה הממוצע החשבוני הוא מדד טוב לנטייה מרכזית? מהו השכר הממוצע של כל הסטודנטים פה ושל ביל גייטס?

10. מדדים לנטייה מרכזית הטייה (Skew) skew<-function(data){ m3<-sum((data-mean(data))^3)/length(data) s3<-sqrt(var(data))^3 m3/s3} skew(data) sdskew<-function(x) sqrt(6/length(x)) שגיאת תקן של הטיה:

11. מדדים לנטייה מרכזית • ממוצע • שונות • הטייה (Skew) • קורטוזיס (Kurtosis)

12. מדדים לנטייה מרכזית • קורטוזיס (Kurtosis) kurtosis<-function(x){ m4<-sum((x-mean(x))^4)/length(x) s4<-var(x)^2 m4/s4-3 } kurtosis(x) sdkurtosis<-function(x) sqrt(24/length(x)) שגיאת תקן של קורטוזיס:

13. התפלגות נורמאלית יכולה לקבל כל ערך של ממוצע ושונות, אבל ה skewness והקורטוזיס שלה צריכים להיות שווים לאפס לערכי skew וקורטוסיס יש שונות משלהם – ואפס צריך להיות מחוץ לרווח הסמך שלהם כדי שהם יהיו שונים במובהק מאפס

14. Residuals כשאנחנו עושים סטטיסטיקה אנחנו יוצרים מודלים של המציאות אחד המודלים הפשוטים ביותר הוא הממוצע: הגובה הממוצע של אזרחי ישראל הוא (נגיד) 173 ס"מ השכר הממוצע הוא 9302 ₪ (למ"ס, נתוני מרץ 2013) והשירות הצבאי הממוצע הוא (אולי) 24 חודשים 2.05 מטר דוב ליאור שירת בצה"ל חודש אחד 46,699 ₪ בחודש (הערכת חסר, ללא הטבות) http://www.haaretz.co.il/1.2057452

15. Residuals כשאנחנו עושים סטטיסטיקה אנחנו יוצרים מודלים של המציאות כך שכאן, המודלים שלנו: 173 ס"מ, 9302 ₪, 24 חודשים, אינם מוצלחים במיוחד ה-Residual הוא הכמות בה ערך מסויים רחוק מהניבוי של המודל. כך שליאור אליהו רחוק 32 ס"מ מהמודל "ישראלי = 173", ורחוק 29 ס"מ מהמודל המורכב יותר "גבר ישראלי = 177, אישה ישראלית = 168". Residual = ₪ 37397 Residual = 32 cm Residual = -23 month IDF service

16. Residuals כשאנחנו עושים סטטיסטיקה אנחנו יוצרים מודלים של המציאות model<-lm(size~Species+sex+Latitude+Longitude) out<-model$residuals write.table(out, file = "residuals.txt",sep="\t",col.names=F,row.names=F) #note that residual values are in the order entered (i.e., not alphabetic, not by residual size – first in, first out) Residual = ₪ 37397 Residual = 32 cm Residual = -23 month service

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

18. איזה מבחן נבחר? זה משתנה בהתאם לטבע ה response variable (=המשתנה התלוי, זה שעל ציר ה-y), ובעיקר לפי טבע ה-predictor variables • אם ה-response variableהוא "הצלחה או כשלון", והשערת האפס היא של שיוויון ביניהם, נשתמש במבחן בינומי (binomial) • אם ה response variableשלנו הוא ספירות נשתמש לרוב במחני חי-בריבוע או G (=log-likelihood). • אבל לרוב ה response variableשלנו יהיה רציף (14 מינים, 78 פרטים, 87.5מ"מ, 54 פעימות לדקה, 7.3 ביצים, 23 מעלות)

19. איזה מבחן נבחר? מהו ה-response variable ? רציף (14 מינים, 78 פרטים, 87.5מ"מ, 54 חודשים, 7.3 ביצים, 23 מעלות) ספירות (שכיחויות: 6 זכרים, 9 נקבות) "הצלחה" או "כישלון" (מצא את הגבינה\אידיוט) חי-בריבוע או G (=log-likelihood) מבחן בינומי (binomial) ראה בהמשך

20. Binomial test in R יש להגדיר את מספר ההצלחות מתוך גודל המדגם הכולל. דוגמה 1: 19 מתוך 34 (לא מובהק). דוגמה 2: 19 מתוך 20 (מובהק) binom.test(19,34) Exact binomial test data: 19 and 34 number of successes = 19, number of trials = 34 p-value = 0.6076 alternative hypothesis: true probability of success is not equal to 0.5 95 percent confidence interval: 0.3788576 0.7281498 sample estimates: probability of success 0.5588235 binom.test(19,20) Exact binomial test data: 19 and 20 number of successes = 19, number of trials = 20, p-value = 4.005e-05 alternative hypothesis: true probability of success is not equal to 0.5 95 percent confidence interval: 0.7512672 0.9987349 sample estimates: probability of success 0.95

21. Chi-square test in R • Data: insularity vs. diet: chisq.test M<-as.table(rbind(c(1901,101,269),c(488,43,177))) chisq.test(M) data: M X-squared = 80.0441, df = 2, p-value < 2.2e-16

22. איזה מבחן נבחר? אם ה response variableשלנו הוא רציףנבחר מבחן לפי טבע ה-predictor variables • אם ה-predictor variableבדיד (אתר א', אתר ב', אתר ג'; זכר\נקבה; מין א', מין ב, מין ג'; שטח עירוני\שטח טבעי; טיפול א', טיפול ב', ביקורת) המבחן יהיה ANOVA (analysis of variance) • אם ה predictor variableרציף (מעלות צלסיוס, כמות מזון, קו רוחב, כמות משקעים, מספר מתחרים, אחוז כיסוי) המבחן יהיהרגרסיה

23. t-test in R t.test(x,y) amy<-read.csv("ssd.csv",header=T) names(amy) kadison<-read.csv("ssd2.csv",header=T) names(kadison) attach(kadison) males<-size[Sex=="male"] females<-size[Sex=="female"] t.test(females,males) Welch Two Sample t-test data: females and males t = -2.1541, df = 6866.57, p-value = 0.03127 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -7.5095545 -0.3536548 sample estimates: mean of x mean of y 88.17030 92.10191

24. t-test in R (2) lm(x~y) amy<-read.csv("ssd.csv",header=T) names(amy) kadison<-read.csv("ssd2.csv",header=T) names(kadison) model<-lm(size~Sex,data=kadison) summary(model)

25. Paired t-test in R t.test(x,y,paired=TRUE) amy<-read.csv("ssd.csv",header=T) names(amy) kadison<-read.csv("ssd2.csv",header=T) names(kadison) attach(kadison) males<-size[Sex=="male"] females<-size[Sex=="female"] t.test(females,males,paired=TRUE) Paired t-test data: females and males t = -10.1917, df = 3503, p-value < 2.2e-16 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -4.687950 -3.175259 sample estimates: mean of the differences -3.931605 • tapply(size,Sex,mean)

26. ANOVA in R aov model<-aov(x~y) ochel<-read.table("di.txt",header=T) names(ochel) [1] "species" "diet" "size" model<-aov(size~diet,data=ochel) summary(model)

27. מבחן post-hoc ל-ANOVA ב-R TukeyHSD(model) Fit: aov(formula = size ~ diet, data = ochel) $diet כל ההשוואות מובהקות

28. מבחן post-hoc ל-ANOVA ב-R model<-aov(SVL~Realm) summary(model) Df Sum Sq Mean Sq F value Pr(>F) Realm 11 5.270 0.479 9.3242 < 2.2e-16 *** Residuals 4851 249.270 0.051 TukeyHSD(model) Fit: aov(formula = SVL ~ Realm) $Realm אפס מחוץ לרווח הסמך

29. correlation in R cor.test(x,y) rapoport<-read.table("rang.txt",header=T) names(rapoport) [1] "latitude" "range_size“ • attach(rapoport) • cor.test(range_size,latitude) • Pearson's product-moment correlation • data: range_size and latitude • t = 9.9823, df = 4910, • p-value < 2.2e-16 • cor 0.1410353 המשתנה “cor” הוא מקדם הקורלציה r

30. regression in R lm (=“linear model”): lm (y~x) אותם נתונים כמו בדוגמה הקודמת model<-lm(range_size~latitude,data=rapoport) summary(model) Call: lm(formula = range_size ~ latitude, data = rapoport) Residuals: Min 1Q Median 3Q Max -4.708 -1.774 0.470 1.465 3.725 Residual standard error: 1.844 on 4910 degrees of freedom Multiple R-squared: 0.01989, Adjusted R-squared: 0.01969 F-statistic: 99.65 on 1 and 4910 DF, p-value: < 2.2e-16

31. aov לעומת lm אולי במפתיע אפשר לבחון נתונים המתאימים ל-ANOVA גם במבחן lm. במקרה כזה נקבל את כל המידע שנותן ה-summary של lm במבחן רגרסיה, כולל (חשוב!) parameter estimates, שגיאות תקן, הבדלים בין פקטורים וערכי-p לכל קונטרסט (בין 2 קטגוריות של המשתנה המסביר הקטגוריאלי)

32. aov לעומת lm אולי במפתיע אפשר לבחון נתונים המתאימים ל-ANOVA גם במבחן lm. במקרה כזה נקבל את כל המידע שנותן ה-summary של lm במבחן רגרסיה, כולל (חשוב!) parameter estimates, שגיאות תקן, הבדלים בין פקטורים וערכי-p לכל קונטרסט (בין 2 קטגוריות של המשתנה המסביר הקטגוריאלי) model2<-lm(size~diet,data=ochel) summary(model2) • Residual standard error: 0.6949 on 2959 degrees of freedom • Multiple R-squared: 0.1258, F-statistic: 212.9 on 2 & 2959 DF, p-value: < 2.2e-16 עוד בהמשך

33. Multiple predictors מה לעשות, החיים מסובכים. לפעמים מה שמעניין אותנו מושפע מיותר מגורם אחד! קצב ליבם של זוחלים, למשל, מושפע גם מגודל גופם וגם מטמפרטורת הסביבה, וגם מהזמן והמהירות בהם נעו לאחרונה Smith, R. J. 1999. Statistics of sexual size dimorphism. Journal of Human Evolution 36: 423-459.

34. Multiple predictors מה לעשות, החיים מסובכים. לפעמים מה שמעניין אותנו מושפע מיותר מגורם אחד! קצב ליבם של זוחלים, למשל, מושפע גם מגודל גופם וגם מטמפרטורת הסביבה, וגם מהזמן והמהירות בהם נעו לאחרונה ניתן להסביר אם כן את המשתנה המעניין (קצב לב) אם יש לנו מידע על כל המשתנים המסבירים. ההנחה היא שכאשר אנו מכניסים את שלושתם למשוואה אנו רואים את השפעתו של כל אחד כאשר שני האחרים "מוחזקים קבועים" (held constant) Smith, R. J. 1999. Statistics of sexual size dimorphism. Journal of Human Evolution 36: 423-459.

35. Multiple predictors ניתן להסביר אם כן את המשתנה המעניין (קצב לב) אם יש לנו מידע על כל המשתנים המסבירים. ההנחה היא שכאשר אנו מכניסים את שלושתם למשוואה אנו רואים את השפעתו של כל אחד כאשר שני האחרים "מוחזקים קבועים" (held constant) הנחה זו נכונה כאשר אין מתאם (גבוה) בין המשתנים המסבירים לבין עצמם Smith, R. J. 1999. Statistics of sexual size dimorphism. Journal of Human Evolution 36: 423-459.

36. איזה מבחן נבחר? אם יש כמה predictor variables(נאמר ארבעה) וכולם קטגוריאלים המבחן יהיה ANOVA(נאמר 4-way ANOVA) אם יש כמה predictor variables(נאמר שבעה) וכולם רציפים המבחן יהיה Multiple Regression

37. איך כותבים מבחן עם כמה משתנים מסבירים? lm(y~a+b+c) משתמשים בפלוס (+) בין המשתנים המסבירים. • למשל במודל שמנסה לחזות ציוני קורסים לפי כמה למדנו, גיל המרצה, כמה התפללנו והאם מעתיקים: model<-lm(Grade~days_studied+age_of_professor+number_ofprayers+sent_sms_texts,data=grades) summary(model)

38. איזה מבחן נבחר? אם יש כמה predictor variables(לפחות 2) – חלקם (לפחות 1) בדידים וחלקם (לפחות 1) רציפים המבחן יהיה ANCOVA (analysis of co-variance)

39. איך זה נראה גרפית? דוגמא: מדדתי אורך שלוש שיניים בשועלים מצויים, זכרים ונקבות, מכל תחום תפוצתם Vulpes vulpes ANOVA שני הגורמים מובהקים: יש הבדל בין השיניים ובין הזוויגים Regression ANCOVA ניב ניב עליון שן שסע תחתונה שן שסע תחתונה שן שסע עליונה שן שסע עליונה

40. איך זה נראה גרפית? דוגמא: אורך שלוש שיניים בשועלים מצויים, כפונקציה של קו הרוחב R-squared: 0.015, F = 32.83, 1 & 2161 DF; p < 0.0001 ANOVA Regression ANCOVA יש כלל ברגמן אבל קל לראות שהמודל* זוועתי: ניבים קטנים יותר משיני שסע *קו הרגרסיה הוא מודל ליחס בין ה predictor ל response

41. b. a. c. response response בדיד מובהק רציף לא Null hypothesis רציף מובהק, בדיד לא response Continuous predictor Continuous predictor Continuous predictor d. שניהם מובהקים response Continuous predictor ANCOVA גרפית קל להבין את זה גרפית: בדוגמה משתנה בדיד אחד עם 2 רמות (מקווקו ומלא), ומשתנה רציף אחד (לאורך ציר ה X) ותודה לדניאל על הגרפים

42. איך זה נראה גרפית? דוגמא: אורך שיניים בשועלים מצויים, כפונקציה של קו הרוחב (ציר ה-X) הזוויג (צבע) ובאיזו שן מדובר (צורה) ANOVA Regression ANCOVA כל הגורמים מובהקים המודל הזה מסביר 96.3% מהשונות העברנו פה שישה קווי רגרסיה מקבילים

43. קריאת תוצאות מודל ANCOVA ב-R Response = intercept + a for level 1 of the 1st categorical predictor variableor + b for level 2 of the 1st + c for level 1 of the 2nd categorical predictorord for level 2… +k*(value of the continuous predictor variable) + error דוגמא ללא אינטראקציות למשל, אם נחזור לשועלים Tooth / sex Latitude כך שהאורך של שן P בזכרים בתל אביב (קו רוחב 32) על פי המודל הוא: 12.726 = 4.342+6.617+0.391+32*0.043

44. איך קוראים תוצאות של lm ב-R כשהמשתנה המסביר קטגוריאלי, R משווה את כל הפקטורים ל-intercept של הפקטור הראשון בסדר אלפביתי כאן המשתנה המסביר הוא דיאטה ואלפביתית הקטגוריה הראשונה היא Carnivore ההבדל בין קרניבורים להרביבורים (t = 19.3) ואומניבורים (t = 9.075) מובהק כך שגודלו הממוצע של קרניבור הוא , 1.038 וזה של אומניבור: 1.038+0.329 = 1.367

45. איך קוראים תוצאות של lm ב-R כשהמשתנה המסביר רציף, R מדווח עבורו את השיפוע עם שגיאת התקן שלו וערכי ה-p וה-t שלו • rapoport<-read.table("rang.txt",header=T) • model<-lm(range~mass+latitude,data=rapoport) • summary(model) כאן המשתנים המסבירים (את גודל טווח התפוצה) הם מסה וקו רוחב ההשפעה של שניהם מובהקת (t = 8.82, ו-t = 9.57, p<<0.05) כך שגודלו תחום התפוצה עולה ב-0.02 (היחידות הן לוג קמ"ר) בכל עליה של מעלה, וב-0.326 יחידות כל עליה של יחידת מסה (לוג גרם, כלומר יחידה אחת היא הכפלה בעשר)

46. איך קוראים תוצאות של lm ב-R כשהמשתנה המסביר רציף, R מדווח עבורו את השיפוע עם שגיאת התקן שלו וערכי ה-p וה-t שלו כך שגודל תחום התפוצה של לטאה במסה של 100 גרם (log(100) = 2) בקוסטה ריקה (קו רוחב 10) יהיה: Intercept+slope*mass+slope*latitude 2.967+0.3255(slope)*2(mass)+0.023(slope)*10latitude = 3.84

47. איך קוראים תוצאות של lm ב-R בANCOVA כשיש גם המשתנים מסבירים רציפים וגם קטגוריאלים, R מדווח intercept לאחרונים ושיפוע לראשונים, עם שגיאות התקן שלו וערכי ה-p וה-t המתאימים model<-lm(range~islands+mass+latitude,data=rapoport) summary(model) כאן יש 3 קטגוריות:מיני יבשת (continental) מינים שאנדמיים לאי יחיד (single island endemics) או מכמה איים (multiple island endemics)

48. איך קוראים תוצאות של lm ב-R ANCOVA: גם משתנים מסבירים רציפים וגם קטגוריאלים model<-lm(range~islands+mass+latitude,data=rapoport) summary(model) 3 קטגוריות:יבשת (continental), אי יחיד (SIE) או כמה איים (MIE) Residual standard error: 1.633 on 4887 degrees of freedom Multiple R-squared: 0.232, Adjusted R-squared: 0.2314 F-statistic: 369 on 4 and 4887 DF, p-value: < 2.2e-16 אלה ערכי R בריבוע, F וכו' של המודל בכללותו

49. איך קוראים תוצאות של lm ב-R model<-lm(range~islands+mass+latitude,data=rapoport) summary(model) כיוון שאלפביתית continental<MIE<SIE החיתוך (intercept) שלנו הוא עבור הקטגוריה הראשונה: מיני יבשת. כך שטווח התפוצה של מיני יבשת גדול משל MIE (שימו לב: הבדל שלילי!) ושל SIE – וכמו כן טווח התפוצה עולה עם המסה, אך יורד עם העליה בקו הרוחב (שיפוע שלילי). ושוב – כל הגורמים מובהקים

50. relevel model<-lm(range~islands+mass+latitude,data=rapoport) summary(model) אבל בגורמים הקטגוריאלים יש לנו בעיה: R מחשב רק את ההבדל בין כל גורם לגורם הראשון באלפבית. כאן בין יבשת לשתי קטגוריות האיים. אבל הוא לא מחשב ולא מדווח על ההבדלים בין שתי קטגוריות האיים וחוץ מזה, הוא לא נותן לנו עבורם שגיאות תקן והבדלים מאפס, אלא רק הבדלים מהיבשת ושגיאת תקן של המבחן הזה (לא את שגיאת התקן של המשתנה עצמו במודל)