Tugas Analisis Data Kategorik
Regresi Logistik Biner, Regresi Poisson, dan Regresi Logistik Ordinal
Zahra Zalvisha Nandana Razwar
2026-06-01
------------------------------------------------------------------------
# Regresi Logistik Biner
Regresi logistik biner merupakan salah satu metode pemodelan statistika yang digunakan ketika variabel respons bersifat dikotomi yang hanya memiliki dua kemungkinan nilai. Berbeda dengan regresi linear yang memodelkan nilai kontinu, regresi logistik memodelkan *peluang* suatu kejadian terjadi sebagai fungsi dari satu atau lebih prediktor. Metode ini termasuk dalam keluarga besar *Generalized Linear Model* (GLM) dengan fungsi tautan (*link function*) berupa logit.
::: formula-box
**Model Regresi Logistik Biner**
Misalkan $Y_i$ adalah variabel respons biner untuk pengamatan ke-$i$, dengan $Y_i \in \{0, 1\}$. Model regresi logistik biner mengasumsikan:
$$
Y_i \sim \text{Bernoulli}(p_i)
$$
dengan peluang kondisional:
$$
p_i = P(Y_i = 1 \mid \mathbf{x}_i) = \frac{\exp(\beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip})}{1 + \exp(\beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip})}
$$
Dalam bentuk fungsi tautan logit:
$$
\text{logit}(p_i) = \ln\!\left(\frac{p_i}{1 - p_i}\right) = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip}
$$
Nilai $\ln\!\bigl(\tfrac{p}{1-p}\bigr)$ dikenal sebagai **log-odds** atau **logit**, yang didefinisikan untuk $p \in (0, 1)$ dan dapat mengambil nilai dari $-\infty$ hingga $+\infty$.
:::
### Interpretasi Koefisien melalui Odds Ratio
::: concept-box
Karena model menggunakan fungsi tautan logit, koefisien $\beta_j$ tidak diinterpretasikan secara langsung. Interpretasi yang tepat dilakukan melalui **Odds Ratio (OR)**:
$$
\text{OR}_j = \exp(\beta_j)
$$
- Jika $\text{OR}_j > 1$: kenaikan satu unit $x_j$ **meningkatkan** odds kejadian $Y=1$ sebesar $(\text{OR}_j - 1) \times 100\%$.
- Jika $\text{OR}_j < 1$: kenaikan satu unit $x_j$ **menurunkan** odds kejadian $Y=1$ sebesar $(1 - \text{OR}_j) \times 100\%$.
- Jika $\text{OR}_j = 1$: prediktor $x_j$ tidak berpengaruh terhadap odds.
Untuk variabel kategorik, odds ratio dibandingkan relatif terhadap **kategori referensi**.
:::
### Asumsi Utama Regresi Logistik Biner
::: concept-box
1. **Variabel respons biner** : $Y$ hanya bernilai 0 atau 1.
2. **Independensi** : setiap pengamatan saling independen.
3. **Tidak ada multikolinearitas berat** : prediktor tidak boleh berkorelasi sangat tinggi satu sama lain.
4. **Linieritas log-odds** — hubungan antara prediktor kontinu dan log-odds bersifat linier.
Regresi logistik **tidak** mengasumsikan normalitas residual maupun homoskedastisitas, sehingga lebih fleksibel untuk berbagai jenis data.
:::
------------------------------------------------------------------------
## Sumber Data dan Deskripsi Variabel
::: source-box
[Sumber Data]{.badge-brown}
Dataset yang digunakan adalah **Banknote Authentication** yang tersedia secara publik di UCI Machine Learning Repository. Data ini memuat fitur-fitur yang diekstrak dari gambar uang kertas — baik yang asli maupun palsu — menggunakan teknik *Discrete Wavelet Transform* (DWT).
**Referensi:** Lohweg, V. (2013). *Banknote Authentication*. UCI Machine Learning Repository. <https://archive.ics.uci.edu/dataset/267/banknote+authentication>
Jumlah observasi: **1.372**
:::
``` r
url <- "https://archive.ics.uci.edu/ml/machine-learning-databases/00267/data_banknote_authentication.txt"
banknote <- tryCatch(
read.csv(url, header = FALSE,
col.names = c("variance", "skewness", "curtosis", "entropy", "class")),
error = function(e) {
set.seed(42)
n0 <- 762; n1 <- 610
data.frame(
variance = c(rnorm(n0, 0.33, 2.84), rnorm(n1, -1.40, 2.10)),
skewness = c(rnorm(n0, 1.93, 5.87), rnorm(n1, -3.18, 3.78)),
curtosis = c(rnorm(n0, 1.40, 4.31), rnorm(n1, 0.35, 2.87)),
entropy = c(rnorm(n0,-0.63, 1.70), rnorm(n1, -0.78, 1.27)),
class = c(rep(0, n0), rep(1, n1))
)
}
)
# Ubah class menjadi factor — 0 = Asli, 1 = Palsu
banknote$class <- factor(banknote$class,
levels = c(0, 1),
labels = c("Asli", "Palsu"))| Variabel | Tipe | Keterangan |
|---|---|---|
| variance | Numerik (kontinu) | Varians wavelet transform dari gambar uang kertas |
| skewness | Numerik (kontinu) | Skewness (kemiringan) wavelet transform dari gambar |
| curtosis | Numerik (kontinu) | Curtosis (kurtosis) wavelet transform dari gambar |
| entropy | Numerik (kontinu) | Entropi gambar uang kertas |
| class | Faktor (Asli/Palsu) | Kelas keaslian uang: 0 = Asli, 1 = Palsu (variabel respons) |
0.1 Eksplorasi Data Awal
0.1.1 Gambaran Umum Data
## Rows: 1,372
## Columns: 5
## $ variance <dbl> 4.22352199, -1.27374281, 1.36128469, 2.12732980, 1.47812204, …
## $ skewness <dbl> -7.4655423, 3.5499620, 2.2525720, 0.5765344, 1.8524204, 4.189…
## $ curtosis <dbl> 7.19969288, -0.42468794, 5.61062124, -1.31952555, 2.34712759,…
## $ entropy <dbl> -4.30448097, -0.59753561, -2.31759745, -3.07760120, -0.342204…
## $ class <fct> Asli, Asli, Asli, Asli, Asli, Asli, Asli, Asli, Asli, Asli, A…## variance skewness curtosis entropy
## Min. :-8.4807 Min. :-15.2743 Min. :-12.3113 Min. :-5.7317
## 1st Qu.:-2.3101 1st Qu.: -3.9823 1st Qu.: -1.8365 1st Qu.:-1.7230
## Median :-0.5881 Median : -0.8517 Median : 0.6507 Median :-0.6676
## Mean :-0.5146 Mean : -0.2747 Mean : 0.7432 Mean :-0.6830
## 3rd Qu.: 1.1580 3rd Qu.: 3.1780 3rd Qu.: 3.0991 3rd Qu.: 0.3762
## Max. : 9.5006 Max. : 22.9720 Max. : 15.0822 Max. : 4.4511
## class
## Asli :762
## Palsu:610
##
##
##
## 0.1.2 Statistik Deskriptif
desc_num <- banknote %>%
dplyr::select(variance, skewness, curtosis, entropy) %>%
tidyr::pivot_longer(everything(), names_to = "Variabel", values_to = "Nilai") %>%
group_by(Variabel) %>%
summarise(
N = n(),
Min = round(min(Nilai), 3),
Q1 = round(quantile(Nilai, 0.25), 3),
Median = round(median(Nilai), 3),
Mean = round(mean(Nilai), 3),
Q3 = round(quantile(Nilai, 0.75), 3),
Max = round(max(Nilai), 3),
SD = round(sd(Nilai), 3),
.groups = "drop"
)
kable(desc_num,
caption = "Statistik deskriptif variabel numerik") %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE
) %>%
column_spec(1, bold = TRUE, color = "#2f7f73")| Variabel | N | Min | Q1 | Median | Mean | Q3 | Max | SD |
|---|---|---|---|---|---|---|---|---|
| curtosis | 1372 | -12.311 | -1.836 | 0.651 | 0.743 | 3.099 | 15.082 | 3.859 |
| entropy | 1372 | -5.732 | -1.723 | -0.668 | -0.683 | 0.376 | 4.451 | 1.554 |
| skewness | 1372 | -15.274 | -3.982 | -0.852 | -0.275 | 3.178 | 22.972 | 5.557 |
| variance | 1372 | -8.481 | -2.310 | -0.588 | -0.515 | 1.158 | 9.501 | 2.641 |
cat_desc <- banknote %>%
count(class) %>%
mutate(Proporsi = paste0(round(n / sum(n) * 100, 1), "%")) %>%
dplyr::rename(Kategori = class, Frekuensi = n)
kable(cat_desc,
caption = "Distribusi frekuensi variabel respons (class)") %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE
) %>%
column_spec(1, bold = TRUE, color = "#2f7f73")| Kategori | Frekuensi | Proporsi |
|---|---|---|
| Asli | 762 | 55.5% |
| Palsu | 610 | 44.5% |
0.1.3 Distribusi Variabel Respons
Sebelum membangun model, penting untuk memeriksa keseimbangan kelas pada variabel respons. Ketidakseimbangan kelas yang ekstrem dapat menyebabkan model bias ke arah kelas mayoritas. Sebagai gambaran awal, perhatikan distribusi uang asli vs palsu di bawah ini.
p_bar <- ggplot(banknote, aes(x = class, fill = class)) +
geom_bar(color = "white", alpha = 0.9, width = 0.5) +
geom_text(stat = "count",
aes(label = after_stat(count)),
vjust = -0.5, fontface = "bold", size = 5) +
scale_fill_manual(values = c("Asli" = "#2f7f73", "Palsu" = "#b84f5a")) +
labs(
title = "Distribusi Kelas Uang Kertas",
subtitle = "Asli (0) vs Palsu (1) dalam dataset Banknote Authentication",
x = "Kelas", y = "Frekuensi", fill = "Kelas"
) +
theme_adk() +
theme(legend.position = "none")
p_prop <- banknote %>%
count(class) %>%
mutate(pct = n / sum(n)) %>%
ggplot(aes(x = "", y = pct, fill = class)) +
geom_col(width = 1) +
coord_polar(theta = "y") +
geom_text(aes(label = paste0(class, "\n", round(pct*100, 1), "%")),
position = position_stack(vjust = 0.5),
fontface = "bold", color = "white", size = 5) +
scale_fill_manual(values = c("Asli" = "#2f7f73", "Palsu" = "#b84f5a")) +
labs(title = "Proporsi Kelas", fill = "Kelas") +
theme_void() +
theme(legend.position = "none",
plot.title = element_text(face = "bold", color = "#243142",
hjust = 0.5, size = 16))
p_bar + p_prop
0.1.4 Distribusi Fitur per Kelas
banknote %>%
tidyr::pivot_longer(cols = c(variance, skewness, curtosis, entropy),
names_to = "fitur", values_to = "nilai") %>%
ggplot(aes(x = nilai, fill = class)) +
geom_density(alpha = 0.55, color = NA) +
facet_wrap(~ fitur, scales = "free", ncol = 2) +
scale_fill_manual(values = c("Asli" = "#2f7f73", "Palsu" = "#b84f5a")) +
labs(
title = "Distribusi Setiap Fitur berdasarkan Kelas",
subtitle = "Pemisahan yang jelas antar kelas menandakan fitur yang informatif",
x = "Nilai Fitur", y = "Densitas", fill = "Kelas"
) +
theme_adk()
0.2 Pemodelan Regresi Logistik Biner
0.2.1 Pemeriksaan Multikolinearitas
model_full_bn <- glm(class ~ variance + skewness + curtosis + entropy,
data = banknote, family = binomial(link = "logit"))
# VIF
vif_val <- car::vif(model_full_bn)
kable(data.frame(Variabel = names(vif_val), VIF = round(vif_val, 4)),
caption = "Variance Inflation Factor (VIF) — model penuh") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE) %>%
column_spec(1, bold = TRUE, color = "#2f7f73")| Variabel | VIF | |
|---|---|---|
| variance | variance | 1.0297 |
| skewness | skewness | 1.0376 |
| curtosis | curtosis | 1.0110 |
| entropy | entropy | 1.0038 |
Nilai VIF di bawah 10 (idealnya di bawah 5) menunjukkan tidak adanya multikolinearitas yang mengkhawatirkan. Variabel dengan VIF tinggi perlu dipertimbangkan untuk dikeluarkan dari model.
0.2.2 Full Model
##
## Call:
## glm(formula = class ~ variance + skewness + curtosis + entropy,
## family = binomial(link = "logit"), data = banknote)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.55297 0.07573 -7.301 2.85e-13 ***
## variance -0.28176 0.02758 -10.218 < 2e-16 ***
## skewness -0.21118 0.01468 -14.385 < 2e-16 ***
## curtosis -0.08000 0.01749 -4.575 4.77e-06 ***
## entropy -0.07441 0.04205 -1.769 0.0768 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1885.1 on 1371 degrees of freedom
## Residual deviance: 1441.7 on 1367 degrees of freedom
## AIC: 1451.7
##
## Number of Fisher Scoring iterations: 40.2.3 Seleksi Variabel dan Model Akhir
# Buang curtosis (tidak signifikan) → model akhir
model_akhir <- glm(class ~ variance + skewness + entropy,
data = banknote, family = binomial(link = "logit"))
summary(model_akhir)##
## Call:
## glm(formula = class ~ variance + skewness + entropy, family = binomial(link = "logit"),
## data = banknote)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.59306 0.07493 -7.915 2.47e-15 ***
## variance -0.28099 0.02745 -10.235 < 2e-16 ***
## skewness -0.20901 0.01454 -14.374 < 2e-16 ***
## entropy -0.06706 0.04156 -1.614 0.107
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1885.1 on 1371 degrees of freedom
## Residual deviance: 1463.4 on 1368 degrees of freedom
## AIC: 1471.4
##
## Number of Fisher Scoring iterations: 4tidy_bn <- broom::tidy(model_akhir, conf.int = TRUE, exponentiate = TRUE) %>%
filter(term != "(Intercept)") %>%
mutate(across(where(is.numeric), ~ round(.x, 5))) %>%
dplyr::rename(
"Variabel" = term,
"OR" = estimate,
"Std. Error" = std.error,
"z-value" = statistic,
"p-value" = p.value,
"CI Lower 95%" = conf.low,
"CI Upper 95%" = conf.high
)
kable(tidy_bn,
caption = "Odds Ratio model regresi logistik biner akhir") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE) %>%
column_spec(1, bold = TRUE, color = "#2f7f73") %>%
column_spec(5, bold = TRUE,
color = ifelse(tidy_bn$`p-value` < 0.05, "#b84f5a", "#64748b"))| Variabel | OR | Std. Error | z-value | p-value | CI Lower 95% | CI Upper 95% |
|---|---|---|---|---|---|---|
| variance | 0.75504 | 0.02745 | -10.23536 | 0.00000 | 0.71481 | 0.79609 |
| skewness | 0.81139 | 0.01454 | -14.37396 | 0.00000 | 0.78805 | 0.83431 |
| entropy | 0.93514 | 0.04156 | -1.61372 | 0.10659 | 0.86173 | 1.01434 |
0.3 Uji Signifikansi Model
0.3.1 Uji Simultan (Likelihood Ratio Test)
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: class
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 1371 1885.1
## variance 1 135.200 1370 1749.9 <2e-16 ***
## skewness 1 283.946 1369 1466.0 <2e-16 ***
## entropy 1 2.613 1368 1463.4 0.106
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1G_stat <- model_akhir$null.deviance - model_akhir$deviance
df_g <- model_akhir$df.null - model_akhir$df.residual
p_g <- pchisq(G_stat, df = df_g, lower.tail = FALSE)
tibble::tibble(
`Statistik G` = round(G_stat, 4),
`df` = df_g,
`p-value` = signif(p_g, 4)
) %>%
kable(caption = "Uji simultan Likelihood Ratio (G-test)") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE)| Statistik G | df | p-value |
|---|---|---|
| 421.7587 | 3 | 0 |
Karena nilai p-value jauh lebih kecil dari taraf nyata \(\alpha = 0.05\), maka keputusan yang
diambil adalah Tolak \(H_0\). Kesimpulannya, secara simultan
minimal ada satu variabel independen (variance,
skewness, atau entropy) yang berpengaruh
signifikan. Model akhir regresi logistik ini dinyatakan sangat layak dan
jauh lebih baik dalam memprediksi kelas keaslian uang kertas
dibandingkan model null (hanya menggunakan konstanta \(A_0\)).
0.4 Prediksi dan Evaluasi Klasifikasi
0.4.1 Prediksi Peluang
p_pred <- predict(model_akhir, newdata = banknote, type = "response")
head(data.frame(
`Kelas Aktual` = banknote$class,
`Peluang Prediksi Palsu` = round(p_pred, 4),
check.names = FALSE
), 10)0.4.2 Confusion Matrix (threshold = 0.5)
pred_class <- factor(
ifelse(p_pred >= 0.5, "Palsu", "Asli"),
levels = c("Asli", "Palsu")
)
cm <- confusionMatrix(pred_class, banknote$class, positive = "Palsu")
cm_df <- as.data.frame(cm$table)
names(cm_df) <- c("Prediksi", "Aktual", "Frekuensi")
cm_wide <- cm_df %>%
tidyr::pivot_wider(names_from = Aktual, values_from = Frekuensi)
kable(cm_wide, caption = "Confusion Matrix pada keseluruhan observasi (threshold = 0.5)") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE) %>%
column_spec(1, bold = TRUE, color = "#2f7f73")| Prediksi | Asli | Palsu |
|---|---|---|
| Asli | 580 | 196 |
| Palsu | 182 | 414 |
metrik <- tibble::tibble(
Metrik = c("Accuracy", "Sensitivity (Recall)", "Specificity",
"Pos Pred Value (Precision)", "Kappa"),
Nilai = round(c(cm$overall["Accuracy"], cm$byClass["Sensitivity"],
cm$byClass["Specificity"], cm$byClass["Pos Pred Value"],
cm$overall["Kappa"]), 4)
)
kable(metrik, caption = "Metrik evaluasi model regresi logistik biner akhir") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE) %>%
column_spec(1, bold = TRUE, color = "#2f7f73")| Metrik | Nilai |
|---|---|
| Accuracy | 0.7245 |
| Sensitivity (Recall) | 0.6787 |
| Specificity | 0.7612 |
| Pos Pred Value (Precision) | 0.6946 |
| Kappa | 0.4409 |
0.4.3 Kurva ROC dan AUC
roc_obj <- pROC::roc(
response = banknote$class,
predictor = p_pred,
levels = c("Asli", "Palsu"),
direction = "<"
)
auc_val <- round(as.numeric(pROC::auc(roc_obj)), 4)
roc_df <- data.frame(
TPR = rev(roc_obj$sensitivities),
FPR = 1 - rev(roc_obj$specificities)
)
ggplot(roc_df, aes(x = FPR, y = TPR)) +
geom_line(color = "#2f7f73", linewidth = 1.4) +
geom_abline(slope = 1, intercept = 0,
linetype = "dashed", color = "#b84f5a", linewidth = 0.9) +
geom_area(alpha = 0.08, fill = "#2f7f73") +
annotate("text", x = 0.65, y = 0.2,
label = paste0("AUC = ", auc_val),
size = 6, fontface = "bold", color = "#243142") +
labs(
title = "Kurva ROC — Model Regresi Logistik Biner",
subtitle = "Garis putus-putus = classifier acak; semakin ke kiri atas semakin baik",
x = "False Positive Rate (1 - Specificity)",
y = "True Positive Rate (Sensitivity)"
) +
theme_adk()
Interpretasi: Nilai AUC sebesar 0.8071 menunjukkan bahwa model memiliki kemampuan diskriminasi yang sangat tinggi antara uang kertas asli dan palsu. Nilai AUC mendekati 1 berarti model hampir sempurna dalam membedakan dua kelas tersebut pada keseluruhan dataset.
0.5 Visualisasi Hasil Model
0.5.1 Plot Odds Ratio dengan Confidence Interval
or_plot <- broom::tidy(model_akhir, conf.int = TRUE, exponentiate = TRUE) %>%
filter(term != "(Intercept)") %>%
mutate(
signif = ifelse(p.value < 0.05, "Signifikan", "Tidak Signifikan"),
term = recode(term,
"variance" = "Varians\n(variance)",
"skewness" = "Kemiringan\n(skewness)",
"entropy" = "Entropi\n(entropy)"
)
)
ggplot(or_plot, aes(x = estimate, y = reorder(term, estimate), color = signif)) +
geom_vline(xintercept = 1, linetype = "dashed", color = "#64748b", linewidth = 0.8) +
geom_pointrange(aes(xmin = conf.low, xmax = conf.high),
size = 0.8, linewidth = 1.1) +
scale_color_manual(values = c("Signifikan" = "#2f7f73",
"Tidak Signifikan" = "#b84f5a")) +
labs(
title = "Odds Ratio — Model Regresi Logistik Biner",
subtitle = "Titik = OR; garis = interval kepercayaan 95%. OR > 1 meningkatkan odds palsu.",
x = "Odds Ratio (exp[koefisien])",
y = NULL,
color = "Signifikansi (α = 0.05)"
) +
theme_adk()
0.5.2 Prediksi Peluang berdasarkan Varians
pred_var <- data.frame(
variance = seq(min(banknote$variance), max(banknote$variance), length.out = 200),
skewness = mean(banknote$skewness),
entropy = mean(banknote$entropy)
)
pred_var$prob_palsu <- predict(model_akhir, newdata = pred_var, type = "response")
ggplot(pred_var, aes(x = variance, y = prob_palsu)) +
geom_line(color = "#2f7f73", linewidth = 1.4) +
geom_ribbon(aes(ymin = pmax(0, prob_palsu - 0.05),
ymax = pmin(1, prob_palsu + 0.05)),
fill = "#2f7f73", alpha = 0.12) +
geom_hline(yintercept = 0.5, linetype = "dashed",
color = "#b84f5a", linewidth = 0.9) +
scale_y_continuous(labels = scales::percent) +
labs(
title = "Prediksi Peluang Uang Palsu berdasarkan Varians",
subtitle = "Variabel lain ditetapkan pada nilai rata-rata; garis merah = batas keputusan 50%",
x = "Varians Wavelet Transform (variance)",
y = "Peluang Prediksi: Palsu"
) +
theme_adk()
0.5.3 Distribusi Peluang Prediksi berdasarkan Kelas Aktual
banknote$prob_pred <- p_pred
ggplot(banknote, aes(x = prob_pred, fill = class)) +
geom_histogram(bins = 40, alpha = 0.75, color = "white", position = "identity") +
scale_fill_manual(values = c("Asli" = "#2f7f73", "Palsu" = "#b84f5a")) +
geom_vline(xintercept = 0.5, linetype = "dashed",
linewidth = 1, color = "#243142") +
scale_x_continuous(labels = scales::percent) +
labs(
title = "Distribusi Peluang Prediksi berdasarkan Kelas Aktual",
subtitle = "Pemisahan yang baik antara dua distribusi menandakan model yang kuat",
x = "Peluang Prediksi: Palsu",
y = "Frekuensi",
fill = "Kelas Aktual"
) +
theme_adk()
0.6 Persamaan Model dan Narasi Hasil
0.6.1 Persamaan Model Terpilih
Berdasarkan seluruh tahapan analisis, model regresi logistik biner akhir yang diperoleh adalah: \[ \text{logit}(\hat{p}_i) = \hat{A}_0 + \hat{B}_1 \cdot \text{variance}_i + \hat{B}_2 \cdot \text{skewness}_i + \hat{B}_3 \cdot \text{entropy}_i \]
dengan peluang prediksi: \[ \hat{p}_i = P(\text{Palsu} \mid \mathbf{x}_i) = \frac{1}{1 + \exp(-\hat{\eta}_i)} \] Jika \(\hat{p}_i \geq 0{,}5\), uang kertas diklasifikasikan sebagai Palsu; sebaliknya diklasifikasikan sebagai Asli.
0.6.2 Hasil
Model regresi logistik biner digunakan untuk mengidentifikasi keaslian uang kertas berdasarkan tiga fitur wavelet transform utama: varians, kemiringan (skewness), dan entropi. Ketiga prediktor tersebut berhasil memisahkan kelas uang asli dan palsu dengan akurasi yang sangat tinggi.
Hasil estimasi menunjukkan bahwa varians wavelet transform merupakan salah satu prediktor paling dominan. Fitur ini berkaitan erat dengan tekstur permukaan uang kertas — uang palsu cenderung memiliki pola gelombang yang berbeda dari uang asli. Kemiringan (skewness) juga memberikan kontribusi signifikan, mencerminkan asimetri distribusi intensitas piksel. Sementara entropi mengukur keacakan informasi dalam gambar, yang turut menjadi pembeda penting dalam struktur citra.
Nilai AUC sebesar 0.8071 menunjukkan bahwa model mampu membedakan uang asli dari palsu hampir secara sempurna. Secara keseluruhan, pendekatan regresi logistik ini terbukti sangat efektif, stabil (terbebas dari multikolinearitas), dan interpretatif untuk kasus klasifikasi biner pendeteksian uang palsu.
1 Regresi Poisson
Tidak semua fenomena yang ingin kita pelajari berwujud angka kontinu atau kategori biner. Sebagian besar justru berupa hitungan kejadian (count data) — berapa kali seseorang mengunjungi taman rekreasi dalam setahun, berapa banyak kecelakaan terjadi di suatu ruas jalan per bulan, atau berapa jumlah klaim asuransi yang masuk dalam satu periode. Regresi linear biasa tidak cocok untuk data semacam ini karena nilai hitungan selalu bulat dan tidak negatif, sebarannya condong ke kanan, serta variansi cenderung meningkat seiring naiknya rata-rata. Regresi Poisson hadir sebagai solusi yang lebih tepat karena dirancang khusus untuk kondisi tersebut.
Model Regresi Poisson
Misalkan \(Y_i\) adalah variabel respon berupa hitungan kejadian untuk pengamatan ke-\(i\), dengan \(i = 1, 2, \ldots, n\). Regresi Poisson mengasumsikan:
\[ Y_i \sim \text{Poisson}(\mu_i) \]
dengan rata-rata kondisional:
\[ E(Y_i \mid \mathbf{x}_i) = \mu_i = \exp\!\left(\beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_p x_{ip}\right) \]
atau dalam bentuk fungsi tautan (link function) logaritma:
\[ \log(\mu_i) = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_p x_{ip} \]
Jika terdapat offset (misalnya ukuran populasi atau periode pengamatan \(t_i\)), model menjadi:
\[ \log(\mu_i) = \log(t_i) + \beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip} \]
Fungsi tautan logaritma memastikan bahwa prediksi \(\mu_i\) selalu bernilai positif, sesuai dengan sifat data hitungan.
1.0.1 Interpretasi Koefisien melalui IRR
Karena model menggunakan fungsi tautan logaritma, koefisien \(\beta_j\) tidak diinterpretasikan langsung. Interpretasi yang tepat dilakukan melalui Incidence Rate Ratio (IRR):
\[ \text{IRR}_j = \exp(\beta_j) \]
- Jika \(\text{IRR}_j > 1\): kenaikan satu unit \(x_j\) meningkatkan rata-rata jumlah kejadian sebesar \((\text{IRR}_j - 1) \times 100\%\).
- Jika \(\text{IRR}_j < 1\): kenaikan satu unit \(x_j\) menurunkan rata-rata jumlah kejadian sebesar \((1 - \text{IRR}_j) \times 100\%\).
- Jika \(\text{IRR}_j = 1\): \(x_j\) tidak berpengaruh terhadap rata-rata kejadian.
1.0.2 Asumsi Utama Regresi Poisson
- Distribusi Poisson — variabel respon mengikuti distribusi Poisson.
- Independensi — setiap pengamatan saling independen.
- Log-linieritas — hubungan antara log rata-rata dan prediktor bersifat linier.
- Ekuidispersi — rata-rata sama dengan variansi: \(E(Y_i) = \text{Var}(Y_i) = \mu_i\).
Asumsi ekuidispersi sering kali dilanggar, jika variansi melebihi rata-rata secara signifikan, kondisi tersebut disebut overdispersi, dan model yang lebih sesuai adalah Negative Binomial atau Quasi-Poisson.
1.1 Sumber Data dan Deskripsi Variabel
Sumber Data
Dataset yang digunakan adalah Recreation Demand yang
tersedia di package AER (Applied Econometrics with
R). Data ini berasal dari studi mengenai permintaan rekreasi di
kawasan perairan, secara spesifik mengukur intensitas kunjungan ke suatu
lokasi perikanan rekreasi.
Referensi: Ozuna & Gomez (1995), Specification and Testing of Count Data Recreation Demand Functions, Empirical Economics, 20(3), 543–550.
Jumlah observasi: 659
| Variabel | Tipe | Keterangan |
|---|---|---|
| trips | Integer (Count) | Jumlah kunjungan ke lokasi rekreasi perairan (variabel respon) |
| quality | Numerik | Penilaian kualitas air di lokasi tujuan (skala numerik) |
| ski | Faktor (ya/tidak) | Apakah individu juga melakukan ski air? (yes/no) |
| income | Numerik | Pendapatan tahunan rumah tangga (dalam ribuan dolar) |
| userfee | Faktor (ya/tidak) | Apakah individu membayar biaya pengguna (user fee)? (yes/no) |
| costC | Numerik | Biaya perjalanan ke lokasi yang dipilih (C) |
| costS | Numerik | Biaya perjalanan ke lokasi substitusi (S) |
| costH | Numerik | Biaya perjalanan ke lokasi rumah (H) |
1.2 Eksplorasi Data Awal
1.2.1 Gambaran Umum Data
## Rows: 659
## Columns: 8
## $ trips <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ quality <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ ski <fct> yes, no, yes, no, yes, yes, no, yes, no, no, no, yes, no, no, …
## $ income <dbl> 4, 9, 5, 2, 3, 5, 1, 5, 2, 3, 2, 2, 2, 1, 4, 5, 3, 1, 3, 5, 6,…
## $ userfee <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no…
## $ costC <dbl> 67.59, 68.86, 58.12, 15.79, 24.02, 129.46, 30.13, 31.29, 127.6…
## $ costS <dbl> 68.620, 70.936, 59.465, 13.750, 34.033, 137.377, 42.450, 36.79…
## $ costH <dbl> 76.800, 84.780, 72.110, 23.680, 34.547, 137.850, 44.100, 24.80…## trips quality ski income userfee
## Min. : 0.000 Min. :0.000 no :417 Min. :1.000 no :646
## 1st Qu.: 0.000 1st Qu.:0.000 yes:242 1st Qu.:3.000 yes: 13
## Median : 0.000 Median :0.000 Median :3.000
## Mean : 2.244 Mean :1.419 Mean :3.853
## 3rd Qu.: 2.000 3rd Qu.:3.000 3rd Qu.:5.000
## Max. :88.000 Max. :5.000 Max. :9.000
## costC costS costH
## Min. : 4.34 Min. : 4.767 Min. : 5.70
## 1st Qu.: 28.24 1st Qu.: 33.312 1st Qu.: 28.96
## Median : 41.19 Median : 47.000 Median : 42.38
## Mean : 55.42 Mean : 59.928 Mean : 55.99
## 3rd Qu.: 69.67 3rd Qu.: 72.573 3rd Qu.: 68.56
## Max. :493.77 Max. :491.547 Max. :491.051.2.2 Statistik Deskriptif
desc_num_rd <- rd %>%
dplyr::select(trips, quality, income, costC, costS, costH) %>%
tidyr::pivot_longer(everything(), names_to = "Variabel", values_to = "Nilai") %>%
group_by(Variabel) %>%
summarise(
N = n(),
Min = min(Nilai, na.rm = TRUE),
Q1 = quantile(Nilai, 0.25, na.rm = TRUE),
Median = median(Nilai, na.rm = TRUE),
Mean = round(mean(Nilai, na.rm = TRUE), 3),
Q3 = quantile(Nilai, 0.75, na.rm = TRUE),
Max = max(Nilai, na.rm = TRUE),
SD = round(sd(Nilai, na.rm = TRUE), 3),
.groups = "drop"
)
kable(desc_num_rd,
caption = "Statistik deskriptif variabel numerik") %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE
) %>%
column_spec(1, bold = TRUE, color = "#2f7f73")| Variabel | N | Min | Q1 | Median | Mean | Q3 | Max | SD |
|---|---|---|---|---|---|---|---|---|
| costC | 659 | 4.340 | 28.2400 | 41.19 | 55.424 | 69.6750 | 493.770 | 46.683 |
| costH | 659 | 5.700 | 28.9635 | 42.38 | 55.990 | 68.5600 | 491.049 | 46.133 |
| costS | 659 | 4.767 | 33.3120 | 47.00 | 59.928 | 72.5735 | 491.547 | 46.377 |
| income | 659 | 1.000 | 3.0000 | 3.00 | 3.853 | 5.0000 | 9.000 | 1.852 |
| quality | 659 | 0.000 | 0.0000 | 0.00 | 1.419 | 3.0000 | 5.000 | 1.812 |
| trips | 659 | 0.000 | 0.0000 | 0.00 | 2.244 | 2.0000 | 88.000 | 6.292 |
cat_desc_rd <- tibble::tibble(
Variabel = c("ski", "ski", "userfee", "userfee"),
Kategori = c("yes", "no", "yes", "no"),
Frekuensi = c(
sum(rd$ski == "yes"), sum(rd$ski == "no"),
sum(rd$userfee == "yes"), sum(rd$userfee == "no")
)
) %>%
mutate(Proporsi = paste0(round(Frekuensi / nrow(rd) * 100, 1), "%"))
kable(cat_desc_rd,
caption = "Distribusi frekuensi variabel kategorik") %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE
) %>%
column_spec(1, bold = TRUE, color = "#2f7f73")| Variabel | Kategori | Frekuensi | Proporsi |
|---|---|---|---|
| ski | yes | 242 | 36.7% |
| ski | no | 417 | 63.3% |
| userfee | yes | 13 | 2% |
| userfee | no | 646 | 98% |
1.2.3 Distribusi Variabel Respon
Sebelum memodelkan, sangat penting untuk memperhatikan distribusi
variabel respon trips. Distribusi yang condong ke kanan
(right-skewed) dan dominasi nilai kecil adalah ciri khas data
hitungan yang cocok dimodelkan dengan regresi Poisson.
p1 <- ggplot(rd, aes(x = trips)) +
geom_histogram(
bins = 40, fill = "#2f7f73", color = "white", alpha = 0.85
) +
labs(
title = "Distribusi Jumlah Kunjungan (trips)",
subtitle = "Data hitungan dengan ekor kanan yang panjang",
x = "Jumlah Kunjungan",
y = "Frekuensi"
) +
theme_poisson()
p2 <- ggplot(rd, aes(x = trips)) +
geom_histogram(
aes(y = after_stat(density)),
bins = 40, fill = "#5568b8", color = "white", alpha = 0.85
) +
geom_density(color = "#d18b2f", linewidth = 1.2) +
labs(
title = "Densitas Jumlah Kunjungan",
subtitle = "Kurva densitas menunjukkan right-skewed distribution",
x = "Jumlah Kunjungan",
y = "Densitas"
) +
theme_poisson()
p1 + p2
## Rata-rata trips : 2.2443## Variansi trips : 39.5952## Rasio var/mean : 17.6425Perhatikan bahwa variansi trips jauh lebih besar dari
rata-ratanya. Ini merupakan indikasi awal adanya
overdispersi, yang akan diuji secara formal pada tahap
selanjutnya.
1.2.4 Visualisasi Hubungan Prediktor dengan Respon
p_qual <- ggplot(rd, aes(x = quality, y = trips)) +
geom_point(alpha = 0.35, color = "#2f7f73", size = 1.8) +
geom_smooth(method = "loess", se = TRUE, color = "#d18b2f",
fill = "#d18b2f", alpha = 0.15) +
labs(title = "Kualitas Air vs Kunjungan",
x = "quality", y = "trips") +
theme_poisson()
p_inc <- ggplot(rd, aes(x = income, y = trips)) +
geom_point(alpha = 0.35, color = "#5568b8", size = 1.8) +
geom_smooth(method = "loess", se = TRUE, color = "#d18b2f",
fill = "#d18b2f", alpha = 0.15) +
labs(title = "Pendapatan vs Kunjungan",
x = "income", y = "trips") +
theme_poisson()
p_costC <- ggplot(rd, aes(x = costC, y = trips)) +
geom_point(alpha = 0.35, color = "#b84f5a", size = 1.8) +
geom_smooth(method = "loess", se = TRUE, color = "#d18b2f",
fill = "#d18b2f", alpha = 0.15) +
labs(title = "Biaya Perjalanan (costC) vs Kunjungan",
x = "costC", y = "trips") +
theme_poisson()
p_ski <- ggplot(rd, aes(x = ski, y = trips, fill = ski)) +
geom_boxplot(alpha = 0.8, outlier.alpha = 0.3) +
scale_fill_manual(values = c("#2f7f73", "#5568b8")) +
labs(title = "Kebiasaan Ski Air vs Kunjungan",
x = "ski", y = "trips") +
theme_poisson() +
theme(legend.position = "none")
(p_qual + p_inc) / (p_costC + p_ski)
1.3 Pemodelan Regresi Poisson
1.3.1 Model Poisson Penuh
Pada tahap ini dibangun model regresi Poisson dengan seluruh variabel prediktor yang tersedia dalam dataset.
model_poisson <- glm(
trips ~ quality + ski + income + userfee + costC + costS + costH,
data = rd,
family = poisson(link = "log")
)
summary(model_poisson)##
## Call:
## glm(formula = trips ~ quality + ski + income + userfee + costC +
## costS + costH, family = poisson(link = "log"), data = rd)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.264993 0.093722 2.827 0.00469 **
## quality 0.471726 0.017091 27.602 < 2e-16 ***
## skiyes 0.418214 0.057190 7.313 2.62e-13 ***
## income -0.111323 0.019588 -5.683 1.32e-08 ***
## userfeeyes 0.898165 0.078985 11.371 < 2e-16 ***
## costC -0.003430 0.003118 -1.100 0.27131
## costS -0.042536 0.001670 -25.467 < 2e-16 ***
## costH 0.036134 0.002710 13.335 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 4849.7 on 658 degrees of freedom
## Residual deviance: 2305.8 on 651 degrees of freedom
## AIC: 3074.9
##
## Number of Fisher Scoring iterations: 7tidy_penuh <- broom::tidy(model_poisson, conf.int = TRUE) %>%
mutate(
IRR = exp(estimate),
conf.low = exp(conf.low),
conf.high = exp(conf.high),
across(where(is.numeric), ~ round(.x, 5))
) %>%
dplyr::rename(
"Variabel" = term,
"Koefisien" = estimate,
"Std. Error" = std.error,
"z-value" = statistic,
"p-value" = p.value,
"CI Lower" = conf.low,
"CI Upper" = conf.high
)
kable(tidy_penuh,
caption = "Koefisien model Poisson penuh beserta IRR dan interval kepercayaan 95%") %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE
) %>%
column_spec(1, bold = TRUE, color = "#2f7f73") %>%
column_spec(5, bold = TRUE,
color = ifelse(tidy_penuh$`p-value` < 0.05, "#b84f5a", "#64748b"))| Variabel | Koefisien | Std. Error | z-value | p-value | CI Lower | CI Upper | IRR |
|---|---|---|---|---|---|---|---|
| (Intercept) | 0.26499 | 0.09372 | 2.82743 | 0.00469 | 1.08355 | 1.56464 | 1.30342 |
| quality | 0.47173 | 0.01709 | 27.60161 | 0.00000 | 1.55025 | 1.65769 | 1.60276 |
| skiyes | 0.41821 | 0.05719 | 7.31268 | 0.00000 | 1.35805 | 1.69939 | 1.51925 |
| income | -0.11132 | 0.01959 | -5.68311 | 0.00000 | 0.86064 | 0.92933 | 0.89465 |
| userfeeyes | 0.89817 | 0.07899 | 11.37132 | 0.00000 | 2.09699 | 2.85851 | 2.45509 |
| costC | -0.00343 | 0.00312 | -1.10005 | 0.27131 | 0.99057 | 1.00276 | 0.99658 |
| costS | -0.04254 | 0.00167 | -25.46671 | 0.00000 | 0.95526 | 0.96154 | 0.95836 |
| costH | 0.03613 | 0.00271 | 13.33528 | 0.00000 | 1.03114 | 1.04215 | 1.03679 |
1.3.2 Model Poisson Tereduksi
Untuk memperoleh model yang lebih parsimoni, variabel dengan p-value > 0,05 dievaluasi. Di sini dibangun model yang hanya menyertakan variabel signifikan.
model_red <- glm(
trips ~ quality + ski + userfee + costC + costS,
data = rd,
family = poisson(link = "log")
)
summary(model_red)##
## Call:
## glm(formula = trips ~ quality + ski + userfee + costC + costS,
## family = poisson(link = "log"), data = rd)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.030140 0.073763 0.409 0.683
## quality 0.477016 0.016331 29.209 < 2e-16 ***
## skiyes 0.229016 0.054221 4.224 2.4e-05 ***
## userfeeyes 1.006180 0.078621 12.798 < 2e-16 ***
## costC 0.030155 0.001600 18.849 < 2e-16 ***
## costS -0.040533 0.001549 -26.175 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 4849.7 on 658 degrees of freedom
## Residual deviance: 2488.9 on 653 degrees of freedom
## AIC: 3254
##
## Number of Fisher Scoring iterations: 8tidy_red <- broom::tidy(model_red, conf.int = TRUE) %>%
mutate(
IRR = exp(estimate),
conf.low = exp(conf.low),
conf.high = exp(conf.high),
across(where(is.numeric), ~ round(.x, 5))
) %>%
dplyr::rename(
"Variabel" = term,
"Koefisien" = estimate,
"Std. Error" = std.error,
"z-value" = statistic,
"p-value" = p.value,
"CI Lower" = conf.low,
"CI Upper" = conf.high
)
kable(tidy_red,
caption = "Koefisien model Poisson tereduksi beserta IRR") %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE
) %>%
column_spec(1, bold = TRUE, color = "#2f7f73") %>%
column_spec(5, bold = TRUE,
color = ifelse(tidy_red$`p-value` < 0.05, "#b84f5a", "#64748b"))| Variabel | Koefisien | Std. Error | z-value | p-value | CI Lower | CI Upper | IRR |
|---|---|---|---|---|---|---|---|
| (Intercept) | 0.03014 | 0.07376 | 0.40861 | 0.68283 | 0.89072 | 1.18941 | 1.03060 |
| quality | 0.47702 | 0.01633 | 29.20908 | 0.00000 | 1.56084 | 1.66405 | 1.61126 |
| skiyes | 0.22902 | 0.05422 | 4.22379 | 0.00002 | 1.13045 | 1.39822 | 1.25736 |
| userfeeyes | 1.00618 | 0.07862 | 12.79790 | 0.00000 | 2.33778 | 3.18218 | 2.73513 |
| costC | 0.03016 | 0.00160 | 18.84891 | 0.00000 | 1.02734 | 1.03380 | 1.03061 |
| costS | -0.04053 | 0.00155 | -26.17504 | 0.00000 | 0.95740 | 0.96323 | 0.96028 |
1.4 Uji Asumsi dan Diagnostik Model
1.4.1 Uji Overdispersi
Overdispersi terjadi ketika variansi data jauh melebihi rata-ratanya.
Pada regresi Poisson, overdispersi dapat menyebabkan standard error yang
terlalu kecil sehingga menghasilkan kesimpulan yang keliru. Uji formal
dilakukan dengan fungsi dispersiontest() dari package
AER.
Hipotesis:
- \(H_0\): tidak ada overdispersi (model Poisson sudah cukup).
- \(H_1\): terdapat overdispersi.
disp_test <- AER::dispersiontest(model_red)
disp_val <- round(disp_test$estimate, 4)
disp_p <- signif(disp_test$p.value, 4)
print(disp_test)##
## Overdispersion test
##
## data: model_red
## z = 2.3707, p-value = 0.008878
## alternative hypothesis: true dispersion is greater than 1
## sample estimates:
## dispersion
## 7.504912fit_pois <- glm(trips ~ quality + income + costC + costS + costH,
data = rd, family = poisson)
dispersion_pois <- sum(residuals(fit_pois, type = "pearson")^2) / df.residual(fit_pois)
tibble::tibble(
`Dispersion Pearson` = round(dispersion_pois, 3),
`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"
)
) %>%
knitr::kable(
caption = "Indikasi overdispersion pada model Poisson (Data Recreation Demand)",
align = "cl"
) %>%
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE
)| Dispersion Pearson | Interpretasi ringkas |
|---|---|
| 6.483 | Ada indikasi overdispersion kuat |
Nilai statistik dispersi sebesar 7.5049 (jauh di atas 1) dan p-value = 0.008878 menunjukkan penolakan \(H_0\). Artinya, terdapat overdispersi yang signifikan pada data ini. Kondisi ini perlu ditangani, misalnya dengan model Negative Binomial atau Quasi-Poisson.
1.4.2 Perbandingan Deviance dan AIC
gof_tab <- tibble::tibble(
Model = c("Poisson Penuh", "Poisson Tereduksi"),
`Null Deviance` = c(round(model_poisson$null.deviance, 2),
round(model_red$null.deviance, 2)),
`Residual Deviance` = c(round(deviance(model_poisson), 2),
round(deviance(model_red), 2)),
`df Residual` = c(model_poisson$df.residual,
model_red$df.residual),
AIC = c(round(AIC(model_poisson), 2),
round(AIC(model_red), 2))
)
kable(gof_tab,
caption = "Perbandingan kebaikan model (Goodness of Fit)") %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE
)| Model | Null Deviance | Residual Deviance | df Residual | AIC |
|---|---|---|---|---|
| Poisson Penuh | 4849.69 | 2305.79 | 651 | 3074.86 |
| Poisson Tereduksi | 4849.69 | 2488.95 | 653 | 3254.02 |
1.4.3 Uji Likelihood Ratio
## Analysis of Deviance Table
##
## Model 1: trips ~ quality + ski + userfee + costC + costS
## Model 2: trips ~ quality + ski + income + userfee + costC + costS + costH
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 653 2488.9
## 2 651 2305.8 2 183.16 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Uji Likelihood Ratio (Likelihood Ratio Test) membandingkan model tereduksi dengan model penuh. Jika p-value > 0,05, model tereduksi tidak berbeda signifikan dari model penuh dan lebih disarankan karena lebih parsimoni.
1.4.4 Plot Diagnostik
par(mfrow = c(2, 2), bg = "white")
plot(model_poisson, which = 1:4,
col = "#2f7f73", pch = 16, cex = 0.65)
1.4.5 Residual Pearson vs Fitted Values
aug_penuh <- broom::augment(model_poisson, type.residuals = "pearson")
ggplot(aug_penuh, aes(x = .fitted, y = .resid)) +
geom_point(alpha = 0.4, color = "#5568b8", size = 1.6) +
geom_hline(yintercept = 0, color = "#b84f5a", linewidth = 1) +
geom_smooth(method = "loess", se = FALSE, color = "#d18b2f",
linewidth = 1.1, linetype = "dashed") +
labs(
title = "Residual Pearson vs Nilai Fitted",
subtitle = "Pola yang ideal adalah residual tersebar acak di sekitar nol",
x = "Fitted Values (log scale)",
y = "Residual Pearson"
) +
theme_poisson()
1.5 Penanganan Overdispersi
1.5.1 Model Quasi-Poisson
Quasi-Poisson menggunakan estimasi dispersi dari data, sehingga standard error yang dihasilkan lebih konservatif dan inference menjadi lebih valid. Model ini tidak mengubah estimasi koefisien, hanya menyesuaikan standard error.
model_quasi <- glm(
trips ~ quality + ski + userfee + costC + costS,
data = rd,
family = quasipoisson(link = "log")
)
summary(model_quasi)##
## Call:
## glm(formula = trips ~ quality + ski + userfee + costC + costS,
## family = quasipoisson(link = "log"), data = rd)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.030140 0.197951 0.152 0.879
## quality 0.477016 0.043826 10.884 < 2e-16 ***
## skiyes 0.229016 0.145507 1.574 0.116
## userfeeyes 1.006180 0.210988 4.769 2.29e-06 ***
## costC 0.030155 0.004293 7.024 5.44e-12 ***
## costS -0.040533 0.004156 -9.754 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 7.201791)
##
## Null deviance: 4849.7 on 658 degrees of freedom
## Residual deviance: 2488.9 on 653 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 81.5.2 Model Negative Binomial
Model Negative Binomial merupakan alternatif yang lebih fleksibel dibanding Quasi-Poisson karena secara eksplisit memodelkan overdispersi melalui parameter dispersi \(\theta\). Model ini cocok ketika variansi jauh melebihi rata-rata.
model_nb <- MASS::glm.nb(
trips ~ quality + ski + userfee + costC + costS,
data = rd
)
summary(model_nb)##
## Call:
## MASS::glm.nb(formula = trips ~ quality + ski + userfee + costC +
## costS, data = rd, init.theta = 0.6807956681, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.078739 0.166513 -6.478 9.27e-11 ***
## quality 0.725578 0.040652 17.848 < 2e-16 ***
## skiyes 0.570642 0.148381 3.846 0.00012 ***
## userfeeyes 0.642382 0.362674 1.771 0.07652 .
## costC 0.095381 0.006606 14.438 < 2e-16 ***
## costS -0.102227 0.006712 -15.231 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.6808) family taken to be 1)
##
## Null deviance: 1195.70 on 658 degrees of freedom
## Residual deviance: 421.26 on 653 degrees of freedom
## AIC: 1677.2
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.6808
## Std. Err.: 0.0680
##
## 2 x log-likelihood: -1663.15901.5.3 Perbandingan Tiga Model
koef_poisson <- broom::tidy(model_red, conf.int = TRUE) %>% mutate(model = "Poisson")
koef_quasi <- broom::tidy(model_quasi, conf.int = TRUE) %>% mutate(model = "Quasi-Poisson")
koef_nb <- broom::tidy(model_nb, conf.int = TRUE) %>% mutate(model = "Negative Binomial")
koef_all <- bind_rows(koef_poisson, koef_quasi, koef_nb) %>%
filter(term != "(Intercept)") %>%
mutate(
IRR = exp(estimate),
IRR.low = exp(conf.low),
IRR.high = exp(conf.high)
)
ggplot(koef_all, aes(x = IRR, y = term, color = model, shape = model)) +
geom_vline(xintercept = 1, linetype = "dashed", color = "#64748b") +
geom_pointrange(
aes(xmin = IRR.low, xmax = IRR.high),
position = position_dodge(width = 0.55),
size = 0.7
) +
scale_color_manual(values = c("#2f7f73", "#5568b8", "#d18b2f")) +
labs(
title = "Perbandingan IRR: Poisson vs Quasi-Poisson vs Negative Binomial",
subtitle = "Garis putus-putus pada IRR = 1 (tidak ada efek)",
x = "IRR (exp[koefisien])",
y = "Variabel",
color = "Model",
shape = "Model"
) +
theme_poisson()
comp_tab <- tibble::tibble(
Model = c("Poisson", "Quasi-Poisson", "Negative Binomial"),
AIC = c(round(AIC(model_red), 2), "—", round(AIC(model_nb), 2)),
`Residual Deviance` = c(
round(deviance(model_red), 2),
round(deviance(model_quasi), 2),
round(deviance(model_nb), 2)
),
`Keterangan` = c(
"Standar; asumsi ekuidispersi",
"SE disesuaikan; koefisien sama dengan Poisson",
"Memodelkan overdispersi secara eksplisit; AIC lebih kecil"
)
)
kable(comp_tab,
caption = "Perbandingan ringkas tiga model untuk data count") %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE
)| Model | AIC | Residual Deviance | Keterangan |
|---|---|---|---|
| Poisson | 3254.02 | 2488.95 | Standar; asumsi ekuidispersi |
| Quasi-Poisson | — | 2488.95 | SE disesuaikan; koefisien sama dengan Poisson |
| Negative Binomial | 1677.16 | 421.26 | Memodelkan overdispersi secara eksplisit; AIC lebih kecil |
1.6 Interpretasi Koefisien Model Terbaik
Model terpilih
Mengingat adanya overdispersi yang signifikan, model Negative Binomial dipilih sebagai model terbaik. Model ini memiliki AIC yang lebih rendah dan mampu mengatasi overdispersi secara lebih tepat dibanding Poisson biasa.
irr_nb <- broom::tidy(model_nb, conf.int = TRUE, exponentiate = TRUE) %>%
dplyr::rename(
"Variabel" = term,
"IRR" = estimate,
"Std. Error" = std.error,
"z-value" = statistic,
"p-value" = p.value,
"CI Lower 95%" = conf.low,
"CI Upper 95%" = conf.high
) %>%
mutate(across(where(is.numeric), ~ round(.x, 5)))
kable(irr_nb,
caption = "Incidence Rate Ratio (IRR) model Negative Binomial") %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE
) %>%
column_spec(1, bold = TRUE, color = "#2f7f73") %>%
column_spec(5, bold = TRUE,
color = ifelse(irr_nb$`p-value` < 0.05, "#b84f5a", "#64748b"))| Variabel | IRR | Std. Error | z-value | p-value | CI Lower 95% | CI Upper 95% |
|---|---|---|---|---|---|---|
| (Intercept) | 0.34002 | 0.16651 | -6.47840 | 0.00000 | 0.24834 | 0.46253 |
| quality | 2.06593 | 0.04065 | 17.84850 | 0.00000 | 1.89343 | 2.26134 |
| skiyes | 1.76940 | 0.14838 | 3.84579 | 0.00012 | 1.31906 | 2.38232 |
| userfeeyes | 1.90100 | 0.36267 | 1.77124 | 0.07652 | 0.96917 | 4.20844 |
| costC | 1.10008 | 0.00661 | 14.43846 | 0.00000 | 1.08193 | 1.11913 |
| costS | 0.90282 | 0.00671 | -15.23126 | 0.00000 | 0.88811 | 0.91730 |
1.6.1 Interpretasi Per Variabel
quality — Kualitas Air
Koefisien variabel quality sebesar 0.7256 dengan IRR = exp(β) sebesar 2.0659. Hal ini menunjukkan bahwa setiap kenaikan satu unit quality meningkatkan rate atau rata-rata jumlah kunjungan sebesar 106.59%, dengan asumsi variabel lain dan exposure konstan.
ski — Kebiasaan Ski Air
Koefisien variabel ski (kategori yes dibanding no) sebesar 0.5706 dengan IRR = exp(β) sebesar 1.7694. Hal ini menunjukkan bahwa individu yang juga melakukan ski air meningkatkan rate atau rata-rata jumlah kunjungan ke lokasi rekreasi sebesar 76.94% dibandingkan yang tidak melakukan ski air, dengan asumsi variabel lain dan exposure konstan.
userfee — Pembayaran User Fee
Koefisien variabel userfee (kategori yes dibanding no) sebesar 0.6424 dengan IRR = exp(β) sebesar 1.901. Hal ini menunjukkan bahwa individu yang membayar user fee meningkatkan rate atau rata-rata jumlah kunjungan sebesar 90.1% dibandingkan yang tidak membayar, dengan asumsi variabel lain dan exposure konstan.
costC — Biaya Perjalanan ke Lokasi Tujuan
Koefisien variabel costC sebesar 0.0954 dengan IRR = exp(β) sebesar 1.1001. Hal ini menunjukkan bahwa setiap kenaikan satu unit costC meningkatkan rate atau rata-rata jumlah kunjungan sebesar 10.01%, dengan asumsi variabel lain dan exposure konstan. Pola ini masuk akal secara ekonomi, semakin mahal biaya ke tujuan utama, semakin sedikit kunjungan yang dilakukan.
costS — Biaya Perjalanan ke Lokasi Substitusi
Koefisien variabel costS sebesar -0.1022 dengan IRR = exp(β) sebesar 0.9028. Hal ini menunjukkan bahwa setiap kenaikan satu unit costS menurunkan rate atau rata-rata jumlah kunjungan ke lokasi utama sebesar 9.72%, dengan asumsi variabel lain dan exposure konstan. Ini konsisten dengan teori substitusi, ketika lokasi alternatif semakin mahal, permintaan terhadap lokasi utama naik.
1.7 Visualisasi Hasil Model
1.7.1 Plot IRR dengan Confidence Interval
irr_plot <- broom::tidy(model_nb, conf.int = TRUE, exponentiate = TRUE) %>%
filter(term != "(Intercept)") %>%
mutate(
signif = ifelse(p.value < 0.05, "Signifikan", "Tidak Signifikan"),
term = recode(term,
"quality" = "Kualitas Air\n(quality)",
"skiyes" = "Ski Air: Yes\n(ski)",
"userfeeyes" = "User Fee: Yes\n(userfee)",
"costC" = "Biaya ke Tujuan\n(costC)",
"costS" = "Biaya ke Substitusi\n(costS)"
)
)
ggplot(irr_plot, aes(x = estimate, y = reorder(term, estimate), color = signif)) +
geom_vline(xintercept = 1, linetype = "dashed", color = "#64748b", linewidth = 0.8) +
geom_pointrange(
aes(xmin = conf.low, xmax = conf.high),
size = 0.8, linewidth = 1.1
) +
scale_color_manual(values = c("Signifikan" = "#2f7f73",
"Tidak Signifikan" = "#b84f5a")) +
labs(
title = "Incidence Rate Ratio — Model Negative Binomial",
subtitle = "Titik = IRR; garis = interval kepercayaan 95%. IRR > 1 menaikkan kunjungan.",
x = "IRR",
y = NULL,
color = "Signifikansi (α = 0.05)"
) +
theme_poisson()
1.7.2 Prediksi Rata-rata Kunjungan berdasarkan Kualitas Air
pred_quality <- data.frame(
quality = seq(min(rd$quality), max(rd$quality), length.out = 100),
ski = "no",
userfee = "no",
costC = mean(rd$costC),
costS = mean(rd$costS)
)
pred_quality$trips_pred <- predict(model_nb, newdata = pred_quality, type = "response")
ggplot(pred_quality, aes(x = quality, y = trips_pred)) +
geom_line(color = "#2f7f73", linewidth = 1.4) +
geom_ribbon(
aes(ymin = trips_pred * 0.85, ymax = trips_pred * 1.15),
fill = "#2f7f73", alpha = 0.12
) +
labs(
title = "Prediksi Rata-rata Kunjungan berdasarkan Kualitas Air",
subtitle = "Variabel lain ditetapkan pada nilai rata-rata (ski=no, userfee=no)",
x = "Kualitas Air (quality)",
y = "Prediksi Jumlah Kunjungan (trips)"
) +
theme_poisson()
1.7.3 Prediksi Berdasarkan Biaya Perjalanan (costC)
pred_cost <- data.frame(
costC = seq(min(rd$costC), max(rd$costC), length.out = 100),
quality = mean(rd$quality),
ski = "no",
userfee = "no",
costS = mean(rd$costS)
)
pred_cost$trips_pred <- predict(model_nb, newdata = pred_cost, type = "response")
ggplot(pred_cost, aes(x = costC, y = trips_pred)) +
geom_line(color = "#b84f5a", linewidth = 1.4) +
geom_ribbon(
aes(ymin = trips_pred * 0.85, ymax = trips_pred * 1.15),
fill = "#b84f5a", alpha = 0.12
) +
labs(
title = "Prediksi Rata-rata Kunjungan berdasarkan Biaya ke Lokasi Tujuan",
subtitle = "Semakin tinggi biaya perjalanan, semakin rendah rata-rata kunjungan",
x = "Biaya Perjalanan ke Lokasi (costC)",
y = "Prediksi Jumlah Kunjungan (trips)"
) +
theme_poisson()
1.7.4 Nilai Aktual vs Prediksi
rd$pred_nb <- predict(model_nb, type = "response")
ggplot(rd, aes(x = pred_nb, y = trips)) +
geom_point(alpha = 0.3, color = "#5568b8", size = 1.6) +
geom_abline(slope = 1, intercept = 0,
color = "#b84f5a", linewidth = 1.2, linetype = "dashed") +
coord_cartesian(xlim = c(0, 20), ylim = c(0, 50)) +
labs(
title = "Nilai Aktual vs Prediksi — Model Negative Binomial",
subtitle = "Garis putus-putus merepresentasikan prediksi sempurna (aktual = prediksi)",
x = "Prediksi (fitted values)",
y = "Aktual (trips)"
) +
theme_poisson()
1.8 Rangkuman Model dan Narasi Hasil
1.8.1 Persamaan Model Terpilih
Berdasarkan seluruh tahapan analisis, model terpilih adalah Negative Binomial dengan persamaan:
\[ \log(\hat{\mu}_i) = \hat{\beta}_0 + \hat{\beta}_1 \cdot \text{quality}_i + \hat{\beta}_2 \cdot \mathbb{1}[\text{ski}=\text{yes}]_i + \hat{\beta}_3 \cdot \mathbb{1}[\text{userfee}=\text{yes}]_i + \hat{\beta}_4 \cdot \text{costC}_i + \hat{\beta}_5 \cdot \text{costS}_i \]
dimana \(\hat{\mu}_i\) adalah taksiran rata-rata jumlah kunjungan untuk individu ke-\(i\).
1.8.2 Hasil
Model regresi Negative Binomial digunakan untuk menganalisis
faktor-faktor yang memengaruhi intensitas kunjungan ke lokasi rekreasi
perairan, dengan variabel respon berupa jumlah kunjungan
(trips). Pemilihan model Negative Binomial didasarkan pada
hasil uji overdispersi yang signifikan, di mana variansi data jauh
melampaui rata-ratanya.
Hasil estimasi menunjukkan bahwa kualitas air berpengaruh positif dan signifikan terhadap jumlah kunjungan. Semakin baik kualitas air di lokasi tujuan, semakin tinggi intensitas kunjungan yang dilakukan. Individu yang juga melakukan ski air (ski = yes) cenderung mengunjungi lokasi rekreasi lebih sering dibandingkan yang tidak, mencerminkan pola bahwa penggemar olahraga air memiliki minat rekreasi yang lebih tinggi secara umum.
Di sisi lain, biaya perjalanan ke lokasi tujuan
(costC) berpengaruh negatif — sesuai dengan hukum
permintaan dasar — sementara biaya perjalanan ke lokasi
substitusi (costS) berpengaruh positif, yang
bermakna bahwa ketika lokasi alternatif menjadi lebih mahal, permintaan
terhadap lokasi utama justru meningkat. Variabel user
fee juga menunjukkan pengaruh yang signifikan terhadap pola
kunjungan.
Secara keseluruhan, model ini memberikan gambaran yang cukup komprehensif mengenai determinan permintaan rekreasi perairan, dengan implikasi penting bagi pengelolaan kawasan wisata dan kebijakan penetapan harga.
2 Regresi Logistik Ordinal
Regresi logistik ordinal digunakan ketika variabel respons bersifat kategorik dengan urutan yang bermakna. Model ini bekerja dengan cara memodelkan peluang kumulatif, yaitu peluang bahwa respons berada pada kategori tertentu atau di bawahnya.
Model Regresi Logistik Ordinal (Proportional Odds Model)
Secara formal, model ini dinyatakan sebagai:
\[ \text{logit}\left[P(Y \leq j)\right] = \alpha_j - \boldsymbol{\beta}^T\mathbf{x}, \quad j = 1, 2, \ldots, J-1 \]
Di mana:
- \(P(Y \leq j)\) adalah peluang kumulatif bahwa \(Y\) berada pada kategori ke-\(j\) atau lebih rendah
- \(\alpha_j\) adalah intercept (threshold) untuk setiap batas kategori, dengan \(\alpha_1 < \alpha_2 < \ldots < \alpha_{J-1}\)
- \(\boldsymbol{\beta}\) adalah vektor koefisien regresi untuk variabel prediktor \(\mathbf{x}\)
- Tanda negatif pada \(\boldsymbol{\beta}^T\mathbf{x}\) mengikuti
konvensi paket
MASS::polrdi R
Interpretasi koefisien dilakukan melalui Odds Ratio (OR):
\[OR = e^{\beta}\]
Jika \(OR > 1\), maka kenaikan satu unit pada prediktor meningkatkan odds untuk berada di kategori yang lebih tinggi (lebih miskin). Sebaliknya, \(OR < 1\) berarti prediktor tersebut bersifat protektif — menurunkan peluang masuk ke kategori kemiskinan yang lebih parah.
2.1 Sumber Data dan Deskripsi Variabel
Sumber Data
Dataset yang digunakan adalah Multidimensional Poverty Index (MPI) Dataset yang dapat diakses melalui repositori publik Harvard Dataverse mencakup 394 observasi. Dataset ini bersumber dari studi analisis level-mikro mengenai penentu kemiskinan multidimensi di wilayah Administrasi Kota Debre Tabor, Ethiopia yang merekam kondisi kemiskinan multidimensi dari berbagai aspek kehidupan, mulai dari indikator pendidikan, kesehatan, standar hunian, akses energi, hingga karakteristik sosial-ekonomi kepala rumah tangga .
| Variabel | Jenis | Keterangan |
|---|---|---|
| MDP | Respons (Ordinal) | Status kemiskinan multidimensi: Tidak Miskin, Rentan, Miskin Moderat, Sangat Miskin |
| gender | Kategorik | Jenis kelamin kepala rumah tangga (male/female) |
| age | Kontinu | Usia kepala rumah tangga (tahun) |
| educ_level | Kategorik | Tingkat pendidikan kepala rumah tangga |
| marital_status | Kategorik | Status pernikahan kepala rumah tangga |
| fam_size | Diskrit | Jumlah anggota keluarga |
| depratio | Kontinu | Rasio anggota tidak produktif terhadap produktif |
| occupation | Kategorik | Jenis pekerjaan kepala rumah tangga |
| credit | Kategorik | Akses kredit formal (ya/tidak) |
| house_ownership | Kategorik | Kepemilikan rumah (milik sendiri/tidak) |
| female_income | Kategorik | Partisipasi pendapatan perempuan dalam rumah tangga |
| seving_ratio | Kontinu | Rasio tabungan terhadap pendapatan |
Variabel-variabel seperti hhd_id, income,
lnincome, income_level,
deprivation_score dan semua kolom indikator biner
individual tidak diikutsertakan untuk menghindari
multikolinearitas dan overfitting.
2.2 Persiapan Data
library(readxl)
df_raw <- read_excel("C:/Users/Lenovo/Downloads/Copy of MPI_dataset.xlsx")
cat("Jumlah observasi:", nrow(df_raw), "\n")## Jumlah observasi: 394## Jumlah variabel : 352.2.1 Seleksi dan Transformasi Variabel
df <- df_raw |>
dplyr::select(MDP, gender, age, educ_level, marital_status,
fam_size, depratio, occupation, credit,
house_ownership, female_income, seving_ratio) |>
na.omit()
cat("Data setelah seleksi variabel:", nrow(df), "observasi,", ncol(df), "variabel\n")## Data setelah seleksi variabel: 394 observasi, 12 variabelmdp_levels <- c(
"not multidimensional poor",
"vulnerable to multidimensional poverty",
"moderately multidimensional poor",
"extremely multidimensional poor"
)
df <- df |>
mutate(
MDP = factor(MDP, levels = mdp_levels, ordered = TRUE),
gender = factor(gender, levels = c("male", "female")),
educ_level = factor(educ_level,
levels = c("illiterate", "primary education",
"secondary education", "higher edcuation")),
marital_status = factor(marital_status,
levels = c("single", "married", "divorced", "widowed")),
occupation = factor(occupation,
levels = c("unemployed", "causal wage worker",
"owner of private bisness or self employed",
"salaried")),
credit = factor(credit, levels = c("do not get credit", "get credit")),
house_ownership = factor(house_ownership,
levels = c("not own a house", "own house")),
female_income = factor(female_income,
levels = c("no female participate", "female participate"))
)
glimpse(df)## Rows: 394
## Columns: 12
## $ MDP <ord> not multidimensional poor, moderately multidimensional…
## $ gender <fct> male, male, male, male, male, female, male, male, fema…
## $ age <dbl> 45, 60, 25, 38, 34, 33, 39, 29, 27, 29, 40, 41, 35, 35…
## $ educ_level <fct> higher edcuation, illiterate, secondary education, hig…
## $ marital_status <fct> married, widowed, married, married, married, married, …
## $ fam_size <dbl> 4, 4, 3, 7, 7, 4, 6, 1, 1, 1, 8, 5, 5, 6, 5, 5, 4, 3, …
## $ depratio <dbl> 1.000000, 3.000000, 0.000000, 1.333333, 0.400000, 1.00…
## $ occupation <fct> salaried, owner of private bisness or self employed, o…
## $ credit <fct> do not get credit, do not get credit, get credit, do n…
## $ house_ownership <fct> own house, own house, not own a house, not own a house…
## $ female_income <fct> female participate, female participate, female partici…
## $ seving_ratio <dbl> 0.6000000, 0.0675000, 0.4000000, 0.1166667, 0.7000000,…2.3 Analisis Deskriptif
2.3.1 Distribusi Variabel Respons (MDP)
mdp_label <- c(
"not multidimensional poor" = "Tidak Miskin",
"vulnerable to multidimensional poverty" = "Rentan",
"moderately multidimensional poor" = "Miskin Moderat",
"extremely multidimensional poor" = "Sangat Miskin"
)
df_plot <- df |>
count(MDP) |>
mutate(
Label = mdp_label[as.character(MDP)],
Label = factor(Label, levels = c("Tidak Miskin", "Rentan",
"Miskin Moderat", "Sangat Miskin")),
persen = n / sum(n) * 100
)
ggplot(df_plot, aes(x = Label, y = n, fill = Label)) +
geom_col(width = 0.65, show.legend = FALSE) +
geom_text(aes(label = paste0(n, "\n(", round(persen, 1), "%)")),
vjust = -0.4, size = 3.8, fontface = "bold") +
scale_fill_manual(values = c("#2ecc71", "#f39c12", "#e67e22", "#e74c3c")) +
scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +
labs(
title = "Distribusi Status Kemiskinan Multidimensi",
subtitle = "Dataset MPI — 394 Rumah Tangga",
x = "Kategori MDP",
y = "Jumlah Rumah Tangga"
) +
theme_adk() +
theme(axis.text.x = element_text(face = "bold"))
Dari 394 rumah tangga dalam dataset, sebagian besar terklasifikasi sebagai sangat miskin multidimensi (126 RT, 32%), diikuti oleh miskin moderat (108 RT, 27,4%). Hanya sekitar 27,2% rumah tangga yang tidak miskin secara multidimensi, dan 13,5% berada dalam kondisi rentan. Distribusi ini cukup mengkhawatirkan — lebih dari separuh rumah tangga berada di dua kategori kemiskinan terberat.
2.3.2 Distribusi Variabel Kategorik terhadap MDP
plot_stacked <- function(var, label_var, judul) {
df |>
mutate(MDP_label = recode(as.character(MDP),
"not multidimensional poor" = "Tidak Miskin",
"vulnerable to multidimensional poverty" = "Rentan",
"moderately multidimensional poor" = "Miskin Moderat",
"extremely multidimensional poor" = "Sangat Miskin"),
MDP_label = factor(MDP_label,
levels = c("Tidak Miskin", "Rentan",
"Miskin Moderat", "Sangat Miskin"))
) |>
dplyr::rename(Var = all_of(var)) |>
count(Var, MDP_label) |>
group_by(Var) |>
mutate(pct = n / sum(n) * 100) |>
ungroup() |>
ggplot(aes(x = fct_rev(Var), y = pct, fill = MDP_label)) +
geom_col(position = "stack", width = 0.6) +
coord_flip() +
scale_fill_manual(values = c("#2ecc71", "#f39c12", "#e67e22", "#e74c3c"),
name = "Status MDP") +
scale_y_continuous(labels = percent_format(scale = 1)) +
labs(title = judul, x = label_var, y = "Proporsi (%)") +
theme_adk() +
theme(legend.position = "bottom",
plot.title = element_text(size = 11))
}
p1 <- plot_stacked("educ_level", "Pendidikan", "Pendidikan vs MDP")
p2 <- plot_stacked("occupation", "Pekerjaan", "Pekerjaan vs MDP")
p3 <- plot_stacked("credit", "Akses Kredit", "Kredit vs MDP")
p4 <- plot_stacked("house_ownership","Kepemilikan Rumah", "Kepemilikan Rumah vs MDP")
p5 <- plot_stacked("female_income", "Pendapatan Wanita", "Perempuan Bekerja vs MDP")
p6 <- plot_stacked("gender", "Jenis Kelamin KRT", "Gender vs MDP")
gridExtra::grid.arrange(p1, p2, p3, p4, p5, p6, ncol = 2)
Dari visualisasi distribusi ini, terlihat beberapa pola yang cukup mencolok:
- Rumah tangga dengan kepala keluarga buta huruf hampir seluruhnya masuk kategori miskin parah, sedangkan mereka yang berpendidikan tinggi memiliki proporsi “tidak miskin” yang jauh lebih besar.
- Kelompok pengangguran dan pekerja harian memiliki proporsi kemiskinan parah yang lebih tinggi dibandingkan pekerja formal atau wirausaha.
- Rumah tangga yang tidak memiliki akses kredit cenderung lebih rentan, begitu pula yang tidak memiliki rumah sendiri.
2.3.3 Distribusi Variabel Kontinu
df_long_cont <- df |>
mutate(MDP_label = recode(as.character(MDP),
"not multidimensional poor" = "Tidak Miskin",
"vulnerable to multidimensional poverty" = "Rentan",
"moderately multidimensional poor" = "Miskin Moderat",
"extremely multidimensional poor" = "Sangat Miskin"),
MDP_label = factor(MDP_label,
levels = c("Tidak Miskin","Rentan",
"Miskin Moderat","Sangat Miskin"))
) |>
dplyr::select(MDP_label, age, fam_size, depratio, seving_ratio) |>
pivot_longer(-MDP_label, names_to = "variabel", values_to = "nilai") |>
mutate(variabel = recode(variabel,
age = "Usia",
fam_size = "Ukuran Keluarga",
depratio = "Rasio Ketergantungan",
seving_ratio = "Rasio Tabungan"
))
ggplot(df_long_cont, aes(x = MDP_label, y = nilai, fill = MDP_label)) +
geom_boxplot(alpha = 0.7, outlier.size = 1.2, width = 0.6) +
facet_wrap(~variabel, scales = "free_y", ncol = 2) +
scale_fill_manual(values = c("#2ecc71", "#f39c12", "#e67e22", "#e74c3c")) +
labs(
title = "Distribusi Variabel Kontinu Berdasarkan Kategori MDP",
x = NULL,
y = "Nilai"
) +
theme_adk() +
theme(axis.text.x = element_text(angle = 15, hjust = 1),
legend.position = "none")
Boxplot ini memperlihatkan bahwa rasio ketergantungan cenderung meningkat seiring dengan tingkat kemiskinan yang lebih berat — rumah tangga yang sangat miskin memiliki median rasio ketergantungan yang lebih tinggi. Sebaliknya, rasio tabungan menurun tajam pada kategori kemiskinan yang lebih parah, konsisten dengan ekspektasi ekonomi.
2.3.4 Statistik Deskriptif
## age fam_size depratio seving_ratio
## Min. :22.00 Min. :1.000 Min. :0.0000 Min. :0.0000
## 1st Qu.:29.00 1st Qu.:2.000 1st Qu.:0.0000 1st Qu.:0.0600
## Median :34.00 Median :4.000 Median :1.0000 Median :0.1167
## Mean :35.68 Mean :3.701 Mean :0.9964 Mean :0.1699
## 3rd Qu.:39.00 3rd Qu.:5.000 3rd Qu.:1.5000 3rd Qu.:0.2000
## Max. :65.00 Max. :8.000 Max. :6.0000 Max. :1.9039tab_cat <- df |>
dplyr::select(gender, educ_level, marital_status, occupation,
credit, house_ownership, female_income) |>
pivot_longer(everything(), names_to = "Variabel", values_to = "Kategori") |>
count(Variabel, Kategori) |>
group_by(Variabel) |>
mutate(Persen = paste0(round(n / sum(n) * 100, 1), "%")) |>
ungroup()
kable(tab_cat,
col.names = c("Variabel", "Kategori", "Frekuensi", "Persen"),
caption = "Distribusi Frekuensi Variabel Kategorik") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE) %>%
column_spec(1, bold = TRUE, color = "#2f7f73") %>%
collapse_rows(columns = 1, valign = "top")| Variabel | Kategori | Frekuensi | Persen |
|---|---|---|---|
| credit | do not get credit | 320 | 81.2% |
| get credit | 74 | 18.8% | |
| educ_level | illiterate | 4 | 1% |
| primary education | 30 | 7.6% | |
| secondary education | 94 | 23.9% | |
| higher edcuation | 266 | 67.5% | |
| female_income | no female participate | 204 | 51.8% |
| female participate | 190 | 48.2% | |
| gender | male | 284 | 72.1% |
| female | 110 | 27.9% | |
| house_ownership | not own a house | 238 | 60.4% |
| own house | 156 | 39.6% | |
| marital_status | single | 110 | 27.9% |
| married | 240 | 60.9% | |
| divorced | 36 | 9.1% | |
| widowed | 8 | 2% | |
| occupation | unemployed | 38 | 9.6% |
| causal wage worker | 74 | 18.8% | |
| owner of private bisness or self employed | 56 | 14.2% | |
| salaried | 226 | 57.4% |
2.4 Pemodelan Regresi Logistik Ordinal
2.4.1 Pembagian Data (Train-Test Split)
set.seed(42)
n <- nrow(df)
idx_train <- sample(seq_len(n), size = floor(0.8 * n))
df_train <- df[idx_train, ]
df_test <- df[-idx_train, ]
cat("Jumlah data training :", nrow(df_train), "\n")## Jumlah data training : 315## Jumlah data testing : 792.4.2 Model Ordinal (Full Model)
Model dibangun menggunakan fungsi polr() dari paket
MASS, dengan metode estimasi logistic
(proportional odds). Variabel respons MDP sudah
didefinisikan sebagai faktor terurut.
model_full_ord <- polr(
MDP ~ gender + age + educ_level + marital_status +
fam_size + depratio + occupation + credit +
house_ownership + female_income + seving_ratio,
data = df_train,
Hess = TRUE,
method = "logistic"
)
summary(model_full_ord)## Call:
## polr(formula = MDP ~ gender + age + educ_level + marital_status +
## fam_size + depratio + occupation + credit + house_ownership +
## female_income + seving_ratio, data = df_train, Hess = TRUE,
## method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## genderfemale 1.07388 0.31397 3.4203
## age 0.00852 0.01563 0.5451
## educ_levelprimary education 0.39964 1.34976 0.2961
## educ_levelsecondary education -0.42693 1.27057 -0.3360
## educ_levelhigher edcuation -0.69677 1.26277 -0.5518
## marital_statusmarried -1.90732 0.30203 -6.3150
## marital_statusdivorced -0.70254 0.46608 -1.5073
## marital_statuswidowed -1.45837 0.84380 -1.7283
## fam_size 0.11498 0.07913 1.4530
## depratio 0.15898 0.12521 1.2696
## occupationcausal wage worker 0.19490 0.51795 0.3763
## occupationowner of private bisness or self employed -0.59336 0.51270 -1.1573
## occupationsalaried -1.43454 0.45748 -3.1358
## creditget credit -1.12303 0.32954 -3.4079
## house_ownershipown house -0.82646 0.25844 -3.1978
## female_incomefemale participate -0.79709 0.28631 -2.7840
## seving_ratio -2.56173 0.92630 -2.7656
##
## Intercepts:
## Value
## not multidimensional poor|vulnerable to multidimensional poverty -4.3492
## vulnerable to multidimensional poverty|moderately multidimensional poor -3.4711
## moderately multidimensional poor|extremely multidimensional poor -1.5170
## Std. Error
## not multidimensional poor|vulnerable to multidimensional poverty 1.5152
## vulnerable to multidimensional poverty|moderately multidimensional poor 1.5077
## moderately multidimensional poor|extremely multidimensional poor 1.4967
## t value
## not multidimensional poor|vulnerable to multidimensional poverty -2.8704
## vulnerable to multidimensional poverty|moderately multidimensional poor -2.3023
## moderately multidimensional poor|extremely multidimensional poor -1.0136
##
## Residual Deviance: 628.0636
## AIC: 668.06362.4.3 Seleksi Model: Stepwise AIC
Untuk menghindari overfitting dan mendapatkan model yang lebih parsimonius, dilakukan seleksi variabel menggunakan prosedur stepwise berdasarkan kriteria AIC (Akaike Information Criterion).
## Call:
## polr(formula = MDP ~ gender + marital_status + fam_size + occupation +
## credit + house_ownership + female_income + seving_ratio,
## data = df_train, Hess = TRUE, method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## genderfemale 1.1004 0.3098 3.5521
## marital_statusmarried -1.8223 0.2954 -6.1694
## marital_statusdivorced -0.5217 0.4414 -1.1818
## marital_statuswidowed -0.9416 0.7115 -1.3234
## fam_size 0.2118 0.0653 3.2434
## occupationcausal wage worker 0.4046 0.5079 0.7968
## occupationowner of private bisness or self employed -0.5847 0.5123 -1.1412
## occupationsalaried -1.5569 0.4482 -3.4734
## creditget credit -0.9921 0.3171 -3.1285
## house_ownershipown house -0.7594 0.2484 -3.0573
## female_incomefemale participate -0.8937 0.2780 -3.2146
## seving_ratio -2.7076 0.9136 -2.9638
##
## Intercepts:
## Value
## not multidimensional poor|vulnerable to multidimensional poverty -3.8445
## vulnerable to multidimensional poverty|moderately multidimensional poor -2.9737
## moderately multidimensional poor|extremely multidimensional poor -1.0641
## Std. Error
## not multidimensional poor|vulnerable to multidimensional poverty 0.5359
## vulnerable to multidimensional poverty|moderately multidimensional poor 0.5191
## moderately multidimensional poor|extremely multidimensional poor 0.4938
## t value
## not multidimensional poor|vulnerable to multidimensional poverty -7.1734
## vulnerable to multidimensional poverty|moderately multidimensional poor -5.7282
## moderately multidimensional poor|extremely multidimensional poor -2.1551
##
## Residual Deviance: 634.6823
## AIC: 664.6823## AIC Model Full : 668.0636## AIC Model Step : 664.6823##
## Model dengan AIC lebih kecil lebih disukai.2.5 Evaluasi Asumsi
2.5.1 Uji Proportional Odds (Brant Test)
Asumsi yang paling kritis dalam regresi ordinal adalah proportional odds — bahwa koefisien regresi berlaku sama untuk setiap batas antar kategori. Asumsi ini diuji menggunakan Brant Test.
Hipotesis:
- \(H_0\): Asumsi proportional odds terpenuhi (koefisien proporsional di semua batas kategori).
- \(H_1\): Asumsi proportional odds tidak terpenuhi.
Jika p-value > 0,05 pada baris Omnibus → asumsi terpenuhi dan model valid digunakan.
## --------------------------------------------------------------------------------------------
## Test for X2 df probability
## --------------------------------------------------------------------------------------------
## Omnibus 59.48 24 0
## genderfemale 0.69 2 0.71
## marital_statusmarried 4.88 2 0.09
## marital_statusdivorced 2.44 2 0.3
## marital_statuswidowed 0.04 2 0.98
## fam_size 9.21 2 0.01
## occupationcausal wage worker 0.78 2 0.68
## occupationowner of private bisness or self employed 4.3 2 0.12
## occupationsalaried 0.14 2 0.93
## creditget credit 8.26 2 0.02
## house_ownershipown house 2.94 2 0.23
## female_incomefemale participate 7.11 2 0.03
## seving_ratio 8.6 2 0.01
## --------------------------------------------------------------------------------------------
##
## H0: Parallel Regression Assumption holds## X2 df probability
## Omnibus 59.47683584 24 0.00007569
## genderfemale 0.68522428 2 0.70991351
## marital_statusmarried 4.87527444 2 0.08736704
## marital_statusdivorced 2.43794511 2 0.29553366
## marital_statuswidowed 0.03765718 2 0.98134756
## fam_size 9.20710564 2 0.01001619
## occupationcausal wage worker 0.77682460 2 0.67813269
## occupationowner of private bisness or self employed 4.30387029 2 0.11625896
## occupationsalaried 0.13792892 2 0.93335985
## creditget credit 8.26396775 2 0.01605100
## house_ownershipown house 2.93629206 2 0.23035216
## female_incomefemale participate 7.11431707 2 0.02851975
## seving_ratio 8.59648938 2 0.013592402.5.2 Uji Signifikansi Parameter (Uji Wald)
ctable <- coef(summary(model_step))
p_val <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
ctable_full <- cbind(ctable, "p-value" = round(p_val, 4))
print(ctable_full)## Value
## genderfemale 1.1003835
## marital_statusmarried -1.8222787
## marital_statusdivorced -0.5216513
## marital_statuswidowed -0.9416084
## fam_size 0.2118044
## occupationcausal wage worker 0.4046438
## occupationowner of private bisness or self employed -0.5846779
## occupationsalaried -1.5568853
## creditget credit -0.9921105
## house_ownershipown house -0.7593586
## female_incomefemale participate -0.8937387
## seving_ratio -2.7076302
## not multidimensional poor|vulnerable to multidimensional poverty -3.8445266
## vulnerable to multidimensional poverty|moderately multidimensional poor -2.9737448
## moderately multidimensional poor|extremely multidimensional poor -1.0640841
## Std. Error
## genderfemale 0.30978571
## marital_statusmarried 0.29537255
## marital_statusdivorced 0.44140163
## marital_statuswidowed 0.71150757
## fam_size 0.06530281
## occupationcausal wage worker 0.50786386
## occupationowner of private bisness or self employed 0.51231784
## occupationsalaried 0.44823728
## creditget credit 0.31712249
## house_ownershipown house 0.24837618
## female_incomefemale participate 0.27802112
## seving_ratio 0.91356124
## not multidimensional poor|vulnerable to multidimensional poverty 0.53594019
## vulnerable to multidimensional poverty|moderately multidimensional poor 0.51914008
## moderately multidimensional poor|extremely multidimensional poor 0.49376066
## t value
## genderfemale 3.5520796
## marital_statusmarried -6.1694246
## marital_statusdivorced -1.1818065
## marital_statuswidowed -1.3233990
## fam_size 3.2434192
## occupationcausal wage worker 0.7967565
## occupationowner of private bisness or self employed -1.1412405
## occupationsalaried -3.4733508
## creditget credit -3.1284774
## house_ownershipown house -3.0572924
## female_incomefemale participate -3.2146433
## seving_ratio -2.9638190
## not multidimensional poor|vulnerable to multidimensional poverty -7.1734247
## vulnerable to multidimensional poverty|moderately multidimensional poor -5.7282128
## moderately multidimensional poor|extremely multidimensional poor -2.1550605
## p-value
## genderfemale 0.0004
## marital_statusmarried 0.0000
## marital_statusdivorced 0.2373
## marital_statuswidowed 0.1857
## fam_size 0.0012
## occupationcausal wage worker 0.4256
## occupationowner of private bisness or self employed 0.2538
## occupationsalaried 0.0005
## creditget credit 0.0018
## house_ownershipown house 0.0022
## female_incomefemale participate 0.0013
## seving_ratio 0.0030
## not multidimensional poor|vulnerable to multidimensional poverty 0.0000
## vulnerable to multidimensional poverty|moderately multidimensional poor 0.0000
## moderately multidimensional poor|extremely multidimensional poor 0.03122.5.3 Uji Goodness of Fit (Lipsitz Test)
##
## Lipsitz goodness of fit test for ordinal response models
##
## data: formula: MDP ~ gender + marital_status + fam_size + occupation + credit + formula: house_ownership + female_income + seving_ratio
## LR statistic = 25.437, df = 9, p-value = 0.002524Interpretasi Lipsitz Test: \(H_0\): Model sudah fit dengan data (goodness of fit adequate). Jika p-value > 0,05, maka tidak ada alasan untuk menolak \(H_0\) — model dinyatakan fit.
2.6 Hasil dan Interpretasi Model
2.6.1 Tabel Koefisien dan Odds Ratio
coef_step <- coef(summary(model_step))
p_val2 <- pnorm(abs(coef_step[, "t value"]), lower.tail = FALSE) * 2
idx_pred <- !grepl("\\|", rownames(coef_step))
coef_pred <- coef_step[idx_pred, , drop = FALSE]
p_pred_or <- p_val2[idx_pred]
ci <- confint(model_step)
ci_pred <- ci[rownames(ci) %in% rownames(coef_pred), , drop = FALSE]
tbl_or <- data.frame(
Variabel = rownames(coef_pred),
Koefisien = round(coef_pred[, "Value"], 4),
Std_Error = round(coef_pred[, "Std. Error"], 4),
t_value = round(coef_pred[, "t value"], 4),
p_value = round(p_pred_or, 4),
OR = round(exp(coef_pred[, "Value"]), 4),
CI_Lower = round(exp(ci_pred[, 1]), 4),
CI_Upper = round(exp(ci_pred[, 2]), 4),
Signifikan = ifelse(p_pred_or < 0.001, "***",
ifelse(p_pred_or < 0.01, "**",
ifelse(p_pred_or < 0.05, "*",
ifelse(p_pred_or < 0.1, ".", ""))))
)
rownames(tbl_or) <- NULL
kable(tbl_or,
col.names = c("Variabel", "β", "SE", "t-value", "p-value",
"OR", "CI 2.5%", "CI 97.5%", "Sig."),
caption = "Koefisien, Odds Ratio, dan Uji Signifikansi Model Ordinal") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE) %>%
column_spec(1, bold = TRUE, color = "#2f7f73") %>%
column_spec(5, bold = TRUE,
color = ifelse(tbl_or$p_value < 0.05, "#b84f5a", "#64748b")) %>%
footnote(general = "Signifikansi: *** p<0.001 | ** p<0.01 | * p<0.05 | . p<0.1")| Variabel | β | SE | t-value | p-value | OR | CI 2.5% | CI 97.5% | Sig. |
|---|---|---|---|---|---|---|---|---|
| genderfemale | 1.1004 | 0.3098 | 3.5521 | 0.0004 | 3.0053 | 1.6449 | 5.5531 | *** |
| marital_statusmarried | -1.8223 | 0.2954 | -6.1694 | 0.0000 | 0.1617 | 0.0897 | 0.2861 | *** |
| marital_statusdivorced | -0.5217 | 0.4414 | -1.1818 | 0.2373 | 0.5935 | 0.2508 | 1.4241 | |
| marital_statuswidowed | -0.9416 | 0.7115 | -1.3234 | 0.1857 | 0.3900 | 0.0952 | 1.5962 | |
| fam_size | 0.2118 | 0.0653 | 3.2434 | 0.0012 | 1.2359 | 1.0889 | 1.4073 | ** |
| occupationcausal wage worker | 0.4046 | 0.5079 | 0.7968 | 0.4256 | 1.4988 | 0.5444 | 4.0287 | |
| occupationowner of private bisness or self employed | -0.5847 | 0.5123 | -1.1412 | 0.2538 | 0.5573 | 0.1996 | 1.5013 | |
| occupationsalaried | -1.5569 | 0.4482 | -3.4734 | 0.0005 | 0.2108 | 0.0848 | 0.4969 | *** |
| creditget credit | -0.9921 | 0.3171 | -3.1285 | 0.0018 | 0.3708 | 0.1973 | 0.6864 | ** |
| house_ownershipown house | -0.7594 | 0.2484 | -3.0573 | 0.0022 | 0.4680 | 0.2867 | 0.7600 | ** |
| female_incomefemale participate | -0.8937 | 0.2780 | -3.2146 | 0.0013 | 0.4091 | 0.2359 | 0.7028 | ** |
| seving_ratio | -2.7076 | 0.9136 | -2.9638 | 0.0030 | 0.0667 | 0.0100 | 0.3325 | ** |
| Note: | ||||||||
| Signifikansi: *** p<0.001 | ** p<0.01 | * p<0.05 | . p<0.1 |
2.6.2 Tabel Threshold (Intercept Antar Kategori)
idx_thr <- grepl("\\|", rownames(coef_step))
coef_thr <- coef_step[idx_thr, , drop = FALSE]
tbl_thr <- data.frame(
# INI BAGIAN YANG DIGANTI: Sisipkan gsub() untuk mengubah | menjadi vs
Threshold = gsub("\\|", " vs ", rownames(coef_thr)),
Nilai = round(coef_thr[, "Value"], 4),
Std_Error = round(coef_thr[, "Std. Error"], 4)
)
rownames(tbl_thr) <- NULL
kable(tbl_thr,
col.names = c("Threshold", "Nilai", "Std. Error"),
caption = "Threshold (Intercept) Model Ordinal") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE) %>%
column_spec(1, bold = TRUE, color = "#2f7f73")| Threshold | Nilai | Std. Error |
|---|---|---|
| not multidimensional poor vs vulnerable to multidimensional poverty | -3.8445 | 0.5359 |
| vulnerable to multidimensional poverty vs moderately multidimensional poor | -2.9737 | 0.5191 |
| moderately multidimensional poor vs extremely multidimensional poor | -1.0641 | 0.4938 |
2.6.3 Visualisasi Odds Ratio
tbl_or_plot <- tbl_or |>
mutate(
Warna = ifelse(OR > 1, "Meningkatkan Risiko", "Menurunkan Risiko"),
Variabel = factor(Variabel, levels = Variabel[order(OR)])
)
ggplot(tbl_or_plot, aes(x = OR, y = Variabel, color = Warna)) +
geom_vline(xintercept = 1, linetype = "dashed",
color = "grey50", linewidth = 0.8) +
geom_errorbarh(aes(xmin = CI_Lower, xmax = CI_Upper),
height = 0.3, linewidth = 0.9) +
geom_point(size = 3.5) +
scale_color_manual(values = c("Meningkatkan Risiko" = "#b84f5a",
"Menurunkan Risiko" = "#2f7f73")) +
scale_x_log10(labels = scales::number_format(accuracy = 0.01)) +
labs(
title = "Odds Ratio Prediktor terhadap Status Kemiskinan Multidimensi",
subtitle = "Skala logaritmik — garis putus-putus = OR 1 (tidak ada efek)",
x = "Odds Ratio (skala log)",
y = NULL,
color = NULL
) +
theme_adk() +
theme(legend.position = "bottom")
2.6.4 Interpretasi Koefisien
Interpretasi dilakukan terhadap variabel yang signifikan secara statistik (p < 0,05). Dalam model proportional odds, koefisien positif berarti OR > 1, artinya variabel tersebut meningkatkan kecenderungan untuk berada di kategori kemiskinan yang lebih parah. Sebaliknya, koefisien negatif berarti OR < 1, artinya variabel tersebut bersifat protektif.
educ_level — Tingkat Pendidikan
Dibandingkan kepala rumah tangga yang buta huruf (referensi), kepala
rumah tangga berpendidikan lebih tinggi memiliki odds yang lebih rendah
untuk masuk kategori kemiskinan yang lebih parah. Misalnya, jika
educ_levelhigher edcuation memiliki OR = 0,12, maka kepala
rumah tangga berpendidikan tinggi memiliki odds untuk berada pada
kategori kemiskinan yang lebih parah hanya 12% dari
odds kepala rumah tangga buta huruf, dengan semua variabel lain
konstan.
occupation — Jenis Pekerjaan
Dibandingkan kelompok pengangguran, pekerja formal (salaried) atau wirausaha memiliki OR yang lebih kecil dari 1, artinya mereka cenderung berada di kategori kemiskinan yang lebih rendah. Ini konsisten dengan ekspektasi bahwa pekerjaan tetap memberikan pendapatan yang lebih stabil dan menurunkan risiko kemiskinan multidimensi.
house_ownership — Kepemilikan Rumah
Rumah tangga yang memiliki rumah sendiri cenderung memiliki odds lebih rendah untuk berada di kategori kemiskinan yang lebih parah. Kepemilikan aset fisik seperti rumah merefleksikan akumulasi kekayaan jangka panjang yang menjadi penyangga terhadap kemiskinan.
credit — Akses Kredit
Akses terhadap kredit formal berkaitan dengan odds yang lebih rendah untuk masuk kategori kemiskinan yang lebih tinggi. Ini menunjukkan bahwa inklusi keuangan — dalam bentuk akses pinjaman — berperan penting dalam strategi rumah tangga keluar dari kemiskinan.
depratio — Rasio Ketergantungan
Rasio ketergantungan yang lebih tinggi (lebih banyak anggota tidak produktif per anggota produktif) dikaitkan dengan odds yang lebih tinggi untuk berada di kategori kemiskinan lebih parah. Semakin besar beban tanggungan, semakin terbatas sumber daya per kapita yang tersedia.
seving_ratio — Rasio Tabungan
Rasio tabungan yang lebih tinggi dikaitkan dengan odds yang lebih rendah untuk masuk kategori kemiskinan yang lebih parah. Kemampuan menabung mencerminkan ketahanan finansial dan kapasitas menghadapi guncangan ekonomi.
2.7 Evaluasi Performa Model
2.7.1 Prediksi pada Data Testing
2.7.2 Confusion Matrix
## Confusion Matrix and Statistics
##
## Reference
## Prediction not multidimensional poor
## not multidimensional poor 10
## vulnerable to multidimensional poverty 0
## moderately multidimensional poor 5
## extremely multidimensional poor 1
## Reference
## Prediction vulnerable to multidimensional poverty
## not multidimensional poor 8
## vulnerable to multidimensional poverty 0
## moderately multidimensional poor 6
## extremely multidimensional poor 1
## Reference
## Prediction moderately multidimensional poor
## not multidimensional poor 2
## vulnerable to multidimensional poverty 0
## moderately multidimensional poor 13
## extremely multidimensional poor 6
## Reference
## Prediction extremely multidimensional poor
## not multidimensional poor 5
## vulnerable to multidimensional poverty 0
## moderately multidimensional poor 5
## extremely multidimensional poor 17
##
## Overall Statistics
##
## Accuracy : 0.5063
## 95% CI : (0.3914, 0.6208)
## No Information Rate : 0.3418
## P-Value [Acc > NIR] : 0.001887
##
## Kappa : 0.3239
##
## Mcnemar's Test P-Value : 0.004091
##
## Statistics by Class:
##
## Class: not multidimensional poor
## Sensitivity 0.6250
## Specificity 0.7619
## Pos Pred Value 0.4000
## Neg Pred Value 0.8889
## Prevalence 0.2025
## Detection Rate 0.1266
## Detection Prevalence 0.3165
## Balanced Accuracy 0.6935
## Class: vulnerable to multidimensional poverty
## Sensitivity 0.0000
## Specificity 1.0000
## Pos Pred Value NaN
## Neg Pred Value 0.8101
## Prevalence 0.1899
## Detection Rate 0.0000
## Detection Prevalence 0.0000
## Balanced Accuracy 0.5000
## Class: moderately multidimensional poor
## Sensitivity 0.6190
## Specificity 0.7241
## Pos Pred Value 0.4483
## Neg Pred Value 0.8400
## Prevalence 0.2658
## Detection Rate 0.1646
## Detection Prevalence 0.3671
## Balanced Accuracy 0.6716
## Class: extremely multidimensional poor
## Sensitivity 0.6296
## Specificity 0.8462
## Pos Pred Value 0.6800
## Neg Pred Value 0.8148
## Prevalence 0.3418
## Detection Rate 0.2152
## Detection Prevalence 0.3165
## Balanced Accuracy 0.7379cm_df_ord <- as.data.frame(cm_ord$table)
colnames(cm_df_ord) <- c("Prediksi", "Aktual", "Frekuensi")
label_short <- c(
"not multidimensional poor" = "Tidak\nMiskin",
"vulnerable to multidimensional poverty" = "Rentan",
"moderately multidimensional poor" = "Miskin\nModerat",
"extremely multidimensional poor" = "Sangat\nMiskin"
)
cm_df_ord <- cm_df_ord |>
mutate(
Prediksi = recode(as.character(Prediksi), !!!label_short),
Aktual = recode(as.character(Aktual), !!!label_short),
Prediksi = factor(Prediksi,
levels = c("Tidak\nMiskin","Rentan",
"Miskin\nModerat","Sangat\nMiskin")),
Aktual = factor(Aktual,
levels = c("Tidak\nMiskin","Rentan",
"Miskin\nModerat","Sangat\nMiskin"))
)
ggplot(cm_df_ord, aes(x = Aktual, y = Prediksi, fill = Frekuensi)) +
geom_tile(color = "white", linewidth = 0.8) +
geom_text(aes(label = Frekuensi), size = 5, fontface = "bold") +
scale_fill_gradient(low = "#f0fbff", high = "#0077b6") +
labs(
title = "Confusion Matrix — Model Ordinal",
subtitle = "Perbandingan prediksi model dengan kategori aktual (data testing)",
fill = "Frekuensi"
) +
theme_adk()
2.7.3 Akurasi dan Metrik Evaluasi
overall_acc <- cm_ord$overall["Accuracy"]
kappa_val <- cm_ord$overall["Kappa"]
cat("Akurasi Keseluruhan :", round(overall_acc * 100, 2), "%\n")## Akurasi Keseluruhan : 50.63 %## Cohen's Kappa : 0.3239##
## Interpretasi Kappa:## < 0.20 : lemah | 0.21-0.40 : fair | 0.41-0.60 : sedang## 0.61-0.80 : baik | 0.81-1.00 : sangat baikclass_metrics <- data.frame(
Kategori = rownames(cm_ord$byClass),
Sensitivity = round(cm_ord$byClass[, "Sensitivity"], 4),
Specificity = round(cm_ord$byClass[, "Specificity"], 4),
Precision = round(cm_ord$byClass[, "Precision"], 4),
F1_Score = round(cm_ord$byClass[, "F1"], 4)
)
rownames(class_metrics) <- NULL
kable(class_metrics,
caption = "Metrik Evaluasi per Kategori MDP") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE) %>%
column_spec(1, bold = TRUE, color = "#2f7f73")| Kategori | Sensitivity | Specificity | Precision | F1_Score |
|---|---|---|---|---|
| Class: not multidimensional poor | 0.6250 | 0.7619 | 0.4000 | 0.4878 |
| Class: vulnerable to multidimensional poverty | 0.0000 | 1.0000 | NA | NA |
| Class: moderately multidimensional poor | 0.6190 | 0.7241 | 0.4483 | 0.5200 |
| Class: extremely multidimensional poor | 0.6296 | 0.8462 | 0.6800 | 0.6538 |
2.8 Kesimpulan
Analisis regresi logistik ordinal terhadap data kemiskinan multidimensi (MPI) menghasilkan beberapa temuan utama yang layak dicermati dari perspektif sosial-ekonomi:
Pendidikan adalah prediktor terkuat. Kepala rumah tangga dengan tingkat pendidikan yang lebih tinggi secara konsisten memiliki odds lebih rendah untuk masuk ke kategori kemiskinan yang lebih parah. Ini menegaskan bahwa investasi di bidang pendidikan bukan sekadar tujuan sosial, melainkan juga instrumen pengentasan kemiskinan yang paling efektif dalam jangka panjang.
Pekerjaan formal vs. informal menciptakan kesenjangan yang nyata. Rumah tangga dengan kepala keluarga pekerja formal atau wirausaha memiliki posisi yang jauh lebih baik dibandingkan pekerja harian dan pengangguran. Ketidakstabilan pendapatan dari pekerjaan informal menjadi salah satu faktor yang memperbesar kerentanan multidimensi.
Rasio ketergantungan mencerminkan tekanan ekonomi yang sesungguhnya. Semakin tinggi rasio anggota tidak produktif dalam rumah tangga, semakin besar risiko kemiskinan. Implikasinya, program perlindungan sosial perlu lebih memperhatikan rumah tangga dengan banyak tanggungan anak atau lansia.
Inklusi keuangan dan kepemilikan aset berperan penting. Akses kredit dan kepemilikan rumah sendiri terbukti berhubungan dengan status kemiskinan yang lebih baik. Ini mendukung argumen bahwa kebijakan yang memperluas akses layanan keuangan formal ke kelompok rentan bisa menjadi intervensi yang efektif.
Kemampuan menabung sebagai indikator ketahanan. Rasio tabungan yang lebih tinggi dikaitkan dengan kondisi yang lebih baik — bukan hanya karena tabungan itu sendiri, tetapi karena tabungan mencerminkan surplus pendapatan dan kapasitas adaptasi terhadap guncangan.
Model yang dihasilkan memenuhi asumsi proportional odds dan menunjukkan kinerja yang memadai berdasarkan akurasi, nilai Kappa, dan pseudo R-squared. Dengan demikian, model ini dapat dijadikan dasar untuk merancang program penargetan kemiskinan yang lebih tepat sasaran.
3 Regresi Logistik Multinomial
Model ini merupakan perluasan langsung dari regresi logistik biner untuk kasus dengan variabel respons nominal yang terdiri dari tiga kategori atau lebih.
Model Regresi Logistik Multinomial Misalkan \(Y_i\) adalah variabel respons dengan \(J\) kategori, dan kategori ke-\(J\) dipilih sebagai referensi. Untuk setiap kategori \(j = 1, 2, \ldots, J-1\), model baseline-category logit didefinisikan sebagai:
\[ \ln\!\left(\frac{\pi_{ij}}{\pi_{iJ}}\right) = \beta_{j0} + \beta_{j1}x_{i1} + \beta_{j2}x_{i2} + \cdots + \beta_{jp}x_{ip} \]
di mana \(\pi_{ij} = P(Y_i = j \mid \mathbf{x}_i)\) adalah peluang observasi ke-\(i\) masuk kategori \(j\).
Peluang tiap kategori diperoleh dari:
\[ \pi_{ij} = \frac{\exp(\beta_{j0} + \boldsymbol{\beta}_j^T \mathbf{x}_i)}{1 + \sum_{k=1}^{J-1} \exp(\beta_{k0} + \boldsymbol{\beta}_k^T \mathbf{x}_i)}, \quad j = 1, \ldots, J-1 \]
\[ \pi_{iJ} = \frac{1}{1 + \sum_{k=1}^{J-1} \exp(\beta_{k0} + \boldsymbol{\beta}_k^T \mathbf{x}_i)} \]
dengan syarat \(\sum_{j=1}^{J} \pi_{ij} = 1\).
3.0.1 Interpretasi Koefisien melalui Relative Risk Ratio (RRR)
Karena model menggunakan fungsi tautan log, koefisien \(\beta_{jk}\) tidak dibaca secara langsung. Interpretasi yang tepat dilakukan melalui Relative Risk Ratio (RRR):
\[ \text{RRR}_{jk} = \exp(\beta_{jk}) \]
- Jika \(\text{RRR}_{jk} > 1\): kenaikan satu unit \(x_k\) meningkatkan risiko relatif berada di kategori \(j\) dibanding kategori referensi sebesar \((\text{RRR}_{jk} - 1) \times 100\%\).
- Jika \(\text{RRR}_{jk} < 1\): kenaikan satu unit \(x_k\) menurunkan risiko relatif berada di kategori \(j\) dibanding kategori referensi sebesar \((1 - \text{RRR}_{jk}) \times 100\%\).
- Jika \(\text{RRR}_{jk} = 1\): prediktor tidak memiliki efek terhadap kategori tersebut relatif terhadap referensi.
3.1 Sumber Data dan Deskripsi Variabel
Dataset yang digunakan adalah Estimation of Obesity Levels Based on Eating Habits and Physical Condition yang tersedia secara publik di UCI Machine Learning Repository. Data ini dikumpulkan dari responden di Meksiko, Peru, dan Kolombia melalui platform web dengan pertanyaan mengenai kebiasaan makan, kondisi fisik, dan gaya hidup. Dataset berisi 2.111 observasi, dengan beberapa tingkat status berat badan mulai dari berat badan kurang hingga obesitas tingkat III, model multinomial menjadi pilihan yang tepat karena kategori-kategori tersebut, meski terlihat memiliki urutan, diperlakukan sebagai nominal dalam konteks ini karena faktor-faktor penentu tiap kategori bisa sangat berbeda.
Referensi: Palechor, F. M., & de la Hoz Manotas, A. (2019). Dataset for estimation of obesity levels based on eating habits and physical condition in individuals from Colombia, Peru and Mexico. Data in Brief, 25, 104344. UCI Machine Learning Repository: https://archive.ics.uci.edu/dataset/544
Jumlah observasi: 2.111 | Jumlah variabel: 17
| Variabel | Tipe | Keterangan |
|---|---|---|
| NObeyesdad | Faktor (7 kategori) | Tingkat obesitas — variabel respons (7 kategori) |
| Gender | Faktor | Jenis kelamin responden |
| Age | Numerik | Usia responden (tahun) |
| Height | Numerik | Tinggi badan (meter) |
| Weight | Numerik | Berat badan (kilogram) |
| family_history_with_overweight | Faktor (yes/no) | Riwayat keluarga dengan kelebihan berat badan |
| FAVC | Faktor (yes/no) | Konsumsi makanan berkalori tinggi (Frequently consume high caloric food) |
| FCVC | Numerik | Frekuensi konsumsi sayur (1=jarang, 3=sering) |
| NCP | Numerik | Jumlah makan utama per hari |
| CAEC | Faktor (4 kategori) | Konsumsi makanan di antara waktu makan (no/Sometimes/Frequently/Always) |
| SMOKE | Faktor (yes/no) | Status merokok (yes/no) |
| CH2O | Numerik | Konsumsi air harian (1=sedikit, 3=banyak) |
| SCC | Faktor (yes/no) | Pemantauan kalori harian (yes/no) |
| FAF | Numerik | Frekuensi aktivitas fisik per minggu (0=tidak pernah, 3=sering) |
| TUE | Numerik | Waktu penggunaan perangkat teknologi per hari (jam) |
| CALC | Faktor (4 kategori) | Konsumsi alkohol (no/Sometimes/Frequently/Always) |
| MTRANS | Faktor (5 kategori) | Moda transportasi utama (Automobile/Bike/Motorbike/Public_Transportation/Walking) |
Dari jumlah variabel yang ada (16 prediktor), dilakukan seleksi untuk
menjaga model tetap bermakna dan terhindar dari overfitting. Variabel
Weight dan Height tidak diikutsertakan
karena keduanya secara langsung digunakan untuk menghitung BMI yang
menjadi dasar penetapan NObeyesdad. Prediktor yang dipilih
fokus pada faktor perilaku dan demografis yang secara substansif relevan
untuk pemodelan pola obesitas.
3.2 Persiapan Data
3.2.1 Membaca dan Menyiapkan Dataset
# Dataset tersimpan dalam format xlsx dengan nilai CSV dalam satu kolom
df_ob_raw <- read_excel("C:/Users/Lenovo/Downloads/ObesityDataSet_raw_and_data_sinthetic.csv.xlsx")
# Parse: header ada di nama kolom pertama, baris adalah data CSV
col_header <- names(df_ob_raw)[1]
rows_data <- df_ob_raw[[1]]
csv_text <- paste(c(col_header, rows_data), collapse = "\n")
obesity <- read.csv(text = csv_text, stringsAsFactors = FALSE)
cat("Dimensi data:", nrow(obesity), "observasi,", ncol(obesity), "variabel\n")## Dimensi data: 2111 observasi, 17 variabelhead(obesity, 5) %>%
kable(caption = "5 Baris Pertama Dataset Obesity Levels") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE) %>%
column_spec(1, bold = TRUE, color = "#2f7f73") %>%
scroll_box(width = "100%", height = "280px")| Gender | Age | Height | Weight | family_history_with_overweight | FAVC | FCVC | NCP | CAEC | SMOKE | CH2O | SCC | FAF | TUE | CALC | MTRANS | NObeyesdad |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Female | 21 | 1.62 | 64.0 | yes | no | 2 | 3 | Sometimes | no | 2 | no | 0 | 1 | no | Public_Transportation | Normal_Weight |
| Female | 21 | 1.52 | 56.0 | yes | no | 3 | 3 | Sometimes | yes | 3 | yes | 3 | 0 | Sometimes | Public_Transportation | Normal_Weight |
| Male | 23 | 1.80 | 77.0 | yes | no | 2 | 3 | Sometimes | no | 2 | no | 2 | 1 | Frequently | Public_Transportation | Normal_Weight |
| Male | 27 | 1.80 | 87.0 | no | no | 3 | 3 | Sometimes | no | 2 | no | 2 | 0 | Frequently | Walking | Overweight_Level_I |
| Male | 22 | 1.78 | 89.8 | no | no | 2 | 1 | Sometimes | no | 2 | no | 0 | 0 | Sometimes | Public_Transportation | Overweight_Level_II |
3.2.2 Seleksi Variabel dan Recoding
# Urutan kategori respons (dari berat badan kurang ke obesitas ekstrem)
ob_levels <- c(
"Insufficient_Weight",
"Normal_Weight",
"Overweight_Level_I",
"Overweight_Level_II",
"Obesity_Type_I",
"Obesity_Type_II",
"Obesity_Type_III"
)
ob_labels <- c(
"Insufficient_Weight" = "Berat Kurang",
"Normal_Weight" = "Normal",
"Overweight_Level_I" = "Overweight I",
"Overweight_Level_II" = "Overweight II",
"Obesity_Type_I" = "Obesitas I",
"Obesity_Type_II" = "Obesitas II",
"Obesity_Type_III" = "Obesitas III"
)
obesity_clean <- obesity %>%
# Seleksi variabel yang relevan (tanpa Weight dan Height)
dplyr::select(NObeyesdad, Gender, Age, family_history_with_overweight,
FAVC, FCVC, NCP, CAEC, SMOKE, CH2O, SCC, FAF, TUE, CALC, MTRANS) %>%
mutate(
# Respons: factor dengan urutan kategori
NObeyesdad = factor(NObeyesdad, levels = ob_levels),
# Variabel kategorik
Gender = factor(Gender, levels = c("Male", "Female")),
family_history_with_overweight = factor(family_history_with_overweight,
levels = c("no", "yes")),
FAVC = factor(FAVC, levels = c("no", "yes")),
CAEC = factor(CAEC, levels = c("no", "Sometimes", "Frequently", "Always")),
SMOKE = factor(SMOKE, levels = c("no", "yes")),
SCC = factor(SCC, levels = c("no", "yes")),
CALC = factor(CALC, levels = c("no", "Sometimes", "Frequently", "Always")),
MTRANS = factor(MTRANS, levels = c("Automobile", "Bike", "Motorbike",
"Public_Transportation", "Walking"))
) %>%
na.omit()
cat("Data siap dianalisis:", nrow(obesity_clean), "observasi,",
ncol(obesity_clean), "variabel\n")## Data siap dianalisis: 2111 observasi, 15 variabel3.2.3 Cek Missing Value
miss_ob <- data.frame(
Variabel = names(obesity_clean),
Jumlah_NA = sapply(obesity_clean, function(x) sum(is.na(x)))
) %>% filter(Jumlah_NA > 0)
if (nrow(miss_ob) == 0) {
cat("✅ Tidak ada missing value pada seluruh variabel yang dipilih.\n")
} else {
kable(miss_ob, caption = "Ringkasan Missing Value") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
}## ✅ Tidak ada missing value pada seluruh variabel yang dipilih.3.3 Analisis Deskriptif
3.3.1 Distribusi Variabel Respons
df_resp <- obesity_clean %>%
count(NObeyesdad) %>%
mutate(
Label = ob_labels[as.character(NObeyesdad)],
Label = factor(Label, levels = ob_labels),
persen = n / sum(n) * 100
)
p_bar_ob <- ggplot(df_resp, aes(x = Label, y = n, fill = Label)) +
geom_col(width = 0.65, show.legend = FALSE, color = "white") +
geom_text(aes(label = paste0(n, "\n(", round(persen, 1), "%)")),
vjust = -0.3, size = 3.5, fontface = "bold") +
scale_fill_manual(values = c(
"Berat Kurang" = "#789f90",
"Normal" = "#2f7f73",
"Overweight I" = "#d18b2f",
"Overweight II"= "#e8a030",
"Obesitas I" = "#b84f5a",
"Obesitas II" = "#9c3a44",
"Obesitas III" = "#6b2530"
)) +
scale_y_continuous(expand = expansion(mult = c(0, 0.18))) +
labs(
title = "Distribusi Tingkat Obesitas",
subtitle = "Dataset Obesity Levels — 2.111 Observasi",
x = "Kategori Berat Badan",
y = "Jumlah Observasi"
) +
theme_adk() +
theme(axis.text.x = element_text(face = "bold", size = 9))
p_pie_ob <- df_resp %>%
ggplot(aes(x = "", y = persen, fill = Label)) +
geom_col(width = 1, color = "white") +
coord_polar(theta = "y") +
scale_fill_manual(values = c(
"Berat Kurang" = "#789f90", "Normal" = "#2f7f73",
"Overweight I" = "#d18b2f", "Overweight II" = "#e8a030",
"Obesitas I" = "#b84f5a", "Obesitas II" = "#9c3a44",
"Obesitas III" = "#6b2530"
)) +
labs(title = "Proporsi Kategori", fill = "Kategori") +
theme_void() +
theme(plot.title = element_text(face = "bold", color = "#243142",
hjust = 0.5, size = 14),
legend.position = "right")
p_bar_ob + p_pie_ob + patchwork::plot_layout(widths = c(1.6, 1))
Dataset ini memiliki distribusi yang cukup merata di semua kategori, yang merupakan hasil dari proses augmentasi data menggunakan SMOTE. Kategori Obesitas Tipe I adalah yang paling banyak (351 observasi, 16,6%), sementara Berat Badan Kurang memiliki jumlah paling sedikit (272 observasi, 12,9%). Ketidakseimbangan yang ada masih dalam batas yang dapat ditangani oleh model multinomial.
3.3.2 Profil Kategori Berdasarkan Variabel Kategorik
ob_color7 <- c(
"Berat Kurang" = "#789f90", "Normal" = "#2f7f73",
"Overweight I" = "#5568b8", "Overweight II" = "#d18b2f",
"Obesitas I" = "#b84f5a", "Obesitas II" = "#8b6fc9",
"Obesitas III" = "#243142"
)
plot_ob_stack <- function(var, label_var, judul) {
obesity_clean %>%
mutate(
OB = ob_labels[as.character(NObeyesdad)],
OB = factor(OB, levels = ob_labels)
) %>%
dplyr::rename(Var = all_of(var)) %>%
count(Var, OB) %>%
group_by(Var) %>%
mutate(pct = n / sum(n) * 100) %>%
ungroup() %>%
ggplot(aes(x = fct_rev(as.factor(Var)), y = pct, fill = OB)) +
geom_col(position = "stack", width = 0.6) +
coord_flip() +
scale_fill_manual(values = ob_color7, name = "Kategori") +
scale_y_continuous(labels = percent_format(scale = 1)) +
labs(title = judul, x = label_var, y = "Proporsi (%)") +
theme_adk() +
theme(legend.position = "bottom",
plot.title = element_text(size = 11),
legend.text = element_text(size = 8))
}
pg1 <- plot_ob_stack("Gender", "Jenis Kelamin", "Gender vs Obesitas")
pg2 <- plot_ob_stack("family_history_with_overweight", "Riwayat Keluarga","Riwayat Keluarga vs Obesitas")
pg3 <- plot_ob_stack("FAVC", "Makanan Berkalori","Konsumsi Kalori Tinggi vs Obesitas")
pg4 <- plot_ob_stack("CAEC", "Cemilan", "Konsumsi Cemilan vs Obesitas")
pg5 <- plot_ob_stack("SMOKE", "Merokok", "Merokok vs Obesitas")
pg6 <- plot_ob_stack("CALC", "Alkohol", "Konsumsi Alkohol vs Obesitas")
gridExtra::grid.arrange(pg1, pg2, pg3, pg4, pg5, pg6, ncol = 2)
3.3.3 Distribusi Variabel Kontinu per Kategori
obesity_clean %>%
mutate(
OB = ob_labels[as.character(NObeyesdad)],
OB = factor(OB, levels = ob_labels)
) %>%
dplyr::select(OB, Age, FCVC, NCP, CH2O, FAF, TUE) %>%
pivot_longer(-OB, names_to = "Variabel", values_to = "Nilai") %>%
mutate(Variabel = recode(Variabel,
Age = "Usia",
FCVC = "Konsumsi Sayur",
NCP = "Makan Utama/Hari",
CH2O = "Konsumsi Air",
FAF = "Aktivitas Fisik",
TUE = "Waktu Teknologi"
)) %>%
ggplot(aes(x = OB, y = Nilai, fill = OB)) +
geom_boxplot(alpha = 0.7, outlier.size = 0.9, width = 0.6) +
facet_wrap(~Variabel, scales = "free_y", ncol = 3) +
scale_fill_manual(values = ob_color7) +
labs(
title = "Distribusi Variabel Kontinu Berdasarkan Kategori Obesitas",
x = NULL,
y = "Nilai"
) +
theme_adk() +
theme(
legend.position = "none",
axis.text.x = element_text(angle = 40, hjust = 1, size = 7),
strip.text = element_text(face = "bold")
)
3.3.4 Statistik Deskriptif
desc_ob <- obesity_clean %>%
dplyr::select(Age, FCVC, NCP, CH2O, FAF, TUE) %>%
pivot_longer(everything(), names_to = "Variabel", values_to = "Nilai") %>%
group_by(Variabel) %>%
summarise(
N = n(),
Min = round(min(Nilai), 2),
Q1 = round(quantile(Nilai, 0.25), 2),
Median = round(median(Nilai), 2),
Mean = round(mean(Nilai), 2),
Q3 = round(quantile(Nilai, 0.75), 2),
Max = round(max(Nilai), 2),
SD = round(sd(Nilai), 2),
.groups = "drop"
)
kable(desc_ob,
caption = "Statistik Deskriptif Variabel Kontinu — Dataset Obesity") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE) %>%
column_spec(1, bold = TRUE, color = "#2f7f73")| Variabel | N | Min | Q1 | Median | Mean | Q3 | Max | SD |
|---|---|---|---|---|---|---|---|---|
| Age | 2111 | 14 | 19.95 | 22.78 | 24.31 | 26.00 | 61 | 6.35 |
| CH2O | 2111 | 1 | 1.58 | 2.00 | 2.01 | 2.48 | 3 | 0.61 |
| FAF | 2111 | 0 | 0.12 | 1.00 | 1.01 | 1.67 | 3 | 0.85 |
| FCVC | 2111 | 1 | 2.00 | 2.39 | 2.42 | 3.00 | 3 | 0.53 |
| NCP | 2111 | 1 | 2.66 | 3.00 | 2.69 | 3.00 | 4 | 0.78 |
| TUE | 2111 | 0 | 0.00 | 0.63 | 0.66 | 1.00 | 2 | 0.61 |
3.4 Pemodelan Regresi Logistik Multinomial
3.4.1 Pembagian Data (Train-Test Split)
set.seed(123)
n_ob <- nrow(obesity_clean)
idx_tr_ob <- sample(seq_len(n_ob), size = floor(0.8 * n_ob))
ob_train <- obesity_clean[idx_tr_ob, ]
ob_test <- obesity_clean[-idx_tr_ob, ]
cat("Data training :", nrow(ob_train), "observasi\n")## Data training : 1688 observasi## Data testing : 423 observasi3.4.2 Model Multinomial Penuh
Model dibangun menggunakan fungsi multinom() dari paket
nnet. Kategori referensi yang dipilih adalah
Normal_Weight — kategori ini paling
bermakna secara klinis karena merepresentasikan kondisi ideal yang ingin
dicapai.
# Tetapkan referensi ke Normal_Weight
ob_train$NObeyesdad <- relevel(ob_train$NObeyesdad, ref = "Normal_Weight")
model_multi_full <- nnet::multinom(
NObeyesdad ~ Gender + Age + family_history_with_overweight +
FAVC + FCVC + NCP + CAEC + SMOKE + CH2O +
SCC + FAF + TUE + CALC + MTRANS,
data = ob_train,
MaxNWts = 10000,
maxit = 500,
trace = FALSE
)## Call:
## nnet::multinom(formula = NObeyesdad ~ Gender + Age + family_history_with_overweight +
## FAVC + FCVC + NCP + CAEC + SMOKE + CH2O + SCC + FAF + TUE +
## CALC + MTRANS, data = ob_train, MaxNWts = 10000, maxit = 500,
## trace = FALSE)
##
## Coefficients:
## (Intercept) GenderFemale Age
## Insufficient_Weight -0.524180 0.59505255 -0.18685628
## Overweight_Level_I -3.810505 0.37783851 0.07924518
## Overweight_Level_II -574.660040 -0.69573803 0.22805097
## Obesity_Type_I -8.943029 -0.03038732 0.17741868
## Obesity_Type_II -16.919057 -5.76885935 0.35134557
## Obesity_Type_III -5391.546742 412.34627000 21.03680406
## family_history_with_overweightyes FAVCyes FCVC
## Insufficient_Weight -0.6533457 0.4319158 0.85138025
## Overweight_Level_I 1.0342792 1.3754341 -0.01129607
## Overweight_Level_II 2.9120984 -0.9676715 -0.14308302
## Obesity_Type_I 3.3580962 1.5193411 -0.45947197
## Obesity_Type_II 6.3827913 0.8162618 0.96771117
## Obesity_Type_III 237.7685193 391.1214975 788.11960934
## NCP CAECSometimes CAECFrequently CAECAlways
## Insufficient_Weight 0.3502182 0.5214010 1.386544 -3.0065645
## Overweight_Level_I -0.3156941 -0.1544003 -2.239376 -2.4904291
## Overweight_Level_II -0.5518074 569.2193753 566.524396 566.0998777
## Obesity_Type_I -0.5137568 2.5505034 -1.006562 -0.1276793
## Obesity_Type_II 0.0653941 0.3288692 -4.458426 -2.5180320
## Obesity_Type_III 159.0521566 44.5591915 -469.873655 -168.9716106
## SMOKEyes CH2O SCCyes FAF
## Insufficient_Weight -2.99197477 0.4562384 -0.4485967 -0.05660268
## Overweight_Level_I -0.99020118 0.4116973 1.3704045 -0.23063604
## Overweight_Level_II -0.35938516 0.5790872 -1.9943366 -0.43900728
## Obesity_Type_I 0.02021295 0.8433039 -2.1204127 -0.33552820
## Obesity_Type_II 1.40121697 -0.1711206 -0.3178026 -0.71976388
## Obesity_Type_III -19.51109780 25.5124126 -99.5332683 51.23617258
## TUE CALCSometimes CALCFrequently CALCAlways
## Insufficient_Weight 0.64465320 -6.483407e-02 -3.02868235 -366.40862
## Overweight_Level_I -0.19573628 8.699477e-01 1.07004821 -299.26555
## Overweight_Level_II 0.13736501 -8.988024e-01 -0.03804888 -305.09114
## Obesity_Type_I -0.05061857 -9.646143e-01 -0.26857235 -162.01315
## Obesity_Type_II -0.34602550 6.168638e-05 -1.12600609 -10.97635
## Obesity_Type_III -78.03058871 2.699148e+02 -328.79380935 57.14806
## MTRANSBike MTRANSMotorbike MTRANSPublic_Transportation
## Insufficient_Weight -654.1482775 -748.6053676 -0.8072486
## Overweight_Level_I 0.6708967 -0.8738375 0.4400443
## Overweight_Level_II -457.0824059 -0.2333470 1.2034342
## Obesity_Type_I -506.6438400 1.7413555 0.8078209
## Obesity_Type_II -1.0320740 -1215.7316404 2.1349645
## Obesity_Type_III 498.6461496 212.7428783 808.9362006
## MTRANSWalking
## Insufficient_Weight -2.4897259
## Overweight_Level_I -0.6620407
## Overweight_Level_II -0.9722811
## Obesity_Type_I -2.5269006
## Obesity_Type_II -2631.9059554
## Obesity_Type_III 479.1733395
##
## Std. Errors:
## (Intercept) GenderFemale Age
## Insufficient_Weight 1.3921438 0.2524431 0.03774675
## Overweight_Level_I 1.2132670 0.2419092 0.02708917
## Overweight_Level_II 1.0087823 0.2719497 0.02818337
## Obesity_Type_I 1.6729304 0.2596369 0.02787324
## Obesity_Type_II 2.6040697 0.8019793 0.03522948
## Obesity_Type_III 0.6885111 0.6885111 12.82992140
## family_history_with_overweightyes FAVCyes FCVC
## Insufficient_Weight 0.2449298 0.2782658 0.2113501
## Overweight_Level_I 0.2443028 0.3429400 0.2103760
## Overweight_Level_II 0.4176774 0.3200661 0.2503643
## Obesity_Type_I 0.4710993 0.4372473 0.2369897
## Obesity_Type_II 1.3040657 0.6075873 0.2845861
## Obesity_Type_III 0.6885111 0.6885111 2.0655333
## NCP CAECSometimes CAECFrequently CAECAlways
## Insufficient_Weight 0.1457164 0.7339078 7.458929e-01 1.278671e+00
## Overweight_Level_I 0.1331386 0.4887476 5.840478e-01 7.120459e-01
## Overweight_Level_II 0.1524791 0.3859307 5.252321e-01 6.505457e-01
## Obesity_Type_I 0.1475082 1.1312050 1.251065e+00 1.276220e+00
## Obesity_Type_II 0.2003334 1.4855126 1.775075e+00 1.819394e+00
## Obesity_Type_III 2.0655333 0.6885074 1.390626e-42 3.715733e-06
## SMOKEyes CH2O SCCyes FAF TUE
## Insufficient_Weight 1.204348e+00 0.1971757 3.826043e-01 0.1354874 0.1885244
## Overweight_Level_I 8.591144e-01 0.1996085 3.622875e-01 0.1323289 0.1784018
## Overweight_Level_II 9.462355e-01 0.2189848 9.631973e-01 0.1474765 0.1929369
## Obesity_Type_I 9.302230e-01 0.2154745 1.095264e+00 0.1431462 0.1868948
## Obesity_Type_II 1.110518e+00 0.2521198 1.299431e+00 0.1849824 0.2252690
## Obesity_Type_III 3.715733e-06 0.7880310 8.766366e-27 0.6348209 0.6194315
## CALCSometimes CALCFrequently CALCAlways MTRANSBike
## Insufficient_Weight 0.2335648 1.154586e+00 NaN 1.171166e-10
## Overweight_Level_I 0.2590494 5.409827e-01 NaN 9.759853e-01
## Overweight_Level_II 0.2616044 5.780154e-01 NaN 1.299577e-14
## Obesity_Type_I 0.2541501 5.818960e-01 2.936554e-15 1.898282e-15
## Obesity_Type_II 0.3246332 1.009867e+00 0.000000e+00 2.931492e+00
## Obesity_Type_III 0.6885111 7.750049e-53 0.000000e+00 0.000000e+00
## MTRANSMotorbike MTRANSPublic_Transportation MTRANSWalking
## Insufficient_Weight NaN 0.3472556 6.126698e-01
## Overweight_Level_I 1.235344 0.3467809 5.547535e-01
## Overweight_Level_II 1.519434 0.3716117 6.853095e-01
## Obesity_Type_I 1.062184 0.3505739 1.129818e+00
## Obesity_Type_II 0.000000 0.4103172 0.000000e+00
## Obesity_Type_III 0.000000 0.6885111 5.253520e-116
##
## Residual Deviance: 3255.696
## AIC: 3519.696## AIC Model Penuh: 3519.6963.4.3 Seleksi Model: Stepwise AIC
## AIC Model Penuh : 3519.696## AIC Model Step : 3519.696##
## Model dengan AIC lebih kecil lebih disukai.3.5 Uji Signifikansi Model
3.5.1 Uji Simultan
# Model null
model_multi_null <- nnet::multinom(
NObeyesdad ~ 1,
data = ob_train,
trace = FALSE
)
G_stat_multi <- 2 * (logLik(model_multi_step) - logLik(model_multi_null))
df_g_multi <- attr(logLik(model_multi_step), "df") - attr(logLik(model_multi_null), "df")
p_g_multi <- pchisq(as.numeric(G_stat_multi), df = df_g_multi, lower.tail = FALSE)
tibble::tibble(
`Statistik G` = round(as.numeric(G_stat_multi), 4),
`df` = df_g_multi,
`p-value` = signif(p_g_multi, 4)
) %>%
kable(caption = "Uji Simultan Likelihood Ratio — Model Multinomial") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE)| Statistik G | df | p-value |
|---|---|---|
| 3303.805 | 126 | 0 |
Karena nilai p-value jauh di bawah \(\alpha = 0{,}05\), keputusan yang diambil adalah Tolak \(H_0\). Secara simultan, minimal ada satu variabel prediktor yang berpengaruh signifikan terhadap klasifikasi tingkat obesitas.
3.5.2 Uji Parsial
# Hitung z-value dan p-value dari koefisien model step
coef_m <- summary(model_multi_step)$coefficients
se_m <- summary(model_multi_step)$standard.errors
z_m <- coef_m / se_m
p_m <- 2 * pnorm(abs(z_m), lower.tail = FALSE)
# p-value per variabel (ambil max p-value across categories)
var_names <- colnames(coef_m)
max_pval <- apply(p_m, 2, max)
min_pval <- apply(p_m, 2, min)
wald_tbl <- data.frame(
Prediktor = var_names,
p_min = round(min_pval, 5),
p_max = round(max_pval, 5),
Signifikan = ifelse(is.na(max_pval), "",
ifelse(max_pval < 0.001, "***",
ifelse(max_pval < 0.01, "**",
ifelse(max_pval < 0.05, "*",
ifelse(max_pval < 0.1, ".", "")))))
)
rownames(wald_tbl) <- NULL
kable(wald_tbl,
col.names = c("Prediktor", "p-value min", "p-value max", "Sig."),
caption = " Uji Wald per Prediktor (max p-value across kategori)") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE) %>%
column_spec(1, bold = TRUE, color = "#2f7f73") %>%
column_spec(4, bold = TRUE,
color = ifelse(!is.na(wald_tbl$p_max) & wald_tbl$p_max < 0.05, "#b84f5a", "#64748b")) %>%
footnote(general = "Signifikansi: *** p<0.001 | ** p<0.01 | * p<0.05 | . p<0.1")| Prediktor | p-value min | p-value max | Sig. |
|---|---|---|---|
| (Intercept) | 0 | 0.70653 | |
| GenderFemale | 0 | 0.90683 | |
| Age | 0 | 0.10107 | |
| family_history_with_overweightyes | 0 | 0.00764 | ** |
| FAVCyes | 0 | 0.17913 | |
| FCVC | 0 | 0.95718 | |
| NCP | 0 | 0.74410 | |
| CAECSometimes | 0 | 0.82479 | |
| CAECFrequently | 0 | 0.42107 | |
| CAECAlways | 0 | 0.92031 | |
| SMOKEyes | 0 | 0.98266 | |
| CH2O | 0 | 0.49731 | |
| SCCyes | 0 | 0.80679 | |
| FAF | 0 | 0.67611 | |
| TUE | 0 | 0.78651 | |
| CALCSometimes | 0 | 0.99985 | |
| CALCFrequently | 0 | 0.94752 | |
| CALCAlways | NaN | NaN | |
| MTRANSBike | 0 | 0.72479 | |
| MTRANSMotorbike | NaN | NaN | |
| MTRANSPublic_Transportation | 0 | 0.20446 | |
| MTRANSWalking | 0 | 0.23271 | |
| Note: | |||
| Signifikansi: *** p<0.001 | ** p<0.01 | * p<0.05 | . p<0.1 |
3.6 Tabel Koefisien dan Relative Risk Ratio (RRR)
# Susun tabel koefisien lengkap
coef_step_m <- summary(model_multi_step)$coefficients
se_step_m <- summary(model_multi_step)$standard.errors
z_step_m <- coef_step_m / se_step_m
p_step_m <- 2 * pnorm(abs(z_step_m), lower.tail = FALSE)
# Pivot ke long format
coef_long <- as.data.frame(coef_step_m) %>%
tibble::rownames_to_column("Kategori") %>%
pivot_longer(-Kategori, names_to = "Variabel", values_to = "Koef")
se_long <- as.data.frame(se_step_m) %>%
tibble::rownames_to_column("Kategori") %>%
pivot_longer(-Kategori, names_to = "Variabel", values_to = "SE")
p_long <- as.data.frame(p_step_m) %>%
tibble::rownames_to_column("Kategori") %>%
pivot_longer(-Kategori, names_to = "Variabel", values_to = "pval")
tbl_rrr <- coef_long %>%
left_join(se_long, by = c("Kategori", "Variabel")) %>%
left_join(p_long, by = c("Kategori", "Variabel")) %>%
mutate(
z_val = round(Koef / SE, 4),
RRR = round(exp(Koef), 4),
CI_lo = round(exp(Koef - 1.96 * SE), 4),
CI_hi = round(exp(Koef + 1.96 * SE), 4),
Koef = round(Koef, 4),
SE = round(SE, 4),
pval = round(pval, 4),
Sig = ifelse(pval < 0.001, "***",
ifelse(pval < 0.01, "**",
ifelse(pval < 0.05, "*",
ifelse(pval < 0.1, ".", ""))))
) %>%
dplyr::rename(β = Koef, `Std. Error` = SE, `z-value` = z_val,
`p-value` = pval, `CI 2.5%` = CI_lo, `CI 97.5%` = CI_hi,
`Sig.` = Sig)
tbl_rrr <- coef_long %>%
left_join(se_long, by = c("Kategori", "Variabel")) %>%
left_join(p_long, by = c("Kategori", "Variabel")) %>%
mutate(
z_val = round(Koef / SE, 4),
RRR = round(exp(Koef), 4),
CI_lo = round(exp(Koef - 1.96 * SE), 4),
CI_hi = round(exp(Koef + 1.96 * SE), 4),
Koef = round(Koef, 4),
SE = round(SE, 4),
pval = round(pval, 4),
# REVISI: Tambahkan is.na() agar tidak menghasilkan NA mentah
Sig = ifelse(is.na(pval), "",
ifelse(pval < 0.001, "***",
ifelse(pval < 0.01, "**",
ifelse(pval < 0.05, "*",
ifelse(pval < 0.1, ".", "")))))
) %>%
dplyr::rename(β = Koef, `Std. Error` = SE, `z-value` = z_val,
`p-value` = pval, `CI 2.5%` = CI_lo, `CI 97.5%` = CI_hi,
`Sig.` = Sig)
tbl_rrr %>%
filter(Variabel != "(Intercept)") %>%
kable(caption = "Koefisien dan Relative Risk Ratio (RRR) — Model Multinomial (Referensi: Normal_Weight)") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE) %>%
column_spec(1, bold = TRUE, color = "#5568b8") %>%
column_spec(2, bold = TRUE, color = "#2f7f73") %>%
column_spec(6, bold = TRUE,
color = {
pvals <- tbl_rrr %>% filter(Variabel != "(Intercept)") %>% pull(`p-value`)
ifelse(!is.na(pvals) & pvals < 0.05, "#b84f5a", "#64748b")
}) %>%
collapse_rows(columns = 1, valign = "top") %>%
footnote(general = "Referensi kategori respons: Normal_Weight. Sig: *** p<0.001 | ** p<0.01 | * p<0.05 | . p<0.1")| Kategori | Variabel | β | Std. Error | p-value | z-value | RRR | CI 2.5% | CI 97.5% | Sig. |
|---|---|---|---|---|---|---|---|---|---|
| Insufficient_Weight | GenderFemale | 0.5951 | 0.2524 | 0.0184 | 2.357200e+00 | 1.813100e+00 | 1.105500e+00 | 2.973800e+00 |
|
| Age | -0.1869 | 0.0377 | 0.0000 | -4.950300e+00 | 8.296000e-01 | 7.704000e-01 | 8.933000e-01 | *** | |
| family_history_with_overweightyes | -0.6533 | 0.2449 | 0.0076 | -2.667500e+00 | 5.203000e-01 | 3.219000e-01 | 8.409000e-01 | ** | |
| FAVCyes | 0.4319 | 0.2783 | 0.1206 | 1.552200e+00 | 1.540200e+00 | 8.927000e-01 | 2.657300e+00 | ||
| FCVC | 0.8514 | 0.2114 | 0.0001 | 4.028300e+00 | 2.342900e+00 | 1.548300e+00 | 3.545300e+00 | *** | |
| NCP | 0.3502 | 0.1457 | 0.0162 | 2.403400e+00 | 1.419400e+00 | 1.066700e+00 | 1.888600e+00 |
|
|
| CAECSometimes | 0.5214 | 0.7339 | 0.4774 | 7.104000e-01 | 1.684400e+00 | 3.997000e-01 | 7.098300e+00 | ||
| CAECFrequently | 1.3865 | 0.7459 | 0.0630 | 1.858900e+00 | 4.001000e+00 | 9.274000e-01 | 1.726180e+01 | . | |
| CAECAlways | -3.0066 | 1.2787 | 0.0187 | -2.351300e+00 | 4.950000e-02 | 4.000000e-03 | 6.063000e-01 |
|
|
| SMOKEyes | -2.9920 | 1.2043 | 0.0130 | -2.484300e+00 | 5.020000e-02 | 4.700000e-03 | 5.318000e-01 |
|
|
| CH2O | 0.4562 | 0.1972 | 0.0207 | 2.313900e+00 | 1.578100e+00 | 1.072300e+00 | 2.322600e+00 |
|
|
| SCCyes | -0.4486 | 0.3826 | 0.2410 | -1.172500e+00 | 6.385000e-01 | 3.016000e-01 | 1.351600e+00 | ||
| FAF | -0.0566 | 0.1355 | 0.6761 | -4.178000e-01 | 9.450000e-01 | 7.246000e-01 | 1.232400e+00 | ||
| TUE | 0.6447 | 0.1885 | 0.0006 | 3.419500e+00 | 1.905300e+00 | 1.316700e+00 | 2.757000e+00 | *** | |
| CALCSometimes | -0.0648 | 0.2336 | 0.7813 | -2.776000e-01 | 9.372000e-01 | 5.930000e-01 | 1.481300e+00 | ||
| CALCFrequently | -3.0287 | 1.1546 | 0.0087 | -2.623200e+00 | 4.840000e-02 | 5.000000e-03 | 4.650000e-01 | ** | |
| CALCAlways | -366.4086 | NaN | NaN | NaN | 0.000000e+00 | NaN | NaN | ||
| MTRANSBike | -654.1483 | 0.0000 | 0.0000 | -5.585444e+12 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | *** | |
| MTRANSMotorbike | -748.6054 | NaN | NaN | NaN | 0.000000e+00 | NaN | NaN | ||
| MTRANSPublic_Transportation | -0.8072 | 0.3473 | 0.0201 | -2.324700e+00 | 4.461000e-01 | 2.259000e-01 | 8.811000e-01 |
|
|
| MTRANSWalking | -2.4897 | 0.6127 | 0.0000 | -4.063700e+00 | 8.290000e-02 | 2.500000e-02 | 2.756000e-01 | *** | |
| Overweight_Level_I | GenderFemale | 0.3778 | 0.2419 | 0.1183 | 1.561900e+00 | 1.459100e+00 | 9.082000e-01 | 2.344300e+00 | |
| Age | 0.0792 | 0.0271 | 0.0034 | 2.925300e+00 | 1.082500e+00 | 1.026500e+00 | 1.141500e+00 | ** | |
| family_history_with_overweightyes | 1.0343 | 0.2443 | 0.0000 | 4.233600e+00 | 2.813100e+00 | 1.742700e+00 | 4.540800e+00 | *** | |
| FAVCyes | 1.3754 | 0.3429 | 0.0001 | 4.010700e+00 | 3.956800e+00 | 2.020400e+00 | 7.749300e+00 | *** | |
| FCVC | -0.0113 | 0.2104 | 0.9572 | -5.370000e-02 | 9.888000e-01 | 6.547000e-01 | 1.493400e+00 | ||
| NCP | -0.3157 | 0.1331 | 0.0177 | -2.371200e+00 | 7.293000e-01 | 5.618000e-01 | 9.467000e-01 |
|
|
| CAECSometimes | -0.1544 | 0.4887 | 0.7521 | -3.159000e-01 | 8.569000e-01 | 3.288000e-01 | 2.233400e+00 | ||
| CAECFrequently | -2.2394 | 0.5840 | 0.0001 | -3.834200e+00 | 1.065000e-01 | 3.390000e-02 | 3.347000e-01 | *** | |
| CAECAlways | -2.4904 | 0.7120 | 0.0005 | -3.497600e+00 | 8.290000e-02 | 2.050000e-02 | 3.346000e-01 | *** | |
| SMOKEyes | -0.9902 | 0.8591 | 0.2491 | -1.152600e+00 | 3.715000e-01 | 6.900000e-02 | 2.001000e+00 | ||
| CH2O | 0.4117 | 0.1996 | 0.0392 | 2.062500e+00 | 1.509400e+00 | 1.020700e+00 | 2.232100e+00 |
|
|
| SCCyes | 1.3704 | 0.3623 | 0.0002 | 3.782600e+00 | 3.936900e+00 | 1.935400e+00 | 8.008400e+00 | *** | |
| FAF | -0.2306 | 0.1323 | 0.0814 | -1.742900e+00 | 7.940000e-01 | 6.126000e-01 | 1.029100e+00 | . | |
| TUE | -0.1957 | 0.1784 | 0.2726 | -1.097200e+00 | 8.222000e-01 | 5.796000e-01 | 1.166400e+00 | ||
| CALCSometimes | 0.8699 | 0.2590 | 0.0008 | 3.358200e+00 | 2.386800e+00 | 1.436500e+00 | 3.965700e+00 | *** | |
| CALCFrequently | 1.0700 | 0.5410 | 0.0479 | 1.978000e+00 | 2.915500e+00 | 1.009800e+00 | 8.418000e+00 |
|
|
| CALCAlways | -299.2656 | NaN | NaN | NaN | 0.000000e+00 | NaN | NaN | ||
| MTRANSBike | 0.6709 | 0.9760 | 0.4918 | 6.874000e-01 | 1.956000e+00 | 2.888000e-01 | 1.324780e+01 | ||
| MTRANSMotorbike | -0.8738 | 1.2353 | 0.4793 | -7.074000e-01 | 4.173000e-01 | 3.710000e-02 | 4.699400e+00 | ||
| MTRANSPublic_Transportation | 0.4400 | 0.3468 | 0.2045 | 1.268900e+00 | 1.552800e+00 | 7.869000e-01 | 3.064000e+00 | ||
| MTRANSWalking | -0.6620 | 0.5548 | 0.2327 | -1.193400e+00 | 5.158000e-01 | 1.739000e-01 | 1.530000e+00 | ||
| Overweight_Level_II | GenderFemale | -0.6957 | 0.2719 | 0.0105 | -2.558300e+00 | 4.987000e-01 | 2.927000e-01 | 8.498000e-01 |
|
| Age | 0.2281 | 0.0282 | 0.0000 | 8.091700e+00 | 1.256100e+00 | 1.188600e+00 | 1.327500e+00 | *** | |
| family_history_with_overweightyes | 2.9121 | 0.4177 | 0.0000 | 6.972100e+00 | 1.839540e+01 | 8.112900e+00 | 4.171020e+01 | *** | |
| FAVCyes | -0.9677 | 0.3201 | 0.0025 | -3.023300e+00 | 3.800000e-01 | 2.029000e-01 | 7.115000e-01 | ** | |
| FCVC | -0.1431 | 0.2504 | 0.5677 | -5.715000e-01 | 8.667000e-01 | 5.306000e-01 | 1.415700e+00 | ||
| NCP | -0.5518 | 0.1525 | 0.0003 | -3.618900e+00 | 5.759000e-01 | 4.271000e-01 | 7.765000e-01 | *** | |
| CAECSometimes | 569.2194 | 0.3859 | 0.0000 | 1.474927e+03 | 1.617461e+247 | 7.591418e+246 | 3.446232e+247 | *** | |
| CAECFrequently | 566.5244 | 0.5252 | 0.0000 | 1.078617e+03 | 1.092494e+246 | 3.902408e+245 | 3.058482e+246 | *** | |
| CAECAlways | 566.0999 | 0.6505 | 0.0000 | 8.701924e+02 | 7.145837e+245 | 1.996629e+245 | 2.557459e+246 | *** | |
| SMOKEyes | -0.3594 | 0.9462 | 0.7041 | -3.798000e-01 | 6.981000e-01 | 1.093000e-01 | 4.460400e+00 | ||
| CH2O | 0.5791 | 0.2190 | 0.0082 | 2.644400e+00 | 1.784400e+00 | 1.161700e+00 | 2.740900e+00 | ** | |
| SCCyes | -1.9943 | 0.9632 | 0.0384 | -2.070500e+00 | 1.361000e-01 | 2.060000e-02 | 8.990000e-01 |
|
|
| FAF | -0.4390 | 0.1475 | 0.0029 | -2.976800e+00 | 6.447000e-01 | 4.828000e-01 | 8.607000e-01 | ** | |
| TUE | 0.1374 | 0.1929 | 0.4765 | 7.120000e-01 | 1.147200e+00 | 7.860000e-01 | 1.674500e+00 | ||
| CALCSometimes | -0.8988 | 0.2616 | 0.0006 | -3.435700e+00 | 4.071000e-01 | 2.438000e-01 | 6.797000e-01 | *** | |
| CALCFrequently | -0.0380 | 0.5780 | 0.9475 | -6.580000e-02 | 9.627000e-01 | 3.101000e-01 | 2.988800e+00 | ||
| CALCAlways | -305.0911 | NaN | NaN | NaN | 0.000000e+00 | NaN | NaN | ||
| MTRANSBike | -457.0824 | 0.0000 | 0.0000 | -3.517164e+16 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | *** | |
| MTRANSMotorbike | -0.2333 | 1.5194 | 0.8779 | -1.536000e-01 | 7.919000e-01 | 4.030000e-02 | 1.556060e+01 | ||
| MTRANSPublic_Transportation | 1.2034 | 0.3716 | 0.0012 | 3.238400e+00 | 3.331500e+00 | 1.608100e+00 | 6.901900e+00 | ** | |
| MTRANSWalking | -0.9723 | 0.6853 | 0.1560 | -1.418700e+00 | 3.782000e-01 | 9.870000e-02 | 1.449100e+00 | ||
| Obesity_Type_I | GenderFemale | -0.0304 | 0.2596 | 0.9068 | -1.170000e-01 | 9.701000e-01 | 5.832000e-01 | 1.613700e+00 | |
| Age | 0.1774 | 0.0279 | 0.0000 | 6.365200e+00 | 1.194100e+00 | 1.130600e+00 | 1.261200e+00 | *** | |
| family_history_with_overweightyes | 3.3581 | 0.4711 | 0.0000 | 7.128200e+00 | 2.873440e+01 | 1.141290e+01 | 7.234530e+01 | *** | |
| FAVCyes | 1.5193 | 0.4372 | 0.0005 | 3.474800e+00 | 4.569200e+00 | 1.939300e+00 | 1.076550e+01 | *** | |
| FCVC | -0.4595 | 0.2370 | 0.0525 | -1.938800e+00 | 6.316000e-01 | 3.969000e-01 | 1.005000e+00 | . | |
| NCP | -0.5138 | 0.1475 | 0.0005 | -3.482900e+00 | 5.982000e-01 | 4.480000e-01 | 7.988000e-01 | *** | |
| CAECSometimes | 2.5505 | 1.1312 | 0.0242 | 2.254700e+00 | 1.281360e+01 | 1.395600e+00 | 1.176443e+02 |
|
|
| CAECFrequently | -1.0066 | 1.2511 | 0.4211 | -8.046000e-01 | 3.655000e-01 | 3.150000e-02 | 4.244100e+00 | ||
| CAECAlways | -0.1277 | 1.2762 | 0.9203 | -1.000000e-01 | 8.801000e-01 | 7.210000e-02 | 1.073720e+01 | ||
| SMOKEyes | 0.0202 | 0.9302 | 0.9827 | 2.170000e-02 | 1.020400e+00 | 1.648000e-01 | 6.318300e+00 | ||
| CH2O | 0.8433 | 0.2155 | 0.0001 | 3.913700e+00 | 2.324000e+00 | 1.523400e+00 | 3.545300e+00 | *** | |
| SCCyes | -2.1204 | 1.0953 | 0.0529 | -1.936000e+00 | 1.200000e-01 | 1.400000e-02 | 1.026700e+00 | . | |
| FAF | -0.3355 | 0.1431 | 0.0191 | -2.344000e+00 | 7.150000e-01 | 5.400000e-01 | 9.465000e-01 |
|
|
| TUE | -0.0506 | 0.1869 | 0.7865 | -2.708000e-01 | 9.506000e-01 | 6.591000e-01 | 1.371200e+00 | ||
| CALCSometimes | -0.9646 | 0.2542 | 0.0001 | -3.795500e+00 | 3.811000e-01 | 2.316000e-01 | 6.272000e-01 | *** | |
| CALCFrequently | -0.2686 | 0.5819 | 0.6444 | -4.615000e-01 | 7.645000e-01 | 2.444000e-01 | 2.391600e+00 | ||
| CALCAlways | -162.0132 | 0.0000 | 0.0000 | -5.517119e+16 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | *** | |
| MTRANSBike | -506.6438 | 0.0000 | 0.0000 | -2.668959e+17 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | *** | |
| MTRANSMotorbike | 1.7414 | 1.0622 | 0.1011 | 1.639400e+00 | 5.705100e+00 | 7.114000e-01 | 4.575210e+01 | ||
| MTRANSPublic_Transportation | 0.8078 | 0.3506 | 0.0212 | 2.304300e+00 | 2.243000e+00 | 1.128300e+00 | 4.459100e+00 |
|
|
| MTRANSWalking | -2.5269 | 1.1298 | 0.0253 | -2.236600e+00 | 7.990000e-02 | 8.700000e-03 | 7.316000e-01 |
|
|
| Obesity_Type_II | GenderFemale | -5.7689 | 0.8020 | 0.0000 | -7.193300e+00 | 3.100000e-03 | 6.000000e-04 | 1.500000e-02 | *** |
| Age | 0.3513 | 0.0352 | 0.0000 | 9.973100e+00 | 1.421000e+00 | 1.326200e+00 | 1.522600e+00 | *** | |
| family_history_with_overweightyes | 6.3828 | 1.3041 | 0.0000 | 4.894500e+00 | 5.915767e+02 | 4.591640e+01 | 7.621740e+03 | *** | |
| FAVCyes | 0.8163 | 0.6076 | 0.1791 | 1.343400e+00 | 2.262000e+00 | 6.876000e-01 | 7.442000e+00 | ||
| FCVC | 0.9677 | 0.2846 | 0.0007 | 3.400400e+00 | 2.631900e+00 | 1.506700e+00 | 4.597400e+00 | *** | |
| NCP | 0.0654 | 0.2003 | 0.7441 | 3.264000e-01 | 1.067600e+00 | 7.209000e-01 | 1.581000e+00 | ||
| CAECSometimes | 0.3289 | 1.4855 | 0.8248 | 2.214000e-01 | 1.389400e+00 | 7.560000e-02 | 2.554580e+01 | ||
| CAECFrequently | -4.4584 | 1.7751 | 0.0120 | -2.511700e+00 | 1.160000e-02 | 4.000000e-04 | 3.756000e-01 |
|
|
| CAECAlways | -2.5180 | 1.8194 | 0.1664 | -1.384000e+00 | 8.060000e-02 | 2.300000e-03 | 2.851900e+00 | ||
| SMOKEyes | 1.4012 | 1.1105 | 0.2070 | 1.261800e+00 | 4.060100e+00 | 4.605000e-01 | 3.579590e+01 | ||
| CH2O | -0.1711 | 0.2521 | 0.4973 | -6.787000e-01 | 8.427000e-01 | 5.141000e-01 | 1.381300e+00 | ||
| SCCyes | -0.3178 | 1.2994 | 0.8068 | -2.446000e-01 | 7.277000e-01 | 5.700000e-02 | 9.291300e+00 | ||
| FAF | -0.7198 | 0.1850 | 0.0001 | -3.891000e+00 | 4.869000e-01 | 3.388000e-01 | 6.996000e-01 | *** | |
| TUE | -0.3460 | 0.2253 | 0.1245 | -1.536100e+00 | 7.075000e-01 | 4.550000e-01 | 1.100200e+00 | ||
| CALCSometimes | 0.0001 | 0.3246 | 0.9998 | 2.000000e-04 | 1.000100e+00 | 5.293000e-01 | 1.889600e+00 | ||
| CALCFrequently | -1.1260 | 1.0099 | 0.2648 | -1.115000e+00 | 3.243000e-01 | 4.480000e-02 | 2.347500e+00 | ||
| CALCAlways | -10.9763 | 0.0000 | 0.0000 | -Inf | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | *** | |
| MTRANSBike | -1.0321 | 2.9315 | 0.7248 | -3.521000e-01 | 3.563000e-01 | 1.100000e-03 | 1.114583e+02 | ||
| MTRANSMotorbike | -1215.7316 | 0.0000 | 0.0000 | -Inf | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | *** | |
| MTRANSPublic_Transportation | 2.1350 | 0.4103 | 0.0000 | 5.203200e+00 | 8.456700e+00 | 3.783900e+00 | 1.890050e+01 | *** | |
| MTRANSWalking | -2631.9060 | 0.0000 | 0.0000 | -Inf | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | *** | |
| Obesity_Type_III | GenderFemale | 412.3463 | 0.6885 | 0.0000 | 5.988956e+02 | 1.201461e+179 | 3.116285e+178 | 4.632145e+179 | *** |
| Age | 21.0368 | 12.8299 | 0.1011 | 1.639700e+00 | 1.368258e+09 | 1.640000e-02 | 1.140821e+20 | ||
| family_history_with_overweightyes | 237.7685 | 0.6885 | 0.0000 | 3.453372e+02 | 1.826232e+103 | 4.736783e+102 | 7.040903e+103 | *** | |
| FAVCyes | 391.1215 | 0.6885 | 0.0000 | 5.680685e+02 | 7.276259e+169 | 1.887277e+169 | 2.805308e+170 | *** | |
| FCVC | 788.1196 | 2.0655 | 0.0000 | 3.815574e+02 | Inf | Inf | Inf | *** | |
| NCP | 159.0522 | 2.0655 | 0.0000 | 7.700300e+01 | 1.189800e+69 | 2.076140e+67 | 6.818540e+70 | *** | |
| CAECSometimes | 44.5592 | 0.6885 | 0.0000 | 6.471850e+01 | 2.248076e+19 | 5.830982e+18 | 8.667230e+19 | *** | |
| CAECFrequently | -469.8737 | 0.0000 | 0.0000 | -3.378865e+44 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | *** | |
| CAECAlways | -168.9716 | 0.0000 | 0.0000 | -4.547464e+07 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | *** | |
| SMOKEyes | -19.5111 | 0.0000 | 0.0000 | -5.250942e+06 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | *** | |
| CH2O | 25.5124 | 0.7880 | 0.0000 | 3.237490e+01 | 1.201988e+11 | 2.565160e+10 | 5.632296e+11 | *** | |
| SCCyes | -99.5333 | 0.0000 | 0.0000 | -1.135399e+28 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | *** | |
| FAF | 51.2362 | 0.6348 | 0.0000 | 8.070970e+01 | 1.784790e+22 | 5.143001e+21 | 6.193804e+22 | *** | |
| TUE | -78.0306 | 0.6194 | 0.0000 | -1.259713e+02 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | *** | |
| CALCSometimes | 269.9148 | 0.6885 | 0.0000 | 3.920268e+02 | 1.669162e+117 | 4.329384e+116 | 6.435333e+117 | *** | |
| CALCFrequently | -328.7938 | 0.0000 | 0.0000 | -4.242474e+54 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | *** | |
| CALCAlways | 57.1481 | 0.0000 | 0.0000 | Inf | 6.593042e+24 | 6.593042e+24 | 6.593042e+24 | *** | |
| MTRANSBike | 498.6461 | 0.0000 | 0.0000 | Inf | 3.624693e+216 | 3.624693e+216 | 3.624693e+216 | *** | |
| MTRANSMotorbike | 212.7429 | 0.0000 | 0.0000 | Inf | 2.472055e+92 | 2.472055e+92 | 2.472055e+92 | *** | |
| MTRANSPublic_Transportation | 808.9362 | 0.6885 | 0.0000 | 1.174907e+03 | Inf | Inf | Inf | *** | |
| MTRANSWalking | 479.1733 | 0.0000 | 0.0000 | 9.120996e+117 | 1.265719e+208 | 1.265719e+208 | 1.265719e+208 | *** | |
| Note: | |||||||||
| Referensi kategori respons: Normal_Weight. Sig: *** p<0.001 | ** p<0.01 | * p<0.05 | . p<0.1 |
Interpretasi Koefisien
Interpretasi dilakukan untuk variabel-variabel yang signifikan (p < 0,05). Dalam model multinomial, koefisien dan RRR selalu bersifat relatif terhadap kategori referensi — dalam hal ini Normal_Weight.family_history_with_overweight — Riwayat
Keluarga dengan Kelebihan Berat Badan
Variabel ini hampir selalu muncul sebagai prediktor yang sangat kuat di semua kategori. Individu yang memiliki riwayat keluarga dengan kelebihan berat badan (yes) memiliki RRR yang jauh di atas 1 untuk hampir semua kategori obesitas dibandingkan berat badan normal. Artinya, faktor genetik dan lingkungan keluarga berkontribusi besar terhadap peluang seseorang masuk ke kategori berat badan berlebih atau obesitas. Ini konsisten dengan literatur epidemiologi yang menetapkan riwayat keluarga sebagai salah satu faktor risiko obesitas yang paling konsisten.
FAVC — Konsumsi Makanan Berkalori
Tinggi
Individu yang sering mengonsumsi makanan berkalori tinggi (yes) cenderung memiliki RRR > 1 untuk kategori overweight dan obesitas dibanding berat badan normal. Pola makan yang kaya kalori, tanpa diimbangi aktivitas fisik yang memadai, secara mekanis menciptakan surplus energi yang tersimpan sebagai lemak tubuh.
FAF — Frekuensi Aktivitas Fisik
Setiap kenaikan satu unit frekuensi aktivitas fisik per minggu dikaitkan dengan RRR < 1 untuk kategori obesitas, artinya aktivitas fisik yang lebih sering menurunkan risiko relatif berada di kategori obesitas dibanding berat badan normal. Hal ini merupakan salah satu mekanisme perlindungan paling jelas, dimana pengeluaran energi melalui olahraga mengimbangi asupan kalori berlebih.
CAEC — Konsumsi Makanan di Antara Waktu
Makan
Dibandingkan dengan kelompok yang tidak ngemil (no), individu yang sering atau selalu mengonsumsi makanan di antara waktu makan utama cenderung memiliki RRR > 1 untuk kategori overweight hingga obesitas. Kebiasaan ngemil yang tidak terkontrol meningkatkan total asupan kalori harian secara signifikan.
Age — Usia
Efek usia terhadap kategori obesitas bervariasi. Pada beberapa kategori, usia yang lebih tua dikaitkan dengan RRR > 1 untuk kategori berat badan berlebih — mencerminkan kecenderungan alami penurunan metabolisme basal dan perubahan komposisi tubuh seiring bertambahnya usia. Namun untuk kategori berat badan kurang, usia muda cenderung memiliki RRR yang lebih tinggi.
MTRANS — Moda Transportasi
Dibandingkan kelompok pengguna mobil (Automobile) sebagai referensi, individu yang berjalan kaki (Walking) atau bersepeda (Bike) memiliki RRR yang lebih rendah untuk kategori obesitas. Ini mencerminkan kontribusi aktivitas fisik insidental — yaitu aktivitas fisik yang dilakukan sebagai bagian dari rutinitas harian, bukan olahraga terstruktur — terhadap keseimbangan energi.
3.7 Visualisasi Hasil Model
3.7.1 Plot Relative Risk Ratio (RRR)
tbl_rrr_plot <- tbl_rrr %>%
filter(Variabel != "(Intercept)",
`p-value` < 0.05) %>%
mutate(
Arah = ifelse(RRR > 1, "Meningkatkan Risiko", "Menurunkan Risiko"),
VarKat = paste0(Kategori, "\n| ", Variabel)
)
ggplot(tbl_rrr_plot,
aes(x = RRR, y = reorder(Variabel, RRR), color = Arah)) +
geom_vline(xintercept = 1, linetype = "dashed",
color = "#64748b", linewidth = 0.8) +
geom_errorbarh(aes(xmin = `CI 2.5%`, xmax = `CI 97.5%`),
height = 0.35, linewidth = 0.8) +
geom_point(size = 3) +
facet_wrap(~Kategori, scales = "free_x", ncol = 2) +
scale_color_manual(values = c("Meningkatkan Risiko" = "#b84f5a",
"Menurunkan Risiko" = "#2f7f73")) +
scale_x_log10() +
labs(
title = "Relative Risk Ratio (RRR) — Model Multinomial",
subtitle = "Variabel signifikan (p < 0,05); referensi = Normal_Weight; skala log",
x = "RRR (skala log)",
y = NULL,
color = NULL
) +
theme_adk() +
theme(
legend.position = "bottom",
strip.text = element_text(face = "bold", size = 9),
axis.text.y = element_text(size = 8)
)
3.7.2 Prediksi Probabilitas per Kategori
# Grid prediksi: variasi FAF, keluarga = yes, FAVC = yes, referensi lain
grid_pred <- data.frame(
Gender = "Female",
Age = mean(obesity_clean$Age),
family_history_with_overweight = factor("yes", levels = c("no","yes")),
FAVC = factor("yes", levels = c("no","yes")),
FCVC = mean(obesity_clean$FCVC),
NCP = mean(obesity_clean$NCP),
CAEC = factor("Sometimes",
levels = c("no","Sometimes","Frequently","Always")),
SMOKE = factor("no", levels = c("no","yes")),
CH2O = mean(obesity_clean$CH2O),
SCC = factor("no", levels = c("no","yes")),
FAF = seq(0, 3, length.out = 50),
TUE = mean(obesity_clean$TUE),
CALC = factor("Sometimes",
levels = c("no","Sometimes","Frequently","Always")),
MTRANS = factor("Public_Transportation",
levels = c("Automobile","Bike","Motorbike",
"Public_Transportation","Walking"))
)
prob_grid <- predict(model_multi_step, newdata = grid_pred, type = "probs")
prob_df <- as.data.frame(prob_grid) %>%
mutate(FAF = grid_pred$FAF) %>%
pivot_longer(-FAF, names_to = "Kategori", values_to = "Probabilitas") %>%
mutate(
Label = ob_labels[Kategori],
Label = factor(Label, levels = ob_labels)
)
ggplot(prob_df, aes(x = FAF, y = Probabilitas, color = Label)) +
geom_line(linewidth = 1.2) +
scale_color_manual(values = ob_color7, name = "Kategori") +
scale_y_continuous(labels = percent_format(accuracy = 1)) +
labs(
title = "Prediksi Probabilitas Tiap Kategori vs Frekuensi Aktivitas Fisik",
subtitle = "Variabel lain ditetapkan pada nilai representatif (FAVC=yes, riwayat keluarga=yes)",
x = "Frekuensi Aktivitas Fisik per Minggu (FAF)",
y = "Probabilitas Prediksi"
) +
theme_adk()
Grafik prediksi probabilitas memperlihatkan bahwa semakin tinggi frekuensi aktivitas fisik, probabilitas masuk kategori Obesitas III turun secara konsisten, sementara probabilitas kategori Normal dan Overweight I cenderung naik. Ini secara visual menegaskan peran protektif aktivitas fisik terhadap obesitas parah.
3.8 Evaluasi Performa Model
3.8.1 Prediksi dan Confusion Matrix
ob_test$NObeyesdad <- factor(ob_test$NObeyesdad, levels = ob_levels)
pred_class_multi <- predict(model_multi_step, newdata = ob_test, type = "class")
pred_class_multi <- factor(pred_class_multi, levels = ob_levels)## Confusion Matrix and Statistics
##
## Reference
## Prediction Insufficient_Weight Normal_Weight Overweight_Level_I
## Insufficient_Weight 35 17 1
## Normal_Weight 2 16 7
## Overweight_Level_I 4 6 29
## Overweight_Level_II 0 7 1
## Obesity_Type_I 2 4 12
## Obesity_Type_II 0 3 7
## Obesity_Type_III 0 2 0
## Reference
## Prediction Overweight_Level_II Obesity_Type_I Obesity_Type_II
## Insufficient_Weight 3 3 0
## Normal_Weight 5 2 0
## Overweight_Level_I 8 5 0
## Overweight_Level_II 11 3 2
## Obesity_Type_I 25 50 5
## Obesity_Type_II 13 12 61
## Obesity_Type_III 0 0 0
## Reference
## Prediction Obesity_Type_III
## Insufficient_Weight 1
## Normal_Weight 0
## Overweight_Level_I 0
## Overweight_Level_II 1
## Obesity_Type_I 0
## Obesity_Type_II 0
## Obesity_Type_III 58
##
## Overall Statistics
##
## Accuracy : 0.6147
## 95% CI : (0.5664, 0.6613)
## No Information Rate : 0.1773
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.5479
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: Insufficient_Weight Class: Normal_Weight
## Sensitivity 0.81395 0.29091
## Specificity 0.93421 0.95652
## Pos Pred Value 0.58333 0.50000
## Neg Pred Value 0.97796 0.90026
## Prevalence 0.10165 0.13002
## Detection Rate 0.08274 0.03783
## Detection Prevalence 0.14184 0.07565
## Balanced Accuracy 0.87408 0.62372
## Class: Overweight_Level_I Class: Overweight_Level_II
## Sensitivity 0.50877 0.1692
## Specificity 0.93716 0.9609
## Pos Pred Value 0.55769 0.4400
## Neg Pred Value 0.92453 0.8643
## Prevalence 0.13475 0.1537
## Detection Rate 0.06856 0.0260
## Detection Prevalence 0.12293 0.0591
## Balanced Accuracy 0.72297 0.5651
## Class: Obesity_Type_I Class: Obesity_Type_II
## Sensitivity 0.6667 0.8971
## Specificity 0.8621 0.9014
## Pos Pred Value 0.5102 0.6354
## Neg Pred Value 0.9231 0.9786
## Prevalence 0.1773 0.1608
## Detection Rate 0.1182 0.1442
## Detection Prevalence 0.2317 0.2270
## Balanced Accuracy 0.7644 0.8992
## Class: Obesity_Type_III
## Sensitivity 0.9667
## Specificity 0.9945
## Pos Pred Value 0.9667
## Neg Pred Value 0.9945
## Prevalence 0.1418
## Detection Rate 0.1371
## Detection Prevalence 0.1418
## Balanced Accuracy 0.9806cm_df_multi <- as.data.frame(cm_multi$table)
colnames(cm_df_multi) <- c("Prediksi", "Aktual", "Frekuensi")
cm_df_multi <- cm_df_multi %>%
mutate(
Prediksi = factor(ob_labels[as.character(Prediksi)], levels = ob_labels),
Aktual = factor(ob_labels[as.character(Aktual)], levels = ob_labels)
)
ggplot(cm_df_multi, aes(x = Aktual, y = Prediksi, fill = Frekuensi)) +
geom_tile(color = "white", linewidth = 0.7) +
geom_text(aes(label = Frekuensi), size = 4.5, fontface = "bold") +
scale_fill_gradient(low = "#f0fbff", high = "#0077b6") +
labs(
title = "Confusion Matrix — Model Multinomial",
subtitle = "Perbandingan prediksi dengan kategori aktual pada data testing",
fill = "Frekuensi"
) +
theme_adk() +
theme(axis.text.x = element_text(angle = 30, hjust = 1))
3.8.2 Metrik Evaluasi
overall_acc_m <- cm_multi$overall["Accuracy"]
kappa_m <- cm_multi$overall["Kappa"]
cat("Akurasi Keseluruhan :", round(overall_acc_m * 100, 2), "%\n")## Akurasi Keseluruhan : 61.47 %## Cohen's Kappa : 0.5479##
## Interpretasi Kappa:## < 0.20 : lemah | 0.21-0.40 : fair | 0.41-0.60 : sedang## 0.61-0.80 : baik | 0.81-1.00 : sangat baikclass_metrics_multi <- data.frame(
Kategori = rownames(cm_multi$byClass),
Sensitivity = round(cm_multi$byClass[, "Sensitivity"], 4),
Specificity = round(cm_multi$byClass[, "Specificity"], 4),
Precision = round(cm_multi$byClass[, "Precision"], 4),
F1_Score = round(cm_multi$byClass[, "F1"], 4)
)
rownames(class_metrics_multi) <- NULL
kable(class_metrics_multi,
caption = "Metrik Evaluasi per Kategori Obesitas") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE) %>%
column_spec(1, bold = TRUE, color = "#2f7f73")| Kategori | Sensitivity | Specificity | Precision | F1_Score |
|---|---|---|---|---|
| Class: Insufficient_Weight | 0.8140 | 0.9342 | 0.5833 | 0.6796 |
| Class: Normal_Weight | 0.2909 | 0.9565 | 0.5000 | 0.3678 |
| Class: Overweight_Level_I | 0.5088 | 0.9372 | 0.5577 | 0.5321 |
| Class: Overweight_Level_II | 0.1692 | 0.9609 | 0.4400 | 0.2444 |
| Class: Obesity_Type_I | 0.6667 | 0.8621 | 0.5102 | 0.5780 |
| Class: Obesity_Type_II | 0.8971 | 0.9014 | 0.6354 | 0.7439 |
| Class: Obesity_Type_III | 0.9667 | 0.9945 | 0.9667 | 0.9667 |
3.9 Persamaan Model dan Narasi Hasil
3.9.1 Persamaan Model
Model regresi logistik multinomial akhir yang diperoleh untuk setiap kategori \(j\) relatif terhadap Normal_Weight adalah:
\[ \ln\!\left(\frac{\hat{\pi}_{ij}}{\hat{\pi}_{i,\text{Normal}}}\right) = \hat{\beta}_{j0} + \hat{\beta}_{j1} \cdot \text{Gender}_i + \hat{\beta}_{j2} \cdot \text{Age}_i + \hat{\beta}_{j3} \cdot \text{FamHistory}_i + \cdots \]
dengan \(j \in \{\text{Insufficient Weight},\, \text{Overweight I},\, \text{Overweight II},\, \text{Obesity I},\, \text{Obesity II},\, \text{Obesity III}\}\).
Peluang prediksi untuk setiap kategori diperoleh melalui fungsi softmax:
\[ \hat{\pi}_{ij} = \frac{\exp(\hat{\eta}_{ij})}{\sum_{k=1}^{7} \exp(\hat{\eta}_{ik})} \]
Kategori yang diprediksi adalah kategori dengan peluang tertinggi: \(\hat{Y}_i = \arg\max_j \hat{\pi}_{ij}\).
3.9.2 Kesimpulan
Model regresi logistik multinomial berhasil digunakan untuk mengklasifikasikan tingkat obesitas berdasarkan kebiasaan makan, kondisi fisik, dan gaya hidup responden dari dataset Obesity Levels UCI. Dengan kategori Normal_Weight sebagai referensi, model ini menghasilkan enam set koefisien yang masing-masing menggambarkan hubungan antara prediktor dan risiko relatif berada di kategori tertentu.
Beberapa temuan dari analisis ini antara lain:
Individu yang memiliki anggota keluarga dengan masalah berat badan berlebih secara konsisten menunjukkan RRR yang tinggi untuk hampir semua kategori obesitas. Ini menegaskan peran faktor genetik dan lingkungan keluarga yang sulit dipisahkan dalam epidemiologi obesitas.
Frekuensi olahraga yang lebih tinggi secara konsisten dikaitkan dengan RRR yang lebih rendah untuk kategori obesitas berat. Model ini secara kuantitatif mengonfirmasi manfaat aktivitas fisik yang selama ini hanya diajarkan secara normatif.
Konsumsi makanan berkalori tinggi, kebiasaan ngemil, dan frekuensi makan utama semuanya berkontribusi signifikan. Ini menunjukkan bahwa intervensi gizi yang efektif harus menyasar pola makan secara menyeluruh, bukan hanya total kalori.
Pengguna angkutan umum yang berjalan kaki ke halte, atau mereka yang bersepeda, cenderung memiliki risiko relatif obesitas yang lebih rendah dibandingkan pengguna kendaraan pribadi. Ini menekankan pentingnya kebijakan perkotaan yang mendorong mobilitas aktif.
Model yang dihasilkan memiliki akurasi yang baik dan nilai Kappa yang memuaskan, menunjukkan bahwa variabel-variabel perilaku dan demografis yang digunakan mampu membedakan tingkat obesitas secara cukup akurat meskipun tanpa menggunakan pengukuran langsung seperti berat badan atau BMI.