Pendekatan Regresi Non-Linear

Dalam analisis data, tidak semua variabel dependen (\(Y\)) berbentuk kontinu atau mengikuti distribusi normal. Seringkali kita dihadapkan pada data kategorikal atau data cacah (counts).

Materi ini akan membahas tiga jenis regresi spesifik untuk menangani karakteristik data tersebut:

  1. Regresi Logistik Binari: Untuk variabel dependen kategorikal nominal (2 kategori)

  2. Regresi Logistik Multinomial: Untuk variabel dependen kategorikal nominal (> 2 kategori).

  3. Regresi Logistik Ordinal: Untuk variabel dependen kategorikal bertingkat/memiliki urutan.

  4. Regresi Poisson: Untuk variabel dependen berupa data cacah (count data) dalam rentang waktu/ruang tertentu.

Regresi Logistik Biner

Konsep Dasar

Regresi Logistik Biner adalah metode statistik yang digunakan ketika variabel dependen (\(Y\)) bersifat kategorikal dikotomi (hanya memiliki dua kemungkinan hasil, seperti “Ya” atau “Tidak”, “Sukses” atau “Gagal”, “Default” atau “Lancar”).

Berbeda dengan Regresi Linear Klasik yang memprediksi nilai kontinu, Regresi Logistik memprediksi probabilitas terjadinya suatu peristiwa (\(P(Y=1)\)). Karena nilai probabilitas harus berada di rentang 0 hingga 1, model ini menggunakan fungsi logistik atau fungsi sigmoid untuk mentransformasikan persamaan linear.

Formula Dasarnya:

\[\ln\left(\frac{P}{1-P}\right)=\beta_0+\beta_1+\beta_2+\dots+\beta_pX_p\]

Di mana \(\frac{P}{1-P}\) disebut sebagai Odds (rasio probabilitas kejadian sukses dibanding gagal).

Contoh Kasus di R

Seorang peneliti ingin membuat model prediksi hujan/tidak hujan di kota Bandung. Model tersebut menggunakan variabel dependen, yaitu curah hujan. Sedangkan, variabel independen sebagai berikut:

  • Arah angin saat kecepatan maksimum

  • Temperatur minimum

  • Lamanya penyinaran matahari 

  • Kelembapan rata-rata

Sumber: https://dataonline.bmkg.go.id/data-harian

Estimasi Parameter Model Regresi Logistik Binari

library(readxl)
Binom <- read_excel("C:/A MAIN FOLDER/Download/laporan_iklim_harian-260531173443.xlsx")
colnames(Binom) <- c("X1","Y","X2","X3","X4")
Binom$Kategori <- ifelse(Binom$Y==0,0,1) 
#Menjadikan Tidak hujan (0) sebagai baseline

#Pemodelan
model_binom <- glm(Kategori~X1+X2+X3+X4,data=Binom)
summary(model_binom)
## 
## Call:
## glm(formula = Kategori ~ X1 + X2 + X3 + X4, data = Binom)
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  1.8062445  2.0975110   0.861  0.39735   
## X1          -0.0307491  0.1072554  -0.287  0.77671   
## X2           0.0077418  0.0591584   0.131  0.89693   
## X3          -0.2308821  0.0826126  -2.795  0.00983 **
## X4           0.0005431  0.0009362   0.580  0.56704   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.1356958)
## 
##     Null deviance: 4.8000  on 29  degrees of freedom
## Residual deviance: 3.3924  on 25  degrees of freedom
## AIC: 31.746
## 
## Number of Fisher Scoring iterations: 2
#Pemodelan stepwise
stepwise_binom <- step(model_binom,direction = "both")
## Start:  AIC=31.75
## Kategori ~ X1 + X2 + X3 + X4
## 
##        Df Deviance    AIC
## - X2    1   3.3947 29.767
## - X1    1   3.4035 29.845
## - X4    1   3.4381 30.148
## <none>      3.3924 31.746
## - X3    1   4.4523 37.903
## 
## Step:  AIC=29.77
## Kategori ~ X1 + X3 + X4
## 
##        Df Deviance    AIC
## - X1    1   3.4039 27.848
## - X4    1   3.4381 28.148
## <none>      3.3947 29.767
## + X2    1   3.3924 31.746
## - X3    1   4.6732 37.356
## 
## Step:  AIC=27.85
## Kategori ~ X3 + X4
## 
##        Df Deviance    AIC
## - X4    1   3.4417 26.179
## <none>      3.4039 27.848
## + X1    1   3.3947 29.767
## + X2    1   3.4035 29.845
## - X3    1   4.7801 36.034
## 
## Step:  AIC=26.18
## Kategori ~ X3
## 
##        Df Deviance    AIC
## <none>      3.4417 26.179
## + X4    1   3.4039 27.848
## + X1    1   3.4381 28.148
## + X2    1   3.4412 28.175
## - X3    1   4.8000 34.159
summary(stepwise_binom)
## 
## Call:
## glm(formula = Kategori ~ X3, data = Binom)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.3039     0.1645   7.924 1.25e-08 ***
## X3           -0.2191     0.0659  -3.324  0.00248 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.1229177)
## 
##     Null deviance: 4.8000  on 29  degrees of freedom
## Residual deviance: 3.4417  on 28  degrees of freedom
## AIC: 26.179
## 
## Number of Fisher Scoring iterations: 2

Koefisien untuk \(X_3\)bernilai negatif, yaitu -0.2191, dengan tingkat signifikansi yang sangat tinggi (p-value = 0.00248). Artinya: Terdapat hubungan berkebalikan (negatif) yang signifikan antara kecepatan angin maksimum dengan curah hujan.

Penentuan Threshold Optimal

library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# Probabilitas prediksi untuk semua observasi
p_all <- predict(stepwise_binom, type = "response")

# ==============================
# Fungsi ROC dan AUC
# ==============================
safe_div <- function(x, y) ifelse(y == 0, 0, x / y)

roc_points <- function(actual, prob) {
  thresholds <- c(Inf, sort(unique(prob), decreasing = TRUE), -Inf)
  
  out <- lapply(thresholds, function(th) {
    pred <- as.integer(prob >= th)
    tp <- sum(pred == 1 & actual == 1)
    tn <- sum(pred == 0 & actual == 0)
    fp <- sum(pred == 1 & actual == 0)
    fn <- sum(pred == 0 & actual == 1)
    sensitivity <- safe_div(tp, tp + fn)
    specificity <- safe_div(tn, tn + fp)
    data.frame(
      threshold = th,
      sensitivity = sensitivity,
      specificity = specificity,
      fpr = 1 - specificity,
      youden = sensitivity + specificity - 1   # <-- tambahkan ini
    )
  })
  bind_rows(out)
}

auc_value <- function(roc_df) {
  roc_df <- roc_df[order(roc_df$fpr, roc_df$sensitivity), ]
  sum(diff(roc_df$fpr) * (head(roc_df$sensitivity, -1) + tail(roc_df$sensitivity, -1)) / 2)
}

# ==============================
# Hitung ROC dan AUC untuk seluruh data
# ==============================
roc_all <- roc_points(Binom$Kategori, p_all)
auc_all <- auc_value(roc_all)

cat("AUC seluruh data:", round(auc_all, 3), "\n")
## AUC seluruh data: 0.868
# ==============================
# Plot kurva ROC
# ==============================
ggplot(roc_all, aes(x = fpr, y = sensitivity)) +
  geom_line(color = "blue", size = 1) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", alpha = 0.5) +
  labs(x = "False Positive Rate (1 - Spesifisitas)",
       y = "True Positive Rate (Sensitivitas)",
       title = paste("Kurva ROC - Seluruh Data (AUC =", round(auc_all, 3), ")")) +
  theme_minimal()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Asumsi: Binom dan stepwise_binom sudah ada
p_all <- predict(stepwise_binom, type = "response")

# Hitung ROC points
roc_all <- roc_points(Binom$Kategori, p_all)

# Ambil threshold optimal (maksimum Youden)
threshold_opt <- roc_all %>%
  filter(is.finite(threshold)) %>%
  arrange(desc(youden)) %>%
  slice(1) %>%
  pull(threshold)

# Tampilkan
cat("Threshold optimal (Youden) =", round(threshold_opt, 4))
## Threshold optimal (Youden) = 0.8657


Nilai threshold optimal untuk model ini adalah 0.8657, sehingga nilai threshold ini digunakan untuk pemodelan.

Evaluasi Model Logistik Binari

# ==============================
# 1. Persiapan (pastikan model stepwise_binom sudah ada)
# ==============================
# Asumsi: Binom sudah diproses dan stepwise_binom sudah di-fit pada seluruh data
# Prediksi probabilitas untuk seluruh data
p_all <- predict(stepwise_binom, type = "response")

# ==============================
# 2. Tentukan threshold custom (ubah sesuai keinginan)
# ==============================
threshold_custom <- threshold_opt   # contoh: 0.5, bisa diganti 0.3, 0.7, dll.

# ==============================
# 3. Buat prediksi kelas berdasarkan threshold
# ==============================
pred_class <- ifelse(p_all >= threshold_custom, 1, 0)

# ==============================
# 4. Confusion matrix
# ==============================
conf_matrix <- table(Actual = Binom$Kategori, Predicted = pred_class)
print("Confusion Matrix:")
## [1] "Confusion Matrix:"
print(conf_matrix)
##       Predicted
## Actual  0  1
##      0  6  0
##      1  8 16
# ==============================
# 5. Hitung metrik evaluasi
# ==============================
tp <- conf_matrix[2,2]  # benar positif (actual=1, pred=1)
tn <- conf_matrix[1,1]  # benar negatif (actual=0, pred=0)
fp <- conf_matrix[1,2]  # false positive (actual=0, pred=1)
fn <- conf_matrix[2,1]  # false negative (actual=1, pred=0)

akurasi <- (tp + tn) / (tp + tn + fp + fn)
sensitivitas <- tp / (tp + fn)  # recall
spesifisitas <- tn / (tn + fp)
presisi <- tp / (tp + fp)
f1_score <- 2 * (presisi * sensitivitas) / (presisi + sensitivitas)

# Tampilkan metrik
cat("\n=== Metrik Evaluasi dengan Threshold =", threshold_custom, "===\n")
## 
## === Metrik Evaluasi dengan Threshold = 0.8657244 ===
cat("Akurasi        :", round(akurasi, 4), "\n")
## Akurasi        : 0.7333
cat("Sensitivitas   :", round(sensitivitas, 4), "\n")
## Sensitivitas   : 0.6667
cat("Spesifisitas   :", round(spesifisitas, 4), "\n")
## Spesifisitas   : 1
cat("Presisi        :", round(presisi, 4), "\n")
## Presisi        : 1
cat("F1-Score       :", round(f1_score, 4), "\n")
## F1-Score       : 0.8

Dengan threshold optimal sebesar 0,8657, model menunjukkan akurasi 73,33%, sensitivitas 66,67%, serta spesifisitas dan presisi sempurna (100%). Hal ini berarti model sangat konservatif: tidak ada kesalahan prediksi positif (false positive), namun sekitar sepertiga kejadian positif terlewat (false negative). Threshold ini cocok diterapkan jika biaya false positive jauh lebih tinggi daripada false negative, tetapi tidak ideal untuk situasi yang membutuhkan deteksi maksimal terhadap kasus positif.

Regresi Logistik Multinomial

Konsep Dasar

Regresi Multinomial digunakan ketika variabel dependen berbentuk kategorikal kualitatif (nominal) dengan lebih dari dua tingkatan, di mana urutan antar kategori tidak penting. Contohnya: Pemilihan moda transportasi (Bus, Kereta, Mobil Pribadi) atau pemilihan jurusan kuliah.

Model ini membandingkan setiap kategori dengan satu kategori acuan (Reference Cell). Jika ada \(K\) kategori, maka akan terbentuk \(K-1\) persamaan logit.

Formula dasarnya:

\[ln(\frac{P(Y=j)}{P(Y=Reference)})=\beta_0+\beta_1X_1+...\beta_nX_n\]

Contoh Kasus di R

Dataset Iris mengukur karakteristik morfologi dari 150 kuntum bunga Iris. Data ini dibagi rata menjadi 3 spesies bunga Iris, yaitu:

  • Iris setosa
  • Iris versicolor
  • Iris virginica

di mana masing-masing spesies memiliki tepat 50 sampel.

Kemudian, tiap bunga diukur berdasarkan beberapa parameter, yaitu:

  • Sepal Length: Panjang kelopak pelindung (daun kelopak).

  • Sepal Width: Lebar kelopak pelindung.

  • Petal Length: Panjang mahkota bunga.

  • Petal Width: Lebar mahkota bunga.

untuk dilakukan pemodelan karakteristik antar bunga, sehingga model dapat mengidentifikasi spesies yang paling cocok untuk bunga tersebut.

Sumber: Fisher, R. (1936). Iris [Dataset]. UCI Machine Learning Repository. https://doi.org/10.24432/C56C76.

Estimasi Parameter Model Regresi Logistik Multinomial

Menggunakan package nnet dengan fungsi multinom(). Kemudian, model akan dicari versi terbaiknya menggunakan metode stepwise

#Menyiapkan Library
library(readxl)
library(nnet)
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
#Persiapan Pemodelan
Multinomial <- read_excel("C:/A MAIN FOLDER/Download/IrisDataset.xlsx")
colnames(Multinomial) <- c("X1","X2","X3","X4","Y")
Multinomial$Kategori <- ifelse(Multinomial$Y=="Iris-setosa",0,
                               ifelse(Multinomial$Y=="Iris-versicolor",1,2))

#Pemodelan
model_multi <- nnet::multinom(Kategori~X1+X2,data = Multinomial)
## # weights:  12 (6 variable)
## initial  value 164.791843 
## iter  10 value 97.897928
## iter  20 value 78.633919
## iter  30 value 76.280985
## iter  40 value 72.686817
## iter  50 value 71.916797
## iter  60 value 71.182614
## iter  70 value 71.079609
## iter  80 value 71.005860
## iter  90 value 70.945999
## iter 100 value 70.819731
## final  value 70.819731 
## stopped after 100 iterations
summary(model_multi)
## Call:
## nnet::multinom(formula = Kategori ~ X1 + X2, data = Multinomial)
## 
## Coefficients:
##   (Intercept)       X1        X2
## 1   -10.68210 155.8685 -186.2485
## 2   -20.78525 190.0916 -169.5363
## 
## Std. Errors:
##   (Intercept)       X1       X2
## 1    6.924874 64.13583 64.15745
## 2    7.315896 64.71516 64.80320
## 
## Residual Deviance: 141.6395 
## AIC: 153.6395
#Stepwise Method
multi_stepwise <- step(model_multi, direction = "both")
## Start:  AIC=153.64
## Kategori ~ X1 + X2
## 
## trying - X1 
## # weights:  9 (4 variable)
## initial  value 164.791843 
## iter  10 value 137.994299
## iter  20 value 133.735491
## iter  30 value 133.165322
## iter  40 value 133.087049
## iter  50 value 133.065998
## iter  60 value 133.056589
## iter  70 value 133.053271
## iter  80 value 133.051382
## iter  90 value 133.049843
## final  value 133.049832 
## converged
## trying - X2 
## # weights:  9 (4 variable)
## initial  value 164.791843 
## iter  10 value 111.210847
## iter  20 value 106.674908
## iter  30 value 105.631820
## iter  40 value 105.239072
## iter  50 value 105.128764
## iter  60 value 105.081199
## iter  70 value 105.029927
## iter  80 value 105.016157
## iter  90 value 104.994082
## final  value 104.987040 
## converged
##        Df      AIC
## <none>  6 153.6395
## - X2    4 217.9741
## - X1    4 274.0997
summary(multi_stepwise)
## Call:
## nnet::multinom(formula = Kategori ~ X1 + X2, data = Multinomial)
## 
## Coefficients:
##   (Intercept)       X1        X2
## 1   -10.68210 155.8685 -186.2485
## 2   -20.78525 190.0916 -169.5363
## 
## Std. Errors:
##   (Intercept)       X1       X2
## 1    6.924874 64.13583 64.15745
## 2    7.315896 64.71516 64.80320
## 
## Residual Deviance: 141.6395 
## AIC: 153.6395

Interpretasi Model Regresi Logistik

\(X_1\) (Sepal Length)

  • Kategori 1 vs Kategori 0: Koefisiennya adalah 155.87 (Positif). Setiap kenaikan satu satuan Sepal Length, maka nilai log-odds bunga tersebut untuk masuk ke Kategori 1 (bukan Kategori 0) akan meningkat sebesar 155.87. Artinya, semakin panjang sepal bunga, semakin besar peluangnya merupakan Kategori 1 dibandingkan Kategori 0.

  • Kategori 2 vs Kategori 0: Koefisiennya adalah 190.09 (Positif). Polanya sama dan bahkan lebih kuat. Semakin panjang sepal bunga, peluang bunga tersebut untuk diklasifikasikan sebagai Kategori 2 dibandingkan Kategori 0 akan meningkat sangat drastis.

\(X_2\) (Sepal Width)

  • Kategori 1 vs Kategori 0: Koefisiennya adalah -186.25 (Negatif). Setiap kenaikan satu satuan Sepal Width, nilai log-odds bunga tersebut untuk masuk ke Kategori 1 dibandingkan Kategori 0 akan menurun sebesar 186.25. Artinya, bunga dengan sepal yang lebar cenderung merupakan Kategori 0.

  • Kategori 2 vs Kategori 0: Koefisiennya adalah -169.54 (Negatif). Peningkatan Sepal Width juga secara signifikan menurunkan peluang bunga masuk ke Kategori 2 dibandingkan Kategori 0.

Kesimpulan Karakteristik

Berdasarkan model terbaik pilihan stepwise ini, kita bisa menyimpulkan sifat fisik dari masing-masing kategori bunga Iris Anda:

  • Kategori 0 (Referensi): Cenderung memiliki Sepal yang Pendek (\(X_1\)) tetapi Lebar (\(X_2\)). (Karakteristik ini sangat identik dengan spesies Setosa pada dataset Iris asli).

  • Kategori 1: Cenderung memiliki Sepal yang Panjang (\(X_1\)) tetapi Ramping/Sempit (\(X_2\)).

  • Kategori 2: Cenderung memiliki Sepal yang Sangat Panjang (\(X_1\)) tetapi Ramping/Sempit (\(X_2\)). (Karakteristik Sepal Length yang paling ekstrem di antara ketiganya).

Evaluasi Akurasi Model Logistik Multinomial

Multinomial$Prediksi <- predict(multi_stepwise,type="class")
conf_multi <- table(Aktual = Multinomial$Kategori,
                    Prediksi = Multinomial$Prediksi)
conf_multi
##       Prediksi
## Aktual  0  1  2
##      0 50  0  0
##      1  6 30 14
##      2  1 11 38
accuracy_multi <- sum(diag(conf_multi))/sum(conf_multi)
accuracy_multi
## [1] 0.7866667

Dengan akurasi mencapai 78.67% hanya menggunakan 2 variabel saja (\(X_1\) dan \(X_2\)) setelah proses stepwise, model ini sudah tergolong cukup baik dan efisien untuk klasifikasi data biologi.

Regresi Logistik Ordinal

Konsep Dasar

Regresi Ordinal (sering menggunakan model Proportional Odds) digunakan apabila variabel dependen berbentuk kategorikal bertingkat/memiliki urutan, namun jarak antar tingkatan tidak dapat diukur secara kuantitatif. Contohnya: Tingkat kepuasan (Puas, Cukup, Tidak Puas) atau tingkat keparahan penyakit.

Formulasi Model:

\[P(Y \le j) = P(Y=1) + P(Y=2) + \dots+P(Y=j)\]

Contoh Kasus di R

Misalkan kita ingin melihat pengaruh variabel independen, sebagai berikut:

  • pared (Apakah minimal salah satu orang tua lulus kuliah: 0 = Tidak, 1 = Ya)

  • public (Jenis sekolah menengah: 0 = Swasta, 1 = Negeri)

  • gpa (Indeks Prestasi Kumulatif / IPK) terhadap

Terhadap apply (Kemungkinan mendaftar, memiliki 3 urutan: 0 = Unlikely, 1 = Somewhat Likely, 2 = Very Likely)

Sumber: https://stats.oarc.ucla.edu/r/dae/ordinal-logistic-regression/

Estimasi Parameter Model Regresi Logistik Ordinal

library(foreign)
Ordinal <- foreign::read.dta("C:/A MAIN FOLDER/Download/ologit.dta")
Ordinal$Kategori <- ifelse(Ordinal$apply==0,"unlikely",
                           ifelse(Ordinal$apply==1,"somewhat likely","very likely"))
colnames(Ordinal) <- c("Kategori","X1","X2","X3","Y")

#Pemodelan
library(MASS)
model_ordinal <- polr(as.factor(Kategori)~X1+X2+X3,data=Ordinal,
                      method = "logistic",
                      Hess = TRUE)
summary(model_ordinal)
## Call:
## polr(formula = as.factor(Kategori) ~ X1 + X2 + X3, data = Ordinal, 
##     Hess = TRUE, method = "logistic")
## 
## Coefficients:
##       Value Std. Error t value
## X1  1.04769     0.2658  3.9418
## X2 -0.05879     0.2979 -0.1974
## X3  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

Interpretasi Variabel Independen (Dampak terhadap Peluang Mendaftar)

Untuk mempermudah interpretasi regresi logistik, koefisien log-odds pada output dapat diubah menjadi nilai Odds Ratio (OR) dengan rumus \(OR = e^{\text{Value}}\).

\(X_1\) - pared (Pendidikan Orang Tua)

  • Koefisien: \(1.04769\) (Positif dan Signifikan, \(t = 3.94\))

  • Interpretasi Kontekstual: Siswa yang memiliki minimal salah satu orang tua lulusan kuliah (pared = 1) memiliki kecenderungan yang jauh lebih tinggi untuk mendaftar ke jenjang pendidikan berikutnya ke kategori yang lebih tinggi (Somewhat Likely atau Very Likely) dibandingkan siswa yang orang tuanya tidak lulus kuliah (pared = 0).

  • Nilai Odds Ratio: \(e^{1.04769} \approx 2.85\). Artinya, siswa dengan orang tua lulusan kuliah memiliki peluang (odds) \(2.85\) kali lebih besar untuk memilih kategori pendaftaran yang lebih tinggi dibandingkan yang tidak, dengan asumsi variabel lain konstan.

\(X_2\) - public (Jenis Sekolah)

  • Koefisien: \(-0.05879\) (Negatif dan Tidak Signifikan, \(t = -0.1974\))

  • Interpretasi Kontekstual: Jenis sekolah menengah (Negeri vs Swasta) tidak memiliki pengaruh yang nyata atau signifikan terhadap keputusan siswa untuk mendaftar. Perbedaan apakah siswa tersebut lulusan sekolah negeri (public = 1) atau swasta (public = 0) tidak mengubah kecenderungan mereka dalam memilih tingkat apply.

\(X_3\) - gpa (IPK / Nilai Akademik)

  • Koefisien: \(0.61594\) (Positif dan Signifikan, \(t = 2.36\))

  • Interpretasi Kontekstual: IPK atau prestasi akademik berpengaruh positif dan signifikan terhadap keputusan mendaftar. Semakin tinggi IPK seorang siswa, semakin tinggi pula kemungkinannya untuk memilih kategori apply yang lebih tinggi.

  • Nilai Odds Ratio: \(e^{0.61594} \approx 1.85\). Artinya, untuk setiap kenaikan \(1.0\) poin pada IPK (GPA), peluang (odds) siswa untuk memilih kategori pendaftaran yang lebih tinggi meningkat sebesar \(1.85\) kali lipat.

Nilai Ambang Batas (Intercepts / Cutpoints)

Nilai intercepts memotong fungsi distribusi kumulatif laten menjadi 3 bagian kategori respon \(Y\):

  • Titik Potong 1 (unlikely | somewhat likely = \(2.2039\)): Menunjukkan batas psikologis atau akademis di mana siswa bergeser dari sikap “Pasti Tidak Mendaftar” (Unlikely) menjadi “Agak Mungkin Mendaftar” (Somewhat Likely).

  • Titik Potong 2 (somewhat likely | very likely = \(4.2994\)): Menunjukkan batas yang lebih tinggi di mana siswa bergeser dari status “Agak Mungkin” menjadi sangat optimis dan yakin untuk mendaftar (Very Likely).

Evaluasi Akurasi Model Logistik Ordinal

Ordinal$Prediksi <- predict(model_ordinal,type="class")
conf_ord <- table(Aktual = Ordinal$Kategori,
                    Prediksi = Ordinal$Prediksi)
conf_ord
##                  Prediksi
## Aktual            unlikely somewhat likely very likely
##   unlikely             201              19           0
##   somewhat likely      110              30           0
##   very likely           27              13           0
accuracy_ord <- sum(diag(conf_ord))/sum(conf_ord)
accuracy_ord
## [1] 0.5775

Secara keseluruhan, model Anda mampu menebak dengan benar status pendaftaran siswa sebesar 57.75%. Angka ini sebenarnya berada di atas tebakan acak (33.3% untuk 3 kelas), namun masih memiliki kelemahan besar.

Asumsi Parallel Lines Pada Model Ordinal

# Pastikan model_ordinal sudah ada
# Hitung linear predictor (Xb) dari model
Xb <- model_ordinal$lp  # nilai linear predictor untuk setiap observasi

# Ambil cutpoints (zeta)
zeta <- model_ordinal$zeta

# Buat data frame untuk plotting cumulative logit
# Untuk setiap batas (j), hitung P(Y <= j) dan logit-nya
cum_logit_data <- data.frame(
  Xb = sort(Xb),
  logit_cat1 = qlogis(plogis(sort(Xb) - zeta[1])),  # logit P(Y <= cat1)
  logit_cat2 = qlogis(plogis(sort(Xb) - zeta[2]))   # logit P(Y <= cat2)
)

# Reshape untuk ggplot
library(tidyr)
cum_logit_long <- pivot_longer(cum_logit_data, 
                               cols = c(logit_cat1, logit_cat2),
                               names_to = "Batas",
                               values_to = "logit_kumulatif")

# Plot
library(ggplot2)
ggplot(cum_logit_long, aes(x = Xb, y = logit_kumulatif, color = Batas)) +
  geom_line(linewidth = 1.2) +
  scale_color_manual(values = c("#2f7f73", "#d18b2f"),
                     labels = c("P(Y ≤ unlikely)", "P(Y ≤ somewhat likely)")) +
  labs(title = "Uji Asumsi Parallel Lines (Proportional Odds)",
       subtitle = "Garis cumulative logit terhadap linear predictor harus sejajar",
       x = "Linear predictor (Xβ)",
       y = "Cumulative logit",
       color = "Batas kumulatif") +
  theme_minimal() +
  theme(legend.position = "bottom")

Secara visual, kedua garis kumulatif logit tersebut memiliki kemiringan (slope) yang identik atau sejajar. Ketika nilai Linear predictor meningkat, kedua kurva merespons dengan tingkat kenaikan yang sama.

Regresi Poisson

Konsep Dasar

Regresi Poisson digunakan untuk data cacah (count data), yaitu data berupa frekuensi terjadinya suatu peristiwa dalam interval waktu atau ruang tertentu. Nilai data ini selalu berupa bilangan bulat non-negatif (\(0, 1, 2, \dots\)).

Syarat utama Regresi Poisson adalah nilai Mean harus sama dengan Variance (\(\mu = \sigma^2\)), kondisi ini disebut dengan Equidispersion.

Formulasi Model:

\[\ln(\lambda)=\beta_0+\beta_1X_1+\dots+\beta_pX_p\]

Di mana \(\lambda\) adalah rata-rata jumlah kejadian (expected count).

Contoh Kasus di R

Seorang peneliti ingin memodelkan berapa banyak penghargaan yang dimenangkan oleh siswa di sebuah sekolah berdasarkan beberapa faktor, yaitu :

  • prog (Jenis program studi: Vocational, General, Academic)

  • math (Nilai ujian matematika)

Sumber Data: https://stats.oarc.ucla.edu/stata/dae/poisson-regression/

Estimasi Parameter Model Regresi Logistik Poisson

Poisson <- read.csv("C:/A MAIN FOLDER/Download/poisson_sim.csv")
Poisson <- Poisson[,2:4]
colnames(Poisson) <- c("Y","X1","X2")
Poisson$X1_cat <- ifelse(Poisson$X1==1,"vocational",
                     ifelse(Poisson$X1==2,"general","academic")) 
#Baseline = Academic

#Pemodelan
model_poisson <- glm(Y~as.factor(X1_cat)+X2,
                     data=Poisson,
                     family = poisson(link="log"))
summary(model_poisson)
## 
## Call:
## glm(formula = Y ~ as.factor(X1_cat) + X2, family = poisson(link = "log"), 
##     data = Poisson)
## 
## Coefficients:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                 -4.87732    0.62818  -7.764 8.21e-15 ***
## as.factor(X1_cat)general     0.71405    0.32001   2.231   0.0257 *  
## as.factor(X1_cat)vocational -0.36981    0.44107  -0.838   0.4018    
## X2                           0.07015    0.01060   6.619 3.63e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 287.67  on 199  degrees of freedom
## Residual deviance: 189.45  on 196  degrees of freedom
## AIC: 373.5
## 
## Number of Fisher Scoring iterations: 6

Interpretasi Model Regresi Poisson

\(X_1\) (Jenis Program Studi)

(Menggunakan Program Academic sebagai Kategori Referensi/Basis)

  • Kategori General vs Kategori Academic: Koefisiennya adalah 0.71405 (Positif). Siswa dari program studi General memiliki nilai log-expected count jumlah penghargaan yang lebih tinggi sebesar 0.71405 dibandingkan siswa dari program Academic. Nilai p-value (\(0.0257 < 0.05\)) menunjukkan bahwa perbedaan ini signifikan secara statistik. Artinya, siswa program General cenderung memenangkan penghargaan lebih banyak daripada siswa program Academic pada tingkat nilai matematika yang sama.

  • Kategori Vocational vs Kategori Academic: Koefisiennya adalah -0.36981 (Negatif). Siswa dari program studi Vocational memiliki nilai log-expected count jumlah penghargaan yang lebih rendah sebesar 0.36981 dibandingkan siswa dari program Academic. Namun, karena p-value besar (\(0.4018 > 0.05\)), perbedaan antara program Vocational dan Academic ini tidak signifikan secara statistik.

\(X_2\) (Nilai Ujian Matematika)
  • Pengaruh Nilai Matematika: Koefisiennya adalah 0.07015 (Positif). Setiap kenaikan satu satuan nilai ujian matematika, maka nilai log-expected count jumlah penghargaan siswa akan meningkat sebesar 0.07015. Nilai p-value yang sangat kecil (\(3.63 \times 10^{-11}\)) menunjukkan bahwa pengaruh ini sangat signifikan. Artinya, tanpa memandang jenis program studinya, semakin tinggi nilai matematika seorang siswa, peluang untuk mendapatkan jumlah penghargaan yang lebih banyak akan meningkat secara nyata.
Kesimpulan Karakteristik

Berdasarkan model terbaik dari hasil analisis regresi Poisson ini, kita bisa menyimpulkan karakteristik siswa dalam memenangkan penghargaan di sekolah Anda:

  • Kategori Program Vocational: Cenderung memiliki jumlah perolehan penghargaan yang paling rendah di antara ketiga program studi (ditunjukkan oleh koefisien negatif), dan jumlah penghargaan mereka sangat bergantung pada peningkatan nilai matematika mereka sendiri.

  • Kategori Program Academic: Memiliki tingkat perolehan penghargaan yang berada di posisi menengah. Jumlah penghargaan mereka lebih baik daripada siswa Vocational, namun masih berada di bawah siswa program General.

  • Kategori Program General: Cenderung menjadi program studi yang paling dominan dalam memenangkan penghargaan (memiliki koefisien positif terbesar dan signifikan). Siswa dari program General dengan Nilai Ujian Matematika (\(X_2\)) yang tinggi merupakan karakteristik kelompok siswa yang paling ekstrem dalam mengumpulkan penghargaan di sekolah tersebut.

Pemeriksaan Asumsi Regresi Poisson

library(dplyr)
library(knitr)
library(kableExtra)
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
dispersion_pois <- sum(residuals(model_poisson, type = "pearson")^2) / df.residual(model_poisson)

tibble::tibble(
  `Dispersion Pearson` = dispersion_pois,
  `Interpretasi ringkas` = dplyr::case_when(
    dispersion_pois < 1.5 ~ "Tidak ada indikasi overdispersion berat",
    dispersion_pois < 2.5 ~ "Ada indikasi overdispersion sedang",
    TRUE ~ "Ada indikasi overdispersion kuat"
  )
) %>%
  mutate(`Dispersion Pearson` = round(`Dispersion Pearson`, 3)) %>%
  kable(
    caption = "Indikasi overdispersion pada model Poisson"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )
Indikasi overdispersion pada model Poisson
Dispersion Pearson Interpretasi ringkas
1.082 Tidak ada indikasi overdispersion berat

Karena nilai dispersi Pearson mendekati 1, tidak ada indikasi overdispersion. Hal ini menandakan bahwa variasi data jumlah penghargaan terkontrol dengan baik oleh model Poisson dan standar error yang dihasilkan sudah akurat.