Transcript Model selection - TAU R Workshop 2015

‫סטטיסטיקה בסיסית והסקה‬
‫סטטיסטית ב‪R-‬‬
‫!!!‪Everything differs‬‬
‫"צפויים להימצא הבדלים בין ‪ x‬ל‪ "y-‬היא‬
‫אמירה טריוויאלית‬
‫הסטטיסטיקאי שבכם שואל "האם ההבדלים שנמצאו גדולים מהצפוי באקראי"‬
‫הביולוג שבכם שואל "למה ההבדלים הם לכיוון ובדרגה שמצאתי"‬
*‫מדדים לנטייה מרכזית‬
‫ממוצע‬.1
‫ הרמוני‬,‫ גיאומטרי‬,‫חשבוני‬
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
:‫דוגמא‬
Geometric mean: (x1*x2*…*xn)1/n
data<-c(2,3,4,5,6,7,8)
exp(mean(log(data)))
[1] 4.549163
:‫דוגמא‬
‫ ממוצע גיאומטרי‬.2
*‫מדדים לנטייה מרכזית‬
)mean( ‫ ממוצע‬.‫א‬.1
)median( ‫ חציון‬.‫ב‬
)mode( ‫ שכיח‬.‫ג‬
* = Moments of central tendency
data<-c(2,3,4,5,6,7,8)
median(data)
[1] 5
:‫דוגמא‬
‫מדדים לנטייה מרכזית‬
‫‪http://www.statmethods.net/management/functions.html‬‬
‫‪.1‬ממוצע‬
‫‪.2‬שונות*‬
‫הממוצע הוא מדד טוב יותר למה שקורה‬
‫באוכלוסיה\מדגם כשהשונות קטנה או‬
‫גדולה?‬
‫‪*Variance = Σ(xi-μ)2 / n‬‬
‫דוגמא‪:‬‬
‫)‪data<-c(2,3,4,5,6,7,8‬‬
‫)‪var(data‬‬
‫‪[1] 4.666667‬‬
‫מדדים לנטייה מרכזית‬
‫‪.1‬ממוצע‬
‫‪.2‬שונות‬
‫המומנט השני של נטייה מרכזית הוא מדד‬
‫לפיזור הנתונים מסביב למומנט הראשון‬
‫דוגמאות למומנט השני הן‪ ,‬למשל‪ ,‬השונות‪ ,‬סטיית התקן‪ ,‬שגיאת התקן‪ ,‬ה‪-‬‬
‫‪ ,coefficient of variation‬ורווח הסמך (של ‪ 99% ,95 ,90‬או מה שלא יהיה)‬
‫מדדים לנטייה מרכזית‬
#for:
data<-c(2,3,4,5,6,7,8)
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
:coefficient of variation
‫מדדים לנטייה מרכזית‬
‫‪.1‬ממוצע‬
‫‪.2‬שונות‬
‫‪.3‬הטייה (‪)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))
:‫שגיאת תקן של הטיה‬
‫מדדים לנטייה מרכזית‬
‫‪.1‬ממוצע‬
‫‪.2‬שונות‬
‫‪.3‬הטייה (‪)Skew‬‬
‫‪.4‬קורטוזיס (‪)Kurtosis‬‬
‫מדדים לנטייה מרכזית‬
)Kurtosis( ‫קורטוזיס‬.1
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‬ס"מ‬
‫השכר הממוצע הוא ‪( ₪ 9302‬למ"ס‪ ,‬נתוני מרץ ‪)2013‬‬
‫והשירות הצבאי הממוצע הוא (אולי) ‪ 24‬חודשים‬
‫‪ 2.05‬מטר‬
‫דוב ליאור‬
‫שירת בצה"ל חודש אחד‬
‫‪ ₪ 46,699‬בחודש‬
‫(הערכת חסר‪ ,‬ללא הטבות)‬
‫‪http://www.haaretz.co.il/1.2057452‬‬
‫‪Residuals‬‬
‫כשאנחנו עושים סטטיסטיקה אנחנו יוצרים מודלים של המציאות‬
‫כך שכאן‪ ,‬המודלים שלנו‪ 173 :‬ס"מ‪ 24 ,₪ 9302 ,‬חודשים‪ ,‬אינם‬
‫מוצלחים במיוחד‬
‫ה‪ Residual-‬הוא הכמות בה ערך מסויים רחוק מהניבוי של המודל‪ .‬כך‬
‫שליאור אליהו רחוק ‪ 32‬ס"מ מהמודל "ישראלי = ‪ ,"173‬ורחוק ‪ 29‬ס"מ‬
‫מהמודל המורכב יותר "גבר ישראלי = ‪ ,177‬אישה ישראלית = ‪."168‬‬
‫‪Residual = ₪37397‬‬
‫‪Residual = 32 cm‬‬
‫‪Residual = -23 month‬‬
‫‪IDF service‬‬
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
‫סטטיסטיקה תיאורית והסקה סטטיסטית‬
‫כשיש לנו נתונים בשלב הראשון כדאי שנתאר אותם‪ :‬נצייר‬
‫גרפים‪ ,‬נחשב ממוצע וכדומה‬
‫בהסקה סטטיסטית אנו בוחנים את התנהגותם של הנתונים שלנו‬
‫אל מול השערה (היפותזה) מסויימת‬
‫את ההשערה שלנו אנו יכולים להציג כמודל סטטיסטי‬
‫למשל‪:‬‬
‫• התפלגות גבהים היא נורמאלית‬
‫• מספר המינים הולך ועולה עם העליה בשטח‬
‫• מספר המינים הולך ועולה עם העליה בשטח על פי ‪power‬‬
‫‪ function‬שלו מעריך חזקה השווה ל‪0.25-‬‬
‫איזה מבחן נבחר?‬
‫זה משתנה בהתאם לטבע ה ‪=( response variable‬המשתנה‬
‫התלוי‪ ,‬זה שעל ציר ה‪ ,)y-‬ובעיקר לפי טבע ה‪predictor variables-‬‬
‫•אם ה‪ response variable-‬הוא "הצלחה או כשלון"‪ ,‬והשערת האפס‬
‫היא של שיוויון ביניהם‪ ,‬נשתמש במבחן בינומי (‪)binomial‬‬
‫•אם ה ‪ response variable‬שלנו הוא ספירות נשתמש לרוב במחני‬
‫חי‪-‬בריבוע או ‪.)log-likelihood=( G‬‬
‫• אבל לרוב ה ‪ response variable‬שלנו יהיה רציף (‪ 14‬מינים‪78 ,‬‬
‫פרטים‪87.5 ,‬מ"מ‪ 54 ,‬פעימות לדקה‪ 7.3 ,‬ביצים‪ 23 ,‬מעלות)‬
‫איזה מבחן נבחר?‬
‫מהו ה‪? response variable-‬‬
‫"הצלחה" או "כישלון"‬
‫(מצא את הגבינה\אידיוט)‬
‫מבחן בינומי‬
‫(‪)binomial‬‬
‫ספירות‬
‫(שכיחויות‪6 :‬‬
‫זכרים‪ 9 ,‬נקבות)‬
‫רציף‬
‫(‪ 14‬מינים‪ 78 ,‬פרטים‪,‬‬
‫‪87.5‬מ"מ‪ 54 ,‬חודשים‪,‬‬
‫‪ 7.3‬ביצים‪ 23 ,‬מעלות)‬
‫חי‪-‬בריבוע או ‪G‬‬
‫(=‪)log-likelihood‬‬
‫ראה בהמשך‬
Binomial test in R
.‫יש להגדיר את מספר ההצלחות מתוך גודל המדגם הכולל‬
)‫ (מובהק‬20 ‫ מתוך‬19 :2 ‫דוגמה‬
.)‫ (לא מובהק‬34 ‫ מתוך‬19 :1 ‫דוגמה‬
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
chisq.test
Data: insularity vs. diet:
habitat
island
island
island
mainland
mainland
mainland
diet
carnivore
herbivore
omnivore
carnivore
herbivore
omnivore
species#
488
43
177
1901
101
269
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
‫איזה מבחן נבחר?‬
‫אם ה ‪ response variable‬שלנו הוא רציף‬
‫נבחר מבחן לפי טבע ה‪predictor variables-‬‬
‫• אם ה‪ predictor variable-‬בדיד (אתר א'‪ ,‬אתר ב'‪ ,‬אתר ג';‬
‫זכר\נקבה; מין א'‪ ,‬מין ב‪ ,‬מין ג'; שטח עירוני\שטח טבעי;‬
‫טיפול א'‪ ,‬טיפול ב'‪ ,‬ביקורת) המבחן יהיה‬
‫)‪(analysis of variance‬‬
‫‪ANOVA‬‬
‫• אם ה ‪ predictor variable‬רציף (מעלות צלסיוס‪ ,‬כמות‬
‫מזון‪ ,‬קו רוחב‪ ,‬כמות משקעים‪ ,‬מספר מתחרים‪ ,‬אחוז כיסוי)‬
‫המבחן יהיה רגרסיה‬
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)
Species
Xenagama_zonura
Xenagama_zonura
Xenosaurus_grandis
Xenosaurus_grandis
Xenosaurus_newmanorum
Xenosaurus_newmanorum
Xenosaurus_penai
Xenosaurus_penai
Xenosaurus_platyceps
Xenosaurus_platyceps
Xenosaurus_rectocollaris
Xenosaurus_rectocollaris
Zonosaurus_anelanelany
Zonosaurus_anelanelany
Zootoca_vivipara
Zootoca_vivipara
Zygaspis_nigra
Zygaspis_nigra
Zygaspis_quadrifrons
Zygaspis_quadrifrons
size
79.7
85
120
133.0
118
126.0
105.8
112
106
121.0
95
111.0
86
93.0
65
75.0
230
240.0
195
227.0
Sex
female
male
male
female
male
female
female
male
male
female
male
female
male
female
male
female
male
female
male
female
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)
Species
Xenagama_zonura
Xenagama_zonura
Xenosaurus_grandis
Xenosaurus_grandis
Xenosaurus_newmanorum
Xenosaurus_newmanorum
Xenosaurus_penai
Xenosaurus_penai
Xenosaurus_platyceps
Xenosaurus_platyceps
Xenosaurus_rectocollaris
Xenosaurus_rectocollaris
Zonosaurus_anelanelany
Zonosaurus_anelanelany
Zootoca_vivipara
Zootoca_vivipara
Zygaspis_nigra
Zygaspis_nigra
Zygaspis_quadrifrons
Zygaspis_quadrifrons
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)
(Intercept)
Sexmale
Estimate
88.17
3.932
standard error
1.291
1.825
t
68.32
2.154
p value
<2e-16 ***
0.031 *
size
79.7
85
120
133.0
118
126.0
105.8
112
106
121.0
95
111.0
86
93.0
65
75.0
230
240.0
195
227.0
Sex
female
male
male
female
male
female
female
male
male
female
male
female
male
female
male
female
male
female
male
female
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)
Species
Xenagama_zonura
Xenagama_zonura
Xenosaurus_grandis
Xenosaurus_grandis
Xenosaurus_newmanorum
Xenosaurus_newmanorum
Xenosaurus_penai
Xenosaurus_penai
Xenosaurus_platyceps
Xenosaurus_platyceps
Xenosaurus_rectocollaris
Xenosaurus_rectocollaris
Zonosaurus_anelanelany
Zonosaurus_anelanelany
Zootoca_vivipara
Zootoca_vivipara
Zygaspis_nigra
Zygaspis_nigra
Zygaspis_quadrifrons
Zygaspis_quadrifrons
size
79.7
85
120
133.0
118
126.0
105.8
112
106
121.0
95
111.0
86
93.0
65
75.0
230
240.0
195
227.0
Sex
female
male
male
female
male
female
female
male
male
female
male
female
male
female
male
female
male
female
male
female
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)
female
88.17
male
92.10
ANOVA in R
aov
model<-aov(x~y)
species
Oligosoma_otagense
Oligosoma_polychroma
Oligosoma_smithi
Oligosoma_striatum
Oligosoma_suteri
Oligosoma_waimatense
Oligosoma_whitakeri
Oligosoma_zelandicum
Omanosaura_jayakari
Ophidiocephalus_taeniatus
ochel<-read.table("di.txt",header=T)
names(ochel)
[1] "species" "diet" "size"
model<-aov(size~diet,data=ochel)
summary(model)
diet
Residuals
Df
2
2959
Sum sq.
205.6
1428.9
Mean Sq. F value Pr(>F)
102.81 212.9 <2e-16 ***
0.48
diet
Omnivorous
Omnivorous
Omnivorous
Herbivorous
Carnivorous
Omnivorous
Herbivorous
Carnivorous
Omnivorous
Carnivorous
size
1.732917
1.002438
1.020078
0.948147
1.657096
1.645922
1.346954
0.89167
1.969343
0.743285
R-‫ ב‬ANOVA-‫ ל‬post-hoc ‫מבחן‬
TukeyHSD(model)
Fit: aov(formula = size ~ diet, data = ochel)
$diet
Herbivorous-Carnivorous
Omnivorous-Carnivorous
Omnivorous-Herbivorous
diff
lwr
upr
p adj
1.150905 1.011075 1.290736 0
0.329414 0.244294 0.414535 0
-0.82149 -0.97824 -0.66474 0
‫כל ההשוואות מובהקות‬
R-‫ ב‬ANOVA-‫ ל‬post-hoc ‫מבחן‬
model<-aov(SVL~Realm)
summary(model)
Df Sum Sq Mean Sq F value
Realm
11 5.270
0.479
9.3242
Pr(>F)
< 2.2e-16 ***
Residuals 4851 249.270 0.051
TukeyHSD(model)
Fit: aov(formula = SVL ~ Realm)
$Realm
diff
Asia-Africa
0.0006
Australia-Africa
0.0191
Caribbean-Africa
-0.0706
centralAmerica-Africa 0.0069
Europe-Africa
0.0686
lwr
-0.0351
-0.0205
-0.1184
-0.0401
-0.0306
upr
0.0363
0.0586
-0.0228
0.0539
0.1679
p adj
1.0000
0.9172
0.0001
1.0000
0.5051
‫אפס מחוץ‬
‫לרווח הסמך‬
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 ‫הקורלציה‬
latitude
35.91883
43.85467
8.622623
5.930
7.17859
19.85003
1.066133
35.46224
47.001
1.375278
28.86541
range_size
-0.50856
-0.40545
-0.21363
-0.09748
0.143959
0.167992
0.275788
0.302834
0.334081
0.347326
0.348311
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
Coefficients:
(Intercept)
latitude
Estimate Std. Error
3.26
0.0517
0.02
0.0024
t value
63.134
9.982
Pr(>|t|)
<2e-16
<2e-16
***
***
Residual standard error: 1.844 on 4910 degrees of freedom
Multiple R-squared: 0.01989, Adjusted R-squared: 0.01969 Fstatistic: 99.65 on 1 and 4910 DF, p-value: < 2.2e-16
‫‪ aov‬לעומת ‪lm‬‬
‫אולי במפתיע אפשר לבחון נתונים המתאימים ל‪ ANOVA-‬גם במבחן ‪.lm‬‬
‫במקרה כזה נקבל את כל המידע שנותן ה‪ summary-‬של ‪ lm‬במבחן רגרסיה‪,‬‬
‫כולל (חשוב!) ‪ ,parameter estimates‬שגיאות תקן‪ ,‬הבדלים בין פקטורים‬
‫וערכי‪ p-‬לכל קונטרסט (בין ‪ 2‬קטגוריות של המשתנה המסביר הקטגוריאלי)‬
lm ‫ לעומת‬aov
.lm ‫ גם במבחן‬ANOVA-‫אולי במפתיע אפשר לבחון נתונים המתאימים ל‬
)!‫ כולל (חשוב‬,‫ במבחן רגרסיה‬lm ‫ של‬summary-‫במקרה כזה נקבל את כל המידע שנותן ה‬
2 ‫ לכל קונטרסט (בין‬p-‫ הבדלים בין פקטורים וערכי‬,‫ שגיאות תקן‬,parameter estimates
)‫קטגוריות של המשתנה המסביר הקטגוריאלי‬
model2<-lm(size~diet,data=ochel)
summary(model2)
(Intercept)
dietHerbivorous
dietOmnivorous
Estimate
1.03832
1.15091
0.32941
Std. Error
0.01423
0.05963
0.0363
t value
72.97
19.3
9.075
Pr(>|t|)
<2e-16 ***
<2e-16 ***
<2e-16 ***
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
‫עוד בהמשך‬
‫‪Multiple predictors‬‬
‫מה לעשות‪ ,‬החיים מסובכים‪ .‬לפעמים מה‬
‫שמעניין אותנו מושפע מיותר מגורם אחד!‬
‫קצב ליבם של זוחלים‪ ,‬למשל‪ ,‬מושפע גם מגודל גופם‬
‫וגם מטמפרטורת הסביבה‪ ,‬וגם מהזמן והמהירות בהם‬
‫נעו לאחרונה‬
‫‪Smith, R. J. 1999. Statistics of sexual size dimorphism. Journal of Human Evolution 36: 423-459.‬‬
‫‪Multiple predictors‬‬
‫מה לעשות‪ ,‬החיים מסובכים‪ .‬לפעמים מה‬
‫שמעניין אותנו מושפע מיותר מגורם אחד!‬
‫קצב ליבם של זוחלים‪ ,‬למשל‪ ,‬מושפע גם מגודל גופם וגם‬
‫מטמפרטורת הסביבה‪ ,‬וגם מהזמן והמהירות בהם נעו לאחרונה‬
‫ניתן להסביר אם כן את המשתנה המעניין (קצב לב) אם יש לנו מידע על כל‬
‫המשתנים המסבירים‪.‬‬
‫ההנחה היא שכאשר אנו מכניסים את שלושתם למשוואה אנו רואים את‬
‫השפעתו של כל אחד כאשר שני האחרים "מוחזקים קבועים" (‪)held constant‬‬
‫‪Smith, R. J. 1999. Statistics of sexual size dimorphism. Journal of Human Evolution 36: 423-459.‬‬
‫‪Multiple predictors‬‬
‫ניתן להסביר אם כן את המשתנה המעניין (קצב לב) אם יש לנו מידע על כל‬
‫המשתנים המסבירים‪.‬‬
‫ההנחה היא שכאשר אנו מכניסים את שלושתם למשוואה אנו רואים את‬
‫השפעתו של כל אחד כאשר שני האחרים "מוחזקים קבועים" (‪)held constant‬‬
‫הנחה זו נכונה כאשר אין מתאם (גבוה) בין המשתנים‬
‫המסבירים לבין עצמם‬
‫‪Smith, R. J. 1999. Statistics of sexual size dimorphism. Journal of Human Evolution 36: 423-459.‬‬
‫איזה מבחן נבחר?‬
‫אם יש כמה ‪( predictor variables‬נאמר ארבעה) וכולם‬
‫קטגוריאלים המבחן יהיה ‪ANOVA‬‬
‫(נאמר ‪)4-way ANOVA‬‬
‫אם יש כמה ‪( predictor variables‬נאמר שבעה) וכולם רציפים‬
‫המבחן יהיה ‪Multiple Regression‬‬
‫איך כותבים מבחן עם כמה משתנים מסבירים?‬
‫משתמשים בפלוס (‪ )+‬בין המשתנים המסבירים‪.‬‬
‫)‪lm(y~a+b+c‬‬
‫למשל במודל שמנסה לחזות ציוני קורסים לפי כמה‬
‫למדנו‪ ,‬גיל המרצה‪ ,‬כמה התפללנו והאם מעתיקים‪:‬‬
‫‪model<-lm(Grade~days_studied+age_of_professor+number_of‬‬
‫)‪prayers+sent_sms_texts,data=grades‬‬
‫)‪summary(model‬‬
‫איזה מבחן נבחר?‬
‫אם יש כמה ‪( predictor variables‬לפחות ‪ – (2‬חלקם‬
‫(לפחות ‪ )1‬בדידים וחלקם (לפחות ‪ )1‬רציפים המבחן יהיה‬
‫‪ANCOVA‬‬
‫)‪(analysis of co-variance‬‬
‫איך זה נראה גרפית?‬
‫דוגמא‪ :‬מדדתי אורך שלוש שיניים בשועלים מצויים‪,‬‬
‫זכרים ונקבות‪ ,‬מכל תחום תפוצתם‬
‫‪ANOVA‬‬
‫‪Regression‬‬
‫‪ANCOVA‬‬
‫שני הגורמים‬
‫מובהקים‪ :‬יש הבדל‬
‫בין השיניים ובין‬
‫הזוויגים‬
‫‪Vulpes vulpes‬‬
‫‪p‬‬
‫‪F‬‬
‫‪MS‬‬
‫‪DF‬‬
‫‪SS‬‬
‫‪0.00‬‬
‫‪417468.9‬‬
‫‪304435.2‬‬
‫‪1‬‬
‫‪304435.2‬‬
‫‪Intercept‬‬
‫‪0.00‬‬
‫‪142.3‬‬
‫‪103.8‬‬
‫‪1‬‬
‫‪103.8‬‬
‫‪sex‬‬
‫‪0.00‬‬
‫‪20704.9‬‬
‫‪15098.9‬‬
‫‪2‬‬
‫‪30197.7‬‬
‫‪tooth‬‬
‫‪0.7‬‬
‫‪2276‬‬
‫‪1659.8‬‬
‫‪Error‬‬
‫ניב‬
‫ניב‬
‫עליון‬
‫שן שסע‬
‫תחתונה‬
‫שן שסע‬
‫עליונה‬
‫שן שסע‬
‫עליונה‬
‫שן שסע‬
‫תחתונה‬
‫איך זה נראה גרפית?‬
‫דוגמא‪ :‬אורך שלוש שיניים בשועלים‬
‫מצויים‪ ,‬כפונקציה של קו הרוחב‬
‫‪ANOVA‬‬
‫‪R-squared: 0.015, F = 32.83,‬‬
‫‪1 & 2161 DF; p < 0.0001‬‬
‫‪Regression‬‬
‫‪ANCOVA‬‬
‫יש כלל ברגמן‬
‫אבל קל לראות שהמודל* זוועתי‪ :‬ניבים‬
‫קטנים יותר משיני שסע‬
‫*קו הרגרסיה הוא מודל ליחס בין ה ‪ predictor‬ל ‪response‬‬
‫‪P‬‬
‫‪t‬‬
‫‪Std.‬‬
‫‪Error‬‬
‫‪Estimate‬‬
‫> ‪0.0001‬‬
‫‪25.25‬‬
‫‪0.378‬‬
‫‪9.554‬‬
‫‪Intercept‬‬
‫> ‪0.0001‬‬
‫‪5.73‬‬
‫‪0.008‬‬
‫‪0.044‬‬
‫‪Latitude‬‬
‫‪ ANCOVA‬גרפית‬
‫קל להבין את זה גרפית‪ :‬בדוגמה משתנה בדיד אחד עם ‪2‬‬
‫רמות (מקווקו ומלא)‪ ,‬ומשתנה רציף אחד (לאורך ציר ה ‪)X‬‬
‫‪response‬‬
‫‪Continuous predictor‬‬
‫‪response‬‬
‫רציף מובהק‪ ,‬בדיד‬
‫לא‬
‫בדיד מובהק‬
‫רציף לא‬
‫‪response‬‬
‫‪c.‬‬
‫‪b.‬‬
‫‪a.‬‬
‫‪Null hypothesis‬‬
‫‪Continuous predictor‬‬
‫‪Continuous predictor‬‬
‫‪d.‬‬
‫שניהם מובהקים‬
‫‪response‬‬
‫‪Continuous predictor‬‬
‫ותודה לדניאל על הגרפים‬
‫איך זה נראה גרפית?‬
‫דוגמא‪ :‬אורך שיניים בשועלים מצויים‪ ,‬כפונקציה של קו‬
‫הרוחב (ציר ה‪ )X-‬הזוויג (צבע) ובאיזו שן מדובר (צורה)‬
‫‪ANOVA‬‬
‫‪Regression‬‬
‫‪p‬‬
‫‪F‬‬
‫‪MS‬‬
‫‪SS‬‬
‫‪Df‬‬
‫‪factor‬‬
‫> ‪0.0001‬‬
‫‪191.7‬‬
‫‪99‬‬
‫‪99‬‬
‫‪1‬‬
‫‪sex‬‬
‫‪2‬‬
‫‪tooth‬‬
‫‪436.1‬‬
‫‪436.1‬‬
‫‪1‬‬
‫‪Latitude‬‬
‫‪0.5‬‬
‫‪1114.2‬‬
‫‪2158‬‬
‫‪Residuals‬‬
‫> ‪28665.3 14332.7 27758.79 0.0001‬‬
‫> ‪0.0001‬‬
‫‪844.69‬‬
‫‪ANCOVA‬‬
‫המודל הזה‬
‫מסביר ‪96.3%‬‬
‫מהשונות‬
‫העברנו פה שישה‬
‫קווי רגרסיה‬
‫מקבילים‬
‫כל הגורמים מובהקים‬
R-‫ ב‬ANCOVA ‫קריאת תוצאות מודל‬
‫דוגמא ללא‬
‫אינטראקציות‬
Response = intercept + a for level 1 of the 1st categorical
predictor variable or + b for level 2 of the 1st + c for level 1 of
the 2nd categorical predictor or d for level 2… +k*(value of the
continuous predictor variable) + error
Tooth / sex
‫ אם נחזור לשועלים‬,‫למשל‬
factor
Estimate
Std.
Error
Intercept
(tooth c)
4.342
0.078
55.92
tooth_m
8.485
0.038
224.21 0.0001>
tooth_p
6.617
0.038
174.84 0.0001>
sex_male
0.391
0.031
12.56
0.0001>
Latitude
0.043
0.001
29.06
0.0001>
Latitude
t
p
0.0001>
:‫) על פי המודל הוא‬32 ‫ בזכרים בתל אביב (קו רוחב‬P ‫כך שהאורך של שן‬
12.726 = 4.342+6.617+0.391+32*0.043
‫איך קוראים תוצאות של ‪ lm‬ב‪R-‬‬
‫כשהמשתנה המסביר קטגוריאלי‪ R ,‬משווה את כל הפקטורים‬
‫ל‪ intercept-‬של הפקטור הראשון בסדר אלפביתי‬
‫)|‪Pr(>|t‬‬
‫*** ‪<2e-16‬‬
‫*** ‪<2e-16‬‬
‫*** ‪<2e-16‬‬
‫‪t value‬‬
‫‪72.97‬‬
‫‪19.3‬‬
‫‪9.075‬‬
‫‪Std. Error‬‬
‫‪0.01423‬‬
‫‪0.05963‬‬
‫‪0.0363‬‬
‫‪Estimate‬‬
‫‪1.03832‬‬
‫‪1.15091‬‬
‫‪0.32941‬‬
‫)‪(Intercept‬‬
‫‪dietHerbivorous‬‬
‫‪dietOmnivorous‬‬
‫כאן המשתנה המסביר הוא דיאטה ואלפביתית הקטגוריה הראשונה היא ‪Carnivore‬‬
‫ההבדל בין קרניבורים להרביבורים (‪ )t = 19.3‬ואומניבורים (‪ )t = 9.075‬מובהק‬
‫כך שגודלו הממוצע של קרניבור הוא ‪ 1.038 ,‬וזה של אומניבור‪1.367 = 1.038+0.329 :‬‬
‫איך קוראים תוצאות של ‪ lm‬ב‪R-‬‬
‫כשהמשתנה המסביר רציף‪ R ,‬מדווח עבורו את השיפוע עם‬
‫שגיאת התקן שלו וערכי ה‪ p-‬וה‪ t-‬שלו‬
‫)‪rapoport<-read.table("rang.txt",header=T‬‬
‫)‪model<-lm(range~mass+latitude,data=rapoport‬‬
‫)‪summary(model‬‬
‫)|‪Pr(>|t‬‬
‫*** ‪<2e-16‬‬
‫*** ‪<2e-16‬‬
‫*** ‪<2e-16‬‬
‫‪t value‬‬
‫‪47.985‬‬
‫‪8.817‬‬
‫‪9.57‬‬
‫‪Std. Error‬‬
‫‪0.061833‬‬
‫‪0.036929‬‬
‫‪0.002365‬‬
‫‪Estimate‬‬
‫‪2.967075‬‬
‫‪0.325585‬‬
‫‪0.02263‬‬
‫)‪(Intercept‬‬
‫‪mass‬‬
‫‪latitude‬‬
‫כאן המשתנים המסבירים (את גודל טווח התפוצה) הם מסה וקו רוחב‬
‫ההשפעה של שניהם מובהקת (‪ ,t = 8.82‬ו‪)p<<0.05 ,t = 9.57-‬‬
‫כך שגודלו תחום התפוצה עולה ב‪( 0.02-‬היחידות הן לוג קמ"ר) בכל עליה של מעלה‪ ,‬וב‪-‬‬
‫‪ 0.326‬יחידות כל עליה של יחידת מסה (לוג גרם‪ ,‬כלומר יחידה אחת היא הכפלה בעשר)‬
‫איך קוראים תוצאות של ‪ lm‬ב‪R-‬‬
‫כשהמשתנה המסביר רציף‪ R ,‬מדווח עבורו את השיפוע עם‬
‫שגיאת התקן שלו וערכי ה‪ p-‬וה‪ t-‬שלו‬
‫)|‪Pr(>|t‬‬
‫*** ‪<2e-16‬‬
‫*** ‪<2e-16‬‬
‫*** ‪<2e-16‬‬
‫‪t value‬‬
‫‪47.985‬‬
‫‪8.817‬‬
‫‪9.57‬‬
‫‪Std. Error‬‬
‫‪0.061833‬‬
‫‪0.036929‬‬
‫‪0.002365‬‬
‫‪Estimate‬‬
‫‪2.967075‬‬
‫‪0.325585‬‬
‫‪0.02263‬‬
‫)‪(Intercept‬‬
‫‪mass‬‬
‫‪latitude‬‬
‫כך שגודל תחום התפוצה של לטאה במסה של ‪ 100‬גרם (‪ )log(100) = 2‬בקוסטה ריקה‬
‫(קו רוחב ‪ )10‬יהיה‪:‬‬
‫‪Intercept+slope*mass+slope*latitude‬‬
‫‪2.967+0.3255(slope)*2(mass)+0.023(slope)*10latitude = 3.84‬‬
R-‫ ב‬lm ‫איך קוראים תוצאות של‬
‫ כשיש גם המשתנים מסבירים רציפים וגם‬ANCOVA‫ב‬
‫ לאחרונים ושיפוע‬intercept ‫ מדווח‬R ,‫קטגוריאלים‬
‫ המתאימים‬t-‫ וה‬p-‫ עם שגיאות התקן שלו וערכי ה‬,‫לראשונים‬
model<-lm(range~islands+mass+latitude,data=rapoport)
summary(model)
(Intercept)
islandsMIE
islandsSIE
mass
latitude
Estimate
3.777624
-1.232342
-2.030331
0.266193
-0.010585
Std. Error
0.059828
0.082848
0.059402
0.033047
0.002139
t value
63.142
-14.875
-34.179
8.055
-4.95
Pr(>|t|)
< 2e-16 ***
< 2e-16 ***
< 2e-16 ***
9.91E-16 ***
7.68E-07 ***
‫) מינים שאנדמיים לאי יחיד‬continental( ‫ מיני יבשת‬:‫ קטגוריות‬3 ‫כאן יש‬
)multiple island endemics( ‫) או מכמה איים‬single island endemics(
R-‫ ב‬lm ‫איך קוראים תוצאות של‬
‫ גם משתנים מסבירים רציפים וגם קטגוריאלים‬:ANCOVA
model<-lm(range~islands+mass+latitude,data=rapoport)
summary(model)
(Intercept)
islandsMIE
islandsSIE
mass
latitude
Estimate
3.777624
-1.232342
-2.030331
0.266193
-0.010585
Std. Error
0.059828
0.082848
0.059402
0.033047
0.002139
t value
63.142
-14.875
-34.179
8.055
-4.95
Pr(>|t|)
< 2e-16 ***
< 2e-16 ***
< 2e-16 ***
9.91E-16 ***
7.68E-07 ***
)MIE( ‫) או כמה איים‬SIE( ‫ אי יחיד‬,)continental( ‫ יבשת‬:‫ קטגוריות‬3
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
‫ וכו' של המודל בכללותו‬F ,‫ בריבוע‬R ‫אלה ערכי‬
‫איך קוראים תוצאות של ‪ lm‬ב‪R-‬‬
‫)‪model<-lm(range~islands+mass+latitude,data=rapoport‬‬
‫)‪summary(model‬‬
‫)|‪Pr(>|t‬‬
‫*** ‪< 2e-16‬‬
‫*** ‪< 2e-16‬‬
‫*** ‪< 2e-16‬‬
‫*** ‪9.91E-16‬‬
‫*** ‪7.68E-07‬‬
‫‪t value‬‬
‫‪63.142‬‬
‫‪-14.875‬‬
‫‪-34.179‬‬
‫‪8.055‬‬
‫‪-4.95‬‬
‫‪Std. Error‬‬
‫‪0.059828‬‬
‫‪0.082848‬‬
‫‪0.059402‬‬
‫‪0.033047‬‬
‫‪0.002139‬‬
‫‪Estimate‬‬
‫‪3.777624‬‬
‫‪-1.232342‬‬
‫‪-2.030331‬‬
‫‪0.266193‬‬
‫‪-0.010585‬‬
‫)‪(Intercept‬‬
‫‪islandsMIE‬‬
‫‪islandsSIE‬‬
‫‪mass‬‬
‫‪latitude‬‬
‫כיוון שאלפביתית ‪ continental<MIE<SIE‬החיתוך (‪ )intercept‬שלנו הוא‬
‫עבור הקטגוריה הראשונה‪ :‬מיני יבשת‪ .‬כך שטווח התפוצה של מיני יבשת גדול‬
‫משל ‪( MIE‬שימו לב‪ :‬הבדל שלילי!) ושל ‪ – SIE‬וכמו כן טווח התפוצה עולה‬
‫עם המסה‪ ,‬אך יורד עם העליה בקו הרוחב (שיפוע שלילי)‪.‬‬
‫ושוב – כל הגורמים מובהקים‬
‫‪relevel‬‬
‫)‪model<-lm(range~islands+mass+latitude,data=rapoport‬‬
‫)‪summary(model‬‬
‫)|‪Pr(>|t‬‬
‫*** ‪< 2e-16‬‬
‫*** ‪< 2e-16‬‬
‫*** ‪< 2e-16‬‬
‫*** ‪9.91E-16‬‬
‫*** ‪7.68E-07‬‬
‫‪t value‬‬
‫‪63.142‬‬
‫‪-14.875‬‬
‫‪-34.179‬‬
‫‪8.055‬‬
‫‪-4.95‬‬
‫‪Std. Error‬‬
‫‪0.059828‬‬
‫‪0.082848‬‬
‫‪0.059402‬‬
‫‪0.033047‬‬
‫‪0.002139‬‬
‫‪Estimate‬‬
‫‪3.777624‬‬
‫‪-1.232342‬‬
‫‪-2.030331‬‬
‫‪0.266193‬‬
‫‪-0.010585‬‬
‫)‪(Intercept‬‬
‫‪islandsMIE‬‬
‫‪islandsSIE‬‬
‫‪mass‬‬
‫‪latitude‬‬
‫אבל בגורמים הקטגוריאלים יש לנו בעיה‪ R :‬מחשב רק את ההבדל בין כל גורם‬
‫לגורם הראשון באלפבית‪ .‬כאן בין יבשת לשתי קטגוריות האיים‪ .‬אבל הוא לא‬
‫מחשב ולא מדווח על ההבדלים בין שתי קטגוריות האיים‬
‫וחוץ מזה‪ ,‬הוא לא נותן לנו עבורם שגיאות תקן והבדלים מאפס‪ ,‬אלא רק‬
‫הבדלים מהיבשת ושגיאת תקן של המבחן הזה (לא את שגיאת התקן של‬
‫המשתנה עצמו במודל)‬
Relevel (2)
model<-lm(range~islands+mass+latitude,data=rapoport)
summary(model)
‫ אליו הוא ישווה את‬,‫ מי יהיה הגורם הראשון‬R-‫ להגדיר ל‬:‫אפשר להתחכם‬
:relevel ‫ על ידי הפקודה‬,‫האחרים‬
model2<lm(range~relevel(islands,”SIE”)+mass+latitude,data=r
apoport)
summary(model2)
‫או‬
model3<lm(range~relevel(islands,”MIE”)+mass+latitude,data=r
apoport)
summary(model3)
Relevel (3)
‫ אליו הוא ישווה את‬,‫ מי יהיה הגורם הראשון‬R-‫ להגדיר ל‬:‫אפשר להתחכם‬
:relevel ‫ על ידי הפקודה‬,‫האחרים‬
model3<lm(range~relevel(islands,”MIE”)+mass+latitude,data=rapoport)
summary(model3)
Estimate Std. Error
t value
Pr(>|t|)
2.545
0.092
27.74
<2e-16
***
relevel(islands, MIE)continental 1.232
relevel(islands, MIE)SIE
-0.8
0.083
14.88
<2e-16
***
-8.54
8.055
<2e-16
9.91E-16
***
***
-4.95
7.68E-07
***
(Intercept)
mass
0.266
0.093
0.033
latitude
-0.01
0.002
:‫שימו לב שהפרמטרים של המודל הכללי נשארו זהים‬
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
‫בחירת מבחן מתאים‬
:‫) מתקיימים‬residuals-‫ התפלגות נורמלית של ה‬,‫אם הנחותיהם של מבחנים פרמטריים (שונות דומה‬
Predictor
Response
test
In R
Categorical
Success/failure
binomial
binom.test
Categorical
Counts
Chi-square/G
chisq.test
Categorical
continuous
ANOVA*
aov
continuous
continuous
Regression/correlation lm
continuous
Categorical/counts
Chi-square/ANOVA
lm
Categorical, multiple
predictors
continuous
Multi-way ANOVA
aov
continuous, multiple
predictors
continuous
Multiple regression
lm
multiple predictors,
Both categorical &
continuous
continuous
ANCOVA
lm
*t-test if there are only 2 categories
‫אינטראקציות‬
‫קל להבין את זה גרפית‪ :‬בדוגמה משתנה בדיד אחד עם ‪ 2‬רמות‬
‫(מקווקו ומלא)‪ ,‬ומשתנה רציף אחד (לאורך ציר ה ‪)X‬‬
‫‪response‬‬
‫‪Continuous predictor‬‬
‫‪d.‬‬
‫‪response‬‬
‫‪response‬‬
‫‪Continuous predictor‬‬
‫‪b.‬‬
‫שניהם מובהקים‪,‬‬
‫אין אינטראקציה‬
‫‪Continuous predictor‬‬
‫‪response‬‬
‫‪Continuous predictor‬‬
‫שניהם מובהקים‪,‬‬
‫יש אינטראקציה‬
‫‪Null hypothesis‬‬
‫‪Continuous predictor‬‬
‫‪Continuous predictor‬‬
‫‪f.‬‬
‫רציף מובהק בדיד‬
‫לא‪ ,‬אין אינטראקציה‬
‫‪response‬‬
‫רציף מובהק בדיד‬
‫לא‪ ,‬יש אינטראקציה‬
‫בדיד מובהק‬
‫רציף לא‬
‫‪response‬‬
‫‪e.‬‬
‫‪c.‬‬
‫‪a.‬‬
‫אינטראקציות ב‪R-‬‬
‫)‪lm(y~a*b‬‬
‫)‪lm(y~a+b+c+a:b‬‬
‫משתמשים בפלוס (‪ )+‬בין המשתנים המסבירים‪ .‬עבור אינטראקציה‬
‫משתמשים בנקודותיים‪ .‬אם רוצים לבחון גם ‪ main effect‬וגם‬
‫אינטראקציה משתמשים בכוכבית למשל במודל שמנסה לחזות ציוני‬
‫קורסים לפי כמה למדנו‪ ,‬גיל המרצה‪ ,‬כמה התפללנו והאם מעתיקים‪:‬‬
‫<‪model‬‬‫_‪lm(Grade~days_studied+age_of_professor+number_of_prayers*sent‬‬
‫)‪sms_texts+age_of_professor:number_of_prayers‬‬
‫כאן ביקשנו גם שתי אינטראקציות‪ :‬בין תפילות להעתקות וגיל לתפילה‬
‫חשוב‪:‬‬
‫הנחת היסוד של מבחנים מרובי ‪ predictors‬היא‬
‫שאין מתאם בין ‪ predictors‬שונים‬
‫מתאם גבוה בין זוג ‪ predictor variables‬קרוי גם ‪,multi-co-linearity‬‬
‫ולעיתים מבוטא על ידי ‪ )1-R2( tolerance‬או על ידי הרציפרוקלי שלו‬
‫‪)VIF = 1/tolerance( Variance Inflation Factors‬‬
‫אם ‪ multicollinearity‬חזק קיים אזי המודל לא יהיה יציב‪ ,‬והערכת‬
‫הפרמטרים עשויה להיות לא נכונה‬
‫אל תכניסו משתנים מנבאים המצויים‬
‫במתאם גבוה בינם לבין עצמם!‬
‫‪Model selection‬‬
‫תמיד‪ ,‬תמיד‪ ,‬תמיד‪ ,‬ככל שנוסיף יותר ‪predictor‬‬
‫‪ variables‬נסביר יותר מהשונות‬
‫היחס יהיה מונוטוני – וטריוויאלי‪ :‬במקרה הרע ביותר‬
‫ה‪ Parameter estimate -‬של משתנים נוספים יהיו אפס‬
‫(למשל מספר מינים = ‪*87+22.5‬קו הרוחב ‪*0 +‬מספר‬
‫המנדטים של המפלגות הדתיות באותו אזור)‬
‫אבל ה ‪ parameter estimate‬לא יהיה אף פעם בדיוק אפס – הוא‬
‫פשוט יהיה נמוך מאוד – נאמר נוסף מין על כל ‪ 5120‬מנדטים‬
‫שנוספים לש"ס‪ ,‬או נגרע מין על כל ‪ 974‬מנדטים הנוספים ליהדות‬
‫התורה‬
‫ה ‪ R2‬שלנו יעלה מ ‪ 0.45‬ל ‪ – 0.45007‬זה שווה לנו?‬
‫‪Model selection‬‬
‫למעשה בכל שאלה סטטיסטית אפשר להסביר‬
‫‪ 100%‬מהשונות באמצעות מספר משתנים השווה‬
‫למספר התצפיות‬
‫רוצים דוגמא?‬
‫מה הגובה שלכם?‬
‫אבל‪ ,‬מהי יכולת הניבוי של‬
‫המודל הזה לנתון הבא???‬
‫‪http://en.wikipedia.org/wiki/Overfitting‬‬
‫‪Model selection‬‬
‫ככל שנוסיף יותר ‪ predictor variables‬נסביר‬
‫יותר מהשונות‬
‫המטרה שלנו כמדענים היא להסביר את מקסימום‬
‫התופעות בעזרת מינימום משתנים‬
‫על תער אוקאם שמעתם?‬
‫כך שאם יש לנו משתנים מסבירים רבים מאוד נרצה‬
‫לדעת אילו מהם מוסיפים כל כך מעט שונות מוסברת‪,‬‬
‫שלא שווה לסבך בגללם את החיים‬
‫‪Model selection‬‬
‫ניתן פשוט לבחון אילו משתנים במודל מובהקים‬
‫‪Backwards (stepwise) elimination‬‬
‫‪ .1‬על פי ערכי ‪p‬‬
‫נתחיל במודל המורכב ביותר‪ ,‬ולמחוק כל פעם את המשתנה‬
‫(או האינטראקציה) שלו מיוחס ה‪ p-‬הגבוה ביותר – עד שכל‬
‫ערכי ה ‪ p‬קטנים מ ‪( 0.05‬או ערך סף אחר)‪ .‬המודל איתו‬
‫נשארנו ייקרא ‪MAM = minimum adequate model‬‬
‫דוגמא‪ :‬מנסים להסביר גודל תטולה בודדת בלטאות (‪)response variable: clutch‬‬
‫באמצעות נתונים על גודל גופן (‪ ,)SVL‬טמפרטורת סביבה המועדפת עליהן (‪,)temp‬‬
‫הגובה החציוני מעל פני הים בו הן חיות (‪ ,)elevation‬ומספר התטולות שהן‬
‫מטילות בשנה (‪)broods‬‬
‫‪Model selection‬‬
‫ניתן פשוט לבחון אילו משתנים במודל מובהקים‬
‫נתחיל במודל המורכב ביותר‪:‬‬
‫)‪model1<-lm(clutch~SVL+temp+elevation+Broods‬‬
‫‪p‬‬
‫‪t‬‬
‫‪Estimate se‬‬
‫‪<0.0001‬‬
‫‪-6.272‬‬
‫‪3.878‬‬
‫‪-24.32‬‬
‫‪Intercept‬‬
‫‪<0.0001‬‬
‫‪10.465‬‬
‫‪1.34‬‬
‫‪14.02‬‬
‫‪SVL‬‬
‫‪0.808‬‬
‫‪0.244‬‬
‫‪0.0809‬‬
‫‪0.0197‬‬
‫‪temp‬‬
‫‪<0.0001‬‬
‫‪4.103‬‬
‫‪0.0005‬‬
‫‪0.0021‬‬
‫‪elevation‬‬
‫‪0.362‬‬
‫‪-0.913‬‬
‫‪0.1069‬‬
‫‪-0.0976‬‬
‫‪Broods‬‬
‫דוגמא‪ :‬מנסים להסביר גודל תטולה בודדת בלטאות באמצעות נתונים על גודל גופן (‪ ,)SVL‬טמפרטורת סביבה המועדפת‬
‫עליהן (‪ ,)temp‬הגובה החציוני מעל פני הים בו הן חיות (‪ ,)elevation‬ומספר התטולות שהן מטילות בשנה (‪)broods‬‬
‫‪Model selection‬‬
‫ניתן פשוט לבחון אילו משתנים במודל מובהקים‬
‫נמחק את הטמפרטורה ונחשב מחדש‪:‬‬
‫)‪model2<-lm(clutch~SVL+elevation+Broods‬‬
‫‪p‬‬
‫‪t‬‬
‫‪Estimate se‬‬
‫‪<0.0001‬‬
‫‪-8.038‬‬
‫‪2.95‬‬
‫‪-23.71‬‬
‫‪Intercept‬‬
‫‪<0.0001‬‬
‫‪10.521‬‬
‫‪1.335‬‬
‫‪14.04‬‬
‫‪SVL‬‬
‫‪<0.0001‬‬
‫‪4.103‬‬
‫‪0.0005‬‬
‫‪0.0021‬‬
‫‪elevation‬‬
‫‪0.341‬‬
‫‪-0.953‬‬
‫‪0.1059‬‬
‫‪-0.1009‬‬
‫‪Broods‬‬
‫דוגמא‪ :‬מנסים להסביר גודל תטולה בודדת בלטאות באמצעות נתונים על גודל גופן (‪ ,)SVL‬טמפרטורת סביבה המועדפת‬
‫עליהן (‪ ,)temp‬הגובה החציוני מעל פני הים בו הן חיות (‪ ,)elevation‬ומספר התטולות שהן מטילות בשנה (‪)broods‬‬
‫‪Model selection‬‬
‫ניתן פשוט לבחון אילו משתנים במודל מובהקים‬
‫נמחק את מספר התטולות ונחשב מחדש‪:‬‬
‫)‪model3<-lm(clutch~SVL+elevation‬‬
‫כל‬
‫המשתנים‬
‫מובהקים‪,‬‬
‫עצור!‬
‫‪t‬‬
‫‪Estimate se‬‬
‫‪p‬‬
‫‪-12.433 <0.0001‬‬
‫‪-23.420 1.884‬‬
‫‪Intercept‬‬
‫‪<0.0001‬‬
‫‪15.151‬‬
‫‪0.9019‬‬
‫‪13.660‬‬
‫‪SVL‬‬
‫‪<0.0001‬‬
‫‪5.491‬‬
‫‪0.0003‬‬
‫‪0.0018‬‬
‫‪elevation‬‬
‫‪Model3 = MAM‬‬
‫גודל תטולה הוא פונקציה של גודל גוף ושל הגובה מעל פני הים‪ ,‬וזהו‬
‫‪Model selection‬‬
‫ניתן פשוט לבחון אילו משתנים במודל מובהקים‬
‫‪forward addition‬‬
‫אפשר להתחיל במודל הפשוט ביותר‪ ,‬ולהוסיף כל פעם משתנה‬
‫נוסף‪ ,‬ולהשאיר אותו אם ה‪ p-‬עבורו קטן מ‪( 0.05-‬או ערך סף אחר)‪.‬‬
‫)‪model1a<-lm(clutch~SVL‬‬
‫)‪Model2a<- lm(clutch~SVL+elevation‬‬
‫)‪Model3a<- lm(clutch~SVL+elevation+Broods‬‬
‫ברגע שנגיע למודל שמכיל גורמים לא מובהקים (‪ model3a‬בדוגמה שלנו)‬
‫נעצור ונבחר את המודל הקודם (‪ model2a‬בדוגמה שלנו) כ‪MAM-‬‬
‫‪Model selection‬‬
‫ניתן פשוט לבחון אילו משתנים במודל מובהקים‬
‫‪forward addition‬‬
‫מתחילים במודל הפשוט ביותר‪ ,‬ומוסיפים כל פעם משתנה נוסף‪,‬‬
‫ולהשאיר אותו אם ה‪ p-‬עבורו קטן מ‪( 0.05-‬או ערך סף אחר)‪.‬‬
‫שימו לב‪ :‬לא כל הצירופים האפשריים בין פרמטרים (והאינטראקציות ביניהם)‬
‫משמשים לא ב‪ forward addition-‬ולא ב‪ backwards elimination-‬כאן‪ ,‬ויכול‬
‫להיות שהצירוף "הכי נכון" לא נבדק‪.‬‬
‫מצד שני מספר המודלים האפשרי עולה בחזקה של מספר הפרמטרים בהם‬
‫משתמשים‪ ,‬כך שבחינת כל המודלים האפשריים לא מעשית אם יש לנו הרבה מאוד‬
‫פרמטרים – אלא אם כן אנו יודעים לתכנת היטב ויש לנו מחשב חזק והרבה זמן‪...‬‬
‫‪Model selection‬‬
‫‪Akaike Information Criterion‬‬
‫‪Hirotsugu Akaike‬‬
‫דרך חלופית לבחירת מודלים‬
‫השוואה בין מודלים על פי ‪ 2‬פרמטרים‪ :‬כמה המודל "טוב" (רמת הדיוק בה‬
‫מתוארת המציאות) לעומת כמה הוא מורכב (כמה פרמטרים הערכנו)‬
‫)‪AIC = 2k-2ln(L‬‬
‫כאשר ‪ K‬הוא מספר הפרמטרים ו‪ L-‬הוא ה ‪ maximum likelihood‬של‬
‫המודל (שבלי להכנס לפירוט מיותר מבטא במקרה הזה את ה ‪residual‬‬
‫‪ – sum of squares‬ככל שהוא קטן יותר המודל טוב יותר‪ ,‬וניתן לכתוב‬
‫[)‪)AIC = 2k+n[ln(RSS‬‬
‫ככל שערך ה ‪ AIC‬נמוך יותר המודל טוב יותר‬
‫‪http://en.wikipedia.org/wiki/Residual_sum_of_squares‬‬
‫‪http://en.wikipedia.org/wiki/Akaike_information_criterion‬‬
Model selection
Akaike Information Criterion
Hirotsugu Akaike
AIC = 2k-2ln(L)
‫שימו לב שהתמיכה במודל חזקה יותר ככל‬
:‫ נמוך יותר‬AIC ‫שערך ה‬
AIC rewards descriptive accuracy via the maximum
likelihood (High L), and penalizes lack of parsimony
according to the number of free parameters (high K)
AIC(model1,model2,model3)
‫‪Model selection‬‬
‫נחזור ללטאות‬
‫)‪model1<-lm(clutch~SVL+temp+elevation+Broods‬‬
‫)‪model2<-lm(clutch~SVL+elevation+Broods‬‬
‫)‪model3<-lm(clutch~SVL+elevation‬‬
‫)‪AIC(model1,model2,model3‬‬
‫‪AIC‬‬
‫‪df‬‬
‫‪1773.50‬‬
‫‪model1 6‬‬
‫‪1771.56‬‬
‫‪model2 5‬‬
‫‪1770.48‬‬
‫‪model3 4‬‬
‫שוב מודל ‪ 3‬הוא הטוב‬
‫ביותר (יש לו ה ‪AIC‬‬
‫‪ score‬הנמוך ביותר)‬
‫דוגמא‪ :‬מנסים להסביר גודל תטולה בודדת בלטאות באמצעות נתונים על גודל גופן (‪ ,)SVL‬טמפרטורת סביבה המועדפת עליהן‬
‫(‪ ,)temp‬הגובה החציוני מעל פני הים בו הן חיות (‪ ,)elevation‬ומספר התטולות שהן מטילות בשנה (‪)broods‬‬
‫‪Akaike Information Criterion‬‬
‫)‪AIC(model1,model2,model3‬‬
‫לא מאפשר לבחון כמה טוב הוא מבחן בודד‪ ,‬אלא רק להשוות‬
‫בין מספר מבחנים הנשענים על אותם נתונים‬
‫התוצאה (‪ )score‬של ‪ AIC‬חסרת משמעות בפני עצמה‪ :‬לא‬
‫ניתן להשוות ‪ AIC‬בין מבחנים של שאלות שונות או‬
‫שמתבססים על נתונים שונים כמו שניתן להגיד ש‬
‫‪0.05=0.05=0.05‬‬
‫בנוסף‪ ,‬כלל האצבע אומר שלא ניתן לומר שהבדלים‬
‫בערכי ‪ AIC‬של פחות מ‪ 2-‬אינם מאפשרים לומר‬
‫איזה מודל טוב יותר‬
‫‪Akaike Information Criterion‬‬
‫כלל האצבע אומר שלא ניתן לומר שהבדלים בערכי ‪ AIC‬של‬
‫פחות מ‪ 2-‬אינם מאפשרים לומר איזה מודל טוב יותר‬
‫)‪AIC(model1,model2,model3‬‬
‫‪∆AIC‬‬
‫‪AIC‬‬
‫‪0‬‬
‫‪1770.48‬‬
‫‪model3 4‬‬
‫‪1.12‬‬
‫‪1771.56‬‬
‫‪model2 5‬‬
‫‪3.02‬‬
‫‪1773.50‬‬
‫‪model1 6‬‬
‫‪df‬‬
‫נסדר את המודלים על פי ערכי ‪ AIC‬מהנמוך (הכי טוב) לגבוה‪ ,‬ונחשב‬
‫עבור כל אחד את הפרש ה ‪ AIC score‬מהמודל עם ה ‪ score‬הנמוך‬
‫ביותר – לקבלת ערך ה ‪ ∆AIC‬של כל מודל‪ .‬בדוגמה הזו לא ניתן לומר‬
‫שמודל ‪ 3‬עדיף על מודל ‪ 2‬כיוון שההפרש ב ‪ AIC‬ביניהם קטן מ ‪2‬‬
‫ וריאציות‬- AIC
AICc = AIC+(2k*[k+1]/[n-k-1])
AICc .1
‫ עבור מדגמים קטנים‬AIC ‫תיקון ל‬
)http://cran.r-project.org/web/packages/AICcmodavg/AICcmodavg.pdf R ‫(חבילת‬
AIC weights .2
“Akaike weights are used in model averaging. They represent the relative
likelihood of a model. To calculate them, for each model first calculate the
relative likelihood of the model, which is just exp(-0.5 * ∆AIC score for
that model). The Akaike weight for a model is this value divided by the sum
of these values across all models.Ӡ
Baysian Information Criterion, BIC .3
‫ למספרים‬AIC ‫ מ‬1‫סובלני פחות‬
††‫גדולים של פרמטרים‬
† http://www.brianomeara.info/tutorials/aic
†† Wagenmakers & Farrell 2004
BIC: -2*ln L + k*ln(n)
‫ גודל המדגם לא משחק‬AIC-‫ ב‬,!‫ מוכפל בגודל המדגם‬k :‫ שימו לב‬.1
‫)‪Generalized linear models (GLM‬‬
‫משמשים ב‪ GLM-‬כשה ‪ response variable‬אינו רציף (ספירות‪,‬‬
‫פרופורציות‪ ,‬בינארי וכו') ‪ -‬או כשהנחות המבחנים הפרמטרים‬
‫(התפלגות נורמלית‪ ,‬שיוויון שונויות) אינן מתקיימות‬
‫‪ GLM‬מורכב שלושה חלקים‪:‬‬
‫‪1. linear predictor; 2. link function; 3. error distribution‬‬
‫הראשון הוא ערך הפרמטר‪ ,‬השני מדבר על טרנספורמציה (למשל הוא ”‪“identity‬‬
‫כשאין טרנספורמציה ו ”‪ “log‬עבור טרנספורמציה לוגריתמית) והשלישי אומר מה‬
‫התפלגות ה ‪ – residuals‬למשל גאמא‪ ,‬פואסון או נורמלית‪ .‬במקרה הפרטי בו‬
‫‪ link=identity‬ו ‪ error=normal‬ה ‪ GLM‬זהה למודל לינארי "רגיל"‬
Generalized linear models (GLM)
:‫מורכב שלושה חלקים‬
1. linear predictor; 2. link function; 3. error distribution
“identity” ‫ השני מדבר על טרנספורמציה (למשל הוא‬,‫הראשון הוא ערך הפרמטר‬
‫“ עבור טרנספורמציה לוגריתמית) והשלישי אומר מה‬log” ‫כשאין טרנספורמציה ו‬
‫ במקרה הפרטי בו‬.‫ פואסון או נורמלית‬,‫ – למשל גאמא‬residuals ‫התפלגות ה‬
"‫ זהה למודל לינארי "רגיל‬GLM-‫ ה‬error=normal-‫ ו‬link=identity
model4<-glm(clutch~log10(SVL)+asin(elevation),family=Gamma)
Log link
arcsin link
Gamma errors
‫‪Non-linear models‬‬
‫לעיתים ברור שהיחס בין ה ‪predictor‬‬
‫ל ‪ response‬אינו לינארי‪.‬‬
‫))‪model2<-lm(y~x+I(x^2‬‬
‫ניתן לבחון מודלים שיודעים להתמודד עם מבנה כזה – למשל לעיתים‬
‫קרובות משלבים משוואה ריבועית (‪ )quadratic‬במודל‪:‬‬
‫‪Response = a(predictor)2+b(predictor)+c‬‬
‫וכרגיל אפשר לבדוק אם המודל הקוואדרטי טוב מהמודל‬
‫הלינארי על ידי ‪ AIC‬או על ידי ‪anova‬‬
Non-linear models
.‫ אינו לינארי‬response ‫ ל‬predictor ‫לעיתים ברור שהיחס בין ה‬
breakpoint ( ‫ניתן לבחון מודלים‬
‫) בהם יש משוואות‬regression
‫לינאריות שונות לערכים שונים של ה‬
predictor
Y = A1.x + K1
for x < breakpoint
Y = A2.x + K2
for x > breakpoint
Losos & Schluter 2000. Analysis of an evolutionary species-area relationship.
Nature 408: 847-850.
:‫ בכל עבודה‬,‫זיכרו‬
"No statistical procedure can substitute for
serious thinking about alternative
evolutionary scenarios and their credibility"
Westoby, Leishman & Lord 1995. On misinterpreting 'phylogenetic correction. J. of Ecology 83: 531-534.