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

cat("Missing values per kolom:\n")

## Missing values per kolom:

colSums(is.na(caesar))

##     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")

Distribusi variabel dependen
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 : 41

1.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))
}

tab_bivariat("religion", "Agama")

Proporsi Caesar menurut Agama
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.

tab_bivariat("age", "Kelompok Usia")

Proporsi Caesar menurut Kelompok Usia
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.

tab_bivariat("location", "Lokasi")

Proporsi Caesar menurut Lokasi
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.

tab_bivariat("period_birth", "Periode Kelahiran")

Proporsi Caesar menurut Periode Kelahiran
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%.

tab_bivariat("weight_birth", "Berat Bayi Lahir")

Proporsi Caesar menurut Berat Bayi Lahir
Kategori	n_Caesar	Proporsi_Caesar
Rendah (<2.5kg)	41	0.436
Normal (2.5-3.9kg)	67	0.372
Tinggi (>=4kg)	39	0.333

tab_bivariat("facility", "Tipe Fasilitas")

Proporsi Caesar menurut Tipe Fasilitas
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

tab_bivariat("marital", "Status Pernikahan")

Proporsi Caesar menurut Status Pernikahan
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

tab_bivariat("ethnic", "Kelompok Etnis")

Proporsi Caesar menurut Kelompok Etnis
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

tab_bivariat("occ_mother", "Pekerjaan Ibu")

Proporsi Caesar menurut Pekerjaan Ibu
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

tab_bivariat("occ_father", "Pekerjaan Ayah")

Proporsi Caesar menurut Pekerjaan Ayah
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")

Distribusi kelas pada training dan testing
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")

Ringkasan kecocokan model penuh
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")

Odds ratio dan p-value — Model Penuh
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_birth

ringkasan_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")

Perbandingan ringkasan Model Penuh vs Reduksi
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).

cat("== VIF Model Penuh ==\n")

## == VIF Model Penuh ==

vif(heart_fit)

##                  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

cat("\n== VIF Model Reduksi ==\n")

## 
## == VIF Model Reduksi ==

vif(heart_fitB)

##                  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 ==

lrtest(null_model, heart_fit)

cat("\n== LRT: Model Reduksi vs Null ==\n")

## 
## == LRT: Model Reduksi vs Null ==

lrtest(null_model, heart_fitB)

cat("\n== LRT: Model Reduksi vs Model Penuh ==\n")

## 
## == LRT: Model Reduksi vs Model Penuh ==

lrtest(heart_fitB, heart_fit)

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.

cat("== Hosmer-Lemeshow — Model Penuh ==\n")

## == Hosmer-Lemeshow — Model Penuh ==

hoslem_penuh <- hoslem.test(train_data$caesar_bin, fitted(heart_fit), g = 10)
print(hoslem_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

cat("\n== Hosmer-Lemeshow — Model Reduksi ==\n")

## 
## == 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"
)

Nagelkerke R² kedua model
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)
}

or_tabel(heart_fitB, "Odds Ratio — Model Reduksi / Model Final (stepAIC)")

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)")

Contoh peluang prediksi pada data testing (Model Reduksi)
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"
)

Nilai AUC kedua model
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)"
)

Threshold optimal Model Reduksi (Indeks Youden)
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")

Confusion Matrix — Model Reduksi, Threshold 0,50
	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), ")"))

Confusion Matrix — Model Reduksi, Threshold Optimal (0.401)
	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")

Perbandingan metrik evaluasi Model Reduksi — Testing
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)")

Perbandingan metrik Model Penuh vs Reduksi (threshold 0,5)
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, dan period_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

cat("Jumlah baris:", nrow(data_multi), "\n")

## Jumlah baris: 1998

glimpse(data_multi)

## 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)

Distribusi jumlah percobaan bunuh diri (variabel respons multinomial)
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")

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)

Uji VIF — Variabel Numerik (Multinomial)
	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)

Uji Box-Tidwell — Multinomial
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.662

2.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)

Ringkasan hasil regresi logistik multinomial (RRR = Relative Risk Ratio)
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 = Ya pada 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)

Confusion Matrix — Model Multinomial
	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)

Ringkasan Pemeriksaan Asumsi — Multinomial
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

glimpse(data_ordinal)

## 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)

Distribusi tingkat kepuasan mahasiswa terhadap fasilitas kampus
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)

Ringkasan skor kepuasan per fasilitas (skala 1–10)
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")

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) ==

print(ordinal::nominal_test(fit_clm))

## 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)

Matriks Cramér’s V antar prediktor (< 0.5 = tidak bermasalah)
	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.5492

3.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)

Ringkasan hasil regresi logistik ordinal (MASS::polr)
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::polr menggunakan 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 polr harus 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 levelDegree memiliki β = 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)

Likelihood Ratio Test dan Pseudo R² — Model Ordinal
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)

Confusion Matrix — Model Ordinal
	Sedang	Tinggi	Sangat Tinggi	Sum
Rendah	2	14	0	16
Sedang	14	56	3	73
Tinggi	6	98	6	110
Sangat Tinggi	2	68	11	81
Sum	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))

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

glimpse(data_ins)

## 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)

Distribusi jumlah anak yang ditanggung asuransi
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)

Statistik deskriptif variabel numerik
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()

4.4 Pemodelan

4.4.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)

Ringkasan hasil regresi Poisson — Model Penuh
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.4.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)

Perbandingan Model Penuh vs Reduksi
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.5 Uji Asumsi

4.5.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)

Uji Overdispersi: Pearson Dispersion Statistic — Model Reduksi
Dispersion Pearson	Interpretasi
1.3246	Tidak ada indikasi overdispersion berat — equidispersion terpenuhi

cat("\n== Formal Overdispersion Test (AER) ==\n")

## 
## == Formal Overdispersion Test (AER) ==

AER::dispersiontest(fit_pois_red, trafo = 1)

## 
##  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.5.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)

Pemeriksaan proporsi nilai nol
Proporsi Y = 0	Catatan
42.9%	Proporsi nol wajar (< 50%) — model Poisson standar sesuai

4.5.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)

Uji Goodness-of-Fit: Deviance Test
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.5.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)

VIF Model Penuh Poisson
	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

cat("\n== VIF Model Reduksi ==\n")

## 
## == 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)

VIF Model Reduksi Poisson
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.5.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 --

lrtest(null_pois, fit_pois)

cat("\n-- LRT: Null vs Model Reduksi --\n")

## 
## -- LRT: Null vs Model Reduksi --

lrtest(null_pois, fit_pois_red)

cat("\n-- LRT: Model Reduksi vs Model Penuh --\n")

## 
## -- LRT: Model Reduksi vs Model Penuh --

lrtest(fit_pois_red, fit_pois)

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.6 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)

IRR dan Interval Kepercayaan 95% — Model Reduksi (Model Final)
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.7 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()

4.8 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)

Perbandingan komprehensif Model Penuh vs Model Reduksi
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%")

Ringkasan komparatif keempat model regresi logistik
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::polr menggunakan 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.

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

Analisis Regresi Logistik: Biner, Multinomial, Ordinal, dan Poisson

Studi Kasus dari Empat Dataset Mendeley

Nama Anda

2026-06-01

1 Regresi Logistik Biner

1.1 Deskripsi Data

1.2 Load Library dan Data

1.3 Persiapan Data

1.3.1 Preprocessing dan Pemberian Label

1.3.2 Pengecekan Data

1.4 Statistik Deskriptif

1.4.1 Distribusi Variabel Prediktor

1.4.2 Tabulasi Silang: Proporsi Caesar per Kategori

1.5 Eksplorasi Visual

1.6 Pembagian Training dan Testing

1.7 Model Regresi Logistik Biner

1.7.1 Rumus Model Logistik

1.7.1.1 Interpretasi Odds Ratio

1.7.2 Model Penuh

1.8 Seleksi Variabel: Stepwise Backward (AIC)

1.9 Uji Asumsi

1.9.1 Multikolinearitas — VIF

1.9.2 Likelihood Ratio Test (LRT)

1.9.3 Hosmer-Lemeshow Goodness-of-Fit

1.9.4 Linearitas Log-Odds (Box-Tidwell)

1.9.5 Nagelkerke \(R^2\)

1.10 Odds Ratio Model Final (Model Reduksi)

1.11 Prediksi Probabilitas

1.12 Kurva ROC dan Threshold Optimal

1.12.1 Rumus Kurva ROC dan AUC

1.13 Confusion Matrix dan Metrik Evaluasi

1.13.1 Rumus Metrik Klasifikasi

1.13.2 Confusion Matrix — Threshold 0,50

1.13.3 Confusion Matrix — Threshold Optimal

1.13.4 Perbandingan Metrik

1.13.5 Perbandingan Model Penuh vs Reduksi

1.14 Distribusi Peluang Prediksi

1.15 Kesimpulan

1.16 Catatan Metodologis

2 Regresi Logistik Multinomial

2.1 Deskripsi Data

2.2 Distribusi Variabel Respons

2.3 Uji Asumsi

2.3.1 Asumsi 1 — Respons Bersifat Nominal

2.3.2 Asumsi 2 — Observasi Independen

2.3.3 Asumsi 3 — Tidak Ada Multikolinearitas Berat (VIF)

2.3.4 Asumsi 4 — Tidak Ada Perfect Separation

2.3.5 Asumsi 5 — Linearitas pada Skala Logit (Box-Tidwell)

2.3.6 Asumsi 6 — Independence of Irrelevant Alternatives (IIA)

2.4 Estimasi Model

2.5 Tabel Koefisien, RRR, dan Interval Kepercayaan

2.6 Prediksi dan Confusion Matrix

2.7 Visualisasi Probabilitas Prediksi

2.8 Ringkasan Asumsi

3 Regresi Logistik Ordinal

3.1 Deskripsi Data

3.2 Statistik Deskriptif

3.3 Visualisasi Deskriptif

3.4 Uji Asumsi

3.4.1 Asumsi 1 — Asumsi Proportional Odds (Uji Nominal)

3.4.2 Asumsi 2 — Tidak Ada Multikolinearitas Berat (Cramér’s V)

3.5 Estimasi Model

3.6 Tabel Koefisien, Odds Ratio, dan Interval Kepercayaan

3.7 Uji Model (Likelihood Ratio Test dan Pseudo R²)

3.8 Prediksi dan Confusion Matrix

3.9 Visualisasi

4 Regresi Poisson

4.1 Deskripsi Data

4.2 Statistik Deskriptif

4.3 Visualisasi Awal

4.4 Pemodelan

4.4.1 Model Penuh

4.4.2 Model Reduksi (Model Final)

4.5 Uji Asumsi

4.5.1 Asumsi 1 — Overdispersi (Pearson Dispersion)

4.5.2 Asumsi 2 — Proporsi Nilai Nol

4.5.3 Asumsi 3 — Goodness-of-Fit (Deviance Test)

4.5.4 Asumsi 4 — Multikolinearitas (VIF)

4.5.5 Asumsi 5 — Likelihood Ratio Test

4.6 IRR Model Reduksi (Model Final)