Logistilise regressiooni mudelid

Tänases praktikumis

Mitme tunnusega regressioonimudel
Mudeldamisstrateegiad ja mudelite võrdlemine
Mudeli mõjude visualiseerimine
Interaktsioonide lisamine

Eelmisel korral vaatasime keeleteaduse andmestikku ade_peal.csv ning püüdsime seletada, kuidas sõltub eesti murretes kohasuhete väljendamisel peal-konstruktsiooni (vs. alaleütleva käände) kasutamine sellest, kas koht, mille peal miski on, on liigutatav või staatiline. Saime teada, et peal-konstruktsiooni kasutamise šanss on koos staatiliste kohtadega väiksem kui siis, kui koht on vähemalt mõtteliselt hõlpsasti liigutatav. Niisiis esineks fraas heinamaa peal vähem tõenäoliselt kui fraas laua peal.

Vaatame sel korral sama andmestiku tasakaalus versiooni ja lisame mudelisse nüüd ka teisi seletavaid tunnuseid.

Esmalt loeme algse andmestiku sisse.

library(tidyverse) # läheb vaja pisut hiljem
ade_peal <- read.delim("data/ade_peal.csv", encoding = "UTF-8", stringsAsFactors = TRUE)

Tuletame meelde, mis tunnused andmestikus on.

summary(ade_peal)

##           TR_liik       tegusõna         tegusõna_klass    LM_lemma   
##  asesõna      :516   olema  :865   liikumisverb :297    maa    : 125  
##  nimisõnafraas:855   käima  :175   olemisverb   :868    koht   : 119  
##  verbifraas   :661   tegema : 82   paiknemisverb: 61    laud   :  86  
##                      panema : 41   tegevusverb  :618    põld   :  85  
##                      kasvama: 38   verbita      :188    meri   :  78  
##                      (Other):644                        põrand :  55  
##                      NA's   :187                        (Other):1484  
##      LM_mobiilsus   LM_komplekssus   LM_pikkus     LM_pikkus_log    sõnajärg   
##  liigutatav: 729   lihtsõna:1833   Min.   :1.000   Min.   :0.000   lm_tr: 613  
##  staatiline:1303   liitsõna: 199   1st Qu.:2.000   1st Qu.:1.000   tr_lm:1419  
##                                    Median :2.000   Median :1.000               
##                                    Mean   :2.655   Mean   :1.261               
##                                    3rd Qu.:3.000   3rd Qu.:1.585               
##                                    Max.   :9.000   Max.   :3.170               
##                                                                                
##     cx           murre       murderühm  
##  ade : 722   Saarte :457   Kesk   :880  
##  peal:1310   Lääne  :380   Lõuna  :430  
##              Kesk   :379   Ranniku:265  
##              Ranna  :161   Saarte :457  
##              Võru   :154                
##              Ida    :121                
##              (Other):380                
##                                     kõneleja   
##  KJ_JOE_Jyri_Maisa_synt-kontrollitud.xml:  27  
##  KJ_PLT_Liina_Roosimaa_EMH1427_lihts.xml:  26  
##  KJ_PIL_Helene_Pikk_synt-parandamata.xml:  25  
##  KJ_SJN_Marie_Kees_EMH16_synt.xml       :  25  
##  KJ_HMD_Aleksander_Idarand_synt.xml     :  24  
##  KJ_KAI_Leena_Lepamaa_synt.xml          :  23  
##  (Other)                                :1882

Tasakaalustamata ja tasakaalus valim

Praegune andmestik on uuritava tunnuse suhtes tasakaalustamata. See tähendab, et uuritava tunnuse tasemed (ade ja peal) ei esine võrdse sagedusega.
Tasakaalustamata valimite kasutamine aitab keele- või mis tahes muu nähtuse kirjeldamisel võtta arvesse ka selle nähtuse avaldumisvariantide loomulikku jaotumist populatsioonis. Antud juhul võib tasakaalustamata valimi põhjal seega oletada, et kui mitte võtta arvesse igasugust muud infot keelelise ja mittekeelelise konteksti kohta, oleks kõnelejate jaoks konkreetsete kohatähenduste väljendamisel peal-konstruktsiooni kasutamine vaikimisi oluliselt loomulikum. Tasakaalustamata valimi põhjal tehtud mudel teaks seda ning ennustaks seega ilmselt paremini järgmist konstruktsiooni, mida kohasuhte väljendamisel kasutataks.
Tasakaalustamata valimi probleem on aga see, et me ei pruugi saada piisavalt infot harvemini esineva variandi kasutamist tingivate piirangute kohta, sest mudel on sagedase variandi poole liialt kaldu ning harvemini esinev variant ei panusta seega ka oluliselt mudeli ennustustäpsusesse.

Teisalt võib meil olla aga ka lihtsalt teoreetiline eeldus, et mingite alternatsioonide korral on mõlemad variandid võrdselt võimalikud, st nende esinemistõenäosus ilma konteksti arvestamata peaks olema 0.5 ja 0.5 ning alles kontekstuaalsed tunnused on need, mis valikut ühe või teise variandi kasuks kallutavad. Tasakaalus valim lähtub sellisest teoreetilisest eeldusest ning võimaldab paremini mõõta erinevate faktorite mõju ka harvemini esineva variandi kasutamisele.
Tasakaalus valimi miinuseks on aga see, et selline teoreetiline eeldus ei pruugi olla tegelikult adekvaatne ning uuritava tunnuse tasemete tegelikku jaotust arvestamata võime mudeliga üle hinnata harvemini esineva variandi kasutust.

Tasakaalus valimi saamiseks on kaks lihtsat moodust (on ka teisi, keerulisemaid):

harvemini esineva variandi üleesindamine (over-sampling), st esinevate vaatluste juhuslik kordamine,
sagedamini esineva variandi alaesindamine (under-sampling), st esinevate vaatluste juhuslik väljajätt.

Esimese meetodiga võib tekkida probleem, et muudame harvemini esineva variandi kasutamise kunstlikult liiga homogeenseks (andmeid on küll rohkem, aga variatiivsus on sama) ning mudel hakkaks mõnel teisel valimil ülesobitama. Teise meetodiga võib jällegi tekkida probleem, et kaotame mingit olulist infot, mudelil on vähem, millest õppida, ning mudel on seetõttu oma ennustustes ebatäpne.

Ehkki üht õiget vastust ega meetodit ei ole, võtame siin nimetatutest teise riski ning tasakaalustame valimi alaesindamise teel nõnda, et võtame andmestiku 1310 peal-konstruktsioonist juhuslikult 722 kasutusjuhtu ehk sama palju, kui andmestikus on käändekonstruktsiooni esinemisjuhtusid.

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

summary(dat)

##           TR_liik       tegusõna         tegusõna_klass    LM_lemma   
##  asesõna      :351   olema  :618   liikumisverb :245    maa    : 103  
##  nimisõnafraas:613   käima  :155   olemisverb   :620    koht   :  81  
##  verbifraas   :480   tegema : 55   paiknemisverb: 44    meri   :  69  
##                      istuma : 27   tegevusverb  :407    põld   :  61  
##                      panema : 25   verbita      :128    laud   :  55  
##                      (Other):436                        põrand :  46  
##                      NA's   :128                        (Other):1029  
##      LM_mobiilsus  LM_komplekssus   LM_pikkus     LM_pikkus_log    sõnajärg   
##  liigutatav:469   lihtsõna:1285   Min.   :1.000   Min.   :0.000   lm_tr: 430  
##  staatiline:975   liitsõna: 159   1st Qu.:2.000   1st Qu.:1.000   tr_lm:1014  
##                                   Median :2.000   Median :1.000               
##                                   Mean   :2.636   Mean   :1.252               
##                                   3rd Qu.:3.000   3rd Qu.:1.585               
##                                   Max.   :9.000   Max.   :3.170               
##                                                                               
##     cx          murre       murderühm  
##  ade :722   Saarte :324   Kesk   :611  
##  peal:722   Kesk   :266   Lõuna  :294  
##             Lääne  :257   Ranniku:215  
##             Ranna  :128   Saarte :324  
##             Võru   :106                
##             Ida    : 88                
##             (Other):275                
##                                     kõneleja   
##  KJ_HMD_Aleksander_Idarand_synt.xml     :  21  
##  KJ_PIL_Helene_Pikk_synt-parandamata.xml:  21  
##  KJ_JOE_Jyri_Maisa_synt-kontrollitud.xml:  19  
##  KJ_JOE_Leena_Toomar_synt-parandatud.xml:  19  
##  KJ_PLT_Liina_Roosimaa_EMH1427_lihts.xml:  19  
##  KJ_SJN_Marie_Kees_EMH16_synt.xml       :  19  
##  (Other)                                :1326

Lihtsam viis juhuvalimit võtta oleks kasutada tidyverse’i paketist funktsiooni sample_n(). NB! Kui jooksutada andmestiku saamiseks allolevat koodi, võib andmestik olla pisut erinev sellest, mida ülemise koodiga saab.

set.seed(1)
dat <- ade_peal %>%
    group_by(cx) %>% # grupeeri vaatlused tunnuse cx järgi 
    sample_n(size = 722, replace = F) %>% # võta kummastki grupist 722 juhuslikku vaatlust
    ungroup() # kaota grupeering

Ülevaade andmetest

Enne mudeli tegemist võiksime vaadata pisut uuritava tunnuse jaotumist muude tunnuse põhjal, mida tahaksime mudelisse kaasata.

names(dat)

##  [1] "TR_liik"        "tegusõna"       "tegusõna_klass" "LM_lemma"      
##  [5] "LM_mobiilsus"   "LM_komplekssus" "LM_pikkus"      "LM_pikkus_log" 
##  [9] "sõnajärg"       "cx"             "murre"          "murderühm"     
## [13] "kõneleja"

# Kategoriaalsed tunnused
seletavad_kat <- c("TR_liik", "tegusõna_klass", "LM_mobiilsus", "LM_komplekssus", "sõnajärg", "murderühm")

# Teeme korraga uuritava tunnuse risttabelid kõikide kat. seletavate tunnustega
# lapply() funktsioon võimaldab anda ette andmestiku kõik kategoriaalsete seletavate tunnuste tulbad: dat[,seletavad_kat]
# ning iga tulba puhul: function(seletav_kat)
# teha sellest tulbast ning uuritava tunnuse tulbast risttabel: table(seletav_kat, dat$cx)
# seletav_kat on siinkohal täiesti suvaline nimi, mida kasutame tulpade kohatäitena, 
# põhimõtteliselt võib see olla ka x, y, trallallaa või midagi muud
# Tulemuseks saame listi, mille iga element on omaette risttabel.
lapply(dat[,seletavad_kat], function(seletav_kat) table(seletav_kat, dat$cx))

## $TR_liik
##                
## seletav_kat     ade peal
##   asesõna       166  185
##   nimisõnafraas 295  318
##   verbifraas    261  219
## 
## $tegusõna_klass
##                
## seletav_kat     ade peal
##   liikumisverb  175   70
##   olemisverb    311  309
##   paiknemisverb  18   26
##   tegevusverb   164  243
##   verbita        54   74
## 
## $LM_mobiilsus
##             
## seletav_kat  ade peal
##   liigutatav 137  332
##   staatiline 585  390
## 
## $LM_komplekssus
##            
## seletav_kat ade peal
##    lihtsõna 602  683
##    liitsõna 120   39
## 
## $sõnajärg
##            
## seletav_kat ade peal
##       lm_tr 216  214
##       tr_lm 506  508
## 
## $murderühm
##            
## seletav_kat ade peal
##     Kesk    273  338
##     Lõuna   108  186
##     Ranniku 157   58
##     Saarte  184  140

# Teeme korraga uuritava tunnuse joonised kõikide kat. seletavate tunnustega
dat %>%
  select(cx, all_of(seletavad_kat)) %>% # vali uuritava ja kat. seletavate tunnuste tulbad
  pivot_longer(., 2:length(.), # tee seletavate tunnuste põhjal (al 2. tulbast) laiast tabelist pikk
               names_to = "tunnused", # kus tunnuste tulbanimed lähevad ühte tulpa
               values_to = "väärtused") %>% # ja nende tasemed/väärtused teise
  ggplot(aes(y = väärtused, # pane tulpdiagrammi y-teljele tunnuste väärtused
               fill = cx)) + # kuva täitevärviga uuritavat tunnust
  geom_bar(position = "fill") + # tee joonis suhteliste sagedustega
  stat_count(geom = "text", # lisa tekstina ka absoluutarvud
             aes(label = after_stat(count)), # tekstisildi jaoks loe kategooriate sagedused kokku
             position = position_fill(vjust = 0.5), # aseta sildid tulpade keskele
             colour = "grey40") + # värvi tekst halliks
  facet_wrap("tunnused", scales = "free") # tee iga tunnusega oma paneel

Sarnase asja tegemiseks võime kasutada ka olemasolevaid ggploti laiendusi, nt funktsiooni ggbivariate() paketist GGally.

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

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

ggbivariate(dat[,c("cx", seletavad_kat)], outcome = "cx")

Näeme, et uuritava tunnuse varieerumine toimub kõikide seletavate tunnuste kõikidel tasemetel. See tähendab, et ühelgi väärtusel ei ole ade või peal suhtes kategoorilist, teist varianti välistavat eelistust. Kui oleks, võiks see regressioonimudelis osutuda probleemiks.
Pane tähele, et uuritava tunnuse jaotumine mingite seletavate tunnuste põhjal on võinud tasakaalus valimis võrreldes algandmestikuga muutuda (nt konstruktsioonide kasutus murderühmades, kus Saarte murderühmas oli algselt rohkem peal-konstruktsioone, nüüd aga ade-konstruktsioone).

Vaatame ka arvulist seletavat tunnust.

# Kohasõna silpide arvu jagunemine käände- ja peal-konstruktsioonides
tapply(dat$LM_pikkus, dat$cx, summary)

## $ade
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   2.000   2.000   2.608   3.000   9.000 
## 
## $peal
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   2.000   2.000   2.665   3.000   9.000

tapply(dat$LM_pikkus_log, dat$cx, summary)

## $ade
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   1.000   1.000   1.247   1.585   3.170 
## 
## $peal
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   1.000   1.000   1.257   1.585   3.170

# Vaatame silpide arvu jagunemist karpdiagrammil
ggplot(dat, aes(x = LM_pikkus_log, y = cx)) +
  geom_boxplot()

# Teisendame kat. uuritava tunnuse 1-deks ja 0-deks
# ja teeme logistilise regressioonijoonega joonise
ggplot(dat, aes(x = LM_pikkus_log, y = ifelse(cx == "peal", 1, 0))) +
  geom_point(alpha = 0.1) +
  geom_smooth(method = "glm",
              method.args = list(family = "binomial"))

## `geom_smooth()` using formula = 'y ~ x'

Andmetest kiire esmase ülevaate saamiseks on R-is erinevaid pakette, nt DataExplorer, GGally või Hmisc, mis hõlbustavad oluliselt andmeanalüütiku igapäevatööd.

Samuti saab esmast kirjeldavat analüüsi küllalt kiiresti teha paketiga sjPlot, mida kasutame hiljem ka mudeli efektide visualiseerimiseks.

#install.packages("sjPlot", dependencies = T)
library(sjPlot)

plot_grpfrq(dat$LM_mobiilsus, dat$cx, drop.empty = TRUE)

dat %>%
  group_by(cx) %>%
  select(LM_mobiilsus) %>%
  plot_stackfrq(show.n = T)

## [[1]]
## [[1]][[1]]
## NULL
## 
## [[1]][[2]]

## 
## 
## [[2]]

Ühe tunnusega mudel tasakaalus valimi põhjal

Kõigepealt teeme mudeli, milles oli ainsa seletava tunnusena kohasõna liigutatavus (LM_mobiilsus).

m1.glm <- glm(cx ~ LM_mobiilsus, data = dat, family = "binomial")
summary(m1.glm)

## 
## Call:
## glm(formula = cx ~ LM_mobiilsus, family = "binomial", data = dat)
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              0.8852     0.1015   8.717   <2e-16 ***
## LM_mobiilsusstaatiline  -1.2906     0.1208 -10.687   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2001.8  on 1443  degrees of freedom
## Residual deviance: 1879.0  on 1442  degrees of freedom
## AIC: 1883
## 
## Number of Fisher Scoring iterations: 4

Kordamiseks:

Mudel ennustab peal-konstruktsiooni esinemise tõenäosust seletavate tunnuste määratud kontekstis, kuna peal on tähestikus vaikimisi tagapool kui ade.
Rida (Intercept) väljendab logaritmitud šansi kaudu peal-konstruktsiooni esinemise tõenäosust seletavate tunnuste baastasemete kontekstis. Mudeli hetkel ainsa seletava tunnuse, LM_mobiilsus baastase on vaikimisi tähestiku järjekorras esimene tase liigutatav. Nullhüpotees selle konteksti kohta on, et peal-konstruktsioonil on seal võrdsed võimalused esineda ja mitte-esineda. Kui nullhüpotees kehtiks, oleksid peal-konstruktsiooni esinemise ja mitte-esinemise (= ade-konstruktsiooni esinemise) tõenäosused kumbki 0.5, peal-konstruktsiooni esinemise šanss seega 0.5/0.5=1 ning selle šansi logaritm tulbas Estimate log(1)=0. Kuna aga Estimate tulbas on hoopis nullist suurem väärtus 0.885154, on peal-konstruktsiooni esinemise šansid liigutatavate kohtadega järelikult suuremad kui 1, täpsemalt exp(m1.glm$coefficients["(Intercept)"])=2.423358, ja peal-konstruktsiooni tõenäosus seega plogis(m1.glm$coefficients["(Intercept)"])=0.7078891. P-väärtus real näitab seda, milline on tõenäosus, et saaksime sellise nullist erineva koefitsiendi juhul, kui nullhüpotees ikkagi kehtiks ja tegelikult alalütleva ja peal-konstruktsiooni kasutuses liigutatavate kohtadega erinevust ei ole. Kui p-väärtus on väiksem kui 0.05, on see tõenäosus väga väike ning saamegi öelda, et erinevus kahe kohakonstruktsiooni kasutuses on vabaliikme kontekstis statistiliselt oluline.
Rida LM_mobiilsusstaatiline väljendab logaritmitud šansisuhte kaudu muutust peal-konstruktsiooni esinemise tõenäosuses, kui muudame seletavate tunnuste kontekstis reanimes täpsustatud tunnuse taseme ära (antud juhul vaatame liigutatavate kohtade asemel hoopis staatilisi kohti). Nullhüpotees sellel real on, et peal-konstruktsiooni võimalused esineda selles uues kontekstis ei erine tema võimalustest esineda vanas, baastasemete kontekstis. Kui nullhüpotees kehtiks, oleksid peal-konstruktsiooni esinemise šansid nii liigutatavate kui ka staatiliste kohtadega sama suured, nende šansside jagatis seega 1 ja selle jagatise logaritm tulbas Estimate log(1)=0. Kuna aga Estimate tulbas on hoopis nullist väiksem väärtus -1.2906192, on peal-konstruktsiooni esinemise šansid staatiliste kohtadega järelikult väiksemad kui liigutatavate kohtadega, täpsemalt exp(m1.glm$coefficients["(Intercept)"]+m1.glm$coefficients["LM_mobiilsusstaatiline"])=0.6666667, ja peal-konstruktsiooni esinemise tõenäosus seega plogis(m1.glm$coefficients["(Intercept)"]+m1.glm$coefficients["LM_mobiilsusstaatiline"])=0.4.

Kui lineaarses regressioonis näitab mudeli väljund kohe ka seda, kas mudel tervikuna on statistiliselt oluline või mitte, siis glm-mudelite pere (Generalized Linear Models), kuhu kuulub ka logistiline regressioon, mudelile ise p-väärtust ei omista. Küll aga on võimalik see p-väärtus välja arvutada n-ö “käsitsi” funktsiooniga pchisq(), kasutades nullmudeli (ilma seletavate tunnusteta mudeli) ning tunnustega mudeli hälbimuse näitajaid (deviance) ja vabadusastmete arvu (df).

pchisq(q = m1.glm$null.deviance-m1.glm$deviance, 
       df = m1.glm$df.null-m1.glm$df.residual, 
       lower.tail = FALSE)

## [1] 1.499986e-28

Saame üldjoontes sama tulemuse, mille tasakaalustamata valimiski: võrreldes liigutatavate kohtadega (nt laud) on staatiliste kohtadega (nt heinamaa) peal-konstruktsiooni kasutamine oluliselt vähem tõenäoline. Erinevus eelmises praktikumis tehtud tasakaalustamata valimiga tehtud mudelist on aga selles, et kui seal oli peal-konstruktsiooni kasutamine mõlemal juhul ikkagi tõenäolisem kui käändekonstruktsiooni kasutamine, siis siin, tasakaalus valimis, see nii enam ei ole. Saame seega tasakaalus valimiga tehtud mudeli põhjal öelda, millised seletavad tunnused mõjutavad uuritava tunnuse varieerumist ja millises suunas, aga šansside ja tõenäosuste väljaarvutamine pole võib-olla ülemäära mõttekas, kuna ei vasta niikuinii päriselt tegelikule keelekasutusele.

Logistilise regressiooni vabaliikme (Intercept) kohta käiv nullhüpotees aga ei erine tasakaalus ja tasakaalustamata valimi puhul: hinnatakse ikkagi tõenäosust, et mõlemad uuritava tunnuse variandid on võrdselt võimalikud.

# Mudel algse, tasakaalustamata andmestikuga,
# kus pole ühtki seletavat tunnust
# ja kus "peal"-konstruktsioon on palju sagedasem kui kääne

m0_alg.glm <- glm(cx ~ 1, ade_peal, family = "binomial")
summary(m0_alg.glm)

## 
## Call:
## glm(formula = cx ~ 1, family = "binomial", data = ade_peal)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.59576    0.04635   12.85   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2644.3  on 2031  degrees of freedom
## Residual deviance: 2644.3  on 2031  degrees of freedom
## AIC: 2646.3
## 
## Number of Fisher Scoring iterations: 4

exp(m0_alg.glm$coefficients["(Intercept)"])

## (Intercept) 
##    1.814404

table(ade_peal$cx)

## 
##  ade peal 
##  722 1310

1310/722

## [1] 1.814404

Kuna nullmudeli vabaliikme koefitsiendi p-väärtus on väike, peame nullhüpoteesi hülgama ning tõdema, et juhul, kui me ei tea mitte midagi kohakonstruktsioonide esinemiskonteksti kohta, ei ole peal- konstruktsiooni ja alalütleva käände konstruktsiooni kasutuse tõenäosused võrdsed.

Mitme seletava tunnusega mudel

Lisame nüüd mudelisse veel tunnuseid.

Mudelite tegemise puhul on kaks põhilist lähenemist:

eksploratiivne lähenemine, kus meil puuduvad tugevad teoreetilised eeldused efektide olulisuse, suuna ja suuruse kohta, püüame leida andmetest uusi seoseid ning optimaalset mudelit, mis ennustaks uuritava tunnuse varieerumist võimalikult hästi;
konfirmatoorne lähenemine, kus kontrollime varasemat teooriat või püstitatud hüpoteese ning püüame leida, millistel tunnustel on uuritava tunnuse varieerumisele meie valimi põhjal statistiliselt oluline mõju ja millistel mitte.

Konfirmatoorne lähenemine

Konfirmatoorse lähenemise puhul võime kõik oma andmestiku tunnused (ja nendevahelised interaktsioonid) panna korraga mudelisse ning raporteerida lihtsalt, millised neist olid olulised, millised mitte, ja mis suunalised oluliste tunnuste mõjud olid.

mmax.glm <- glm(cx ~ TR_liik + tegusõna_klass + LM_mobiilsus + LM_komplekssus + LM_pikkus_log + sõnajärg + murderühm, data = dat, family = "binomial")

summary(mmax.glm)

## 
## Call:
## glm(formula = cx ~ TR_liik + tegusõna_klass + LM_mobiilsus + 
##     LM_komplekssus + LM_pikkus_log + sõnajärg + murderühm, 
##     family = "binomial", data = dat)
## 
## Coefficients:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  0.63892    0.28548   2.238 0.025219 *  
## TR_liiknimisõnafraas        -0.23292    0.15150  -1.537 0.124194    
## TR_liikverbifraas           -0.31230    0.16484  -1.895 0.058145 .  
## tegusõna_klassolemisverb     0.58245    0.18262   3.189 0.001425 ** 
## tegusõna_klasspaiknemisverb  0.99535    0.36951   2.694 0.007066 ** 
## tegusõna_klasstegevusverb    0.99429    0.18914   5.257 1.46e-07 ***
## tegusõna_klassverbita        0.83611    0.25032   3.340 0.000837 ***
## LM_mobiilsusstaatiline      -1.24274    0.12971  -9.581  < 2e-16 ***
## LM_komplekssusliitsõna      -1.20587    0.21419  -5.630 1.80e-08 ***
## LM_pikkus_log                0.13668    0.09111   1.500 0.133566    
## sõnajärgtr_lm               -0.03367    0.13287  -0.253 0.799948    
## murderühmLõuna               0.26866    0.15637   1.718 0.085781 .  
## murderühmRanniku            -1.16511    0.18810  -6.194 5.86e-10 ***
## murderühmSaarte             -0.63889    0.15189  -4.206 2.60e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2001.8  on 1443  degrees of freedom
## Residual deviance: 1722.3  on 1430  degrees of freedom
## AIC: 1750.3
## 
## Number of Fisher Scoring iterations: 4

Baastasemete kontekst on selline, kus

TR_liik == “asesõna”,
tegusõna_klass == “liikumisverb”,
LM_mobiilsus == “liigutatav”,
LM_komplekssus == “lihtsõna”,
LM_pikkus_log == 0,
sõnajärg == “lm_tr”,
murderühm == “Kesk”.

Näiteks võiks baastasemete konteksti moodustada selline andmestikus tegelikult mitte-esinev lause, nagu käe peal kõnnib see, mille on öelnud Kesk-murderühma kõneleja.

Sellises kontekstis on peal-konstruktsiooni esinemise šansid suuremad kui selle mitte-esinemise (= käändekonstruktsiooni esinemise) šansid, sest (Intercept)-real on koefitsient > 0. Iga järgnev rida väljundis hindab nüüd seda, kas peal-konstruktsiooni esinemise šansid kasvavad või kahanevad, kui panna baastasemete konteksti selle mingi konkreetse tunnuse väärtus (nt TR_liik==“asesõna” -> TR_liik==“nimisõnafraas”), aga jätta ülejäänud kontekst samaks (nt käe peal kõnnib sipelgas).

Koondtabelist nähtub, et baastasemete kontekstiga võrreldes ei mõjuta peal-konstruktsiooni esinemise šansse oluliselt

trajektori (ehk selle, mis kuskil peal on) liigi muutmine mis tahes muuks liigiks,
orientiiri (ehk selle, mille peal miski on) silpide arvu (logaritmi) kasvatamine,
sõnajärje muutmine vastupidiseks (millegi peal on miski -> miski on millegi peal),
Kesk murderühma asemel Lõuna murderühma keelekasutuse vaatamine.

Baastasemete kontekstiga võrreldes kasvatab peal-konstruktsiooni esinemise šansse oluliselt

mis tahes teise tähendusklassi kuuluva tegusõna kasutamine või tegusõna väljajätt.

Baastasemete kontekstiga võrreldes kahandavad peal-konstruktsiooni esinemise šansse oluliselt

staatilise kohasõna kasutamine liigutatava kohasõna asemel (nt laud -> heinamaa),
liitsõnalise kohasõna kasutamine lihtsõnalise kohasõna asemel (nt laud -> heina+maa),
Ranniku või Saarte murderühma kuuluva keelekasutuse vaatamine.

library(report)
report(mmax.glm)

## We fitted a logistic model (estimated using ML) to predict cx with TR_liik,
## tegusõna_klass, LM_mobiilsus, LM_komplekssus, LM_pikkus_log, sõnajärg and
## murderühm (formula: cx ~ TR_liik + tegusõna_klass + LM_mobiilsus +
## LM_komplekssus + LM_pikkus_log + sõnajärg + murderühm). The model's explanatory
## power is moderate (Tjur's R2 = 0.18). The model's intercept, corresponding to
## TR_liik = asesõna, tegusõna_klass = liikumisverb, LM_mobiilsus = liigutatav,
## LM_komplekssus = lihtsõna, LM_pikkus_log = 0, sõnajärg = lm_tr and murderühm =
## Kesk, is at 0.64 (95% CI [0.08, 1.20], p = 0.025). Within this model:
## 
##   - The effect of TR liik [nimisõnafraas] is statistically non-significant and
## negative (beta = -0.23, 95% CI [-0.53, 0.06], p = 0.124; Std. beta = -0.23, 95%
## CI [-0.53, 0.06])
##   - The effect of TR liik [verbifraas] is statistically non-significant and
## negative (beta = -0.31, 95% CI [-0.64, 0.01], p = 0.058; Std. beta = -0.31, 95%
## CI [-0.64, 0.01])
##   - The effect of tegusõna klass [olemisverb] is statistically significant and
## positive (beta = 0.58, 95% CI [0.23, 0.94], p = 0.001; Std. beta = 0.58, 95% CI
## [0.23, 0.94])
##   - The effect of tegusõna klass [paiknemisverb] is statistically significant and
## positive (beta = 1.00, 95% CI [0.28, 1.73], p = 0.007; Std. beta = 1.00, 95% CI
## [0.28, 1.73])
##   - The effect of tegusõna klass [tegevusverb] is statistically significant and
## positive (beta = 0.99, 95% CI [0.63, 1.37], p < .001; Std. beta = 0.99, 95% CI
## [0.63, 1.37])
##   - The effect of tegusõna klass [verbita] is statistically significant and
## positive (beta = 0.84, 95% CI [0.35, 1.33], p < .001; Std. beta = 0.84, 95% CI
## [0.35, 1.33])
##   - The effect of LM mobiilsus [staatiline] is statistically significant and
## negative (beta = -1.24, 95% CI [-1.50, -0.99], p < .001; Std. beta = -1.24, 95%
## CI [-1.50, -0.99])
##   - The effect of LM komplekssus [liitsõna] is statistically significant and
## negative (beta = -1.21, 95% CI [-1.64, -0.79], p < .001; Std. beta = -1.21, 95%
## CI [-1.64, -0.79])
##   - The effect of LM pikkus log is statistically non-significant and positive
## (beta = 0.14, 95% CI [-0.04, 0.32], p = 0.134; Std. beta = 0.09, 95% CI [-0.03,
## 0.21])
##   - The effect of sõnajärg [tr_lm] is statistically non-significant and negative
## (beta = -0.03, 95% CI [-0.29, 0.23], p = 0.800; Std. beta = -0.03, 95% CI
## [-0.29, 0.23])
##   - The effect of murderühm [Lõuna] is statistically non-significant and positive
## (beta = 0.27, 95% CI [-0.04, 0.58], p = 0.086; Std. beta = 0.27, 95% CI [-0.04,
## 0.58])
##   - The effect of murderühm [Ranniku] is statistically significant and negative
## (beta = -1.17, 95% CI [-1.54, -0.80], p < .001; Std. beta = -1.17, 95% CI
## [-1.54, -0.80])
##   - The effect of murderühm [Saarte] is statistically significant and negative
## (beta = -0.64, 95% CI [-0.94, -0.34], p < .001; Std. beta = -0.64, 95% CI
## [-0.94, -0.34])
## 
## 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.

Kuna aga iga taseme mõju hinnatakse ainult üldise baastasemete konteksti suhtes, ei saa me siit teada, kas tunnustel üldiselt on mingi mõju või mitte. Samuti ei saa me mudeli väljundist teada 3 või enama tasemega kategoriaalsete seletavate tunnuste omavahelisi erinevusi.

Kui kasutaksime nüüd tunnuste mõju üldiseks hindamiseks lihtsalt anova() käsku, võime saada natuke erinevaid tulemusi vastavalt sellele, mis järjekorras tunnuseid mudelisse lisaksime, sest see funktsioon teeb nn I tüüpi ANOVA testi.

anova(mmax.glm, test = "Chisq")

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: cx
## 
## Terms added sequentially (first to last)
## 
## 
##                Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                            1443     2001.8              
## TR_liik         2    5.572      1441     1996.2   0.06167 .  
## tegusõna_klass  4   67.162      1437     1929.1 9.010e-14 ***
## LM_mobiilsus    1  102.663      1436     1826.4 < 2.2e-16 ***
## LM_komplekssus  1   33.561      1435     1792.8 6.906e-09 ***
## LM_pikkus_log   1    2.344      1434     1790.5   0.12575    
## sõnajärg        1    0.320      1433     1790.2   0.57141    
## murderühm       3   67.894      1430     1722.3 1.205e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Muudame järjekorda
anova(glm(cx ~ murderühm + sõnajärg + LM_pikkus_log + LM_komplekssus + LM_mobiilsus + tegusõna_klass + TR_liik, data = dat, family = "binomial"), test = "Chisq")

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: cx
## 
## Terms added sequentially (first to last)
## 
## 
##                Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                            1443     2001.8              
## murderühm       3   81.217      1440     1920.6 < 2.2e-16 ***
## sõnajärg        1    0.229      1439     1920.4    0.6323    
## LM_pikkus_log   1    0.111      1438     1920.2    0.7388    
## LM_komplekssus  1   52.248      1437     1868.0 4.891e-13 ***
## LM_mobiilsus    1  111.784      1436     1756.2 < 2.2e-16 ***
## tegusõna_klass  4   30.015      1432     1726.2 4.861e-06 ***
## TR_liik         2    3.913      1430     1722.3    0.1413    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Muudame järjekorda
anova(glm(cx ~ LM_pikkus_log + sõnajärg + TR_liik + LM_komplekssus + murderühm + tegusõna_klass + LM_mobiilsus, data = dat, family = "binomial"), test = "Chisq")

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: cx
## 
## Terms added sequentially (first to last)
## 
## 
##                Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                            1443     2001.8              
## LM_pikkus_log   1    0.077      1442     2001.7   0.78178    
## sõnajärg        1    0.013      1441     2001.7   0.90776    
## TR_liik         2    5.631      1439     1996.1   0.05989 .  
## LM_komplekssus  1   51.853      1438     1944.2 5.982e-13 ***
## murderühm       3   78.656      1435     1865.6 < 2.2e-16 ***
## tegusõna_klass  4   45.043      1431     1820.5 3.896e-09 ***
## LM_mobiilsus    1   98.244      1430     1722.3 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Kui tahaksime nt määrata ka, milline tunnus on kõige olulisem ja kasutada selleks Deviance tulpa, mis näitab, kui palju vähendab mingi tunnus seletamata varieerumist uuritavas tunnuses, siis sõltub tulemus sellest, mis järjekorras oleme tunnuseid mudelisse lisanud.

Seepärast on konfirmatoorses analüüsis parem kasutada II tüüpi ANOVAt, mis võrdleb iga tunnuse mõju koos kõikide teiste lisatavate tunnustega. Selleks on funktsioon Anova() paketist car või funktsioon drop1. III tüüpi ANOVAt, millest on varem ka juttu olnud, on kasulikum kasutada siis, kui mudelis on ka seletavate tunnuste vahelisi interaktsioone ehk koosmõjusid.

car::Anova(mmax.glm)

## Analysis of Deviance Table (Type II tests)
## 
## Response: cx
##                LR Chisq Df Pr(>Chisq)    
## TR_liik           3.913  2     0.1413    
## tegusõna_klass   31.393  4  2.545e-06 ***
## LM_mobiilsus     98.244  1  < 2.2e-16 ***
## LM_komplekssus   35.015  1  3.272e-09 ***
## LM_pikkus_log     2.257  1     0.1330    
## sõnajärg          0.064  1     0.7999    
## murderühm        67.894  3  1.205e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# või
drop1(mmax.glm, test = "Chisq")

## Single term deletions
## 
## Model:
## cx ~ TR_liik + tegusõna_klass + LM_mobiilsus + LM_komplekssus + 
##     LM_pikkus_log + sõnajärg + murderühm
##                Df Deviance    AIC    LRT  Pr(>Chi)    
## <none>              1722.3 1750.3                     
## TR_liik         2   1726.2 1750.2  3.913    0.1413    
## tegusõna_klass  4   1753.7 1773.7 31.393 2.545e-06 ***
## LM_mobiilsus    1   1820.5 1846.5 98.244 < 2.2e-16 ***
## LM_komplekssus  1   1757.3 1783.3 35.015 3.272e-09 ***
## LM_pikkus_log   1   1724.5 1750.5  2.257    0.1330    
## sõnajärg        1   1722.4 1748.4  0.064    0.7999    
## murderühm       3   1790.2 1812.2 67.894 1.205e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Võrdle ka
car::Anova(mmax.glm)
car::Anova(glm(cx ~ murderühm + sõnajärg + LM_pikkus_log + LM_komplekssus + LM_mobiilsus + tegusõna_klass + TR_liik, data = dat, family = "binomial"))
car::Anova(glm(cx ~ LM_pikkus_log + sõnajärg + TR_liik + LM_komplekssus + murderühm + tegusõna_klass + LM_mobiilsus, data = dat, family = "binomial"))

Tulba LRT ehk likelihood-ratio test (või LR Chisq) põhjal saab nüüd otsustada, et kõige olulisemad tunnused peal-konstruktsiooni ja käändekonstruktsiooni valikul on kohasõna liigutatavus, murderühm, kohasõna komplekssus ning tegusõna klass.

3 või enama tasemega faktorite tasemete omavaheliseks võrdluseks võib teha nüüd samuti mitu mudelit, kus igal pool muudame faktori baastaset, ning jagada p-väärtused läbi võrdluste arvuga (nn Bonferroni parandus).
Baastasemeteks võib aga määrata ka mingid tunnuste tüüpilised väärtused, nt iga faktori kõige sagedama väärtuse. Nii hindame mudelis parameetreid mingis oletatavasti kõige tüüpilisemas kontekstis.

# Mis tase igas kategoriaalses seletavas tunnuses on kõige sagedam?
lapply(dat[,seletavad_kat], table)

## $TR_liik
## 
##       asesõna nimisõnafraas    verbifraas 
##           351           613           480 
## 
## $tegusõna_klass
## 
##  liikumisverb    olemisverb paiknemisverb   tegevusverb       verbita 
##           245           620            44           407           128 
## 
## $LM_mobiilsus
## 
## liigutatav staatiline 
##        469        975 
## 
## $LM_komplekssus
## 
## lihtsõna liitsõna 
##     1285      159 
## 
## $sõnajärg
## 
## lm_tr tr_lm 
##   430  1014 
## 
## $murderühm
## 
##    Kesk   Lõuna Ranniku  Saarte 
##     611     294     215     324

# Muudame faktorite tasemeid ja salvestame uue andmestiku eraldi objekti dat2
dat %>%
  mutate(TR_liik = relevel(TR_liik, ref = "nimisõnafraas"),
         tegusõna_klass = relevel(tegusõna_klass, ref = "olemisverb"),
         LM_mobiilsus = relevel(LM_mobiilsus, ref = "staatiline"),
         sõnajärg = relevel(sõnajärg, ref = "tr_lm")) -> dat2

mmax2.glm <- glm(cx ~ TR_liik + tegusõna_klass + LM_mobiilsus + LM_komplekssus + LM_pikkus_log + sõnajärg + murderühm, data = dat2, family = "binomial")

summary(mmax2.glm)

## 
## Call:
## glm(formula = cx ~ TR_liik + tegusõna_klass + LM_mobiilsus + 
##     LM_komplekssus + LM_pikkus_log + sõnajärg + murderühm, 
##     family = "binomial", data = dat2)
## 
## Coefficients:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                 -0.28796    0.17574  -1.639  0.10131    
## TR_liikasesõna               0.23292    0.15150   1.537  0.12419    
## TR_liikverbifraas           -0.07938    0.14688  -0.540  0.58888    
## tegusõna_klassliikumisverb  -0.58245    0.18262  -3.189  0.00143 ** 
## tegusõna_klasspaiknemisverb  0.41290    0.35092   1.177  0.23934    
## tegusõna_klasstegevusverb    0.41184    0.14837   2.776  0.00551 ** 
## tegusõna_klassverbita        0.25366    0.21372   1.187  0.23527    
## LM_mobiilsusliigutatav       1.24274    0.12971   9.581  < 2e-16 ***
## LM_komplekssusliitsõna      -1.20587    0.21419  -5.630 1.80e-08 ***
## LM_pikkus_log                0.13668    0.09111   1.500  0.13357    
## sõnajärglm_tr                0.03367    0.13287   0.253  0.79995    
## murderühmLõuna               0.26866    0.15637   1.718  0.08578 .  
## murderühmRanniku            -1.16511    0.18810  -6.194 5.86e-10 ***
## murderühmSaarte             -0.63889    0.15189  -4.206 2.60e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2001.8  on 1443  degrees of freedom
## Residual deviance: 1722.3  on 1430  degrees of freedom
## AIC: 1750.3
## 
## Number of Fisher Scoring iterations: 4

car::Anova(mmax2.glm)

## Analysis of Deviance Table (Type II tests)
## 
## Response: cx
##                LR Chisq Df Pr(>Chisq)    
## TR_liik           3.913  2     0.1413    
## tegusõna_klass   31.393  4  2.545e-06 ***
## LM_mobiilsus     98.244  1  < 2.2e-16 ***
## LM_komplekssus   35.015  1  3.272e-09 ***
## LM_pikkus_log     2.257  1     0.1330    
## sõnajärg          0.064  1     0.7999    
## murderühm        67.894  3  1.205e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Nüüd näeme, et ehkki eelmises mudelis erinesid liikumisverbid kõikidest teistest verbitüüpidest, siis näiteks paiknemisverbid olemisverbidest peal-konstruktsiooni kasutamise tõenäosuses ei erine.

Eksploratiivne lähenemine

Eksploratiivse lähenemise puhul, kus sooviksime andmetest alles seoseid leida ning selle põhjalt hüpoteese luua, tahaksime jõuda võimalikult lihtsa, aga maksimaalse ennustusjõuga mudelini. Selleks võime kasutada inkrementaalset mudeliehitamist, mida juba lineaarse regressiooni puhul vaatasime.

Logistilise regressiooni mudeleid ei saa võrrelda R-ruudu alusel nagu lineaarses regressioonis, aga saab võrrelda Akaike informatsioonikriteeriumit AIC, mis väljendab seletamata varieerumise hulka andmetes. Teisest perspektiivist võib AIC-d kirjeldada ka kui seda informatsiooni hulka, mis läheb kaotsi, kui algandmete varieeruvust mudeli kaudu seletada: mudel ei ole kunagi täiesti täpne representatsioon algandmetest. AIC arvestab valimi suuruse ja mudelisse kaasatud parameetrite arvuga (nagu kohandatud R-ruut lineaarse regressiooni mudelis) ja mida väiksem AIC väärtus on, seda parem. AIC-d võib võrrelda erinevate mudelite vahel, mis põhinevad samal andmestikul ja millel on sama uuritav tunnus, ilma et nendes mudelites peaks tingimata olema jagatud seletavaid tunnuseid.

Nn nested ehk hierarhiliste mudelite võrdlemiseks, kus igas komplekssemas mudelis sisalduvad ka lihtsamates mudelites olevad mõjud, saab kasutada ka logistilise regressiooni mudelite puhul anova()-funktsiooni, mis teeb mudelite jääkidel esimest tüüpi ANOVA. Kui tahame väljundis näha ka jääkide erinevuse statistilist olulisust, tuleb logistiliste mudelite puhul täpsustada ka, et testi liik on hii-ruut-test (test = "Chisq").

m0.glm <- glm(cx ~ 1, data = dat, family = "binomial")

m1.glm <- glm(cx ~ LM_mobiilsus, data = dat, family = "binomial")
anova(m0.glm, m1.glm, test = "Chisq") # teeb paremaks

## Analysis of Deviance Table
## 
## Model 1: cx ~ 1
## Model 2: cx ~ LM_mobiilsus
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1      1443     2001.8                          
## 2      1442     1879.0  1   122.86 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Tulbad Resid.Dev ja Resid.Df on vastavalt jääkhälbimuse näitaja ja sellega seotud vabadusastmete arv, mida näeb ka mudelist summary() käsuga. Df on vabadusastmete arv mudelis, mis arvuliste tunnuste puhul on alati 1 ja kategoriaalsete tunnuste puhul tasemete arv -1. Deviance näitab siin kahe järjestikuse mudeli jääkhälbimuste ümardatud vahet, 122.86 on 2001.8-1879.0. Mida suurem on mudeli Deviance, seda rohkem varieerumisest see mudel seletab, võrreldes sellele eelneva mudeliga.

m2.glm <- glm(cx ~ LM_mobiilsus + LM_komplekssus, data = dat, family = "binomial")
anova(m1.glm, m2.glm, test = "Chisq") # teeb paremaks

## Analysis of Deviance Table
## 
## Model 1: cx ~ LM_mobiilsus
## Model 2: cx ~ LM_mobiilsus + LM_komplekssus
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1      1442     1879.0                          
## 2      1441     1840.8  1   38.178 6.459e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m3.glm <- glm(cx ~ LM_mobiilsus + LM_komplekssus + tegusõna_klass, data = dat, family = "binomial")
anova(m2.glm, m3.glm, test = "Chisq") # teeb paremaks

## Analysis of Deviance Table
## 
## Model 1: cx ~ LM_mobiilsus + LM_komplekssus
## Model 2: cx ~ LM_mobiilsus + LM_komplekssus + tegusõna_klass
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1      1441     1840.8                          
## 2      1437     1796.8  4   44.027 6.332e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m4.glm <- glm(cx ~ LM_mobiilsus + LM_komplekssus + tegusõna_klass + murderühm, data = dat, family = "binomial")
anova(m3.glm, m4.glm, test = "Chisq") # teeb paremaks

## Analysis of Deviance Table
## 
## Model 1: cx ~ LM_mobiilsus + LM_komplekssus + tegusõna_klass
## Model 2: cx ~ LM_mobiilsus + LM_komplekssus + tegusõna_klass + murderühm
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1      1437     1796.8                          
## 2      1434     1728.8  3   67.993 1.148e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m5.glm <- glm(cx ~ LM_mobiilsus + LM_komplekssus + tegusõna_klass + murderühm + TR_liik, data = dat, family = "binomial")
anova(m4.glm, m5.glm, test = "Chisq") # ei tee paremaks

## Analysis of Deviance Table
## 
## Model 1: cx ~ LM_mobiilsus + LM_komplekssus + tegusõna_klass + murderühm
## Model 2: cx ~ LM_mobiilsus + LM_komplekssus + tegusõna_klass + murderühm + 
##     TR_liik
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1      1434     1728.8                     
## 2      1432     1724.6  2   4.1175   0.1276

Ja nii edasi.

# Võrdleme ka Akaike informatsioonikriteeriumit
AIC(m0.glm, m1.glm, m2.glm, m3.glm, m4.glm, m5.glm)

##        df      AIC
## m0.glm  1 2003.809
## m1.glm  2 1882.954
## m2.glm  3 1846.776
## m3.glm  7 1810.749
## m4.glm 10 1748.756
## m5.glm 12 1748.638

m5.glm-mudeli AIC on kõige väiksem, aga arvestades kulu, mis läheb trajektori liigi mõju hindamiseks, ei tasu selle tunnuse lisamine ning mudeli komplekssemaks tegemine ära ning peaks eelistama lihtsamat, ilma trajektori liigita mudelit.

anova() ei ütle iseenesest millist mudelit peaksime eelistama, vaid lihtsalt seda, kas erinevus kahe või enama mudeli vahel on statistiliselt oluline, st seda, kas väiksem jääkhälbimus (Resid. Dev, mis on lähedane AIC-ga) võib olla saadud juhuslikult või mitte.

Saaksime kasutada ka logistilise regressiooni puhul step-protseduuri, kus laseme R-il meie eest inkrementaalse modelleerimise töö ära teha.

mopt.glm <- step(glm(cx ~ 1, data = dat, family = "binomial"), scope = cx ~ TR_liik * tegusõna_klass * LM_mobiilsus * LM_komplekssus * LM_pikkus_log * sõnajärg * murderühm, direction = "forward")

## Start:  AIC=2003.81
## cx ~ 1
## 
##                  Df Deviance    AIC
## + LM_mobiilsus    1   1879.0 1883.0
## + murderühm       3   1920.6 1928.6
## + tegusõna_klass  4   1935.3 1945.3
## + LM_komplekssus  1   1953.4 1957.4
## + TR_liik         2   1996.2 2002.2
## <none>                2001.8 2003.8
## + LM_pikkus_log   1   2001.7 2005.7
## + sõnajärg        1   2001.8 2005.8
## 
## Step:  AIC=1882.95
## cx ~ LM_mobiilsus
## 
##                  Df Deviance    AIC
## + murderühm       3   1794.9 1804.9
## + tegusõna_klass  4   1831.0 1843.0
## + LM_komplekssus  1   1840.8 1846.8
## <none>                1879.0 1883.0
## + TR_liik         2   1875.8 1883.8
## + sõnajärg        1   1878.5 1884.5
## + LM_pikkus_log   1   1878.9 1884.9
## 
## Step:  AIC=1804.88
## cx ~ LM_mobiilsus + murderühm
## 
##                          Df Deviance    AIC
## + LM_komplekssus          1   1758.5 1770.5
## + tegusõna_klass          4   1762.1 1780.1
## <none>                        1794.9 1804.9
## + TR_liik                 2   1792.0 1806.0
## + sõnajärg                1   1794.8 1806.8
## + LM_pikkus_log           1   1794.8 1806.8
## + LM_mobiilsus:murderühm  3   1791.2 1807.2
## 
## Step:  AIC=1770.5
## cx ~ LM_mobiilsus + murderühm + LM_komplekssus
## 
##                               Df Deviance    AIC
## + tegusõna_klass               4   1728.8 1748.8
## + LM_pikkus_log                1   1756.2 1770.2
## <none>                             1758.5 1770.5
## + LM_komplekssus:murderühm     3   1754.0 1772.0
## + TR_liik                      2   1756.0 1772.0
## + LM_mobiilsus:LM_komplekssus  1   1758.5 1772.5
## + sõnajärg                     1   1758.5 1772.5
## + LM_mobiilsus:murderühm       3   1754.9 1772.9
## 
## Step:  AIC=1748.76
## cx ~ LM_mobiilsus + murderühm + LM_komplekssus + tegusõna_klass
## 
##                                 Df Deviance    AIC
## + LM_pikkus_log                  1   1726.2 1748.2
## + TR_liik                        2   1724.6 1748.6
## <none>                               1728.8 1748.8
## + tegusõna_klass:LM_mobiilsus    4   1720.9 1748.9
## + LM_komplekssus:murderühm       3   1724.2 1750.2
## + LM_mobiilsus:LM_komplekssus    1   1728.7 1750.7
## + sõnajärg                       1   1728.7 1750.7
## + LM_mobiilsus:murderühm         3   1725.7 1751.7
## + tegusõna_klass:murderühm      12   1709.5 1753.5
## + tegusõna_klass:LM_komplekssus  4   1726.0 1754.0
## 
## Step:  AIC=1748.24
## cx ~ LM_mobiilsus + murderühm + LM_komplekssus + tegusõna_klass + 
##     LM_pikkus_log
## 
##                                 Df Deviance    AIC
## + LM_pikkus_log:murderühm        3   1714.4 1742.4
## <none>                               1726.2 1748.2
## + LM_komplekssus:LM_pikkus_log   1   1724.3 1748.3
## + TR_liik                        2   1722.4 1748.4
## + tegusõna_klass:LM_mobiilsus    4   1718.9 1748.9
## + LM_komplekssus:murderühm       3   1721.8 1749.8
## + tegusõna_klass:LM_pikkus_log   4   1720.0 1750.0
## + LM_mobiilsus:LM_pikkus_log     1   1726.0 1750.0
## + LM_mobiilsus:LM_komplekssus    1   1726.2 1750.2
## + sõnajärg                       1   1726.2 1750.2
## + LM_mobiilsus:murderühm         3   1723.5 1751.5
## + tegusõna_klass:LM_komplekssus  4   1723.5 1753.5
## + tegusõna_klass:murderühm      12   1707.6 1753.6
## 
## Step:  AIC=1742.43
## cx ~ LM_mobiilsus + murderühm + LM_komplekssus + tegusõna_klass + 
##     LM_pikkus_log + murderühm:LM_pikkus_log
## 
##                                 Df Deviance    AIC
## + LM_komplekssus:LM_pikkus_log   1   1711.8 1741.8
## + TR_liik                        2   1710.3 1742.3
## + tegusõna_klass:LM_mobiilsus    4   1706.4 1742.4
## <none>                               1714.4 1742.4
## + LM_mobiilsus:LM_pikkus_log     1   1714.0 1744.0
## + LM_mobiilsus:LM_komplekssus    1   1714.4 1744.4
## + sõnajärg                       1   1714.4 1744.4
## + tegusõna_klass:LM_pikkus_log   4   1708.4 1744.4
## + LM_mobiilsus:murderühm         3   1712.0 1746.0
## + LM_komplekssus:murderühm       3   1712.2 1746.2
## + tegusõna_klass:murderühm      12   1695.5 1747.5
## + tegusõna_klass:LM_komplekssus  4   1711.9 1747.9
## 
## Step:  AIC=1741.8
## cx ~ LM_mobiilsus + murderühm + LM_komplekssus + tegusõna_klass + 
##     LM_pikkus_log + murderühm:LM_pikkus_log + LM_komplekssus:LM_pikkus_log
## 
##                                 Df Deviance    AIC
## + tegusõna_klass:LM_mobiilsus    4   1703.7 1741.7
## + TR_liik                        2   1707.7 1741.7
## <none>                               1711.8 1741.8
## + LM_mobiilsus:LM_pikkus_log     1   1711.1 1743.1
## + LM_mobiilsus:LM_komplekssus    1   1711.7 1743.7
## + tegusõna_klass:LM_pikkus_log   4   1705.8 1743.8
## + sõnajärg                       1   1711.8 1743.8
## + LM_mobiilsus:murderühm         3   1709.4 1745.4
## + LM_komplekssus:murderühm       3   1709.6 1745.6
## + tegusõna_klass:LM_komplekssus  4   1709.1 1747.1
## + tegusõna_klass:murderühm      12   1693.2 1747.2
## 
## Step:  AIC=1741.65
## cx ~ LM_mobiilsus + murderühm + LM_komplekssus + tegusõna_klass + 
##     LM_pikkus_log + murderühm:LM_pikkus_log + LM_komplekssus:LM_pikkus_log + 
##     LM_mobiilsus:tegusõna_klass
## 
##                                 Df Deviance    AIC
## <none>                               1703.7 1741.7
## + tegusõna_klass:LM_pikkus_log   4   1695.9 1741.9
## + TR_liik                        2   1700.0 1742.0
## + LM_mobiilsus:LM_pikkus_log     1   1702.7 1742.7
## + LM_mobiilsus:LM_komplekssus    1   1703.3 1743.3
## + sõnajärg                       1   1703.6 1743.6
## + LM_komplekssus:murderühm       3   1701.0 1745.0
## + LM_mobiilsus:murderühm         3   1701.1 1745.1
## + tegusõna_klass:LM_komplekssus  4   1701.0 1747.0
## + tegusõna_klass:murderühm      12   1685.9 1747.9

car::Anova(mopt.glm)

## Analysis of Deviance Table (Type II tests)
## 
## Response: cx
##                              LR Chisq Df Pr(>Chisq)    
## LM_mobiilsus                   89.725  1  < 2.2e-16 ***
## murderühm                      67.203  3  1.694e-14 ***
## LM_komplekssus                 36.433  1  1.580e-09 ***
## tegusõna_klass                 30.980  4  3.090e-06 ***
## LM_pikkus_log                   1.990  1   0.158358    
## murderühm:LM_pikkus_log        13.214  3   0.004196 ** 
## LM_komplekssus:LM_pikkus_log    2.777  1   0.095645 .  
## LM_mobiilsus:tegusõna_klass     8.151  4   0.086213 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Tundub, et orientiiri pikkus võiks mingit mõju omada koos murderühmaga. Vaatame seda natuke hiljem.

Mudeli mõjude visualiseerimine

On oluline, et oskaksime regressioonimudeli väljundeid lugeda ja tõlgendada. Tihtipeale on aga oma tulemuste teistele esitamiseks mugavam kasutada eri visualiseerimisviise.

Vaatame siin paketti sjPlot, mis on loodud sotsiaalteaduste statistiliste analüüside visualiseerimiseks. Pakett kasutab tegelikult sama ggeffects paketi funktsioone, mida eelmistel kordadel oleme samuti põgusalt vaadanud.

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

Funktsioon tab_model() kuvab mudeli väljundi html-ina, ilusti vormistatud tabelis.

tab_model(mmax.glm, use.viewer = F)

	cx
Predictors	Odds Ratios	CI	p
(Intercept)	1.89	1.08 – 3.32	0.025
TR liik [nimisõnafraas]	0.79	0.59 – 1.07	0.124
TR liik [verbifraas]	0.73	0.53 – 1.01	0.058
tegusõna klass [olemisverb]	1.79	1.26 – 2.57	0.001
tegusõna klass [paiknemisverb]	2.71	1.32 – 5.64	0.007
tegusõna klass [tegevusverb]	2.70	1.87 – 3.93	<0.001
tegusõna klass [verbita]	2.31	1.42 – 3.78	0.001
LM mobiilsus [staatiline]	0.29	0.22 – 0.37	<0.001
LM komplekssus [liitsõna]	0.30	0.19 – 0.45	<0.001
LM pikkus log	1.15	0.96 – 1.37	0.134
sõnajärg [tr_lm]	0.97	0.75 – 1.25	0.800
murderühm [Lõuna]	1.31	0.96 – 1.78	0.086
murderühm [Ranniku]	0.31	0.21 – 0.45	<0.001
murderühm [Saarte]	0.53	0.39 – 0.71	<0.001
Observations	1444
R² Tjur	0.178

Vaikimisi on logistilises regressioonis kõik koefitsiendid tabelis juba uuesti šansisuheteks teisendatud. Tabelit saab aga kõiksugu viisidel kohandada, näidata tavalisi mudeli väljundeid, lisada näidatavaid tulpasid, järjestada neid ümber jne.

tab_model(mmax.glm, 
          show.reflvl = T, # näita tunnuste baastasemeid
          show.ci = F, # ära näita usaldusvahemikke
          show.se = T, # näita koefitsiendi standardviga
          show.aic = T, # näita mudeli AIC-d
          title = "Adessiivi ja peal-kaassõna kasutust mõjutavad tegurid Eesti murretes", # tabeli pealkiri
          string.pred = "Tunnused", # tunnuste tulbanimi
          string.est = "Šansisuhe", # koefitsientide tulbanimi
          string.se = "Standardviga", # standardvea tulbanimi
          string.intercept = "Vabaliige", # vabaliikme nimi
          p.adjust = "holm", # Holmi-Bonferroni parandustega p-väärtused
          p.style = "numeric_stars", # näita ka stat.olulisuse tärnikesi
          use.viewer = F) # ära näita tabelit Vieweris, vaid eraldi brauseris

Adessiivi ja peal-kaassõna kasutust mõjutavad tegurid Eesti murretes
	cx
Tunnused	Šansisuhe	Standardviga	p
asesõna	Reference
nimisõnafraas	0.79	0.12	0.373
verbifraas	0.73	0.12	0.291
liikumisverb	Reference
LM_pikkus_log	1.15	0.10	0.373
olemisverb	1.79 ^*	0.33	0.011
paiknemisverb	2.71 ^*	1.00	0.049
tegevusverb	2.70 ^***	0.51	<0.001
verbita	2.31 ^**	0.58	0.008
liigutatav	Reference
staatiline	0.29 ^***	0.04	<0.001
lihtsõna	Reference
liitsõna	0.30 ^***	0.06	<0.001
lm_tr	Reference
tr_lm	0.97	0.13	0.800
Kesk	Reference
Lõuna	1.31	0.20	0.343
Ranniku	0.31 ^***	0.06	<0.001
Saarte	0.53 ^***	0.08	<0.001
Vabaliige	1.89	0.54	0.151
Observations	1444
R² Tjur	0.178
AIC	1750.292
p<0.05 p<0.01 * p<0.001

Windowsi operatsioonisüsteemiga võib olla vajalik täpitähtede kuvamiseks määrata kodeeringuks encoding = "windows-1257". Muudes operatsioonisüsteemides peaks toimima vaikimisi säte encoding = "UTF-8".
Tabeli saab salvestada ka näiteks eraldi docx-faili, kui lisada funktsiooni veel argument file = "failinimi.docx".

Mudeli väljundite kuvamiseks joonisel on aga funktsioon plot_model().

plot_model(mmax.glm)

Vaikimisi kuvatakse jälle juba teisendatud šansisuhteid (vahemikus 0 kuni lõpmatus), aga võib näidata ka logaritmitud šansisuhteid (transform = NULL) või tõenäosuseid (transform = "plogis"). Tunnused, mille usaldusvahemik ei kata 1, on statistiliselt olulised ning mõju suunda näitab see, kas punkt on sinine ning 1st suurem või punane ja 1st väiksem.

Ka joonisel võime aga igasugu asju muuta.

plot_model(mmax.glm, 
           sort.est = T, # sorteeri koefitsiendid väärtuste järgi
           show.values = T, # näita koefitsientide väärtusi
           value.offset = 0.3, # tõsta neid natuke ülespoole
           title = "peal-konstruktsiooni šansside muutus võrreldes baastasemete kontekstiga", # lisa pealkiri
           p.adjust = "holm") # Holmi-Bonferroni parandustega p-väärtused

Samuti võime kuvada ainult kindlaid tunnuseid.

plot_model(mmax.glm, 
           terms = c("tegusõna_klassolemisverb",
                     "tegusõna_klasspaiknemisverb",
                     "tegusõna_klasstegevusverb", 
                     "tegusõna_klassverbita"))

Lisaks üldiste mudeli väljundite kuvamisele saab joonisel näidata ka eraldi konkreetsete tunnuste mõjusid (nn marginal effects). type = "pred" kuvab näiteks peamõjusid, kusjuures y-teljel on nüüd protsentideks arvutatud esinemistõenäosused teiste seletavate tunnuste baastasemete kontekstis (aga arvuliste tunnuste keskväärtuse korral!). Võid kontrollida konteksti käsuga ggeffects::ggpredict(mmax.glm, terms = "tegusõna_klass").

plot_model(mmax.glm, 
           type = "pred", 
           terms = "tegusõna_klass", 
           show.data = T, # lisa ka tegelikud andmepunktid
           jitter = 0.1) # hajuta neid pisut

Vaatame ka meie ainsat arvulist seletavat tunnust ja kasutame ggplot-paketi omadusi, et graafiku kujundust muuta.

plot_model(mmax.glm,
           type = "pred",
           terms = "LM_pikkus_log",
           show.data = T,
           jitter = 0.1) +
  theme(plot.title = element_text(hjust = 0.5)) # pealkiri joondatud keskele

Interaktsioonidega mudel

Nagu lineaarses regressioonis, võib ka logistilise regressiooni mudelis hinnata tunnustevahelisi interaktsioone. See tähendab, et ühe seletava tunnuse mõju sõltub teise seletava tunnuse väärtustest. Näiteks hüpoteetiliselt:

Gripihooajal grippi haigestumise (vs. mittehaigestumise) tõenäosust vähendab vaktsineeritud olemine ainult juhul, kui vaktsineeritud on vahemikus oktoobrist detsembrini. Muul juhul statistiliselt oluline mõju puudub;
Teietamise (vs. sinatamise) tõenäosus kasvab vastavalt vestluspartneri vanuse tõusule ainult juhul, kui vestluspartner on võõras. Tuttavate vestluspartneritega vanusel statistiliselt oluline mõju puudub;
Leksikaalses otsustuskatses sõna tõenäosus saada klassifitseerituks päriselt olemasolevaks sõnaks (vs. päriselt mitte olemasolevaks sõnaks) kasvab vastavalt täishäälikute hulga suurenemisele sõnas juhul, kui sõna pikkus on alla 8 tähemärgi, ning kahaneb juhul, kui sõna pikkus on üle 8 tähemärgi.

Vaatame nüüd, kas võiksime oma mudelisse panna ka mõne interaktsiooni, näiteks orientiiri pikkuse ja murderühma vahel, mida enne step-meetodiga avastasime.

ggplot(dat, aes(x = LM_pikkus_log, y = ifelse(cx == "peal", 1, 0))) +
  geom_point(alpha = 0.1) +
  geom_smooth(method = "glm", method.args = list(family = "binomial")) +
  facet_wrap("murderühm")

## `geom_smooth()` using formula = 'y ~ x'

Tundub, et Saarte ja Lõuna murderühmas võiks pikkusel olla positiivne mõju, Kesk murderühmas negatiivne ja Ranniku rühmas üldse mitte mingit mõju.

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

AIC(m4.glm, m5.glm)

##        df      AIC
## m4.glm 10 1748.756
## m5.glm 14 1742.431

anova(m4.glm, m5.glm, test = "Chisq") # teeb paremaks

## Analysis of Deviance Table
## 
## Model 1: cx ~ LM_mobiilsus + LM_komplekssus + tegusõna_klass + murderühm
## Model 2: cx ~ LM_mobiilsus + LM_komplekssus + tegusõna_klass + LM_pikkus_log * 
##     murderühm
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)   
## 1      1434     1728.8                        
## 2      1430     1714.4  4   14.324 0.006329 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova(m5.glm, type = "III")

## Analysis of Deviance Table (Type III tests)
## 
## Response: cx
##                         LR Chisq Df Pr(>Chisq)    
## LM_mobiilsus              93.505  1  < 2.2e-16 ***
## LM_komplekssus            35.596  1  2.427e-09 ***
## tegusõna_klass            31.058  4  2.979e-06 ***
## LM_pikkus_log              1.094  1   0.295560    
## murderühm                 24.811  3  1.691e-05 ***
## LM_pikkus_log:murderühm   11.811  3   0.008058 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Kui mudelis on olulisi interaktsioone, võiks kasutada III tüüpi ANOVAt. Loe lähemalt siit.

tab_model(m5.glm, 
          transform = NULL, # näita logaritmitud šansisuhteid
          show.reflvl = T, # näita tunnuste baastasemeid
          show.ci = F, # ära näita usaldusvahemikke
          show.se = T, # näita koefitsiendi standardviga
          show.aic = T, # näita mudeli AIC-d
          use.viewer = F) # ära näita tabelit Vieweris, vaid eraldi brauseris

	cx
Predictors	Log-Odds	std. Error	p
(Intercept)	0.75	0.26	0.004
liigutatav	Reference
staatiline	-1.22	0.13	<0.001
lihtsõna	Reference
liitsõna	-1.21	0.21	<0.001
LM_pikkus_log	-0.14	0.13	0.296
LM_pikkus_log:murderühmLõuna	0.61	0.23	0.008
LM_pikkus_log:murderühmRanniku	0.12	0.30	0.688
LM_pikkus_log:murderühmSaarte	0.65	0.23	0.006
liikumisverb	Reference
olemisverb	0.60	0.18	0.001
paiknemisverb	0.95	0.37	0.010
tegevusverb	0.99	0.19	<0.001
verbita	0.87	0.25	<0.001
Kesk	Reference
Lõuna	-0.49	0.32	0.124
Ranniku	-1.34	0.41	0.001
Saarte	-1.45	0.33	<0.001
Observations	1444
R² Tjur	0.183
AIC	1742.431

Kui mudelis on mingite tunnuste vahel oluline interaktsioon, ei saa mudeli peamõjusid enam üksinda tõlgendada. Nii tuleks pikkuse mõju tõlgendada mudelis m5.glm ainult koos murde mõjuga. Pikkusel (LM_pikkus_log) ei ole niisiis olulist mõju, juhul kui tegemist on Kesk-murderühmaga. Lõuna-murderühm (murderühmLõuna) ei erine oluliselt Kesk-murderühmast, juhul kui pikkuse logaritm on 0, aga silpide lisandudes muutub erinevus oluliseks (LM_pikkus_log:murderühmLõuna). Ranniku-murderühmaga (murderühmRanniku) on vastupidi: erinevus Kesk-murderühmaga on statistiliselt oluline juhul kui pikkuse logaritm on 0, aga silpide arvu kasvades midagi oluliselt ei muutu. Saarte-murderühm omakorda erineb Kesk-murderühmast mõlemal juhul: murded erinevad väga lühikeste sõnade puhul (mille pikkuse logaritm on 0), samuti erineb murretes pikkuse mõju, sest Kesk-murderühmas pikkuse lisandumine peal-konstruktsiooni esinemise šansse eriti ei mõjuta, Saarte-murderühmas aga muutub pikemate sõnadega peal-konstruktsiooni esinemine järjest tõenäolisemaks.

plot_model(m5.glm, type = "int")

# või
plot_model(m5.glm, type = "pred", terms = c("LM_pikkus_log", "murderühm"))

Ülesanne

Proovime nüüd konfirmatoorse lähenemisega, millised tegurid osutuksid oluliseks siis, kui kasutame andmete tasakaalustamiseks harvemini esineva variandi (ade) üleesindamist.

ade <- ade_peal[ade_peal$cx == "ade",] # ainult alalütleva käände andmed
peal <- 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(ade), # võta ridadelt 1 kuni 722 juhuslik valim,
                 size = nrow(peal), # mille suurus oleks sama, mis peal-andmestikul
                 replace = T) # korda vaatlusi

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

Visualiseerime.

seletavad_kat <- c("TR_liik", "tegusõna_klass", "LM_mobiilsus", "LM_komplekssus", "sõnajärg", "murderühm")

dat.uus %>%
  select(cx, seletavad_kat) %>% 
  pivot_longer(., 2:length(.), 
               names_to = "tunnused", 
               values_to = "väärtused") %>%
  ggplot(aes(x = väärtused, 
               fill = cx)) + 
  geom_bar(position = "fill") +
  stat_count(geom = "text", 
             aes(label = stat(count)), 
             position = position_fill(vjust = 0.5), 
             colour = "grey40") + 
  facet_wrap("tunnused", scales = "free") +
  coord_flip()

ggplot(dat.uus, aes(x = LM_pikkus_log, y = cx)) +
  geom_boxplot()

Tee mudel, kuhu on kaasatud samad tunnused, mis mudelisse mmax.glm. Millised erinevused kahe mudeli vahel on?

Kordamisküsimused

1. Milline väide on õige?

Logistilise regressiooni vabaliige näitab šanssi, kui mudelis on 1 seletav tunnus, ja šansisuhet, kui mudelis on mitu seletavat tunnust.
Mitme seletava tunnusega mudel on alati parem kui ühe seletava tunnusega mudel.
Logistilise regressioonimudeli headust ei saa hinnata R-ruudu järgi samamoodi nagu lineaarse regressioonimudeli puhul.

2. Milliseid mudeli näitajaid näeb logistilise regressioonimudeli väljundist käsuga `summary()`?

AIC-d.
Hälbimust.
Mudeli p-väärtust.
Koefitsientide p-väärtusi.

3. Teeme kursuse küsimustiku andmete pealt mudeli, mis ennustab vastajate lemmikjooki sünniaasta, ülikoolis õpitud aastate arvu, kvantmeetodite kogemuse, lemmiklooma ja küsimustikule vastamise aasta järgi. Millised peamõjud tuleks mudelisse jätta?

load("data/kysimustik_2024.RData")
kysimustik$aasta <- as.numeric(as.character(kysimustik$aasta))
m <- glm(lemmikjook ~ synniaasta + kaua_opid + kogemused_kvant + lemmikloom + aasta, data = kysimustik, family = "binomial")
car::Anova(m)

## Analysis of Deviance Table (Type II tests)
## 
## Response: lemmikjook
##                 LR Chisq Df Pr(>Chisq)  
## synniaasta        2.5030  1    0.11363  
## kaua_opid         0.0036  1    0.95215  
## kogemused_kvant   2.9600  1    0.08535 .
## lemmikloom        1.4113  3    0.70289  
## aasta             4.9644  1    0.02587 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mitte ühtegi
kogemused_kvant ja vastamise aasta
vastamise aasta
kogemused_kvant, vastamise aasta ja lemmikloom
kõik

Järgmisel korral

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

Funktsioonide sõnastik

car::Anova() - leia seletavate tunnuste statistiline olulisus mudelis (type = “II”, kui interaktsioone pole, ja type == “III”, kui mudelis on ka interaktsioonid)
GGally::ggbivariate() - anna ülevaade uuritava tunnuse jaotumisest kõikides seletavate tunnuste kategooriates
lapply() - rakenda mingit funktsiooni igale listi liikmele (nt igale tulbale andmetabelis)
pchisq() - arvuta terve logistilise regresiooni mudeli statistiline olulisus
rbind() - liida mitu andmestikku, mil on ühesugused tulbanimed, ridupidi kokku üheks andmestikuks
sample() - võta vektorist kindla suurusega juhuslik valim
sjPlot::plot_model() - näita logistilise regressiooni mudeli koefitsiente või seletavate tunnuste mõjusid joonistel
sjPlot::tab_model() - näita logistilise regressiooni mudeli koefitsiente ilusas HTML-vormingus tabelis
tidyverse::sample_n() - võta andmestikust kindla suurusega juhuslik valim

15. praktikum

15.04.2024

Logistilise regressiooni mudelid

Tänases praktikumis

Tasakaalustamata ja tasakaalus valim

Ülevaade andmetest

Ühe tunnusega mudel tasakaalus valimi põhjal

Mitme seletava tunnusega mudel

Konfirmatoorne lähenemine

Eksploratiivne lähenemine

Mudeli mõjude visualiseerimine

Interaktsioonidega mudel

Ülesanne

Kordamisküsimused

1. Milline väide on õige?

2. Milliseid mudeli näitajaid näeb logistilise regressioonimudeli väljundist käsuga `summary()`?

3. Teeme kursuse küsimustiku andmete pealt mudeli, mis ennustab vastajate lemmikjooki sünniaasta, ülikoolis õpitud aastate arvu, kvantmeetodite kogemuse, lemmiklooma ja küsimustikule vastamise aasta järgi. Millised peamõjud tuleks mudelisse jätta?

Järgmisel korral

Funktsioonide sõnastik

15. praktikum

15.04.2024

Logistilise regressiooni mudelid

Tänases praktikumis

Tasakaalustamata ja tasakaalus valim

Ülevaade andmetest

Ühe tunnusega mudel tasakaalus valimi põhjal

Mitme seletava tunnusega mudel

Konfirmatoorne lähenemine

Eksploratiivne lähenemine

Mudeli mõjude visualiseerimine

Interaktsioonidega mudel

Ülesanne

Kordamisküsimused

1. Milline väide on õige?

2. Milliseid mudeli näitajaid näeb logistilise regressioonimudeli väljundist käsuga summary()?

3. Teeme kursuse küsimustiku andmete pealt mudeli, mis ennustab vastajate lemmikjooki sünniaasta, ülikoolis õpitud aastate arvu, kvantmeetodite kogemuse, lemmiklooma ja küsimustikule vastamise aasta järgi. Millised peamõjud tuleks mudelisse jätta?

Järgmisel korral

Funktsioonide sõnastik

2. Milliseid mudeli näitajaid näeb logistilise regressioonimudeli väljundist käsuga `summary()`?