Regresi logistik di R: glm() dengan family = binomial

Saat hasilnya ya/tidak

Regresi linear memprediksi angka. Tapi banyak pertanyaan yang layak dimodelkan bersifat biner: apakah pelanggan berhenti berlangganan, apakah pasien sembuh, apakah email-nya diklik. Memasang garis lurus pada hasil 0/1 langsung rusak: garis itu dengan senang hati memprediksi probabilitas −0,3 atau 1,4, yang tidak masuk akal.

Regresi logistik memperbaikinya dengan memodelkan probabilitas hasilnya melalui transformasi log-odds (logit): log(p / (1 − p)) = intercept + slope × x. Skala log-odds membentang di seluruh garis bilangan, jadi persamaan linear cocok secara alami di sana, dan memetakannya kembali menekan setiap prediksi ke dalam (0, 1) di sepanjang kurva berbentuk S yang sudah dikenal. Harga dari trik ini: koefisiennya berada pada skala log-odds, dan seluruh permainan membaca regresi logistik adalah menerjemahkannya kembali menjadi sesuatu yang dipahami manusia.

Membuat model: glm() dengan family = binomial

glm() (generalized linear model) adalah kakak dari lm(); family = binomial memilih regresi logistik. Di mtcars, am mencatat jenis transmisi (1 = manual, 0 = otomatis). Apakah mobil yang irit bahan bakar cenderung bertransmisi manual?

fit <- glm(am ~ mpg, data = mtcars, family = binomial)
summary(fit)

Dua hal sebelum membaca output-nya. Pertama, family = binomial tidak opsional: kalau dihilangkan, glm() diam-diam memasang kuadrat terkecil biasa. Kedua, hasilnya harus biner: 0/1, logical, atau factor dengan dua level (R memodelkan probabilitas level kedua).

Sekarang summary-nya, blok demi blok:

• Coefficients: Estimate untuk mpg sekitar 0,31, dan itu adalah kemiringan log-odds: setiap tambahan satu mpg menambah 0,31 pada log-odds bertransmisi manual. Positif berarti "menaikkan probabilitas", negatif berarti "menurunkan"; di luar tandanya, intuisi siapa pun tidak bekerja pada skala ini, itulah sebabnya bagian berikutnya ada.
• z value dan Pr(>|z|): logikanya sama dengan uji t pada regresi (Estimate ÷ Std. Error, lalu p-value untuk "mungkinkah ini nol?"), hanya saja memakai pendekatan normal, sehingga z bukan t. Di sini p ≈ 0,011: hubungan antara mpg dan jenis transmisi kemungkinan besar bukan derau.
• Null deviance vs Residual deviance: deviance adalah ukuran buruknya kecocokan di dunia glm (makin kecil makin baik). Null deviance (43,2 pada 31 df) adalah model yang hanya berisi intercept; residual deviance (29,7 pada 30 df) adalah modelmu. Penurunan sekitar 13,6 untuk satu df prediktor adalah padanan glm dari "R-squared naik."
• AIC: skor pembanding model yang menyeimbangkan kecocokan dengan kerumitan; yang lebih rendah menang. Tidak bermakna sendirian, tapi berguna untuk membandingkan model-model kandidat pada data yang sama.

Dari log-odds ke odds ratio: exp(coef())

Memangkatkan dengan exp memindahkan koefisien dari skala log-odds yang aditif ke skala odds yang multiplikatif:

fit <- glm(am ~ mpg, data = mtcars, family = binomial)

print(exp(coef(fit)))
print(exp(confint(fit)))

exp(0.307) ≈ 1.36, dan inilah pola kalimat jujur yang perlu dihafal: "setiap tambahan satu mpg mengalikan odds transmisi manual sekitar 1,36." Odds ratio di atas 1 menaikkan odds, di bawah 1 menurunkannya, tepat 1 berarti tidak ada efek. Itulah sebabnya keputusan interval kepercayaan untuk odds ratio adalah "apakah intervalnya tidak memuat 1?" (bukan nol; nol adalah batasnya di skala log-odds).

Perhatikan bahasanya: odds bukan probabilitas. Odds = p / (1 − p), jadi probabilitas 0,75 sama dengan odds 3. Mengalikan odds dengan 1,36 tidak sama dengan mengalikan probabilitas dengan 1,36, dan saat hasilnya umum terjadi, selisihnya besar. Odds ratio sebesar 2 untuk hasil yang langka berperilaku seperti "kira-kira dua kali lipat risikonya"; untuk hasil yang terjadi 50% dari waktu, jelas tidak begitu. Jangan pernah melaporkan odds ratio dengan kata-kata risk ratio ("1,36 kali lebih mungkin") kecuali hasilnya langka.

Probabilitas prediksi: jebakan type = "response"

Bug regresi logistik yang paling umum di dunia nyata:

fit <- glm(am ~ mpg, data = mtcars, family = binomial)
new_cars <- data.frame(mpg = c(15, 20, 25, 30))

print(predict(fit, newdata = new_cars))                      # log-odds!
print(predict(fit, newdata = new_cars, type = "response"))   # probabilities

Pemanggilan pertama mengembalikan default type = "link": prediksi pada skala log-odds, lengkap dengan nilai negatif. Pemanggilan kedua mengembalikan probabilitas yang sebenarnya. Kalau "probabilitas"-mu pernah keluar negatif atau di atas 1, inilah alasannya. Jalankan blok itu: mobil 15 mpg hampir tidak punya peluang bertransmisi manual, mobil 30 mpg sangat mungkin manual, dan kurva S melengkung di antaranya.

Klasifikasi: ambang batas dan tabel confusion

Probabilitas menjadi kelas prediksi dengan memilih ambang batas (0,5 adalah pilihan default-nya), dan kartu skor yang jujur adalah tabel prediksi vs aktual:

fit <- glm(am ~ mpg, data = mtcars, family = binomial)

probs <- predict(fit, type = "response")   # fitted probabilities, one per car
predicted <- probs > 0.5

print(table(actual = mtcars$am, predicted = predicted))

Sel-sel diagonal adalah tebakan yang benar; dua sel di luar diagonal adalah dua kesalahan yang berbeda (memprediksi manual untuk mobil otomatis, dan sebaliknya). Akurasi keseluruhan saja bisa membuat model terlihat jauh lebih baik: kalau 95% pelanggan tidak berhenti berlangganan, "prediksi tidak ada yang berhenti" mendapat skor 95% sambil tidak menangkap satu pun yang berhenti. Jadi selalu lihat kedua jenis error. Dan 0,5 adalah konvensi, bukan hukum: saat dua kesalahan itu punya biaya yang berbeda, geser ambang batasnya.

Satu catatan kejujuran: tabel ini menilai model pada data yang sama dengan yang dipakai untuk membuatnya, sehingga terlihat lebih baik. Evaluasi yang sebenarnya menyisihkan data yang tidak pernah dilihat model.

Banyak prediktor

Persis seperti lm(): tambahkan suku dengan +, dan setiap penafsiran mendapat keterangan tambahan "dengan yang lain dijaga tetap":

fit2 <- glm(am ~ mpg + qsec, data = mtcars, family = binomial)
summary(fit2)
print(exp(coef(fit2)))

Setiap koefisien yang dipangkatkan dengan exp sekarang adalah pengali odds untuk kenaikan satu unit pada prediktor itu di antara mobil-mobil yang sama pada prediktor lainnya. Mesinnya bisa dikembangkan, tapi begitu pula peringatan dari regresi linear: prediktor yang saling berkorelasi mengacak koefisien satu sama lain.

Peringatan

• Separasi sempurna. Kalau sebuah prediktor membagi hasilnya dengan sempurna (setiap mobil di atas mpg tertentu manual, setiap mobil di bawahnya otomatis), koefisien maximum likelihood ingin bernilai tak hingga. R memberi peringatan, glm.fit: fitted probabilities numerically 0 or 1 occurred, dan melaporkan koefisien yang sangat besar dengan standard error yang tidak masuk akal. Jangan memakai angka-angka itu; sederhanakan modelnya, kumpulkan lebih banyak data, atau pakai metode berpenalti (package brglm2 atau logistf).
• Kejadian yang cukup. Batasan yang menentukan adalah jumlah hasil yang lebih langka, bukan total baris. Aturan praktis lama menginginkan sekitar 10-15 kejadian per prediktor; contoh 32 mobil di sini untuk mengajarkan mekanismenya, bukan templat untuk ukuran sampel yang layak dipublikasikan.
• Odds ratio bukan risk ratio saat hasilnya umum terjadi. Sudah dibahas di atas, diulang karena penelaah akan menangkapnya meskipun kamu tidak.

Yang kamu bawa pulang

• Hasil biner → glm(y ~ x, data = df, family = binomial); jangan pernah lupa family-nya.
• Koefisien mentah adalah log-odds; exp(coef(fit)) memberikan odds ratio, dan nilai nol untuk intervalnya adalah 1.
• Pola kalimatnya: "setiap kenaikan satu unit pada x mengalikan odds hasilnya dengan exp(b)."
• predict(..., type = "response") untuk probabilitas; default-nya mengembalikan log-odds, kebingungan nomor satu.
• Klasifikasikan dengan ambang batas dan nilai dengan tabel confusion; akurasi saja bisa menipu.
• Waspadai peringatan separasi, hitung jumlah kejadianmu, dan jangan mendandani odds ratio sebagai risk ratio.

Berikutnya: mesin di balik setiap interval yang sudah kamu lihat sejauh ini, yaitu interval kepercayaan dengan t.test(), confint(), dan prop.test().