Lineaarse regressiooni mudelid 2

Tänases praktikumis

Mitme tunnusega lineaarne regressioon
Tunnuste vahelised interaktsioonid
Parima mudeli valimine

Andmestikud endiselt:

ldt.csv
kandidaadid2019.RData
fonkorp2 (loeme otse DataDOIst, ei pea alla laadima)

Kaks arvulist tunnust

Näide 1: Kas lisaks sõna pikkusele võiks reaktsiooniaega mõjutada sõna sagedus?

Loeme sisse ldt.csv andmestiku, kust eelmisel korral vaatasime reaktsiooniaja seost sõna pikkusega. See mäletatavasti oli statistiliselt oluline seos, mis kirjeldas umbes 35% reaktsiooniaja varieerumisest.

ldt <- read.csv("data/ldt.csv")
# Jätame alles ainult read, kus Mean_RT < 1200 ms
ldt2 <- ldt[ldt$Mean_RT < 1200,]

Nüüd proovime, kas mudel läheb paremaks, kui me lisame sinna seletavaid tunnuseid. Üks sõnade äratundmiseks kuluvat reaktsiooniaega mõjutav tunnus võiks lisaks sõna pikkusele olla sõna sagedus: võiks oletada, sagedasemaid sõnu tunneme me ära kiiremini kui vähemsagedasi.

Enne mudelisse lisamist vaatame selle faktori jaotust.

hist(ldt2$Freq)

Sagedus on tugevalt paremale kaldu, seda on juba histogrammilt näha. Proovime logaritmida.

Kuna siin on 0-sagedusega sõnu ja logaritm nullist on miinus lõpmatus,

log(0)

## [1] -Inf

siis siin saab kasutada käsku log1p(), mis liidab algsele väärtusele ühe ja siis logaritmib, vältides nii -Inf (ehk miinus lõpmatus) väärtuseid.

# Võrdle!
log(1:3); log1p(0:2)

## [1] 0.0000000 0.6931472 1.0986123

## [1] 0.0000000 0.6931472 1.0986123

Testime normaaljaotust logaritmitud väärtustega:

qqnorm(log1p(ldt2$Freq))
qqline(log1p(ldt2$Freq))

shapiro.test(log1p(ldt2$Freq))

## 
##  Shapiro-Wilk normality test
## 
## data:  log1p(ldt2$Freq)
## W = 0.9794, p-value = 0.1313

Logaritmimine aitas, nüüd on tunnus normaalajotusega.

Läheme edasi mudeliga. Mudelis lisatakse sõltumatud tunnused +-märgiga. Kui enne ühe seletava tunnusega mudeli valem oli y ~ x, siis nüüd on see y ~ x1 + x2.

ldt2.mod1 <- lm(Mean_RT ~ Length + log1p(Freq), 
               data = ldt2)

Kuidas mudelit lugeda?

summary(ldt2.mod1)

## 
## Call:
## lm(formula = Mean_RT ~ Length + log1p(Freq), data = ldt2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -219.823  -58.871   -2.198   44.140  243.543 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  762.232     50.251  15.168  < 2e-16 ***
## Length        18.895      4.265   4.430 2.54e-05 ***
## log1p(Freq)  -21.348      3.894  -5.483 3.51e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 87.47 on 94 degrees of freedom
## Multiple R-squared:  0.5048, Adjusted R-squared:  0.4943 
## F-statistic: 47.91 on 2 and 94 DF,  p-value: 4.509e-15

Samamoodi nagu ühe seletava tunnusega mudeli puhul tuleks vaadata kõige pealt väljundi viimast rida, et kas mudel on oluline. See mudel on, p<0.001
Teiseks tuleks vaadata eelviimast rida, R-ruutu. See mudel seletab 50% reaktsiooniaja varieerumisest, seega on see päris hea mudel ja faktori lisamine tegi asja paremaks, sest ainult sõna pikkust arvestav mudel, mida testisime eelmises praktikumis, seletas 35%. Mudelite võrdlemisest tuleb natuke põhjalikumalt juttu peagi.

Vaatame nüüd koefitsientide tabelit (Coefficients).

Siin mudelis on kaks üldefekti, mis mõjuvad uuritavale tunnusele üksteisest sõltumata. Vabaliige on reaktsiooniaeg juhul, kui kõigi seletavate tunnuste väärtus on 0. Iga seletava tunnuse efekt on uuritava tunnuse väärtuse muutus siis, kui see seletav tunnus kasvab ühe ühiku võrra. Kui seletava tunnuse efekt on negatiivne, siis seos on vastupidine: seletav tunnus kasvab ja uuritav kahaneb.
Mudeli võime kirja panna valemina
Reaktsiooniaeg = 762 + 19 * sõna pikkus - 21 * log1p(sõna sagedus)
Kuna sagedus on mudelis logaritmitud, siis seda peame ka mudeli prognooside arvutamisel arvestama. Selle kordajaga peaks korrutama logaritmitud sagedusi.

Me võime mudelit testida, arvutades välja prognoosi seletavate tunnuste keskmiste väärtuste korral ning võrrelda seda reaalse andmestiku keskmisega. Tegelikud keskmised väärtused on:

apply(ldt2[,-1], 2, mean)

##      Length        Freq     Mean_RT 
##    8.092784 3452.103093  792.117938

paneme need valemisse
RT = 762 + 19 * pikkus - 21 * log1p(sagedus)

762.232 + 18.895 * 8.09 + -21.348 * log1p(3452.10)

## [1] 741.1698

Mudeli prognoose võime muidugi ka otse mudelist lugeda ja valemites rakendada. Vaatame lm() käsuga saadud objekti, et saada natuke aimu, mis struktuuriga objekt see on:

str(ldt2.mod1); names(ldt2.mod1$coefficients)

Käsk str() aitab ehk kõige paremini lm objekti lahti muukida: tegemist on listilaadse objektiga, kus element “coefficients” on vektor mudeli prognooside väärtustega. Kasutame nüüd neid koefitsentide väärtusi, et arvutada mudeliprognoose. Nii saame arvutada, mis võiks olla reaktsiooniaeg keskmise pikkusega keskmise sagedusega sõna korral:

ldt2.mod1$coefficients["(Intercept)"] + 
  mean(ldt2$Length) * ldt2.mod1$coefficients["Length"] + 
  log1p(mean(ldt2$Freq)) * ldt2.mod1$coefficients["log1p(Freq)"]

## (Intercept) 
##    741.2168

Joonistame mudeli põhjal regressioonijooned nii, et teise tunnuse väärtus oleks keskmine, mitte 0.

par(mfcol = c(1,2)) # Kuva järgmised joonised ühes reas, kahes tulbas
plot(Mean_RT ~ Length, data = ldt2, 
     xlab = "Sõna pikkus", ylab = "Keskmine reaktsioonikiirus", 
     main = "Sõna pikkuse mõju\nkeskmise sagedusega sõnadel")
lines(x = 3:14, 
      y = (ldt2.mod1$coefficients["(Intercept)"] + 
               3:14 * ldt2.mod1$coefficients["Length"] +
               log1p(mean(ldt2$Freq)) * ldt2.mod1$coefficients["log1p(Freq)"]), 
      col = "red")

plot(Mean_RT ~ log1p(Freq), data = ldt2, 
     xlab = "Sõna sagedus", ylab = "Keskmine reaktsioonikiirus", 
     main = "Sõna sageduse mõju\nkeskmise pikkusega sõnadel", 
     axes=F)
lines(x = 0:12, 
      y = (ldt2.mod1$coefficients["(Intercept)"] + 
               mean(ldt2$Length) * ldt2.mod1$coefficients["Length"] + 
               0:12 * ldt2.mod1$coefficients["log1p(Freq)"]), 
      col = "red")
# Samuti arvestame sellega, et sagedus oli logaritmitud, aga joonisel tahaks päris väärtusi näha.
box(); axis(side = 2); axis(side = 1, at = 0:12, labels = round(exp(0:12)-1))

par(mfcol = c(1, 1))

Sama ggplotiga:

library(ggplot2)
library(gridExtra) # pakett, millega saab erinevaid jooniseid paneelidena ühe joonise peale kokku panna
paneel1 <- ggplot(data = ldt2, aes(x = Length, y = Mean_RT)) +
    geom_point(shape = 1) +
    geom_abline(intercept = ldt2.mod1$coefficients["(Intercept)"] + 
                    log1p(mean(ldt2$Freq)) * ldt2.mod1$coefficients["log1p(Freq)"],
                slope = ldt2.mod1$coefficients["Length"],
                color = "red") +
    labs(x = "Sõna pikkus", y = "Keskmine reaktsioonikiirus",
         title = "Sõna pikkuse mõju keskmise sagedusega sõnadel")


paneel2 <- ggplot(data = ldt2, aes(x = log1p(Freq), y = Mean_RT)) +
    geom_point(shape = 1) +
    geom_abline(intercept = ldt2.mod1$coefficients["(Intercept)"] + 
                    mean(ldt2$Length) * ldt2.mod1$coefficients["Length"], 
                slope = ldt2.mod1$coefficients["log1p(Freq)"],
                color = "red") +
    scale_x_continuous(breaks = seq(0,12,1),  
                       labels = round(exp(seq(0,12,1)))-1) + # Kuva tavalised sagedused
    labs(x = "Sõna sagedus", y = "Keskmine reaktsioonikiirus",
         title = "Sõna sageduse mõju keskmise pikkusega sõnadel")
grid.arrange(paneel1, paneel2, ncol=2)

Kui tahame sageduse skaalat joonisele lineaarsena (st x-telje intervallid on ühesuurused), siis tuleb regressioonijoon joonistada kõverana.

plot(Mean_RT  ~Freq, data = ldt2, 
     xlab="Sõna sagedus", ylab="Keskmine reaktsioonikiirus", main="Sõna sagedus", 
     xlim=c(0, 5000))

lines(x = 0:5000, y = (ldt2.mod1$coefficients["(Intercept)"] + 
                           mean(ldt2$Length) * ldt2.mod1$coefficients["Length"] + 
                           log1p(0:5000) * ldt2.mod1$coefficients["log1p(Freq)"]), 
      col="red")
box(); axis(side=2); axis(side=1, at=0:12, labels = round(exp(0:12)-1))

Mudeleid saab visualiseerida ka ggploti laiendavate pakettide abil.

# install.packages("ggeffects")
library(ggeffects)
# plot(ggpredict(model = ldt2.mod1), facets = T, add.data = T) + labs(x = "") # see andis veateadet millegi pärast, paketis on vist midagi muutunud?
paneel1 = plot(ggpredict(model = ldt2.mod1), add.data = T)$Length
paneel2 = plot(ggpredict(model = ldt2.mod1), add.data = T)$Freq
grid.arrange(paneel1, paneel2, ncol=2)

# install.packages("sjPlot")
library(sjPlot)
plot_model(ldt2.mod1, type = "pred", grid = T) + labs(x = "") #[Maarja-Liisa, see annab veateadet]

Väike lõbus vahepala: interaktiivsed joonised

Paketiga plotly saab ggploti joonised teha interaktiivseks, nii et RStuudios Viewer aknas või RMarkdownist genereeritud html-dokumendis saab joonisel hiirega üle sõites näha punktide väärtuseid ning saab sisse ja välja zoomida:

# install.packages("plotly")
library(plotly)
ggplotly(paneel1)

Ja kaks paneeli kõrvuti:

subplot(ggplotly(paneel1), ggplotly(paneel2))

Peamõju ja interaktsioon

Peamõju ehk üldefekt toimib mitme tunnusega mudelis teistest tunnustest sõltumatuna. Näiteks kui reaktsiooniaeg sõltub sõna pikkusest ja sagedusest, nagu viimases näites, siis peaks sõna pikkuse efekt olema alati sama suur sõltumata sõna sagedusest ja sõna sageduse efekt alati sama suur sõltumata sõna pikkusest.

Kuid seletavate tunnuste vahel võib ka esineda koosmõju ehk interaktsioon. See tähendab, et ühe tunnuse mõju mõjutab teise tunnuse mõju. Näiteks pikemaid sõnu tuntakse ära aeglasemalt ja sagedasemaid sõnu kiiremini, aga a) pikki sagedasemaid sõnu tuntakse ära sama kiiresti kui lühikesi sagedasi sõnu, või b) lühikesi sagedasi sõnu tuntakse ära eriti kiiresti. Sellisel juhul on lisaks kummagi tunnuse peamõjule ka tunnuste vahel interaktsioon.

Peamõju ja interaktsioon lm() funktsioonis:

Mudeli süntaksis eraldatakse erinevad mõjud plussiga: y ~ x1 + x2
Interaktsiooni märkimiseks koolon: y ~ x1 + x2 + x1:x2
Tärniga saab tähistada lühidalt faktorite vahelist peamõju ja interaktsiooni korraga, tulemus on sama, mis eelmisel real: y ~ x1 * x2

Kõik peamõjude kombinatsioonid võivad anda interaktsiooni:

Peamõjusid võib olla ka rohkem kui kaks. Ja sel juhul võvad interaktsioonid olla ka kõikide kombinatsioonide vahel:

2 faktorit, 1 kahene interaktsioon: y ~ x1 + x2 + x1:x2
3 faktorit, 3 kahest ja üks kolmene interaktsioon: y ~ x1 + x2 + x3 + x1:x2 + x1:x3 + x2:x3 + x1:x2:x3
4 faktorit, 6 kahest, 4 kolmest ja 1 neljane interaktsioon:
y ~ x1 + x2 + x3 + x4 +
x1:x2 + x1:x3 + x1:x4 + x2:x3 + x2:x4 + x3:x4 +
x1:x2:x3 + x1:x2:x4 + x1:x3:x4 + x2:x3:x4 +
x1:x2:x3:x4
jne. Mida rohkem faktoreid, seda rohkem kombinatsioone ja seda keerukam mudel.

Kaks arvulist tunnust interaktsiooniga

Näide 2: Kas sõna pikkus ja sagedus mõjutavad reaktsiooniaega?

Kui enne vaatasime ainult peamõjudega sõna pikkuse ja sageduse efekti reaktsiooniajale, siis nüüd lisame ka interaktsiooni. R-i lm() süntaksis võiks selleks kasutada lihtsalt tärni (*), aga kirjutame siin praegu pikalt välja:

ldt2.mod2 <- lm(Mean_RT ~ Length + log1p(Freq) + Length:log1p(Freq), 
               data=ldt2)
summary(ldt2.mod2)

## 
## Call:
## lm(formula = Mean_RT ~ Length + log1p(Freq) + Length:log1p(Freq), 
##     data = ldt2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -220.071  -59.205   -0.926   44.653  243.970 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        771.2282    94.2795   8.180 1.42e-12 ***
## Length              17.8945     9.8342   1.820   0.0720 .  
## log1p(Freq)        -22.6969    12.5576  -1.807   0.0739 .  
## Length:log1p(Freq)   0.1575     1.3934   0.113   0.9103    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 87.93 on 93 degrees of freedom
## Multiple R-squared:  0.5049, Adjusted R-squared:  0.4889 
## F-statistic: 31.61 on 3 and 93 DF,  p-value: 3.545e-14

Mudeli väljundisse lisandub rida Length:log1p(Freq). Antud juhul see interaktsiooni efekt ei ole oluline, p on suurem kui 0.05.

Kuidas interaktsioone tõlgendada?

Kuigi siin näites ei osutunud interaktsioon oluliseks, siis selle mudeli põhjal peaks reaktsiooni aja prognoose arvutama nii:

RT = 771 + 18*pikkus + -23*sagedus + 0.16*pikkus*sagedus

Ja see tähendaks, et reaktsiooniaeg on:

pikem pikemate sõnadega (kasvab 18 ms sõna igatähemärgi kohta),
lühem sagedate sõnadega (kahaneb 23 ms logaritmitud sagedusühiku kohta), aga
natuke pikem pikkade sagedate sõnadega (täiendavalt pikeneb 0.16 ms iga pikkuse ja sageduse ühiku kohta).

Peamõjud ilma interaktsioonita:

Kui faktorite vahel ei ole interaktsioone, siis on ühe tunnuse regressioonijooned teise tunnuse erinevate tasemete puhul paralleelsed.

Peamõjud interaktsiooniga

Kuna interaktsioon mudelis ldt2.mod2 ei olnud oluline, siis selleks et interaktsiooni efekti illustreerida, toome paar fiktiivset näidet: Esiteks oletame, et pikkuse ja sageduse vahel on positiivne interaktsioon, näiteks:

RT = 700 + 15*pikkus + -25*sagedus + 3*pikkus*sagedus

Niisiis, positiivse interaktsiooni mõjul on madala sagedusega sõnadel pikkusel väiksem efekt (must joon vasakul), aga sagedasemate sõnade korral on pikkuse efekt tugevam (järsem kalle, helesinine joon). Sõna sagedusele, mille peaefekt oli negatiivne, mõjub interaktsioon vastupidiselt.

Ja kui interaktsioon oleks negatiivne:

RT = 1000 + 15*pikkus + -25*sagedus + -3*pikkus*sagedus

Aga veelkord, need on fiktiivsed näited, mistõttu ei lähe jooned joonistel ka tegelike andmepunktidega päriselt kokku :)

Mudelite võrdlemine

Occami habemenoa printsiip

Occam’s razor ehk __Occami habemenoa printsiip (https://et.wikipedia.org/wiki/Ockhami_habemenuga)__ ütleb, et parim seletus on kõige lihtsam seletus. Occami habemenoa printsiibile toetudes tuleks alati mudelist välja visata kõik, mis ei ole oluline, et mudel oleks võimalikult lihtne.

Mudelite võrdlemiseks on mitu üksteist täiendavat meetodit: võib lihtsalt vaadata, et ei oleks ebaolulisi tunnuseid mudelis, võib vaadata mudeli seletusvõimet R-ruudu näol, võib võrrelda erinevusi mudeli jääkides või teha spetsiaalselt mudelite headuse võrdlemiseks mõeldud teste (mis üldiselt toetuvad ka mudeli jääkide võrdlemisele).

Võrdleme mudeleid: R-ruut, ANOVA ja AIC

R-ruut ehk determinatsioonikordaja näitab protsentuaalselt, kui palju uuritava tunnuse varieerumisest mudel kirjeldab. Parem mudel on see, mis rohkem kirjeldab, mille R-ruut on suurem.

R-ruut läheb alati natuke paremaks, kui rohkem seletavaid tunnuseid lisada. Ka siis, kui seletavad tunnused on marginaalse efektiga. Occami habemenoa printsiipi aitab järgida kohandatud R-ruut, mis vähendab R-ruudu väärtust iga seletava tunnuse lisamisel: täiendava tunnuse lisamisest saadav efekt mudeli seletusvõimele peab olema suurem, kui “karistus” selle lisamise eest.

Esimese, ilma interaktsioonideta mudeli R-ruut on 0.50481 ja kohandatud R-ruut 0.49428.
Teise, interaktsioonidega mudeli R-ruut on 0.50488 ja kohandatud R-ruut 0.48891.

Seega interaktsiooni lisamisest üldine mudeli seletusvõime kasvas küll 0,007%, aga see on liiga vähe, et ühte lisafaktorit mudelisse lisada, mistõttu kohandatud R-ruut langes. Parem mudel on seega ilma interaktsioonita mudel.

Kohandatud R-ruutu vaatame siis, kui mudeleid võrdleme. Siis, kui oleme parima mudeli välja valinud ja raporteerime, kui hea kirjeldusvõimega mudelit esitame, siis peaks vaatama tavalist R-ruutu.

ANOVA test mudeli jääkidele

Anova testiga saab ka võrrelda mudeleid, täpsemalt mudeli jääkide (residuals) jaotust. Niimoodi peaks aga võrdlema samal andmestikul põhinevaid ühe faktori võrra erinevaid mudeleid.

anova(ldt2.mod1, ldt2.mod2)

## Analysis of Variance Table
## 
## Model 1: Mean_RT ~ Length + log1p(Freq)
## Model 2: Mean_RT ~ Length + log1p(Freq) + Length:log1p(Freq)
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     94 719165                           
## 2     93 719066  1    98.759 0.0128 0.9103

Kui kahe mudeli jääkides ei ole olulist erinevust, siis peaks valima lihtsama mudeli. Kui erinevus on oluline, tähendab see, et mudeli seletusvõime muutus faktori lisamisest oluliselt ja peaks valima keerukama mudeli. Siin võrreldud mudelite ANOVA tulemus on p=0.9 ehk olulist erinevust mudelite vahel pole ja peaksime neist valima lihtsama.

Akaike informatsioonikriteerium

AIC ehk Akaike informatsioonikriteerium on alternatiiv R-ruudule hinnata mudeli headust. Erinevalt R-ruudust on tulemus lihtsalt üks number, mis ei paigutu mingile skaalale ja on võrreldav ainult sama andmestiku pealt tehtud mudelite piires. Mida väiksem AIC, seda parem mudel.

AIC(ldt2.mod1, ldt2.mod2)

##           df      AIC
## ldt2.mod1  4 1147.654
## ldt2.mod2  5 1149.641

Inkrementaalne modelleerimine

Lähtudes printsiibist, et parim mudel on võimalikult lihtne, peaks mudelist kõik eba- või väheolulise välja viskama. Seda peaks tegema sammhaaval: valima kõige keerukama interaktsiooni ja/või väiksema efektiga tunnuse ning vaatama, kas seda saab välja visata. Toimingut tuleks korrata seni, kuni jõuame optimaalse mudelini.

Sageli kasutatakse ka hoopis inkrementaalset meetodit, st alustatakse 0-mudelist, kus ei ole ühtegi seletavat faktorit, ning hakatakse ükshaaval lisama tunnuseid ja nende interaktsioone seni, kuni saavutatakse optimaalne mudel. Mõlemal juhul on eesmärk leida võimalikult lihtne, aga samas kõiki olulisi faktoreid sisaldav mudel.

Näide 3: täiskasvanud kõnelejate kõnetempot mõjutavad tegurid

Kasutame Lippus, Pilvik, Lõo & Lindström 2024 kõnetempo repositooriumist täiskasvanute spontaanse kõne andmestikku fonkorp2:

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

Arvutame kõneleja ja vestluspartneri kõnetempo:

fonkorp2$konetempo <- fonkorp2$silpide_arv/fonkorp2$kestus
fonkorp2$kaaskoneleja_konetempo <- fonkorp2$kaaskoneleja_silpide_arv/fonkorp2$kaaskoneleja_kestus

Testime kõnetempo varieerumist (uuritava tunnusena). Alustame 0-mudelist ja lisame kõigepealt kõneleja vanuse, soo ja nende interaktsiooni.

mod0 <- lm(konetempo ~ 1, data = fonkorp2)
mod1a <- lm(konetempo ~ vanus, data = fonkorp2)
mod1b <- lm(konetempo ~ sugu, data = fonkorp2)

Testime 0-mudelit ühe seletava tunnusega muutuja vastu:

anova(mod0, mod1a)

## Analysis of Variance Table
## 
## Model 1: konetempo ~ 1
## Model 2: konetempo ~ vanus
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1    163 82.776                           
## 2    162 81.938  1   0.83762 1.6561    0.2

anova(mod0, mod1b)

## Analysis of Variance Table
## 
## Model 1: konetempo ~ 1
## Model 2: konetempo ~ sugu
##   Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
## 1    163 82.776                                
## 2    162 77.631  1    5.1445 10.736 0.001286 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova testi põhjal võiks öelda praegu, et sugu on oluline kõnetempo mõjutaja, vanus mitte. Aga vaatame edasi:

mod2 <- lm(konetempo ~ sugu+vanus, data = fonkorp2)
mod3 <- lm(konetempo ~ sugu*vanus, data = fonkorp2)

anova(mod1b, mod2)

## Analysis of Variance Table
## 
## Model 1: konetempo ~ sugu
## Model 2: konetempo ~ sugu + vanus
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1    162 77.631                           
## 2    161 76.605  1    1.0266 2.1576 0.1438

anova(mod2, mod3)

## Analysis of Variance Table
## 
## Model 1: konetempo ~ sugu + vanus
## Model 2: konetempo ~ sugu * vanus
##   Res.Df    RSS Df Sum of Sq     F  Pr(>F)  
## 1    161 76.605                             
## 2    160 75.210  1    1.3947 2.967 0.08691 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ANOVA testi põhjal jääks vanus ikkagi välja. Aga vaatame ka Aikaike informatsioonikriteeriumit. AIC käsule võime sisse panna kõik (samal andmestikul põhinevad) mudelid korraga. Võidab see, kellel on kõige väiksem AIC väärtus.

AIC(mod0, mod1a, mod1b, mod2, mod3)

##       df      AIC
## mod0   2 357.2804
## mod1a  3 357.6124
## mod1b  3 348.7573
## mod2   4 348.5740
## mod3   5 347.5608

AIC põhjal on parim soo ja vanuse interaktsiooni sisaldav mudel. Vaatame nüüd selle väljundit:

summary(mod3)

## 
## Call:
## lm(formula = konetempo ~ sugu * vanus, data = fonkorp2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.72781 -0.45981  0.02138  0.47721  2.00232 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.070240   0.245400  20.661   <2e-16 ***
## suguN       -0.166092   0.325352  -0.510   0.6104    
## vanus       -0.013951   0.006147  -2.269   0.0246 *  
## suguN:vanus  0.013830   0.008029   1.722   0.0869 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6856 on 160 degrees of freedom
## Multiple R-squared:  0.0914, Adjusted R-squared:  0.07436 
## F-statistic: 5.365 on 3 and 160 DF,  p-value: 0.001516

Mudeli väljundist näeme, et soo efekt nüüd on ebaoluline, aga on vanuse efekt on oluline ning peaaegu-et-oluline on ka soo ja vanuse interaktsioon.

Kuna sugu on faktor, siis Intercept ja teiste tunnuste peaefektid kehtivad soo baastaseme ehk M korral. Seega vanuse efekt meestel on -0.014 silp/sek aasta kohta. Naistel peaks vanuse efekti ja interaktsiooni kokku liitma ja -0.014 + 0.014 = 0

joonis <- ggplot(fonkorp2, aes(y = konetempo, x = vanus, col = sugu))+
    geom_point()+
    geom_smooth(method = "lm")
ggplotly(joonis)

Multikollineaarsus

Kui meil on rohkem sõltumatuid tunnuseid, mida mudelis testida, on oluline kontrollida, et sõltumatud tunnused ei oleks omavahel liiga tugevasti seotud, sest mitu väga sarnaselt käituvat faktorit teevad mudeli kehvemaks.

Kui andmestikus on mitu sama tunnuse teisendit (nt sünniaasta ja vanus), siis peaks valima neist ainult ühe.
Sarnaselt käituvate tunnuste puhul tuleks samuti valida üks, mis kõige paremini seletab.
Võib testida tunnustevahelist korrelatsiooni.
VIF ehk Variance Inflation Factor (nt car::vif(model1))

cov2cor(cov(fonkorp2[,c("kaaskoneleja_konetempo","vanus","kaaskoneleja_vanus", "vanusevahe")]))

##                        kaaskoneleja_konetempo      vanus kaaskoneleja_vanus
## kaaskoneleja_konetempo              1.0000000  0.1122187         -0.1005938
## vanus                               0.1122187  1.0000000          0.2277180
## kaaskoneleja_vanus                 -0.1005938  0.2277180          1.0000000
## vanusevahe                         -0.1712356 -0.6214024          0.6214024
##                        vanusevahe
## kaaskoneleja_konetempo -0.1712356
## vanus                  -0.6214024
## kaaskoneleja_vanus      0.6214024
## vanusevahe              1.0000000

Vanusevahe on nii kõneleja kui vestluspartneri vanusega tugevas korrelatsioonis, seega mudelisse võiks võtta kas vanusevahe või kummagi kõneleja vanused. Muud tunnused on nõrgalt korreleeruvad ja võivad koos mudelis olla.

Stepwise protseduur

Selle asemel, et ise ükshaaval mudelisse faktoreid juurde lisada või ära võtta ja vahesamme testida, võib seda lasta ka R-il automaatselt teha.

Käsk step() kasutab Akaike informatsioonikriteeriumit (AIC), et ehitada samm-sammuliselt faktoreid kas lisades või eemaldades optimaalne mudel. Proovime leida faktoreid lisades parima mudeli täiskasvanud kõnelejate kõnetempo seletamiseks.

Alustame 0-mudelist ja liigume keerukama suunas. Nullmudeli valemiks on y ~ 1, see läheb käsu step() objektiks.
Teine argument on scope, kuhu praegusel juhul paneme maksimaalse mudeli, st kõikide faktorite ja nendevaheliste interaktsioonidega.
Kolmas argument on direction, mis määrab faktorite lisamise suuna. Variandid on “both”, “backward”, “forward”, siit praegu valime viimase, mis tähendab, et alustatakse kõige väiksemast mudelist ja lisatakse sellele faktoreid seni, kuni mudel ei lähe enam paremaks. “Backward” tähendaks, et alustatakse maksimaalsest ja faktoreid võetakse vähemaks seni, kuni mudel ei lähe enam paremaks, ja nii võib tulemus tulla veidi erinev.

Tulemuseks saadud mudelis ei pruugi kõik faktorid olla olulised p-väärtuse poolest, sest mudeleid võrreldakse ainult AIC põhjal.

Testime seletavate tunnustena kõneleja vanust ja sugu ning vestluspartneri vanust, sugu ja kõnetempot ning kõigi tunnuste vahelist interaktsiooni.

mod_step <- step(lm(konetempo ~ 1, data = fonkorp2), scope = konetempo ~ vanus * sugu * kaaskoneleja_konetempo * kaaskoneleja_sugu * kaaskoneleja_vanus, direction = "forward")

## Start:  AIC=-110.13
## konetempo ~ 1
## 
##                          Df Sum of Sq    RSS     AIC
## + kaaskoneleja_konetempo  1   17.6947 65.081 -147.57
## + sugu                    1    5.1445 77.631 -118.66
## + kaaskoneleja_vanus      1    1.0424 81.734 -110.21
## <none>                                82.776 -110.13
## + vanus                   1    0.8376 81.938 -109.80
## + kaaskoneleja_sugu       1    0.6237 82.152 -109.37
## 
## Step:  AIC=-147.57
## konetempo ~ kaaskoneleja_konetempo
## 
##                      Df Sum of Sq    RSS     AIC
## + sugu                1    3.6490 61.432 -155.04
## + kaaskoneleja_vanus  1    2.1068 62.974 -150.97
## + vanus               1    1.9490 63.132 -150.56
## <none>                            65.081 -147.57
## + kaaskoneleja_sugu   1    0.0715 65.010 -145.75
## 
## Step:  AIC=-155.04
## konetempo ~ kaaskoneleja_konetempo + sugu
## 
##                               Df Sum of Sq    RSS     AIC
## + vanus                        1   2.13360 59.299 -158.83
## + kaaskoneleja_vanus           1   1.64084 59.791 -157.48
## <none>                                     61.432 -155.04
## + kaaskoneleja_sugu            1   0.20370 61.229 -153.58
## + sugu:kaaskoneleja_konetempo  1   0.15303 61.279 -153.45
## 
## Step:  AIC=-158.83
## konetempo ~ kaaskoneleja_konetempo + sugu + vanus
## 
##                                Df Sum of Sq    RSS     AIC
## + kaaskoneleja_vanus            1   2.82343 56.475 -164.83
## + vanus:sugu                    1   1.29091 58.008 -160.44
## <none>                                      59.299 -158.83
## + vanus:kaaskoneleja_konetempo  1   0.64797 58.651 -158.64
## + kaaskoneleja_sugu             1   0.13930 59.159 -157.22
## + sugu:kaaskoneleja_konetempo   1   0.12842 59.170 -157.19
## 
## Step:  AIC=-164.83
## konetempo ~ kaaskoneleja_konetempo + sugu + vanus + kaaskoneleja_vanus
## 
##                                             Df Sum of Sq    RSS     AIC
## + vanus:sugu                                 1   1.26531 55.210 -166.55
## <none>                                                   56.475 -164.83
## + vanus:kaaskoneleja_konetempo               1   0.47129 56.004 -164.21
## + kaaskoneleja_sugu                          1   0.21095 56.264 -163.45
## + kaaskoneleja_konetempo:kaaskoneleja_vanus  1   0.11514 56.360 -163.17
## + sugu:kaaskoneleja_vanus                    1   0.07597 56.399 -163.06
## + vanus:kaaskoneleja_vanus                   1   0.03302 56.442 -162.93
## + sugu:kaaskoneleja_konetempo                1   0.00148 56.474 -162.84
## 
## Step:  AIC=-166.55
## konetempo ~ kaaskoneleja_konetempo + sugu + vanus + kaaskoneleja_vanus + 
##     sugu:vanus
## 
##                                             Df Sum of Sq    RSS     AIC
## <none>                                                   55.210 -166.55
## + vanus:kaaskoneleja_konetempo               1   0.53426 54.676 -166.15
## + kaaskoneleja_sugu                          1   0.11512 55.095 -164.89
## + kaaskoneleja_konetempo:kaaskoneleja_vanus  1   0.09467 55.115 -164.83
## + vanus:kaaskoneleja_vanus                   1   0.03828 55.172 -164.66
## + sugu:kaaskoneleja_konetempo                1   0.02604 55.184 -164.63
## + sugu:kaaskoneleja_vanus                    1   0.00024 55.210 -164.55

summary(mod_step)

## 
## Call:
## lm(formula = konetempo ~ kaaskoneleja_konetempo + sugu + vanus + 
##     kaaskoneleja_vanus + sugu:vanus, data = fonkorp2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.79108 -0.37142  0.02127  0.35300  1.67235 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             2.602218   0.389497   6.681 3.84e-10 ***
## kaaskoneleja_konetempo  0.486138   0.066242   7.339 1.07e-11 ***
## suguN                  -0.217945   0.280627  -0.777 0.438536    
## vanus                  -0.018587   0.005373  -3.459 0.000696 ***
## kaaskoneleja_vanus      0.010022   0.003542   2.830 0.005265 ** 
## suguN:vanus             0.013174   0.006923   1.903 0.058871 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5911 on 158 degrees of freedom
## Multiple R-squared:  0.333,  Adjusted R-squared:  0.3119 
## F-statistic: 15.78 on 5 and 158 DF,  p-value: 1.374e-12

Kuidas mudelit tõlgendada? Alustama peaks olulistest peaefektidest.

kõnetempot mõjutab kaaskõneleja kõnetempo: kui vestluspartner on 1 silp/sek kiirem, siis on ka kõneleja keskmiselt 0.486 silpi/sek kiirem.
vanusel on negatiivne efekt: kui vanus kasvab 1 aasta võrra, on kõneleja vanus keskmiselt -0.019 aeglasem, aga siin on ka (peaaegu-et-oluline) interaktsioon: kui kõneleja on naine, siis vanuse kasvades 1 aasta võrra muutub kõnetempo ainult -0.019 + 0.013 = -0.005 võrra.
lisaks on vestluspartneri vanuse efekt: kui vestluspartneri vanus kasvab 1 aasta võrra, kasvab kõneleja kõnetempo 0.01 võrra.

Ja lõpetuseks teeme paar rehkendust selle mudeli põhjal:

Mis võiks olla minu kõnetempo, kui olen 44-aastane mees ja räägin 25-aastase kiire kõnelejaga (3.5 silpi/sek)?

 2.60221810 + # vabaliige on alati sama
    3.5 * 0.48613842 + # vestluspartneri tempo 3.5
    0 * -0.21794464 + # minu sugu M ehk M == N on FALSE ehk 0
    44 * -0.01858738 + # minu vanus 44
    25 * 0.01002157 + # kaasvestleja vanus 25
    0 * 44 * 0.01317365 # on kokku 0, interaktsioon ei rakendu, sest sugu pole N

## [1] 3.736397

Mis võiks olla minu kõnetempo, kui olen 80-aastane naine ja räägin 44-aastase aeglase kõnelejaga (5.5 silpi/sek)?

 2.60221810 + # vabaliige
    5.5 * 0.48613842 + # vestluspartneri tempo 5.5
    1 * -0.21794464 + # minu sugu N ehk N == N on TRUE ehk 1
    80 * -0.01858738 + # minu vanus 80
    44 * 0.01002157 + # kaasvestleja vanus 44
    1 * 80 * 0.01317365 # N == N on TRUE ehk 1 * minu vanus 80

## [1] 5.065885

Näide 4: Mis mõjutab valimistulemust?

Loeme valimistulemuste andmestiku R-i.

load("data/kandidaadid2019.RData")

Vaatasime eelmisel korral eraldi soo ja vanuse efekti häältesaagile ja leidsime, et mõlemal faktoril on väike, aga oluline efekt. Nüüd võiksime ühes mudelis testida kõiki võimalikke faktoreid ja välja selgitada, millised kirjeldavad häältesaaki kõige paremini. Andmestikus on sellised faktorid, mille seosed on andmestikku lisatud eeldusel, et neil on ehk hääletustulemusega mingi seos:

colnames(kandidaadid2019)

##  [1] "ringkond"               "nimekiri"               "nimi"                  
##  [4] "sünniaeg"               "sugu"                   "tähtkuju"              
##  [7] "haridus"                "erakond"                "amet"                  
## [10] "kontakt"                "kandidaadi_link"        "eesnimi"               
## [13] "nr"                     "hääli_kokku"            "e.hääli"               
## [16] "hääli_välisriigist"     "kokku_ringkonnas"       "hääli_kokku_prots_ring"

Ettevalmistused suuremaks mudeliks

Enne suurema mudeli juurde asumist tuleks natuke ettevalmistusi teha ja võimalikud seletavad tunnused üle vaadata.

Multikollineaarsus

Kui andmestikus on mitu sama tunnuse teisendit (nt sünniaasta ja vanus), siis peaks valima neist ainult ühe. Sarnaselt käituvate tunnuste puhul tuleks samuti valida üks, mis kõige paremini seletab (nt ei ole mõtet kandidaatide mudelisse lisada nii parteilist kuuluvust kui valimisnimekirja).

cov2cor(cov(kandidaadid2019[,c("e.hääli", "hääli_välisriigist", "kokku_ringkonnas","nr")]))

##                        e.hääli hääli_välisriigist kokku_ringkonnas          nr
## e.hääli             1.00000000         0.57115821       0.06475214 -0.06060796
## hääli_välisriigist  0.57115821         1.00000000       0.01377868  0.02979023
## kokku_ringkonnas    0.06475214         0.01377868       1.00000000 -0.04607272
## nr                 -0.06060796         0.02979023      -0.04607272  1.00000000

Vaatame enne mudeldamist üle ülejäänud tunnused

Kandidaatide tabelis on erakondlik_kuuluvus ja nimekiri_voi_uksikkandidaat väga sarnased, sest enamik kandidaate riigikogu valimistel kandideerib oma erakonna nimekirjas, seega tasuks valida neist üks:

table(kandidaadid2019$erakond) # üks ühene rühm

## 
##                      Eesti Keskerakond      Eesti Konservatiivne Rahvaerakond 
##                                    116                                    117 
##                   Eesti Reformierakond Eesti Vabaduspartei - Põllumeeste Kogu 
##                                    125                                      1 
##                      Eesti Vabaerakond         Eestimaa Ühendatud Vasakpartei 
##                                     87                                     11 
##                     Elurikkuse Erakond                      Erakond Eesti 200 
##                                     70                                     79 
##             Erakond Eestimaa Rohelised                         Isamaa Erakond 
##                                     59                                    112 
##           Sotsiaaldemokraatlik Erakond 
##                                    102

table(kandidaadid2019$nimekiri) # valime selle

## 
##                 Eesti Keskerakond Eesti Konservatiivne Rahvaerakond 
##                               125                               125 
##              Eesti Reformierakond                 Eesti Vabaerakond 
##                               125                               125 
##    Eestimaa Ühendatud Vasakpartei                Elurikkuse Erakond 
##                                11                                73 
##                 Erakond Eesti 200        Erakond Eestimaa Rohelised 
##                               125                               125 
##                    Isamaa Erakond      Sotsiaaldemokraatlik Erakond 
##                               125                               125 
##                  Üksikkandidaadid 
##                                15

table(kandidaadid2019$haridus) # eemalda ühene rühm
table(kandidaadid2019$tähtkuju)

Kandidaatide töökohtadest oleks vist liiga keeruline mingit selgete piiridega faktorit kujundada ilma käsitsi sekkumata :(

head(kandidaadid2019$amet)

Kandidaatide kontaktanmded ei moodusta ka faktoriseeritavat tunnust, siis võiks ehk ainult vaadata, kas meiliaadress on või ei ole ja vaadata, kas häältesaak võiks seostuda sellega, kas kandidaat kasutab e-posti või mitte.

head(kandidaadid2019$kontakt)

## [1] "5056540  andres.herkel@gmail.com"                                          
## [2] "56691750"                                                                  
## [3] "5012319  jaak.prozes@fennougria.ee"                                        
## [4] "5093067  jaanika.klopets@gmail.com"                                        
## [5] "5028390  jaanus.ojangu@gmail.com"                                          
## [6] "56354614   juri.pino@gmail.com  \n                      Pardi 3-9, Tallinn"

length(grep("@", kandidaadid2019$kontakt))

## [1] 1051

Teeme muutuja kasutab_meili binaarse faktorina selle põhjal, kas aadressis esineb “@” sümbol:

kandidaadid2019$kasutab_meili <- "ei"; kandidaadid2019$kasutab_meili[grep("@", kandidaadid2019$kontakt)] <- "jah"
kandidaadid2019$kasutab_meili <- factor(kandidaadid2019$kasutab_meili, levels = c("jah","ei"))

Et lisada veel võimalikke kirjeldavaid tunnuseid, tekitame kandidaadi nime põhjal ühe tunnuse, mis võiks kirjeldada nime keerukust: kas nime pikkus (~keerukus) mõjutab häältesaaki :)

head(kandidaadid2019$nimi) 
summary(nchar(kandidaadid2019$nimi))

kandidaadid2019$nime_pikkus <- nchar(kandidaadid2019$nimi)

Ja lõpetuseks: sünniaeg on date-formaadis, mis võib osutuda tülikaks, selle võiks teisendada vanuseks. Valimised toimusid 2019, nii et arvutame vanuse valimiste ajal.

is(kandidaadid2019$sünniaeg)

## [1] "Date"     "oldClass"

head(kandidaadid2019$sünniaeg)

## [1] "1962-08-14" "1984-09-07" "1965-04-04" "1972-12-10" "1966-09-09"
## [6] "1970-07-03"

kandidaadid2019$vanus <- 2019 -as.numeric(substr(kandidaadid2019$sünniaeg, 1,4))

Modelleerimise sammud

Null-mudel ilma ühegi seletava tunnuseta

Kui mudelil ei ole ühtegi seletavat tunnust, uuritava tunnuse varieerumist seletab tema varieerumine ise, siis sellise mudeli valemiks on y ~ 1.

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

## 
## Call:
## lm(formula = log(hääli_kokku) ~ 1, data = kandidaadid2019, 
##     subset = haridus != "Algharidus")
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.0081 -1.1162  0.0026  1.0334  4.8991 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.00805    0.04732   105.8   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.568 on 1097 degrees of freedom

Sõltumatute tunnustena võiks testida:
nimekiri_voi_uksikkandidaat, ringkond, haridus, sugu, vanus

ja nalja pärast ka:
sodiaak, kasutab_meili, nime_pikkus

Aga kuna me eelmisel korral juba vaatasime sugu ja haridust, siis alustame neist.

Esimene samm: lisame ühe faktori

Esimeseks faktoriks lisame hariduse ja võrdleme seda ANOVA testi abil nullmudeliga.

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

## Analysis of Variance Table
## 
## Model 1: log(hääli_kokku) ~ 1
## Model 2: log(hääli_kokku) ~ haridus
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1   1097 2697.3                                  
## 2   1095 2500.1  2    197.28 43.202 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Faktor haridus muudab mudelit oluliselt, nii et jääb sisse. Vaatame mudeli väljundit:

summary(m1)

## 
## 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

anova(m1)

## 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

Teine samm

Lisame teise faktorina soo ning nüüd võrdleme seda mitte nullmudeli vaid eelmise sammuga:

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

## Analysis of Variance Table
## 
## Model 1: log(hääli_kokku) ~ haridus
## Model 2: log(hääli_kokku) ~ haridus + sugu
##   Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
## 1   1095 2500.1                                
## 2   1094 2478.5  1    21.547 9.5107 0.002094 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Jällegi on mõju oluline, sugu peaks ka mudelisse jääma. Vaatame mudeli kokkuvõtet ja ANOVAt ka.

summary(m2)

## 
## Call:
## lm(formula = log(hääli_kokku) ~ haridus + sugu, data = kandidaadid2019, 
##     subset = haridus != "Algharidus")
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.4398 -1.0472  0.0236  0.9863  4.8537 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         4.43984    0.09602  46.237  < 2e-16 ***
## haridusKõrgharidus  0.91185    0.10669   8.547  < 2e-16 ***
## haridusPõhiharidus -1.00316    0.41286  -2.430  0.01527 *  
## sugunaine          -0.29826    0.09671  -3.084  0.00209 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.505 on 1094 degrees of freedom
## Multiple R-squared:  0.08112,    Adjusted R-squared:  0.07861 
## F-statistic:  32.2 on 3 and 1094 DF,  p-value: < 2.2e-16

anova(m2)

## 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.5377 < 2.2e-16 ***
## sugu         1   21.55  21.547  9.5107  0.002094 ** 
## Residuals 1094 2478.52   2.266                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Kolmas samm

Lisame kahe olulise faktori interaktsiooni ja võrdleme seda ainult peamõjusid sisaldava mudeliga:

m3 <- lm(log(hääli_kokku) ~ haridus + sugu + haridus:sugu, data = kandidaadid2019, subset = haridus != "Algharidus")
anova(m2, m3)

## Analysis of Variance Table
## 
## Model 1: log(hääli_kokku) ~ haridus + sugu
## Model 2: log(hääli_kokku) ~ haridus + sugu + haridus:sugu
##   Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
## 1   1094 2478.5                              
## 2   1092 2463.7  2    14.794 3.2787 0.03805 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ka see on oluline! Vaatame jälle ka mudeli väljundi üle.

summary(m3)

## 
## Call:
## lm(formula = log(hääli_kokku) ~ haridus + sugu + haridus:sugu, 
##     data = kandidaadid2019, subset = haridus != "Algharidus")
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5352 -1.0322  0.0279  0.9414  4.8504 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.5352     0.1078  42.054  < 2e-16 ***
## haridusKõrgharidus             0.7737     0.1260   6.138 1.16e-09 ***
## haridusPõhiharidus            -0.6330     0.5122  -1.236  0.21672    
## sugunaine                     -0.6541     0.2083  -3.140  0.00174 ** 
## haridusKõrgharidus:sugunaine   0.4772     0.2355   2.026  0.04300 *  
## haridusPõhiharidus:sugunaine  -0.9474     0.8633  -1.097  0.27269    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.502 on 1092 degrees of freedom
## Multiple R-squared:  0.08661,    Adjusted R-squared:  0.08243 
## F-statistic: 20.71 on 5 and 1092 DF,  p-value: < 2.2e-16

Kuidas seda interaktsiooniga mudelit tõlgendama peaks?

Kuna seletavad tunnused (sugu ja haridus) on nominaalsed, siis mõlemal tunnusel on baastase ja mudeli väljundis toodud efektid eeldavad teiste tunnuste baastaset.
Lõikepunkt ehk intercept on siis olukord, kui haridus = kesk ja sugu = mees, siis on hääli exp(4.5352) = 93.2421619.
Järgmisel real kõrghariduse efekt kehtib siis, kui sugu = mees, seega kõrgharidusega meestel on hääli exp(4.5352 + 0.7737) = 202.1277655.
Edasi, põhihariduse efekt kehtib endiselt siis, kui sugu = mees, põhiharidusega meestel on hääli exp(4.5352 - 0.6330) = 49.5112541
Nüüd järgmiseks tuleb soo efekt, mis kehtib siis, kui haridus on baastasemel, ehk exp(4.5352 - 0.6541) = 48.477511
Ja siis läheb põnevaks: soo ja hariduse vahel on interaktsioon. See tähendab, et kui me tahame teada, palju kõrgharidusega naiskandidaadid hääli said, siis me peame kokku liitma kõik vastavad peaefektid ja interaktsioonid, mis neid puudutavad: baastase + haridus:kõrg + sugu:naine + interaktsioon_kõrg_naine: exp(4.5352 + 0.7737 - 0.6541 + 0.4772) = 169.3554905. Ja põhiharidusega naise hääled siis vastavalt exp(4.5352 - 0.6330 - 0.6541 - 0.9474) = 9.9811668.

library(ggplot2)
df <- ggeffects::ggpredict(m3, terms = c("haridus", "sugu"))

str(df)

## Classes 'ggeffects' and 'data.frame':    6 obs. of  6 variables:
##  $ x        : Factor w/ 3 levels "Keskharidus (sh keskeriharidus)",..: 1 1 2 2 3 3
##  $ predicted: num  93.2 48.5 202.1 169.3 49.5 ...
##  $ std.error: num  0.1078 0.1783 0.0652 0.0884 0.5007 ...
##  $ conf.low : num  75.5 34.2 177.8 142.4 18.5 ...
##  $ conf.high: num  115.2 68.8 229.7 201.4 132.2 ...
##  $ group    : Factor w/ 2 levels "mees","naine": 1 2 1 2 1 2
##  - attr(*, "legend.labels")= chr [1:2] "mees" "naine"
##  - attr(*, "x.is.factor")= chr "1"
##  - attr(*, "averaged_predictions")= logi FALSE
##  - attr(*, "continuous.group")= logi FALSE
##  - attr(*, "rawdata")='data.frame':  1098 obs. of  5 variables:
##   ..$ response: int [1:1098] 222 43 35 23 15 36 9 5 18 12 ...
##   ..$ x       : num [1:1098] 2 2 2 2 2 2 1 2 1 1 ...
##   ..$ group   : Factor w/ 2 levels "mees","naine": 1 2 1 2 1 1 1 1 2 1 ...
##   ..$ facet   : Factor w/ 1 level "1": 1 1 1 1 1 1 1 1 1 1 ...
##   ..$ rowname : chr [1:1098] "1" "2" "3" "4" ...
##  - attr(*, "title")= chr "Predicted values of hääli_kokku"
##  - attr(*, "x.title")= chr "haridus"
##  - attr(*, "y.title")= chr "hääli_kokku"
##  - attr(*, "legend.title")= chr "sugu"
##  - attr(*, "x.axis.labels")= chr [1:3] "Keskharidus (sh keskeriharidus)" "Kõrgharidus" "Põhiharidus"
##  - attr(*, "constant.values")= Named list()
##  - attr(*, "terms")= chr [1:2] "haridus" "sugu"
##  - attr(*, "original.terms")= chr [1:2] "haridus" "sugu"
##  - attr(*, "at.list")=List of 2
##   ..$ haridus: chr [1:3] "Keskharidus (sh keskeriharidus)" "Kõrgharidus" "Põhiharidus"
##   ..$ sugu   : chr [1:2] "mees" "naine"
##  - attr(*, "ci.lvl")= num 0.95
##  - attr(*, "type")= chr "fe"
##  - attr(*, "response.name")= chr "hääli_kokku"
##  - attr(*, "back.transform")= logi TRUE
##  - attr(*, "response.transform")= chr "log(hääli_kokku)"
##  - attr(*, "untransformed.predictions")= num [1:6] 4.54 3.88 5.31 5.13 3.9 ...
##  - attr(*, "family")= chr "gaussian"
##  - attr(*, "link")= chr "identity"
##  - attr(*, "logistic")= chr "0"
##  - attr(*, "is.trial")= chr "0"
##  - attr(*, "link_inverse")=function (eta)  
##  - attr(*, "link_function")=function (mu)  
##  - attr(*, "fitfun")= chr "lm"
##  - attr(*, "verbose")= logi TRUE
##  - attr(*, "model.name")= chr "m3"

ggplot(df, aes(x = x, y = predicted, group = group, color = group)) + 
  geom_point() + 
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.1)

library(sjPlot)
plot_model(m3, type = "int")

Neljas samm

Lisame ühe täiendava faktori, näiteks vanuse:

m4 <- lm(log(hääli_kokku) ~ haridus + sugu + haridus:sugu + vanus, data = kandidaadid2019, subset = haridus != "Algharidus")
anova(m3, m4)

## Analysis of Variance Table
## 
## Model 1: log(hääli_kokku) ~ haridus + sugu + haridus:sugu
## Model 2: log(hääli_kokku) ~ haridus + sugu + haridus:sugu + vanus
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1   1092 2463.7                           
## 2   1091 2463.7  1  0.023058 0.0102 0.9195

Mudel ei läinud paremaks, nii et see jääb välja, aga vaatame enne mudeli väljundit:

summary(m4)

## 
## Call:
## lm(formula = log(hääli_kokku) ~ haridus + sugu + haridus:sugu + 
##     vanus, data = kandidaadid2019, subset = haridus != "Algharidus")
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5350 -1.0300  0.0265  0.9444  4.8484 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   4.5179243  0.2020855  22.356  < 2e-16 ***
## haridusKõrgharidus            0.7718122  0.1274698   6.055 1.93e-09 ***
## haridusPõhiharidus           -0.6297054  0.5134616  -1.226  0.22032    
## sugunaine                    -0.6544883  0.2084666  -3.140  0.00174 ** 
## vanus                         0.0003709  0.0036707   0.101  0.91953    
## haridusKõrgharidus:sugunaine  0.4780132  0.2357734   2.027  0.04286 *  
## haridusPõhiharidus:sugunaine -0.9450990  0.8640236  -1.094  0.27427    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.503 on 1091 degrees of freedom
## Multiple R-squared:  0.08662,    Adjusted R-squared:  0.0816 
## F-statistic: 17.24 on 6 and 1091 DF,  p-value: < 2.2e-16

Selle mudeli p-väärtus on endiselt oluline, st ebaoluline faktor ei tee mudelit otseselt halvemaks, aga kuna ta ei tee seda ka paremaks, siis ei ole seda tarvis sinna panna. Nii et vanus jääb välja, selle asemel proovime midagi muud.

Viies samm

Lisame vanuse asemel näiteks nimekirja, kus kandidaat kandideeris. Kuna sel tulbal on tabelis üsna pikk nimi ja mudeli väljundis esineb see tunnus tulba nimega, siis selle võiks välja vahetada millegi lühema vastu nt “nimek”.

m4 <- lm(log(hääli_kokku) ~ haridus + sugu + haridus:sugu + nimekiri, data = kandidaadid2019, subset = haridus != "Algharidus")
anova(m3, m4)

## Analysis of Variance Table
## 
## Model 1: log(hääli_kokku) ~ haridus + sugu + haridus:sugu
## Model 2: log(hääli_kokku) ~ haridus + sugu + haridus:sugu + nimekiri
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1   1092 2463.7                                  
## 2   1082 1359.8 10      1104 87.847 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

See on jälle väga oluline faktor, jääb sisse!

summary(m4)

## 
## Call:
## lm(formula = log(hääli_kokku) ~ haridus + sugu + haridus:sugu + 
##     nimekiri, data = kandidaadid2019, subset = haridus != "Algharidus")
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.8567 -0.7391 -0.1469  0.6311  3.4952 
## 
## Coefficients:
##                                           Estimate Std. Error t value Pr(>|t|)
## (Intercept)                                5.74918    0.12857  44.718  < 2e-16
## haridusKõrgharidus                         0.50538    0.09623   5.252 1.81e-07
## haridusPõhiharidus                         0.02772    0.38488   0.072   0.9426
## sugunaine                                 -0.06665    0.15743  -0.423   0.6721
## nimekiriEesti Konservatiivne Rahvaerakond -0.25788    0.14364  -1.795   0.0729
## nimekiriEesti Reformierakond               0.16144    0.14181   1.138   0.2552
## nimekiriEesti Vabaerakond                 -2.80746    0.14416 -19.475  < 2e-16
## nimekiriEestimaa Ühendatud Vasakpartei    -2.57966    0.35284  -7.311 5.14e-13
## nimekiriElurikkuse Erakond                -2.14443    0.16593 -12.924  < 2e-16
## nimekiriErakond Eesti 200                 -1.55219    0.14251 -10.891  < 2e-16
## nimekiriErakond Eestimaa Rohelised        -2.31718    0.14474 -16.009  < 2e-16
## nimekiriIsamaa Erakond                    -0.58272    0.14200  -4.104 4.37e-05
## nimekiriSotsiaaldemokraatlik Erakond      -0.75939    0.14204  -5.346 1.10e-07
## nimekiriÜksikkandidaadid                  -1.72573    0.30809  -5.601 2.70e-08
## haridusKõrgharidus:sugunaine               0.08283    0.17662   0.469   0.6392
## haridusPõhiharidus:sugunaine              -0.79832    0.64521  -1.237   0.2162
##                                              
## (Intercept)                               ***
## haridusKõrgharidus                        ***
## haridusPõhiharidus                           
## sugunaine                                    
## nimekiriEesti Konservatiivne Rahvaerakond .  
## nimekiriEesti Reformierakond                 
## nimekiriEesti Vabaerakond                 ***
## nimekiriEestimaa Ühendatud Vasakpartei    ***
## nimekiriElurikkuse Erakond                ***
## nimekiriErakond Eesti 200                 ***
## nimekiriErakond Eestimaa Rohelised        ***
## nimekiriIsamaa Erakond                    ***
## nimekiriSotsiaaldemokraatlik Erakond      ***
## nimekiriÜksikkandidaadid                  ***
## haridusKõrgharidus:sugunaine                 
## haridusPõhiharidus:sugunaine                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.121 on 1082 degrees of freedom
## Multiple R-squared:  0.4959, Adjusted R-squared:  0.4889 
## F-statistic: 70.96 on 15 and 1082 DF,  p-value: < 2.2e-16

Mudeli väljundist näeme, et valimisnimekirja lisamine kasvatas väga palju mudeli R-ruutu, enne oli see kuskil 8% kandis, nüüd on 50%.

Järgmised sammud…

Edasi võiks sammhaaval testida kõik võimalikud kombinatsioonid läbi. Mis on oluline, jääb sisse, mis ei ole, see välja.

Alternatiiv oleks alustada maksimaalsest mudelist, mis oleks:

log(hääli_kokku)~ nimekiri * ringkond * haridus * sugu * vanus * tähtkuju * kasutab_meili * nime_pikkus

See sisaldab 8 faktori peamõju, 28 kahe faktori interaktsiooni, 56 kolme faktori interaktsiooni 70 neljast, 56 viiest, 28 kuuest 8 seitsmest ja 1 kaheksa faktori interakstiooni :O

Tõenäoliselt enamik neist interaktsioonidest ja osa peamõjudest ei ole olulised, aga nende läbitestimine on üsna aeganõudev ning nii keeruka mudeli sisukas seletamine võib olla raske.

ANOVA tüüp III

R-i baaspaketi käsk anova() arvutab ainult esimest tüüpi ANOVA-t, mis põhineb lineaarsel mudelil, kus faktori üks tase on baastasemeks ja kõiki teisi tasemeid võrreldakse baastaseme väärtusega (ja tulemused sõltuvad pisut sellest, millises järjekorras tunnuseid mudelisse lisada). Kui meil aga ei ole üheselt defineeritavat baastaset ja me tahame rühmi võrrelda rühmaülese keskmise suhtes, siis sobib selleks paremini kolmandat tüüpi ANOVA. Paketis car on käsk Anova() (NB! suure algustähega), millel on rohkem valikuid kui baaspaketis.

# kui sul pole see pakett alla laetud, siis kõigepealt installi
# install.packages("car")
library(car)

## Loading required package: carData

Anova(m3, type="III")

## Anova Table (Type III tests)
## 
## Response: log(hääli_kokku)
##              Sum Sq   Df   F value    Pr(>F)    
## (Intercept)  3990.2    1 1768.5709 < 2.2e-16 ***
## haridus        97.8    2   21.6741 5.876e-10 ***
## sugu           22.2    1    9.8577  0.001737 ** 
## haridus:sugu   14.8    2    3.2787  0.038049 *  
## Residuals    2463.7 1092                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Post-hoc test keerukama mudeli puhul

Vaatame Tukey HSD tulemust mudeliga, kus on kahe faktori peamõjud ja nende vaheline interaktsioon:

TukeyHSD(aov(log(hääli_kokku) ~ haridus + sugu + haridus:sugu, data = kandidaadid2019, subset = haridus != "Algharidus"))

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(hääli_kokku) ~ haridus + sugu + haridus:sugu, data = kandidaadid2019, subset = haridus != "Algharidus")
## 
## $haridus
##                                                   diff        lwr         upr
## Kõrgharidus-Keskharidus (sh keskeriharidus)  0.8865119  0.6373791  1.13564458
## Põhiharidus-Keskharidus (sh keskeriharidus) -1.0297743 -1.9964842 -0.06306436
## Põhiharidus-Kõrgharidus                     -1.9162861 -2.8664479 -0.96612435
##                                                 p adj
## Kõrgharidus-Keskharidus (sh keskeriharidus) 0.0000000
## Põhiharidus-Keskharidus (sh keskeriharidus) 0.0335883
## Põhiharidus-Kõrgharidus                     0.0000075
## 
## $sugu
##                  diff        lwr        upr     p adj
## naine-mees -0.2964824 -0.4852892 -0.1076756 0.0021137
## 
## $`haridus:sugu`
##                                                                                   diff
## Kõrgharidus:mees-Keskharidus (sh keskeriharidus):mees                       0.77369590
## Põhiharidus:mees-Keskharidus (sh keskeriharidus):mees                      -0.63304216
## Keskharidus (sh keskeriharidus):naine-Keskharidus (sh keskeriharidus):mees -0.65413385
## Kõrgharidus:naine-Keskharidus (sh keskeriharidus):mees                      0.59673935
## Põhiharidus:naine-Keskharidus (sh keskeriharidus):mees                     -2.23461586
## Põhiharidus:mees-Kõrgharidus:mees                                          -1.40673806
## Keskharidus (sh keskeriharidus):naine-Kõrgharidus:mees                     -1.42782975
## Kõrgharidus:naine-Kõrgharidus:mees                                         -0.17695655
## Põhiharidus:naine-Kõrgharidus:mees                                         -3.00831176
## Keskharidus (sh keskeriharidus):naine-Põhiharidus:mees                     -0.02109169
## Kõrgharidus:naine-Põhiharidus:mees                                          1.22978152
## Põhiharidus:naine-Põhiharidus:mees                                         -1.60157370
## Kõrgharidus:naine-Keskharidus (sh keskeriharidus):naine                     1.25087321
## Põhiharidus:naine-Keskharidus (sh keskeriharidus):naine                    -1.58048201
## Põhiharidus:naine-Kõrgharidus:naine                                        -2.83135521
##                                                                                   lwr
## Kõrgharidus:mees-Keskharidus (sh keskeriharidus):mees                       0.4138731
## Põhiharidus:mees-Keskharidus (sh keskeriharidus):mees                      -2.0951656
## Keskharidus (sh keskeriharidus):naine-Keskharidus (sh keskeriharidus):mees -1.2489066
## Kõrgharidus:naine-Keskharidus (sh keskeriharidus):mees                      0.1987404
## Põhiharidus:naine-Keskharidus (sh keskeriharidus):mees                     -4.1768376
## Põhiharidus:mees-Kõrgharidus:mees                                          -2.8481673
## Keskharidus (sh keskeriharidus):naine-Kõrgharidus:mees                     -1.9697411
## Kõrgharidus:naine-Kõrgharidus:mees                                         -0.4905112
## Põhiharidus:naine-Kõrgharidus:mees                                         -4.9350029
## Keskharidus (sh keskeriharidus):naine-Põhiharidus:mees                     -1.5383262
## Kõrgharidus:naine-Põhiharidus:mees                                         -0.2216485
## Põhiharidus:naine-Põhiharidus:mees                                         -3.9933244
## Kõrgharidus:naine-Keskharidus (sh keskeriharidus):naine                     0.6828955
## Põhiharidus:naine-Keskharidus (sh keskeriharidus):naine                    -3.5645236
## Põhiharidus:naine-Kõrgharidus:naine                                        -4.7655396
##                                                                                    upr
## Kõrgharidus:mees-Keskharidus (sh keskeriharidus):mees                       1.13351866
## Põhiharidus:mees-Keskharidus (sh keskeriharidus):mees                       0.82908128
## Keskharidus (sh keskeriharidus):naine-Keskharidus (sh keskeriharidus):mees -0.05936113
## Kõrgharidus:naine-Keskharidus (sh keskeriharidus):mees                      0.99473830
## Põhiharidus:naine-Keskharidus (sh keskeriharidus):mees                     -0.29239411
## Põhiharidus:mees-Kõrgharidus:mees                                           0.03469120
## Keskharidus (sh keskeriharidus):naine-Kõrgharidus:mees                     -0.88591838
## Kõrgharidus:naine-Kõrgharidus:mees                                          0.13659810
## Põhiharidus:naine-Kõrgharidus:mees                                         -1.08162064
## Keskharidus (sh keskeriharidus):naine-Põhiharidus:mees                      1.49614280
## Kõrgharidus:naine-Põhiharidus:mees                                          2.68121153
## Põhiharidus:naine-Põhiharidus:mees                                          0.79017697
## Kõrgharidus:naine-Keskharidus (sh keskeriharidus):naine                     1.81885092
## Põhiharidus:naine-Keskharidus (sh keskeriharidus):naine                     0.40355954
## Põhiharidus:naine-Kõrgharidus:naine                                        -0.89717078
##                                                                                p adj
## Kõrgharidus:mees-Keskharidus (sh keskeriharidus):mees                      0.0000000
## Põhiharidus:mees-Keskharidus (sh keskeriharidus):mees                      0.8190969
## Keskharidus (sh keskeriharidus):naine-Keskharidus (sh keskeriharidus):mees 0.0214345
## Kõrgharidus:naine-Keskharidus (sh keskeriharidus):mees                     0.0002924
## Põhiharidus:naine-Keskharidus (sh keskeriharidus):mees                     0.0134362
## Põhiharidus:mees-Kõrgharidus:mees                                          0.0604392
## Keskharidus (sh keskeriharidus):naine-Kõrgharidus:mees                     0.0000000
## Kõrgharidus:naine-Kõrgharidus:mees                                         0.5913273
## Põhiharidus:naine-Kõrgharidus:mees                                         0.0001331
## Keskharidus (sh keskeriharidus):naine-Põhiharidus:mees                     1.0000000
## Kõrgharidus:naine-Põhiharidus:mees                                         0.1505630
## Põhiharidus:naine-Põhiharidus:mees                                         0.3954551
## Kõrgharidus:naine-Keskharidus (sh keskeriharidus):naine                    0.0000000
## Põhiharidus:naine-Keskharidus (sh keskeriharidus):naine                    0.2055641
## Põhiharidus:naine-Kõrgharidus:naine                                        0.0004525

Tukey HSD test annab kõikide faktori tasemete vaheliste paaride võrdlused. Keerulisema ja rohketasemelise faktori korral võib see väljund olla üsna pikk ja kirju. Tavaliselt post-hoc testi tulemusi ei raporteerita tabelina vaid tehakse sellest teksti sees üldistusi või korjatakse sealt välja võrdlused, mis on sisuliselt huvitavad.

Kordamisküsimused

Kordamisküsimuste testi saad ka moodlis teha, seal näed pärast ka õigeid vastuseid ja kommentaare.

1. Kui mudelis on oluline kahe tunnuse interaktsioon, siis

peamõjud enam ei loe
peamõjude efekt muutub vastupidiseks
peamõjude vahel on täiendav seos
peamõjud toimivad ühte moodi

2. Ockhami habemenoa printsiipi järgides peaks valima mudeli, millel on

kõrge R-ruut ja väike p-väärtus
madal kohandatud R-ruut ja palju seletavaid tunnuseid
võimalikult vähe seletavaid tunnuseid
võimalikult kõrge AIC

3. Parem mudel on see, millel on

kõrge kohandatud R-ruut ja madal AIC
madal kohandatud R-ruut ja madal AIC
kõrge kohandatud R-ruut ja kõrge AIC
madal kohandatud R-ruut ja kõrge AIC

4. Lineaarses mudelis peaks vältima multikolineaarsust. See tähendab, et:

seletavad tunnused ei tohi omavahel korreleeruda
seletavad tunnused peavad omavahel korreleeruma
uuritava ja seletava tunnuse vahel peab olema lineaarne seos
seletavate tunnusete vahel peab olema interaktsioon

5. Kolmandat tüüpi ANOVA (Type III) tähendab, et

võrreldakse kolme rühma
rühmasid võrreldakse valimi keskmisega
testitakse kahte peamõju ja interaktsiooni
rühmasid võrreldakse baastasemega

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

##                       Estimate Std. Error    t value     Pr(>|t|)
## (Intercept)        771.2282027  94.279541  8.1802286 1.422446e-12
## Length              17.8944667   9.834225  1.8196112 7.203563e-02
## log1p(Freq)        -22.6968893  12.557574 -1.8074263 7.392896e-02
## Length:log1p(Freq)   0.1574769   1.393383  0.1130177 9.102602e-01

7. Mis võiks olla reaktsiooniaeg, kui sõnas on 26 tähemärki ja sõna on keskmise sagedusega?

771.23 + 26 * 17.89 + -22.70 + 26 * 0.16 = 1217.83 millisekundit
771.23 + 26 * 17.89 + 5.68 * -22.70 + 26 * 5.68 * 0.16 = 1131.0628 millisekundit
771.23 + 26 * 17.89 = 1236.37 millisekundit
keskmise sagedusega sõnad ei ole kunagi nii pikad

Järgmisel korral

Lineaarse regressiooni segamudelid

Funktsioonide sõnastik

lm() – lineaarne mudel
anova() – anova test (vaikimisi I tüüpi) -AIC() – Aikaike informatsioonikriteerium
step() – stepwise protseduur parima mudeli valimiseks
car::Anova(lm(), type="III") – kolmandat tüüpi anova
cov2cor(cov()) – tunnuste vahelise konrrelatsiooni testimine
car::vif(model1) – VIF ehk Variance Inflation Factor

12. praktikum

03.04.2024

Lineaarse regressiooni mudelid 2

Tänases praktikumis

Kaks arvulist tunnust

Näide 1: Kas lisaks sõna pikkusele võiks reaktsiooniaega mõjutada sõna sagedus?

Kuidas mudelit lugeda?

Väike lõbus vahepala: interaktiivsed joonised

Peamõju ja interaktsioon

Kaks arvulist tunnust interaktsiooniga

Näide 2: Kas sõna pikkus ja sagedus mõjutavad reaktsiooniaega?

Kuidas interaktsioone tõlgendada?

Mudelite võrdlemine

Occami habemenoa printsiip

Võrdleme mudeleid: R-ruut, ANOVA ja AIC

ANOVA test mudeli jääkidele

Akaike informatsioonikriteerium

Inkrementaalne modelleerimine

Näide 3: täiskasvanud kõnelejate kõnetempot mõjutavad tegurid

Multikollineaarsus

Stepwise protseduur

Näide 4: Mis mõjutab valimistulemust?

Ettevalmistused suuremaks mudeliks

Multikollineaarsus

Vaatame enne mudeldamist üle ülejäänud tunnused

Modelleerimise sammud

Null-mudel ilma ühegi seletava tunnuseta

Esimene samm: lisame ühe faktori

Teine samm

Kolmas samm

Neljas samm

Viies samm

Järgmised sammud…

ANOVA tüüp III

Post-hoc test keerukama mudeli puhul

Kordamisküsimused

1. Kui mudelis on oluline kahe tunnuse interaktsioon, siis

2. Ockhami habemenoa printsiipi järgides peaks valima mudeli, millel on

3. Parem mudel on see, millel on

4. Lineaarses mudelis peaks vältima multikolineaarsust. See tähendab, et:

5. Kolmandat tüüpi ANOVA (Type III) tähendab, et

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

7. Mis võiks olla reaktsiooniaeg, kui sõnas on 26 tähemärki ja sõna on keskmise sagedusega?

Järgmisel korral

Funktsioonide sõnastik