Ezen a laboron gyakorlatias szemszögből fogjuk megközelíteni illetve átismételni a lineáris regressziót.
Ebben a fejezetben a Learning Statistics with R jegyzet alapján gyorsan átismételjük a lineáris regresszió lényegét. Ehhez először használjuk a parenthood
nevű adatokat, amit innen lehet letölteni a Data sets-re kattintva.
load("./parenthood.Rdata")
parenthood[1:5, ]
## dan.sleep baby.sleep dan.grump day
## 1 7.59 10.18 56 1
## 2 7.91 11.66 60 2
## 3 5.14 7.92 82 3
## 4 7.71 9.61 55 4
## 5 6.68 9.75 67 5
plot(parenthood$dan.sleep, parenthood$dan.grump,
xlab="Mennyit aludtam (órák)", ylab="Mennyire vagyok grumpy (0-100)",
pch=19)
Látjuk, hogy lineáris összefüggés van a két változó között, tehát \(Y_i = a_1 X_i + a_0 + \epsilon_i\), ahol \(Y\) azt jelöli, hogy mennyire vagyok rosszkedvű egy nap, \(X\) pedig, hogy hány órát aludtam aznap éjszaka Ha az \(Y\)-t előre szeretnénk jelezni \(X\) ismeretében, vagy le szeretnénk írni egy formulával \(X\) és \(Y\) összefüggését, akkor például végezhetünk egy lineáris regressziót, azaz illeszthetünk ezekre az adatokra egy egyenest, amit a következő formula ír le:
\[\hat{Y}_i = \hat{a}_1 X_i + \hat{a}_0\] Ahol \(X_i\) a magyarázó változó \(i\)-edik megfigyelése, amihez \(Y_i\) célváltozó tartozik, valamint \(\hat{Y}_i\) az \(Y_i\) becslése (a predikció amit a regresszió végez). Az \(\hat{a}_1\) jelöli az illesztett egyenes meredekségét (slope), \(\hat{a}_0\) pedig megmutatja, hogy hol metszi az egyenes az \(y\)-tengelyt (intercept). A különbséget egy valós megfigyelés és a becslése között reziduálisnak hívjuk (residual):
\[r_i = Y_i - \hat{Y}_i\]
azaz a teljes lineáris regresszió modell:
\[Y_i = \hat{a}_1 X_i + \hat{a}_0 + r_i\]
A legjobban illeszkedő egyenes az, ahol a reziduálisak “kicsik”. Pontosabban a becsült \(\hat{a}_0, \,\hat{a}_1\) együtthatók azok lesznek, amik minimalizálják a reziduálisok négyzetösszegét (sum of quared residuals): \(\sum_i (Y_i - \hat{Y}_i)^2\).
Megj.: Ezen a kurzuson a legkisebb négyzetek módszerével (ordinary least squares (OLS) regression) közelítjük az együtthatókat, végezzük el a regressziót; de nem csak ez az egy optimalizáló, eltérések/reziduálisok négyzetösszegét minimalizáló eljárás létezik.
Ismét.: R-ben az lm
függvénnyel (linear model) lehet lineáris regressziót végezni, aminek meg kell adnunk egy formulát és egy adathalmazt.
regression.1 <- lm(formula = dan.grump ~ dan.sleep,
data = parenthood)
Tehát a formulában a dan.grump
(rosszkedvűség) változót szeretnénk becsülni (kimeneti változó) a dan.sleep
magyarázó változó segítségével.
regression.1
##
## Call:
## lm(formula = dan.grump ~ dan.sleep, data = parenthood)
##
## Coefficients:
## (Intercept) dan.sleep
## 125.956 -8.937
Látjuk hogy az intercept \(\hat{a}_0 = 125.96\) és a meredekség \(\hat{a}_1 = -8.94\). Azaz az alábbi ábrán lévő egyenesnek az egyenlete a következő: \[ \hat{Y}_i = -8.94 X_i + 125.96\]
plot(parenthood$dan.sleep, parenthood$dan.grump,
xlab="Mennyit aludtam (órák)", ylab="Mennyire vagyok grumpy (0-100)",
pch=19)
abline(regression.1, col="red")
Modell értelmezése
A legfontosabb, hogy tudjuk hogyan kell értelmezni az együtthatókat. Kezdjük az \(\hat{a}_1\) együtthatóval, ami így 1 dimenzióban a meredekség: az \(\hat{a}_1 = -8.49\) azt jelenti, hogy ha növelem az \(X\) értékét 1-gyel, akkor az \(Y\) csökken 8.49-cel, azaz ha 1 órával többet alszom, akkor növelem a közérzetem, és a rosszkedvűségem 8.94 grumpiness ponttal kevesebb lesz. Mit jelent az \(\hat{a}_0\)? Definíció szerint \(\hat{a}_0\) az \(Y_i\) lineáris közelítése, ha \(X_i = 0\), azaz ha nem alszom semmit, akkor a rosszkedvűségem már olyan nagy lesz (\(Y_i = 125.96\)) hogy a skálán már nem is lehet mérni (mérhetetlenül nagy lesz).
A többváltozós lineáris regresszió csak annyiban különbözik az előbbitől, hogy több tagot (magyarázó változót) adunk az egyenlethez. Azaz úgy akarjuk előrejelezni a dan.grump
változót hogy egyszerre használjuk hozzá a dan.sleep
és a baby.sleep
változókat, ahol az utóbbi azt jelöli hogy Danielle Navarro fia mennyit aludt az éjszaka (a többi pedig Danielle rosszkedvűségét és alvását jelöli). Tehát legyen megint \(Y_i\) az \(i\)-edik napon a rosszkedvűség, \(X_{i 1}\) az hogy mennyit aludt Danielle az \(i\)-edik nap éjszaka (dan.sleep
) és \(X_{i 2}\) pedig jelölje azt hogy Danielle kisfia mennyit aludt az \(i\)-edik éjszakán. Tehát az elméleti regresszió modell: \(Y_i = a_2 X_{i 2} + a_1 X_{i 1} + a_0 + \epsilon_i\), az adatokra illesztett lineáris regresszió pedig:
\[Y_i = \hat{a}_2 X_{i 2} + \hat{a}_1 X_{i 1} + \hat{a}_0 + r_i \]
Habár most több együtthatónk van, az \(\hat{a}_0,\, \hat{a}_1\) és \(\hat{a}_2\) becslések számolása ugyanúgy történik: azok az együtthatók amik minimalizálják a reziduálisok négyzetösszegét.
Az R kódban csak a formulát kell kiegészítenünk az új magyarázóváltozóval:
regression.2 <- lm(dan.grump ~ dan.sleep + baby.sleep, data = parenthood)
regression.2
##
## Call:
## lm(formula = dan.grump ~ dan.sleep + baby.sleep, data = parenthood)
##
## Coefficients:
## (Intercept) dan.sleep baby.sleep
## 125.96557 -8.95025 0.01052
Látjuk, hogy a dan.sleep
változónak még mindig nagy az együtthatója, azaz minél több órát alszik Danielle, annál kevésbé lesz grumpy, viszont a baby.sleep
változó elég kicsit, ami az sugallja hogy nem befolyásolja annyira Danielle rosszkedvűségét az ha többet vagy kevesebbet aludt a fia (de ettől még lehet hogy szignifikánsan nem nulla ez az együttható).
plot(parenthood$baby.sleep, parenthood$dan.grump,
xlab="Mennyit aludtam (órák)", ylab="Mennyire vagyok grumpy (0-100)",
pch=19)
Két magyarázóváltozó esetén nem egy egyenest kapunk, hanem egy síkot. Magasabb dimenzióban pedig értelemszerűen hipersíkot.
Egyenest illeszteni bármire lehet de lehet, hogy amit kapunk az haszontalan. A regresszió csak egy \(\hat{Y}_i\) becslést ad \(Y_i\)-re, nyilván ha ez a kettő minél közelebb van egymáshoz, annál jobb a modellünk.
Az \(R^2\) érték
Az \(R^2\) definiálásához szükségünk van a reziduálisok négyzetösszegére (sum of the squared residuals: \(SS_{\text{res}}\)), amiről reméljük hogy minél kisebb.
\[SS_{\text{res}} = \sum_i (Y_i - \hat{Y}_i)^2\]
Pontosabban azt szeretnénk hogy a célváltozó varianciájához (total sum of squares: \(SS_{\text{tot}}\)) képest legyen kicsi,
\[SS_{\text{tot}} = \sum_i (Y_i - \overline{Y})^2\]
ahol \(\overline{Y}\) az \(Y_i\)-k átlaga. Az alábbi ábrán a piros négyzetek jelölik az átlagtól való eltérések négyzetét, a kék négyzetek pedig a reziduálisok négyzetét. A piros négyzetek területének összege az \(SS_{\text{tot}}\) a kék négyzetek területének összege pedig az \(SS_{\text{res}}\)
Az \(R^2\) értéket a következőképpen definiáljuk:
\[R^2 = 1 - \frac{SS_{\text{res}}}{SS_{\text{tot}}} \]
Tehát minél jobban illeszkedik a lineáris regresszió (\(f\) a jobb oldali képen) az átlaghoz képest, annál közelebb lesz az \(R^2\) értéke az 1-hez. Az \(R^2\) tehát azt mondja meg, hogy a magyarázóváltozókkal a célváltózó varianciájának mekkora része magyarázható meg.
Forrás: wikipédia
Számoljuk ki kézzel a regression.1
R négyzetét:
X <- parenthood$dan.sleep # the predictor / magyarázó változó
Y <- parenthood$dan.grump # the outcome / célváltozó
Az illesztett egyenes egyenlete a következő volt: \(\hat{Y} = -8.937 X + 125.956\), hozzuk létre az \(\hat{Y}\) változót is:
Y.pred <- -8.937 * X + 125.956
Most számoljuk ki az \(SS_{\text{res}}\) és \(SS_{\text{tot}}\) összegeket:
SS.res <- sum( (Y - Y.pred)^2 )
SS.tot <- sum( (Y - mean(Y))^2 )
print(SS.res)
## [1] 1838.714
print(SS.tot)
## [1] 9998.59
Ez a két nagy szám így nem sokat mond még nekünk, de már ki tudjuk számolni az \(R^2\) értéket:
R.squared <- 1 - (SS.res / SS.tot)
R.squared
## [1] 0.8161027
Látjuk, hogy az \(R^2 = .816\), ami azt jelenti, hogy a magyarázóváltozónk (dan.sleep
) megmagyarázza a célváltozó (dan.grump
) varianciájának 81.6%-át. Persze nem muszáj ezeket a kalkulációkat mindig végig csinálni, elég csak a summary()
függvényt használni.
summary(regression.1)
##
## Call:
## lm(formula = dan.grump ~ dan.sleep, data = parenthood)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.025 -2.213 -0.399 2.681 11.750
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 125.9563 3.0161 41.76 <2e-16 ***
## dan.sleep -8.9368 0.4285 -20.85 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.332 on 98 degrees of freedom
## Multiple R-squared: 0.8161, Adjusted R-squared: 0.8142
## F-statistic: 434.9 on 1 and 98 DF, p-value: < 2.2e-16
summary(regression.1)$r.squared
## [1] 0.8161027
Regresszió és korreláció kapcsolata
Tudjuk hogy a Pearson-féle korreláció lineáris kapcsolatot mér. Továbbá azt is tudjuk, hogy ez \(r\) betűvel szoktuk jelölni. Véletlen lenne? A korreláció négyzete \(r^2\) megegyezik az \(R^2\)-tel egyváltozós lineáris regresszió esetén:
r <- cor(X, Y) # X és Y korrelációja
r^2
## [1] 0.8161027
Látjuk, hogy a két szám megegyezik. Tehát más szavakkal, a Pearson-féle korreláció többé-kevésbé ekvivalens egy egydimenziós lineáris regresszióval.
Adjusted \(R^2\) (\(R_*^2\))
Geometriai megfontolásból könnyű meggondolni, hogy ha minél több magyarázóváltozót adunk a modellhez, annál nagyobb lesz az \(R^2\) (vagy legalábbis nem csökken). Az adjusted R négyzet (\(R_*^2\)) az R négyzetnek egy olyan módosítása, ami ezt a problémát kiküszöböli és figyelembe veszi hogy mennyi megfigyelésünk, illetve mennyi magyarázóváltozónk van. Például, ha csak 1 változónk van és három megfigyelésünk akkor az \((X, Y)\) pontok általános esetben nem esnek egy egyenesre (tehát \(R^2 <1\)), de ha hozzáadnánk még egy változót és marad ez a 3 megfigyelés, akkor ez a három megfigyelés már meghatároz egy síkot, azaz \(R^2 = 1\) lenne.
A lineáris regresszió, amit a tárgy keretein belül tárgyalunk rengeteg feltevésen alapszik:
Több módszert is lehet alkalmazni:
shapiro.test()
).hist(x = residuals( regression.2 ), xlab="Reziduális értéke", main="Reziduális hisztogram", breaks = 14)
shapiro.test(residuals(regression.2))
##
## Shapiro-Wilk normality test
##
## data: residuals(regression.2)
## W = 0.99228, p-value = 0.8414
Szerencsére meggyőződhetünk arról, hogy itt nincs veszélyben a normalitás feltevés.
Fontos, hogy ellenőrizzük, hogy valóban lineáris kapcsolat áll fenn a magyarázóváltozók és a célváltozó között. Itt, amit tehetünk, hogy megvizsgáljuk a becsült \(\hat{Y}_i\) és a valós \(Y_i\) kapcsolatát, azaz ábrázoljuk ezeket a pontokat egy pontdiagramon és reméljük hogy egy egyenes körül összpontosulnak. Az \(\hat{Y}\) kiszámolásához használhatjuk a fitted.values()
függvényt.
yhat <- fitted.values( regression.2)
plot( x = yhat,
y = parenthood$dan.grump,
xlab = "Fitted values",
ylab = "Observed values")
Egy másik, talán informatívabb ábra, ha nem a becsült pontok és a megfigyeléseket ábrázoljuk, hanem a reziduálisok és becsült értékek kapcsolatát (residual analysis). Ehhez az ábrához nem is kell külön kiszámolnunk a reziduálisokat, elég csak simán a plot()
függvény which
paraméterét átállítani.
Megj.: A which=2
beállítással lehet Q–Q plotot készíteni.
plot(x=regression.2, which=1)
Az ábrán nem csak a reziduálisokat látjuk az illesztett egyenes függvényében, hanem egy piros vonalat is ami mutatja ennek a két változónak a kapcsolatát. Ideális esetben egy egyenes horizontális vonalat kellene látnunk. Látjuk, hogy van egy kis görbület. Ennyiből még nem világos hogy kellene emiatt aggódnunk vagy sem, de könnyű meggondolni hogy ha valójában négyzetes összefüggés van \(X\) és \(Y\) között, akkor egy U alakú ábrát kellene kapnunk:
x <- seq(-1, 4, 0.1)
y <- x^2 + rnorm(length(x)) # Legyen y = x^2 + random zaj
plot(x, y)
abline(lm(y ~ x), col="blue")
A kék vonal a lineáris regresszió eredménye, nézzük ilyenkor hogy néz ki az előbbi ábra:
plot( lm(y ~ x), which=1 )
Ezzel az ábrával a homoszkedaszticitást is ellenőrizhetjük heurisztikusan (de erre is vannak szofisztikáltabb módszerek, pl. non-constant variance test: ncv.test
).
Gyakorlatban, ha azt sejtjük, hogy négyzetes összefüggés (is) állhat fenn, akkor hozzáadjuk a változó négyzetét (is) a lineáris modellhez:
x2 <- x^2
regression.3 <- lm(y ~ x + x2)
summary(regression.3)
##
## Call:
## lm(formula = y ~ x + x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.85308 -0.58520 -0.04232 0.47481 2.72563
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.008026 0.194350 -0.041 0.967
## x -0.148052 0.230131 -0.643 0.523
## x2 1.043730 0.070250 14.857 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9717 on 48 degrees of freedom
## Multiple R-squared: 0.9634, Adjusted R-squared: 0.9618
## F-statistic: 631.1 on 2 and 48 DF, p-value: < 2.2e-16
plot(regression.3 , which=1 )
Látjuk, hogy az \(x\) együtthatója kicsi, de ami fontosabb hogy az ehhez tartozó \(p\)-value nagy, azaz nem különbözik szignifikánsan a nullától az \(x\) együtthatója. Viszont az \(x^2\) nagyonis szignifikáns. Ezt az új mesterséges változót hozzáadva a modellhez már az ábra is jobban hasonlít az ideálishoz.
A kollinearitás vizsgálatra is léteznek kifinomult módszerek, pl. variance inflation factor (VIF), de ez túlmutat a tárgy keretein. Gyakorlatban sokszor csak a változók korrelációját ellenőrizzük például egy ú.n. hőtérkép/heatmap-en.
A példa kedvéért tegyük fel, hogy a parenthood
adatok day
változóját szeretnénk előrejelezni. Nézzük a korrelációkat:
cor( parenthood )
## dan.sleep baby.sleep dan.grump day
## dan.sleep 1.00000000 0.62794934 -0.90338404 -0.09840768
## baby.sleep 0.62794934 1.00000000 -0.56596373 -0.01043394
## dan.grump -0.90338404 -0.56596373 1.00000000 0.07647926
## day -0.09840768 -0.01043394 0.07647926 1.00000000
Látjuk hogy elég erős korrelációk vannak a magyarázóváltozók között! Nézzük hogy néz ki ugyanez egy heatmap-en:
cormat <- round(cor(parenthood), 2) # korrelációk kerekítve
library(reshape2)
melted_cormat <- melt(cormat) # erre azért van szükség mert a ggplot majd így tudja szépen ábrázolni
head(melted_cormat)
## Var1 Var2 value
## 1 dan.sleep dan.sleep 1.00
## 2 baby.sleep dan.sleep 0.63
## 3 dan.grump dan.sleep -0.90
## 4 day dan.sleep -0.10
## 5 dan.sleep baby.sleep 0.63
## 6 baby.sleep baby.sleep 1.00
library(ggplot2)
ggplot(data = melted_cormat, aes(x=Var1, y=Var2, fill=value)) +
geom_tile() +
scale_fill_gradient2(low = "blue", high = "orangered",
midpoint = 0, limit = c(-1,1),
name="Pearson\nCorrelation")
Miért okozhat gondot a modell értelmezésénél a multikollinearitás?
Tegyük fel, hogy előre szeretnénk jelezni a műszaki szakokra jelentkező diákok első féléves tanulmányi átlagát, és a magyarázó változók a érettségi tárgyakon elért pontszámok. Általában aki jó fizikából, az jó szokott lenni matematikából is és fordítva, tehát a fizika és matematika érettségi jegyek erősen korrelálnak. Arra számítunk, hogy mikor elvégezzük a lineáris regressziót, akkor mindkét tárgynak az együtthatója szignifikánsan nem nulla lesz, de egy ilyen vegyes modellben azt tapasztalhatjuk, hogy csak az egyikük lesz szignifikáns, viszont egyváltozós regresszióban mindkettő szignifikáns változó. Miért? Nem azért, mert az egyik egyáltalán nem fontos, hanem azért mert mivel erősen korrelálnak, ha az egyiket ismerjük, akkor az kölcsönösen szinte minden információt tartalmaz a másik változóról, azaz egyik ismeretében a másik már nem ad több információt.
Az eddig tárgyalt többváltozós lineáris modellben a marginális hatások a többi változó nagyságától függetlenül állandóak voltak, azaz ugyanúgy lehet értelmezni a többváltozós modellt mint az egyváltozósat: ha egy egységgel növelem az \(X_2\) változót akkor \(\hat{a}_2\) egységgel változik a célváltozó, ha minden mást változatlanul hagyok.
Ha két változó között interakció van, akkor az egyik marginális hatásának nagyságát befolyásolja a másik változó nagysága. Például tegyük fel, hogy meg akarjuk vizsgálni egy koleszterinszint csökkentő szer hatását, és arra is kíváncsiak vagyunk, hogy másképp reagálnak-e a férfiak a tablettára mint a nők. Tehát legyen az \(Y\) változó az, hogy mennyivel csökkent a koleszterinszint, \(X\) az hogy mekkora dózisban kapott a páciens koleszterint csökkentőt, és legyen \(G\) bináris változó, ami 0 értéket vesz fel, ha a páciens férfi, és 1 az értéke ha a páciens nő. Ekkor az interakciós taggal a modell a következő: \[ Y_i = a_1 X_i + a_2 X_i \cdot G_i + \epsilon_i \] Ekkor az együtthatókat a következőképpen lehet értelmezni: a férfiak esetén a második tag nulla (\(G_i = 0\)), ezért náluk ha a dózist 1 egységgel növeljük, akkor \(a_1\) egységgel változik a koleszterinszint. A nők esetében a második tag már nem nulla (\(G_i = 1\)), ekkor ha a dózist egy egységgel növeljük, akkor a koleszterinszint már \(a_1 + a_2\) egységgel változik. Például, ha az \(a_2\) negatív, akkor azt lehet mondani, hogy nőknél kevésbé csökkenti koleszterinszintet mint a férfiak esetében, ha pozitív, akkor pedig fordítva.
R példa
Bővítsük a múltórai modellünket interakciós tagokkal. R-ben a :
karakterrel lehet interakciós tagokat létrehozni, ami kezelni tudja a kategorikus és numerikus változókat is vegyesen.
Megj.: Interakciós tagok kezelésére van a *
operátor is, de annak kicsit bonyolultabb az “algebrája” és a működése.
library(readxl)
carprices<-read_excel("kuiper.xlsx")
Az alapmodell, amit kibővítünk legyen a következő:
regression.4 <- lm(Price ~ Mileage + Cylinder, carprices)
summary(regression.4)
##
## Call:
## lm(formula = Price ~ Mileage + Cylinder, data = carprices)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10264 -5121 -2838 3102 35477
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3145.75032 1325.93436 2.372 0.0179 *
## Mileage -0.15243 0.03464 -4.401 1.22e-05 ***
## Cylinder 4027.67463 204.61180 19.684 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8042 on 801 degrees of freedom
## Multiple R-squared: 0.3398, Adjusted R-squared: 0.3382
## F-statistic: 206.2 on 2 and 801 DF, p-value: < 2.2e-16
Látjuk, hogy mindkét változónk szignifikáns, és az együtthatókat a következőképp lehet értelmezni: ha egy egységgel növelem a Mileage
változót, akkor csökken a Price
0.15 egységgel, azaz minél több kilométert/mérföldet teszünk meg egy autóval, annál kisebb a piaci értéke. Viszont a Cylinder
együtthatója pozitív, tehát minél több hengeres egy autó motorja, annál drágább az autó.
Nézzük, mi történik akkor, ha nem numerikus a Cylinder
hanem kategorikus:
regression.5 <- lm(Price ~ Mileage + as.factor(Cylinder), carprices)
summary(regression.5)
##
## Call:
## lm(formula = Price ~ Mileage + as.factor(Cylinder), data = carprices)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16843 -4460 -2559 2555 28763
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.106e+04 7.153e+02 29.448 < 2e-16 ***
## Mileage -1.592e-01 3.075e-02 -5.178 2.84e-07 ***
## as.factor(Cylinder)6 2.132e+03 5.422e+02 3.932 9.14e-05 ***
## as.factor(Cylinder)8 2.102e+04 7.995e+02 26.293 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7138 on 800 degrees of freedom
## Multiple R-squared: 0.4804, Adjusted R-squared: 0.4785
## F-statistic: 246.6 on 3 and 800 DF, p-value: < 2.2e-16
Hasonló eredményt kapunk, mint az előbb, viszont ezzel a modell specifikációval nagyobb \(R^2\) értéket értünk el, úgy tűnik, hogy a hengerek számát több értelme van kategorikusan kezelni. Amikor numerikusan vettük be a modellbe, azzal úgy kezeltük a változót mintha a 4 és 6 hengeres motor között ugyanakkora lenne a különbség árban, mint a 6 és 8 hengeres motor között. Azzal, hogy faktorrá alakítjuk megengedjük kvázi azt, hogy ehhez a két különbséghez különböző értéket rendeljen. Az együtthatókból látszik is, hogy mekkora különbség van a modell becslésében: a 4 hengeres motorhoz képest, amit a kategorizálás az alapnak vett, a 8 henger ~21 ezerrel emelt az áron, míg a 6 henger csak ~2.1 ezerrel.
Nézzük meg a kategorikus Cylinder
és a Mileage
interakcióját:
regression.6 <- lm(Price ~ Mileage + as.factor(Cylinder) + Mileage:as.factor(Cylinder), carprices)
summary(regression.6)
##
## Call:
## lm(formula = Price ~ Mileage + as.factor(Cylinder) + Mileage:as.factor(Cylinder),
## data = carprices)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16788 -4488 -2431 2589 25190
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.944e+04 9.414e+02 20.648 < 2e-16 ***
## Mileage -7.839e-02 4.333e-02 -1.809 0.070778 .
## as.factor(Cylinder)6 3.812e+03 1.430e+03 2.666 0.007831 **
## as.factor(Cylinder)8 2.787e+04 1.955e+03 14.256 < 2e-16 ***
## Mileage:as.factor(Cylinder)6 -8.363e-02 6.694e-02 -1.249 0.211919
## Mileage:as.factor(Cylinder)8 -3.475e-01 9.070e-02 -3.831 0.000138 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7082 on 798 degrees of freedom
## Multiple R-squared: 0.4898, Adjusted R-squared: 0.4866
## F-statistic: 153.2 on 5 and 798 DF, p-value: < 2.2e-16
Látjuk, hogy egy kicsivel nagyobb \(R^2\)-et értünk el, nézzük az együtthatók értelmezését:
Zárás
Zárásként lássunk egy példát egy sokkal pontosabb modellre, ugyanezen a problémán, amiből pl. látszik, hogy a márka és a forma sokkal jobb becslést eredményez a használtautó árra, hiszen az \(R^2\) sokkal közelebb van az 1-hez, mint eddig.
regression.7 <- lm(Price ~ Make + Type + Mileage + as.factor(Cylinder) + Liter, carprices)
summary(regression.7)
##
## Call:
## lm(formula = Price ~ Make + Type + Mileage + as.factor(Cylinder) +
## Liter, data = carprices)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8666.7 -1242.1 -89.2 1100.1 14973.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.885e+04 8.924e+02 21.123 < 2e-16 ***
## MakeCadillac 1.450e+04 5.180e+02 28.002 < 2e-16 ***
## MakeChevrolet -2.271e+03 3.560e+02 -6.379 3.03e-10 ***
## MakePontiac -2.355e+03 3.639e+02 -6.473 1.69e-10 ***
## MakeSAAB 9.905e+03 4.502e+02 22.001 < 2e-16 ***
## MakeSaturn -2.090e+03 4.708e+02 -4.440 1.03e-05 ***
## TypeCoupe -1.164e+04 4.647e+02 -25.045 < 2e-16 ***
## TypeHatchback -1.173e+04 5.454e+02 -21.501 < 2e-16 ***
## TypeSedan -1.179e+04 4.111e+02 -28.670 < 2e-16 ***
## TypeWagon -8.157e+03 5.006e+02 -16.292 < 2e-16 ***
## Mileage -1.862e-01 1.064e-02 -17.492 < 2e-16 ***
## as.factor(Cylinder)6 -3.313e+03 6.200e+02 -5.343 1.20e-07 ***
## as.factor(Cylinder)8 -3.673e+03 1.246e+03 -2.947 0.0033 **
## Liter 5.697e+03 3.427e+02 16.624 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2463 on 790 degrees of freedom
## Multiple R-squared: 0.9389, Adjusted R-squared: 0.9379
## F-statistic: 934.1 on 13 and 790 DF, p-value: < 2.2e-16
A reziduálisok is rendben vannak “grafikusan”:
hist(x = residuals( regression.7 ), xlab="Reziduális értéke", main="Reziduális hisztogram", breaks = 100)
A varianciaanalízis (angol rövidítéssel ANOVA) több egymással rokon statisztikai módszer összefoglaló neve. Ezek a módszerek a különböző t-próbák (2-nél) több csoportra vett általánosításaiként is tekinthetőek. Az elméleti alapokat a varianciák felbonthatósága adja. A Stat1 kurzus kereteibe az egyszempontú varianciaanalízis fér bele.
Megj.: lásd a Bolla-Krámli könyv 266. oldalától
Az egyszempontos varianciaanalízis a független kétmintás t-próba általánosítása, amikor több mint két minta átlagát szeretnénk egymáshoz hasonlítani. Valamilyen szempont alapján \(k\) db külön csoportban végezzük a megfigyeléseket, így van \(k\) db többnyire nem egyező méretű minta: \(n=\sum_{i=1}^k n_i\). A minta összmérete továbbra is legyen \(n\). A modell feltevései:
A \(b_i\)-k helyett a \(b=a_i+m=a_i +\frac1n \sum_{i=1}^k n_i b_i\) felbontást használjuk, ahol az \(m\) a várható értékek súlyozott átlaga.
F-próba építhető a következő hipotézisre: \[ H_0:\ a_1=\cdots=a_k=0 \quad\text{vs.}\quad H_1: \exists i\text{, hogy}\ a_i\neq 0. \] Az ehhez tartozó próbastatisztika: \[ F=\frac{s_a^2}{s_e^2}=\frac{Q_a}{Q_e}\cdot \frac{n-k}{k-1}\sim \mathcal{F}_{(k-1,n-k)}, \] ahol:
amelyek a teljes mintára (\(n\)) vett négyzetes eltérések összegéből (\(Q\)) származnak: \(Q=Q_a+Q_e\). Részletekért lásd a tankönyv 270. oldalát és az itteni táblázatot.
Lineáris modellt feltételezünk: \(X_{ij}=m+a_i+\epsilon_{ij}\), és ez alapján becsüljük az \(a_i\)-ket.
Próbáljuk ki a fentieket a már vizsgált iris
adaton.
head(iris)
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
summary(iris)
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100 setosa :50
## 1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300 versicolor:50
## Median :5.800 Median :3.000 Median :4.350 Median :1.300 virginica :50
## Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199
## 3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800
## Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500
Vizsgáljuk meg, hogy a faj (Species
) szempontjából van-e szignifikáns eltérés a csészelevelek szélességében várható értékben.
Először nézzünk egy Box-Whiskers diagramot.
library(ggplot2)
ggplot(iris, aes(x=as.factor(Species), y=Sepal.Width))+
geom_boxplot()
Ez alapján mondható, hogy van különbség, de a 3*50-es elemszám lehet túl kevés ahhoz, hogy kép alapján döntsünk. A változók normalitását most tegyük fel, és az ábra alapján a szórások is kb. egyezőnek tűnnek.
Meg kellene győződni arról, hogy a fenti képletben \(m\)-mel jelölt paraméter (a csoportbeli várható értékek súlyozott átlag) nem nulla-e. Ez általában feltehető, hiszen többnyire valamilyen nemnegatív változóval dolgozunk, mint itt is. Egyébként előtesztelhető sima egymintás t-próbával:
t.test(iris$Sepal.Width)
##
## One Sample t-test
##
## data: iris$Sepal.Width
## t = 85.908, df = 149, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 2.987010 3.127656
## sample estimates:
## mean of x
## 3.057333
A csoportbeli szórások egyezőségét csak ránézésre fogadtuk el eddig, de ez is előtesztelhető az ú.n. Bartlett-teszttel, amelynek hipotézise: \[ H_0:\ \sigma_1=\cdots=\sigma_k \quad\text{vs.}\quad H_1: \exists i,j\text{, hogy}\ \sigma_i\neq \sigma_j. \] A tesztstatisztika részletei a tankönyv 273. oldalán. Ennek a próbának szintén feltétele a normalitás.
bartlett.test(Sepal.Width ~ Species, data=iris)
##
## Bartlett test of homogeneity of variances
##
## data: Sepal.Width by Species
## Bartlett's K-squared = 2.0911, df = 2, p-value = 0.3515
A \(p\)-érték alapján elfogadhatjuk a nullhipotézist, hogy a szórások megegyeznek a 3 kategóriánkra.
Ez a fajta előtesztelés a kétmintás esetben annak felelt meg, amikor a varianciára nézett F-próba után végeztünk független kétmintás t-próbát. A kurzus keretein túlmutat, milyen problémákat szülhet ez. Kétmintás esetben a Welch-féle t-próba volt az alternatíva, amit tudunk használni különböző varianciák esetén is. Az ANOVA esetén is van ilyen, erre később visszatérünk.
Most nézzük végre az egyszempontú ANOVA-nkat a beépített aov()
fv.-nyel. Ez a lineáris regresszióhoz használt lm()
fv.-hez hasonlóan formulát vár bemenetként, de máshogy is működik. A létrejövő objektum summary()
fv.-e adja vissza a fenn taglalt F-próba eredményeit.
anov<-aov(Sepal.Width ~ Species, data = iris)
print(anov)
## Call:
## aov(formula = Sepal.Width ~ Species, data = iris)
##
## Terms:
## Species Residuals
## Sum of Squares 11.34493 16.96200
## Deg. of Freedom 2 147
##
## Residual standard error: 0.3396877
## Estimated effects may be unbalanced
summary(anov)
## Df Sum Sq Mean Sq F value Pr(>F)
## Species 2 11.35 5.672 49.16 <2e-16 ***
## Residuals 147 16.96 0.115
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A \(p\)-érték alapján szignifikáns különbség van a 3 faj esetén a csészelevelek szélességében.
Az aov()
fv. valójában csak egy wrapper a következőre:
anova(lm(Sepal.Width ~ Species, data = iris))
## Analysis of Variance Table
##
## Response: Sepal.Width
## Df Sum Sq Mean Sq F value Pr(>F)
## Species 2 11.345 5.6725 49.16 < 2.2e-16 ***
## Residuals 147 16.962 0.1154
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Érdemes azt végiggondolni (vagy elolvasni a tankönyvben), hogy valójában, amit eddig csináltunk ANOVA címszóval az nem más, mint a lineáris modellben az együtthatók becsléseire együtt nézett F-próba, amit a summary()
fv. már a múltkor is kiadott nekünk (lásd utolsó sor):
summary(lm(Sepal.Width ~ Species, data = iris))
##
## Call:
## lm(formula = Sepal.Width ~ Species, data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.128 -0.228 0.026 0.226 0.972
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.42800 0.04804 71.359 < 2e-16 ***
## Speciesversicolor -0.65800 0.06794 -9.685 < 2e-16 ***
## Speciesvirginica -0.45400 0.06794 -6.683 4.54e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3397 on 147 degrees of freedom
## Multiple R-squared: 0.4008, Adjusted R-squared: 0.3926
## F-statistic: 49.16 on 2 and 147 DF, p-value: < 2.2e-16
Most vizsgáljuk azt, hogy a használtató árak logaritmusa (logPrice
) szignifikánsan eltér-e márkánként (Make
). Azért vizsgáljuk most az árak logaritmusát, mert az eredeti árak nagyon nem tűntek számunkra normálisnak, és a varianciában lévő különbségek is nagyon szembetűnőek voltak. Ilyenkor egy log-transzformáció tud segíteni valamennyit a helyzeten, de csodákra nyilván ez sem képes. Ezzel végülis azt feltételeztük most az ANOVA modelljében, hogy az igazi ár háttérváltozója log-normális eloszlású. Ez a valósághoz valamennyivel közelebb állhat: árak vagy kereseti adatok eloszlásai többnyire nem szimmetrikusak, ezek modellezésére valamilyen vastagabb farkú eloszlás többnyire értelmesebb feltételezés.
carprices$logPrice <- log(carprices$Price)
ggplot(carprices, aes(x=as.factor(Make), y=logPrice))+
geom_boxplot()
A Box-Whisker diagram alapján elég nagy különbség van az árakban. Ez nem annyira meglepő annak, aki tisztában van ezekkel a márkákkal. Meglepő lehet esetleg, hogy még a log-skálán is mennyi kiugró érték van több márkánál is.
Sajnos a varianciák sem tűnnek nagyon megegyezőnek, ezt a Bartlett-teszt is megerősíti:
bartlett.test(logPrice ~ Make, data=carprices)
##
## Bartlett test of homogeneity of variances
##
## data: logPrice by Make
## Bartlett's K-squared = 221.99, df = 5, p-value < 2.2e-16
Ez alapján a szórások szignifikánsan különböznek a csoportonként. Azt viszont meg kell jegyezni, hogy a Bartlett-teszt érzékeny a normalitásra, ami itt az ábra alapján is sérülni látszik. Egy nemparaméteres alternatíva, amelyik nem feltételez normalitást, ugyanerre a hipotézisre Fligner-Killeen teszt:
fligner.test(logPrice ~ Make, data=carprices)
##
## Fligner-Killeen test of homogeneity of variances
##
## data: logPrice by Make
## Fligner-Killeen:med chi-squared = 103.1, df = 5, p-value < 2.2e-16
Ez alapján sem egyeznek meg a varianciák.
Általában az egyszempontú ANOVA nem annyira érzékeny a normalitás sérülésére, de a varianciák különbözősége viszont problémás lehet. Emiatt a sima egyszempontú ANOVA a fenti F-próbával nem megfelelő modell itt.
Két népszerűbb módosítás is van:
oneway.test()
kruskal.test()
welch<-oneway.test(logPrice ~ Make, data = carprices)
welch
##
## One-way analysis of means (not assuming equal variances)
##
## data: logPrice and Make
## F = 500.59, num df = 5.00, denom df = 271.15, p-value < 2.2e-16
kruwal<-kruskal.test(logPrice ~ Make, data = carprices)
kruwal
##
## Kruskal-Wallis rank sum test
##
## data: logPrice by Make
## Kruskal-Wallis chi-squared = 476.5, df = 5, p-value < 2.2e-16
A két fenti próba feltételei jobban illenek erre az adatra. Ezek alapján magabiztosabban mondhatjuk, hogy szignifikáns különbség van a használt autó árak logaritmusában márkánként.