Analisis Regresi Logistik

1. Multinomial Logistic Regression

Penjelasan Analisis Regresi Logistik Multinomial

Tujuan Utama Analisis

Analisis regresi logistik mutinomial pada dataset yang kami gunakan bertujuan untuk menganalisis faktor-faktor yang mempengaruhi pilihan program pendidikan siswa, dengan beberapa fokus pada 2 variabel berikut:

1. Variabel Prediktor:
   • Status Sosial Ekonomi (SES)
   • Skor Menulis (Write)
2. Variabel Dependen:
   • Program Pendidikan (3 kategori):
     • Program Akademik
     • Program Umum
     • Program Vokasi

library(foreign)
library(nnet)
library(ggplot2)
library(reshape2)

data <- "/Users/nabilanasywazen/Downloads/hsbdemo.dta"
data <- read.dta(data)
data

with(data, table(ses, prog))

        prog
ses      general academic vocation
  low         16       19       12
  middle      20       44       31
  high         9       42        7

Menunjukkan distribusi program berdasarkan status sosial ekonomi.

with(data, do.call(rbind, tapply(write, prog, function(x) c(Mean = mean(x), Std = sd(x)))))

             Mean      Std
general  51.33333 9.397775
academic 56.25714 7.943343
vocation 46.76000 9.318754

Menghitung mean dan standar deviasi skor menulis berdasarkan kategori program.

Siswa dalam program “academic” cenderung memiliki skor tertinggi dengan variasi nilai yang lebih konsisten, siswa dalam program general memiliki nilai di tengah-tengah, sementara siswa dalam program “vocation” memiliki skor terrendah.

data$prog2 <- relevel(data$prog, ref = "academic")
model_test <- multinom(prog2 ~ ses + write, data = data)

# weights:  15 (8 variable)
initial  value 219.722458 
iter  10 value 179.982880
final  value 179.981726 
converged

model_summary <- summary(model_test)
print(model_summary)

Call:
multinom(formula = prog2 ~ ses + write, data = data)

Coefficients:
         (Intercept)  sesmiddle    seshigh      write
general     2.852198 -0.5332810 -1.1628226 -0.0579287
vocation    5.218260  0.2913859 -0.9826649 -0.1136037

Std. Errors:
         (Intercept) sesmiddle   seshigh      write
general     1.166441 0.4437323 0.5142196 0.02141097
vocation    1.163552 0.4763739 0.5955665 0.02221996

Residual Deviance: 359.9635 
AIC: 375.9635 

z_values <- model_summary$coefficients / model_summary$standard.errors
p_values <- 2 * (1 - pnorm(abs(z_values)))

z_values

         (Intercept)  sesmiddle   seshigh     write
general     2.445214 -1.2018081 -2.261334 -2.705562
vocation    4.484769  0.6116747 -1.649967 -5.112689

p_values

          (Intercept) sesmiddle    seshigh        write
general  0.0144766100 0.2294379 0.02373856 6.818902e-03
vocation 0.0000072993 0.5407530 0.09894976 3.176045e-07

• Nilai write signifikan untuk kedua kategori (general dan vocational vs academic)
• SES menengah dan tinggi memiliki perbedaan signifikan terhadap SES rendah dalam mempengaruhi pemilihan program

odds_ratios <- exp(coef(model_test))
odds_ratios

         (Intercept) sesmiddle   seshigh     write
general     17.32582 0.5866769 0.3126026 0.9437172
vocation   184.61262 1.3382809 0.3743123 0.8926116

• Odds ratio < 1: variabel mengurangi peluang kategori dibandingkan referensi
• Odds ratio > 1: variabel meningkatkan peluang kategori dibandingkan referensi

predicted_probabilities <- predict(model_test, type = "probs")
head(predicted_probabilities)

   academic   general  vocation
1 0.1482764 0.3382454 0.5134781
2 0.1202017 0.1806283 0.6991700
3 0.4186747 0.2368082 0.3445171
4 0.1726885 0.3508384 0.4764731
5 0.1001231 0.1689374 0.7309395
6 0.3533566 0.2377976 0.4088458

eg_ses <- data.frame(
  ses = c("low", "middle", "high"),
  write = mean(data$write)
)

predicted_rpogram <- predict(model_test, newdata = eg_ses, type = "probs")

predicted_rpogram

   academic   general  vocation
1 0.4396845 0.3581917 0.2021238
2 0.4777488 0.2283353 0.2939159
3 0.7009007 0.1784939 0.1206054

data_write <- data.frame(
  ses = rep(c("low", "middle", "high"), each = 41),
  write = rep(30:70, 3)
)

predicted_output <- cbind(data_write,
  predict(model_test, newdata = data_write, type = "probs", se = TRUE)
)

mean_probability <- by(predicted_output[, 3:5], predicted_output$ses, colMeans)

mean_probability

predicted_output$ses: high
 academic   general  vocation 
0.6164315 0.1808037 0.2027648 
------------------------------------------------------------------------ 
predicted_output$ses: low
 academic   general  vocation 
0.3972977 0.3278174 0.2748849 
------------------------------------------------------------------------ 
predicted_output$ses: middle
 academic   general  vocation 
0.4256198 0.2010864 0.3732938 

data_long <- melt(
  pp.write, 
  id.vars = c("ses", "write"), 
  value.name = "probability"
)

head(data_long)

ggplot(data_long, aes(x = write, y = probability, colour = ses)) + 
  geom_line() + 
  facet_grid(variable ~ ., scales = "free") + 
  scale_colour_manual(values = c("low" = "blue", "middle" = "green", "high" = "red")) +
  labs(
    title = "Probabilitas Program Berdasarkan Skor Menulis dan SES",
    x = "Skor Menulis",
    y = "Probabilitas",
    colour = "Status Sosial Ekonomi"
  )

Kesimpulan Analisis:

Risiko Relative Risk Ratio

Rasio Risiko Relatif untuk Skor Menulis (write):

Rasio risiko relatif untuk peningkatan satu unit dalam variabel menulis (write) adalah 0.9437 untuk berada di program umum (general) dibandingkan dengan program akademik (academic). Ini berarti bahwa setiap peningkatan satu unit pada skor menulis menurunkan peluang untuk memilih program akademik dan meningkatkan peluang untuk memilih program umum.

Rasio Risiko Relatif untuk Status Sosial Ekonomi (ses):

Rasio risiko relatif untuk peralihan dari ses = 1 (rendah) ke ses = 3 (tinggi) adalah 0.3126 untuk berada di program umum dibandingkan dengan program akademik. Ini berarti bahwa individu dengan ses tinggi (ses = 3) memiliki kemungkinan lebih kecil untuk memilih program umum dibandingkan dengan program akademik, dibandingkan dengan individu yang memiliki ses rendah (ses = 1).

Hasil Prediksi Probabilitas:

1. Pengaruh Status Sosial Ekonomi (SES) terhadap Pemilihan Program

Individu yang memiliki status sosial ekonomi tinggi (SES = tinggi) cenderung memilih program akademik dengan persentase mencapai 61. 64%. Hal ini berbeda dengan individu yang memiliki status sosial ekonomi rendah (SES = rendah), yang menampilkan kecenderungan lebih merata dalam memilih antara program general (32. 78%) dan vokasi (27. 49%), serta memiliki kemungkinan yang lebih rendah untuk memilih program akademik (39. 73%). Sementara itu, individu dengan status sosial ekonomi menengah (SES = menengah) menunjukkan peluang yang hampir seimbang antara memilih program akademik (42. 56%) dan vokasi (37. 33%), meskipun mereka masih lebih memilih program akademik dibandingkan dengan program general (20. 11%).

2. Perbandingan antar Program

Program Akademik: Program ini menjadi pilihan utama bagi individu dengan status sosial ekonomi tinggi, dengan probabilitas pemilihan yang melebihi 60%. Sebaliknya, individu dengan status sosial ekonomi rendah atau menengah menunjukkan kecenderungan yang lebih rendah untuk memilih program ini.

Program General: Program ini paling banyak diminati oleh individu dengan status sosial ekonomi rendah, di mana 32. 78% dari mereka memilihnya. Namun, probabilitas pemilihan program ini menurun pada individu dari status sosial ekonomi tinggi, hanya mencapai 18. 08%.

Program Vokasi: Program ini menunjukkan popularitas yang relatif tinggi di kalangan individu dengan status sosial ekonomi rendah dan menengah, dengan angka tertinggi terdapat pada individu dengan SES rendah (27. 49%) dan SES menengah (37. 33%).

---

2. Ordinal Logistic Regression

Penjelasan Analisis Ordinal Logistik Multinomial

Model ini adalah ordinal logistic regression untuk memodelkan peluang (Probability) seseorang memilih untuk melanjutkan pendidikan (apply: unlikely, somewhat likely, very likely) berdasarkan:

• pared: apakah orang tua pernah kuliah (1 = ya, 0 = tidak),
• public: apakah siswa berasal dari sekolah negeri (1 = ya, 0 = tidak),
• gpa: nilai rata-rata siswa.

library(foreign)
library(ggplot2)
library(MASS)
library(Hmisc)
library(reshape2)

# Baca dataset
dat <- read.dta("/Users/nabilanasywazen/Downloads/ologit.dta")

# Lihat data awal
head(dat)

# Cek tabel untuk setiap variabel kategorik
lapply(dat[, c("apply", "pared", "public")], table)

$apply

       unlikely somewhat likely     very likely 
            220             140              40 

$pared

  0   1 
337  63 

$public

  0   1 
343  57 

# Tabel silang tiga variabel
ftable(xtabs(~ public + apply + pared, data = dat))

                       pared   0   1
public apply                        
0      unlikely              175  14
       somewhat likely        98  26
       very likely            20  10
1      unlikely               25   6
       somewhat likely        12   4
       very likely             7   3

# Rangkuman statistik GPA
summary(dat$gpa)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.900   2.720   2.990   2.999   3.270   4.000 

# Standar deviasi GPA
sd(dat$gpa)

[1] 0.3979409

# Boxplot GPA berdasarkan apply, pared, dan public
ggplot(dat, aes(x = apply, y = gpa)) +
  geom_boxplot() +
  geom_jitter(alpha = 0.5) +
  facet_grid(pared ~ public, margins = TRUE) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Fit model ordered logit
m <- polr(apply ~ pared + public + gpa, data = dat, Hess = TRUE)

# Lihat hasil ringkasan model
summary(m)

Call:
polr(formula = apply ~ pared + public + gpa, data = dat, Hess = TRUE)

Coefficients:
          Value Std. Error t value
pared   1.04769     0.2658  3.9418
public -0.05879     0.2979 -0.1974
gpa     0.61594     0.2606  2.3632

Intercepts:
                            Value   Std. Error t value
unlikely|somewhat likely     2.2039  0.7795     2.8272
somewhat likely|very likely  4.2994  0.8043     5.3453

Residual Deviance: 717.0249 
AIC: 727.0249 

# Ambil tabel koefisien
(ctable <- coef(summary(m)))

                                  Value Std. Error    t value
pared                        1.04769010  0.2657894  3.9418050
public                      -0.05878572  0.2978614 -0.1973593
gpa                          0.61594057  0.2606340  2.3632399
unlikely|somewhat likely     2.20391473  0.7795455  2.8271792
somewhat likely|very likely  4.29936315  0.8043267  5.3452947

# Hitung p-value
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2

# Gabungkan tabel koefisien dan p-value
(ctable <- cbind(ctable, "p value" = p))

                                  Value Std. Error    t value      p value
pared                        1.04769010  0.2657894  3.9418050 8.087072e-05
public                      -0.05878572  0.2978614 -0.1973593 8.435464e-01
gpa                          0.61594057  0.2606340  2.3632399 1.811594e-02
unlikely|somewhat likely     2.20391473  0.7795455  2.8271792 4.696004e-03
somewhat likely|very likely  4.29936315  0.8043267  5.3452947 9.027008e-08

# Confidence interval (CI)
confint(m)  # metode profile likelihood

Waiting for profiling to be done...

            2.5 %    97.5 %
pared   0.5281768 1.5721750
public -0.6522060 0.5191384
gpa     0.1076202 1.1309148

# Confidence interval (asumsi normalitas)
confint.default(m)

            2.5 %    97.5 %
pared   0.5267524 1.5686278
public -0.6425833 0.5250119
gpa     0.1051074 1.1267737

# Hitung odds ratio
exp(coef(m))

    pared    public       gpa 
2.8510579 0.9429088 1.8513972 

# Odds Ratio dan CI
exp(cbind(OR = coef(m), confint(m)))

Waiting for profiling to be done...

              OR     2.5 %   97.5 %
pared  2.8510579 1.6958376 4.817114
public 0.9429088 0.5208954 1.680579
gpa    1.8513972 1.1136247 3.098490

# Fungsi membuat ringkasan logits
sf <- function(y) {
  c('Y>=1' = qlogis(mean(y >= 1)),
    'Y>=2' = qlogis(mean(y >= 2)),
    'Y>=3' = qlogis(mean(y >= 3)))
}

# Terapkan fungsi ke data
(s <- with(dat, summary(as.numeric(apply) ~ pared + public + gpa, fun = sf)))

as.numeric(apply)      N= 400  

+-------+-----------+---+----+-----------+---------+
|       |           |  N|Y>=1|       Y>=2|     Y>=3|
+-------+-----------+---+----+-----------+---------+
|  pared|         No|337| Inf|-0.37833644|-2.440735|
|       |        Yes| 63| Inf| 0.76546784|-1.347074|
+-------+-----------+---+----+-----------+---------+
| public|         No|343| Inf|-0.20479441|-2.345006|
|       |        Yes| 57| Inf|-0.17589067|-1.547563|
+-------+-----------+---+----+-----------+---------+
|    gpa|[1.90,2.73)|102| Inf|-0.39730180|-2.772589|
|       |[2.73,3.00)| 99| Inf|-0.26415158|-2.302585|
|       |[3.00,3.28)|100| Inf|-0.20067070|-2.090741|
|       |[3.28,4.00]| 99| Inf| 0.06062462|-1.803594|
+-------+-----------+---+----+-----------+---------+
|Overall|           |400| Inf|-0.20067070|-2.197225|
+-------+-----------+---+----+-----------+---------+

# Model logit untuk Y>=2
glm(I(as.numeric(apply) >= 2) ~ pared, family = "binomial", data = dat)

Call:  glm(formula = I(as.numeric(apply) >= 2) ~ pared, family = "binomial", 
    data = dat)

Coefficients:
(Intercept)        pared  
    -0.3783       1.1438  

Degrees of Freedom: 399 Total (i.e. Null);  398 Residual
Null Deviance:      550.5 
Residual Deviance: 534.1    AIC: 538.1

# Model logit untuk Y>=3
glm(I(as.numeric(apply) >= 3) ~ pared, family = "binomial", data = dat)

Call:  glm(formula = I(as.numeric(apply) >= 3) ~ pared, family = "binomial", 
    data = dat)

Coefficients:
(Intercept)        pared  
     -2.441        1.094  

Degrees of Freedom: 399 Total (i.e. Null);  398 Residual
Null Deviance:      260.1 
Residual Deviance: 252.2    AIC: 256.2

# Koreksi hasil logits
s[, 4] <- s[, 4] - s[, 3]
s[, 3] <- s[, 3] - s[, 3]

# Print hasilnya
s

as.numeric(apply)      N= 400  

+-------+-----------+---+----+----+---------+
|       |           |  N|Y>=1|Y>=2|     Y>=3|
+-------+-----------+---+----+----+---------+
|  pared|         No|337| Inf|   0|-2.062399|
|       |        Yes| 63| Inf|   0|-2.112541|
+-------+-----------+---+----+----+---------+
| public|         No|343| Inf|   0|-2.140211|
|       |        Yes| 57| Inf|   0|-1.371672|
+-------+-----------+---+----+----+---------+
|    gpa|[1.90,2.73)|102| Inf|   0|-2.375287|
|       |[2.73,3.00)| 99| Inf|   0|-2.038434|
|       |[3.00,3.28)|100| Inf|   0|-1.890070|
|       |[3.28,4.00]| 99| Inf|   0|-1.864219|
+-------+-----------+---+----+----+---------+
|Overall|           |400| Inf|   0|-1.996554|
+-------+-----------+---+----+----+---------+

# Plot hasil logits
plot(s, which = 1:3, pch = 1:3, xlab = 'logit', main = ' ', xlim = range(s[,3:4]))

# Buat data baru untuk prediksi
newdat <- data.frame(
  pared = rep(0:1, 200),
  public = rep(0:1, each = 200),
  gpa = rep(seq(1.9, 4, length.out = 100), 4)
)

# Prediksi probabilitas
newdat <- cbind(newdat, predict(m, newdat, type = "probs"))

# Lihat hasil prediksi
head(newdat)

NA

# Ubah data prediksi menjadi long format
lnewdat <- melt(newdat, id.vars = c("pared", "public", "gpa"),
                variable.name = "Level", value.name = "Probability")

# Lihat hasil
head(lnewdat)

# Plot probabilitas prediksi
ggplot(lnewdat, aes(x = gpa, y = Probability, colour = Level)) +
  geom_line() +
  facet_grid(pared ~ public, labeller = "label_both")

Interpretasi Hasil Model

Berikut penjelasan per variabelnya:

1. Pendidikan Orang Tua Koefisien: 1.048 (signifikan, t = 3.94) OR (Odds Ratio): 2.85 Kesimpulan: Jika orang tua pernah kuliah, peluang siswa untuk lebih “likely” melanjutkan pendidikan meningkat 2.85 kali lipat dibandingkan yang orang tuanya tidak kuliah.
2. Jenis Sekolah: Koefisien: -0.059 (tidak signifikan, t = -0.20) OR: 0.94 Kesimpulan: Koefisien negatif kecil sehingga asal dari sekolah negeri tidak berpengaruh signifikan terhadap keputusan apply.
3. GPA: Koefisien: 0.616 (signifikan, t = 2.36) OR: 1.85 Kesimpulan: Setiap kenaikan 1 poin di IPK siswa, peluang untuk melamar akan meningkat 85% lebih besar.

Interpretasi Grafik

Grafik menunjukkan hubungan antara GPA dan probabilitas level apply pada kombinasi kondisi public dan pared.

• Semakin tinggi GPA, probabilitas untuk kategori “unlikely” menurun, sedangkan untuk “somewhat likely” dan “very likely” meningkat.
• Siswa dengan orang tua yang kuliah (pared = 1) memiliki probabilitas lebih tinggi untuk “somewhat likely” dan “very likely”, dibandingkan siswa dengan orang tua tidak kuliah (pared = 0).
• Status sekolah (public) tidak mengubah pola secara drastis, sesuai dengan hasil model bahwa variabel ini tidak signifikan.

Kesimpulan:

1. Pendidikan orang tua (pared) dan nilai akademik (gpa) merupakan faktor penting yang mempengaruhi keputusan siswa untuk melanjutkan pendidikan.
2. Status sekolah (public) tidak berpengaruh signifikan terhadap keputusan tersebut. Kebijakan intervensi untuk meningkatkan angka 3. lanjut studi sebaiknya fokus pada: Meningkatkan prestasi akademik (GPA), Memberikan dukungan ekstra kepada siswa yang orang tuanya tidak berpendidikan tinggi. ```