Analisis Regresi Logistik: Biner, Multinomial, Ordinal, dan Poisson
Studi Kasus dari Empat Dataset Mendeley
Salsabila Najwa Effianti
2026-06-01
suppressPackageStartupMessages({
library(readxl)
library(MASS) # polr — load sebelum tidyverse
library(ordinal) # clm, nominal_test
library(nnet) # multinom
library(broom)
library(knitr)
library(kableExtra)
library(scales)
library(tibble)
library(purrr)
library(car) # vif
library(AER) # dispersiontest
library(lmtest) # lrtest
library(ResourceSelection) # Hosmer-Lemeshow
library(tidyverse) # load paling akhir
})theme_analisis <- function(base_size = 12) {
theme_minimal(base_size = base_size) +
theme(
plot.title = element_text(face = "bold", size = base_size + 4, color = "#243142"),
plot.subtitle = element_text(size = base_size - 1, color = "#64748b"),
axis.title = element_text(face = "bold", color = "#243142"),
axis.text = element_text(color = "#243142"),
panel.grid.minor = element_blank(),
panel.grid.major = element_line(color = "#d9e2ec", linewidth = 0.35),
legend.position = "bottom",
legend.title = element_text(face = "bold"),
strip.text = element_text(face = "bold", color = "#243142"),
plot.background = element_rect(fill = "#ffffff", color = NA),
panel.background = element_rect(fill = "#ffffff", color = NA)
)
}
pal <- c("#2f7f73", "#5568b8", "#d18b2f", "#b84f5a", "#789f90", "#243142", "#8b6fc9")1 Regresi Logistik Biner
1.1 Deskripsi Data
Dataset ini diperoleh dari Mendeley Data (DOI: 10.17632/mz3c929xwr.1) yang dikumpulkan oleh Boateng & Abaye (2019). Data mencakup 395 ibu hamil dari sepuluh wilayah Ghana yang mengunjungi klinik antenatal dan postnatal selama periode 2012–2016 (data dari Ghana Health Service). Variabel respons adalah jenis persalinan (cesar atau bukan), yang merupakan variabel biner, sehingga regresi logistik biner merupakan metode yang tepat. Prediktor meliputi faktor non-medis seperti tingkat pendidikan, paritas, berat lahir bayi, riwayat persalinan cesar sebelumnya, lokasi, usia ibu, dan periode dalam tahun kelahiran. Analisis ini bertujuan memodelkan probabilitas riwayat persalinan caesar pada ibu bersalin menggunakan regresi logistik biner. Dataset yang digunakan adalah Boateng EY et al. Caesarean Data 2019 yang bersumber dari Mendeley Data.
Variabel dependen:
PreviousCaesarean Delivery
- 1 = Tidak Caesar (dikode ulang menjadi 0)
- 2 = Ya Caesar (dikode ulang menjadi 1)
Variabel prediktor (semuanya kategorik):
| Variabel | Keterangan |
|---|---|
| Religion | Agama (Catholic, Protestant, Islam, Traditional/Other) |
| Age | Kelompok usia ibu |
| Location | Lokasi tempat tinggal |
| PeriodofChildBirth | Periode kelahiran anak |
| WeightofBirth | Berat bayi lahir |
| FacilityType | Tipe fasilitas kesehatan |
| MaritalStatus | Status pernikahan |
| EthnicGroup | Kelompok etnis |
| Mother’s Occupation | Pekerjaan ibu |
| Father’s Occupation | Pekerjaan ayah |
1.2 Load Library dan Data
raw_caesar <- read_excel(
"C:/Users/Salsabila Najwa/Downloads/Boateng EY et al Caesarean Data 2019.xlsx"
)
glimpse(raw_caesar)## Rows: 395
## Columns: 11
## $ Religion <dbl> 1, 2, 3, 4, 3, 2, 3, 2, 3, 4, 2, 3, 3, 4,…
## $ `PreviousCaesarean Delivery` <dbl> 1, 2, 1, 1, 1, 1, 2, 1, 1, 2, 1, 1, 1, 2,…
## $ Age <dbl> 3, 1, 4, 3, 4, 3, 1, 4, 3, 1, 4, 4, 3, 2,…
## $ Location <dbl> 1, 2, 1, 3, 3, 1, 1, 2, 3, 2, 3, 2, 3, 1,…
## $ PeriodofChildBirth <dbl> 2, 4, 3, 2, 2, 3, 1, 2, 3, 4, 2, 1, 2, 4,…
## $ WeightofBirth <dbl> 2, 3, 2, 3, 2, 1, 3, 1, 2, 3, 2, 3, 1, 2,…
## $ FacilityType <dbl> 1, 2, 1, 3, 1, 4, 2, 3, 4, 1, 2, 3, 4, 5,…
## $ MaritalStatus <dbl> 1, 2, 3, 4, 5, 3, 4, 5, 2, 1, 2, 3, 2, 1,…
## $ EthnicGroup <dbl> 5, 2, 3, 4, 5, 4, 2, 5, 1, 3, 5, 2, 3, 2,…
## $ `Mother'sOcuupation` <dbl> 3, 2, 4, 5, 1, 2, 3, 4, 3, 4, 2, 3, 2, 1,…
## $ `Father'sOccupation` <dbl> 4, 2, 1, 5, 2, 3, 4, 1, 4, 5, 2, 3, 1, 4,…📊 Interpretasi output: Dataset terdiri dari 395 observasi dan 11 variabel. Seluruh variabel tersimpan dalam tipe numerik (
<dbl>) karena menggunakan kode angka, sehingga perlu diubah menjadi faktor berlabel agar analisis dan interpretasi lebih mudah dipahami.
1.3 Persiapan Data
1.3.1 Preprocessing dan Pemberian Label
caesar <- raw_caesar %>%
rename(
religion = Religion,
prev_caesar = `PreviousCaesarean Delivery`,
age = Age,
location = Location,
period_birth = PeriodofChildBirth,
weight_birth = WeightofBirth,
facility = FacilityType,
marital = MaritalStatus,
ethnic = EthnicGroup,
occ_mother = `Mother'sOcuupation`,
occ_father = `Father'sOccupation`
) %>%
mutate(
caesar_bin = ifelse(prev_caesar == 2, 1, 0),
religion = factor(religion,
levels = 1:4,
labels = c("Catholic", "Protestant",
"Islam", "Traditional/Other")),
age = factor(age,
levels = 1:4,
labels = c("<20 tahun", "20-29 tahun",
"30-39 tahun", ">=40 tahun")),
location = factor(location,
levels = 1:3,
labels = c("Urban", "Rural", "Periurban")),
period_birth = factor(period_birth,
levels = 1:4,
labels = c("<1 tahun", "1-3 tahun",
"3-5 tahun", ">5 tahun")),
weight_birth = factor(weight_birth,
levels = 1:3,
labels = c("Rendah (<2.5kg)",
"Normal (2.5-3.9kg)",
"Tinggi (>=4kg)")),
facility = factor(facility,
levels = 1:6,
labels = c("Teaching", "Regional", "Distrik",
"Health Centre", "Swasta", "Lainnya")),
marital = factor(marital,
levels = 1:5,
labels = c("Lajang", "Menikah", "Kohabitasi",
"Cerai", "Janda/Duda")),
ethnic = factor(ethnic,
levels = 1:6,
labels = c("Akan", "Mole-Dagbani", "Ewe",
"Ga-Dangme", "Guan", "Lainnya")),
occ_mother = factor(occ_mother,
levels = 1:5,
labels = c("Petani", "Pedagang", "PNS",
"Pelajar", "Lainnya")),
occ_father = factor(occ_father,
levels = 1:5,
labels = c("Petani", "Pedagang", "PNS",
"Pelajar", "Lainnya"))
) %>%
# Hapus kategori langka (n=2) untuk menghindari quasi-separation
filter(!(facility == "Lainnya"), !(ethnic == "Lainnya")) %>%
droplevels()1.3.2 Pengecekan Data
## Missing values per kolom:## religion prev_caesar age location period_birth weight_birth
## 0 0 0 0 0 0
## facility marital ethnic occ_mother occ_father caesar_bin
## 0 0 0 0 0 0📊 Interpretasi output: Tidak terdapat missing value pada seluruh variabel. Dataset sudah bersih dan siap digunakan untuk analisis tanpa memerlukan imputasi data.
dist_dep <- data.frame(
Status = c("Tidak Caesar (0)", "Caesar (1)"),
Frekuensi = as.integer(table(caesar$caesar_bin)),
Proporsi = round(prop.table(table(caesar$caesar_bin)), 4)
)
knitr::kable(dist_dep, caption = "Distribusi variabel dependen") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Status | Frekuensi | Proporsi.Var1 | Proporsi.Freq |
|---|---|---|---|
| Tidak Caesar (0) | 244 | 0 | 0.624 |
| Caesar (1) | 147 | 1 | 0.376 |
📊 Interpretasi output: Dari 391 ibu yang dianalisis (setelah membuang 4 observasi langka), sebanyak 244 ibu (62,4%) tidak memiliki riwayat persalinan caesar, dan 147 ibu (37,6%) memiliki riwayat persalinan caesar. Proporsi kelas tidak seimbang sempurna, namun proporsi caesar masih cukup besar untuk dipelajari oleh model.
1.4 Statistik Deskriptif
1.4.1 Distribusi Variabel Prediktor
caesar %>%
dplyr::select(religion, age, location, period_birth, weight_birth,
facility, marital, ethnic, occ_mother, occ_father) %>%
summary()## religion age location period_birth
## Catholic : 38 <20 tahun : 75 Urban :191 <1 tahun : 80
## Protestant :109 20-29 tahun:117 Rural :140 1-3 tahun:105
## Islam :126 30-39 tahun:107 Periurban: 60 3-5 tahun:120
## Traditional/Other:118 >=40 tahun : 92 >5 tahun : 86
##
## weight_birth facility marital ethnic
## Rendah (<2.5kg) : 94 Teaching :83 Lajang :63 Akan :77
## Normal (2.5-3.9kg):180 Regional :84 Menikah :89 Mole-Dagbani:97
## Tinggi (>=4kg) :117 Distrik :77 Kohabitasi:98 Ewe :99
## Health Centre:72 Cerai :92 Ga-Dangme :65
## Swasta :75 Janda/Duda:49 Guan :53
## occ_mother occ_father
## Petani : 68 Petani : 65
## Pedagang: 85 Pedagang: 89
## PNS :100 PNS :104
## Pelajar : 89 Pelajar : 92
## Lainnya : 49 Lainnya : 411.4.2 Tabulasi Silang: Proporsi Caesar per Kategori
tab_bivariat <- function(var, label) {
caesar %>%
count(!!sym(var), caesar_bin) %>%
group_by(!!sym(var)) %>%
mutate(Proporsi_Caesar = round(n / sum(n), 3)) %>%
filter(caesar_bin == 1) %>%
dplyr::select(Kategori = !!sym(var), n_Caesar = n, Proporsi_Caesar) %>%
knitr::kable(caption = paste("Proporsi Caesar menurut", label)) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)
}| Kategori | n_Caesar | Proporsi_Caesar |
|---|---|---|
| Catholic | 9 | 0.237 |
| Protestant | 37 | 0.339 |
| Islam | 56 | 0.444 |
| Traditional/Other | 45 | 0.381 |
💡 Interpretasi: Proporsi caesar tertinggi pada kelompok Islam (44,4%) dan terendah pada kelompok Catholic (23,7%). Perbedaan ini mengindikasikan adanya kemungkinan pengaruh kelompok agama terhadap riwayat persalinan caesar.
| Kategori | n_Caesar | Proporsi_Caesar |
|---|---|---|
| <20 tahun | 43 | 0.573 |
| 20-29 tahun | 57 | 0.487 |
| 30-39 tahun | 35 | 0.327 |
| >=40 tahun | 12 | 0.130 |
💡 Interpretasi: Proporsi caesar cenderung menurun seiring bertambahnya usia. Kelompok usia termuda (<20 tahun) memiliki proporsi caesar tertinggi (57,3%), sementara kelompok ≥40 tahun memiliki proporsi terendah (13,0%). Pola ini menunjukkan indikasi awal bahwa variabel usia berpotensi menjadi prediktor penting.
| Kategori | n_Caesar | Proporsi_Caesar |
|---|---|---|
| Urban | 40 | 0.209 |
| Rural | 98 | 0.700 |
| Periurban | 9 | 0.150 |
💡 Interpretasi: Ibu yang tinggal di daerah Rural memiliki proporsi caesar sangat tinggi (70,0%), jauh di atas Urban (20,9%) dan Periurban (15,0%). Kesenjangan ini sangat mencolok dan mengisyaratkan bahwa lokasi merupakan prediktor yang kuat.
| Kategori | n_Caesar | Proporsi_Caesar |
|---|---|---|
| <1 tahun | 19 | 0.238 |
| 1-3 tahun | 35 | 0.333 |
| 3-5 tahun | 48 | 0.400 |
| >5 tahun | 45 | 0.523 |
💡 Interpretasi: Semakin lama jarak sejak kelahiran anak terakhir, semakin tinggi proporsi caesar. Periode >5 tahun menunjukkan proporsi tertinggi (52,3%), dibandingkan <1 tahun yang hanya 23,8%.
| Kategori | n_Caesar | Proporsi_Caesar |
|---|---|---|
| Rendah (<2.5kg) | 41 | 0.436 |
| Normal (2.5-3.9kg) | 67 | 0.372 |
| Tinggi (>=4kg) | 39 | 0.333 |
| Kategori | n_Caesar | Proporsi_Caesar |
|---|---|---|
| Teaching | 35 | 0.422 |
| Regional | 35 | 0.417 |
| Distrik | 27 | 0.351 |
| Health Centre | 28 | 0.389 |
| Swasta | 22 | 0.293 |
| Kategori | n_Caesar | Proporsi_Caesar |
|---|---|---|
| Lajang | 23 | 0.365 |
| Menikah | 29 | 0.326 |
| Kohabitasi | 37 | 0.378 |
| Cerai | 31 | 0.337 |
| Janda/Duda | 27 | 0.551 |
| Kategori | n_Caesar | Proporsi_Caesar |
|---|---|---|
| Akan | 21 | 0.273 |
| Mole-Dagbani | 36 | 0.371 |
| Ewe | 39 | 0.394 |
| Ga-Dangme | 33 | 0.508 |
| Guan | 18 | 0.340 |
| Kategori | n_Caesar | Proporsi_Caesar |
|---|---|---|
| Petani | 27 | 0.397 |
| Pedagang | 30 | 0.353 |
| PNS | 41 | 0.410 |
| Pelajar | 32 | 0.360 |
| Lainnya | 17 | 0.347 |
| Kategori | n_Caesar | Proporsi_Caesar |
|---|---|---|
| Petani | 24 | 0.369 |
| Pedagang | 40 | 0.449 |
| PNS | 41 | 0.394 |
| Pelajar | 32 | 0.348 |
| Lainnya | 10 | 0.244 |
1.5 Eksplorasi Visual
ggplot(caesar, aes(x = age, fill = factor(caesar_bin))) +
geom_bar(position = "fill", alpha = 0.85) +
scale_fill_manual(values = c("0" = "#2196F3", "1" = "#F44336"),
labels = c("Tidak Caesar", "Caesar")) +
scale_y_continuous(labels = percent_format()) +
labs(title = "Proporsi Caesar berdasarkan Kelompok Usia",
x = "Kelompok Usia", y = "Proporsi", fill = "Status") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 15, hjust = 1))
Proporsi Caesar berdasarkan Kelompok Usia
💡 Interpretasi: Grafik memperlihatkan penurunan proporsi caesar yang konsisten seiring bertambahnya usia. Kelompok <20 tahun memiliki lebih dari separuh kasus caesar, sedangkan kelompok ≥40 tahun justru memiliki proporsi sangat rendah.
ggplot(caesar, aes(x = location, fill = factor(caesar_bin))) +
geom_bar(position = "fill", alpha = 0.85) +
scale_fill_manual(values = c("0" = "#2196F3", "1" = "#F44336"),
labels = c("Tidak Caesar", "Caesar")) +
scale_y_continuous(labels = percent_format()) +
labs(title = "Proporsi Caesar berdasarkan Lokasi",
x = "Lokasi", y = "Proporsi", fill = "Status") +
theme_minimal()
Proporsi Caesar berdasarkan Lokasi
💡 Interpretasi: Daerah Rural mendominasi kasus caesar dengan proporsi mencapai 70%, sangat kontras dengan Urban dan Periurban. Kondisi ini kemungkinan terkait perbedaan akses dan jenis fasilitas kesehatan yang tersedia.
ggplot(caesar, aes(x = facility, fill = factor(caesar_bin))) +
geom_bar(position = "fill", alpha = 0.85) +
scale_fill_manual(values = c("0" = "#2196F3", "1" = "#F44336"),
labels = c("Tidak Caesar", "Caesar")) +
scale_y_continuous(labels = percent_format()) +
labs(title = "Proporsi Caesar berdasarkan Tipe Fasilitas",
x = "Tipe Fasilitas", y = "Proporsi", fill = "Status") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 15, hjust = 1))
Proporsi Caesar berdasarkan Tipe Fasilitas
1.6 Pembagian Training dan Testing
Data dibagi menjadi 80% training dan 20% testing secara stratified agar proporsi caesar dan tidak caesar tetap seimbang pada kedua subset.
stratified_split <- function(y, prop = 0.8) {
idx_by_class <- split(seq_along(y), y)
train_idx <- lapply(
idx_by_class,
function(idx) sample(idx, size = floor(length(idx) * prop))
)
unlist(train_idx, use.names = FALSE)
}
set.seed(42)
train_id <- stratified_split(caesar$caesar_bin, prop = 0.8)
train_data <- caesar[train_id, ]
test_data <- caesar[-train_id, ]
split_summary <- bind_rows(
train_data %>% count(caesar_bin) %>% mutate(data = "Training"),
test_data %>% count(caesar_bin) %>% mutate(data = "Testing")
) %>%
group_by(data) %>%
mutate(Proporsi = scales::percent(n / sum(n), accuracy = 0.1)) %>%
ungroup() %>%
mutate(Status = ifelse(caesar_bin == 1, "Caesar", "Tidak Caesar")) %>%
dplyr::select(Data = data, Status, Jumlah = n, Proporsi)
knitr::kable(split_summary,
caption = "Distribusi kelas pada training dan testing") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Data | Status | Jumlah | Proporsi |
|---|---|---|---|
| Training | Tidak Caesar | 195 | 62.5% |
| Training | Caesar | 117 | 37.5% |
| Testing | Tidak Caesar | 49 | 62.0% |
| Testing | Caesar | 30 | 38.0% |
📊 Interpretasi output: Data training berisi 312 observasi dan data testing berisi 79 observasi. Karena split dilakukan secara stratified, proporsi caesar pada training (37,5%) dan testing (38,0%) mendekati proporsi data keseluruhan. Ini penting agar evaluasi testing tidak bias akibat perbedaan komposisi kelas.
1.7 Model Regresi Logistik Biner
1.7.1 Rumus Model Logistik
Misalkan \(Y_i\) adalah status persalinan caesar ibu ke-\(i\), dengan:
\[ Y_i = \begin{cases} 1, & \text{jika memiliki riwayat persalinan caesar,}\\ 0, & \text{jika tidak.} \end{cases} \]
Peluang persalinan caesar dinotasikan sebagai:
\[ p_i = P(Y_i = 1 \mid X_i). \]
Pada regresi logistik biner, peluang tersebut dimodelkan melalui fungsi logit:
\[ \text{logit}(p_i) = \ln\left(\frac{p_i}{1-p_i}\right) = \eta_i = \beta_0 + \beta_1X_{1i} + \beta_2X_{2i} + \cdots + \beta_kX_{ki}. \]
Bentuk peluang prediksi diperoleh dengan transformasi balik:
\[ \hat{p}_i = \frac{\exp(\hat{\eta}_i)}{1+\exp(\hat{\eta}_i)} = \frac{1}{1+\exp(-\hat{\eta}_i)}. \]
1.7.1.1 Interpretasi Odds Ratio
Untuk prediktor kategorik, odds ratio (OR) dihitung sebagai:
\[ OR_j = \exp(\hat{\beta}_j). \]
- Jika \(OR_j > 1\): kategori
tersebut berkaitan dengan peningkatan odds persalinan caesar
dibandingkan kategori referensi.
- Jika \(OR_j < 1\): kategori tersebut berkaitan dengan penurunan odds persalinan caesar dibandingkan kategori referensi.
1.7.2 Model Penuh
heart_fit <- glm(
caesar_bin ~ religion + age + location + period_birth +
weight_birth + facility + marital + ethnic +
occ_mother + occ_father,
data = train_data,
family = binomial(link = "logit")
)
ringkasan_model <- data.frame(
Keterangan = c("Jumlah observasi training", "Null deviance",
"Residual deviance", "Derajat bebas residual", "AIC"),
Nilai = c(
nobs(heart_fit),
round(heart_fit$null.deviance, 3),
round(heart_fit$deviance, 3),
heart_fit$df.residual,
round(AIC(heart_fit), 3)
)
)
knitr::kable(ringkasan_model,
caption = "Ringkasan kecocokan model penuh") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Keterangan | Nilai |
|---|---|
| Jumlah observasi training | 312.000 |
| Null deviance | 412.815 |
| Residual deviance | 259.290 |
| Derajat bebas residual | 278.000 |
| AIC | 327.290 |
📊 Interpretasi output: Model diestimasi menggunakan 312 observasi training. Nilai residual deviance (259,290) lebih kecil dari null deviance (412,815), menunjukkan bahwa prediktor yang dimasukkan memberikan perbaikan yang berarti dibanding model kosong yang hanya menggunakan intercept. Nilai AIC sebesar 327,290 menjadi acuan untuk dibandingkan dengan model reduksi; semakin kecil AIC, semakin baik keseimbangan antara kecocokan model dan kompleksitas parameter.
coef_table <- broom::tidy(heart_fit) %>%
filter(term != "(Intercept)") %>%
mutate(
odds_ratio = exp(estimate),
ci_low = exp(estimate - 1.96 * std.error),
ci_high = exp(estimate + 1.96 * std.error)
) %>%
arrange(p.value) %>%
transmute(
Variabel = term,
`Odds Ratio` = round(odds_ratio, 3),
`IK 95%` = paste0(round(ci_low, 3), " - ", round(ci_high, 3)),
`p-value` = signif(p.value, 3)
)
knitr::kable(coef_table,
caption = "Odds ratio dan p-value — Model Penuh") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed","responsive"), full_width = TRUE)| Variabel | Odds Ratio | IK 95% | p-value |
|---|---|---|---|
| locationRural | 13.478 | 6.24 - 29.11 | 0.00e+00 |
| age>=40 tahun | 0.053 | 0.015 - 0.182 | 3.10e-06 |
| period_birth>5 tahun | 6.119 | 2.055 - 18.223 | 1.14e-03 |
| religionIslam | 6.150 | 1.774 - 21.317 | 4.18e-03 |
| religionTraditional/Other | 5.511 | 1.504 - 20.201 | 1.00e-02 |
| occ_fatherLainnya | 0.185 | 0.049 - 0.694 | 1.24e-02 |
| age30-39 tahun | 0.329 | 0.129 - 0.841 | 2.02e-02 |
| period_birth3-5 tahun | 3.246 | 1.15 - 9.162 | 2.61e-02 |
| facilitySwasta | 0.299 | 0.097 - 0.927 | 3.64e-02 |
| period_birth1-3 tahun | 2.681 | 0.954 - 7.536 | 6.14e-02 |
| maritalJanda/Duda | 3.345 | 0.891 - 12.557 | 7.35e-02 |
| facilityHealth Centre | 0.463 | 0.166 - 1.291 | 1.41e-01 |
| occ_fatherPelajar | 0.448 | 0.153 - 1.317 | 1.45e-01 |
| occ_motherPelajar | 0.519 | 0.185 - 1.459 | 2.14e-01 |
| religionProtestant | 2.154 | 0.627 - 7.396 | 2.23e-01 |
| ethnicGa-Dangme | 1.934 | 0.626 - 5.976 | 2.52e-01 |
| facilityDistrik | 0.605 | 0.228 - 1.606 | 3.13e-01 |
| maritalCerai | 1.712 | 0.589 - 4.976 | 3.23e-01 |
| weight_birthTinggi (>=4kg) | 0.627 | 0.245 - 1.602 | 3.29e-01 |
| maritalKohabitasi | 1.618 | 0.561 - 4.663 | 3.73e-01 |
| occ_fatherPNS | 0.711 | 0.268 - 1.885 | 4.93e-01 |
| ethnicEwe | 0.699 | 0.251 - 1.949 | 4.94e-01 |
| ethnicMole-Dagbani | 0.706 | 0.257 - 1.945 | 5.01e-01 |
| age20-29 tahun | 0.746 | 0.309 - 1.804 | 5.16e-01 |
| maritalMenikah | 0.703 | 0.236 - 2.096 | 5.27e-01 |
| weight_birthNormal (2.5-3.9kg) | 0.766 | 0.333 - 1.763 | 5.31e-01 |
| occ_fatherPedagang | 0.783 | 0.273 - 2.248 | 6.49e-01 |
| facilityRegional | 0.827 | 0.308 - 2.225 | 7.07e-01 |
| ethnicGuan | 0.812 | 0.269 - 2.452 | 7.12e-01 |
| locationPeriurban | 1.135 | 0.384 - 3.353 | 8.19e-01 |
| occ_motherLainnya | 0.870 | 0.255 - 2.97 | 8.24e-01 |
| occ_motherPNS | 1.072 | 0.39 - 2.949 | 8.93e-01 |
| occ_motherPedagang | 0.952 | 0.334 - 2.713 | 9.27e-01 |
📊 Interpretasi output: Koefisien diurutkan berdasarkan p-value terkecil. Variabel dengan pengaruh paling signifikan adalah
locationRural(OR = 13,478; p < 0,001), yang berarti ibu yang tinggal di daerah Rural memiliki odds persalinan caesar 13,5 kali lebih tinggi dibandingkan ibu di daerah Urban. Kelompok usia ≥40 tahun memiliki OR = 0,053 (p < 0,001), artinya odds caesar pada kelompok ini hanya sekitar 5% dari kelompok usia <20 tahun (referensi) — jauh lebih rendah. Interval kepercayaan yang tidak melewati angka 1 dan p-value kecil memberikan bukti statistik yang lebih kuat bahwa prediktor tersebut relevan dalam model.
1.8 Seleksi Variabel: Stepwise Backward (AIC)
Pemilihan variabel dilakukan menggunakan metode stepwise
backward berbasis AIC (stepAIC). Dimulai dari
model penuh, variabel dihapus satu per satu selama penghapusannya
menurunkan AIC. Proses berhenti otomatis ketika tidak ada variabel yang
bila dihapus akan memperbaiki AIC.
cat("== Stepwise Backward Selection (AIC) ==\n")
heart_fitB <- stepAIC(heart_fit,
direction = "backward",
trace = TRUE)
cat("\n== Formula Model Reduksi (hasil stepAIC) ==\n")
print(formula(heart_fitB))## == Stepwise Backward Selection (AIC) ==
## Start: AIC=327.29
## caesar_bin ~ religion + age + location + period_birth + weight_birth +
## facility + marital + ethnic + occ_mother + occ_father
##
## Df Deviance AIC
## - occ_mother 4 262.14 322.14
## - ethnic 4 263.68 323.68
## - weight_birth 2 260.25 324.25
## - facility 4 265.12 325.12
## - marital 4 267.15 327.15
## <none> 259.29 327.29
## - occ_father 4 267.83 327.83
## - period_birth 3 271.15 333.15
## - religion 3 272.03 334.03
## - age 3 293.84 355.84
## - location 2 321.52 385.52
##
## Step: AIC=322.14
## caesar_bin ~ religion + age + location + period_birth + weight_birth +
## facility + marital + ethnic + occ_father
##
## Df Deviance AIC
## - ethnic 4 266.14 318.14
## - weight_birth 2 263.40 319.40
## - facility 4 267.49 319.49
## - marital 4 269.52 321.52
## <none> 262.14 322.14
## - occ_father 4 270.62 322.62
## - period_birth 3 272.60 326.60
## - religion 3 274.25 328.25
## - age 3 296.10 350.10
## - location 2 326.57 382.57
##
## Step: AIC=318.14
## caesar_bin ~ religion + age + location + period_birth + weight_birth +
## facility + marital + occ_father
##
## Df Deviance AIC
## - facility 4 271.14 315.14
## - weight_birth 2 267.80 315.80
## - marital 4 273.93 317.93
## - occ_father 4 273.97 317.97
## <none> 266.14 318.14
## - period_birth 3 275.61 321.61
## - religion 3 278.33 324.33
## - age 3 300.52 346.52
## - location 2 333.63 381.63
##
## Step: AIC=315.14
## caesar_bin ~ religion + age + location + period_birth + weight_birth +
## marital + occ_father
##
## Df Deviance AIC
## - weight_birth 2 272.21 312.21
## - occ_father 4 277.86 313.86
## - marital 4 278.94 314.94
## <none> 271.14 315.14
## - period_birth 3 280.32 318.32
## - religion 3 282.31 320.31
## - age 3 304.57 342.57
## - location 2 341.14 381.14
##
## Step: AIC=312.21
## caesar_bin ~ religion + age + location + period_birth + marital +
## occ_father
##
## Df Deviance AIC
## - occ_father 4 279.23 311.23
## - marital 4 279.76 311.76
## <none> 272.21 312.21
## - period_birth 3 281.51 315.51
## - religion 3 283.20 317.20
## - age 3 305.36 339.36
## - location 2 341.80 377.80
##
## Step: AIC=311.23
## caesar_bin ~ religion + age + location + period_birth + marital
##
## Df Deviance AIC
## - marital 4 286.41 310.41
## <none> 279.23 311.23
## - period_birth 3 289.36 315.36
## - religion 3 290.66 316.66
## - age 3 310.25 336.25
## - location 2 348.14 376.14
##
## Step: AIC=310.41
## caesar_bin ~ religion + age + location + period_birth
##
## Df Deviance AIC
## <none> 286.41 310.41
## - period_birth 3 295.84 313.84
## - religion 3 296.16 314.16
## - age 3 315.06 333.06
## - location 2 357.10 377.10
##
## == Formula Model Reduksi (hasil stepAIC) ==
## caesar_bin ~ religion + age + location + period_birthringkasan_2model <- data.frame(
Keterangan = c("Null deviance", "Residual deviance",
"Derajat bebas residual", "AIC"),
`Model Penuh` = c(
round(heart_fit$null.deviance, 3),
round(heart_fit$deviance, 3),
heart_fit$df.residual,
round(AIC(heart_fit), 3)
),
`Model Reduksi` = c(
round(heart_fitB$null.deviance, 3),
round(heart_fitB$deviance, 3),
heart_fitB$df.residual,
round(AIC(heart_fitB), 3)
),
check.names = FALSE
)
knitr::kable(ringkasan_2model,
caption = "Perbandingan ringkasan Model Penuh vs Reduksi") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Keterangan | Model Penuh | Model Reduksi |
|---|---|---|
| Null deviance | 412.815 | 412.815 |
| Residual deviance | 259.290 | 286.406 |
| Derajat bebas residual | 278.000 | 300.000 |
| AIC | 327.290 | 310.406 |
📊 Interpretasi output: Proses stepwise backward memilih model reduksi dengan prediktor
religion + age + location + period_birth. AIC model reduksi (310,406) lebih kecil dari model penuh (327,290), yang berarti model reduksi lebih parsimoni — menghasilkan kecocokan yang lebih baik relatif terhadap jumlah parameter yang digunakan.
1.9 Uji Asumsi
1.9.1 Multikolinearitas — VIF
Multikolinearitas diperiksa menggunakan Variance Inflation Factor (VIF). Untuk variabel kategorik digunakan Generalized VIF (GVIF).
## == VIF Model Penuh ==## GVIF Df GVIF^(1/(2*Df))
## religion 1.948777 3 1.117619
## age 1.687728 3 1.091148
## location 1.695294 2 1.141067
## period_birth 1.798096 3 1.102729
## weight_birth 1.428211 2 1.093196
## facility 1.819808 4 1.077713
## marital 1.791369 4 1.075593
## ethnic 1.936754 4 1.086136
## occ_mother 1.809523 4 1.076950
## occ_father 1.791945 4 1.075637##
## == VIF Model Reduksi ==## GVIF Df GVIF^(1/(2*Df))
## religion 1.244396 3 1.037114
## age 1.189696 3 1.029373
## location 1.212979 2 1.049454
## period_birth 1.112396 3 1.017911📊 Interpretasi output: Seluruh nilai GVIF^(1/(2·Df)) pada model penuh maupun model reduksi berada jauh di bawah 2. Karena nilai VIF < 5 (umumnya < 10 masih diterima), tidak ada indikasi multikolinearitas yang mengkhawatirkan. Artinya, setiap variabel prediktor memberikan informasi yang relatif independen satu sama lain.
1.9.2 Likelihood Ratio Test (LRT)
LRT digunakan untuk menguji apakah penambahan prediktor dalam model secara signifikan memperbaiki kecocokan.
\[ G^2 = -2\ln\left(\frac{L_{\text{model sederhana}}}{L_{\text{model kompleks}}}\right) = -2(\ell_{\text{sederhana}} - \ell_{\text{kompleks}}) \]
\(G^2\) mengikuti distribusi \(\chi^2\) dengan derajat bebas = selisih jumlah parameter antar model. Jika \(p\text{-value} < 0{,}05\), model kompleks secara signifikan lebih baik.
null_model <- glm(caesar_bin ~ 1, data = train_data, family = binomial)
cat("== LRT: Model Penuh vs Null ==\n")## == LRT: Model Penuh vs Null ==##
## == LRT: Model Reduksi vs Null ==##
## == LRT: Model Reduksi vs Model Penuh ==📊 Interpretasi output:
- LRT model penuh vs null: \(\chi^2\) = 153,53, p < 0,001 → model penuh secara signifikan lebih baik dari model null.
- LRT model reduksi vs null: \(\chi^2\) = 126,41, p < 0,001 →
model reduksi juga secara signifikan lebih baik dari model null.
- LRT model reduksi vs penuh: \(\chi^2\) = 27,116, p = 0,207 → tidak signifikan. Artinya, variabel-variabel yang dihapus oleh stepAIC tidak memberikan kontribusi berarti secara statistik. Model reduksi cukup dan dapat dipilih atas dasar kesederhanaan (parsimony).
1.9.3 Hosmer-Lemeshow Goodness-of-Fit
Uji Hosmer-Lemeshow mengevaluasi kecocokan model dengan membandingkan frekuensi observasi dan prediksi pada kelompok-kelompok peluang. Hipotesis:
- \(H_0\): Tidak ada perbedaan signifikan antara nilai yang diobservasi dan diprediksi model (model fit).
- \(H_1\): Ada perbedaan signifikan (model tidak fit).
Jika \(p > 0{,}05\), model dianggap fit.
## == Hosmer-Lemeshow — Model Penuh ==##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: train_data$caesar_bin, fitted(heart_fit)
## X-squared = 5.8796, df = 8, p-value = 0.6607##
## == Hosmer-Lemeshow — Model Reduksi ==hoslem_reduksi <- hoslem.test(train_data$caesar_bin, fitted(heart_fitB), g = 10)
print(hoslem_reduksi)##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: train_data$caesar_bin, fitted(heart_fitB)
## X-squared = 4.8908, df = 8, p-value = 0.7692📊 Interpretasi output:
- Model Penuh: \(\chi^2\) = 5,880, p = 0,661 →
\(p > 0{,}05\), model penuh
fit.
- Model Reduksi: \(\chi^2\) = 4,891, p = 0,769 → \(p > 0{,}05\), model reduksi juga fit.
Nilai p yang besar pada kedua model menunjukkan tidak ada perbedaan sistematis antara probabilitas yang diobservasi dan diprediksi.
1.9.4 Linearitas Log-Odds (Box-Tidwell)
Uji Box-Tidwell hanya relevan untuk prediktor kontinu. Karena seluruh prediktor dalam dataset ini bersifat kategorik, uji Box-Tidwell tidak diperlukan dan asumsi linearitas log-odds dianggap terpenuhi secara otomatis.
1.9.5 Nagelkerke \(R^2\)
Nagelkerke \(R^2\) merupakan ukuran pseudo-R² yang menunjukkan proporsi variasi yang dapat dijelaskan model, analog dengan \(R^2\) pada regresi linier.
\[ R^2_N = \frac{1 - \exp\!\left[\tfrac{2}{n}(\ell_{\text{null}} - \ell_{\text{model}})\right]}{1 - \exp\!\left[\tfrac{2}{n}\,\ell_{\text{null}}\right]} \]
n <- nrow(train_data)
nagelkerke <- function(fit, null) {
round((1 - exp((logLik(null) - logLik(fit)) * (2/n))) /
(1 - exp(logLik(null) * (2/n))), 4)
}
knitr::kable(
data.frame(
Model = c("Model Penuh", "Model Reduksi"),
`Nagelkerke R²` = c(
nagelkerke(heart_fit, null_model),
nagelkerke(heart_fitB, null_model)
),
check.names = FALSE
),
caption = "Nagelkerke R² kedua model"
) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Model | Nagelkerke R² |
|---|---|
| Model Penuh | 0.5297 |
| Model Reduksi | 0.4540 |
📊 Interpretasi output: Model penuh menjelaskan sekitar 52,97% variasi dalam status caesar, sementara model reduksi menjelaskan 45,40%. Penurunan \(R^2_N\) sebesar ~7,5 poin persentase ini sebanding dengan penghematan 22 parameter — konsisten dengan kesimpulan LRT bahwa model reduksi sudah cukup.
1.10 Odds Ratio Model Final (Model Reduksi)
or_tabel <- function(fit, judul) {
broom::tidy(fit) %>%
filter(term != "(Intercept)") %>%
mutate(
odds_ratio = exp(estimate),
ci_low = exp(estimate - 1.96 * std.error),
ci_high = exp(estimate + 1.96 * std.error)
) %>%
arrange(p.value) %>%
transmute(
Variabel = term,
`Odds Ratio` = round(odds_ratio, 3),
`IK 95%` = paste0(round(ci_low, 3), " - ", round(ci_high, 3)),
`p-value` = signif(p.value, 3)
) %>%
knitr::kable(caption = judul) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed","responsive"),
full_width = TRUE)
}
or_tabel(heart_fitB, "Odds Ratio — Model Reduksi / Model Final (stepAIC)")| Variabel | Odds Ratio | IK 95% | p-value |
|---|---|---|---|
| locationRural | 10.822 | 5.63 - 20.803 | 0.00e+00 |
| age>=40 tahun | 0.086 | 0.03 - 0.246 | 5.20e-06 |
| period_birth>5 tahun | 4.164 | 1.602 - 10.82 | 3.42e-03 |
| religionIslam | 4.598 | 1.546 - 13.677 | 6.09e-03 |
| religionTraditional/Other | 3.714 | 1.213 - 11.369 | 2.15e-02 |
| age30-39 tahun | 0.418 | 0.181 - 0.965 | 4.10e-02 |
| period_birth3-5 tahun | 2.480 | 1.015 - 6.061 | 4.63e-02 |
| period_birth1-3 tahun | 1.966 | 0.797 - 4.85 | 1.42e-01 |
| religionProtestant | 2.214 | 0.733 - 6.691 | 1.59e-01 |
| age20-29 tahun | 0.666 | 0.298 - 1.486 | 3.21e-01 |
| locationPeriurban | 1.007 | 0.383 - 2.647 | 9.89e-01 |
📊 Interpretasi output:
locationRural(OR = 10,822; IK 95%: 5,63–20,80; p < 0,001): Ibu di daerah Rural memiliki odds persalinan caesar 10,8 kali lebih tinggi dibandingkan ibu di Urban. Ini adalah prediktor paling kuat dalam model.age>=40 tahun(OR = 0,086; IK 95%: 0,03–0,246; p < 0,001): Ibu berusia ≥40 tahun memiliki odds caesar hanya 8,6% dari ibu berusia <20 tahun — penurunan odds yang sangat besar.age30-39 tahun(OR = 0,418; IK 95%: 0,181–0,965; p = 0,041): Ibu usia 30–39 tahun memiliki odds caesar sekitar 42% dari kelompok <20 tahun.period_birth>5 tahun(OR = 4,164; IK 95%: 1,60–10,82; p = 0,003): Jarak kelahiran >5 tahun meningkatkan odds caesar hingga 4,2 kali dibandingkan jarak <1 tahun.period_birth3-5 tahun(OR = 2,480; IK 95%: 1,015–6,061; p = 0,046): Jarak kelahiran 3–5 tahun meningkatkan odds caesar 2,5 kali dibandingkan jarak <1 tahun.religionIslam(OR = 4,598; IK 95%: 1,546–13,677; p = 0,006): Ibu beragama Islam memiliki odds caesar 4,6 kali lebih tinggi dibandingkan Catholic (referensi).religionTraditional/Other(OR = 3,714; IK 95%: 1,213–11,369; p = 0,021): Ibu dengan agama Traditional/Other memiliki odds caesar 3,7 kali lebih tinggi dari Catholic.
Prediktor locationPeriurban,
age20-29 tahun, period_birth1-3 tahun, dan
religionProtestant tidak signifikan secara statistik
(p > 0,05), namun tetap dipertahankan karena variabelnya
secara keseluruhan dipilih oleh stepAIC sebagai satu kesatuan
kategorik.
1.11 Prediksi Probabilitas
p_trainFull <- predict(heart_fit, newdata = train_data, type = "response")
p_trainB <- predict(heart_fitB, newdata = train_data, type = "response")
p_testFull <- predict(heart_fit, newdata = test_data, type = "response")
p_testB <- predict(heart_fitB, newdata = test_data, type = "response")
contoh_pred <- head(data.frame(
`Status aktual` = test_data$caesar_bin,
`Status (label)` = ifelse(test_data$caesar_bin == 1, "Caesar", "Tidak Caesar"),
`Peluang prediksi caesar` = round(p_testB, 4),
check.names = FALSE
), 8)
knitr::kable(contoh_pred, caption = "Contoh peluang prediksi pada data testing (Model Reduksi)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Status aktual | Status (label) | Peluang prediksi caesar |
|---|---|---|
| 0 | Tidak Caesar | 0.2584 |
| 1 | Caesar | 0.4280 |
| 0 | Tidak Caesar | 0.2153 |
| 1 | Caesar | 0.8120 |
| 0 | Tidak Caesar | 0.0542 |
| 0 | Tidak Caesar | 0.7892 |
| 1 | Caesar | 0.5835 |
| 0 | Tidak Caesar | 0.0226 |
📊 Interpretasi output: Tabel ini memperlihatkan beberapa contoh observasi testing beserta peluang prediksi caesar dari model reduksi. Semakin besar nilai peluang, semakin tinggi risiko ibu diklasifikasikan sebagai memiliki riwayat caesar. Peluang ini belum menjadi kelas akhir sampai dibandingkan dengan threshold.
1.12 Kurva ROC dan Threshold Optimal
1.12.1 Rumus Kurva ROC dan AUC
Kurva ROC (Receiver Operating Characteristic) mengevaluasi trade-off antara sensitivity dan specificity pada berbagai nilai threshold \(c\).
Untuk setiap threshold \(c\), dihitung:
\[ TPR(c) = \text{Sensitivity}(c) = \frac{TP(c)}{TP(c)+FN(c)} \]
\[ FPR(c) = 1 - \text{Specificity}(c) = \frac{FP(c)}{FP(c)+TN(c)} \]
Kurva ROC menggambarkan pasangan \((FPR(c),\, TPR(c))\) untuk berbagai nilai \(c\).
AUC (Area Under the Curve) adalah luas area di bawah kurva ROC:
\[ AUC = \int_0^1 TPR(FPR)\,d(FPR). \]
Semakin dekat AUC ke 1, semakin baik kemampuan model membedakan kelas. AUC = 0,5 setara dengan tebakan acak.
Threshold optimal dipilih dari data training menggunakan Indeks Youden:
\[ J = \text{Sensitivity} + \text{Specificity} - 1 \]
\[ c_{\text{optimal}} = \arg\max_c \left\{\text{Sensitivity}(c) + \text{Specificity}(c) - 1\right\} \]
safe_div <- function(num, den) ifelse(den == 0, NA_real_, num / den)
roc_points <- function(actual, prob) {
thresholds <- c(Inf, sort(unique(prob), decreasing = TRUE), -Inf)
out <- lapply(thresholds, function(th) {
pred <- as.integer(prob >= th)
tp <- sum(pred == 1 & actual == 1)
tn <- sum(pred == 0 & actual == 0)
fp <- sum(pred == 1 & actual == 0)
fn <- sum(pred == 0 & actual == 1)
data.frame(
threshold = th,
sensitivity = safe_div(tp, tp + fn),
specificity = safe_div(tn, tn + fp),
fpr = 1 - safe_div(tn, tn + fp),
youden = safe_div(tp, tp + fn) + safe_div(tn, tn + fp) - 1
)
})
bind_rows(out)
}
auc_value <- function(roc_df) {
roc_sorted <- roc_df %>% arrange(fpr, sensitivity)
sum(diff(roc_sorted$fpr) *
(head(roc_sorted$sensitivity, -1) +
tail(roc_sorted$sensitivity, -1)) / 2)
}roc_train_full <- roc_points(train_data$caesar_bin, p_trainFull) %>%
mutate(model = "Model Penuh - Training")
roc_train_B <- roc_points(train_data$caesar_bin, p_trainB) %>%
mutate(model = "Model Reduksi - Training")
roc_test_full <- roc_points(test_data$caesar_bin, p_testFull) %>%
mutate(model = "Model Penuh - Testing")
roc_test_B <- roc_points(test_data$caesar_bin, p_testB) %>%
mutate(model = "Model Reduksi - Testing")
knitr::kable(
data.frame(
Model = c("Model Penuh — Training", "Model Penuh — Testing",
"Model Reduksi — Training", "Model Reduksi — Testing"),
AUC = round(c(
auc_value(roc_train_full),
auc_value(roc_test_full),
auc_value(roc_train_B),
auc_value(roc_test_B)
), 3)
),
caption = "Nilai AUC kedua model"
) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Model | AUC |
|---|---|
| Model Penuh — Training | 0.880 |
| Model Penuh — Testing | 0.860 |
| Model Reduksi — Training | 0.852 |
| Model Reduksi — Testing | 0.870 |
📊 Interpretasi output: AUC model reduksi pada training (0,852) dan testing (0,870) sangat dekat satu sama lain, menandakan model tidak mengalami overfitting. AUC testing yang sedikit lebih tinggi mengindikasikan model mampu menggeneralisasi dengan baik ke data baru. Nilai AUC > 0,8 umumnya dikategorikan sebagai diskriminasi yang baik.
optimal_B <- roc_train_B %>%
filter(is.finite(threshold)) %>%
arrange(desc(youden), desc(sensitivity)) %>%
slice(1)
threshold_opt <- optimal_B$threshold[1]
knitr::kable(
optimal_B %>%
transmute(
`Threshold optimal` = round(threshold, 3),
Sensitivity = round(sensitivity, 3),
Specificity = round(specificity, 3),
`Indeks Youden` = round(youden, 3)
),
caption = "Threshold optimal Model Reduksi (Indeks Youden)"
) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Threshold optimal | Sensitivity | Specificity | Indeks Youden |
|---|---|---|---|
| 0.401 | 0.709 | 0.831 | 0.54 |
📊 Interpretasi output: Threshold optimal pada data training adalah 0,401, menghasilkan sensitivity 0,709 dan specificity 0,831 dengan indeks Youden = 0,540. Threshold ini dipilih karena memberikan keseimbangan terbaik antara kemampuan mendeteksi ibu caesar dan menjaga ibu tidak caesar tidak terlalu banyak salah terklasifikasi.
bind_rows(roc_test_full, roc_test_B) %>%
ggplot(aes(x = fpr, y = sensitivity, color = model)) +
geom_path(linewidth = 1.1) +
geom_abline(intercept = 0, slope = 1,
linetype = "dashed", color = "#6c757d") +
coord_equal() +
scale_color_manual(values = c(
"Model Penuh - Testing" = "#0077b6",
"Model Reduksi - Testing" = "#e76f51"
)) +
labs(
title = "Kurva ROC — Caesar Delivery (Testing)",
x = "False Positive Rate (1 - Specificity)",
y = "Sensitivity (True Positive Rate)",
color = "Model"
) +
theme_minimal(base_size = 12) +
theme(legend.position = "bottom")
Kurva ROC Model Penuh dan Reduksi pada Data Testing
💡 Interpretasi: Kurva ROC kedua model berada jauh di atas garis diagonal (tebakan acak), mengkonfirmasi kemampuan diskriminasi yang baik. Kurva model penuh dan reduksi pada data testing sangat berdekatan, yang konsisten dengan hasil LRT bahwa variabel yang dihapus tidak memberikan kontribusi signifikan.
1.13 Confusion Matrix dan Metrik Evaluasi
1.13.1 Rumus Metrik Klasifikasi
Untuk klasifikasi biner status caesar, notasi confusion matrix yang digunakan adalah:
- TP (True Positive): caesar aktual,
diprediksi caesar.
- TN (True Negative): tidak caesar aktual,
diprediksi tidak caesar.
- FP (False Positive): tidak caesar aktual,
diprediksi caesar.
- FN (False Negative): caesar aktual, diprediksi tidak caesar.
\[\text{Akurasi} = \frac{TP + TN}{TP + TN + FP + FN}\]
\[\text{Error rate} = 1 - \text{Akurasi}\]
\[\text{Sensitivity} = \text{Recall} = \text{TPR} = \frac{TP}{TP + FN}\]
\[\text{Specificity} = \text{TNR} = \frac{TN}{TN + FP}\]
\[\text{Presisi} = \frac{TP}{TP + FP}\]
\[\text{NPV} = \frac{TN}{TN + FN}\]
\[\text{F1-score} = 2 \times \frac{\text{Presisi} \times \text{Sensitivity}}{\text{Presisi} + \text{Sensitivity}}\]
\[\text{Balanced Accuracy} = \frac{\text{Sensitivity} + \text{Specificity}}{2}\]
Dalam konteks persalinan caesar, sensitivity penting karena menunjukkan kemampuan model mendeteksi ibu yang benar-benar memiliki riwayat caesar. Specificity menunjukkan kemampuan model mengenali ibu tanpa riwayat caesar agar tidak terlalu banyak yang keliru terklasifikasi.
classification_metrics <- function(actual, prob, threshold = 0.5) {
pred <- as.integer(prob >= threshold)
tp <- sum(pred == 1 & actual == 1)
tn <- sum(pred == 0 & actual == 0)
fp <- sum(pred == 1 & actual == 0)
fn <- sum(pred == 0 & actual == 1)
sensitivity <- safe_div(tp, tp + fn)
specificity <- safe_div(tn, tn + fp)
precision <- safe_div(tp, tp + fp)
npv <- safe_div(tn, tn + fn)
accuracy <- safe_div(tp + tn, tp + tn + fp + fn)
data.frame(
threshold = threshold,
accuracy = accuracy,
error_rate = 1 - accuracy,
sensitivity = sensitivity,
specificity = specificity,
precision = precision,
negative_predictive_value = npv,
f1_score = safe_div(2 * precision * sensitivity, precision + sensitivity),
balanced_accuracy = (sensitivity + specificity) / 2,
false_positive_rate = 1 - specificity,
false_negative_rate = 1 - sensitivity
)
}
format_metrics_indonesia <- function(x) {
x %>%
transmute(
Threshold = threshold,
Akurasi = accuracy,
`Error rate` = error_rate,
Sensitivity = sensitivity,
Specificity = specificity,
Presisi = precision,
NPV = negative_predictive_value,
`F1-score` = f1_score,
`Balanced accuracy` = balanced_accuracy,
FPR = false_positive_rate,
FNR = false_negative_rate
)
}
confusion_matrix_tbl <- function(actual, prob, threshold = 0.5) {
pred_label <- factor(
ifelse(prob >= threshold, "Prediksi Caesar", "Prediksi Tidak Caesar"),
levels = c("Prediksi Caesar", "Prediksi Tidak Caesar")
)
actual_label <- factor(
ifelse(actual == 1, "Aktual Caesar", "Aktual Tidak Caesar"),
levels = c("Aktual Caesar", "Aktual Tidak Caesar")
)
addmargins(table(actual_label, pred_label))
}1.13.2 Confusion Matrix — Threshold 0,50
cm_05 <- confusion_matrix_tbl(test_data$caesar_bin, p_testB, 0.5)
knitr::kable(cm_05,
caption = "Confusion Matrix — Model Reduksi, Threshold 0,50") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Prediksi Caesar | Prediksi Tidak Caesar | Sum | |
|---|---|---|---|
| Aktual Caesar | 17 | 13 | 30 |
| Aktual Tidak Caesar | 3 | 46 | 49 |
| Sum | 20 | 59 | 79 |
1.13.3 Confusion Matrix — Threshold Optimal
cm_opt <- confusion_matrix_tbl(test_data$caesar_bin, p_testB, threshold_opt)
knitr::kable(cm_opt,
caption = paste0("Confusion Matrix — Model Reduksi, Threshold Optimal (",
round(threshold_opt, 3), ")")) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Prediksi Caesar | Prediksi Tidak Caesar | Sum | |
|---|---|---|---|
| Aktual Caesar | 19 | 11 | 30 |
| Aktual Tidak Caesar | 7 | 42 | 49 |
| Sum | 26 | 53 | 79 |
1.13.4 Perbandingan Metrik
metrik_compare <- bind_rows(
classification_metrics(test_data$caesar_bin, p_testB, 0.5) %>%
mutate(aturan = "Threshold 0.50"),
classification_metrics(test_data$caesar_bin, p_testB, threshold_opt) %>%
mutate(aturan = "Threshold Optimal ROC")
) %>%
format_metrics_indonesia() %>%
mutate(
`Aturan klasifikasi` = c("Threshold 0.50", "Threshold Optimal ROC"),
across(where(is.numeric), round, 3)
) %>%
dplyr::select(`Aturan klasifikasi`, everything())
knitr::kable(metrik_compare,
caption = "Perbandingan metrik evaluasi Model Reduksi — Testing") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed","responsive"), full_width = TRUE)| Aturan klasifikasi | Threshold | Akurasi | Error rate | Sensitivity | Specificity | Presisi | NPV | F1-score | Balanced accuracy | FPR | FNR |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Threshold 0.50 | 0.500 | 0.797 | 0.203 | 0.567 | 0.939 | 0.850 | 0.780 | 0.680 | 0.753 | 0.061 | 0.433 |
| Threshold Optimal ROC | 0.401 | 0.772 | 0.228 | 0.633 | 0.857 | 0.731 | 0.792 | 0.679 | 0.745 | 0.143 | 0.367 |
📊 Interpretasi output:
- Threshold 0,50: Akurasi 79,7%, sensitivity 56,7%,
specificity 93,9%, F1-score 0,680.
- Threshold Optimal (0,401): Akurasi 77,2%, sensitivity 63,3%, specificity 85,7%, balanced accuracy 74,5%.
Pada threshold optimal, sensitivity meningkat dari 56,7% menjadi 63,3% — artinya lebih banyak ibu yang benar-benar caesar berhasil dideteksi. Konsekuensinya, specificity turun dari 93,9% menjadi 85,7% karena beberapa ibu tanpa riwayat caesar ikut terklasifikasi sebagai caesar. Pemilihan threshold bergantung pada prioritas: jika lebih penting mendeteksi semua riwayat caesar, threshold optimal lebih disarankan.
1.13.5 Perbandingan Model Penuh vs Reduksi
metrik_2model <- bind_rows(
classification_metrics(test_data$caesar_bin, p_testFull, 0.5) %>%
mutate(model = "Model Penuh"),
classification_metrics(test_data$caesar_bin, p_testB, 0.5) %>%
mutate(model = "Model Reduksi")
) %>%
dplyr::select(model, accuracy, sensitivity, specificity, f1_score, balanced_accuracy) %>%
mutate(across(where(is.numeric), round, 3)) %>%
rename(
Model = model,
Akurasi = accuracy,
Sensitivity = sensitivity,
Specificity = specificity,
`F1-score` = f1_score,
`Balanced Acc` = balanced_accuracy
)
knitr::kable(metrik_2model,
caption = "Perbandingan metrik Model Penuh vs Reduksi (threshold 0,5)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| Model | Akurasi | Sensitivity | Specificity | F1-score | Balanced Acc |
|---|---|---|---|---|---|
| Model Penuh | 0.810 | 0.667 | 0.898 | 0.727 | 0.782 |
| Model Reduksi | 0.797 | 0.567 | 0.939 | 0.680 | 0.753 |
📊 Interpretasi output: Model penuh memiliki akurasi (81,0%), sensitivity (66,7%), dan F1-score (0,727) yang sedikit lebih tinggi dari model reduksi. Namun perbedaan ini kecil dan tidak signifikan secara statistik (terkonfirmasi oleh LRT p = 0,207). Mengingat model reduksi jauh lebih sederhana (12 parameter vs 34 parameter) dengan AIC yang lebih rendah dan semua uji asumsi terpenuhi, model reduksi dipilih sebagai model final.
1.14 Distribusi Peluang Prediksi
test_data %>%
mutate(peluang_caesar = p_testB,
Status = ifelse(caesar_bin == 1, "Caesar", "Tidak Caesar")) %>%
ggplot(aes(x = peluang_caesar, fill = Status)) +
geom_density(alpha = 0.55, color = "white", linewidth = 0.7) +
geom_vline(xintercept = threshold_opt,
color = "#fb8500", linewidth = 1.2, linetype = "dashed") +
annotate(
"label",
x = threshold_opt, y = Inf,
label = paste0("threshold optimal = ", round(threshold_opt, 3)),
vjust = 1.4, fill = "#fff3b0", color = "#5f370e"
) +
scale_fill_manual(values = c("Tidak Caesar" = "#2a9d8f",
"Caesar" = "#e76f51")) +
labs(
title = "Distribusi Peluang Prediksi Caesar — Data Testing (Model Reduksi)",
x = "Peluang prediksi caesar",
y = "Kepadatan",
fill = "Status aktual"
) +
theme_minimal(base_size = 12)
Distribusi Peluang Prediksi Caesar pada Data Testing
💡 Interpretasi: Grafik kepadatan memperlihatkan sebaran peluang prediksi untuk ibu caesar (merah) dan tidak caesar (hijau). Kedua distribusi terpisah dengan cukup jelas: distribusi ibu tidak caesar terpusat di peluang rendah, sementara distribusi ibu caesar condong ke peluang lebih tinggi. Area tumpang-tindih antara kedua kurva menunjukkan bagian data yang sulit dibedakan oleh model — observasi di area inilah yang paling rentan salah klasifikasi. Garis vertikal menunjukkan threshold optimal; observasi di sebelah kanan garis diklasifikasikan sebagai caesar.
1.15 Kesimpulan
Regresi logistik biner berhasil dibangun untuk memodelkan probabilitas riwayat persalinan caesar berdasarkan karakteristik ibu dan fasilitas kesehatan. Berikut ringkasan temuan utama:
Pemilihan model: Stepwise backward (AIC) memilih empat prediktor utama —
religion,age,location, danperiod_birth— dengan AIC = 310,406, lebih parsimoni dari model penuh (AIC = 327,290). LRT mengkonfirmasi model reduksi sudah cukup (p = 0,207).Asumsi terpenuhi: Tidak ada multikolinearitas (VIF < 2), uji Hosmer-Lemeshow menunjukkan kedua model fit (p > 0,05), dan Box-Tidwell tidak diperlukan karena semua prediktor kategorik.
Prediktor paling berpengaruh:
- Lokasi Rural merupakan faktor risiko terkuat (OR = 10,8).
- Usia ≥40 tahun merupakan faktor protektif terkuat (OR = 0,086).
- Periode kelahiran >5 tahun meningkatkan odds caesar 4,2 kali.
- Agama Islam meningkatkan odds caesar 4,6 kali dibandingkan Catholic.
Performa model: AUC testing = 0,870 (diskriminasi baik). Pada threshold optimal 0,401, sensitivity = 63,3% dan specificity = 85,7%. Model tidak menunjukkan tanda overfitting (AUC training dan testing sangat dekat).
1.16 Catatan Metodologis
- Semua prediktor kategorik → uji Box-Tidwell tidak diperlukan.
- Hosmer-Lemeshow: \(p > 0{,}05\) → model fit.
- Pemilihan model via stepAIC (backward) lebih objektif daripada seleksi manual berdasarkan p-value per level.
- Baseline (kategori referensi) tiap faktor: Catholic, <20 tahun, Urban, <1 tahun (period_birth).
- Threshold klasifikasi dapat disesuaikan dengan kebutuhan klinis: threshold lebih rendah meningkatkan sensitivity namun menurunkan specificity.
2 Regresi Logistik Multinomial
2.1 Deskripsi Data
Dataset ini berasal dari Mendeley Data (DOI: 10.17632/xppzm3kv9g.2)
oleh Mondol & Kabir (2025). Data merupakan survei kesehatan mental
terstruktur yang memuat 1.998 respons dengan 21 fitur
demografis, gaya hidup, perilaku, dan psikologis. Variabel respons
adalah jumlah percobaan bunuh diri
(Suicide_Attempts) yang dikodekan sebagai: 0 = Tidak
Pernah, 1 = 1 Kali, 2 = 2 Kali, 3 = 3 Kali atau lebih. Kategori ini
bersifat nominal (tidak ada asumsi urutan atau jarak
yang sama antar kategori), menjadikan regresi logistik multinomial
sebagai metode yang sesuai.
data_multi <- read.csv("C:/Users/Salsabila Najwa/Downloads/Mental Health Classification.csv")
names(data_multi) <- trimws(names(data_multi))
data_multi <- data_multi %>%
mutate(
Suicide_Attempts = factor(
Suicide_Attempts,
levels = c(0, 1, 2, 3),
labels = c("Tidak Pernah", "1 Kali", "2 Kali", "3 Kali")
),
Gender = factor(Gender, levels = c(0, 1),
labels = c("Perempuan", "Laki-laki")),
Education_Level = factor(Education_Level,
levels = c(0, 1, 2, 3),
labels = c("SD", "SMP", "SMA", "Perguruan Tinggi")),
Employment_Status = factor(Employment_Status,
levels = c(0, 1, 2, 3, 4),
labels = c("Pengangguran", "Paruh Waktu",
"Penuh Waktu", "Pelajar", "Pensiunan")),
Low_Energy = factor(Low_Energy, levels = c(0, 1, 2),
labels = c("Tidak", "Kadang", "Sering")),
Low_SelfEsteem = factor(Low_SelfEsteem, levels = c(0, 1, 2),
labels = c("Tidak", "Kadang", "Sering")),
Search_Depression_Online = factor(Search_Depression_Online,
levels = c(0, 1),
labels = c("Tidak", "Ya")),
Worsening_Depression = factor(Worsening_Depression,
levels = c(0, 1),
labels = c("Tidak", "Ya")),
Self_Harm = factor(Self_Harm, levels = c(0, 1),
labels = c("Tidak", "Ya")),
Mental_Health_Support = factor(Mental_Health_Support,
levels = c(0, 1),
labels = c("Tidak", "Ya"))
) %>%
na.omit()
data_multi$Suicide_Attempts <- relevel(data_multi$Suicide_Attempts, ref = "Tidak Pernah")
cat("Referensi kategori:", levels(data_multi$Suicide_Attempts)[1], "\n")## Referensi kategori: Tidak Pernah## Jumlah baris: 1998## Rows: 1,998
## Columns: 21
## $ Gender <fct> Laki-laki, Laki-laki, Laki-laki, Perempuan, L…
## $ Age <int> 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 2…
## $ Education_Level <fct> SMA, SMA, SMA, SMA, SMA, SMA, SMA, SMP, SMA, …
## $ Employment_Status <fct> Pelajar, Penuh Waktu, Pelajar, Penuh Waktu, P…
## $ Depression_Type <int> 5, 5, 2, 6, 6, 9, 9, 4, 11, 2, 9, 9, 5, 2, 2,…
## $ Symptoms <int> 11, 0, 5, 7, 5, 1, 5, 3, 14, 5, 5, 1, 13, 5, …
## $ Low_Energy <fct> Kadang, Kadang, Kadang, Tidak, Tidak, Kadang,…
## $ Low_SelfEsteem <fct> Kadang, Kadang, Kadang, Kadang, Tidak, Tidak,…
## $ Search_Depression_Online <fct> Ya, Tidak, Ya, Tidak, Tidak, Tidak, Tidak, Ya…
## $ Worsening_Depression <fct> Ya, Ya, Ya, Ya, Tidak, Ya, Ya, Ya, Tidak, Tid…
## $ Your.overeating.level <int> 4, 4, 4, 2, 8, 9, 10, 11, 8, 8, 7, 11, 8, 8, …
## $ How.many.times.you.eat <int> 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, …
## $ SocialMedia_Hours <int> 10, 8, 10, 4, 3, 10, 8, 4, 10, 1, 5, 2, 5, 8,…
## $ SocialMedia_WhileEating <int> 3, 3, 3, 3, 3, 3, 2, 3, 2, 0, 2, 3, 0, 2, 1, …
## $ Sleep_Hours <int> 10, 4, 4, 3, 7, 7, 8, 7, 5, 5, 4, 6, 8, 5, 7,…
## $ Nervous_Level <int> 10, 10, 10, 10, 1, 5, 7, 8, 4, 2, 8, 2, 10, 5…
## $ Depression_Score <int> 10, 10, 10, 10, 3, 7, 8, 8, 3, 2, 9, 6, 6, 5,…
## $ Coping_Methods <int> 11, 0, 0, 5, 0, 0, 3, 3, 2, 9, 9, 13, 6, 3, 9…
## $ Self_Harm <fct> Tidak, Tidak, Ya, Ya, Tidak, Tidak, Tidak, Ti…
## $ Mental_Health_Support <fct> Tidak, Ya, Ya, Ya, Tidak, Tidak, Tidak, Tidak…
## $ Suicide_Attempts <fct> Tidak Pernah, Tidak Pernah, 1 Kali, 1 Kali, T…2.2 Distribusi Variabel Respons
data_multi %>%
count(Suicide_Attempts) %>%
mutate(Proporsi = percent(n / sum(n), accuracy = 0.1)) %>%
rename(`Percobaan Bunuh Diri` = Suicide_Attempts, Frekuensi = n) %>%
kable(caption = "Distribusi jumlah percobaan bunuh diri (variabel respons multinomial)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Percobaan Bunuh Diri | Frekuensi | Proporsi |
|---|---|---|
| Tidak Pernah | 559 | 28.0% |
| 1 Kali | 495 | 24.8% |
| 2 Kali | 461 | 23.1% |
| 3 Kali | 483 | 24.2% |
data_multi %>%
count(Suicide_Attempts) %>%
mutate(p = n / sum(n)) %>%
ggplot(aes(x = Suicide_Attempts, y = p, fill = Suicide_Attempts)) +
geom_col(width = 0.65, alpha = 0.94) +
geom_text(aes(label = percent(p, accuracy = 0.1)),
vjust = -0.4, fontface = "bold", color = "#243142") +
scale_y_continuous(labels = percent_format(), limits = c(0, 0.40)) +
scale_fill_manual(values = pal) +
labs(title = "Distribusi Jumlah Percobaan Bunuh Diri",
subtitle = "Respons multinomial: kategori tidak memiliki urutan alami.",
x = NULL, y = "Proporsi") +
theme_analisis() + theme(legend.position = "none")
interpretasi:Berdasarkan visualisasi bar chart
di atas, terlihat proporsi distribusi frekuensi percobaan bunuh diri
pada responden. Berikut adalah rinciannya:
- Kategori Mayoritas: Kelompok responden terbanyak adalah mereka yang “Tidak Pernah” melakukan percobaan bunuh diri, dengan proporsi sebesar 28.0%.
- Kategori Lainnya: Bagi responden yang memiliki
riwayat percobaan, distribusinya tersebar secara cukup merata pada
ketiga kategori frekuensi lainnya, yaitu:
- 1 Kali: 24.8%
- 3 Kali: 24.2%
- 2 Kali: 23.1% (proporsi terkecil)
Kesimpulan: Secara keseluruhan, sebaran data antarkategori terbilang cukup seimbang (balanced), di mana masing-masing kelompok berada pada rentang proporsi yang saling berdekatan (antara 23% hingga 28%). Selain itu, variabel dependen ini diperlakukan sebagai respons multinomial. Hal ini berarti dalam pemodelan regresi selanjutnya, kategori-kategori ini dianalisis tanpa mengasumsikan adanya tingkatan atau urutan matematis yang pasti (nominal), melainkan sebagai kelas-kelas independen yang saling lepas.
2.3 Uji Asumsi
2.3.1 Asumsi 1 — Respons Bersifat Nominal
Kategori percobaan bunuh diri (0, 1, 2, 3 kali) diperlakukan sebagai label nominal. Meskipun menggunakan angka, tidak ada asumsi bahwa jarak antara “1 Kali” dan “2 Kali” sama dengan antara “2 Kali” dan “3 Kali” secara klinis. Asumsi ini terpenuhi secara konseptual.
2.3.2 Asumsi 2 — Observasi Independen
Setiap baris mewakili satu individu berbeda tanpa struktur data berulang. Asumsi ini terpenuhi.
2.3.3 Asumsi 3 — Tidak Ada Multikolinearitas Berat (VIF)
num_vars <- c("Age", "Depression_Score", "Symptoms",
"SocialMedia_Hours", "Sleep_Hours", "Nervous_Level")
vif_manual <- sapply(num_vars, function(y_var) {
x_vars <- setdiff(num_vars, y_var)
formula <- as.formula(paste(y_var, "~", paste(x_vars, collapse = " + ")))
r2 <- summary(lm(formula, data = data_multi))$r.squared
1 / (1 - r2)
})
data.frame(
Variabel = names(vif_manual),
VIF = round(vif_manual, 3),
Keterangan = case_when(
vif_manual > 10 ~ "Multikolinearitas Berat",
vif_manual > 5 ~ "Perlu Diperhatikan",
TRUE ~ "Aman"
)
) %>%
arrange(desc(VIF)) %>%
kable(caption = "Uji VIF — Variabel Numerik (Multinomial)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Variabel | VIF | Keterangan | |
|---|---|---|---|
| Symptoms | Symptoms | 1.019 | Aman |
| Nervous_Level | Nervous_Level | 1.017 | Aman |
| Sleep_Hours | Sleep_Hours | 1.003 | Aman |
| Depression_Score | Depression_Score | 1.002 | Aman |
| SocialMedia_Hours | SocialMedia_Hours | 1.002 | Aman |
| Age | Age | 1.001 | Aman |
2.3.4 Asumsi 4 — Tidak Ada Perfect Separation
fit_check_m <- nnet::multinom(
Suicide_Attempts ~ Age + Depression_Score + Symptoms +
SocialMedia_Hours + Sleep_Hours + Nervous_Level +
Gender + Education_Level + Employment_Status +
Low_Energy + Low_SelfEsteem + Search_Depression_Online +
Worsening_Depression + Self_Harm + Mental_Health_Support,
data = data_multi, trace = FALSE, maxit = 500
)
coef_ext <- coef(fit_check_m)
extreme <- which(abs(coef_ext) > 10, arr.ind = TRUE)
if (length(extreme) == 0) {
cat("Tidak ditemukan koefisien ekstrem (|β| > 10). Tidak ada indikasi perfect separation.\n")
} else {
cat("Ditemukan koefisien ekstrem:\n"); print(extreme)
}## Tidak ditemukan koefisien ekstrem (|β| > 10). Tidak ada indikasi perfect separation.2.3.5 Asumsi 5 — Linearitas pada Skala Logit (Box-Tidwell)
cont_vars_m <- c("Age", "Depression_Score", "Symptoms",
"SocialMedia_Hours", "Sleep_Hours", "Nervous_Level")
data_bt_m <- data_multi
for (v in cont_vars_m) {
data_bt_m[[paste0(v, "_log")]] <- data_bt_m[[v]] * log(data_bt_m[[v]] + 1)
}
bt_formula_m <- as.formula(paste(
"Suicide_Attempts ~",
paste(c(cont_vars_m, paste0(cont_vars_m, "_log"),
"Gender", "Education_Level", "Employment_Status",
"Low_Energy", "Low_SelfEsteem", "Search_Depression_Online",
"Worsening_Depression", "Self_Harm", "Mental_Health_Support"),
collapse = " + ")
))
fit_bt_m <- nnet::multinom(bt_formula_m, data = data_bt_m, trace = FALSE, maxit = 500)
bt_sum_m <- summary(fit_bt_m)
bt_coef_m <- as.data.frame(bt_sum_m$coefficients)
bt_se_m <- as.data.frame(bt_sum_m$standard.errors)
bt_long <- bt_coef_m %>%
rownames_to_column("kategori") %>%
pivot_longer(-kategori, names_to = "variabel", values_to = "estimate") %>%
left_join(
bt_se_m %>% rownames_to_column("kategori") %>%
pivot_longer(-kategori, names_to = "variabel", values_to = "se"),
by = c("kategori", "variabel")
) %>%
mutate(z = estimate / se, p = 2 * (1 - pnorm(abs(z)))) %>%
filter(grepl("_log$", variabel))
bt_long %>%
mutate(across(c(estimate, se, z, p), ~round(.x, 4)),
Keterangan = ifelse(p < 0.05, "Tidak Linear (p < 0.05)", "Linear (p ≥ 0.05)")) %>%
rename(Kategori = kategori, Variabel = variabel,
Estimate = estimate, SE = se, `z-value` = z, `p-value` = p) %>%
kable(caption = "Uji Box-Tidwell — Multinomial") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| Kategori | Variabel | Estimate | SE | z-value | p-value | Keterangan |
|---|---|---|---|---|---|---|
| 1 Kali | Age_log | -0.0158 | 0.0383 | -0.4116 | 0.6807 | Linear (p ≥ 0.05) |
| 1 Kali | Depression_Score_log | 0.0271 | 0.0199 | 1.3569 | 0.1748 | Linear (p ≥ 0.05) |
| 1 Kali | Symptoms_log | -0.0067 | 0.0519 | -0.1288 | 0.8975 | Linear (p ≥ 0.05) |
| 1 Kali | SocialMedia_Hours_log | 0.0307 | 0.0491 | 0.6246 | 0.5323 | Linear (p ≥ 0.05) |
| 1 Kali | Sleep_Hours_log | 0.2522 | 0.1756 | 1.4358 | 0.1510 | Linear (p ≥ 0.05) |
| 1 Kali | Nervous_Level_log | -0.0901 | 0.0920 | -0.9791 | 0.3275 | Linear (p ≥ 0.05) |
| 2 Kali | Age_log | 0.0470 | 0.0385 | 1.2203 | 0.2223 | Linear (p ≥ 0.05) |
| 2 Kali | Depression_Score_log | -0.0007 | 0.0211 | -0.0316 | 0.9748 | Linear (p ≥ 0.05) |
| 2 Kali | Symptoms_log | -0.0060 | 0.0532 | -0.1134 | 0.9097 | Linear (p ≥ 0.05) |
| 2 Kali | SocialMedia_Hours_log | 0.0526 | 0.0497 | 1.0579 | 0.2901 | Linear (p ≥ 0.05) |
| 2 Kali | Sleep_Hours_log | 0.1489 | 0.1787 | 0.8332 | 0.4047 | Linear (p ≥ 0.05) |
| 2 Kali | Nervous_Level_log | -0.0913 | 0.0993 | -0.9198 | 0.3577 | Linear (p ≥ 0.05) |
| 3 Kali | Age_log | 0.0350 | 0.0382 | 0.9175 | 0.3589 | Linear (p ≥ 0.05) |
| 3 Kali | Depression_Score_log | 0.0037 | 0.0205 | 0.1815 | 0.8560 | Linear (p ≥ 0.05) |
| 3 Kali | Symptoms_log | 0.0171 | 0.0519 | 0.3296 | 0.7417 | Linear (p ≥ 0.05) |
| 3 Kali | SocialMedia_Hours_log | 0.0893 | 0.0487 | 1.8328 | 0.0668 | Linear (p ≥ 0.05) |
| 3 Kali | Sleep_Hours_log | 0.5281 | 0.1740 | 3.0357 | 0.0024 | Tidak Linear (p < 0.05) |
| 3 Kali | Nervous_Level_log | -0.1062 | 0.0940 | -1.1300 | 0.2585 | Linear (p ≥ 0.05) |
2.3.6 Asumsi 6 — Independence of Irrelevant Alternatives (IIA)
Setiap pasangan kategori percobaan bunuh diri merepresentasikan outcome klinis yang secara substantif berbeda. Tidak ada indikasi bahwa kehadiran atau ketiadaan satu kategori akan mengubah odds relatif antara dua kategori lainnya. Asumsi IIA terpenuhi secara konseptual.
2.4 Estimasi Model
fit_multi <- nnet::multinom(
Suicide_Attempts ~ Age + Depression_Score + Symptoms +
SocialMedia_Hours + Sleep_Hours + Nervous_Level +
Gender + Education_Level + Employment_Status +
Low_Energy + Low_SelfEsteem + Search_Depression_Online +
Worsening_Depression + Self_Harm + Mental_Health_Support,
data = data_multi,
trace = FALSE,
maxit = 500
)
summary(fit_multi)## Call:
## nnet::multinom(formula = Suicide_Attempts ~ Age + Depression_Score +
## Symptoms + SocialMedia_Hours + Sleep_Hours + Nervous_Level +
## Gender + Education_Level + Employment_Status + Low_Energy +
## Low_SelfEsteem + Search_Depression_Online + Worsening_Depression +
## Self_Harm + Mental_Health_Support, data = data_multi, trace = FALSE,
## maxit = 500)
##
## Coefficients:
## (Intercept) Age Depression_Score Symptoms SocialMedia_Hours
## 1 Kali -0.52201844 -0.001039709 -0.001512692 -0.005492664 0.013402636
## 2 Kali -0.02659206 0.002275849 0.010800912 0.020540527 0.003738022
## 3 Kali 0.40741320 -0.002289252 0.002080553 0.009941910 0.007301051
## Sleep_Hours Nervous_Level GenderLaki-laki Education_LevelSMP
## 1 Kali 0.01790473 0.002439378 -0.05593599 0.6922532
## 2 Kali -0.02233086 0.064942771 0.13699636 -0.2023597
## 3 Kali -0.01927725 -0.011582379 0.15368084 0.1912035
## Education_LevelSMA Education_LevelPerguruan Tinggi
## 1 Kali 0.286773893 -0.9418570
## 2 Kali -0.229627637 -0.6105408
## 3 Kali -0.002501825 -0.6110936
## Employment_StatusParuh Waktu Employment_StatusPenuh Waktu
## 1 Kali -0.1470176 -0.04859015
## 2 Kali -0.3315887 -0.72543500
## 3 Kali -0.3212932 -0.49188875
## Employment_StatusPelajar Employment_StatusPensiunan Low_EnergyKadang
## 1 Kali -0.3844465 -0.2718530 0.078752472
## 2 Kali -0.8256868 -0.8146593 0.004060511
## 3 Kali -0.7910615 -0.4707275 -0.025609757
## Low_EnergySering Low_SelfEsteemKadang Low_SelfEsteemSering
## 1 Kali 0.1782386 0.02727663 -0.1490788
## 2 Kali -0.7514141 0.01871303 -2.8733099
## 3 Kali 0.6788065 -0.07268337 -1.1086254
## Search_Depression_OnlineYa Worsening_DepressionYa Self_HarmYa
## 1 Kali 0.17570942 -0.15744146 0.2600122
## 2 Kali 0.08718718 0.11181122 -0.3720007
## 3 Kali -0.03447103 0.02912569 -0.2975029
## Mental_Health_SupportYa
## 1 Kali 0.16909547
## 2 Kali 0.03103914
## 3 Kali 0.18637308
##
## Std. Errors:
## (Intercept) Age Depression_Score Symptoms SocialMedia_Hours
## 1 Kali 0.6789573 0.007859415 0.006983756 0.01850682 0.01671063
## 2 Kali 0.6579760 0.008043611 0.007104825 0.01840619 0.01697956
## 3 Kali 0.6635840 0.007876796 0.007002876 0.01866766 0.01680727
## Sleep_Hours Nervous_Level GenderLaki-laki Education_LevelSMP
## 1 Kali 0.02745599 0.03423221 0.1246034 0.5254511
## 2 Kali 0.02796110 0.03499556 0.1273167 0.5131668
## 3 Kali 0.02760008 0.03475881 0.1254921 0.5128235
## Education_LevelSMA Education_LevelPerguruan Tinggi
## 1 Kali 0.4634991 0.8840247
## 2 Kali 0.4328037 0.7249215
## 3 Kali 0.4439645 0.7326779
## Employment_StatusParuh Waktu Employment_StatusPenuh Waktu
## 1 Kali 0.7029243 0.3721877
## 2 Kali 0.6445642 0.3661353
## 3 Kali 0.6606192 0.3635101
## Employment_StatusPelajar Employment_StatusPensiunan Low_EnergyKadang
## 1 Kali 0.4085262 0.4288994 0.1449799
## 2 Kali 0.4010624 0.4227477 0.1492561
## 3 Kali 0.4022852 0.4145291 0.1469093
## Low_EnergySering Low_SelfEsteemKadang Low_SelfEsteemSering
## 1 Kali 0.9680429 0.1701792 0.6615911
## 2 Kali 0.9234058 0.1789822 1.1630823
## 3 Kali 0.8710105 0.1718529 0.6946210
## Search_Depression_OnlineYa Worsening_DepressionYa Self_HarmYa
## 1 Kali 0.1368546 0.1819007 0.2515295
## 2 Kali 0.1403579 0.1867209 0.2701093
## 3 Kali 0.1406419 0.1851493 0.2717444
## Mental_Health_SupportYa
## 1 Kali 0.1912115
## 2 Kali 0.2014442
## 3 Kali 0.1945859
##
## Residual Deviance: 5452.662
## AIC: 5590.6622.5 Tabel Koefisien, RRR, dan Interval Kepercayaan
multi_sum <- summary(fit_multi)
coef_m <- as.data.frame(multi_sum$coefficients)
se_m <- as.data.frame(multi_sum$standard.errors)
result_multi <- coef_m %>%
rownames_to_column("kategori") %>%
pivot_longer(-kategori, names_to = "variabel", values_to = "estimate") %>%
left_join(
se_m %>% rownames_to_column("kategori") %>%
pivot_longer(-kategori, names_to = "variabel", values_to = "se"),
by = c("kategori", "variabel")
) %>%
mutate(
z = estimate / se,
p_value = 2 * (1 - pnorm(abs(z))),
RRR = exp(estimate),
CI_low = exp(estimate - 1.96 * se),
CI_high = exp(estimate + 1.96 * se)
) %>%
mutate(across(c(estimate, se, z, p_value, RRR, CI_low, CI_high), ~round(.x, 4)))
result_multi %>%
kable(
caption = "Ringkasan hasil regresi logistik multinomial (RRR = Relative Risk Ratio)",
col.names = c("Kategori vs Referensi", "Variabel", "Estimate",
"SE", "z-value", "p-value", "RRR", "CI 2.5%", "CI 97.5%")
) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| Kategori vs Referensi | Variabel | Estimate | SE | z-value | p-value | RRR | CI 2.5% | CI 97.5% |
|---|---|---|---|---|---|---|---|---|
| 1 Kali | (Intercept) | -0.5220 | 0.6790 | -0.7689 | 0.4420 | 0.5933 | 0.1568 | 2.2451 |
| 1 Kali | Age | -0.0010 | 0.0079 | -0.1323 | 0.8948 | 0.9990 | 0.9837 | 1.0145 |
| 1 Kali | Depression_Score | -0.0015 | 0.0070 | -0.2166 | 0.8285 | 0.9985 | 0.9849 | 1.0122 |
| 1 Kali | Symptoms | -0.0055 | 0.0185 | -0.2968 | 0.7666 | 0.9945 | 0.9591 | 1.0313 |
| 1 Kali | SocialMedia_Hours | 0.0134 | 0.0167 | 0.8020 | 0.4225 | 1.0135 | 0.9808 | 1.0472 |
| 1 Kali | Sleep_Hours | 0.0179 | 0.0275 | 0.6521 | 0.5143 | 1.0181 | 0.9647 | 1.0744 |
| 1 Kali | Nervous_Level | 0.0024 | 0.0342 | 0.0713 | 0.9432 | 1.0024 | 0.9374 | 1.0720 |
| 1 Kali | GenderLaki-laki | -0.0559 | 0.1246 | -0.4489 | 0.6535 | 0.9456 | 0.7407 | 1.2072 |
| 1 Kali | Education_LevelSMP | 0.6923 | 0.5255 | 1.3174 | 0.1877 | 1.9982 | 0.7135 | 5.5965 |
| 1 Kali | Education_LevelSMA | 0.2868 | 0.4635 | 0.6187 | 0.5361 | 1.3321 | 0.5370 | 3.3043 |
| 1 Kali | Education_LevelPerguruan Tinggi | -0.9419 | 0.8840 | -1.0654 | 0.2867 | 0.3899 | 0.0689 | 2.2052 |
| 1 Kali | Employment_StatusParuh Waktu | -0.1470 | 0.7029 | -0.2092 | 0.8343 | 0.8633 | 0.2177 | 3.4237 |
| 1 Kali | Employment_StatusPenuh Waktu | -0.0486 | 0.3722 | -0.1306 | 0.8961 | 0.9526 | 0.4593 | 1.9757 |
| 1 Kali | Employment_StatusPelajar | -0.3844 | 0.4085 | -0.9411 | 0.3467 | 0.6808 | 0.3057 | 1.5163 |
| 1 Kali | Employment_StatusPensiunan | -0.2719 | 0.4289 | -0.6338 | 0.5262 | 0.7620 | 0.3287 | 1.7661 |
| 1 Kali | Low_EnergyKadang | 0.0788 | 0.1450 | 0.5432 | 0.5870 | 1.0819 | 0.8143 | 1.4375 |
| 1 Kali | Low_EnergySering | 0.1782 | 0.9680 | 0.1841 | 0.8539 | 1.1951 | 0.1792 | 7.9693 |
| 1 Kali | Low_SelfEsteemKadang | 0.0273 | 0.1702 | 0.1603 | 0.8727 | 1.0277 | 0.7362 | 1.4345 |
| 1 Kali | Low_SelfEsteemSering | -0.1491 | 0.6616 | -0.2253 | 0.8217 | 0.8615 | 0.2356 | 3.1507 |
| 1 Kali | Search_Depression_OnlineYa | 0.1757 | 0.1369 | 1.2839 | 0.1992 | 1.1921 | 0.9116 | 1.5588 |
| 1 Kali | Worsening_DepressionYa | -0.1574 | 0.1819 | -0.8655 | 0.3867 | 0.8543 | 0.5981 | 1.2203 |
| 1 Kali | Self_HarmYa | 0.2600 | 0.2515 | 1.0337 | 0.3013 | 1.2969 | 0.7922 | 2.1234 |
| 1 Kali | Mental_Health_SupportYa | 0.1691 | 0.1912 | 0.8843 | 0.3765 | 1.1842 | 0.8141 | 1.7227 |
| 2 Kali | (Intercept) | -0.0266 | 0.6580 | -0.0404 | 0.9678 | 0.9738 | 0.2681 | 3.5362 |
| 2 Kali | Age | 0.0023 | 0.0080 | 0.2829 | 0.7772 | 1.0023 | 0.9866 | 1.0182 |
| 2 Kali | Depression_Score | 0.0108 | 0.0071 | 1.5202 | 0.1285 | 1.0109 | 0.9969 | 1.0250 |
| 2 Kali | Symptoms | 0.0205 | 0.0184 | 1.1160 | 0.2644 | 1.0208 | 0.9846 | 1.0583 |
| 2 Kali | SocialMedia_Hours | 0.0037 | 0.0170 | 0.2201 | 0.8258 | 1.0037 | 0.9709 | 1.0377 |
| 2 Kali | Sleep_Hours | -0.0223 | 0.0280 | -0.7986 | 0.4245 | 0.9779 | 0.9258 | 1.0330 |
| 2 Kali | Nervous_Level | 0.0649 | 0.0350 | 1.8557 | 0.0635 | 1.0671 | 0.9964 | 1.1429 |
| 2 Kali | GenderLaki-laki | 0.1370 | 0.1273 | 1.0760 | 0.2819 | 1.1468 | 0.8936 | 1.4719 |
| 2 Kali | Education_LevelSMP | -0.2024 | 0.5132 | -0.3943 | 0.6933 | 0.8168 | 0.2987 | 2.2332 |
| 2 Kali | Education_LevelSMA | -0.2296 | 0.4328 | -0.5306 | 0.5957 | 0.7948 | 0.3403 | 1.8565 |
| 2 Kali | Education_LevelPerguruan Tinggi | -0.6105 | 0.7249 | -0.8422 | 0.3997 | 0.5431 | 0.1312 | 2.2486 |
| 2 Kali | Employment_StatusParuh Waktu | -0.3316 | 0.6446 | -0.5144 | 0.6069 | 0.7178 | 0.2029 | 2.5390 |
| 2 Kali | Employment_StatusPenuh Waktu | -0.7254 | 0.3661 | -1.9813 | 0.0476 | 0.4841 | 0.2362 | 0.9922 |
| 2 Kali | Employment_StatusPelajar | -0.8257 | 0.4011 | -2.0587 | 0.0395 | 0.4379 | 0.1995 | 0.9612 |
| 2 Kali | Employment_StatusPensiunan | -0.8147 | 0.4227 | -1.9271 | 0.0540 | 0.4428 | 0.1934 | 1.0140 |
| 2 Kali | Low_EnergyKadang | 0.0041 | 0.1493 | 0.0272 | 0.9783 | 1.0041 | 0.7494 | 1.3453 |
| 2 Kali | Low_EnergySering | -0.7514 | 0.9234 | -0.8137 | 0.4158 | 0.4717 | 0.0772 | 2.8819 |
| 2 Kali | Low_SelfEsteemKadang | 0.0187 | 0.1790 | 0.1046 | 0.9167 | 1.0189 | 0.7174 | 1.4470 |
| 2 Kali | Low_SelfEsteemSering | -2.8733 | 1.1631 | -2.4704 | 0.0135 | 0.0565 | 0.0058 | 0.5523 |
| 2 Kali | Search_Depression_OnlineYa | 0.0872 | 0.1404 | 0.6212 | 0.5345 | 1.0911 | 0.8287 | 1.4366 |
| 2 Kali | Worsening_DepressionYa | 0.1118 | 0.1867 | 0.5988 | 0.5493 | 1.1183 | 0.7756 | 1.6125 |
| 2 Kali | Self_HarmYa | -0.3720 | 0.2701 | -1.3772 | 0.1684 | 0.6894 | 0.4060 | 1.1705 |
| 2 Kali | Mental_Health_SupportYa | 0.0310 | 0.2014 | 0.1541 | 0.8775 | 1.0315 | 0.6950 | 1.5309 |
| 3 Kali | (Intercept) | 0.4074 | 0.6636 | 0.6140 | 0.5392 | 1.5029 | 0.4093 | 5.5181 |
| 3 Kali | Age | -0.0023 | 0.0079 | -0.2906 | 0.7713 | 0.9977 | 0.9824 | 1.0132 |
| 3 Kali | Depression_Score | 0.0021 | 0.0070 | 0.2971 | 0.7664 | 1.0021 | 0.9884 | 1.0159 |
| 3 Kali | Symptoms | 0.0099 | 0.0187 | 0.5326 | 0.5943 | 1.0100 | 0.9737 | 1.0476 |
| 3 Kali | SocialMedia_Hours | 0.0073 | 0.0168 | 0.4344 | 0.6640 | 1.0073 | 0.9747 | 1.0411 |
| 3 Kali | Sleep_Hours | -0.0193 | 0.0276 | -0.6984 | 0.4849 | 0.9809 | 0.9293 | 1.0354 |
| 3 Kali | Nervous_Level | -0.0116 | 0.0348 | -0.3332 | 0.7390 | 0.9885 | 0.9234 | 1.0582 |
| 3 Kali | GenderLaki-laki | 0.1537 | 0.1255 | 1.2246 | 0.2207 | 1.1661 | 0.9118 | 1.4913 |
| 3 Kali | Education_LevelSMP | 0.1912 | 0.5128 | 0.3728 | 0.7093 | 1.2107 | 0.4431 | 3.3080 |
| 3 Kali | Education_LevelSMA | -0.0025 | 0.4440 | -0.0056 | 0.9955 | 0.9975 | 0.4178 | 2.3814 |
| 3 Kali | Education_LevelPerguruan Tinggi | -0.6111 | 0.7327 | -0.8341 | 0.4042 | 0.5428 | 0.1291 | 2.2818 |
| 3 Kali | Employment_StatusParuh Waktu | -0.3213 | 0.6606 | -0.4864 | 0.6267 | 0.7252 | 0.1987 | 2.6472 |
| 3 Kali | Employment_StatusPenuh Waktu | -0.4919 | 0.3635 | -1.3532 | 0.1760 | 0.6115 | 0.2999 | 1.2468 |
| 3 Kali | Employment_StatusPelajar | -0.7911 | 0.4023 | -1.9664 | 0.0493 | 0.4534 | 0.2061 | 0.9974 |
| 3 Kali | Employment_StatusPensiunan | -0.4707 | 0.4145 | -1.1356 | 0.2561 | 0.6245 | 0.2771 | 1.4074 |
| 3 Kali | Low_EnergyKadang | -0.0256 | 0.1469 | -0.1743 | 0.8616 | 0.9747 | 0.7308 | 1.3000 |
| 3 Kali | Low_EnergySering | 0.6788 | 0.8710 | 0.7793 | 0.4358 | 1.9715 | 0.3576 | 10.8698 |
| 3 Kali | Low_SelfEsteemKadang | -0.0727 | 0.1719 | -0.4229 | 0.6723 | 0.9299 | 0.6640 | 1.3023 |
| 3 Kali | Low_SelfEsteemSering | -1.1086 | 0.6946 | -1.5960 | 0.1105 | 0.3300 | 0.0846 | 1.2877 |
| 3 Kali | Search_Depression_OnlineYa | -0.0345 | 0.1406 | -0.2451 | 0.8064 | 0.9661 | 0.7334 | 1.2728 |
| 3 Kali | Worsening_DepressionYa | 0.0291 | 0.1851 | 0.1573 | 0.8750 | 1.0296 | 0.7162 | 1.4800 |
| 3 Kali | Self_HarmYa | -0.2975 | 0.2717 | -1.0948 | 0.2736 | 0.7427 | 0.4360 | 1.2651 |
| 3 Kali | Mental_Health_SupportYa | 0.1864 | 0.1946 | 0.9578 | 0.3382 | 1.2049 | 0.8228 | 1.7643 |
Panduan Interpretasi RRR — Regresi Logistik Multinomial:
Kategori referensi adalah “Tidak Pernah”. Model memiliki \((K-1) = 3\) persamaan log-odds (untuk 1 Kali, 2 Kali, 3 Kali masing-masing vs Tidak Pernah): \[\log\!\left(\frac{P(Y = k)}{P(Y = \text{Referensi})}\right) = \beta_{k0} + \beta_{k1} X_1 + \cdots + \beta_{kp} X_p, \quad k \in \{1\text{ Kali}, 2\text{ Kali}, 3\text{ Kali}\}\]
RRR (Relative Risk Ratio) = exp(β): - RRR > 1: Kenaikan satu satuan prediktor meningkatkan odds relatif masuk ke kategori tersebut dibanding referensi. - RRR < 1: Kenaikan satu satuan prediktor menurunkan odds relatif masuk ke kategori tersebut. - Signifikan jika p-value < 0.05 atau CI 95% tidak mencakup angka 1.
Contoh: Jika
Self_Harm = Yapada baris “3 Kali” memiliki RRR = 4.5 (p < 0.05), maka individu yang pernah menyakiti diri sendiri memiliki odds 4,5 kali lebih besar untuk melakukan percobaan bunuh diri 3 kali dibandingkan tidak pernah sama sekali, dengan variabel lain konstan.
2.6 Prediksi dan Confusion Matrix
pred_prob_m <- predict(fit_multi, type = "probs")
pred_class_m <- predict(fit_multi, type = "class")
data_pred_m <- data_multi %>%
bind_cols(as.data.frame(pred_prob_m)) %>%
mutate(prediksi = pred_class_m)
conf_m <- table(Aktual = data_pred_m$Suicide_Attempts, Prediksi = pred_class_m)
kable(addmargins(conf_m), caption = "Confusion Matrix — Model Multinomial") %>%
kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)| Tidak Pernah | 1 Kali | 2 Kali | 3 Kali | Sum | |
|---|---|---|---|---|---|
| Tidak Pernah | 266 | 127 | 103 | 63 | 559 |
| 1 Kali | 215 | 156 | 82 | 42 | 495 |
| 2 Kali | 211 | 88 | 117 | 45 | 461 |
| 3 Kali | 230 | 106 | 88 | 59 | 483 |
| Sum | 922 | 477 | 390 | 209 | 1998 |
acc_m <- sum(diag(conf_m)) / sum(conf_m)
cat("Akurasi model multinomial:", round(acc_m * 100, 2), "%\n")## Akurasi model multinomial: 29.93 %2.7 Visualisasi Probabilitas Prediksi
grid_m <- expand.grid(
Age = mean(data_multi$Age, na.rm = TRUE),
Depression_Score = seq(min(data_multi$Depression_Score),
max(data_multi$Depression_Score), length.out = 100),
Symptoms = mean(data_multi$Symptoms, na.rm = TRUE),
SocialMedia_Hours = mean(data_multi$SocialMedia_Hours, na.rm = TRUE),
Sleep_Hours = mean(data_multi$Sleep_Hours, na.rm = TRUE),
Nervous_Level = mean(data_multi$Nervous_Level, na.rm = TRUE),
Gender = factor("Perempuan", levels = levels(data_multi$Gender)),
Education_Level = factor("SMA", levels = levels(data_multi$Education_Level)),
Employment_Status = factor("Penuh Waktu", levels = levels(data_multi$Employment_Status)),
Low_Energy = factor("Kadang", levels = levels(data_multi$Low_Energy)),
Low_SelfEsteem = factor("Kadang", levels = levels(data_multi$Low_SelfEsteem)),
Search_Depression_Online = factor("Tidak", levels = levels(data_multi$Search_Depression_Online)),
Worsening_Depression = factor("Tidak", levels = levels(data_multi$Worsening_Depression)),
Self_Harm = factor("Tidak", levels = levels(data_multi$Self_Harm)),
Mental_Health_Support = factor("Tidak", levels = levels(data_multi$Mental_Health_Support))
)
grid_prob_m <- predict(fit_multi, newdata = grid_m, type = "probs")
prob_cols_m <- colnames(grid_prob_m)
grid_m %>%
bind_cols(as.data.frame(grid_prob_m)) %>%
pivot_longer(all_of(prob_cols_m), names_to = "Percobaan", values_to = "probabilitas") %>%
mutate(Percobaan = factor(Percobaan,
levels = c("Tidak Pernah","1 Kali","2 Kali","3 Kali"))) %>%
ggplot(aes(x = Depression_Score, y = probabilitas, color = Percobaan)) +
geom_line(linewidth = 1.35) +
scale_y_continuous(labels = percent_format()) +
scale_color_manual(values = pal) +
labs(
title = "Prediksi Probabilitas Jumlah Percobaan Bunuh Diri",
subtitle = "Variabel lain ditahan pada nilai rata-rata / modus; Gender = Perempuan",
x = "Skor Depresi (Depression_Score)",
y = "Probabilitas Prediksi",
color = "Jumlah Percobaan"
) +
theme_analisis()
💡 Interpretasi: Apabila garis 3 Kali meningkat seiring kenaikan skor depresi, model menunjukkan bahwa individu dengan tingkat depresi lebih tinggi memiliki probabilitas lebih besar untuk melakukan percobaan bunuh diri berulang. Garis Tidak Pernah yang menurun mengindikasikan pergeseran risiko dari tidak pernah mencoba ke percobaan berulang seiring peningkatan skor depresi.
2.8 Ringkasan Asumsi
tibble(
No = 1:6,
Asumsi = c("Respons Bersifat Nominal", "Observasi Independen",
"Tidak Ada Multikolinearitas Berat (VIF)", "Tidak Ada Perfect Separation",
"Linearitas pada Skala Logit (Box-Tidwell)",
"Independence of Irrelevant Alternatives (IIA)"),
Metode = c("Konseptual", "Konseptual / desain data", "VIF",
"Cek konvergensi & koefisien ekstrem",
"Interaksi X ln(X)", "Pendekatan konseptual"),
Keterangan = c("Terpenuhi — kategori tidak berurutan secara alami",
"Terpenuhi — setiap baris adalah individu independen",
"Lihat tabel VIF di atas",
"Lihat output konvergensi di atas",
"Lihat tabel Box-Tidwell di atas",
"Terpenuhi — tiap kategori merepresentasikan outcome klinis berbeda")
) %>%
kable(caption = "Ringkasan Pemeriksaan Asumsi — Multinomial") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| No | Asumsi | Metode | Keterangan |
|---|---|---|---|
| 1 | Respons Bersifat Nominal | Konseptual | Terpenuhi — kategori tidak berurutan secara alami |
| 2 | Observasi Independen | Konseptual / desain data | Terpenuhi — setiap baris adalah individu independen |
| 3 | Tidak Ada Multikolinearitas Berat (VIF) | VIF | Lihat tabel VIF di atas |
| 4 | Tidak Ada Perfect Separation | Cek konvergensi & koefisien ekstrem | Lihat output konvergensi di atas |
| 5 | Linearitas pada Skala Logit (Box-Tidwell) | Interaksi X ln(X) | Lihat tabel Box-Tidwell di atas |
| 6 | Independence of Irrelevant Alternatives (IIA) | Pendekatan konseptual | Terpenuhi — tiap kategori merepresentasikan outcome klinis berbeda |
3 Regresi Logistik Ordinal
3.1 Deskripsi Data
Dataset ini bersumber dari Mendeley Data (DOI: 10.17632/yv7ddpjw9y.1) oleh Khairunesa Isa (2020), Universiti Tun Hussein Onn Malaysia. Data mencakup 280 responden (mahasiswa diploma, sarjana, dan pascasarjana di universitas awam Malaysia) yang menilai kepuasan terhadap 14 fasilitas kampus menggunakan skala Likert 1–10. Skor rata-rata kemudian dikategorikan menjadi empat level kepuasan yang berurutan (Rendah < Sedang < Tinggi < Sangat Tinggi), menjadikan regresi logistik ordinal sebagai metode yang tepat.
raw_data <- read_excel("C:/Users/Salsabila Najwa/Downloads/Rawdata (1).xlsx", sheet = 1)
colnames(raw_data) <- c(
"timestamp","jantina","umur","etnik","agama","semester","level",
"universiti","bursary","health_centre","library","sport_centre",
"islamic_centre","atm","residential","transportation",
"wireless","col18","toilet","cafeteria","cleanliness","welfare"
)
data_ordinal <- raw_data %>%
select(-timestamp, -col18) %>%
mutate(
umur = case_when(
umur == "26 < 300 tahun" ~ "≥ 26 tahun",
TRUE ~ umur
),
universiti = str_to_upper(str_trim(universiti))
) %>%
mutate(across(bursary:welfare, as.numeric)) %>%
drop_na() %>%
mutate(
skor_rata = rowMeans(select(., bursary:welfare), na.rm = TRUE),
kepuasan = cut(
skor_rata,
breaks = c(-Inf, 4.5, 6.0, 7.5, Inf),
labels = c("Rendah","Sedang","Tinggi","Sangat Tinggi"),
ordered_result = TRUE
),
jantina = factor(jantina, levels = c("Lelaki","Perempuan")),
umur = factor(umur, levels = c("≤ 19 tahun","20 < 25 tahun","≥ 26 tahun")),
level = factor(level, levels = c("Diploma","Degree","Siswazah")),
universiti = fct_lump_min(factor(universiti), min = 10,
other_level = "Universiti Lain")
) %>%
drop_na(kepuasan, jantina, umur, level, universiti)
cat("Observasi final:", nrow(data_ordinal), "\n")## Observasi final: 280## Rows: 280
## Columns: 22
## $ jantina <fct> Lelaki, Perempuan, Perempuan, Perempuan, Perempuan, Lel…
## $ umur <fct> 20 < 25 tahun, 20 < 25 tahun, 20 < 25 tahun, 20 < 25 ta…
## $ etnik <chr> "Melayu", "Melayu", "Melayu", "Melayu", "Melayu", "Mela…
## $ agama <chr> "Islam", "Islam", "Islam", "Islam", "Islam", "Islam", "…
## $ semester <chr> "Semester 2 Tahun 2", "Semester 2 Tahun 1", "Semester 2…
## $ level <fct> Degree, Degree, Degree, Degree, Degree, Degree, Degree,…
## $ universiti <fct> UTHM, UTHM, UTHM, UTHM, UTHM, UTHM, UTHM, UTHM, Univers…
## $ bursary <dbl> 8, 3, 4, 9, 8, 7, 8, 7, 6, 2, 7, 8, 6, 8, 6, 9, 10, 4, …
## $ health_centre <dbl> 3, 1, 2, 9, 9, 8, 7, 8, 6, 2, 8, 8, 9, 3, 9, 9, 10, 10,…
## $ library <dbl> 8, 10, 8, 9, 9, 8, 10, 9, 6, 10, 10, 10, 9, 9, 10, 9, 1…
## $ sport_centre <dbl> 7, 5, 8, 9, 7, 8, 8, 7, 6, 6, 10, 8, 6, 7, 7, 6, 10, 1,…
## $ islamic_centre <dbl> 8, 6, 8, 10, 7, 10, 10, 9, 6, 7, 10, 10, 8, 8, 5, 10, 1…
## $ atm <dbl> 7, 3, 8, 9, 8, 6, 9, 8, 6, 7, 10, 8, 9, 7, 3, 5, 6, 10,…
## $ residential <dbl> 7, 3, 8, 10, 9, 7, 10, 8, 6, 7, 9, 10, 4, 8, 7, 9, 10, …
## $ transportation <dbl> 7, 5, 4, 9, 5, 6, 8, 8, 6, 4, 2, 10, 7, 7, 3, 4, 8, 10,…
## $ wireless <dbl> 6, 3, 4, 6, 4, 9, 7, 8, 6, 3, 3, 10, 4, 6, 4, 3, 9, 6, …
## $ toilet <dbl> 5, 3, 8, 7, 9, 4, 9, 8, 7, 5, 6, 3, 6, 8, 5, 6, 8, 10, …
## $ cafeteria <dbl> 7, 3, 6, 5, 7, 9, 9, 8, 6, 6, 2, 10, 8, 7, 4, 7, 7, 10,…
## $ cleanliness <dbl> 9, 3, 8, 6, 8, 6, 9, 8, 6, 8, 7, 9, 8, 8, 8, 7, 6, 10, …
## $ welfare <dbl> 7, 3, 8, 6, 8, 7, 9, 8, 5, 6, 8, 10, 7, 8, 6, 7, 6, 10,…
## $ skor_rata <dbl> 6.846154, 3.923077, 6.461538, 8.000000, 7.538462, 7.307…
## $ kepuasan <ord> Tinggi, Rendah, Tinggi, Sangat Tinggi, Sangat Tinggi, T…3.2 Statistik Deskriptif
# Distribusi variabel dependen
data_ordinal %>%
count(kepuasan) %>%
mutate(Proporsi = percent(n / sum(n), accuracy = 0.1)) %>%
rename(Kategori = kepuasan, Jumlah = n) %>%
kable(caption = "Distribusi tingkat kepuasan mahasiswa terhadap fasilitas kampus") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Kategori | Jumlah | Proporsi |
|---|---|---|
| Rendah | 16 | 5.7% |
| Sedang | 73 | 26.1% |
| Tinggi | 110 | 39.3% |
| Sangat Tinggi | 81 | 28.9% |
# Ringkasan skor per fasilitas
data_ordinal %>%
select(bursary:welfare) %>%
summarise(across(everything(),
list(mean = ~round(mean(.x, na.rm=TRUE), 2),
sd = ~round(sd(.x, na.rm=TRUE), 2)))) %>%
pivot_longer(everything(), names_to = c("fasilitas",".value"),
names_pattern = "^(.+)_(mean|sd)$") %>%
rename(Fasilitas = fasilitas, `Rata-rata` = mean, `Std. Deviasi` = sd) %>%
kable(caption = "Ringkasan skor kepuasan per fasilitas (skala 1–10)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Fasilitas | Rata-rata | Std. Deviasi |
|---|---|---|
| bursary | 6.51 | 2.01 |
| health_centre | 7.15 | 2.07 |
| library | 7.92 | 1.77 |
| sport_centre | 6.72 | 1.89 |
| islamic_centre | 7.46 | 2.02 |
| atm | 6.71 | 2.02 |
| residential | 6.93 | 2.22 |
| transportation | 5.78 | 2.31 |
| wireless | 5.46 | 2.56 |
| toilet | 6.18 | 1.91 |
| cafeteria | 6.74 | 2.07 |
| cleanliness | 6.68 | 2.04 |
| welfare | 6.74 | 1.95 |
3.3 Visualisasi Deskriptif
# Distribusi kepuasan
data_ordinal %>%
count(kepuasan) %>%
mutate(p = n / sum(n)) %>%
ggplot(aes(x = kepuasan, y = p, fill = kepuasan)) +
geom_col(width = 0.65, alpha = 0.94) +
geom_text(aes(label = percent(p, accuracy = 0.1)), vjust = -0.4, fontface = "bold") +
scale_y_continuous(labels = percent_format(), limits = c(0, 0.55)) +
scale_fill_manual(values = pal[1:4]) +
labs(title = "Distribusi Tingkat Kepuasan Mahasiswa",
subtitle = "Respons ordinal: kategori memiliki urutan alami",
x = NULL, y = "Proporsi") +
theme_analisis() + theme(legend.position = "none")
interpretasi: Interpretasi Distribusi Tingkat Kepuasan
Mahasiswa (Variabel Respons)
Berdasarkan bar chart distribusi kepuasan, variabel respons memiliki 4 kategori yang memiliki urutan alami (ordinal), yaitu dari “Rendah” hingga “Sangat Tinggi”. Sebaran datanya adalah sebagai berikut: Tingkat Kepuasan Tertinggi: Mayoritas mahasiswa berada pada kelompok kepuasan “Tinggi” (39.3%) dan “Sangat Tinggi” (28.9%). Jika diakumulasikan, terdapat sekitar 68.2% mahasiswa yang merasa puas pada level atas. * Tingkat Kepuasan Rendah: Hanya sebagian kecil mahasiswa yang merasa kepuasannya berada di tingkat “Rendah” (5.7%), sementara tingkat “Sedang” berada di angka 26.1%. Kesimpulan: Karakteristik data yang memiliki tingkatan logis ini memperkuat alasan penggunaan Regresi Logistik Ordinal dibandingkan logistik multinomial biasa, karena model ordinal akan memanfaatkan informasi urutan tersebut untuk menghasilkan estimasi yang lebih efisien.
3.4 Uji Asumsi
3.4.1 Asumsi 1 — Asumsi Proportional Odds (Uji Nominal)
Asumsi proportional odds menyatakan bahwa efek setiap prediktor
konstan di seluruh threshold kategori. Pengujian
dilakukan dengan ordinal::nominal_test().
data_model_o <- data_ordinal %>%
select(kepuasan, jantina, umur, level, universiti) %>%
drop_na()
fit_clm <- ordinal::clm(
kepuasan ~ jantina + umur + level + universiti,
data = data_model_o, link = "logit"
)
cat("== Uji Nominal (Proportional Odds) ==\n")## == Uji Nominal (Proportional Odds) ==## Tests of nominal effects
##
## formula: kepuasan ~ jantina + umur + level + universiti
## Df logLik AIC LRT Pr(>Chi)
## <none> -337.77 703.55
## jantina 2 -337.43 706.86 0.6895 0.7084
## umur 4 -334.95 705.91 5.6438 0.2274
## level 4 -334.83 705.67 5.8804 0.2083
## universiti 12 -332.41 716.82 10.7263 0.5525💡 Interpretasi: H₀: asumsi proportional odds terpenuhi. Jika p > 0.05 untuk semua prediktor, efek masing-masing prediktor dapat diasumsikan konstan di sepanjang threshold. Apabila ada prediktor dengan p < 0.05, model partial proportional odds atau model generalized perlu dipertimbangkan.
3.4.2 Asumsi 2 — Tidak Ada Multikolinearitas Berat (Cramér’s V)
cramers_v <- function(x, y) {
tbl <- table(x, y)
chi <- suppressWarnings(chisq.test(tbl)$statistic)
n <- sum(tbl); k <- min(nrow(tbl), ncol(tbl))
sqrt(chi / (n * (k - 1)))
}
pred_vars_o <- c("jantina","umur","level","universiti")
cv_matrix <- outer(pred_vars_o, pred_vars_o, Vectorize(function(a, b) {
if (a == b) return(1)
round(cramers_v(data_model_o[[a]], data_model_o[[b]]), 3)
}))
rownames(cv_matrix) <- pred_vars_o; colnames(cv_matrix) <- pred_vars_o
kable(cv_matrix,
caption = "Matriks Cramér's V antar prediktor (< 0.5 = tidak bermasalah)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| jantina | umur | level | universiti | |
|---|---|---|---|---|
| jantina | 1.000 | 0.088 | 0.079 | 0.228 |
| umur | 0.088 | 1.000 | 0.493 | 0.165 |
| level | 0.079 | 0.493 | 1.000 | 0.191 |
| universiti | 0.228 | 0.165 | 0.191 | 1.000 |
3.5 Estimasi Model
fit_ord <- MASS::polr(
kepuasan ~ jantina + umur + level + universiti,
data = data_model_o,
method = "logistic",
Hess = TRUE
)
summary(fit_ord)## Call:
## MASS::polr(formula = kepuasan ~ jantina + umur + level + universiti,
## data = data_model_o, Hess = TRUE, method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## jantinaPerempuan 0.4365 0.2424 1.8011
## umur20 < 25 tahun -1.6724 1.0465 -1.5981
## umur≥ 26 tahun -2.1960 1.2358 -1.7770
## levelDegree 2.0818 0.8333 2.4981
## levelSiswazah 1.5245 1.0118 1.5068
## universitiUMS 0.5699 0.4522 1.2603
## universitiUNIMAS 1.1186 0.6646 1.6831
## universitiUPNM 1.5541 0.7068 2.1989
## universitiUTEM 0.4571 0.5318 0.8595
## universitiUTHM 0.7708 0.4285 1.7990
## universitiUniversiti Lain 0.3822 0.4395 0.8697
##
## Intercepts:
## Value Std. Error t value
## Rendah|Sedang -1.7031 1.0326 -1.6494
## Sedang|Tinggi 0.4111 1.0155 0.4048
## Tinggi|Sangat Tinggi 2.1667 1.0243 2.1153
##
## Residual Deviance: 675.5492
## AIC: 703.54923.6 Tabel Koefisien, Odds Ratio, dan Interval Kepercayaan
ord_coef <- coef(summary(fit_ord))
result_ord <- as.data.frame(ord_coef) %>%
rownames_to_column("parameter") %>%
dplyr::rename(estimate = Value, std_error = `Std. Error`, t_value = `t value`) %>%
mutate(
p_value = 2 * (1 - pnorm(abs(t_value))),
jenis = ifelse(grepl("\\|", parameter), "Cutpoint (α)", "Koefisien (β)"),
OR = ifelse(jenis == "Koefisien (β)", round(exp(-estimate), 4), NA_real_),
CI_low = ifelse(jenis == "Koefisien (β)",
round(exp(-estimate - 1.96 * std_error), 4), NA_real_),
CI_high = ifelse(jenis == "Koefisien (β)",
round(exp(-estimate + 1.96 * std_error), 4), NA_real_)
) %>%
mutate(across(c(estimate, std_error, t_value, p_value), ~round(.x, 4)))
result_ord %>%
kable(
caption = "Ringkasan hasil regresi logistik ordinal (MASS::polr)",
col.names = c("Parameter", "Estimate (β)", "SE", "t-value", "p-value",
"Jenis", "OR = exp(−β)", "CI 2.5% OR", "CI 97.5% OR")
) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| Parameter | Estimate (β) | SE | t-value | p-value | Jenis | OR = exp(−β) | CI 2.5% OR | CI 97.5% OR |
|---|---|---|---|---|---|---|---|---|
| jantinaPerempuan | 0.4365 | 0.2424 | 1.8011 | 0.0717 | Koefisien (β) | 0.6463 | 0.4019 | 1.0393 |
| umur20 < 25 tahun | -1.6724 | 1.0465 | -1.5981 | 0.1100 | Koefisien (β) | 5.3250 | 0.6847 | 41.4098 |
| umur≥ 26 tahun | -2.1960 | 1.2358 | -1.7770 | 0.0756 | Koefisien (β) | 8.9889 | 0.7976 | 101.2999 |
| levelDegree | 2.0818 | 0.8333 | 2.4981 | 0.0125 | Koefisien (β) | 0.1247 | 0.0244 | 0.6386 |
| levelSiswazah | 1.5245 | 1.0118 | 1.5068 | 0.1319 | Koefisien (β) | 0.2177 | 0.0300 | 1.5818 |
| universitiUMS | 0.5699 | 0.4522 | 1.2603 | 0.2076 | Koefisien (β) | 0.5656 | 0.2331 | 1.3722 |
| universitiUNIMAS | 1.1186 | 0.6646 | 1.6831 | 0.0924 | Koefisien (β) | 0.3267 | 0.0888 | 1.2021 |
| universitiUPNM | 1.5541 | 0.7068 | 2.1989 | 0.0279 | Koefisien (β) | 0.2114 | 0.0529 | 0.8446 |
| universitiUTEM | 0.4571 | 0.5318 | 0.8595 | 0.3900 | Koefisien (β) | 0.6331 | 0.2233 | 1.7953 |
| universitiUTHM | 0.7708 | 0.4285 | 1.7990 | 0.0720 | Koefisien (β) | 0.4627 | 0.1998 | 1.0714 |
| universitiUniversiti Lain | 0.3822 | 0.4395 | 0.8697 | 0.3845 | Koefisien (β) | 0.6824 | 0.2884 | 1.6147 |
| Rendah|Sedang | -1.7031 | 1.0326 | -1.6494 | 0.0991 | Cutpoint (α) | NA | NA | NA |
| Sedang|Tinggi | 0.4111 | 1.0155 | 0.4048 | 0.6856 | Cutpoint (α) | NA | NA | NA |
| Tinggi|Sangat Tinggi | 2.1667 | 1.0243 | 2.1153 | 0.0344 | Cutpoint (α) | NA | NA | NA |
⚠️ Penting — Interpretasi Khusus
MASS::polr(Spesifikasi α − Xβ):Package
MASS::polrmenggunakan spesifikasi negatif pada linier prediktor: \[\text{logit}[P(Y \leq j \mid \mathbf{x})] = \alpha_j - \mathbf{x}^{\top}\boldsymbol{\beta}\] di mana \(\alpha_j\) adalah threshold (cutpoint) ke-\(j\) dan \(\boldsymbol{\beta}\) adalah vektor koefisien prediktor.Konsekuensinya, tanda koefisien (β) yang dikeluarkan
polrharus dibalik untuk mendapatkan OR yang intuitif:\[\text{OR untuk } P(Y > j) = \exp(-\hat{\beta})\]
- OR = exp(−β) > 1: Kenaikan prediktor meningkatkan odds berada di kategori yang lebih tinggi (semakin puas).
- OR = exp(−β) < 1: Kenaikan prediktor menurunkan odds berada di kategori yang lebih tinggi.
Berbeda dengan regresi logistik biner dan multinomial yang menggunakan OR = exp(+β).
Contoh: Jika
levelDegreememiliki β = 2.082, maka OR = exp(−2.082) = 0.125. Ini berarti mahasiswa Degree memiliki odds 0,125 kali untuk berada di kategori kepuasan lebih rendah dibanding Diploma — atau equivalennya, mahasiswa Degree lebih cenderung berada di kategori kepuasan lebih tinggi (OR “naik” = 1/0.125 = 8.0 kali lebih besar). Hasil ini signifikan (p = 0.012).
3.7 Uji Model (Likelihood Ratio Test dan Pseudo R²)
null_ord <- MASS::polr(kepuasan ~ 1, data = data_model_o,
method = "logistic", Hess = TRUE)
lr_stat <- 2 * (logLik(fit_ord) - logLik(null_ord))
lr_df <- length(coef(fit_ord))
lr_pval <- pchisq(as.numeric(lr_stat), df = lr_df, lower.tail = FALSE)
n_o <- nrow(data_model_o)
R2_cox_o <- 1 - exp((2/n_o) * (logLik(null_ord) - logLik(fit_ord)))
R2_nagel <- as.numeric(R2_cox_o) / (1 - exp((2/n_o) * as.numeric(logLik(null_ord))))
data.frame(
Keterangan = c("Log-Likelihood model penuh", "Log-Likelihood model null",
"Statistik LR (χ²)", "Derajat bebas", "p-value LR", "Nagelkerke R²"),
Nilai = c(round(as.numeric(logLik(fit_ord)), 3),
round(as.numeric(logLik(null_ord)), 3),
round(as.numeric(lr_stat), 3),
lr_df,
round(lr_pval, 4),
round(R2_nagel, 4))
) %>%
kable(caption = "Likelihood Ratio Test dan Pseudo R² — Model Ordinal") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Keterangan | Nilai |
|---|---|
| Log-Likelihood model penuh | -337.7750 |
| Log-Likelihood model null | -347.1730 |
| Statistik LR (χ²) | 18.7970 |
| Derajat bebas | 11.0000 |
| p-value LR | 0.0648 |
| Nagelkerke R² | 0.0709 |
3.8 Prediksi dan Confusion Matrix
pred_prob_o <- predict(fit_ord, newdata = data_model_o, type = "probs")
pred_class_o <- predict(fit_ord, newdata = data_model_o, type = "class")
conf_o <- table(Aktual = data_model_o$kepuasan, Prediksi = pred_class_o)
kable(addmargins(conf_o), caption = "Confusion Matrix — Model Ordinal") %>%
kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)| Rendah | Sedang | Tinggi | Sangat Tinggi | Sum | |
|---|---|---|---|---|---|
| Rendah | 0 | 2 | 14 | 0 | 16 |
| Sedang | 0 | 14 | 56 | 3 | 73 |
| Tinggi | 0 | 6 | 98 | 6 | 110 |
| Sangat Tinggi | 0 | 2 | 68 | 11 | 81 |
| Sum | 0 | 24 | 236 | 20 | 280 |
acc_o <- sum(diag(conf_o)) / sum(conf_o)
cat("Akurasi model ordinal:", round(acc_o * 100, 2), "%\n")## Akurasi model ordinal: 43.93 %3.9 Visualisasi
# Prediksi probabilitas per level pengajian
grid_level <- expand.grid(
jantina = factor("Perempuan", levels = levels(data_model_o$jantina)),
umur = factor("20 < 25 tahun", levels = levels(data_model_o$umur)),
level = factor(c("Diploma","Degree","Siswazah"), levels = levels(data_model_o$level)),
universiti = factor("UTHM", levels = levels(data_model_o$universiti))
)
predict(fit_ord, newdata = grid_level, type = "probs") %>%
as.data.frame() %>%
bind_cols(grid_level) %>%
pivot_longer(c("Rendah","Sedang","Tinggi","Sangat Tinggi"),
names_to = "kepuasan", values_to = "probabilitas") %>%
mutate(kepuasan = factor(kepuasan,
levels = c("Rendah","Sedang","Tinggi","Sangat Tinggi"))) %>%
ggplot(aes(x = level, y = probabilitas, fill = kepuasan)) +
geom_col(position = "stack", alpha = 0.9, width = 0.6) +
scale_y_continuous(labels = percent_format()) +
scale_fill_manual(values = pal[1:4]) +
labs(title = "Prediksi Probabilitas Kepuasan berdasarkan Level Pengajian",
subtitle = "Jantina = Perempuan, Umur = 20–25 tahun, Universiti = UTHM",
x = "Level Pengajian", y = "Probabilitas", fill = "Kepuasan") +
theme_analisis()
# Parallel lines — bukti visual proportional odds
grid_univ <- data.frame(
jantina = factor("Perempuan", levels = levels(data_model_o$jantina)),
umur = factor("20 < 25 tahun", levels = levels(data_model_o$umur)),
level = factor("Degree", levels = levels(data_model_o$level)),
universiti = factor(c("UTHM","UMS","UMP","UTEM","UNIMAS","UPNM","Universiti Lain"),
levels = levels(data_model_o$universiti))
)
eps <- 1e-6
predict(fit_ord, newdata = grid_univ, type = "probs") %>%
as.data.frame() %>%
bind_cols(grid_univ) %>%
mutate(
`P(Y ≤ Rendah)` = Rendah,
`P(Y ≤ Sedang)` = Rendah + Sedang,
`P(Y ≤ Tinggi)` = Rendah + Sedang + Tinggi
) %>%
pivot_longer(c(`P(Y ≤ Rendah)`,`P(Y ≤ Sedang)`,`P(Y ≤ Tinggi)`),
names_to = "batas", values_to = "prob") %>%
mutate(prob = pmin(pmax(prob, eps), 1 - eps),
logit_kumulatif = qlogis(prob)) %>%
ggplot(aes(x = universiti, y = logit_kumulatif, color = batas, group = batas)) +
geom_point(size = 3) + geom_line(linewidth = 1.1) +
scale_color_manual(values = pal[c(1,3,2)]) +
labs(title = "Parallel Lines — Bukti Asumsi Proportional Odds",
subtitle = "Garis sejajar pada cumulative logit menunjukkan asumsi terpenuhi",
x = "Universiti", y = "Logit kumulatif", color = "Batas kumulatif") +
theme_analisis() +
theme(axis.text.x = element_text(angle = 20, hjust = 1))
interpretasi: bagian 1. Pengujian Asumsi Proportional
Odds (Parallel Lines Plot)
Sebelum menarik kesimpulan dari model logistik ordinal, dilakukan
pengujian terhadap asumsi paling kritisnya, yaitu Asumsi
Proportional Odds (Asumsi Garis Paralel). Asumsi ini
mensyaratkan bahwa pengaruh dari variabel prediktor (dalam hal ini
Universiti) bernilai konstan pada setiap tingkatan logit
kumulatif.
Berdasarkan grafik Parallel Lines di atas, diperoleh
evaluasi sebagai berikut: Grafik menampilkan tiga garis batas kumulatif:
\(P(Y \le \text{Rendah})\), \(P(Y \le \text{Sedang})\), dan \(P(Y \le \text{Tinggi})\). Ketiga kurva
tersebut menunjukkan pola pergerakan yang sangat seirama dan
sejajar (paralel) di seluruh kategori universitas. Sebagai
contoh, saat nilai logit kumulatif menurun di UPNM, ketiga
garis turun bersamaan; dan ketika melonjak di UTEM,
ketiganya juga naik secara proporsional tanpa ada garis yang saling
berpotongan (crossing).
Kesimpulan: Karena visualisasi menunjukkan garis-garis yang paralel, maka Asumsi Proportional Odds terpenuhi. Hal ini menandakan model regresi logistik ordinal yang digunakan sudah valid dan reliabel untuk mengestimasi tingkat kepuasan.
bagian 2. Interpretasi Prediksi Probabilitas Kepuasan berdasarkan Level Pengajian
Grafik stacked bar chart di atas memvisualisasikan prediksi probabilitas tingkat kepuasan mahasiswa berdasarkan Level Pengajian (Diploma, Degree, dan Siswazah). Nilai prediksi ini dikontrol berdasarkan profil responden spesifik, yaitu: Jantina = Perempuan, Umur = 20–25 tahun, dan Universiti = UTHM.
Berikut adalah poin-poin analisis utama dari model ordinal tersebut:
Level Diploma: Mahasiswa Diploma diprediksi memiliki peluang kepuasan yang cenderung lebih rendah dibandingkan jenjang lainnya. Kategori “Sedang” sangat dominan (mendekati 50%), serta memiliki probabilitas merasakan kepuasan “Rendah” yang paling tinggi dibandingkan level Degree dan Siswazah. Level Degree (Sarjana): Pada jenjang ini, terjadi pergeseran probabilitas (probability shift) yang sangat positif ke arah kanan (tingkat kepuasan yang lebih tinggi). Peluang untuk mencapai tingkat kepuasan “Tinggi” dan “Sangat Tinggi” mendominasi secara mutlak (akumulasi di atas 75%), sementara probabilitas tingkat “Rendah” menyusut hingga hampir tidak terlihat. Level Siswazah (Pascasarjana): Kelompok mahasiswa pascasarjana menunjukkan pola yang kembali sedikit bergeser ke tengah. Tingkat kepuasan mereka didominasi oleh kategori “Tinggi” dan “Sedang”, namun peluang untuk merasa “Sangat Tinggi” masih terjaga cukup baik di kisaran 25%.
Kesimpulan: Model logistik ordinal berhasil menunjukkan bahwa Level Pengajian memberikan pengaruh nyata terhadap probabilitas tingkat kepuasan mahasiswa. Pada profil mahasiswi UTHM usia 20–25 tahun, kelompok responden dengan jenjang pendidikan Degree diprediksi memiliki tingkat kepuasan yang paling optimal dibandingkan jenjang Diploma maupun Siswazah.
4 Regresi Poisson
4.1 Deskripsi Data
Dataset ini diperoleh dari Mendeley Data (DOI: 10.17632/69xxkwghy4.1)
oleh Chinnathambi (2022), SRM University. Dataset memuat 1.338
pemegang polis asuransi kesehatan dengan variabel meliputi
usia, jenis kelamin, BMI, jumlah anak yang ditanggung, kelayakan diskon,
wilayah, dan biaya medis. Variabel respons adalah jumlah anak
yang ditanggung (children) — data hitung
(count data) bilangan bulat nonnegatif — menjadikan regresi
Poisson sebagai metode yang sesuai.
raw_ins <- read.csv("C:/Users/Salsabila Najwa/Downloads/medical_insurance.csv")
data_ins <- raw_ins %>%
select(-id, -premium) %>%
mutate(
sex = factor(sex, levels = c("male","female")),
region = factor(region, levels = c("northeast","northwest","southeast","southwest")),
discount_eligibility = factor(discount_eligibility, levels = c("no","yes"))
) %>%
drop_na()
cat("Observasi setelah cleaning:", nrow(data_ins), "\n")## Observasi setelah cleaning: 1338## Rows: 1,338
## Columns: 7
## $ age <int> 19, 18, 28, 33, 32, 31, 46, 37, 37, 60, 25, 62, 2…
## $ sex <fct> female, male, male, male, male, female, female, f…
## $ bmi <dbl> 27.9, 33.8, 33.0, 22.7, 28.9, 25.7, 33.4, 27.7, 2…
## $ children <int> 0, 1, 3, 0, 0, 0, 1, 3, 2, 0, 0, 0, 0, 0, 0, 1, 1…
## $ discount_eligibility <fct> yes, no, no, no, no, no, no, no, no, no, no, yes,…
## $ region <fct> southwest, southeast, southeast, northwest, north…
## $ expenses <dbl> 16884.92, 1725.55, 4449.46, 21984.47, 3866.86, 37…4.2 Statistik Deskriptif
# Distribusi variabel dependen
data_ins %>%
count(children) %>%
mutate(Proporsi = percent(n / sum(n), accuracy = 0.1)) %>%
rename(`Jumlah Anak` = children, Frekuensi = n) %>%
kable(caption = "Distribusi jumlah anak yang ditanggung asuransi") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Jumlah Anak | Frekuensi | Proporsi |
|---|---|---|
| 0 | 574 | 42.9% |
| 1 | 324 | 24.2% |
| 2 | 240 | 17.9% |
| 3 | 157 | 11.7% |
| 4 | 25 | 1.9% |
| 5 | 18 | 1.3% |
# Statistik numerik
data_ins %>%
select(age, bmi, children, expenses) %>%
summarise(across(everything(),
list(mean = ~round(mean(.x), 2), sd = ~round(sd(.x), 2),
min = ~min(.x), max = ~max(.x)))) %>%
pivot_longer(everything(), names_to = c("variabel",".value"),
names_pattern = "^(.+)_(mean|sd|min|max)$") %>%
rename(Variabel = variabel, `Rata-rata` = mean, `Std. Deviasi` = sd,
Min = min, Max = max) %>%
kable(caption = "Statistik deskriptif variabel numerik") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Variabel | Rata-rata | Std. Deviasi | Min | Max |
|---|---|---|---|---|
| age | 39.21 | 14.05 | 18.00 | 64.00 |
| bmi | 30.67 | 6.10 | 16.00 | 53.10 |
| children | 1.09 | 1.21 | 0.00 | 5.00 |
| expenses | 13270.42 | 12110.01 | 1121.87 | 63770.43 |
4.3 Visualisasi Awal
data_ins %>%
count(children) %>%
ggplot(aes(x = children, y = n)) +
geom_col(width = 0.75, fill = pal[1], alpha = 0.92) +
geom_text(aes(label = n), vjust = -0.4, fontface = "bold") +
scale_x_continuous(breaks = 0:5) +
labs(title = "Distribusi Jumlah Anak yang Ditanggung Asuransi",
subtitle = "Respons Poisson: variabel hitung bilangan bulat nonnegatif",
x = "Jumlah anak", y = "Frekuensi") +
theme_analisis()
ggplot(data_ins, aes(x = age, y = children, color = sex)) +
geom_jitter(alpha = 0.4, height = 0.2, width = 0.3) +
geom_smooth(method = "loess", se = FALSE, linewidth = 1.2) +
scale_color_manual(values = pal[c(2,4)]) +
labs(title = "Hubungan Usia dan Jumlah Anak berdasarkan Jenis Kelamin",
x = "Usia", y = "Jumlah anak", color = "Jenis kelamin") +
theme_analisis()
interpretasi: interpretasi Distribusi jumlah anak yang
ditanggung asuransi menunjukkan bahwa sebagian besar responden tidak
memiliki anak yang ditanggung (574 responden).
Frekuensi kemudian menurun seiring bertambahnya jumlah anak, yaitu
324 responden memiliki 1 anak, 240
responden memiliki 2 anak, dan 157 responden
memiliki 3 anak. Jumlah responden dengan 4–5 anak relatif sedikit. Pola
ini menunjukkan distribusi yang condong ke kanan (miring ke
kanan), sehingga sebagian besar responden memiliki jumlah anak yang
ditanggung relatif rendah. Interpretasi Grafik menunjukkan bahwa
hubungan antara usia dan jumlah anak membentuk pola nonlinier. Jumlah
anak cenderung meningkat seiring pertambahan usia hingga sekitar usia
40–45 tahun, kemudian menurun pada kelompok usia yang lebih tua. Pola
ini terlihat baik pada laki-laki maupun perempuan. Selain itu, tidak
terdapat perbedaan yang mencolok antara laki-laki dan perempuan dalam
jumlah anak yang dimiliki, karena kedua kurva tren menunjukkan pola yang
hampir sama. Sebagian besar responden memiliki 0–2 anak pada berbagai
kelompok usia. ## Pemodelan
4.3.1 Model Penuh
fit_pois <- glm(
children ~ age + sex + bmi + discount_eligibility + region + expenses,
data = data_ins,
family = poisson(link = "log")
)
summary(fit_pois)##
## Call:
## glm(formula = children ~ age + sex + bmi + discount_eligibility +
## region + expenses, family = poisson(link = "log"), data = data_ins)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.690e-02 1.643e-01 0.285 0.7753
## age -9.278e-04 2.144e-03 -0.433 0.6652
## sexfemale -4.022e-02 5.253e-02 -0.766 0.4439
## bmi -3.122e-03 4.751e-03 -0.657 0.5110
## discount_eligibilityyes -3.645e-01 1.201e-01 -3.036 0.0024 **
## regionnorthwest 9.844e-02 7.507e-02 1.311 0.1897
## regionsoutheast 8.615e-03 7.711e-02 0.112 0.9110
## regionsouthwest 9.959e-02 7.552e-02 1.319 0.1873
## expenses 1.603e-05 4.057e-06 3.950 7.8e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2001.6 on 1337 degrees of freedom
## Residual deviance: 1979.5 on 1329 degrees of freedom
## AIC: 3886.8
##
## Number of Fisher Scoring iterations: 5# Tabel IRR model penuh
pois_coef <- as.data.frame(coef(summary(fit_pois))) %>%
rownames_to_column("parameter") %>%
dplyr::rename(estimate = Estimate, std_error = `Std. Error`,
z_value = `z value`, p_value = `Pr(>|z|)`) %>%
mutate(
IRR = round(exp(estimate), 4),
CI_low = round(exp(estimate - 1.96 * std_error), 4),
CI_high = round(exp(estimate + 1.96 * std_error), 4),
perubahan_persen = round(100 * (exp(estimate) - 1), 3)
) %>%
mutate(across(c(estimate, std_error, z_value, p_value), ~round(.x, 4)))
pois_coef %>%
kable(
caption = "Ringkasan hasil regresi Poisson — Model Penuh",
col.names = c("Parameter","Estimate","SE","z-value","p-value",
"IRR","CI 2.5%","CI 97.5%","Perubahan (%)")
) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| Parameter | Estimate | SE | z-value | p-value | IRR | CI 2.5% | CI 97.5% | Perubahan (%) |
|---|---|---|---|---|---|---|---|---|
| (Intercept) | 0.0469 | 0.1643 | 0.2855 | 0.7753 | 1.0480 | 0.7595 | 1.4461 | 4.801 |
| age | -0.0009 | 0.0021 | -0.4327 | 0.6652 | 0.9991 | 0.9949 | 1.0033 | -0.093 |
| sexfemale | -0.0402 | 0.0525 | -0.7657 | 0.4439 | 0.9606 | 0.8666 | 1.0647 | -3.942 |
| bmi | -0.0031 | 0.0048 | -0.6572 | 0.5110 | 0.9969 | 0.9876 | 1.0062 | -0.312 |
| discount_eligibilityyes | -0.3645 | 0.1201 | -3.0357 | 0.0024 | 0.6946 | 0.5489 | 0.8788 | -30.543 |
| regionnorthwest | 0.0984 | 0.0751 | 1.3113 | 0.1897 | 1.1034 | 0.9525 | 1.2783 | 10.345 |
| regionsoutheast | 0.0086 | 0.0771 | 0.1117 | 0.9110 | 1.0087 | 0.8672 | 1.1732 | 0.865 |
| regionsouthwest | 0.0996 | 0.0755 | 1.3188 | 0.1873 | 1.1047 | 0.9527 | 1.2810 | 10.472 |
| expenses | 0.0000 | 0.0000 | 3.9505 | 0.0001 | 1.0000 | 1.0000 | 1.0000 | 0.002 |
4.3.2 Model Reduksi (Model Final)
Berdasarkan model penuh, hanya discount_eligibility (p =
0.002) dan expenses (p < 0.001) yang signifikan.
Prediktor lain (age, sex, bmi, region) tidak signifikan dan
dikeluarkan.
fit_pois_red <- glm(
children ~ discount_eligibility + expenses,
data = data_ins,
family = poisson(link = "log")
)
summary(fit_pois_red)##
## Call:
## glm(formula = children ~ discount_eligibility + expenses, family = poisson(link = "log"),
## data = data_ins)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.706e-02 4.213e-02 -0.880 0.3790
## discount_eligibilityyes -3.239e-01 1.058e-01 -3.061 0.0022 **
## expenses 1.419e-05 3.365e-06 4.217 2.48e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2001.6 on 1337 degrees of freedom
## Residual deviance: 1984.1 on 1335 degrees of freedom
## AIC: 3879.4
##
## Number of Fisher Scoring iterations: 5# Perbandingan dua model
data.frame(
Keterangan = c("Jumlah prediktor","Null deviance",
"Residual deviance","df residual","AIC"),
`Model Penuh` = c(8, round(fit_pois$null.deviance, 3),
round(fit_pois$deviance, 3),
fit_pois$df.residual, round(AIC(fit_pois), 3)),
`Model Reduksi` = c(2, round(fit_pois_red$null.deviance, 3),
round(fit_pois_red$deviance, 3),
fit_pois_red$df.residual, round(AIC(fit_pois_red), 3)),
check.names = FALSE
) %>%
kable(caption = "Perbandingan Model Penuh vs Reduksi") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Keterangan | Model Penuh | Model Reduksi |
|---|---|---|
| Jumlah prediktor | 8.000 | 2.000 |
| Null deviance | 2001.581 | 2001.581 |
| Residual deviance | 1979.481 | 1984.111 |
| df residual | 1329.000 | 1335.000 |
| AIC | 3886.785 | 3879.414 |
4.4 Uji Asumsi
4.4.1 Asumsi 1 — Overdispersi (Pearson Dispersion)
disp_pearson <- sum(residuals(fit_pois_red, type = "pearson")^2) /
df.residual(fit_pois_red)
data.frame(
`Dispersion Pearson` = round(disp_pearson, 4),
Interpretasi = case_when(
disp_pearson < 1.5 ~ "Tidak ada indikasi overdispersion berat — equidispersion terpenuhi",
disp_pearson < 2.5 ~ "Overdispersion sedang — pertimbangkan Negative Binomial",
TRUE ~ "Overdispersion kuat — Negative Binomial lebih sesuai"
),
check.names = FALSE
) %>%
kable(caption = "Uji Overdispersi: Pearson Dispersion Statistic — Model Reduksi") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| Dispersion Pearson | Interpretasi |
|---|---|
| 1.3246 | Tidak ada indikasi overdispersion berat — equidispersion terpenuhi |
##
## == Formal Overdispersion Test (AER) ==##
## Overdispersion test
##
## data: fit_pois_red
## z = 7.5083, p-value = 2.995e-14
## alternative hypothesis: true alpha is greater than 0
## sample estimates:
## alpha
## 0.3222182💡 Interpretasi: Pearson dispersion = 1.32 (< 1.5) → tidak ada indikasi overdispersi berat. Meskipun
dispersiontest()formal signifikan (p < 0.05), hal ini wajar pada sampel besar (n = 1338) karena uji sangat sensitif terhadap ukuran sampel. Model Poisson masih valid sebagai pendekatan dalam konteks ini.
4.4.2 Asumsi 2 — Proporsi Nilai Nol
zero_share <- mean(data_ins$children == 0)
data.frame(
`Proporsi Y = 0` = percent(zero_share, accuracy = 0.1),
Catatan = ifelse(zero_share > 0.5,
"Proporsi nol > 50% — pertimbangkan Zero-Inflated Poisson",
"Proporsi nol wajar (< 50%) — model Poisson standar sesuai"),
check.names = FALSE
) %>%
kable(caption = "Pemeriksaan proporsi nilai nol") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Proporsi Y = 0 | Catatan |
|---|---|
| 42.9% | Proporsi nol wajar (< 50%) — model Poisson standar sesuai |
4.4.3 Asumsi 3 — Goodness-of-Fit (Deviance Test)
dev_pval_red <- pchisq(deviance(fit_pois_red), df.residual(fit_pois_red),
lower.tail = FALSE)
dev_pval_ful <- pchisq(deviance(fit_pois), df.residual(fit_pois),
lower.tail = FALSE)
data.frame(
Keterangan = c("Residual deviance","Derajat bebas","p-value"),
`Model Penuh` = c(round(fit_pois$deviance, 3), fit_pois$df.residual,
round(dev_pval_ful, 4)),
`Model Reduksi` = c(round(fit_pois_red$deviance, 3), fit_pois_red$df.residual,
round(dev_pval_red, 4)),
check.names = FALSE
) %>%
kable(caption = "Uji Goodness-of-Fit: Deviance Test") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Keterangan | Model Penuh | Model Reduksi |
|---|---|---|
| Residual deviance | 1979.481 | 1984.111 |
| Derajat bebas | 1329.000 | 1335.000 |
| p-value | 0.000 | 0.000 |
💡 Interpretasi: Kedua model menghasilkan p ≈ 0.000. Pada n = 1338, deviance test sangat sensitif sehingga bahkan model yang cukup baik dapat menghasilkan p < 0.05. Evaluasi model sebaiknya juga mempertimbangkan AIC, IRR, dan logika substansi.
4.4.4 Asumsi 4 — Multikolinearitas (VIF)
# Helper: ekstrak nilai VIF yang bisa dibandingkan dari output car::vif()
# Untuk variabel kategorik multi-level, car::vif() mengembalikan GVIF matrix
extract_vif <- function(vif_out) {
if (is.matrix(vif_out)) {
# GVIF^(1/(2*Df)) dikuadratkan agar setara dengan VIF skalar
adj <- vif_out[, "GVIF^(1/(2*Df))"]^2
data.frame(Variabel = rownames(vif_out),
GVIF = round(vif_out[, "GVIF"], 3),
Df = as.integer(vif_out[, "Df"]),
`VIF_adj` = round(adj, 3),
Keterangan = case_when(adj > 10 ~ "Berat", adj > 5 ~ "Sedang", TRUE ~ "Aman"),
check.names = FALSE)
} else {
data.frame(Variabel = names(vif_out),
VIF = round(as.numeric(vif_out), 3),
Keterangan = case_when(
as.numeric(vif_out) > 10 ~ "Berat",
as.numeric(vif_out) > 5 ~ "Sedang",
TRUE ~ "Aman"))
}
}
cat("== VIF Model Penuh ==\n")## == VIF Model Penuh ==extract_vif(car::vif(fit_pois)) %>%
kable(caption = "VIF Model Penuh Poisson") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Variabel | GVIF | Df | VIF_adj | Keterangan | |
|---|---|---|---|---|---|
| age | age | 1.338 | 1 | 1.338 | Aman |
| sex | sex | 1.010 | 1 | 1.010 | Aman |
| bmi | bmi | 1.212 | 1 | 1.212 | Aman |
| discount_eligibility | discount_eligibility | 3.481 | 1 | 3.481 | Aman |
| region | region | 1.100 | 3 | 1.032 | Aman |
| expenses | expenses | 3.942 | 1 | 3.942 | Aman |
##
## == VIF Model Reduksi ==extract_vif(car::vif(fit_pois_red)) %>%
kable(caption = "VIF Model Reduksi Poisson") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Variabel | VIF | Keterangan |
|---|---|---|
| discount_eligibility | 2.704 | Aman |
| expenses | 2.704 | Aman |
💡 Interpretasi: Semua VIF < 5 (model penuh) dan VIF = 2.704 (model reduksi, di bawah threshold 5). Tidak ada masalah multikolinearitas.
4.4.5 Asumsi 5 — Likelihood Ratio Test
null_pois <- glm(children ~ 1, data = data_ins, family = poisson(link = "log"))
cat("-- LRT: Null vs Model Penuh --\n")## -- LRT: Null vs Model Penuh --##
## -- LRT: Null vs Model Reduksi --##
## -- LRT: Model Reduksi vs Model Penuh --💡 Interpretasi: - Null vs Model Reduksi: p < 0.001 → model reduksi signifikan lebih baik dari model null. - Model Reduksi vs Penuh: p = 0.592 (> 0.05) → model reduksi tidak berbeda signifikan dari model penuh. Model reduksi lebih dipilih karena lebih parsimoni (prinsip Occam’s razor) dengan AIC lebih rendah (3879.4 < 3886.8).
4.5 IRR Model Reduksi (Model Final)
irr_red <- as.data.frame(coef(summary(fit_pois_red))) %>%
rownames_to_column("parameter") %>%
dplyr::rename(estimate = Estimate, std_error = `Std. Error`,
z_value = `z value`, p_value = `Pr(>|z|)`) %>%
mutate(
IRR = round(exp(estimate), 4),
CI_low = round(exp(estimate - 1.96 * std_error), 4),
CI_high = round(exp(estimate + 1.96 * std_error), 4),
perubahan_persen = round(100 * (exp(estimate) - 1), 3)
) %>%
mutate(across(c(estimate, std_error, z_value, p_value), ~round(.x, 4)))
irr_red %>%
kable(
caption = "IRR dan Interval Kepercayaan 95% — Model Reduksi (Model Final)",
col.names = c("Parameter","Estimate","SE","z-value","p-value",
"IRR","CI 2.5%","CI 97.5%","Perubahan (%)")
) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| Parameter | Estimate | SE | z-value | p-value | IRR | CI 2.5% | CI 97.5% | Perubahan (%) |
|---|---|---|---|---|---|---|---|---|
| (Intercept) | -0.0371 | 0.0421 | -0.8798 | 0.3790 | 0.9636 | 0.8872 | 1.0466 | -3.638 |
| discount_eligibilityyes | -0.3239 | 0.1058 | -3.0612 | 0.0022 | 0.7233 | 0.5878 | 0.8900 | -27.669 |
| expenses | 0.0000 | 0.0000 | 4.2170 | 0.0000 | 1.0000 | 1.0000 | 1.0000 | 0.001 |
Panduan Interpretasi IRR — Regresi Poisson:
Model regresi Poisson menggunakan fungsi log sebagai link: \[\log(\mu) = \beta_0 + \beta_1 X_1 + \cdots + \beta_k X_k\] di mana \(\mu = E(Y) = \lambda\) adalah rata-rata jumlah kejadian yang diharapkan.
IRR (Incidence Rate Ratio) = exp(β): - IRR > 1: Kenaikan satu satuan prediktor meningkatkan rata-rata jumlah anak yang diprediksi. - IRR < 1: Kenaikan satu satuan prediktor menurunkan rata-rata jumlah anak yang diprediksi. - Signifikan jika p-value < 0.05 atau CI 95% tidak mencakup angka 1.
Interpretasi prediktor signifikan:
discount_eligibility = yes(IRR = 0.723, CI 95%: 0.587–0.890, p = 0.002): Pemegang asuransi yang memenuhi syarat diskon memiliki prediksi rata-rata jumlah anak 0,723 kali dibandingkan yang tidak memenuhi syarat diskon, dengan biaya medis konstan. Artinya, kelayakan diskon dikaitkan dengan jumlah anak yang lebih sedikit (~27,7% lebih rendah). Secara kontekstual, diskon kemungkinan diberikan kepada pemegang polis dengan tanggungan lebih ringan sehingga risiko klaim lebih rendah.
expenses(IRR = 1.0000 per 1 unit, p < 0.001): Pengaruh per 1 unit expenses sangat kecil (~0.001%), namun pada skala praktis:
- Kenaikan expenses 10.000 unit → IRR kumulatif ≈ 1.153 → jumlah anak naik ~15.3%.
- Kenaikan expenses 50.000 unit → IRR kumulatif ≈ 2.034 → jumlah anak naik ~103%. Temuan ini konsisten: keluarga dengan lebih banyak anak cenderung memiliki total biaya medis lebih besar.
4.6 Visualisasi Model Reduksi
# Prediksi vs expenses per kelompok diskon
grid_exp <- expand.grid(
discount_eligibility = factor(c("no","yes"), levels = c("no","yes")),
expenses = seq(min(data_ins$expenses), max(data_ins$expenses),
length.out = 200)
)
pred_exp <- predict(fit_pois_red, newdata = grid_exp, type = "link", se.fit = TRUE)
grid_exp %>%
mutate(
fit = exp(pred_exp$fit),
fit_low = exp(pred_exp$fit - 1.96 * pred_exp$se.fit),
fit_high = exp(pred_exp$fit + 1.96 * pred_exp$se.fit)
) %>%
ggplot(aes(x = expenses, y = fit,
color = discount_eligibility, fill = discount_eligibility)) +
geom_ribbon(aes(ymin = fit_low, ymax = fit_high), alpha = 0.15, color = NA) +
geom_line(linewidth = 1.3) +
scale_color_manual(values = pal[c(2,4)],
labels = c("Tidak eligible diskon","Eligible diskon")) +
scale_fill_manual(values = pal[c(2,4)],
labels = c("Tidak eligible diskon","Eligible diskon")) +
labs(title = "Prediksi Jumlah Anak berdasarkan Expenses dan Kelayakan Diskon",
subtitle = "Model Reduksi Poisson: children ~ discount_eligibility + expenses",
x = "Biaya medis (expenses)", y = "Prediksi rata-rata jumlah anak",
color = "Kelayakan diskon", fill = "Kelayakan diskon") +
theme_analisis()
# Bar chart IRR
irr_red %>%
filter(parameter != "(Intercept)") %>%
mutate(
signif = ifelse(p_value < 0.05, "Signifikan (p<0.05)", "Tidak signifikan"),
parameter = case_when(
parameter == "discount_eligibilityyes" ~ "Diskon: Yes vs No",
parameter == "expenses" ~ "Expenses",
TRUE ~ parameter
)
) %>%
ggplot(aes(x = reorder(parameter, IRR), y = IRR, fill = signif)) +
geom_col(alpha = 0.85, width = 0.5) +
geom_errorbar(aes(ymin = CI_low, ymax = CI_high), width = 0.15, linewidth = 0.8) +
geom_hline(yintercept = 1, linetype = "dashed", color = "#6c757d", linewidth = 0.8) +
scale_fill_manual(values = c("Signifikan (p<0.05)" = pal[1], "Tidak signifikan" = "#adb5bd")) +
coord_flip() +
labs(title = "Incidence Rate Ratio (IRR) — Model Reduksi",
x = "Prediktor",
y = "IRR (garis putus-putus = IRR 1.0 / tidak ada efek)",
fill = "Signifikansi") +
theme_analisis()
**Interpretasi*: 1. Hasil Prediksi Model Poisson
Berdasarkan grafik hasil prediksi model regresi Poisson, variabel expenses menunjukkan hubungan positif terhadap jumlah anak (children). Semakin besar biaya medis yang dikeluarkan, semakin besar pula nilai ekspektasi jumlah anak yang diprediksi oleh model.
Variabel discount eligibility juga menunjukkan perbedaan prediksi jumlah anak. Individu yang memenuhi syarat diskon (eligible discount) cenderung memiliki nilai ekspektasi jumlah anak yang lebih rendah dibandingkan individu yang tidak memenuhi syarat diskon (not eligible discount) pada tingkat biaya medis yang sama.
Dengan demikian, model mengindikasikan bahwa biaya medis berasosiasi positif dengan jumlah anak, sementara status kelayakan diskon berasosiasi dengan perbedaan tingkat jumlah anak yang diprediksi.
2.Interpretasi Incidence Rate Ratio (IRR)
Gambar di atas menampilkan nilai Incidence Rate Ratio (IRR) dari model regresi Poisson yang telah direduksi. Garis putus-putus pada IRR = 1 menunjukkan tidak adanya pengaruh variabel prediktor terhadap jumlah anak.
Variabel expenses memiliki nilai IRR yang sangat mendekati 1. Hal ini menunjukkan bahwa setiap peningkatan satu satuan biaya medis hanya memberikan perubahan yang sangat kecil terhadap ekspektasi jumlah anak. Meskipun demikian, berdasarkan hasil pengujian model, variabel ini tetap berkontribusi dalam menjelaskan variasi jumlah anak.
Sementara itu, variabel discount eligibility memiliki nilai IRR sekitar 0,72 dengan interval kepercayaan yang seluruhnya berada di bawah nilai 1. Hasil ini menunjukkan bahwa individu yang memenuhi syarat diskon (eligible discount) memiliki ekspektasi jumlah anak yang lebih rendah dibandingkan individu yang tidak memenuhi syarat diskon (not eligible discount).
Secara kuantitatif, nilai IRR sebesar 0,72 mengindikasikan bahwa ekspektasi jumlah anak pada kelompok yang memenuhi syarat diskon adalah sekitar 0,72 kali dibandingkan kelompok yang tidak memenuhi syarat diskon, atau setara dengan penurunan sekitar 28% \((1 - 0,72) \times 100\%\).
Karena interval kepercayaan tidak memotong nilai 1 dan hasil pengujian menunjukkan nilai p-value kurang dari 0,05, maka pengaruh variabel discount eligibility terhadap jumlah anak dapat dikatakan signifikan secara statistik.
4.7 Ringkasan Perbandingan Model
data.frame(
Kriteria = c("Jumlah prediktor","AIC","Residual deviance",
"Prediktor signifikan","VIF multikolinearitas",
"LRT vs null","LRT reduksi vs penuh"),
`Model Penuh` = c("8", round(AIC(fit_pois),1), round(fit_pois$deviance,1),
"2 dari 8", "Semua GVIF^(1/(2Df)) < 2 (aman)",
"p < 0.001 (signifikan)", "—"),
`Model Reduksi` = c("2", round(AIC(fit_pois_red),1),
round(fit_pois_red$deviance,1), "2 dari 2",
"2.704 < 5 (aman)", "p < 0.001 (signifikan)",
"p = 0.592 (tidak berbeda)"),
check.names = FALSE
) %>%
kable(caption = "Perbandingan komprehensif Model Penuh vs Model Reduksi") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| Kriteria | Model Penuh | Model Reduksi |
|---|---|---|
| Jumlah prediktor | 8 | 2 |
| AIC | 3886.8 | 3879.4 |
| Residual deviance | 1979.5 | 1984.1 |
| Prediktor signifikan | 2 dari 8 | 2 dari 2 |
| VIF multikolinearitas | Semua GVIF^(1/(2Df)) < 2 (aman) | 2.704 < 5 (aman) |
| LRT vs null | p < 0.001 (signifikan) | p < 0.001 (signifikan) |
| LRT reduksi vs penuh | — | p = 0.592 (tidak berbeda) |
5 Ringkasan Komparatif Keempat Model
tibble(
Model = c("Logistik Biner","Logistik Multinomial",
"Logistik Ordinal","Poisson"),
Dataset = c("Caesarean Delivery (Ghana Health Service)",
"Depression & Mental Health Classification",
"Student Satisfaction on Facilities (Malaysia)",
"Medical Insurance Dataset"),
`Variabel Respons` = c("Jenis persalinan (0/1)",
"Jumlah percobaan bunuh diri (0–3)",
"Tingkat kepuasan (Rendah–Sangat Tinggi)",
"Jumlah anak yang ditanggung (0–5)"),
`Skala Respons` = c("Biner (nominal 2 kategori)",
"Nominal (4 kategori tidak berurutan)",
"Ordinal (4 kategori berurutan)",
"Count (bilangan bulat nonnegatif)"),
`Ukuran Efek` = c("Odds Ratio (OR) = exp(β)",
"Relative Risk Ratio (RRR) = exp(β)",
"OR = exp(−β) karena MASS::polr: α_j − Xβ",
"Incidence Rate Ratio (IRR) = exp(β)"),
`Package R` = c("stats::glm(family=binomial)",
"nnet::multinom()",
"MASS::polr() / ordinal::clm()",
"stats::glm(family=poisson)")
) %>%
kable(caption = "Ringkasan komparatif keempat model regresi logistik") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE) %>%
column_spec(5, width = "22%")| Model | Dataset | Variabel Respons | Skala Respons | Ukuran Efek | Package R |
|---|---|---|---|---|---|
| Logistik Biner | Caesarean Delivery (Ghana Health Service) | Jenis persalinan (0/1) | Biner (nominal 2 kategori) | Odds Ratio (OR) = exp(β) | stats::glm(family=binomial) |
| Logistik Multinomial | Depression & Mental Health Classification | Jumlah percobaan bunuh diri (0–3) | Nominal (4 kategori tidak berurutan) | Relative Risk Ratio (RRR) = exp(β) | nnet::multinom() |
| Logistik Ordinal | Student Satisfaction on Facilities (Malaysia) | Tingkat kepuasan (Rendah–Sangat Tinggi) | Ordinal (4 kategori berurutan) | OR = exp(−β) karena MASS::polr: α_j − Xβ | MASS::polr() / ordinal::clm() |
| Poisson | Medical Insurance Dataset | Jumlah anak yang ditanggung (0–5) | Count (bilangan bulat nonnegatif) | Incidence Rate Ratio (IRR) = exp(β) | stats::glm(family=poisson) |
Catatan Kritis — Perbedaan Tanda pada
MASS::polr:
Model Spesifikasi Link OR / IRR Arah Interpretasi Biner (glm) logit(P) = Xβ exp(+β) β > 0 → OR > 1 → event lebih mungkin Multinomial (nnet) log(P_k/P_ref) = Xβ exp(+β) β > 0 → RRR > 1 → lebih mungkin ke kategori k Ordinal (MASS::polr) logit[P(Y≤j)] = α_j − Xβ exp(−β) β > 0 → lebih mungkin ke kategori TINGGI Poisson (glm) log(λ) = Xβ exp(+β) β > 0 → IRR > 1 → rata-rata count naik Karena
MASS::polrmenggunakan tanda negatif pada prediktor, koefisien positif (β > 0) justru berarti prediktor meningkatkan odds berada di kategori yang lebih tinggi — berlawanan dengan intuisi biasa dari regresi logistik.