Lineaarse regressiooni mudelid 1

Tänases praktikumis

Lineaarne regressioon
- arvulise seletava tunnusega
- kategoriaalse seletava tunnusega
ANOVA
Post-hoc test (Tukey HSD)

Kasutame andmestikke:

ldt.csv
teke2 (loeme otse dataDOI repositooriumist, ei pea faili alla laadima)
kandidaadid2019.RData

Seletav tunnus arvuline

Lineaarne regressioon seletab ühe arvulise tunnuse muutumist sõltuvalt teis(t)e tunnus(t)e muutumisest.
Kahe arvulise tunnuse puhul on see sama mis korrelatsioon, aga kui korrelatsiooni puhul me ei ütle, kumb tunnus teise muutuse põhjustab, siis lineaarses mudelis on sõltuv tunnus y ja seletav tunnus x.
Sõltuv tunnus peab olema arvuline tunnus.
Seletav tunnus võib olla arvuline või ka nominaalne kategoriaalne tunnus (faktor)

Mudeli eeldused

Nagu t-testigi puhul, on lineaarse regressiooni eelduseks mõõtmiste sõltumatus ja normaaljaotus, aga nüüd peaks testima pigem mudeli jääkide jaotust, mitte iga üksiku tunnuse oma. Ehk siis:

mõõtmised on üksteisest sõltumatud
mudeli jäägid on normaaljaotusega
mudeli jääkide dispersioon on homoskedastiline
mudeli jääkides ei ole autokorrelatsiooni
ja lisaks, et seos on lineaarne

Näide 1: Kas lühemaid sõnu tuntakse kiiremini ära?

Kasutame üleeelmisest korrast tuttavat Levshina õpiku andmestikku ldt. Eelmisel korral leidsime, et sõnade äratundmise keskmise reaktsiooniaja ja sõna pikkuse vahel oli positiivne korrelatsioon. Proovime sama lineaarse regressioonimudeliga testida.

ldt <- read.csv("data/ldt.csv")
plot(Mean_RT~Length, data=ldt, xlab="Sõna pikkus (tähemärkides)", ylab="Keskmine reaktsiooniaeg (ms)")
m <- lm(Mean_RT~Length, data=ldt, subset=Mean_RT < 1200)
abline(m, col="red")

Kas mudeli eeldused on täidetud?

Iga rida andmestikus käib ühe katses esinenud sõna kohta, mida testis esineb ühe korra. Katses osales rohkem kui üks katseisik, kes igale sõnale vastasid, aga andmestikus on iga sõna kohta rühma keskmised reaktsiooniajad.
Kas tunnused on normaaljaotusega? Testime!
- sõna pikkus on
- reaktsiooniaeg on, kui logaritmida
- reaktsiooniaeg on normaaljaotusega ka siis, kui välja jätta 3 eriti suurt väärtust :)

Sõltuv ja seletav tunnus

Kui korrelatsiooni puhul me testisime seda, kas kahe arvulise tunnuse vahel on seos, siis nüüd otsime põhjuslikku seost. Kahe arvulise tunnuse puhul on see tegelikult sõnastuslik küsimus, aga siiski tasuks läbi mõelda, kumb tunnus võiks põhjustada teise muutumist:

Kas see, kui kiiresti sõna ära tuntakse võib mõjutada seda, mitu tähemärki sõnas on?
Kas sõna pikkus võib mõjutada seda, kui kiiresti see ära tuntakse?

Lineaarse mudeli saab käsuga lm(), mida kasutasime juba regressioonijoone joonistamiseks.

ldt2 <- ldt[ldt$Mean_RT < 1200,]
m <- lm(Mean_RT~Length, data=ldt2)

Mudelist ülevaate saamiseks kasutame käsku summary()

summary(m)

Kas mudeli jäägid on normaaljaotusega?

Mudeli jäägid (residuals) on see osa uuritava tunnuse väärtustest, mis jääb mudelis kirjeldavate tunnuste kaudu seletamata. See on siis iga mõõtmispunkti kaugus regressioonijoonest.

plot(Mean_RT~Length, data=ldt2);abline(lm(m), col = "red")
for (i in 1:nrow(ldt2)) {
  lines(rep(ldt2$Length[i], 
            length(ldt2$Mean_RT[i]:m$fitted[i])),
        ldt2$Mean_RT[i]:m$fitted[i],
        lty = "dashed")
}

Mudeli väljund on oma olemuselt list (vt mode(m) ja str(m)). Mudeli jäägid on vektorina ..$residuals all ja selle vektori iga väärtus vastab ühele andmestiku mõõtmispunktile.

Lineaarse mudeli eeldus on, et mudeli jäägid on normaaljaotusega. Testime:

hist(residuals(m))

shapiro.test(residuals(m))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(m)
## W = 0.98574, p-value = 0.3801

Kuidas mudelit lugeda?

summary(m)

## 
## Call:
## lm(formula = Mean_RT ~ Length, data = ldt2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -259.247  -63.097   -0.795   53.320  224.378 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  547.379     35.950  15.226  < 2e-16 ***
## Length        30.242      4.261   7.096 2.29e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 99.96 on 95 degrees of freedom
## Multiple R-squared:  0.3465, Adjusted R-squared:  0.3396 
## F-statistic: 50.36 on 1 and 95 DF,  p-value: 2.289e-10

Kas mudel on oluline? Mudeli väljundi viimane rida annab mudeli üldise p-väärtuse, mis on 2.289e-10 ehk 2.289 * 10 ^ -10 ehk 0.0000000002289. Niisiis mudel on oluline tasemel p<0.001***. See tähendab, et tõenäosus juhuslikult selline tulemus saada on väga väike ja võime hüljata 0-hüpoteesi.
Kui hea mudel on? Eelviimane rida, R² (R-squared) näitab protsentuaalselt, kui palju uuritava tunnuse varieerumisest mudel kirjeldab. Antud juhul 0.3464513 ehk võib öelda, et sõna pikkus kirjeldab 35% reaktsiooniaja varieerumisest.
Estimate ehk mudeli prognoos näitab, kui palju kirjeldava tunnuse muutudes uuritav tunnus muutub.
- Esimene rida (Intercept) ehk vabaliige ütleb, mis on uuritava tunnuse väärtus siis, kui seletava tunnuse väärtus on 0.
- Teine rida (slope) ehk kalle ütleb, kui palju uuritav tunnus muutub siis, kui seletava tunnuse väärtus muutub ühe ühiku võrra.

Vaatame joonist uuesti natuke lähemalt:

Vabaliige on y telje väärtus, kui x on 0, kalle on y väärtuse muutus, kui x on 1. See, et tegelikult andmed x-teljel 0 ja 1 punkte ei läbi (ei saagi läbida, sest ühegi sõna pikkus ei saa olla 0), ei oma tähtsust. Kuna tegemist on lineaarse regressiooniga, siis me võime vabaliikme ja kalde väärtuste põhjal tuletada ükskõik millise punkti väärtused.

Kuna suhe on lineaarne, siis eeldame, et kirjeldava tunnuse kasvades näiteks 5 ühiku võrra kasvab sõltuv tunnus samuti 5*regressioonikorda võrra.
Mudeli võib kirja panna y = A + Bx
Reaktsiooniaeg = 547 + 30 * sõna pikkus

Harjutus

Harjutus e Näide 2: Kas kõnetempo sõltub vanusest?

Kasutame kõnetempo andmestikku artiklist Lippus, Pilvik, Lõo & Lindström 2024 (ilmub aprillis Rakenduslingvistika Ühingu aastaraamatus) ja vaatame, kas kõnetempo sõltub vanusest.

Loeme DataDOI repositooriumist TeKE korpusel põhineva teismeliste kõnelejate spontaanse kõne andmestiku.

load(url("https://datadoi.ee/bitstream/handle/33/592/teke_globaalne_konetempo.Rda?sequence=21&isAllowed=y"))

Arvutame kõnetempo, mis peaks olema kõneleja silpide arv jagatud tema kõnevoorude kestusega:

teke2$konetempo <- teke2$silpide_arv/teke2$kestus

Kontrollime mudeli eeldusi:

Sõltumatute tunnuste eeldus: iga vestluse iga kõneleja kohta on andmestikus küll üks rida, aga on suurem osa kõnelejaid osalesid rohkem kui ühes vestluses. Õigem oleks kasutada segamudelit.
Mis jaotusega on kõnetempo ja vanus? Kontrollime qqnorm ja Shapiro testiga.

Mis võiks olla uuritav tunnus, mis kirjeldav? Sõnastame hüpoteesid.
Teeme lineaarse mudeli.
Kas mudeli jäägid on normaaljaotusega?
Vaatame mudeli väljundit.

Kas võime 0-hüpoteesi ümber lükata? Kas mudel on oluline?
Kui palju lehekülgede arvu varieeruvusest väitekirja kaitsmise aasta seletab?
Kuidas tõlgendada mudeli prognoositavaid väärtusi?

Teeme joonise ka

Vaatame veel korraks mudeli koefitsiente.

## (Intercept)       vanus 
##   3.1210763   0.1631468

Mis on mudeli vabaliige (Intercept) väärtus? Mudel arvestab seda lõikumisena 0-punktis ehk et see väärtus oleks uuritaval tunnusel siis, kui seletava tunnuse vanus oleks 0. Tegelikult meil algavad vaatlused alles väärtusest 10. Kui eeldada lineaarset suhet, siis mudeli järgi oli vanuses 0 doktoritööde maht 3 silpi/sek. Eeldades lineaarset kasvu, hindaks see mudel, et kasvades iga eluaastaga 0.163 silpi/sek võrra, jõudis tööde maht 10. eluaastaks 5 silp/sek ja 18. eluaastaks 6 silp/sek.

Eeldades, et seos on lineaarne, saame me selle kehtimist laiendada näiteks inimese elueale 0-90 eluaastani

plot(konetempo~vanus, data=teke2, main="Kõnetempo muutumine sõltuvalt vanusest", xlim=c(0, 90), ylim=c(3,18))
abline(mudel)

Tegelikult ei peaks muidugi vaatama mudelit tegelike andmete muutumispiirkonnast väljaspool. Mõningatel juhtudel on kasulik teisendada andmeid nii, et 0 oleks näiteks esimene vaadeldav aasta, kust meil on mõõtmisi (st siin 10-aastased), sellisel juhul on vabaliige ja kirjeldavate tunnuste efektid paremini interpreteeritavad. Siin andmestikus väga suurt nihet sellega pole, aga kui me näiteks vanuse asemel hindaks sünniaastat ja need algaksid umbes 2000. aastast, siis lineaarne kasv 0.163 silpi/sek võrra aastas eeldaks, et 0 aastal oli kõnetempo -323 silpi sekundis, mis ei meigi mingit senssi.

Seletav tunnus nominaalne

Lineaarse regressiooni mudeli seletavaks tunnuseks võib olla ka nominaalne tunnus ehk faktor.
Faktori puhul valitakse üks tase baastasemeks.
Teiste tasemete väärtust hinnatakse baastaseme suhtes.
Baastase on mudeli vabaliige ehk intercept ning iga teise taseme kohta arvutatakse erinevus baastasemest.
Erinevus arvulise kirjeldava tunnusega mudelist on see, et kui arvulise tunnuse puhul on üks kalde kordaja (slope), siis nominaalsel tunnusel on iga baastasemest erineva taseme kohta eraldi kordaja. Arvulise tunnuse puhul saab kalde kordaja tunnuse väärtusega läbi korrutades välja arvutada kalde väärtuse igas skaalalpunktis. Nominaalse tunnuse puhul rakendatakse vastava taseme kalde kordajat binaarselt TRUE/FALSE väärtustega.

Kahe tasemega faktor

Näide 3: Kas kandidaadi sugu mõjutab häältesaaki (2019 riigikogu valmisite andmetes)?

Võtame valimistulemuste andmestiku ja vaatame alustuseks kahe tasemega faktorit: kas sugu mõjutab häältesaaki? Me oleme juba t-testiga teinud kindlaks, et meeskandidaadid said keskmiselt rohkem hääli kui naiskandidaadid, kuigi efekti suurus oli väga väike.

load("data/kandidaadid2019.RData")

Kas mudeli eeldused on täidetud?

Tunnus hääli_kokku ei ole normaaljaotusega, aga kui see logaritmida, siis vastab normaaljaotuse tunnustele. Kuna kandidaate oli rohkem kui 30, siis võiks ka normaaljaotuse eeldust ignoreerida, aga kui logaritmimine aitab normaalajotust saavutada, siis võiks seda ikkagi teha.
Mõõtmised on sõltumatud: iga kandidaat on individuaalne isik ja esindatud andmestikus ühe reaga.

Sõltuva ja sõltmatu tunnuse määramisel ei ole siin ilmselt vaja pikalt mõelda, sest juba tehniliselt kui me tahame kasutada lineaarset mudelit, siis uuritav tunnus peab olema arvuline. Ja kui oleks võimalus, et kandidaadi sugu võiks mõjutada see, kui palju ta hääli sai, siis peaksime valima hoopis logistilise regressioonimudeli, millest tuleb juttu umbes üle-ülejärgmises praktikumis :)

Milline tase valida baastasemeks? Mõnel juhul on lihtne ja selge, et üks väärtus on kuidagi moodi vaikimisi väärtus ja teised väärtused on markeeritud. Aga ka siis, kui meil pole mingit sisulist või ideoloogilist eelistust, peame valima ühe taseme baastasemeks.

Kahetasemelise faktori puhul ei ole mudeli väljundi tõlgenduse seisukohast väga suurt vahet, mis on baastase, aga sellest sõltub, mis on vabaliige ehk intercept. Vaikimisi on baastase faktori esimene tase ja kui me pole nominaalse tunnuse tasemete järjekorda faktorina eksplitsiitselt määranud, siis enamasti on see lihtsalt tähestikulises järjekorras.

sugu.lm <- lm(log(hääli_kokku)~sugu, data=kandidaadid2019)
summary(sugu.lm)

## 
## Call:
## lm(formula = log(hääli_kokku) ~ sugu, data = kandidaadid2019)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.0875 -1.1053  0.0254  1.0231  5.0573 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.08749    0.05773  88.120   <2e-16 ***
## sugunaine   -0.23767    0.10018  -2.372   0.0178 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.564 on 1097 degrees of freedom
## Multiple R-squared:  0.005104,   Adjusted R-squared:  0.004197 
## F-statistic: 5.628 on 1 and 1097 DF,  p-value: 0.01785

Muus osas on nominaalse faktoriga mudeli väljund sama nagu arvuilse sõltumatu tunnuse puhul:

mudeli üldist olulisust näitab p-väärtus viimasel real, mis ütleb, kui tõenäone on sellist tulemust saada juhuslikult. p = 0.01785 < ɑ, seega võib mudelit oluliseks pidada.
mudeli üldist headust nätab R-ruut eelviimasel real, mis ütleb, mitu protsenti uuritava tunnuse varieerumisest seletav(ad) tunnused kirjeldab. Seos on võrdlemisi nõrk, R-ruut on 0.005 ehk sugu kirjeldab ainult 0.5% häältesaagi varieerumisest, 99,5% ulatuses sõltub see millestki muust.

AGA nüüd on vabaliikmeks ehk mõtteliseks 0-punktiks faktori baastase. St esimene rida (Intercept) ütleb, mis on häältesaagi väärtus juhul, kui sugu on baastasemel ehk “mees” ning teine rida ehk slope ütleb seda, kui palju see muutub siis, kui faktori tase on “naine”.

Kuidas tulemust tõlgendada?

summary(sugu.lm)

## 
## Call:
## lm(formula = log(hääli_kokku) ~ sugu, data = kandidaadid2019)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.0875 -1.1053  0.0254  1.0231  5.0573 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.08749    0.05773  88.120   <2e-16 ***
## sugunaine   -0.23767    0.10018  -2.372   0.0178 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.564 on 1097 degrees of freedom
## Multiple R-squared:  0.005104,   Adjusted R-squared:  0.004197 
## F-statistic: 5.628 on 1 and 1097 DF,  p-value: 0.01785

Mudel annab prognoositud väärtused: meeste häälesaak on 5.08749 ühikut ja naiste oma on meeste omast -0.23767 võrra erinev. Ühikud on üldiselt samad mudeli sisendiks olnud ühikud, ehk hääled, AGA siin peab meeles pidama, et kui normaliseerimise eesmärgil sai uuritav tunnus logaritmitud, siis mudeli prognoositud väärtustest algse ühiku kätte saamiseks peame tegema logaritmimise pöördtehte ehk astendama:

# Meeste häältesaak on vabaliikme väärtus
exp(sugu.lm$coefficients[1])

## (Intercept) 
##    161.9826

# Naiste väärtuse saamiseks tuleb koefitsiendid kokku liita
exp(sugu.lm$coefficients[1]+sugu.lm$coefficients[2])

## (Intercept) 
##    127.7179

Teeme selguse mõttes sama mudeli ka ilma logaritmimata uuritava tunnusega, see küll eirab mudeli normaaljaotuse eeldust, aga on oluliselt lihtsam tõlgendada:

sugu.lm2 <- lm(hääli_kokku~sugu, data=kandidaadid2019)
summary(sugu.lm2)

## 
## Call:
## lm(formula = hääli_kokku ~ sugu, data = kandidaadid2019)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
##  -542.2  -442.7  -344.1   -91.2 19627.9 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   543.16      47.28  11.488   <2e-16 ***
## sugunaine     -98.08      82.04  -1.195    0.232    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1281 on 1097 degrees of freedom
## Multiple R-squared:  0.001301,   Adjusted R-squared:  0.0003907 
## F-statistic: 1.429 on 1 and 1097 DF,  p-value: 0.2322

Nüüd võime koefitsientide põhjal öelda, et mehed said keskmiselt 543 häält ja naised 98 häält vähem ehk keskmiselt 445 häält. Ja saame mudeli prognoositud väärtusi otse võrrelda aritmeetiliste keskmiste väärtustega:

tapply(kandidaadid2019$hääli_kokku, kandidaadid2019$sugu, mean)

##     mees    naine 
## 543.1567 445.0795

Mitme tasemega faktor

Näide 4: Kas haridus mõjutab häältesaaki?

Proovime nüüd testida seletavat tunnust, millel oleks rohkem kui kaks taset. Näiteks kandidaatide haridus:

table(kandidaadid2019$haridus)

## 
##                      Algharidus Keskharidus (sh keskeriharidus) 
##                               1                             265 
##                     Kõrgharidus                     Põhiharidus 
##                             819                              14

Haridusel on kandidaatide andmestikus 4 taset, aga kuna algharidusega kandidaate on ainult 1, siis ilmselt oleks mõistlik jätta see kandidaat analüüsist välja, sest üks mõõtmispunkt on liiga vähe, et seost kirjeldada.

haridus.lm <- lm(log(hääli_kokku)~haridus, data=kandidaadid2019, subset = haridus != "Algharidus")
summary(haridus.lm)

## 
## Call:
## lm(formula = log(hääli_kokku) ~ haridus, data = kandidaadid2019, 
##     subset = haridus != "Algharidus")
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.3599 -1.0568  0.0468  0.9530  4.6607 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         4.35993    0.09282  46.971   <2e-16 ***
## haridusKõrgharidus  0.88651    0.10679   8.302    3e-16 ***
## haridusPõhiharidus -1.02977    0.41437  -2.485   0.0131 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.511 on 1095 degrees of freedom
## Multiple R-squared:  0.07314,    Adjusted R-squared:  0.07144 
## F-statistic:  43.2 on 2 and 1095 DF,  p-value: < 2.2e-16

Mudeli väljund ütleb meile et:

mudel on erinev juhuslikust, oluline, p<0.001
seos on endiselt võrdlemisi nõrk, kuigi tugevam, kui soo puhul: mudel kirjeldab 7% häältesaagi varieerumisest.
Mudeli prognoosi osa, koefitsientide tabel toob välja kirjeldava tunnuse võrdlused baastaseme ja teiste tasemete vahel. Siin võiks vaadata nii esimest tulpa Estimate ehk mudeli hinnang väärtuse suurusele ja viimast tulpa p-väärtus ehk hinnang taseme hinnangu olulisusele.
- Esimesel real on intercept ehk baastase, siin mudelis keskharidus, mille väärtus on 4.36 logartimitud häält ehk 78.25 häält. Siin real p <0.001*** ei ole eriti tähenduslik (tähendab, et on nullist erinev);
- Teisel real on kõrgharidusega kandidaatide erinevus keskharidusega kandidaatide häältesaagist ehk 4.36. Kui see oleks algsetes ühikutes, siis võiksime öelda, et see rühm sai nii palju rohkem hääli, aga kuna me logaritmisime väärtused enne modelleerimist, siis peame kõigepealt prognoositud tulemuse kokku arvutama ja siis astendama (käsuga exp()): kõrgharitute logaritmitud häältesaak oli 4.36 + 0.89 = 5.25 ehk 189.89. Ja see erinevus keskharidusega kandidaatidest on oluliselt erinev p<0.001.
- Kolmandal real on põhiharidusega kandidaatide erinevus Keskharidusega kandidaatidest, mis on jällegi 4.36 + -1.03 = 3.33 ehk 27.94 häält. See erinevus on ka oluliselt erinev juhuslikust tasemel p<0.05.
- Mudel otseselt mittebaastasemete erinevust ei võrdle, sest kõigi tasemete hinnatavad väärtused on esitatud ja selle kaudu on teada, kas tasemete hulgas leidub mõni, mis on teistest erinev, mille põhjal otsustada, kas faktoril üldiselt on mõju või ei ole.

ANOVA

ANalysis Of VAriance ehk dispersioonanalüüs

ANOVA võrdleb rühmade hajuvust valimi üldise hajuvusega. R-is on ANOVA seotud lineaarse mudeliga, selle saab, kui lineaarne mudel anda käsu anova() sisendiks. Mõnes mõttes see on lihtsalt alternatiivne viis seost kirjeldada, selle kaudu saab anda ühe üldhinnangu kategoriaalse faktori olulisusele, samas kui lineaarne mudel annab iga taseme võrdluse baastasemega.

Võimaldab hinnata faktori üldist mõju uuritavale tunnusele
Kui faktoril on kaks taset, siis on tulemus sarnane t-testile
Kui tasemeid on rohkem, siis osutab faktori üldisele mõjule, aga mitte sellele, milliste faktori tasemete vahel erinevus on
- ANOVA ütleb seda, kas faktori tasemete hulgas leidub mõni tase, mis on teistest oluliselt erinev, aga ei ülte, mitu või milline,
- ja kui kolme- või enamatasemeline faktor on oluline, peame lisaks tegema post-hoc testi, kui tahame välja selgitada, milliste tasemete vahel tegelikult erinevused on ja milliste tasemete vahel ei ole.

ANOVA, faktor 2 tasemega

Näide 3b: Sugu ja häältesaak ANOVAga

Anova testi tegemiseks paneme lihtsalt lineaarse mudeli väljundi summary() käsu asemel anova() käsu sisse:

sugu.lm <- lm(log(hääli_kokku)~sugu, data=kandidaadid2019)
anova(sugu.lm)

## Analysis of Variance Table
## 
## Response: log(hääli_kokku)
##             Df  Sum Sq Mean Sq F value  Pr(>F)  
## sugu         1   13.77 13.7696  5.6281 0.01785 *
## Residuals 1097 2683.88  2.4466                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Kuidas väljundit lugeda?

Tabelis on rida iga faktori kohta (siin üks faktor: sugu) ja lõpuks rida mudeli jääkide kohta.
Iga faktori puhul hinnatakse rühmadevahelist dispersiooni valimi üldise populatsiooniga, millest kokku tuleb ANOVA statistik F.
F väärtuse statistiline olulisus sõltub vabadusastmetest. Faktori vabadusasmed on tasemete arv miinus üks, ehk soo puhul 2-1=1; Valimi vabadusastmete arv on mõõtmiste arv miinus võrreldavate konditsioonide arv ehk 1099 (kandidaati) - 2 (mehed/naised) = 1097.

Kahe tasemega faktori puhul on tulemus võrdlemisi sarnane t-testiga.

t.test(log(hääli_kokku)~sugu, data=kandidaadid2019)

## 
##  Welch Two Sample t-test
## 
## data:  log(hääli_kokku) by sugu
## t = 2.4084, df = 756.57, p-value = 0.01626
## alternative hypothesis: true difference in means between group mees and group naine is not equal to 0
## 95 percent confidence interval:
##  0.04394483 0.43138548
## sample estimates:
##  mean in group mees mean in group naine 
##            5.087489            4.849824

ANOVA, faktor 3 tasemega

Näide 4b: Haridus ja häältesaak

Nüüd näide rohkem kui kahe tasemega faktorist: kandidaatide andmestikus haridus on 4 taset, või kui ainult ühe mõõtmisega põhiharidus välja visata, siis on 3 taset. Ainult üks mõõtmine on endiselt liiga vähe. Aga kolme rühmaga enam t-testi teha ei saa.

haridus.lm <- lm(log(hääli_kokku)~haridus, data=kandidaadid2019, subset = haridus != "Algharidus")
anova(haridus.lm)

## Analysis of Variance Table
## 
## Response: log(hääli_kokku)
##             Df  Sum Sq Mean Sq F value    Pr(>F)    
## haridus      2  197.27  98.637  43.202 < 2.2e-16 ***
## Residuals 1095 2500.07   2.283                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova tulemuse põhjal võime öelda, et haridusel on oluline üldefekt.
See tähendab, et “leidub vähemalt üks rühm, mille keskmine on teistest oluliselt erinev”.
Artiklis saab ANOVA tulemuse teksti sees kirja panna nii:
F(2, 1095) = 43.202, p < 0.001
Kuna siin on rohkem kui 2 rühma, siis selleks, et otsustada, kes kellest erines, peaks tegema post-hoc testi.

Post-hoc test

Iga kord, kui me teeme mõne statistilise testi, on teatav võimalus, et me saame olulise tulemuse juhuslikult. Sama andmestiku peal veidi erinevate kombinatsioonidega sama testide tegemisel kõik juhuslikud efektid kasvavad. Samas, kui oleme ANOVA testiga teada saanud, et rohkem kui kahe tasemega faktoril on oluline üldefekt, tahaks me teada, milliste faktori tasemete vahel erinevus oli. Selleks peaks faktori tasemeid paari kaupa testima, aga tulemuste tõlgendamisel arvestama kordustega. Seda nimetatakse post-hoc testimiseks.

Bonferroni parandus

Üks võimalus oleks korrata lineaarset mudelit nii mitu korda, nagu on tarvis kõigite tasemete võrdlemiseks, muutes iga testi korral vastavalt faktori baastaset. Seejärel tulemuste tõlgendamisel jagatakse usaldusnivoo korduste arvuga. Seda nimetatakse Bonferroni paranduseks (Bonferroni correction). Näiteks kui meil on 3 faktori taset, siis kõigi võrdluste saamiseks peame kordama testi 2 korda ja langetama ɑ taset vastavalt 0.05/2 = 0.025.

Näites 4 saime teada, et keskharidusega (baastase) kandidaatide häältesaak oli oluliselt erinev kõrgharidusega (p<0.001) ja põhiharidusega (p<0.05) kandidaatidest:

summary(haridus.lm)

ning ANOVA testi põhjal võime öelda, et haridus on oluline faktor [F(2, 1095) = 43.202, p < 0.001], sest leidub vähemalt üks haridustase, mille puhul keskmine häältesaak on teistest erinev:

anova(haridus.lm)

siis me ei tea, kas kõrghariduse ja põhiharidusega kandidaatide häältesaak omavahel ka erinevad on.

Kordame lihtsalt lineaarset mudelit muudetud baasväärtusega. Sellks võib andmestikus muuta ära baasväärtuse käsuga relevel() (või lihtsalt käsuga factor()), aga kui me ei taha andmestikku muuta, võib ainult lm() käsu sees muuta faktori baastaset käsuga C(). Lisaks, kuna ühe haridustaseme jätsime välja, tuleks kasutada käsku droplevels(), mis eemaldab faktoritasemed, mida ei esine.

haridus.lm1 <- lm(log(hääli_kokku)~C(droplevels(haridus), base=1), data=kandidaadid2019[ kandidaadid2019$haridus != "Algharidus",])
haridus.lm2 <- lm(log(hääli_kokku)~C(droplevels(haridus), base=3), data=kandidaadid2019[ kandidaadid2019$haridus != "Algharidus",])
summary(haridus.lm1)

## 
## Call:
## lm(formula = log(hääli_kokku) ~ C(droplevels(haridus), base = 1), 
##     data = kandidaadid2019[kandidaadid2019$haridus != "Algharidus", 
##         ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.3599 -1.0568  0.0468  0.9530  4.6607 
## 
## Coefficients:
##                                             Estimate Std. Error t value
## (Intercept)                                  4.35993    0.09282  46.971
## C(droplevels(haridus), base = 1)Kõrgharidus  0.88651    0.10679   8.302
## C(droplevels(haridus), base = 1)Põhiharidus -1.02977    0.41437  -2.485
##                                             Pr(>|t|)    
## (Intercept)                                   <2e-16 ***
## C(droplevels(haridus), base = 1)Kõrgharidus    3e-16 ***
## C(droplevels(haridus), base = 1)Põhiharidus   0.0131 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.511 on 1095 degrees of freedom
## Multiple R-squared:  0.07314,    Adjusted R-squared:  0.07144 
## F-statistic:  43.2 on 2 and 1095 DF,  p-value: < 2.2e-16

summary(haridus.lm2)

## 
## Call:
## lm(formula = log(hääli_kokku) ~ C(droplevels(haridus), base = 3), 
##     data = kandidaadid2019[kandidaadid2019$haridus != "Algharidus", 
##         ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.3599 -1.0568  0.0468  0.9530  4.6607 
## 
## Coefficients:
##                                                                 Estimate
## (Intercept)                                                       3.3302
## C(droplevels(haridus), base = 3)Keskharidus (sh keskeriharidus)   1.0298
## C(droplevels(haridus), base = 3)Kõrgharidus                       1.9163
##                                                                 Std. Error
## (Intercept)                                                         0.4038
## C(droplevels(haridus), base = 3)Keskharidus (sh keskeriharidus)     0.4144
## C(droplevels(haridus), base = 3)Kõrgharidus                         0.4073
##                                                                 t value
## (Intercept)                                                       8.246
## C(droplevels(haridus), base = 3)Keskharidus (sh keskeriharidus)   2.485
## C(droplevels(haridus), base = 3)Kõrgharidus                       4.705
##                                                                 Pr(>|t|)    
## (Intercept)                                                     4.64e-16 ***
## C(droplevels(haridus), base = 3)Keskharidus (sh keskeriharidus)   0.0131 *  
## C(droplevels(haridus), base = 3)Kõrgharidus                     2.86e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.511 on 1095 degrees of freedom
## Multiple R-squared:  0.07314,    Adjusted R-squared:  0.07144 
## F-statistic:  43.2 on 2 and 1095 DF,  p-value: < 2.2e-16

Nüüd saame mudeli väljunditest välja noppida kõikide paaride võrdluste p-väärtused. Kuna kordasime testi 2 korda, siis peame usaldusnivood langetama poole võrra, mis tähendab seda, et harjumuspäraste p-väärtuste saamiseks hoopis korutame need korduste arvuga:

summary(haridus.lm1)$coefficients[2:3, "Pr(>|t|)"] *2

## C(droplevels(haridus), base = 1)Kõrgharidus 
##                                5.993082e-16 
## C(droplevels(haridus), base = 1)Põhiharidus 
##                                2.619409e-02

summary(haridus.lm2)$coefficients[2:3, "Pr(>|t|)"] *2

## C(droplevels(haridus), base = 3)Keskharidus (sh keskeriharidus) 
##                                                    2.619409e-02 
##                     C(droplevels(haridus), base = 3)Kõrgharidus 
##                                                    5.721540e-06

Ja lõpetuseks veel, kui e astmetega komakohtade arvutamine tundub tüütu, võib R-i sundida teisendama:

format(summary(haridus.lm1)$coefficients[2:3, "Pr(>|t|)"] *2, scientific = F)

## C(droplevels(haridus), base = 1)Kõrgharidus 
##                  "0.0000000000000005993082" 
## C(droplevels(haridus), base = 1)Põhiharidus 
##                  "0.0261940889199223717054"

format(summary(haridus.lm2)$coefficients[2:3, "Pr(>|t|)"] *2, scientific = F)

## C(droplevels(haridus), base = 3)Keskharidus (sh keskeriharidus) 
##                                                 "0.02619408892" 
##                     C(droplevels(haridus), base = 3)Kõrgharidus 
##                                                 "0.00000572154"

Siit nüüd saame välja lugeda, et erinevused on:

1 (kesk) ja 2 (kõrg) p = 0.0000000000000005993082
1 (kesk) ja 3 (põhi) p = 0.02619408892
2 (kõrg) ja 3 (põhi) p = 0.00000572154

Või siis selle asemel et lihtsalt korduste arvuga p-väärtuseid läbi korrutada, võib kasutada käsku p.adjust()

p.adjust(p = summary(haridus.lm1)$coefficients[2:3, "Pr(>|t|)"], n=2, method="bonferroni")

## C(droplevels(haridus), base = 1)Kõrgharidus 
##                                5.993082e-16 
## C(droplevels(haridus), base = 1)Põhiharidus 
##                                2.619409e-02

Ja teisendame tavaliseks komakohtadega arvuks:

format(p.adjust(p = summary(haridus.lm1)$coefficients[2:3, "Pr(>|t|)"], n=2, method="bonferroni"), scientific = F)

## C(droplevels(haridus), base = 1)Kõrgharidus 
##                  "0.0000000000000005993082" 
## C(droplevels(haridus), base = 1)Põhiharidus 
##                  "0.0261940889199223717054"

Tukey HSD

Teine võimalus on kasutada Tukey Honest Significant Differences testi, mis teeb kõigi kombinatsioonide võrdlused ja korrigeerib p-väärtuseid vastavalt.

R-is on Tukey testi sisendiks anova test. Anova testi tegemiseks on R-is mitu eri süntaksit, üks variant on panna lineaarne mudel anova käsu sisse, nagu me varem tegimegi.

anova(lm(y~x))
summary(aov(y~x))

Need annavad täpselt sama tulemuse. Kuid Tukey testi sisendiks peab kasutama käsku aov().

haridus.anova <- aov(log(hääli_kokku)~haridus, data=kandidaadid2019, subset = haridus != "Algharidus")
summary(haridus.anova)

##               Df Sum Sq Mean Sq F value Pr(>F)    
## haridus        2  197.3   98.64    43.2 <2e-16 ***
## Residuals   1095 2500.1    2.28                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(haridus.anova)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(hääli_kokku) ~ haridus, data = kandidaadid2019, subset = haridus != "Algharidus")
## 
## $haridus
##                                                   diff        lwr         upr
## Kõrgharidus-Keskharidus (sh keskeriharidus)  0.8865119  0.6358934  1.13713034
## Põhiharidus-Keskharidus (sh keskeriharidus) -1.0297743 -2.0022493 -0.05729919
## Põhiharidus-Kõrgharidus                     -1.9162861 -2.8721144 -0.96045787
##                                                 p adj
## Kõrgharidus-Keskharidus (sh keskeriharidus) 0.0000000
## Põhiharidus-Keskharidus (sh keskeriharidus) 0.0349434
## Põhiharidus-Kõrgharidus                     0.0000085

Tukey testi väljund annab paarikaupa võrdluses:

kahe taseme vahelise vahelise erinevuse
usaldusvahemiku (lwr ja upr)
korrigeeritud p-väärtuse

Kordamisküsimused

Kordamisküsimuste testi saad teha moodlis, seal näed peale vastamist ka õigeid vastuseid ja kommentaare. Testi võib teha mitu korda. Küsimustik on ainult kordamiseks mõeldud ja vastuseid ei hinnata ega arvestata kursuse soorituse hindamisel.

1. Lineaarse mudeli eeldused on:

mõõtmised on üksteisest sõltumatud
mõõtmisi peab olema vähemalt 30
väärtused peavad olema logaritmitud
mudeli jäägid on normaaljaotusega

2. Vaata Näites 1 esitatud lineaarset mudelit sõnade äratundmise reaktsiooniaja ja sõna pikkuse seosest:

##             Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) 547.3792  35.949589 15.226300 3.130610e-27
## Length       30.2416   4.261484  7.096495 2.289494e-10

Estimate tulbas on mudeli hinnang reaktsiooniaja kestuse kohta

esimesel real väärtus siis, kui sõna pikkus on 0 tähemärki ja teisel real siis, kui tähemärke on 1
esimesel real keskmine ja teisel real see, kui palju sõna pikkus mõjutab
esimesel real keskmine ja teisel real siis, kui tähemärke on 1
esimesel real väärtus siis, kui sõna pikkus on 0 tähemärki ja teisel real ühe tähemärgi efekt

3. Vaata veel mudeli väljundit. Mis võiks olla reaktsiooniaeg, kui sõnas on 26 tähemärki

547.38 + 26 * 30.24 = 1333.66 millisekundit
26 * 547.38 + 26 * 30.24 = 1.501814^{4} millisekundit
547.38 + 30.24 = 577.62 millisekundit
nii pikka sõna ei ole

4. Nominaalse seletava tunnusega lineaarse mudeli lõikepunkt (Intercept) on

uuritava tunnuse väärtus siis, kui seletava tunnuse väärtus on 0
uuritava tunnuse väärtus, kui ükski seletava tunnuse tasemetest ei kehti
uuritava tunnuse väärus kõige sagesasema seletava tunnuse taseme korral
uuritava tunnuse väärtus seletava tunnuse baastaseme korral

5. Mudeli headust näitab

p-väärtus
R-ruut
seletava tunnuse p-väärtus
seletava tunnuse efekti suurus

6. Artiklis, kus raporteeritakse ANOVA testi tulemusi, on tekstis kirjas: F(2, 167)=27.532; p<0.001. Kuidas seda tõlgendada?

Andmestikus oli 170 mõõtmist, seletaval tunnusel 3 rühma, vähemalt üks oli teistest erinev
Andmestikus oli 170 mõõtmist, seletaval tunnusel 3 rühma, mis üksteisest erinesid
Andmestikus oli 167 mõõtmist, seletaval tunnusel 2 rühma, mis üksteisest erinesid
P-väärtus on väike, aga seal ei ole rohkem sisukat informatsiooni

7. Post-hoc test on:

alternatiiv t-testile, kui on rohkem kui kaks rühma
kordustestid, mida peab alati tegema, kui on rohkem kui kaks rühma
kordustestid kõikide faktori tasemete võrdlemiseks
kordustestid rohkem kui kahe tasemega faktori olulisuse kontrollimiseks

Järgmisel korral

Mitme tunnusega lineaarsed mudelid
Interaktsioonid
Mudelite võrdlemine ja optimaalse mudeli leidmine

Funktsioonide sõnastik

lm(y ~ x) – lineaarne mudel
summary(lm()) – lineaarse mudeli väljund
anova(lm()) – ANOVA test
summary(aov()) – ANOVA test (sama mis eelmine)
TukeyHSD(aov()) – Tukey test (post-hoc ANOVAle)
p.adjust() – p-väärtuste korrigeerimine (Bonferroni, Holm …)
C(mingi_faktor, base = taseme_number) – vahetab faktori baastaset. NB! suure tähega C() ei ole sama käsk mis väikse tähega c().
droplevels(faktor) – jätab faktori tasemetest välja andmestikus esindamata tasemed.

11. praktikum

01.04.2024

Lineaarse regressiooni mudelid 1

Tänases praktikumis

Seletav tunnus arvuline

Mudeli eeldused

Näide 1: Kas lühemaid sõnu tuntakse kiiremini ära?

Kas mudeli eeldused on täidetud?

Sõltuv ja seletav tunnus

Kas mudeli jäägid on normaaljaotusega?

Kuidas mudelit lugeda?

Harjutus

Harjutus e Näide 2: Kas kõnetempo sõltub vanusest?

Seletav tunnus nominaalne

Kahe tasemega faktor

Näide 3: Kas kandidaadi sugu mõjutab häältesaaki (2019 riigikogu valmisite andmetes)?

Kuidas tulemust tõlgendada?

Mitme tasemega faktor

Näide 4: Kas haridus mõjutab häältesaaki?

ANOVA

ANOVA, faktor 2 tasemega

Näide 3b: Sugu ja häältesaak ANOVAga

ANOVA, faktor 3 tasemega

Näide 4b: Haridus ja häältesaak

Post-hoc test

Bonferroni parandus

Tukey HSD

Kordamisküsimused

1. Lineaarse mudeli eeldused on:

2. Vaata Näites 1 esitatud lineaarset mudelit sõnade äratundmise reaktsiooniaja ja sõna pikkuse seosest:

Estimate tulbas on mudeli hinnang reaktsiooniaja kestuse kohta

3. Vaata veel mudeli väljundit. Mis võiks olla reaktsiooniaeg, kui sõnas on 26 tähemärki

4. Nominaalse seletava tunnusega lineaarse mudeli lõikepunkt (Intercept) on

5. Mudeli headust näitab

6. Artiklis, kus raporteeritakse ANOVA testi tulemusi, on tekstis kirjas: F(2, 167)=27.532; p<0.001. Kuidas seda tõlgendada?

7. Post-hoc test on:

Järgmisel korral

Funktsioonide sõnastik