Lineaarse regressiooni mudelid 1

Tänases praktikumis

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

Kasutame andmestikke:

sõnapikkused & reaktsiooniajad ldt.csv
kõnetempo teke2 (loeme otse dataDOI repositooriumist, ei pea faili alla laadima)
välteandmestik (loeme otse dataDOI repositooriumist)

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. Niisiis, lineaarse regressioonimudeli eeldused on:

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 üle-eelmisest korrast tuttavat Levshina õpiku andmestikku ldt. Eelmisel korral testiseme kahe arvulise tunnuse vahelist seost ning 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 vastas, kuid 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 testisime me seda, kas kahe arvulise tunnuse vahel on seos, siis nüüd otsime põhjuslikku seost. Kahe arvulise tunnuse puhul võib see olla vaid 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 R-is on objektitüübi poolest list (vt mode(m) ja str(m)). Mudeli jäägid on vektorina selle listi elemendis ..$residuals 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, siin siis 547 ms. Kalle on y väärtuse muutus, kui x muutub 1 ühiku võrra. Kui sõna pikeneb ühe tähemärgi võrra, kasvab reaktsiooniaeg 30 millisekundi võrra. 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 teismeliste kõnetempo andmestikku ja vaatame, kas kõnetempo sõltub vanusest. Andmestik on pärit artiklist Lippus, Pärtel, Maarja-Liisa Pilvik, Kaidi Lõo & Liina Lindström. 2024. Kõnetempo ja -soravuse varieerumine eesti keeles. Eesti Rakenduslingvistika Ühingu aastaraamat = Estonian Papers in Applied Linguistics 20. 149–163. https://doi.org/10.5128/ERYa20.09.

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

Kuigi artikkel, kus seda andmestikku on kasutatud, uurib kõnetempo varieerumist, ei ole tabelis see valmis kujul olemas vaid tuleb arvutada. Kõnetempona mõõdetakse, mitu silpi lausub kõneleja ühe sekundi jooksul. Andmestikus on tulp kõneleja silpide arvuga ja teine tulp tema kõnevoorude kogukestusega. Arvutame kõnetempo, jagades silpide arvu 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.

Välteandmestik

Kasutame andmestikku, mis on kogutud eesti keele spontaanse kõne foneetilisest korpusest, kasutatud artiklis Lippus, P., Asu, E. L., Teras, P., & Tuisk, T. (2013). Quantity-related variation of duration, pitch and vowel quality in spontaneous Estonian. Journal of Phonetics, 41(1), 17–28. https://doi.org/10.1016/j.wocn.2012.09.005

Selles töös uuriti eesti keele välte akustilisi tunnuseid. Eesti keeles on 3 väldet (nt sada : saada : saada). Töös uuriti häälikute kestus, põhitooni liikumist ja vokaalikvaliteeti sõltuvalt vältes ja sõna lauserõhust.

Andmestik on saadaval DataDOI repositooriumis: http://dx.doi.org/10.15155/repo-16 Loeme sisse otse veebilingilt.

dat <- read.delim("http://datadoi.ut.ee/bitstream/handle/33/51/Lippus_etal_JPhon2013_dataset.txt?sequence=4&isAllowed=y")

Andmestikus on info 1714 kahesilbilise sõna kohta. Tegemist on korpusuuringuga, mis tüüpiliselt rikub tavalise lineaarse mudeli mõõtmiste sõltumatuse nõuet: meil on ühelt kõnelejalt korjatud kokku kõik uurimiskriteeriumitele vastavad sõnad ja ühelt kõnelejalt on andmestikus rohkem kui üks rida; samuti võivad andmestikus esineda korduvad sõnad.

# ühelt kõnelejalt on vähemalt 23 ja kuni 203 sõna
quantile(table(dat$SP))

##   0%  25%  50%  75% 100% 
##   23   41   64  108  203

# ainult 179 sõna esineb ühe korra, mõned sõnad esinevad 50-100 korda
table(table(dat$Word))

## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  20  21 
## 179  63  28  18  13   8   7   8   5   4   4   1   1   1   3   4   1   2   2   1 
##  23  24  25  26  31  36  38  40  50  54  56  93 117 
##   2   1   1   1   1   1   1   1   1   1   1   1   1

Et andmed oleksid tavalise lineaarse mudeliga analüüsitavad ja ei rikuks sõltumatute mõõtmiste nõuet, keskmistame andmed iga kõneleja iga uuritava tunnuse kohta:

# Kodeerime numbrina sisestatud välte ümber faktoriteks:
dat$Quantity <- factor(dat$Quantity); levels(dat$Quantity) <- paste0("Q", levels(dat$Quantity))

library(tidyverse)
# keskmistatud andmestik
dat2 <- dat %>% 
    dplyr::filter(Struct == "cvcv") %>% 
    group_by(SP, Gender, Quantity, Intonation) %>% 
    summarise(C1 = mean(C1, na.rm=T), V1 = mean(V1), C2 = mean(C3), V2 = mean(V2), 
              TP.time = mean(TP.time, na.rm=T), S1.f0.end = mean(S1.f0.end, na.rm=T),
              S1.f0.range = mean(S1.f0.range, na.rm=T), S2.f0.range = mean(S2.f0.range, na.rm=T),
              Euk1 =mean(Euk1), Euk2 = mean(Euk2))

Kahe tasemega faktor

Näide 3: Vokaali kestus ja lauserõhk

Näitena kahetasemega seletavast faktorist võiksime küsida, kas lauserõhulisus mõjutab sõna rõhulise vokaali kestust? Rõhulise vokaali kestus on tulbas “V1”, see on mõõdetud millisekundites, ja lauserõhulisus on tulbas Intonation, tasemed “deaccented” on rõhutu ja “H*L” aktsentueeritud (kõrge rõhulise silbi põhitoonile järgneb põhitooni langus).

boxplot(V1 ~ Intonation, data = dat2)

Kas mudeli eeldused on täidetud?

Mõõtmiste sõltumatusega oleme juba tegelenud.
Vokaali kestus ei ole normaaljaotusega ja on paremale kaldu, seega võiks logaritmida. Kuna mõõtmisi oli rohkem kui 30, siis võiks ka normaaljaotuse eeldust ignoreerida, aga kui logaritmimine aitab normaalajotust saavutada, siis võiks seda ikkagi teha.

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. Kui eeldaksime, et lauserõhulisus sõltub kestuses, võiks seda testida logistilise regressioonimudeliga, millest tuleb juttu umbes üle-ülejärgmises praktikumis :)

Milline tase valida 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.

inton.lm <- lm(log(V1)~Intonation, data=dat2)
summary(inton.lm)

## 
## Call:
## lm(formula = log(V1) ~ Intonation, data = dat2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.12689 -0.42499  0.07334  0.34445  0.85198 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    4.42273    0.06114  72.336   <2e-16 ***
## IntonationH*L  0.18561    0.08225   2.257    0.026 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4366 on 112 degrees of freedom
## Multiple R-squared:  0.0435, Adjusted R-squared:  0.03496 
## F-statistic: 5.093 on 1 and 112 DF,  p-value: 0.02596

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.02596 < ɑ, seega võib mudelit oluliseks pidada.
mudeli üldist headust nätab R-ruut eelviimasel real, mis ütleb, mitu protsenti uuritava tunnuse varieerumisest mudel kirjeldab. Seos on võrdlemisi nõrk, R-ruut on 0.0435 ehk lauserõhk kirjeldab ainult 4.35% vokaali kestuse varieerumisest. Seega 95,65% 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 vokaali kestus juhul, kui lauserõhk on baastasemel ehk “deaccented” ning teine rida ehk slope ütleb seda, kui palju see muutub siis, kui faktori tase on “H*L”.

Kuidas tulemust tõlgendada?

summary(inton.lm)

## 
## Call:
## lm(formula = log(V1) ~ Intonation, data = dat2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.12689 -0.42499  0.07334  0.34445  0.85198 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    4.42273    0.06114  72.336   <2e-16 ***
## IntonationH*L  0.18561    0.08225   2.257    0.026 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4366 on 112 degrees of freedom
## Multiple R-squared:  0.0435, Adjusted R-squared:  0.03496 
## F-statistic: 5.093 on 1 and 112 DF,  p-value: 0.02596

Mudel annab prognoositud väärtused: lauserõhutute vokaalide kestus on 4.42273 ühikut ja rõhuliste oma on 0.18561 võrra suurem. Ühikud on üldiselt samad mudeli sisendiks olnud ühikud, ehk hääliku kestus, 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:

# lauserõhutute vokaalide kestus on vabaliikme väärtus
exp(inton.lm$coefficients[1])

## (Intercept) 
##     83.3231

# aktsentueeritud häälikute kestuse saamiseks tuleb koefitsiendid kokku liita
exp(inton.lm$coefficients[1]+inton.lm$coefficients[2])

## (Intercept) 
##    100.3173

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

inton.lm2 <- lm(V1~Intonation, data=dat2)
summary(inton.lm2)

## 
## Call:
## lm(formula = V1 ~ Intonation, data = dat2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -64.725 -38.099  -1.825  30.522 110.848 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     91.725      5.871  15.625   <2e-16 ***
## IntonationH*L   17.427      7.897   2.207   0.0294 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 41.92 on 112 degrees of freedom
## Multiple R-squared:  0.04167,    Adjusted R-squared:  0.03311 
## F-statistic:  4.87 on 1 and 112 DF,  p-value: 0.02937

Nüüd võime koefitsientide põhjal öelda, et deaktsentueeritud vokaalid on 92 millisekundit ja aktsentueeritud on 17 millisekundit rohkem ehk 109 ms. Ja saame mudeli prognoositud väärtusi otse võrrelda aritmeetiliste keskmiste väärtustega:

tapply(dat2$V1, dat2$Intonation, mean)

## deaccented        H*L 
##   91.72461  109.15181

Mitme tasemega faktor

Näide 4: Vokaali kestus ja välde

Proovime nüüd testida seletavat tunnust, millel oleks rohkem kui kaks taset. Näiteks sobib selleks eesti keele välde, millel on kolm taset.

boxplot(V1 ~ Quantity, data = dat2)

valde.lm <- lm(log(V1)~Quantity, data=dat2)
summary(valde.lm)

## 
## Call:
## lm(formula = log(V1) ~ Quantity, data = dat2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.89277 -0.12681  0.00443  0.15610  0.60252 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.02934    0.04177   96.46   <2e-16 ***
## QuantityQ2   0.62961    0.05947   10.59   <2e-16 ***
## QuantityQ3   0.85241    0.05870   14.52   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2575 on 111 degrees of freedom
## Multiple R-squared:  0.6703, Adjusted R-squared:  0.6643 
## F-statistic: 112.8 on 2 and 111 DF,  p-value: < 2.2e-16

Mudeli väljund ütleb meile et:

mudel on erinev juhuslikust, oluline, p<0.001
seos on üsna tugev, palju tugevam, kui lauserõhu puhul: mudel kirjeldab 67% vokaali kestuse 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 esimene välde (Q1), mille väärtus on 4.03 logartimitud millisekundit ehk 56.22 millisekundit. Siin real p <0.001*** ei ole eriti tähenduslik (tähendab, et on nullist erinev);
- Teisel real on teise välte (Q2) erinevus esimese välte (Q1) kestusest ehk 4.03. Kui see oleks algsetes ühikutes, siis võiksime öelda, et see teise välte vokaal on nii mitu millisekundit pikem, aga kuna me logaritmisime väärtused enne modelleerimist, siis peame kõigepealt prognoositud tulemuse kokku arvutama ja siis astendama (käsuga exp()): teise välte logaritmitud kestus oli 4.03 + 0.63 = 4.66 ehk 105.53 millisekundit. Ja see erinevus esimese ja teise välte vahel on oluliselt erinev p<0.001.
- Kolmandal real on kolmanda välte (Q3) kestuse erinevus jällegi esimesest vältes, mis on 4.03 + 0.85 = 4.88 ehk 131.86 millisekundit. See erinevus on ka oluliselt erinev juhuslikust tasemel p <0.001.
- 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: Kestus ja lauserõhk ANOVAga

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

inton.lm <- lm(log(V1)~Intonation, data=dat2)
anova(inton.lm)

## Analysis of Variance Table
## 
## Response: log(V1)
##             Df Sum Sq Mean Sq F value  Pr(>F)  
## Intonation   1  0.971 0.97100  5.0931 0.02596 *
## Residuals  112 21.353 0.19065                  
## ---
## 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: Intonation) ja lõpuks rida mudeli jääkide kohta.
Iga faktori puhul hinnatakse rühmadevahelist dispersiooni valimi üldise dispersiooniga, millest kokku tuleb ANOVA statistik F.
F väärtuse statistiline olulisus sõltub vabadusastmetest. Faktori vabadusasmed on tasemete arv miinus üks, ehk lauserõhu puhul 2-1=1; Valimi vabadusastmete arv on mõõtmiste arv miinus võrreldavate konditsioonide arv ehk 114 (rida tabelis dat2) - 2 (deaccented/H*L) = 112.

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

t.test(log(V1)~Intonation, data=dat2)

## 
##  Welch Two Sample t-test
## 
## data:  log(V1) by Intonation
## t = -2.241, df = 103.92, p-value = 0.02715
## alternative hypothesis: true difference in means between group deaccented and group H*L is not equal to 0
## 95 percent confidence interval:
##  -0.34985792 -0.02136616
## sample estimates:
## mean in group deaccented        mean in group H*L 
##                 4.422726                 4.608338

ANOVA, faktor 3 tasemega

Näide 4b: Kestus ja välde

Nüüd näide rohkem kui kahe tasemega faktorist: vokaali kestus ja välde.

valde.lm <- lm(log(V1)~Quantity, data=dat2)
anova(valde.lm)

## Analysis of Variance Table
## 
## Response: log(V1)
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## Quantity    2 14.9633  7.4816  112.83 < 2.2e-16 ***
## Residuals 111  7.3606  0.0663                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova tulemuse põhjal võime öelda, et vältel 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, 111) = 112.83, p < 0.001
Kuna siin on rohkem kui 2 rühma, siis selleks, et otsustada, kes milline rühm millisest 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 esimese välte vokaalist oli teise välte vokaal pikem (p<0.001) ja esimese välte vokaalist oli ka kolmanda välte vokaal pikem (p<0.001):

summary(valde.lm)

ning ANOVA testi põhjal võime öelda, et välde on oluline faktor [F(2, 111) = 112.83, p < 0.001], sest leidub vähemalt üks välteaste, mille puhul vokaali kestus on erinev:

anova(valde.lm)

AGA me ei tea, kas teise ja kolmanda välte vokaalid omavahel ka erinevad on.

Post-hoc testimiseks Bonferroni korrektsiooniga võime lihtsalt korrata 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().

valde.lm1 <- lm(log(V1)~C(droplevels(Quantity), base=1), data=dat2)
valde.lm2 <- lm(log(V1)~C(droplevels(Quantity), base=3), data=dat2)
summary(valde.lm1)

## 
## Call:
## lm(formula = log(V1) ~ C(droplevels(Quantity), base = 1), data = dat2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.89277 -0.12681  0.00443  0.15610  0.60252 
## 
## Coefficients:
##                                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                          4.02934    0.04177   96.46   <2e-16 ***
## C(droplevels(Quantity), base = 1)Q2  0.62961    0.05947   10.59   <2e-16 ***
## C(droplevels(Quantity), base = 1)Q3  0.85241    0.05870   14.52   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2575 on 111 degrees of freedom
## Multiple R-squared:  0.6703, Adjusted R-squared:  0.6643 
## F-statistic: 112.8 on 2 and 111 DF,  p-value: < 2.2e-16

summary(valde.lm2)

## 
## Call:
## lm(formula = log(V1) ~ C(droplevels(Quantity), base = 3), data = dat2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.89277 -0.12681  0.00443  0.15610  0.60252 
## 
## Coefficients:
##                                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                          4.88175    0.04123  118.39  < 2e-16 ***
## C(droplevels(Quantity), base = 3)Q1 -0.85241    0.05870  -14.52  < 2e-16 ***
## C(droplevels(Quantity), base = 3)Q2 -0.22280    0.05910   -3.77 0.000263 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2575 on 111 degrees of freedom
## Multiple R-squared:  0.6703, Adjusted R-squared:  0.6643 
## F-statistic: 112.8 on 2 and 111 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(valde.lm1)$coefficients[2:3, "Pr(>|t|)"] *2

## C(droplevels(Quantity), base = 1)Q2 C(droplevels(Quantity), base = 1)Q3 
##                        3.179586e-18                        4.039967e-27

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

## C(droplevels(Quantity), base = 3)Q1 C(droplevels(Quantity), base = 3)Q2 
##                        4.039967e-27                        5.265602e-04

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

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

##   C(droplevels(Quantity), base = 1)Q2   C(droplevels(Quantity), base = 1)Q3 
## "0.000000000000000003179585974469365" "0.000000000000000000000000004039967"

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

##   C(droplevels(Quantity), base = 3)Q1   C(droplevels(Quantity), base = 3)Q2 
## "0.000000000000000000000000004039967" "0.000526560230228536240923431499539"

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

Q1 (esimene) ja Q2 (teine) p = 0.000000000000000003179585974469365
Q1 (esimene) ja Q3 (kolmas) p = 0.000000000000000000000000004039967
Q2 (teine) ja Q3 (kolmas) p = 0.0005

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(valde.lm1)$coefficients[2:3, "Pr(>|t|)"], n=2, method="bonferroni")

## C(droplevels(Quantity), base = 1)Q2 C(droplevels(Quantity), base = 1)Q3 
##                        3.179586e-18                        4.039967e-27

Ja teisendame tavaliseks komakohtadega arvuks:

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

##   C(droplevels(Quantity), base = 1)Q2   C(droplevels(Quantity), base = 1)Q3 
## "0.000000000000000003179585974469365" "0.000000000000000000000000004039967"

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

valde.anova <- aov(log(V1)~Quantity, data=dat2)
summary(valde.anova)

##              Df Sum Sq Mean Sq F value Pr(>F)    
## Quantity      2 14.963   7.482   112.8 <2e-16 ***
## Residuals   111  7.361   0.066                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(valde.anova)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(V1) ~ Quantity, data = dat2)
## 
## $Quantity
##            diff        lwr       upr     p adj
## Q2-Q1 0.6296095 0.48832342 0.7708956 0.0000000
## Q3-Q1 0.8524131 0.71297461 0.9918516 0.0000000
## Q3-Q2 0.2228036 0.08241398 0.3631932 0.0007646

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

31.03.2025

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

Välteandmestik

Kahe tasemega faktor

Näide 3: Vokaali kestus ja lauserõhk

Kuidas tulemust tõlgendada?

Mitme tasemega faktor

Näide 4: Vokaali kestus ja välde

ANOVA

ANOVA, faktor 2 tasemega

Näide 3b: Kestus ja lauserõhk ANOVAga

ANOVA, faktor 3 tasemega

Näide 4b: Kestus ja välde

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