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

1 BAB 1. Regresi Logistik Biner

Bab ini membahas pemodelan peluang asteroid tergolong Hazardous menggunakan regresi logistik biner.

1.1 Gambaran Umum Analisis

Dokumen ini menyajikan analisis regresi logistik biner untuk memodelkan peluang suatu asteroid tergolong berbahaya atau tidak berbahaya berdasarkan karakteristik orbit dan ukuran asteroid.

Respons yang digunakan adalah:

\[ Y = \begin{cases} 1, & \text{jika asteroid tergolong Hazardous} \\ 0, & \text{jika asteroid tidak tergolong Hazardous} \end{cases} \]

Prediktor yang digunakan adalah:

1. Absolute Magnitude
2. Est Dia in KM(max)
3. Velocity km per sec
4. Miss Dist.(kilometers)
5. Eccentricity

1.2 Pendahuluan

Asteroid dekat bumi memiliki karakteristik fisik dan orbital yang dapat digunakan untuk menilai potensi bahayanya. Salah satu bentuk klasifikasi yang umum digunakan adalah apakah asteroid tersebut termasuk Hazardous atau Non-Hazardous.

Karena variabel respons hanya terdiri dari dua kategori, metode yang sesuai adalah regresi logistik biner. Model ini tidak langsung memodelkan nilai kategori, tetapi memodelkan peluang suatu asteroid masuk ke kategori berbahaya.

Tujuan analisis:
Mengetahui apakah variabel Absolute Magnitude, Est Dia in KM(max), Velocity km per sec, Miss Dist.(kilometers), dan Eccentricity dapat digunakan untuk menjelaskan peluang asteroid tergolong Hazardous.

1.3 Dasar Teori Regresi Logistik Biner

Regresi logistik biner digunakan ketika variabel respons memiliki dua kategori, misalnya:

\[ Y = 1 \quad \text{atau} \quad Y = 0 \]

Misalkan:

\[ \pi(x) = P(Y = 1 \mid X = x) \]

Model regresi logistik biner ditulis sebagai:

\[ \log \left(\frac{\pi(x)}{1-\pi(x)}\right) = \beta_0 + \beta_1X_1 + \beta_2X_2 + \cdots + \beta_pX_p \]

Bagian:

\[ \frac{\pi(x)}{1-\pi(x)} \]

disebut odds, yaitu perbandingan antara peluang kejadian sukses dengan peluang tidak sukses.

Jika koefisien \(\beta_j\) dieksponensialkan, maka diperoleh:

\[ OR_j = e^{\beta_j} \]

Nilai \(OR_j\) disebut odds ratio.

Cara membaca odds ratio:

Jika \(OR > 1\), maka kenaikan prediktor meningkatkan odds asteroid tergolong Hazardous.
Jika \(OR < 1\), maka kenaikan prediktor menurunkan odds asteroid tergolong Hazardous.
Jika \(OR = 1\), maka prediktor tidak mengubah odds secara berarti.

1.4 Import Data

data_raw <- readxl::read_excel("~/Salwa/SMT 4/ADK 1/tugas regresi logostik/nasa_logistik_biner.xlsx")

dim(data_raw)

## [1] 4687    6

names(data_raw)

## [1] "Hazardous"                    "Absolute Magnitude"          
## [3] "Est Dia in KM(max)"           "Relative Velocity km per sec"
## [5] "Miss Dist.(kilometers)"       "Eccentricity"

Dataset berhasil diimpor ke R. Pada tahap ini, nama kolom diperiksa terlebih dahulu agar variabel yang digunakan dalam model sesuai dengan nama kolom pada file Excel.

1.5 Pemeriksaan dan Penyesuaian Nama Variabel

Pada file Excel, variabel kecepatan dapat memiliki nama lengkap Relative Velocity km per sec. Dalam laporan ini, variabel tersebut tetap ditulis sebagai Velocity km per sec agar sesuai dengan tujuan analisis.

find_col <- function(data, candidates) {
  nm <- names(data)
  nm_lower <- tolower(trimws(nm))
  cand_lower <- tolower(trimws(candidates))

  exact <- match(cand_lower, nm_lower)
  exact <- exact[!is.na(exact)]
  if (length(exact) > 0) return(nm[exact[1]])

  for (cand in cand_lower) {
    idx <- grep(cand, nm_lower, fixed = TRUE)
    if (length(idx) > 0) return(nm[idx[1]])
  }

  stop(paste("Kolom tidak ditemukan. Kandidat:", paste(candidates, collapse = ", ")))
}

col_hazard <- find_col(data_raw, c("Hazardous"))
col_abs_mag <- find_col(data_raw, c("Absolute Magnitude"))
col_dia <- find_col(data_raw, c("Est Dia in KM(max)", "Est Dia in KM"))
col_vel <- find_col(data_raw, c("Velocity km per sec", "Relative Velocity km per sec"))
col_miss <- find_col(data_raw, c("Miss Dist.(kilometers)", "Miss Dist"))
col_ecc <- find_col(data_raw, c("Eccentricity"))

kolom_digunakan <- tibble::tibble(
  Variabel_Analisis = c(
    "Hazardous",
    "Absolute Magnitude",
    "Est Dia in KM(max)",
    "Velocity km per sec",
    "Miss Dist.(kilometers)",
    "Eccentricity"
  ),
  Nama_Kolom_di_Excel = c(
    col_hazard,
    col_abs_mag,
    col_dia,
    col_vel,
    col_miss,
    col_ecc
  )
)

kolom_digunakan %>%
  knitr::kable(caption = "Penyesuaian nama variabel analisis dengan nama kolom pada Excel") %>%
  kableExtra::kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Penyesuaian nama variabel analisis dengan nama kolom pada Excel

Variabel_Analisis	Nama_Kolom_di_Excel
Hazardous	Hazardous
Absolute Magnitude	Absolute Magnitude
Est Dia in KM(max)	Est Dia in KM(max)
Velocity km per sec	Relative Velocity km per sec
Miss Dist.(kilometers)	Miss Dist.(kilometers)
Eccentricity	Eccentricity

1.6 Cleaning Data

data_bin <- data_raw %>%
  dplyr::transmute(
    Hazardous_raw = .data[[col_hazard]],
    Absolute_Magnitude = suppressWarnings(as.numeric(.data[[col_abs_mag]])),
    Est_Dia_KM_max = suppressWarnings(as.numeric(.data[[col_dia]])),
    Velocity_km_per_sec = suppressWarnings(as.numeric(.data[[col_vel]])),
    Miss_Dist_kilometers = suppressWarnings(as.numeric(.data[[col_miss]])),
    Eccentricity = suppressWarnings(as.numeric(.data[[col_ecc]]))
  ) %>%
  dplyr::mutate(
    Hazardous_num = dplyr::case_when(
      Hazardous_raw %in% c(TRUE, "TRUE", "True", "true", "Yes", "YES", "1", 1) ~ 1,
      Hazardous_raw %in% c(FALSE, "FALSE", "False", "false", "No", "NO", "0", 0) ~ 0,
      TRUE ~ NA_real_
    ),
    Hazardous = factor(
      Hazardous_num,
      levels = c(0, 1),
      labels = c("Tidak Hazardous", "Hazardous")
    )
  ) %>%
  dplyr::filter(
    !is.na(Hazardous_num),
    !is.na(Absolute_Magnitude),
    !is.na(Est_Dia_KM_max),
    !is.na(Velocity_km_per_sec),
    !is.na(Miss_Dist_kilometers),
    !is.na(Eccentricity)
  )

glimpse(data_bin)

## Rows: 4,687
## Columns: 8
## $ Hazardous_raw        <chr> "TRUE", "FALSE", "TRUE", "FALSE", "TRUE", "FALSE"…
## $ Absolute_Magnitude   <dbl> 21.6, 21.3, 20.3, 27.4, 21.6, 19.6, 19.6, 19.2, 1…
## $ Est_Dia_KM_max       <dbl> 0.28447230, 0.32661790, 0.51765448, 0.01968067, 0…
## $ Velocity_km_per_sec  <dbl> 6.115834, 18.113985, 7.590711, 11.173874, 9.84083…
## $ Miss_Dist_kilometers <dbl> 62753692, 57298148, 7622912, 42683616, 61010824, …
## $ Eccentricity         <dbl> 0.4255491, 0.3516743, 0.3482483, 0.2165783, 0.210…
## $ Hazardous_num        <dbl> 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0…
## $ Hazardous            <fct> Hazardous, Tidak Hazardous, Hazardous, Tidak Haza…

Data dibersihkan dengan cara memilih variabel yang diperlukan, mengubah variabel numerik ke format angka, serta mengubah respons Hazardous menjadi kode biner:
\(Y = 1\) untuk asteroid berbahaya dan \(Y = 0\) untuk asteroid tidak berbahaya.

1.7 Statistik Deskriptif

1.7.1 Distribusi Variabel Respons

tab_respon <- data_bin %>%
  dplyr::count(Hazardous) %>%
  dplyr::mutate(
    Proporsi = n / sum(n),
    Persentase = scales::percent(Proporsi, accuracy = 0.1)
  )

tab_respon %>%
  knitr::kable(caption = "Distribusi variabel respons Hazardous") %>%
  kableExtra::kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Distribusi variabel respons Hazardous

Hazardous	n	Proporsi	Persentase
Tidak Hazardous	3932	0.8389162	83.9%
Hazardous	755	0.1610838	16.1%

ggplot(tab_respon, aes(x = Hazardous, y = Proporsi, fill = Hazardous)) +
  geom_col(width = 0.65, alpha = 0.92) +
  geom_text(aes(label = Persentase), vjust = -0.4, fontface = "bold") +
  scale_y_continuous(labels = scales::percent_format(), limits = c(0, max(tab_respon$Proporsi) + 0.12)) +
  scale_fill_manual(values = c("#2563eb", "#dc2626")) +
  labs(
    title = "Distribusi Status Hazardous",
    subtitle = "Proporsi asteroid berdasarkan kategori berbahaya dan tidak berbahaya",
    x = NULL,
    y = "Proporsi"
  ) +
  theme_minimal(base_size = 13) +
  theme(legend.position = "none")

Tabel dan grafik di atas menunjukkan komposisi asteroid berdasarkan status Hazardous. Pemeriksaan ini penting karena regresi logistik biner membutuhkan dua kategori respons yang jelas. Jika salah satu kategori terlalu sedikit, model dapat menjadi kurang stabil.

1.7.2 Ringkasan Prediktor Numerik

summary_prediktor <- data_bin %>%
  dplyr::select(
    Absolute_Magnitude,
    Est_Dia_KM_max,
    Velocity_km_per_sec,
    Miss_Dist_kilometers,
    Eccentricity
  ) %>%
  dplyr::summarise(
    dplyr::across(
      dplyr::everything(),
      list(
        Min = ~min(.x, na.rm = TRUE),
        Mean = ~mean(.x, na.rm = TRUE),
        Median = ~median(.x, na.rm = TRUE),
        Max = ~max(.x, na.rm = TRUE),
        SD = ~sd(.x, na.rm = TRUE)
      )
    )
  ) %>%
  tidyr::pivot_longer(
    cols = dplyr::everything(),
    names_to = c("Variabel", ".value"),
    names_pattern = "(.+)_(Min|Mean|Median|Max|SD)"
  )

summary_prediktor %>%
  dplyr::mutate(
    dplyr::across(dplyr::where(is.numeric), ~round(.x, 4))
  ) %>%
  knitr::kable(caption = "Statistik deskriptif prediktor numerik") %>%
  kableExtra::kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = TRUE
  )

Statistik deskriptif prediktor numerik

Variabel	Min	Mean	Median	Max	SD
Absolute_Magnitude	11.1600	2.226790e+01	2.190000e+01	3.21000e+01	2.89100e+00
Est_Dia_KM_max	0.0023	4.575000e-01	2.478000e-01	3.48369e+01	8.26400e-01
Velocity_km_per_sec	0.3355	1.397080e+01	1.291790e+01	4.46337e+01	7.29320e+00
Miss_Dist_kilometers	26609.8867	3.841347e+07	3.964771e+07	7.47816e+07	2.18111e+07
Eccentricity	0.0075	3.826000e-01	3.725000e-01	9.60300e-01	1.80400e-01

Ringkasan statistik digunakan untuk melihat rentang nilai tiap prediktor. Variabel seperti Miss Dist.(kilometers) biasanya memiliki skala yang sangat besar, sehingga pada tahap pemodelan variabel numerik akan distandarisasi agar estimasi model lebih stabil dan mudah dibandingkan.

1.8 Visualisasi Eksploratif

data_long <- data_bin %>%
  dplyr::select(
    Hazardous,
    Absolute_Magnitude,
    Est_Dia_KM_max,
    Velocity_km_per_sec,
    Miss_Dist_kilometers,
    Eccentricity
  ) %>%
  dplyr::mutate(
    Miss_Dist_kilometers = log10(Miss_Dist_kilometers + 1),
    Est_Dia_KM_max = log10(Est_Dia_KM_max + 1)
  ) %>%
  tidyr::pivot_longer(
    cols = -Hazardous,
    names_to = "Prediktor",
    values_to = "Nilai"
  )

ggplot(data_long, aes(x = Hazardous, y = Nilai, fill = Hazardous)) +
  geom_boxplot(alpha = 0.82, outlier.alpha = 0.35) +
  facet_wrap(~ Prediktor, scales = "free_y") +
  scale_fill_manual(values = c("#2563eb", "#dc2626")) +
  labs(
    title = "Perbandingan Prediktor Berdasarkan Status Hazardous",
    subtitle = "Beberapa variabel diskalakan log agar pola lebih mudah dibaca",
    x = NULL,
    y = "Nilai prediktor"
  ) +
  theme_minimal(base_size = 12.5) +
  theme(legend.position = "none")

Visualisasi eksploratif membantu melihat apakah distribusi prediktor berbeda antara asteroid Hazardous dan Tidak Hazardous. Jika pola distribusi antar-kategori terlihat berbeda, prediktor tersebut berpotensi membantu klasifikasi dalam model regresi logistik biner.

1.9 Pembagian Data Training dan Testing

Data dibagi menjadi data training dan testing. Data training digunakan untuk membangun model, sedangkan data testing digunakan untuk mengevaluasi kemampuan prediksi model pada data yang tidak digunakan saat estimasi.

set.seed(123)

idx_train <- sample(seq_len(nrow(data_bin)), size = floor(0.8 * nrow(data_bin)))

data_train <- data_bin[idx_train, ]
data_test <- data_bin[-idx_train, ]

dim(data_train)

## [1] 3749    8

dim(data_test)

## [1] 938   8

Pembagian data dilakukan dengan proporsi sekitar 80% untuk training dan 20% untuk testing. Tujuannya agar evaluasi model tidak hanya dilakukan pada data yang sama dengan data pembentukan model.

1.10 Standardisasi Prediktor

Karena prediktor memiliki skala yang sangat berbeda, dilakukan standardisasi pada variabel numerik.

\[ Z = \frac{X-\bar{X}}{s_X} \]

Standardisasi membuat setiap prediktor berada pada skala yang sebanding.

prediktor_num <- c(
  "Absolute_Magnitude",
  "Est_Dia_KM_max",
  "Velocity_km_per_sec",
  "Miss_Dist_kilometers",
  "Eccentricity"
)

center_train <- sapply(data_train[prediktor_num], mean, na.rm = TRUE)
scale_train <- sapply(data_train[prediktor_num], sd, na.rm = TRUE)

scale_train[scale_train == 0 | is.na(scale_train)] <- 1

standardize_data <- function(data, center, scale) {
  data_out <- data
  for (v in names(center)) {
    data_out[[paste0(v, "_z")]] <- (data_out[[v]] - center[[v]]) / scale[[v]]
  }
  data_out
}

data_train_z <- standardize_data(data_train, center_train, scale_train)
data_test_z <- standardize_data(data_test, center_train, scale_train)

head(data_train_z)

Standardisasi tidak mengubah hubungan dasar antarvariabel, tetapi membantu model lebih stabil karena seluruh prediktor berada pada skala yang relatif sama. Interpretasi koefisien setelah standardisasi dibaca sebagai perubahan odds akibat kenaikan 1 simpangan baku pada prediktor.

1.11 Estimasi Model Regresi Logistik Biner

Model yang digunakan adalah:

\[ \log \left(\frac{\pi(x)}{1-\pi(x)}\right) = \beta_0 + \beta_1Z_1 + \beta_2Z_2 + \beta_3Z_3 + \beta_4Z_4 + \beta_5Z_5 \]

dengan:

\[ \pi(x)=P(\text{Hazardous}=1 \mid X) \]

fit_logit <- tryCatch(
  glm(
    Hazardous_num ~ Absolute_Magnitude_z + Est_Dia_KM_max_z +
      Velocity_km_per_sec_z + Miss_Dist_kilometers_z + Eccentricity_z,
    data = data_train_z,
    family = binomial(link = "logit")
  ),
  error = function(e) {
    message("Model gagal diestimasi: ", e$message)
    NULL
  }
)

if (!is.null(fit_logit)) {
  summary(fit_logit)
} else {
  cat("Model tidak dapat diestimasi. Periksa kembali data dan kategori respons.")
}

## 
## Call:
## glm(formula = Hazardous_num ~ Absolute_Magnitude_z + Est_Dia_KM_max_z + 
##     Velocity_km_per_sec_z + Miss_Dist_kilometers_z + Eccentricity_z, 
##     family = binomial(link = "logit"), data = data_train_z)
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            -2.62562    0.09877 -26.582  < 2e-16 ***
## Absolute_Magnitude_z   -2.80502    0.18242 -15.376  < 2e-16 ***
## Est_Dia_KM_max_z       -1.74710    0.17126 -10.201  < 2e-16 ***
## Velocity_km_per_sec_z   0.24210    0.05583   4.336 1.45e-05 ***
## Miss_Dist_kilometers_z -0.43915    0.05899  -7.445 9.73e-14 ***
## Eccentricity_z          0.02368    0.05763   0.411    0.681    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 3283.9  on 3748  degrees of freedom
## Residual deviance: 2586.7  on 3743  degrees of freedom
## AIC: 2598.7
## 
## Number of Fisher Scoring iterations: 7

Model regresi logistik biner diestimasi menggunakan fungsi glm() dengan family binomial(link = "logit"). Model ini digunakan untuk memprediksi peluang asteroid tergolong Hazardous berdasarkan lima prediktor yang telah distandarisasi.

1.12 Ringkasan Koefisien, Odds Ratio, dan Signifikansi

if (!is.null(fit_logit)) {
  coef_tab <- as.data.frame(coef(summary(fit_logit))) %>%
    tibble::rownames_to_column("Parameter") %>%
    dplyr::rename(
      Estimate = Estimate,
      Std_Error = `Std. Error`,
      z_value = `z value`,
      p_value = `Pr(>|z|)`
    ) %>%
    dplyr::mutate(
      Odds_Ratio = exp(Estimate),
      CI_low = exp(Estimate - 1.96 * Std_Error),
      CI_high = exp(Estimate + 1.96 * Std_Error),
      Perubahan_Odds_Persen = 100 * (Odds_Ratio - 1),
      Keputusan = dplyr::case_when(
        p_value < 0.001 ~ "Signifikan pada 0.1%",
        p_value < 0.01 ~ "Signifikan pada 1%",
        p_value < 0.05 ~ "Signifikan pada 5%",
        p_value < 0.10 ~ "Signifikan pada 10%",
        TRUE ~ "Tidak signifikan"
      )
    )

  coef_tab %>%
    dplyr::mutate(
      dplyr::across(
        c(Estimate, Std_Error, z_value, p_value, Odds_Ratio, CI_low, CI_high, Perubahan_Odds_Persen),
        ~round(.x, 4)
      )
    ) %>%
    knitr::kable(
      caption = "Ringkasan koefisien regresi logistik biner",
      col.names = c(
        "Parameter", "Estimate", "SE", "z-value", "p-value",
        "Odds Ratio", "CI 2.5%", "CI 97.5%", "Perubahan Odds (%)", "Keputusan"
      )
    ) %>%
    kableExtra::kable_styling(
      bootstrap_options = c("striped", "hover", "condensed", "responsive"),
      full_width = TRUE
    )
}

Ringkasan koefisien regresi logistik biner

Parameter	Estimate	SE	z-value	p-value	Odds Ratio	CI 2.5%	CI 97.5%	Perubahan Odds (%)	Keputusan
(Intercept)	-2.6256	0.0988	-26.5821	0.0000	0.0724	0.0597	0.0879	-92.7605	Signifikan pada 0.1%
Absolute_Magnitude_z	-2.8050	0.1824	-15.3765	0.0000	0.0605	0.0423	0.0865	-93.9494	Signifikan pada 0.1%
Est_Dia_KM_max_z	-1.7471	0.1713	-10.2015	0.0000	0.1743	0.1246	0.2438	-82.5722	Signifikan pada 0.1%
Velocity_km_per_sec_z	0.2421	0.0558	4.3363	0.0000	1.2739	1.1419	1.4212	27.3927	Signifikan pada 0.1%
Miss_Dist_kilometers_z	-0.4392	0.0590	-7.4446	0.0000	0.6446	0.5742	0.7236	-35.5418	Signifikan pada 0.1%
Eccentricity_z	0.0237	0.0576	0.4108	0.6812	1.0240	0.9146	1.1464	2.3959	Tidak signifikan

Tabel di atas menampilkan koefisien model, nilai signifikansi, dan odds ratio. Karena prediktor telah distandarisasi, maka odds ratio menunjukkan perubahan odds asteroid menjadi Hazardous untuk setiap kenaikan 1 simpangan baku pada prediktor, dengan prediktor lain dianggap konstan.

1.13 Interpretasi Variabel Secara Otomatis

if (!is.null(fit_logit)) {
  coef_interpretasi <- coef_tab %>%
    dplyr::filter(Parameter != "(Intercept)") %>%
    dplyr::mutate(
      Nama = dplyr::case_when(
        Parameter == "Absolute_Magnitude_z" ~ "Absolute Magnitude",
        Parameter == "Est_Dia_KM_max_z" ~ "Est Dia in KM(max)",
        Parameter == "Velocity_km_per_sec_z" ~ "Velocity km per sec",
        Parameter == "Miss_Dist_kilometers_z" ~ "Miss Dist.(kilometers)",
        Parameter == "Eccentricity_z" ~ "Eccentricity",
        TRUE ~ Parameter
      ),
      Arah = ifelse(Odds_Ratio > 1, "meningkatkan", "menurunkan"),
      Persen = abs(round(Perubahan_Odds_Persen, 2)),
      Signif = ifelse(p_value < 0.05, "signifikan", "tidak signifikan")
    )

  cat("<div class='interpretasi'>")
  cat("<strong>Interpretasi hasil model:</strong><br><br>")

  for (i in seq_len(nrow(coef_interpretasi))) {
    cat(
      paste0(
        "<p><strong>", coef_interpretasi$Nama[i], "</strong>: ",
        "kenaikan 1 simpangan baku pada variabel ini diperkirakan ",
        coef_interpretasi$Arah[i], " odds asteroid tergolong Hazardous sekitar ",
        coef_interpretasi$Persen[i], "%, dengan variabel lain konstan. ",
        "Secara statistik, pengaruh variabel ini tergolong <strong>",
        coef_interpretasi$Signif[i], "</strong> pada taraf 5%.</p>"
      )
    )
  }

  cat("</div>")
}

Interpretasi hasil model:

Absolute Magnitude: kenaikan 1 simpangan baku pada variabel ini diperkirakan menurunkan odds asteroid tergolong Hazardous sekitar 93.95%, dengan variabel lain konstan. Secara statistik, pengaruh variabel ini tergolong signifikan pada taraf 5%.

Est Dia in KM(max): kenaikan 1 simpangan baku pada variabel ini diperkirakan menurunkan odds asteroid tergolong Hazardous sekitar 82.57%, dengan variabel lain konstan. Secara statistik, pengaruh variabel ini tergolong signifikan pada taraf 5%.

Velocity km per sec: kenaikan 1 simpangan baku pada variabel ini diperkirakan meningkatkan odds asteroid tergolong Hazardous sekitar 27.39%, dengan variabel lain konstan. Secara statistik, pengaruh variabel ini tergolong signifikan pada taraf 5%.

Miss Dist.(kilometers): kenaikan 1 simpangan baku pada variabel ini diperkirakan menurunkan odds asteroid tergolong Hazardous sekitar 35.54%, dengan variabel lain konstan. Secara statistik, pengaruh variabel ini tergolong signifikan pada taraf 5%.

Eccentricity: kenaikan 1 simpangan baku pada variabel ini diperkirakan meningkatkan odds asteroid tergolong Hazardous sekitar 2.4%, dengan variabel lain konstan. Secara statistik, pengaruh variabel ini tergolong tidak signifikan pada taraf 5%.

1.14 Uji Kesesuaian Model Secara Keseluruhan

Model penuh dibandingkan dengan model null. Model null hanya berisi intercept tanpa prediktor.

Hipotesis:

\[ H_0 : \beta_1 = \beta_2 = \cdots = \beta_p = 0 \]

\[ H_1 : \text{minimal ada satu } \beta_j \neq 0 \]

if (!is.null(fit_logit)) {
  fit_null <- glm(
    Hazardous_num ~ 1,
    data = data_train_z,
    family = binomial(link = "logit")
  )

  lr_test <- anova(fit_null, fit_logit, test = "Chisq")

  lr_test %>%
    knitr::kable(caption = "Likelihood Ratio Test antara model null dan model penuh") %>%
    kableExtra::kable_styling(
      bootstrap_options = c("striped", "hover", "condensed", "responsive"),
      full_width = FALSE
    )
}

Likelihood Ratio Test antara model null dan model penuh

Resid. Df	Resid. Dev	Df	Deviance	Pr(>Chi)
3748	3283.880	NA	NA	NA
3743	2586.706	5	697.1745	0

Likelihood Ratio Test digunakan untuk melihat apakah model dengan prediktor lebih baik dibandingkan model tanpa prediktor. Jika p-value kecil, maka secara keseluruhan prediktor yang digunakan mampu memberikan tambahan informasi dalam menjelaskan status Hazardous.

1.15 Pseudo R-Squared

Pada regresi logistik, ukuran \(R^2\) tidak sama seperti regresi linear. Salah satu ukuran yang umum digunakan adalah McFadden Pseudo \(R^2\):

\[ R^2_{McFadden} = 1-\frac{\log L_{\text{model}}}{\log L_{\text{null}}} \]

if (!is.null(fit_logit)) {
  loglik_model <- as.numeric(logLik(fit_logit))
  loglik_null <- as.numeric(logLik(fit_null))
  pseudo_r2 <- 1 - (loglik_model / loglik_null)

  tibble::tibble(
    `Log Likelihood Model` = loglik_model,
    `Log Likelihood Null` = loglik_null,
    `McFadden Pseudo R2` = pseudo_r2
  ) %>%
    dplyr::mutate(dplyr::across(dplyr::everything(), ~round(.x, 4))) %>%
    knitr::kable(caption = "Nilai McFadden Pseudo R-Squared") %>%
    kableExtra::kable_styling(
      bootstrap_options = c("striped", "hover", "condensed", "responsive"),
      full_width = FALSE
    )
}

Nilai McFadden Pseudo R-Squared

Log Likelihood Model	Log Likelihood Null	McFadden Pseudo R2
-1293.353	-1641.94	0.2123

Pseudo \(R^2\) digunakan sebagai ukuran ringkas kemampuan model dalam memperbaiki model null. Nilainya tidak harus dibaca seperti \(R^2\) pada regresi linear, tetapi semakin besar nilainya menunjukkan model lebih informatif dibandingkan model tanpa prediktor.

1.16 Pemeriksaan Multikolinearitas Sederhana

Multikolinearitas terjadi ketika antar-prediktor memiliki hubungan linear yang kuat. Salah satu ukuran yang sering digunakan adalah Variance Inflation Factor (VIF).

\[ VIF_j = \frac{1}{1-R_j^2} \]

calc_vif <- function(data, vars) {
  out <- lapply(vars, function(v) {
    others <- setdiff(vars, v)
    form <- as.formula(paste(v, "~", paste(others, collapse = " + ")))
    fit <- lm(form, data = data)
    r2 <- summary(fit)$r.squared
    vif <- 1 / (1 - r2)
    data.frame(Variabel = v, VIF = vif)
  })
  dplyr::bind_rows(out)
}

vif_tab <- calc_vif(
  data_train_z,
  c(
    "Absolute_Magnitude_z",
    "Est_Dia_KM_max_z",
    "Velocity_km_per_sec_z",
    "Miss_Dist_kilometers_z",
    "Eccentricity_z"
  )
) %>%
  dplyr::mutate(
    Variabel = dplyr::case_when(
      Variabel == "Absolute_Magnitude_z" ~ "Absolute Magnitude",
      Variabel == "Est_Dia_KM_max_z" ~ "Est Dia in KM(max)",
      Variabel == "Velocity_km_per_sec_z" ~ "Velocity km per sec",
      Variabel == "Miss_Dist_kilometers_z" ~ "Miss Dist.(kilometers)",
      Variabel == "Eccentricity_z" ~ "Eccentricity",
      TRUE ~ Variabel
    ),
    VIF = round(VIF, 4),
    Keterangan = dplyr::case_when(
      VIF < 5 ~ "Aman",
      VIF < 10 ~ "Perlu diperhatikan",
      TRUE ~ "Indikasi multikolinearitas tinggi"
    )
  )

vif_tab %>%
  knitr::kable(caption = "Pemeriksaan multikolinearitas menggunakan VIF") %>%
  kableExtra::kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Pemeriksaan multikolinearitas menggunakan VIF

Variabel	VIF	Keterangan
Absolute Magnitude	1.8961	Aman
Est Dia in KM(max)	1.5513	Aman
Velocity km per sec	1.5882	Aman
Miss Dist.(kilometers)	1.3063	Aman
Eccentricity	1.5028	Aman

VIF digunakan untuk melihat apakah ada prediktor yang terlalu kuat berkorelasi dengan prediktor lain. Secara praktis, VIF di bawah 5 biasanya masih dianggap aman, sedangkan nilai yang sangat tinggi menandakan potensi multikolinearitas.

1.17 Prediksi Peluang pada Data Testing

if (!is.null(fit_logit)) {
  data_test_z$prob_hazardous <- predict(
    fit_logit,
    newdata = data_test_z,
    type = "response"
  )

  data_test_z$pred_class <- ifelse(
    data_test_z$prob_hazardous >= 0.5,
    1,
    0
  )

  data_test_z$pred_label <- factor(
    data_test_z$pred_class,
    levels = c(0, 1),
    labels = c("Tidak Hazardous", "Hazardous")
  )

  data_test_z %>%
    dplyr::select(Hazardous, prob_hazardous, pred_label) %>%
    head(10) %>%
    dplyr::mutate(prob_hazardous = round(prob_hazardous, 4)) %>%
    knitr::kable(caption = "Contoh hasil prediksi peluang pada data testing") %>%
    kableExtra::kable_styling(
      bootstrap_options = c("striped", "hover", "condensed", "responsive"),
      full_width = FALSE
    )
}

Contoh hasil prediksi peluang pada data testing

Hazardous	prob_hazardous	pred_label
Tidak Hazardous	0.2623	Tidak Hazardous
Hazardous	0.3759	Tidak Hazardous
Tidak Hazardous	0.0543	Tidak Hazardous
Tidak Hazardous	0.1389	Tidak Hazardous
Tidak Hazardous	0.2352	Tidak Hazardous
Tidak Hazardous	0.0005	Tidak Hazardous
Tidak Hazardous	0.0077	Tidak Hazardous
Hazardous	0.2517	Tidak Hazardous
Tidak Hazardous	0.3843	Tidak Hazardous
Tidak Hazardous	0.2561	Tidak Hazardous

Nilai prob_hazardous menunjukkan peluang prediksi asteroid tergolong Hazardous. Jika peluang prediksi minimal 0.5, asteroid diklasifikasikan sebagai Hazardous. Jika kurang dari 0.5, asteroid diklasifikasikan sebagai Tidak Hazardous.

1.18 Confusion Matrix dan Ukuran Evaluasi

if (!is.null(fit_logit)) {
  conf_mat <- table(
    Aktual = data_test_z$Hazardous,
    Prediksi = data_test_z$pred_label
  )

  conf_mat

  TP <- sum(data_test_z$Hazardous_num == 1 & data_test_z$pred_class == 1)
  TN <- sum(data_test_z$Hazardous_num == 0 & data_test_z$pred_class == 0)
  FP <- sum(data_test_z$Hazardous_num == 0 & data_test_z$pred_class == 1)
  FN <- sum(data_test_z$Hazardous_num == 1 & data_test_z$pred_class == 0)

  accuracy <- (TP + TN) / (TP + TN + FP + FN)
  sensitivity <- ifelse((TP + FN) == 0, NA, TP / (TP + FN))
  specificity <- ifelse((TN + FP) == 0, NA, TN / (TN + FP))
  precision <- ifelse((TP + FP) == 0, NA, TP / (TP + FP))

  eval_tab <- tibble::tibble(
    Ukuran = c("Accuracy", "Sensitivity", "Specificity", "Precision"),
    Nilai = c(accuracy, sensitivity, specificity, precision)
  )

  eval_tab %>%
    dplyr::mutate(Nilai = scales::percent(Nilai, accuracy = 0.1)) %>%
    knitr::kable(caption = "Ukuran evaluasi klasifikasi pada data testing") %>%
    kableExtra::kable_styling(
      bootstrap_options = c("striped", "hover", "condensed", "responsive"),
      full_width = FALSE
    )
}

Ukuran evaluasi klasifikasi pada data testing

Ukuran	Nilai
Accuracy	83.5%
Sensitivity	6.9%
Specificity	99.1%
Precision	61.1%

Confusion matrix membandingkan kategori aktual dengan kategori hasil prediksi.
Accuracy menunjukkan proporsi prediksi yang benar secara keseluruhan.
Sensitivity menunjukkan kemampuan model mengenali asteroid yang benar-benar Hazardous.
Specificity menunjukkan kemampuan model mengenali asteroid yang tidak Hazardous.

1.19 Kurva ROC dan AUC

AUC digunakan untuk mengukur kemampuan model membedakan kategori \(Y=1\) dan \(Y=0\). Nilai AUC berada antara 0.5 sampai 1.

if (!is.null(fit_logit)) {
  roc_data <- data.frame(
    y = data_test_z$Hazardous_num,
    prob = data_test_z$prob_hazardous
  ) %>%
    dplyr::arrange(dplyr::desc(prob))

  thresholds <- sort(unique(roc_data$prob), decreasing = TRUE)

  roc_points <- lapply(thresholds, function(th) {
    pred <- ifelse(roc_data$prob >= th, 1, 0)
    TP <- sum(roc_data$y == 1 & pred == 1)
    TN <- sum(roc_data$y == 0 & pred == 0)
    FP <- sum(roc_data$y == 0 & pred == 1)
    FN <- sum(roc_data$y == 1 & pred == 0)

    data.frame(
      threshold = th,
      TPR = ifelse((TP + FN) == 0, 0, TP / (TP + FN)),
      FPR = ifelse((FP + TN) == 0, 0, FP / (FP + TN))
    )
  }) %>%
    dplyr::bind_rows()

  roc_points <- dplyr::bind_rows(
    data.frame(threshold = Inf, TPR = 0, FPR = 0),
    roc_points,
    data.frame(threshold = -Inf, TPR = 1, FPR = 1)
  ) %>%
    dplyr::arrange(FPR, TPR)

  auc_value <- sum(diff(roc_points$FPR) *
                     (head(roc_points$TPR, -1) + tail(roc_points$TPR, -1)) / 2)

  ggplot(roc_points, aes(x = FPR, y = TPR)) +
    geom_line(color = "#2563eb", linewidth = 1.25) +
    geom_abline(linetype = "dashed", color = "gray45") +
    coord_equal() +
    labs(
      title = "Kurva ROC",
      subtitle = paste0("AUC = ", round(auc_value, 4)),
      x = "False Positive Rate (1 - Specificity)",
      y = "True Positive Rate (Sensitivity)"
    ) +
    theme_minimal(base_size = 13)
}

Kurva ROC menunjukkan performa model pada berbagai kemungkinan threshold. Nilai AUC yang semakin mendekati 1 menunjukkan kemampuan diskriminasi model yang semakin baik dalam membedakan asteroid Hazardous dan Tidak Hazardous.

1.20 Visualisasi Odds Ratio

if (!is.null(fit_logit)) {
  coef_tab %>%
    dplyr::filter(Parameter != "(Intercept)") %>%
    dplyr::mutate(
      Variabel = dplyr::case_when(
        Parameter == "Absolute_Magnitude_z" ~ "Absolute Magnitude",
        Parameter == "Est_Dia_KM_max_z" ~ "Est Dia in KM(max)",
        Parameter == "Velocity_km_per_sec_z" ~ "Velocity km per sec",
        Parameter == "Miss_Dist_kilometers_z" ~ "Miss Dist.(kilometers)",
        Parameter == "Eccentricity_z" ~ "Eccentricity",
        TRUE ~ Parameter
      )
    ) %>%
    ggplot(aes(x = reorder(Variabel, Odds_Ratio), y = Odds_Ratio)) +
    geom_hline(yintercept = 1, linetype = "dashed", color = "gray45") +
    geom_pointrange(aes(ymin = CI_low, ymax = CI_high), color = "#2563eb", linewidth = 0.7) +
    coord_flip() +
    labs(
      title = "Odds Ratio dan Interval Kepercayaan 95%",
      subtitle = "Prediktor telah distandarisasi",
      x = NULL,
      y = "Odds Ratio"
    ) +
    theme_minimal(base_size = 13)
}

Grafik odds ratio memudahkan pembaca melihat arah dan kekuatan pengaruh prediktor. Titik di atas 1 menunjukkan peningkatan odds asteroid menjadi Hazardous, sedangkan titik di bawah 1 menunjukkan penurunan odds. Garis interval yang melewati 1 menandakan efeknya belum kuat secara statistik pada taraf 5%.

1.21 Contoh Prediksi Individual

Misalkan ingin memprediksi peluang Hazardous untuk asteroid dengan karakteristik tertentu, kita dapat menggunakan data baru.

if (!is.null(fit_logit)) {
  data_baru <- data.frame(
    Absolute_Magnitude = median(data_bin$Absolute_Magnitude, na.rm = TRUE),
    Est_Dia_KM_max = median(data_bin$Est_Dia_KM_max, na.rm = TRUE),
    Velocity_km_per_sec = median(data_bin$Velocity_km_per_sec, na.rm = TRUE),
    Miss_Dist_kilometers = median(data_bin$Miss_Dist_kilometers, na.rm = TRUE),
    Eccentricity = median(data_bin$Eccentricity, na.rm = TRUE)
  )

  data_baru_z <- standardize_data(data_baru, center_train, scale_train)

  peluang_baru <- predict(fit_logit, newdata = data_baru_z, type = "response")

  data_baru %>%
    dplyr::mutate(
      Peluang_Hazardous = round(peluang_baru, 4),
      Klasifikasi = ifelse(peluang_baru >= 0.5, "Hazardous", "Tidak Hazardous")
    ) %>%
    knitr::kable(caption = "Contoh prediksi peluang Hazardous untuk data baru") %>%
    kableExtra::kable_styling(
      bootstrap_options = c("striped", "hover", "condensed", "responsive"),
      full_width = TRUE
    )
}

Contoh prediksi peluang Hazardous untuk data baru

Absolute_Magnitude	Est_Dia_KM_max	Velocity_km_per_sec	Miss_Dist_kilometers	Eccentricity	Peluang_Hazardous	Klasifikasi
21.9	0.247765	12.91789	39647712	0.3724502	0.1318	Tidak Hazardous

Contoh prediksi individual menunjukkan cara menggunakan model untuk menghitung peluang asteroid baru tergolong Hazardous. Nilai peluang ini dapat digunakan sebagai dasar klasifikasi dengan threshold tertentu, misalnya 0.5.

1.22 Kesimpulan

Berdasarkan analisis regresi logistik biner, model digunakan untuk menjelaskan peluang asteroid tergolong Hazardous berdasarkan lima prediktor numerik. Prediktor yang memiliki odds ratio di atas 1 cenderung meningkatkan peluang asteroid masuk kategori berbahaya, sedangkan prediktor dengan odds ratio di bawah 1 cenderung menurunkan peluang tersebut.

Secara umum, regresi logistik biner sesuai digunakan pada kasus ini karena variabel respons hanya memiliki dua kategori, yaitu Hazardous dan Tidak Hazardous. Hasil model dapat dibaca melalui koefisien, p-value, odds ratio, confusion matrix, serta AUC. Dengan demikian, analisis ini tidak hanya memberikan prediksi, tetapi juga membantu menjelaskan faktor-faktor yang berkaitan dengan status bahaya asteroid.

2 BAB 2. Regresi Logistik Multinomial

Bab ini membahas klasifikasi genre musik Spotify menggunakan regresi logistik multinomial.

2.1 Pendahuluan

Analisis ini bertujuan untuk memodelkan genre lagu Spotify menggunakan regresi logistik multinomial. Model ini digunakan karena variabel respon memiliki lebih dari dua kategori dan kategori tersebut bersifat nominal atau tidak memiliki urutan alami.

Variabel respon yang digunakan adalah track_genre, yang terdiri dari lima kategori:

• pop
• rock
• hip-hop
• classical
• r-n-b

Variabel prediktor yang digunakan adalah:

• danceability
• energy
• acousticness
• tempo
• valence

2.2 Membaca Data Excel

data <- read_excel("~/Salwa/SMT 4/ADK 1/tugas regresi logostik/spotify_filtered.xlsx")

glimpse(data)

## Rows: 5,000
## Columns: 21
## $ X                <dbl> 16000, 16001, 16002, 16003, 16004, 16005, 16006, 1600…
## $ track_id         <chr> "7wrYBASu0OoxoDErd4Edxd", "72HdutlIHBZJ7WT1xVAAZT", "…
## $ artists          <chr> "Bombay Jayashri", "Shankar;Ehsaan;Loy;Alisha Chinai;…
## $ album_name       <chr> "Rehnaa Hai Terre Dil Mein", "Bunty Aur Babli", "Hind…
## $ track_name       <chr> "Zara Zara", "Kajra Re", "Zara Zara - Lofi", "Vaseega…
## $ popularity       <dbl> 58, 59, 54, 68, 59, 62, 59, 0, 2, 1, 0, 1, 0, 2, 0, 1…
## $ duration_ms      <dbl> 298266, 482586, 219437, 299146, 387716, 298266, 30662…
## $ explicit         <chr> "False", "False", "False", "False", "False", "False",…
## $ danceability     <dbl> 0.6430, 0.4840, 0.6080, 0.6950, 0.5830, 0.6430, 0.642…
## $ energy           <dbl> 0.2680, 0.8980, 0.6380, 0.2930, 0.3080, 0.2680, 0.562…
## $ key              <dbl> 11, 0, 11, 11, 7, 11, 9, 0, 9, 9, 2, 2, 7, 2, 2, 0, 4…
## $ loudness         <dbl> -15.073, -4.132, -6.008, -16.278, -18.303, -15.073, -…
## $ mode             <dbl> 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1,…
## $ speechiness      <dbl> 0.0900, 0.1640, 0.0292, 0.0431, 0.0465, 0.0900, 0.060…
## $ acousticness     <dbl> 0.593, 0.365, 0.581, 0.596, 0.581, 0.593, 0.584, 0.96…
## $ instrumentalness <dbl> 2.06e-06, 0.00e+00, 1.72e-02, 1.58e-02, 1.06e-02, 2.0…
## $ liveness         <dbl> 0.3160, 0.0910, 0.4480, 0.1320, 0.2570, 0.3160, 0.104…
## $ valence          <dbl> 0.6200, 0.6800, 0.4390, 0.6370, 0.2410, 0.6200, 0.671…
## $ tempo            <dbl> 143.813, 91.975, 140.109, 143.804, 118.226, 143.813, …
## $ time_signature   <dbl> 4, 4, 4, 4, 4, 4, 4, 3, 4, 4, 4, 4, 3, 4, 3, 3, 4, 4,…
## $ track_genre      <chr> "classical", "classical", "classical", "classical", "…

Dataset yang digunakan merupakan data lagu Spotify. Variabel respon adalah genre lagu, sedangkan variabel prediktor merupakan karakteristik audio lagu seperti danceability, energy, acousticness, tempo, dan valence.

2.3 Memilih Variabel Penelitian

data_model <- data %>%
  dplyr::select(
    genre = track_genre,
    danceability,
    energy,
    acousticness,
    tempo,
    valence
  ) %>%
  dplyr::filter(
    genre %in% c("pop", "rock", "hip-hop", "classical", "r-n-b")
  ) %>%
  dplyr::mutate(
    genre = factor(
      genre,
      levels = c("classical", "pop", "rock", "hip-hop", "r-n-b")
    )
  ) %>%
  tidyr::drop_na()

glimpse(data_model)

## Rows: 5,000
## Columns: 6
## $ genre        <fct> classical, classical, classical, classical, classical, cl…
## $ danceability <dbl> 0.6430, 0.4840, 0.6080, 0.6950, 0.5830, 0.6430, 0.6420, 0…
## $ energy       <dbl> 0.2680, 0.8980, 0.6380, 0.2930, 0.3080, 0.2680, 0.5620, 0…
## $ acousticness <dbl> 0.593, 0.365, 0.581, 0.596, 0.581, 0.593, 0.584, 0.965, 0…
## $ tempo        <dbl> 143.813, 91.975, 140.109, 143.804, 118.226, 143.813, 149.…
## $ valence      <dbl> 0.6200, 0.6800, 0.4390, 0.6370, 0.2410, 0.6200, 0.6710, 0…

Kategori acuan yang digunakan adalah classical. Artinya, model akan membandingkan peluang genre pop, rock, hip-hop, dan r-n-b terhadap genre classical.

2.4 Statistik Deskriptif

data_model %>%
  dplyr::summarise(
    Jumlah_Data = n(),
    Rata_Rata_Danceability = mean(danceability),
    Rata_Rata_Energy = mean(energy),
    Rata_Rata_Acousticness = mean(acousticness),
    Rata_Rata_Tempo = mean(tempo),
    Rata_Rata_Valence = mean(valence)
  ) %>%
  kable(
    digits = 3,
    caption = "Statistik Deskriptif Variabel Numerik"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Statistik Deskriptif Variabel Numerik

Jumlah_Data	Rata_Rata_Danceability	Rata_Rata_Energy	Rata_Rata_Acousticness	Rata_Rata_Tempo	Rata_Rata_Valence
5000	0.581	0.559	0.408	119.201	0.522

data_model %>%
  dplyr::group_by(genre) %>%
  dplyr::summarise(
    Jumlah = n(),
    Mean_Danceability = mean(danceability),
    Mean_Energy = mean(energy),
    Mean_Acousticness = mean(acousticness),
    Mean_Tempo = mean(tempo),
    Mean_Valence = mean(valence),
    .groups = "drop"
  ) %>%
  kable(
    digits = 3,
    caption = "Rata-rata Karakteristik Audio Berdasarkan Genre"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = TRUE
  )

Rata-rata Karakteristik Audio Berdasarkan Genre

genre	Jumlah	Mean_Danceability	Mean_Energy	Mean_Acousticness	Mean_Tempo	Mean_Valence
classical	1000	0.382	0.190	0.920	107.950	0.381
pop	1000	0.630	0.606	0.344	120.927	0.506
rock	1000	0.544	0.679	0.209	126.321	0.539
hip-hop	1000	0.736	0.683	0.194	116.770	0.551
r-n-b	1000	0.614	0.638	0.371	124.039	0.633

2.5 Distribusi Variabel Respon

data_model %>%
  dplyr::count(genre) %>%
  dplyr::mutate(Proporsi = n / sum(n)) %>%
  kable(
    digits = 3,
    caption = "Distribusi Genre Lagu"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Distribusi Genre Lagu

genre	n	Proporsi
classical	1000	0.2
pop	1000	0.2
rock	1000	0.2
hip-hop	1000	0.2
r-n-b	1000	0.2

data_model %>%
  dplyr::count(genre) %>%
  dplyr::mutate(Proporsi = n / sum(n)) %>%
  ggplot(aes(x = genre, y = Proporsi, fill = genre)) +
  geom_col(width = 0.65, alpha = 0.95) +
  geom_text(
    aes(label = percent(Proporsi, accuracy = 0.1)),
    vjust = -0.4,
    fontface = "bold",
    color = "#243142"
  ) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
  scale_fill_manual(values = c(
    "classical" = "#5568b8",
    "pop" = "#d18b2f",
    "rock" = "#b84f5a",
    "hip-hop" = "#2f7f73",
    "r-n-b" = "#7c3aed"
  )) +
  labs(
    title = "Distribusi Genre Lagu Spotify",
    subtitle = "Variabel respon terdiri dari lima kategori nominal.",
    x = "Genre",
    y = "Proporsi",
    fill = "Genre"
  ) +
  theme_profesional()

2.6 Visualisasi Awal Karakteristik Audio

data_long <- data_model %>%
  tidyr::pivot_longer(
    cols = c(danceability, energy, acousticness, tempo, valence),
    names_to = "Variabel",
    values_to = "Nilai"
  )

ggplot(data_long, aes(x = genre, y = Nilai, fill = genre)) +
  geom_boxplot(alpha = 0.85, outlier.alpha = 0.25) +
  facet_wrap(~ Variabel, scales = "free_y") +
  scale_fill_manual(values = c(
    "classical" = "#5568b8",
    "pop" = "#d18b2f",
    "rock" = "#b84f5a",
    "hip-hop" = "#2f7f73",
    "r-n-b" = "#7c3aed"
  )) +
  labs(
    title = "Sebaran Karakteristik Audio Berdasarkan Genre",
    subtitle = "Boxplot digunakan untuk melihat perbedaan pola audio pada masing-masing genre.",
    x = "Genre",
    y = "Nilai",
    fill = "Genre"
  ) +
  theme_profesional() +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))

2.7 Rumusan Model Regresi Logistik Multinomial

Regresi logistik multinomial digunakan ketika variabel respon terdiri dari lebih dari dua kategori nominal. Misalkan kategori acuan adalah classical.

Untuk setiap kategori \(j\) selain kategori acuan, model dituliskan sebagai:

\[ \log\left(\frac{\pi_j}{\pi_{classical}}\right) = \beta_{0j} + \beta_{1j}danceability + \beta_{2j}energy + \beta_{3j}acousticness + \beta_{4j}tempo + \beta_{5j}valence \]

dengan:

\[ \pi_j = P(Y = j) \]

dan:

\[ j \in \{pop, rock, hip\text{-}hop, r\text{-}n\text{-}b\} \]

Secara lengkap, model membentuk empat persamaan logit:

\[ \log\left(\frac{\pi_{pop}}{\pi_{classical}}\right) = \beta_{0,pop} + \beta_{1,pop}danceability + \beta_{2,pop}energy + \beta_{3,pop}acousticness + \beta_{4,pop}tempo + \beta_{5,pop}valence \]

\[ \log\left(\frac{\pi_{rock}}{\pi_{classical}}\right) = \beta_{0,rock} + \beta_{1,rock}danceability + \beta_{2,rock}energy + \beta_{3,rock}acousticness + \beta_{4,rock}tempo + \beta_{5,rock}valence \]

\[ \log\left(\frac{\pi_{hip-hop}}{\pi_{classical}}\right) = \beta_{0,hip-hop} + \beta_{1,hip-hop}danceability + \beta_{2,hip-hop}energy + \beta_{3,hip-hop}acousticness + \beta_{4,hip-hop}tempo + \beta_{5,hip-hop}valence \]

\[ \log\left(\frac{\pi_{r-n-b}}{\pi_{classical}}\right) = \beta_{0,r-n-b} + \beta_{1,r-n-b}danceability + \beta_{2,r-n-b}energy + \beta_{3,r-n-b}acousticness + \beta_{4,r-n-b}tempo + \beta_{5,r-n-b}valence \]

Koefisien positif menunjukkan bahwa peningkatan variabel prediktor meningkatkan log-odds suatu genre dibandingkan genre classical. Koefisien negatif menunjukkan bahwa peningkatan variabel prediktor menurunkan log-odds genre tersebut dibandingkan genre classical.

2.8 Membagi Data Latih dan Data Uji

set.seed(123)

indeks_latih <- data_model %>%
  dplyr::mutate(id = dplyr::row_number()) %>%
  dplyr::group_by(genre) %>%
  dplyr::slice_sample(prop = 0.8) %>%
  dplyr::pull(id)

data_latih <- data_model[indeks_latih, ]
data_uji <- data_model[-indeks_latih, ]

dim(data_latih)

## [1] 4000    6

dim(data_uji)

## [1] 1000    6

Data dibagi menjadi 80% data latih dan 20% data uji. Pembagian dilakukan secara proporsional berdasarkan genre agar setiap kategori tetap terwakili pada data latih dan data uji.

2.9 Membentuk Model Regresi Logistik Multinomial

model_multinom <- nnet::multinom(
  genre ~ danceability + energy + acousticness + tempo + valence,
  data = data_latih,
  trace = FALSE,
  maxit = 1000
)

summary(model_multinom)

## Call:
## nnet::multinom(formula = genre ~ danceability + energy + acousticness + 
##     tempo + valence, data = data_latih, trace = FALSE, maxit = 1000)
## 
## Coefficients:
##         (Intercept) danceability   energy acousticness       tempo   valence
## pop       -1.208437    10.673277 4.633309    -6.252155 0.005873803 -4.739093
## rock       2.598124     4.285420 4.229088    -8.479957 0.002703960 -2.429884
## hip-hop   -7.118033    18.295148 7.072579    -7.404278 0.004620625 -5.636935
## r-n-b     -1.486820     7.557334 5.173339    -5.744653 0.005057915 -1.606806
## 
## Std. Errors:
##         (Intercept) danceability    energy acousticness       tempo   valence
## pop       0.7688838    0.7874017 0.5879389    0.5754970 0.002990815 0.4720482
## rock      0.7605027    0.7800742 0.5996610    0.5862316 0.002995161 0.4774287
## hip-hop   0.8787138    0.8859196 0.6494834    0.6039369 0.003301820 0.5042145
## r-n-b     0.7561353    0.7685413 0.5721998    0.5666061 0.002938021 0.4586678
## 
## Residual Deviance: 8528.958 
## AIC: 8576.958

2.10 Koefisien Model

koefisien <- coef(model_multinom)

round(koefisien, 4)

##         (Intercept) danceability energy acousticness  tempo valence
## pop         -1.2084      10.6733 4.6333      -6.2522 0.0059 -4.7391
## rock         2.5981       4.2854 4.2291      -8.4800 0.0027 -2.4299
## hip-hop     -7.1180      18.2951 7.0726      -7.4043 0.0046 -5.6369
## r-n-b       -1.4868       7.5573 5.1733      -5.7447 0.0051 -1.6068

Baris pada tabel koefisien menunjukkan genre yang dibandingkan terhadap kategori acuan classical. Kolom menunjukkan parameter model. Misalnya, koefisien danceability pada baris pop menunjukkan pengaruh danceability terhadap log-odds genre pop dibandingkan classical.

2.11 Uji Signifikansi Parameter

hasil_summary <- summary(model_multinom)

z_value <- hasil_summary$coefficients / hasil_summary$standard.errors
p_value <- 2 * (1 - pnorm(abs(z_value)))

tabel_uji <- data.frame(
  Genre = rep(rownames(hasil_summary$coefficients),
              each = ncol(hasil_summary$coefficients)),
  Variabel = rep(colnames(hasil_summary$coefficients),
                 times = nrow(hasil_summary$coefficients)),
  Koefisien = as.vector(t(hasil_summary$coefficients)),
  Std_Error = as.vector(t(hasil_summary$standard.errors)),
  Z_Value = as.vector(t(z_value)),
  P_Value = as.vector(t(p_value))
) %>%
  dplyr::mutate(
    Kesimpulan = ifelse(P_Value < 0.05, "Signifikan", "Tidak signifikan")
  )

tabel_uji %>%
  dplyr::mutate(
    dplyr::across(c(Koefisien, Std_Error, Z_Value, P_Value), ~ round(.x, 4))
  ) %>%
  kable(
    caption = "Uji Signifikansi Koefisien Model Multinomial"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = TRUE
  ) %>%
  row_spec(
    which(tabel_uji$P_Value < 0.05),
    bold = TRUE,
    color = "white",
    background = "#16a34a"
  )

Uji Signifikansi Koefisien Model Multinomial

Genre	Variabel	Koefisien	Std_Error	Z_Value	P_Value	Kesimpulan
pop	(Intercept)	-1.2084	0.7689	-1.5717	0.1160	Tidak signifikan
pop	danceability	10.6733	0.7874	13.5551	0.0000	Signifikan
pop	energy	4.6333	0.5879	7.8806	0.0000	Signifikan
pop	acousticness	-6.2522	0.5755	-10.8639	0.0000	Signifikan
pop	tempo	0.0059	0.0030	1.9639	0.0495	Signifikan
pop	valence	-4.7391	0.4720	-10.0394	0.0000	Signifikan
rock	(Intercept)	2.5981	0.7605	3.4163	0.0006	Signifikan
rock	danceability	4.2854	0.7801	5.4936	0.0000	Signifikan
rock	energy	4.2291	0.5997	7.0525	0.0000	Signifikan
rock	acousticness	-8.4800	0.5862	-14.4652	0.0000	Signifikan
rock	tempo	0.0027	0.0030	0.9028	0.3666	Tidak signifikan
rock	valence	-2.4299	0.4774	-5.0895	0.0000	Signifikan
hip-hop	(Intercept)	-7.1180	0.8787	-8.1005	0.0000	Signifikan
hip-hop	danceability	18.2951	0.8859	20.6510	0.0000	Signifikan
hip-hop	energy	7.0726	0.6495	10.8895	0.0000	Signifikan
hip-hop	acousticness	-7.4043	0.6039	-12.2600	0.0000	Signifikan
hip-hop	tempo	0.0046	0.0033	1.3994	0.1617	Tidak signifikan
hip-hop	valence	-5.6369	0.5042	-11.1796	0.0000	Signifikan
r-n-b	(Intercept)	-1.4868	0.7561	-1.9663	0.0493	Signifikan
r-n-b	danceability	7.5573	0.7685	9.8333	0.0000	Signifikan
r-n-b	energy	5.1733	0.5722	9.0411	0.0000	Signifikan
r-n-b	acousticness	-5.7447	0.5666	-10.1387	0.0000	Signifikan
r-n-b	tempo	0.0051	0.0029	1.7215	0.0852	Tidak signifikan
r-n-b	valence	-1.6068	0.4587	-3.5032	0.0005	Signifikan

2.12 Odds Ratio

Odds ratio diperoleh dengan mengeksponensialkan koefisien:

\[ OR = e^{\beta} \]

odds_ratio <- exp(coef(model_multinom))

tabel_or <- data.frame(
  Genre = rep(rownames(odds_ratio), each = ncol(odds_ratio)),
  Variabel = rep(colnames(odds_ratio), times = nrow(odds_ratio)),
  Odds_Ratio = as.vector(t(odds_ratio))
)

tabel_or %>%
  dplyr::mutate(Odds_Ratio = round(Odds_Ratio, 4)) %>%
  kable(
    caption = "Odds Ratio Model Regresi Logistik Multinomial"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = TRUE
  )

Odds Ratio Model Regresi Logistik Multinomial

Genre	Variabel	Odds_Ratio
pop	(Intercept)	2.987000e-01
pop	danceability	4.318623e+04
pop	energy	1.028539e+02
pop	acousticness	1.900000e-03
pop	tempo	1.005900e+00
pop	valence	8.700000e-03
rock	(Intercept)	1.343850e+01
rock	danceability	7.263310e+01
rock	energy	6.865460e+01
rock	acousticness	2.000000e-04
rock	tempo	1.002700e+00
rock	valence	8.800000e-02
hip-hop	(Intercept)	8.000000e-04
hip-hop	danceability	8.820272e+07
hip-hop	energy	1.179185e+03
hip-hop	acousticness	6.000000e-04
hip-hop	tempo	1.004600e+00
hip-hop	valence	3.600000e-03
r-n-b	(Intercept)	2.261000e-01
r-n-b	danceability	1.914733e+03
r-n-b	energy	1.765032e+02
r-n-b	acousticness	3.200000e-03
r-n-b	tempo	1.005100e+00
r-n-b	valence	2.005000e-01

Jika odds ratio lebih dari 1, maka peningkatan variabel prediktor meningkatkan odds suatu genre dibandingkan genre classical. Jika odds ratio kurang dari 1, maka peningkatan variabel prediktor menurunkan odds suatu genre dibandingkan genre classical.

2.13 Interpretasi Otomatis Odds Ratio

tabel_interpretasi <- tabel_or %>%
  dplyr::filter(Variabel != "(Intercept)") %>%
  dplyr::mutate(
    Interpretasi = dplyr::case_when(
      Odds_Ratio > 1 ~ paste0(
        "Pada genre ", Genre, ", variabel ", Variabel,
        " memiliki OR = ", round(Odds_Ratio, 4),
        ". Artinya, kenaikan variabel ini meningkatkan odds genre ",
        Genre, " dibandingkan classical sebesar ",
        round((Odds_Ratio - 1) * 100, 2),
        "%, dengan asumsi variabel lain konstan."
      ),
      Odds_Ratio < 1 ~ paste0(
        "Pada genre ", Genre, ", variabel ", Variabel,
        " memiliki OR = ", round(Odds_Ratio, 4),
        ". Artinya, kenaikan variabel ini menurunkan odds genre ",
        Genre, " dibandingkan classical sebesar ",
        round((1 - Odds_Ratio) * 100, 2),
        "%, dengan asumsi variabel lain konstan."
      ),
      TRUE ~ paste0(
        "Pada genre ", Genre, ", variabel ", Variabel,
        " memiliki OR = 1 sehingga tidak mengubah odds terhadap classical."
      )
    )
  ) %>%
  dplyr::left_join(
    tabel_uji %>% dplyr::select(Genre, Variabel, P_Value, Kesimpulan),
    by = c("Genre", "Variabel")
  )

tabel_interpretasi %>%
  dplyr::mutate(
    Odds_Ratio = round(Odds_Ratio, 4),
    P_Value = round(P_Value, 4)
  ) %>%
  dplyr::select(Genre, Variabel, Odds_Ratio, P_Value, Kesimpulan, Interpretasi) %>%
  kable(
    caption = "Interpretasi Odds Ratio Berdasarkan Genre"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = TRUE
  )

Interpretasi Odds Ratio Berdasarkan Genre

Genre	Variabel	Odds_Ratio	P_Value	Kesimpulan	Interpretasi
pop	danceability	4.318623e+04	0.0000	Signifikan	Pada genre pop, variabel danceability memiliki OR = 43186.2305. Artinya, kenaikan variabel ini meningkatkan odds genre pop dibandingkan classical sebesar 4318523.05%, dengan asumsi variabel lain konstan.
pop	energy	1.028539e+02	0.0000	Signifikan	Pada genre pop, variabel energy memiliki OR = 102.8539. Artinya, kenaikan variabel ini meningkatkan odds genre pop dibandingkan classical sebesar 10185.39%, dengan asumsi variabel lain konstan.
pop	acousticness	1.900000e-03	0.0000	Signifikan	Pada genre pop, variabel acousticness memiliki OR = 0.0019. Artinya, kenaikan variabel ini menurunkan odds genre pop dibandingkan classical sebesar 99.81%, dengan asumsi variabel lain konstan.
pop	tempo	1.005900e+00	0.0495	Signifikan	Pada genre pop, variabel tempo memiliki OR = 1.0059. Artinya, kenaikan variabel ini meningkatkan odds genre pop dibandingkan classical sebesar 0.59%, dengan asumsi variabel lain konstan.
pop	valence	8.700000e-03	0.0000	Signifikan	Pada genre pop, variabel valence memiliki OR = 0.0087. Artinya, kenaikan variabel ini menurunkan odds genre pop dibandingkan classical sebesar 99.13%, dengan asumsi variabel lain konstan.
rock	danceability	7.263310e+01	0.0000	Signifikan	Pada genre rock, variabel danceability memiliki OR = 72.6331. Artinya, kenaikan variabel ini meningkatkan odds genre rock dibandingkan classical sebesar 7163.31%, dengan asumsi variabel lain konstan.
rock	energy	6.865460e+01	0.0000	Signifikan	Pada genre rock, variabel energy memiliki OR = 68.6546. Artinya, kenaikan variabel ini meningkatkan odds genre rock dibandingkan classical sebesar 6765.46%, dengan asumsi variabel lain konstan.
rock	acousticness	2.000000e-04	0.0000	Signifikan	Pada genre rock, variabel acousticness memiliki OR = 2e-04. Artinya, kenaikan variabel ini menurunkan odds genre rock dibandingkan classical sebesar 99.98%, dengan asumsi variabel lain konstan.
rock	tempo	1.002700e+00	0.3666	Tidak signifikan	Pada genre rock, variabel tempo memiliki OR = 1.0027. Artinya, kenaikan variabel ini meningkatkan odds genre rock dibandingkan classical sebesar 0.27%, dengan asumsi variabel lain konstan.
rock	valence	8.800000e-02	0.0000	Signifikan	Pada genre rock, variabel valence memiliki OR = 0.088. Artinya, kenaikan variabel ini menurunkan odds genre rock dibandingkan classical sebesar 91.2%, dengan asumsi variabel lain konstan.
hip-hop	danceability	8.820272e+07	0.0000	Signifikan	Pada genre hip-hop, variabel danceability memiliki OR = 88202716.6462. Artinya, kenaikan variabel ini meningkatkan odds genre hip-hop dibandingkan classical sebesar 8820271564.62%, dengan asumsi variabel lain konstan.
hip-hop	energy	1.179185e+03	0.0000	Signifikan	Pada genre hip-hop, variabel energy memiliki OR = 1179.1851. Artinya, kenaikan variabel ini meningkatkan odds genre hip-hop dibandingkan classical sebesar 117818.51%, dengan asumsi variabel lain konstan.
hip-hop	acousticness	6.000000e-04	0.0000	Signifikan	Pada genre hip-hop, variabel acousticness memiliki OR = 6e-04. Artinya, kenaikan variabel ini menurunkan odds genre hip-hop dibandingkan classical sebesar 99.94%, dengan asumsi variabel lain konstan.
hip-hop	tempo	1.004600e+00	0.1617	Tidak signifikan	Pada genre hip-hop, variabel tempo memiliki OR = 1.0046. Artinya, kenaikan variabel ini meningkatkan odds genre hip-hop dibandingkan classical sebesar 0.46%, dengan asumsi variabel lain konstan.
hip-hop	valence	3.600000e-03	0.0000	Signifikan	Pada genre hip-hop, variabel valence memiliki OR = 0.0036. Artinya, kenaikan variabel ini menurunkan odds genre hip-hop dibandingkan classical sebesar 99.64%, dengan asumsi variabel lain konstan.
r-n-b	danceability	1.914733e+03	0.0000	Signifikan	Pada genre r-n-b, variabel danceability memiliki OR = 1914.7334. Artinya, kenaikan variabel ini meningkatkan odds genre r-n-b dibandingkan classical sebesar 191373.34%, dengan asumsi variabel lain konstan.
r-n-b	energy	1.765032e+02	0.0000	Signifikan	Pada genre r-n-b, variabel energy memiliki OR = 176.5032. Artinya, kenaikan variabel ini meningkatkan odds genre r-n-b dibandingkan classical sebesar 17550.32%, dengan asumsi variabel lain konstan.
r-n-b	acousticness	3.200000e-03	0.0000	Signifikan	Pada genre r-n-b, variabel acousticness memiliki OR = 0.0032. Artinya, kenaikan variabel ini menurunkan odds genre r-n-b dibandingkan classical sebesar 99.68%, dengan asumsi variabel lain konstan.
r-n-b	tempo	1.005100e+00	0.0852	Tidak signifikan	Pada genre r-n-b, variabel tempo memiliki OR = 1.0051. Artinya, kenaikan variabel ini meningkatkan odds genre r-n-b dibandingkan classical sebesar 0.51%, dengan asumsi variabel lain konstan.
r-n-b	valence	2.005000e-01	0.0005	Signifikan	Pada genre r-n-b, variabel valence memiliki OR = 0.2005. Artinya, kenaikan variabel ini menurunkan odds genre r-n-b dibandingkan classical sebesar 79.95%, dengan asumsi variabel lain konstan.

2.14 Uji Kelayakan Model dengan Likelihood Ratio Test

model_null <- nnet::multinom(
  genre ~ 1,
  data = data_latih,
  trace = FALSE,
  maxit = 1000
)

LR_stat <- 2 * (as.numeric(logLik(model_multinom)) - as.numeric(logLik(model_null)))
df_LR <- attr(logLik(model_multinom), "df") - attr(logLik(model_null), "df")
p_LR <- pchisq(LR_stat, df = df_LR, lower.tail = FALSE)

tibble(
  Statistik_LR = LR_stat,
  Derajat_Bebas = df_LR,
  P_Value = p_LR
) %>%
  kable(
    digits = 4,
    caption = "Likelihood Ratio Test Model Multinomial"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Likelihood Ratio Test Model Multinomial

Statistik_LR	Derajat_Bebas	P_Value
4346.545	20	0

Jika p-value pada Likelihood Ratio Test kurang dari 0,05, maka model dengan prediktor danceability, energy, acousticness, tempo, dan valence lebih baik dibandingkan model tanpa prediktor.

2.15 Prediksi pada Data Uji

prediksi_kelas <- predict(
  model_multinom,
  newdata = data_uji,
  type = "class"
)

prediksi_prob <- predict(
  model_multinom,
  newdata = data_uji,
  type = "probs"
)

head(prediksi_prob)

##    classical        pop       rock      hip-hop       r-n-b
## 1 0.11289667 0.33752637 0.15873562 5.719927e-02 0.333642076
## 2 0.01392852 0.39580204 0.13714403 1.517983e-01 0.301327121
## 3 0.03091195 0.30725951 0.14541834 1.014093e-01 0.415000898
## 4 0.91336123 0.02594349 0.01500451 8.293724e-04 0.044861399
## 5 0.96560849 0.01001403 0.01662610 3.675918e-05 0.007714622
## 6 0.96493229 0.01268953 0.01264323 7.302781e-05 0.009661922

data_hasil_prediksi <- data_uji %>%
  dplyr::bind_cols(as.data.frame(prediksi_prob)) %>%
  dplyr::mutate(Prediksi = prediksi_kelas)

head(data_hasil_prediksi)

2.16 Confusion Matrix dan Akurasi

conf_matrix <- table(
  Aktual = data_uji$genre,
  Prediksi = prediksi_kelas
)

conf_matrix

##            Prediksi
## Aktual      classical pop rock hip-hop r-n-b
##   classical       174  11    5       0    10
##   pop              17  45   35      69    34
##   rock              4  34  113      35    14
##   hip-hop           0  22   22     134    22
##   r-n-b            15  42   45      36    62

akurasi <- sum(diag(conf_matrix)) / sum(conf_matrix)

tibble(
  Akurasi = akurasi
) %>%
  dplyr::mutate(Akurasi = percent(Akurasi, accuracy = 0.01)) %>%
  kable(
    caption = "Akurasi Model Regresi Logistik Multinomial"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Akurasi Model Regresi Logistik Multinomial

Akurasi
52.80%

Confusion matrix menunjukkan perbandingan antara genre aktual dan genre hasil prediksi. Akurasi menunjukkan proporsi lagu yang berhasil diklasifikasikan dengan benar oleh model.

2.17 Contoh Prediksi Lagu Baru

Misalkan terdapat lagu baru dengan karakteristik:

• danceability = 0.70
• energy = 0.65
• acousticness = 0.20
• tempo = 120
• valence = 0.55

lagu_baru <- tibble(
  danceability = 0.70,
  energy = 0.65,
  acousticness = 0.20,
  tempo = 120,
  valence = 0.55
)

prob_lagu_baru <- predict(
  model_multinom,
  newdata = lagu_baru,
  type = "probs"
)

prob_lagu_baru

##    classical          pop         rock      hip-hop        r-n-b 
## 0.0006361376 0.2900611302 0.1788474300 0.3325239966 0.1979313055

tibble(
  Genre = names(prob_lagu_baru),
  Probabilitas = as.numeric(prob_lagu_baru)
) %>%
  dplyr::mutate(Probabilitas = percent(Probabilitas, accuracy = 0.01)) %>%
  kable(
    caption = "Probabilitas Prediksi Genre Lagu Baru"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Probabilitas Prediksi Genre Lagu Baru

Genre	Probabilitas
classical	0.06%
pop	29.01%
rock	17.88%
hip-hop	33.25%
r-n-b	19.79%

Genre dengan probabilitas terbesar merupakan genre yang paling mungkin untuk lagu baru tersebut berdasarkan karakteristik audio yang diberikan.

2.18 Visualisasi Pengaruh Danceability terhadap Probabilitas Genre

grid_dance <- tibble(
  danceability = seq(
    min(data_model$danceability, na.rm = TRUE),
    max(data_model$danceability, na.rm = TRUE),
    length.out = 200
  ),
  energy = mean(data_model$energy, na.rm = TRUE),
  acousticness = mean(data_model$acousticness, na.rm = TRUE),
  tempo = mean(data_model$tempo, na.rm = TRUE),
  valence = mean(data_model$valence, na.rm = TRUE)
)

prob_dance <- predict(
  model_multinom,
  newdata = grid_dance,
  type = "probs"
)

plot_dance <- grid_dance %>%
  dplyr::bind_cols(as.data.frame(prob_dance)) %>%
  tidyr::pivot_longer(
    cols = c("classical", "pop", "rock", "hip-hop", "r-n-b"),
    names_to = "Genre",
    values_to = "Probabilitas"
  )

ggplot(plot_dance, aes(x = danceability, y = Probabilitas, color = Genre)) +
  geom_line(linewidth = 1.25) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
  scale_color_manual(values = c(
    "classical" = "#5568b8",
    "pop" = "#d18b2f",
    "rock" = "#b84f5a",
    "hip-hop" = "#2f7f73",
    "r-n-b" = "#7c3aed"
  )) +
  labs(
    title = "Prediksi Probabilitas Genre Berdasarkan Danceability",
    subtitle = "Variabel lain dibuat tetap pada nilai rata-rata.",
    x = "Danceability",
    y = "Probabilitas Prediksi",
    color = "Genre"
  ) +
  theme_profesional()

2.19 Visualisasi Pengaruh Energy terhadap Probabilitas Genre

grid_energy <- tibble(
  danceability = mean(data_model$danceability, na.rm = TRUE),
  energy = seq(
    min(data_model$energy, na.rm = TRUE),
    max(data_model$energy, na.rm = TRUE),
    length.out = 200
  ),
  acousticness = mean(data_model$acousticness, na.rm = TRUE),
  tempo = mean(data_model$tempo, na.rm = TRUE),
  valence = mean(data_model$valence, na.rm = TRUE)
)

prob_energy <- predict(
  model_multinom,
  newdata = grid_energy,
  type = "probs"
)

plot_energy <- grid_energy %>%
  dplyr::bind_cols(as.data.frame(prob_energy)) %>%
  tidyr::pivot_longer(
    cols = c("classical", "pop", "rock", "hip-hop", "r-n-b"),
    names_to = "Genre",
    values_to = "Probabilitas"
  )

ggplot(plot_energy, aes(x = energy, y = Probabilitas, color = Genre)) +
  geom_line(linewidth = 1.25) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
  scale_color_manual(values = c(
    "classical" = "#5568b8",
    "pop" = "#d18b2f",
    "rock" = "#b84f5a",
    "hip-hop" = "#2f7f73",
    "r-n-b" = "#7c3aed"
  )) +
  labs(
    title = "Prediksi Probabilitas Genre Berdasarkan Energy",
    subtitle = "Variabel lain dibuat tetap pada nilai rata-rata.",
    x = "Energy",
    y = "Probabilitas Prediksi",
    color = "Genre"
  ) +
  theme_profesional()

2.20 Visualisasi Pengaruh Acousticness terhadap Probabilitas Genre

grid_acoustic <- tibble(
  danceability = mean(data_model$danceability, na.rm = TRUE),
  energy = mean(data_model$energy, na.rm = TRUE),
  acousticness = seq(
    min(data_model$acousticness, na.rm = TRUE),
    max(data_model$acousticness, na.rm = TRUE),
    length.out = 200
  ),
  tempo = mean(data_model$tempo, na.rm = TRUE),
  valence = mean(data_model$valence, na.rm = TRUE)
)

prob_acoustic <- predict(
  model_multinom,
  newdata = grid_acoustic,
  type = "probs"
)

plot_acoustic <- grid_acoustic %>%
  dplyr::bind_cols(as.data.frame(prob_acoustic)) %>%
  tidyr::pivot_longer(
    cols = c("classical", "pop", "rock", "hip-hop", "r-n-b"),
    names_to = "Genre",
    values_to = "Probabilitas"
  )

ggplot(plot_acoustic, aes(x = acousticness, y = Probabilitas, color = Genre)) +
  geom_line(linewidth = 1.25) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
  scale_color_manual(values = c(
    "classical" = "#5568b8",
    "pop" = "#d18b2f",
    "rock" = "#b84f5a",
    "hip-hop" = "#2f7f73",
    "r-n-b" = "#7c3aed"
  )) +
  labs(
    title = "Prediksi Probabilitas Genre Berdasarkan Acousticness",
    subtitle = "Variabel lain dibuat tetap pada nilai rata-rata.",
    x = "Acousticness",
    y = "Probabilitas Prediksi",
    color = "Genre"
  ) +
  theme_profesional()

2.21 Kesimpulan

Berdasarkan analisis regresi logistik multinomial, genre lagu Spotify dapat dimodelkan menggunakan karakteristik audio berupa danceability, energy, acousticness, tempo, dan valence. Model ini sesuai karena variabel respon terdiri dari lima kategori nominal, yaitu pop, rock, hip-hop, classical, dan r-n-b.

Hasil model dapat digunakan untuk melihat karakteristik audio yang membedakan genre lagu serta memprediksi probabilitas sebuah lagu termasuk ke dalam masing-masing genre.

3 BAB 3. Regresi Logistik Ordinal

Bab ini membahas pemodelan respons bertingkat kualitas kopi menggunakan regresi logistik ordinal.

3.1 Pendahuluan

Analisis ini menggunakan regresi logistik ordinal untuk memodelkan variabel respons kualitas_kopi. Variabel respons ini bersifat bertingkat, sehingga kategori kualitas kopi tidak diperlakukan sebagai kategori nominal biasa, tetapi sebagai kategori yang memiliki urutan alami.

Tujuan analisis ini adalah mengetahui pengaruh Altitude_clean, Processing Method, Variety, Moisture Percentage, dan Country of Origin terhadap kecenderungan kopi berada pada kategori kualitas yang lebih tinggi.

3.2 Dasar Teori

Regresi logistik ordinal digunakan ketika variabel respons memiliki lebih dari dua kategori dan kategori tersebut memiliki urutan. Model yang digunakan adalah proportional odds model atau cumulative logit model.

Jika respons ordinal memiliki kategori:

\[ Y_1 < Y_2 < \cdots < Y_J \]

maka model cumulative logit dapat ditulis sebagai:

\[ \log\left(\frac{P(Y \leq j)}{P(Y > j)}\right) = \alpha_j - \beta_1X_1 - \beta_2X_2 - \cdots - \beta_pX_p \]

untuk:

\[ j = 1,2,\ldots,J-1 \]

Keterangan:

• \(Y\) adalah variabel respons ordinal.
• \(\alpha_j\) adalah cutpoint atau batas antar kategori.
• \(\beta\) adalah koefisien prediktor.
• \(X\) adalah variabel prediktor.

Pada fungsi MASS::polr(), koefisien positif menunjukkan bahwa kenaikan nilai prediktor cenderung meningkatkan peluang respon berada pada kategori ordinal yang lebih tinggi, dengan asumsi prediktor lain konstan.

3.3 Import Data

data_raw <- readxl::read_excel("~/Salwa/SMT 4/ADK 1/tugas regresi logostik/data_kopi_ordinal.xlsx")

dim_data <- dim(data_raw)
kolom_data <- names(data_raw)

data_raw %>%
  head(10) %>%
  kable(caption = "Preview Data Awal") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = TRUE
  )

Preview Data Awal

kualitas_kopi	Altitude_clean	Processing Method	Variety	Moisture Percentage	Country of Origin
Tinggi	1700	Double Anaerobic Washed	Castillo	11.8	Colombia
Tinggi	1200	Washed / Wet	Gesha	10.5	Taiwan
Tinggi	1300	Semi Washed	Java	10.4	Laos
Tinggi	1900	Washed / Wet	Gesha	11.8	Costa Rica
Tinggi	1850	Honey,Mossto	Red Bourbon	11.6	Colombia
Tinggi	1668	Washed / Wet	Gesha	10.7	Guatemala
Tinggi	1250	Washed / Wet	Gesha	9.1	Taiwan
Tinggi	1200	Natural / Dry	Sl34+Gesha	10.0	Taiwan
Tinggi	1250	Washed / Wet	SL34	10.8	Taiwan
Tinggi	1400	Washed / Wet	Bourbon	11.0	Tanzania, United Republic Of

cat("Jumlah baris data:", nrow(data_raw), "\n")

## Jumlah baris data: 197

cat("Jumlah kolom data:", ncol(data_raw), "\n")

## Jumlah kolom data: 6

cat("Nama kolom:\n")

## Nama kolom:

print(names(data_raw))

## [1] "kualitas_kopi"       "Altitude_clean"      "Processing Method"  
## [4] "Variety"             "Moisture Percentage" "Country of Origin"

3.4 Validasi Data

Sebelum model dijalankan, data perlu dicek terlebih dahulu. Regresi logistik ordinal dengan MASS::polr() membutuhkan:

1. data memiliki baris observasi,
2. variabel respons tersedia,
3. variabel prediktor tersedia,
4. respons memiliki minimal 3 kategori setelah data kosong dibuang,
5. prediktor numerik tidak bernilai NA, NaN, atau Inf.

required_cols <- c(
  "kualitas_kopi",
  "Altitude_clean",
  "Processing Method",
  "Variety",
  "Moisture Percentage",
  "Country of Origin"
)

missing_cols <- setdiff(required_cols, names(data_raw))

validasi <- tibble::tibble(
  Pemeriksaan = c(
    "Jumlah baris lebih dari 0",
    "Semua kolom yang dibutuhkan tersedia",
    "Respons kualitas_kopi tersedia",
    "Respons memiliki minimal 3 kategori"
  ),
  Status = c(
    nrow(data_raw) > 0,
    length(missing_cols) == 0,
    "kualitas_kopi" %in% names(data_raw),
    ifelse("kualitas_kopi" %in% names(data_raw),
           dplyr::n_distinct(data_raw$kualitas_kopi, na.rm = TRUE) >= 3,
           FALSE)
  )
)

validasi %>%
  dplyr::mutate(Status = ifelse(Status, "Memenuhi", "Belum memenuhi")) %>%
  kable(caption = "Validasi Awal Data") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Validasi Awal Data

Pemeriksaan	Status
Jumlah baris lebih dari 0	Memenuhi
Semua kolom yang dibutuhkan tersedia	Memenuhi
Respons kualitas_kopi tersedia	Memenuhi
Respons memiliki minimal 3 kategori	Memenuhi

if (length(missing_cols) > 0) {
  cat("Kolom yang belum ada dalam file Excel:\n")
  print(missing_cols)
}

kategori_respons <- if ("kualitas_kopi" %in% names(data_raw)) {
  dplyr::n_distinct(data_raw$kualitas_kopi, na.rm = TRUE)
} else {
  0
}

bisa_model <- nrow(data_raw) > 0 && length(missing_cols) == 0 && kategori_respons >= 3

if (!bisa_model) {
  cat('<div class="danger-box">
')
  cat('<b>Model belum bisa dijalankan.</b><br>')
  cat('File Excel yang terbaca belum memenuhi syarat untuk regresi logistik ordinal. ')
  cat('Pastikan file Excel berisi data observasi lengkap, bukan hanya header kolom. ')
  cat('Minimal harus ada 3 kategori pada variabel <code>kualitas_kopi</code>.')
  cat('</div>')
}

3.5 Persiapan Data

if (bisa_model) {
  data_kopi <- data_raw %>%
    dplyr::select(
      kualitas_kopi,
      Altitude_clean,
      `Processing Method`,
      Variety,
      `Moisture Percentage`,
      `Country of Origin`
    ) %>%
    dplyr::mutate(
      kualitas_kopi = as.character(kualitas_kopi),
      kualitas_kopi = stringr::str_squish(kualitas_kopi),
      kualitas_kopi = factor(kualitas_kopi, ordered = TRUE),
      Altitude_clean = as.numeric(Altitude_clean),
      `Moisture Percentage` = as.numeric(`Moisture Percentage`),
      `Processing Method` = as.factor(`Processing Method`),
      Variety = as.factor(Variety),
      `Country of Origin` = as.factor(`Country of Origin`)
    ) %>%
    dplyr::filter(
      !is.na(kualitas_kopi),
      !is.na(Altitude_clean),
      !is.na(`Moisture Percentage`),
      !is.na(`Processing Method`),
      !is.na(Variety),
      !is.na(`Country of Origin`)
    )

  data_kopi <- data_kopi %>%
    dplyr::mutate(
      `Processing Method` = forcats::fct_lump_min(`Processing Method`, min = 5, other_level = "Other"),
      Variety = forcats::fct_lump_min(Variety, min = 5, other_level = "Other"),
      `Country of Origin` = forcats::fct_lump_min(`Country of Origin`, min = 5, other_level = "Other"),
      Altitude_z = as.numeric(scale(Altitude_clean)),
      Moisture_z = as.numeric(scale(`Moisture Percentage`))
    ) %>%
    dplyr::filter(is.finite(Altitude_z), is.finite(Moisture_z))

  data_kopi$kualitas_kopi <- droplevels(data_kopi$kualitas_kopi)

  data_kopi %>%
    glimpse()
}

## Rows: 197
## Columns: 8
## $ kualitas_kopi         <ord> Tinggi, Tinggi, Tinggi, Tinggi, Tinggi, Tinggi, …
## $ Altitude_clean        <dbl> 1700, 1200, 1300, 1900, 1850, 1668, 1250, 1200, …
## $ `Processing Method`   <fct> Other, Washed / Wet, Other, Washed / Wet, Other,…
## $ Variety               <fct> Other, Gesha, Other, Gesha, Other, Gesha, Gesha,…
## $ `Moisture Percentage` <dbl> 11.8, 10.5, 10.4, 11.8, 11.6, 10.7, 9.1, 10.0, 1…
## $ `Country of Origin`   <fct> "Colombia", "Taiwan", "Other", "Costa Rica", "Co…
## $ Altitude_z            <dbl> 0.56444850, -0.12640089, 0.01176899, 0.84078826,…
## $ Moisture_z            <dbl> 0.8752429267, -0.1581812943, -0.2376754651, 0.87…

Variabel numerik Altitude_clean dan Moisture Percentage distandarisasi menjadi Altitude_z dan Moisture_z. Tujuannya adalah membuat proses estimasi lebih stabil, terutama jika skala antar variabel sangat berbeda.

3.6 Distribusi Respons Ordinal

data_kopi %>%
  dplyr::count(kualitas_kopi) %>%
  dplyr::mutate(proporsi = n / sum(n)) %>%
  kable(
    digits = 3,
    caption = "Distribusi Kategori Kualitas Kopi"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Distribusi Kategori Kualitas Kopi

kualitas_kopi	n	proporsi
Rendah	3	0.015
Sedang	149	0.756
Tinggi	45	0.228

data_kopi %>%
  dplyr::count(kualitas_kopi) %>%
  dplyr::mutate(proporsi = n / sum(n)) %>%
  ggplot(aes(x = kualitas_kopi, y = proporsi, fill = kualitas_kopi)) +
  geom_col(width = 0.65, alpha = 0.92) +
  geom_text(aes(label = scales::percent(proporsi, accuracy = 0.1)),
            vjust = -0.35, fontface = "bold") +
  scale_y_continuous(labels = scales::percent_format(), limits = c(0, NA)) +
  labs(
    title = "Distribusi Kualitas Kopi",
    subtitle = "Respons ordinal: kategori kualitas memiliki urutan alami.",
    x = "Kualitas Kopi",
    y = "Proporsi"
  ) +
  theme_minimal(base_size = 13) +
  theme(legend.position = "none")

Distribusi respons digunakan untuk melihat apakah kategori kualitas kopi cukup terwakili. Jika salah satu kategori memiliki jumlah data sangat sedikit, model ordinal dapat menjadi tidak stabil.

3.7 Estimasi Model Regresi Logistik Ordinal

fit_sukses <- FALSE
if (bisa_model) {
  n_level_model <- nlevels(data_kopi$kualitas_kopi)
  cat("Jumlah level respons setelah cleaning:", n_level_model, "\n")

  bisa_fit <- nrow(data_kopi) > 0 && n_level_model >= 3
} else {
  bisa_fit <- FALSE
}

## Jumlah level respons setelah cleaning: 3

## Model utama dengan semua prediktor sesuai permintaan.
## tryCatch dipakai supaya proses knit tidak berhenti jika polr gagal optimasi.
fit_formula_full <- kualitas_kopi ~ Altitude_z + `Processing Method` + Variety + Moisture_z + `Country of Origin`

fit_ord <- tryCatch(
  MASS::polr(
    formula = fit_formula_full,
    data = data_kopi,
    method = "logistic",
    Hess = TRUE,
    control = list(maxit = 1000)
  ),
  error = function(e) e
)

fit_sukses <- inherits(fit_ord, "polr")

if (fit_sukses) {
  summary(fit_ord)
} else {
  cat("Model utama gagal dijalankan. Pesan error dari polr:\n")
  cat(conditionMessage(fit_ord), "\n")
}

## Model utama gagal dijalankan. Pesan error dari polr:
## initial value in 'vmmin' is not finite

## Jika model penuh gagal, jalankan model cadangan yang lebih stabil
## dengan prediktor numerik saja. Ini menjaga dokumen tetap bisa di-knit.
if (!exists("fit_sukses")) fit_sukses <- FALSE

if (!fit_sukses) {
  fit_formula_simple <- kualitas_kopi ~ Altitude_z + Moisture_z

  fit_ord_simple <- tryCatch(
    MASS::polr(
      formula = fit_formula_simple,
      data = data_kopi,
      method = "logistic",
      Hess = TRUE,
      control = list(maxit = 1000)
    ),
    error = function(e) e
  )

  fit_simple_sukses <- inherits(fit_ord_simple, "polr")

  if (fit_simple_sukses) {
    fit_ord <- fit_ord_simple
    fit_sukses <- TRUE
    model_yang_dipakai <- "Model cadangan: kualitas_kopi ~ Altitude_z + Moisture_z"
    summary(fit_ord)
  } else {
    model_yang_dipakai <- "Tidak ada model yang berhasil diestimasi"
    cat("Model cadangan juga gagal dijalankan. Pesan error:\n")
    cat(conditionMessage(fit_ord_simple), "\n")
  }
} else {
  model_yang_dipakai <- "Model utama: semua prediktor"
}

## Call:
## MASS::polr(formula = fit_formula_simple, data = data_kopi, control = list(maxit = 1000), 
##     Hess = TRUE, method = "logistic")
## 
## Coefficients:
##              Value Std. Error t value
## Altitude_z  0.1829     0.1564  1.1697
## Moisture_z -0.1323     0.1587 -0.8335
## 
## Intercepts:
##               Value   Std. Error t value
## Rendah|Sedang -4.1884  0.5823    -7.1924
## Sedang|Tinggi  1.2289  0.1712     7.1788
## 
## Residual Deviance: 239.4656 
## AIC: 247.4656

if (exists("bisa_fit") && !bisa_fit) {
  cat('<div class="danger-box">
')
  cat('<b>Model tidak dijalankan.</b><br>')
  cat('Setelah proses cleaning, variabel <code>kualitas_kopi</code> memiliki kurang dari 3 kategori atau data observasi tidak tersedia. ')
  cat('Regresi logistik ordinal dengan <code>MASS::polr()</code> membutuhkan minimal 3 level respons.')
  cat('</div>')
}

if (exists("bisa_fit") && bisa_fit && exists("fit_sukses") && !fit_sukses) {
  cat('<div class="danger-box">
')
  cat('<b>Model gagal diestimasi walaupun data lolos validasi awal.</b><br>')
  cat('Penyebab paling umum adalah kategori respons atau prediktor yang terlalu jarang, separation, atau kombinasi kategori yang terlalu banyak. ')
  cat('Solusi analisis: gabungkan kategori yang sangat jarang, kurangi jumlah level pada prediktor kategorik, atau gunakan model ordinal alternatif.')
  cat('</div>')
}

if (exists("fit_sukses") && fit_sukses) {
  cat('<div class="interpretasi">
')
  cat('<b>Model berhasil diestimasi.</b><br>')
  cat('Model yang dipakai: ', model_yang_dipakai, '.')
  cat('</div>')
}

Model berhasil diestimasi.
Model yang dipakai: Model cadangan: kualitas_kopi ~ Altitude_z + Moisture_z .

3.8 Ringkasan Koefisien dan Odds Ratio

ord_coef <- coef(summary(fit_ord))

result_ord <- as.data.frame(ord_coef) %>%
  tibble::rownames_to_column("parameter") %>%
  dplyr::rename(
    estimate = Value,
    std_error = `Std. Error`,
    t_value = `t value`
  ) %>%
  dplyr::mutate(
    p_value = 2 * (1 - pnorm(abs(t_value))),
    jenis = ifelse(grepl("\\|", parameter), "Cutpoint", "Koefisien"),
    OR = ifelse(jenis == "Koefisien", exp(estimate), NA_real_),
    CI_low = ifelse(jenis == "Koefisien", exp(estimate - 1.96 * std_error), NA_real_),
    CI_high = ifelse(jenis == "Koefisien", exp(estimate + 1.96 * std_error), NA_real_)
  )

result_ord %>%
  dplyr::mutate(
    dplyr::across(c(estimate, std_error, t_value, p_value, OR, CI_low, CI_high), ~ round(.x, 4))
  ) %>%
  kable(
    caption = "Ringkasan Hasil Regresi Logistik Ordinal",
    col.names = c("Parameter", "Estimate", "SE", "z/t-value", "p-value", "Jenis", "Odds Ratio", "CI 2.5%", "CI 97.5%")
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = TRUE
  )

Ringkasan Hasil Regresi Logistik Ordinal

Parameter	Estimate	SE	z/t-value	p-value	Jenis	Odds Ratio	CI 2.5%	CI 97.5%
Altitude_z	0.1829	0.1564	1.1697	0.2421	Koefisien	1.2007	0.8838	1.6313
Moisture_z	-0.1323	0.1587	-0.8335	0.4046	Koefisien	0.8761	0.6419	1.1958
Rendah|Sedang	-4.1884	0.5823	-7.1924	0.0000	Cutpoint	NA	NA	NA
Sedang|Tinggi	1.2289	0.1712	7.1788	0.0000	Cutpoint	NA	NA	NA

Interpretasi umum: jika koefisien bernilai positif, maka kenaikan prediktor tersebut cenderung meningkatkan peluang kopi berada pada kategori kualitas yang lebih tinggi. Jika koefisien bernilai negatif, maka kenaikan prediktor tersebut cenderung menurunkan peluang kopi berada pada kategori kualitas yang lebih tinggi.

3.9 Interpretasi Otomatis Variabel Signifikan

interpretasi <- result_ord %>%
  dplyr::filter(jenis == "Koefisien") %>%
  dplyr::mutate(
    keputusan = ifelse(p_value < 0.05, "signifikan", "tidak signifikan"),
    arah = ifelse(estimate > 0, "meningkatkan", "menurunkan"),
    perubahan_odds = round((OR - 1) * 100, 2)
  )

cat('<div class="interpretasi">
')

cat('<b>Interpretasi:</b><br>')

Interpretasi:

for (i in seq_len(nrow(interpretasi))) {
  cat('- Variabel <b>', interpretasi$parameter[i], '</b> ', interpretasi$keputusan[i],
      ' pada taraf 5%. Nilai koefisiennya ', round(interpretasi$estimate[i], 4),
      ', sehingga variabel ini cenderung ', interpretasi$arah[i],
      ' peluang kopi berada pada kategori kualitas yang lebih tinggi. ', sep = '')
  cat('Odds ratio = ', round(interpretasi$OR[i], 4), '.', '<br>', sep = '')
}

• Variabel Altitude_z tidak signifikan pada taraf 5%. Nilai koefisiennya 0.1829, sehingga variabel ini cenderung meningkatkan peluang kopi berada pada kategori kualitas yang lebih tinggi. Odds ratio = 1.2007.
  - Variabel Moisture_z tidak signifikan pada taraf 5%. Nilai koefisiennya -0.1323, sehingga variabel ini cenderung menurunkan peluang kopi berada pada kategori kualitas yang lebih tinggi. Odds ratio = 0.8761.
  

cat('</div>')

3.10 Prediksi Probabilitas

pred_prob_ord <- predict(fit_ord, type = "probs")

head(pred_prob_ord) %>%
  as.data.frame() %>%
  kable(digits = 4, caption = "Contoh Prediksi Probabilitas Tiap Kategori") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = TRUE
  )

Contoh Prediksi Probabilitas Tiap Kategori

Rendah	Sedang	Tinggi
0.0151	0.7607	0.2242
0.0150	0.7590	0.2260
0.0145	0.7532	0.2323
0.0144	0.7525	0.2331
0.0143	0.7511	0.2346
0.0136	0.7429	0.2434

data_kopi_pred <- data_kopi %>%
  dplyr::bind_cols(as.data.frame(pred_prob_ord)) %>%
  dplyr::mutate(prediksi = predict(fit_ord, type = "class"))

head(data_kopi_pred, 10) %>%
  kable(caption = "Data dengan Hasil Prediksi") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = TRUE
  )

Data dengan Hasil Prediksi

kualitas_kopi	Altitude_clean	Processing Method	Variety	Moisture Percentage	Country of Origin	Altitude_z	Moisture_z	Rendah	Sedang	Tinggi	prediksi
Tinggi	1700	Other	Other	11.8	Colombia	0.5644485	0.8752429	0.0151295	0.7606878	0.2241827	Sedang
Tinggi	1200	Washed / Wet	Gesha	10.5	Taiwan	-0.1264009	-0.1581813	0.0149762	0.7590376	0.2259862	Sedang
Tinggi	1300	Other	Other	10.4	Other	0.0117690	-0.2376755	0.0144573	0.7532355	0.2323072	Sedang
Tinggi	1900	Washed / Wet	Gesha	11.8	Costa Rica	0.8407883	0.8752429	0.0143945	0.7525096	0.2330958	Sedang
Tinggi	1850	Other	Other	11.6	Colombia	0.7717033	0.7162546	0.0142759	0.7511243	0.2345998	Sedang
Tinggi	1668	Washed / Wet	Gesha	10.7	Guatemala	0.5202341	0.0008070	0.0136073	0.7429460	0.2434467	Sedang
Tinggi	1250	Washed / Wet	Gesha	9.1	Taiwan	-0.0573160	-1.2710997	0.0127921	0.7320457	0.2551622	Sedang
Tinggi	1200	Natural / Dry	Other	10.0	Taiwan	-0.1264009	-0.5556521	0.0142200	0.7504651	0.2353149	Sedang
Tinggi	1250	Washed / Wet	SL34	10.8	Taiwan	-0.0573160	0.0803012	0.0152577	0.7620466	0.2226957	Sedang
Tinggi	1400	Washed / Wet	Bourbon	11.0	Tanzania, United Republic Of	0.1499389	0.2392896	0.0150062	0.7593630	0.2256308	Sedang

3.11 Confusion Matrix dan Akurasi

conf_ord <- table(
  Aktual = data_kopi_pred$kualitas_kopi,
  Prediksi = data_kopi_pred$prediksi
)

conf_ord

##         Prediksi
## Aktual   Rendah Sedang Tinggi
##   Rendah      0      3      0
##   Sedang      0    148      1
##   Tinggi      0     45      0

accuracy_ord <- sum(diag(conf_ord)) / sum(conf_ord)
accuracy_ord

## [1] 0.751269

Akurasi pada respons ordinal perlu dibaca hati-hati. Kesalahan prediksi antar kategori yang berdekatan lebih ringan dibanding kesalahan prediksi yang lompat jauh. Karena itu, confusion matrix hanya menjadi evaluasi awal, bukan satu-satunya ukuran kualitas model.

3.12 Visualisasi Prediksi Probabilitas

if ("Altitude_z" %in% names(data_kopi)) {
  grid_ord <- data.frame(
    Altitude_z = seq(min(data_kopi$Altitude_z), max(data_kopi$Altitude_z), length.out = 120),
    Moisture_z = mean(data_kopi$Moisture_z),
    `Processing Method` = names(sort(table(data_kopi$`Processing Method`), decreasing = TRUE))[1],
    Variety = names(sort(table(data_kopi$Variety), decreasing = TRUE))[1],
    `Country of Origin` = names(sort(table(data_kopi$`Country of Origin`), decreasing = TRUE))[1],
    check.names = FALSE
  )

  grid_prob_ord <- predict(fit_ord, newdata = grid_ord, type = "probs")

  grid_ord_plot <- grid_ord %>%
    dplyr::bind_cols(as.data.frame(grid_prob_ord)) %>%
    pivot_longer(
      cols = dplyr::all_of(levels(data_kopi$kualitas_kopi)),
      names_to = "kualitas_kopi",
      values_to = "probabilitas"
    )

  ggplot(grid_ord_plot, aes(x = Altitude_z, y = probabilitas, color = kualitas_kopi)) +
    geom_line(linewidth = 1.2) +
    scale_y_continuous(labels = scales::percent_format()) +
    labs(
      title = "Prediksi Probabilitas Kualitas Kopi Berdasarkan Altitude",
      subtitle = "Variabel lain ditahan pada nilai rata-rata atau kategori terbanyak.",
      x = "Altitude_clean terstandarisasi",
      y = "Probabilitas prediksi",
      color = "Kualitas Kopi"
    ) +
    theme_minimal(base_size = 13)
}

Grafik probabilitas menunjukkan perubahan peluang setiap kategori kualitas kopi ketika altitude berubah. Pada model ordinal, yang diasumsikan sejajar adalah cumulative logit, bukan kurva probabilitas kategori tunggal.

3.13 Pemeriksaan Asumsi Proportional Odds

Asumsi utama regresi logistik ordinal adalah proportional odds atau parallel lines assumption. Artinya, pengaruh prediktor dianggap sama pada setiap batas kategori kumulatif.

Secara matematis:

\[ \log\left(\frac{P(Y \leq j)}{P(Y > j)}\right) = \alpha_j - \mathbf{x}^T\boldsymbol{\beta} \]

Slope \(\beta\) sama untuk semua batas \(j\), sedangkan cutpoint \(\alpha_j\) berbeda antar batas kategori.

if (requireNamespace("ordinal", quietly = TRUE)) {
  library(ordinal)

  fit_clm <- ordinal::clm(
    kualitas_kopi ~ Altitude_z + `Processing Method` + Variety + Moisture_z + `Country of Origin`,
    data = data_kopi,
    link = "logit"
  )

  ordinal::nominal_test(fit_clm)
} else {
  cat("Paket 'ordinal' belum terpasang. Bagian nominal_test dilewati agar proses knit tetap berjalan.")
}

Jika hasil nominal_test() memiliki p-value kecil pada suatu variabel, hal tersebut dapat menjadi indikasi bahwa asumsi proportional odds untuk variabel tersebut kurang terpenuhi. Namun keputusan akhir tetap perlu mempertimbangkan teori, ukuran sampel, dan tujuan analisis.

3.14 Kesimpulan

if (exists("fit_sukses") && fit_sukses) {
  cat('<div class="interpretasi">
')
  cat('Berdasarkan model regresi logistik ordinal, variabel Altitude_clean, Processing Method, Variety, Moisture Percentage, dan Country of Origin dapat digunakan untuk menjelaskan kecenderungan kategori kualitas kopi. Variabel dengan p-value kurang dari 0,05 dapat dianggap memiliki pengaruh yang signifikan terhadap peluang kopi berada pada kategori kualitas yang lebih tinggi.')
  cat('</div>')
} else {
  cat('<div class="danger-box">
')
  cat('Kesimpulan model belum dapat dibuat karena data pada file Excel belum memenuhi syarat analisis. File yang digunakan harus memuat observasi lengkap dan variabel kualitas_kopi harus memiliki minimal 3 kategori.')
  cat('</div>')
}

Berdasarkan model regresi logistik ordinal, variabel Altitude_clean, Processing Method, Variety, Moisture Percentage, dan Country of Origin dapat digunakan untuk menjelaskan kecenderungan kategori kualitas kopi. Variabel dengan p-value kurang dari 0,05 dapat dianggap memiliki pengaruh yang signifikan terhadap peluang kopi berada pada kategori kualitas yang lebih tinggi.

4 BAB 4. Regresi Poisson

Bab ini membahas pemodelan data hitungan win_count pada candy menggunakan regresi Poisson.

4.1 Pendahuluan

Analisis ini bertujuan untuk memodelkan jumlah kemenangan candy dalam head-to-head matchup menggunakan regresi Poisson. Variabel respon yang digunakan adalah win_count, sedangkan variabel prediktor yang digunakan adalah:

1. chocolate
2. caramel
3. peanutyalmondy
4. sugarpercent
5. pricepercent

Catatan

Tujuan analisis: mengetahui apakah kandungan cokelat, karamel, kacang/almond, tingkat gula, dan tingkat harga berpengaruh terhadap jumlah kemenangan candy. Karena variabel respon berupa data hitungan nonnegatif, regresi Poisson menjadi salah satu model yang sesuai untuk digunakan.

4.2 Dasar Teori Regresi Poisson

Regresi Poisson digunakan ketika variabel respon berbentuk data hitung, yaitu bilangan bulat nonnegatif seperti 0, 1, 2, 3, dan seterusnya. Pada penelitian ini, win_count menunjukkan banyaknya kemenangan suatu candy, sehingga secara konsep dapat dimodelkan sebagai variabel count.

Misalkan:

\[ Y_i = \text{jumlah kemenangan candy ke-}i \]

Dalam model Poisson, diasumsikan:

\[ Y_i \sim \text{Poisson}(\mu_i) \]

dengan:

\[ E(Y_i \mid X_i) = \mu_i \]

Model regresi Poisson menggunakan fungsi penghubung log, sehingga:

\[ \log(\mu_i) = \beta_0 + \beta_1X_{1i}+\beta_2X_{2i}+\cdots+\beta_pX_{pi} \]

atau setara dengan:

\[ \mu_i = \exp(\beta_0 + \beta_1X_{1i}+\beta_2X_{2i}+\cdots+\beta_pX_{pi}) \]

Untuk penelitian ini, model yang digunakan adalah:

\[ \log(\mu_i) = \beta_0 + \beta_1\text{chocolate}_i + \beta_2\text{caramel}_i + \beta_3\text{peanutyalmondy}_i + \beta_4\text{sugarpercent}_i + \beta_5\text{pricepercent}_i \]

Interpretasi

Interpretasi utama: pada regresi Poisson, nilai exp(koefisien) disebut Incidence Rate Ratio (IRR). Jika IRR > 1, maka prediktor meningkatkan rata-rata jumlah kemenangan. Jika IRR < 1, maka prediktor menurunkan rata-rata jumlah kemenangan, dengan asumsi prediktor lain konstan.

4.3 Import Data

data_raw <- readxl::read_excel("~/Salwa/SMT 4/ADK 1/tugas regresi logostik/data_candy_poissonn.xlsx")

glimpse(data_raw)

## Rows: 85
## Columns: 6
## $ win_count      <dbl> 67, 68, 32, 46, 52, 50, 57, 23, 38, 35, 39, 36, 25, 42,…
## $ chocolate      <dbl> 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0…
## $ caramel        <dbl> 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ peanutyalmondy <dbl> 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ sugarpercent   <dbl> 0.732, 0.604, 0.011, 0.011, 0.906, 0.465, 0.604, 0.313,…
## $ pricepercent   <dbl> 0.860, 0.511, 0.116, 0.511, 0.511, 0.767, 0.767, 0.511,…

4.4 Pemeriksaan Awal Data

required_cols <- c(
  "win_count",
  "chocolate",
  "caramel",
  "peanutyalmondy",
  "sugarpercent",
  "pricepercent"
)

missing_cols <- setdiff(required_cols, names(data_raw))

if(length(missing_cols) > 0){
  stop(paste("Kolom berikut tidak ditemukan di dataset:", paste(missing_cols, collapse = ", ")))
}

cat("Semua kolom yang dibutuhkan tersedia di dataset.")

## Semua kolom yang dibutuhkan tersedia di dataset.

data_candy <- data_raw %>%
  dplyr::select(dplyr::all_of(required_cols)) %>%
  dplyr::mutate(
    win_count = as.numeric(win_count),
    chocolate = factor(chocolate, levels = c(0, 1), labels = c("Tidak", "Ya")),
    caramel = factor(caramel, levels = c(0, 1), labels = c("Tidak", "Ya")),
    peanutyalmondy = factor(peanutyalmondy, levels = c(0, 1), labels = c("Tidak", "Ya")),
    sugarpercent = as.numeric(sugarpercent),
    pricepercent = as.numeric(pricepercent)
  ) %>%
  dplyr::filter(
    !is.na(win_count),
    !is.na(chocolate),
    !is.na(caramel),
    !is.na(peanutyalmondy),
    !is.na(sugarpercent),
    !is.na(pricepercent)
  ) %>%
  dplyr::filter(win_count >= 0) %>%
  dplyr::mutate(
    win_count = round(win_count),
    sugarpercent_10 = sugarpercent / 0.10,
    pricepercent_10 = pricepercent / 0.10
  ) %>%
  droplevels()

summary(data_candy)

##    win_count     chocolate   caramel   peanutyalmondy  sugarpercent   
##  Min.   :22.00   Tidak:48   Tidak:71   Tidak:71       Min.   :0.0110  
##  1st Qu.:39.00   Ya   :37   Ya   :14   Ya   :14       1st Qu.:0.2200  
##  Median :48.00                                        Median :0.4650  
##  Mean   :50.26                                        Mean   :0.4786  
##  3rd Qu.:60.00                                        3rd Qu.:0.7320  
##  Max.   :84.00                                        Max.   :0.9880  
##   pricepercent    sugarpercent_10 pricepercent_10
##  Min.   :0.0110   Min.   :0.110   Min.   :0.110  
##  1st Qu.:0.2550   1st Qu.:2.200   1st Qu.:2.550  
##  Median :0.4650   Median :4.650   Median :4.650  
##  Mean   :0.4689   Mean   :4.786   Mean   :4.689  
##  3rd Qu.:0.6510   3rd Qu.:7.320   3rd Qu.:6.510  
##  Max.   :0.9760   Max.   :9.880   Max.   :9.760

Catatan

Catatan cleaning: variabel biner chocolate, caramel, dan peanutyalmondy diperlakukan sebagai faktor dengan kategori referensi Tidak. Variabel sugarpercent dan pricepercent tetap digunakan dalam bentuk proporsi. Untuk interpretasi tambahan, dibuat juga versi per 10 percentage point, yaitu sugarpercent_10 dan pricepercent_10.

4.5 Eksplorasi Data

4.5.1 Ringkasan Deskriptif

data_candy %>%
  dplyr::summarise(
    n = n(),
    mean_win_count = mean(win_count),
    median_win_count = median(win_count),
    min_win_count = min(win_count),
    max_win_count = max(win_count),
    mean_sugarpercent = mean(sugarpercent),
    mean_pricepercent = mean(pricepercent)
  ) %>%
  dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3))) %>%
  kable(
    caption = "Ringkasan deskriptif variabel numerik"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Ringkasan deskriptif variabel numerik

n	mean_win_count	median_win_count	min_win_count	max_win_count	mean_sugarpercent	mean_pricepercent
85	50.259	48	22	84	0.479	0.469

4.5.2 Distribusi Variabel Respon

data_candy %>%
  dplyr::count(win_count) %>%
  dplyr::mutate(proporsi = n / sum(n)) %>%
  dplyr::mutate(proporsi = percent(proporsi, accuracy = 0.1)) %>%
  kable(
    caption = "Distribusi jumlah kemenangan candy"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Distribusi jumlah kemenangan candy

win_count	n	proporsi
22	1	1.2%
23	1	1.2%
25	1	1.2%
27	1	1.2%
28	1	1.2%
30	1	1.2%
32	2	2.4%
33	1	1.2%
34	1	1.2%
35	4	4.7%
36	1	1.2%
37	1	1.2%
38	3	3.5%
39	6	7.1%
41	2	2.4%
42	3	3.5%
43	3	3.5%
44	1	1.2%
45	1	1.2%
46	5	5.9%
47	2	2.4%
48	1	1.2%
49	1	1.2%
50	3	3.5%
51	1	1.2%
52	1	1.2%
53	2	2.4%
55	6	7.1%
56	1	1.2%
57	3	3.5%
59	1	1.2%
60	2	2.4%
61	1	1.2%
62	1	1.2%
63	1	1.2%
64	1	1.2%
66	2	2.4%
67	3	3.5%
68	1	1.2%
69	1	1.2%
71	2	2.4%
73	3	3.5%
77	2	2.4%
82	2	2.4%
84	1	1.2%

ggplot(data_candy, aes(x = win_count)) +
  geom_histogram(binwidth = 5, fill = "#2F7F73", color = "white", alpha = 0.92) +
  labs(
    title = "Distribusi Jumlah Kemenangan Candy",
    subtitle = "Variabel respon win_count berupa data hitungan nonnegatif.",
    x = "Jumlah kemenangan",
    y = "Frekuensi"
  ) +
  theme_candy()

Catatan

Interpretasi: grafik distribusi digunakan untuk melihat pola awal variabel win_count. Karena nilai respon berupa bilangan hitungan, regresi Poisson relevan sebagai pendekatan awal. Namun, tetap perlu diperiksa apakah terjadi overdispersion setelah model diestimasi.

4.5.3 Distribusi Prediktor Biner

data_candy %>%
  dplyr::select(chocolate, caramel, peanutyalmondy) %>%
  pivot_longer(cols = dplyr::everything(), names_to = "prediktor", values_to = "kategori") %>%
  dplyr::count(prediktor, kategori) %>%
  group_by(prediktor) %>%
  dplyr::mutate(proporsi = n / sum(n)) %>%
  ungroup() %>%
  dplyr::mutate(proporsi = percent(proporsi, accuracy = 0.1)) %>%
  kable(
    caption = "Distribusi prediktor biner"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Distribusi prediktor biner

prediktor	kategori	n	proporsi
caramel	Tidak	71	83.5%
caramel	Ya	14	16.5%
chocolate	Tidak	48	56.5%
chocolate	Ya	37	43.5%
peanutyalmondy	Tidak	71	83.5%
peanutyalmondy	Ya	14	16.5%

data_candy %>%
  dplyr::select(chocolate, caramel, peanutyalmondy) %>%
  pivot_longer(cols = dplyr::everything(), names_to = "prediktor", values_to = "kategori") %>%
  ggplot(aes(x = kategori, fill = kategori)) +
  geom_bar(alpha = 0.92, width = 0.65) +
  facet_wrap(~ prediktor, scales = "free_y") +
  scale_fill_manual(values = c("Tidak" = "#B84F5A", "Ya" = "#2F7F73")) +
  labs(
    title = "Distribusi Prediktor Biner",
    subtitle = "Kategori referensi pada model adalah Tidak.",
    x = NULL,
    y = "Frekuensi",
    fill = "Kategori"
  ) +
  theme_candy()

4.6 Estimasi Model Regresi Poisson

Model utama yang digunakan adalah:

model_formula <- win_count ~ chocolate + caramel + peanutyalmondy + sugarpercent + pricepercent

fit_pois <- tryCatch(
  glm(
    formula = model_formula,
    data = data_candy,
    family = poisson(link = "log")
  ),
  error = function(e){
    message("Model Poisson gagal diestimasi: ", e$message)
    NULL
  }
)

model_success <- !is.null(fit_pois)

if(model_success){
  summary(fit_pois)
} else {
  cat("Model gagal diestimasi. Periksa kembali struktur data dan distribusi variabel.")
}

## 
## Call:
## glm(formula = model_formula, family = poisson(link = "log"), 
##     data = data_candy)
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       3.67261    0.03730  98.465  < 2e-16 ***
## chocolateYa       0.33956    0.03775   8.995  < 2e-16 ***
## caramelYa         0.02775    0.04172   0.665  0.50585    
## peanutyalmondyYa  0.12761    0.04128   3.092  0.00199 ** 
## sugarpercent      0.17937    0.06006   2.987  0.00282 ** 
## pricepercent     -0.07263    0.06805  -1.067  0.28581    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 364.31  on 84  degrees of freedom
## Residual deviance: 196.83  on 79  degrees of freedom
## AIC: 694.56
## 
## Number of Fisher Scoring iterations: 4

Interpretasi

Model yang diestimasi: model Poisson dengan link log. Artinya, model tidak memprediksi win_count secara langsung dalam bentuk linear, melainkan memprediksi log dari rata-rata jumlah kemenangan. Hasil akhirnya dapat dikembalikan ke skala asli menggunakan fungsi eksponensial.

4.7 Ringkasan Koefisien dan IRR

if(model_success){
  pois_coef <- as.data.frame(coef(summary(fit_pois))) %>%
    tibble::rownames_to_column("parameter") %>%
    dplyr::rename(
      estimate = Estimate,
      std_error = `Std. Error`,
      z_value = `z value`,
      p_value = `Pr(>|z|)`
    ) %>%
    dplyr::mutate(
      IRR = exp(estimate),
      CI_low = exp(estimate - 1.96 * std_error),
      CI_high = exp(estimate + 1.96 * std_error),
      perubahan_persen = 100 * (IRR - 1),
      keputusan = ifelse(p_value < 0.05, "Signifikan", "Tidak signifikan")
    )

  pois_coef %>%
    dplyr::mutate(
      dplyr::across(
        c(estimate, std_error, z_value, p_value, IRR, CI_low, CI_high, perubahan_persen),
        ~ round(.x, 4)
      )
    ) %>%
    kable(
      caption = "Ringkasan hasil regresi Poisson",
      col.names = c(
        "Parameter", "Estimate", "SE", "z-value", "p-value",
        "IRR", "CI 2.5%", "CI 97.5%", "Perubahan (%)", "Keputusan"
      )
    ) %>%
    kable_styling(
      bootstrap_options = c("striped", "hover", "condensed", "responsive"),
      full_width = TRUE
    ) %>%
    row_spec(which(pois_coef$p_value < 0.05), bold = TRUE, color = "#1F5E4D")
} else {
  box_warning("Tabel koefisien tidak ditampilkan karena model gagal diestimasi.")
}

Ringkasan hasil regresi Poisson

Parameter	Estimate	SE	z-value	p-value	IRR	CI 2.5%	CI 97.5%	Perubahan (%)	Keputusan
(Intercept)	3.6726	0.0373	98.4648	0.0000	39.3544	36.5800	42.3392	3835.4398	Signifikan
chocolateYa	0.3396	0.0378	8.9948	0.0000	1.4043	1.3042	1.5122	40.4337	Signifikan
caramelYa	0.0278	0.0417	0.6653	0.5059	1.0281	0.9474	1.1157	2.8143	Tidak signifikan
peanutyalmondyYa	0.1276	0.0413	3.0917	0.0020	1.1361	1.0478	1.2318	13.6113	Signifikan
sugarpercent	0.1794	0.0601	2.9868	0.0028	1.1965	1.0636	1.3459	19.6469	Signifikan
pricepercent	-0.0726	0.0681	-1.0674	0.2858	0.9299	0.8138	1.0626	-7.0059	Tidak signifikan

4.8 Interpretasi Koefisien

Interpretasi

Variabel yang signifikan pada taraf 5%: chocolateYa (IRR = 1.404, p-value = 0); peanutyalmondyYa (IRR = 1.136, p-value = 0.002); sugarpercent (IRR = 1.196, p-value = 0.0028). Artinya, variabel-variabel tersebut memiliki hubungan yang bermakna secara statistik terhadap rata-rata jumlah kemenangan candy, setelah mengendalikan prediktor lain dalam model.

Secara umum, interpretasi koefisien regresi Poisson adalah sebagai berikut:

• Jika koefisien bernilai positif, maka nilai IRR lebih dari 1. Artinya, prediktor tersebut meningkatkan rata-rata jumlah kemenangan.
• Jika koefisien bernilai negatif, maka nilai IRR kurang dari 1. Artinya, prediktor tersebut menurunkan rata-rata jumlah kemenangan.
• Jika p-value kurang dari 0,05, maka prediktor tersebut signifikan secara statistik pada taraf 5%.

4.8.1 Interpretasi Khusus Prediktor

Catatan

Chocolate: Candy dengan kategori Ya memiliki rata-rata jumlah kemenangan sekitar 40.43% lebih tinggi dibanding candy dengan kategori Tidak, dengan asumsi prediktor lain konstan. Hubungan ini signifikan pada taraf 5%.

Caramel: Candy dengan kategori Ya memiliki rata-rata jumlah kemenangan sekitar 2.81% lebih tinggi dibanding candy dengan kategori Tidak, dengan asumsi prediktor lain konstan. Hubungan ini tidak signifikan pada taraf 5%.

Peanutyalmondy: Candy dengan kategori Ya memiliki rata-rata jumlah kemenangan sekitar 13.61% lebih tinggi dibanding candy dengan kategori Tidak, dengan asumsi prediktor lain konstan. Hubungan ini signifikan pada taraf 5%.

Sugarpercent: Setiap kenaikan 10 percentage point pada Sugarpercent diperkirakan membuat rata-rata jumlah kemenangan meningkat sekitar 1.81%, dengan asumsi prediktor lain konstan. Hubungan ini signifikan pada taraf 5%.

Pricepercent: Setiap kenaikan 10 percentage point pada Pricepercent diperkirakan membuat rata-rata jumlah kemenangan menurun sekitar 0.72%, dengan asumsi prediktor lain konstan. Hubungan ini tidak signifikan pada taraf 5%.

4.9 Uji Simultan Model

Uji simultan dilakukan untuk melihat apakah seluruh prediktor secara bersama-sama berperan dalam menjelaskan variasi win_count.

Hipotesis yang digunakan adalah:

\[ H_0: \beta_1 = \beta_2 = \cdots = \beta_p = 0 \]

\[ H_1: \text{minimal ada satu } \beta_j \neq 0 \]

Statistik uji likelihood ratio dapat ditulis sebagai:

\[ G^2 = -2(\ell_0 - \ell_1) \]

if(model_success){
  fit_null <- glm(
    win_count ~ 1,
    data = data_candy,
    family = poisson(link = "log")
  )

  lr_test <- anova(fit_null, fit_pois, test = "Chisq")
  lr_test %>%
    kable(
      caption = "Uji simultan likelihood ratio"
    ) %>%
    kable_styling(
      bootstrap_options = c("striped", "hover", "condensed", "responsive"),
      full_width = FALSE
    )
}

Uji simultan likelihood ratio

Resid. Df	Resid. Dev	Df	Deviance	Pr(>Chi)
84	364.3122	NA	NA	NA
79	196.8349	5	167.4773	0

Interpretasi

Interpretasi uji simultan: p-value = 0 < 0,05, sehingga H0 ditolak. Artinya, secara bersama-sama prediktor dalam model berpengaruh signifikan terhadap rata-rata jumlah kemenangan candy.

4.10 Pemeriksaan Overdispersion

Salah satu asumsi penting regresi Poisson adalah equidispersion, yaitu:

\[ E(Y_i \mid X_i) = Var(Y_i \mid X_i) \]

Dalam praktik, data count sering mengalami overdispersion, yaitu varians lebih besar daripada mean. Indikasi overdispersion dapat diperiksa menggunakan statistik:

\[ \hat{\phi} = \frac{\sum r_{P,i}^2}{df} \]

dengan \(r_{P,i}\) adalah residual Pearson dan \(df\) adalah derajat bebas residual.

if(model_success){
  dispersion_pois <- sum(residuals(fit_pois, type = "pearson")^2) / df.residual(fit_pois)

  dispersion_table <- tibble::tibble(
    `Dispersion Pearson` = dispersion_pois,
    `Interpretasi ringkas` = dplyr::case_when(
      dispersion_pois < 1.5 ~ "Tidak ada indikasi overdispersion berat",
      dispersion_pois < 2.5 ~ "Ada indikasi overdispersion sedang",
      TRUE ~ "Ada indikasi overdispersion kuat"
    )
  ) %>%
    dplyr::mutate(`Dispersion Pearson` = round(`Dispersion Pearson`, 3))

  dispersion_table %>%
    kable(
      caption = "Indikasi overdispersion pada model Poisson"
    ) %>%
    kable_styling(
      bootstrap_options = c("striped", "hover", "condensed", "responsive"),
      full_width = FALSE
    )
}

Indikasi overdispersion pada model Poisson

Dispersion Pearson	Interpretasi ringkas
2.486	Ada indikasi overdispersion sedang

Perhatian

Interpretasi overdispersion: terdapat indikasi overdispersion sedang. Model Poisson masih dapat dilaporkan, tetapi sebaiknya dibandingkan dengan quasi-Poisson atau negative binomial untuk analisis lanjutan.

4.11 Prediksi Nilai Rata-Rata Win Count

if(model_success){
  data_pred <- data_candy %>%
    dplyr::mutate(
      prediksi_mean = predict(fit_pois, type = "response"),
      residual_pearson = residuals(fit_pois, type = "pearson")
    )

  head(data_pred, 10) %>%
    dplyr::select(win_count, prediksi_mean, residual_pearson, chocolate, caramel, peanutyalmondy, sugarpercent, pricepercent) %>%
    dplyr::mutate(
      prediksi_mean = round(prediksi_mean, 3),
      residual_pearson = round(residual_pearson, 3),
      sugarpercent = round(sugarpercent, 3),
      pricepercent = round(pricepercent, 3)
    ) %>%
    kable(
      caption = "Contoh hasil prediksi model Poisson"
    ) %>%
    kable_styling(
      bootstrap_options = c("striped", "hover", "condensed", "responsive"),
      full_width = TRUE
    )
}

Contoh hasil prediksi model Poisson

win_count	prediksi_mean	residual_pearson	chocolate	caramel	peanutyalmondy	sugarpercent	pricepercent
67	60.871	0.786	Ya	Ya	Tidak	0.732	0.860
68	59.347	1.123	Ya	Tidak	Tidak	0.604	0.511
32	39.101	-1.136	Tidak	Tidak	Tidak	0.011	0.116
46	37.995	1.299	Tidak	Tidak	Tidak	0.011	0.511
52	44.612	1.106	Tidak	Tidak	Tidak	0.906	0.511
50	64.553	-1.811	Ya	Tidak	Ya	0.465	0.767
57	68.045	-1.339	Ya	Ya	Ya	0.604	0.767
23	45.570	-3.343	Tidak	Tidak	Ya	0.313	0.511
38	45.219	-1.074	Tidak	Tidak	Tidak	0.906	0.325
35	44.040	-1.362	Tidak	Ya	Tidak	0.604	0.325

if(model_success){
  ggplot(data_pred, aes(x = prediksi_mean, y = win_count)) +
    geom_point(color = "#2F7F73", alpha = 0.75, size = 2.4) +
    geom_smooth(method = "lm", se = FALSE, color = "#B84F5A", linewidth = 1.1) +
    labs(
      title = "Aktual vs Prediksi Rata-Rata Win Count",
      subtitle = "Titik yang semakin dekat dengan pola naik menunjukkan prediksi model semakin masuk akal.",
      x = "Prediksi rata-rata win_count",
      y = "win_count aktual"
    ) +
    theme_candy()
}

Catatan

Interpretasi prediksi: nilai prediksi pada regresi Poisson adalah rata-rata jumlah kemenangan yang diharapkan berdasarkan karakteristik candy. Prediksi ini tidak harus berupa bilangan bulat karena yang diprediksi adalah nilai harapan atau mean dari distribusi Poisson.

4.12 Visualisasi Efek Sugarpercent dan Pricepercent

4.12.1 Efek Sugarpercent terhadap Prediksi Win Count

if(model_success){
  grid_sugar <- expand.grid(
    chocolate = factor("Ya", levels = levels(data_candy$chocolate)),
    caramel = factor("Tidak", levels = levels(data_candy$caramel)),
    peanutyalmondy = factor("Tidak", levels = levels(data_candy$peanutyalmondy)),
    sugarpercent = seq(min(data_candy$sugarpercent), max(data_candy$sugarpercent), length.out = 120),
    pricepercent = mean(data_candy$pricepercent)
  )

  pred_sugar <- predict(fit_pois, newdata = grid_sugar, type = "link", se.fit = TRUE)

  grid_sugar_plot <- grid_sugar %>%
    dplyr::mutate(
      fit_link = pred_sugar$fit,
      se_link = pred_sugar$se.fit,
      pred_mean = exp(fit_link),
      low = exp(fit_link - 1.96 * se_link),
      high = exp(fit_link + 1.96 * se_link)
    )

  ggplot(grid_sugar_plot, aes(x = sugarpercent, y = pred_mean)) +
    geom_ribbon(aes(ymin = low, ymax = high), fill = "#2F7F73", alpha = 0.16) +
    geom_line(color = "#2F7F73", linewidth = 1.3) +
    scale_x_continuous(labels = percent_format(accuracy = 1)) +
    labs(
      title = "Prediksi Win Count berdasarkan Sugarpercent",
      subtitle = "Variabel lain ditahan pada nilai tertentu/rata-rata.",
      x = "Sugarpercent",
      y = "Prediksi rata-rata win_count"
    ) +
    theme_candy()
}

4.12.2 Efek Pricepercent terhadap Prediksi Win Count

if(model_success){
  grid_price <- expand.grid(
    chocolate = factor("Ya", levels = levels(data_candy$chocolate)),
    caramel = factor("Tidak", levels = levels(data_candy$caramel)),
    peanutyalmondy = factor("Tidak", levels = levels(data_candy$peanutyalmondy)),
    sugarpercent = mean(data_candy$sugarpercent),
    pricepercent = seq(min(data_candy$pricepercent), max(data_candy$pricepercent), length.out = 120)
  )

  pred_price <- predict(fit_pois, newdata = grid_price, type = "link", se.fit = TRUE)

  grid_price_plot <- grid_price %>%
    dplyr::mutate(
      fit_link = pred_price$fit,
      se_link = pred_price$se.fit,
      pred_mean = exp(fit_link),
      low = exp(fit_link - 1.96 * se_link),
      high = exp(fit_link + 1.96 * se_link)
    )

  ggplot(grid_price_plot, aes(x = pricepercent, y = pred_mean)) +
    geom_ribbon(aes(ymin = low, ymax = high), fill = "#5568B8", alpha = 0.16) +
    geom_line(color = "#5568B8", linewidth = 1.3) +
    scale_x_continuous(labels = percent_format(accuracy = 1)) +
    labs(
      title = "Prediksi Win Count berdasarkan Pricepercent",
      subtitle = "Variabel lain ditahan pada nilai tertentu/rata-rata.",
      x = "Pricepercent",
      y = "Prediksi rata-rata win_count"
    ) +
    theme_candy()
}

4.13 Perbandingan Model Poisson dan Quasi-Poisson

Jika terdapat indikasi overdispersion, model quasi-Poisson dapat digunakan untuk memperbaiki standard error. Model quasi-Poisson tidak mengubah bentuk mean, tetapi mengizinkan varians lebih fleksibel:

\[ Var(Y_i \mid X_i) = \phi \mu_i \]

if(model_success){
  fit_quasi <- glm(
    formula = model_formula,
    data = data_candy,
    family = quasipoisson(link = "log")
  )

  quasi_coef <- as.data.frame(coef(summary(fit_quasi))) %>%
    tibble::rownames_to_column("parameter") %>%
    dplyr::rename(
      estimate = Estimate,
      std_error = `Std. Error`,
      t_value = `t value`,
      p_value = `Pr(>|t|)`
    ) %>%
    dplyr::mutate(
      IRR = exp(estimate),
      perubahan_persen = 100 * (IRR - 1),
      keputusan = ifelse(p_value < 0.05, "Signifikan", "Tidak signifikan")
    )

  quasi_coef %>%
    dplyr::mutate(dplyr::across(c(estimate, std_error, t_value, p_value, IRR, perubahan_persen), ~ round(.x, 4))) %>%
    kable(
      caption = "Ringkasan model quasi-Poisson sebagai pembanding",
      col.names = c("Parameter", "Estimate", "SE", "t-value", "p-value", "IRR", "Perubahan (%)", "Keputusan")
    ) %>%
    kable_styling(
      bootstrap_options = c("striped", "hover", "condensed", "responsive"),
      full_width = TRUE
    )
}

Ringkasan model quasi-Poisson sebagai pembanding

Parameter	Estimate	SE	t-value	p-value	IRR	Perubahan (%)	Keputusan
(Intercept)	3.6726	0.0588	62.4435	0.0000	39.3544	3835.4398	Signifikan
chocolateYa	0.3396	0.0595	5.7042	0.0000	1.4043	40.4337	Signifikan
caramelYa	0.0278	0.0658	0.4219	0.6742	1.0281	2.8143	Tidak signifikan
peanutyalmondyYa	0.1276	0.0651	1.9607	0.0534	1.1361	13.6113	Tidak signifikan
sugarpercent	0.1794	0.0947	1.8942	0.0619	1.1965	19.6469	Tidak signifikan
pricepercent	-0.0726	0.1073	-0.6769	0.5005	0.9299	-7.0059	Tidak signifikan

Catatan

Interpretasi pembanding: quasi-Poisson terutama digunakan ketika terdapat overdispersion. Koefisien dan IRR biasanya sama dengan model Poisson, tetapi standard error dan p-value dapat berubah karena quasi-Poisson mengoreksi varians.

4.14 Kesimpulan

Interpretasi

Regresi Poisson digunakan untuk memodelkan win_count sebagai data hitungan. Prediktor yang signifikan pada taraf 5% adalah chocolateYa, peanutyalmondyYa, sugarpercent. Nilai IRR digunakan untuk membaca arah dan besar pengaruh prediktor terhadap rata-rata jumlah kemenangan candy. Pemeriksaan overdispersion perlu diperhatikan karena jika nilai dispersion terlalu besar, quasi-Poisson atau negative binomial lebih disarankan sebagai model lanjutan.

5 Kesimpulan Umum

Regresi logistik biner digunakan untuk respons dua kategori, regresi logistik multinomial digunakan untuk respons lebih dari dua kategori tanpa urutan, regresi logistik ordinal digunakan untuk respons yang memiliki tingkatan, sedangkan regresi Poisson digunakan untuk respons berbentuk hitungan.