Lineaarse regressiooni segamudelid

Tänases praktikumis

  • Segamudelid
    • juhusliku vabaliikmega
    • juhusliku kaldega
    • ANOVA ja post-hoc testid segamudelitega

Andmestikud:

  • fonkorp2 (loeme otse DataDOIst, ei pea alla laadima)
  • välteandmestik (loeme otse DataDOI-st)

Segamudelid

Segamudelid: fikseeritud ja juhuslikud faktorid

Kui seni oleme tegelenud mudelitega, mis eeldavad sõltumatuid mõõtmisi, siis segamudel arvestab sellega, et osa mõõtmisi võib tulla korduvatelt subjektidelt.

  • Fikseeritud faktoril on teada kindel hulk tasemeid, mille mõju kaudu me otseselt seletame uuritava tunnuse varieerumist.
  • Juhuslikul faktoril võib olla määramata hulk tasemeid ja meie valimis on nendest juhuslik valik, millel võib ka olla mõju uuritava tunnuse varieerumisele ja millega tuleb seetõttu arvestada, ent mis ei seleta ise otseselt uuritava tunnuse varieerumist ega paku uurijale omaette huvi.
  • Juhuslik faktor on enamasti “müra”, mis tuleks välja filtreerida. Juhuslikud faktorid võimaldavad korrigeerida regressioonimudeli fikseeritud mõjude koefitsiente, arvestades sellega, et iga üksik subjekt võib olla järjekindlalt teistest pisut erinev.

Juhuslik vabaliige ja kalle

Juhuslik vabaliige (random intercept) tähendab seda, et lisaks keskmisele vabaliikmele korrigeeritakse seda iga juhusliku faktori taseme jaoks (nt iga indiviidi jaoks algaks koefitsientide võrdlemine natuke erinevalt baastasemelt).

Juhuslik kalle (random slope) tähendaks, et lisaks fikseeritud faktori efektile korrigeeritakse seda efekti vastavalt iga juhusliku faktori tasemele (nt iga indiviidi kohta).

Näide 1: Millest sõltub kõneleja aktiivsus vestluses?

See näide lähtub artiklis Aare, K., & Lippus, P. (2017). Some gender patterns in Estonian dyadic conversations. J. E. Abrahamsen, J. Koreman, & W. van Dommelen (Toim), Nordic Prosody. Proceedings of the XIIth Conference, Trondheim 2016 (lk 29–38). Peter Lang. https://doi.org/10.13140/RG.2.2.17105.94568

Selles uurimuses on vaadatud mitmeid üldiseid hääle kõrgust ja häälelaadi kirjeldavaid parameetreid seoses kõneleja soo ja vanuse ning vestluspartneri vastavate tunnustega. Korpus, mille põhjal on uurimus tehtud, sisaldab kahe kõneleja vaheliste vestluste salvestusi. Iga vestlus kestab umbes pool tundi ja iga vestluse iga osaleja kohta on salvestuse põhjal mõõdetud tema hääle põhitooni keskmist kõrgust (hertsides), põhitooni ulatust (pooltoonides), kui palju ta selle vestluse jooksul rääkis (% salvestuse ajast) ja kui palju tema hääl sel ajal kärises (% kõneldud ajast). Lisaks häälelaadile vaadati selles artiklis ka seda, kui palju kõneleja üldse vestluse jooksul rääkis. Leiti, et kõnes osalemine sõltus kõneleja soost (mehed räägivad rohkem), vestluspartneri soost (meestega räägitakse vähem) ja kõnelejate vanusevahest (endast vanematega räägitaske vähem).

“Aare Lippus 2017”
“Aare Lippus 2017”

Täiskasvanud kõnelejate kõnetempo andmestik, mida kasutasime eelmises praktikumis, on natuke uuem variant sellest samast andmestikust.

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

Kui 2016. aastal olid andmestikus andmed 88 eri kõnelejalt, siis nüüd on neid 139. Proovime selle mudeli uuenenud andmete peal läbi mängida, loodetavasti ei ole tulemused kardinaalselt muutunud.

Arvutame proportsionaalse kõneaja (kõneleja kõnevoorude kestus jagatud vestluse kogukestusega):

fonkorp2$kestus_prop <- fonkorp2$kestus/fonkorp2$faili_kogukestus
  1. kõneleja sugu (meie andmestikus 2 taset),
  2. kõneleja vanus (pidev tunnus),
  3. vestluspartneri sugu (2 taset),
  4. vestluspartneri vanus (pidev tunnus),
  5. kõnelejate vanusevahe (pidev tunnus).

Kuna Aare ja Lippus (2017) oli soo baastasemeks naine ja fonkorp2 andmestikus on mees, siis vahetame parema võrreldavuse eesmärgil siin ka tasemete järjekorra ära:

fonkorp2$sugu <- factor(fonkorp2$sugu, levels = c("N", "M"))
fonkorp2$kaaskoneleja_sugu <- factor(fonkorp2$kaaskoneleja_sugu, levels = c("N", "M"))

Kas lineaarse mudeli eeldused on täidetud?

  • Uuritav tunnus on üllataval kombel isegi Shapiro testi põhjal normaaljaotusega. Kirjeldavad arvulised tunnused (vanus ja vanusevahe) ei ole normaaljaotusega, siin ei aita ka logaritmimine. Aga õigupoolest huvitab meid rohkem mudeli jääkide jaotus ja iga üksik tunnus ei peagi normaaljaotusega olema. Tihtipeale võib aga juhtuda, et kui arvulised tunnused on väga ühele või teisele poole kaldu, ei ole ka mudeli jääkide jaotus normaalne. Mudeli jääkide jaotust saab aga kontrollida pärast mudeli tegemist.
  • Seletavate tunnuste vahel ei tohiks olla kolineaarsust. Eelmisest korrast mäletame, et kummagi kõneleja vanused olid vanusevahega tugevas korrelatsioonis. Seega kasutame see kord (nagu Aare ja Lippus 2017) vanusevahet.
  • Uuritava tunnuse ja arvulise seletava tunnuse vahel peaks olema lineaarne seos.
library(tidyverse)
library(plotly)
joonis1 <- ggplot(data = fonkorp2, aes(x = vanusevahe, y = kestus_prop, col = kaaskoneleja_sugu)) + geom_point() + geom_smooth(method = "lm") + facet_wrap("sugu")
ggplotly(joonis1)
  • Iga vestluse iga osaleja kohta on andmestikus üks mõõtmispunkt, mis justkui tähendaks, et mõõtmised on sõltumatud. Üldiselt on igas vestluses osalejad erinevad, aga on üksikuid koruvaid osalejaid. Kuna neid ei ole väga palju, siis me võiksime kasutada tavalist lineaarset mudelit ja kui tahame väga hoolega järgida sõltumatute mõõtmiste nõuet, siis võiksime valida korduvate osalejate salvestustest ühe ja ülejäänud kordused andmestikust eemaldada. Vaatame üle, palju on korduvaid osalejaid:
table(table(fonkorp2$koneleja))
## 
##   0   1   2   3   4   5 
##  11 125   8   2   3   1

125 osalejat esineb andmestikus ainult ühe korra, aga 14 osalejat osales mitu korda, üks neist lausa viis korda. (11 null-osalejat on sellest, et välja on jäätud outlierid, näiteks üle 70-aastased osalejad, keda oli väga vähe).

Miks see korduvate mõõtmispunktide küsimus on üldse oluline? Selle pärast, et lisaks fikseeritud faktoritele (sugu, vanus) võib uuritavat tunnust mõjutada (õigemini: mõjutab) tema individuaalne eripära. Kuna meil on andmestikus mõned korduvad kõnelejad, siis võib nende individuaalne efekt kallutada fikseeritud faktorite mõju.

Kõneleja isikut me ei saa fikseeritud faktorina mudelisse lisada, sest meil on valimis juhuslikud 139 indiviidi (8 miljardist inimesest või 1 miljonist eestlasest). Pealegi me nii väga ei tahagi konkreetse indiviidi efekti määratleda, me lihtsalt tahame selle olemasoluga arvestada.

Testime kõigepealt tavalise lineaarse mudeliga

lm0 <- lm(kestus_prop ~ 1, data = fonkorp2)
# lisame 1. tunnusena kõneleja soo
lm1 <- lm(kestus_prop ~ sugu, data = fonkorp2)
anova(lm0, lm1)
# see oli oluline, lisame teise tunnusena vestluspartneri soo
lm2 <- lm(kestus_prop ~ sugu + kaaskoneleja_sugu, data = fonkorp2)
anova(lm1, lm2)
# see oli ka oluline, nüüd lisame soo ja vestluspartneri soo interaktsiooni
lm3 <- lm(kestus_prop ~ sugu * kaaskoneleja_sugu, data = fonkorp2)
anova(lm2, lm3)
# see pole oluline, proovime selle asemel vanusevahe peaefekti
lm4 <-lm(kestus_prop ~ sugu + kaaskoneleja_sugu + vanusevahe, data = fonkorp2)
anova(lm2, lm4) # võrdleme 2 peaefektiga mudelit kolme peafektiga mudeliga
# see on oluline; võiks edasi proovida interaktsioone
lm5a <-lm(kestus_prop ~ sugu + kaaskoneleja_sugu + vanusevahe + sugu:vanusevahe, data = fonkorp2)
lm5b <-lm(kestus_prop ~ sugu + kaaskoneleja_sugu + vanusevahe + kaaskoneleja_sugu:vanusevahe, data = fonkorp2)
anova(lm4, lm5a)
anova(lm4, lm5b)
# pole oluline, parim mudel on lm4 kolme faktori peamõjuga ja ilma interaktsioonideta

Parim mudel on niisiis see:

lm4 <-lm(kestus_prop ~ sugu + kaaskoneleja_sugu + vanusevahe, data = fonkorp2)
summary(lm4)
## 
## Call:
## lm(formula = kestus_prop ~ sugu + kaaskoneleja_sugu + vanusevahe, 
##     data = fonkorp2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.33730 -0.13112  0.00642  0.12387  0.41342 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         0.475428   0.020862  22.789  < 2e-16 ***
## suguM               0.070865   0.026830   2.641  0.00908 ** 
## kaaskoneleja_suguM -0.083173   0.026830  -3.100  0.00229 ** 
## vanusevahe         -0.004382   0.000788  -5.560 1.11e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1698 on 160 degrees of freedom
## Multiple R-squared:  0.2307, Adjusted R-squared:  0.2163 
## F-statistic:    16 on 3 and 160 DF,  p-value: 3.824e-09
ggeffects::ggpredict(lm4)
sjPlot::plot_model(lm4, type = "pred", show.data = T, terms = c("vanusevahe", "kaaskoneleja_sugu","sugu"))

Kontrollime jääkide normaalsust.

hist(lm4$residuals)

shapiro.test(lm4$residuals) # eeldus on täidetud
## 
##  Shapiro-Wilk normality test
## 
## data:  lm4$residuals
## W = 0.98552, p-value = 0.08653

Näide 1a: Sama lineaarse segamudeliga (lme4 pakett)

Installime ja aktiveerime paketi lme4, kus tuleb segamudeli käsk lmer(). Lisaks oleks hea installida ka paket lmerTest, mis lisab mudeli väljundile p-väärtused (vaikimis on segamudeli väljund ilma p-väärtusteta).

# install.packages("lme4", "lmerTest")
library(lme4)
library(lmerTest)

Segamudelis peab olema vähemalt üks juhuslik faktor. Lisame juhusliku kõneleja efekti vabaliikmele (random intercept). lmer() mudelis tähistatakse juhuslikku vabaliiget (1|juhuslik_faktor). Muus osas on mudeli sisend täpselt sama, nagu tavalise lineaarse mudeli puhul.

sm4 <- lmer(kestus_prop ~ sugu + kaaskoneleja_sugu + vanusevahe + (1|koneleja), data = fonkorp2)

Kas mudel muutus juhusliku faktori lisamisest?

Võrdleme ainult fikseeritud faktoritega mudelit ja lisatud juhusliku mõjuga mudelit:

anova(sm4, lm4)
## refitting model(s) with ML (instead of REML)
## Data: fonkorp2
## Models:
## lm4: kestus_prop ~ sugu + kaaskoneleja_sugu + vanusevahe
## sm4: kestus_prop ~ sugu + kaaskoneleja_sugu + vanusevahe + (1 | koneleja)
##     npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)  
## lm4    5 -110.30 -94.804 60.152  -120.30                       
## sm4    6 -113.41 -94.810 62.705  -125.41 5.1063  1    0.02384 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Segamudel on oluliselt erinev. Vaatame nüüd mudeli väljundit käsuga summary():

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: kestus_prop ~ sugu + kaaskoneleja_sugu + vanusevahe + (1 | koneleja)
##    Data: fonkorp2
## 
## REML criterion at convergence: -95.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.06916 -0.67141  0.00286  0.59425  1.99382 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  koneleja (Intercept) 0.008335 0.0913  
##  Residual             0.020119 0.1418  
## Number of obs: 164, groups:  koneleja, 139
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)         4.859e-01  2.139e-02  1.480e+02  22.718  < 2e-16 ***
## suguM               6.319e-02  2.784e-02  1.414e+02   2.269  0.02476 *  
## kaaskoneleja_suguM -8.131e-02  2.632e-02  1.586e+02  -3.090  0.00237 ** 
## vanusevahe         -4.585e-03  7.927e-04  1.586e+02  -5.784 3.78e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) suguM  kskn_M
## suguM       -0.541              
## kasknlj_sgM -0.470 -0.122       
## vanusevahe   0.003  0.039 -0.052

Kuna siin mudelis on seletavate tunnustena kõneleja enda omadused, mis ühe kõneleja puhul korduda ei saa, siis siia juhuslikku kallet lisada ei saa.

Mudeli väljund on üldiselt väga sarnane tavalise lineaarse mudeli väljundile. Mudeli üldist olulisust ei ole küll eraldi real välja toodud, aga on mudeli jääkide jaotus (residuals). Et nullhüpotees hüljata, peaks mudeli jäägid läbima nulli.

Väljundi lõpus on fikseeritud efektide korrelatsioonimaatriks, mille jälgimine võiks aidata multikollineaarsust vältida.

Kõige huvitavam osa on segamudeli väljundis fikseeritud efektide koefitsientide tabel. Selle lugemine ja tõlgendamine on põhimõtteliselt sama nagu tavalise lineaarse mudeli puhul. Ilma lmerTest paketi kasutamata ei esitata fikseeritud efektide tabelis p-väärtusi. Sel juhul peaks hindama faktori olulisust statistiku t väärtuse alusel: rusikareegel on see, et kui selle absoluutväärtus on suurem kui kaks |t| > 2, siis on faktor oluline.

Kui varasemates näidetes meil on sageli olnud logaritmitud tunnuseid, siis siin on tõlgendamine seda lihtsam, et kõik tunnused on algsetes ühikutes. Võib-olla lihtsam on mudeli efekte hinnata siis, kui teisendame proportsionaalse kestuse protsentideks:

 summary(sm4)$coefficients[,"Estimate"] *100
##        (Intercept)              suguM kaaskoneleja_suguM         vanusevahe 
##         48.5908982          6.3190324         -8.1313863         -0.4585191
  • vabaliige (intercept) on siis keskmine põhitoon siis, kui kõigi seletavade tunnuste väärtus on 0 või nad on baastasemel. See tähendab siis, et kui vestlesid kaks samavanust naist, rääkis kumbki 48% ajast.
  • Kõneleja soo efekt tähendab seda, et kui mees vestles naisega, siis mees rääkis 48+6 = 54% ajast.
  • Vestluspartneri soo efekt tähendab seda et kui naine vestles mehega, siis naine rääkis 48-8 = 40% ajast.
  • Vanusevahe on andmestikus arvutatud nii, et positiivne vanusevahe tähendab kõnelejast vanemat, negatiivne nooremat vestluspartnerit. Seega endast vanemaga rääkides väheneb kõnelemise aeg 0.45% iga vanusevahe aasta kohta, nii et kui naine räägib endast 30 aastat vanema naisega, räägib ta 48 + 30*-0.45 = 35% ajast, aga kui näiteks mees räägib endast 30 aastat noorema naisega, siis ta räägib 48 + 6 - 30*-0.45 = 68% ajast.
sjPlot::plot_model(sm4, type = "pred", show.data = T, terms = c("vanusevahe", "kaaskoneleja_sugu","sugu"))

Testime ka mudeli jääke

Lineaarsele segamudelile on endiselt eeldus, et mudeli jäägid on normaaljaotusega:

par(mfcol=c(1,2))
hist(residuals(sm4))
qqnorm(residuals(sm4))
qqline(residuals(sm4))

par(mfcol=c(1,1))
shapiro.test(residuals(sm4))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(sm4)
## W = 0.98909, p-value = 0.2363

Tundub, et on ilus normaaljaotus.

Vaatame juhuslikke lõikepunkte

Segamudelis on juhuslik faktor, siin mudelis juhuslik lõikepunkt iga kõneleja kohta. Kui meie mudel hindab seda, kui palju kõneleja vestluses osaleb, siis üldine mudeli fikseeritud lõikepunkt:

 summary(sm4)$coefficients["(Intercept)","Estimate"]
## [1] 0.485909

tähendab, et keskmiselt osalesid kõnelejad vestluses 48% ajast, aga juhuslik lõikepunkt kõneleja kohta tähendab seda, et selle võrra erines iga üksik indiviid keskmisest (võttes arvesse tema korduvaid osalemisi erinevate vestluspartneritega). Juhuslikke faktoreid saab mudeli objektist käsuga ranef().

head(ranef(sm4)$koneleja)
##      (Intercept)
## KJ3 -0.004377745
## KJ4  0.091351304
## KJ1 -0.082068129
## KJ5 -0.091768838
## KJ2 -0.041298905
## KJ6  0.032933419

Siit näeme siis seda, et KJ1 oli keskmisest 8% vähem jutukas, KJ2 4% ja KJ3 0.4% vähem jutukad, aga KJ4 rääkis keskmisest 9% rohkem.

Näide 1b: Sama GAM-mudeliga (mgcv pakett)

lme4 on üks populaarsemaid R-i pakette lineaarsete segamudelite jaoks. Teine väga populaarne pakett, millega saab nii lineaarseid segamudeleid kui ka mittelineaarsete seostega üldistatud aditiivseid segamudeleid (genralised aditive mixed models, GAMM), on pakett mgcv. Proovime korra sama mudelit teise paketiga:

#install.packages("mgcv")
library(mgcv)

Käsu gam() süntaks on üldiselt sama nagu eelmistel mudelitel, aga kui lmer() juhuslik vabaliige tähistatakse (1|tunnus), siis sama asi gam-is on s(tunnus, bs="re").

gam4 <- gam(kestus_prop ~ sugu + kaaskoneleja_sugu + vanusevahe + s(koneleja, bs = "re"), data = fonkorp2)
summary(gam4)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## kestus_prop ~ sugu + kaaskoneleja_sugu + vanusevahe + s(koneleja, 
##     bs = "re")
## 
## Parametric coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         0.4844236  0.0211617  22.891  < 2e-16 ***
## suguM               0.0643661  0.0274997   2.341   0.0208 *  
## kaaskoneleja_suguM -0.0813726  0.0262607  -3.099   0.0024 ** 
## vanusevahe         -0.0045566  0.0007876  -5.785 5.48e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##               edf Ref.df     F p-value   
## s(koneleja) 35.12    137 0.397 0.00756 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.419   Deviance explained = 55.5%
## GCV = 0.02806  Scale est. = 0.021366  n = 164

GAM-mudelisse on võimalik lisada smooth-efektidena mittelineaarseid seoseid. Praegu siin mudelis on fikseeritud efektid lineaarsed, aga juhuslik kõneleja efekt on smooth efektina. Proovime korra ka seda, kas vanusevahe efekt võiks olla mittelineaarne:

gam5 <- gam(kestus_prop ~ sugu + kaaskoneleja_sugu + s(vanusevahe) + s(koneleja, bs = "re"), data = fonkorp2)
summary(gam5)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## kestus_prop ~ sugu + kaaskoneleja_sugu + s(vanusevahe) + s(koneleja, 
##     bs = "re")
## 
## Parametric coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         0.48442    0.02116  22.891   <2e-16 ***
## suguM               0.06437    0.02750   2.341   0.0208 *  
## kaaskoneleja_suguM -0.08137    0.02626  -3.099   0.0024 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                 edf Ref.df      F p-value    
## s(vanusevahe)  1.00      1 33.472 < 2e-16 ***
## s(koneleja)   35.12    137  0.397 0.00756 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.419   Deviance explained = 55.5%
## GCV = 0.02806  Scale est. = 0.021366  n = 164

Mittelineaarsete efektide puhul peaks vaatama edf kordajat: kui see on suurem kui 2, siis on seoses mittelineaarsust, siin mudlis vanusevahe edf = 1, seega tegemist on siiski lineaarse seosega.

Näide 2: vokaali kestus ja välde

Kasutame teist andmestikku, mis on kogutud samast 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. Andmestikus on info 1714 kahesilbilise sõna kohta, mõõdetud on häälikute kestust, põhitooni liikumist, vokaalikvaliteeti jms. Samuti on märgitud sõnade välde, fookusrõhk jms.

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

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

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 rohkem kui üks rida andmestikus; samuti võivad andmestikus esineda korduvad sõnad.

# ühelt kõnelejalt on vähemalt 23 ja kuni 203 sõna
quantile(table(dat2$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(dat2$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

Eesti keeles on 3 väldet (nt sada : saada : saada) ja välde kõige rohkem mõjutab rõhulise vokaali kestust. Kui uurida vokaali kestust, siis fikseeritud faktorina võiks kestuse varieerumist seletada:

  1. välde,
  2. lauserõhk (fookus)
  3. sõnaliik.

Vokaalide kestust võib mõjutada ka kõnetempo, mis võib olla kõnelejati erinev. Samuti võib mõjutada vokaali häälikuline kontekst ehk sõna, milles ta esineb. Nii nagu kõnelejatest on meil juhuslik valim, on meil ka keeles esinevatest sõnadest juhuslik valim, mistõttu me ei saa nende fikseeritud efekti arvestada.

Kodeerime mõned numbritena sisestatud faktorid ümber faktoriteks:

dat2$Quantity <- factor(dat2$Quantity)
dat2$Focus <- factor(dat2$Focus); levels(dat2$Focus) <- c("no","yes")

Kas rõhulise silbi vokaali kestus sõltub vältest?

Jätame välja kinnise esisilbiga (nn metsa-tüüpi) sõnad, sest neis on esimene vokaal fonoloogiliselt lühike, sest väldet kannab silbilõpukonsonant.

boxplot(V1 ~ Quantity, data = dat2, subset = Struct == "cvcv", 
        ylab = "V1 kestus (ms)", xlab = "Välde")

Kui me tahaks tavalist anovat või lineaarset mudelit ilma juhuslike efektideta

Kui me tahaks rakendada selle andmestiku peal varasematest praktikumidest tuttavaid meetodeid ehk ilma juhuslike efektideta anovat või lineaarset mudelit, siis selle kuidagi saama kooskõlla eeldusega, et tegemist on sõltumatute mõõtmistega. Üks võimalus selleks on keskmistada nii, et näiteks iga kõneleja kohta oleks iga välte kohta üks mõõtmine. (Umbes nii on ka Lippus jt 2013 artiklis analüüs läbi viidud.)

dat2 %>% 
    filter(Struct == "cvcv") %>% 
    group_by(SP, Quantity) %>% 
    summarise(V1 = mean(V1)) %>% 
    lm(V1 ~ Quantity, data = .) %>% 
    summary()
## 
## Call:
## lm(formula = V1 ~ Quantity, data = .)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -38.466 -11.619  -3.632   5.832  77.201 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   58.366      4.865  11.996  < 2e-16 ***
## Quantity2     53.005      6.881   7.703 1.57e-10 ***
## Quantity3     84.433      6.881  12.271  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22.3 on 60 degrees of freedom
## Multiple R-squared:  0.7194, Adjusted R-squared:  0.7101 
## F-statistic: 76.93 on 2 and 60 DF,  p-value: < 2.2e-16

Juhuslik kõneleja vabaliige

Laeme lineaarsete segamudelite paketi lme4 ja lisapaketi lmerTest, mis lisab lineaarsele segamudelile p-väärtused.

# install.packages("lme4", "lmerTest")
library(lme4)
library(lmerTest)

Kui teeme 0-mudeli, siis lmer’i eeldus on, et mudelis on vähemalt üks juhuslik efekt. Kui ühtegi fikseeritud efekti ei ole, siis saab olla selleks ainult juhuslik vabaliige.

Juhuslikku vabaliiget tähistab (1|juhuslik_faktor)

Lisame esimeseks vabaliikmeks kõneleja, seda võiks tõlgendada kui individuaalset kõnetempo efekti.

m0 <- lmer(V1 ~ (1|SP), data = dat2, subset = Struct == "cvcv", REML = F)

Ja siis lisame ka ühe fikseeritud faktori, selleks võiks olla välde:

m1 <- lmer(V1 ~ Quantity + (1|SP), data = dat2, subset = Struct == "cvcv", REML = F)
anova(m0, m1)
## Data: dat2
## Subset: Struct == "cvcv"
## Models:
## m0: V1 ~ (1 | SP)
## m1: V1 ~ Quantity + (1 | SP)
##    npar   AIC   BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## m0    3 13826 13841 -6909.8    13820                         
## m1    5 12637 12663 -6313.5    12627 1192.8  2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m1)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: V1 ~ Quantity + (1 | SP)
##    Data: dat2
##  Subset: Struct == "cvcv"
## 
##      AIC      BIC   logLik deviance df.resid 
##  12636.9  12662.9  -6313.5  12626.9     1340 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1132 -0.5804 -0.0664  0.4995  6.4701 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  SP       (Intercept) 154.0    12.41   
##  Residual             672.3    25.93   
## Number of obs: 1345, groups:  SP, 21
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   57.534      2.916   21.688   19.73 2.43e-15 ***
## Quantity2     60.672      2.107 1344.435   28.79  < 2e-16 ***
## Quantity3     81.405      2.117 1336.609   38.46  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) Qntty2
## Quantity2 -0.159       
## Quantity3 -0.146  0.216

Välte efekt oli oluline. Esimeses vältes oli vokaali pikkus 58 millisekundit, teises 118 ms ja kolmandas vältes 139 ms.

Vaatame nüüd kõnelejatepõhiseid vabaliikmeid:

ranef(m1)
## $SP
##      (Intercept)
## SP1  -18.8245525
## SP10  22.8715306
## SP11  -3.4002725
## SP12  -4.2580884
## SP13  10.6962704
## SP14  14.9482950
## SP15   1.0772266
## SP16  31.0334271
## SP17  -5.7978398
## SP18  -2.7243061
## SP19  -6.2832400
## SP2  -12.4845180
## SP20  -0.3906256
## SP21  -2.9826901
## SP3  -10.8790722
## SP4   -7.1449626
## SP5   -7.5051669
## SP6    0.1493260
## SP7   -7.2053662
## SP8   -3.4094710
## SP9   12.5140962
## 
## with conditional variances for "SP"

SP1 kõnetempo on keskmiselt 19 ms kiirem (Q1 58-19=39, Q2 118-19=99, Q3 139-19=120), SP10 on jällegi 23 ms aeglasem, jne. Individuaalne vabaliige tähendab seda, et see efekt liitub kõigele üldiselt.

sjPlot::plot_model(m1, type = "re", sort.est = T, grid = F)

Siin joonisel on kast-vurrud diagrammina vokaalikestuse jaotused kolmes vältes ja punase joonega on lisatud mudeli hinnang välte efektile:

# baasgraafikaga
plot(V1 ~ Quantity, data = dat2, subset = Struct=="cvcv")
coefs <- summary(m1)$coefficients[,"Estimate"]
lines(y = coefs["(Intercept)"] + c(0, coefs[c("Quantity2", "Quantity3")]), 
      x = 1:3, lwd = 3, col = "red")
#ggplotiga sama
coefs <- summary(m1)$coefficients[,"Estimate"]
ggplot() +
    geom_boxplot(data = dat2 %>% filter(Struct=="cvcv"), aes(x = Quantity, y = V1)) +
    geom_line(aes(x = 1:3,
                  y = c(coefs["(Intercept)"],
                        coefs["(Intercept)"] + coefs["Quantity2"],
                        coefs["(Intercept)"] + coefs["Quantity3"]),
                  group = 1),
              color = "red", size = 1)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Siin joonisel on eri värviga iga kõneleja individuaalne joon, mida on mudeli juhusliku vabaliikme põhjal korrigeeritud:

# baasgraafikaga
plot(V1 ~ Quantity, data = dat2, subset = Struct == "cvcv")
coefs <- summary(m1)$coefficients[,"Estimate"]
for(i in 1:nrow(ranef(m1)$SP)) lines(y = ranef(m1)$SP[i,] + coefs["(Intercept)"] + c(0,coefs[c("Quantity2", "Quantity3")]), x = 1:3, col = i)
# ggplotiga sama
coefs <- summary(m1)$coefficients[,"Estimate"]

ranef(m1)$SP %>%
    mutate(koneleja = rownames(.),
           Q1 = coefs["(Intercept)"] + `(Intercept)`,
           Q2 = Q1 + coefs["Quantity2"],
           Q3 = Q1 + coefs["Quantity3"]) %>%
    pivot_longer(., cols = starts_with("Q"), names_prefix = "Q",
                 names_to = "välde", values_to = "kõnetempo") -> df

ggplot() +
    geom_boxplot(data = dat2 %>% filter(Struct == "cvcv"), 
                 aes(x = Quantity, y = V1)) +
    geom_line(data = df,
              aes(x = välde, y = kõnetempo, color = koneleja, group = koneleja),
              show.legend = F)

Juhuslik kalle

Lisaks individuaalsele üldisele kõnetempole on välte mõju ka kõnelejati veidi erinev. See tähendab, et aeglasema kõnetempoga võib-olla aeglustatakse pikemaid häälikuid rohkem, lühemaid vähem, või vastupidi.

Juhuslikku kallet (random slope) tähistab (1+fikseeritud_faktor|juhuslik_faktor). Juhuslik kalle on seotud mingi fikseeritud faktoriga, see see, mille võrra fikseeritud efekti mõju korrigeeritakse.

m2 <- lmer(V1 ~ Quantity + (1+Quantity|SP), data = dat2, subset = Struct=="cvcv", REML = F)
## boundary (singular) fit: see help('isSingular')
anova(m1, m2)
## Data: dat2
## Subset: Struct == "cvcv"
## Models:
## m1: V1 ~ Quantity + (1 | SP)
## m2: V1 ~ Quantity + (1 + Quantity | SP)
##    npar   AIC   BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## m1    5 12637 12663 -6313.5    12627                         
## m2   10 12559 12611 -6269.6    12539 87.635  5  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Juhuslik kalle on oluline. Vaatame mudeli väljundit

summary(m2)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: V1 ~ Quantity + (1 + Quantity | SP)
##    Data: dat2
##  Subset: Struct == "cvcv"
## 
##      AIC      BIC   logLik deviance df.resid 
##  12559.3  12611.3  -6269.6  12539.3     1335 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8905 -0.6038 -0.0848  0.5112  5.8316 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr     
##  SP       (Intercept)  60.51    7.779            
##           Quantity2   236.33   15.373   0.87     
##           Quantity3   260.67   16.145   0.37 0.78
##  Residual             619.61   24.892            
## Number of obs: 1345, groups:  SP, 21
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)   57.465      1.947 22.172   29.51  < 2e-16 ***
## Quantity2     56.704      3.981 23.510   14.24 4.70e-13 ***
## Quantity3     82.215      4.218 16.897   19.49 5.12e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) Qntty2
## Quantity2 0.551        
## Quantity3 0.156  0.616 
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

Ja vaatame juhuslikke efekte:

ranef(m2)
## $SP
##      (Intercept)  Quantity2   Quantity3
## SP1  -10.3494542 -19.125377 -10.5251008
## SP10  10.5354382  20.922136  13.5982954
## SP11  -3.4889768  -6.642542  -3.9353321
## SP12  -3.4230339  -8.623142  -8.0412354
## SP13  10.0191260  11.496970  -3.7401192
## SP14   5.8237724  19.720315  23.7028965
## SP15  -2.1819277   6.616913  18.9172514
## SP16  19.8104599  37.861742  22.6332461
## SP17  -4.9076920 -12.353419 -11.5094722
## SP18  -1.4339647  -2.162829  -0.4915445
## SP19  -3.0321872 -10.818933 -13.4355273
## SP2   -9.3539776 -11.872582   1.2313868
## SP20   2.0350326  -6.018883 -17.3408825
## SP21  -1.0760245  -2.736578  -2.5791713
## SP3   -7.4199823 -17.402320 -14.8708023
## SP4   -6.1005736  -9.375062  -2.4358574
## SP5   -3.9831694  -5.155547   0.3260958
## SP6   -0.1277132   6.150779  12.5466269
## SP7   -1.5409553  -5.727450  -7.2829918
## SP8    1.9903220  -8.726964 -22.5973654
## SP9    8.2054812  23.972774  25.8296033
## 
## with conditional variances for "SP"
sjPlot::plot_model(m2, type = "re")

Juhuslike kalletega joonis:

# baasgraafikaga
plot(V1 ~ Quantity, data = dat2, subset = Struct == "cvcv")
coefs <- summary(m2)$coefficients[,"Estimate"]
for(i in 1:nrow(ranef(m2)$SP)) lines(y = coefs["(Intercept)"] + 
                                         c(ranef(m2)$SP[i,"(Intercept)"], 
                                           ranef(m2)$SP[i,"Quantity2"] + coefs["Quantity2"], 
                                           ranef(m2)$SP[i,"Quantity3"] + coefs["Quantity3"]), 
                                     x = 1:3, 
                                     col = i)
# ggplotiga sama
ranef(m2)$SP %>%
    mutate(koneleja = rownames(.),
           Q1 = coefs["(Intercept)"] + `(Intercept)`,
           Q2 = coefs["(Intercept)"] + coefs["Quantity2"] + Quantity2,
           Q3 = coefs["(Intercept)"] + coefs["Quantity3"] + Quantity3) %>%
    select(koneleja, Q1, Q2, Q3) %>%
    pivot_longer(., cols = starts_with("Q"), names_to = "välde", values_to = "kõnetempo", names_prefix = "Q") -> df2

ggplot() +
    geom_boxplot(data = dat2 %>% filter(Struct == "cvcv"), aes(x = Quantity, y = V1)) +
    geom_line(data = df2, aes(x = välde, y = kõnetempo, color = koneleja, group = koneleja), 
              show.legend = F)

Mitu juhuslikku faktorit

Lisaks keelejuhi efektile võib korduval sõnal olla efekt.

m3 <- lmer(V1 ~ Quantity + (1|SP) + (1|Word), data = dat2, subset = Struct == "cvcv", REML=F)

Võrdleme seda mudeliga, kus oli ainult üks juhuslik vabaliige (1|SP)

anova(m1, m3)
## Data: dat2
## Subset: Struct == "cvcv"
## Models:
## m1: V1 ~ Quantity + (1 | SP)
## m3: V1 ~ Quantity + (1 | SP) + (1 | Word)
##    npar   AIC   BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## m1    5 12637 12663 -6313.5    12627                         
## m3    6 12420 12452 -6204.1    12408 218.63  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sõna juhuslik efekt on oluline. Mida see tähendab? Siin mudelis on siis nüüd 2 juhuslikku vabaliiget, st et vabaliiget kohandatakse iga korduva sõna JA iga korduva kõneleja suhtes. Sisuliselt see tähendab, et nii nagu igal kõnelejal võib olla individuaalne kõnetempol võib mõni sõna (näteks sageduse või mingi häälikukombinatsiooni koartikulatsiooni tõttu) olla regulaarselt lühem või pikem kui teised.

summary(m3)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: V1 ~ Quantity + (1 | SP) + (1 | Word)
##    Data: dat2
##  Subset: Struct == "cvcv"
## 
##      AIC      BIC   logLik deviance df.resid 
##  12420.3  12451.5  -6204.1  12408.3     1339 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9578 -0.5067 -0.0435  0.4552  7.3923 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Word     (Intercept) 269.7    16.42   
##  SP       (Intercept) 141.3    11.89   
##  Residual             472.1    21.73   
## Number of obs: 1345, groups:  Word, 262; SP, 21
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   61.084      3.181  36.290   19.20   <2e-16 ***
## Quantity2     56.312      3.160 299.892   17.82   <2e-16 ***
## Quantity3     87.331      3.263 322.127   26.77   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) Qntty2
## Quantity2 -0.314       
## Quantity3 -0.294  0.489

Käsuga ranef() saab juhuslike faktorite efekte vaadata:

ranef(m3)
sjPlot::plot_model(m3, type = "re", sort.est = T, grid = F)
## [[1]]

## 
## [[2]]

Muidugi ühes mudelis võib olla ka mitu juhuslikku faktorit, millest mõnel on ka juhuslikud kalded ja teisel ei ole, aga peame ka arvestama, et iga juhusliku efekti lisamine võtab fikseeritud efektidelt seletusjõudu vähemaks.

ANOVA segamudeliga

Segamudelist saab ka ANOVA väljundi, et kirjeldada faktori üldefekti olulisust.

anova(m3)
## Type III Analysis of Variance Table with Satterthwaite's method
##          Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Quantity 352128  176064     2 365.49  372.97 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Siit võib siis väita, et vältel on vokaali kestusele oluline efekt F(2, 365.49) = 372.97, p < 0.001.

Post-hoc test

Ja samuti peaks tegema post-hoc testi segamudeli puhul, kui on faktor, millel on rohkem kui kaks taset ja tahame võrrelda ka mittebaastasemetevahelisi erinevusi. Tukey HSD testi tegemiseks saab kasutada näiteks multcomp paketti.

#install.packages("multcomp")
library(multcomp)
summary(glht(m3, linfct = mcp(Quantity = "Tukey")), test = adjusted("holm"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = V1 ~ Quantity + (1 | SP) + (1 | Word), data = dat2, 
##     REML = F, subset = Struct == "cvcv")
## 
## Linear Hypotheses:
##            Estimate Std. Error z value Pr(>|z|)    
## 2 - 1 == 0   56.312      3.160  17.819   <2e-16 ***
## 3 - 1 == 0   87.331      3.263  26.766   <2e-16 ***
## 3 - 2 == 0   31.018      3.249   9.547   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- holm method)

Post-hoc testi tulemusel võib öelda, et vokaali kestus on oluliselt erinev kõigis kolmes vältes p<0.001.

R-ile ja pakettidele viitamine teadustöös

Kui kasutad oma töös R-i ja selle pakette, siis peaks neile ka viitama, esiteks selle pärast, et meetodi kirjelduses peab täpselt raporteerima, mis vahenditega analüüs on läbi viidud, teiseks selle pärast, et enamik R-i pakette on valminud teadustöö tulemusena ja nende autoreid hinnatakse bibliomeetria ehk viitamiste arvu järgi. Enamasti on paketi metainfos öeldud, kuidas sellele paketile võiks viidata ja R-is saab viite kerge vaevaga kätte käsuga citation(). Näiteks:

# baaspakett
citation()
## To cite R in publications use:
## 
##   R Core Team (2023). _R: A Language and Environment for Statistical
##   Computing_. R Foundation for Statistical Computing, Vienna, Austria.
##   <https://www.R-project.org/>.
## 
## A BibTeX entry for LaTeX users is
## 
##   @Manual{,
##     title = {R: A Language and Environment for Statistical Computing},
##     author = {{R Core Team}},
##     organization = {R Foundation for Statistical Computing},
##     address = {Vienna, Austria},
##     year = {2023},
##     url = {https://www.R-project.org/},
##   }
## 
## We have invested a lot of time and effort in creating R, please cite it
## when using it for data analysis. See also 'citation("pkgname")' for
## citing R packages.
# LME4 
citation("lme4")
## To cite lme4 in publications use:
## 
##   Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker (2015).
##   Fitting Linear Mixed-Effects Models Using lme4. Journal of
##   Statistical Software, 67(1), 1-48. doi:10.18637/jss.v067.i01.
## 
## A BibTeX entry for LaTeX users is
## 
##   @Article{,
##     title = {Fitting Linear Mixed-Effects Models Using {lme4}},
##     author = {Douglas Bates and Martin M{\"a}chler and Ben Bolker and Steve Walker},
##     journal = {Journal of Statistical Software},
##     year = {2015},
##     volume = {67},
##     number = {1},
##     pages = {1--48},
##     doi = {10.18637/jss.v067.i01},
##   }
citation("lmerTest")
## To cite lmerTest in publications use:
## 
##   Kuznetsova A, Brockhoff PB, Christensen RHB (2017). "lmerTest
##   Package: Tests in Linear Mixed Effects Models." _Journal of
##   Statistical Software_, *82*(13), 1-26. doi:10.18637/jss.v082.i13
##   <https://doi.org/10.18637/jss.v082.i13>.
## 
## A BibTeX entry for LaTeX users is
## 
##   @Article{,
##     title = {{lmerTest} Package: Tests in Linear Mixed Effects Models},
##     author = {Alexandra Kuznetsova and Per B. Brockhoff and Rune H. B. Christensen},
##     journal = {Journal of Statistical Software},
##     year = {2017},
##     volume = {82},
##     number = {13},
##     pages = {1--26},
##     doi = {10.18637/jss.v082.i13},
##   }
# mgcv
citation("mgcv")
## 2011 for generalized additive model method; 2016 for beyond exponential
## family; 2004 for strictly additive GCV based model method and basics of
## gamm; 2017 for overview; 2003 for thin plate regression splines.
## 
##   Wood, S.N. (2011) Fast stable restricted maximum likelihood and
##   marginal likelihood estimation of semiparametric generalized linear
##   models. Journal of the Royal Statistical Society (B) 73(1):3-36
## 
##   Wood S.N., N. Pya and B. Saefken (2016) Smoothing parameter and model
##   selection for general smooth models (with discussion). Journal of the
##   American Statistical Association 111:1548-1575.
## 
##   Wood, S.N. (2004) Stable and efficient multiple smoothing parameter
##   estimation for generalized additive models. Journal of the American
##   Statistical Association. 99:673-686.
## 
##   Wood, S.N. (2017) Generalized Additive Models: An Introduction with R
##   (2nd edition). Chapman and Hall/CRC.
## 
##   Wood, S.N. (2003) Thin-plate regression splines. Journal of the Royal
##   Statistical Society (B) 65(1):95-114.
## 
## To see these entries in BibTeX format, use 'print(<citation>,
## bibtex=TRUE)', 'toBibtex(.)', or set
## 'options(citation.bibtex.max=999)'.
# sjPlot
citation("sjPlot")
## To cite package 'sjPlot' in publications use:
## 
##   Lüdecke D (2023). _sjPlot: Data Visualization for Statistics in
##   Social Science_. R package version 2.8.15,
##   <https://CRAN.R-project.org/package=sjPlot>.
## 
## A BibTeX entry for LaTeX users is
## 
##   @Manual{,
##     title = {sjPlot: Data Visualization for Statistics in Social Science},
##     author = {Daniel Lüdecke},
##     year = {2023},
##     note = {R package version 2.8.15},
##     url = {https://CRAN.R-project.org/package=sjPlot},
##   }

Kordamisküsimused

Kordamisküsimuste testi saad teha ka Moodle’is, seal näed peale vastamist ka õigeid vastuseid ja selgitusi.

1. Juhuslik faktor on juhuslik selle poolest, et

  1. selle mõju uuritavale tunnusele ei ole oluline
  2. meil on selle mõjust suva
  3. meil on valimis juhuslik valik populatsiooni esindajatest
  4. selle mõju ei ole võimalik kirjeldada

2. Varasemates praktikumides kasutatud andmesti ldt.csv on töödeldud nii, et seal on sõltumatud mõõtmised, sest iga sõna kohta on arvutatud katseisikute keskmine reaktsiooniaeg. Eelmises praktikumis leidsime, et reaktsiooniaega mõjutab sõna pikkus ja sagedus, aga nende vahel ei ole olulist interaktsiooni. Kui meil oleks kasutada katseandmestik sellisel kujul, et iga rida tabelis on ühe katsisiku reaktsiooniaeg vastavale sõnale, siis peaks segamudelis testima ka

  1. peaks valima, kas katseisiku mõju või sõna mõju
  2. sõna mõju juhusliku pikkuse ja sageduse kaldena, katseisikut vabaliikmena
  3. nii sõna kui katseisiku mõju juhusliku faktorina
  4. sõna juhuliku vabaliikmena ja sõna juhuliku kaldena

3. Oletame, et reaktsiooniaega mõjutab sõna pikkus ja sagedus, mille vahel ei ole interaktsiooni. Mudel arvestab sõna mõjuga reaktsioooniajale ning katseisikute reaktsoonikiirus on erinev sõltuvalt sõna sagedusest. Milline selle mudeli valem võiks olla?

  1. RT ~ pikkus * sagedus + (1|katsesiik) + (1|sõna)
  2. RT ~ pikkus + sagedus + (1|katsesiik) + (1|sõna)
  3. RT ~ pikkus + sagedus + (1+sagedus|katsesiik) + (1+pikkus|sõna)
  4. RT ~ pikkus + sagedus + (1+sagedus|katsesiik) + (1|sõna)

4. Kuidas hinnata lineaarse segamudeli fikseeritud faktori olulisust?

  1. tuleb võrrelda ANOVA testiga mudelit, kus faktor on sees, ja mudelit, kust see on välja jäetud
  2. kui faktori p > 0.05, siis on oluline
  3. kui faktori |t| > 2, siis on oluline
  4. kui faktori Estimate > 0.2, siis on oluline

5. Kui meil on vaja võrrelda rohkem kui 2 tasemega faktorite kõiki tasemeid

  1. peame tegema post-hoc testi
  2. teeme lihtsalt testi mitu korda, vahetame baastaset
  3. post-hoc testi tegemine on keeruline ja seda tuleks vältida
  4. segamudeliga ei peagi post-hoc testi tegema

Järgmisel korral

  • Logistiline regressioon e uuritav tunnus on nominaalne (kahe tasemega).

Funktsioonide sõnastik

  • lme4::lmer() – lineaarne segamudel paketiga lme4
  • lmer::ranef() – segamudeli juhuslikud efektid
  • residuals() – mudeli jäägid (sama mis mudel$residuals)
  • mgcv::gam() – segamudel mgcv paketiga
  • sjPlot::plot_model – mudeli efektide visualiseerimine