Membuat Data Simulasi Model Regresi Logistik

Data simulasi adalah data yang dihasilkan secara artifisial menggunakan metode statistik atau algoritma tertentu, yang meniru karakteristik dari data dunia nyata. Proses simulasi data melibatkan pembuatan sampel data dengan karakteristik yang diinginkan, seperti distribusi tertentu, korelasi antar variabel, serta pola yang relevan dengan studi yang dilakukan.

Terdapat banyak alasan mengapa peneliti memerlukan data simulasi, beberapa alasan yang mendasarinya antara lain sebagai berikut:

• Pembelajaran: Data simulasi adalah alat yang penting untuk tujuan pembelajaran dan pendidikan. Ini memungkinkan peneliti untuk bereksperimen dengan berbagai karakteristik data dan model tanpa harus mengumpulkan data dunia nyata yang mungkin sulit atau mahal untuk diperoleh.
• Pengujian Metode Statistik: Data simulasi memungkinkan para peneliti untuk menguji dan mengevaluasi metode statistik atau algoritma machine learning dalam kondisi yang terkontrol. Hal ini penting untuk memastikan bahwa metode tersebut bekerja dengan baik sebelum diterapkan pada data dunia nyata.
• Pemahaman Teoretis: Dengan menggunakan data simulasi, kita dapat memahami perilaku teoretis dari model statistik. Ini membantu dalam memahami bagaimana model berfungsi, bagaimana variabel mempengaruhi hasil, dan bagaimana model bereaksi terhadap perubahan data.
• Validasi dan Verifikasi: Simulasi data memungkinkan validasi dan verifikasi model, memastikan bahwa model memberikan hasil yang konsisten dan akurat. Ini juga membantu dalam mendeteksi bias atau kesalahan dalam model.

Secara garis besar terdapat 5 proses utama yang dilakukan untuk membangkitkan data simulasi pada pemodelan regresi logistik.

• Menentukan kriteria data: penentuan ini merupakan proses pengaturan karakteristik data yang akan dibangkitkan. Hal ini mencakup banyaknya data dan jumlah peubah prediktor, bagaimana nilai tengahnya, bagaimana keragamannya dan bagaimana keterkaitannya dengan peubah prediktor lainnya.
• Membangkitkan peubah prediktor: setelah skenario ditetapkan, maka kita bangkitkan peubah-peubah prediktor sesuai skenario tersebut.
• Mendefinisikan koefisien model regresi logistik: menentukan koefisien yang akan digunakan dalam model regresi logistik, termasuk penambahan noise.
• Menghitung nilai respon linear (log odds): berdasarkan koefisien yang didefinisikan sebelumnya, kita dapat menghitung nilai peubah respon dari model linear.
• Menghitung nilai peluang: dari nilai log odds, kita dapat menghitung nilai peluangnya menggunakan fungsi logistik.
• Menentukan Kelas respon: penetapan kate dilakukan dengan mengatur titik cut-off tertentu pada nilai probabilitas (misalkan 0,5 di mana data dengan nilai probabilitas kurang dari 0,5 ditetapkan sebagai kelas A dan selebihnya sebagai kelas B).

Tahapan-tahapan di atas akan dijabarkan lebih rinci pada bagian berikutnya beserta implementasinya menggunakan bahasa R.

Misalkan kita ingin membangkitkan sebanyak 100 data dengan 5 peubah prediktor (X1, X2, X3, X4 dan X5) berdasarkan skenario sebagai berikut:

Peubah $\text{X} \sim \mathcal{N}(\boldsymbol{\mu}, \Sigma)$ (Multivariate Normal) dengan vektor nilai tengah $\boldsymbol{\mu}$ dan matriks kovarian $\Sigma$ adalah:

Pada skenario di atas, kita mengatur X1 dan X2 memiliki hubungan linier yang rendah dengan peubah lainnya (0,05, 0,10, 0,15). Sementara itu untuk X3, X4 dan X5 diatur sedemikian rupa sehingga memiliki korelasi yang cukup tinggi (0,65, 0,80 dan 0,90).

Penggunaan sebaran multivariat normal dengan $\boldsymbol{\mu}=0$ dan $\Sigma=1$, akan mempermudah kita mengatur nilai koefisien korelasi antar peubah prediktor yang diharapkan. Pada data yang menyebar dengan nilai tengah 0 dan varians 1, maka nilai-nilai koefisien korelasi akan sama dengan nilai varians-kovarian nya. Dengan demikian, saat implementasi menggunakan R nantinya, kita tidak rumit menghitung nilai matriks kovarian karena cukup berdasarkan nilai koefisien korelasi yang diinginkan.

Meskipun demikian, pembangkitan ini tidak terbatas untuk sebaran seperti itu saja. Setelah memperoleh data simulasi awal, jika diinginkan kita dapat melakukan berbagai operasi matematika pada hasil simulasi untuk memperoleh data simulasi baru dengan nilai rata-rata dan varians yang berbeda.

Misalnya untuk simulasi data di mana terdapat perbedaan satuan dan keragaman data, maka kita dapat merubah skala data serta variansnya menjadi lebih kecil dan lebih besar, dengan tetap mempertahankan nilai koefisien korelasinya.

Pada contoh ini, kita akan merubah nilai rata-rata dan keragaman data sehingga memiliki kriteria sebagai berikut:

• X1 akan dibuat memiliki rata-rata 10
• X2 dengan rata-rata 20 dan varians 16
• X3 dengan rata-rata 1, varians 0.04 dan koefisen korelasi menjadi negatif.
• X4 dan X5 tetap

Pembangkitan data simulasi yang menyebar multivariat normal dapat dilakukan menggunakan fungsi mvnorm dari paket MASS. Tiga parameter utama yang diperlukan yaitu:

• n: jumlah sampel data yang akan dibangkitkan
• mu: vektor yang berisi nilai rata-rata sejumlah variabel yang akan dibangkitkan
• sigma: matriks kovarian

Secara default pembangkitan menggunakan mvnorm tidak akan menghasilkan nilai matriks kovarian yang sama persis, melainkan hanya mendekati saja. Jika ingin tepat sama dengan skenario, maka kita dapat mengatur parameter empirical=T.

Berikut implementasinya dalam bahasa R:

library(MASS)

# mengatur seed untuk reproduksi
set.seed(123)

# Vektor mean
X_mu <- rep(0, 5)

# matriks var-cov
X_sig <- matrix(c(1.00, 0.10, 0.10, 0.15, 0.05,
                  0.10, 1.00, 0.15, 0.05, 0.20,
                  0.10, 0.15, 1.00, 0.80, 0.65,
                  0.15, 0.05, 0.80, 1.00, 0.90,
                  0.05, 0.20, 0.65, 0.90, 1.00), 
                  nrow = 5, byrow = TRUE) 

# membangkitkan data
X <- mvrnorm(100, mu = X_mu, Sigma = X_sig,  empirical = T)

head(X, 5)

cat("\n\nNilai rata-rata peubah X:\n")
apply(X, 2, function(x){round(mean(x))})

cat("\n\nMatriks Var-Cov peubah X:\n")
round(cov(X), 2)

cat("\n\nMatriks Korelasi peubah X:\n")
round(cor(X), 2)

# OUTPUT

           [,1]       [,2]        [,3]       [,4]       [,5]
[1,]  1.2392191 -0.9211449 -1.49929229 -0.4731595 -0.3383906
[2,]  1.6702834 -0.1974900  0.02616934  0.7128291  0.4707361
[3,]  0.7245258  1.6863432  0.42332041  0.1950561  0.5670124
[4,]  0.8224963  0.2697408 -0.56225192 -0.4553403 -0.4489560
[5,] -0.5334866  0.5064252 -0.72456483 -0.8823165 -0.7326693

Nilai rata-rata peubah X:
[1] 0 0 0 0 0

Matriks Var-Cov peubah X:
     [,1] [,2] [,3] [,4] [,5]
[1,] 1.00 0.10 0.10 0.15 0.05
[2,] 0.10 1.00 0.15 0.05 0.20
[3,] 0.10 0.15 1.00 0.80 0.65
[4,] 0.15 0.05 0.80 1.00 0.90
[5,] 0.05 0.20 0.65 0.90 1.00

Matriks Korelasi peubah X:
     [,1] [,2] [,3] [,4] [,5]
[1,] 1.00 0.10 0.10 0.15 0.05
[2,] 0.10 1.00 0.15 0.05 0.20
[3,] 0.10 0.15 1.00 0.80 0.65
[4,] 0.15 0.05 0.80 1.00 0.90
[5,] 0.05 0.20 0.65 0.90 1.00

Langkah ini bersifat opsional. Jika diperlukan, kita dapat merubah nilai tengah serta varians dari data yang sudah dibangkitkan dan tetap menjaga nilai koefisien korelasinya. Seluruh operasi yang bersifat linier tidak akan mempengaruhi nilai koefisien korelasi antar peubah tersebut (kecuali tanda + atau – saja).

Sesuai skenario yang sudah ditentukan maka kita dapat membuat X1, X2 dan X3 yang baru dengan kode berikut ini:

# membuat copy variabel X ke X_new
X_new <- X

# transformasi data
X_new[,1] <- X[,1] + 10      # mean = 10, vars tetap
X_new[,2] <- X[,2]*4 + 20    # mean = 20, vars (4^2) = 16
X_new[,3] <- 1 - X[,3]/5     # mean = 1, vars (1/5)^2 = 0.04, arah korelasi negatif

# menampilkan data X_new
head(X_new, 10)

# menampilkan nilai rata-rata
cat("\n\nNilai rata-rata peubah X_new:\n")
apply(X_new, 2, function(x){round(mean(x))})

# menampilkan matriks kovarian
cat("\n\nMatriks Var-Cov peubah X_new:\n")
round(cov(X_new), 2)

# menampilkan matriks korelasi
cat("\n\nMatriks Korelasi peubah X_new:\n")
round(cor(X_new), 2)

# OUTPUT

           [,1]     [,2]      [,3]       [,4]       [,5]
 [1,] 11.239219 16.31542 1.2998585 -0.4731595 -0.3383906
 [2,] 11.670283 19.21004 0.9947661  0.7128291  0.4707361
 [3,] 10.724526 26.74537 0.9153359  0.1950561  0.5670124
 [4,] 10.822496 21.07896 1.1124504 -0.4553403 -0.4489560
 [5,]  9.466513 22.02570 1.1449130 -0.8823165 -0.7326693
 [6,] 10.306783 26.83223 0.7181908  1.6249666  1.9576109
 [7,] 10.809229 27.05050 1.0839079 -0.7415477 -0.8355571
 [8,]  8.601573 18.93217 1.2724415 -0.9705854 -1.1817385
 [9,] 10.028374 14.49141 1.0763083  0.8879220  0.9207116
[10,]  8.878925 14.60134 0.7596142  1.1781028  0.9127282

Nilai rata-rata peubah X_new:
[1] 10 20  1  0  0

Matriks Var-Cov peubah X_new:
      [,1]  [,2]  [,3]  [,4]  [,5]
[1,]  1.00  0.40 -0.02  0.15  0.05
[2,]  0.40 16.00 -0.12  0.20  0.80
[3,] -0.02 -0.12  0.04 -0.16 -0.13
[4,]  0.15  0.20 -0.16  1.00  0.90
[5,]  0.05  0.80 -0.13  0.90  1.00

Matriks Korelasi peubah X_new:
      [,1]  [,2]  [,3]  [,4]  [,5]
[1,]  1.00  0.10 -0.10  0.15  0.05
[2,]  0.10  1.00 -0.15  0.05  0.20
[3,] -0.10 -0.15  1.00 -0.80 -0.65
[4,]  0.15  0.05 -0.80  1.00  0.90
[5,]  0.05  0.20 -0.65  0.90  1.00

Berdasarkan hasil di atas, data simulasi yang baru memiliki karakteristik yang lebih beragam. Namun hubungan antar variabel tetap sama seperti sebelumnya. Hal ini dapat dilihat dari matriks korelasi yang dihasilkan tetap sama (kecuali tanda negatif untuk peubah X3).

Pada bagian ini kita akan menentukan koefisien model untuk membangkitkan nilai peubah respon bagi model regresi logistik. Penentuannya juga tentunya bersifat subjektif dan tergatung kebutuhan peneliti, menentukan peubah mana yang akan diberikan pengaruh besar, kecil, positif atau negatif.

Di sini, kita akan menggunakan pengaturan sebagai berikut:

Sehingga diperoleh persamaan linier dari model yaitu:

Setelah koefisien $\beta$ ditentukan, maka kita dapat menghitung nilai log odds untuk setiap baris data yang sudah dibangkitkan. Pada contoh di bawah ini, kita membuat fungsi untuk menghitung log odds berdasarkan data dan parameter $\beta$ yang diberikan. Perhitungan kita buat dalam bentuk fungsi agar dapat digunakan berulang kali dengan pengaturan parameter yang berbeda.

set.seed(42)

# fungsi menghitung log-odds
calc_log_odds <- function(dt, pars, noise) {
  
  const <- matrix(rep(1, nrow(dt)), ncol = 1)     # membuat kolom untuk konstanta b0
  dt_full <- cbind(const, dt)                     # menggabunngkan nilai konstanta ke data
  log_odds <- dt_full %*% pars + noise            # menghitung nilai log-odds semua baris
  
  return (log_odds)
}

# menentukan beta_0, beta_1, ...beta 5
params <- c(0.1, -1.0, 0.5, 0.25, -0.1, 0.01)

# membangkitkan data noise
noise <- rnorm(nrow(X_new), 0, 0.25)

# menghitung log odds data X_new berdasarkan parameter params
result_log_odds <- calc_log_odds(X_new, params, noise)

cat("Nilai Log Odds:\n")
head(result_log_odds, 10)

# OUTPUT

Nilai Log Odds:
            [,1]
 [1,] -2.2698727
 [2,] -1.9243220
 [3,]  3.0539411
 [4,]  0.2943579
 [5,]  2.1145374
 [6,]  3.2194270
 [7,]  3.5306797
 [8,]  1.3442005
 [9,] -1.9885713
[10,] -1.4127121

Nilai peluang dapat dihitung dari nilai log odds yang dihasilkan pada langkah sebelumnya melalui fungsi logistik. Adapun persamaannya adalah sebagai berikut:

Di bawah ini merupakan implementasi persamaan di atas menggunakan fungsi R. Selanjutnya fungsi tersebut kita gunakan untuk menghitung nilai peluang dari log odds yang telah diperoleh sebelumnya.

# fungsi menghitung nilai peluang berdasarkan log-odds
get_probs <- function(log_odds) {
  return (1/(1 + exp(-log_odds)))
}

result_probs <- get_probs(result_log_odds)

cat("Nilai peluang:\n")
head(result_probs, 10)

# OUTPUT

Nilai peluang:
            [,1]
 [1,] 0.09364901
 [2,] 0.12738038
 [3,] 0.95495237
 [4,] 0.57306269
 [5,] 0.89230813
 [6,] 0.96155884
 [7,] 0.97154821
 [8,] 0.79317987
 [9,] 0.12040810
[10,] 0.19580664

Langkah terakhir yang perlu dilakukan yaitu menentukan kelas peubah respon untuk masing-masing data. Penentuan kelas didasari dengan nilai cut-off. Nilai yang umum digunakan adalah 0.5. Misalkan, menggunakan cut-off tersebut, maka data dengan nilai peluang lebih dari 0,5 akan dikategorikan sebagai kelas A dan data dengan peluang kurang dari itu dimasukkan ke dalam kelas B. Tentunya kita dapat menggunakan nilai cut-off lainnya sesuai kebutuhan.

Tidak hanya terbatas untuk 2 kelas, kita dapat juga menentukan batas cut-off sebanyak 2 titik sehingga kita dapat menghasilkan peubah respon dengan 3 kelas.

Pada contoh di bawah ini, kita akan mencoba menerapkan beberapa cut-off yang berbeda pada data yang kita miliki.

# contoh 1 (cut-off 0.5)
output1 <- sapply(result_probs, function(x) ifelse(x>0.5, "Kat. A", "Kat. B"))

# contoh 2 (cut-off 0.7)
output2 <- sapply(result_probs, function(x) ifelse(x>0.7, "Kat. A", "Kat. B"))


cutoff_fun <- function(dt, cutoff1, cutoff2) {
  # asumsikan cutoff1 selalu lebih kecil dari cutoff2
  if(dt < cutoff1) "Kat. A"
  else if (dt < cutoff2) "Kat. B"
  else "Kat. C"
}

# contoh 3 (cut-off 0,33 dan 0,67)
output3 <- sapply(result_probs, function(x) cutoff_fun(x, 0.33, 0.67))


# menggabungkan data X dan output menjadi dataframe
data_complete <- data.frame(X_new, Y1=output1, Y2=output2, Y3=output3)

print(data_complete)

# OUTPUT

            X1        X2         X3          X4          X5      Y1      Y2      Y3
  1   11.23922  16.31542  1.2998585  -0.4731595  -0.3383906  Kat. B  Kat. B  Kat. A
  2   11.67028  19.21004  0.9947661   0.7128291   0.4707361  Kat. B  Kat. B  Kat. A
  3   10.72453  26.74537  0.9153359   0.1950561   0.5670124  Kat. A  Kat. A  Kat. C
  4    0.82250  21.07896  1.1124504  -0.4553403  -0.4489560  Kat. A  Kat. B  Kat. B
  .        ...       ...        ...         ...         ...     ...     ...     ...
 98    8.36084  26.41011  0.9457376   0.0855687   0.7503880  Kat. A  Kat. A  Kat. C
 99   10.43854  21.49085  1.1660664  -0.9405982  -1.0522992  Kat. A  Kat. B  Kat. C
100   10.10606  16.31728  1.3464233  -0.9455142  -0.9359788  Kat. B  Kat. B  Kat. A

Selanjutnya data simulasi yang sudah dihasilkan dapat disimpan ke dalam file dan digunakan lebih lanjut untuk melakukan analisis.

note: untuk pembangkitan data simulasi dengan lebih dari 2 kelas, terdapat cara lain yang lebih ‘elegan’ yaitu mengikuti kerangka kerja model regresi logistik multinomial atau regresi logistik ordinal. Misalkan pada regresi logistik ordinal dengan 3 kelas respon, maka kita dapat menentukan cut-off (berupa 2 nilai intecept yang berbeda) untuk untuk setiap kelasnya.

• mvrnorm function – RDocumentation

• Pembangkitan Bilangan Acak pada Bahasa R
• Regresi Logistik untuk Klasifikasi Biner