1 Pendahuluan

  • Sebagai awalan, coretan ini berfokus ke pembelajaran (sehingga ukuran metrik evaluasi yang mungkin tidak terlalu bagus bisa diabaikan, makna diabaikan disini jangan terlalu terfokus pada nilai metrik evaluasi untuk coretan ini melainkan bagaimana step by step bagaimana mendapatkan interpretasi hingga metrik evaluasi apabila kita menggunakan regresi logistik).

  • Oh ya, coretan ini sebagai salah satu sarana untuk berbagi dan mungkin tak lupt dari kesalahan, sehingga kalau ada perbaikan atau kurang sesuai, silahkan bernada pada email kami: “”

  • Analisis Regresi Logistik adalah metode statistika untuk memahami hubungan antara variabel independen/ prediktor/ bebas/ feature dan variabel dependen/ respon/ terikat/ target biner.

  • Makna biner dalam konteks analisis regresi logistik adalah bahwa variabel dependen hanya memiliki dua kemungkinan kategori atau nilai.

  • Contoh:

    1. Ya (1) / Tidak (0) → Misalnya, apakah seseorang lulus ujian atau tidak.

    2. Sakit (1) / Sehat (0) → Misalnya, apakah seseorang terkena penyakit tertentu atau tidak.

    3. Membeli (1) / Tidak membeli (0) → Misalnya, apakah seseorang akan membeli produk atau tidak.

  • Pada regresi logistik, kita akan mengenal istilah kategori referensi. Kategori referensi merupakan kategori dasar. Cara menentukan kategori referensi diantaranya:

    1. Kategori yang dianggap baseline standar dalam penelitian

    2. Kategori paling umum

    3. Kategori dengan makna khusus (misalnya tidak mendapatkan perlakukan).

  • Dalam regresi logistik biner, model digunakan untuk memperkirakan probabilitas suatu kejadian terjadi berdasarkan variabel independen. Variabel dependen biner biasanya memiliki dua nilai, seperti 0 dan 1.

  • Proses analisis regresi logistik melibatkan pembuatan model matematika yang memetakan nilai-nilai variabel independent ke probabilitas variabel dependent yang berada dalam rentang antara 0 dan 1. Salah satu fungsi yang digunakan pada model ini adalah logit atau log-odds untuk memodelkan hubungan.

  • Selain fungsi logit, pembaca bisa menggunakan fungsi hubung yang lain seperti probit, dll (silahkan membaca Agresti (link)).

  • Fungsi logit mengubah probabilitas ke dalam rentang nilai yang tidak terbatas. Tujuan fungsi logit adalah memungkinkan penggunaan metode statistika konvensional.

  • Metode yang digunakan adalah maximum likelihood estimation untuk menentukan parameter-model yang paling sesuai dengan data.

  • Distribusi yang diikuti oleh variabel dependen biner adalah distribusi Bernoulli.

  • Fungsi kepekatan peluang (PDF) dari distribusi Bernoulli digunakan dalam analisis regresi logistik.

\[f(y_i)=\pi (x_i)^{y_i} (1-\pi(x_i))^{1-y_i}, y_i = 0,1\] Keterangan

\(y_i\): Variabel respon ke-\(i\) (dengan nilai 0 atau 1)

\(\pi_i\): Peluang kejadian ke-\(i\).

Secara umum, model regresi logistik dapat ditulis

\[\mu = E[y] = \log \left[ \frac{\pi (x)}{1-\pi (x)} \right] = \beta_0 + \beta_1 X + \beta_2 X+ \cdots + \beta_p X_p \]

\[ \pi (x) = P(Y=1 | \pi(x)) \] dan

\[1 - \pi (x) = 1 - P(Y=1 | \pi(x)) = P (Y=0 | \pi(x))\]

atau dapat juga ditulis

\[\pi(x_i) = \frac{\exp(\beta_0 + \beta_1x_{1i} + \beta_2x_{2i}+\cdots + \beta_px_{pi})}{1+\exp(\beta_0 + \beta_1x_{1i} + \beta_2x_{2i}+\cdots + \beta_px_{pi}}\] Keterangan:

\(\pi(x_i)\): Peluang terjadinya kategori variabel respon.

\(x_{ji}\): variabel prediktor ke-\(j\).

\(p\): Banyaknya variabel prediktor.

\(\beta_0\): intercept.

\(\beta_0, \beta_1, \beta_2, \cdots, \beta_p\): koefisien regresi untuk setiap variabel prediktor.

  • Pemodelan Regresi Logistik Sederhana dengan variabel prediktor biner memerlukan reference category/ kategori pembanding/ kategori baseline. Kategori pembanding dikodekan sebagai \(0\), sedangkan kategori lainnya dikodekan sebagai \(1\).

  • Proses pengkodean ini merupakan transformasi variabel biner menjadi variabel numerik.

2 Install Package

rm(list=ls()) # Menghapus semua dataset dan variabel
graphics.off() # Menutup semua grafik

library(tidyverse)
library(broom)
library(broom.helpers)
library(skimr)
library(performance)
library(DataExplorer)
library(dplyr)
library(rsample)
library(yardstick)

3 Data

  • Data yang digunakan adalah data stillbirth di kabupaten Aceh Timur.

  • Stillbirth merupakan kelahiran bayi dalam keadaan tidak bernyawa yang telah mencapai umur kehamilan 20 minggu atau sebelum masa persalinan. Terdapat beberapa penyebab stillbirth yaitu infeksi selama kehamilan, kelainan atau cacat bawaan, kondisi kesehatan ibu, usia ibu, dan sebagainya.

  • Data diperoleh dari artikel dengan judul “Analisis Regresi Logistik Biner Untuk Menentukan Faktor Stillbirth Di Kabupaten Aceh Timur” (link artikel).

  • Data tersebut sudah kami siapkan dan dapat anda akses pada link berikut: (data).

data <- read.csv("data.csv", header = TRUE, sep = ",")
data <- data %>% select(-no)  # Menghapus kolom "
head(data)

Keterangan:

\(Y\) : Status kejadian stillbirth (0: tidak terjadi stillbirth, 1: terjadi stillbirth)

\(X_1\): Infeksi selama kehamilan (0 = tidak terjadinya infeksi kehamilan, 1 = terjadinya infeksi kehamilan)

\(X_2\): Kelainan atau cacat bawaan (0 = tidak mengalami kelainan atau cacat bawaan, 1 = mengalami kelainan atau cacat bawaan)

\(X_3\); Kondisi Ibu (0 = preklamsia, 1 = diabetes)

\(X_4\): Usia Ibu

  • Karena fungsi glm untuk membutuhkan semua variabel (kolom) kategorik (biner atau nominal) dalam tipe data factor maka kita akan konversi terlebih dahulu.
glimpse(data)
## Rows: 55
## Columns: 5
## $ Y  <int> 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1, 0…
## $ X1 <int> 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 0…
## $ X2 <int> 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1…
## $ X3 <int> 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0…
## $ X4 <int> 33, 38, 31, 32, 24, 24, 26, 36, 30, 25, 40, 31, 39, 25, 26, 38, 31,…
data <- data %>%
  mutate(across(c(Y, X1, X2,X3), as.factor))  # Mengubah hanya kolom Y, X1, dan X2 menjadi faktor
glimpse(data)  # Cek perubahan tipe data
## Rows: 55
## Columns: 5
## $ Y  <fct> 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1, 0…
## $ X1 <fct> 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 0…
## $ X2 <fct> 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1…
## $ X3 <fct> 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0…
## $ X4 <int> 33, 38, 31, 32, 24, 24, 26, 36, 30, 25, 40, 31, 39, 25, 26, 38, 31,…

3.1 Eksplorasi Umum

skim_without_charts(data = data)
Data summary
Name data
Number of rows 55
Number of columns 5
_______________________
Column type frequency:
factor 4
numeric 1
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
Y 0 1 FALSE 2 1: 30, 0: 25
X1 0 1 FALSE 2 1: 30, 0: 25
X2 0 1 FALSE 2 0: 37, 1: 18
X3 0 1 FALSE 2 0: 41, 1: 14

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100
X4 0 1 30.47 7.52 2 26 31 36 40

3.2 Korelasi

  1. Korelasi variabel prediktor (numerik) dengan variabel response
plot_boxplot(data = data,by = "Y",
             geom_boxplot_args = list(fill="#03A9F4"),
             ggtheme = theme_classic())

  1. Korelasi variabel prediktok (kategorik) dengan variabel response
plot_bar(data = data,by = "Y",ggtheme = theme_classic(),ncol = 2)

4 Splitting Data

  • Metode Holdout Sample digunakan untuk membagi dataset menjadi training data dan testing data.

  • Training data/ data latih digunakan untuk melatih model, sedangkan testing data/ data uji digunakan untuk mengevaluasi performa prediksi model.

  • Training data harus lebih banyak dibandingkan dengan testing data.

  • Pembagian dataset biasanya dinyatakan dalam bentuk proporsi atau persentase (misalnya \(0.9\) atau \(90\%\) untuk training dan \(0.1 \text{ atau} 10\%\) untuk testing).

  • Pembagian data dilakukan secara acak menggunakan random sampling.

  • Jika variabel respon bersifat kategorik, disarankan menggunakan stratified random sampling agar distribusi kategori tetap seimbang. Berikut sintaksnya.

library(dplyr)
set.seed(12)
intrain <- initial_split(data = data,prop = 0.9,strata =Y)
data_train <- training(intrain)
data_train
data_test <- testing(intrain)
data_test
print("Level, angka pertama adalah kelas positif/ kategori referensinya/ baselinnya")
## [1] "Level, angka pertama adalah kelas positif/ kategori referensinya/ baselinnya"
data$Y %>% levels
## [1] "0" "1"
  • Di R, kategori pembanding biasanya merupakan kategori pertama yang muncul dalam levels(). Ilustatrasinya adalah sebagai berikut.
levels(data$Y)
## [1] "0" "1"
levels(data$X1)
## [1] "0" "1"
levels(data$X2)
## [1] "0" "1"
levels(data$X3)
## [1] "0" "1"

Berdasarkan output diatas, penjelasannya adalah sebagai berikut.

  1. Pada variabel \(Y\), Kategori pembanding-nya/ kategori referensinya adalah “0” karena berada di urutan pertama (0: tidak terjadi stillbirth, 1: terjadi stillbirth)

  2. Pada variabel \(X_1\), Kategori pembanding-nya/ kategori referensinya adalah “0” karena berada di urutan pertama (0 = tidak terjadinya infeksi kehamilan, 1 = terjadinya infeksi kehamilan).

  3. Pada variabel \(X_2\), Kategori pembanding-nya/ kategori referensinya adalah “0” karena berada di urutan pertama (0 = tidak mengalami kelainan atau cacat bawaan, 1 = mengalami kelainan atau cacat bawaan).

  4. Pada variabel \(X_3\), Kategori pembanding-nya/ kategori referensinya adalah “0” karena berada di urutan pertama (0 = preklamsia, 1 = diabetes).

kita bisa merubah Kategori pembanding dengan fungsi relevel dan mengisi argumen ref dengan kategori yang akan dijadikan sebagai kategori pembanding. Berikut ilustrasinya

# coba <- relevel(x = data$Y,ref = "1")
# levels(coba)

5 Regresi Logistik

5.1 Pembentukan Model

reglog <- glm(Y=="1"~X1+X2+X3+X4, family = binomial(link = "logit"), data = data_train)
summary(reglog)
## 
## Call:
## glm(formula = Y == "1" ~ X1 + X2 + X3 + X4, family = binomial(link = "logit"), 
##     data = data_train)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -2.248445   1.774764  -1.267   0.2052  
## X11          2.921707   1.366034   2.139   0.0324 *
## X21          2.118471   1.408886   1.504   0.1327  
## X31          2.675139   1.256575   2.129   0.0333 *
## X4          -0.009075   0.044503  -0.204   0.8384  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 67.417  on 48  degrees of freedom
## Residual deviance: 51.821  on 44  degrees of freedom
## AIC: 61.821
## 
## Number of Fisher Scoring iterations: 5
  • Berdasarkan summary diatas model logit yang terbentuk adalah sebagai berikut.

\[Logit [\hat{\pi} (X)] = -2.248445 + 2.921707 X_{1(1)} + 2.118471X_{2(1)} + 2.675139 X_{3(1)} + 0.01942 X_4\]

  • sintaks Y==“1” berarti kita menenetapkan kelas positifnya adalah “1”atau \(P(Y=1|\pi(x))\), sedangkan kelas negatifnya adalah 0 atau \(P(Y=0|\pi(x))\)

  • Jika kita tidak mendefinisikan dengan sintaks Y==“1” maka secara default yang akan menjadi kelas positif adalah kategori pertama yang ditunjukkan dari output fungsi levels. Dalam kasus variabel respon Y, secara default yang menjadi kelas positif adalah 0 Berikut cara mengecek kategori pertama dari levels.

levels(data$Y)
## [1] "0" "1"

  • sintaks ~X1+X2+X3+X4 menandakan bahwa X1, X2, X3, dan X4 adalah variabel prediktor.

  • sintaks family(link = "logit") digunakan untuk menandakan kita menggunakan regresi logistik dengan link function adalah logit.

  • untuk mengeluarkan nilai dugaan koefisien, Selang kepercayaan (Confidence Interval) dan p-value dapat menggunakan fungsi parameters dari package parameters.

tidy_parameters(reglog,conf.int = TRUE,conf.level = 0.95,ci_method="wald")

  • Intercept umumnya dapat diinterpretasikan dengan melakukan transformasi terlebih dahulu ke dalam bentuk peluang dengan mengatur nilai peubah prediktor sama dengan 0. Perhitungan ini bisa menggunakan fungsi predict dan argumen dari predict yaitu type=“response”.

  • Koefisien dari prediktor X1 [1] (\(\beta_1\)) adalah 2.921706941 (positif) berarti saat seorang ibu memiliki infeksi kehamilan (X1) maka kemungkinan bayi akan terkena stillbirth semakin besar. \(p-val = 0.03244989 < 0.05 = \alpha\) untuk prediktor X1 [1] yang berarti sudah cukup bukti untuk menolak \(H_0:\beta_1 = 0\) yang berarti X1 [1] berpengaruh terhadap \(Y\). 95% CI untuk (\(\beta_1\)) dapat diinterpretasikan jika sudah kita transformasi ke dalam Odd-Ratio

  • Koefisien dari prediktor X2 [1] (\(\beta_2\)) adalah 2.118471137 (positif) berarti saat seorang ibu memiliki kelainan/ cacat bawaan maka kemungkinan bayi akan terkena stillbirth semakin besar. \(p-val = 0.13267149 > 0.05 = \alpha\) untuk prediktor X2 [1] yang berarti tidak cukup bukti untuk menolak \(H_0:\beta_2 = 0\) yang berarti X2 [1] tidak berpengaruh signifikan terhadap \(Y\).

  • Koefisien dari prediktor X3 [1] (\(\beta_3\)) adalah 2.675139441 (positif) berarti saat seorang ibu memiliki diabetes maka kemungkinan bayi akan terkena stillbirth semakin besar. \(p-val = 0.03326147 < 0.05 = \alpha\) untuk prediktor X3 [1] yang berarti cukup bukti untuk menolak \(H_0:\beta_3 = 0\) yang berarti X3 [1] berpengaruh signifikan terhadap \(Y\).

  • Koefisien dari prediktor X4 (\(\beta_4\)) adalah -0.009075333 (negatif) berarti saat seorang ibu semakin berumumr maka kemungkinan bayi akan terkena stillbirth semakin kecil, kendati demikian \(p-val = 0.83841197 > 0.05 = \alpha\) untuk prediktor X4 yang berarti tidak cukup bukti untuk menolak \(H_0:\beta_4 = 0\) yang berarti X4 tidak berpengaruh signifikan terhadap \(Y\).

  • Untuk pembelajaran ini variabel yang tidak signifikan tetap digunakan dengan pertimbangan variabel tersebut penting dalam penelitian dan sebagai bahan pembelajaran. Sementara untuk penelitian riilnya, anda dapat mempertimbangakan untuk tetap menggunakan atau dilakukan eliminasi menggunakan stepwise.

  • Untuk mendapatkan nilai Odds-Ratio dari koefisien dan intercept kita dapat menerapkan rumus berikut

\[\text{Odds-Ratio} (\beta) = \exp (\beta) \] sehingga

beta <- coef(reglog)
OR_beta <- exp(beta)
OR_beta
## (Intercept)         X11         X21         X31          X4 
##   0.1055633  18.5729634   8.3184101  14.5143736   0.9909657

atau lebih mudah jika kita menggunakan argumen exponentiate = TRUE pada fungsi parameters. Selain nilai koefisien dan intercept berubah menjadi Odds-Ratio, nilai standard error dan selang kepercayaan (CI) koefisien juga berubah menyesuaikan. Namun, nilai statistik z dan p-value tidak berubah.

tidy_parameters(reglog,conf.int = TRUE,conf.level = 0.95,ci_method="wald",exponentiate = TRUE) %>% 
  mutate(across(where(is.numeric),function(x) round(x,2) ))

Penjelasan

  • Odds-Ratio dari prediktor X1 [1] (\(\beta_1\)) adalah 18.57 berarti seorang ibu yang memiliki infeksi kehamilan mempunyai odds 1757 % (\([18.57-1]\times 100\% = 17.57 \times 100\% = 1757 \%\)) lebih tinggi dibandingkan dengan seorang ibu yang tidak memiliki infeksi kehamilan dalam potensi terjadinya stillbirth.

  • Odds-Ratio dari prediktor X2 [1] (\(\beta_2\)) adalah 8.32 berarti seorang ibu yang memiliki cacat bawaan mempunyai odds 732 % (\([8.32-1]\times 100\% = 7.32 \times 100 \% = 732 \%\)) lebih tinggi dibandingkan dengan seorang ibu yang tidak memiliki cacat bawaan dalam potensi terjadinya stillbirth.

  • Odds-Ratio dari prediktor X3 [1] (\(\beta_3\)) adalah 14.51 berarti seorang ibu yang memiliki diabetes mempunyai odds 578 % (\([6.78-1]\times 100\% = 5.14 \times 100 \% = 578 \%\)) lebih tinggi dibandingkan dengan seorang ibu yang tidak memiliki diabetes dalam potensi terjadinya stillbirth.

  • Odds-Ratio dari prediktor X4 (\(\beta_4\)) adalah 0.99 berarti semakin usia seorang ibu hampir tidak ada pengaruh. Jika dihitung nilai odds nya sebesar \(1 \% (1-0.99) \times 100 \% = 0.01 \times 100 \% = 1 \%\) lebih rendah dengan peningkatan umur.

  • Selanjutnya, misal kita hendak Peluang terjadinya stillbirth

    • Contoh 1: jika nilai \(X_1=0\) (tidak terjadi infeksi kehamilan), \(X_2=0\) (tidak memiliki cacat bawaan), \(X_3=0\) (tidak memiliki diabetes), dan berumur 40 tahun.
predict(object = reglog, newdata = data.frame(X1 = "0", X2 = "0", X3 = "0", X4 = 40),type = "response")
##          1 
## 0.06840462

Berdasarkan hasil diatas, diperoleh nilai peluang terjadinya stillbirth sebesar 0.06840462. Secara matematis dapat dihitung sebagai berikut. Persamaan probabilitas kejadian \(Y = 1\) dalam regresi logistik adalah:

\[\begin{align} \pi(x) &= \frac{\exp(\beta_0 + \beta_1 X_1 + \beta_2 X_2 + \beta_3 X_3 + \beta_4 X_4 )}{1 + \exp(\beta_0 + \beta_1 X_1 + \beta_2 X_2 + \beta_3 X_3 + \beta_4 X_4} \\ &= \frac{\exp(-2.248445 + 2.921707[0] + 2.118471 [0] + 2.675139[0] - 0.009075[40] )}{1 + \exp(-2.248445 + 2.921707[0] + 2.118471 [0] + 2.675139[0] - 0.009075[40] )} \\ &= 0.068405462 \end{align}\]

  • Contoh 2: jika nilai \(X_1=1\) (terjadi infeksi kehamilan), \(X_2=1\) (memiliki cacat bawaan), \(X_3=1\) ( memiliki diabetes), dan berumur 40 tahun.
predict(object = reglog, newdata = data.frame(X1 = "1", X2 = "1", X3 = "1", X4 = 40),type = "response")
##         1 
## 0.9939634

Berdasarkan hasil diatas, diperoleh nilai peluang terjadinya stillbirth sebesar 0.09833882. Secara matematis dapat dihitung sebagai berikut.

Persamaan probabilitas kejadian \(Y = 1\) dalam regresi logistik adalah:

\[\begin{align} \pi(x) &= \frac{\exp(\beta_0 + \beta_1 X_1 + \beta_2 X_2 + \beta_3 X_3 + \beta_4 X_4 )}{1 + \exp(\beta_0 + \beta_1 X_1 + \beta_2 X_2 + \beta_3 X_3 + \beta_4 X_4} \\ &= \frac{\exp(-2.248445 + 2.921707[1] + 2.118471 [1] + 2.675139[1] - 0.009075[40] )}{1 + \exp(-2.248445 + 2.921707[1] + 2.118471 [1] + 2.675139[1] - 0.009075[40] )} \\ &= 0.993963475 \end{align}\]

5.2 Goodness-of-fit Model

Goodness-of-fits Model biasa menggunakan ukuran-ukuran atau metrik-metrik tertentu yang didasarkan dari log-likehood. Ukuran AICc, BIC, RMSE biasanya digunakan untuk menentukan model terbaik dari beberapa model yang tersedia. Jika hanya satu model saja yang tersedia ukuran ini tidak bisa digunakan. Semakin kecil nila dari AICc, BIC, RMSE maka semakin baik pula modelnya.

Selain itu ada juga ukuran \(R^2\) yang bisa digunakan jika model yang tersedia hanya satu. Nilai berkisar antara 0 dan 1. Semakin mendekati 1 nilai semakin baik modelnya. Namun, dalam regresi logistik, \(R^2\) tidak didefinisikan secara langsung seperti pada regresi linear karena modelnya berbasis probabilitas (likelihood), bukan selisih kuadrat (sum of squares). Sebagai gantinya, beberapa pseudo-R² dikembangkan untuk regresi logistik, salah satunya adalah Tjur’s R².

Tjur’s R² (juga disebut sebagai Coefficient of Discrimination) mengukur perbedaan rata-rata probabilitas yang diprediksi antara kelas positif dan kelas negatif.

Formula \(Tjur's R^2\)

\[R^2_{Tjur}=\bar{\pi}_1-\bar{\pi}_0\]

\(\bar{\pi}_1\): rata-rata probabilitas prediksi untuk kelas sukses (\(Y=1\)).

\(\bar{\pi}_0\): rata-rata probabilitas prediksi untuk kelas \(Y=0\)

Semakin tinggi nilai \(R^2_{Tjur}\), semakin baik model dalam membedakan antara kelas 0 dan kelas 1. Jika \(R^2_{Tjur}\) = 0, maka model tidak memiliki daya diskriminasi. Jika \(R^2_{Tjur}\) mendekati 1, maka model sangat baik dalam memisahkan kelas.

Alternatif \(pseudo-R^2\) dalam regresi logistik

Metode Interpretasi
McFadden’s R² Menggunakan perbandingan likelihood antara model penuh dan model kosong. Nilai 0.2 - 0.4 dianggap baik.
Cox & Snell’s R² Berdasarkan likelihood ratio, tetapi nilainya tidak mencapai 1.
Nagelkerke’s R² Koreksi dari Cox & Snell’s R², sehingga bisa mencapai maksimum 1. 
performance(model = reglog)

Kemudian, Goodness-of-fits juga bisa menggunakan uji hipotesis statistik yaitu Likelihood Ratio Test. Secara umum, uji ini membandingkan dua model regresi dan menyimpulkan model mana yang terbaik.

Jika kita hanya memliki satu model saja maka model pembandingnya adalah null-model atau model tanpa variabel prediktor dan hanya memuat intercept saja.

Hipotesis dari Likelihood Ratio Test adalah sebagai berikut:

\[\begin{align} H_0 &: \beta_1 = \beta_2 = \beta_3 = \beta_4 \text{ null model fit (cocok untuk data)} \\ H_1 &: \exists j, \beta_j \neq 0, j=1,2,3,4 \text{ model yang diusulkan fit (cocok untuk data)} \end{align}\]

Likelihood Ratio Test Statistics \(G\) didefinisikan

\[G=-2 \cdot \log \left( \frac{L_0}{L_1}\right)\]

Di bawah hipotesis nol, \(G\) mengikuti distribusi chi-squared dengan derajat kebebasan yang sama dengan perbedaan jumlah parameter yang diduga pada kedua model. \(L_0\) adalah likelihood dari data pada model nol dan \(L_1\) adalah likelihood dari data pada model yang diusulkan.

Jika likelihood ratio test statistic \(G\) cukup besar, maka kita menolak hipotesis nol dan memilih hipotesis alternatif, yang menunjukkan bahwa parameter tambahan dalam model yang diusulkan secara signifikan meningkatkan fit (kecocokan untuk data).

Berikut adalah sintaksnya.

test_likelihoodratio(reglog)
## Only one model was provided, however, at least two are required for
##   comparison.
##   Fitting a null-model as reference now.

Berdasarkan hasil diatas, \(H_0\) ditolak atau model yang diusulkan fit (cocok untuk data).

5.3 Evaluasi Performa Prediksi

Berikut adalah sintaks untuk mendapatkan prediksi regresi logistik menggunakan testing data

pred_reglog0 <- predict(reglog,newdata = data_test,type = "response")
head(pred_reglog0)
##         1         2         3         4         5         6 
## 0.9555136 0.9555136 0.3877723 0.3791901 0.9558978 0.4117206

Hasil prediksi dari regresi logistik berupa nilai peluang, semakin besar nilai peluangnya maka semakin besar kemungkinan hasil prediksinya adalah kelas positifnya atau dalam hal ini adalah terjadi stillbirth. Kemudian kita akan satukan hasil prediksi ini dengan kolom \(Y\) pada data testing dalam bentuk data.frame. Kita definsikan kolom truth adalah nilai aktual dari data testing variabel \(Y\) data yang kita miliki.

pred_reglog <- data.frame(estimate_prob=pred_reglog0,truth=data_test$Y)
pred_reglog

untuk melakukan evaluasi performa prediksi, maka nilai peluang dari hasil prediksi regresi logistik harus kita konversi menjadi nilai kelas asal yaitu “1” dan “0”. Caranya adalah dengan menentukan suatu threshold tertentu pada nilai peluangnya. Sehingga threshold itulah yang menentukan hasil prediksi mana yang menjadi kelas “1” atau “0”. Berikut adalah sintaksnya.

Umumnya nilai threshold yang digunakan adalah \(0.5\), nilai threshold ini bisa kita ganti berdasarkan subjektivitas pengguna ataupun dengan metode-motode optimalisasi threshold.

threshold = 0.5
pred_reglog <- pred_reglog %>% 
  mutate(estimate=if_else(estimate_prob>threshold,"1","0")) %>% #kategorisasi
  mutate(estimate=factor(estimate,levels=c("1","0")))   %>%  # ubah ke factor
  # karena positive class dataset adalah 0 kita perlu rubah dulu ke 1
  mutate(truth = factor(truth,levels=c("1","0"))) # ubah kelas positif menjadi 1
pred_reglog

5.4 Confusion Matrix

Hasil prediksi ini biasanya kita evaluasi dengan suatu nilai metrik tertentu. Sebelum ke metrics tersebut, kita akan mengenal confusion matrix yang menjadi dasar perhitungan metrik-metrik tersebut. Berikut adalah komponen dari confusion matrix.

Berikut sintaks untuk menghitung CM.

confussion_matrix <- pred_reglog %>% conf_mat(truth=truth,estimate=estimate)

autoplot(confussion_matrix,type = "heatmap")+
  scale_fill_viridis_c(direction = -1,option = "mako",alpha = 0.4)
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

Dalam Confusion Matriks, terdapat 4 komponen, yakni sebagai berikut.

  1. TRUE POSITIVE (TP): Banyaknya amatan kelas positif yang hasil prediksi modelnya juga kelas positif

  2. TRUE NEGATIVE (TN): Banyaknya amatan kelas negatif yang hasil prediksi modelnya juga kelas negatif

  3. FALSE NEGATIVE (FN): Banyaknya amatan kelas negatif yang hasil prediksi modelnya kelas positif

  4. FALSE POSITIVE (FP): Banyaknya amatan kelas positif yang hasil prediksi modelnya kelas negatif.

Dari 4 komponen tersebut, dapat dihitung beberapa ukuran/ metrik evaluasi yakni

  1. Akurasi

  • Accuracy(Akurasi) adalah ukuran banyaknya hasil prediksi model yang benar dari keseluruhan hasil prediksi model.

  • Rumus

\[\text{Accuracy} = \frac{TP+TN}{TP+TN+FP+FN}\]

  1. Precision (disebut juga positive predictive value)

  • Precision (Presisi) adalah proporsi prediksi kelas positif yang benar dari total kelas positif yang diprediksi.

  • Rumus

\[\text{Precision} = \frac{TP}{TP+FP}\]

  1. Recall (juga disebut Sensitivity, Hit Rate, atau True Positive Rate)

  • Recall adalah proporsi kelas positif yang diprediksi benar dari total kelas positif yang sebenarnya.

  • Rumus

\[\text{Recall} = \frac{TP}{TP+FN}\]

  1. F1-Score

  • F1-Score adalah rata-rata harmonik dari presisi dan recall. Metric ini mempertimbangkan keseimbangan antara precision dan recall.

  • Rumus

\[\text{F1} = 2 \times \frac{\text{Precision}\times \text{ Recall} } {\text{Precision} + \text{ Recall}}\]

  1. Specificity (juga disebut True Negative Rate)

  • Specificity mengukur proporsi kelas negatif yang diprediksi benar dari total kelas negatif yang sebenarnya.

  • Rumus

\[\text{Specificity} = \frac{TN}{TN+FP}\]

  1. False Positive Rate (FPR)

  • FPR mengukur proporsi FALSE POSITIVE dari total kelas negatif sebenarnya.

  • Rumus

\[\text{FPR} = 1 - \text{Specificity} = \frac{FP}{TN+FP}\]

  1. False Negative Rate (FNR)

  • FNR mengukur proporsi FALSE NEGATIVE dari total kelas positif sebenarnya.

  • Rumus

\[\text{FNR} = \frac{FN}{TP+FN}\]

  1. Koefisien Korelasi Matthews (MCC)

  • MCC adalah ukuran kualitas klasifikasi biner, mempertimbangkan keempat elemen confusion matrix.

  • Rumus

\[\text{MCC} = \frac{TP \times TN - FP \times FN}{\sqrt{(TP+FP)(TP+FN)(TN+FP)(TN+FN)}}\]

  1. Balanced Accuracy

  • Balanced Accuracy adalah metrik evaluasi yang berguna ketika dataset memiliki ketimpangan (imbalance) di antara kelas-kelas yang berbeda. Metrik ini memberikan akurasi yang seimbang antara kelas mayoritas dan minoritas.

  • Rumus

\[\text{Balanced Accuracy} = \frac{\text{Specificity} + \text{Sensitivity}}{2}\]

Kesemua metriks tersebut dapat dikeluarkan menggunakan sintaks berikut.

confussion_matrix %>% summary()

atau bisa juga ditampilkan salah satu saja seperti berikut.

pred_reglog %>% accuracy(estimate = estimate,truth = truth)

atau bisa juga lebih dari satu

multi_metric <- metric_set(accuracy,sensitivity,specificity,bal_accuracy,mcc,f_meas)
pred_reglog %>% 
  multi_metric(estimate = estimate,truth = truth)

5.5 Kurva ROC

Evaluasi performa prediksi juga bisa didapatkan dengan menggunakan nilai peluang hasil prediksi yang didapatkan dari regresi logistik. Nilai peluang tersebut pada langkah sebelumnya disimpan dalam kolom estimate_prob.

pred_reglog
kurva_roc <- pred_reglog %>% roc_curve(truth = truth,estimate_prob)
kurva_roc 
autoplot(kurva_roc)

  • Berdasarkan dua output diatas, kita bisa melihat bahwa kurva ROC dibentuk berdasarkan metrik sensitivity dan specificity yang dihitung dari berbagai macam kemungkinan threshold. Oleh karena itu, kurva ROC ini juga sering digunakan untuk mendapatkan threshold optimal berdasarkan sensitivity dan specificity.

  • Kurva ROC ini biasanya digunakan untuk menentukan model dengan performa prediksi terbaik. Namun, pada ilustrasi ini kita hanya menggunakan satu model sehingga kurva ROC ini kurang bisa digunakan. Tetapi, terdapat metrik yang berbasiskan kurva roc yang bisa digunakan tanpa membandingkan dengan model lain yaitu AREA Under Curve (AUC).

pred_reglog %>% 
  roc_auc(truth = truth,estimate_prob)

5.6 Treshold Optimal

Penentuan Treshold Optimal sifatnya adalah optional, umumnya Treshold yang digunakan sebesar 0.5. Pada kali ini kita akan mencoba mendapatkan threshold optimal berdasarkan sensitivity dan specificity menggunakan Youden’s J statistic yang didefinisikan sebagai berikut.

\[\text{Youden's J Statistics = sensitivity + specificity - 1}\] Semakin mendekati 1 nilainya maka semakin baik thresholdnya dan semakin nilainya mendekati 0 semakin buruk thresholdnya. Jika kita terapkan pada sintaks adalah sebagai berikut.

kurva_roc %>% 
  mutate(youden_index=sensitivity+specificity-1)

kemudian karena kita mencari Youden’s J statistic yang mendekati satu berarti sama saja kita mencari nilai threshold yang Youden’s J statistic paling maksimum. Kita bisa mendapatkan nilai maksimum ini dengan bantuan fungsi slice_max

optimal_thresh <- kurva_roc %>% 
  mutate(youden_index=sensitivity+specificity-1) %>% 
  slice_max(order_by = youden_index)
optimal_thresh

Kemudian hasil prediksi dengan threshold optimal bisa didapatkan dengan sintaks berikut

pred_reglog2 <- data.frame(estimate_prob=pred_reglog0,truth=data_test$Y)
threshold2 = optimal_thresh %>% pull(.threshold)
pred_reglog2 <- pred_reglog2 %>% 
  mutate(estimate=if_else(estimate_prob>threshold2,"1","0")) %>% 
  mutate(estimate=factor(estimate,levels=c("1","0"))) %>% 
  # karena positive class dataset adalah 0 kita perlu rubah dulu ke 1
  mutate(truth = factor(truth,levels=c("1","0")))
pred_reglog2

Kemudian jika kita bandingkan antara threshold 0.5 dengan threshold optimal pada beberapa metric maka hasilnya adalah sebagai berikut.

res_metric2 <- pred_reglog2 %>% 
  multi_metric(truth = truth,estimate = estimate) %>% 
  mutate(threshold=threshold2) %>% 
  relocate(threshold,.after = .metric)

res_metric <- pred_reglog %>% 
  multi_metric(truth = truth,estimate = estimate) %>% 
  mutate(threshold=threshold) %>% 
  relocate(threshold,.after = .metric)
res_metric2 %>% 
  bind_rows(res_metric) %>% 
  select(-.estimator) %>% 
  pivot_wider(id_cols = .metric,names_from = threshold,
              names_prefix ="thres=", 
              values_from = .estimate)

6 Prediksi Data Baru

untuk mengilustrasikan prediksi pada data baru, kita akan membuat suatu data baru “buatan” dengan mengambil sampel acak dari data yang kita miliki. Data baru bisanya tidak memuat variabel response sehingga kita akan hilangkan kolom \(Y\).

set.seed(2025)
data_baru <- data %>% slice_sample(n = 10) 
data_baru
set.seed(2025)
data_baru <- data %>% slice_sample(n = 10) %>% select(-Y)
data_baru
pred_data_baru <- predict(reglog,newdata = data_baru,type = "response")
pred_data_baru
##          1          2          3          4          5          6          7 
## 0.38132885 0.07379667 0.59237109 0.61193519 0.59237109 0.40952425 0.60977789 
##          8          9         10 
## 0.38992906 0.41392044 0.39642476
pred_data_baru_df <- data.frame(estimate=pred_data_baru) %>% 
                      mutate(estimate=if_else(estimate>threshold2,"1","0")) %>% 
                      mutate(estimate=factor(estimate,levels=c("0","1")))
pred_data_baru_df

Alhamdulillah, selesai

7 Daftar Pustaka

  1. Agresti (link).

  2. Geery Dito (link ).

  3. Analisis Regresi Logistik Biner Untuk Menentukan Faktor Stillbirth Di Kabupaten Aceh Timur (link).

  4. Analisis Data dengan Regresi Logistik di R Studio ( link).

  5. Jimy Badrus (link ).