Regressioonimudelite puhul on võimalik hinnata ühe seletava tunnuse taseme/ühiku muutuse mõju uuritava tunnuse varieerumisele. Selles mõttes on regressioon lähim meetod kontrollitud katse tegemisele, kui tegelikult on kasutada ainult vaatlusandmed. Regressioonimudelid on parameetrilised mudelid, mis tähendab, et uute vaatluste väärtuste ennustamiseks piisab kindlate parameetrite hinnangute (nt regressioonikordaja ja selle standardvea) teadmisest. Näiteks lineaarses regressioonis, kus on 1 seletav tunnus, piisab vabaliikmest ja seletava tunnuse koefitsiendist selleks, et mingi kindlusega ennustada kõiki võimalikke väärtusi. Seda aga juhul, kui andmed vastavad kindlale jaotusele ning uuritavat tunnust saab näiteks ka päriselt ennustada lineaarse funktsiooni kaudu. Kui tegelikult ei ole uuritav tunnus kirjeldatav seletavate tunnuste lineaarse (st sirgjoonelise) kombinatsiooni kaudu, ei ennusta meie mudel tegelikult uuritava tunnuse väärtuseid populatsioonis kuigi hästi.
Puudused:
Otsustuspuud on mitteparameetriliste algoritmide klass. Need ei sobita andmetele kindlat funktsiooni, vaid jõuavad sobiva funktsioonini andmete enda kaudu. Mitteparameetrilised mudelid on seega paindlikumad ning suudavad paremini arvestada korraga eri tunnuste mõjuga erinevates tingimustes. Seepärast sobivad need paremini kirjeldama komplekssemaid uurimisküsimusi, mille puhul praktikas ei pruugi nt lihtne lineaarne ennustusmudel hästi töötada. Mitteparameetrilised mudelid õpivad andmetest, seega mida rohkem on vaatlusandmeid, seda tugevam on mitteparameetrilise mudeli ennustusvõime.
Otsustuspuud
Otsustuspuu tööpõhimõtet ja ka nimetust motiveerib kujutlus puu kasvatamisest altpoolt ülespoole.
https://medium.com/analytics-vidhya/mathematics-behind-decision-tree-73ee2ef82164/
Otsustuspuu (decision tree) algoritm põhineb vaatluste binaarsel rekursiivsel ositamisel (binary recursive partitioning). See tähendab, et
Rekursiooni olemust näitlikustatakse sageli Sierpinski kolmnurga kaudu, kus üht kolmnurka saab jagada lõputult väiksemateks kolmnurkadeks.
Sierpinski kolmnurk
Otsustuspuid on eriti viimase 10 aasta jooksul kasutatud palju ka keeleteaduses. Näiteks uurisid Ruutma jt (2016) otsustuspuude abil, millest võiks sõltuda eesti murretes see, kas kaassõnad läbi, mööda, vastu, üle ja ümber paiknevad nimisõna ees (nt ümber järve) või järel (kahe kilo ümber).
Nad leidsid näiteks, et kui suhe, mida kaassõnafraas väljendab,
märgib mingit mõõdet ehk kvantiteeti või ligikaudsust (ta
oli siis juba seitsmekümne aasta ümber),
kasutatakse kaassõna ümber eeskätt nimisõna järel, tagasõna ehk
postpositsioonina.
Kui aga kaassõnafraas väljendab konkreetset ruumisuhet, sõltub kaassõna
asend sellest, kas nimisõna, millega see koos esineb, viitab elusale või
elutule referendile. Elusatega kasutatakse samuti pigem tagasõna (nt
loomad jooksid mu ümber), elututega pigem
eessõna (nt need pandi ümber vankri).
Otsustuspuude algoritmides on üht-teist erinevat, nt selles osas,
mida kasutatakse vaatluste jagamise kriteerumina.
Tuntuim algoritm on ilmselt CART (Classification
and Regression Tree). CART-mudelid jagunevad
1-(0.5^2 + 0.5^2)=0.5
. Kolme klassi puhul oleks siis Gini
indeks 1-((1/3)^2+(1/3)^2+(1/3)^2)=0.67
, nelja klassi puhul
1-((1/4)^2+(1/4)^2+(1/4)^2+(1/4)^2)=0.75
jne. Kui sõlme
jääb näiteks kahest uuritava tunnuse klassist ainult üks, siis on Gini
indeks 1-(1^2+0^2)=0
.Näites seletatakse USA 1990. aasta andmetel seda, mitu miili ühe galloni kütusega sõita saab (arvuline tunnus) vastavalt auto hinnale, tüübile, valmistamisriigile, usaldusväärsusele. Keskmine näitaja on 24.58 miili galloni kohta (kasutatud 60 auto andmeid). Kõrgeim väärtus (32.08) on ennustatud nt odavamatele kui 9446 dollarit maksvatele autodele (neid on 12), kallimatel autodel (neid on 48) hakkab omakorda mängima see, mis tüüpi auto on: väikeste autodega saab rohkem sõita kui suurtega.
Selles näites omakorda seletatakse küfoosi (lülisamba kõverdumine) olemasolu lapspatsientidel pärast korrigeerivat operatsiooni sõltuvalt patsiendi vanusest operatsiooni tegemise ajal (kuudes), selgroolülide arvust, mida operatsioon puudutas, ning sellest, mitmendast selgroolülist alates operatsiooni tehti. Patsientidest 17-l oli küfoos ka pärast operatsiooni ja 64-l ei olnud (seepärast ennustatakse kogu valimile pigem klassi absent). Kui kõige ülemise lüli järjekorranumber, mida operatsioon puudutas, oli 8 või väiksem (mudel käsitleb seda pideva tunnusena ja seepärast lubab ka väärtust 8.5), jäi küfoos pigem alles (8-l lapsel puudu, 11-l olemas). Kui alustati sellest altpoolt (eriti kui altpoolt kui 14. lüli), siis saadi haigusest pigem lahti, ehkki rolli mängis ka see, mis vanuses operatsioon tehti.
https://www.statmethods.net/advstats/cart.html
Vaata andmestikke install.packages("rpart")
,
rpart::cu.summary
ja
rpart::kyphosis
.
Probleeme klassikaliste otsustuspuu algoritmidega:
R-is on puumudelite tegemiseks erinevate algoritmidega paketid
rpart
, tree
, maptree
,
party
, partykit
jt.
Meie vaatame lähemalt natuke uuemat tingimusliku otsustuspuu algoritmi (conditional inference tree), mis kasutab sõlmedes ositamise alusena p-väärtusi (mitte nt mõnd informatsioonikriteeriumit, nagu Gini-indeks). Seletavate tunnuste valik ei ole kallutatud ja puud ei ole ülesobitamise hirmus vaja trimmida. Edaspidi viitame otsustuspuu nimetusega just tingimuslikule otsustuspuule.
Kasutame siin paketti partykit
ja selle funktsiooni
ctree()
(loe dokumentatsiooni siit
või siit).
partykit
on paketi party
edasiarendus koos
lisavõimalustega visualiseerimiseks.
Võtame sama andmestiku ade_peal.csv
, millel tegime
logistilise regressiooni mudelit ning mudeldame nüüd andmestikku
otsustuspuude kaudu.
ade_peal <- read.delim("data/ade_peal.csv", encoding = "UTF-8", stringsAsFactors = TRUE)
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(peal), # võta ridadelt 1 kuni 1310 juhuslik valim,
size = nrow(ade), # mille suurus oleks sama, mis ade-andmestikul
replace = F) # ära korda ühtegi vaatlust
peal_juhuvalim <- peal[juh_id,] # võta peal-andmestikust ainult juhuslikult valitud read
dat <- droplevels(rbind(ade, peal_juhuvalim)) # liida kaks andmestikku ridupidi uuesti kokku
#library(partykit)
# Teeme puumudeli
puu1 <- ctree(cx ~ TR_liik + tegusõna_klass + LM_mobiilsus + LM_komplekssus + LM_pikkus_log + sõnajärg + murderühm,
data = dat)
# Visualiseerime mudelit
plot(puu1)
Vaikimisi sätetega tehakse üsna keerukas puu. Saame paremini sõlmi eraldi vaadata.
Võime aga lihtsama tõlgendamise eesmärgil siin puud ka veidi
trimmida. Sellega vähendame mudeli täpsust, aga samas vähendame ka
ülesobitamise võimalust. Lisame funktsiooni argumendi
minbucket = 200
, millega ütleme, et igas lehes peab olema
vähemalt 200 vaatlust (vaikimisi on minbucket
7). Ütleme ka
argumendiga alpha = 0.01
, et 1. tüüpi vea tegemise võimalus
sõlmedes (tunnuse alusel jagatakse andmeid edasi, ehkki tegelikult on
seletava ja uuritava tunnuse vaheline seos juhuslik) võiks olla ainult
1%, ehk olulisusnivoo α = 0,05 asemel kasutame nivood α = 0,01.
puu2 <- ctree(cx ~ TR_liik + tegusõna_klass + LM_mobiilsus + LM_komplekssus + LM_pikkus_log + sõnajärg + murderühm,
data = dat,
control = ctree_control(minbucket = 200, alpha = 0.01))
plot(puu2)
LM_mobiilsus
, mis jagab vaatlused kahte
klassi selle põhjal, kas kohasõna on liigutatav või staatiline. Enamasti
on esimese jagunemise jaoks valitud tunnus üks uuritava tunnusega kõige
tugevamalt seotud tunnustest, samas mitte tingimata ise see ainuke
tugevasti seotud tunnus.TR_liik
(asesõna ja verbifraas eristuvad
peal-konstruktsiooni kasutuses liigutatavate kohasõnadega
nimisõnafraasist) ning grupis “staatiline” tunnus
murderühm
(eristuvad vaatlused, mis
tulevad Kesk- või Lõuna-murderühma kõnelejatelt, ja vaatlused, mis
tulevad Ranniku- või Saarte-murderühma kõnelejatelt).Node 3
),
ennustatakse peal-konstruktsiooni esinemise tõenäosuseks ~0.8
ja mitte-esinemise tõenäosuseks ~0.2. Siia vaatluste gruppi kuulub 250
vaatlust (n = 250
). Kuna peal esinemise tõenäosus
on suurem kui selle mitte-esinemise tõenäosus, on nendele vaatlustele
ennustatud cx
klass “peal”.Kokku on puul 7 sõlme (lehed kaasa arvatud). Nende info on mudelis listide kujul.
## Classes 'constparty', 'party' hidden list of 8
## $ node :Class 'partynode' hidden list of 5
## $ data :'data.frame': 1444 obs. of 8 variables:
## ..- attr(*, "terms")=Classes 'terms', 'formula' language cx ~ TR_liik + tegusõna_klass + LM_mobiilsus + LM_komplekssus + LM_pikkus_log + sõnajärg + murderühm
## .. .. ..- attr(*, "variables")= language list(cx, TR_liik, tegusõna_klass, LM_mobiilsus, LM_komplekssus, LM_pikkus_log, sõnajärg, murderühm)
## .. .. ..- attr(*, "factors")= int [1:8, 1:7] 0 1 0 0 0 0 0 0 0 0 ...
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. ..- attr(*, "term.labels")= chr [1:7] "TR_liik" "tegusõna_klass" "LM_mobiilsus" "LM_komplekssus" ...
## .. .. ..- attr(*, "order")= int [1:7] 1 1 1 1 1 1 1
## .. .. ..- attr(*, "intercept")= int 1
## .. .. ..- attr(*, "response")= int 1
## .. .. ..- attr(*, ".Environment")=<environment: 0x00000283ab2f2818>
## .. .. ..- attr(*, "predvars")= language list(cx, TR_liik, tegusõna_klass, LM_mobiilsus, LM_komplekssus, LM_pikkus_log, sõnajärg, murderühm)
## .. .. ..- attr(*, "dataClasses")= Named chr [1:8] "factor" "factor" "factor" "factor" ...
## .. .. .. ..- attr(*, "names")= chr [1:8] "cx" "TR_liik" "tegusõna_klass" "LM_mobiilsus" ...
## $ fitted:'data.frame': 1444 obs. of 3 variables:
## $ terms :Classes 'terms', 'formula' language cx ~ TR_liik + tegusõna_klass + LM_mobiilsus + LM_komplekssus + LM_pikkus_log + sõnajärg + murderühm
## .. ..- attr(*, "variables")= language list(cx, TR_liik, tegusõna_klass, LM_mobiilsus, LM_komplekssus, LM_pikkus_log, sõnajärg, murderühm)
## .. ..- attr(*, "factors")= int [1:8, 1:7] 0 1 0 0 0 0 0 0 0 0 ...
## .. .. ..- attr(*, "dimnames")=List of 2
## .. ..- attr(*, "term.labels")= chr [1:7] "TR_liik" "tegusõna_klass" "LM_mobiilsus" "LM_komplekssus" ...
## .. ..- attr(*, "order")= int [1:7] 1 1 1 1 1 1 1
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: 0x00000283ab2f2818>
## $ names : NULL
## $ info :List of 2
## $ update:function (subset, weights, control, doFit = TRUE)
## $ trafo :function (subset, weights, info, estfun, object, ...)
# 1. sõlme info: LM_mobiilsus = "liigutatav" ja LM_mobiilsus = "staatiline"
# See näitab terve puu jagunemisi
# Sama, mis lihtsalt puu2
puu2[[1]]
##
## Model formula:
## cx ~ TR_liik + tegusõna_klass + LM_mobiilsus + LM_komplekssus +
## LM_pikkus_log + sõnajärg + murderühm
##
## Fitted party:
## [1] root
## | [2] LM_mobiilsus in liigutatav
## | | [3] TR_liik in asesõna, verbifraas: peal (n = 250, err = 21.6%)
## | | [4] TR_liik in nimisõnafraas: peal (n = 219, err = 37.9%)
## | [5] LM_mobiilsus in staatiline
## | | [6] murderühm in Kesk, Lõuna: ade (n = 616, err = 49.4%)
## | | [7] murderühm in Ranniku, Saarte: ade (n = 359, err = 24.0%)
##
## Number of inner nodes: 3
## Number of terminal nodes: 4
##
## Model formula:
## cx ~ TR_liik + tegusõna_klass + LM_mobiilsus + LM_komplekssus +
## LM_pikkus_log + sõnajärg + murderühm
##
## Fitted party:
## [7] root: ade (n = 359, err = 24.0%)
##
## Number of inner nodes: 0
## Number of terminal nodes: 1
##
## Model formula:
## cx ~ TR_liik + tegusõna_klass + LM_mobiilsus + LM_komplekssus +
## LM_pikkus_log + sõnajärg + murderühm
##
## Fitted party:
## [6] root: ade (n = 616, err = 49.4%)
##
## Number of inner nodes: 0
## Number of terminal nodes: 1
Iga sõlme/jagunemise kohta näeme niisiis
n
)err
) ehk nende
vaatluste osakaalu, mis jäävad vähemusse. Ehk kui nt 6. sõlme (lõpulehe)
ennustatud klass on “ade” ning err = 49.4%, on järelikult sellesse lehte
jagatud tegelike vaatluste hulgas “peal” osakaal 0.494 ja “ade” osakaal
seega 0.506. See leht ei ole niisiis kuigi stabiilne.Kui vaatame nüüd uuesti üldpilti, siis näeme, et liigutatavate
kohtade hulgas (LM_mobiilsus = "liigutatav"
) on
peal-konstruktsiooni esinemise tõenäosus keskmiselt suurem kui
staatiliste kohtade hulgas (LM_mobiilsus = "staatiline"
),
ent peal-konstruktsiooni esinemist mõjutavad kummaski rühmas
veel muud tegurid, mis esinemise tõenäosust tõstavad või langetavad.
Vaatame lähemalt näiteks 1) konteksti, kus
peal-konstruktsiooni esinemise tõenäosus on kõige suurem
(Node 3
), ning 2) konteksti, kus see on kõige väiksem
(Node 7
).
1) peal esinemise tõenäosus on kõige suurem (79.4%) kontekstis, kus
##
## Model formula:
## cx ~ TR_liik + tegusõna_klass + LM_mobiilsus + LM_komplekssus +
## LM_pikkus_log + sõnajärg + murderühm
##
## Fitted party:
## [3] root: peal (n = 250, err = 21.6%)
##
## Number of inner nodes: 0
## Number of terminal nodes: 1
Vaatame, kas andmestikus nendele tunnustele vastavate vaatluste arv klapib.
library(tidyverse)
dat %>%
filter(LM_mobiilsus == "liigutatav",
TR_liik %in% c("asesõna", "verbifraas")) %>%
nrow()
## [1] 250
2) peal esinemise tõenäosus on kõige väiksem (24%) kontekstis, kus
##
## Model formula:
## cx ~ TR_liik + tegusõna_klass + LM_mobiilsus + LM_komplekssus +
## LM_pikkus_log + sõnajärg + murderühm
##
## Fitted party:
## [7] root: ade (n = 359, err = 24.0%)
##
## Number of inner nodes: 0
## Number of terminal nodes: 1
Vaatame, kas andmestikus nendele tunnustele vastavate vaatluste arv klapib.
## [1] 359
Igas sõlmes võistleb jagunemise tekitamiseks omavahel mitu seletavat
tunnust, mille hulgast valitakse p-väärtuse alusel ainult 1.
Seda, mis tunnused veel sõlmedes kandideerisid, saab vaadata käsuga
sctest.constparty()
.
# Kõikide sõlmede jagunemised
# sctest.constparty(puu2)
# 1. sõlme jagunemine
sctest.constparty(puu2)[1]
## $`1`
## TR_liik tegusõna_klass LM_mobiilsus LM_komplekssus LM_pikkus_log
## statistic 5.5626041 6.487519e+01 1.199936e+02 4.633788e+01 0.0766731
## p.value 0.3609178 1.913658e-12 4.442550e-27 6.966501e-11 0.9999765
## sõnajärg murderühm
## statistic 0.01323792 7.91153e+01
## p.value 0.99999995 3.32574e-16
Selline info võib olla informatiivne siis, kui sõlme valitud seletav
tunnus oli ainult õige natuke parem kui mõni teine tunnus, mida
jagunemiseks ei valitud. Esimese sõlme jagunemises näeme, et kõik
tunnused peale trajektori liigi (!), kohafraasi pikkuse ja sõnajärje on
uuritava tunnusega statistiliselt oluliselt seotud, ent kohasõna
mobiilsuse p-väärtus oli mitu korda väiksem kui teiste tunnuste
p-väärtused, seega valiti see esimese jagunemise aluseks. Mida
allapoole puus aga liigume, seda väiksemaks jääb jagunemiseks sobivate
tunnuste valik.
Seda, mis statistikut tunnustevahelise seose hindamiseks kasutatakse,
sõltub sellest, mis tüüpi (arvuline/kategoriaalne) on seletav tunnus ja
mis tüüpi on uuritav tunnus. Loe soovi korral rohkem nt siit.
Miks siis see TR_liik
ikkagi siin puus 2. sõlmes
oluliseks osutus, kui ta uuritava tunnuse varieerumisega 1. sõlmes
individuaalselt üldse seotud ei ole ning osutus ebaoluliseks ka
regressioonis?
## $`2`
## TR_liik tegusõna_klass LM_komplekssus LM_pikkus_log sõnajärg
## statistic 1.895034e+01 21.43852387 9.4776945 1.8361491 1.0913917
## p.value 4.603138e-04 0.00155401 0.0124144 0.6856231 0.8784283
## murderühm
## statistic 13.11021221
## p.value 0.02613611
Puumudelid töötavad n-ö lokaalselt, st iga järgmist olulist tunnust otsitakse eelnevate jagunemiste põhjal saadud gruppides, mitte enam tervest andmestikust. Regressioon jällegi hindab mõjusid n-ö globaalselt, fikseerides küll mingi kindla võrdluskonteksti, ent hinnates kõiki parameetreid sellesama konteksti suhtes. Nõnda on puumudelid peamõjude tuvastamisel pisut kehvemad kui regressioonimudelid (näiteks ei saa me siit kuskilt teada, et tegelikult on ka kohafraasi komplekssus tugevalt peal- ja käändekonstruktsiooni varieerumisega seotud). Samas võimaldavad need hõlpsamalt leida ning visualiseerida interaktsioone ning tähele tuleb siinjuures panna ka, milliste parameetritega puu oleme ehitanud. Kui lubame näiteks, et lõpusõlmes võib olla ka ainult 50 vaatlust, et minna arvesse kui uuritava tunnuse suhtes teistest erinev vaatluste klass, leiame jooniselt jälle ka kohasõna komplekssuse, mis mängib olulisemat rolli Kesk- ja Lõuna-murderühma staatiliste kohasõnadega.
Ehkki otsustuspuude kasutamine on tavalisem klassifitseerimisprobleemide puhul, saab seda küllalt edukalt kasutada ka regressiooniülesannetes, st mingi arvulise uuritava tunnuse keskväärtuse ennustamise jaoks.
Uurime näiteks puumudelitega uuesti kõigepealt seda, kas sõna äratundmiseks kuluv aeg sõltub sõna pikkusest ja sagedusest.
Regressioonipuu visuaalne väljund (vaikimisi karpdiagramm) annab iga lõpliku sõlme kohta uuritava tunnuse ennustatud kvartiiljaotuse selles sõlmes.
Näeme, et reaktsiooniaega mõjutab eeskätt sõna pikkus: kui sõna on pikem kui 10 tähemärki, siis on reaktsiooniaeg pikem kui siis, kui sõna pikkus on kuni 10 tähemärki. Lühemate sõnade puhul jällegi tekib erinevus lühikeste (<= 7 tm) ja keskmise pikkusega (>7 tm) sõnade vahel. Nendes gruppides omakorda hakkab rolli mängima sõna sagedus: mõlemal juhul on sagedamate sõnade äratundmiseks kuluv aeg lühem.
Mudelist saame kätte ka iga vaatluste grupi ennustatud keskväärtuse.
##
## Model formula:
## Mean_RT ~ Length + Freq
##
## Fitted party:
## [9] root: 989.824 (n = 18, err = 665437.8)
##
## Number of inner nodes: 0
## Number of terminal nodes: 1
Arvuliste uuritavate tunnustega näitab iga lehe n endiselt sellesse rühma kuuluvate vaatluste arvu, ent err sedapuhku jääkide ruutude summat, mis mäletatavasti oli seotud üksikute vaatluste erinevusega grupi keskmisest.
## [1] 1125.42 948.33 816.94 931.13 1054.77 931.04 620.79 978.81 942.86
## [10] 1314.33 768.41 915.96 1089.00 904.47 1216.81 923.74 875.27 1458.75
## [1] 989.8239
## [1] 665437.8
Mida suurem on jääkide ruutude summa, seda erinevamad vaatlused sellesse lehte kuuluvad. Näiteks joonise järgi peaks olema sõlmes 8 jääkide ruutude summa väiksem, kuna karpdiagrammi kuju viitab sellele, et vaatlused on koondunud tihedamalt keskväärtuse ümber.
##
## Model formula:
## Mean_RT ~ Length + Freq
##
## Fitted party:
## [8] root: 746.839 (n = 16, err = 61399.5)
##
## Number of inner nodes: 0
## Number of terminal nodes: 1
Konkreetne arv mõistagi sõltub sellest, mis ühikutes on uuritav
tunnus. Näiteks kui teisendaksime reaktsiooniaja millisekunditest
minutitesse, oleks sõlmede err
väärtused peaaegu 0.
##
## Model formula:
## Mean_RT_min ~ Length + Freq
##
## Fitted party:
## [1] root
## | [2] Length <= 10
## | | [3] Length <= 7
## | | | [4] Freq <= 702: 0.013 (n = 14, err = 0.0)
## | | | [5] Freq > 702: 0.011 (n = 25, err = 0.0)
## | | [6] Length > 7
## | | | [7] Freq <= 393: 0.014 (n = 27, err = 0.0)
## | | | [8] Freq > 393: 0.012 (n = 16, err = 0.0)
## | [9] Length > 10: 0.016 (n = 18, err = 0.0)
##
## Number of inner nodes: 4
## Number of terminal nodes: 5
Vaatame nüüd uuesti ka 2019. a valimiste andmeid.
load("data/kandidaadid2019.RData")
puu4 <- ctree(hääli_kokku ~ sugu + haridus + tähtkuju,
data = kandidaadid2019,
subset = haridus != "Algharidus")
plot(puu4)
Kuna häältesaak on tugevalt madalamate väärtuste poole kaldu, siis on parem ka puumudelis kasutada logaritmitud uuritavat tunnust, et vähendada erandlikult palju hääli saanud kandidaatide mõju.
puu5 <- ctree(log(hääli_kokku) ~ sugu + haridus + tähtkuju,
data = kandidaadid2019,
subset = haridus != "Algharidus")
plot(puu5)
Näeme, et kõrgharidusega kandidaatide seas statistiliselt olulist erinevust mees- ja naiskandidaatide häältesaagis ei olnud, aga kesk- ja põhiharidusega kandidaatide seas said naised keskmiselt pisut vähem hääli. Näeme siiski ka, et see erinevus ei ole väga suur. Kandidaadi tähtkuju häältesaagi puhul mingit rolli ei mängi.
##
## Model formula:
## log(hääli_kokku) ~ sugu + haridus + tähtkuju
##
## Fitted party:
## [1] root
## | [2] haridus in Keskharidus (sh keskeriharidus), Põhiharidus
## | | [3] sugu in mees: 4.507 (n = 203, err = 487.2)
## | | [4] sugu in naine: 3.777 (n = 76, err = 152.4)
## | [5] haridus in Kõrgharidus: 5.246 (n = 819, err = 1845.1)
##
## Number of inner nodes: 2
## Number of terminal nodes: 3
Teeme mõned seletavad tunnused juurde ja proovime neid ka mudelisse panna.
# Kas kandidaat kasutab meili?
kandidaadid2019$kasutab_meili <- as.factor(ifelse(grepl("@", kandidaadid2019$kontakt), "jah", "ei"))
# Kui pikk on kandidaadi nimi
kandidaadid2019$nime_pikkus <- nchar(kandidaadid2019$nimi)
# Kui vana kandidaat on?
kandidaadid2019$vanus <- 2019 -as.numeric(substr(kandidaadid2019$sünniaeg, 1, 4))
puu6 <- ctree(log(hääli_kokku) ~ sugu + haridus + tähtkuju + kasutab_meili +
nime_pikkus + vanus,
data = kandidaadid2019,
subset = haridus != "Algharidus")
plot(puu6)
Ühelgi lisatud tunnusel ei ole kandidaadi saadud häälesaagiga puumudelis olulist seost.
Ülesanne: kontrolli, kas mõni vaadeldud tunnustest mõjutas logaritmitud häältesaaki 2023. aasta valimistel?
# loe kategoriaalsed tunnused sisse faktoritena
valimised <- read.delim("data/kandidaadid_2023.csv", stringsAsFactors = T)
valimised$kasutab_meili <- as.factor(ifelse(grepl("@", valimised$kontakt), "jah", "ei"))
valimised$nime_pikkus <- nchar(as.character(valimised$nimi))
names(valimised)
## [1] "ringkond" "nimekiri"
## [3] "nimi" "sünniaeg"
## [5] "vanus" "sugu"
## [7] "tähtkuju" "haridus"
## [9] "erakond" "amet"
## [11] "kontakt" "hääli_kokku"
## [13] "e.hääli" "hääli_välisriigist"
## [15] "hääli_väljaspool_ringkonda" "kasutab_meili"
## [17] "nime_pikkus"
Vaatame veel, mis saab siis, kui ühelgi seletaval tunnusel ei ole uuritava tunnuse seletamisel olulist rolli.
load("data/kysimustik_2024.RData")
# Kategoriaalne uuritav tunnus
puu7 <- ctree(lemmikloom ~ ., data = kysimustik)
plot(puu7)
# Arvuline uuritav tunnus
puu8 <- ctree(kaua_opid ~ kogemused_kvant + lemmikjook + kursuse_labimine,
data = kysimustik)
plot(puu8)
Sel juhul saame lihtsalt uuritava tunnuse jaotuse andmetes.
Pakett partykit
pakub rea võimalusi otsustuspuude
visualiseerimiseks. Vaatame neist paari. Põhilised elemendid, mida
graafikul saame modifitseerida, on
inner_panel
: puu sõlmed, kust jagatakse vaatlused mingi
tunnuse alusel kaheks;terminal_panel
: puu lehed ehk mudeli lõpusõlmed, kus on
klassid, tõenäosused jms.Neid elemente saavad omakorda muuta erinevad
node
-funktsioonid, nt
node_inner()
: üldine sisemiste sõlmede modifitseerimise
funktsioonnode_terminal()
: üldine lõpusõlmede modifitseerimise
funktsioonnode_barplot()
: tulpdiagrammesitus kategoriaalse
uuritava tunnuse klassidelenode_boxplot()
: karpdiagrammesitus arvulisele
uuritavale tunnusele# Teeme visualiseerimise näitlikustamiseks
# ühe veidi rohkem trimmitud puu, kus
# alpha = 0.01 ja
# minimaalne vaatluste arv klassides on 200.
puu9 <- ctree(cx ~ TR_liik + tegusõna_klass + LM_mobiilsus + LM_komplekssus + LM_pikkus_log + sõnajärg + murderühm,
data = dat,
control = ctree_control(alpha = 0.01,
minbucket = 200))
plot(puu9)
# Kuva puu lehed ilma tulpdiagrammita
# Tee tekst graafikul väiksemaks ja halliks
# Kuva sisemised sõlmed ilma numbri ja p-väärtuseta
plot(puu9,
type = "simple",
gp = gpar(fontsize = 9,
col = "darkslategrey"),
inner_panel = node_inner(obj = puu9,
id = FALSE,
pval = FALSE))
# Tee tekst graafikul väiksemaks ja halliks
# Kuva sisemised sõlmed ilma numbri ja p-väärtuseta
# Kuva lõpulehed tulpdiagrammina, kus
# tulbad on üksteise kõrval ja vastupidises järjekorras
# ning värvitud oranžiks ja halliks.
plot(puu9,
gp = gpar(fontsize = 8,
col = "darkslategrey"),
inner_panel = node_inner(obj = puu9,
id = FALSE,
pval = FALSE),
terminal_panel = node_barplot(obj = puu9,
beside = TRUE,
reverse = TRUE,
fill = c("orange2",
"grey35")))
Lisaks partykit
-paketi enda võimalustele on olemas ka
pakett ggparty
, mis on tehtud ggplot2
laiendusena partykit
paketi objektide visualiseerimiseks.
Kui rohkem huvi, tasub selle kohta uurida nt siit
või siit.
# install.packages("ggparty")
library(ggparty)
ggparty(puu9) +
geom_edge() + # "oksad"
geom_edge_label(max_length = 5) + # okste tekst lühendada
geom_node_label( # sõlmede tekst
ids = "inner", # puu sisemistele sõlmedele
line_list = list(
aes(label = splitvar), # ühel real tunnuse nimi
aes(label = paste("N = ", # teisel real kleebi kokku "N = " ja
nodesize))), # vaatluste arv
line_gpar = list( # teksti graafilised parameetrid
list(size = 13), # tunnuse nime suurus
list(size = 10))) + # vaatluste arvu suurus
geom_node_label( # sõlmede tekst
ids = "terminal", # puu lõpusõlmedele ("lehtedele")
aes(label = paste0("Sõlm ", # kleebi kokku "Sõlm ",
id, # sõlme number
", N = ", # tekst " , N = " ja
nodesize)), # vaatluste arv
nudge_y = 0.015, # liiguta teksti natuke kõrgemale
nudge_x = 0.02) + # liiguta teksti natuke paremale
geom_node_plot( # lõpusõlmede ehk lehtede joonis
gglist = list( # siia sisse tulevad kõik ggploti tavalised funktsioonid
geom_bar(aes(x = "",
fill = cx),
position = "fill",
color = "grey50"),
theme_minimal(),
scale_fill_manual(values = c("grey35",
"orange2")),
labs(x = "", y = ""),
scale_y_continuous(breaks = c(0,0.5,1))))
Nii nagu logistilises regressiooniski ennustab otsustuspuu sisuliselt igale teatud omadustega vaatlusele mingi uuritava tunnuse klassi (kat. uuritava tunnusega) või keskväärtuse (arv. uuritava tunnusega), mille põhjalt saab hinnata mudeli sobivust päris andmetega.
Nagu logistilises regressioonis, saame ka kategoriaalse uuritava tunnusega puumudeli puhul hinnata mudeli headust klassifitseerimistäpsuse abil.
# Mudeli ennustatud klassid
# Puumudeli puhul saab funktsiooniga predict()
# kategoriaalse uuritava tunnuse puhul
# vaikimisi kohe ennustatud klassid, mitte tõenäosused.
# Klassi kohal olev number näitab selle lehe/lõpusõlme
# numbrit, kuhu vaatlus on klassifitseeritud.
enn <- predict(puu2)
head(enn)
## 7 7 7 3 4 4
## ade ade ade peal peal peal
## Levels: ade peal
## [1] ade ade ade ade ade ade
## Levels: ade peal
## enn
## teg ade peal
## ade 585 137
## peal 390 332
# Klassifitseerimistäpsus
(klass <- (t["peal", "peal"]+t["ade", "ade"])/sum(t)) # klassifitseerimistäpsus
## [1] 0.6350416
Meie puumudel ennustab ~63% juhtudest vaatlusele õige uuritava tunnuse klassi. See on pisut parem kui ilma seletavate tunnusteta kõikidele vaatlustele lihtsalt sagedamat või suvaliselt üht klassi ennustades (50%).
##
## ade peal
## 0.5 0.5
Binaarse uuritava tunnuse puhul saab leida ka
C-indeksi/AUC, mis hindab mudeli eristusvõimet mitte
ennustatud klasside, vaid tõenäosuste põhjal. Kuna nüüd on meil vaja
tõenäosusi, lisame funktsiooni predict()
argumendi
type = "prob"
. See annab meile iga vaatluse jaoks
ennustatud “ade” ja “peal” klassi tõenäosused (kokku annavad 1).
## ade peal
## 7 0.7604457 0.2395543
## 7 0.7604457 0.2395543
## 7 0.7604457 0.2395543
## 3 0.2160000 0.7840000
## 4 0.3789954 0.6210046
## 4 0.3789954 0.6210046
C-indeksi arvutamiseks vajame sellest tabelist ainult ühe klassi tõenäosusi (“peal”).
# Leiame puumudeli ennustatud tõenäosused
enn_probs <- predict(puu2, type = "prob")[,"peal"]
head(enn_probs)
## 7 7 7 3 4 4
## 0.2395543 0.2395543 0.2395543 0.7840000 0.6210046 0.6210046
## Area under the curve: 0.6975
# Või
# Teisendame uuritava tunnuse tegelikud väärtused
# samuti tõenäosusteks 0 ("ei") ja 1 ("jah")
teg_probs <- ifelse(dat$cx == "peal", 1, 0)
library(Hmisc)
# C-indeks
somers2(enn_probs, teg_probs)
## C Dxy n Missing
## 0.6974682 0.3949363 1444.0000000 0.0000000
Kuna meie C väärtus on veidi väiksem kui 0.7, võime paigutada mudeli eristusvõime kehva ja rahuldava piirimaile. Ilmselt aitaks, kui oleksime oma p-väärtuste ja lõpusõlmede vaatluste arvuga natuke heldemad ega trimmiks puud nii rangelt. Teisalt võib see viia selleni, et meie puumudel hakkab ülesobitama ega tööta teistel valimitel eriti hästi.
Kui uuritaval tunnusel on rohkem kui 2 klassi, on AUC kasutamine
pisut keerukam. On pakutud, et sellisel juhul võib leida eraldi
AUC-väärtused kõikide tunnuse klasside paaride vahel ning seejärel võtta
nende aritmeetiline keskmine. Selleks on R-is eraldi funktsiooon multiclass.roc()
paketist pROC
.
Arvulise uuritava tunnusega mudeli jaoks vaatame uuesti reaktsiooniaegade andmestikku.
Üheks variandiks hinnata mudeli ligikaudset headust on leida mudeli ennustatud keskväärtuste ja tegelike vaatluste väärtuste vaheline korrelatsioon.
## 7 7 4 4 9 9
## 855.6022 855.6022 786.9471 786.9471 989.8239 989.8239
## [1] 819.19 977.63 908.22 766.30 1125.42 948.33
## [1] 0.7026489
Vaikimisi kasutab funktsioon cor()
mäletatavasti
Pearsoni korrelatsioonikordajat, mille kasutamise eelduseks on see, et
tunnused on normaaljaotusega ning tunnuste vahel on lineaarne seos.
Ehkki puumudelite hindamise puhul seda sageli ignoreeritakse ning
lähtutakse Pearsoni kordajast kui ligikaudsest näitajast, võiksime
siiski selle kasutamise eeldused üle kontrollida.
# Lineaarne seos?
ggplot(data = data.frame(enn, teg),
aes(x = teg, y = enn)) +
geom_point(alpha = 0.5) +
geom_smooth() +
ggtitle("Ennustatud ja tegelikud väärtused")
##
## Shapiro-Wilk normality test
##
## data: teg
## W = 0.92006, p-value = 1.418e-05
##
## Shapiro-Wilk normality test
##
## data: enn
## W = 0.87098, p-value = 7.677e-08
Kuna meil ei ole tegemist päris lineaarse seosega, siis oleks hea kontrollida mudeli headuse hindamisel seost ka Spearmani korrelatsioonikordajaga.
## [1] 0.7352883
Selleks, et puumudeli (ligikaudset) headust võiks võrrelda ka nt lineaarse regressioonimudeli omaga, võib Pearsoni korrelatsioonikordaja põhjal leida mudeli poolt seletatud varieerumise protsendi R2.
## [1] 0.4937155
Teine variant jõuda sellesama Pearsoni R2-ni on lahutada 1st koguruutude summa (Total Sum of Squares) ja jääkide ruutude summa (Residual Sum of Squares) jagatis.
# koguruutude summa (Total Sum of Squares):
# seletamata variatsioon andmetes
TSS <- sum((teg - mean(teg))^2)
# jääkide ruutude summa (Residual Sum of Squares):
# seletamata variatsioon pärast seletavate tunnuste kasutamist
RSS <- sum((teg - enn)^2)
# R-ruut
(R2 <- 1-RSS/TSS)
## [1] 0.4937155
Tingimuslikud otsustuspuud on väga head tunnustevaheliste interaktsioonide visuaalseks vaatlemiseks, ent need ei ole kõige stabiilsemad mudelid ning võivad olla tundlikud sisendandmete suhtes. Kõikumised andmestikus võivad anda natuke teistsuguse lõpptulemuse, eriti väiksemate andmestike puhul.
Selleks, et ühe puu põhjalt mitte liiga palju järeldada, tuleks
mudelit valideerida. Valideerimiseks võib andmestiku jagada jällegi
kaheks: treening- ja testandmeteks. Treeningandmetel tehtud mudeli
ennustusi saab testida testandmetel funktsiooniga
predict(puumudel, newdata = testandmestik)
ning arvutada
nende põhjalt mudeli täpsuse näitajad.
Teeme seda siin ilma puud trimmimata (st jätame määramata argumendid
minbucket
ja alpha
).
# Võtame treeningandmestikku 70% originaalandmestiku vaatlustest
set.seed(100)
trenn_id <- sample(1:nrow(dat), round(0.7*nrow(dat)), replace = FALSE)
trenn <- dat[trenn_id,]
test <- dat[-trenn_id,]
puu_trenn <- ctree(cx ~ TR_liik + tegusõna_klass + LM_mobiilsus + LM_komplekssus + LM_pikkus_log + sõnajärg + murderühm,
data = trenn)
Mudeli täpsus treeningandmestikul:
enn_trenn <- predict(puu_trenn)
teg_trenn <- trenn$cx
t_trenn <- table(teg_trenn, enn_trenn)
# klassifitseerimistäpsus
(klass_trenn <- (t_trenn["peal","peal"]+t_trenn["ade", "ade"])/sum(t_trenn))
## [1] 0.6567755
# library(Hmisc)
enn_probs_trenn <- predict(puu_trenn, type = "prob")[,"peal"]
teg_probs_trenn <- ifelse(trenn$cx == "peal", 1, 0)
somers2(enn_probs_trenn, teg_probs_trenn)["C"]
## C
## 0.7153017
Mudeli täpsus testandmestikul:
enn_test <- predict(puu_trenn, newdata = test)
teg_test <- test$cx
t_test <- table(teg_test, enn_test)
# klassifitseerimistäpsus
(klass_test <- (t_test["peal","peal"]+t_test["ade","ade"])/sum(t_test))
## [1] 0.6605081
# library(Hmisc)
enn_probs <- predict(puu_trenn, newdata = test, type = "prob")[,"peal"]
teg_probs <- ifelse(test$cx == "peal", 1, 0)
somers2(enn_probs, teg_probs)["C"]
## C
## 0.7075913
Antud juhul sobib mudel kirjeldama ka testandmestikku, ent kui
headuse näitajad test- ja treeningandmestiku vahel on väga erinevad,
siis võib proovida mudelit erinevate parameetrite abil trimmida (vt
control = ctree_control()
).
Teine variant ennustustes kindlamaks saamiseks on kasutada nn ansamblimeetodeid (ensemble methods) nagu juhumetsad, millest räägime järgmises praktikumis.
ggparty::ggparty()
- joonista puumudel
ggplot2
võimalusi kasutadesHmisc::somers2()
- leia binaarse kategoriaalse uuritava
tunnusega puumudeli C-indeks/AUCpartykit::ctree()
- tee tingimusliku otsustuspuu
mudelpartykit::sctest.constparty()
- küsi jagunemiseks
proovitud seletavate tunnuste statistikuid ja p-väärtusi iga jagunemise
sõlme kohtapredict()
- leia puumudeli ennustused igale
konkreetsele andmestiku vaatluselepROC::auc()
- leia binaarse kategoriaalse uuritava
tunnusega puumudeli C-indeks/AUC