Logistilise regressiooni mudelid ja segamudelid

Tänases praktikumis

  • Mudelite headuse näitajad
  • Logistilise regressiooni eelduste kontroll
  • (Mudeli valideerimine)
  • Segamudelid

Mudelite headuse näitajad

Eelmisel korral lõpetasime tasakaalus alalütleva ja peal-kaassõna andmestiku peal tehtud mudeliga cx ~ LM_mobiilsus + LM_komplekssus + tegusõna_klass + LM_pikkus_log * murderühm.

library(tidyverse)
library(sjPlot) # mudeli visualiseerimiseks
ade_peal <- read.delim("data/ade_peal.csv", encoding = "UTF-8", stringsAsFactors = TRUE)
ade <- droplevels(ade_peal[ade_peal$cx == "ade",]) # ainult alalütleva käände andmed
peal <- droplevels(ade_peal[ade_peal$cx == "peal",]) # ainult peal-kaassõna andmed

set.seed(1) # seeme juhusliku valiku tegemise korratavuseks
juh_id <- sample(x = 1:nrow(peal), # võta ridadelt 1 kuni 1310 juhuslik valim,
                 size = nrow(ade), # mille suurus oleks sama, mis ade-andmestikul
                 replace = F) # ära korda ühtegi vaatlust

peal_juhuvalim <- peal[juh_id,] # võta peal-andmestikust ainult juhuslikult valitud read
dat <- rbind(ade, peal_juhuvalim) # liida kaks andmestikku ridupidi uuesti kokku

Teeme sel korral sama mudeli, aga määrame kõigepealt kategoriaalsete tunnuste baastasemeteks nende kõige tüüpilisema (= sagedama) väärtuse ehk moodi.

sapply(dat[,c("LM_mobiilsus", "LM_komplekssus", "tegusõna_klass", "murderühm")], # iga tulba/tunnuse kohta neljast
       function(x) names(which.max(table(x)))) # leia selle tunnuse mood
##   LM_mobiilsus LM_komplekssus tegusõna_klass      murderühm 
##   "staatiline"     "lihtsõna"   "olemisverb"         "Kesk"

Tüüpiline kontekst oleks niisiis midagi sellist nagu auto on teel / tee peal.

dat %>%
  mutate(LM_mobiilsus = relevel(LM_mobiilsus, ref = "staatiline"),
         LM_komplekssus = relevel(LM_komplekssus, ref = "lihtsõna"),
         tegusõna_klass = relevel(tegusõna_klass, ref = "olemisverb"),
         murderühm = relevel(murderühm, ref = "Kesk")) -> dat2
m5.glm <- glm(cx ~ LM_mobiilsus + LM_komplekssus + tegusõna_klass + LM_pikkus_log*murderühm, data = dat2, family = "binomial")
plot_model(m5.glm, 
           show.values = T, # näita koefitsientide väärtuseid
           value.offset = 0.3, # punktist veidi ülevalpool
           sort.est = T, # sorteeri koefitsiendid
           p.adjust = "holm", # korrigeeri p-väärtuseid enam kui kahe tasemega kat. tunnuste võrdluseks
           title = "Baaskontekst:\nLM mobiilsus [staatiline], LM komplekssus [lihtsõna],\nLM pikkus log [0 ehk 1 silp], murderühm [Kesk]",
           wrap.title = 150) + # luba pealkirjas ühele reale kuni 150 tähemärki
  theme(plot.title = element_text(size = 12)) # muuda pealkirja teksti suurust ggplotiga 

# Visualiseerime peamõjusid
plot_model(m5.glm, type = "pred", terms = c("LM_mobiilsus"), axis.title = c("LM mobiilsus", "'peal' ennustatud tõenäosus"))

plot_model(m5.glm, type = "pred", terms = c("LM_komplekssus"), axis.title = c("LM komplekssus", "'peal' ennustatud tõenäosus"))

plot_model(m5.glm, type = "pred", terms = c("tegusõna_klass"), axis.title = c("Tegusõna klass", "'peal' ennustatud tõenäosus"))

Võime paremaks võrdluseks määrata kõikidel ka y-telje 0 kuni 100 protsendini.

plot_model(m5.glm, type = "pred", terms = c("LM_mobiilsus"), axis.title = c("LM mobiilsus", "'peal' ennustatud tõenäosus")) +
  scale_y_continuous(limits = c(0, 1), labels = scales::percent)
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

plot_model(m5.glm, type = "pred", terms = c("LM_komplekssus"), axis.title = c("LM komplekssus", "'peal' ennustatud tõenäosus")) +
  scale_y_continuous(limits = c(0, 1), labels = scales::percent)
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

plot_model(m5.glm, type = "pred", terms = c("tegusõna_klass"), axis.title = c("Tegusõna klass", "'peal' ennustatud tõenäosus")) +
  scale_y_continuous(limits = c(0, 1), labels = scales::percent)
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

# Visualiseerime interaktsiooni
plot_model(m5.glm, type = "int", grid = T)

plot_model(m5.glm, type = "int", grid = T) +
  scale_y_continuous(limits = c(0, 1), labels = scales::percent)
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

AIC ja anova aitasid meil leida olulisima mudeli, st mudeli, kus iga tunnus panustab piisaval määral uuritava tunnuse seletamisse selleks, et seda mudelis sees hoida. Nüüd vaatame ka, kui hästi selline mudel uuritava tunnuse varieerumist seletada suudab.
Logistilises regressioonis ei pakuta meile R-ruutu, mille põhjal mudeli headust ja täpsust hinnata. Küll aga kasutatakse mõnikord erinevaid nn pseudo-R2 variante, nagu nt Nagelkerke R2, Cox-Snelli R2, McFaddeni R2, või Tjuri R2 (mida näidatakse ka nt sjPlot::tab_model() või report::report_text() funktsioonide väljundis). Neid võib pidada ligikaudseteks mudeli headuse näitajateks, ent need on üldjuhul oluliselt madalamad kui lineaarsete mudelite R-ruudud, samuti ei ole nende tõlgendamine nii ühemõtteline, kuna eri pseudo-R-ruudud on erinevalt arvutatud. Sellegipoolest võivad need olla kasulikud inkrementaalsel modelleerimisel samal andmestikul tehtud mudelite võrdlemiseks.

Publitseerimisel raporteeritakse pseudo-R-ruudust enam klassifitseerimistäpsust ja C-indeksit/AUC-d. Nendest järgmiseks pikemalt.

Klassifitseerimistäpsus

Iga vaatluse kohta andmestikus ennustab mudel mingi tõenäosuse, millega peal-konstruktsioon mudelisse kaasatud seletavate tunnuste põhjal esineda võiks. Sisuliselt on mudelis peal ja ade kodeeritud vastavalt kui 1 (100% peal-konstruktsioon) ja 0 (0% peal-konstruktsioon). Kui meil on andmestikus rida, kus cx = "peal" ja mudel ennustab sellele vaatlusele peal-konstruktsiooni esinemise tõenäosuse, mille väärtus on > 0.5, siis ennustab mudel põhimõtteliselt õigesti, sest peab sellises kontekstis peal-konstruktsiooni esinemist tõenäolisemaks kui käändekonstruktsiooni esinemist. Kui ennustab väärtuse, mis on < 0.5, siis valesti. Mudeli klassifitseerimistäpsus näitab nende juhtude osakaalu, kus mudel ennustab vaatlusele õige klassi/taseme.

Ennustatud tõenäosused saab mudelist kätte käsuga mudel$fitted.values.

# Vaatame ennustusi esimesele 6 vaatlusele/reale
head(m5.glm$fitted.values) # tõenäosused
##         6         7         8         9        11        13 
## 0.2593372 0.3959104 0.1963169 0.6000458 0.6693924 0.7146624
head(m5.glm$linear.predictors) # vastavad šansi logaritmid
##          6          7          8          9         11         13 
## -1.0494165 -0.4225345 -1.4094748  0.4056559  0.7054383  0.9181372
# Võrdleme seda andmestiku tulba "cx" tegelike väärtustega
head(dat2$cx) # esimesed 6 vaatlust on kõik "ade"
## [1] ade ade ade ade ade ade
## Levels: ade peal
# Tegelikud väärtused (kujul 0 ja 1) on ka mudelis olemas (y viitab uuritavale tunnusele)
head(m5.glm$y) # esimesed 6 vaatlust on kõik 0
##  6  7  8  9 11 13 
##  0  0  0  0  0  0

Näeme, et esimesed 6 vaatlust on käändekonstruktsiooni kasutusjuhud. Kuna panime oma tasakaalus valimi kokku nii, et kõigepealt ade andmestik ja siis peal andmestik, on tegelikult kõik 722 esimest vaatlust käändekonstruktsioonid.
Esimese 6 puhul ennustab mudel 3 juhul õigesti (vaatlused 1-3) ja 3 juhul valesti (vaatlused 4-6), sest viimastele ennustatud peal-konstruktsiooni esinemise tõenäosus on suurem kui 0.5, ehkki tegelikult on andmestiku järgi seal kontekstis “ade”.

Võib kasutada ka funktsiooni predict(), mis töötab igasuguste mudelitega ning võimaldab küsida siin nii ennustatud tõenäosusi kui ka šansi logaritme. predict()-funktsioon on eriti kasulik siis, kui meil on eraldi treeningandmestik, mille põhjal mudelit ehitada, ning testandmestik, mille peal mudelit hinnata. Vaatame seda mudeli valideerimise juures.

# Ennustatud tõenäosused
head(predict(m5.glm, type = "response"))
##         6         7         8         9        11        13 
## 0.2593372 0.3959104 0.1963169 0.6000458 0.6693924 0.7146624
# Ennustatud šansi logaritmid
head(predict(m5.glm, type = "link"))
##          6          7          8          9         11         13 
## -1.0494165 -0.4225345 -1.4094748  0.4056559  0.7054383  0.9181372

Kui tahame saada mudeli klassifitseerimistäpsust, on meil vaja teada kolme arvu:

  • Kõikide vaatluste arvu.
  • Nende vaatluste hulka, kus ennustatud tõenäosus on > 0.5 JA tegelik väärtus andmestikus on cx = "peal". Neid nimetatakse ka tõeliselt positiivseteks juhtudeks (true positives, TP).
  • Nende vaatluste hulka, kus ennustatud tõenäosus on < 0.5 JA tegelik väärtus andmestikus on cx = "ade". Neid nimetatakse tõeliselt negatiivseteks juhtudeks (true negatives, TN).

Vastavaid valesid ennustusi nimetatakse seega valepositiivseteks (false positives, FP - mudel ennustab > 0.5, tegelik vaatlus on “ade” ehk 0) ja valenegatiivseteks (false negatives, FN - mudel ennustab < 0.5, tegelik vaatlus on “peal” ehk 1).

Klassifitseerimistäpsuse leidmiseks on erinevaid viise. Näiteks

# vaatlustele ennustatud tõenäosused
enn <- m5.glm$fitted.values
head(enn)
##         6         7         8         9        11        13 
## 0.2593372 0.3959104 0.1963169 0.6000458 0.6693924 0.7146624
# "cx" tulp andmestikus
teg <- as.character(dat2$cx)
head(teg)
## [1] "ade" "ade" "ade" "ade" "ade" "ade"
# tee uus objekt, milles on sama palju ridu
# kui dat2 andmestikus vaatlusi/ridu ning
# milles esineb "peal", 
# kui vastav väärtus enn-objektis on > 0.5, 
# ja "ade", kui vastav väärtus enn-objektis on < 0.5
enn_bin <- ifelse(enn > 0.5, "peal", "ade")
head(enn_bin)
##      6      7      8      9     11     13 
##  "ade"  "ade"  "ade" "peal" "peal" "peal"
# tee tegelikest ja ennustatud binaarsetest 
# väärtustest sagedustabel, kus on
# tõeliselt positiivsed, valepositiivsed,
# tõeliselt negatiivsed ja valenegatiivsed
# juhud
(klass_tabel <- table(teg, enn_bin))
##       enn_bin
## teg    ade peal
##   ade  489  233
##   peal 220  502
# leia klassifitseerimistäpsus
# selleks liida kokku juhud, kus tegelik vaatlus oli "ade"
# ja ka mudel ennustas "ade", ning juhud, kus tegelik
# vaatlus oli "peal" ja ka mudel ennustas "peal" 
# (= õigete ennustuste hulk), 
# ning jaga tulemus vaatluste arvuga
(klass_acc <- klass_tabel["ade","ade"]+klass_tabel["peal","peal"])/sum(klass_tabel)
## [1] 0.6862881
# või
# sum(diag(klass_tabel))/sum(klass_tabel)

Mudel suudab seega klassifitseerida õigesti 68.6% andmestikus olevatest vaatlustest.

Oletame, et otsustame, et meil ei ole mingeid seletavaid tunnuseid vaja ja saame hakkama ka sedasi, et ennustame uuritava tunnuse vaatluste klasse, pakkudes alati lihtsalt sagedamat varianti. Kuna praegusel juhul on andmestik tasakaalus, siis võiksime ennustada kogu aeg ükskõik kas ainult “peal”-konstruktsiooni või ainult “ade”-konstruktsiooni, mõlematpidi oleks meie naiivne klassifitseerimistäpsus 0.5 ehk 50%.

Seletavate tunnustega mudel ennustab seega ~18.6% paremini kui ühegi tunnuseta mudel. See ei ole väga hea tulemus, aga midagi siiski.

TP, FP, TN ja FN risttabelit (klass_tabel) nimetatakse mudelite hindamisel sageli nimega eksimismaatriks (confusion matrix) ja selle põhjal saab tegelikult arvutada veel mitmeid klassifitseerimisheaduse näitajaid, vt rohkem siit.

Näiteks võib veel joonistada nn ROC-ruumi (ROC space). ROC-lühend tuleb fraasist receiver operating characteristic, see illustreerib binaarse klassifitseerimismudeli ennustusvõimet ning seda saab leida TPR-i (true positive rate) ja FPR-i (false positive rate) järgi. TPR on “peal” ennustuste osakaal kõikidest “peal” tegelikest vaatlustest ning FPR on “peal” ennustuste osakaal kõikidest tegelikest “ade” vaatlustest. TPR on inglise keeles ka sensitivity või recall, eesti keeles saagis.

klass_tabel
##       enn_bin
## teg    ade peal
##   ade  489  233
##   peal 220  502
# TPR ja FPR
TPR <- klass_tabel["peal", "peal"]/sum(klass_tabel["peal",]) # kui palju peal-konstruktsioone mudel õigesti tuvastas?
FPR <- klass_tabel["ade", "peal"]/sum(klass_tabel["ade",]) # kui palju ade-konstruktsioonidest mudel peal-konstruktsioonideks pidas?

plot(FPR, TPR, xlim = c(0,1), ylim = c(0,1), main = "ROC-ruum")
lines(x = 0:1, y = 0:1, lty = "dashed", col = "red")

Punane katkendjoon on nn ROC-kurv, millel asuvad näitajad siis, kui mudel ennustab juhuslikult “1”/“peal” või “0”/“ade”, sest kui võtta mudelist välja kõik “peal” vaatlustele ennustatud tõenäosused ning kõik “ade” vaatlustele ennustatud tõenäosused, siis need tõenäosusjaotused kattuvad täielikult. Head näitajad on need, mis jäävad ROC-ruumis punasest diagonaalsest joonest ülespoole, sest need näitavad paremat klassifitseerimisvõimet kui juhuslik arvamine. Mida lähemal on täpp x-teljel 0-le ja y-teljel 1-le, seda paremini mudel vaatlusi klassifitseerib. Mida lähemal on täpp x-teljel 1-le ja y-teljel 0-le, seda halvemini mudel vaatlusi klassifitseerib (mudel ennustab alati vaatlusele “1” väärtust “0” ja vastupidi).

C-indeks ehk AUC

C ehk konkordantsiindeks (ka Area Under the ROC Curve ehk AUC. NB! mitte ajada segi AIC-ga!) näitab logistilise mudeli eristusvõimet: kui hästi suudab mudel uuritava tunnuse klasse eristada. Teisisõnu näitab C nende kordade proportsiooni, mil kõikide võimalike vaatluste paaride korral ennustab meie mudel suuremat sündmuse “1” esinemise tõenäosust vaatlusele, millel on ka päriselt andmetes sündmus “1”, ja suuremat sündmuse “0” esinemise tõenäosust vaatlusele, millel on ka päriselt andmetes sündmus “0”.

Näiteks kahe vaatluse puhul, millest ühe tegelik väärtus on peal ja teise oma ade, võib mudel küll mõlemale ennustada peal esinemisele väiksema tõenäosuse kui 0.5, ent kui esimese puhul on see tõenäosus suurem kui teise puhul, eristab mudel neid siiski õigesti.

C-d peetakse niisiis täpsemaks ja paindlikumaks mudeli eristusvõime hindajaks kui klassifitseerimistäpsust, kuna see arvestab ennustusi tõenäosustena ehk pidevana, samas kui klassifitseerimistäpsus jagab ennustused lihtsalt 0.5 juurest kahte klassi, hoolimata sellest, kas ennustatud tõenäosus on nt 0.1 või 0.49.

Indeksi skaalat võib tõlgendada nii:

  • C = 0.5: mudel ei erista, 50-50 tõenäosus täppi panna
  • 0.5 < C <= 0.7: kehv eristusvõime
  • 0.7 < C <= 0.8: rahuldav eristusvõime
  • 0.8 < C <= 0.9: väga hea eristusvõime
  • C > 0.9: suurepärane eristusvõime

Indeksi arvutamiseks võib kasutada näiteks pakette pROC või Hmisc.

# install.packages(c("pROC", "Hmisc"))
library(pROC)
auc(dat2$cx, # tegelik uuritava tunnuse tulp andmestikus
    m5.glm$fitted.values) # igale vaatlusele ennustatud "peal" tõenäosused
## Area under the curve: 0.746
library(Hmisc)
somers2(m5.glm$fitted.values, # igale vaatlusele ennustatud "peal" tõenäosused
        as.numeric(dat2$cx)-1) # 0ks ja 1ks teisendatud tegelik uuritava tunnuse tulp andmestikus
##            C          Dxy            n      Missing 
##    0.7459629    0.4919257 1444.0000000    0.0000000

Kuna meie mudeli C on 0.746, võime pidada mudeli eristusvõimet rahuldavaks.

Kui mudelil on ideaalne eristusvõime ja see suudab täiuslikult eristada 1sid ja 0e (“peal”-konstruktsioone ja “ade”-konstruktsioone), siis kahe klassi tõenäosusjaotused ei kattu ja ROC-kurv on täisnurkne. Kui mudelil ei ole mingit eristusvõimet (kahe klassi tõenäosusjaotused kattuvad täielikult), siis on ROC-kurv diagonaalne sirge ROC-ruumis. Kui mudel ennustab tegelikele vaatlustele täpselt vastupidiseid klasse, siis on ROC-kurv teistpidi täisnurkne (kahe klassi tõenäosusjaotused ei kattu, aga need on vastupidised).

Meie mudeli uuritava tunnuse klassidele ennustatud tõenäosuste jaotused näevad välja umbes nii:

data.frame(teg = dat2$cx, enn = m5.glm$fitted.values) %>%
ggplot(aes(x = enn, color = teg)) +
  geom_density(alpha = 0.5, linewidth = 1) +
  scale_color_manual(values = c("red", "darkgreen")) +
  geom_vline(xintercept = 0.5, linetype = "dashed") +
  labs(x = "Ennustatud tõenäosused", 
       y = "Tõenäosustihedus", 
       title = "Ennustatud tõenäosuste jaotus",
       color = "Tegelik klass")

Ja ROC-kurv nii:

library(ROCR)
prediction(m5.glm$fitted.values, dat2$cx) %>%
  performance(., "tpr", "fpr") -> perf

data.frame(x = unlist(perf@x.values), y = unlist(perf@y.values)) %>%
ggplot(aes(x = x, y = y)) +
  geom_line(color = "red") +
  labs(x = perf@x.name, y = perf@y.name, title = "ROC-kurv")

AUC ehk Area Under the ROC-curve näitab niisiis selle ala suurust, mis jääb joonest alla ehk x-telje poole.

Mudeli headuse näitajate kohta loe veel siit.

Mudeli eelduste kontroll ja diagnostika

Ehkki logistiline regressioon ei eelda nt jääkide normaaljaotust, on tegemist siiski parameetrilise mudeliga (milles hinnatavad parameetrid on mudeli väljundi koefitsiendid ehk regressioonikordajad).

Logistilisel regressioonil on mõned eeldused:

  • vaatlused on üksteisest sõltumatud,
  • suhe logaritmitud šansi ja arvuliste seletavate tunnuste vahel on lineaarne,
  • seletavate tunnuste vahel ei esine multikollineaarsust.

Vaatluste sõltumatus

Vaatluste sõltumatus tähendaks, et meil on nt ühe kõneleja kohta andmestikus üksainus vaatlus/rida, kus ta siis kas on kasutanud “peal”- või käändekonstruktsiooni. Keeleteaduse andmestikud on harva sellised ning tihtipeale kasutatakse ära kõik vaatlused/lausungid, mida saada on. Meil on andmestikus tunnus kõneleja, mis identifitseerib unikaalsed kõnelejad. Selles nimes sisaldub failinimi, selles ka lindistatud kõneleja nimi. Harvadel juhtudel osaleb ühes lindistuses ka mitu kõnelejat, mispuhul samas failis kõnelejad on eristatud koodidega “KJ1” ja “KJ2” vmt. Seda tunnust võib pidada niisiis heaks indikaatoriks selle kohta, kui palju mingi kõneleja vaatlusi on panustanud.

head(sort(table(dat$kõneleja), decreasing = TRUE)) # kõige produktiivsemad kõnelejad
## 
##      KJ_HMD_Aleksander_Idarand_synt.xml KJ_PIL_Helene_Pikk_synt-parandamata.xml 
##                                      21                                      21 
## KJ_JOE_Jyri_Maisa_synt-kontrollitud.xml KJ_JOE_Leena_Toomar_synt-parandatud.xml 
##                                      19                                      19 
## KJ_PLT_Liina_Roosimaa_EMH1427_lihts.xml        KJ_SJN_Marie_Kees_EMH16_synt.xml 
##                                      19                                      19
summary(as.numeric(table(dat$kõneleja))) # kui palju kõnelejad üldiselt on vaatlusi panustanud?
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   2.000   4.000   5.049   7.000  21.000

Kõige produktiivsematelt keelejuhtidelt on andmestikus tervelt paarkond vaatlust, samas kui keskmiselt on kõneleja kohta 4-5 vaatlust. Seetõttu on esimene eeldus tegelikult rikutud.

Kui arvuliste tunnuste puhul saame kasutada sellisel juhul ka näiteks tavalist regressioonimudelit keskmistatud andmetega, siis kategoriaalsete tunnustega oleks kohane kasutada segamõjudega mudelit, millest räägime veidi aja pärast.

Lineaarne suhe šansi logaritmi ja arvuliste tunnuste vahel

Meil on mudelis üks arvuline tunnus, LM_pikkus_log, kusjuures tunnus ei ole peamõjuna, vaid osaleb murderühmaga interaktsioonis. Lineaarsuse kontrollimiseks sobib näiteks visuaalne vaatlus.

data.frame(logit = m5.glm$linear.predictors, 
           pikkus = dat$LM_pikkus_log, 
           ryhm = dat$murderühm) %>%
  ggplot(aes(x = pikkus, y = logit)) +
  geom_point(alpha = 0.1) +
  facet_wrap("ryhm") +
  geom_smooth(method = "loess")

Näeme, et suhe kohafraasi logaritmitud pikkuse ja šansi logaritmi vahel ei ole üheski murderühmas päris lineaarne, sest 2-, 3- või 4-silbiliste kohafraaside kohal toimuvad väikesed jõnksud. Samas ei ole jooned ka liiga kurvilised, nii et võime teatud ettevaatusega arvulise tunnuse interaktsiooni mudelisse alles jätta.

Seletavate tunnuste sõltumatus

Nagu lineaarse regressiooni mudelite puhul juba nägime, võib seletavate tunnuste enda vahel olla tugev korrelatsioon, mis tähendab, et tunnused seletavad (vähemalt osaliselt) sama osa uuritava tunnuse varieerumisest. Seda nimetatakse multikollineaarsuseks.

Multikollineaarsust saab kontrollida funktsiooniga vif() paketist car. Lühend VIF ehk Variation Inflation Factor hindab, kui palju mingi regressioonikoefitsiendi (= parameetri) varieeruvus kasvab, kui tunnused on omavahel seotud. Kui VIF on 1, siis ei korreleeru tunnus ühegi teise seletava tunnusega mudelis. VIF, mis on väiksem kui 5, näitab, et multikollineaarsust olulisel määral ei esine. 5-10 näitab juba sellist seotust, mis võib olla problemaatiline. Kui VIF on suurem kui 10, siis võib eeldada, et regressioonikoefitsientide hinnangud on tugevalt kallutatud, kuna esineb multikollineaarsust. Kui ühel tunnusel on kõrge VIF-skoor, peab olema veel vähemalt üks tunnus, millel on samamoodi kõrge skoor.

car::vif(m5.glm)
##                               GVIF Df GVIF^(1/(2*Df))
## LM_mobiilsus              1.040257  1        1.019930
## LM_komplekssus            1.083548  1        1.040936
## tegusõna_klass            1.045885  4        1.005624
## LM_pikkus_log             2.223639  1        1.491187
## murderühm                99.020461  3        2.150903
## LM_pikkus_log:murderühm 115.125115  3        2.205608

GVIF tähistab fraasi generalized VIF, mis sobitab VIF-skoorid ka üldistatud lineaarsetele mudelitele ning mida näidatakse siis, kui tunnuste vabadusastmete arv on > 1. Multikollineaarsus peamõjude ja interaktsioonide vahel on ootuspärane ja üldiselt ka vältimatu, kuna VIF-skoorid mõõdavad seda, kas kaks tunnust seletavad (osaliselt) sama osa varieeruvusest, aga interaktsioon juba oma olemuselt sisaldabki ka peamõjusid.
Siin väljastatakse ka GVIF-skoorid, mis kohandavad skoore vabadusastmete arvuga. Neid peaksimegi pigem vaatama.
Näeme, et siin tugevat multikollineaarsust ei esine. Kui esineks, oleks ohutum ühest või teisest seletavast tunnusest loobuda.

Saame kollineaarsust hinnata ka funktsiooni check_collinearity() abil paketist performance, aga see näitab meile ainult tavalisi VIF-skoore.

# install.packages("performance")
library(performance)
check_collinearity(m5.glm)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##            Term  VIF       VIF 95% CI Increased SE Tolerance Tolerance 95% CI
##    LM_mobiilsus 1.04 [  1.01,   1.16]         1.02      0.96     [0.86, 0.99]
##  LM_komplekssus 1.08 [  1.04,   1.17]         1.04      0.92     [0.85, 0.96]
##  tegusõna_klass 1.05 [  1.01,   1.15]         1.02      0.96     [0.87, 0.99]
##   LM_pikkus_log 2.22 [  2.06,   2.41]         1.49      0.45     [0.41, 0.49]
## 
## High Correlation
## 
##                     Term    VIF       VIF 95% CI Increased SE Tolerance
##                murderühm  99.02 [ 89.47, 109.60]         9.95      0.01
##  LM_pikkus_log:murderühm 115.13 [104.01, 127.43]        10.73  8.69e-03
##  Tolerance 95% CI
##      [0.01, 0.01]
##      [0.01, 0.01]

performance pakett pakub veel mitmeid funktsioone mudeli diagnostika tegemiseks, täpsuse leidmiseks jpm-ks ning seda saab kasutada mitmesuguste regressioonimudelite puhul.

Jäägid

Logistilise regressiooni mudeli jäägid (deviance residuals) näitavad erinevust vaatlustele ennustatud tõenäosuse (0 kuni 1) ja nende tõelise väärtuse (1 või 0) vahel. Jääkide arvutamiseks kasutatakse valemit sqrt(-2*log(ennustatud)), kui vaatluse väärtus on 1 (siin: “peale”) ja valemit -sqrt(-2*log(1-ennustatud)), kui vaatluse väärtus on 0 (siin: “ade”). Näiteks:

(t <- m5.glm$y[1]) # esimese rea vaatlus (adessiiv ehk 0)
## 6 
## 0
(e <- m5.glm$fitted.values[1]) # esimesele vaatlusele ennustatud tõenäosus olla "peal"-konstruktsioon
##         6 
## 0.2593372
residuals(m5.glm, type = "deviance")[1] # esimese vaatluse deviance residual
##          6 
## -0.7748674
# on sama mis
-sqrt(-2*log(1-e))
##          6 
## -0.7748674

Tahaksime, et mudeli jäägid ei oleks väga suured. Näiteks kui käändekonstruktsiooni (“ade”) vaatlusele ennustaks mudel 0.99 tõenäosust olla “peal”-konstruktsioon, oleks selle jääk -sqrt(-2*log(1-0.99))=-3.034854, ja kui peal-konstruktsiooni vaatlusele ennustaks mudel ainult 0.01 tõenäosust olla päris “peal”-konstruktsioon, oleks selle jääk sqrt(-2*log(0.01))=3.034854. Ka summary(m5.glm) näitab meile jääkide jaotuse ära. Meie mudeli kõige väiksem jääk on -2.14635, kõige suurem 2.25620 ja mediaan 0.07362.

# Kõige enam valesti ennustatud adesiivi vaatlus
dat2[which.min(residuals(m5.glm, "deviance")),]
##        TR_liik tegusõna tegusõna_klass LM_lemma LM_mobiilsus LM_komplekssus
## 461 verbifraas   tegema    tegevusverb     külg   liigutatav       lihtsõna
##     LM_pikkus LM_pikkus_log sõnajärg  cx murre murderühm
## 461         4             2    lm_tr ade Tartu     Lõuna
##                         kõneleja
## 461 KJ1_VON_Ann_Plaks_EMH502.xml
# Kõige enam valesti ennustatud peal-konstruktsiooni vaatlus
dat2[which.max(residuals(m5.glm, "deviance")),]
##      TR_liik tegusõna tegusõna_klass LM_lemma LM_mobiilsus LM_komplekssus
## 1137 asesõna    olema     olemisverb  maalapp   staatiline       liitsõna
##      LM_pikkus LM_pikkus_log sõnajärg   cx murre murderühm
## 1137         5      2.321928    tr_lm peal Kirde   Ranniku
##                                         kõneleja
## 1137 KJ_JOH_Heinrich_Kaart_synt-kontrollitud.xml

Logistilise regressiooni jäägid aga ei pea olema normaaljaotusega! Seetõttu on need kasulikud üldiselt pigem ebatavaliste vaatluste tuvastamiseks ja üldise ennustuste homogeensuse ja headuse hindamiseks.

Ebatavalisi vaatlusi võime tuvastada ka käsuga influencePlot() paketist car. Joonisel on näidatud

  • y-teljel jäägid. Need, mille absoluutväärtus on suurem kui 2, võivad olla probleemsed;
  • x-teljel nn hat-väärtused, mis näitavad, kui palju mingi vaatlus on oma seletavate tunnuste väärtuste poolest teistest erinev. Mida suurem väärtus on punktil x-teljel, seda “imelikum” vaatlus teistest on;
  • punktide suurusega nn Cooki kaugust, mis näitab vaatluse mõjukust (“kui palju mudel muutuks, kui selle vaatluse eemaldaksime?”). Mida mõjukam on vaatlus, seda suurema punktina seda joonisel kujutatakse.
car::influencePlot(m5.glm)

##        StudRes         Hat       CookD
## 461  -2.158860 0.005943756 0.003870346
## 1089 -1.186948 0.066234105 0.005183414
## 1840  1.799979 0.057865804 0.016453229
## 1137  2.285507 0.011202760 0.009613434

Kontrollida võiks andmestikust niisiis neid ridu, kus on korraga nt kas suur jääk ja suur Cooki kaugus või suur hat-väärtus ja suur Cooki kaugus, ja need andmestikust välja jätta, kui see on põhjendatud (alati ei ole).

dat2 %>%
  mutate(enn = m5.glm$fitted.values) %>%
  filter(rownames(.) %in% c("461", "1089", "1840", "1137"))
##         TR_liik tegusõna tegusõna_klass LM_lemma LM_mobiilsus LM_komplekssus
## 1    verbifraas   tegema    tegevusverb     külg   liigutatav       lihtsõna
## 2 nimisõnafraas    olema     olemisverb    talje   liigutatav       lihtsõna
## 3 nimisõnafraas    olema     olemisverb käsivars   liigutatav       liitsõna
## 4       asesõna    olema     olemisverb  maalapp   staatiline       liitsõna
##   LM_pikkus LM_pikkus_log sõnajärg   cx murre murderühm
## 1         4      2.000000    lm_tr  ade Tartu     Lõuna
## 2         8      3.000000    lm_tr  ade Kirde   Ranniku
## 3         9      3.169925    tr_lm peal Ranna   Ranniku
## 4         5      2.321928    tr_lm peal Kirde   Ranniku
##                                      kõneleja        enn
## 1                KJ1_VON_Ann_Plaks_EMH502.xml 0.90008208
## 2                 KJ_JOH_Liina_Laur_synt1.xml 0.48856865
## 3    KJ_HLJ_Johannes_Karineeme_synt-kontr.xml 0.22058898
## 4 KJ_JOH_Heinrich_Kaart_synt-kontrollitud.xml 0.07845498

LISA: Mudeli valideerimine

Eeskätt eksploratiivse mudeldamise käigus võib juhtuda, et meie mudel n-ö ülesobitab (ingl overfitting). See tähendab, et mudel sobitub küll meie andmetega, st andmetega, millel teda on treenitud, aga ei suuda uue valimi puhul enam kuigi palju seletada. Ülesobitamist võib juhtuda eriti näiteks siis, kui meil on vähe vaatlusi/ridu, aga palju seletavaid tunnuseid.

Ülesobitamise kontrollimiseks jagatakse piisavalt suure andmestiku olemasolul enamasti vaatlused kaheks: treeningandmestikuks ja testandmestikuks. Esimese peal ehitatakse mudelit nii, nagu meie seda oleme teinud, ja teise peal kontrollitakse, kas mudeli tulemused testandmetel on üldjoontes samad, mis treeningandmetel. Konfirmatoorse mudeldamise puhul põhimõtteliselt kontrollime juba olemasolevat teoreetilist mudelit enda andmetel, seega toimib selle strateegia puhul meie valim ainult testandmestikuna.

Alati ei ole andmestike ositamine aga võimalik, sest andmestikus on vähe vaatlusi, seal on kategoriaalseid seletavaid tunnuseid, millel on mitu taset ja mõnd taset esineb oluliselt harvem kui teisi vm. Sellisel juhul saab kasutada bootstrap-meetodit. See tähendab, et algsest andmestikust võetakse hulk juhuslikke valimeid, kuhu mõned algandmestiku vaatlused võivad sattuda mitu korda ja mõned hoopis välja jääda. Kui selliseid valimeid võtta piisavalt palju ja nende peal treenitud mudeleid testida, siis saab lõpuks kokku üsna usaldusväärsed üldistused.

Katsetame mõlemad variandid läbi, aga liigume tagantpoolt ettepoole ehk vaatame kõigepealt bootstrap-meetodiga valideerimist.

Bootstrap

Kuna meie interaktsiooniga mudel on praegu treenitud ja testitud samal andmestikul, siis valideerime mudelit bootstrap-meetodiga. Saame kasutada rms paketist funktsiooni lrm, et teha samasugune logistilise regressiooni mudel, nagu meil on juba glm-mudel, ning funktsiooni validate(), et oma mudelit valideerida. validate() funktsiooniga võetakse meie andmestikust hulk juhuslikke sama suuri valimeid, kus mõned algandmestiku read võivad korduda ja mõned read visatakse välja (tekitatakse n-ö juhuslikku varieeruvust). Nendel andmestikel treenitakse mudelid, mille üldistusi testitakse siis terve meie algandmestiku peal.

# install.packages("rms")
library(rms)
m5.lrm <- lrm(cx ~ LM_mobiilsus + LM_komplekssus + tegusõna_klass + LM_pikkus_log*murderühm, data = dat2, x = TRUE, y = TRUE)

# valideerimine võtab natuke aega
# B näitab meie andmestikust juhuslikult
# kordustega võetud valimite arvu
set.seed(1)
rms::validate(m5.lrm, B = 200)
##           index.orig training   test optimism index.corrected   n
## Dxy           0.4919   0.4977 0.4812   0.0165          0.4754 200
## R2            0.2406   0.2481 0.2321   0.0160          0.2246 200
## Intercept     0.0000   0.0000 0.0039  -0.0039          0.0039 200
## Slope         1.0000   1.0000 0.9553   0.0447          0.9553 200
## Emax          0.0000   0.0000 0.0108   0.0108          0.0108 200
## D             0.1983   0.2054 0.1906   0.0148          0.1835 200
## U            -0.0014  -0.0014 0.0003  -0.0016          0.0003 200
## Q             0.1997   0.2068 0.1903   0.0165          0.1832 200
## B             0.2044   0.2027 0.2063  -0.0036          0.2079 200
## g             1.1613   1.1888 1.1328   0.0560          1.1053 200
## gp            0.2467   0.2501 0.2419   0.0082          0.2385 200

Oluline on tulp optimism, mis näitab, kui palju on mudel treeningandmestikel ridade regressioonikordajaid (koefitsiente) üle hinnanud. Meid huvitavad selles tulbas read, kus on mudeli headuse näitajad Dxy (seotud C-indeksiga ehk 2 * (C-0.5)), R2, ja Slope, mis peegeldab seletavate tunnuste koefitsientide hinnanguid. optimism on tulpade training ja test vahe. training tulbas on iga rea statistiku väärtus, mis on üldistatud üle kõikide bootstrap-meetodiga tehtud mudelite. test tulba statistiku väärtus on saadud siis, kui bootstrap-meetoditega saadud mudeleid on testitud terve originaalandmestiku peal. Mida lähemal treening- ja testandmete tulemused üksteisele on, seda väiksem on optimism ja seda kindlamad võime olla selles, et meie mudel on üldistatav ka väljapoole meie valimit ja seega teaduslikult relevantne. Kui nimetatud statistikute puhul jäävad optimism-tulba väärtused alla 0.05, siis peaks kõik korras olema. Loe lähemalt nt siit.

Treening- ja testandmestik

Kirjutamata reegel on, et treening- ja testandmete suhe võiks olla u 7/3 ehk 70% originaalandmestikust võiks minna treeningandmestikuks ning ülejäänud 30% testandmestikuks. Juhusliku valimi võtmiseks saab kasutada funktsiooni sample() (selleks, et saada samasugust juhuvalimit ka teistel kordadel koodi jooksutamisel, määrame ka seemne (seed)).

# 70% vaatluste arvust on
round(nrow(dat2)*0.7)
## [1] 1011
# 1011 juhuslikku ridade indeksit andmestikust
set.seed(1)
trenn_index <- sample(1:nrow(dat2), 
                      round(nrow(dat2)*0.7), 
                      replace = FALSE)

# treeningandmestik
trenn <- dat2[trenn_index,]

# testandmestik
test <- dat2[-trenn_index,]

Teeme treeningandmetel mudeli.

m_trenn.glm <- glm(cx ~ LM_mobiilsus + LM_komplekssus + tegusõna_klass + LM_pikkus_log * murderühm, data = trenn, family = "binomial")

car::Anova(m_trenn.glm, type = "III")
## Analysis of Deviance Table (Type III tests)
## 
## Response: cx
##                         LR Chisq Df Pr(>Chisq)    
## LM_mobiilsus              79.884  1  < 2.2e-16 ***
## LM_komplekssus            25.177  1  5.231e-07 ***
## tegusõna_klass            27.653  4  1.466e-05 ***
## LM_pikkus_log              0.517  1   0.472328    
## murderühm                 12.136  3   0.006932 ** 
## LM_pikkus_log:murderühm    8.964  3   0.029770 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Arvutame nüüd treenitud mudeli headuse näitajad testandmestiku peal. C-indeksi (AUC) jaoks peame nüüd kasutama funktsioone paketist ROCR või Hmisc, sest me ei saa vaadata treeningmudeli enda C-d, vaid peame vaatama testandmestiku ennustuste pealt arvutatud näitajat.

# Klassiftseerimistäpsus
enn <- predict(m_trenn.glm, newdata = test)
teg <- as.character(test[,"cx"])
enn_bin <- ifelse(enn > 0.5, "peal", "ade")
(klass_acc <- sum(diag(table(teg, enn_bin)))/sum(table(teg, enn_bin)))
## [1] 0.6258661
# C-indeks/AUC
somers2(enn, as.numeric(test$cx)-1)["C"]
##         C 
## 0.7150214

Mudeli raporteerimine

Mudeli esitamisel on hea tava esitada koefitsientide tabel, mõni mudeli headuse näitaja (nt AUC) ning erinevus nullmudeli ja seletavate tunnustega mudeli hälbimuse (deviance) vahel.

tab_model(m5.glm, 
          transform = NULL, 
          show.ci = F, 
          show.se = T, 
          show.r2 = F, 
          show.loglik = T, 
          show.stat = T, 
          show.dev = T, 
          show.reflvl = T, 
          p.adjust = "holm", 
          title = "(AUC = 0.746, null deviance = 2001.809)", 
          string.stat = "z")
(AUC = 0.746, null deviance = 2001.809)
  cx
Predictors Log-Odds std. Error z p
(Intercept) 0.13 0.21 0.62 1.000
staatiline Reference
liigutatav 1.22 0.13 9.37 <0.001
lihtsõna Reference
liitsõna -1.21 0.21 -5.68 <0.001
LM_pikkus_log -0.14 0.13 -1.04 1.000
LM_pikkus_log:murderühmLõuna 0.61 0.23 2.64 0.058
LM_pikkus_log:murderühmRanniku 0.12 0.30 0.40 1.000
LM_pikkus_log:murderühmSaarte 0.65 0.23 2.77 0.050
olemisverb Reference
liikumisverb -0.60 0.18 -3.40 0.007
paiknemisverb 0.35 0.35 1.01 1.000
tegevusverb 0.39 0.14 2.76 0.050
verbita 0.27 0.21 1.28 0.996
Kesk Reference
Lõuna -0.49 0.32 -1.54 0.742
Ranniku -1.34 0.41 -3.25 0.012
Saarte -1.45 0.33 -4.35 <0.001
Observations 1444
Deviance 1714.431
log-Likelihood -857.216

Samuti võiks lisada mõned joonised seletavate tunnuste mõjude kohta.

Segamudelid

Nagu 13. praktikumis räägitud, siis siis segamudel arvestab sellega, et osa mõõtmisi võib tulla korduvatelt subjektidelt ning võimaldab mudelisse kaasata nii fikseeritud faktoreid kui ka juhuslikke faktoreid. 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 (enamasti) uurijale omaette huvi.

Tegeleme siin juhuslike vabaliikmetega (random intercept), mis korrigeeritakse iga juhusliku faktori taseme jaoks (nt iga indiviidi jaoks algaks koefitsientide võrdlemine natuke erinevalt tasemelt, sest igal kõnelejal on oma tõenäosuslik keelemudel, mille põhjal ta käändekonstruktsiooni ja peal-konstruktsiooni vahel valib).
Mudelisse võib aga lisada ka juhuslikke kaldeid (random slope), mis korrigeerivad valitud fikseeritud tunnuse efekti vastavalt iga juhusliku faktori tasemele (nt iga indiviidi kohta on mingi tunnuse mõju erinev, mõni kõneleja on näiteks tegusõna klassi mõju suhtes tundlikum kui mõni teine kõneleja).

Lisame oma regressioonimudelisse individuaalse kõneleja kui juhusliku faktori. Saame kasutada tunnust kõneleja.

library(lme4)
# mudeli tegemine võib võtta pisut aega
m1.glmer <- glmer(cx ~ LM_mobiilsus + LM_komplekssus + tegusõna_klass + LM_pikkus_log * murderühm + (1|kõneleja), data = dat2, family = "binomial")

Kui saame teate, et Model failed to converge, tähendab see seda, et mudeli algoritm ei suutnud mudeli parameetreid optimaalselt hinnata, mudeli tulemusi ei saa usaldada ja neid ei peaks nt publikatsioonides raporteerima. See võib juhtuda näiteks siis, kui mudel on andmete hulka arvestades liiga kompleksne või kui andmestikus on mingeid erindeid, mis parameetrite hindamist oluliselt raskendavad. Ent vahel võib olukorda parandada see, kui vahetada mudeli optimeerimise algoritmi. Sobiva leidmiseks võib kasutada näiteks funktsiooni allFit(). Funktsioon tuleb samuti paketist lme4, ent lisaks võidakse nõuda, et installitakse ka paketid optimx ja dfoptim.

#install.packages(c("optimx", "dfoptim"))
allFit(m1.glmer, verbose = F) -> allfits
summary(allfits)$which.OK
##                        bobyqa                   Nelder_Mead 
##                          TRUE                          TRUE 
##                    nlminbwrap                         nmkbw 
##                          TRUE                          TRUE 
##               optimx.L-BFGS-B nloptwrap.NLOPT_LN_NELDERMEAD 
##                          TRUE                          TRUE 
##     nloptwrap.NLOPT_LN_BOBYQA 
##                          TRUE
summary(allfits)$msgs
## $bobyqa
## NULL
## 
## $Nelder_Mead
## NULL
## 
## $nlminbwrap
## NULL
## 
## $nmkbw
## [1] "Model failed to converge with max|grad| = 0.00433066 (tol = 0.002, component 1)"
## 
## $`optimx.L-BFGS-B`
## [1] "Model failed to converge with max|grad| = 0.00335709 (tol = 0.002, component 1)"
## 
## $nloptwrap.NLOPT_LN_NELDERMEAD
## [1] "Model failed to converge with max|grad| = 0.0042745 (tol = 0.002, component 1)"
## 
## $nloptwrap.NLOPT_LN_BOBYQA
## [1] "Model failed to converge with max|grad| = 0.0111781 (tol = 0.002, component 1)"

Sagedasti kasutatakse nt bobyqa (Bound Optimization BY Quadratic Approximation) optimeerijat.

m1.glmer <- glmer(cx ~ LM_mobiilsus + LM_komplekssus + tegusõna_klass + LM_pikkus_log * murderühm + (1|kõneleja), 
                  data = dat2, family = "binomial", 
                  control = glmerControl(optimizer= "bobyqa"))
summary(m1.glmer)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## cx ~ LM_mobiilsus + LM_komplekssus + tegusõna_klass + LM_pikkus_log *  
##     murderühm + (1 | kõneleja)
##    Data: dat2
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   1718.3   1797.4   -844.1   1688.3     1429 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7744 -0.7213  0.0333  0.7220  3.7499 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  kõneleja (Intercept) 0.5285   0.7269  
## Number of obs: 1444, groups:  kõneleja, 272
## 
## Fixed effects:
##                                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                     0.1815362  0.2373131   0.765  0.44429    
## LM_mobiilsusliigutatav          1.3834638  0.1502642   9.207  < 2e-16 ***
## LM_komplekssusliitsõna         -1.3640095  0.2342015  -5.824 5.74e-09 ***
## tegusõna_klassliikumisverb     -0.6286121  0.1911790  -3.288  0.00101 ** 
## tegusõna_klasspaiknemisverb     0.5416890  0.3779737   1.433  0.15182    
## tegusõna_klasstegevusverb       0.3923725  0.1545602   2.539  0.01113 *  
## tegusõna_klassverbita           0.3389669  0.2331886   1.454  0.14605    
## LM_pikkus_log                  -0.1359389  0.1430258  -0.950  0.34188    
## murderühmLõuna                 -0.5412748  0.3669294  -1.475  0.14017    
## murderühmRanniku               -1.4194984  0.4790704  -2.963  0.00305 ** 
## murderühmSaarte                -1.5983041  0.3856839  -4.144 3.41e-05 ***
## LM_pikkus_log:murderühmLõuna    0.6238778  0.2539599   2.457  0.01403 *  
## LM_pikkus_log:murderühmRanniku -0.0007955  0.3266860  -0.002  0.99806    
## LM_pikkus_log:murderühmSaarte   0.6183636  0.2529755   2.444  0.01451 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 14 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

Mudeli väljundist näeme juba paljusid tuttavaid asju. Lisaks on seal tabel juhuslike mõjude kohta (Random effects), mis näitab siin vabaliikme hinnangulist varieerumist juhusliku faktori väärtuste vahel (dispersiooni ja standardhälvet). Mida väiksem varieeruvus on, seda sarnasemad on meie juhusliku mõju väärtused üksteisega (antud juhul kõnelejad) ning seda väiksem mõju juhuslikul faktoril mudelis on.

Mudeli tõlgendus on sarnane tavalise logistilise regressiooni mudeli tõlgendamisega, ent koefitsientide tõlgendamisel tuleb nüüd lisaks teiste juhuslike mõjude fikseerituna (st nende baastasemetel või nullis) hoidmisele hoida ka juhuslik mõju fikseerituna. See tähendab, et näiteks LM_mobiilsusliigutatav koefitsient näitab šansside logaritmitud muutust JUHUL, KUI kohasõna on lihtsõna, tegusõna väljendab olemist, kõneleja on Kesk-murderühma kõneleja, kohafraasi pikkus on 1 silp (log2(1) = 0) JA tegu on ühe ja sama kõnelejaga või teise kõnelejaga, kellel on täpselt samasugune juhuslik mõju. Aga kui mitte konkreetsetest šansisuhetest rääkida, siis võib tõlgendamisse suhtuda veidi vabamalt.

# Võrdleme mudeleid
tab_model(m5.glm, m1.glmer, show.aic = T)
  cx cx
Predictors Odds Ratios CI p Odds Ratios CI p
(Intercept) 1.14 0.76 – 1.71 0.532 1.20 0.75 – 1.91 0.444
LM mobiilsus [liigutatav] 3.37 2.62 – 4.36 <0.001 3.99 2.97 – 5.35 <0.001
LM komplekssus [liitsõna] 0.30 0.19 – 0.45 <0.001 0.26 0.16 – 0.40 <0.001
tegusõna klass
[liikumisverb] 		0.55 0.39 – 0.77 0.001 0.53 0.37 – 0.78 0.001
tegusõna klass
[paiknemisverb] 		1.42 0.72 – 2.83 0.313 1.72 0.82 – 3.61 0.152
tegusõna klass
[tegevusverb] 		1.47 1.12 – 1.94 0.006 1.48 1.09 – 2.00 0.011
tegusõna klass [verbita] 1.31 0.87 – 2.00 0.199 1.40 0.89 – 2.22 0.146
LM pikkus log 0.87 0.67 – 1.13 0.296 0.87 0.66 – 1.16 0.342
murderühm [Lõuna] 0.61 0.33 – 1.14 0.124 0.58 0.28 – 1.19 0.140
murderühm [Ranniku] 0.26 0.12 – 0.58 0.001 0.24 0.09 – 0.62 0.003
murderühm [Saarte] 0.23 0.12 – 0.45 <0.001 0.20 0.09 – 0.43 <0.001
LM pikkus log × murderühm
[Lõuna] 		1.85 1.17 – 2.92 0.008 1.87 1.13 – 3.07 0.014
LM pikkus log × murderühm
[Ranniku] 		1.13 0.63 – 2.02 0.688 1.00 0.53 – 1.90 0.998
LM pikkus log × murderühm
[Saarte] 		1.91 1.22 – 3.05 0.006 1.86 1.13 – 3.05 0.015
Random Effects
σ2   3.29
τ00   0.53 kõneleja
ICC   0.14
N   272 kõneleja
Observations 1444 1444
R2 Tjur 0.183 0.258 / 0.361
AIC 1742.431 1718.271

Tundub, et AIC läheb väiksemaks ning parameetrite hinnangutes ja statistilises olulises erilisi muutuseid ei toimu. Võib-olla võiks kaaluda kohafraasi pikkuse ja murderühma interaktsiooni väljajätmist, kuna seal kaoksid pärast Bonferroni või Holmi-Bonferroni parandusi statistiliselt olulised erinevused. Jätame selle praegu siiski mudelisse alles.

Vaatame ka mudeli headuse näitajaid.

Klassifitseerimistäpsus

enn <- fitted(m1.glmer)
teg <- as.character(dat2[,"cx"])
enn_bin <- ifelse(enn > 0.5, "peal", "ade")
klass_tabel <- table(teg, enn_bin)
(sum(diag(klass_tabel))/sum(klass_tabel))
## [1] 0.7472299

C-indeks (AUC)

library(Hmisc)
somers2(enn, as.numeric(dat2$cx)-1)["C"]
##         C 
## 0.8261562

Mudeli headuse näitajad on läinud oluliselt paremaks!

report::report_text(m1.glmer)
## We fitted a logistic mixed model (estimated using ML and BOBYQA optimizer) to
## predict cx with LM_mobiilsus, LM_komplekssus, tegusõna_klass, LM_pikkus_log and
## murderühm (formula: cx ~ LM_mobiilsus + LM_komplekssus + tegusõna_klass +
## LM_pikkus_log * murderühm). The model included kõneleja as random effect
## (formula: ~1 | kõneleja). The model's total explanatory power is substantial
## (conditional R2 = 0.36) and the part related to the fixed effects alone
## (marginal R2) is of 0.26. The model's intercept, corresponding to LM_mobiilsus
## = staatiline, LM_komplekssus = lihtsõna, tegusõna_klass = olemisverb,
## LM_pikkus_log = 0 and murderühm = Kesk, is at 0.18 (95% CI [-0.28, 0.65], p =
## 0.444). Within this model:
## 
##   - The effect of LM mobiilsus [liigutatav] is statistically significant and
## positive (beta = 1.38, 95% CI [1.09, 1.68], p < .001; Std. beta = 1.38, 95% CI
## [1.09, 1.68])
##   - The effect of LM komplekssus [liitsõna] is statistically significant and
## negative (beta = -1.36, 95% CI [-1.82, -0.90], p < .001; Std. beta = -1.36, 95%
## CI [-1.82, -0.90])
##   - The effect of tegusõna klass [liikumisverb] is statistically significant and
## negative (beta = -0.63, 95% CI [-1.00, -0.25], p = 0.001; Std. beta = -0.63,
## 95% CI [-1.00, -0.25])
##   - The effect of tegusõna klass [paiknemisverb] is statistically non-significant
## and positive (beta = 0.54, 95% CI [-0.20, 1.28], p = 0.152; Std. beta = 0.54,
## 95% CI [-0.20, 1.28])
##   - The effect of tegusõna klass [tegevusverb] is statistically significant and
## positive (beta = 0.39, 95% CI [0.09, 0.70], p = 0.011; Std. beta = 0.39, 95% CI
## [0.09, 0.70])
##   - The effect of tegusõna klass [verbita] is statistically non-significant and
## positive (beta = 0.34, 95% CI [-0.12, 0.80], p = 0.146; Std. beta = 0.34, 95%
## CI [-0.12, 0.80])
##   - The effect of LM pikkus log is statistically non-significant and negative
## (beta = -0.14, 95% CI [-0.42, 0.14], p = 0.342; Std. beta = -0.09, 95% CI
## [-0.27, 0.09])
##   - The effect of murderühm [Lõuna] is statistically non-significant and negative
## (beta = -0.54, 95% CI [-1.26, 0.18], p = 0.140; Std. beta = 0.24, 95% CI
## [-0.17, 0.65])
##   - The effect of murderühm [Ranniku] is statistically significant and negative
## (beta = -1.42, 95% CI [-2.36, -0.48], p = 0.003; Std. beta = -1.42, 95% CI
## [-1.94, -0.91])
##   - The effect of murderühm [Saarte] is statistically significant and negative
## (beta = -1.60, 95% CI [-2.35, -0.84], p < .001; Std. beta = -0.82, 95% CI
## [-1.24, -0.40])
##   - The effect of LM pikkus log × murderühm [Lõuna] is statistically significant
## and positive (beta = 0.62, 95% CI [0.13, 1.12], p = 0.014; Std. beta = 0.41,
## 95% CI [0.08, 0.73])
##   - The effect of LM pikkus log × murderühm [Ranniku] is statistically
## non-significant and negative (beta = -7.96e-04, 95% CI [-0.64, 0.64], p =
## 0.998; Std. beta = -5.21e-04, 95% CI [-0.42, 0.42])
##   - The effect of LM pikkus log × murderühm [Saarte] is statistically significant
## and positive (beta = 0.62, 95% CI [0.12, 1.11], p = 0.015; Std. beta = 0.40,
## 95% CI [0.08, 0.73])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald z-distribution approximation.

Ka siin on nimetatud mudeli headuse näitajatena pseudo-R-ruutusid, mille arvutamiseks leiab valemi nt siit. Tingimuslik (conditional) R-ruut arvestab seletatud varieerumise hulga puhul nii fikseeritud kui ka juhuslikke mõjusid, tavaline R-ruut ainult fikseeritud mõjusid.

Vaatame ka juhuslike mõjude jaotumist funktsiooniga ranef(). Tegemist on juhuslike vabaliikmete kohandustega ehk vahega mudeli üldistatud vabaliikme ning iga indiviidi jaoks ennustatud vabaliikme vahel. Teisisõnu näitavad need väärtused, kui palju mingi kõneleja peal-konstruktsiooni kasutamises keskmisest peal-kasutamisest erineb.

juh <- ranef(m1.glmer, condVar = TRUE)
head(juh$kõneleja)
##                                         (Intercept)
## KJ_AKS_Jaan_Koort_synt-kontrollitud.xml -0.15309378
## KJ_AMB_Marie_Raude_synt.xml             -0.01684003
## KJ_ANS_Ann_Ratas_F104.xml                0.33500414
## KJ_ANS_Ann_Ratas_F11.xml                -0.60454228
## KJ_ANS_Ann_Usin_F14_synt.xml            -0.51187475
## KJ_ANS_Ann_Usin_F19_synt.xml            -0.27189041
mean(juh$kõneleja$`(Intercept)`) # keskmine on alati 0 lähedal
## [1] -0.001099424
median(juh$kõneleja$`(Intercept)`) # mediaan on alati 0 lähedal
## [1] 0.01700568
hist(juh$kõneleja$`(Intercept)`)

Võime teha vabaliikmete kohanduste jaotumise joonise.

plot_model(m1.glmer, type = "re", sort.est = T, grid = F) # võib tahta, et installiksite enne paketi "glmmTMB"

Või ka ggplotiga.

# Teeme juhuslike mõjude objektist andmetabeli
data.frame(juh) %>%
  ggplot(aes(y = grp, x = condval)) + # y-teljel kõneleja kood, x-teljel kohandus vabaliikmes
  geom_point(alpha = 0.3) +
  geom_errorbar(aes(xmin = condval - 2*condsd, xmax = condval + 2*condsd), alpha = 0.3) + # lisa usaldusvahemikud
  geom_vline(xintercept = 0, color = "red", linetype = "dashed") + # tõmba 0-punkti punane joon
  labs(x = "Kohandus vabaliikmes",
       y = "Informandid")

Vaatame ainult kõnelejaid, kelle puhul usaldusvahemik ei kata nulli.

data.frame(juh) %>%
  filter(condval - 2*condsd > 0 | condval + 2*condsd < 0) %>%
  ggplot(aes(y = grp, x = condval)) +
  geom_point(alpha = 0.3) +
  geom_errorbar(aes(xmin = condval - 2*condsd, xmax = condval + 2*condsd), alpha = 0.3) +
  geom_vline(xintercept = 0, color = "red", linetype = "dashed") # tõmba 0-punkti punane joon

Vaatame näiteks vabaliikme kohandusi ka kõneleja murderühma järgi.

data.frame(juh) %>%
  merge(., dat2[,c("kõneleja", "murderühm")], by.x = "grp", by.y = "kõneleja") %>%
  ggplot(aes(y = grp, x = condval)) +
  facet_grid("murderühm", scales = "free", space = "free") +
  geom_point(alpha = 0.3) +
  geom_errorbar(aes(xmin = condval -2*condsd, xmax = condval + 2*condsd), alpha = 0.3) +
  geom_vline(xintercept = 0, color = "red", linetype = "dashed")

Võime hinnata ka seda, kui suure osa juhuslik mõju mudelis uuritava tunnuse varieerumisest ära seletab. Selleks kasutame funktsiooni icc() paketist performance.

performance::icc(m1.glmer)
## # Intraclass Correlation Coefficient
## 
##     Adjusted ICC: 0.138
##   Unadjusted ICC: 0.103

ICC väärtus ulatub 0st 1ni. Kui juhuslik mõju ei toimi andmetes üldse mingi grupeeriva tunnusena, siis on indeksi väärtus 0, ja kui juhuslik mõju kirjeldab ära kogu uuritava tunnuse varieerumise, on indeksi väärtus 1. Kui gruppe ehk juhusliku faktori tasemeid on vähe, võib tavaline ICC juhusliku faktori mõju ülehinnata. Kui gruppe on palju (nagu meil), võib tavaline ICC juhusliku faktori mõju jällegi alahinnata. Kohandatud ICC (Adjusted ICC) võtab arvesse aga ka seda, kui palju gruppe ehk juhusliku faktori tasemeid on, ning on seega üldjuhul täpsem mõõdik. Juhusliku faktori poolt kirjeldatud osa varieerumisest on siin ~0.14 ehk 14%.

Segamudelitele kehtivad samad eeldused, mis ainult fikseeritud mõjudega mudelitele, välja arvatud vaatluste sõltumatuse nõue (kui juhuslik faktor grupeerib mingeid vaatlusi). Seega tuleks ka segamudeli puhul hinnata šansi logaritmi ja arvuliste tunnuste vahelise suhte lineaarsust ning vaadata, kas mudeli tunnuste vahel esineb multikollineaarsust.

data.frame(logit = predict(m1.glmer, type = "link"), 
           pikkus = dat2$LM_pikkus_log, 
           ryhm = dat2$murderühm) %>%
  ggplot(aes(x = pikkus, y = logit)) +
  geom_point(alpha = 0.1) +
  facet_wrap("ryhm") +
  geom_smooth(method = "loess")

car::vif(m1.glmer)
##                              GVIF Df GVIF^(1/(2*Df))
## LM_mobiilsus             1.043260  1        1.021401
## LM_komplekssus           1.083210  1        1.040774
## tegusõna_klass           1.034016  4        1.004190
## LM_pikkus_log            2.209743  1        1.486521
## murderühm               35.936626  3        1.816587
## LM_pikkus_log:murderühm 45.804045  3        1.891548

Siin võiksime vaadata parempoolses tulbas olevaid, standardiseeritud ja võrreldavaid GVIF-skoore. Need on < 5, seega ei pea tunnuste omavahelise seotuse pärast muretsema.

Juhuslikke kaldeid tunnuse kõneleja põhjalt on keeruline lisada, kuna meil on tunnusel kõneleja väga palju tasemeid, mudelis mitu mitme tasemega fikseeritud tunnust ning seega esineb väga harvu tunnusekombinatsioone.

Küll aga võiksime proovida lisada mudelisse juhusliku vabaliikmena ka kohasõna lemma LM_lemma, kuna käände- või peal-konstruktsiooni kasutus võib sõltuda osaliselt ka sellest, mis sõnad on meie andmestikku sattunud.

m2.glmer <- glmer(cx ~ LM_mobiilsus + LM_komplekssus + tegusõna_klass + LM_pikkus_log * murderühm + (1|kõneleja) + (1|LM_lemma), data = dat2, family = "binomial", control = glmerControl(optimizer = "bobyqa"))
summary(m2.glmer)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## cx ~ LM_mobiilsus + LM_komplekssus + tegusõna_klass + LM_pikkus_log *  
##     murderühm + (1 | kõneleja) + (1 | LM_lemma)
##    Data: dat2
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   1521.5   1605.9   -744.7   1489.5     1428 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4666 -0.4747  0.0122  0.4267  5.2983 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  LM_lemma (Intercept) 3.9322   1.9830  
##  kõneleja (Intercept) 0.8606   0.9277  
## Number of obs: 1444, groups:  LM_lemma, 338; kõneleja, 272
## 
## Fixed effects:
##                                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                      0.7615     0.4038   1.886 0.059303 .  
## LM_mobiilsusliigutatav           1.2021     0.3330   3.610 0.000306 ***
## LM_komplekssusliitsõna          -2.0527     0.4880  -4.207 2.59e-05 ***
## tegusõna_klassliikumisverb      -0.5956     0.2585  -2.304 0.021232 *  
## tegusõna_klasspaiknemisverb      0.7024     0.5247   1.339 0.180682    
## tegusõna_klasstegevusverb        0.5003     0.2060   2.428 0.015172 *  
## tegusõna_klassverbita            0.4245     0.3060   1.387 0.165412    
## LM_pikkus_log                    0.1396     0.2002   0.697 0.485510    
## murderühmLõuna                  -1.0224     0.4714  -2.169 0.030083 *  
## murderühmRanniku                -1.7099     0.6063  -2.820 0.004797 ** 
## murderühmSaarte                 -1.8558     0.4894  -3.792 0.000149 ***
## LM_pikkus_log:murderühmLõuna     0.8059     0.3310   2.434 0.014919 *  
## LM_pikkus_log:murderühmRanniku  -0.2263     0.4130  -0.548 0.583729    
## LM_pikkus_log:murderühmSaarte    0.7240     0.3258   2.222 0.026281 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 14 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
anova(m1.glmer, m2.glmer, test = "Chisq")
## Data: dat2
## Models:
## m1.glmer: cx ~ LM_mobiilsus + LM_komplekssus + tegusõna_klass + LM_pikkus_log * murderühm + (1 | kõneleja)
## m2.glmer: cx ~ LM_mobiilsus + LM_komplekssus + tegusõna_klass + LM_pikkus_log * murderühm + (1 | kõneleja) + (1 | LM_lemma)
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## m1.glmer   15 1718.3 1797.4 -844.14   1688.3                         
## m2.glmer   16 1521.5 1605.9 -744.73   1489.5 198.81  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Saame parema mudeli küll, ehkki ka jääkide ulatus läheb suuremaks!

Kordamisküsimused

1. Mida näitab klassifitseerimistäpsus?

  1. Klassifitseerimistäpsus näitab nende kordade arvu, mil mudel ennustab tegelikule vaatlusele õige klassi.
  2. Klassifitseerimistäpsus näitab nende kordade osakaalu, mil mudel ennustab konkreetsele vaatlusele õige klassi.
  3. Klassifitseerimistäpsus näitab, kui oluline mudel on.
  4. Klassifitseerimistäpsus näitab, kui suur osa ennustatud tõenäosustest on > 0.5.

2. Mille poolest erineb C-indeks (AUC) klassifitseerimistäpsusest?

  1. Esitatakse alati katkendliku punase joonega.
  2. Võib saada ka 1st suuremaid väärtusi.
  3. On paindlikum mudeli headuse näitaja.
  4. Hindab mudeli võimet eristada uuritava tunnuse klasse tõenäosuste kaudu.

3. Milliseid eeldusi tuleb ainult fikseeritud mõjudega logistilise regressiooni puhul kontrollida?

  1. Seda, et ühelt subjektilt ei oleks andmestikus mitut vaatlust.
  2. Seda, et uuritav tunnus oleks normaaljaotusega.
  3. Seda, et arvulised seletavad tunnused oleksid normaaljaotusega.
  4. Seda, et mudeli ennustatud šansi logaritmi ja arvulise seletava tunnuse vaheline suhe oleks lineaarne.
  5. Seda, et seletavad tunnused ei seletaks uuritava tunnuse varieerumisest sama asja.
  6. Enne mudeldamist hii-ruut-testiga seda, et kõik mudelisse kaasatavad kategoriaalsed tunnused ikka on uuritava tunnusega oluliselt seotud.

4. Mida tähendab bootstrapping statistikas?

  1. Sellise mudeli tegemist, kus andmepunktidele ennustatud väärtuste jaotus on saapakujuline.
  2. See on valikumeetod, kus üldkogumi ehk populatsiooni mingi parameetri hindamiseks jagatakse andmestik käsitsi test- ja treeningandmeteks, tehakse kummagi peal eraldi mudeleid ja võrreldakse siis mudelite näitajaid.
  3. See on valikumeetod, kus üldkogumi ehk populatsiooni mingi parameetri hindamiseks võetakse olemasolevast andmestikust hulk sama suuri juhuslikke valimeid, milles osa vaatlusi võib korduda ja osa vaatlusi jäetakse välja.

Järgmiseks kasutame ära R-i n-ö sisseehitatud andmestikke ning teeme logistilise regressiooni mudeli andmestikust, mis sisaldab andmeid Torontos marihuaana omamise eest arreteeritud isikute kohta (vt täpsemalt siit). Tunnuste väärtused on sellised:

  • Released: kas vahistatu lasti kohtukutse vastu vabaks või mitte (= vangistati).
    • Tasemed: Yes, No.
  • Colour: vahistatu nahavärv.
    • Tasemed: Black, White.
  • Year: vahistamise aasta.
    • Väärtused: 1997-2002
  • Age: vahistatu vanus vahistamise ajal
    • Väärtused: 12-66
  • Sex: vahistatu sugu.
    • Tasemed: Male, Female.
  • Employed: kas vahistatul on töökoht või mitte.
    • Tasemed: Yes, No.
  • Citizen: kas vahistatu on Kanada kodanik või mitte.
    • Tasemed: Yes, No.
  • Checks: kui mitu korda vahistatu nimi varasemalt politsei andmebaasis esines.
    • Väärtused: 0-6.
# install.packages("carData")
# Laadime andmestiku
Arrests <- carData::Arrests

5. Ennusta seda, kas arreteeritu lasti kohtukutsega vabadusse või vangistati (tunnus released, vastavalt tasemed Yes ja No). Lisa mudelisse fikseeritud peamõjudena nahavärv, vanus ja sugu. Milline tunnustest ei tee mudelit oluliselt paremaks?

m1.glm <- glm(released ~ colour + age + sex, data = Arrests, family = "binomial")
  1. nahavärv
  2. vanus
  3. sugu
  4. ükski ei tee

6. Milline nende kolme tunnuse interaktsioonist teeb mudelit paremaks?

  1. colour*age
  2. colour*sex
  3. age*sex
  4. ükski ei tee

7. Tee nüüd allolev mudel. Mis on selle mudeli klassifitseerimistäpsus?

m2.glm <- glm(released ~ colour + age + employed + citizen + checks + colour*age, data = Arrests, family = "binomial")
  1. 0.01
  2. 0.19
  3. 0.82
  4. 0.83

8. Mis on selle mudeli C/AUC?

  1. pole võimalik arvutada
  2. 0.73
  3. 0.83

9. Kas mudelis esineb multikollineaarsust?

  1. Ei.
  2. Jah, mingil määral, aga ainult interaktsioonis osalevate tunnuste vahel.
  3. Jah, olulisel määral.

Järgmisel korral

  • Tingimuslikud otsustuspuud

Funktsioonide sõnastik

  • car::influencePlot() - tuvasta andmestikust mudeli seisukohast ebatavalised vaatlused
  • car::vif() - kontrolli seletavate tunnuste vahelist multikollineaarsust
  • Hmisc::somers2() - leia binaarse kategoriaalse uuritava tunnusega mudeli C-indeks/AUC
  • lme4::allFit() - leia sobiv segamudeli optimeerimise algoritm
  • lme4::glmer() - tee logistilise regressiooni segamudel
  • lme4::ranef() - küsi segamudelist välja juhuslikud mõjud
  • performance::check_collinearity() - kontrolli seletavate tunnuste vahelist multikollineaarsust
  • performance::icc() - leia juhuslike mõjude rühmasisese korrelatsiooni koefitsient
  • pROC::auc() - leia binaarse kategoriaalse uuritava tunnusega mudeli C-indeks/AUC
  • rms::lrm() - tee logistilise regressiooni mudel, mis annab väljundis rohkem infot kui glm()
  • rms::validate() - valideeri logistilise regressiooni mudelit bootstrap-meetodiga
  • sapply() - rakenda igale listi liikmele (nt igale andmetabeli tulbale) mingit funktsiooni ja tagasta vastused vektorina