Laporan ini menganalisis berbagai aspek prediktif perilaku pelanggan, kinerja akademik, perputaran karyawan, dan kunjungan layanan kesehatan menggunakan empat dataset sekunder yang berbeda. Setiap dataset diaplikasikan pada satu jenis model statistik yang spesifik karena struktur dan bentuk variabel responnya yang berbeda, mulai dari data biner, ordinal, multinomial, hingga data diskret (cacah).
Empat model yang digunakan adalah:
Churn = Yes) atau bertahan (Churn = No),
berdasarkan profil layanan dan biaya bulanan mereka.G3). Variabel respon ini dikelompokkan ke dalam skala
ordinal (Rendah < Sedang < Tinggi) untuk melihat pengaruh latar
belakang keluarga, konsumsi alkohol, dan waktu belajar terhadap
pencapaian siswa.JobRole) di
perusahaan. Variabel respon memiliki banyak kategori tanpa urutan,
sehingga setiap kategori peran akan dibandingkan terhadap satu kategori
acuan (baseline).Number of Doctors Visited). Model ini sangat cocok karena
datanya berupa bilangan bulat non-negatif (\(0, 1, 2, \dots\)) yang dipengaruhi oleh
faktor kesehatan fisik, mental, dan gangguan tidur pasien.| Istilah | Penjelasan dalam laporan ini |
|---|---|
| Churn / Retensi | Kondisi ketika pelanggan memutuskan untuk berhenti berlangganan layanan (Telco Customer Churn). Dalam data biner, 1 (Yes) berarti churn (berhenti) dan 0 (No) berarti loyal (bertahan). |
| Performa G3 | Nilai akhir siswa pada mata pelajaran bahasa Portugis (student-por.csv). Dalam model ordinal, nilai ini dikelompokkan menjadi tiga tingkatan: Rendah, Sedang, dan Tinggi. |
| JobRole | Peran jabatan karyawan dalam dataset HR Employee Attrition. Digunakan sebagai respon multinomial karena memiliki banyak kategori tanpa urutan. |
| Data Cacah (Count Data) | Data yang nilainya berupa bilangan bulat non-negatif hasil dari penghitungan. Dalam laporan ini, contohnya adalah jumlah dokter yang dikunjungi oleh pasien lansia. |
| Odds Ratio (OR) | Ukuran untuk membaca pengaruh variabel prediktor pada model logistik (biner dan ordinal). Untuk model ordinal dengan fungsi MASS::polr di R, nilai \(\text{OR} = \exp(\beta)\) diinterpretasikan sebagai perubahan peluang (odds) untuk berada pada kategori performa/risiko yang lebih tinggi, karena spesifikasi modelnya menggunakan formula:\[\text{logit}(P(Y \le j)) = \alpha_j - x\beta\] |
| Relative risk ratio / RRR | Ukuran untuk membaca pengaruh prediktor pada regresi multinomial. RRR membandingkan satu kategori respon terhadap kategori acuan. |
| Incidence rate ratio / IRR | Ukuran untuk membaca pengaruh prediktor pada regresi Poisson. IRR lebih dari 1 berarti rata-rata jumlah kejadian lebih tinggi. |
| Data leakage | Kondisi ketika model diberi variabel yang sebenarnya sudah mengandung jawaban. Variabel seperti ini harus dihindari agar model tidak terlihat bagus secara palsu. |
Dataset yang digunakan bersumber dari repositori publik/sekunder: IBM Business Analytics Community untuk data perputaran pelanggan telekomunikasi (Telco Customer Churn) dan perputaran karyawan (HR Employee Attrition), UCI Machine Learning Repository untuk data performa akademik siswa bahasa Portugis (student-por), serta National Poll on Healthy Aging (NPHA) untuk data kunjungan dokter pasien lansia.
| Model | Dataset | Variabel respon | Bentuk respon | Tujuan analisis |
|---|---|---|---|---|
| Logistik biner | Telco Customer Churn | Churn |
0/1 | Menjelaskan peluang pelanggan berhenti berlangganan |
| Ordinal | Student Performance (Portugis) | Performa_G3 |
Rendah < Sedang < Tinggi | Menjelaskan faktor penentu tingkat capaian akademik siswa |
| Multinomial | HR Employee Attrition | JobRole |
Banyak kategori tanpa urutan | Menjelaskan kecenderungan penempatan atau peran jabatan karyawan |
| Poisson | NPHA Doctor Visits | Number of Doctors Visited |
Hitungan/count | Menjelaskan frekuensi jumlah kunjungan pasien ke dokter |
Regresi logistik biner dipakai ketika variabel respon hanya memiliki
dua kategori (dikotomi). Dalam laporan ini, responnya adalah apakah
pelanggan berhenti berlangganan atau tidak (Churn = Yes
atau Churn = No).
\[ \log\left(\frac{P(Y=1)}{P(Y=0)}\right) = \beta_0 + \beta_1X_1 + \cdots + \beta_pX_p \]
Regresi ordinal dipakai ketika variabel respon memiliki tingkatan atau urutan yang jelas. Dalam laporan ini, responnya adalah performa nilai akhir siswa (G3) yang dikelompokkan menjadi kategori Rendah < Sedang < Tinggi.
Pada analisis ini model diestimasi menggunakan fungsi
MASS::polr. Spesifikasi model polr adalah:
\[ \log\left(\frac{P(Y \leq j)}{P(Y > j)}\right) = \alpha_j - x\beta \]
Artinya, tanda koefisien harus dibaca dengan hati-hati. Jika koefisien \(\beta\) bernilai positif, maka nilai \(\alpha_j - x\beta\) menjadi lebih kecil. Akibatnya, peluang kumulatif untuk berada pada kategori rendah \(P(Y \le j)\) menurun, sehingga peluang untuk berada pada kategori yang lebih tinggi meningkat. Dengan urutan tingkat nilai Rendah < Sedang < Tinggi, koefisien positif berarti variabel prediktor tersebut mendorong kecenderungan siswa untuk bergeser ke tingkat performa akademik yang lebih tinggi/baik.
Regresi multinomial dipakai ketika variabel respon memiliki lebih
dari dua kategori dan kategorinya tidak memiliki urutan atau tingkatan
(nominal). Dalam laporan ini, respon yang digunakan adalah jenis peran
jabatan karyawan (JobRole), di mana setiap kategori peran
akan dibandingkan terhadap satu kategori acuan (baseline / reference
cell).
Regresi Poisson dipakai ketika variabel respon berupa data cacah atau jumlah frekuensi kejadian (count data). Dalam laporan ini, responnya adalah jumlah dokter yang dikunjungi oleh pasien lansia. Karena data berbentuk individual-level (satu baris per pasien), model Poisson digunakan tanpa offset.
\[ Y_i \sim Poisson(\mu_i) \]
\[ \log(\mu_i) = \beta_0 + \beta_1X_{1i} + \cdots + \beta_pX_{pi} \]
Berbeda dengan data kelompok/agregat, data tingkat individu ini menganalisis rata-rata nilai harapan (\(\mu_i\)) jumlah kunjungan dokter secara langsung berdasarkan kondisi kesehatan fisik, mental, dan gangguan tidur pasien tanpa memerlukan variabel offset.
required_packages <- c(
"readxl", "readr", "dplyr", "tidyr", "ggplot2",
"MASS", "nnet", "broom", "knitr", "scales", "stringr", "forcats", "tibble", "purrr",
"pROC", "ordinal"
)
missing_packages <- required_packages[!(required_packages %in% rownames(installed.packages()))]
if (length(missing_packages) > 0) {
install.packages(missing_packages, repos = "https://cloud.r-project.org")
}
invisible(lapply(required_packages, library, character.only = TRUE))
fmt_num <- function(x, digits = 3) {
ifelse(is.na(x), NA_character_, format(round(x, digits), nsmall = digits, big.mark = ","))
}
fmt_p <- function(x) {
ifelse(is.na(x), NA_character_, ifelse(x < 0.001, "< 0.001", fmt_num(x, 3)))
}
ratio_phrase <- function(ratio, measure, up_text, down_text) {
if (is.na(ratio)) return("nilai rasio tidak tersedia")
if (ratio > 1) {
paste0(measure, " = ", fmt_num(ratio, 3), ", sehingga ", up_text)
} else if (ratio < 1) {
paste0(measure, " = ", fmt_num(ratio, 3), ", sehingga ", down_text)
} else {
paste0(measure, " = 1, sehingga tidak menunjukkan perubahan kecenderungan")
}
}
modus <- function(x) {
ux <- na.omit(unique(x))
ux[which.max(tabulate(match(x, ux)))]
}
# Label untuk model biner (Telco Customer Churn)
label_biner <- function(term) {
dplyr::case_when(
term == "tenure" ~ "durasi berlangganan (bulan)",
term == "MonthlyCharges" ~ "biaya bulanan pelanggan",
term == "TotalCharges" ~ "total biaya kumulatif pelanggan",
term == "SeniorCitizen" ~ "status lansia (1 = Ya, 0 = Tidak)",
stringr::str_starts(term, "gender") ~ paste0("jenis kelamin: ", stringr::str_remove(term, "gender"), " dibanding kategori acuan"),
stringr::str_starts(term, "Contract") ~ paste0("jenis kontrak: ", stringr::str_remove(term, "Contract"), " dibanding kategori acuan"),
stringr::str_starts(term, "InternetService")~ paste0("layanan internet: ", stringr::str_remove(term, "InternetService"), " dibanding kategori acuan"),
stringr::str_starts(term, "PaymentMethod") ~ paste0("metode pembayaran: ", stringr::str_remove(term, "PaymentMethod"), " dibanding kategori acuan"),
stringr::str_starts(term, "OnlineSecurity") ~ paste0("keamanan online: ", stringr::str_remove(term, "OnlineSecurity"), " dibanding kategori acuan"),
stringr::str_starts(term, "OnlineBackup") ~ paste0("pencadangan online: ", stringr::str_remove(term, "OnlineBackup"), " dibanding kategori acuan"),
stringr::str_starts(term, "DeviceProtection")~ paste0("perlindungan perangkat: ", stringr::str_remove(term, "DeviceProtection"), " dibanding kategori acuan"),
stringr::str_starts(term, "TechSupport") ~ paste0("dukungan teknis: ", stringr::str_remove(term, "TechSupport"), " dibanding kategori acuan"),
stringr::str_starts(term, "StreamingTV") ~ paste0("layanan streaming TV: ", stringr::str_remove(term, "StreamingTV"), " dibanding kategori acuan"),
stringr::str_starts(term, "StreamingMovies")~ paste0("layanan streaming film: ", stringr::str_remove(term, "StreamingMovies"), " dibanding kategori acuan"),
stringr::str_starts(term, "MultipleLines") ~ paste0("multiple lines: ", stringr::str_remove(term, "MultipleLines"), " dibanding kategori acuan"),
stringr::str_starts(term, "PhoneService") ~ paste0("layanan telepon: ", stringr::str_remove(term, "PhoneService"), " dibanding kategori acuan"),
stringr::str_starts(term, "Partner") ~ paste0("status pasangan: ", stringr::str_remove(term, "Partner"), " dibanding kategori acuan"),
stringr::str_starts(term, "Dependents") ~ paste0("status tanggungan: ", stringr::str_remove(term, "Dependents"), " dibanding kategori acuan"),
stringr::str_starts(term, "PaperlessBilling")~ paste0("tagihan tanpa kertas: ", stringr::str_remove(term, "PaperlessBilling"), " dibanding kategori acuan"),
TRUE ~ term
)
}
# Label untuk model ordinal (Student Performance)
label_ordinal <- function(term) {
dplyr::case_when(
term == "studytime" ~ "waktu belajar mingguan (skala 1-4)",
term == "failures" ~ "jumlah kegagalan kelas sebelumnya",
term == "Dalc" ~ "konsumsi alkohol hari kerja (skala 1-5)",
term == "Walc" ~ "konsumsi alkohol akhir pekan (skala 1-5)",
term == "age" ~ "usia siswa",
term == "absences" ~ "jumlah absensi/ketidakhadiran sekolah",
term == "traveltime" ~ "waktu perjalanan ke sekolah",
term == "famrel" ~ "kualitas hubungan keluarga",
term == "freetime" ~ "waktu luang setelah sekolah",
term == "goout" ~ "frekuensi pergi bermain bersama teman",
term == "health" ~ "status kesehatan saat ini",
term == "Medu" ~ "tingkat pendidikan ibu",
term == "Fedu" ~ "tingkat pendidikan ayah",
stringr::str_starts(term, "school") ~ paste0("sekolah: ", stringr::str_remove(term, "school"), " dibanding kategori acuan"),
stringr::str_starts(term, "sex") ~ paste0("jenis kelamin siswa: ", stringr::str_remove(term, "sex"), " dibanding kategori acuan"),
stringr::str_starts(term, "address") ~ paste0("tipe alamat tinggal: ", stringr::str_remove(term, "address"), " dibanding kategori acuan"),
stringr::str_starts(term, "Mjob") ~ paste0("pekerjaan ibu: ", stringr::str_remove(term, "Mjob"), " dibanding kategori acuan"),
stringr::str_starts(term, "Fjob") ~ paste0("pekerjaan ayah: ", stringr::str_remove(term, "Fjob"), " dibanding kategori acuan"),
stringr::str_starts(term, "internet") ~ paste0("akses internet di rumah: ", stringr::str_remove(term, "internet"), " dibanding acuan"),
stringr::str_starts(term, "romantic") ~ paste0("memiliki hubungan romantis: ", stringr::str_remove(term, "romantic"), " dibanding acuan"),
stringr::str_starts(term, "famsup") ~ paste0("dukungan pendidikan keluarga: ", stringr::str_remove(term, "famsup"), " dibanding acuan"),
stringr::str_starts(term, "schoolsup")~ paste0("dukungan pendidikan sekolah: ", stringr::str_remove(term, "schoolsup"), " dibanding acuan"),
stringr::str_starts(term, "paid") ~ paste0("les berbayar: ", stringr::str_remove(term, "paid"), " dibanding acuan"),
stringr::str_starts(term, "activities")~paste0("kegiatan ekstrakulikuler: ", stringr::str_remove(term, "activities"), " dibanding acuan"),
stringr::str_starts(term, "higher") ~ paste0("ingin melanjutkan pendidikan tinggi: ", stringr::str_remove(term, "higher"), " dibanding acuan"),
stringr::str_starts(term, "nursery") ~ paste0("pernah ikut taman kanak-kanak: ", stringr::str_remove(term, "nursery"), " dibanding acuan"),
stringr::str_starts(term, "famsize") ~ paste0("ukuran keluarga: ", stringr::str_remove(term, "famsize"), " dibanding acuan"),
stringr::str_starts(term, "Pstatus") ~ paste0("status orang tua: ", stringr::str_remove(term, "Pstatus"), " dibanding acuan"),
stringr::str_starts(term, "reason") ~ paste0("alasan pilih sekolah: ", stringr::str_remove(term, "reason"), " dibanding acuan"),
stringr::str_starts(term, "guardian") ~ paste0("wali siswa: ", stringr::str_remove(term, "guardian"), " dibanding acuan"),
TRUE ~ term
)
}
# Label untuk model multinomial (HR Employee Attrition - JobRole)
label_multinom <- function(term) {
dplyr::case_when(
term == "Age" ~ "usia karyawan",
term == "MonthlyIncome" ~ "pendapatan bulanan karyawan",
term == "DistanceFromHome" ~ "jarak dari rumah ke kantor (km)",
term == "PercentSalaryHike" ~ "persentase kenaikan gaji",
term == "TotalWorkingYears" ~ "total masa kerja (tahun)",
term == "YearsAtCompany" ~ "lama bekerja di perusahaan saat ini",
term == "YearsInCurrentRole" ~ "lama memegang peran jabatan saat ini",
term == "JobSatisfaction" ~ "tingkat kepuasan kerja (skala 1-4)",
term == "WorkLifeBalance" ~ "keseimbangan hidup-kerja (skala 1-4)",
term == "EnvironmentSatisfaction"~ "kepuasan lingkungan kerja (skala 1-4)",
term == "JobInvolvement" ~ "keterlibatan dalam pekerjaan (skala 1-4)",
term == "JobLevel" ~ "level jabatan karyawan",
term == "NumCompaniesWorked" ~ "jumlah perusahaan sebelumnya",
term == "TrainingTimesLastYear" ~ "frekuensi pelatihan tahun lalu",
term == "StockOptionLevel" ~ "level opsi saham karyawan",
term == "YearsSinceLastPromotion"~ "tahun sejak promosi terakhir",
term == "YearsWithCurrManager" ~ "lama dengan manajer saat ini",
term == "DailyRate" ~ "tarif harian karyawan",
term == "HourlyRate" ~ "tarif per jam karyawan",
term == "MonthlyRate" ~ "tarif bulanan karyawan",
stringr::str_starts(term, "Attrition") ~ paste0("status keluar karyawan: ", stringr::str_remove(term, "Attrition"), " dibanding acuan"),
stringr::str_starts(term, "BusinessTravel") ~ paste0("frekuensi perjalanan bisnis: ", stringr::str_remove(term, "BusinessTravel"), " dibanding acuan"),
stringr::str_starts(term, "Department") ~ paste0("departemen kerja: ", stringr::str_remove(term, "Department"), " dibanding acuan"),
stringr::str_starts(term, "EducationField") ~ paste0("bidang pendidikan: ", stringr::str_remove(term, "EducationField"), " dibanding acuan"),
stringr::str_starts(term, "Gender") ~ paste0("jenis kelamin karyawan: ", stringr::str_remove(term, "Gender"), " dibanding acuan"),
stringr::str_starts(term, "MaritalStatus") ~ paste0("status pernikahan: ", stringr::str_remove(term, "MaritalStatus"), " dibanding acuan"),
stringr::str_starts(term, "OverTime") ~ paste0("bekerja lembur: ", stringr::str_remove(term, "OverTime"), " dibanding acuan"),
TRUE ~ term
)
}
# Label untuk model Poisson (NPHA Doctor Visits)
label_poisson <- function(term) {
dplyr::case_when(
term == "Age" ~ "kelompok usia pasien",
term == "Phyiscal.Health" ~ "status kesehatan fisik pasien",
term == "Mental.Health" ~ "status kesehatan mental pasien",
term == "Dental.Health" ~ "status kesehatan gigi pasien",
term == "Employment" ~ "status pekerjaan pasien",
term == "Trouble.Sleeping" ~ "tingkat kesulitan tidur umum",
term == "Stress.Keeps.Patient.from.Sleeping" ~ "gangguan tidur akibat stres",
term == "Medication.Keeps.Patient.from.Sleeping" ~ "gangguan tidur akibat efek obat",
term == "Pain.Keeps.Patient.from.Sleeping" ~ "gangguan tidur akibat rasa nyeri fisik",
term == "Bathroom.Needs.Keeps.Patient.from.Sleeping" ~ "gangguan tidur akibat kebutuhan ke toilet",
term == "Uknown.Keeps.Patient.from.Sleeping" ~ "gangguan tidur akibat penyebab tidak diketahui",
term == "Prescription.Sleep.Medication" ~ "penggunaan obat tidur resep dokter",
term == "Race" ~ "kategori ras pasien",
term == "Gender" ~ "jenis kelamin pasien",
TRUE ~ term
)
}data_dir <- "."
file_biner <- file.path(data_dir, "WA_Fn-UseC_-Telco-Customer-Churn.csv")
file_ordinal <- file.path(data_dir, "student-por.csv")
file_multinom<- file.path(data_dir, "WA_Fn-UseC_-HR-Employee-Attrition.csv")
file_poisson <- file.path(data_dir, "NPHA-doctor-visits.csv")
check_file <- function(path) {
if (!file.exists(path)) {
stop(paste("File tidak ditemukan:", path,
"\nPastikan file Rmd dan dataset berada di folder yang sama atau ubah data_dir."))
}
}
invisible(lapply(c(file_biner, file_ordinal, file_multinom, file_poisson), check_file))
data_biner_raw <- readr::read_csv(file_biner, show_col_types = FALSE)
data_ordinal_raw <- readr::read_csv(file_ordinal, show_col_types = FALSE)
data_multinom_raw<- readr::read_csv(file_multinom, show_col_types = FALSE)
data_poisson_raw <- readr::read_csv(file_poisson, show_col_types = FALSE)dimensi_data <- tibble::tibble(
Dataset = c("Telco Customer Churn", "Student Performance (Portugis)",
"HR Employee Attrition", "NPHA Doctor Visits"),
Jumlah_Baris = c(nrow(data_biner_raw), nrow(data_ordinal_raw),
nrow(data_multinom_raw), nrow(data_poisson_raw)),
Jumlah_Kolom = c(ncol(data_biner_raw), ncol(data_ordinal_raw),
ncol(data_multinom_raw), ncol(data_poisson_raw))
)
knitr::kable(dimensi_data, caption = "Dimensi Data Awal")| Dataset | Jumlah_Baris | Jumlah_Kolom |
|---|---|---|
| Telco Customer Churn | 7043 | 21 |
| Student Performance (Portugis) | 649 | 33 |
| HR Employee Attrition | 1470 | 35 |
| NPHA Doctor Visits | 714 | 15 |
Tujuan analisis ini adalah mengetahui faktor-faktor yang berhubungan dengan peluang pelanggan telekomunikasi melakukan churn (berhenti berlangganan). Churn berarti pelanggan memutuskan untuk tidak melanjutkan layanan dari perusahaan.
Variabel respon adalah Churn:
Nilai Churn |
Arti sederhana |
|---|---|
| No (0) | Pelanggan tidak berhenti berlangganan (loyal) |
| Yes (1) | Pelanggan berhenti berlangganan (churn) |
| Variabel | Peran | Alasan_dipakai |
|---|---|---|
| Churn | Respon | Mewakili status churn pelanggan (biner: Yes/No) |
| gender | Prediktor | Karakteristik demografi pelanggan |
| SeniorCitizen | Prediktor | Status lansia dapat mempengaruhi loyalitas pelanggan |
| Partner | Prediktor | Pelanggan dengan pasangan mungkin lebih stabil |
| Dependents | Prediktor | Tanggungan keluarga dapat mempengaruhi keputusan churn |
| tenure | Prediktor | Semakin lama berlangganan, biasanya semakin loyal |
| PhoneService | Prediktor | Penggunaan layanan telepon |
| MultipleLines | Prediktor | Penggunaan lebih dari satu jalur telepon |
| InternetService | Prediktor | Jenis layanan internet yang digunakan |
| OnlineSecurity | Prediktor | Penggunaan layanan keamanan online |
| OnlineBackup | Prediktor | Penggunaan layanan pencadangan online |
| DeviceProtection | Prediktor | Penggunaan layanan perlindungan perangkat |
| TechSupport | Prediktor | Penggunaan layanan dukungan teknis |
| StreamingTV | Prediktor | Penggunaan layanan streaming TV |
| StreamingMovies | Prediktor | Penggunaan layanan streaming film |
| Contract | Prediktor | Jenis kontrak berpengaruh kuat terhadap churn |
| PaperlessBilling | Prediktor | Metode tagihan |
| PaymentMethod | Prediktor | Metode pembayaran yang digunakan |
| MonthlyCharges | Prediktor | Biaya bulanan yang dibayar pelanggan |
| TotalCharges | Prediktor | Total biaya kumulatif selama berlangganan |
Variabel customerID tidak digunakan karena hanya
merupakan identitas unik pelanggan dan tidak memiliki nilai prediktif.
TotalCharges perlu diperiksa tipe datanya karena terkadang
terbaca sebagai karakter akibat nilai kosong.
data_biner <- data_biner_raw %>%
dplyr::mutate(
TotalCharges = as.numeric(TotalCharges),
Churn_num = ifelse(Churn == "Yes", 1L, 0L),
Churn_label = factor(Churn, levels = c("No", "Yes")),
gender = factor(gender),
Partner = factor(Partner),
Dependents = factor(Dependents),
PhoneService = factor(PhoneService),
MultipleLines = factor(MultipleLines),
InternetService = factor(InternetService),
OnlineSecurity = factor(OnlineSecurity),
OnlineBackup = factor(OnlineBackup),
DeviceProtection = factor(DeviceProtection),
TechSupport = factor(TechSupport),
StreamingTV = factor(StreamingTV),
StreamingMovies = factor(StreamingMovies),
Contract = factor(Contract),
PaperlessBilling = factor(PaperlessBilling),
PaymentMethod = factor(PaymentMethod)
) %>%
tidyr::drop_na()
knitr::kable(head(data_biner %>%
dplyr::select(customerID, gender, SeniorCitizen, tenure, Contract,
MonthlyCharges, TotalCharges, Churn_label), 6),
caption = "Cuplikan Data Biner Setelah Cleaning")| customerID | gender | SeniorCitizen | tenure | Contract | MonthlyCharges | TotalCharges | Churn_label |
|---|---|---|---|---|---|---|---|
| 7590-VHVEG | Female | 0 | 1 | Month-to-month | 29.85 | 29.85 | No |
| 5575-GNVDE | Male | 0 | 34 | One year | 56.95 | 1889.50 | No |
| 3668-QPYBK | Male | 0 | 2 | Month-to-month | 53.85 | 108.15 | Yes |
| 7795-CFOCW | Male | 0 | 45 | One year | 42.30 | 1840.75 | No |
| 9237-HQITU | Female | 0 | 2 | Month-to-month | 70.70 | 151.65 | Yes |
| 9305-CDSKC | Female | 0 | 8 | Month-to-month | 99.65 | 820.50 | Yes |
Variabel TotalCharges dikonversi ke numerik. Baris
dengan nilai kosong dihapus menggunakan drop_na(). Variabel
kategorik seperti gender, Contract,
InternetService, dan PaymentMethod diubah
menjadi factor agar R dapat membuat dummy variable secara otomatis.
dist_biner <- data_biner %>%
dplyr::count(Churn_label) %>%
dplyr::mutate(Persen = scales::percent(n / sum(n)))
knitr::kable(dist_biner, caption = "Distribusi Status Churn Pelanggan")| Churn_label | n | Persen |
|---|---|---|
| No | 5163 | 73% |
| Yes | 1869 | 27% |
summary_biner <- data_biner %>%
dplyr::summarise(
Rata_rata_tenure = mean(tenure),
Median_tenure = median(tenure),
Rata_rata_MonthlyCharges = mean(MonthlyCharges),
Median_MonthlyCharges = median(MonthlyCharges)
)
knitr::kable(summary_biner, caption = "Ringkasan Prediktor Utama Data Biner")| Rata_rata_tenure | Median_tenure | Rata_rata_MonthlyCharges | Median_MonthlyCharges |
|---|---|---|---|
| 32.42179 | 29 | 64.79821 | 70.35 |
ggplot(data_biner, aes(x = Churn_label)) +
geom_bar(fill = c("#2f7f73", "#b84f5a")) +
labs(
title = "Distribusi Status Churn Pelanggan Telekomunikasi",
x = "Status pelanggan",
y = "Jumlah pelanggan"
) +
theme_minimal()Model biner dibentuk untuk menjelaskan peluang
Churn = Yes, yaitu pelanggan berhenti berlangganan.
\[ \log\left(\frac{P(\text{Churn}=1)}{P(\text{Churn}=0)}\right) = \beta_0 + \beta_1\,\text{gender} + \beta_2\,\text{SeniorCitizen} + \cdots + \beta_p X_p \]
model_biner <- glm(
Churn_num ~ gender + SeniorCitizen + Partner + Dependents + tenure +
PhoneService + MultipleLines + InternetService +
OnlineSecurity + OnlineBackup + DeviceProtection + TechSupport +
StreamingTV + StreamingMovies + Contract + PaperlessBilling +
PaymentMethod + MonthlyCharges + TotalCharges,
data = data_biner,
family = binomial(link = "logit")
)
summary(model_biner)##
## Call:
## glm(formula = Churn_num ~ gender + SeniorCitizen + Partner +
## Dependents + tenure + PhoneService + MultipleLines + InternetService +
## OnlineSecurity + OnlineBackup + DeviceProtection + TechSupport +
## StreamingTV + StreamingMovies + Contract + PaperlessBilling +
## PaymentMethod + MonthlyCharges + TotalCharges, family = binomial(link = "logit"),
## data = data_biner)
##
## Coefficients: (7 not defined because of singularities)
## Estimate Std. Error z value
## (Intercept) 1.16528747 0.81513548 1.430
## genderMale -0.02183273 0.06480439 -0.337
## SeniorCitizen 0.21677504 0.08453049 2.564
## PartnerYes -0.00038400 0.07782976 -0.005
## DependentsYes -0.14848782 0.08973115 -1.655
## tenure -0.06058758 0.00623568 -9.716
## PhoneServiceYes 0.17146779 0.64869215 0.264
## MultipleLinesNo phone service NA NA NA
## MultipleLinesYes 0.44839544 0.17725585 2.530
## InternetServiceFiber optic 1.74747491 0.79807951 2.190
## InternetServiceNo -1.78629472 0.80726809 -2.213
## OnlineSecurityNo internet service NA NA NA
## OnlineSecurityYes -0.20542004 0.17868789 -1.150
## OnlineBackupNo internet service NA NA NA
## OnlineBackupYes 0.02604185 0.17540119 0.148
## DeviceProtectionNo internet service NA NA NA
## DeviceProtectionYes 0.14737500 0.17637397 0.836
## TechSupportNo internet service NA NA NA
## TechSupportYes -0.18049681 0.18060228 -0.999
## StreamingTVNo internet service NA NA NA
## StreamingTVYes 0.59050741 0.32630948 1.810
## StreamingMoviesNo internet service NA NA NA
## StreamingMoviesYes 0.59929571 0.32668370 1.834
## ContractOne year -0.66079529 0.10758532 -6.142
## ContractTwo year -1.35710624 0.17644517 -7.691
## PaperlessBillingYes 0.34235364 0.07449538 4.596
## PaymentMethodCredit card (automatic) -0.08779182 0.11407927 -0.770
## PaymentMethodElectronic check 0.30446727 0.09449652 3.222
## PaymentMethodMailed check -0.05758719 0.11491136 -0.501
## MonthlyCharges -0.04034353 0.03175504 -1.270
## TotalCharges 0.00032894 0.00007063 4.657
## Pr(>|z|)
## (Intercept) 0.15284
## genderMale 0.73619
## SeniorCitizen 0.01033 *
## PartnerYes 0.99606
## DependentsYes 0.09796 .
## tenure < 0.0000000000000002 ***
## PhoneServiceYes 0.79153
## MultipleLinesNo phone service NA
## MultipleLinesYes 0.01142 *
## InternetServiceFiber optic 0.02855 *
## InternetServiceNo 0.02691 *
## OnlineSecurityNo internet service NA
## OnlineSecurityYes 0.25031
## OnlineBackupNo internet service NA
## OnlineBackupYes 0.88197
## DeviceProtectionNo internet service NA
## DeviceProtectionYes 0.40339
## TechSupportNo internet service NA
## TechSupportYes 0.31759
## StreamingTVNo internet service NA
## StreamingTVYes 0.07035 .
## StreamingMoviesNo internet service NA
## StreamingMoviesYes 0.06658 .
## ContractOne year 0.0000000008145911 ***
## ContractTwo year 0.0000000000000146 ***
## PaperlessBillingYes 0.0000043143131817 ***
## PaymentMethodCredit card (automatic) 0.44156
## PaymentMethodElectronic check 0.00127 **
## PaymentMethodMailed check 0.61627
## MonthlyCharges 0.20392
## TotalCharges 0.0000032026251598 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 8143.4 on 7031 degrees of freedom
## Residual deviance: 5826.3 on 7008 degrees of freedom
## AIC: 5874.3
##
## Number of Fisher Scoring iterations: 6
Chunk model-biner di atas sudah membangun model
regresi logistik biner menggunakan fungsi glm()
dengan family = binomial(link = "logit"). Persamaan umum
hasil estimasi adalah:
\[ \widehat{\eta_i} = \log\left(\frac{\widehat{P}(\text{Churn}_i=1)}{\widehat{P}(\text{Churn}_i=0)}\right) = \hat{\beta}_0 + \hat{\beta}_1 X_{1i} + \cdots + \hat{\beta}_p X_{pi} \]
Dengan:
\[ \widehat{P}(\text{Churn}_i=1) = \frac{\exp(\widehat{\eta_i})}{1+\exp(\widehat{\eta_i})} \]
| Variabel_Mudah | term | estimate | Odds_Ratio |
|---|---|---|---|
| (Intercept) | (Intercept) | 1.1653 | 3.2068 |
| jenis kelamin: Male dibanding kategori acuan | genderMale | -0.0218 | 0.9784 |
| status lansia (1 = Ya, 0 = Tidak) | SeniorCitizen | 0.2168 | 1.2421 |
| status pasangan: Yes dibanding kategori acuan | PartnerYes | -0.0004 | 0.9996 |
| status tanggungan: Yes dibanding kategori acuan | DependentsYes | -0.1485 | 0.8620 |
| durasi berlangganan (bulan) | tenure | -0.0606 | 0.9412 |
| layanan telepon: Yes dibanding kategori acuan | PhoneServiceYes | 0.1715 | 1.1870 |
| multiple lines: No phone service dibanding kategori acuan | MultipleLinesNo phone service | NA | NA |
| multiple lines: Yes dibanding kategori acuan | MultipleLinesYes | 0.4484 | 1.5658 |
| layanan internet: Fiber optic dibanding kategori acuan | InternetServiceFiber optic | 1.7475 | 5.7401 |
| layanan internet: No dibanding kategori acuan | InternetServiceNo | -1.7863 | 0.1676 |
| keamanan online: No internet service dibanding kategori acuan | OnlineSecurityNo internet service | NA | NA |
| keamanan online: Yes dibanding kategori acuan | OnlineSecurityYes | -0.2054 | 0.8143 |
| pencadangan online: No internet service dibanding kategori acuan | OnlineBackupNo internet service | NA | NA |
| pencadangan online: Yes dibanding kategori acuan | OnlineBackupYes | 0.0260 | 1.0264 |
| perlindungan perangkat: No internet service dibanding kategori acuan | DeviceProtectionNo internet service | NA | NA |
| perlindungan perangkat: Yes dibanding kategori acuan | DeviceProtectionYes | 0.1474 | 1.1588 |
| dukungan teknis: No internet service dibanding kategori acuan | TechSupportNo internet service | NA | NA |
| dukungan teknis: Yes dibanding kategori acuan | TechSupportYes | -0.1805 | 0.8349 |
| layanan streaming TV: No internet service dibanding kategori acuan | StreamingTVNo internet service | NA | NA |
| layanan streaming TV: Yes dibanding kategori acuan | StreamingTVYes | 0.5905 | 1.8049 |
| layanan streaming film: No internet service dibanding kategori acuan | StreamingMoviesNo internet service | NA | NA |
| layanan streaming film: Yes dibanding kategori acuan | StreamingMoviesYes | 0.5993 | 1.8208 |
| jenis kontrak: One year dibanding kategori acuan | ContractOne year | -0.6608 | 0.5164 |
| jenis kontrak: Two year dibanding kategori acuan | ContractTwo year | -1.3571 | 0.2574 |
| tagihan tanpa kertas: Yes dibanding kategori acuan | PaperlessBillingYes | 0.3424 | 1.4083 |
| metode pembayaran: Credit card (automatic) dibanding kategori acuan | PaymentMethodCredit card (automatic) | -0.0878 | 0.9160 |
| metode pembayaran: Electronic check dibanding kategori acuan | PaymentMethodElectronic check | 0.3045 | 1.3559 |
| metode pembayaran: Mailed check dibanding kategori acuan | PaymentMethodMailed check | -0.0576 | 0.9440 |
| biaya bulanan pelanggan | MonthlyCharges | -0.0403 | 0.9605 |
tabel_biner <- broom::tidy(model_biner) %>%
dplyr::mutate(
Odds_Ratio = exp(estimate),
Variabel_Mudah = purrr::map_chr(term, label_biner),
Keputusan_5persen= ifelse(p.value < 0.05, "Signifikan", "Tidak signifikan")
) %>%
dplyr::arrange(p.value)
knitr::kable(
tabel_biner %>%
dplyr::select(Variabel_Mudah, term, estimate, std.error,
statistic, p.value, Odds_Ratio, Keputusan_5persen) %>%
head(25),
digits = 4,
caption = "Koefisien dan Odds Ratio Regresi Logistik Biner"
)| Variabel_Mudah | term | estimate | std.error | statistic | p.value | Odds_Ratio | Keputusan_5persen |
|---|---|---|---|---|---|---|---|
| durasi berlangganan (bulan) | tenure | -0.0606 | 0.0062 | -9.7163 | 0.0000 | 0.9412 | Signifikan |
| jenis kontrak: Two year dibanding kategori acuan | ContractTwo year | -1.3571 | 0.1764 | -7.6914 | 0.0000 | 0.2574 | Signifikan |
| jenis kontrak: One year dibanding kategori acuan | ContractOne year | -0.6608 | 0.1076 | -6.1421 | 0.0000 | 0.5164 | Signifikan |
| total biaya kumulatif pelanggan | TotalCharges | 0.0003 | 0.0001 | 4.6574 | 0.0000 | 1.0003 | Signifikan |
| tagihan tanpa kertas: Yes dibanding kategori acuan | PaperlessBillingYes | 0.3424 | 0.0745 | 4.5956 | 0.0000 | 1.4083 | Signifikan |
| metode pembayaran: Electronic check dibanding kategori acuan | PaymentMethodElectronic check | 0.3045 | 0.0945 | 3.2220 | 0.0013 | 1.3559 | Signifikan |
| status lansia (1 = Ya, 0 = Tidak) | SeniorCitizen | 0.2168 | 0.0845 | 2.5645 | 0.0103 | 1.2421 | Signifikan |
| multiple lines: Yes dibanding kategori acuan | MultipleLinesYes | 0.4484 | 0.1773 | 2.5297 | 0.0114 | 1.5658 | Signifikan |
| layanan internet: No dibanding kategori acuan | InternetServiceNo | -1.7863 | 0.8073 | -2.2128 | 0.0269 | 0.1676 | Signifikan |
| layanan internet: Fiber optic dibanding kategori acuan | InternetServiceFiber optic | 1.7475 | 0.7981 | 2.1896 | 0.0286 | 5.7401 | Signifikan |
| layanan streaming film: Yes dibanding kategori acuan | StreamingMoviesYes | 0.5993 | 0.3267 | 1.8345 | 0.0666 | 1.8208 | Tidak signifikan |
| layanan streaming TV: Yes dibanding kategori acuan | StreamingTVYes | 0.5905 | 0.3263 | 1.8097 | 0.0703 | 1.8049 | Tidak signifikan |
| status tanggungan: Yes dibanding kategori acuan | DependentsYes | -0.1485 | 0.0897 | -1.6548 | 0.0980 | 0.8620 | Tidak signifikan |
| (Intercept) | (Intercept) | 1.1653 | 0.8151 | 1.4296 | 0.1528 | 3.2068 | Tidak signifikan |
| biaya bulanan pelanggan | MonthlyCharges | -0.0403 | 0.0318 | -1.2705 | 0.2039 | 0.9605 | Tidak signifikan |
| keamanan online: Yes dibanding kategori acuan | OnlineSecurityYes | -0.2054 | 0.1787 | -1.1496 | 0.2503 | 0.8143 | Tidak signifikan |
| dukungan teknis: Yes dibanding kategori acuan | TechSupportYes | -0.1805 | 0.1806 | -0.9994 | 0.3176 | 0.8349 | Tidak signifikan |
| perlindungan perangkat: Yes dibanding kategori acuan | DeviceProtectionYes | 0.1474 | 0.1764 | 0.8356 | 0.4034 | 1.1588 | Tidak signifikan |
| metode pembayaran: Credit card (automatic) dibanding kategori acuan | PaymentMethodCredit card (automatic) | -0.0878 | 0.1141 | -0.7696 | 0.4416 | 0.9160 | Tidak signifikan |
| metode pembayaran: Mailed check dibanding kategori acuan | PaymentMethodMailed check | -0.0576 | 0.1149 | -0.5011 | 0.6163 | 0.9440 | Tidak signifikan |
| jenis kelamin: Male dibanding kategori acuan | genderMale | -0.0218 | 0.0648 | -0.3369 | 0.7362 | 0.9784 | Tidak signifikan |
| layanan telepon: Yes dibanding kategori acuan | PhoneServiceYes | 0.1715 | 0.6487 | 0.2643 | 0.7915 | 1.1870 | Tidak signifikan |
| pencadangan online: Yes dibanding kategori acuan | OnlineBackupYes | 0.0260 | 0.1754 | 0.1485 | 0.8820 | 1.0264 | Tidak signifikan |
| status pasangan: Yes dibanding kategori acuan | PartnerYes | -0.0004 | 0.0778 | -0.0049 | 0.9961 | 0.9996 | Tidak signifikan |
| multiple lines: No phone service dibanding kategori acuan | MultipleLinesNo phone service | NA | NA | NA | NA | NA | NA |
data_biner$prob_churn <- predict(model_biner, type = "response")
data_biner$pred_churn <- ifelse(data_biner$prob_churn >= 0.5, 1, 0)
tabel_klasifikasi_biner <- table(
Prediksi = factor(data_biner$pred_churn, levels = c(0, 1),
labels = c("Tidak churn", "Churn")),
Aktual = data_biner$Churn_label
)
tabel_klasifikasi_biner## Aktual
## Prediksi No Yes
## Tidak churn 4639 833
## Churn 524 1036
## [1] 0.807025
## [1] 5874.272
Visualisasi ini digunakan untuk melihat sebaran probabilitas churn yang dihasilkan oleh model. Probabilitas yang lebih tinggi menunjukkan pelanggan yang diprediksi memiliki risiko churn lebih besar.
ggplot(data_biner, aes(x = prob_churn)) +
geom_histogram(bins = 40, fill = "#2f7f73", color = "white") +
labs(
title = "Distribusi Probabilitas Prediksi Churn",
x = "Probabilitas prediksi churn",
y = "Jumlah pelanggan"
) +
theme_minimal()ggplot(data_biner, aes(x = Churn_label, y = prob_churn, fill = Churn_label)) +
geom_boxplot() +
scale_fill_manual(values = c("No" = "#2f7f73", "Yes" = "#b84f5a")) +
labs(
title = "Probabilitas Prediksi Churn menurut Status Aktual",
x = "Status aktual",
y = "Probabilitas prediksi churn",
fill = "Status churn"
) +
theme_minimal()Grafik pertama menunjukkan sebaran probabilitas churn hasil model. Grafik kedua membandingkan probabilitas prediksi antara pelanggan yang benar-benar churn dan tidak churn. Jika model cukup membedakan risiko, maka kelompok yang aktual churn seharusnya memiliki probabilitas prediksi yang lebih tinggi.
ROC Curve digunakan untuk melihat kemampuan model membedakan pelanggan churn dan tidak churn pada berbagai threshold. AUC mendekati 1 menunjukkan kemampuan pemisahan yang lebih baik, sedangkan AUC mendekati 0,5 menunjukkan kemampuan model mendekati tebakan acak.
roc_biner <- pROC::roc(
response = data_biner$Churn_num,
predictor = data_biner$prob_churn,
quiet = TRUE
)
auc_biner <- pROC::auc(roc_biner)
auc_biner## Area under the curve: 0.8482
AUC menunjukkan kemampuan model membedakan pelanggan yang churn dan tidak churn. Semakin besar AUC, semakin baik kemampuan diskriminasi model. Namun AUC tidak menjelaskan arah pengaruh variabel; arah pengaruh tetap dibaca melalui koefisien dan odds ratio.
Regresi logistik biner digunakan karena variabel respon hanya memiliki dua kategori, yaitu churn dan tidak churn. Model ini mengasumsikan bahwa observasi bersifat independen, tidak terdapat multikolinearitas berat antar prediktor, tidak terjadi perfect separation, dan prediktor numerik memiliki hubungan yang cukup linear terhadap log-odds churn.
Dalam analisis churn, jumlah pelanggan yang churn sering kali lebih kecil daripada pelanggan yang tidak churn (kelas tidak seimbang). Oleh karena itu, akurasi tidak boleh menjadi satu-satunya ukuran evaluasi. Confusion matrix, probabilitas prediksi, dan AUC perlu dibaca bersama dengan interpretasi odds ratio.
sig_biner <- tabel_biner %>%
dplyr::filter(term != "(Intercept)", p.value < 0.05) %>%
dplyr::arrange(p.value)
cat("Pada taraf signifikansi 5%, model biner menemukan ", nrow(sig_biner),
" prediktor yang berhubungan signifikan dengan peluang churn pelanggan.\n\n", sep = "")Pada taraf signifikansi 5%, model biner menemukan 10 prediktor yang berhubungan signifikan dengan peluang churn pelanggan.
if (nrow(sig_biner) > 0) {
cat("Prediktor paling penting berdasarkan p-value terkecil adalah:\n\n")
top_biner <- head(sig_biner, 8)
for (i in seq_len(nrow(top_biner))) {
row <- top_biner[i, ]
arah <- ifelse(row$Odds_Ratio > 1,
"meningkatkan odds churn pelanggan",
"menurunkan odds churn pelanggan")
cat("- **", row$Variabel_Mudah, "** memiliki OR = ", fmt_num(row$Odds_Ratio, 3),
" dan p-value = ", fmt_p(row$p.value), ", sehingga variabel ini ", arah,
", dengan asumsi variabel lain tetap.\n", sep = "")
}
cat("\n")
}Prediktor paling penting berdasarkan p-value terkecil adalah:
cat("Akurasi klasifikasi model biner adalah **", scales::percent(akurasi_biner),
"**. AIC model adalah **", fmt_num(AIC_biner, 2), "**.\n\n", sep = "")Akurasi klasifikasi model biner adalah 81%. AIC model adalah 5,874.27.
cat("**Kesimpulan substantif:** jenis kontrak, lama berlangganan, dan jenis layanan internet merupakan faktor utama yang menjelaskan peluang churn pelanggan telekomunikasi. Pelanggan dengan kontrak bulanan dan berlangganan lebih singkat cenderung memiliki risiko churn yang lebih tinggi.\n")Kesimpulan substantif: jenis kontrak, lama berlangganan, dan jenis layanan internet merupakan faktor utama yang menjelaskan peluang churn pelanggan telekomunikasi. Pelanggan dengan kontrak bulanan dan berlangganan lebih singkat cenderung memiliki risiko churn yang lebih tinggi.
Tujuan analisis ini adalah mengetahui faktor-faktor yang berhubungan
dengan tingkat performa akademik siswa yang diukur
melalui nilai akhir (G3) pada mata pelajaran bahasa
Portugis. Nilai G3 dikelompokkan menjadi tiga kategori ordinal: Rendah,
Sedang, dan Tinggi.
Variabel respon adalah Performa_G3 dengan urutan:
\[ \text{Rendah} < \text{Sedang} < \text{Tinggi} \]
| Variabel | Peran | Alasan_dipakai |
|---|---|---|
| Performa_G3 | Respon | Nilai akhir G3 dikelompokkan menjadi Rendah/Sedang/Tinggi |
| studytime | Prediktor | Waktu belajar berkaitan langsung dengan capaian nilai |
| failures | Prediktor | Riwayat kegagalan akademik sebelumnya |
| Dalc | Prediktor | Konsumsi alkohol hari kerja (skala 1-5) |
| Walc | Prediktor | Konsumsi alkohol akhir pekan (skala 1-5) |
| absences | Prediktor | Banyaknya ketidakhadiran dapat menurunkan capaian |
| famrel | Prediktor | Hubungan keluarga yang baik mendukung belajar |
| goout | Prediktor | Frekuensi bermain bersama teman |
| health | Prediktor | Status kesehatan siswa |
| Medu | Prediktor | Pendidikan ibu berkaitan dengan dukungan belajar di rumah |
| Fedu | Prediktor | Pendidikan ayah berkaitan dengan dukungan belajar di rumah |
| age | Prediktor | Usia siswa |
| internet | Prediktor | Akses internet di rumah |
| higher | Prediktor | Motivasi melanjutkan pendidikan tinggi |
| traveltime | Prediktor | Waktu perjalanan ke sekolah |
Variabel G1 dan G2 tidak digunakan dalam
model karena merupakan nilai antara yang sangat berkorelasi dengan
G3. Penggunaan keduanya akan menyebabkan data leakage
karena nilai akhir G3 sangat ditentukan oleh G1 dan G2. Selain itu,
variabel school, sex, dan variabel kategorik
lainnya digunakan sebagai prediktor tambahan.
data_ordinal <- data_ordinal_raw %>%
dplyr::mutate(
Performa_G3 = dplyr::case_when(
G3 <= 9 ~ "Rendah",
G3 <= 14 ~ "Sedang",
G3 >= 15 ~ "Tinggi"
),
Performa_G3 = ordered(Performa_G3, levels = c("Rendah", "Sedang", "Tinggi")),
school = factor(school),
sex = factor(sex),
address = factor(address),
famsize = factor(famsize),
Pstatus = factor(Pstatus),
Mjob = factor(Mjob),
Fjob = factor(Fjob),
reason = factor(reason),
guardian = factor(guardian),
schoolsup = factor(schoolsup),
famsup = factor(famsup),
paid = factor(paid),
activities = factor(activities),
nursery = factor(nursery),
higher = factor(higher),
internet = factor(internet),
romantic = factor(romantic)
) %>%
tidyr::drop_na()
knitr::kable(
head(data_ordinal %>%
dplyr::select(school, sex, age, studytime, failures, Dalc, Walc,
absences, G3, Performa_G3), 6),
caption = "Cuplikan Data Ordinal Setelah Cleaning"
)| school | sex | age | studytime | failures | Dalc | Walc | absences | G3 | Performa_G3 |
|---|---|---|---|---|---|---|---|---|---|
| GP | F | 18 | 2 | 0 | 1 | 1 | 4 | 11 | Sedang |
| GP | F | 17 | 2 | 0 | 1 | 1 | 2 | 11 | Sedang |
| GP | F | 15 | 2 | 0 | 2 | 3 | 6 | 12 | Sedang |
| GP | F | 15 | 3 | 0 | 1 | 1 | 0 | 14 | Sedang |
| GP | F | 16 | 2 | 0 | 1 | 2 | 0 | 13 | Sedang |
| GP | M | 16 | 2 | 0 | 1 | 2 | 6 | 13 | Sedang |
Nilai G3 dikelompokkan menjadi tiga kategori ordinal:
Rendah (0–9), Sedang (10–14), dan Tinggi (15–20). Pengelompokan ini
mengacu pada konvensi sistem penilaian Portugis dan bertujuan agar model
ordinal dapat diaplikasikan dengan tiga tingkatan yang bermakna secara
substantif.
dist_ordinal <- data_ordinal %>%
dplyr::count(Performa_G3) %>%
dplyr::mutate(Persen = scales::percent(n / sum(n)))
knitr::kable(dist_ordinal, caption = "Distribusi Performa G3 Siswa")| Performa_G3 | n | Persen |
|---|---|---|
| Rendah | 100 | 15.4% |
| Sedang | 418 | 64.4% |
| Tinggi | 131 | 20.2% |
ringkasan_ordinal <- data_ordinal %>%
dplyr::summarise(
Rata_rata_G3 = mean(G3),
Median_G3 = median(G3),
Rata_rata_studytime= mean(studytime),
Rata_rata_absences = mean(absences),
Rata_rata_failures = mean(failures)
)
knitr::kable(ringkasan_ordinal, caption = "Ringkasan Prediktor Utama Data Ordinal")| Rata_rata_G3 | Median_G3 | Rata_rata_studytime | Rata_rata_absences | Rata_rata_failures |
|---|---|---|---|---|
| 11.90601 | 12 | 1.930663 | 3.659476 | 0.2218798 |
ggplot(data_ordinal, aes(x = Performa_G3, fill = Performa_G3)) +
geom_bar() +
scale_fill_manual(values = c("Rendah" = "#b84f5a", "Sedang" = "#d18b2f", "Tinggi" = "#2f7f73")) +
labs(
title = "Distribusi Performa Akademik Siswa (G3)",
x = "Kategori performa",
y = "Jumlah siswa",
fill = "Performa"
) +
theme_minimal()Model ordinal dibentuk untuk menjelaskan faktor-faktor yang
berhubungan dengan tingkat performa akademik siswa. Model ini
menggunakan spesifikasi MASS::polr() sebagai berikut:
\[ \log\left(\frac{P(Y \leq j)}{P(Y > j)}\right) = \alpha_j - x\beta \]
Dalam penelitian ini, variabel respon adalah
Performa_G3, dengan urutan:
\[ \text{Rendah} < \text{Sedang} < \text{Tinggi} \]
Karena model menggunakan bentuk \(\alpha_j - x\beta\), maka koefisien positif menunjukkan penurunan peluang kumulatif berada pada kategori rendah. Dengan kata lain, koefisien positif menunjukkan peningkatan kecenderungan siswa untuk bergeser ke tingkat performa yang lebih tinggi.
model_ordinal <- MASS::polr(
Performa_G3 ~ school + sex + age + address + famsize + Pstatus +
Medu + Fedu + Mjob + Fjob + reason + guardian +
traveltime + studytime + failures + schoolsup + famsup + paid +
activities + nursery + higher + internet + romantic +
famrel + freetime + goout + Dalc + Walc + health + absences,
data = data_ordinal,
Hess = TRUE
)
summary(model_ordinal)## Call:
## MASS::polr(formula = Performa_G3 ~ school + sex + age + address +
## famsize + Pstatus + Medu + Fedu + Mjob + Fjob + reason +
## guardian + traveltime + studytime + failures + schoolsup +
## famsup + paid + activities + nursery + higher + internet +
## romantic + famrel + freetime + goout + Dalc + Walc + health +
## absences, data = data_ordinal, Hess = TRUE)
##
## Coefficients:
## Value Std. Error t value
## schoolMS -1.135536 0.23216 -4.89120
## sexM -0.446014 0.20402 -2.18617
## age 0.280442 0.08380 3.34650
## addressU 0.009603 0.21598 0.04447
## famsizeLE3 0.295601 0.20013 1.47705
## PstatusT -0.014614 0.28619 -0.05107
## Medu 0.172790 0.12616 1.36960
## Fedu 0.066272 0.11440 0.57928
## Mjobhealth 0.441884 0.43667 1.01195
## Mjobother 0.026157 0.24728 0.10578
## Mjobservices 0.247434 0.30451 0.81257
## Mjobteacher 0.218389 0.40863 0.53444
## Fjobhealth -0.603004 0.62190 -0.96962
## Fjobother -0.328684 0.37268 -0.88194
## Fjobservices -0.595829 0.39149 -1.52194
## Fjobteacher 0.134257 0.54458 0.24653
## reasonhome 0.087483 0.23226 0.37666
## reasonother -0.097900 0.30474 -0.32126
## reasonreputation 0.318822 0.24219 1.31639
## guardianmother -0.342715 0.21594 -1.58709
## guardianother -0.228469 0.43264 -0.52809
## traveltime 0.127330 0.13257 0.96044
## studytime 0.255165 0.11373 2.24355
## failures -1.024736 0.18080 -5.66772
## schoolsupyes -1.076297 0.30649 -3.51163
## famsupyes -0.087918 0.18617 -0.47223
## paidyes -0.619844 0.38686 -1.60223
## activitiesyes 0.087654 0.18229 0.48085
## nurseryyes -0.189096 0.22327 -0.84693
## higheryes 1.492243 0.31799 4.69273
## internetyes 0.099099 0.22815 0.43437
## romanticyes -0.119833 0.18764 -0.63865
## famrel 0.098758 0.09529 1.03643
## freetime -0.026837 0.09290 -0.28889
## goout -0.055981 0.08845 -0.63289
## Dalc -0.142805 0.12732 -1.12161
## Walc -0.039912 0.09733 -0.41009
## health -0.110034 0.06371 -1.72701
## absences -0.059224 0.02117 -2.79714
##
## Intercepts:
## Value Std. Error t value
## Rendah|Sedang 2.7716 1.6248 1.7059
## Sedang|Tinggi 6.9768 1.6527 4.2214
##
## Residual Deviance: 927.9076
## AIC: 1009.908
Chunk model-ordinal di atas sudah membangun
model regresi logistik ordinal menggunakan
MASS::polr(). Bentuk hasil estimasinya adalah:
\[ \log\left(\frac{\widehat{P}(\text{Performa}_i \leq j)}{\widehat{P}(\text{Performa}_i > j)}\right) = \hat{\alpha}_j - \widehat{\eta_i} \]
Dengan:
\[ \widehat{\eta_i} = \hat{\beta}_1X_{1i} + \hat{\beta}_2X_{2i} + \cdots + \hat{\beta}_pX_{pi} \]
Karena MASS::polr() memakai bentuk
alpha_j - x beta, maka exp(beta) dibaca
sebagai odds menuju performa yang lebih tinggi.
| Parameter | Estimate | Std_Error | t_value | Jenis | Odds_Ratio_Higher_Perf |
|---|---|---|---|---|---|
| schoolMS | -1.1355 | 0.2322 | -4.8912 | Koefisien prediktor | 0.3212 |
| sexM | -0.4460 | 0.2040 | -2.1862 | Koefisien prediktor | 0.6402 |
| age | 0.2804 | 0.0838 | 3.3465 | Koefisien prediktor | 1.3237 |
| addressU | 0.0096 | 0.2160 | 0.0445 | Koefisien prediktor | 1.0096 |
| famsizeLE3 | 0.2956 | 0.2001 | 1.4770 | Koefisien prediktor | 1.3439 |
| PstatusT | -0.0146 | 0.2862 | -0.0511 | Koefisien prediktor | 0.9855 |
| Medu | 0.1728 | 0.1262 | 1.3696 | Koefisien prediktor | 1.1886 |
| Fedu | 0.0663 | 0.1144 | 0.5793 | Koefisien prediktor | 1.0685 |
| Mjobhealth | 0.4419 | 0.4367 | 1.0120 | Koefisien prediktor | 1.5556 |
| Mjobother | 0.0262 | 0.2473 | 0.1058 | Koefisien prediktor | 1.0265 |
| Mjobservices | 0.2474 | 0.3045 | 0.8126 | Koefisien prediktor | 1.2807 |
| Mjobteacher | 0.2184 | 0.4086 | 0.5344 | Koefisien prediktor | 1.2441 |
| Fjobhealth | -0.6030 | 0.6219 | -0.9696 | Koefisien prediktor | 0.5472 |
| Fjobother | -0.3287 | 0.3727 | -0.8819 | Koefisien prediktor | 0.7199 |
| Fjobservices | -0.5958 | 0.3915 | -1.5219 | Koefisien prediktor | 0.5511 |
| Fjobteacher | 0.1343 | 0.5446 | 0.2465 | Koefisien prediktor | 1.1437 |
| reasonhome | 0.0875 | 0.2323 | 0.3767 | Koefisien prediktor | 1.0914 |
| reasonother | -0.0979 | 0.3047 | -0.3213 | Koefisien prediktor | 0.9067 |
| reasonreputation | 0.3188 | 0.2422 | 1.3164 | Koefisien prediktor | 1.3755 |
| guardianmother | -0.3427 | 0.2159 | -1.5871 | Koefisien prediktor | 0.7098 |
| guardianother | -0.2285 | 0.4326 | -0.5281 | Koefisien prediktor | 0.7958 |
| traveltime | 0.1273 | 0.1326 | 0.9604 | Koefisien prediktor | 1.1358 |
| studytime | 0.2552 | 0.1137 | 2.2435 | Koefisien prediktor | 1.2907 |
| failures | -1.0247 | 0.1808 | -5.6677 | Koefisien prediktor | 0.3589 |
| schoolsupyes | -1.0763 | 0.3065 | -3.5116 | Koefisien prediktor | 0.3409 |
| famsupyes | -0.0879 | 0.1862 | -0.4722 | Koefisien prediktor | 0.9158 |
| paidyes | -0.6198 | 0.3869 | -1.6022 | Koefisien prediktor | 0.5380 |
| activitiesyes | 0.0877 | 0.1823 | 0.4808 | Koefisien prediktor | 1.0916 |
| nurseryyes | -0.1891 | 0.2233 | -0.8469 | Koefisien prediktor | 0.8277 |
| higheryes | 1.4922 | 0.3180 | 4.6927 | Koefisien prediktor | 4.4471 |
| internetyes | 0.0991 | 0.2281 | 0.4344 | Koefisien prediktor | 1.1042 |
| romanticyes | -0.1198 | 0.1876 | -0.6386 | Koefisien prediktor | 0.8871 |
| famrel | 0.0988 | 0.0953 | 1.0364 | Koefisien prediktor | 1.1038 |
| freetime | -0.0268 | 0.0929 | -0.2889 | Koefisien prediktor | 0.9735 |
| goout | -0.0560 | 0.0885 | -0.6329 | Koefisien prediktor | 0.9456 |
| Dalc | -0.1428 | 0.1273 | -1.1216 | Koefisien prediktor | 0.8669 |
| Walc | -0.0399 | 0.0973 | -0.4101 | Koefisien prediktor | 0.9609 |
| health | -0.1100 | 0.0637 | -1.7270 | Koefisien prediktor | 0.8958 |
| absences | -0.0592 | 0.0212 | -2.7971 | Koefisien prediktor | 0.9425 |
| Rendah|Sedang | 2.7716 | 1.6248 | 1.7059 | Cutpoint | NA |
| Sedang|Tinggi | 6.9768 | 1.6527 | 4.2214 | Cutpoint | NA |
Pada tabel berikut,
Odds_Ratio_Higher_Perf = exp(estimate) dibaca sebagai odds
untuk berada pada kategori performa yang lebih tinggi. Ini sesuai dengan
parameterisasi MASS::polr() yang berbentuk:
\[ \operatorname{logit}\{P(Y \leq j)\} = \alpha_j - x\beta \]
coef_ordinal <- coef(summary(model_ordinal)) %>%
as.data.frame()
coef_ordinal$term <- rownames(coef_ordinal)
names(coef_ordinal)[1:3] <- c("estimate", "std_error", "t_value")
coef_ordinal <- coef_ordinal %>%
dplyr::mutate(
p_value = 2 * pnorm(abs(t_value), lower.tail = FALSE),
Odds_Ratio_Higher_Perf = exp(estimate),
Odds_Ratio_Lower_Cumul = exp(-estimate),
Jenis = ifelse(stringr::str_detect(term, "\\|"), "Cutpoint", "Prediktor"),
Variabel_Mudah = purrr::map_chr(term, label_ordinal),
Keputusan_5persen = ifelse(p_value < 0.05, "Signifikan", "Tidak signifikan")
) %>%
dplyr::relocate(term, Jenis)
knitr::kable(
coef_ordinal %>%
dplyr::filter(Jenis == "Prediktor") %>%
dplyr::arrange(p_value) %>%
dplyr::select(Variabel_Mudah, term, estimate, std_error, t_value,
p_value, Odds_Ratio_Higher_Perf, Odds_Ratio_Lower_Cumul,
Keputusan_5persen) %>%
head(25),
digits = 4,
caption = "Koefisien dan Odds Ratio Regresi Ordinal MASS::polr"
)| Variabel_Mudah | term | estimate | std_error | t_value | p_value | Odds_Ratio_Higher_Perf | Odds_Ratio_Lower_Cumul | Keputusan_5persen | |
|---|---|---|---|---|---|---|---|---|---|
| failures | jumlah kegagalan kelas sebelumnya | failures | -1.0247 | 0.1808 | -5.6677 | 0.0000 | 0.3589 | 2.7864 | Signifikan |
| schoolMS | sekolah: MS dibanding kategori acuan | schoolMS | -1.1355 | 0.2322 | -4.8912 | 0.0000 | 0.3212 | 3.1128 | Signifikan |
| higheryes | ingin melanjutkan pendidikan tinggi: yes dibanding acuan | higheryes | 1.4922 | 0.3180 | 4.6927 | 0.0000 | 4.4471 | 0.2249 | Signifikan |
| schoolsupyes | sekolah: supyes dibanding kategori acuan | schoolsupyes | -1.0763 | 0.3065 | -3.5116 | 0.0004 | 0.3409 | 2.9338 | Signifikan |
| age | usia siswa | age | 0.2804 | 0.0838 | 3.3465 | 0.0008 | 1.3237 | 0.7554 | Signifikan |
| absences | jumlah absensi/ketidakhadiran sekolah | absences | -0.0592 | 0.0212 | -2.7971 | 0.0052 | 0.9425 | 1.0610 | Signifikan |
| studytime | waktu belajar mingguan (skala 1-4) | studytime | 0.2552 | 0.1137 | 2.2435 | 0.0249 | 1.2907 | 0.7748 | Signifikan |
| sexM | jenis kelamin siswa: M dibanding kategori acuan | sexM | -0.4460 | 0.2040 | -2.1862 | 0.0288 | 0.6402 | 1.5621 | Signifikan |
| health | status kesehatan saat ini | health | -0.1100 | 0.0637 | -1.7270 | 0.0842 | 0.8958 | 1.1163 | Tidak signifikan |
| paidyes | les berbayar: yes dibanding acuan | paidyes | -0.6198 | 0.3869 | -1.6022 | 0.1091 | 0.5380 | 1.8586 | Tidak signifikan |
| guardianmother | wali siswa: mother dibanding acuan | guardianmother | -0.3427 | 0.2159 | -1.5871 | 0.1125 | 0.7098 | 1.4088 | Tidak signifikan |
| Fjobservices | pekerjaan ayah: services dibanding kategori acuan | Fjobservices | -0.5958 | 0.3915 | -1.5219 | 0.1280 | 0.5511 | 1.8145 | Tidak signifikan |
| famsizeLE3 | ukuran keluarga: LE3 dibanding acuan | famsizeLE3 | 0.2956 | 0.2001 | 1.4770 | 0.1397 | 1.3439 | 0.7441 | Tidak signifikan |
| Medu | tingkat pendidikan ibu | Medu | 0.1728 | 0.1262 | 1.3696 | 0.1708 | 1.1886 | 0.8413 | Tidak signifikan |
| reasonreputation | alasan pilih sekolah: reputation dibanding acuan | reasonreputation | 0.3188 | 0.2422 | 1.3164 | 0.1880 | 1.3755 | 0.7270 | Tidak signifikan |
| Dalc | konsumsi alkohol hari kerja (skala 1-5) | Dalc | -0.1428 | 0.1273 | -1.1216 | 0.2620 | 0.8669 | 1.1535 | Tidak signifikan |
| famrel | kualitas hubungan keluarga | famrel | 0.0988 | 0.0953 | 1.0364 | 0.3000 | 1.1038 | 0.9060 | Tidak signifikan |
| Mjobhealth | pekerjaan ibu: health dibanding kategori acuan | Mjobhealth | 0.4419 | 0.4367 | 1.0120 | 0.3116 | 1.5556 | 0.6428 | Tidak signifikan |
| Fjobhealth | pekerjaan ayah: health dibanding kategori acuan | Fjobhealth | -0.6030 | 0.6219 | -0.9696 | 0.3322 | 0.5472 | 1.8276 | Tidak signifikan |
| traveltime | waktu perjalanan ke sekolah | traveltime | 0.1273 | 0.1326 | 0.9604 | 0.3368 | 1.1358 | 0.8804 | Tidak signifikan |
| Fjobother | pekerjaan ayah: other dibanding kategori acuan | Fjobother | -0.3287 | 0.3727 | -0.8819 | 0.3778 | 0.7199 | 1.3891 | Tidak signifikan |
| nurseryyes | pernah ikut taman kanak-kanak: yes dibanding acuan | nurseryyes | -0.1891 | 0.2233 | -0.8469 | 0.3970 | 0.8277 | 1.2082 | Tidak signifikan |
| Mjobservices | pekerjaan ibu: services dibanding kategori acuan | Mjobservices | 0.2474 | 0.3045 | 0.8126 | 0.4165 | 1.2807 | 0.7808 | Tidak signifikan |
| romanticyes | memiliki hubungan romantis: yes dibanding acuan | romanticyes | -0.1198 | 0.1876 | -0.6386 | 0.5231 | 0.8871 | 1.1273 | Tidak signifikan |
| goout | frekuensi pergi bermain bersama teman | goout | -0.0560 | 0.0885 | -0.6329 | 0.5268 | 0.9456 | 1.0576 | Tidak signifikan |
Pada model ordinal dengan MASS::polr(), nilai
Odds_Ratio_Higher_Perf = exp(beta) digunakan untuk membaca
kecenderungan menuju performa yang lebih tinggi. Karena performa
diurutkan sebagai Rendah < Sedang < Tinggi, maka OR di atas 1
berarti variabel tersebut berkaitan dengan performa akademik yang lebih
baik, sedangkan OR di bawah 1 berarti berkaitan dengan performa yang
lebih rendah.
pred_ordinal <- predict(model_ordinal, newdata = data_ordinal, type = "class")
perf_levels <- levels(data_ordinal$Performa_G3)
pred_ordinal_f <- factor(as.character(pred_ordinal), levels = perf_levels, ordered = TRUE)
aktual_ordinal_f <- factor(as.character(data_ordinal$Performa_G3), levels = perf_levels, ordered = TRUE)
conf_ordinal <- table(Prediksi = pred_ordinal_f, Aktual = aktual_ordinal_f)
conf_ordinal## Aktual
## Prediksi Rendah Sedang Tinggi
## Rendah 30 19 0
## Sedang 69 383 105
## Tinggi 1 16 26
akurasi_ordinal <- mean(as.character(pred_ordinal_f) == as.character(aktual_ordinal_f), na.rm = TRUE)
akurasi_ordinal## [1] 0.6764253
## [1] 1009.908
prob_ordinal <- predict(model_ordinal, type = "probs")
knitr::kable(
head(prob_ordinal, 6),
digits = 4,
caption = "Cuplikan Prediksi Probabilitas untuk Model Ordinal"
)| Rendah | Sedang | Tinggi |
|---|---|---|
| 0.0334 | 0.6648 | 0.3018 |
| 0.0252 | 0.6087 | 0.3661 |
| 0.2191 | 0.7304 | 0.0505 |
| 0.0318 | 0.6556 | 0.3126 |
| 0.0294 | 0.6409 | 0.3296 |
| 0.0284 | 0.6340 | 0.3375 |
Visualisasi ini digunakan untuk melihat bagaimana peluang
masing-masing kategori performa berubah ketika absences
(jumlah absensi) meningkat. Variabel lain ditahan pada nilai median atau
kategori paling umum.
grid_ordinal <- data.frame(
absences = seq(
quantile(data_ordinal$absences, 0.05, na.rm = TRUE),
quantile(data_ordinal$absences, 0.95, na.rm = TRUE),
length.out = 100
),
school = factor(modus(data_ordinal$school), levels = levels(data_ordinal$school)),
sex = factor(modus(data_ordinal$sex), levels = levels(data_ordinal$sex)),
age = median(data_ordinal$age, na.rm = TRUE),
address = factor(modus(data_ordinal$address), levels = levels(data_ordinal$address)),
famsize = factor(modus(data_ordinal$famsize), levels = levels(data_ordinal$famsize)),
Pstatus = factor(modus(data_ordinal$Pstatus), levels = levels(data_ordinal$Pstatus)),
Medu = median(data_ordinal$Medu, na.rm = TRUE),
Fedu = median(data_ordinal$Fedu, na.rm = TRUE),
Mjob = factor(modus(data_ordinal$Mjob), levels = levels(data_ordinal$Mjob)),
Fjob = factor(modus(data_ordinal$Fjob), levels = levels(data_ordinal$Fjob)),
reason = factor(modus(data_ordinal$reason), levels = levels(data_ordinal$reason)),
guardian = factor(modus(data_ordinal$guardian),levels = levels(data_ordinal$guardian)),
traveltime = median(data_ordinal$traveltime, na.rm = TRUE),
studytime = median(data_ordinal$studytime, na.rm = TRUE),
failures = 0,
schoolsup = factor(modus(data_ordinal$schoolsup), levels = levels(data_ordinal$schoolsup)),
famsup = factor(modus(data_ordinal$famsup), levels = levels(data_ordinal$famsup)),
paid = factor(modus(data_ordinal$paid), levels = levels(data_ordinal$paid)),
activities = factor(modus(data_ordinal$activities),levels = levels(data_ordinal$activities)),
nursery = factor(modus(data_ordinal$nursery), levels = levels(data_ordinal$nursery)),
higher = factor(modus(data_ordinal$higher), levels = levels(data_ordinal$higher)),
internet = factor(modus(data_ordinal$internet), levels = levels(data_ordinal$internet)),
romantic = factor(modus(data_ordinal$romantic), levels = levels(data_ordinal$romantic)),
famrel = median(data_ordinal$famrel, na.rm = TRUE),
freetime = median(data_ordinal$freetime,na.rm = TRUE),
goout = median(data_ordinal$goout, na.rm = TRUE),
Dalc = median(data_ordinal$Dalc, na.rm = TRUE),
Walc = median(data_ordinal$Walc, na.rm = TRUE),
health = median(data_ordinal$health, na.rm = TRUE)
)
prob_grid_ordinal <- predict(model_ordinal, newdata = grid_ordinal, type = "probs")
plot_prob_ordinal <- grid_ordinal %>%
dplyr::bind_cols(as.data.frame(prob_grid_ordinal)) %>%
tidyr::pivot_longer(
cols = levels(data_ordinal$Performa_G3),
names_to = "Performa_G3",
values_to = "Probabilitas"
) %>%
dplyr::mutate(
Performa_G3 = factor(Performa_G3, levels = levels(data_ordinal$Performa_G3))
)
ggplot(
plot_prob_ordinal,
aes(x = absences, y = Probabilitas, color = Performa_G3)
) +
geom_line(linewidth = 1.1) +
scale_color_manual(values = c("Rendah" = "#b84f5a", "Sedang" = "#d18b2f", "Tinggi" = "#2f7f73")) +
scale_y_continuous(labels = scales::percent_format()) +
labs(
title = "Prediksi Probabilitas Performa G3 Berdasarkan Jumlah Absensi",
subtitle = "Variabel lain ditahan pada median atau kategori paling umum",
x = "Jumlah absensi",
y = "Probabilitas prediksi",
color = "Performa"
) +
theme_minimal()Grafik ini membantu melihat perubahan peluang kategori performa ketika jumlah absensi meningkat. Jika probabilitas kategori Rendah meningkat ketika absensi naik, maka model menunjukkan bahwa ketidakhadiran yang lebih banyak berkaitan dengan performa akademik yang lebih rendah.
Asumsi utama regresi logistik ordinal adalah proportional odds atau parallel lines. Artinya, pengaruh prediktor diasumsikan sama pada setiap batas kumulatif performa.
cek_proportional_odds <- tryCatch({
model_ordinal_clm <- ordinal::clm(
Performa_G3 ~ school + sex + age + address + famsize + Pstatus +
Medu + Fedu + Mjob + Fjob + reason + guardian +
traveltime + studytime + failures + schoolsup + famsup + paid +
activities + nursery + higher + internet + romantic +
famrel + freetime + goout + Dalc + Walc + health + absences,
data = data_ordinal,
link = "logit"
)
print(summary(model_ordinal_clm))
ordinal::nominal_test(model_ordinal_clm)
}, error = function(e) {
message("Pemeriksaan proportional odds dengan nominal_test() tidak berhasil dijalankan: ", e$message)
NULL
})## formula:
## Performa_G3 ~ school + sex + age + address + famsize + Pstatus + Medu + Fedu + Mjob + Fjob + reason + guardian + traveltime + studytime + failures + schoolsup + famsup + paid + activities + nursery + higher + internet + romantic + famrel + freetime + goout + Dalc + Walc + health + absences
## data: data_ordinal
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 649 -463.95 1009.91 6(0) 3.96e-07 2.8e+05
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## schoolMS -1.135542 0.232158 -4.891 0.0000010020 ***
## sexM -0.446020 0.204016 -2.186 0.028801 *
## age 0.280465 0.083800 3.347 0.000817 ***
## addressU 0.009618 0.215975 0.045 0.964478
## famsizeLE3 0.295612 0.200129 1.477 0.139647
## PstatusT -0.014650 0.286185 -0.051 0.959174
## Medu 0.172778 0.126160 1.370 0.170840
## Fedu 0.066286 0.114404 0.579 0.562314
## Mjobhealth 0.441937 0.436664 1.012 0.311503
## Mjobother 0.026205 0.247282 0.106 0.915605
## Mjobservices 0.247434 0.304506 0.813 0.416462
## Mjobteacher 0.218447 0.408628 0.535 0.592936
## Fjobhealth -0.603218 0.621898 -0.970 0.332065
## Fjobother -0.328801 0.372680 -0.882 0.377635
## Fjobservices -0.595932 0.391492 -1.522 0.127958
## Fjobteacher 0.134110 0.544577 0.246 0.805477
## reasonhome 0.087477 0.232255 0.377 0.706441
## reasonother -0.097894 0.304740 -0.321 0.748029
## reasonreputation 0.318791 0.242194 1.316 0.188087
## guardianmother -0.342690 0.215939 -1.587 0.112518
## guardianother -0.228479 0.432635 -0.528 0.597423
## traveltime 0.127357 0.132574 0.961 0.336729
## studytime 0.255167 0.113733 2.244 0.024860 *
## failures -1.024713 0.180801 -5.668 0.0000000145 ***
## schoolsupyes -1.076327 0.306494 -3.512 0.000445 ***
## famsupyes -0.087886 0.186174 -0.472 0.636880
## paidyes -0.619860 0.386861 -1.602 0.109094
## activitiesyes 0.087664 0.182290 0.481 0.630582
## nurseryyes -0.189084 0.223271 -0.847 0.397061
## higheryes 1.492255 0.317990 4.693 0.0000026953 ***
## internetyes 0.099105 0.228145 0.434 0.664003
## romanticyes -0.119837 0.187635 -0.639 0.523036
## famrel 0.098765 0.095286 1.037 0.299968
## freetime -0.026826 0.092897 -0.289 0.772752
## goout -0.055985 0.088454 -0.633 0.526778
## Dalc -0.142799 0.127321 -1.122 0.262046
## Walc -0.039907 0.097326 -0.410 0.681777
## health -0.110030 0.063713 -1.727 0.084176 .
## absences -0.059222 0.021173 -2.797 0.005158 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## Rendah|Sedang 2.772 1.625 1.706
## Sedang|Tinggi 6.977 1.653 4.222
Jika hasil nominal_test() menunjukkan p-value kecil pada
prediktor tertentu, maka terdapat indikasi bahwa efek prediktor tersebut
tidak sepenuhnya memenuhi asumsi proportional odds. Namun hasil ini
tidak boleh dibaca secara mekanis; interpretasi tetap perlu
mempertimbangkan tujuan analisis, teori pendidikan, ukuran sampel, dan
kesederhanaan model.
Jika asumsi proportional odds tidak terpenuhi, maka model ordinal
sederhana dapat menjadi terlalu kaku. Alternatif yang dapat
dipertimbangkan adalah partial proportional odds model, generalized
ordered logit model, atau regresi multinomial sebagai pembanding
non-ordinal. Dalam laporan ini, model MASS::polr() tetap
digunakan sebagai model utama karena tujuan analisis adalah memahami
penerapan regresi logistik ordinal pada data performa akademik.
sig_ordinal <- coef_ordinal %>%
dplyr::filter(Jenis == "Prediktor", p_value < 0.05) %>%
dplyr::arrange(p_value)
cat("Pada taraf signifikansi 5%, model ordinal menemukan ", nrow(sig_ordinal),
" prediktor yang berhubungan signifikan dengan performa akademik siswa.\n\n", sep = "")Pada taraf signifikansi 5%, model ordinal menemukan 8 prediktor yang berhubungan signifikan dengan performa akademik siswa.
if (nrow(sig_ordinal) > 0) {
cat("Prediktor paling penting berdasarkan p-value terkecil adalah:\n\n")
top_ordinal <- head(sig_ordinal, 8)
for (i in seq_len(nrow(top_ordinal))) {
row <- top_ordinal[i, ]
arah <- ifelse(row$Odds_Ratio_Higher_Perf > 1,
"meningkatkan kecenderungan siswa masuk ke performa yang lebih tinggi",
"menurunkan kecenderungan siswa masuk ke performa yang lebih tinggi")
cat("- **", row$Variabel_Mudah, "** memiliki OR kategori lebih tinggi = ",
fmt_num(row$Odds_Ratio_Higher_Perf, 3),
" dan p-value = ", fmt_p(row$p_value), ", sehingga variabel ini ", arah,
", dengan asumsi variabel lain tetap. Dalam parameterisasi `polr`, ini berasal dari model `logit(P(Y <= j)) = alpha_j - x beta`.\n", sep = "")
}
cat("\n")
}Prediktor paling penting berdasarkan p-value terkecil adalah:
polr, ini berasal dari model
logit(P(Y <= j)) = alpha_j - x beta.polr, ini berasal dari model
logit(P(Y <= j)) = alpha_j - x beta.polr, ini berasal dari model
logit(P(Y <= j)) = alpha_j - x beta.polr, ini berasal dari model
logit(P(Y <= j)) = alpha_j - x beta.polr, ini berasal
dari model logit(P(Y <= j)) = alpha_j - x beta.polr, ini berasal dari model
logit(P(Y <= j)) = alpha_j - x beta.polr, ini berasal dari model
logit(P(Y <= j)) = alpha_j - x beta.polr, ini berasal dari model
logit(P(Y <= j)) = alpha_j - x beta.cat("Akurasi klasifikasi model ordinal adalah **", scales::percent(akurasi_ordinal),
"**. AIC model adalah **", fmt_num(AIC_ordinal, 2), "**.\n\n", sep = "")Akurasi klasifikasi model ordinal adalah 68%. AIC model adalah 1,009.91.
cat("**Kesimpulan substantif:** model ordinal menunjukkan faktor-faktor yang membedakan siswa berkinerja rendah dan tinggi. Variabel seperti jumlah kegagalan kelas sebelumnya, waktu belajar, dan konsumsi alkohol menjadi faktor penting dalam menentukan tingkat capaian akademik. Karena model menggunakan `MASS::polr`, OR di atas 1 berarti variabel tersebut mendorong siswa ke performa yang lebih tinggi.\n")Kesimpulan substantif: model ordinal menunjukkan
faktor-faktor yang membedakan siswa berkinerja rendah dan tinggi.
Variabel seperti jumlah kegagalan kelas sebelumnya, waktu belajar, dan
konsumsi alkohol menjadi faktor penting dalam menentukan tingkat capaian
akademik. Karena model menggunakan MASS::polr, OR di atas 1
berarti variabel tersebut mendorong siswa ke performa yang lebih
tinggi.
Tujuan analisis ini adalah mengetahui faktor-faktor yang berhubungan
dengan peran jabatan karyawan (JobRole) di
perusahaan. Peran jabatan tidak diperlakukan sebagai urutan, karena
kategori seperti Sales Executive,
Research Scientist, Manager, dan lainnya
adalah jenis pekerjaan yang berbeda.
Variabel respon adalah JobRole:
| Kategori JobRole | Arti sederhana |
|---|---|
| Healthcare Representative | Perwakilan layanan kesehatan |
| Human Resources | Staf SDM/HR |
| Laboratory Technician | Teknisi laboratorium |
| Manager | Manajer |
| Manufacturing Director | Direktur manufaktur |
| Research Director | Direktur riset |
| Research Scientist | Ilmuwan riset |
| Sales Executive | Eksekutif penjualan (acuan) |
| Sales Representative | Perwakilan penjualan |
| Variabel | Peran | Alasan_dipakai |
|---|---|---|
| JobRole | Respon | Mewakili peran jabatan karyawan (multinomial tanpa urutan) |
| Age | Prediktor | Usia karyawan berpengaruh pada jenis peran yang diemban |
| MonthlyIncome | Prediktor | Pendapatan bulanan mencerminkan level jabatan |
| TotalWorkingYears | Prediktor | Total pengalaman kerja mempengaruhi penempatan jabatan |
| YearsAtCompany | Prediktor | Lama bekerja di perusahaan |
| YearsInCurrentRole | Prediktor | Lama memegang peran jabatan saat ini |
| Department | Prediktor | Departemen kerja berkaitan erat dengan peran jabatan |
| EducationField | Prediktor | Bidang pendidikan mempengaruhi penempatan jabatan |
| JobLevel | Prediktor | Level jabatan mencerminkan hierarki dan peran |
| Attrition | Prediktor | Status keluar karyawan mungkin berkaitan dengan peran |
| BusinessTravel | Prediktor | Frekuensi perjalanan bisnis bervariasi antar peran |
| OverTime | Prediktor | Bekerja lembur dapat bervariasi antar jenis jabatan |
| Gender | Prediktor | Jenis kelamin sebagai faktor demografi |
| MaritalStatus | Prediktor | Status pernikahan sebagai faktor demografi |
Variabel EmployeeNumber, EmployeeCount,
Over18, dan StandardHours tidak digunakan
karena bersifat konstan atau hanya merupakan identifier. Variabel
JobSatisfaction, EnvironmentSatisfaction,
WorkLifeBalance, dan JobInvolvement adalah
ukuran kepuasan yang mungkin dipengaruhi oleh peran jabatan, bukan
sebaliknya, sehingga perlu berhati-hati dalam interpretasi
kausalitas.
data_multinom <- data_multinom_raw %>%
dplyr::mutate(
JobRole = relevel(factor(JobRole), ref = "Sales Executive"),
Attrition = factor(Attrition),
BusinessTravel = factor(BusinessTravel),
Department = factor(Department),
EducationField = factor(EducationField),
Gender = factor(Gender),
MaritalStatus = factor(MaritalStatus),
OverTime = factor(OverTime)
) %>%
dplyr::select(
JobRole, Age, MonthlyIncome, TotalWorkingYears, YearsAtCompany,
YearsInCurrentRole, Department, EducationField, JobLevel,
Attrition, BusinessTravel, OverTime, Gender, MaritalStatus,
JobSatisfaction, WorkLifeBalance, EnvironmentSatisfaction,
DistanceFromHome, PercentSalaryHike, NumCompaniesWorked,
TrainingTimesLastYear, YearsSinceLastPromotion, YearsWithCurrManager
) %>%
tidyr::drop_na()
knitr::kable(
head(data_multinom %>%
dplyr::select(JobRole, Age, MonthlyIncome, Department,
EducationField, JobLevel, Attrition, OverTime), 6),
caption = "Cuplikan Data Multinomial Setelah Cleaning"
)| JobRole | Age | MonthlyIncome | Department | EducationField | JobLevel | Attrition | OverTime |
|---|---|---|---|---|---|---|---|
| Sales Executive | 41 | 5993 | Sales | Life Sciences | 2 | Yes | Yes |
| Research Scientist | 49 | 5130 | Research & Development | Life Sciences | 2 | No | No |
| Laboratory Technician | 37 | 2090 | Research & Development | Other | 1 | Yes | Yes |
| Research Scientist | 33 | 2909 | Research & Development | Life Sciences | 1 | No | Yes |
| Laboratory Technician | 27 | 3468 | Research & Development | Medical | 1 | No | No |
| Laboratory Technician | 32 | 3068 | Research & Development | Life Sciences | 1 | No | No |
Sales Executive digunakan sebagai kategori acuan. Jadi
hasil model akan membaca peran jabatan lain dibandingkan dengan Sales
Executive.
dist_multinom <- data_multinom %>%
dplyr::count(JobRole) %>%
dplyr::mutate(Persen = scales::percent(n / sum(n))) %>%
dplyr::arrange(desc(n))
knitr::kable(dist_multinom, caption = "Distribusi Peran Jabatan Karyawan (JobRole)")| JobRole | n | Persen |
|---|---|---|
| Sales Executive | 326 | 22.18% |
| Research Scientist | 292 | 19.86% |
| Laboratory Technician | 259 | 17.62% |
| Manufacturing Director | 145 | 9.86% |
| Healthcare Representative | 131 | 8.91% |
| Manager | 102 | 6.94% |
| Sales Representative | 83 | 5.65% |
| Research Director | 80 | 5.44% |
| Human Resources | 52 | 3.54% |
ringkasan_multinom <- data_multinom %>%
dplyr::summarise(
Rata_rata_Age = mean(Age),
Median_Age = median(Age),
Rata_rata_Income = mean(MonthlyIncome),
Median_Income = median(MonthlyIncome),
Rata_rata_WorkingYears = mean(TotalWorkingYears)
)
knitr::kable(ringkasan_multinom, caption = "Ringkasan Prediktor Utama Data Multinomial")| Rata_rata_Age | Median_Age | Rata_rata_Income | Median_Income | Rata_rata_WorkingYears |
|---|---|---|---|---|
| 36.92381 | 36 | 6502.931 | 4919 | 11.27959 |
ggplot(data_multinom, aes(x = forcats::fct_infreq(JobRole))) +
geom_bar(fill = "#2f7f73") +
coord_flip() +
labs(
title = "Distribusi Peran Jabatan Karyawan (JobRole)",
x = "Job Role",
y = "Jumlah karyawan"
) +
theme_minimal()Model multinomial membandingkan setiap kategori JobRole
terhadap kategori acuan Sales Executive. Model membangun
persamaan logit untuk setiap kategori sekaligus.
Untuk setiap kategori jabatan \(k\)
selain Sales Executive, bentuk modelnya adalah:
\[ \log\left( \frac{P(Y_i = k)}{P(Y_i = \text{Sales Executive})} \right) = \beta_{0k} + \beta_{1k}X_{1i} + \cdots + \beta_{pk}X_{pi} \]
model_multinom <- nnet::multinom(
JobRole ~ Age + MonthlyIncome + TotalWorkingYears + YearsAtCompany +
YearsInCurrentRole + Department + EducationField + JobLevel +
Attrition + BusinessTravel + OverTime + Gender + MaritalStatus +
JobSatisfaction + WorkLifeBalance + EnvironmentSatisfaction +
DistanceFromHome + PercentSalaryHike + NumCompaniesWorked +
TrainingTimesLastYear + YearsSinceLastPromotion + YearsWithCurrManager,
data = data_multinom,
trace = FALSE,
maxit = 300,
MaxNWts= 15000
)
summary(model_multinom)## Call:
## nnet::multinom(formula = JobRole ~ Age + MonthlyIncome + TotalWorkingYears +
## YearsAtCompany + YearsInCurrentRole + Department + EducationField +
## JobLevel + Attrition + BusinessTravel + OverTime + Gender +
## MaritalStatus + JobSatisfaction + WorkLifeBalance + EnvironmentSatisfaction +
## DistanceFromHome + PercentSalaryHike + NumCompaniesWorked +
## TrainingTimesLastYear + YearsSinceLastPromotion + YearsWithCurrManager,
## data = data_multinom, trace = FALSE, maxit = 300, MaxNWts = 15000)
##
## Coefficients:
## (Intercept) Age MonthlyIncome
## Healthcare Representative -193.7347 8.9739256 -0.1375353
## Human Resources 1234.7691 2.6355715 -0.1999099
## Laboratory Technician -485.0004 8.9783870 -0.1385331
## Manager -656.0151 -13.5323613 0.3189753
## Manufacturing Director -276.9259 8.9579585 -0.1376054
## Research Director -1170.9908 -13.5773081 0.3188088
## Research Scientist -286.0277 8.9775036 -0.1385099
## Sales Representative 1032.0450 0.1515729 -0.2168137
## TotalWorkingYears YearsAtCompany YearsInCurrentRole
## Healthcare Representative 2.459191 7.424926 -10.194318
## Human Resources 3.452041 -11.170283 17.246048
## Laboratory Technician 2.577644 7.519987 -10.094700
## Manager -11.984791 -12.395691 -3.372370
## Manufacturing Director 2.451742 7.359127 -10.130969
## Research Director -11.928672 -12.435679 -3.257568
## Research Scientist 2.586693 7.544561 -10.096615
## Sales Representative 36.454061 -3.872958 -23.728702
## DepartmentResearch & Development DepartmentSales
## Healthcare Representative 795.29932 -430.0472
## Human Resources -891.43044 -2018.5754
## Laboratory Technician 1137.37243 -1001.6746
## Manager -743.23921 -1165.0631
## Manufacturing Director 915.03478 -463.8049
## Research Director -68.55588 -553.8204
## Research Scientist 1001.24888 -1434.9438
## Sales Representative 186.13608 336.0655
## EducationFieldLife Sciences EducationFieldMarketing
## Healthcare Representative -70.699573 -123.41860
## Human Resources 5.413813 209.81586
## Laboratory Technician -112.121333 -261.14018
## Manager -152.376589 34.05217
## Manufacturing Director -107.988059 28.97191
## Research Director -308.389921 -224.70220
## Research Scientist -174.381050 -271.15342
## Sales Representative -28.628892 -127.10677
## EducationFieldMedical EducationFieldOther
## Healthcare Representative -68.16040 -50.99249
## Human Resources -61.60781 -204.35700
## Laboratory Technician -109.97361 -91.98464
## Manager -111.96348 -95.97247
## Manufacturing Director -105.44327 -88.85988
## Research Director -268.21754 -253.59824
## Research Scientist -172.29094 -154.55773
## Sales Representative 61.02722 -112.47609
## EducationFieldTechnical Degree JobLevel AttritionYes
## Healthcare Representative 264.86632 250.3995 -58.47558
## Human Resources 251.35120 351.2665 -163.20997
## Laboratory Technician 222.86984 246.2751 -57.05670
## Manager 81.70314 -354.3686 -269.12643
## Manufacturing Director 227.37777 250.8466 -58.19880
## Research Director -73.69057 -354.3405 -271.02778
## Research Scientist 161.48105 245.8634 -58.11508
## Sales Representative 590.51387 -573.2649 -257.84277
## BusinessTravelTravel_Frequently
## Healthcare Representative 31.66026
## Human Resources 51.25161
## Laboratory Technician 32.34446
## Manager -193.06969
## Manufacturing Director 32.06820
## Research Director -193.13917
## Research Scientist 32.45998
## Sales Representative 193.24043
## BusinessTravelTravel_Rarely OverTimeYes GenderMale
## Healthcare Representative 78.97647 -22.28982 -31.89363
## Human Resources 66.81188 11.24389 -49.84380
## Laboratory Technician 79.25857 -22.45033 -31.58060
## Manager -99.05467 -66.01441 -107.02564
## Manufacturing Director 79.30070 -22.31203 -32.36521
## Research Director -98.08181 -65.09154 -106.29600
## Research Scientist 79.48062 -21.66776 -31.87298
## Sales Representative 329.33854 -116.30310 -16.47816
## MaritalStatusMarried MaritalStatusSingle
## Healthcare Representative 224.92532 191.99427
## Human Resources -10.40213 80.91409
## Laboratory Technician 225.16249 191.98459
## Manager 74.28209 39.81838
## Manufacturing Director 224.99214 191.99838
## Research Director 73.90523 39.80470
## Research Scientist 225.03891 192.22898
## Sales Representative 75.43468 184.90523
## JobSatisfaction WorkLifeBalance
## Healthcare Representative -13.34660 59.31338
## Human Resources -49.94713 137.59045
## Laboratory Technician -13.39397 59.59713
## Manager 87.82190 -14.02239
## Manufacturing Director -13.39672 59.47529
## Research Director 87.83321 -14.02890
## Research Scientist -13.32346 59.50167
## Sales Representative 33.94052 -32.11259
## EnvironmentSatisfaction DistanceFromHome
## Healthcare Representative -54.55074 -7.0213057
## Human Resources -58.08572 -4.3066132
## Laboratory Technician -54.57459 -7.0190072
## Manager -83.74942 -3.1874636
## Manufacturing Director -54.39896 -7.0271722
## Research Director -84.06256 -3.1434318
## Research Scientist -54.64904 -7.0238500
## Sales Representative -65.51621 0.3767013
## PercentSalaryHike NumCompaniesWorked
## Healthcare Representative 18.96011 -11.061296
## Human Resources 10.82848 -29.396468
## Laboratory Technician 18.93470 -10.932305
## Manager 15.68289 -6.823188
## Manufacturing Director 18.97728 -11.072696
## Research Director 15.66871 -6.623813
## Research Scientist 18.97797 -10.949900
## Sales Representative 17.49796 -24.087189
## TrainingTimesLastYear YearsSinceLastPromotion
## Healthcare Representative -19.736312 -7.350464
## Human Resources -14.094838 7.578287
## Laboratory Technician -19.646913 -7.453880
## Manager 5.955875 22.271297
## Manufacturing Director -19.732722 -7.447505
## Research Director 5.695241 22.118682
## Research Scientist -19.827761 -7.416186
## Sales Representative -6.406207 -0.646007
## YearsWithCurrManager
## Healthcare Representative -28.76618
## Human Resources -13.22622
## Laboratory Technician -28.76197
## Manager -24.60483
## Manufacturing Director -28.66626
## Research Director -24.54373
## Research Scientist -28.80312
## Sales Representative -19.70450
##
## Std. Errors:
## (Intercept) Age MonthlyIncome
## Healthcare Representative 0.001636025530333 0.0119006104465 0.00007059469
## Human Resources NaN NaN NaN
## Laboratory Technician 0.002070922882813 0.0105008697927 0.00007596846
## Manager 0.000587173384428 0.0183056492320 0.00004132232
## Manufacturing Director 0.001700838859934 0.0117261630003 0.00006914167
## Research Director 0.000587173311970 0.0183056440108 0.00004147625
## Research Scientist 0.002102577166302 0.0105565178262 0.00007585756
## Sales Representative 0.000000004943403 0.0000002323387 0.00001471157
## TotalWorkingYears YearsAtCompany
## Healthcare Representative 0.02387362661423726023502 0.02537669001318218
## Human Resources 0.00000000000000002074766 NaN
## Laboratory Technician 0.02418860892282232583073 0.03122426517284106
## Manager 0.02472947126574341278338 0.02276377874454226
## Manufacturing Director 0.02382080419751222283287 0.02493787863774681
## Research Director 0.02472947256476230870614 0.02276377839226153
## Research Scientist 0.02393260740736621178248 0.03074529372671153
## Sales Representative 0.00000002471866503520047 0.00000000001675842
## YearsInCurrentRole
## Healthcare Representative 0.03052644707925256453707
## Human Resources 0.00000000000000001865721
## Laboratory Technician 0.03156598688980761374667
## Manager 0.03357632803474904642504
## Manufacturing Director 0.03015770760929062030420
## Research Director 0.03357632818128254009737
## Research Scientist 0.03126323567955484328262
## Sales Representative 0.00000000000651716210437
## DepartmentResearch & Development
## Healthcare Representative 0.001636029183405088
## Human Resources NaN
## Laboratory Technician 0.002070922882813129
## Manager 0.000587173311970338
## Manufacturing Director 0.001700838859483119
## Research Director 0.000587173311970375
## Research Scientist 0.002102577166302162
## Sales Representative 0.000000000000931023
## DepartmentSales
## Healthcare Representative 0.00000000638827162086904348390833941
## Human Resources NaN
## Laboratory Technician 0.00000000000000000000020519128635730
## Manager 0.00000000010732756995771270353088461
## Manufacturing Director 0.00000000000263920206070341407727717
## Research Director 0.00000000000000000000000000001813052
## Research Scientist 0.00000000000000000000000000000000000
## Sales Representative 0.00000000494327645416542970930456891
## EducationFieldLife Sciences
## Healthcare Representative 0.0008845428188583855
## Human Resources NaN
## Laboratory Technician 0.0025405656761900621
## Manager 0.0012534732492635056
## Manufacturing Director 0.0010222772619373137
## Research Director 0.0012534732492634698
## Research Scientist 0.0026764714699485884
## Sales Representative 0.0000000000009310232
## EducationFieldMarketing
## Healthcare Representative 0.00000000001094151366904666682396823906842087126278784126043319702148437500000000
## Human Resources 0.00000000000000000477413122486281363537202748759114001586567610502243041992187500
## Laboratory Technician NaN
## Manager 0.00000000000000000000001241511571628525825109617963803998463845346122980117797852
## Manufacturing Director 0.00000000000263920206070341407727716798170547463087132200598716735839843750000000
## Research Director 0.00000000000000000000000000000000000000000000000000000000000000000000000003822273
## Research Scientist 0.00000000000000000000000000000000000000000000000000000000000000000000000000000000
## Sales Representative 0.00000000000000000000265061710532166948601009992536603476764867082238197326660156
## EducationFieldMedical
## Healthcare Representative 0.001560525154756
## Human Resources NaN
## Laboratory Technician 0.001867098412612
## Manager 0.001817337771521
## Manufacturing Director 0.001658248709124
## Research Director 0.001817337771521
## Research Scientist 0.001904897865753
## Sales Representative 0.000000004943276
## EducationFieldOther
## Healthcare Representative 0.000504691717260775983902587604745804128469899296760559082031250000000000000000000000000000000
## Human Resources NaN
## Laboratory Technician 0.001063163402480520047799217309147934429347515106201171875000000000000000000000000000000000000
## Manager 0.000742123304667742229694116229410383311915211379528045654296875000000000000000000000000000000
## Manufacturing Director 0.000415907051641658371231302648851624326198361814022064208984375000000000000000000000000000000
## Research Director 0.000742123304667742446534550726511270113405771553516387939453125000000000000000000000000000000
## Research Scientist 0.001070603243568866297374952978316287044435739517211914062500000000000000000000000000000000000
## Sales Representative 0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000002517525
## EducationFieldTechnical Degree
## Healthcare Representative 0.00064060429485782621923662
## Human Resources NaN
## Laboratory Technician 0.00134278618064070085803829
## Manager 0.00092682365934712357055658
## Manufacturing Director 0.00059318261718578808600749
## Research Director 0.00092682360980200762228670
## Research Scientist 0.00144881331639602740045392
## Sales Representative 0.00000000000000000004873736
## JobLevel
## Healthcare Representative 0.00205699073121571500769
## Human Resources 0.00000000000000002914509
## Laboratory Technician 0.00151841061368537789088
## Manager 0.00070455836363053875142
## Manufacturing Director 0.00215055573440868248031
## Research Director 0.00070455824250549830678
## Research Scientist 0.00162570040323980363230
## Sales Representative 0.00000000988667945900040
## AttritionYes
## Healthcare Representative 0.00075810334407961472383969
## Human Resources NaN
## Laboratory Technician 0.00199288243452218610113347
## Manager 0.00044450665026605817907804
## Manufacturing Director 0.00075284637701391839744586
## Research Director 0.00044450665026605817907804
## Research Scientist 0.00176686046629414617893661
## Sales Representative 0.00000000000000000000276749
## BusinessTravelTravel_Frequently
## Healthcare Representative 0.001016306935012
## Human Resources NaN
## Laboratory Technician 0.001493608623050
## Manager 0.000916612847806
## Manufacturing Director 0.000874829397026
## Research Director 0.000916612847806
## Research Scientist 0.001410495145240
## Sales Representative 0.000000004943276
## BusinessTravelTravel_Rarely
## Healthcare Representative 0.001628706137962928788762
## Human Resources 0.000000000000000001088456
## Laboratory Technician 0.002052301019080409860001
## Manager 0.000568828934115224145798
## Manufacturing Director 0.001635578296240736313749
## Research Director 0.000568828850159129898942
## Research Scientist 0.002004366972167395243715
## Sales Representative 0.000000000000931023052908
## OverTimeYes
## Healthcare Representative 0.0006797489045432548440336951
## Human Resources NaN
## Laboratory Technician 0.0021428841740646622712751501
## Manager 0.0019207375782153575903554721
## Manufacturing Director 0.0005678841930602837421993012
## Research Director 0.0019207374879906787337735530
## Research Scientist 0.0022334226167275263510991756
## Sales Representative 0.0000000000000000000003393096
## GenderMale
## Healthcare Representative 0.0010257715808876994613
## Human Resources 0.0000000000000000477711
## Laboratory Technician 0.0041382449293419498737
## Manager 0.0012302687527956453380
## Manufacturing Director 0.0008377975329476091554
## Research Director 0.0012302687573596303786
## Research Scientist 0.0041032135366639117466
## Sales Representative 0.0000000049432764535056
## MaritalStatusMarried
## Healthcare Representative 0.0015205226581476033193701669
## Human Resources NaN
## Laboratory Technician 0.0016468321399982083219626361
## Manager 0.0018976766897374022259603610
## Manufacturing Director 0.0015410703304662427652604872
## Research Director 0.0018976766927394060364281403
## Research Scientist 0.0019902173909439784330044354
## Sales Representative 0.0000000000000000000003393096
## MaritalStatusSingle JobSatisfaction
## Healthcare Representative 0.001411449394486838761556 0.00893059191295
## Human Resources 0.000000000000000009002845 NaN
## Laboratory Technician 0.002684128858123998109808 0.03791829075527
## Manager 0.001200400435271526625863 0.00504277495354
## Manufacturing Director 0.001664901327257130780532 0.00874554171666
## Research Director 0.001200400435271549177269 0.00504277486374
## Research Scientist 0.002861877633423594909784 0.03782582266990
## Sales Representative 0.000000004943276453077128 0.00000001482996
## WorkLifeBalance EnvironmentSatisfaction
## Healthcare Representative 0.00547548034186 0.01139910817194
## Human Resources NaN NaN
## Laboratory Technician 0.00845676629116 0.03699610529655
## Manager 0.00224101470282 0.00207440467826
## Manufacturing Director 0.00571756094418 0.01166427156443
## Research Director 0.00224101470289 0.00207440478203
## Research Scientist 0.00862005858322 0.03700484664013
## Sales Representative 0.00000001483008 0.00000000988706
## DistanceFromHome
## Healthcare Representative 0.0110799574626983237529698
## Human Resources 0.0000000000000000009644243
## Laboratory Technician 0.0098012203090086895407307
## Manager 0.0139714346613242224037776
## Manufacturing Director 0.0109086601371174169577971
## Research Director 0.0139714346914975948849014
## Research Scientist 0.0098925655351133935655472
## Sales Representative 0.0000000197748811811796246
## PercentSalaryHike NumCompaniesWorked
## Healthcare Representative 0.022842101846402788206669 0.04092498436839
## Human Resources 0.000000000000000007724598 NaN
## Laboratory Technician 0.019848978120084289494818 0.03579188579651
## Manager 0.023626888478305234364552 0.04508776702955
## Manufacturing Director 0.021688991152221732960914 0.03993711255418
## Research Director 0.023626888821036241133955 0.04508776708183
## Research Scientist 0.019642981068175021514000 0.03582250349854
## Sales Representative 0.000000093923771402119989 0.00000001482996
## TrainingTimesLastYear YearsSinceLastPromotion
## Healthcare Representative 0.03941278442076 0.03856968680186983
## Human Resources NaN NaN
## Laboratory Technician 0.03289107784399 0.03863004946199855
## Manager 0.00460650650370 0.03232816839753523
## Manufacturing Director 0.04056147162006 0.03898866456054612
## Research Director 0.00460650644860 0.03232816839753531
## Research Scientist 0.03238756166211 0.03820453001122903
## Sales Representative 0.00000001483008 0.00000000001117228
## YearsWithCurrManager
## Healthcare Representative 0.03144712159051799
## Human Resources NaN
## Laboratory Technician 0.03141118633364032
## Manager 0.03645186601656455
## Manufacturing Director 0.03074123837897624
## Research Director 0.03645186611473998
## Research Scientist 0.03172568501362495
## Sales Representative 0.00000000001582739
##
## Residual Deviance: 1541.73
## AIC: 2021.73
Chunk model-multinomial di atas sudah membangun
model regresi logistik multinomial menggunakan
nnet::multinom(). Persamaan logit yang dibentuk adalah
untuk setiap kategori JobRole dibandingkan Sales Executive,
misalnya:
\[ \log\left(\frac{P(Y=\mathrm{Manager})}{P(Y=\mathrm{Sales\ Executive})}\right) = \beta_{0,\mathrm{Manager}} + \beta_{1,\mathrm{Manager}}\,\mathrm{Age} + \cdots + \beta_{p,\mathrm{Manager}}\,X_p \]
summary_multinom <- summary(model_multinom)
coef_multinom <- summary_multinom$coefficients
se_multinom <- summary_multinom$standard.errors
multinom_table <- purrr::map_dfr(rownames(coef_multinom), function(kelas) {
tibble::tibble(
Kategori = kelas,
term = colnames(coef_multinom),
estimate = as.numeric(coef_multinom[kelas, ]),
std_error = as.numeric(se_multinom[kelas, ])
)
}) %>%
dplyr::mutate(
z_value = estimate / std_error,
p_value = 2 * pnorm(abs(z_value), lower.tail = FALSE),
Relative_Risk_Ratio = exp(estimate),
Variabel_Mudah = purrr::map_chr(term, label_multinom),
Keputusan_5persen = ifelse(p_value < 0.05, "Signifikan", "Tidak signifikan")
) %>%
dplyr::arrange(Kategori, p_value)
knitr::kable(
multinom_table %>%
dplyr::select(Kategori, Variabel_Mudah, term, estimate, std_error,
z_value, p_value, Relative_Risk_Ratio, Keputusan_5persen) %>%
head(40),
digits = 4,
caption = "Koefisien dan Relative Risk Ratio Regresi Multinomial"
)| Kategori | Variabel_Mudah | term | estimate | std_error | z_value | p_value | Relative_Risk_Ratio | Keputusan_5persen |
|---|---|---|---|---|---|---|---|---|
| Healthcare Representative | (Intercept) | (Intercept) | -193.7347 | 0.0016 | -118417.8912 | 0 | 0.0000 | Signifikan |
| Healthcare Representative | usia karyawan | Age | 8.9739 | 0.0119 | 754.0727 | 0 | 7894.5316 | Signifikan |
| Healthcare Representative | pendapatan bulanan karyawan | MonthlyIncome | -0.1375 | 0.0001 | -1948.2385 | 0 | 0.8715 | Signifikan |
| Healthcare Representative | total masa kerja (tahun) | TotalWorkingYears | 2.4592 | 0.0239 | 103.0087 | 0 | 11.6954 | Signifikan |
| Healthcare Representative | lama bekerja di perusahaan saat ini | YearsAtCompany | 7.4249 | 0.0254 | 292.5884 | 0 | 1677.2749 | Signifikan |
| Healthcare Representative | lama memegang peran jabatan saat ini | YearsInCurrentRole | -10.1943 | 0.0305 | -333.9504 | 0 | 0.0000 | Signifikan |
| Healthcare Representative | departemen kerja: Research & Development dibanding acuan | DepartmentResearch & Development | 795.2993 | 0.0016 | 486115.6060 | 0 | Inf | Signifikan |
| Healthcare Representative | departemen kerja: Sales dibanding acuan | DepartmentSales | -430.0472 | 0.0000 | -67318234280.1656 | 0 | 0.0000 | Signifikan |
| Healthcare Representative | bidang pendidikan: Life Sciences dibanding acuan | EducationFieldLife Sciences | -70.6996 | 0.0009 | -79927.8130 | 0 | 0.0000 | Signifikan |
| Healthcare Representative | bidang pendidikan: Marketing dibanding acuan | EducationFieldMarketing | -123.4186 | 0.0000 | -11279847252144.3848 | 0 | 0.0000 | Signifikan |
| Healthcare Representative | bidang pendidikan: Medical dibanding acuan | EducationFieldMedical | -68.1604 | 0.0016 | -43677.8606 | 0 | 0.0000 | Signifikan |
| Healthcare Representative | bidang pendidikan: Other dibanding acuan | EducationFieldOther | -50.9925 | 0.0005 | -101036.9016 | 0 | 0.0000 | Signifikan |
| Healthcare Representative | bidang pendidikan: Technical Degree dibanding acuan | EducationFieldTechnical Degree | 264.8663 | 0.0006 | 413463.2277 | 0 | 10714723834898563653422080068024864426208660466464480806626082860882400662820204660080062440040420860606468484480606.0000 | Signifikan |
| Healthcare Representative | level jabatan karyawan | JobLevel | 250.3995 | 0.0021 | 121730.9990 | 0 | 5586464935741923457206606602042888460440244608880482264006488082844064864246288682468880042606224002004862202.0000 | Signifikan |
| Healthcare Representative | status keluar karyawan: Yes dibanding acuan | AttritionYes | -58.4756 | 0.0008 | -77134.0468 | 0 | 0.0000 | Signifikan |
| Healthcare Representative | frekuensi perjalanan bisnis: Travel_Frequently dibanding acuan | BusinessTravelTravel_Frequently | 31.6603 | 0.0010 | 31152.2620 | 0 | 56218101026456.5156 | Signifikan |
| Healthcare Representative | frekuensi perjalanan bisnis: Travel_Rarely dibanding acuan | BusinessTravelTravel_Rarely | 78.9765 | 0.0016 | 48490.3152 | 0 | 19908882855787660328482624008408666.0000 | Signifikan |
| Healthcare Representative | bekerja lembur: Yes dibanding acuan | OverTimeYes | -22.2898 | 0.0007 | -32791.2601 | 0 | 0.0000 | Signifikan |
| Healthcare Representative | jenis kelamin karyawan: Male dibanding acuan | GenderMale | -31.8936 | 0.0010 | -31092.3280 | 0 | 0.0000 | Signifikan |
| Healthcare Representative | status pernikahan: Married dibanding acuan | MaritalStatusMarried | 224.9253 | 0.0015 | 147926.3188 | 0 | 48286430097978824912026266200400440486828886444864062862208846680668886802408260848040248406084860.0000 | Signifikan |
| Healthcare Representative | status pernikahan: Single dibanding acuan | MaritalStatusSingle | 191.9943 | 0.0014 | 136026.3191 | 0 | 241018371863804662582200004664686062268446842024062804624648808426848004600880282084.0000 | Signifikan |
| Healthcare Representative | tingkat kepuasan kerja (skala 1-4) | JobSatisfaction | -13.3466 | 0.0089 | -1494.4811 | 0 | 0.0000 | Signifikan |
| Healthcare Representative | keseimbangan hidup-kerja (skala 1-4) | WorkLifeBalance | 59.3134 | 0.0055 | 10832.5441 | 0 | 57474425417065050072440480.0000 | Signifikan |
| Healthcare Representative | kepuasan lingkungan kerja (skala 1-4) | EnvironmentSatisfaction | -54.5507 | 0.0114 | -4785.5269 | 0 | 0.0000 | Signifikan |
| Healthcare Representative | jarak dari rumah ke kantor (km) | DistanceFromHome | -7.0213 | 0.0111 | -633.6943 | 0 | 0.0009 | Signifikan |
| Healthcare Representative | persentase kenaikan gaji | PercentSalaryHike | 18.9601 | 0.0228 | 830.0509 | 0 | 171502285.1237 | Signifikan |
| Healthcare Representative | jumlah perusahaan sebelumnya | NumCompaniesWorked | -11.0613 | 0.0409 | -270.2822 | 0 | 0.0000 | Signifikan |
| Healthcare Representative | frekuensi pelatihan tahun lalu | TrainingTimesLastYear | -19.7363 | 0.0394 | -500.7592 | 0 | 0.0000 | Signifikan |
| Healthcare Representative | tahun sejak promosi terakhir | YearsSinceLastPromotion | -7.3505 | 0.0386 | -190.5762 | 0 | 0.0006 | Signifikan |
| Healthcare Representative | lama dengan manajer saat ini | YearsWithCurrManager | -28.7662 | 0.0314 | -914.7478 | 0 | 0.0000 | Signifikan |
| Human Resources | total masa kerja (tahun) | TotalWorkingYears | 3.4520 | 0.0000 | 166382212338542496.0000 | 0 | 31.5648 | Signifikan |
| Human Resources | lama memegang peran jabatan saat ini | YearsInCurrentRole | 17.2460 | 0.0000 | 924363539835938176.0000 | 0 | 30893240.0541 | Signifikan |
| Human Resources | bidang pendidikan: Marketing dibanding acuan | EducationFieldMarketing | 209.8159 | 0.0000 | 43948490017981497344.0000 | 0 | 13239438112204088689824262462024004266822026442288022684220462664422686806864624240464426282.0000 | Signifikan |
| Human Resources | level jabatan karyawan | JobLevel | 351.2665 | 0.0000 | 12052341435678699520.0000 | 0 | 357368189835990547028220622822246644268288082062806224040860606264640640042686446800226604084864026648284442686446428000060040028028806282624664840648288.0000 | Signifikan |
| Human Resources | frekuensi perjalanan bisnis: Travel_Rarely dibanding acuan | BusinessTravelTravel_Rarely | 66.8119 | 0.0000 | 61382254040884273152.0000 | 0 | 103759853300141322560488404226.0000 | Signifikan |
| Human Resources | jenis kelamin karyawan: Male dibanding acuan | GenderMale | -49.8438 | 0.0000 | -1043388109321112064.0000 | 0 | 0.0000 | Signifikan |
| Human Resources | status pernikahan: Single dibanding acuan | MaritalStatusSingle | 80.9141 | 0.0000 | 8987613563558907904.0000 | 0 | 138210813838098451120288882046268028.0000 | Signifikan |
| Human Resources | jarak dari rumah ke kantor (km) | DistanceFromHome | -4.3066 | 0.0000 | -4465475619272303616.0000 | 0 | 0.0135 | Signifikan |
| Human Resources | persentase kenaikan gaji | PercentSalaryHike | 10.8285 | 0.0000 | 1401818318499469312.0000 | 0 | 50437.1094 | Signifikan |
| Human Resources | (Intercept) | (Intercept) | 1234.7691 | NaN | NaN | NaN | Inf | NA |
pred_multinom <- predict(model_multinom, type = "class")
conf_multinom <- table(Prediksi = pred_multinom, Aktual = data_multinom$JobRole)
conf_multinom## Aktual
## Prediksi Sales Executive Healthcare Representative
## Sales Executive 326 0
## Healthcare Representative 0 65
## Human Resources 0 0
## Laboratory Technician 0 8
## Manager 0 0
## Manufacturing Director 0 51
## Research Director 0 0
## Research Scientist 0 7
## Sales Representative 0 0
## Aktual
## Prediksi Human Resources Laboratory Technician Manager
## Sales Executive 0 0 0
## Healthcare Representative 0 12 0
## Human Resources 52 0 0
## Laboratory Technician 0 132 0
## Manager 0 0 78
## Manufacturing Director 0 10 0
## Research Director 0 0 24
## Research Scientist 0 105 0
## Sales Representative 0 0 0
## Aktual
## Prediksi Manufacturing Director Research Director
## Sales Executive 0 0
## Healthcare Representative 37 0
## Human Resources 0 0
## Laboratory Technician 9 0
## Manager 0 14
## Manufacturing Director 89 0
## Research Director 0 66
## Research Scientist 10 0
## Sales Representative 0 0
## Aktual
## Prediksi Research Scientist Sales Representative
## Sales Executive 0 0
## Healthcare Representative 5 0
## Human Resources 0 0
## Laboratory Technician 75 0
## Manager 0 0
## Manufacturing Director 15 0
## Research Director 0 0
## Research Scientist 197 0
## Sales Representative 0 83
akurasi_multinom <- mean(as.character(pred_multinom) == as.character(data_multinom$JobRole))
akurasi_multinom## [1] 0.7401361
## [1] 2021.73
prob_multinom <- predict(model_multinom, type = "probs")
knitr::kable(
head(prob_multinom, 6),
digits = 4,
caption = "Cuplikan Prediksi Probabilitas untuk Model Multinomial"
)| Sales Executive | Healthcare Representative | Human Resources | Laboratory Technician | Manager | Manufacturing Director | Research Director | Research Scientist | Sales Representative |
|---|---|---|---|---|---|---|---|---|
| 1 | 0.0000 | 0 | 0.0000 | 0 | 0.0000 | 0 | 0.0000 | 0 |
| 0 | 0.1600 | 0 | 0.3373 | 0 | 0.2671 | 0 | 0.2357 | 0 |
| 0 | 0.0003 | 0 | 0.6314 | 0 | 0.0002 | 0 | 0.3681 | 0 |
| 0 | 0.0007 | 0 | 0.2911 | 0 | 0.0007 | 0 | 0.7075 | 0 |
| 0 | 0.0034 | 0 | 0.5798 | 0 | 0.0022 | 0 | 0.4145 | 0 |
| 0 | 0.0017 | 0 | 0.4071 | 0 | 0.0019 | 0 | 0.5893 | 0 |
Visualisasi probabilitas prediksi digunakan untuk melihat bagaimana
peluang masing-masing peran jabatan berubah ketika
MonthlyIncome berubah. Variabel lain ditahan pada nilai
yang umum.
grid_multinom <- data.frame(
MonthlyIncome = seq(
quantile(data_multinom$MonthlyIncome, 0.05, na.rm = TRUE),
quantile(data_multinom$MonthlyIncome, 0.95, na.rm = TRUE),
length.out = 100
),
Age = median(data_multinom$Age, na.rm = TRUE),
TotalWorkingYears = median(data_multinom$TotalWorkingYears, na.rm = TRUE),
YearsAtCompany = median(data_multinom$YearsAtCompany, na.rm = TRUE),
YearsInCurrentRole = median(data_multinom$YearsInCurrentRole, na.rm = TRUE),
Department = factor(modus(data_multinom$Department), levels = levels(data_multinom$Department)),
EducationField = factor(modus(data_multinom$EducationField), levels = levels(data_multinom$EducationField)),
JobLevel = median(data_multinom$JobLevel, na.rm = TRUE),
Attrition = factor(modus(data_multinom$Attrition), levels = levels(data_multinom$Attrition)),
BusinessTravel = factor(modus(data_multinom$BusinessTravel), levels = levels(data_multinom$BusinessTravel)),
OverTime = factor(modus(data_multinom$OverTime), levels = levels(data_multinom$OverTime)),
Gender = factor(modus(data_multinom$Gender), levels = levels(data_multinom$Gender)),
MaritalStatus = factor(modus(data_multinom$MaritalStatus), levels = levels(data_multinom$MaritalStatus)),
JobSatisfaction = median(data_multinom$JobSatisfaction, na.rm = TRUE),
WorkLifeBalance = median(data_multinom$WorkLifeBalance, na.rm = TRUE),
EnvironmentSatisfaction= median(data_multinom$EnvironmentSatisfaction, na.rm = TRUE),
DistanceFromHome = median(data_multinom$DistanceFromHome, na.rm = TRUE),
PercentSalaryHike = median(data_multinom$PercentSalaryHike, na.rm = TRUE),
NumCompaniesWorked = median(data_multinom$NumCompaniesWorked, na.rm = TRUE),
TrainingTimesLastYear = median(data_multinom$TrainingTimesLastYear, na.rm = TRUE),
YearsSinceLastPromotion= median(data_multinom$YearsSinceLastPromotion, na.rm = TRUE),
YearsWithCurrManager = median(data_multinom$YearsWithCurrManager, na.rm = TRUE)
)
prob_grid_multinom <- predict(model_multinom, newdata = grid_multinom, type = "probs")
plot_prob_multinom <- grid_multinom %>%
dplyr::bind_cols(as.data.frame(prob_grid_multinom)) %>%
tidyr::pivot_longer(
cols = levels(data_multinom$JobRole),
names_to = "JobRole",
values_to = "Probabilitas"
)
ggplot(
plot_prob_multinom,
aes(x = MonthlyIncome, y = Probabilitas, color = JobRole)
) +
geom_line(linewidth = 0.9) +
scale_y_continuous(labels = scales::percent_format()) +
labs(
title = "Prediksi Probabilitas Peran Jabatan Berdasarkan Pendapatan Bulanan",
subtitle = "Variabel lain ditahan pada median atau kategori paling umum",
x = "Pendapatan bulanan",
y = "Probabilitas prediksi",
color = "Job Role"
) +
theme_minimal() +
theme(legend.position = "bottom")Grafik ini membantu membaca hasil regresi multinomial secara lebih
substantif. Jika probabilitas kategori Manager atau
Research Director meningkat ketika
MonthlyIncome naik, maka model menunjukkan bahwa karyawan
dengan pendapatan lebih tinggi cenderung menduduki jabatan manajerial
atau direktoral. Sebaliknya, jabatan yang probabilitasnya menurun
seiring kenaikan pendapatan cenderung dikaitkan dengan level pendapatan
yang lebih rendah.
Pada analisis ini, uji goodness-of-fit devians \(G^2\) tidak digunakan sebagai evaluasi utama karena data HR berbentuk data individu dengan beberapa prediktor kontinu. Uji \(G^2\) lebih cocok digunakan pada data berbentuk grouped data atau tabel frekuensi. Evaluasi model multinomial dalam laporan ini difokuskan pada koefisien model, relative risk ratio, AIC, prediksi probabilitas, dan confusion matrix.
Regresi logistik multinomial memiliki asumsi Independence of
Irrelevant Alternatives (IIA). Dalam konteks peran jabatan, asumsi ini
menyatakan bahwa perbandingan peluang antara dua kategori
JobRole tidak berubah secara drastis hanya karena kategori
lain ditambahkan atau dihapus. Asumsi IIA dicatat sebagai keterbatasan
dalam laporan ini.
sig_multi <- multinom_table %>%
dplyr::filter(term != "(Intercept)", p_value < 0.05) %>%
dplyr::arrange(Kategori, p_value)
cat("Pada taraf signifikansi 5%, model multinomial menemukan ", nrow(sig_multi),
" kombinasi prediktor-kategori yang signifikan. Semua interpretasi dibandingkan terhadap kategori acuan `Sales Executive`.\n\n", sep = "")Pada taraf signifikansi 5%, model multinomial menemukan 211 kombinasi
prediktor-kategori yang signifikan. Semua interpretasi dibandingkan
terhadap kategori acuan Sales Executive.
if (nrow(sig_multi) > 0) {
for (kelas in unique(sig_multi$Kategori)) {
cat("**JobRole ", kelas, " dibanding Sales Executive**\n\n", sep = "")
top_kelas <- sig_multi %>% dplyr::filter(Kategori == kelas) %>% head(5)
for (i in seq_len(nrow(top_kelas))) {
row <- top_kelas[i, ]
arah <- ifelse(row$Relative_Risk_Ratio > 1,
paste0("meningkatkan kecenderungan jabatan ", kelas, " dibanding Sales Executive"),
paste0("menurunkan kecenderungan jabatan ", kelas, " dibanding Sales Executive"))
cat("- **", row$Variabel_Mudah, "** memiliki RRR = ", fmt_num(row$Relative_Risk_Ratio, 3),
" dan p-value = ", fmt_p(row$p_value), ", sehingga variabel ini ", arah,
", dengan asumsi variabel lain tetap.\n", sep = "")
}
cat("\n")
}
}JobRole Healthcare Representative dibanding Sales Executive
JobRole Human Resources dibanding Sales Executive
JobRole Laboratory Technician dibanding Sales Executive
JobRole Manager dibanding Sales Executive
JobRole Manufacturing Director dibanding Sales Executive
JobRole Research Director dibanding Sales Executive
JobRole Research Scientist dibanding Sales Executive
JobRole Sales Representative dibanding Sales Executive
cat("Akurasi klasifikasi model multinomial adalah **", scales::percent(akurasi_multinom),
"**. AIC model adalah **", fmt_num(AIC_multinom, 2), "**.\n\n", sep = "")Akurasi klasifikasi model multinomial adalah 74%. AIC model adalah 2,021.73.
cat("**Kesimpulan substantif:** model multinomial menunjukkan karakteristik karyawan yang membedakan peran jabatan mereka. Faktor seperti pendapatan bulanan, departemen, total pengalaman kerja, dan bidang pendidikan merupakan pembeda utama antar kategori `JobRole`. Faktor signifikan dengan RRR lebih dari 1 meningkatkan kecenderungan suatu jabatan dibanding `Sales Executive`, sedangkan RRR kurang dari 1 menurunkannya.\n")Kesimpulan substantif: model multinomial menunjukkan
karakteristik karyawan yang membedakan peran jabatan mereka. Faktor
seperti pendapatan bulanan, departemen, total pengalaman kerja, dan
bidang pendidikan merupakan pembeda utama antar kategori
JobRole. Faktor signifikan dengan RRR lebih dari 1
meningkatkan kecenderungan suatu jabatan dibanding
Sales Executive, sedangkan RRR kurang dari 1
menurunkannya.
Tujuan analisis ini adalah mengetahui faktor-faktor yang berhubungan dengan jumlah dokter yang dikunjungi oleh pasien lansia. Data bersumber dari National Poll on Healthy Aging (NPHA). Variabel respon berupa data cacah (bilangan bulat non-negatif) sehingga Regresi Poisson merupakan pilihan model yang tepat.
| Variabel | Peran | Alasan_dipakai |
|---|---|---|
| Number of Doctors Visited | Respon | Jumlah dokter yang dikunjungi pasien lansia (data cacah) |
| Age | Prediktor | Kelompok usia pasien dapat memengaruhi frekuensi kunjungan |
| Phyiscal Health | Prediktor | Status kesehatan fisik berkaitan dengan kebutuhan kunjungan |
| Mental Health | Prediktor | Status kesehatan mental dapat mendorong atau menghambat kunjungan |
| Dental Health | Prediktor | Status kesehatan gigi sebagai dimensi kesehatan lain |
| Employment | Prediktor | Status pekerjaan mempengaruhi akses dan waktu kunjungan |
| Trouble Sleeping | Prediktor | Gangguan tidur umum mungkin berkaitan dengan jumlah kondisi medis |
| Stress Keeps Patient from Sleeping | Prediktor | Gangguan tidur akibat stres |
| Medication Keeps Patient from Sleeping | Prediktor | Gangguan tidur akibat efek obat (indikator penggunaan obat) |
| Pain Keeps Patient from Sleeping | Prediktor | Gangguan tidur akibat nyeri fisik |
| Bathroom Needs Keeps Patient from Sleeping | Prediktor | Gangguan tidur akibat kebutuhan ke toilet |
| Uknown Keeps Patient from Sleeping | Prediktor | Gangguan tidur penyebab tidak diketahui |
| Prescription Sleep Medication | Prediktor | Penggunaan obat tidur resep (indikator konsumsi obat) |
| Race | Prediktor | Kategori ras pasien sebagai faktor sosio-demografis |
| Gender | Prediktor | Jenis kelamin pasien |
# Bersihkan nama kolom agar mudah digunakan dalam R
data_poisson <- data_poisson_raw %>%
dplyr::rename(
Num_Doctors_Visited = `Number of Doctors Visited`,
Age = Age,
Physical_Health = `Phyiscal Health`,
Mental_Health = `Mental Health`,
Dental_Health = `Dental Health`,
Employment = Employment,
Stress_Sleep = `Stress Keeps Patient from Sleeping`,
Medication_Sleep = `Medication Keeps Patient from Sleeping`,
Pain_Sleep = `Pain Keeps Patient from Sleeping`,
Bathroom_Sleep = `Bathroom Needs Keeps Patient from Sleeping`,
Unknown_Sleep = `Uknown Keeps Patient from Sleeping`,
Trouble_Sleeping = `Trouble Sleeping`,
Prescription_Sleep_Med = `Prescription Sleep Medication`,
Race = Race,
Gender = Gender
) %>%
dplyr::mutate(
Num_Doctors_Visited = as.integer(Num_Doctors_Visited),
Age = as.integer(Age),
Physical_Health = as.integer(Physical_Health),
Mental_Health = as.integer(Mental_Health),
Dental_Health = as.integer(Dental_Health),
Employment = as.integer(Employment),
Stress_Sleep = as.integer(Stress_Sleep),
Medication_Sleep = as.integer(Medication_Sleep),
Pain_Sleep = as.integer(Pain_Sleep),
Bathroom_Sleep = as.integer(Bathroom_Sleep),
Unknown_Sleep = as.integer(Unknown_Sleep),
Trouble_Sleeping = as.integer(Trouble_Sleeping),
Prescription_Sleep_Med= as.integer(Prescription_Sleep_Med),
Race = as.integer(Race),
Gender = as.integer(Gender)
) %>%
tidyr::drop_na()
knitr::kable(
head(data_poisson, 8),
caption = "Cuplikan Data Poisson Setelah Cleaning"
)| Num_Doctors_Visited | Age | Physical_Health | Mental_Health | Dental_Health | Employment | Stress_Sleep | Medication_Sleep | Pain_Sleep | Bathroom_Sleep | Unknown_Sleep | Trouble_Sleeping | Prescription_Sleep_Med | Race | Gender |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 3 | 2 | 4 | 3 | 3 | 3 | 0 | 0 | 0 | 0 | 1 | 2 | 3 | 1 | 2 |
| 2 | 2 | 4 | 2 | 3 | 3 | 1 | 0 | 0 | 1 | 0 | 3 | 3 | 1 | 1 |
| 3 | 2 | 3 | 2 | 3 | 3 | 0 | 0 | 0 | 0 | 1 | 3 | 3 | 4 | 1 |
| 1 | 2 | 3 | 2 | 3 | 3 | 0 | 0 | 0 | 1 | 0 | 3 | 3 | 4 | 2 |
| 3 | 2 | 3 | 3 | 3 | 3 | 1 | 0 | 0 | 0 | 0 | 2 | 3 | 1 | 2 |
| 2 | 2 | 3 | 2 | 4 | 3 | 0 | 0 | 0 | 1 | 0 | 3 | 3 | 1 | 1 |
| 3 | 2 | 4 | 1 | 1 | 3 | 0 | 0 | 1 | 1 | 0 | 2 | 1 | 1 | 1 |
| 2 | 2 | 3 | 2 | 6 | 3 | 1 | 0 | 0 | 0 | 0 | 3 | 3 | 1 | 1 |
ringkasan_poisson <- data_poisson %>%
dplyr::summarise(
Jumlah_observasi = dplyr::n(),
Rata_rata_kunjungan = mean(Num_Doctors_Visited),
Median_kunjungan = median(Num_Doctors_Visited),
Variansi_kunjungan = var(Num_Doctors_Visited),
Min_kunjungan = min(Num_Doctors_Visited),
Max_kunjungan = max(Num_Doctors_Visited)
)
knitr::kable(ringkasan_poisson, caption = "Ringkasan Data Poisson (Jumlah Kunjungan Dokter)",
digits = 3)| Jumlah_observasi | Rata_rata_kunjungan | Median_kunjungan | Variansi_kunjungan | Min_kunjungan | Max_kunjungan |
|---|---|---|---|---|---|
| 714 | 2.112 | 2 | 0.467 | 1 | 3 |
dist_poisson <- data_poisson %>%
dplyr::count(Num_Doctors_Visited) %>%
dplyr::mutate(Persen = scales::percent(n / sum(n)))
knitr::kable(dist_poisson, caption = "Distribusi Jumlah Dokter yang Dikunjungi")| Num_Doctors_Visited | n | Persen |
|---|---|---|
| 1 | 131 | 18% |
| 2 | 372 | 52% |
| 3 | 211 | 30% |
ggplot(data_poisson, aes(x = Num_Doctors_Visited)) +
geom_histogram(binwidth = 1, fill = "#2f7f73", color = "white") +
labs(
title = "Distribusi Jumlah Dokter yang Dikunjungi Pasien Lansia",
x = "Jumlah dokter dikunjungi",
y = "Jumlah pasien"
) +
theme_minimal()Jika rata-rata dan variansi jumlah kunjungan relatif dekat, maka asumsi Poisson (mean = variansi) terpenuhi dengan baik. Jika variansi jauh lebih besar dari rata-rata, perlu dipertimbangkan quasi-Poisson atau Negative Binomial sebagai alternatif.
Model Poisson yang digunakan untuk data individual tanpa offset:
\[ \log(\mu_i) = \beta_0 + \beta_1\,\text{Age}_i + \beta_2\,\text{Physical\_Health}_i + \cdots + \beta_p X_{pi} \]
Karena data berbentuk individual (satu baris per pasien), model Poisson tanpa offset digunakan untuk mengestimasi rata-rata jumlah dokter yang dikunjungi (\(\mu_i\)) secara langsung.
model_poisson <- glm(
Num_Doctors_Visited ~ Age + Physical_Health + Mental_Health + Dental_Health +
Employment + Stress_Sleep + Medication_Sleep + Pain_Sleep +
Bathroom_Sleep + Unknown_Sleep + Trouble_Sleeping +
Prescription_Sleep_Med + Race + Gender,
data = data_poisson,
family = poisson(link = "log")
)
summary(model_poisson)##
## Call:
## glm(formula = Num_Doctors_Visited ~ Age + Physical_Health + Mental_Health +
## Dental_Health + Employment + Stress_Sleep + Medication_Sleep +
## Pain_Sleep + Bathroom_Sleep + Unknown_Sleep + Trouble_Sleeping +
## Prescription_Sleep_Med + Race + Gender, family = poisson(link = "log"),
## data = data_poisson)
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.6157139 0.2493428 2.469 0.0135 *
## Age NA NA NA NA
## Physical_Health 0.0644736 0.0348218 1.852 0.0641 .
## Mental_Health -0.0148086 0.0305598 -0.485 0.6280
## Dental_Health -0.0126340 0.0209388 -0.603 0.5463
## Employment 0.0362067 0.0465152 0.778 0.4363
## Stress_Sleep 0.0515882 0.0652498 0.791 0.4292
## Medication_Sleep 0.1333316 0.1086001 1.228 0.2195
## Pain_Sleep 0.0061158 0.0679182 0.090 0.9282
## Bathroom_Sleep 0.0416547 0.0566048 0.736 0.4618
## Unknown_Sleep 0.0290452 0.0614102 0.473 0.6362
## Trouble_Sleeping 0.0130477 0.0422874 0.309 0.7577
## Prescription_Sleep_Med -0.0511115 0.0468408 -1.091 0.2752
## Race -0.0200849 0.0264055 -0.761 0.4469
## Gender 0.0009173 0.0528779 0.017 0.9862
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 167.31 on 713 degrees of freedom
## Residual deviance: 157.05 on 700 degrees of freedom
## AIC: 2050.6
##
## Number of Fisher Scoring iterations: 4
Chunk model-poisson di atas sudah membangun
model regresi Poisson menggunakan glm() dengan
family = poisson(link = "log"). Bentuk hasil estimasinya
adalah:
\[ \log(\hat{\mu_i}) = \hat{\beta}_0 + \hat{\beta}_1 X_{1i} + \cdots + \hat{\beta}_p X_{pi} \]
Koefisien yang dieksponensialkan menjadi
Incidence Rate Ratio (IRR). IRR > 1 berarti variabel
tersebut berkaitan dengan rata-rata jumlah kunjungan dokter yang lebih
tinggi.
| Variabel_Mudah | term | estimate | IRR |
|---|---|---|---|
| (Intercept) | (Intercept) | 0.6157 | 1.8510 |
| kelompok usia pasien | Age | NA | NA |
| Physical_Health | Physical_Health | 0.0645 | 1.0666 |
| Mental_Health | Mental_Health | -0.0148 | 0.9853 |
| Dental_Health | Dental_Health | -0.0126 | 0.9874 |
| status pekerjaan pasien | Employment | 0.0362 | 1.0369 |
| Stress_Sleep | Stress_Sleep | 0.0516 | 1.0529 |
| Medication_Sleep | Medication_Sleep | 0.1333 | 1.1426 |
| Pain_Sleep | Pain_Sleep | 0.0061 | 1.0061 |
| Bathroom_Sleep | Bathroom_Sleep | 0.0417 | 1.0425 |
| Unknown_Sleep | Unknown_Sleep | 0.0290 | 1.0295 |
| Trouble_Sleeping | Trouble_Sleeping | 0.0130 | 1.0131 |
| Prescription_Sleep_Med | Prescription_Sleep_Med | -0.0511 | 0.9502 |
| kategori ras pasien | Race | -0.0201 | 0.9801 |
| jenis kelamin pasien | Gender | 0.0009 | 1.0009 |
tabel_poisson <- broom::tidy(model_poisson) %>%
dplyr::mutate(
Incidence_Rate_Ratio = exp(estimate),
Perubahan_Persen = 100 * (Incidence_Rate_Ratio - 1),
Variabel_Mudah = purrr::map_chr(term, label_poisson),
Keputusan_5persen = ifelse(p.value < 0.05, "Signifikan", "Tidak signifikan")
) %>%
dplyr::arrange(p.value)
knitr::kable(
tabel_poisson %>%
dplyr::select(Variabel_Mudah, term, estimate, std.error, statistic,
p.value, Incidence_Rate_Ratio, Perubahan_Persen, Keputusan_5persen),
digits = 4,
caption = "Koefisien dan Incidence Rate Ratio Regresi Poisson"
)| Variabel_Mudah | term | estimate | std.error | statistic | p.value | Incidence_Rate_Ratio | Perubahan_Persen | Keputusan_5persen |
|---|---|---|---|---|---|---|---|---|
| (Intercept) | (Intercept) | 0.6157 | 0.2493 | 2.4693 | 0.0135 | 1.8510 | 85.0978 | Signifikan |
| Physical_Health | Physical_Health | 0.0645 | 0.0348 | 1.8515 | 0.0641 | 1.0666 | 6.6597 | Tidak signifikan |
| Medication_Sleep | Medication_Sleep | 0.1333 | 0.1086 | 1.2277 | 0.2195 | 1.1426 | 14.2629 | Tidak signifikan |
| Prescription_Sleep_Med | Prescription_Sleep_Med | -0.0511 | 0.0468 | -1.0912 | 0.2752 | 0.9502 | -4.9827 | Tidak signifikan |
| Stress_Sleep | Stress_Sleep | 0.0516 | 0.0652 | 0.7906 | 0.4292 | 1.0529 | 5.2942 | Tidak signifikan |
| status pekerjaan pasien | Employment | 0.0362 | 0.0465 | 0.7784 | 0.4363 | 1.0369 | 3.6870 | Tidak signifikan |
| kategori ras pasien | Race | -0.0201 | 0.0264 | -0.7606 | 0.4469 | 0.9801 | -1.9885 | Tidak signifikan |
| Bathroom_Sleep | Bathroom_Sleep | 0.0417 | 0.0566 | 0.7359 | 0.4618 | 1.0425 | 4.2534 | Tidak signifikan |
| Dental_Health | Dental_Health | -0.0126 | 0.0209 | -0.6034 | 0.5463 | 0.9874 | -1.2554 | Tidak signifikan |
| Mental_Health | Mental_Health | -0.0148 | 0.0306 | -0.4846 | 0.6280 | 0.9853 | -1.4699 | Tidak signifikan |
| Unknown_Sleep | Unknown_Sleep | 0.0290 | 0.0614 | 0.4730 | 0.6362 | 1.0295 | 2.9471 | Tidak signifikan |
| Trouble_Sleeping | Trouble_Sleeping | 0.0130 | 0.0423 | 0.3085 | 0.7577 | 1.0131 | 1.3133 | Tidak signifikan |
| Pain_Sleep | Pain_Sleep | 0.0061 | 0.0679 | 0.0900 | 0.9282 | 1.0061 | 0.6135 | Tidak signifikan |
| jenis kelamin pasien | Gender | 0.0009 | 0.0529 | 0.0173 | 0.9862 | 1.0009 | 0.0918 | Tidak signifikan |
| kelompok usia pasien | Age | NA | NA | NA | NA | NA | NA | NA |
## [1] 2050.627
## [1] 0.2243558
if (dispersion_poisson > 1.5) {
cat("Terdapat indikasi overdispersion (rasio deviance/df =",
round(dispersion_poisson, 3),
"). Inferensi p-value sebaiknya dilihat juga menggunakan quasi-Poisson.\n")
} else {
cat("Tidak ada indikasi overdispersion yang kuat berdasarkan rasio deviance/df residual =",
round(dispersion_poisson, 3), ".\n")
}## Tidak ada indikasi overdispersion yang kuat berdasarkan rasio deviance/df residual = 0.224 .
Selain rasio deviance terhadap derajat bebas residual, overdispersion juga dapat diperiksa menggunakan Pearson dispersion.
\[ \hat{\phi} = \frac{\sum r_{P,i}^2}{df} \]
Jika nilai dispersion jauh lebih besar dari 1, maka terdapat indikasi overdispersion.
pearson_dispersion <- sum(residuals(model_poisson, type = "pearson")^2) /
df.residual(model_poisson)
tibble::tibble(
Ukuran = "Pearson dispersion",
Nilai = pearson_dispersion,
Interpretasi = dplyr::case_when(
pearson_dispersion < 1.5 ~ "Tidak ada indikasi overdispersion berat",
pearson_dispersion < 2.5 ~ "Ada indikasi overdispersion sedang",
TRUE ~ "Ada indikasi overdispersion kuat"
)
) %>%
dplyr::mutate(Nilai = round(Nilai, 3)) %>%
knitr::kable(caption = "Pemeriksaan Overdispersion dengan Pearson Dispersion")| Ukuran | Nilai | Interpretasi |
|---|---|---|
| Pearson dispersion | 0.212 | Tidak ada indikasi overdispersion berat |
Jika nilai Pearson dispersion jauh lebih besar dari 1, maka variasi jumlah kunjungan dokter lebih besar daripada yang diasumsikan oleh model Poisson. Dalam kondisi tersebut, hasil signifikansi dari model Poisson standar perlu dibaca hati-hati dan model quasi-Poisson lebih aman digunakan sebagai pembanding.
model_quasipoisson <- glm(
Num_Doctors_Visited ~ Age + Physical_Health + Mental_Health + Dental_Health +
Employment + Stress_Sleep + Medication_Sleep + Pain_Sleep +
Bathroom_Sleep + Unknown_Sleep + Trouble_Sleeping +
Prescription_Sleep_Med + Race + Gender,
data = data_poisson,
family = quasipoisson(link = "log")
)
tabel_quasipoisson <- broom::tidy(model_quasipoisson) %>%
dplyr::mutate(
Incidence_Rate_Ratio = exp(estimate),
Perubahan_Persen = 100 * (Incidence_Rate_Ratio - 1),
Variabel_Mudah = purrr::map_chr(term, label_poisson),
Keputusan_5persen = ifelse(p.value < 0.05, "Signifikan", "Tidak signifikan")
) %>%
dplyr::arrange(p.value)
knitr::kable(
tabel_quasipoisson %>%
dplyr::select(Variabel_Mudah, term, estimate, std.error, statistic,
p.value, Incidence_Rate_Ratio, Perubahan_Persen, Keputusan_5persen),
digits = 4,
caption = "Koefisien dan IRR Quasi-Poisson"
)| Variabel_Mudah | term | estimate | std.error | statistic | p.value | Incidence_Rate_Ratio | Perubahan_Persen | Keputusan_5persen |
|---|---|---|---|---|---|---|---|---|
| (Intercept) | (Intercept) | 0.6157 | 0.1149 | 5.3599 | 0.0000 | 1.8510 | 85.0978 | Signifikan |
| Physical_Health | Physical_Health | 0.0645 | 0.0160 | 4.0189 | 0.0001 | 1.0666 | 6.6597 | Signifikan |
| Medication_Sleep | Medication_Sleep | 0.1333 | 0.0500 | 2.6649 | 0.0079 | 1.1426 | 14.2629 | Signifikan |
| Prescription_Sleep_Med | Prescription_Sleep_Med | -0.0511 | 0.0216 | -2.3685 | 0.0181 | 0.9502 | -4.9827 | Signifikan |
| Stress_Sleep | Stress_Sleep | 0.0516 | 0.0301 | 1.7161 | 0.0866 | 1.0529 | 5.2942 | Tidak signifikan |
| status pekerjaan pasien | Employment | 0.0362 | 0.0214 | 1.6895 | 0.0916 | 1.0369 | 3.6870 | Tidak signifikan |
| kategori ras pasien | Race | -0.0201 | 0.0122 | -1.6510 | 0.0992 | 0.9801 | -1.9885 | Tidak signifikan |
| Bathroom_Sleep | Bathroom_Sleep | 0.0417 | 0.0261 | 1.5973 | 0.1106 | 1.0425 | 4.2534 | Tidak signifikan |
| Dental_Health | Dental_Health | -0.0126 | 0.0096 | -1.3097 | 0.1907 | 0.9874 | -1.2554 | Tidak signifikan |
| Mental_Health | Mental_Health | -0.0148 | 0.0141 | -1.0518 | 0.2932 | 0.9853 | -1.4699 | Tidak signifikan |
| Unknown_Sleep | Unknown_Sleep | 0.0290 | 0.0283 | 1.0266 | 0.3050 | 1.0295 | 2.9471 | Tidak signifikan |
| Trouble_Sleeping | Trouble_Sleeping | 0.0130 | 0.0195 | 0.6697 | 0.5033 | 1.0131 | 1.3133 | Tidak signifikan |
| Pain_Sleep | Pain_Sleep | 0.0061 | 0.0313 | 0.1955 | 0.8451 | 1.0061 | 0.6135 | Tidak signifikan |
| jenis kelamin pasien | Gender | 0.0009 | 0.0244 | 0.0377 | 0.9700 | 1.0009 | 0.0918 | Tidak signifikan |
| kelompok usia pasien | Age | NA | NA | NA | NA | NA | NA | NA |
data_poisson$fitted_poisson <- fitted(model_poisson)
ggplot(data_poisson, aes(x = fitted_poisson, y = Num_Doctors_Visited)) +
geom_point(alpha = 0.35, color = "#2f7f73") +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "#b84f5a") +
labs(
title = "Observed vs Fitted Jumlah Dokter yang Dikunjungi",
x = "Nilai fitted model",
y = "Jumlah kunjungan dokter aktual"
) +
theme_minimal()Visualisasi ini digunakan untuk melihat prediksi rata-rata kunjungan dokter berdasarkan status kesehatan fisik dan mental pasien.
grid_poisson <- expand.grid(
Physical_Health = 1:5,
Mental_Health = c(1, 3, 5),
Age = median(data_poisson$Age, na.rm = TRUE),
Dental_Health = median(data_poisson$Dental_Health, na.rm = TRUE),
Employment = modus(data_poisson$Employment),
Stress_Sleep = modus(data_poisson$Stress_Sleep),
Medication_Sleep = modus(data_poisson$Medication_Sleep),
Pain_Sleep = modus(data_poisson$Pain_Sleep),
Bathroom_Sleep = modus(data_poisson$Bathroom_Sleep),
Unknown_Sleep = modus(data_poisson$Unknown_Sleep),
Trouble_Sleeping = modus(data_poisson$Trouble_Sleeping),
Prescription_Sleep_Med= modus(data_poisson$Prescription_Sleep_Med),
Race = modus(data_poisson$Race),
Gender = modus(data_poisson$Gender)
)
pred_rate_poisson <- predict(
model_poisson,
newdata = grid_poisson,
type = "link",
se.fit = TRUE
)
grid_poisson_plot <- grid_poisson %>%
dplyr::mutate(
fit_link = pred_rate_poisson$fit,
se_link = pred_rate_poisson$se.fit,
Predicted_Count = exp(fit_link),
Count_Lower = exp(fit_link - 1.96 * se_link),
Count_Upper = exp(fit_link + 1.96 * se_link),
Mental_Health = factor(Mental_Health,
labels = c("Mental baik (1)", "Mental sedang (3)", "Mental buruk (5)"))
)
ggplot(
grid_poisson_plot,
aes(x = Physical_Health, y = Predicted_Count, color = Mental_Health, group = Mental_Health)
) +
geom_line(linewidth = 1.1) +
geom_ribbon(aes(ymin = Count_Lower, ymax = Count_Upper, fill = Mental_Health),
alpha = 0.15, color = NA) +
labs(
title = "Prediksi Rata-rata Jumlah Kunjungan Dokter",
subtitle = "Berdasarkan status kesehatan fisik dan mental pasien",
x = "Status kesehatan fisik (1 = sangat baik, 5 = sangat buruk)",
y = "Prediksi rata-rata jumlah kunjungan dokter",
color = "Status mental",
fill = "Status mental"
) +
theme_minimal()Grafik ini menunjukkan prediksi rata-rata jumlah kunjungan dokter menurut status kesehatan fisik dan mental pasien, dengan variabel lain ditahan pada nilai paling umum. Jika garis cenderung naik ketika status kesehatan fisik memburuk (nilai lebih tinggi), maka model menunjukkan bahwa kondisi fisik yang lebih buruk berkaitan dengan lebih banyak kunjungan dokter.
if (dispersion_poisson > 1.5) {
tabel_pakai_poisson <- tabel_quasipoisson
metode_pakai <- "quasi-Poisson"
} else {
tabel_pakai_poisson <- tabel_poisson
metode_pakai <- "Poisson"
}
sig_poisson <- tabel_pakai_poisson %>%
dplyr::filter(term != "(Intercept)", p.value < 0.05) %>%
dplyr::arrange(p.value)
cat("Rasio deviance/df residual adalah **", fmt_num(dispersion_poisson, 3), "**. ", sep = "")Rasio deviance/df residual adalah 0.224.
if (dispersion_poisson > 1.5) {
cat("Nilai ini menunjukkan overdispersion, sehingga interpretasi signifikansi lebih aman mengacu pada model quasi-Poisson.\n\n")
} else {
cat("Nilai ini tidak menunjukkan overdispersion yang kuat, sehingga model Poisson standar masih memadai.\n\n")
}Nilai ini tidak menunjukkan overdispersion yang kuat, sehingga model Poisson standar masih memadai.
cat("Dengan acuan inferensi ", metode_pakai, ", terdapat ", nrow(sig_poisson),
" prediktor signifikan pada taraf 5%.\n\n", sep = "")Dengan acuan inferensi Poisson, terdapat 0 prediktor signifikan pada taraf 5%.
if (nrow(sig_poisson) > 0) {
cat("Prediktor paling penting berdasarkan p-value terkecil adalah:\n\n")
top_poisson <- head(sig_poisson, 10)
for (i in seq_len(nrow(top_poisson))) {
row <- top_poisson[i, ]
arah <- ifelse(row$Incidence_Rate_Ratio > 1,
"berkaitan dengan rata-rata kunjungan dokter yang lebih tinggi",
"berkaitan dengan rata-rata kunjungan dokter yang lebih rendah")
cat("- **", row$Variabel_Mudah, "** memiliki IRR = ", fmt_num(row$Incidence_Rate_Ratio, 3),
" dengan perubahan sekitar ", fmt_num(row$Perubahan_Persen, 2), "%",
" dan p-value = ", fmt_p(row$p.value), ", sehingga variabel ini ", arah,
", dengan asumsi variabel lain tetap.\n", sep = "")
}
cat("\n")
}
cat("AIC model Poisson adalah **", fmt_num(AIC_poisson, 2), "**.\n\n", sep = "")AIC model Poisson adalah 2,050.63.
cat("**Kesimpulan substantif:** jumlah dokter yang dikunjungi pasien lansia tidak hanya dipengaruhi oleh status kesehatan fisik dan mental, tetapi juga oleh gangguan tidur, status pekerjaan, dan karakteristik demografi. Karena data berbentuk individual tanpa offset, model Poisson mengestimasi rata-rata jumlah kunjungan dokter secara langsung berdasarkan kondisi tiap pasien.\n")Kesimpulan substantif: jumlah dokter yang dikunjungi pasien lansia tidak hanya dipengaruhi oleh status kesehatan fisik dan mental, tetapi juga oleh gangguan tidur, status pekerjaan, dan karakteristik demografi. Karena data berbentuk individual tanpa offset, model Poisson mengestimasi rata-rata jumlah kunjungan dokter secara langsung berdasarkan kondisi tiap pasien.
| Model | Dataset | Respon | Fungsi link | Ukuran interpretasi | Evaluasi utama |
|---|---|---|---|---|---|
| Logistik biner | Telco Customer Churn | Churn 0/1 | Logit | Odds ratio | Confusion matrix, AIC, ROC/AUC |
| Ordinal | Student Performance (Portugis) | Performa G3 (Rendah/Sedang/Tinggi) | Cumulative logit MASS::polr:
logit(P(Y <= j)) = alpha_j - x beta |
exp(beta) sebagai OR menuju performa lebih tinggi |
Confusion matrix, AIC, proportional odds check |
| Multinomial | HR Employee Attrition | JobRole (9 kategori) | Baseline-category logit vs Sales Executive | Relative risk ratio | Confusion matrix, AIC, prediksi probabilitas, catatan IIA |
| Poisson | NPHA Doctor Visits | Jumlah kunjungan dokter | Log (tanpa offset, data individu) | Incidence rate ratio dan perubahan persen | AIC, deviance dispersion, Pearson dispersion, overdispersion |
cat("Berdasarkan empat analisis yang dilakukan, setiap model menjawab pertanyaan substantif yang berbeda dari empat konteks data yang beragam.\n\n")Berdasarkan empat analisis yang dilakukan, setiap model menjawab pertanyaan substantif yang berbeda dari empat konteks data yang beragam.
cat("- **Model biner** menjelaskan peluang pelanggan telekomunikasi melakukan churn (berhenti berlangganan). Model ini cocok karena respon hanya terdiri dari dua kategori: churn dan tidak churn. Faktor seperti jenis kontrak, lama berlangganan, dan layanan internet menjadi penentu utama.\n")cat("- **Model ordinal** menjelaskan tingkat performa akademik siswa berdasarkan nilai akhir G3 yang dikelompokkan menjadi Rendah, Sedang, dan Tinggi. Model ini cocok karena performa memiliki urutan yang bermakna. Faktor seperti kegagalan kelas sebelumnya, waktu belajar, dan konsumsi alkohol menjadi prediktor penting.\n")cat("- **Model multinomial** menjelaskan peran jabatan karyawan (JobRole) di perusahaan. Model ini cocok karena jabatan memiliki beberapa kategori yang tidak diperlakukan sebagai urutan. Pendapatan bulanan, departemen, dan total pengalaman kerja membedakan antar peran jabatan.\n")cat("- **Model Poisson** menjelaskan jumlah dokter yang dikunjungi oleh pasien lansia. Model ini cocok karena respon berupa hitungan kejadian yang bersifat diskret dan non-negatif. Status kesehatan fisik, mental, dan gangguan tidur menjadi faktor utama penjelas.\n\n")cat("Secara keseluruhan, analisis ini menunjukkan bahwa pemilihan model harus disesuaikan dengan skala pengukuran variabel respon. Model biner untuk respon dikotomi, model ordinal untuk respon bertingkat, model multinomial untuk respon nominal banyak kategori, dan model Poisson untuk data cacah. Keempat model ini merupakan bagian penting dari keluarga generalized linear model (GLM).\n")Secara keseluruhan, analisis ini menunjukkan bahwa pemilihan model harus disesuaikan dengan skala pengukuran variabel respon. Model biner untuk respon dikotomi, model ordinal untuk respon bertingkat, model multinomial untuk respon nominal banyak kategori, dan model Poisson untuk data cacah. Keempat model ini merupakan bagian penting dari keluarga generalized linear model (GLM).