1. Stack Overflow adatbázis

Töltsük be a Moodle-n található StackOverflowHungary2020.xlsx adattáblát, ami Stack Overflow programozói közösségi oldal 2020-as felmérése a világ amatőr és profi programozóiról 60 változó szerint. A teljes adatbázis (és a korábbi+újabb évek felmérései) erről a linkről elérhető.
A mi Moodle-n található Excel fájlunkban csak a 2020-as felmérés 210 magyar kitöltőjének válaszai szerepelnek az alábbi 9 változó szerint :

  • Age: A válaszadó életkora (év)
  • Age1stCode: A válaszadó életkora első sor programkódjának megírásakor (év)
  • YearsCodePro: Programozási tapasztalat a tanulmányokat nem beleszámítva (év)
  • MonthlyHuf: Havi bruttó fizetés Forintban
  • Gender: Válaszadó neme
  • EdLevel: Legmagasabb befejezett iskolai végzettség
  • Employment: Foglalkoztatási státusz (teljes munkidő; részmunkaidő; egyéni vállalkozó)
  • JobSat: Elégedettség a jelenlegi munkahelyen
  • OpSys: Használt operációs rendszer (Windows; Linux; MacOS)

Állítsuk be a Working Directory-t, majd olvassuk be az Excel fájlt egy ProgData névre hallgató data frame-be:

setwd("~/Oktatás 2021221/Öko. I/R jegyzet")

library(readxl)
ProgData <- readxl::read_excel("StackOverflowHungary2020.xlsx")

str(ProgData)
## tibble [210 x 9] (S3: tbl_df/tbl/data.frame)
##  $ Age         : num [1:210] 31 27 22 31 24 33 43 23 38 22 ...
##  $ Age1stCode  : num [1:210] 12 13 11 14 15 24 10 20 12 9 ...
##  $ YearsCodePro: num [1:210] 8 6 5 3 2 1 20 0.5 15 3 ...
##  $ MonthlyHuf  : num [1:210] 788080 1023642 1702347 511981 337703 ...
##  $ Gender      : chr [1:210] "Man" "Man" "Man" "Man" ...
##  $ EdLevel     : chr [1:210] "Unfinished university" "Master’s degree" "Bachelor’s degree" "Master’s degree" ...
##  $ Employment  : chr [1:210] "Employed full-time" "Employed full-time" "Independent contractor, freelancer, or self-employed" "Employed full-time" ...
##  $ JobSat      : chr [1:210] "Slightly dissatisfied" "Slightly satisfied" "Very satisfied" "Slightly satisfied" ...
##  $ OpSys       : chr [1:210] "Windows" "Windows" "Windows" "Linux-based" ...

Láthatjuk, hogy meg is van mind a \(210\) megfigyelésünk és mind a \(9\) változónk. Azt vehetjük észre rögtön az str eredményéből, hogy az utolsó 5 oszlopban lévő változót érdemes factor adattípusúvá konvertálni, mivel nominális, azaz szöveges jelentéstartalmú változók ezek, relatíve kevés lehetséges értékkel.
Ha “egy kupacban”, azaz egymás mellett lévő oszlopokon szeretnénk ugyan azt a műveletet végrehajtani egy R data frame-ben, akkor sokat könnyít az életünkön a lapply függvény. Ez a függvény egy data frame kijelölt oszlopain hajtja végre ugyan azt a műveletet. Jelen esetben a factorrá konvertálást a ProgData data frame utolsó 5, azaz \(5-9.\) oszlopain:

ProgData[,5:9] <- lapply(ProgData[,5:9], as.factor)

str(ProgData)
## tibble [210 x 9] (S3: tbl_df/tbl/data.frame)
##  $ Age         : num [1:210] 31 27 22 31 24 33 43 23 38 22 ...
##  $ Age1stCode  : num [1:210] 12 13 11 14 15 24 10 20 12 9 ...
##  $ YearsCodePro: num [1:210] 8 6 5 3 2 1 20 0.5 15 3 ...
##  $ MonthlyHuf  : num [1:210] 788080 1023642 1702347 511981 337703 ...
##  $ Gender      : Factor w/ 3 levels "Man","NoAnswer",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ EdLevel     : Factor w/ 6 levels "Bachelor’s degree",..: 6 3 1 3 5 1 3 4 2 5 ...
##  $ Employment  : Factor w/ 3 levels "Employed full-time",..: 1 1 3 1 1 3 3 1 1 3 ...
##  $ JobSat      : Factor w/ 5 levels "Neither satisfied nor dissatisfied",..: 2 3 5 3 3 3 3 3 3 3 ...
##  $ OpSys       : Factor w/ 3 levels "Linux-based",..: 3 3 3 1 3 1 3 3 3 2 ...

Voilá, meg is vagyunk a konvertálásokkal, és nem kellett mind az 5 oszlopra végigpötyögni ugyan azt a kódot. :)

2. Multinomiális Logisztikus regresszió alapfeladata

Ebből a ProgData adatbázisból nagyon sok érdekes modellt lehetne építeni. Adná magát egy sima OLS regresszió a MonthlyHuf eredményváltozóra pl. De most nem ebbe az irányba indulunk el, hanem a JobSat-ra, mint eredményváltozóra akarunk egy regressziós modellt építeni!
Tehát meg akarjuk tudni, hogy mitől függ az, hogy egy programozó mennyire elégedett a munkájával.

A problémáink ott kezdődnek, hogy mint látjuk a fenti str eredményből, a JobSat változó egy nominális, azaz szövegesen megadott változó, 5 lehetséges értékkel. A levels függvénnyel meg is nézhetjük, hogy mik ezek az értékek:

levels(ProgData$JobSat)
## [1] "Neither satisfied nor dissatisfied" "Slightly dissatisfied"             
## [3] "Slightly satisfied"                 "Very dissatisfied"                 
## [5] "Very satisfied"

Szóval, az eredményváltozónk nominális, én nem csak 2 lehetséges értéke van, hanem 5. Szóval, az eddig ismert, sima bináris logisztikus regressziót elfelejthetjük. Jön helyette a multinomiális logisztikus regresszió! Éljen! :)

De mindjárt látjuk, hogy nem akkora spanyolviasz ez a multinomiális logit modell. Lássuk a matematikai modellt!

Szóval, az eredményváltozó, azaz \(Y\) nomiális mérési skálájú (szöveges adat), több mint 2 lehetséges értékkel.

Esetünkben \(Y=JobSat\), melynek értékkészlete:

  • \(1 =\) Very dissatisfied
  • \(2 =\) Slightly dissatisfied
  • \(3 =\) Neither satisfied nor dissatisfied
  • \(4 =\) Slightly satisfied
  • \(5 =\) Very satisfied

Ekkor azt mondjuk, hogy \(K=5\), azaz \(K\) az \(Y\) értékkészletét jelöli.

Magyarázóváltozóknak első körben tekintsünk csak 3 mennyiségi változót:

  • \(X_1 =\) Age
  • \(X_1 =\) YearsCodePro
  • \(X_1 =\) MonthlyHuf

Ekkor a magyarázóváltozók száma \(k=3\).

A jelölések bevezetése után jön a TRÜKK! Az \(Y\) eredményváltozóra vezessünk be egy dummy kódolást referencia kategóriával!
Azaz:

És minden dummy változóra építsünk külön-külön bináris logisztikus regressziót. Ezzel \(5-1=4\) db regressziónk lesz, ahol a bináris logit \(0\) szintje az eredményváltozó referencia kategóriája:

  • \(Logit(Y_1)=\ln(\frac{P(Y=1)}{P(Y=5)})=\beta_{11} \times X_1 +\beta_{12} \times X_2 + \beta_{13} \times X_3 + \beta_{10}\)
  • \(Logit(Y_2)=\ln(\frac{P(Y=2)}{P(Y=5)})=\beta_{21} \times X_1 +\beta_{22} \times X_2 + \beta_{23} \times X_3 + \beta_{20}\)
  • \(Logit(Y_3)=\ln(\frac{P(Y=3)}{P(Y=5)})=\beta_{31} \times X_1 +\beta_{32} \times X_2 + \beta_{33} \times X_3 + \beta_{30}\)
  • \(Logit(Y_4)=\ln(\frac{P(Y=4)}{P(Y=5)})=\beta_{41} \times X_1 +\beta_{42} \times X_2 + \beta_{43} \times X_3 + \beta_{40}\)

Ezeket az egyenleteket emeljük \(e\)-re (\(\exp\)), majd szorozzuk fel \(P(Y=5)\)-tel:

  • \(P(Y=1)=P(Y=5) \times \exp(\beta_{11} \times X_1 +\beta_{12} \times X_2 + \beta_{13} \times X_3 + \beta_{10})\)
  • \(P(Y=2)=P(Y=5) \times \exp(\beta_{21} \times X_1 +\beta_{22} \times X_2 + \beta_{23} \times X_3 + \beta_{20})\)
  • \(P(Y=3)=P(Y=5) \times \exp(\beta_{31} \times X_1 +\beta_{32} \times X_2 + \beta_{33} \times X_3 + \beta_{30})\)
  • \(P(Y=4)=P(Y=5) \times \exp(\beta_{41} \times X_1 +\beta_{42} \times X_2 + \beta_{43} \times X_3 + \beta_{40})\)

Egyszerűsítsük egy kicsit a fenti egyenleteket azzal, hogy amiket a regressziós egyenletekből visszakapunk, az valójában mindig az aktuális, \(i\)-edik dummy változónk logitja. Azaz: \[\exp(Logit(Y_i))=\exp(\beta_{i1} \times X_1 +\beta_{i2} \times X_2 + \beta_{i3} \times X_3 + \beta_{i0}), \forall i\]

Ezzel a 4 regressziós egyenletünk:

  • \(P(Y=1)=P(Y=5) \times \exp(Logit(Y_1))\)
  • \(P(Y=2)=P(Y=5) \times \exp(Logit(Y_2))\)
  • \(P(Y=3)=P(Y=5) \times \exp(Logit(Y_3))\)
  • \(P(Y=4)=P(Y=5) \times \exp(Logit(Y_4))\)

Tehát, 4 db regressziós egyenletünk van. A 4 db dummy-val leírható eredményváltozó kategória oddsának logaritmusát (logitját) becslik meg közvetlenül a regressziók. Ha pedig tudnánk, hogy mi a modell szerint a referencia kategóriába esés valószínűsége, azaz esetünkben \(P(Y=5)\), akkor a fenti modellegyenletekből minden egyéb \(P(Y=i)\) valószínűséget is meg tudnánk becsülni!

Gondolkozzunk akkor! Nyilván igaz a valószínűségekre, hogy: \[P(Y=1) + P(Y=2) + P(Y=3) + P(Y=4) + P(Y=5) = 1\]

Ebből pedig következik, hogy: \[P(Y=5) = 1 - P(Y=1) - P(Y=2) - P(Y=3) - P(Y=4)\]

Ezen a ponton pedig felhasználhatjuk a modellegyenletekből a \(P(Y=i)\) értékekre kapott \(P(Y=i)=P(Y=5) \times \exp(Logit(Y_i))\) összefüggéseket és ezeket behelyettesíthetjük a fenti \(P(Y=5)\)-öt megadó formulába: \[P(Y=5) = 1 - [P(Y=5) \times \exp(Logit(Y_1))] - [P(Y=5) \times \exp(Logit(Y_2))] - … - [P(Y=5) \times \exp(Logit(Y_4))]\]

Azaz, egyszerűbb alakban, \(-P(Y=5)\)-öt kiemelve: \[P(Y=5)=1-P(Y=5) \times \sum_{i=1}^4{\exp(Logit(Y_i))}\]

Ha ebből az egyenletből pedig kifejezzük \(P(Y=5)\)-öt, akkor már egy egészen kulturált formulára jutunk. \[P(Y=5)=\frac{1}{1+\sum_i{\exp(Logit(Y_i))}}\]

Ezt az összefüggést felhasználva megvan a többi valószínűség becslése is az \(Y_i, i=1,2,3,4\) dummy változókra felírt modellegyenletek \(P(Y=i)=P(Y=5) \times \exp(Logit(Y_i))\) alapján:

  • \(P(Y=1) = \frac{\exp(Logit(Y_1))}{1 + \sum_i{\exp(Logit(Yi))}}\)
  • \(P(Y=2) = \frac{\exp(Logit(Y_2))}{1 + \sum_i{\exp(Logit(Yi))}}\)
  • \(P(Y=3) = \frac{\exp(Logit(Y_3))}{1 + \sum_i{\exp(Logit(Yi))}}\)
  • \(P(Y=4) = \frac{\exp(Logit(Y_4))}{1 + \sum_i{\exp(Logit(Yi))}}\)

Végül pedig azt az értéket becsüljük a regresszióval \(Y\)-ra, aminek a legnagyobb a \(P(Y=i),i=1,2,3,4,5\) bekövetezési valószínűsége!

Ez alapján megvan az is hogyan lehet a multinomiális logit modellegyenlet együtthatóit, azaz \(\beta_{ij}\)-it megbecsülni maximum likelihood elven!

Hiszen a bináris logit modellben a \(\beta\) együtthatók változtatásával a \(-\ln(L)=-(\sum_{y=1}{ln(P(Y=1))}+\sum_{y=0}{ln(P(Y=0))})\) negatív log-likelihood értéket minimalizáltunk.

Most is ez a helyzet, csak kicsit több tagú lesz a negatív log-likelihood. Hiszen most az \(Y\) eredményváltozónak nem csak 2 db \(0\) és \(1\) kategóriája van, hanem \(5\) db kategóriát kell figyelembe venni: \[-\ln(L)=-(\sum_{y=1}{ln(P(Y=1))}+...+\sum_{y=5}{ln(P(Y=5))})\]

Tehát, a likelihood érték meghatározása egy mintaelemre most is ugyan azon az elven megy, mint a bináris logit esetében. Azt mondjuk, hogy a valóságban \(y=i\) elégedettségi szintbe tartozó programozók a modell szerint az egyenletekből kapott \(P(Y=i)\) valószínűséggel fordulnak elő a világban. Mivel most id lépésben feltesszük, hogy a 210 elemű mintánk FÜGGETLEN megfigyelésekből áll, a teljes minta bekövetkezési valószínűsége az egyéni Likelihood értékek szorzata. Azaz logaritmus szinten a teljes minta log-likelihoodja a logaritmizált egyéni likelihoodok összege!

És kész, ezzel a \(-\ln(L)\) minimalizálással meg is vannak minden egyenletben a \(\beta_{ij}\) együtthatók maximum likelihood becslései! :)

Ez azért is jó hír, mert ha a modell egészének az illeszkedése leírható egy \(-\ln(L)\)-el, mint a bináris logit modell esetében, akkor az azt jelenti, hogy a multinomiális esetben az \(IC\)-k és a pszeudo R-négyzet mutató és a \(\chi^2\)-féle illeszkedésvizsgálat ugyan úgy számíthatók, mint a bináris logit esetében!

3. Multinomiális Logisztikus regresszió R-ben

Na, akkor ezek után csináljunk egy ilyen multinomiális logit modellt R-ben végre!

De mielőtt tényleg belevágunk be kéne állítani egy értelmes referencia kategóriát a nominális eredményvátozónknak. Ez legyen most a legalacsonyabb elégedettségi szint, a Very dissatisfied. Használhatjuk itt is a relevel függvényt, mint a nominális magyarázóváltozóknál tettük:

ProgData$JobSat <- relevel(ProgData$JobSat, ref = "Very dissatisfied")

A függvény, ami megalkotja nekünk a multinomiális logisztikus regressziót multinom névre hallgat és a nnet csomagban lakik. Szóval, telepítsük a csomagot és hivatkozzuk is be R-be:

install.packages("nnet")
library(nnet) # Szokásos módon ne törődjünk az esetleges Warningokkal! :)
## Warning: package 'nnet' was built under R version 4.1.2

A multinom függvény igazából teljesen úgy használható, mint a sima lm függvény. Mivel csak multinomiális logit modelleket csinál ez a függvény, nem kell benne szórakozni a family paraméterrel, mint a glm függvényben kellett.

Első körben magyarázóváltozónk az elégedettségi szinthez (JobSat) legyen az a 3 db, amit az elején is mondtunk: Age, YearsCodePro és MonthlyHuf.

multi_logit <- multinom(JobSat ~ Age + YearsCodePro + MonthlyHuf, data = ProgData)
## # weights:  25 (16 variable)
## initial  value 337.981962 
## iter  10 value 325.590915
## iter  20 value 301.285247
## final  value 301.237278 
## converged
summary(multi_logit)
## Call:
## multinom(formula = JobSat ~ Age + YearsCodePro + MonthlyHuf, 
##     data = ProgData)
## 
## Coefficients:
##                                    (Intercept)          Age YearsCodePro
## Neither satisfied nor dissatisfied  -0.7592164  0.056689359  -0.06215057
## Slightly dissatisfied                0.2416452  0.054695642  -0.05323145
## Slightly satisfied                   1.5794347  0.003601398   0.01897727
## Very satisfied                       2.3996360 -0.056026755   0.09288114
##                                       MonthlyHuf
## Neither satisfied nor dissatisfied -1.389351e-07
## Slightly dissatisfied              -7.399280e-07
## Slightly satisfied                 -3.679717e-07
## Very satisfied                      2.350005e-10
## 
## Std. Errors:
##                                     (Intercept)          Age YearsCodePro
## Neither satisfied nor dissatisfied 1.963987e-13 6.525743e-12 1.816675e-12
## Slightly dissatisfied              2.820433e-13 9.256965e-12 2.478016e-12
## Slightly satisfied                 1.734406e-13 5.810523e-12 1.846768e-12
## Very satisfied                     8.636479e-14 2.930356e-12 9.909756e-13
##                                      MonthlyHuf
## Neither satisfied nor dissatisfied 2.838243e-07
## Slightly dissatisfied              2.877590e-07
## Slightly satisfied                 2.426048e-07
## Very satisfied                     2.276008e-07
## 
## Residual Deviance: 602.4746 
## AIC: 634.4746

No, hát most a summary egy jó nagy táblát dobott vissza, hiszen \(5-1=4\) db egyenletünk van. És így sok helyünk nincs, szóval a magyarázóváltozók parciális teszteléséhez, azaz a \(H_0: \beta_{ij}=0\)-hoz tartozó p-értékeket nem is adja meg a summary, csak az egyes együtthatók standard hibáit.

Ezen kívül leolvasható még az eredményből a modell \(-2\ln(L)\) értéke, mint Residual Deviance a sima bináris logithoz hasonló módon, és az ebből számított \(AIC\).

Egyelőre elemezzük, hogy mit mutatnak meg nekünk az együtthatók, azaz a \(§\beta_{ij}\) értékek! Mivel minden egyenlet egy \(Logit(Y_i)\) transzformációra lett felírva, így az értelmezéshez érdemes lesz őket \(e\)-re emelni az exp függvénnyel, mint a sima bináris esetben tettük. Először kiolvassuk az együtthatókat a coef függvénnyel, aztán jöhet az \(e\)-adik hatvány:

exp(coef(multi_logit))
##                                    (Intercept)       Age YearsCodePro
## Neither satisfied nor dissatisfied    0.468033 1.0583270    0.9397414
## Slightly dissatisfied                 1.273342 1.0562191    0.9481605
## Slightly satisfied                    4.852212 1.0036079    1.0191585
## Very satisfied                       11.019164 0.9455138    1.0973313
##                                    MonthlyHuf
## Neither satisfied nor dissatisfied  0.9999999
## Slightly dissatisfied               0.9999993
## Slightly satisfied                  0.9999996
## Very satisfied                      1.0000000

Amit így kapunk azok megintcsak az egyes \(Y_i\) kategóriák ODDS-aira vett marginális hatásai a magyarázóváltozóknak! Csak most ezekhez az oddsokhoz hozzátartozik, hogy mindig az egyenletekből kihagyott referencia kategóriához képest értendők!!

Nézzük végig mit jelentenek ez alapján az Age \(e\)-ra emelt együtthatói az egyes egyenletekben:

  • \(\beta_{Age1}=1.058\) –> Minden más magyarázóváltozó változatlansága mellett, ha a programozó kora 1 évvel nő, akkor a “Neither satisfied nor dissatisfied” kategóriába tartozás oddsa a “Very dissatisfied” kategóriához nézve \(1.058\)-szorosra nő, vagyis \(5.8\%\)-kal nő.
  • \(\beta_{Age2}=1.056\) –> Minden más magyarázóváltozó változatlansága mellett, ha a programozó kora 1 évvel nő, akkor a “Slightly dissatisfied” kategóriába tartozás oddsa a “Very dissatisfied” kategóriához nézve \(1.056\)-szorosra nő, vagyis \(5.6\%\)-kal nő.
  • \(\beta_{Age3}=1.004\) –> Minden más magyarázóváltozó változatlansága mellett, ha a programozó kora 1 évvel nő, akkor a “Slightly satisfied” kategóriába tartozás oddsa a “Very dissatisfied” kategóriához nézve \(1.004\)-szorosra nő, vagyis \(0.4\%\)-kal nő.
  • \(\beta_{Age4}=0.946\) –> Minden más magyarázóváltozó változatlansága mellett, ha a programozó kora 1 évvel nő, akkor a “Very satisfied” kategóriába tartozás oddsa a “Very dissatisfied” kategóriához nézve \(0.946\)-szorosra csökken, vagyis \(5.4\%\)-kal csökken.

Összességében tehát elmondható, hogy a kor növekedésével (változatlan programozási tapasztalat és fizetés mellett), várhatóan csökken a programozók elégedettsége a munkájukkal. Mivel a “legpozitív hangulatú”, Very satisfied kategória oddsaira hat negatívan a magasabb kor, és a második legpozitívabb Slightly satisfied kategória oddsain lényegében nem változtat a referencia, Very dissatisfied kategóriáoz képest.

Ezek alapján pedig azt tudjuk elmondani a YearsCodePro változóra nézve, hogy minél több éve foglalkozik a programozó profi szinten programozással (változatlan életkor és fizetés mellett), várhatóan annál elégedettebb a munkájával. Ugyanis, a két “legpozitív hangulatú”, Very satisfied és Slightly satisfied kategória oddsaira hat pozitívan (1-nél nagyobb \(e\)-re emelt együtthatókkal) a magas profi programozási tapasztalat a Very dissatisfied kategóriához képest.

Érdekes módon a változatlan életkor és profi programozási tapasztalat mellett úgy tűnik, mintha a munkával vett elégedettségre nem hatna érdemben a fizetés, mivel minden elégedettségi kategóriára nézve az \(e\)-re emelt együtthatója közel \(1\). Azaz, a ceteris paribus magasabb fizetés várhatóan érdemben nem növeli és nem csökkenti az egyéb elégedettségi kategóriák oddsait a Very dissatisfied kategóriához képest.

Szumma szummárum, a fiatal, de még így is nagy programozási tapasztalattal bíró programozók a legelégedettebbek a munkájukkal, fizetéstől függetlenül. Ezt mondja jelenleg nekünk a multinomiális logit modell.

NODE, a magyarázóváltozók közül melyek a szignifikánsak és hogyan tudunk közöttük ez alapján egy fontossági sorrendet felállítani az elégedettségi kategóriák előrejelzését nézve?

3. A magyarázóváltozók parciális tesztelése

Mivel egy magyarázóváltozónak jelenleg mind a négy egeynletben szerepelnek együtthatói, így azt a \(H_0\)-t, ami szerint az adott magyarázóváltozó nem szignifikáns, úgy lehet leafordítani a multinomiális logit modellünk nyelvére, hogy a magyarázóváltozónak mind a négy egyenletben egyszerre tekinthető az együtthatója 0-nak a mintán kívüli világban! Ezzel az alábbi lesz a konkrétan megvizsgált \(H_0\) és \(H_1\) párunk:

  • \(H_0:\) A változó nem szignifikáns = A változó együtthatója minden modellegyenletben = 0
  • \(H_1:\) A változó szignifikáns = A változónak van legalább egy együtthatója valamelyik egyenletben, ami NEM 0

A fenti \(H_0\) és \(H_1\) pároshoz próbafüggvényt úgy kapunk, hogy legyártjuk azt a modellt, amiből a tesztelt változót kihagyjuk, és összenézzük a két modell \(-2\ln(L)\) értékét. A bővebb modell \(-2\ln(L)\) értékét jelöljük \(-2\ln(L_{H_1})\)-nek, míg a nullhipotézisben feltételezett szűkebb modellét \(-2\ln(L_{H_0})\)-nak. A próbafüggvény ezek után e két érték különbsége, ami igaz \(H_0\) esetén \(\chi^2_m\) eloszlású, pont mint a sima bináris logisztikus regresszió \(\chi^2\)-féle illeszkedésvizsgálata esetében:

A különbség annyi, hogy itt \(-2\ln(L_{H_0})\) nem a nullmodell, hanem az egy változóval szűkebb modell \(-2\ln(L)\) értéke. Illetve, az \(m\) szabadságfok a kihagyott \(\beta\)-k száma, azaz gyakorlatilag az egyenletek száma. Hiszen egy magyarázóváltozó minden egyenletben csak egyszer, egy \(\beta\)-val szerepel. Kivéve persze a nominális magyarázóváltozókat, hiszen ők minden egyenletben annyiszor szerepelnek, ahány lehetséges értékük van, mínusz 1. Azaz ők az \(m\) szabadságfokhoz dummy-k száma * egyenletek száma értéket adnak hozzá :)

A \(\chi^2\)-féle illeszkedésvizsgálat most bemutatott alkalmazását a multinomiális logit modellekben az lmtest csomag lrtest függvénye valósítja meg. Az lr elnevezés a likelihood-ratio rövidítésből jön, mivel a log-likelihoodok különbsége az eredeti (nem logaritmált) likelihood-ok hányadosa. :)

Az lrtest függvénybe egyszerűen csak bele kell rakni a vizsgált modell R objektumát, és a változó nevét, amit vizsgálunk. Pl. az Age változó szignifikanciája a modellünkben:

library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
lrtest(multi_logit, "Age")
## # weights:  20 (12 variable)
## initial  value 337.981962 
## iter  10 value 326.036197
## iter  20 value 303.806372
## final  value 303.806030 
## converged
## Likelihood ratio test
## 
## Model 1: JobSat ~ Age + YearsCodePro + MonthlyHuf
## Model 2: JobSat ~ YearsCodePro + MonthlyHuf
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  16 -301.24                     
## 2  12 -303.81 -4 5.1375     0.2735

A p-érték \(27.35\%\), ami elég magas, az Age változó nem szignifikánsnak vehető minden szokásos szignifikanci-szinten. Lássuk akkor a többi változót:

lrtest(multi_logit, "YearsCodePro")
## # weights:  20 (12 variable)
## initial  value 337.981962 
## iter  10 value 308.824226
## iter  20 value 305.462539
## final  value 305.462515 
## converged
## Likelihood ratio test
## 
## Model 1: JobSat ~ Age + YearsCodePro + MonthlyHuf
## Model 2: JobSat ~ Age + MonthlyHuf
##   #Df  LogLik Df  Chisq Pr(>Chisq)  
## 1  16 -301.24                       
## 2  12 -305.46 -4 8.4505     0.0764 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lrtest(multi_logit, "MonthlyHuf")
## # weights:  20 (12 variable)
## initial  value 337.981962 
## iter  10 value 304.669638
## final  value 303.238079 
## converged
## Likelihood ratio test
## 
## Model 1: JobSat ~ Age + YearsCodePro + MonthlyHuf
## Model 2: JobSat ~ Age + YearsCodePro
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  16 -301.24                     
## 2  12 -303.24 -4 4.0016     0.4058

Hm, úgy tűnik, hogy csak a YearsCodePro változó mondható szignifikánsnak, és az is csak \(\alpha=10\%\)-on. Ez szomorú, de hát van még minekünk több lehetséges magyarázóváltozónk, nézzük meg, hogy szuperál a modell, ha minden lehetséges magyarázóváltozót beleveszünk!

3.1. A Teljes Modell magyarázóváltozóinak szignifikanciája

Mielőtt ész nélkül megcsinálnánk a minden magyarázóváltozót használó teljes modellt, nézzünk azért rá a változóink egyszerű leíró statisztikáira a summary-vel:

summary(ProgData)
##       Age          Age1stCode     YearsCodePro      MonthlyHuf     
##  Min.   :21.00   Min.   : 3.00   Min.   : 0.500   Min.   :  61284  
##  1st Qu.:26.00   1st Qu.:10.25   1st Qu.: 3.000   1st Qu.: 511981  
##  Median :31.00   Median :13.00   Median : 6.000   Median : 767652  
##  Mean   :32.17   Mean   :13.92   Mean   : 8.688   Mean   : 888295  
##  3rd Qu.:38.00   3rd Qu.:16.00   3rd Qu.:14.000   3rd Qu.:1023682  
##  Max.   :54.00   Max.   :35.00   Max.   :33.000   Max.   :6649792  
##       Gender                     EdLevel  
##  Man     :200   Bachelor’s degree    :71  
##  NoAnswer:  6   Doctoral degree      :11  
##  Woman   :  4   Master’s degree      :54  
##                 Other                :14  
##                 Secondary school     :19  
##                 Unfinished university:41  
##                                                 Employment 
##  Employed full-time                                  :167  
##  Employed part-time                                  : 10  
##  Independent contractor, freelancer, or self-employed: 33  
##                                                            
##                                                            
##                                                            
##                                 JobSat           OpSys    
##  Very dissatisfied                 :15   Linux-based: 58  
##  Neither satisfied nor dissatisfied:23   MacOS      : 27  
##  Slightly dissatisfied             :39   Windows    :125  
##  Slightly satisfied                :69                    
##  Very satisfied                    :64                    
##

Amit itt rögtön ki tudunk szúrni az a Gender változó. Ezt nem érdemes magyarázóváltozóként alkalmazni a logit modellben, mert lényegében csak egyféle lehetséges értéke van, a férfi. A többi kategória gyakorisága nagyon ritka. Ez nem meglepő amúgy, programozókról beszélünk… :)

A többi lehetséges magyarázóváltozó első szemmelverésre rendben lévőnek tűnik, nem látszik semmi durva outlier pl. a felső kvartilisek és a maximumok viszonyából sem a numerikus változókban.

Akkor futtassuk a modellt, ami a JobSat magyarázatára minden egyéb magyarázóváltozót használ, csak a Gendert nem, majd nézzük meg szépen minden magyarázóváltozó szignifikanciáját az LR-teszttel:

multi_logit <- multinom(JobSat ~ .-Gender, data = ProgData)
## # weights:  75 (56 variable)
## initial  value 337.981962 
## iter  10 value 321.589734
## iter  20 value 284.983944
## iter  30 value 277.192749
## iter  40 value 276.509416
## iter  50 value 276.417469
## iter  60 value 276.358376
## iter  70 value 276.357310
## iter  70 value 276.357307
## iter  70 value 276.357307
## final  value 276.357307 
## converged
lrtest(multi_logit, "Age")
## # weights:  70 (52 variable)
## initial  value 337.981962 
## iter  10 value 322.149564
## iter  20 value 280.906386
## iter  30 value 279.003929
## iter  40 value 278.684730
## iter  50 value 278.553336
## iter  60 value 278.549269
## final  value 278.549244 
## converged
## Likelihood ratio test
## 
## Model 1: JobSat ~ (Age + Age1stCode + YearsCodePro + MonthlyHuf + Gender + 
##     EdLevel + Employment + OpSys) - Gender
## Model 2: JobSat ~ Age1stCode + YearsCodePro + MonthlyHuf + EdLevel + Employment + 
##     OpSys
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  56 -276.36                     
## 2  52 -278.55 -4 4.3839     0.3565
lrtest(multi_logit, "Age1stCode")
## # weights:  70 (52 variable)
## initial  value 337.981962 
## iter  10 value 325.208447
## iter  20 value 279.944675
## iter  30 value 278.308986
## iter  40 value 277.979728
## iter  50 value 277.892274
## iter  60 value 277.873682
## final  value 277.873437 
## converged
## Likelihood ratio test
## 
## Model 1: JobSat ~ (Age + Age1stCode + YearsCodePro + MonthlyHuf + Gender + 
##     EdLevel + Employment + OpSys) - Gender
## Model 2: JobSat ~ Age + YearsCodePro + MonthlyHuf + EdLevel + Employment + 
##     OpSys
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  56 -276.36                     
## 2  52 -277.87 -4 3.0323     0.5524
lrtest(multi_logit, "YearsCodePro")
## # weights:  70 (52 variable)
## initial  value 337.981962 
## iter  10 value 326.104558
## iter  20 value 282.181743
## iter  30 value 280.232862
## iter  40 value 279.824473
## iter  50 value 279.734431
## iter  60 value 279.724666
## final  value 279.724527 
## converged
## Likelihood ratio test
## 
## Model 1: JobSat ~ (Age + Age1stCode + YearsCodePro + MonthlyHuf + Gender + 
##     EdLevel + Employment + OpSys) - Gender
## Model 2: JobSat ~ Age + Age1stCode + MonthlyHuf + EdLevel + Employment + 
##     OpSys
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  56 -276.36                     
## 2  52 -279.73 -4 6.7344     0.1506
lrtest(multi_logit, "MonthlyHuf")
## # weights:  70 (52 variable)
## initial  value 337.981962 
## iter  10 value 304.164320
## iter  20 value 279.270860
## iter  30 value 278.010536
## iter  40 value 277.760426
## iter  50 value 277.687197
## iter  60 value 277.680882
## final  value 277.680831 
## converged
## Likelihood ratio test
## 
## Model 1: JobSat ~ (Age + Age1stCode + YearsCodePro + MonthlyHuf + Gender + 
##     EdLevel + Employment + OpSys) - Gender
## Model 2: JobSat ~ Age + Age1stCode + YearsCodePro + EdLevel + Employment + 
##     OpSys
##   #Df  LogLik Df Chisq Pr(>Chisq)
## 1  56 -276.36                    
## 2  52 -277.68 -4 2.647     0.6185
lrtest(multi_logit, "EdLevel")
## # weights:  50 (36 variable)
## initial  value 337.981962 
## iter  10 value 322.020791
## iter  20 value 298.505987
## iter  30 value 290.090623
## iter  40 value 289.881605
## iter  50 value 289.850645
## final  value 289.850460 
## converged
## Likelihood ratio test
## 
## Model 1: JobSat ~ (Age + Age1stCode + YearsCodePro + MonthlyHuf + Gender + 
##     EdLevel + Employment + OpSys) - Gender
## Model 2: JobSat ~ Age + Age1stCode + YearsCodePro + MonthlyHuf + Employment + 
##     OpSys
##   #Df  LogLik  Df  Chisq Pr(>Chisq)
## 1  56 -276.36                      
## 2  36 -289.85 -20 26.986     0.1356
lrtest(multi_logit, "Employment")
## # weights:  65 (48 variable)
## initial  value 337.981962 
## iter  10 value 321.721269
## iter  20 value 289.500207
## iter  30 value 281.355972
## iter  40 value 280.726635
## iter  50 value 280.665703
## iter  60 value 280.641506
## final  value 280.641226 
## converged
## Likelihood ratio test
## 
## Model 1: JobSat ~ (Age + Age1stCode + YearsCodePro + MonthlyHuf + Gender + 
##     EdLevel + Employment + OpSys) - Gender
## Model 2: JobSat ~ Age + Age1stCode + YearsCodePro + MonthlyHuf + EdLevel + 
##     OpSys
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  56 -276.36                     
## 2  48 -280.64 -8 8.5678     0.3801
lrtest(multi_logit, "OpSys")
## # weights:  65 (48 variable)
## initial  value 337.981962 
## iter  10 value 321.734240
## iter  20 value 285.693568
## iter  30 value 281.354767
## iter  40 value 280.980166
## iter  50 value 280.901543
## iter  60 value 280.898815
## final  value 280.898803 
## converged
## Likelihood ratio test
## 
## Model 1: JobSat ~ (Age + Age1stCode + YearsCodePro + MonthlyHuf + Gender + 
##     EdLevel + Employment + OpSys) - Gender
## Model 2: JobSat ~ Age + Age1stCode + YearsCodePro + MonthlyHuf + EdLevel + 
##     Employment
##   #Df  LogLik Df Chisq Pr(>Chisq)
## 1  56 -276.36                    
## 2  48 -280.90 -8 9.083     0.3353

Ajjaj az a helyzet, hogy egy magyarázóváltozónk sem szignifikáns! Ez azt sugallja, hogy az aktuális modellünk NEM tér el szignifikánsan a nullmodelltől.

Gyártsuk is le a nullmodellt, és győződjünk meg erről az LR-teszt azon verziójával, amit akkor kapunk, ha azt gondoljuk a \(H_0\)-ban, hogy minden \(\beta\) minden egyenletben nullának vehető a mintán kívüli világban \(\forall\beta_{ij}=0\).
Ez a folyamat úgy történik pont, mint a sima bináris logit modellnél történt:

# Nullmodell = Magyarázóváltpzók nélküli modell építése
multi_logit_null <- multinom(JobSat ~ 1, data = ProgData)
## # weights:  10 (4 variable)
## initial  value 337.981962 
## iter  10 value 308.954691
## iter  10 value 308.954691
## final  value 308.954691 
## converged
# Khi-négyzet próba a nullmodellre és az aktuális modellre
anova(multi_logit, multi_logit_null)
##                                                                                             Model
## 1                                                                                               1
## 2 (Age + Age1stCode + YearsCodePro + MonthlyHuf + Gender + EdLevel + Employment + OpSys) - Gender
##   Resid. df Resid. Dev   Test    Df LR stat.   Pr(Chi)
## 1       836   617.9094           NA       NA        NA
## 2       784   552.7146 1 vs 2    52 65.19477 0.1034081

A p-érték \(10.34\%\), tehát éppenhogy, de még \(\alpha=10\%\)-on sem lehet elvetni azt a \(H_0\)-t, hogy a modellünkben minden \(\beta_{ij}=0\).

Na, ilyet aztán még nem láttunk! Most akkor mivan?

Nos a megfejtésünk az, hogy túl sok paraméterünk, azaz túl sok \(\beta\)-nk van a modellben a megfigyelések számához képest!!

Általános ajánlás, hogy egy regressziós modellben az elemszám, azaz \(n\) legyen nagyobb, mint a paraméterek számának ötszöröse: \(n > p \times 5\) vagy \(\frac{n}{p} > 5\).

Ez azért van, mert minél több \(\beta\)-ja, azaz paramétere van a regressziónknak annál jobban tud ráilleszkedni a megfigyelt adatokra.
Hiszen annál több lehetősége van a negatív log-likelihoodot, tehát \(-\ln(L)\)-t vagy az \(SSE\)-t minimalizálni. Azaz a modell túlilleszkedik.

Ez az a jelenség, ami ellen védekezünk modellszelekciónál a Wald-próbával és az \(IC\)-k használatával, mivel ezek az eszközök büntetik a feleslegesen bevont extra magyarázóváltozókat, azaz \(\beta\)-kat. De a túlillesztés egy modellben önmagában is jelentkezhet gondot, ha az elemszámhoz képest túl sok a \(\beta\).

Ez a jelenség, hogy egy modell önmagában is túlillesztett legyen az elemszámhoz képest a sima lineáris és bináris logisztikus regresszióknál nem gyakran áll fenn, mivel nagyon sok magyarázóváltozó kell hozzá. Pl. a mostani \(210\)-es elemszámunkhoz \(\frac{n}{5}=\frac{210}{5}=42\) paraméter kell, hogy a határt elérjük.

Ellenben most a multinomiális regressziónál van \(4\) egyenletünk, ami alapból minden \(\beta\)-t négyszerez, és van ugye nominális magyarázóváltozó is, ami sok-sok dummy-t jelent…Nézzük meg hány \(\beta\) is van most a modellben:

length(coef(multi_logit))
## [1] 56

Ínnye, \(56>42\)… :( Tehát a modellünkben túl sok paraméter van a megfigyelések számához képest. Ilyenkor pedig minden hipotézisvizsgálat elkezd “túl óvatos lenni”, mondván egy ennyire túlillesztett modellben semmi változó hatása nem lehet szignifikáns a mintán kívüli világan!

Hogyan tegyük helyre a kérdést? Hát csökkentsük az egyenletek számát! Ezt hogy lehet elérni? A ritka eredményváltozó kategóriák összevonásával!
Ezekben a ritka gyakoriságú kategóriákban valószínűleg egyébként is elég béna a modell magyarázóereje, szóval valószínűleg nem kár az összevonásukért. :)

4. Az eredményváltozó ritka kategóriáinak összevonása

Az összevonandó kategóriák kijelöléséhez lássuk először az eredményváltozónk egyes értékeinek gyakoriságait!

summary(ProgData$JobSat)
##                  Very dissatisfied Neither satisfied nor dissatisfied 
##                                 15                                 23 
##              Slightly dissatisfied                 Slightly satisfied 
##                                 39                                 69 
##                     Very satisfied 
##                                 64

Úgy látszik, hogy a 3 legritkább kategória tök logikusan egybevonható, mivel ezek a “nem satisfied” elégedettségi szintek.
Tegyük is meg ezt a lépést!

Először character adattípusúvá kell konvertálnunk a JobSat változót, hogy az értékkészletét meg tudjuk buherálni, majd elvégezzük a 3 kérdéses kategória átnevezését Not satisfied-á. Végül újra factor adattípusra lépünk vissza, hogy továbbra is tudjon a multinomiális modellünk mit kezdeni a változóval, és beállítjuk az új, “Not satisfied” kategóriát referenciának.
Ez a szép hosszadalmas lépéssorozat R kód szinten a következőképp néz ki:

ProgData$JobSat <- as.character(ProgData$JobSat)
ProgData$JobSat[ProgData$JobSat == "Neither satisfied nor dissatisfied"] <- "Not satisfied"
ProgData$JobSat[ProgData$JobSat == "Slightly dissatisfied"] <- "Not satisfied"
ProgData$JobSat[ProgData$JobSat == "Very dissatisfied"] <- "Not satisfied"
ProgData$JobSat <- as.factor(ProgData$JobSat)

ProgData$JobSat <- relevel(ProgData$JobSat, ref = "Not satisfied")

summary(ProgData$JobSat)
##      Not satisfied Slightly satisfied     Very satisfied 
##                 77                 69                 64

Szuper, a summary eredménye alapján, tök jól összevonásra került a 3 ritka kategória. Most már az eredményváltozónak csak 3, többi kevésbé hasonló gyakoriságú csoportja van, így a multinom modellnek csak \(3-1=2\) egyenlete lesz mindösszesen.

Lássuk is akkor az új modellt, és paramétereinek (azaz \(\beta\)-inak) számát!

multi_logit <- multinom(JobSat ~ .-Gender, data = ProgData)
## # weights:  45 (28 variable)
## initial  value 230.708581 
## iter  10 value 217.156107
## iter  20 value 208.119464
## iter  30 value 208.107978
## iter  30 value 208.107977
## iter  30 value 208.107977
## final  value 208.107977 
## converged
length(coef(multi_logit))
## [1] 28

Na már csak \(28\) paraméterünk van, ami kisebb, mint \(\frac{n}{5}=\frac{210}{5}=42\). Hurrá!

Akkor lásszuk van-e már szignifikáns változónk az LR-tesztek szerint:

lrtest(multi_logit, "Age")
## # weights:  42 (26 variable)
## initial  value 230.708581 
## iter  10 value 213.243280
## iter  20 value 210.237437
## final  value 210.231796 
## converged
## Likelihood ratio test
## 
## Model 1: JobSat ~ (Age + Age1stCode + YearsCodePro + MonthlyHuf + Gender + 
##     EdLevel + Employment + OpSys) - Gender
## Model 2: JobSat ~ Age1stCode + YearsCodePro + MonthlyHuf + EdLevel + Employment + 
##     OpSys
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  28 -208.11                     
## 2  26 -210.23 -2 4.2476     0.1196
lrtest(multi_logit, "Age1stCode")
## # weights:  42 (26 variable)
## initial  value 230.708581 
## iter  10 value 212.536379
## iter  20 value 209.360065
## final  value 209.353919 
## converged
## Likelihood ratio test
## 
## Model 1: JobSat ~ (Age + Age1stCode + YearsCodePro + MonthlyHuf + Gender + 
##     EdLevel + Employment + OpSys) - Gender
## Model 2: JobSat ~ Age + YearsCodePro + MonthlyHuf + EdLevel + Employment + 
##     OpSys
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  28 -208.11                     
## 2  26 -209.35 -2 2.4919     0.2877
lrtest(multi_logit, "YearsCodePro")
## # weights:  42 (26 variable)
## initial  value 230.708581 
## iter  10 value 214.245861
## iter  20 value 211.517476
## final  value 211.506235 
## converged
## Likelihood ratio test
## 
## Model 1: JobSat ~ (Age + Age1stCode + YearsCodePro + MonthlyHuf + Gender + 
##     EdLevel + Employment + OpSys) - Gender
## Model 2: JobSat ~ Age + Age1stCode + MonthlyHuf + EdLevel + Employment + 
##     OpSys
##   #Df  LogLik Df  Chisq Pr(>Chisq)  
## 1  28 -208.11                       
## 2  26 -211.51 -2 6.7965    0.03343 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lrtest(multi_logit, "MonthlyHuf")
## # weights:  42 (26 variable)
## initial  value 230.708581 
## iter  10 value 211.419515
## iter  20 value 208.515045
## final  value 208.507922 
## converged
## Likelihood ratio test
## 
## Model 1: JobSat ~ (Age + Age1stCode + YearsCodePro + MonthlyHuf + Gender + 
##     EdLevel + Employment + OpSys) - Gender
## Model 2: JobSat ~ Age + Age1stCode + YearsCodePro + EdLevel + Employment + 
##     OpSys
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  28 -208.11                     
## 2  26 -208.51 -2 0.7999     0.6704
lrtest(multi_logit, "EdLevel")
## # weights:  30 (18 variable)
## initial  value 230.708581 
## iter  10 value 222.483208
## iter  20 value 216.669008
## final  value 216.668572 
## converged
## Likelihood ratio test
## 
## Model 1: JobSat ~ (Age + Age1stCode + YearsCodePro + MonthlyHuf + Gender + 
##     EdLevel + Employment + OpSys) - Gender
## Model 2: JobSat ~ Age + Age1stCode + YearsCodePro + MonthlyHuf + Employment + 
##     OpSys
##   #Df  LogLik  Df  Chisq Pr(>Chisq)  
## 1  28 -208.11                        
## 2  18 -216.67 -10 17.121    0.07172 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lrtest(multi_logit, "Employment")
## # weights:  39 (24 variable)
## initial  value 230.708581 
## iter  10 value 219.610428
## iter  20 value 210.379812
## final  value 210.371184 
## converged
## Likelihood ratio test
## 
## Model 1: JobSat ~ (Age + Age1stCode + YearsCodePro + MonthlyHuf + Gender + 
##     EdLevel + Employment + OpSys) - Gender
## Model 2: JobSat ~ Age + Age1stCode + YearsCodePro + MonthlyHuf + EdLevel + 
##     OpSys
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  28 -208.11                     
## 2  24 -210.37 -4 4.5264     0.3394
lrtest(multi_logit, "OpSys")
## # weights:  39 (24 variable)
## initial  value 230.708581 
## iter  10 value 218.451457
## iter  20 value 211.273584
## final  value 211.271181 
## converged
## Likelihood ratio test
## 
## Model 1: JobSat ~ (Age + Age1stCode + YearsCodePro + MonthlyHuf + Gender + 
##     EdLevel + Employment + OpSys) - Gender
## Model 2: JobSat ~ Age + Age1stCode + YearsCodePro + MonthlyHuf + EdLevel + 
##     Employment
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  28 -208.11                     
## 2  24 -211.27 -4 6.3264     0.1761

No, egész tűrhető a dolog! A YearsCodePro szignifikáns változó lett \(\alpha=5\%\)-on, míg az EdLevel \(\alpha=5\%\)-on. A két igazán nem szignifikáns változónk \(20\%\) feltti p-értékekkel az Age1stCode, az Employment és a MonthlyHuf.

Vegyük ki ezt a 3 elég haszontalannak tűnő változót a modellből, és nézzik meg az \(IC\)-k, és az LR-teszt szerint javult-e ezzel a modell. Most az LR-teszt \(\chi^2\)-próbájához a p-értéket az anova függvénnyel számoljuk R-ben, mivel ezzel a függvénnyel lehet úgy használni a tesztet, hogy két (akár több, mint egy változóban eltérő modellt) hasonlítunk össze. Ugye az lrtest függvényben mindig csak egy adott változó \(\beta\)-inak kinullázását vizsgáltuk.
Ne feledjük, hogy az LR-tesztben a \(H_0\) továbbra is az, hogy a szűkített modell a preferált!

multi_logit_szuk <- multinom(JobSat ~ Age + YearsCodePro + EdLevel + OpSys, data = ProgData)
## # weights:  33 (20 variable)
## initial  value 230.708581 
## iter  10 value 212.437840
## iter  20 value 211.910561
## final  value 211.910444 
## converged
# minden eszköz a szűkített modellt preferálja
AIC(multi_logit, multi_logit_szuk)
##                  df      AIC
## multi_logit      28 472.2160
## multi_logit_szuk 20 463.8209
BIC(multi_logit, multi_logit_szuk)
##                  df     BIC
## multi_logit      28 565.935
## multi_logit_szuk 20 530.763
anova(multi_logit, multi_logit_szuk)
##                                                                                             Model
## 1                                                            Age + YearsCodePro + EdLevel + OpSys
## 2 (Age + Age1stCode + YearsCodePro + MonthlyHuf + Gender + EdLevel + Employment + OpSys) - Gender
##   Resid. df Resid. Dev   Test    Df LR stat.   Pr(Chi)
## 1       400   423.8209           NA       NA        NA
## 2       392   416.2160 1 vs 2     8 7.604933 0.4729804

Oké: mindkét \(IC\) a szűkített modellnél kisebb, és az LR-teszt p-értéke is jó magas, azaz minden szokásos \(\alpha\)-n elfogadhatjuk \(H_0\)-t. Szumma, szummárum, a szűkített modellünk jobban illeszkedik az adatokra!!

Lássunk hát pár együttható értelmezést a szűkített modellből, természetesen \(e\) hatványra emelve!

exp(coef(multi_logit_szuk))
##                    (Intercept)       Age YearsCodePro EdLevelDoctoral degree
## Slightly satisfied   0.8905682 0.9835542     1.051245              0.3674308
## Very satisfied      19.6831061 0.8724643     1.184561              3.5242892
##                    EdLevelMaster’s degree EdLevelOther EdLevelSecondary school
## Slightly satisfied               1.228248     1.638280               2.1607138
## Very satisfied                   2.012958     0.817426               0.7179254
##                    EdLevelUnfinished university OpSysMacOS OpSysWindows
## Slightly satisfied                     3.814164  0.5008606    0.7802951
## Very satisfied                         2.742998  0.4183552    0.3739768

Lássuk, a két relatíve szignifikánsnak vehető változót: a YearsCodePro-t és az EdLevelt!

  • A YearsCodePro esetében láthatjuk, hogy minden egyéb változó változatlansága mellett +1 év programozói tapasztalat \(5\%\)-kal emeli a Slightly satisfied, míg \(18\%\)-kal emeli a Very satisfied kategóriába kerülés oddsait a Not satisfied szinthez képest!
  • Az EdLevel tekintetében pl. azt láthatjuk, hogy az egyetemet nem befejező programozók a legelégedettebbek munkájukkal, hiszen minden egyéb változó változatlansága mellett, ha a programozó nem végzett egyetemet, akkor a két satidfied kategória oddsai \(3.81\)-szeresre, illetve \(2.74\)-szeresre nőnek a referencia kategóriában szereplő Bachelor diplomával rendelkezőkhöz képest.
  • Bár a legmagasabb EdLevel szerinti oddsok a Very satisfied kategóriában a doktorival bíró programozóknak vannak a Bachelor diplomásokhoz képest. Ugyanakkor, a doktorisok \(1\) alatti együtthatója a Slightly satisfied kategóriában arra utal, hogy valószínűbb, hogy ők Not satisfiedok, mint Slightly satisfiedok a Bachelor diplomásokhoz képest. Szóval a doktorisok vagy nagyon elégedettek vagy nem elégedettek. Mintha nem lenne náluk középút. :)

5. Klasszifikációs Mátrix és Pszeudo R-négyzet

No, akkor próbáljuk megtudni, hogy úgy összegészében mennyire jó előrejelzéseket készít a modellünk a programozók elégedettségére!

Ehhez kérjük le a modellünk becsléseit mind a \(210\) megfigyelésünk JobSat értékére. Mivel a modell a predikciók során csak annyit csinál, hogy mindenkire kiszámolja mindhárom kategória bekövetkezési valószínűségét, és kiválasztja a legnagyobb valószínűségű kategóriát, így az R predict függvényével simán le lehet kérni közvetlenül a becslést az eredményváltozó legvalószínűbb kategóriájára. Ezt a predict függvényben egy type="class beállítással érhetjük el.
Az első két paraméter ugyan az a függvényben, amit már megszoktunk a bináris logitnál is: a modell, amiből előrejelzéseket kérünk és a data frame, aminek a soraira kellenek az előrejelzések:

ProgData$JobSatPredicted <- predict(multi_logit_szuk, newdata = ProgData, type = "class")
str(ProgData)
## tibble [210 x 10] (S3: tbl_df/tbl/data.frame)
##  $ Age            : num [1:210] 31 27 22 31 24 33 43 23 38 22 ...
##  $ Age1stCode     : num [1:210] 12 13 11 14 15 24 10 20 12 9 ...
##  $ YearsCodePro   : num [1:210] 8 6 5 3 2 1 20 0.5 15 3 ...
##  $ MonthlyHuf     : num [1:210] 788080 1023642 1702347 511981 337703 ...
##  $ Gender         : Factor w/ 3 levels "Man","NoAnswer",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ EdLevel        : Factor w/ 6 levels "Bachelor’s degree",..: 6 3 1 3 5 1 3 4 2 5 ...
##  $ Employment     : Factor w/ 3 levels "Employed full-time",..: 1 1 3 1 1 3 3 1 1 3 ...
##  $ JobSat         : Factor w/ 3 levels "Not satisfied",..: 1 2 3 2 2 2 2 2 2 2 ...
##  $ OpSys          : Factor w/ 3 levels "Linux-based",..: 3 3 3 1 3 1 3 3 3 2 ...
##  $ JobSatPredicted: Factor w/ 3 levels "Not satisfied",..: 2 3 1 1 2 1 3 1 3 1 ...

nagyon jó, megvannak a JobSat becslések! Ezekből már a bináris esehet hasonló módon jöhet a klasszifikációs mátrix a table függvénnyel. Most mentsük is el ezt a klasszifikációs mátrixot egy külön R objektumba, mielőtt megnéznénk az értékeit!

klassz_matrix <- table(ProgData[,c("JobSat", "JobSatPredicted")])
klassz_matrix
##                     JobSatPredicted
## JobSat               Not satisfied Slightly satisfied Very satisfied
##   Not satisfied                 46                 18             13
##   Slightly satisfied            29                 21             19
##   Very satisfied                21                 11             32

Nagyon jó! Látjuk az eredményből, hogy van pl. \(46\) Not satisfied programozónk, akiket helyesen azonosít a modell, míg van \(13\) valóságban Not satisfied programozónk, akiket a modell tévesen Very satisfiednak prediktált.

A klasszifikációs mátrixból ugyan úgy számolhatunk Accuracy, Precision és Recall mutatókat, mint bináris esetben. Annyi lesz az extra, hogy Precisionból és Recallból értelemszerűen \(3\) érték lesz a \(3\) kategóriára és nem csak \(2\).

Lássuk az Accuracy-t! Itt ugye a mátrix átlójában szereplő helyes találatok számát kell összeadni, és elosztani a teljes elemszámmal, azaz a mátrix elemeinek teljes összegével. Ezt a műveletet könnyen el tudjuk végezni a klassz_matrix objektumon. Hiszen be tudjuk vetni az átlóban lakó elemek kilvasására a diag függvényt:

round((sum(diag(klassz_matrix))/sum(klassz_matrix))*100,2) # szorzunk 100-al és 2 tizedesre is kerekítünk
## [1] 47.14

Azt mondja az eredmény, hogy a megfigyelések majdnem felét, \(47.14\%\)-át helyesen sorolta be a modell. Tűrhető. “Nem jó, de nem is tragikus.”

Lássuk a Precision mutatókat! Ebben a mátrixban most az oszlopokban vannak a prediktált kategóriák. Azaz minden kategóriában a megfelelő átló értéket kell osztani az adott kategória oszlopösszegével. Tehát, pl. a “Slightly satisfied” kategóriában a Precision a \(\frac{21}{18+21+11}\) műveletet jelenti. Ezt R-ben gy tudjuk könnyen kicsikarni, hogy az átlóból kiválasztjuk a 2. elemet a számálóba, majd a nevezőben összeadjuk a 2. oszlop értékeit. A 2. oszlopát a mátrixnak pedig klassz_matrix[,2] módon lehet kiválasztani.

Alkalmazzuk a fenti elvet mindhárom eredményváltozó kategóriánkra!

# Not satisfied
round((sum(diag(klassz_matrix)[1])/sum(klassz_matrix[,1]))*100,2)
## [1] 47.92
# Slightly satisfied
round((sum(diag(klassz_matrix)[2])/sum(klassz_matrix[,2]))*100,2)
## [1] 42
# Very satisfied
round((sum(diag(klassz_matrix)[3])/sum(klassz_matrix[,3]))*100,2)
## [1] 50

Ezek szerint a legpontosabban a modellünk a Very satisfied kategóriát jelzi előre: a Very satisfied kategóriájú becslései az esetek felében, \(50\%\)-ában bejönnek!

Végül nézzük a Recallokat! Most a valós eredményváltozó kategóriák a mátrix soraiban vannak, így csak annyir kell változtatni a Precision esetében használt képleteken, hogy a helyes találatok számát az adott kategóriában (ami az átlóban van) nem az adott kategória oszlopösszegéhez, hanem sorösszegéhez kell viszonyítani a nevezőben! A klasszifikációs mátrixunk \(i\)-edik sorát pedig klassz_matrix[i,] módon lehet kiválasztani.

Alkalmazzuk a fenti elvet mindhárom eredményváltozó kategóriánkra!

# Not satisfied
round((sum(diag(klassz_matrix)[1])/sum(klassz_matrix[1,]))*100,2)
## [1] 59.74
# Slightly satisfied
round((sum(diag(klassz_matrix)[2])/sum(klassz_matrix[2,]))*100,2)
## [1] 30.43
# Very satisfied
round((sum(diag(klassz_matrix)[3])/sum(klassz_matrix[3,]))*100,2)
## [1] 50

Az eredmények alapján a modellünk igen jó a Not satisfied kategóriák felderítésében, a valóságban Not satisfied programozók majdnem \(60\%\)-át helyesen azonosította.

Mivel most nincs cut-off értékünk, hiszen egyszerűen a legnagyobb valószínűséggel bekövetkező kategóriát prediktálja a modell, így a klasszifikációs mátrixnak ez az egy verziója van multinomiális esetben. Nincs lehetőségünk most arra, hogy a Precison/Recall értékekkel játszunk a modell megváltoztatása nélkül.

5.1. A pszeudo R-négyzet multinomiális esetben

A modell becslési pontosságát meg tudjuk mérni multinomiális esetben is a McFadden-féle pszeudo R-négyzet mutatóval. Hiszen láttuk a 2. fejezetben, hogy a multinomiális modell egészének az illeszkedése is leírható egy \(-\ln(L)\)-el, mint a bináris logit modell esetében. Ez nyilván igaz a nullmodellre is.

Szóvak a pszeudo \(R^2\) alábbi képlete a multinomiális esetben is alkalmazható:

A dolog szépséghibája, hogy a bináris esetben használt függvény a a pszeudo \(R^2\)-re a DescTools csomagból multinomiális esetben nem működik, így manuálisan kell kiszámítanunk az értéket.

De ezt meg is tudjuk tenni minden további nélkül, mivel a deviance függvénnyel ki tudjuk szedni a multinomiális modellekből is a \(-2\ln(L)\) értékeket. A kettes szorzó a fenti képlet szempontjából pedig nem számít, hiszen kiegyszerűsödik az osztás során.

Számoljuk is akkor ki ezt a csudiszép pszeudo \(R^2\)-et!

# -2*ln(L) a viszgálttárgymodellre
neg_lnL_targy <- deviance(multi_logit_szuk)

# nullmodell legyártása, és a -2*ln(L) értékének elmentése
multi_null <- multinom(JobSat ~ 1, data = ProgData)
## # weights:  6 (2 variable)
## initial  value 230.708581 
## final  value 230.097698 
## converged
neg_lnL_null <- deviance(multi_null)

# pszeudo R-négyzet számolása
(neg_lnL_null-neg_lnL_targy)/neg_lnL_null
## [1] 0.07904144

Az eredményünk alapján a szűkített modellünk \(7.9\%\)-kal ad jobb becsléseket a JobSat értékekre, mint a magyarázóváltozók nélküli nullmodell. Ez egy \(10\%\) alatti R-négyzet, szóval a modellünk magyarázóereje ez alapján gyengének minősíthető. De legalább több, mint a semmi. :)
Plusz, azt láttuk a klasszifikációs mtárixból, hogy pl. a valóságban Not satisfied programozók azonosításában egész jól szuperál a modell. Szóval, gyenge, de valamire mégiscsak jó. :)

6. Marginális hatások vizualizációja Multinomiális Logit modellben

Záróakkordként nézzük meg hogyan lehet a multinomiális modellben a magyarázóváltozók valószínűségekre gyakorolt marginális hatásait vizualizálni.
Mivel lássuk be, több egyenlet együtthatóiból leolvasgatni az ODDSokra gyakorolt marginális hatásokat kellően agymegsüllyesztő vállalkozás. :)

Lássunk valami emészthetőbbet!

Először nézzük meg hogy az egyes JobSat kategóriáknak milyen bekövetkezési valószínűséget ad a szűkített multinomiális modellünk. Első körben egy külön, BecsultVszek névre hallgató data frame-be rakjuk el ezeket a valószínűségeket. Előcsalni őket a multinomiális modellből a predict függvénnyel lehet type = "probs" beállítással:

BecsultVszek <- predict(multi_logit_szuk, newdata = ProgData, type = "probs")
str(BecsultVszek)
##  num [1:210, 1:3] 0.222 0.362 0.404 0.368 0.418 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : chr [1:210] "1" "2" "3" "4" ...
##   ..$ : chr [1:3] "Not satisfied" "Slightly satisfied" "Very satisfied"

Nagyon jó, egy \(3\) oszlopos data frame lett az eredmény a \(3\) eredményváltozó kategória bekövetkezési valószínűségeivel mind a \(210\) megfigyelt programozóra.

Kössük össze ezeket a valószínűségeket az eredeti ProgData data frame-el egy új, ProgDataWithProbs nevű data frame-ben régi barátunkkal, a cbind függvénnyel:

ProgDataWithProbs <- cbind(ProgData, BecsultVszek)
str(ProgDataWithProbs)
## 'data.frame':    210 obs. of  13 variables:
##  $ Age               : num  31 27 22 31 24 33 43 23 38 22 ...
##  $ Age1stCode        : num  12 13 11 14 15 24 10 20 12 9 ...
##  $ YearsCodePro      : num  8 6 5 3 2 1 20 0.5 15 3 ...
##  $ MonthlyHuf        : num  788080 1023642 1702347 511981 337703 ...
##  $ Gender            : Factor w/ 3 levels "Man","NoAnswer",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ EdLevel           : Factor w/ 6 levels "Bachelor’s degree",..: 6 3 1 3 5 1 3 4 2 5 ...
##  $ Employment        : Factor w/ 3 levels "Employed full-time",..: 1 1 3 1 1 3 3 1 1 3 ...
##  $ JobSat            : Factor w/ 3 levels "Not satisfied",..: 1 2 3 2 2 2 2 2 2 2 ...
##  $ OpSys             : Factor w/ 3 levels "Linux-based",..: 3 3 3 1 3 1 3 3 3 2 ...
##  $ JobSatPredicted   : Factor w/ 3 levels "Not satisfied",..: 2 3 1 1 2 1 3 1 3 1 ...
##  $ Not satisfied     : num  0.222 0.362 0.404 0.368 0.418 ...
##  $ Slightly satisfied: num  0.525 0.266 0.25 0.28 0.465 ...
##  $ Very satisfied    : num  0.253 0.372 0.345 0.353 0.117 ...

Oké, akkor mevannak a becsült valószínűségekkel kiegészített adatok is! Erre az egészre azért van szükség, mert a marginális hatások ábrázolásához ezt a data frame-t kicsit átdolgozzuk, és nem akarjuk, hogy véletlen az eredeti tábla elromoljon.

Első körben mondjuk két legszignifikánsabbnak vehető magyarázóváltozó valószínűségekre gyakorolt hatásait vizsgáljuk meg: a YearsCodePro és az EdLevel hatásait.
Szűkítsük le a ProgDataWithProbs tábla oszlopait egy új, temp_tabla objektumban úgy, hogy csak ezt a két változót és a három valószínűséget jelölő oszlopot tartalmazza:

temp_tabla <- ProgDataWithProbs[,c("YearsCodePro", "EdLevel",
                                   "Not satisfied","Slightly satisfied",
                                   "Very satisfied")]

# nézzünk bele az eredményül kapott temp_tabla első 6 sorába
head(temp_tabla)
##   YearsCodePro               EdLevel Not satisfied Slightly satisfied
## 1            8 Unfinished university     0.2220232          0.5249303
## 2            6       Master’s degree     0.3616731          0.2662595
## 3            5     Bachelor’s degree     0.4043842          0.2504966
## 4            3       Master’s degree     0.3677947          0.2795227
## 5            2      Secondary school     0.4175159          0.4653321
## 6            1     Bachelor’s degree     0.5555360          0.3009031
##   Very satisfied
## 1      0.2530466
## 2      0.3720674
## 3      0.3451191
## 4      0.3526827
## 5      0.1171520
## 6      0.1435609

Ok, úgy néz ki jók vagyunk! Csak az a \(2+3=5\) oszlop szerepel itt, amit akartunk!

Ezen a ponton alakítsuk át a becsült valószínűségekkel kiegészített táblát úgy, hogy az egyes elégedettségi kategóriák valószínűségei ne oszlopokban, hanem külön sorokban szerepeljenek.
Tehát így minden programozónak \(3\) sora lesz a \(3\) elégedettségi kategória miatt. Így az új táblának összesen \(210 \times 3=630\) sora lesz.

A művelethez a reshape2 csomag melt függvényét vetjük be. Telepítsük is a csomagot és hivatkozzuk be a függvényeit az R-be:

install.packages("reshape2")
library(reshape2) # Szokásos módon ne törődjünk az esetleges Warningokkal! :)
## Warning: package 'reshape2' was built under R version 4.1.2

A melt függvényt magát úgy alkalmazzuk, hogy először megadjuk a data frame-t, amin dolgozunk, ez most nekünk a temp_tabla. Utána megadjuk az id.vars paraméterben, hogy mely változókat NE törje be több sorba. Ez a két vizsgált magyarázóváltozónk lesz: YearsCodePro és EdLevel. Végül pedig egy value.name paraméterben megadjuk, hogy mi legyen annak az oszlopnak a neve, amiben a \(3\) kategória valószínűségei egységesen tárolásra kerülnek. Erre logikus név most nekünk a Probability.

temp_tabla <- melt(temp_tabla, id.vars = c("YearsCodePro", "EdLevel"), value.name = "Probability")

# nézzük meg hány sora van az új temp_tabla-nak
nrow(temp_tabla)
## [1] 630
# nézzünk bele az eredményül kapott temp_tabla első 6 sorába
head(temp_tabla)
##   YearsCodePro               EdLevel      variable Probability
## 1            8 Unfinished university Not satisfied   0.2220232
## 2            6       Master’s degree Not satisfied   0.3616731
## 3            5     Bachelor’s degree Not satisfied   0.4043842
## 4            3       Master’s degree Not satisfied   0.3677947
## 5            2      Secondary school Not satisfied   0.4175159
## 6            1     Bachelor’s degree Not satisfied   0.5555360

Oké, a várakozásainknak megfelelően megvan a \(210 \times 3=630\) soros új temp_tabla, amiben külön soron vannak a \(3\) eredményváltozó kategória multinomiális modell szerinti valószínűségei minden megfigyelésre.

Ebből a táblából pedig már simán lehet ábrázolni a 3 elégedettségi kategória valószínűségeit a YearsCodePro változó függvényében vonaldiagramon! Csak ggplot2 kell hozzá!

Az ábrán x koordináták a YearsCodePro értékek, y koordináták a valószínűségek, és külön színt kapnak az egyes eredményváltozó, azaz \(JobSat\) kategóriák, amiket a variable oszlopban tárolunk a temp_tabla data frame-ben:

library(ggplot2)

ggplot(temp_tabla, aes(x = YearsCodePro, y = Probability, colour = variable)) +
  geom_line(size = 1) # vonal vastagítása a size-al

Láthatjuk, hogy kb \(15\) év munkatapasztalat után egyértelműen magasabban futnak a kék és zöld vonalak, azaz a Very satisfied és Slightly satisfied kategóriák bekövetkezései valószínűségei. Tehát, hasonló dolgot látunk valószínűségek szintjén is, amit a YearsCodePro együtthatói is elmondtak az oddsokról: a tapasztalat növekedésével nőnek a pozitív hangulatú kategóriák (Very satisfied és Slightly satisfied) oddsai a Not satisfied kategóriához képest.

Az ábrát fel is turbózhatjuk! Egy facet_grid nevű extra réteg alkalmazásával megnézetjük a YearsCodePro marginális hatásait EdLevel kategóriánként külön-külön.
A lentó kódban az EdLevel ~ . felírás technikai szükségszerűség, míg a scales="free" beállítás azt eredményezi, hogy minden EdLevel esetén külön skálázása lesz a diagram tengelyeinek.

ggplot(temp_tabla, aes(x = YearsCodePro, y = Probability, colour = variable)) + geom_line(size = 1) +
  facet_grid(EdLevel ~ ., scales = "free")

Érdekes, hogy a valószínűségekre a YearsCodePro-nak alapvetően a Bachelor’s degree és Master’s degree végzettségek esetén van hatása. A többi iskolai végzettség esetén nem annyira.
Pl. a Unfinished university végzettség esetén kb. minden YearsCodePro érték mellett a két elégedett kategória (Very satisfied és Slightly satisfied) a legnagyobb valószínűségű.
Esetleg még a doktori vézettség esetén látványos, hogy a tapasztalat növekedésével csökken folyamatosan a Not satisfied kategória valószínűsége.

Félév végi ujjgyakorlatként elkészíthetjük még a fenti ábrát egy lépésben OpSys szerinti megbontásban is!
Ne feledjük, hogy ha az OpSys változót is használni akarjuk, akkor a temp_tabla data frame-t is újra kell generálni!

# a tábla szűkítést és a melt függvény alkalmazást lehet egyben is az alább látható módon
temp_tabla <- melt(ProgDataWithProbs[,c("YearsCodePro", "OpSys", "Not satisfied",
                                           "Slightly satisfied", "Very satisfied")],
                      id.vars = c("YearsCodePro", "OpSys"), value.name = "Probability")
head(temp_tabla)
##   YearsCodePro       OpSys      variable Probability
## 1            8     Windows Not satisfied   0.2220232
## 2            6     Windows Not satisfied   0.3616731
## 3            5     Windows Not satisfied   0.4043842
## 4            3 Linux-based Not satisfied   0.3677947
## 5            2     Windows Not satisfied   0.4175159
## 6            1 Linux-based Not satisfied   0.5555360
ggplot(temp_tabla, aes(x = YearsCodePro, y = Probability, colour = variable)) + geom_line(size = 1) +
  facet_grid(OpSys ~ ., scales = "free")

Úgy látszik a MacOS-en dolgozók eléggé sokáig Not satisfied-ok legnagyobb valószínűséggel. De \(20\) év munkatapasztalat után már legnagyobb valószínűségel Very Satisfied-ok…a Mac-hez úgy látszik idő kell. :)
Ellenben a Linuxosok már \(10\) év munkatapasztalat után eléggé Very Satisfied-nak festenek legnagyobb valószínűséggel.