# Generalized Linear Models
#---Multinomial Logistic Regression----
# This code is from: https://bookdown.org/sarahwerth2024/CategoricalBook/multinomial-logit-regression-r.html#lab-overview-13
df_MAR <- read.csv("D:/SEMESTER 4/Analisis Multivariat/hsbdemo.csv", header = TRUE, sep = ",")
#df_MAR
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
df_MAR <- df_MAR %>%
mutate(
prog = factor(prog,
levels = c("general", "academic", "vocation")),
ses = factor(ses,
levels = c("low", "middle", "high")),
female = factor(female,
levels = c("female", "male")),
schtyp = factor(schtyp,
levels = c("public", "private")),
)
str(df_MAR)
## 'data.frame': 200 obs. of 13 variables:
## $ id : int 45 108 15 67 153 51 164 133 2 53 ...
## $ female : Factor w/ 2 levels "female","male": 1 2 2 2 2 1 2 2 1 2 ...
## $ ses : Factor w/ 3 levels "low","middle",..: 1 2 3 1 2 3 2 2 2 2 ...
## $ schtyp : Factor w/ 2 levels "public","private": 1 1 1 1 1 1 1 1 1 1 ...
## $ prog : Factor w/ 3 levels "general","academic",..: 3 1 3 3 3 1 3 3 3 3 ...
## $ read : int 34 34 39 37 39 42 31 50 39 34 ...
## $ write : int 35 33 39 37 31 36 36 31 41 37 ...
## $ math : int 41 41 44 42 40 42 46 40 33 46 ...
## $ science: int 29 36 26 33 39 31 39 34 42 39 ...
## $ socst : int 26 36 42 32 51 39 46 31 41 31 ...
## $ honors : chr "not enrolled" "not enrolled" "not enrolled" "not enrolled" ...
## $ awards : int 0 0 0 0 0 0 0 0 0 0 ...
## $ cid : int 1 1 1 1 1 1 1 1 1 1 ...
#Multinomial with 1 IV
library("nnet")
fit_basic <- multinom(prog ~ ses, data = df_MAR)
## # weights: 12 (6 variable)
## initial value 219.722458
## iter 10 value 195.705517
## final value 195.705188
## converged
summary(fit_basic)
## Call:
## multinom(formula = prog ~ ses, data = df_MAR)
##
## Coefficients:
## (Intercept) sesmiddle seshigh
## academic 0.1718494 0.6165982 1.36861772
## vocation -0.2876865 0.7259513 0.03638183
##
## Std. Errors:
## (Intercept) sesmiddle seshigh
## academic 0.3393103 0.4334269 0.5000539
## vocation 0.3818815 0.4775889 0.6323011
##
## Residual Deviance: 391.4104
## AIC: 403.4104
#or
library("broom")
library("kableExtra")
## Warning: package 'kableExtra' was built under R version 4.4.3
##
## Attaching package: 'kableExtra'
##
## The following object is masked from 'package:dplyr':
##
## group_rows
tidy(fit_basic, conf.int = TRUE) %>%
kable() %>%
kable_styling("basic", full_width = FALSE)
|
y.level
|
term
|
estimate
|
std.error
|
statistic
|
p.value
|
conf.low
|
conf.high
|
|
academic
|
(Intercept)
|
0.1718494
|
0.3393103
|
0.5064669
|
0.6125289
|
-0.4931865
|
0.8368853
|
|
academic
|
sesmiddle
|
0.6165982
|
0.4334269
|
1.4226119
|
0.1548487
|
-0.2329028
|
1.4660993
|
|
academic
|
seshigh
|
1.3686177
|
0.5000539
|
2.7369405
|
0.0062014
|
0.3885301
|
2.3487053
|
|
vocation
|
(Intercept)
|
-0.2876865
|
0.3818815
|
-0.7533398
|
0.4512457
|
-1.0361606
|
0.4607875
|
|
vocation
|
sesmiddle
|
0.7259513
|
0.4775889
|
1.5200341
|
0.1285024
|
-0.2101056
|
1.6620083
|
|
vocation
|
seshigh
|
0.0363818
|
0.6323011
|
0.0575388
|
0.9541160
|
-1.2029056
|
1.2756693
|
### INTERPRETASI LOG ODDS SES ONLY###
#>pada program academic, anak dari SES "middle" memiliki log-odds 0.6166, yang artinya memiliki probabilitas lebih besar daripada anak dari SES "low" (0.1718). anak dari SES "high" memiliki log-odds 1.3686, yang artinya memiliki probabilitas lebih besar daripada anak SES "low" (0,1718).
#>pada program vocation, anak dari SES "middle" memiliki log-odds 0,729, yang artinya memiliki probabilitas lebih besar daripada anak dari SES "low" (-0,2877). anak dari SES "high" memiliki log-odds 0.0364, yang artinya memiliki probabilitas lebih besar daripada anak dari SES "low" (-0,2877)
#>log-odds anak dari SES "low" -0,2877 pada vocation, berarti probabilitas mereka memilih vocation lebih rendah daripada general. sebaliknya, log-odds anak dari SES "low" pada academic 0,1718 pada academic, berarti probabilitas mereka memilih academic lebih besar daripada general.
##---Multinomial with Multiple IV----
fit_full <- multinom(prog ~ ses + female + schtyp, data = df_MAR)
## # weights: 18 (10 variable)
## initial value 219.722458
## iter 10 value 190.175247
## final value 190.148334
## converged
##---1. Interpret for the relative log odds----
summary(fit_full)
## Call:
## multinom(formula = prog ~ ses + female + schtyp, data = df_MAR)
##
## Coefficients:
## (Intercept) sesmiddle seshigh femalemale schtypprivate
## academic 0.2105078 0.5506420 1.3414584 -0.2065183 0.5233744
## vocation -0.2186660 0.8887273 0.1479293 -0.1125838 -1.5231430
##
## Std. Errors:
## (Intercept) sesmiddle seshigh femalemale schtypprivate
## academic 0.3619400 0.4512813 0.5112042 0.3734168 0.5157763
## vocation 0.4089818 0.4922409 0.6421229 0.4249998 0.8560186
##
## Residual Deviance: 380.2967
## AIC: 400.2967
#or
tidy(fit_full, conf.int = TRUE) %>%
kable() %>%
kable_styling("basic", full_width = FALSE)
|
y.level
|
term
|
estimate
|
std.error
|
statistic
|
p.value
|
conf.low
|
conf.high
|
|
academic
|
(Intercept)
|
0.2105078
|
0.3619400
|
0.5816098
|
0.5608295
|
-0.4988815
|
0.9198972
|
|
academic
|
sesmiddle
|
0.5506420
|
0.4512813
|
1.2201746
|
0.2223987
|
-0.3338531
|
1.4351372
|
|
academic
|
seshigh
|
1.3414584
|
0.5112042
|
2.6241145
|
0.0086875
|
0.3395166
|
2.3434002
|
|
academic
|
femalemale
|
-0.2065183
|
0.3734168
|
-0.5530503
|
0.5802290
|
-0.9384019
|
0.5253653
|
|
academic
|
schtypprivate
|
0.5233744
|
0.5157763
|
1.0147314
|
0.3102339
|
-0.4875285
|
1.5342774
|
|
vocation
|
(Intercept)
|
-0.2186660
|
0.4089818
|
-0.5346595
|
0.5928853
|
-1.0202557
|
0.5829236
|
|
vocation
|
sesmiddle
|
0.8887273
|
0.4922409
|
1.8054723
|
0.0710008
|
-0.0760471
|
1.8535016
|
|
vocation
|
seshigh
|
0.1479293
|
0.6421229
|
0.2303754
|
0.8178001
|
-1.1106085
|
1.4064671
|
|
vocation
|
femalemale
|
-0.1125838
|
0.4249998
|
-0.2649032
|
0.7910841
|
-0.9455682
|
0.7204006
|
|
vocation
|
schtypprivate
|
-1.5231430
|
0.8560186
|
-1.7793340
|
0.0751850
|
-3.2009086
|
0.1546226
|
### MODEL KATEGORI AKADEMIK DIKOMPARE DENGAN KATEGORI GENERAL ###
#0.210 + 0.550X1 + 1.341X2 - 0.206X3 + 0.523X4
### MODEL KATEGORI VOKASI DIKOMPARE DENGAN KATEGORI GENERAL ###
#-0.219 + 0.889X1 + 0.148X2 - 0.113X3 - 1.523X4
### INTERPRETASI LOG ODDS ###
#X= status sosial ekonomi (ses)
#Kategori = academic
#>Kelompok yang status sosial ekonominya menengah (dibandingkan dengan kelompok yang status sosial ekonominya rendah) dikaitkan dengan peningkatan 0,550 dalam log-odds kelompok yang menginginkan jenis program pendidikan akademis dibandingkan dengan yang memiliki jenis program pendidikan umum, secara statistik tidak signifikan.
#>Kelompok yang status sosial ekonominya tinggi (dibandingkan dengan kelompok yang status sosial ekonominya rendah) dikaitkan dengan peningkatan 1,341 dalam log-odds kelompok yang menginginkan jenis program pendidikan akademis dibandingkan dengan yang memiliki jenis program pendidikan umum, secara statistik signifikan.
#Kategori = vocation
#>Kelompok yang status sosial ekonominya menengah (dibandingkan dengan kelompok yang status sosial ekonominya rendah) dikaitkan dengan peningkatan 0,889 dalam log-odds kelompok yang menginginkan jenis program pendidikan kejuruan dibandingkan dengan yang memiliki jenis program pendidikan umum, secara statistik tidak signifikan.
#>Kelompok yang berstatus sosial ekonomi tinggi (dibandingkan kelompok yang berstatus sosial ekonomi rendah) dikaitkan dengan peningkatan 1,148 dalam log-odds kelompok yang menginginkan jenis program pendidikan kejuruan dibandingkan dengan kelompok yang memiliki jenis program pendidikan umum, hal ini tidak signifikan secara statistik.
#X= jenis kelamin (perempuan)
#Kategori= academic
#>Kelompok yang berjenis kelamin laki-laki (bandingkan dengan kelompok yang berjenis kelamin perempuan) dikaitkan dengan penurunan 0,207 dalam log-odds kelompok yang menginginkan jenis program pendidikan akademis dibandingkan dengan yang memiliki jenis program pendidikan umum, hal ini tidak signifikan secara statistik.
#Kategori= vocation
#>Kelompok yang berjenis kelamin laki-laki (bandingkan dengan kelompok yang berjenis kelamin perempuan) dikaitkan dengan penurunan 0,113 dalam log-odds kelompok yang menginginkan jenis program pendidikan kejuruan dibandingkan dengan yang memiliki jenis program pendidikan umum, hal ini tidak signifikan secara statistik.
#X= jenis sekolah (schtyp)
#Kategori= academic
#>Kelompok yang jenis sekolahnya swasta (bandingkan dengan kelompok yang jenis sekolahnya negeri) dikaitkan dengan peningkatan 0,523 dalam log-odds kelompok yang menginginkan jenis program pendidikan akademis dibandingkan dengan yang memiliki jenis program pendidikan umum, hal ini tidak signifikan secara statistik.
#Kategori= vocation
#>Kelompok yang jenis sekolahnya swasta (bandingkan dengan kelompok yang jenis sekolahnya negeri) dikaitkan dengan penurunan 1,523 dalam log-odds kelompok yang menginginkan jenis program pendidikan kejuruan dibandingkan dengan yang memiliki jenis program pendidikan umum, hal ini tidak signifikan secara statistik.
##---2. Interpret for the relative risk ratios----
exp(coef(fit_full))
## (Intercept) sesmiddle seshigh femalemale schtypprivate
## academic 1.2343047 1.734366 3.824617 0.8134114 1.6877131
## vocation 0.8035901 2.432032 1.159431 0.8935225 0.2180256
#or
tidy(fit_full, conf.int = TRUE, exponentiate = TRUE) %>%
kable() %>%
kable_styling("basic", full_width = FALSE)
|
y.level
|
term
|
estimate
|
std.error
|
statistic
|
p.value
|
conf.low
|
conf.high
|
|
academic
|
(Intercept)
|
1.2343047
|
0.3619400
|
0.5816098
|
0.5608295
|
0.6072094
|
2.509032
|
|
academic
|
sesmiddle
|
1.7343661
|
0.4512813
|
1.2201746
|
0.2223987
|
0.7161590
|
4.200221
|
|
academic
|
seshigh
|
3.8246172
|
0.5112042
|
2.6241145
|
0.0086875
|
1.4042685
|
10.416595
|
|
academic
|
femalemale
|
0.8134114
|
0.3734168
|
-0.5530503
|
0.5802290
|
0.3912526
|
1.691076
|
|
academic
|
schtypprivate
|
1.6877131
|
0.5157763
|
1.0147314
|
0.3102339
|
0.6141423
|
4.637973
|
|
vocation
|
(Intercept)
|
0.8035901
|
0.4089818
|
-0.5346595
|
0.5928853
|
0.3605028
|
1.791268
|
|
vocation
|
sesmiddle
|
2.4320324
|
0.4922409
|
1.8054723
|
0.0710008
|
0.9267725
|
6.382128
|
|
vocation
|
seshigh
|
1.1594309
|
0.6421229
|
0.2303754
|
0.8178001
|
0.3293585
|
4.081510
|
|
vocation
|
femalemale
|
0.8935225
|
0.4249998
|
-0.2649032
|
0.7910841
|
0.3884588
|
2.055256
|
|
vocation
|
schtypprivate
|
0.2180256
|
0.8560186
|
-1.7793340
|
0.0751850
|
0.0407252
|
1.167217
|
### INTERPRETASI RELATIVE RISK RASIO ###
#X= status sosial ekonomi (ses)
#Kategori = academic
#>Kelompok yang persentase status sosial ekonomi anggotanya sedang (dibandingkan dengan kelompok yang status sosial ekonominya rendah) melipatgandakan peluang untuk menginginkan jenis program pendidikan akademis dibandingkan dengan memiliki jenis program pendidikan umum sebesar 1,73 (73%), secara statistik tidak signifikan.
#>Kelompok yang persentase status sosial ekonomi anggotanya tinggi (dibandingkan dengan kelompok yang status sosial ekonominya rendah) melipatgandakan peluang untuk menginginkan jenis program pendidikan akademis dibandingkan dengan memiliki jenis program pendidikan umum sebesar 3,82 (282%), secara statistik signifikan.
#Kategori = vocation
#>suatu kelompok yang persentase anggotanya berstatus sosial ekonomi sedang (dibandingkan dengan kelompok yang status sosial ekonominya rendah) melipatgandakan peluang untuk menginginkan jenis program pendidikan kejuruan dibandingkan dengan yang memiliki jenis program pendidikan umum sebesar 2,43 (143%), secara statistik tidak signifikan.
#>suatu kelompok yang persentase anggotanya berstatus sosial ekonomi tinggi (dibandingkan dengan kelompok yang status sosial ekonominya rendah) melipatgandakan peluang untuk menginginkan jenis program pendidikan kejuruan dibandingkan dengan yang memiliki jenis program pendidikan umum sebesar 1,16 (16%), tetapi tidak signifikan secara statistik.
#X= jenis kelamin (female)
#Kategori= academic
#>Kelompok yang persentase jenis kelamin anggotanya adalah laki-laki (dibandingkan dengan kelompok yang jenis kelaminnya adalah perempuan) peluangnya berkelipatan pada menginginkan jenis program pendidikan akademis dibandingkan dengan memiliki jenis program pendidikan umum sebesar 0,81 (-18%), secara statistik tidak signifikan.
#Kategori= vocation
#>Kelompok yang persentase jenis kelamin anggotanya adalah laki-laki (dibandingkan dengan kelompok yang jenis kelaminnya adalah perempuan) peluangnya berkelipatan pada menginginkan jenis program pendidikan kejuruan dibandingkan dengan memiliki jenis program pendidikan umum sebesar 0,89 (-11%), secara statistik tidak signifikan.
#X= jenis sekolah (schtyp)
#Kategori= academic
#>persentase anggota kelompok yang jenis sekolahnya swasta (dibandingkan dengan kelompok yang jenis sekolahnya negeri) memiliki kelipatan peluang menginginkan jenis program pendidikan akademis dibandingkan dengan memiliki jenis program pendidikan umum sebesar 1,69 (69%), secara statistik tidak signifikan.
#Kategori= vocation
#>persentase anggota kelompok yang jenis sekolahnya swasta (dibandingkan dengan kelompok yang jenis sekolahnya negeri) memiliki kelipatan peluang menginginkan jenis program pendidikan kejuruan dibandingkan dengan memiliki jenis program pendidikan umum sebesar 0,21 (-0,79%), secara statistik tidak signifikan.
##---3. Interpret for the marginal effect----
library("marginaleffects")
## Warning: package 'marginaleffects' was built under R version 4.4.3
mfx_schtyp <- avg_comparisons(fit_full, variables = "schtyp", type = "probs")
mfx_schtyp
##
## Group Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## general -0.0149 0.0824 -0.181 0.85674 0.2 -0.1763 0.147
## academic 0.2417 0.0893 2.706 0.00681 7.2 0.0666 0.417
## vocation -0.2268 0.0538 -4.218 < 0.001 15.3 -0.3322 -0.121
##
## Term: schtyp
## Type: probs
## Comparison: private - public
mfx_ses <- avg_comparisons(fit_full, variables = "ses", type = "probs")
mfx_ses
##
## Group Contrast Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## general high - low -0.1882 0.0864 -2.177 0.02949 5.1 -0.3576 -0.0187
## general middle - low -0.1342 0.0837 -1.603 0.10889 3.2 -0.2983 0.0299
## academic high - low 0.2949 0.0942 3.130 0.00175 9.2 0.1102 0.4796
## academic middle - low 0.0267 0.0894 0.299 0.76511 0.4 -0.1485 0.2019
## vocation high - low -0.1067 0.0731 -1.460 0.14418 2.8 -0.2499 0.0365
## vocation middle - low 0.1075 0.0760 1.415 0.15720 2.7 -0.0415 0.2565
##
## Term: ses
## Type: probs
### INTERPRETASI MARGINAL EFFECT ###
#schtyp
#>Sekolah "private" meningkatkan kemungkinan anak memilih program "academic" sebesar 24.17% dan menurunkan "vocation" sebesar 22.68%, tidak signifikan terhadap program general.
#SES
#>pada program general, anak dari SES "high" memiliki probabilitas lebih rendah (18,82%) secara signifikan untuk memilih program "general" dibanding anak dari SES "low"
#>pada program academic, anak dari SES "high" memiliki probabilitas lebih tinggi (29,49%) untuk memilih program "academic" dibanding anak dari SES "low"
#>tidak ada perbedaan signifikan antara anak dari SES "middle" dan "low" dalam program apapun. begitupun dengan anak dari SES "high" dan "low" dalam program vocation.