Ajánlott irodalom, segédanyag, DataCamp, játék

Lineáris regresszió

Ezen a laboron gyakorlatias szemszögből fogjuk megközelíteni illetve átismételni a lineáris regressziót.

Ismétlés a Learning Statistics with R jegyzet alapján

Egyváltozós lineáris regresszió

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).

Többváltozós lineáris regresszió (multiple linear regression)

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.

Regresszió illeszkedés jóságának mérése

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ó feltevései

A lineáris regresszió, amit a tárgy keretein belül tárgyalunk rengeteg feltevésen alapszik:

  • Normalitás. Mint a statisztikai modellek nagy részének, a lineáris regressziónak is feltétele valamilyen normalitás: itt ez az \(\epsilon_i\) hibák normális eloszlását jelenti, (tehát az \(X\) és \(Y\) lehet nem-normális amíg a reziduálisok normálisak).
  • Linearitás. Alapvető feltevés hogy az \(X\) és \(Y\) között valóban lineáris kapcsolat legyen. Mindegy hogy egy- vagy többváltozós regresszióról van szó.
  • Homoszkedaszticitás (homogén szórás). A regressziós modell felteszi, hogy minden \(\epsilon_i\) nulla várhatóértékű \(\sigma\) szórású háttéreloszlásból származik és a \(\sigma\) ugyanaz mindegyik reziduálisra (variancia koordinátánként azonos). Gyakorlatban nem tudjuk ellenőrizni hogy mindegyik azonos eloszlású, ehelyett azt vizsgáljuk hogy a szórása a reziduálisnak állandó.
  • Korrelálatlan magyarázóváltozók. Itt arról van szó, hogy nem akarjuk, hogy a magyarázóváltozók erősen korreláltak legyenek egymással (kollinearitás/multikollinearitás). A változók kollinearitása a modell értelmezésénél/kiértékelésénél okoz gondot. (Ennek a gyakorlati feltevésnek a design mátrixra vonatkozó teljes rang feltétel tekinthető valamelyest az elméleti megfelelőjének.)
  • Outlier mentesség. Ez is inkább egy gyakorlati feltevés, és az előbbiekből lényegében következik, de fontos megjegyezni, hogy ahogy a korreláció is, a lineáris regresszió is nagyon érzékeny a kiugró adatokra, outlierekre. Ha vannak vagy lehetnek outlierek az adatban, akkor a modell megbízhatósága megkérdőjelezhető. Lásd a példákat és ábrákat a Learning statistics with R 475. oldalától.

Eddig nem vizsgált gyakorlati dolgok

Reziduálisok normalitásának vizsgálata

Több módszert is lehet alkalmazni:

  • Első körben nem árthat egy hisztogramon ábrázolni a reziduálisokat.
  • Sharpio-Wilk tesztet (vagy más normalitástesztet) lehet futtatni (shapiro.test()).
  • A gyakorlatban gyakran használják a QQ-plotot.
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.

Linearitás vizsgálata

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.

Kollinearitás vizsgálata

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.

Interakciós tagok

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:

  • Továbbra is, minél nagyobb távolságot tettek meg az autóval, annál kevesebbet ér
  • Hasonlóan a 6 vagy 8 hengeres autók többe kerülnek mint a 4 hengeresek.
  • Az interakciós tagokat így értelmezhetjük itt:
    • Ha az autónak 6 vagy 8 hengere van, akkor több egységgel csökken az ára +1 mérföld megtétele után mint egy 4 hengeres autónak.
    • Az előzővel együtt tekintve, a több hengeres autók drágábbak, de megtett mérföldönként több egységgel csökken az áruk.
    • Megjegyezzük, hogy csak az egyik interakciós tag lett szignifikáns, itt közrejátszhat a multikollinearitás jelensége is.

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)

Varianciaanalízis - ANOVA

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.

Egyszempontú varianciaanalízis

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:

  • Az \(i\). csoportban \(X_{ij}\sim \mathcal{N}(b_i,\sigma^2),\ j=1\dots n_i\), és ugye \(i=1,\dots,k\).
  • A mintaelemek csoporton belül és csoportok között is függetlenek.
  • A szórások minden csoportban azonosak, de a várható értékek különbözhetnek.

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:

  • \(s_a\) a csoportok közötti empirikus szórás,
  • \(s_e\) a csoportokon belüli empirikus szórás,
  • \(Q_a\) a csoportok közötti négyzetes eltérések összege,
  • \(Q_b\) a csoportokon belüli négyzetes eltérések összege,

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.

R példa

Íriszek

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

Használtautó árak

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:

  • Welch-féle ANOVA (Welch F-test)
    • ez a normalitást még mindig felteszi, de a varianciák lehetnek különbözőek;
    • a kétmintás Welch-teszthez hasonló korrekcióval számol;
    • R-ben oneway.test()
  • Kruskal-Wallis rangösszeg teszt
    • ez a normalitást sem feltételezi;
    • Wilcoxon rangösszeg próbához hasonlóan számol;
    • R-ben 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.