1.11k likes | 1.49k Views
סטטיסטיקה בסיסית והסקה סטטיסטית ב- R. Everything differs !!!. "צפויים להימצא הבדלים בין x ל- y " היא אמירה טריוויאלית. הסטטיסטיקאי שבכם שואל "האם ההבדלים שנמצאו גדולים מהצפוי באקראי " הביולוג שבכם שואל "למה ההבדלים הם לכיוון ובדרגה שמצאתי ". מדדים לנטייה מרכזית *. ממוצע.
E N D
Everything differs!!! "צפויים להימצא הבדלים בין x ל-y" היא אמירה טריוויאלית הסטטיסטיקאי שבכם שואל "האם ההבדלים שנמצאו גדולים מהצפוי באקראי" הביולוג שבכם שואל "למה ההבדלים הם לכיוון ובדרגה שמצאתי"
מדדים לנטייה מרכזית* • ממוצע חשבוני, גיאומטרי, הרמוני Arithmetic mean: Σxi/n Geometric mean: (x1*x2*…*xn)1/n Harmonic mean: * = Moments of central tendency
מדדים לנטייה מרכזית ב-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 דוגמא: ניתן גם לעבוד מהקובץ: • dat<-read.csv("island_type_final2.csv") • Attach(dat) • mean(lat) • [1] 17.40439 data<-c(2,3,4,5,6,7,8) exp(mean(log(data))) [1] 4.549163
מדדים לנטייה מרכזית* • א. ממוצע (mean) • ב. חציון (median) • ג. שכיח (mode) * = Moments of central tendency דוגמא: דוגמאמהקובץ: data<-c(2,3,4,5,6,7,8) median(data) • [1] 5 median(mass) [1] 0.69
מדדים לנטייה מרכזית http://www.statmethods.net/management/functions.html • ממוצע • שונות* הממוצע הוא מדד טוב יותר למה שקורה באוכלוסיה\מדגם כשהשונות קטנה או גדולה? דוגמא: • data<-c(2,3,4,5,6,7,8) • var(data) • [1] 4.666667 • Var(lat) • [1] 89.20388 *Variance = Σ(xi-μ)2 / n
מדדים לנטייה מרכזית • ממוצע • שונות המומנט השני של נטייה מרכזית הוא מדד לפיזור הנתונים מסביב למומנט הראשון דוגמאות למומנט השני הן, למשל, השונות, סטיית התקן, שגיאת התקן, ה-coefficient of variation, ורווח הסמך (של 90, 95, 99% או מה שלא יהיה).
מדדים לנטייה מרכזית #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
מדדים לנטייה מרכזית • ממוצע • שונות • הטייה (Skew) התפלגות שכיחויות מוטה אינה סימטרית! האם בהתפלגות שכיחויות מוטה הממוצע החשבוני הוא מדד טוב לנטייה מרכזית? מהו השכר הממוצע של כל הסטודנטים פה ושל ביל גייטס?
מדדים לנטייה מרכזית הטייה (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)) שגיאת תקן של הטיה:
מדדים לנטייה מרכזית • ממוצע • שונות • הטייה (Skew) • קורטוזיס (Kurtosis)
מדדים לנטייה מרכזית • קורטוזיס (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)) שגיאת תקן של קורטוזיס:
התפלגות נורמאלית יכולה לקבל כל ערך של ממוצע ושונות, אבל ה-skewness והקורטוזיס שלה צריכים להיות שווים לאפס לערכי skew וקורטוסיס יש שונות משלהם – ואפס צריך להיות מחוץ לרווח הסמך שלהם כדי שהם יהיו שונים במובהק מאפס
Residuals כשאנחנו עושים סטטיסטיקה אנחנו יוצרים מודלים של המציאות אחד המודלים הפשוטים ביותר הוא הממוצע: הגובה הממוצע של אזרחי ישראל הוא (נגיד) 173 ס"מ השכר הממוצע הוא 9271 ₪ (למ"ס, נתוני אפריל 2014) והשירות הצבאי הממוצע הוא (אולי) 24 חודשים 2.06 מטר דוב ליאור שירת בצה"ל חודש אחד 46,699 ₪ בחודש (הערכת חסר, ללא הטבות) http://www.haaretz.co.il/1.2057452
Residuals כשאנחנו עושים סטטיסטיקה אנחנו יוצרים מודלים של המציאות כך שכאן, המודלים שלנו: 173 ס"מ, 9271 ₪, 24 חודשים, אינם מוצלחים במיוחד ה-Residual הוא הכמות בה ערך מסויים רחוק מהניבוי של המודל. כך שליאור אליהו רחוק 32 ס"מ מהמודל "ישראלי = 173", ורחוק 29 ס"מ מהמודל המורכב יותר "גבר ישראלי = 177, אישה ישראלית = 168". Residual = ₪ 37428 Residual = 33 cm Residual = -23 month IDF service
Residuals כשאנחנו עושים סטטיסטיקה אנחנו יוצרים מודלים של המציאות dat<-read.csv("island_type_final2.csv") model<-lm(mass~iso+area+age+lat, data=dat) out<-model$residuals out 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 = ₪ 37428 Residual = 33 cm Residual = -23 month service
סטטיסטיקה תיאורית והסקה סטטיסטית כשיש לנו נתונים בשלב הראשון כדאי שנתאר אותם: נצייר גרפים, נחשב ממוצע וכדומה • בהסקה סטטיסטית אנו בוחנים את התנהגותם של הנתונים שלנו אל מול השערה (היפותזה) מסויימת • את ההשערה שלנו אנו יכולים להציג כמודל סטטיסטי • למשל: • התפלגות גבהים היא נורמאלית • מספר המינים הולך ועולה עם העליה בשטח • מספר המינים הולך ועולה עם העליה בשטח על פי power functionשלו מעריך חזקה השווה ל-0.25
התפלגות שכיחויות* כמה תצפיות יש בכל קטגוריה (bin) • dat<-read.csv("island_type_final2.csv") • attach(dat) • names(dat) • Hist(mass) מבטאת את כל מרחב התצפיות * = frequency distribution, in graphic form = “histogram”
התפלגות שכיחויות מה למדנו? • dat<-read.csv("island_type_final2.csv") • attach(dat) • Hist(mass) • אין מסות קטנות מעשירית הגרם או גדולות מ-100 ק"ג. • לטאות במסות של 1-10 גרם נפוצות – גדולות יותר או קטנות יותר נדירות • ההתפלגות היא חד שיאית, ומוטה לימין
התפלגות שכיחויות לא חייבת להיות כל כך מכוערת • dat<-read.csv("island_type_final2.csv") • attach(dat) • hist(mass, col="purple",breaks=25,xlab="log mass (g)",main="masses of island lizards - great data by Maria",cex.axis=1.2,cex.lab=1.5)
הצגת משתנה קטגורי לעומת משתנה רציף אחר • dat<-read.csv("island_type_final2.csv") • attach(dat) • plot(type,brood) תמיד העדיפו Box & Whiskers plotעלbar plot!
הצגת משתנה רציף לעומת משתנה רציף אחר • dat<-read.csv("island_type_final2.csv") • attach(dat) • plot(mass,clutch) • plot(mass,clutch,pch=16,col="blue“)
איזה מבחן נבחר? זה משתנה בהתאם לטבע ה-response variable(=המשתנה התלוי, זה שעל ציר ה-y), ובעיקר לפי טבע ה-predictor variables • אם ה-response variableהוא "הצלחה או כשלון", והשערת האפס היא של שיוויון ביניהם, נשתמש במבחן בינומי (binomial). • אם ה-response variableשלנו הוא ספירות נשתמש לרוב במחני 2χאו G (=log-likelihood). • לרוב ה-response variableשלנו יהיה רציף (14 מינים, 78 פרטים, 87.5מ"מ, 54 פעימות לדקה, 7.3 ביצים, 23 מעלות).
איזה מבחן נבחר? מהו ה-response variable ? רציף (14 מינים, 78 פרטים, 87.5מ"מ, 54 חודשים, 7.3 ביצים, 23 מעלות) ספירות (שכיחויות: 6 זכרים, 9 נקבות) "הצלחה" או "כישלון" (מצא את הגבינה\אידיוט) חי-בריבוע או G (=log-likelihood) מבחן בינומי (binomial) עוד מעט...
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
Chi-square test in R • Data: lizard insularity & diet: chisq.test M<-as.table(rbind(c(1901,101,269),c(488,43,177))) chisq.test(M) data: M χ2 = 80.04, df = 2, p-value < 2.2e-16
Chi-square test in R והפעם עם בסיס הנתונים שלנו: chisq.test • dat<-read.csv("island_type_final2.csv") • install.packages("reshape") • Library(reshape) • cast(dat, type ~ what, length) M<-as.table(rbind(c(7,45,45),c(1,30,14),c(23,110,44))) chisq.test(M) data: M χ2 = 17.568, df = 4, p-value = 0.0015
איזה מבחן נבחר? אם ה response variableשלנו הוא רציףנבחר מבחן לפי טבע ה-predictor variables • אם ה-predictor variableבדיד (אתר א', אתר ב', אתר ג'; זכר\נקבה; מין א', מין ב, מין ג'; שטח עירוני\שטח טבעי; טיפול א', טיפול ב', ביקורת) המבחן יהיה ANOVA(analysis of variance) • אם ה predictor variableרציף (מעלות צלסיוס, כמות מזון, קו רוחב, כמות משקעים, מספר מתחרים, אחוז כיסוי) המבחן יהיהרגרסיה
t-test in R t.test(x,y) dimorphism<-read.csv("ssd.csv",header=T) attach(dimorphism) names(dimorphism) 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
t-test in R (2) lm(x~y) dimorphism<-read.csv("ssd.csv",header=T) attach(dimorphism) names(dimorphism) model<-lm(size~Sex,data=dimorphism) summary(model)
Paired t-test in R t.test(x,y,paired=TRUE) dimorphism<-read.csv("ssd.csv",header=T) attach(dimorphism) names(dimorphism) males<-size[Sex=="male"] females<-size[Sex=="female"] t.test(females,males, paired=TRUE) Paired t-test data: females and males t = -10.192, df = 3503, p-value < 2.2e-16 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -4.688 -3.175 sample estimates: mean of the differences -3.931 • tapply(size,Sex,mean)
ANOVA in R aov model<-aov(x~y) • island<-read.csv("island_type_final2.csv",header=T) names(island) • [1] "species" "what" "family" "insular" "Archipelago" "largest_island" [7] "area" "type" "age" "iso" "lat" "mass" • [13] "clutch" "brood" "hatchling" "productivity“ • model<-aov(clutch~type,data=island) summary(model)
מבחן post-hoc ל-ANOVA ב-R TukeyHSD(model) ההבדלים אינם מובהקים (שימו לב שאפס תמיד ברווח הסמך. ההבדל בין איי מדף יבשת לאיי פלטות טקטוניות קרוב למובהקות, p = 0.056)
correlation in R cor.test(x,y) • island<-read.csv("island_type_final2.csv",header=T) • names(island) • [1] "species" "what" "family" "insular" "Archipelago" "largest_island" • [7] "area" "type" "age" "iso" "lat" "mass" • [13] "clutch" "brood" "hatchling" "productivity“ • attach(island) • cor.test(mass,lat) • Pearson's product-moment correlation • data: mass and lat • t = -1.138, df = 317, p-value = 0.256 • alternative hypothesis: true correlation is not equal to 0 • 95 percent confidence interval: -0.17239 0.04635 • sample estimates: cor -0.06378 המשתנה “cor” הוא מקדם הקורלציה r
regression in R אותם נתונים כמו בדוגמה הקודמת lm (=“linear model”): lm (y~x) model<-lm(mass~lat,data=island) summary(model) • Call: lm(formula = mass ~ lat, data = island) • Residuals: • Min 1Q Median 3Q Max • -4.708 -1.774 0.470 1.465 3.725 • Residual standard error: 0.8206 on 317 degrees of freedom • Multiple R-squared: 0.004069, Adjusted R-squared: 0.0009268 • F-statistic: 1.295 on 1 and 317 DF, p-value: 0.256
aov לעומת lm אולי במפתיע אפשר לבחון נתונים המתאימים ל-ANOVA גם במבחן lm. במקרה כזה נקבל את כל המידע שנותן ה-summary של lm במבחן רגרסיה, כולל (חשוב!) parameter estimates, שגיאות תקן, הבדלים בין פקטורים וערכי-p לכל קונטרסט (בין 2 קטגוריות של המשתנה המסביר הקטגוריאלי)
aov לעומת lm אולי במפתיע אפשר לבחון נתונים המתאימים ל-ANOVA גם במבחן lm. במקרה כזה נקבל את כל המידע שנותן ה-summary של lm במבחן רגרסיה, כולל (חשוב!) parameter estimates, שגיאות תקן, הבדלים בין פקטורים וערכי-p לכל קונטרסט (בין 2 קטגוריות של המשתנה המסביר הקטגוריאלי) • island<-read.csv("island_type_final2.csv",header=T) • model<-aov(clutch~type,data=island) • model2<-lm(clutch~type,data=island) • summary(model) • summary(model2) טבלת aov טבלת lm • Residual standard error: 0.2893 on 289 degrees of freedom • (27 observations deleted due to missingness) • Multiple R-squared: 0.0189, Adjusted R-squared: 0.01211 • F-statistic: 2.784 on 2 and 289 DF, p-value: 0.06346 עוד בהמשך
הנחותיהם של מבחנים סטטיסטיים (כל המבחנים הסטטיסטיים) A non-random, non-independent sample of Israeli people • מדגם אקראי (הנחה של כל מבחן – לא רק פרמטרי) • אי תלות (מרחבית, בזמן, פילוגנטית וכו')
הנחותיהם של מבחנים פרמטריים א. ANOVA בנוסף להנחות של כל מבחן • שיוויון שוניות (Homoscedasticity) • התפלגות נורמלית של ה-residuals "Comments on earlier drafts of this manuscript made it clear that for many readers who analyze data but who are not particularly interested in statistical questions, any discussion of statistical methods becomes uncomfortable when the term ‘‘error variance’’ is introduced.“ Smith, R. J. 2009. Use and misuse of the reduced major axis for line-fitting. American Journal of Physical Anthropology 140: 476-486. Richard Smith & 3 friends קריאה: Sokal & Rohlf 1995. Biometry. 3rd edition. Pages 392-409 (especially 406-407 for normality)
Always look at your data Don’t just rely on the statistics! Anscombe's quartet Summary statistics are the same for all four data sets: • mean (7.5), • standard deviation (4.12) • correlation (0.816) • regression line (y = 3 + 0.5x) Anscombe 1973. Graphs in statistical analysis. The American Statistician27: 17–21. Briggs Henan University 2010
הנחותיהם של מבחנים פרמטריים ב: רגרסיה • שוויון שונויות (Homoscedasticity) Smith, R. J. 2009. Use and misuse of the reduced major axis for line-fitting. American Journal of Physical Anthropology 140: 476-486.
הנחותיהם של מבחנים פרמטריים ב: רגרסיה • שוויון שונויות (Homoscedasticity) • המשתנה המסביר (explanatory variable) נמדד ללא שגיאה Smith, R. J. 2009. Use and misuse of the reduced major axis for line-fitting. American Journal of Physical Anthropology 140: 476-486.
הנחותיהם של מבחנים פרמטריים ב: רגרסיה • שוויון שונויות (Homoscedasticity) • התפלגות נורמלית של ה residuals עבור כל ערך של המשתנה המסביר • התפלגות נורמלית של ה-residualsעבור כל ערך של המשתנה המסביר Smith, R. J. 2009. Use and misuse of the reduced major axis for line-fitting. American Journal of Physical Anthropology 140: 476-486.
הנחותיהם של מבחנים פרמטריים ב: רגרסיה • שוויון שונויות (Homoscedasticity) • המשתנה המסביר (explanatory variable) נמדד ללא שגיאה • התפלגות נורמלית של ה residuals עבור כל ערך של המשתנה המסביר • שוויון שונויות עבור כל ערך של המשתנה המסביר Smith, R. J. 2009. Use and misuse of the reduced major axis for line-fitting. American Journal of Physical Anthropology 140: 476-486.
הנחותיהם של מבחנים פרמטריים ב: רגרסיה • שוויון שונויות (Homoscedasticity) • המשתנה המסביר (explanatory variable) נמדד ללא שגיאה • התפלגות נורמלית של ה residuals עבור כל ערך של המשתנה המסביר • שוויון שונויות עבור כל ערך של המשתנה המסביר • יחס לינארי בין המשתנה המסביר למשתנה המגיב Smith, R. J. 2009. Use and misuse of the reduced major axis for line-fitting. American Journal of Physical Anthropology 140: 476-486.
איך נבדוק אם המודל מתנהג לפי ההנחות? ל-R יש model diagnostic functions שימושיות למדי שמאפשרות לנו להעריך (באופן גרפי, איכותי) עד כמה הנחות מודלים (בדגש על רגרסיה) מתקיימות https://www.youtube.com/watch?v=eTZ4VUZHzxw ראו גם: http://stat.ethz.ch/R-manual/R-patched/library/stats/html/plot.lm.html
מה עושים אם הנתונים לא מתנהגים לפי ההנחות? • מתעלמים ומקווים שהמבחן רובוסטי להפרה של הנחותיו: זה לא תמיד בלתי סביר כמו שזה נשמע • משתמשים במבחנים א-פרמטריים • משתמשים ב-generalized linear models (או בקיצור glm) בגדול זה אומר: • טרנספורמציות (ב-glmזה אומר שינוי link functions) • שינוי ב-error distribution(ב-glm) להתפלגויות לא נורמאליות • משתמשים במבחן לא לינארי • משתמשים ברנדומיזציות (ראו בהרצאות של רועי הולצמן)
מבחנים א-פרמטריים מבחנים א-פרמטריים לא מניחים שוויון שונויות או התפלגות נורמלית. הם מבוססים על rank. חסרונות: לא קיימים מבחנים למודלים מרובי predictors לעיתים קרובות ה-statistical power שלהם נמוך משל מבחן פרמטרי מקביל לא מאפשרים הערכת פרמטרים (שיפועים ונקודות חיתוך)
נראה לי ממש לא בסדר שבמצגת שלמה אין לי תמונות של חיות מבחנים א-פרמטריים מבחנים א-פרמטריים לא מניחים שוויון שונויות או התפלגות נורמלית. הם מבוססים על rank. חסרונות: לא קיימים מבחנים למודלים מרובי predictors לעיתים קרובות ה-statistical power שלהם נמוך משל מבחן פרמטרי מקביל לא מאפשרים הערכת פרמטרים (שיפועים ונקודות חיתוך)
כמה מבחנים א-פרמטריים שימושיים Orycteropusafer למצולם אין קשר להרצאה מבחן χ2הוא מבחן א-פרמטרי קולמוגורוב-סמירנוב הוא מבחן א-פרטרי להשוואת שתי התפלגויות שכיחויות (או השוואת התפלגות "שלנו" להתפלגות ידועה – דוגמת הנורמלית) Mann-Whitney U = Wilcoxon rank sumהוא מבחן א-פרטרי המחליף מבחן t Wilcoxon two-sample (=Wilcoxon signed-rank) testמחליף מבחן paired-t Kruskal-Wallisמחליף מבחן one-way ANOVA מבחןSpearmanומבחןKendall’s-tauמחליפים מבחן קורלציה