Kategorički prediktori
autor: Blagoje Ivanović
# Prvo cemo se baviti binarnim prediktorima, tj. prediktorima koji uzimaju samo
# 2 razlicite vrednosti. Najcesce se u modelu predstavljaju kao indikatori,
# odnosno promenljive koje su 0 ili 1.
# Pogledajmo primer podataka koji pokazuju efekat upotrebe cigareta u trudnoci
# na tezinu novorodjenceta.
# Podatke preuzeti sa:
# http://www.matf.bg.ac.rs/p/files/94-birthsmokers.txt
birthsmokers <- read.table("http://www.matf.bg.ac.rs/p/files/94-birthsmokers.txt", header = TRUE)
# Promenljiva Smoke je ucitana kao factor, sto u R-u oznacava kategorcku
# promenljivu, koja uzima konacan broj vrednosti koje nisu uredjene.
birthsmokers$Smoke## [1] yes no yes no yes yes no yes no no yes no yes no no yes no
## [18] yes no yes no yes no yes yes no no yes no yes yes no
## Levels: no yes# Nacrtajmo podatke
plot(birthsmokers)
# Wgt (tezina) ima linearnu vezu sa Gest (nedelja trudnoce) i to cemo
# modelovati, ali cemo posmatrati i uticaj pusenja na tezinu.
# Pravimo model sa oba prediktora
model <- lm(Wgt ~ Gest + Smoke, birthsmokers)
summary(model)##
## Call:
## lm(formula = Wgt ~ Gest + Smoke, data = birthsmokers)
##
## Residuals:
## Min 1Q Median 3Q Max
## -223.693 -92.063 -9.365 79.663 197.507
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2389.573 349.206 -6.843 1.63e-07 ***
## Gest 143.100 9.128 15.677 1.07e-15 ***
## Smokeyes -244.544 41.982 -5.825 2.58e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 115.5 on 29 degrees of freedom
## Multiple R-squared: 0.8964, Adjusted R-squared: 0.8892
## F-statistic: 125.4 on 2 and 29 DF, p-value: 5.289e-15# Bitno je tumacenje koeficijenata modela. Koeficijent koji je vezan za Smoke je
# oznacen kao "Smokeyes", jer je faktorska promenljiva i R automatski kreira
# indikatorsku promenljivu za nivo "yes" koja je 1 ukoliko je Smoke 'yes' i 0
# ako nije. Za nivo "no" ne pravi posebnu promenljivu jer je jednoznacno
# odredjena sa indikatorom za "yes".
# Kako je kodirana promenljiva Smoke u modelu mozemo da vidimo iz model matrice
model.matrix(model)## (Intercept) Gest Smokeyes
## 1 1 38 1
## 2 1 38 0
## 3 1 36 1
## 4 1 34 0
## 5 1 39 1
## 6 1 35 1
## 7 1 40 0
## 8 1 42 1
## 9 1 37 0
## 10 1 40 0
## 11 1 36 1
## 12 1 39 0
## 13 1 39 1
## 14 1 39 0
## 15 1 35 0
## 16 1 39 1
## 17 1 41 0
## 18 1 42 1
## 19 1 38 0
## 20 1 39 1
## 21 1 42 0
## 22 1 38 1
## 23 1 37 0
## 24 1 42 1
## 25 1 41 1
## 26 1 39 0
## 27 1 40 0
## 28 1 42 1
## 29 1 35 0
## 30 1 41 1
## 31 1 38 1
## 32 1 36 0
## attr(,"assign")
## [1] 0 1 2
## attr(,"contrasts")
## attr(,"contrasts")$Smoke
## [1] "contr.treatment"# ili funkcijom contrasts
contrasts(birthsmokers$Smoke)## yes
## no 0
## yes 1# Interpretacija ovog koeficijenta je: Ako je imamo dva podatka za koja je
# vrednost Gest ista, dok je kod prvog Smoke "yes" a drugog "no", za prvi ce
# vrednost zavisne promenljive Wgt biti za 244.544 niza nego za drugi.
# U sustini, ovime pravimo dve regresione prave, jednu za Smoke "yes" i jednu za
# "no", koje dele koeficijent pravca (uz Gest), tj. paralelne su.
# Ilustrujmo na slici:
# Crtamo Wgt u odnosu na Gest i bojimo tacke u zavisnosti od vrednosti Smoke
plot(Wgt ~ Gest, birthsmokers,
col = ifelse(birthsmokers$Smoke == "yes", "red", "blue"))
# Prava koja bi bila bez razmatranja Smoke je
#abline(lm(Wgt ~ Gest, birthsmokers))
# Prava za Smoke == "no":
abline(model$coef[1], model$coef[2], col = "blue")
# Prava za Smoke == "yes":
abline(model$coef[1] + model$coef[3], model$coef[2], col = "red")
# Prave su paralelne, a Smoke=="yes" prava je ispod Smoke=="no" prave. Ovo znaci
# da je u proseku tezina bebe od majke pusaca 244.5 grama manja od tezine bebe
# majke nepusaca.
# Mozemo testirati da li je razlika izmedju pusaca i nepusaca znacajna
# testiranjem hipoteze H0: b_{smokeyes} = 0, sto vidimo da je znacajno, jer ima
# 3 zvezdice koeficijent uz smokeyes. Mogli smo i da proverimo interval
# poverenja za koeficijente
confint(model)## 2.5 % 97.5 %
## (Intercept) -3103.7795 -1675.3663
## Gest 124.4312 161.7694
## Smokeyes -330.4064 -158.6817# Interval poverenja za Smokeyes ne sadrzi nulu, pa odbacujemo H0.
#
# Mogli smo da promenimo referentni nivo faktora da ne bude Smoke=="no", tj. da
# kodiramo Smoke=="yes" sa 0, a Smoke=="no" sa 1.
birthsmokers$Smoke2 <- relevel(birthsmokers$Smoke, "yes")
# Kodiranje je sad:
contrasts(birthsmokers$Smoke2)## no
## yes 0
## no 1# Kad napravimo model sa referentnim nivoom "yes", dobijamo
model2 <- lm(Wgt ~ Gest + Smoke2, birthsmokers)
summary(model2)##
## Call:
## lm(formula = Wgt ~ Gest + Smoke2, data = birthsmokers)
##
## Residuals:
## Min 1Q Median 3Q Max
## -223.693 -92.063 -9.365 79.663 197.507
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2634.117 358.872 -7.340 4.37e-08 ***
## Gest 143.100 9.128 15.677 1.07e-15 ***
## Smoke2no 244.544 41.982 5.825 2.58e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 115.5 on 29 degrees of freedom
## Multiple R-squared: 0.8964, Adjusted R-squared: 0.8892
## F-statistic: 125.4 on 2 and 29 DF, p-value: 5.289e-15# Sada pise Smoke2no i koeficijent je 244.544, sto je suprotno od Smokeyes u
# prvobitnom modelu, sto je ocekivano.
# Zasto koristiti kategoricku promenljivu umesto dva razlicita modela, kad vec
# dobijamo dve prave?
# Podelimo podatke na pusace i nepusace i napravimo dva modela Wgt ~ Gest
model_smoke <- lm(Wgt ~ Gest, birthsmokers[birthsmokers$Smoke == "yes", ])
model_nosmoke <- lm(Wgt ~ Gest, birthsmokers[birthsmokers$Smoke == "no", ])
summary(model_smoke)##
## Call:
## lm(formula = Wgt ~ Gest, data = birthsmokers[birthsmokers$Smoke ==
## "yes", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -228.53 -64.86 -19.10 93.89 184.53
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2474.56 553.97 -4.467 0.000532 ***
## Gest 139.03 14.11 9.851 1.12e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 126.6 on 14 degrees of freedom
## Multiple R-squared: 0.8739, Adjusted R-squared: 0.8649
## F-statistic: 97.04 on 1 and 14 DF, p-value: 1.125e-07summary(model_nosmoke)##
## Call:
## lm(formula = Wgt ~ Gest, data = birthsmokers[birthsmokers$Smoke ==
## "no", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -171.52 -101.59 23.28 83.63 139.48
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2546.14 457.29 -5.568 6.93e-05 ***
## Gest 147.21 11.97 12.294 6.85e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 106.9 on 14 degrees of freedom
## Multiple R-squared: 0.9152, Adjusted R-squared: 0.9092
## F-statistic: 151.1 on 1 and 14 DF, p-value: 6.852e-09# Nacrtamo li ih...
abline(model_smoke)
abline(model_nosmoke)
# Dobijamo jako slicne prave kao i sa kategorickim, s tim sto su koeficijenti
# pravca za nijansu drugaciji.