Tänases praktikumis

Eelmisel korral

Mitme tunnusega regressioonimudel
Mudelite võrdlemine ja optimaalse mudeli leidmine
Interaktsioonide lisamine

Sel korral

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

Mudelite headuse näitajad

Eelmisel korral lõpetasime esimese isiku suure andmestiku peal tehtud mudeliga Pron ~ Lopp + Aeg + Sonepikkus + Ref_kaugus_num + Verbiklass + Murre + Aeg:Verbiklass.

# Laadime andmestiku
load("isik1suur.RData")

# Kustutame ebavajalikud sümbolid.
# Selleks asendame need mitte millegagi.
isik1suur$Sone <- gsub("['#()\\*\\.`]", "", isik1suur$Sone)

# Arvutame tegusõnade pikkuse tähemärkides
isik1suur$Sonepikkus <- nchar(isik1suur$Sone) 

# Teeme mudeli
m7.glm <- glm(Pron ~ Lopp + Aeg + Sonepikkus + Ref_kaugus_num + Verbiklass + Murre + Verbiklass:Aeg, data = isik1suur, family = "binomial")

Vaatame veel kord seletavate tunnuste mõju asesõna väljendamisele.

summary(m7.glm)
## 
## Call:
## glm(formula = Pron ~ Lopp + Aeg + Sonepikkus + Ref_kaugus_num + 
##     Verbiklass + Murre + Verbiklass:Aeg, family = "binomial", 
##     data = isik1suur)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4790  -1.0172   0.4542   1.0166   2.0954  
## 
## Coefficients:
##                                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                      3.010527   0.326947   9.208  < 2e-16 ***
## Loppjah                         -1.469137   0.128783 -11.408  < 2e-16 ***
## Aegpr                            0.779829   0.287908   2.709  0.00676 ** 
## Sonepikkus                      -0.216372   0.025250  -8.569  < 2e-16 ***
## Ref_kaugus_num                   0.047638   0.005467   8.714  < 2e-16 ***
## Verbiklasskommunikatiivne        1.001958   0.313384   3.197  0.00139 ** 
## Verbiklassmuu                   -0.649403   0.242598  -2.677  0.00743 ** 
## MurreKESK                       -0.636810   0.121661  -5.234 1.66e-07 ***
## MurreRANNA                      -0.777364   0.139647  -5.567 2.60e-08 ***
## MurreSAARTE                     -0.690195   0.142992  -4.827 1.39e-06 ***
## MurreVORU                       -1.810241   0.180786 -10.013  < 2e-16 ***
## Aegpr:Verbiklasskommunikatiivne -1.064950   0.412732  -2.580  0.00987 ** 
## Aegpr:Verbiklassmuu             -0.441221   0.297074  -1.485  0.13748    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 5617.4  on 4065  degrees of freedom
## Residual deviance: 4893.1  on 4053  degrees of freedom
## AIC: 4919.1
## 
## Number of Fisher Scoring iterations: 4

library(effects)
plot(predictorEffects(m7.glm), type = "response", rows = 3, cols = 2)

Roosad vurrud on tõenäosushinnangute 95% usaldusvahemikud. Kui vurrud on väga laiad või kui hinnangute vurrud olulisel määral kattuvad, siis on üldjuhul koefitsiendi z-väärtus pigem nullile lähemal ja p-väärtus pigem suurem, sest kindlus, millega saame öelda, et kaks konteksti üksteisest oluliselt erinevad, on väiksem.

Sõnastame veel kord mudeli tulemused.

Pöördelõpu esinemine tegusõna vormil (nt olen) vähendab oluliselt asesõna mina/ma esinemise tõenäosust. Pöördelõpuga tegusõna vormiga koos asesõna pigem ei esine.
Kognitiivsete tegusõnadega (nt teadma) ja muude (mittekognitiivsete ja mittekommunikatiivsete, nt jooksma) tegusõnadega asesõna kasutamise tõenäosus olevikuvormides kasvab, kommunikatiivsete (nt ütlema) tegusõnadega kahaneb. Sealjuures on nii kognitiivsete kui ka kommunikatiivsete tegusõnadega asesõna esinemise tõenäosus mõlemas ajavormis suurem kui muude tegusõnadega. Kui muu tegusõna on minevikuvormis, siis sellega asesõna tõenäolisemalt koos ei esine.
Tegusõna tähemärkide arvu tõustes asesõna esinemise tõenäosus pigem kahaneb. Seega on tulemus vastupidine sellele, mida hüpoteesina esitasime. Näeme ka, et seos ei ole päris lineaarne (sellest mudeli eelduste kontrolli juures pikemalt).
Mida kaugemal on eelmine viide esimesele isikule, seda tõenäolisem on asesõna kasutus.
Aeg*Verbiklass näitab põhimõtteliselt sama infot, mis Verbiklass*Aeg, aga lihtsalt teisest perspektiivist.
Idamurde alal on asesõna kasutamise tõenäosus kõige kõrgem. Kõikides teistes murretes madalam. Eriti madal on asesõna kasutamise tõenäosus Võru murdes, kus asesõna tõenäolisemalt ei kasutata.

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.

Klassifitseerimistäpsus

Iga vaatluse kohta andmestikus ennustab mudel mingi tõenäosuse, millega asesõna mudelisse kaasatud seletavate tunnuse põhjal koos tegusõnaga esineda võiks. Kui meil on andmestikus rida, kus Pron = "jah" ja mudel ennustab sellele vaatlusele asesõna esinemise tõenäosuse, mille väärtus on > 0.5, siis ennustab mudel õigesti. Kui 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(m7.glm$fitted.values)
##         1         2         3         4         5         6 
## 0.6715597 0.7458081 0.5036226 0.3711365 0.3998122 0.2356891

# Võrdleme seda andmestiku tulba "Pron" tegelike väärtustega
head(isik1suur$Pron)
## [1] ei ei ei ei ei ei
## Levels: ei jah

# Ka tegelikud väärtused (kujul 0 ja 1) on tegelikult ka mudelis olemas
head(m7.glm$y)
## 1 2 3 4 5 6 
## 0 0 0 0 0 0

Näeme, et esimese 6 vaatluse puhul ennustab mudel 3 juhul õigesti (vaatlused 3-6) ja 3 juhul valesti (vaatlused 1-3), sest nendele ennustatud asesõna esinemise tõenäosus on suurem kui 0.5, ehkki tegelikult on andmestiku järgi seal kontekstis tegelikult “ei”.

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 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(m7.glm, newdata = isik1suur, type = "response"))
##         1         2         3         4         5         6 
## 0.6715597 0.7458081 0.5036226 0.3711365 0.3998122 0.2356891

# Ennustatud šansi logaritmid
head(predict(m7.glm, newdata = isik1suur, type = "link"))
##           1           2           3           4           5           6 
##  0.71524800  1.07637907  0.01449054 -0.52734409 -0.40624783 -1.17646133

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 Pron = "jah". 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 Pron = "ei". 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 “ei”) ja valenegatiivseteks (false negatives, FN - mudel ennustab < 0.5, tegelik vaatlus on “jah”).

Klassifitseerimistäpsuse leidmiseks on erinevaid viise. Näiteks

# vaatlustele ennustatud tõenäosused
enn <- m7.glm$fitted.values
head(enn)
##         1         2         3         4         5         6 
## 0.6715597 0.7458081 0.5036226 0.3711365 0.3998122 0.2356891

# "Pron" tulp andmestikus
teg <- as.character(isik1suur[,"Pron"])
head(teg)
## [1] "ei" "ei" "ei" "ei" "ei" "ei"

# tee uus objekt, milles on sama palju ridu
# kui 1. isiku andmestikus vaatlusi/ridu ning
# milles esineb "jah", 
# kui vastav väärtus enn-objektis on > 0.5, 
# ja "ei", kui vastav väärtus enn-objektis on < 0.5
enn_bin <- ifelse(enn > 0.5, "jah", "ei")
head(enn_bin)
##     1     2     3     4     5     6 
## "jah" "jah" "jah"  "ei"  "ei"  "ei"

# tee tegelikest ja ennustatud binaarsetest 
# väärtustest sagedustabel, kus on
# tõeliselt positiivsed, valepositiivsed,
# tõeliselt negatiivsed ja valenegatiivsed
# juhud
klass_tabel <- table(enn_bin, teg)
klass_tabel
##        teg
## enn_bin   ei  jah
##     ei  1301  741
##     jah  592 1432

# leia klassifitseerimistäpsus
# selleks liida kokku juhud, kus mudel ennustas "ei"
# ja ka tegelik vaatlus oli "ei", ning juhud, kus mudel
# ennustas "jah" ja ka tegelik vaatlus oli "jah" (= õigete
# ennustuste hulk), ning jaga tulemus vaatluste arvuga
(klass_acc <- klass_tabel["ei","ei"]+klass_tabel["jah","jah"])/sum(klass_tabel)
## [1] 0.6721594

Mudel suudab seega klassifitseerida õigesti 67% 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.

(t <- table(isik1suur$Pron))
## 
##   ei  jah 
## 1893 2173
max(t)
## [1] 2173

Sellisel juhul oleks meie klassifitseerimistäpsus

max(t)/sum(t)
## [1] 0.5344319

Seletavate tunnustega mudel ennustab seega ~14% 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 confusion matrix ja selle põhjal saab tegelikult arvutada veel mitmeid klassifitseerimisheaduse näitajaid, vt rohkem siit.

Näiteks võib veel joonistatada nn ROC-ruumi (ROC space) TPR-i (true positive rate) ja FPR-i (false positive rate) järgi. TPR on “jah” ennustuste osakaal kõikidest “jah” vaatlustest ning FPR on “jah” ennustuste osakaal kõikidest “ei” vaatlustest.

klass_tabel
##        teg
## enn_bin   ei  jah
##     ei  1301  741
##     jah  592 1432

# TPR ja FPR
TPR <- klass_tabel["jah", "jah"]/sum(klass_tabel[,"jah"])
FPR <- klass_tabel["jah", "ei"]/sum(klass_tabel[,"ei"])

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”/“jah” või “0”/“ei”, sest kui võtta mudelist välja kõik “jah” vaatustele ennustatud tõenäosused ning kõik “ei” 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

Selleks, et saada veelgi rohkem infot selle kohta, kui hästi mudel uuritava tunnuse tasemeid eristada suudab, teeme logistilise regressiooni mudeli funktsiooniga lrm() paketist rms. Selle tulemus on sama, mis glm-mudeliga, ent väljund annab veidi rohkem informatsiooni.

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

# siin ole vaja täpsustada family = "binomial"
m7.lrm <- lrm(Pron ~ Lopp + Aeg + Sonepikkus + Ref_kaugus_num + Verbiklass + Murre + Verbiklass:Aeg, data = isik1suur)

# ei ole vaja ka eraldi `summary()` käsku
m7.lrm
## Logistic Regression Model
##  
##  lrm(formula = Pron ~ Lopp + Aeg + Sonepikkus + Ref_kaugus_num + 
##      Verbiklass + Murre + Verbiklass:Aeg, data = isik1suur)
##  
##                        Model Likelihood     Discrimination    Rank Discrim.    
##                           Ratio Test           Indexes           Indexes       
##  Obs          4066    LR chi2     724.30    R2       0.218    C       0.732    
##   ei          1893    d.f.            12    g        1.100    Dxy     0.465    
##   jah         2173    Pr(> chi2) <0.0001    gr       3.004    gamma   0.465    
##  max |deriv| 3e-11                          gp       0.232    tau-a   0.231    
##                                             Brier    0.208                     
##  
##                                      Coef    S.E.   Wald Z Pr(>|Z|)
##  Intercept                            3.0105 0.3269   9.21 <0.0001 
##  Lopp=jah                            -1.4691 0.1288 -11.41 <0.0001 
##  Aeg=pr                               0.7798 0.2879   2.71 0.0068  
##  Sonepikkus                          -0.2164 0.0252  -8.57 <0.0001 
##  Ref_kaugus_num                       0.0476 0.0055   8.71 <0.0001 
##  Verbiklass=kommunikatiivne           1.0020 0.3134   3.20 0.0014  
##  Verbiklass=muu                      -0.6494 0.2426  -2.68 0.0074  
##  Murre=KESK                          -0.6368 0.1217  -5.23 <0.0001 
##  Murre=RANNA                         -0.7774 0.1396  -5.57 <0.0001 
##  Murre=SAARTE                        -0.6902 0.1430  -4.83 <0.0001 
##  Murre=VORU                          -1.8102 0.1808 -10.01 <0.0001 
##  Aeg=pr * Verbiklass=kommunikatiivne -1.0650 0.4127  -2.58 0.0099  
##  Aeg=pr * Verbiklass=muu             -0.4412 0.2971  -1.49 0.1375  
##

Alumises tabelis näeme sama infot, mida ka glm-mudeli väljundis. Ülemises tabelis on aga rida muid näitajaid, mis annavad meile mudeli headuse kohta kasulikku infot.

Neist sagedamini kasutatavad:

Model Likelihood Ratio Test näitab, kas mudel on statistiliselt oluline (raporteerib statistiku väärtuse, vabadusastmete arvu ja p-väärtuse). Nullhüpotees on, et mudeli hälbimus ehk seletamata variatsioon ei erine ühegi seletava tunnuseta mudeli ehk nullmudeli omast.
R2 on Nagelkerke pseudo-R2 ja viitab logistilises regressioonis mudeli headusele ainult umbkaudselt, aga ei näita, kui palju mudel variatiivsusest kirjeldab nagu lineaarses regressioonis!
C ehk konkordantsiindeks (ka Area Under the ROC Curve ehk AUC. NB! mitte ajada segi AIC-ga!) näitab logistilise mudeli ennustustäpsust: kui hästi suudab mudel uuritava tunnuse klasse eristada. Näitab 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”. C-d peetakse täpsemaks ja paindlikumaks mudeli eristusvõime hindajaks kui klassifitseerimistäpsust, kuna see arvestab ennustusi tõenäosustena ehk pideva muutujana, 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.
- 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
Dxy on 2*(C−0.5).

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

Saame vastavaid näitajaid küsida ka mudelist ilma mudelit ennast välja trükkimata.

# kõik näitajad
m7.lrm$stats
##          Obs    Max Deriv   Model L.R.         d.f.            P            C 
## 4.066000e+03 2.733103e-11 7.242972e+02 1.200000e+01 0.000000e+00 7.324160e-01 
##          Dxy        Gamma        Tau-a           R2        Brier            g 
## 4.648319e-01 4.654339e-01 2.313707e-01 2.179064e-01 2.080628e-01 1.099967e+00 
##           gr           gp 
## 3.004066e+00 2.320769e-01

# C-indeks
m7.lrm$stats["C"]
##        C 
## 0.732416

Kui mudelil on ideaalne eristusvõime ja see suudab täiuslikult eristada “ei”-sid ja “jah”-isid (kahe klassi tõenäosusjaotused ei kattu), siis on ROC-kurv 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:

enn_jah <- enn[which(isik1suur$Pron == "jah")]
enn_ei <- enn[which(isik1suur$Pron == "ei")]

plot(density(enn_jah, adjust = 2), col="darkgreen", xlab="Tõenäosused", ylab="", type="l", lwd=2, cex=2, xlim = c(0, 1), ylim = c(0,3), main = "Kahe klassi ennustatud\ntõenäosuste jaotus")
lines(density(enn_ei, adjust = 2), col="darkred", type="l", lwd=2, cex=2)
legend("topright", legend = c("ei", "jah"), fill = c("darkred", "darkgreen"))

Ja ROC-kurv nii:

library(ROCR)
## Loading required package: gplots
## 
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
## 
##     lowess
enn2 <- prediction(m7.glm$fitted.values, isik1suur$Pron)
perf <- performance(enn2, "tpr", "fpr")
plot(perf, col = "red", lty = "dashed")

Mudeli headuse näitajate kohta loe veel siit.

Mudeli eelduste kontroll

Ehkki logistiline regressioon ei eelda nt mingeid normaaljaotusi, on tegemist siiski parameetrilise mudeliga. Logistilise regressiooni eeldused on:

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

Vaatluste sõltumatus

Vaatluste sõltumatus tähendaks, et meil on nt ühe informandi kohta andmestikus üksainus vaatlus/rida, kus ta siis kas on kasutanud asesõna või mitte. Keeleteaduse andmestikud on harva sellised ning tihtipeale kasutatakse ära kõik vaatlused/lausungid, mida saada on. Meil on andmestikus tunnus Fail, mis näitab lindistuse failinime. Selles failinimes sisaldub ka lindistatud informandi nimi. Harvadel juhtudel osaleb ühes lindistuses ka mitu informanti, ent üldjuhul võib failinime pidada heaks indikaatoriks selle kohta, kui palju mingi informant vaatlusi on panustanud.

head(sort(table(isik1suur$Fail), decreasing = TRUE))
## 
##   PLV_Leena_Kurvits_EMH_0222_synt.xml           KEI_Mihkel_Anion_EMH787.xml 
##                                   118                                   108 
##            MUH_Raissa_Molder_synt.xml HLJ_Johannes_Karineeme_synt-kontr.xml 
##                                   103                                   100 
##   PEE_Maria_Limberg_EMH0123_lihts.xml MUS_Johannes_Lepp_EMH2904_1_lihts.xml 
##                                    91                                    90

Kõige “produktiivsematelt” keelejuhtidelt on andmestikus tervelt sadakond vaatlust, mistõttu on esimene regressiooni eeldus tegelikult rikutud.

Sellisel juhul oleks kohane kasutada segamõjudega mudelit. Räägime sellest veidi aja pärast.

Lineaarne suhe šansi logaritmi ja arvuliste tunnuste vahel

Meil on mudelis kaks arvulist tunnust: Sonapikkus ja Ref_kaugus_num. Lineaarsuse kontrollimiseks sobib visuaalne vaatlus.

Teeme kõigepealt tabeli, kus oleks kõrvuti arvuliste tunnuste tulbad ning mudeliga vastavale vaatlusele ennustatud asesõna esinemise šansi logaritm.

# vaatlustele ennustatud tõenäosused
enn <- m7.glm$fitted.values

# arvuliste tunnuste tulbad eraldi objektis
arvulised <- isik1suur[, c("Sonepikkus", "Ref_kaugus_num")]
head(arvulised)
##   Sonepikkus Ref_kaugus_num
## 1          5             19
## 2          7              1
## 3         10             27
## 4          7              2
## 5          6              0
## 6         10              2

# lisa "arvulised" tabelisse tulp "logit", 
# kus on tõenäosustest arvutatud šansi logaritmid
arvulised$logit <- log(enn/(1-enn))
# võiks teha ka nii
# arvulised$logit <- predict(m7.glm, isik1suur, type = "link")

head(arvulised)
##   Sonepikkus Ref_kaugus_num       logit
## 1          5             19  0.71524800
## 2          7              1  1.07637907
## 3         10             27  0.01449054
## 4          7              2 -0.52734409
## 5          6              0 -0.40624783
## 6         10              2 -1.17646133

Nüüd tuleb teha sellest nn laias formaadis tabelist pikk. See tähendab seda, et kui praegu on meil kaks arvulist tunnust eraldi tulpades ja need “jagavad” ühist logit-väärtust, siis pikas formaadis on meil üks nn võtme (key) tulp, mis sisaldab arvuliste tunnuste tulbanimesid (“Sonepikkus” ja “Ref_kaugus_num”), ning väärtuse (value) tulp, milles on nende tulpade sisu. Logit-väärtuseid korratakse kummagi võtmetunnuse kohta eraldi. Meie tabel läheb niisiis kaks korda pikemaks. Kasutame selleks funktsiooni gather() paketist tidyr (vt ka siia).

# tee laiast formaadist pikk
# install.packages("tidyr")
arvulised_pikk <- tidyr::gather(arvulised, key = "arv_tunnus", value = "arv_vaartus", -logit)

nrow(arvulised)
## [1] 4066
ncol(arvulised)
## [1] 3
nrow(arvulised_pikk)
## [1] 8132
ncol(arvulised_pikk)
## [1] 3
head(arvulised_pikk)
##         logit arv_tunnus arv_vaartus
## 1  0.71524800 Sonepikkus           5
## 2  1.07637907 Sonepikkus           7
## 3  0.01449054 Sonepikkus          10
## 4 -0.52734409 Sonepikkus           7
## 5 -0.40624783 Sonepikkus           6
## 6 -1.17646133 Sonepikkus          10
tail(arvulised_pikk)
##           logit     arv_tunnus arv_vaartus
## 8127 -0.2709709 Ref_kaugus_num          10
## 8128 -0.2927888 Ref_kaugus_num           5
## 8129 -0.1716926 Ref_kaugus_num           3
## 8130  0.7334308 Ref_kaugus_num          22
## 8131  0.2403738 Ref_kaugus_num           0
## 8132  0.6188885 Ref_kaugus_num           0

Lõpuks teeme joonise

library(ggplot2)
ggplot(data = arvulised_pikk, 
       aes(y = logit, x = arv_vaartus))+
  geom_point(size = 0.5, alpha = 0.5) +
  facet_wrap("arv_tunnus") +
  geom_smooth(method = "loess")

Näeme, et suhe viitekauguse ja šansi logaritmi vahel on enam-vähem lineaarne, aga tegusõna pikkuse puhul on lühemate sõnadega sõnapikkuse kasvades šansi logaritmi langus järsem kui keskmise pikkusega sõnadega ning päris pikki sõnu on tegelikult liiga vähe, et nende kohta midagi arvata, nii et ennustuste usaldusvahemik läheb väga laiaks.

Proovime mudelit logaritmitud tegusõna pikkusega.

isik1suur$Sonepikkus_log <- log(isik1suur$Sonepikkus)

# Eemalda mudelist Sonepikkus ja lisa Sonepikkus_log
m8.glm <- update(m7.glm, ~ . - Sonepikkus + Sonepikkus_log)
anova(m8.glm, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Pron
## 
## Terms added sequentially (first to last)
## 
## 
##                Df Deviance Resid. Df Resid. Dev Pr(>Chi)    
## NULL                            4065     5617.4             
## Lopp            1  230.070      4064     5387.3  < 2e-16 ***
## Aeg             1   78.957      4063     5308.3  < 2e-16 ***
## Ref_kaugus_num  1   92.268      4062     5216.1  < 2e-16 ***
## Verbiklass      2  107.843      4060     5108.2  < 2e-16 ***
## Murre           4  129.258      4056     4979.0  < 2e-16 ***
## Sonepikkus_log  1   92.194      4055     4886.8  < 2e-16 ***
## Aeg:Verbiklass  2    6.149      4053     4880.6  0.04621 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(m8.glm, test = "Chisq")
## Single term deletions
## 
## Model:
## Pron ~ Lopp + Aeg + Ref_kaugus_num + Verbiklass + Murre + Sonepikkus_log + 
##     Aeg:Verbiklass
##                Df Deviance    AIC     LRT Pr(>Chi)    
## <none>              4880.6 4906.6                     
## Lopp            1   5005.1 5029.1 124.446  < 2e-16 ***
## Ref_kaugus_num  1   4960.5 4984.5  79.876  < 2e-16 ***
## Murre           4   5018.4 5036.4 137.731  < 2e-16 ***
## Sonepikkus_log  1   4970.0 4994.0  89.352  < 2e-16 ***
## Aeg:Verbiklass  2   4886.8 4908.8   6.149  0.04621 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# 'drop = FALSE' säilitab data.frame'i struktuuri
sonepikkus_log <- isik1suur[,"Sonepikkus_log", drop = FALSE]
sonepikkus_log$logit <- log(enn/(1-enn))

ggplot(data = sonepikkus_log, 
       aes(x = Sonepikkus_log, y = logit)) +
  geom_point(size = 0.5, alpha = 0.5) +
  geom_smooth(method = "loess")

Suhe logaritmitud šansi ja logaritmitud sõnapikkuse vahel on lineaarsele lähemal, ent suhte tõlgendamine keerulisem: iga logaritmitud ühiku kohta sõna pikkuses langeb asesõna esinemise logaritmitud šanss. Lisaks ei näita selle mudeli vabaliige enam konteksti, kus tegusõna pikkus on 0, vaid konteksti, kus tegusõna pikkuse logaritm on 0, ehk konteksti, kus tegusõna pikkus on 1 tähemärk.

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 rms. Lühend VIF ehk Variation Inflation Factor hindab, kui palju mingi regressioonikoefitsiendi varieeruvus kasvab, kui tunnused on omavahel seotud. 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.

# vif() töötab lrm-mudeli peal, seega peame enne tegema logaritmitud sõnapikkusega mudelist lrm-variandi
m8.lrm <- lrm(Pron ~ Lopp + Aeg + Sonepikkus_log + Ref_kaugus_num + Verbiklass + Murre + Verbiklass:Aeg, data = isik1suur)
rms::vif(m8.lrm)
##                            Lopp=jah                              Aeg=pr 
##                            3.270791                           15.587441 
##                      Sonepikkus_log                      Ref_kaugus_num 
##                            1.270974                            1.020250 
##          Verbiklass=kommunikatiivne                      Verbiklass=muu 
##                            4.196000                            5.188075 
##                          Murre=KESK                         Murre=RANNA 
##                            2.519176                            1.917617 
##                        Murre=SAARTE                          Murre=VORU 
##                            3.478855                            4.256265 
## Aeg=pr * Verbiklass=kommunikatiivne             Aeg=pr * Verbiklass=muu 
##                            3.585707                           14.584589

# glm-mudelite jaoks võib kasutada paketti "performance" 
# ja funktsiooni check_collinearity(), mis näitab
# ainult peamõjude skoore
# install.packages("performance")
performance::check_collinearity(m8.glm)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##       Parameter  VIF Increased SE
##            Lopp 3.27         1.81
##  Ref_kaugus_num 1.02         1.01
##           Murre 3.27         1.81
##  Sonepikkus_log 1.27         1.13
## 
## Moderate Correlation
## 
##   Parameter  VIF Increased SE
##  Verbiklass 5.73         2.39
## 
## High Correlation
## 
##       Parameter   VIF Increased SE
##             Aeg 15.59         3.95
##  Aeg:Verbiklass 28.72         5.36

Näeme, et multikollineaarsus mõjutab koefitsienti Aeg=pr, Verbiklass=muu ja interaktsiooni koefitsienti Aeg=pr * Verbiklass=muu. Seega oleks ohutum siinkohal sellise andmete kodeerimisskeemi puhul loobuda aja ja verbiklassi interaktsioonist, kui selle panus mudelisse ei ole väga suur.

# Eelmine glm-mudel
m8.glm$call
## glm(formula = Pron ~ Lopp + Aeg + Ref_kaugus_num + Verbiklass + 
##     Murre + Sonepikkus_log + Aeg:Verbiklass, family = "binomial", 
##     data = isik1suur)
# Eelmine lrm-mudel
m8.lrm$call
## lrm(formula = Pron ~ Lopp + Aeg + Sonepikkus_log + Ref_kaugus_num + 
##     Verbiklass + Murre + Verbiklass:Aeg, data = isik1suur)

# Interaktsioonita mudelid
m9.glm <- update(m8.glm, ~ . -Verbiklass:Aeg)
m9.lrm <- update(m8.lrm, ~ . -Verbiklass:Aeg)

# Vaatame, kui palju interaktsioon mudelit paremaks teeb
anova(m9.glm, m8.glm, test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: Pron ~ Lopp + Aeg + Ref_kaugus_num + Verbiklass + Murre + Sonepikkus_log
## Model 2: Pron ~ Lopp + Aeg + Ref_kaugus_num + Verbiklass + Murre + Sonepikkus_log + 
##     Aeg:Verbiklass
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1      4055     4886.8                       
## 2      4053     4880.6  2   6.1493  0.04621 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# C-indeksite võrdlus
m8.lrm$stats["C"]
##         C 
## 0.7335807
m9.lrm$stats["C"]
##         C 
## 0.7331709

Kuna interaktsioon mudelis on statistilise olulisuse piiri peal (ja kui rakendaksime Bonferroni parandust, oleks see ebaoluline) ning interaktsiooniga mudeli C ei ole tegelikult interaktsioonita mudeli omast suurem, siis on parem see interaktsioon mudelist välja jätta.

# Kontrollime uue mudeli peal multikollineaarsust
rms::vif(m9.lrm)
##                   Lopp=jah                     Aeg=pr 
##                   3.272564                   1.160019 
##             Sonepikkus_log             Ref_kaugus_num 
##                   1.266316                   1.019802 
## Verbiklass=kommunikatiivne             Verbiklass=muu 
##                   1.741677                   1.893261 
##                 Murre=KESK                Murre=RANNA 
##                   2.519351                   1.917714 
##               Murre=SAARTE                 Murre=VORU 
##                   3.480109                   4.262993
performance::check_collinearity(m9.glm)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##       Parameter  VIF Increased SE
##            Lopp 3.27         1.81
##             Aeg 1.16         1.08
##  Ref_kaugus_num 1.02         1.01
##      Verbiklass 1.13         1.06
##           Murre 3.27         1.81
##  Sonepikkus_log 1.27         1.13

Nüüd on multikollineaarsuse näitajad kõik korras.

Seletavate tunnuste peamõjud ei muutu.

plot(predictorEffects(m9.glm), type = "response", cols = 2, rows = 3)

Mudeli valideerimine

Igasuguse 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.

Alati ei ole see 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 mudelite tulemusi üldistada, siis saab lõpuks kokku üsna usaldusväärsed andmed.

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

Enne valideerimist aga proovime veel natuke mudeldada. Katsetame ära võimalikud interaktsioonid Lopp:Sonepikkus_log, Lopp:Murre, Lopp:Aeg, ja Ref_kaugus_num:Verbiklass. Vaatame, kas mõnel neist võiks mingit mõju olla.

m10.glm <- update(m9.glm, ~ . + Sonepikkus_log:Lopp + Lopp:Murre + Lopp:Aeg + Ref_kaugus_num:Verbiklass)

library(visreg)
visreg(m10.glm, xvar = "Sonepikkus_log", by = "Lopp")

visreg(m10.glm, xvar = "Lopp", by = "Murre")

visreg(m10.glm, xvar = "Lopp", by = "Aeg")

visreg(m10.glm, xvar = "Ref_kaugus_num", by = "Verbiklass")


anova(m10.glm, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Pron
## 
## Terms added sequentially (first to last)
## 
## 
##                           Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                       4065     5617.4              
## Lopp                       1  230.070      4064     5387.3 < 2.2e-16 ***
## Aeg                        1   78.957      4063     5308.3 < 2.2e-16 ***
## Ref_kaugus_num             1   92.268      4062     5216.1 < 2.2e-16 ***
## Verbiklass                 2  107.843      4060     5108.2 < 2.2e-16 ***
## Murre                      4  129.258      4056     4979.0 < 2.2e-16 ***
## Sonepikkus_log             1   92.194      4055     4886.8 < 2.2e-16 ***
## Lopp:Sonepikkus_log        1    0.697      4054     4886.1  0.403790    
## Lopp:Murre                 4   12.850      4050     4873.2  0.012030 *  
## Lopp:Aeg                   1   10.359      4049     4862.9  0.001288 ** 
## Ref_kaugus_num:Verbiklass  2    0.725      4047     4862.2  0.696044    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(m10.glm, test = "Chisq")
## Single term deletions
## 
## Model:
## Pron ~ Lopp + Aeg + Ref_kaugus_num + Verbiklass + Murre + Sonepikkus_log + 
##     Lopp:Sonepikkus_log + Lopp:Murre + Lopp:Aeg + Ref_kaugus_num:Verbiklass
##                           Df Deviance    AIC     LRT Pr(>Chi)   
## <none>                         4862.2 4900.2                    
## Lopp:Sonepikkus_log        1   4862.2 4898.2  0.0001 0.992953   
## Lopp:Murre                 4   4876.0 4906.0 13.8121 0.007920 **
## Lopp:Aeg                   1   4872.5 4908.5 10.3105 0.001323 **
## Ref_kaugus_num:Verbiklass  2   4862.9 4896.9  0.7247 0.696044   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Võiksime proovida jätta sisse interaktsioonid Lopp:Murre ja Lopp:Aeg, ehkki need ei tee kumbki mudelit nii palju paremaks kui peamõjud (vt nt anova() väljundis Deviance’i).

m10.glm <- update(m10.glm, ~ . - Sonepikkus_log:Lopp - Ref_kaugus_num:Verbiklass)
anova(m10.glm, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Pron
## 
## Terms added sequentially (first to last)
## 
## 
##                Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                            4065     5617.4              
## Lopp            1  230.070      4064     5387.3 < 2.2e-16 ***
## Aeg             1   78.957      4063     5308.3 < 2.2e-16 ***
## Ref_kaugus_num  1   92.268      4062     5216.1 < 2.2e-16 ***
## Verbiklass      2  107.843      4060     5108.2 < 2.2e-16 ***
## Murre           4  129.258      4056     4979.0 < 2.2e-16 ***
## Sonepikkus_log  1   92.194      4055     4886.8 < 2.2e-16 ***
## Lopp:Murre      4   12.863      4051     4873.9 0.0119655 *  
## Lopp:Aeg        1   11.044      4050     4862.9 0.0008898 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(m10.glm, test = "Chisq")
## Single term deletions
## 
## Model:
## Pron ~ Lopp + Aeg + Ref_kaugus_num + Verbiklass + Murre + Sonepikkus_log + 
##     Lopp:Murre + Lopp:Aeg
##                Df Deviance    AIC     LRT  Pr(>Chi)    
## <none>              4862.9 4894.9                      
## Ref_kaugus_num  1   4941.1 4971.1  78.218 < 2.2e-16 ***
## Verbiklass      2   4998.3 5026.3 135.414 < 2.2e-16 ***
## Sonepikkus_log  1   4950.4 4980.4  87.521 < 2.2e-16 ***
## Lopp:Murre      4   4876.7 4900.7  13.806 0.0079411 ** 
## Lopp:Aeg        1   4873.9 4903.9  11.044 0.0008898 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m9.glm, m10.glm, test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: Pron ~ Lopp + Aeg + Ref_kaugus_num + Verbiklass + Murre + Sonepikkus_log
## Model 2: Pron ~ Lopp + Aeg + Ref_kaugus_num + Verbiklass + Murre + Sonepikkus_log + 
##     Lopp:Murre + Lopp:Aeg
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1      4055     4886.8                          
## 2      4050     4862.9  5   23.907 0.0002263 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Kontrollime multikollineaarsust.

# lrm-mudelid
m10.lrm <- update(m9.lrm, ~ . + Lopp:Murre + Lopp:Aeg)
rms::vif(m10.lrm)
##                   Lopp=jah                     Aeg=pr 
##                 141.459774                   3.440939 
##             Sonepikkus_log             Ref_kaugus_num 
##                   1.260849                   1.020023 
## Verbiklass=kommunikatiivne             Verbiklass=muu 
##                   1.758000                   1.906485 
##                 Murre=KESK                Murre=RANNA 
##                 191.902547                 133.737896 
##               Murre=SAARTE                 Murre=VORU 
##                 121.735717                  93.196776 
##      Lopp=jah * Murre=KESK     Lopp=jah * Murre=RANNA 
##                 193.229034                 134.083255 
##    Lopp=jah * Murre=SAARTE      Lopp=jah * Murre=VORU 
##                  62.597767                   1.487646 
##          Lopp=jah * Aeg=pr 
##                   3.933321

# glm-mudelid
performance::check_collinearity(m10.glm)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##       Parameter  VIF Increased SE
##             Aeg 3.44         1.85
##  Ref_kaugus_num 1.02         1.01
##      Verbiklass 1.13         1.06
##  Sonepikkus_log 1.26         1.12
##        Lopp:Aeg 3.93         1.98
## 
## High Correlation
## 
##   Parameter       VIF Increased SE
##        Lopp    141.46        11.89
##       Murre 293679.09       541.92
##  Lopp:Murre 249371.60       499.37

Murrete ja pöördelõpu skoorid on hiiglaslikud, seega peaksime selle interaktsiooni ka siiski välja jätma. Üldiselt ei ole väga suurt vahet, kas töötada glm- või lrm-mudelitega. Küll aga nõuavad teatud toimingud mudelifunktsioonis lisaargumente. Näiteks ei saa lrm-mudeleid otse kasutada funktsioonides visreg() või allEffects().

m10.lrm <- update(m10.lrm, ~ . - Lopp:Murre)
rms::vif(m10.lrm)
##                   Lopp=jah                     Aeg=pr 
##                   3.674987                   3.390416 
##             Sonepikkus_log             Ref_kaugus_num 
##                   1.258537                   1.019319 
## Verbiklass=kommunikatiivne             Verbiklass=muu 
##                   1.751955                   1.903513 
##                 Murre=KESK                Murre=RANNA 
##                   2.510896                   1.916406 
##               Murre=SAARTE                 Murre=VORU 
##                   3.442134                   4.208010 
##          Lopp=jah * Aeg=pr 
##                   3.889191

# VÕI

m10.glm <- update(m10.glm, ~ . - Lopp:Murre)
performance::check_collinearity(m10.glm)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##       Parameter  VIF Increased SE
##            Lopp 3.67         1.92
##             Aeg 3.39         1.84
##  Ref_kaugus_num 1.02         1.01
##      Verbiklass 1.13         1.06
##           Murre 3.27         1.81
##  Sonepikkus_log 1.26         1.12
##        Lopp:Aeg 3.89         1.97

Nüüd on multikollineaarsusega kõik korras.

plot(predictorEffects(m10.glm), type = "response", cols = 2, rows = 3)

Näeme, et pöördelõpu puudumisel on tegusõna ajavormil märksa suurem mõju asesõna väljendamise tõenäosusele kui pöördelõpuga vormides.

Vaatame ka mudeli headuse näitajaid.

# Klassiftseerimistäpsus
enn <- m10.glm$fitted.values
teg <- as.character(isik1suur[,"Pron"])
enn_bin <- ifelse(enn > 0.5, "jah", "ei")
klass_tabel <- table(enn_bin, teg)
(klass_acc <- klass_tabel["ei","ei"]+klass_tabel["jah","jah"])/sum(klass_tabel)
## [1] 0.6714215

# C-indeks
m10.lrm$stats["C"]
##         C 
## 0.7346992

Bootstrap

Kuna meie mudel on praegu treenitud ja testitud samal andmestikul, siis valideerime mudelit bootstrap-meetodiga. Saame kasutada sellest samast rms paketist funktsiooni validate(), ent esmalt peame oma lrm-mudelisse lisama argumendid x = TRUE ja y = TRUE. 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.

m10.lrm <- update(m10.lrm, x = TRUE, y = TRUE)

# valideerimine võtab natuke aega
# B näitab meie andmestikust juhuslikult
# kordustega võetud valimite arvu
rms::validate(m10.lrm, B = 200)
##           index.orig training   test optimism index.corrected   n
## Dxy           0.4694   0.4740 0.4663   0.0077          0.4616 200
## R2            0.2224   0.2272 0.2197   0.0075          0.2149 200
## Intercept     0.0000   0.0000 0.0036  -0.0036          0.0036 200
## Slope         1.0000   1.0000 0.9769   0.0231          0.9769 200
## Emax          0.0000   0.0000 0.0058   0.0058          0.0058 200
## D             0.1819   0.1863 0.1795   0.0068          0.1752 200
## U            -0.0005  -0.0005 0.0001  -0.0005          0.0001 200
## Q             0.1824   0.1868 0.1795   0.0073          0.1751 200
## B             0.2072   0.2063 0.2078  -0.0015          0.2087 200
## g             1.1171   1.1345 1.1076   0.0270          1.0902 200
## gp            0.2335   0.2361 0.2321   0.0040          0.2295 200

Oluline on tulp optimism, mis näitab, kui palju on mudel treeningandmestikel ridade statistikuid üle hinnanud. Meid huvitavad selles tulbas read, kus on mudeli headuse näitajad Dxy (seotud mäletatavasti C-indeksiga), 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.

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 kasutame funktsiooni sample() (selleks, et saada samasugust juhuvalimit ka teistel kordadel koodi jooksutamisel, määrame ka seemne (seed)).

# 70% vaatluste arvust on
round(nrow(isik1suur)*0.7)
## [1] 2846

# 2846 juhuslikku ridade indeksit andmestikust
set.seed(1)
trenn_index <- sample(1:nrow(isik1suur), 
                      round(nrow(isik1suur)*0.7), 
                      replace = FALSE)

# treeningandmestik
trenn <- isik1suur[trenn_index,]

# testandmestik
test <- isik1suur[-trenn_index,]

Teeme treeningandmetel mudeli.

m_trenn.glm <- glm(Pron ~ Lopp + Aeg + Ref_kaugus_num + Verbiklass + Murre + Sonepikkus_log + Lopp:Aeg, data = isik1suur, family = "binomial")
m_trenn.lrm <- lrm(Pron ~ Lopp + Aeg + Ref_kaugus_num + Verbiklass + Murre + Sonepikkus_log + Lopp:Aeg, data = isik1suur)

anova(m_trenn.glm, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Pron
## 
## Terms added sequentially (first to last)
## 
## 
##                Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                            4065     5617.4              
## Lopp            1  230.070      4064     5387.3 < 2.2e-16 ***
## Aeg             1   78.957      4063     5308.3 < 2.2e-16 ***
## Ref_kaugus_num  1   92.268      4062     5216.1 < 2.2e-16 ***
## Verbiklass      2  107.843      4060     5108.2 < 2.2e-16 ***
## Murre           4  129.258      4056     4979.0 < 2.2e-16 ***
## Sonepikkus_log  1   92.194      4055     4886.8 < 2.2e-16 ***
## Lopp:Aeg        1   10.101      4054     4876.7  0.001482 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(m_trenn.glm, test = "Chisq")
## Single term deletions
## 
## Model:
## Pron ~ Lopp + Aeg + Ref_kaugus_num + Verbiklass + Murre + Sonepikkus_log + 
##     Lopp:Aeg
##                Df Deviance    AIC     LRT  Pr(>Chi)    
## <none>              4876.7 4900.7                      
## Ref_kaugus_num  1   4955.8 4977.8  79.090 < 2.2e-16 ***
## Verbiklass      2   5012.9 5032.9 136.172 < 2.2e-16 ***
## Murre           4   5013.8 5029.8 137.081 < 2.2e-16 ***
## Sonepikkus_log  1   4966.3 4988.3  89.572 < 2.2e-16 ***
## Lopp:Aeg        1   4886.8 4908.8  10.101  0.001482 ** 
## ---
## 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[,"Pron"])
enn_bin <- ifelse(enn > 0.5, "jah", "ei")
(klass_acc <- klass_tabel["ei","ei"]+klass_tabel["jah","jah"])/sum(klass_tabel)
## [1] 0.6714215

# C-indeks/AUC
# install.packages("ROCR")
library(ROCR)
enn2 <- prediction(enn, test$Pron)
auc <- performance(enn2, measure = "auc")
(auc@y.values[[1]])
## [1] 0.7405512

# C-indeksit võib arvutada ka nii
# install.packages("Hmisc")
library(Hmisc)
somers2(enn, as.numeric(test$Pron)-1)["C"]
##         C 
## 0.7405512

Segamudelid

Nagu 10. 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 uurijale omaette huvi.

Lisame oma regressioonimudelisse individuaalse kõneleja kui juhusliku faktori. Kuna üldjuhul on murdekorpuses ühe faili kohta üks kõneleja, siis saame kasutada tunnust Fail.

Tegeleme siin juhuslike vabaliikmetega (random intercept), mis korrigeeritakse iga juhusliku faktori taseme jaoks (nt iga indiviidi jaoks algaks koefitsientide võrdlemine natuke erinevalt tasemelt). Mudelisse võib lisada ka juhuslikke kaldeid (random slope) mis korrigeerivad valitud fikseeritud faktori efekti vastavalt iga juhusliku faktori tasemele (nt iga indiviidi kohta on mingi tunnuse mõju erinev).

library(lme4)
# mudeli tegemine võib võtta pisut aega
# mudelis võiks täpsustada ka argumendi `control = glmerControl(optimizer = "bobyqa"))`,
# vastasel juhul võib saada veateate, et "model failed to converge"
m1.glmer <- glmer(Pron ~ Lopp + Aeg + Sonepikkus_log + Ref_kaugus_num + Verbiklass + Murre + Lopp:Aeg + (1|Fail), data = isik1suur, family = "binomial", control = glmerControl(optimizer= "bobyqa"))
summary(m1.glmer)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Pron ~ Lopp + Aeg + Sonepikkus_log + Ref_kaugus_num + Verbiklass +  
##     Murre + Lopp:Aeg + (1 | Fail)
##    Data: isik1suur
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   4619.1   4701.1  -2296.5   4593.1     4053 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.3644 -0.7231  0.2137  0.7064  3.4781 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Fail   (Intercept) 0.8493   0.9216  
## Number of obs: 4066, groups:  Fail, 161
## 
## Fixed effects:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                4.463256   0.427770  10.434  < 2e-16 ***
## Loppjah                   -1.566935   0.181706  -8.623  < 2e-16 ***
## Aegpr                      0.748670   0.148308   5.048 4.46e-07 ***
## Sonepikkus_log            -1.285424   0.148494  -8.656  < 2e-16 ***
## Ref_kaugus_num             0.045763   0.005843   7.832 4.81e-15 ***
## Verbiklasskommunikatiivne  0.418451   0.216078   1.937 0.052798 .  
## Verbiklassmuu             -0.912629   0.158773  -5.748 9.03e-09 ***
## MurreKESK                 -1.022764   0.300351  -3.405 0.000661 ***
## MurreRANNA                -0.841752   0.350364  -2.403 0.016283 *  
## MurreSAARTE               -0.508185   0.304862  -1.667 0.095527 .  
## MurreVORU                 -2.082659   0.362207  -5.750 8.93e-09 ***
## Loppjah:Aegpr             -0.533622   0.177945  -2.999 0.002710 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Loppjh Aegpr  Snpkk_ Rf_kg_ Vrbklssk Vrbklssm MrKESK MRANNA
## Loppjah     -0.319                                                            
## Aegpr       -0.232  0.176                                                     
## Sonepkks_lg -0.629 -0.119  0.184                                              
## Ref_kags_nm -0.095 -0.013  0.039 -0.006                                       
## Vrbklsskmmn -0.308  0.000  0.075  0.103  0.020                                
## Verbiklassm -0.485 -0.045  0.123  0.251  0.035  0.655                         
## MurreKESK   -0.456  0.021 -0.006  0.001  0.002  0.009    0.003                
## MurreRANNA  -0.387  0.005 -0.003  0.016 -0.027  0.004   -0.012    0.546       
## MurreSAARTE -0.568  0.306  0.011  0.007  0.005  0.006   -0.008    0.622  0.535
## MurreVORU   -0.571  0.440 -0.012  0.032 -0.003 -0.006   -0.006    0.535  0.454
## Loppjh:Agpr  0.067 -0.288 -0.793 -0.012 -0.001 -0.006    0.040   -0.006  0.002
##             MSAART MrVORU
## Loppjah                  
## Aegpr                    
## Sonepkks_lg              
## Ref_kags_nm              
## Vrbklsskmmn              
## Verbiklassm              
## MurreKESK                
## MurreRANNA               
## MurreSAARTE              
## MurreVORU    0.667       
## Loppjh:Agpr -0.011  0.011

Mudeli väljundist näeme juba paljusid tuttavaid asju. Lisaks on seal tabel juhuslike mõjude kohta (Random effects), mis näitab vabaliikme 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 informandid) ning seda väiksem mõju juhuslikul faktoril mudelis on. Allpool esitab glmer-mudel ka fikseeritud mõjude omavahelise korrelatsiooni, mille põhjal saab tuvastada ka multikollineaarsust.

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 pöördelõpu koefitsient näitab šansside logaritmitud muutust JUHUL, KUI ajavorm on minevik, tegusõna on 1 tähemärgi pikkune (ehk log(0)), eelmine viide on 0 sõna kaugusel, tegusõna on kognitiivne, murre on idamurre JA tegu on 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.

Põhimõtteliselt on ka logistilise regressiooni mudelite puhul võimalik Tukey post-hoc testiga võrrelda ühe seletava tunnuse tasemeid omavahel. Näiteks:

library(multcomp)
summary(glht(m1.glmer, mcp(Murre="Tukey")))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: glmer(formula = Pron ~ Lopp + Aeg + Sonepikkus_log + Ref_kaugus_num + 
##     Verbiklass + Murre + Lopp:Aeg + (1 | Fail), data = isik1suur, 
##     family = "binomial", control = glmerControl(optimizer = "bobyqa"))
## 
## Linear Hypotheses:
##                     Estimate Std. Error z value Pr(>|z|)    
## KESK - IDA == 0      -1.0228     0.3004  -3.405  0.00585 ** 
## RANNA - IDA == 0     -0.8418     0.3504  -2.403  0.11168    
## SAARTE - IDA == 0    -0.5082     0.3049  -1.667  0.44794    
## VORU - IDA == 0      -2.0827     0.3622  -5.750  < 0.001 ***
## RANNA - KESK == 0     0.1810     0.3132   0.578  0.97767    
## SAARTE - KESK == 0    0.5146     0.2631   1.956  0.28193    
## VORU - KESK == 0     -1.0599     0.3242  -3.270  0.00923 ** 
## SAARTE - RANNA == 0   0.3336     0.3184   1.048  0.82917    
## VORU - RANNA == 0    -1.2409     0.3724  -3.332  0.00749 ** 
## VORU - SAARTE == 0   -1.5745     0.2772  -5.679  < 0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Üksteisest erinevad oluliselt niisiis

idamurre ja keskmurre,
idamurre ja Võru murre,
Võru murre ja keskmurre,
Võru murre ja rannamurre,
Võru murre ja saarte murre.

Kas juhusliku vabaliikme lisamine mudelisse aitab mudelil paremini uuritava tunnuse varieerumist seletada? Kuna me ei saa anova()-ga otseselt võrrelda ilma juhuslike mõjudeta glm-mudeleid ja juhuslike mõjudega glmer-mudeleid, siis peame natuke trikitama. Selleks tuleb esmalt teha testmudel, milles kõik vaatlused pärineksid justkui ühelt kõnelejalt. Samas peame lisama argumendi check.nlev.gtr.1 = "ignore", mis ütleb, et mudel ignoreeriks seda, et juhuslike mõjudena on lisatud ainult üks subjekt.

idconst <- factor(rep(1, nrow(isik1suur)))
head(idconst)
## [1] 1 1 1 1 1 1
## Levels: 1

m0.glmer <- glmer(Pron ~ Lopp + Aeg + Sonepikkus_log + Ref_kaugus_num + Verbiklass + Murre + Lopp:Aeg + (1|idconst), data = isik1suur, family = "binomial", control = glmerControl(optimizer = "bobyqa", check.nlev.gtr.1 = "ignore"))

anova(m0.glmer, m1.glmer, test = "Chisq")
## Data: isik1suur
## Models:
## m0.glmer: Pron ~ Lopp + Aeg + Sonepikkus_log + Ref_kaugus_num + Verbiklass + 
## m0.glmer:     Murre + Lopp:Aeg + (1 | idconst)
## m1.glmer: Pron ~ Lopp + Aeg + Sonepikkus_log + Ref_kaugus_num + Verbiklass + 
## m1.glmer:     Murre + Lopp:Aeg + (1 | Fail)
##          Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)    
## m0.glmer 13 4902.7 4984.7 -2438.3   4876.7                            
## m1.glmer 13 4619.1 4701.1 -2296.5   4593.1 283.6      0  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Tuleb välja, et juhuslik vabaliige teeb mudelit paremaks ja vähendab AIC-d oluliselt.

Vaatame ka mudeli headuse näitajaid.

Klassifitseerimistäpsus

enn <- fitted(m1.glmer)
teg <- as.character(isik1suur[,"Pron"])
enn2 <- ifelse(enn>0.5, "jah", "ei")
klass_tabel <- table(enn2, teg)
(klass_tabel[1]+klass_tabel[4])/sum(klass_tabel)
## [1] 0.7238072

C-indeks (AUC)

library(Hmisc)
somers2(enn, as.numeric(isik1suur$Pron)-1)["C"]
##         C 
## 0.8096971

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 informant asesõna kasutamises keskmisest asesõna kasutamisest erineb.

juh <- ranef(m1.glmer, condVar = TRUE)
head(juh$Fail)
##                                           (Intercept)
## AKS_Jaan_Koort_synt-kontrollitud.xml       -0.3072526
## AKS_Jaan_Ounap_synt-kontrollitud.xml        1.4173364
## AMB_AnnaMaria_Toome_synt-kontrollitud.xml   0.4697477
## AMB_Juuli_Tomp_synt1.xml                    1.0625399
## AMB_Juuli_Tomp_synt2.xml                    1.4616972
## AMB_Marie_Raude_synt.xml                    0.2743242
mean(juh[[1]]$`(Intercept)`) # keskmine on 0 lähedal
## [1] -0.007220133
median(juh[[1]]$`(Intercept)`) # mediaan on 0 lähedal
## [1] 0.09238185

Võime teha vabaliikmete kohanduste jaotumise joonise, nt paketiga lattice

# install.packages("lattice")
lattice::dotplot(ranef(m1.glmer))
## $Fail

Või ka ggplotiga.

# Teeme juhuslike mõjude objektist andmetabeli
juh_df <- as.data.frame(juh)
head(juh_df)
##   grpvar        term                                       grp    condval
## 1   Fail (Intercept)      AKS_Jaan_Koort_synt-kontrollitud.xml -0.3072526
## 2   Fail (Intercept)      AKS_Jaan_Ounap_synt-kontrollitud.xml  1.4173364
## 3   Fail (Intercept) AMB_AnnaMaria_Toome_synt-kontrollitud.xml  0.4697477
## 4   Fail (Intercept)                  AMB_Juuli_Tomp_synt1.xml  1.0625399
## 5   Fail (Intercept)                  AMB_Juuli_Tomp_synt2.xml  1.4616972
## 6   Fail (Intercept)                  AMB_Marie_Raude_synt.xml  0.2743242
##      condsd
## 1 0.2487469
## 2 0.3502290
## 3 0.3749068
## 4 0.2287077
## 5 0.3244283
## 6 0.4706823

library(ggplot2)
ggplot(data = juh_df, 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") +
  labs(x = "Kohandus vabaliikmes",
       y = "Informandid") +
  theme_minimal() +
  theme(axis.text.y = element_blank(),
        panel.grid = element_blank(),
        panel.border = element_rect(fill = "transparent"))

Võime ka keskenduda nendele informantidele, kes kasutavad asesõnu keskmiselt oluliselt rohkem, või nendele, kes oluliselt vähem.

ggplot(data = juh_df, 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") +
  labs(x = "Kohandus vabaliikmes",
       y = "Informandid",
       title = "Informandid, kes kasutavad asesõnu vähem") +
  theme_minimal() +
  theme(panel.grid = element_blank(),
        panel.border = element_rect(fill = "transparent"),
        plot.title = element_text(hjust = 0.5)) +
  coord_cartesian(ylim = c(1, 6))

ggplot(data = juh_df, 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") +
  labs(x = "Kohandus vabaliikmes",
       y = "Informandid",
       title = "Informandid, kes kasutavad asesõnu rohkem") +
  theme_minimal() +
  theme(panel.grid = element_blank(),
        panel.border = element_rect(fill = "transparent"),
        plot.title = element_text(hjust = 0.5)) +
  coord_cartesian(ylim = c(152, 161))

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.205
##   Conditional ICC: 0.152

Kohandatud ICC-väärtus (Adjusted ICC) on seotud ainult juhuslike mõjudega, tingimuslik ICC (Conditional ICC) aga arvestab ka fikseeritud mõjudega. Kui huvi pakuvad ainult juhuslikud mõjud, siis tuleks vaadata kohandatud ICC-d. ICC väärtus ulatub 0st 1ni. Kui juhuslik mõju ei toimi andmetes üldse mingi grupeeriva tunnusena (= kõik vaatlused on ikkagi erinevad), siis on indeksi väärtus 0, ja kui kõik vaatlused kuuluvad ühte gruppi, on indeksi väärtus 1. Juhusliku faktori poolt kirjeldatud osa varieerumisest on siin ~0.21 ehk 21%.

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.

arvulised <- isik1suur[, c("Sonepikkus_log", "Ref_kaugus_num")]
arvulised$logit <- predict(m1.glmer, isik1suur, type = "link")
arvulised_pikk <- tidyr::gather(arvulised, key = "arv_tunnus", value = "arv_vaartus", -logit)
ggplot(data = arvulised_pikk, 
       aes(y = logit, x = arv_vaartus))+
  geom_point(size = 0.5, alpha = 0.5) +
  facet_wrap("arv_tunnus", scales = "free") +
  geom_smooth(method = "loess")

performance::check_collinearity(m1.glmer)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##       Parameter  VIF Increased SE
##            Lopp 1.69         1.30
##             Aeg 3.09         1.76
##  Sonepikkus_log 1.17         1.08
##  Ref_kaugus_num 1.01         1.00
##      Verbiklass 1.13         1.06
##           Murre 1.53         1.24
##        Lopp:Aeg 3.18         1.78

VIF-skoorid on < 5, seega ei pea tunnuste seotuse pärast muretsema.

Juhuslikke kaldeid tunnuse Fail põhjalt on keeruline lisada, kuna meil on Fail-il väga palju tasemeid, mudelis mitu mitme tasemega fikseeritud tunnust ning esineb väga harvu tunnusekombinatsioone. Võiksime ju näiteks oletada, et kuna informandid on oma murdetaustast ja üldkeelest erineval määral mõjutatud, on ka murdeala mõju keelekasutusele varieeruv. Kui üritame aga seda juhusliku kaldena mudeldada, saame teate, et mudel ei suuda parameetreid hinnata, sest meil pole selleks lihtsalt piisavalt palju andmeid. Mudel küll tehakse ära ja saame selle väljundit vaadata, ent me ei saa olla päris kindlad selle täpsuses.

# mudeli tegemine võib võtta natuke aega
m2.glmer <- glmer(Pron ~ Lopp + Aeg + Sonepikkus_log + Ref_kaugus_num + Verbiklass + Murre + Lopp:Aeg + (1 + Murre|Fail), data = isik1suur, family = "binomial", control = glmerControl(optimizer = "bobyqa"))
## Warning in optwrap(optimizer, devfun, start, rho$lower, control = control, :
## convergence code 1 from bobyqa: bobyqa -- maximum number of function evaluations
## exceeded
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 5 negative eigenvalues
summary(m1.glmer)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Pron ~ Lopp + Aeg + Sonepikkus_log + Ref_kaugus_num + Verbiklass +  
##     Murre + Lopp:Aeg + (1 | Fail)
##    Data: isik1suur
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   4619.1   4701.1  -2296.5   4593.1     4053 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.3644 -0.7231  0.2137  0.7064  3.4781 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Fail   (Intercept) 0.8493   0.9216  
## Number of obs: 4066, groups:  Fail, 161
## 
## Fixed effects:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                4.463256   0.427770  10.434  < 2e-16 ***
## Loppjah                   -1.566935   0.181706  -8.623  < 2e-16 ***
## Aegpr                      0.748670   0.148308   5.048 4.46e-07 ***
## Sonepikkus_log            -1.285424   0.148494  -8.656  < 2e-16 ***
## Ref_kaugus_num             0.045763   0.005843   7.832 4.81e-15 ***
## Verbiklasskommunikatiivne  0.418451   0.216078   1.937 0.052798 .  
## Verbiklassmuu             -0.912629   0.158773  -5.748 9.03e-09 ***
## MurreKESK                 -1.022764   0.300351  -3.405 0.000661 ***
## MurreRANNA                -0.841752   0.350364  -2.403 0.016283 *  
## MurreSAARTE               -0.508185   0.304862  -1.667 0.095527 .  
## MurreVORU                 -2.082659   0.362207  -5.750 8.93e-09 ***
## Loppjah:Aegpr             -0.533622   0.177945  -2.999 0.002710 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Loppjh Aegpr  Snpkk_ Rf_kg_ Vrbklssk Vrbklssm MrKESK MRANNA
## Loppjah     -0.319                                                            
## Aegpr       -0.232  0.176                                                     
## Sonepkks_lg -0.629 -0.119  0.184                                              
## Ref_kags_nm -0.095 -0.013  0.039 -0.006                                       
## Vrbklsskmmn -0.308  0.000  0.075  0.103  0.020                                
## Verbiklassm -0.485 -0.045  0.123  0.251  0.035  0.655                         
## MurreKESK   -0.456  0.021 -0.006  0.001  0.002  0.009    0.003                
## MurreRANNA  -0.387  0.005 -0.003  0.016 -0.027  0.004   -0.012    0.546       
## MurreSAARTE -0.568  0.306  0.011  0.007  0.005  0.006   -0.008    0.622  0.535
## MurreVORU   -0.571  0.440 -0.012  0.032 -0.003 -0.006   -0.006    0.535  0.454
## Loppjh:Agpr  0.067 -0.288 -0.793 -0.012 -0.001 -0.006    0.040   -0.006  0.002
##             MSAART MrVORU
## Loppjah                  
## Aegpr                    
## Sonepkks_lg              
## Ref_kags_nm              
## Vrbklsskmmn              
## Verbiklassm              
## MurreKESK                
## MurreRANNA               
## MurreSAARTE              
## MurreVORU    0.667       
## Loppjh:Agpr -0.011  0.011

Sellisel juhul tuleks proovida mudelist ükshaaval välja visata juhuslikke mõjusid ja fikseeritud mõjusid ja vaadata, milline mudel lõpuks kõige paremini sobitab. Võib juhtuda näiteks, et kui visata kaks fikseeritud tunnust välja ja jätta juhuslik mõju sisse, võib mudel sobitada paremini kui siis, kui visata juhuslik mõju välja.

Kordamisküsimused

1. Mida näitab klassifitseerimistäpsus?

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

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

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

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

Seda, et ühelt subjektilt ei oleks andmestikus mitut vaatlust.
Seda, et uuritav tunnus oleks normaaljaotusega.
Seda, et arvulised seletavad tunnused oleksid normaaljaotusega.
Seda, et mudeli ennustatud šansi logaritmi ja arvulise seletava tunnuse vaheline suhe oleks lineaarne.
Seda, et seletavad tunnused ei seletaks uuritava tunnuse varieerumisest sama asja.
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?

Sellise mudeli tegemist, kus andmepunktidele ennustatud väärtuste jaotus on saapakujuline.
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.
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")

nahavärv
vanus
sugu
ükski ei tee

6. Milline nende kolme tunnuse interaktsioonist teeb mudelit paremaks?

colour*age
colour*sex
age*sex
ü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")

0.01
0.19
0.82
0.83

8. Mis on selle mudeli C/AUC?

pole võimalik arvutada
0.73
0.83

9. Kas mudelis esineb multikollineaarsust?

Ei.
Jah, mingil määral.
Jah, olulisel määral.

Järgmisel korral

Tingimuslikud otsustuspuud

13. Logistilise regressiooni mudelid (3). Segamudelid (2) (01.04.2020)