Analisis Regresi Logistik
Biner | Multinomial | Ordinal | Poisson
Shuffiatul Rahmi
2026-05-31
1 Regresi Logistik Biner
1.1 Pendahuluan & Deskripsi Dataset
Regresi logistik biner merupakan metode statistika yang digunakan untuk memodelkan hubungan antara satu atau lebih variabel prediktor dengan variabel respons dikotomi (dua kategori: 0/1). Metode ini banyak digunakan dalam bidang kesehatan untuk memprediksi kemungkinan seseorang mengalami suatu penyakit berdasarkan faktor-faktor risiko tertentu.
Pada penelitian ini digunakan Diabetes Prediction Dataset untuk memodelkan status diabetes pasien berdasarkan variabel kesehatan dan karakteristik individu. Variabel respons berupa status diabetes dengan dua kategori, yaitu “Diabetes” dan “Tidak Diabetes”.
Sumber data: Diabetes Prediction Dataset – tersedia melalui Mendeley Data:
https://data.mendeley.com/datasets/pd64hfttwy/1
Variabel Penelitian
| Variabel | Tipe | Keterangan | Satuan/Kategori |
|---|---|---|---|
| Outcome | Respons | Status diabetes | 0 = Tidak Diabetes, 1 = Diabetes |
| Pregnancies | Numerik | Jumlah kehamilan | Count |
| Glucose | Numerik | Kadar glukosa | mg/dL |
| BloodPressure | Numerik | Tekanan darah | mmHg |
| SkinThickness | Numerik | Ketebalan kulit | mm |
| Insulin | Numerik | Kadar insulin | mu U/ml |
| BMI | Numerik | Body Mass Index | kg/m^2 |
| DiabetesPedigreeFunction | Numerik | Riwayat diabetes keluarga | Skor |
| Age | Numerik | Umur pasien | Tahun |
1.2 Pemuatan Data
# Ganti path sesuai lokasi file Anda
diabetes_path <- "C:/Users/hp/Downloads/diabetes.csv"
if (!file.exists(diabetes_path)) {
tmp <- tempfile(fileext=".csv")
download.file("https://raw.githubusercontent.com/plotly/datasets/master/diabetes.csv", tmp, quiet=TRUE)
diabetes_path <- tmp
}
diabetes <- read.csv(diabetes_path)
# Transformasi data
diabetes <- diabetes %>%
mutate(
diabetes_status = ifelse(Outcome == 1, "Diabetes", "Tidak Diabetes"),
diabetes_status = factor(diabetes_status,
levels = c("Tidak Diabetes", "Diabetes")),
log_Glucose = log(Glucose + 1),
log_Age = log(Age + 1)
)
str(diabetes)## 'data.frame': 2251 obs. of 12 variables:
## $ Pregnancies : int 2 0 0 0 1 0 4 8 2 2 ...
## $ Glucose : int 138 84 145 135 139 173 99 194 83 89 ...
## $ BloodPressure : int 62 82 0 68 62 78 72 80 65 90 ...
## $ SkinThickness : int 35 31 0 42 41 32 17 0 28 30 ...
## $ Insulin : int 0 125 0 250 480 265 0 0 66 0 ...
## $ BMI : num 33.6 38.2 44.2 42.3 40.7 46.5 25.6 26.1 36.8 33.5 ...
## $ DiabetesPedigreeFunction: num 0.127 0.233 0.63 0.365 0.536 ...
## $ Age : int 47 23 31 24 21 58 28 67 24 42 ...
## $ Outcome : int 1 0 1 1 0 0 0 0 0 0 ...
## $ diabetes_status : Factor w/ 2 levels "Tidak Diabetes",..: 2 1 2 2 1 1 1 1 1 1 ...
## $ log_Glucose : num 4.93 4.44 4.98 4.91 4.94 ...
## $ log_Age : num 3.87 3.18 3.47 3.22 3.09 ...## [1] 2251 121.3 Eksplorasi Data (EDA)
1.3.1 Statistik Deskriptif
## Pregnancies Glucose BloodPressure BMI
## Min. : 0.000 Min. : 0.0 Min. : 0.00 Min. : 0.00
## 1st Qu.: 1.000 1st Qu.: 99.0 1st Qu.: 64.00 1st Qu.:27.50
## Median : 3.000 Median :117.0 Median : 72.00 Median :32.00
## Mean : 3.772 Mean :121.6 Mean : 69.38 Mean :32.15
## 3rd Qu.: 6.000 3rd Qu.:143.0 3rd Qu.: 80.00 3rd Qu.:36.60
## Max. :17.000 Max. :200.0 Max. :122.00 Max. :80.60
## Age Outcome
## Min. :21.00 Min. :0.0000
## 1st Qu.:24.00 1st Qu.:0.0000
## Median :29.00 Median :0.0000
## Mean :33.46 Mean :0.3496
## 3rd Qu.:41.00 3rd Qu.:1.0000
## Max. :81.00 Max. :1.0000Interpretasi: Dataset terdiri dari 2.251 observasi. Rata-rata kadar glukosa 121,6 mg/dL, rata-rata BMI 32,15 kg/m^2, dan rata-rata usia 33,46 tahun. Proporsi pasien diabetes 34,96% (787 pasien) dan tidak diabetes 65,04% (1.464 pasien).
1.3.2 Distribusi Status Diabetes
##
## Tidak Diabetes Diabetes
## 1464 787##
## Tidak Diabetes Diabetes
## 65.04 34.96Interpretasi: Proporsi pasien tidak diabetes (65,04%) lebih besar dibandingkan pasien diabetes (34,96%). Ketidakseimbangan kelas ini perlu diperhatikan karena dapat memengaruhi kemampuan model dalam mendeteksi kasus diabetes secara akurat.
1.3.3 Visualisasi Data
1.3.3.1 Distribusi Status Diabetes
ggplot(diabetes, aes(x = diabetes_status, fill = diabetes_status)) +
geom_bar() +
geom_text(
stat = "count",
aes(label = paste0(after_stat(count), "\n(",
round(after_stat(count) / nrow(diabetes) * 100, 1), "%)")),
vjust = -0.3
) +
labs(title = "Distribusi Status Diabetes",
x = "Status Diabetes", y = "Jumlah") +
theme_minimal() +
theme(legend.position = "none")
Interpretasi: Dari total 2.251 pasien, sebanyak 1.464 pasien (65,04%) berstatus tidak diabetes dan 787 pasien (34,96%) berstatus diabetes.
1.3.3.2 Boxplot Glucose berdasarkan Status Diabetes
ggplot(diabetes, aes(x = diabetes_status, y = Glucose, fill = diabetes_status)) +
geom_boxplot(alpha = 0.8) +
labs(title = "Glucose berdasarkan Status Diabetes",
x = "Status Diabetes", y = "Glucose (mg/dL)") +
theme_minimal() +
theme(legend.position = "none")
Interpretasi: Pasien diabetes memiliki median kadar glukosa yang lebih tinggi, mengindikasikan bahwa kadar glukosa merupakan faktor pembeda yang kuat antara kedua kelompok.
1.3.3.3 Histogram Usia berdasarkan Status Diabetes
ggplot(diabetes, aes(x = Age, fill = diabetes_status)) +
geom_histogram(binwidth = 5, alpha = 0.7, position = "identity") +
labs(title = "Distribusi Umur berdasarkan Status Diabetes",
x = "Usia (Tahun)", y = "Frekuensi") +
theme_minimal()
Interpretasi: Mayoritas pasien berusia 20-35 tahun. Pasien diabetes cenderung memiliki distribusi usia yang lebih menyebar ke kelompok usia lebih tua, mengindikasikan bahwa risiko diabetes meningkat seiring bertambahnya usia.
1.4 Uji Asumsi
1.4.1 Uji Linearitas Logit (Box-Tidwell)
diabetes_bt <- diabetes %>%
mutate(
Glucose_log = Glucose * log(Glucose + 1),
BMI_log = BMI * log(BMI + 1),
Age_log = Age * log(Age + 1)
)
fit_bt <- glm(
Outcome ~ Pregnancies + Glucose + Glucose_log + BloodPressure +
SkinThickness + Insulin + BMI + BMI_log +
DiabetesPedigreeFunction + Age + Age_log,
data = diabetes_bt, family = binomial
)
box_tidwell <- broom::tidy(fit_bt) %>%
dplyr::filter(grepl("log", term)) %>%
dplyr::select(term, std.error, statistic, p.value) %>%
dplyr::mutate(across(where(is.numeric), \(x) round(x, 4)))
knitr::kable(box_tidwell, caption = "Hasil Uji Box-Tidwell")| term | std.error | statistic | p.value |
|---|---|---|---|
| Glucose_log | 0.0053 | 6.0706 | 0.0000 |
| BMI_log | 0.0470 | -1.3747 | 0.1692 |
| Age_log | 0.0329 | -8.7730 | 0.0000 |
Interpretasi: Variabel Glucose (statistik = 6,07; p < 0,001) dan Age (statistik = -8,77; p < 0,001) tidak memenuhi asumsi linearitas logit. Transformasi logaritma diterapkan pada kedua variabel tersebut (menjadi
log_Glucosedanlog_Age). Variabel BMI (p = 0,169) telah memenuhi asumsi.
1.4.2 Uji Multikolinearitas (VIF)
fit_vif <- lm(
Outcome ~ Pregnancies + log_Glucose + BloodPressure + SkinThickness +
Insulin + BMI + DiabetesPedigreeFunction + log_Age,
data = diabetes
)
vif_result <- vif(fit_vif)
vif_df <- data.frame(Variabel = names(vif_result), VIF = round(vif_result, 3))
knitr::kable(vif_df, caption = "Hasil Uji VIF")| Variabel | VIF | |
|---|---|---|
| Pregnancies | Pregnancies | 1.578 |
| log_Glucose | log_Glucose | 1.129 |
| BloodPressure | BloodPressure | 1.202 |
| SkinThickness | SkinThickness | 1.516 |
| Insulin | Insulin | 1.402 |
| BMI | BMI | 1.289 |
| DiabetesPedigreeFunction | DiabetesPedigreeFunction | 1.067 |
| log_Age | log_Age | 1.720 |
Interpretasi: Seluruh variabel memiliki nilai VIF di bawah 10 (tertinggi log_Age = 1,720; terendah DiabetesPedigreeFunction = 1,067). Tidak terdapat masalah multikolinearitas dalam model.
1.5 Pembagian Data
set.seed(123)
idx <- createDataPartition(diabetes$Outcome, p = 0.8, list = FALSE)
train_data <- diabetes[idx, ]
test_data <- diabetes[-idx, ]
cat("Jumlah training:", nrow(train_data), "\n")## Jumlah training: 1801## Jumlah testing : 450Interpretasi: Data dibagi 80:20 menggunakan stratified random sampling. Sebanyak 1.801 observasi untuk training dan 450 untuk testing.
set.seed(123)memastikan hasil reproducible.
1.6 Pembuatan Model
1.6.1 Model Penuh
model_full <- glm(
Outcome ~ Pregnancies + log_Glucose + BloodPressure + SkinThickness +
Insulin + BMI + DiabetesPedigreeFunction + log_Age,
data = train_data, family = binomial(link = "logit")
)
summary(model_full)##
## Call:
## glm(formula = Outcome ~ Pregnancies + log_Glucose + BloodPressure +
## SkinThickness + Insulin + BMI + DiabetesPedigreeFunction +
## log_Age, family = binomial(link = "logit"), data = train_data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.128e+01 1.391e+00 -15.302 < 2e-16 ***
## Pregnancies 1.107e-01 2.174e-02 5.091 3.57e-07 ***
## log_Glucose 2.832e+00 2.775e-01 10.206 < 2e-16 ***
## BloodPressure -9.535e-03 3.613e-03 -2.639 0.00831 **
## SkinThickness 9.225e-04 4.424e-03 0.209 0.83482
## Insulin 3.452e-05 6.036e-04 0.057 0.95440
## BMI 8.428e-02 9.531e-03 8.843 < 2e-16 ***
## DiabetesPedigreeFunction 8.317e-01 1.907e-01 4.362 1.29e-05 ***
## log_Age 1.156e+00 2.379e-01 4.861 1.17e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2334.1 on 1800 degrees of freedom
## Residual deviance: 1796.4 on 1792 degrees of freedom
## AIC: 1814.4
##
## Number of Fisher Scoring iterations: 5Interpretasi: Variabel Pregnancies, log_Glucose, BloodPressure, BMI, DiabetesPedigreeFunction, dan log_Age menunjukkan pengaruh signifikan (p < 0,05). Variabel SkinThickness (p = 0,835) dan Insulin (p = 0,954) tidak signifikan. AIC model penuh = 1.814,4.
1.6.2 Seleksi Variabel Stepwise
##
## Call:
## glm(formula = Outcome ~ Pregnancies + log_Glucose + BloodPressure +
## BMI + DiabetesPedigreeFunction + log_Age, family = binomial(link = "logit"),
## data = train_data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.293109 1.349924 -15.774 < 2e-16 ***
## Pregnancies 0.110517 0.021714 5.090 3.59e-07 ***
## log_Glucose 2.838463 0.260472 10.897 < 2e-16 ***
## BloodPressure -0.009414 0.003577 -2.632 0.0085 **
## BMI 0.085055 0.009044 9.405 < 2e-16 ***
## DiabetesPedigreeFunction 0.838875 0.188155 4.458 8.26e-06 ***
## log_Age 1.148060 0.234161 4.903 9.44e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2334.1 on 1800 degrees of freedom
## Residual deviance: 1796.4 on 1794 degrees of freedom
## AIC: 1810.4
##
## Number of Fisher Scoring iterations: 5Interpretasi: Metode stepwise dua arah mengeluarkan SkinThickness dan Insulin. Model final terdiri dari 6 prediktor: Pregnancies, log_Glucose, BloodPressure, BMI, DiabetesPedigreeFunction, dan log_Age. AIC model final = 1.810,4 – lebih efisien dari model penuh.
1.7 Uji Signifikansi Model
1.7.1 Uji G (Likelihood Ratio Test)
G2 <- model_final$null.deviance - model_final$deviance
df_G2 <- model_final$df.null - model_final$df.residual
p_G2 <- pchisq(G2, df = df_G2, lower.tail = FALSE)
uji_G <- data.frame(G2 = round(G2, 3), df = df_G2, p_value = signif(p_G2, 4))
knitr::kable(uji_G, caption = "Hasil Uji G")| G2 | df | p_value |
|---|---|---|
| 537.691 | 6 | 0 |
Interpretasi: G^2 = 537,691; df = 6; p-value < 0,001. Model secara simultan signifikan – minimal satu prediktor berpengaruh nyata terhadap probabilitas diabetes.
1.7.2 Uji Wald
wald_final <- broom::tidy(model_final) %>%
filter(term != "(Intercept)") %>%
mutate(
W2 = (estimate / std.error)^2,
Keterangan = ifelse(p.value < 0.05, "Signifikan", "Tidak Signifikan")
) %>%
transmute(
Variabel = term,
W2 = round(W2, 3),
p_value = round(p.value, 4),
Keterangan = Keterangan
)
knitr::kable(wald_final, caption = "Hasil Uji Wald")| Variabel | W2 | p_value | Keterangan |
|---|---|---|---|
| Pregnancies | 25.904 | 0.0000 | Signifikan |
| log_Glucose | 118.753 | 0.0000 | Signifikan |
| BloodPressure | 6.926 | 0.0085 | Signifikan |
| BMI | 88.451 | 0.0000 | Signifikan |
| DiabetesPedigreeFunction | 19.878 | 0.0000 | Signifikan |
| log_Age | 24.038 | 0.0000 | Signifikan |
Interpretasi: Seluruh variabel dalam model final signifikan secara parsial (p < 0,05). log_Glucose memiliki W^2 tertinggi (118,753), diikuti BMI (88,451) – keduanya prediktor paling dominan.
1.8 Odds Ratio
or_final <- broom::tidy(model_final) %>%
mutate(
OR = exp(estimate),
CI_Low = exp(estimate - 1.96 * std.error),
CI_High = exp(estimate + 1.96 * std.error)
) %>%
transmute(
Variabel = term,
Odds_Ratio = round(OR, 3),
CI_95 = paste0(round(CI_Low, 3), " - ", round(CI_High, 3)),
p_value = round(p.value, 4)
)
knitr::kable(or_final, caption = "Odds Ratio Model Final")| Variabel | Odds_Ratio | CI_95 | p_value |
|---|---|---|---|
| (Intercept) | 0.000 | 0 - 0 | 0.0000 |
| Pregnancies | 1.117 | 1.07 - 1.165 | 0.0000 |
| log_Glucose | 17.089 | 10.257 - 28.474 | 0.0000 |
| BloodPressure | 0.991 | 0.984 - 0.998 | 0.0085 |
| BMI | 1.089 | 1.07 - 1.108 | 0.0000 |
| DiabetesPedigreeFunction | 2.314 | 1.6 - 3.346 | 0.0000 |
| log_Age | 3.152 | 1.992 - 4.988 | 0.0000 |
Interpretasi: log_Glucose memiliki OR tertinggi (17,089; 95% CI: 10,257-28,474) – setiap kenaikan satu satuan log glukosa meningkatkan peluang diabetes sebesar 17,089 kali. log_Age (OR = 3,152), DiabetesPedigreeFunction (OR = 2,314), Pregnancies (OR = 1,117), dan BMI (OR = 1,089) juga berpengaruh positif. BloodPressure (OR = 0,991) sedikit menurunkan peluang diabetes setelah mengontrol variabel lain.
1.9 Uji Hosmer-Lemeshow
##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: train_data$Outcome, fitted(model_final)
## X-squared = 35.923, df = 8, p-value = 1.814e-05Interpretasi: X^2 = 35,923; df = 8; p-value < 0,05. Secara statistik terdapat perbedaan antara nilai diamati dan diprediksi. Hal ini dapat terjadi karena sampel besar membuat uji sangat sensitif. Evaluasi melalui AUC dan confusion matrix pada data testing tetap menunjukkan hasil yang baik.
1.10 Evaluasi Model
1.10.1 Confusion Matrix
prob_test <- predict(model_final, newdata = test_data, type = "response")
pred_test <- ifelse(prob_test >= 0.5, 1, 0)
cm <- table(Aktual = test_data$Outcome, Prediksi = pred_test)
cm## Prediksi
## Aktual 0 1
## 0 250 45
## 1 60 95Interpretasi: Model memprediksi dengan benar 250 kasus tidak diabetes (TN) dan 95 kasus diabetes (TP). Terdapat 45 False Positive dan 60 False Negative. Kasus False Negative perlu perhatian khusus mengingat implikasi klinis yang lebih serius.
1.10.2 Metrik Evaluasi
TP <- cm[2,2]; TN <- cm[1,1]; FP <- cm[1,2]; FN <- cm[2,1]
accuracy <- (TP + TN) / sum(cm)
sensitivity <- TP / (TP + FN)
specificity <- TN / (TN + FP)
precision <- TP / (TP + FP)
f1_score <- 2 * precision * sensitivity / (precision + sensitivity)
hasil_metric <- data.frame(
Metrik = c("Accuracy", "Sensitivity", "Specificity", "Precision", "F1-Score"),
Nilai = round(c(accuracy, sensitivity, specificity, precision, f1_score), 4)
)
knitr::kable(hasil_metric, caption = "Metrik Evaluasi")| Metrik | Nilai |
|---|---|
| Accuracy | 0.7667 |
| Sensitivity | 0.6129 |
| Specificity | 0.8475 |
| Precision | 0.6786 |
| F1-Score | 0.6441 |
Interpretasi: Akurasi 76,67%; spesifisitas 84,75%; sensitivitas 61,29%; F1-Score 64,41%. Model baik dalam mendeteksi pasien tidak diabetes, namun masih perlu peningkatan sensitivitas untuk kasus diabetes.
1.11 Kurva ROC dan AUC
roc_obj <- roc(test_data$Outcome, prob_test)
plot(roc_obj, col = "blue", lwd = 2,
main = paste("ROC Curve | AUC =", round(auc(roc_obj), 4)))
abline(a = 0, b = 1, lty = 2, col = "gray")
## AUC: 0.8325Interpretasi: AUC = 0,8325 (kategori baik, > 0,8). Model memiliki probabilitas 83,25% untuk memberikan skor prediksi lebih tinggi pada pasien diabetes dibandingkan pasien tidak diabetes yang dipilih secara acak.
1.12 Pseudo R-Squared
## fitting null model for pseudo-r2## llh llhNull G2 McFadden r2ML
## -898.2238537 -1167.0694906 537.6912738 0.2303596 0.2581079
## r2CU
## 0.3553334Interpretasi: McFadden R^2 = 0,2304 (model menjelaskan ~23,04% variasi status diabetes). Nagelkerke R^2 = 0,3553 (~35,53%). Nilai McFadden antara 0,20-0,40 sudah dianggap memadai untuk regresi logistik.
1.13 Kesimpulan Bagian I
- Regresi logistik biner berhasil digunakan untuk memodelkan status diabetes pada 2.251 observasi.
- Tidak terdapat masalah multikolinearitas (VIF tertinggi 1,720).
- Variabel Glucose dan Age ditransformasi logaritma karena tidak memenuhi asumsi linearitas logit.
- Model final (6 prediktor) lebih parsimonious dibanding model penuh, dengan log_Glucose dan BMI sebagai prediktor paling dominan.
- Model secara simultan dan parsial signifikan (Uji G: G^2 = 537,691; p < 0,001).
- Evaluasi testing: akurasi 76,67%, spesifisitas 84,75%, sensitivitas 61,29%, AUC = 0,8325.
2 Regresi Logistik Multinomial
2.1 Pendahuluan & Deskripsi Dataset
Regresi logistik multinomial merupakan perluasan dari regresi logistik biner yang digunakan ketika variabel respons bersifat nominal dengan lebih dari dua kategori tanpa urutan alami. Metode ini memodelkan log-odds setiap kategori relatif terhadap kategori referensi sebagai fungsi linear dari variabel prediktor.
Pada analisis ini digunakan Dhaka Obesity: Dataset based on Eating Behavior and Physical Conditions untuk memodelkan moda transportasi berdasarkan karakteristik demografis, kebiasaan makan, dan kondisi fisik individu.
Sumber data: Dhaka Obesity Dataset – Mendeley Data:
https://data.mendeley.com/datasets/vhnks4n4h7/2
2.1.1 Variabel Penelitian
| Variabel | Tipe | Keterangan | Kategori |
|---|---|---|---|
| Transport.Mode | Respons | Moda transportasi (digabung) | Active, Public, Private |
| Gender | Kategorik | Jenis kelamin | Male, Female |
| Family.Hist.OW | Kategorik | Riwayat keluarga overweight | Yes, No |
| High.Caloric.Food | Kategorik | Konsumsi makanan tinggi kalori | Yes, No |
| Vegetable.Freq | Kategorik | Frekuensi konsumsi sayur | Always, Sometimes, Never |
| Smoking | Kategorik | Kebiasaan merokok | Yes, No |
| Alcohol.Intake | Kategorik | Konsumsi alkohol | Always, Frequently, Sometimes, I do not drink |
| Physical.Exercise | Kategorik | Frekuensi olahraga | I do not have, 1 or 2 days, 2 or 4 days, 4 or 5 days, Almost Everyday |
| Device.Usage | Kategorik | Durasi penggunaan perangkat/hari | 0-2 hours, 3-5 hours, More than 5 hours |
| Obesity.Level | Kategorik | Tingkat obesitas | Insufficient_Weight, Normal_Weight, Overweight, Obesity_Type_I/II/III |
2.2 Pemuatan Data
# Ganti path sesuai lokasi file Anda
dhaka_path <- "C:/Users/hp/Downloads/Dhaka Obesity.csv"
if (!file.exists(dhaka_path)) stop("File 'Dhaka Obesity.csv' tidak ditemukan di: ", dhaka_path)
data_raw <- read.csv(dhaka_path)
cat("Dimensi data mentah:", nrow(data_raw), "baris x", ncol(data_raw), "kolom\n")## Dimensi data mentah: 2182 baris x 17 kolom## Rows: 2,182
## Columns: 17
## $ Gender <chr> "Male", "Female", "Male", "Mal…
## $ Age <int> 29, 25, 23, 22, 19, 24, 30, 16…
## $ Height..m. <dbl> 1.65, 1.65, 1.70, 1.68, 1.75, …
## $ Weight..kg. <dbl> 101, 53, 70, 112, 67, 70, 90, …
## $ Family.history.of.overweight <chr> "Yes", "No", "No", "Yes", "No"…
## $ High.caloric.food.consumption <chr> "Yes", "No", "No", "Yes", "Yes…
## $ Vegetable.consumption.frequency <chr> "Sometimes", "Always", "Always…
## $ Daily.main.meals.frequency <chr> "Three", "Three", "Between 1-2…
## $ Between.meal.food.consumption.frequency <chr> "Frequently", "Always", "Somet…
## $ Smoking <chr> "No", "No", "Yes", "No", "Yes"…
## $ Alcohol.intake <chr> "I do not drink", "I do not dr…
## $ Daily.water.intake <chr> "Between 1 and 2 L", "Between …
## $ Monitor.calories <chr> "No", "No", "No", "No", "No", …
## $ Physical.exercise <chr> "I do not have", "Almost Every…
## $ Daily.device.usage.duration <chr> "More than 5 hours", "More tha…
## $ Mode.of.transportation <chr> "Private Car", "Public Transpo…
## $ Obesity.level <chr> "Obesity_Type_II", "Normal_Wei…2.3 Preprocessing
data_multi <- data_raw %>%
rename(
Transport.Mode = `Mode.of.transportation`,
Family.Hist.OW = `Family.history.of.overweight`,
High.Caloric.Food = `High.caloric.food.consumption`,
Vegetable.Freq = `Vegetable.consumption.frequency`,
Alcohol.Intake = `Alcohol.intake`,
Physical.Exercise = `Physical.exercise`,
Device.Usage = `Daily.device.usage.duration`,
Obesity.Level = `Obesity.level`
) %>%
mutate(
Transport.Mode = case_when(
Transport.Mode %in% c("Walking", "Bike") ~ "Active",
Transport.Mode %in% c("Public Transportation", "Rickshaw") ~ "Public",
Transport.Mode == "Private Car" ~ "Private",
TRUE ~ NA_character_
)
) %>%
dplyr::select(
Transport.Mode, Gender, Family.Hist.OW, High.Caloric.Food,
Vegetable.Freq, Smoking, Alcohol.Intake, Physical.Exercise,
Device.Usage, Obesity.Level
) %>%
drop_na() %>%
mutate(
Transport.Mode = factor(Transport.Mode, levels = c("Private", "Active", "Public")),
Gender = factor(Gender),
Family.Hist.OW = factor(Family.Hist.OW),
High.Caloric.Food = factor(High.Caloric.Food),
Vegetable.Freq = factor(Vegetable.Freq),
Smoking = factor(Smoking),
Alcohol.Intake = factor(Alcohol.Intake),
Physical.Exercise = factor(Physical.Exercise),
Device.Usage = factor(Device.Usage),
Obesity.Level = factor(Obesity.Level)
)
cat("Dimensi data setelah preprocessing:",
nrow(data_multi), "baris x", ncol(data_multi), "kolom\n")## Dimensi data setelah preprocessing: 2182 baris x 10 kolom## Rows: 2,182
## Columns: 10
## $ Transport.Mode <fct> Private, Public, Public, Public, Active, Public, Pri…
## $ Gender <fct> Male, Female, Male, Male, Male, Male, Male, Male, Ma…
## $ Family.Hist.OW <fct> Yes, No, No, Yes, No, Yes, No, Yes, Yes, No, No, Yes…
## $ High.Caloric.Food <fct> Yes, No, No, Yes, Yes, No, Yes, Yes, Yes, No, No, Ye…
## $ Vegetable.Freq <fct> Sometimes, Always, Always, Sometimes, Always, Always…
## $ Smoking <fct> No, No, Yes, No, Yes, No, No, No, No, No, No, No, No…
## $ Alcohol.Intake <fct> I do not drink, I do not drink, Sometimes, I do not …
## $ Physical.Exercise <fct> I do not have, Almost Everyday, I do not have, I do …
## $ Device.Usage <fct> More than 5 hours, More than 5 hours, More than 5 ho…
## $ Obesity.Level <fct> Obesity_Type_II, Normal_Weight, Normal_Weight, Obesi…Interpretasi: Variabel moda transportasi yang semula 5 kategori digabungkan menjadi 3 kategori nominal: Active (Walking + Bike), Public (Public Transportation + Rickshaw), dan Private (Private Car). Kategori Private ditetapkan sebagai referensi.
2.4 Eksplorasi Data (EDA)
2.4.1 Statistik Deskriptif
## Transport.Mode Gender Family.Hist.OW High.Caloric.Food Vegetable.Freq
## Private: 493 Female: 917 No :1354 No : 959 Always : 548
## Active : 372 Male :1265 Yes: 828 Yes:1223 Never : 179
## Public :1317 Sometimes:1455
##
##
##
## Smoking Alcohol.Intake Physical.Exercise
## No :1898 Always : 8 1 or 2 days : 292
## Yes: 284 Frequently : 19 2 or 4 days : 120
## I do not drink:2059 4 or 5 days : 52
## Sometimes : 96 Almost Everyday: 95
## I do not have :1623
##
## Device.Usage Obesity.Level
## 0-2 hours : 127 Insufficient_Weight:378
## 3-5 hours : 776 Normal_Weight :530
## More than 5 hours:1279 Obesity_Type_I :317
## Obesity_Type_II :283
## Obesity_Type_III :251
## Overweight :4232.4.2 Distribusi Variabel Respons
dist_resp <- data_multi %>%
count(Transport.Mode) %>%
mutate(Proporsi = n / sum(n))
dist_resp %>%
kable(
digits = 3,
caption = "Distribusi Moda Transportasi",
col.names = c("Moda Transportasi", "Frekuensi", "Proporsi")
) %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE
)| Moda Transportasi | Frekuensi | Proporsi |
|---|---|---|
| Private | 493 | 0.226 |
| Active | 372 | 0.170 |
| Public | 1317 | 0.604 |
Interpretasi: Tabel menampilkan frekuensi dan proporsi setiap kategori moda transportasi. Ketidakseimbangan kelas dapat memengaruhi performa model.
2.4.3 Visualisasi Data
2.4.3.1 Distribusi Moda Transportasi
data_multi %>%
count(Transport.Mode) %>%
mutate(Proporsi = n / sum(n)) %>%
ggplot(aes(x = Transport.Mode, y = Proporsi, fill = Transport.Mode)) +
geom_col(width = 0.6, alpha = 0.92) +
geom_text(aes(label = percent(Proporsi, accuracy = 0.1)),
vjust = -0.5, fontface = "bold", color = "#0c2340") +
scale_y_continuous(labels = percent_format(), limits = c(0, 0.7)) +
scale_fill_manual(values = c("#0077b6", "#5568b8", "#90e0ef")) +
labs(title = "Distribusi Moda Transportasi",
subtitle = "Variabel respons multinomial: Active, Public, Private",
x = "Moda Transportasi", y = "Proporsi") +
theme_minimal(base_size = 13) +
theme(legend.position = "none")
Interpretasi: Grafik batang menunjukkan distribusi proporsi ketiga kategori moda transportasi. Perbedaan proporsi memberikan gambaran awal komposisi data sebelum pemodelan.
2.4.3.2 Moda Transportasi berdasarkan Jenis Kelamin
data_multi %>%
count(Gender, Transport.Mode) %>%
group_by(Gender) %>%
mutate(Proporsi = n / sum(n)) %>%
ggplot(aes(x = Transport.Mode, y = Proporsi, fill = Transport.Mode)) +
geom_col(width = 0.6, alpha = 0.92) +
geom_text(aes(label = percent(Proporsi, accuracy = 0.1)),
vjust = -0.4, size = 3.5, color = "#0c2340") +
facet_wrap(~Gender) +
scale_y_continuous(labels = percent_format(), limits = c(0, 0.75)) +
scale_fill_manual(values = c("#0077b6", "#5568b8", "#90e0ef")) +
labs(title = "Moda Transportasi berdasarkan Jenis Kelamin",
subtitle = "Perbandingan proporsi antar kategori pada tiap kelompok gender",
x = "Moda Transportasi", y = "Proporsi") +
theme_minimal(base_size = 13) +
theme(legend.position = "none")
Interpretasi: Perbedaan pola distribusi antar gender mengindikasikan apakah gender berpotensi menjadi prediktor yang signifikan.
2.4.3.3 Moda Transportasi berdasarkan Tingkat Obesitas
data_multi %>%
count(Obesity.Level, Transport.Mode) %>%
group_by(Obesity.Level) %>%
mutate(Proporsi = n / sum(n)) %>%
ggplot(aes(x = Transport.Mode, y = Proporsi, fill = Transport.Mode)) +
geom_col(width = 0.6, alpha = 0.92) +
geom_text(aes(label = percent(Proporsi, accuracy = 0.1)),
vjust = -0.4, size = 3, color = "#0c2340") +
facet_wrap(~Obesity.Level, ncol = 3) +
scale_y_continuous(labels = percent_format(), limits = c(0, 0.85)) +
scale_fill_manual(values = c("#0077b6", "#5568b8", "#90e0ef")) +
labs(title = "Moda Transportasi berdasarkan Tingkat Obesitas",
subtitle = "Perbandingan proporsi di setiap kelompok tingkat obesitas",
x = "Moda Transportasi", y = "Proporsi") +
theme_minimal(base_size = 12) +
theme(legend.position = "none",
axis.text.x = element_text(angle = 15, hjust = 1))
Interpretasi: Jika terdapat perbedaan pola proporsi yang nyata antar kelompok obesitas, hal ini mengindikasikan bahwa tingkat obesitas berhubungan dengan pilihan moda transportasi.
2.4.3.4 Moda Transportasi berdasarkan Frekuensi Olahraga
data_multi %>%
count(Physical.Exercise, Transport.Mode) %>%
group_by(Physical.Exercise) %>%
mutate(Proporsi = n / sum(n)) %>%
ggplot(aes(x = Transport.Mode, y = Proporsi, fill = Transport.Mode)) +
geom_col(width = 0.6, alpha = 0.92) +
geom_text(aes(label = percent(Proporsi, accuracy = 0.1)),
vjust = -0.4, size = 3, color = "#0c2340") +
facet_wrap(~Physical.Exercise, ncol = 3) +
scale_y_continuous(labels = percent_format(), limits = c(0, 0.85)) +
scale_fill_manual(values = c("#0077b6", "#5568b8", "#90e0ef")) +
labs(title = "Moda Transportasi berdasarkan Frekuensi Olahraga",
subtitle = "Perbandingan proporsi di setiap kelompok frekuensi olahraga",
x = "Moda Transportasi", y = "Proporsi") +
theme_minimal(base_size = 12) +
theme(legend.position = "none",
axis.text.x = element_text(angle = 15, hjust = 1))
Interpretasi: Individu yang lebih aktif berolahraga mungkin cenderung memilih moda transportasi aktif seperti berjalan atau bersepeda.
2.5 Estimasi Model Regresi Logistik Multinomial
2.5.1 Pembuatan Model
Model diestimasi menggunakan multinom() dari paket
nnet. Kategori referensi: Private
(kendaraan pribadi).
fit_multi <- nnet::multinom(
Transport.Mode ~
Gender + Family.Hist.OW + High.Caloric.Food +
Vegetable.Freq + Smoking + Alcohol.Intake +
Physical.Exercise + Device.Usage + Obesity.Level,
data = data_multi,
trace = FALSE
)
summary(fit_multi)## Call:
## nnet::multinom(formula = Transport.Mode ~ Gender + Family.Hist.OW +
## High.Caloric.Food + Vegetable.Freq + Smoking + Alcohol.Intake +
## Physical.Exercise + Device.Usage + Obesity.Level, data = data_multi,
## trace = FALSE)
##
## Coefficients:
## (Intercept) GenderMale Family.Hist.OWYes High.Caloric.FoodYes
## Active 2.849214 1.3298548 -0.5600626 0.1195473
## Public 2.578307 0.7939367 -0.2835522 0.0764246
## Vegetable.FreqNever Vegetable.FreqSometimes SmokingYes
## Active -0.1917066 -0.4061606 0.8958892
## Public -0.3265678 -0.2318609 0.5148171
## Alcohol.IntakeFrequently Alcohol.IntakeI do not drink
## Active -0.5800169 0.8300789
## Public -1.0593413 0.9392414
## Alcohol.IntakeSometimes Physical.Exercise2 or 4 days
## Active 0.1146742 -0.7054970
## Public -0.1040944 -0.4304671
## Physical.Exercise4 or 5 days Physical.ExerciseAlmost Everyday
## Active 0.25483204 0.04186118
## Public 0.06579829 0.12111409
## Physical.ExerciseI do not have Device.Usage3-5 hours
## Active -0.3302981 -1.2066934
## Public 0.1490701 -0.3568612
## Device.UsageMore than 5 hours Obesity.LevelNormal_Weight
## Active -1.2905785 -2.226003
## Public -0.2829183 -1.814427
## Obesity.LevelObesity_Type_I Obesity.LevelObesity_Type_II
## Active -2.732621 -6.024309
## Public -2.403442 -3.326246
## Obesity.LevelObesity_Type_III Obesity.LevelOverweight
## Active -18.385402 -2.974187
## Public -3.187067 -2.634862
##
## Std. Errors:
## (Intercept) GenderMale Family.Hist.OWYes High.Caloric.FoodYes
## Active 1.237569 0.1726598 0.1780637 0.1887925
## Public 1.054628 0.1228972 0.1204943 0.1369889
## Vegetable.FreqNever Vegetable.FreqSometimes SmokingYes
## Active 0.4221799 0.1901581 0.2672792
## Public 0.2510444 0.1582811 0.2311155
## Alcohol.IntakeFrequently Alcohol.IntakeI do not drink
## Active 1.255766 1.1152051
## Public 1.087594 0.9319366
## Alcohol.IntakeSometimes Physical.Exercise2 or 4 days
## Active 1.1487063 0.3528724
## Public 0.9623817 0.2895111
## Physical.Exercise4 or 5 days Physical.ExerciseAlmost Everyday
## Active 0.5025385 0.4260159
## Public 0.4525628 0.3835619
## Physical.ExerciseI do not have Device.Usage3-5 hours
## Active 0.2258542 0.3357033
## Public 0.1906409 0.2986977
## Device.UsageMore than 5 hours Obesity.LevelNormal_Weight
## Active 0.3276401 0.4209734
## Public 0.2933518 0.3960059
## Obesity.LevelObesity_Type_I Obesity.LevelObesity_Type_II
## Active 0.4466642 0.6744271
## Public 0.4060578 0.4084504
## Obesity.LevelObesity_Type_III Obesity.LevelOverweight
## Active 8.187970e-06 0.4291956
## Public 4.099004e-01 0.3962747
##
## Residual Deviance: 3501.439
## AIC: 3585.439Interpretasi: Output menampilkan dua set koefisien: Active vs Private dan Public vs Private. Koefisien positif meningkatkan peluang berada di kategori tersebut dibandingkan Private, sebaliknya untuk negatif.
2.5.2 Tabel Koefisien, SE, z-value, dan p-value
Karena multinom() tidak menghasilkan p-value secara
langsung, nilai z dan p-value dihitung dari:
\[z = \frac{\hat{\beta}}{SE(\hat{\beta})}, \quad p\text{-value} = 2\{1 - \Phi(|z|)\}\]
multi_sum <- summary(fit_multi)
coef_df <- as.data.frame(multi_sum$coefficients)
se_df <- as.data.frame(multi_sum$standard.errors)
coef_long <- coef_df %>%
rownames_to_column("Kategori") %>%
pivot_longer(-Kategori, names_to = "Variabel", values_to = "Estimate")
se_long <- se_df %>%
rownames_to_column("Kategori") %>%
pivot_longer(-Kategori, names_to = "Variabel", values_to = "SE")
result_multi <- coef_long %>%
left_join(se_long, by = c("Kategori", "Variabel")) %>%
mutate(
z_value = Estimate / SE,
p_value = 2 * (1 - pnorm(abs(z_value))),
RRR = exp(Estimate),
CI_low = exp(Estimate - 1.96 * SE),
CI_high = exp(Estimate + 1.96 * SE),
Signifikan = ifelse(p_value < 0.05, "Ya", "Tidak")
)
result_multi %>%
mutate(across(c(Estimate, SE, z_value, p_value, RRR, CI_low, CI_high),
~ round(.x, 4))) %>%
kable(
caption = "Ringkasan Hasil Regresi Logistik Multinomial",
col.names = c("Kategori", "Variabel", "Estimate", "SE", "z-value",
"p-value", "RRR", "CI 2.5%", "CI 97.5%", "Signifikan")
) %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE
) %>%
column_spec(10, bold = TRUE,
color = ifelse(result_multi$Signifikan == "Ya", "#0077b6", "#b84f5a"))| Kategori | Variabel | Estimate | SE | z-value | p-value | RRR | CI 2.5% | CI 97.5% | Signifikan |
|---|---|---|---|---|---|---|---|---|---|
| Active | (Intercept) | 2.8492 | 1.2376 | 2.3023 | 0.0213 | 17.2742 | 1.5274 | 195.3608 | Ya |
| Active | GenderMale | 1.3299 | 0.1727 | 7.7022 | 0.0000 | 3.7805 | 2.6951 | 5.3030 | Ya |
| Active | Family.Hist.OWYes | -0.5601 | 0.1781 | -3.1453 | 0.0017 | 0.5712 | 0.4029 | 0.8097 | Ya |
| Active | High.Caloric.FoodYes | 0.1195 | 0.1888 | 0.6332 | 0.5266 | 1.1270 | 0.7784 | 1.6316 | Tidak |
| Active | Vegetable.FreqNever | -0.1917 | 0.4222 | -0.4541 | 0.6498 | 0.8255 | 0.3609 | 1.8885 | Tidak |
| Active | Vegetable.FreqSometimes | -0.4062 | 0.1902 | -2.1359 | 0.0327 | 0.6662 | 0.4589 | 0.9671 | Ya |
| Active | SmokingYes | 0.8959 | 0.2673 | 3.3519 | 0.0008 | 2.4495 | 1.4507 | 4.1361 | Ya |
| Active | Alcohol.IntakeFrequently | -0.5800 | 1.2558 | -0.4619 | 0.6442 | 0.5599 | 0.0478 | 6.5619 | Tidak |
| Active | Alcohol.IntakeI do not drink | 0.8301 | 1.1152 | 0.7443 | 0.4567 | 2.2935 | 0.2578 | 20.4071 | Tidak |
| Active | Alcohol.IntakeSometimes | 0.1147 | 1.1487 | 0.0998 | 0.9205 | 1.1215 | 0.1180 | 10.6562 | Tidak |
| Active | Physical.Exercise2 or 4 days | -0.7055 | 0.3529 | -1.9993 | 0.0456 | 0.4939 | 0.2473 | 0.9862 | Ya |
| Active | Physical.Exercise4 or 5 days | 0.2548 | 0.5025 | 0.5071 | 0.6121 | 1.2902 | 0.4818 | 3.4549 | Tidak |
| Active | Physical.ExerciseAlmost Everyday | 0.0419 | 0.4260 | 0.0983 | 0.9217 | 1.0427 | 0.4524 | 2.4033 | Tidak |
| Active | Physical.ExerciseI do not have | -0.3303 | 0.2259 | -1.4624 | 0.1436 | 0.7187 | 0.4616 | 1.1189 | Tidak |
| Active | Device.Usage3-5 hours | -1.2067 | 0.3357 | -3.5945 | 0.0003 | 0.2992 | 0.1549 | 0.5777 | Ya |
| Active | Device.UsageMore than 5 hours | -1.2906 | 0.3276 | -3.9390 | 0.0001 | 0.2751 | 0.1447 | 0.5229 | Ya |
| Active | Obesity.LevelNormal_Weight | -2.2260 | 0.4210 | -5.2878 | 0.0000 | 0.1080 | 0.0473 | 0.2464 | Ya |
| Active | Obesity.LevelObesity_Type_I | -2.7326 | 0.4467 | -6.1178 | 0.0000 | 0.0650 | 0.0271 | 0.1561 | Ya |
| Active | Obesity.LevelObesity_Type_II | -6.0243 | 0.6744 | -8.9325 | 0.0000 | 0.0024 | 0.0006 | 0.0091 | Ya |
| Active | Obesity.LevelObesity_Type_III | -18.3854 | 0.0000 | -2245416.3254 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | Ya |
| Active | Obesity.LevelOverweight | -2.9742 | 0.4292 | -6.9297 | 0.0000 | 0.0511 | 0.0220 | 0.1185 | Ya |
| Public | (Intercept) | 2.5783 | 1.0546 | 2.4448 | 0.0145 | 13.1748 | 1.6674 | 104.1026 | Ya |
| Public | GenderMale | 0.7939 | 0.1229 | 6.4602 | 0.0000 | 2.2121 | 1.7386 | 2.8146 | Ya |
| Public | Family.Hist.OWYes | -0.2836 | 0.1205 | -2.3532 | 0.0186 | 0.7531 | 0.5947 | 0.9537 | Ya |
| Public | High.Caloric.FoodYes | 0.0764 | 0.1370 | 0.5579 | 0.5769 | 1.0794 | 0.8252 | 1.4119 | Tidak |
| Public | Vegetable.FreqNever | -0.3266 | 0.2510 | -1.3008 | 0.1933 | 0.7214 | 0.4410 | 1.1800 | Tidak |
| Public | Vegetable.FreqSometimes | -0.2319 | 0.1583 | -1.4649 | 0.1430 | 0.7931 | 0.5815 | 1.0815 | Tidak |
| Public | SmokingYes | 0.5148 | 0.2311 | 2.2275 | 0.0259 | 1.6733 | 1.0638 | 2.6322 | Ya |
| Public | Alcohol.IntakeFrequently | -1.0593 | 1.0876 | -0.9740 | 0.3300 | 0.3467 | 0.0411 | 2.9222 | Tidak |
| Public | Alcohol.IntakeI do not drink | 0.9392 | 0.9319 | 1.0078 | 0.3135 | 2.5580 | 0.4117 | 15.8923 | Tidak |
| Public | Alcohol.IntakeSometimes | -0.1041 | 0.9624 | -0.1082 | 0.9139 | 0.9011 | 0.1366 | 5.9428 | Tidak |
| Public | Physical.Exercise2 or 4 days | -0.4305 | 0.2895 | -1.4869 | 0.1370 | 0.6502 | 0.3686 | 1.1468 | Tidak |
| Public | Physical.Exercise4 or 5 days | 0.0658 | 0.4526 | 0.1454 | 0.8844 | 1.0680 | 0.4399 | 2.5930 | Tidak |
| Public | Physical.ExerciseAlmost Everyday | 0.1211 | 0.3836 | 0.3158 | 0.7522 | 1.1288 | 0.5322 | 2.3938 | Tidak |
| Public | Physical.ExerciseI do not have | 0.1491 | 0.1906 | 0.7819 | 0.4342 | 1.1608 | 0.7988 | 1.6866 | Tidak |
| Public | Device.Usage3-5 hours | -0.3569 | 0.2987 | -1.1947 | 0.2322 | 0.6999 | 0.3897 | 1.2568 | Tidak |
| Public | Device.UsageMore than 5 hours | -0.2829 | 0.2934 | -0.9644 | 0.3348 | 0.7536 | 0.4241 | 1.3392 | Tidak |
| Public | Obesity.LevelNormal_Weight | -1.8144 | 0.3960 | -4.5818 | 0.0000 | 0.1629 | 0.0750 | 0.3541 | Ya |
| Public | Obesity.LevelObesity_Type_I | -2.4034 | 0.4061 | -5.9190 | 0.0000 | 0.0904 | 0.0408 | 0.2004 | Ya |
| Public | Obesity.LevelObesity_Type_II | -3.3262 | 0.4085 | -8.1436 | 0.0000 | 0.0359 | 0.0161 | 0.0800 | Ya |
| Public | Obesity.LevelObesity_Type_III | -3.1871 | 0.4099 | -7.7752 | 0.0000 | 0.0413 | 0.0185 | 0.0922 | Ya |
| Public | Obesity.LevelOverweight | -2.6349 | 0.3963 | -6.6491 | 0.0000 | 0.0717 | 0.0330 | 0.1560 | Ya |
Interpretasi: Variabel dengan p-value < 0,05 (hijau) berpengaruh signifikan. RRR > 1 berarti peningkatan odds relatif; RRR < 1 berarti penurunan odds relatif dibandingkan kategori referensi Private.
2.6 Interpretasi Koefisien & Relative Risk Ratio (RRR)
Kategori referensi: Private (kendaraan pribadi)
- Baris Active: log-odds memilih Active dibanding Private
- Baris Public: log-odds memilih Public dibanding Private
Jika \(\exp(\hat{\beta}) = 1{,}45\), individu pada kategori tersebut memiliki odds 1,45x lebih besar untuk memilih moda tersebut dibanding kendaraan pribadi, ceteris paribus.
result_multi %>%
filter(Signifikan == "Ya", Variabel != "(Intercept)") %>%
ggplot(aes(x = RRR, y = reorder(Variabel, RRR), color = Kategori)) +
geom_point(size = 3.5) +
geom_errorbarh(aes(xmin = CI_low, xmax = CI_high), height = 0.3) +
geom_vline(xintercept = 1, linetype = "dashed", color = "#b84f5a") +
scale_color_manual(values = c("#5568b8", "#90e0ef")) +
labs(
title = "Relative Risk Ratio (RRR) -- Variabel Signifikan",
subtitle = "Garis merah putus-putus: RRR = 1 (tidak ada efek)",
x = "RRR (dengan 95% CI)",
y = "Variabel",
color = "Kategori vs Private"
) +
theme_minimal(base_size = 12)
Interpretasi: Forest plot menampilkan RRR beserta 95% CI untuk variabel signifikan. CI yang tidak melewati RRR = 1 mengonfirmasi signifikansi statistik.
2.7 Uji Signifikansi Model (Likelihood Ratio Test)
model_null <- nnet::multinom(Transport.Mode ~ 1, data = data_multi, trace = FALSE)
G2 <- 2 * (logLik(fit_multi) - logLik(model_null))
df_G2 <- attr(logLik(fit_multi), "df") - attr(logLik(model_null), "df")
p_G2 <- pchisq(as.numeric(G2), df = df_G2, lower.tail = FALSE)
data.frame(
G2 = round(as.numeric(G2), 3),
df = df_G2,
p_value = signif(p_G2, 4)
) %>%
kable(caption = "Hasil Uji Likelihood Ratio Test (Uji G)") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| G2 | df | p_value |
|---|---|---|
| 611.305 | 40 | 0 |
Interpretasi: Jika p-value < 0,05, H_0 ditolak – minimal satu prediktor berpengaruh nyata terhadap probabilitas pemilihan moda transportasi.
2.8 Prediksi Probabilitas
pred_prob <- predict(fit_multi, type = "probs")
head(as.data.frame(pred_prob)) %>%
kable(digits = 4, caption = "Probabilitas Prediksi untuk 6 Observasi Pertama") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Private | Active | Public |
|---|---|---|
| 0.3935 | 0.0121 | 0.5944 |
| 0.1450 | 0.1779 | 0.6771 |
| 0.0901 | 0.3452 | 0.5646 |
| 0.3935 | 0.0121 | 0.5944 |
| 0.0569 | 0.3752 | 0.5679 |
| 0.1919 | 0.1803 | 0.6278 |
data_pred_multi <- data_multi %>%
bind_cols(as.data.frame(pred_prob)) %>%
mutate(prediksi = predict(fit_multi, type = "class"))
head(data_pred_multi[, c("Transport.Mode", "prediksi", "Active", "Public", "Private")]) %>%
kable(digits = 4, caption = "Perbandingan Aktual vs Prediksi beserta Probabilitas") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Transport.Mode | prediksi | Active | Public | Private |
|---|---|---|---|---|
| Private | Public | 0.0121 | 0.5944 | 0.3935 |
| Public | Public | 0.1779 | 0.6771 | 0.1450 |
| Public | Public | 0.3452 | 0.5646 | 0.0901 |
| Public | Public | 0.0121 | 0.5944 | 0.3935 |
| Active | Public | 0.3752 | 0.5679 | 0.0569 |
| Public | Public | 0.1803 | 0.6278 | 0.1919 |
Interpretasi: Setiap observasi mendapatkan tiga nilai probabilitas yang totalnya = 1. Kelas prediksi ditentukan berdasarkan kategori dengan probabilitas tertinggi.
2.9 Confusion Matrix & Evaluasi Model
2.9.1 Confusion Matrix
conf_multi <- table(Aktual = data_pred_multi$Transport.Mode,
Prediksi = data_pred_multi$prediksi)
conf_multi %>%
kable(caption = "Confusion Matrix Regresi Logistik Multinomial") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Private | Active | Public | |
|---|---|---|---|
| Private | 114 | 5 | 374 |
| Active | 2 | 24 | 346 |
| Public | 108 | 22 | 1187 |
Interpretasi: Elemen diagonal = prediksi benar. Elemen di luar diagonal = kesalahan klasifikasi antar kategori.
2.9.2 Metrik Evaluasi
accuracy_multi <- sum(diag(conf_multi)) / sum(conf_multi)
n_class <- nrow(conf_multi)
precision <- recall <- f1 <- numeric(n_class)
for (i in seq_len(n_class)) {
TP <- conf_multi[i, i]
FP <- sum(conf_multi[, i]) - TP
FN <- sum(conf_multi[i, ]) - TP
precision[i] <- ifelse((TP + FP) > 0, TP / (TP + FP), 0)
recall[i] <- ifelse((TP + FN) > 0, TP / (TP + FN), 0)
f1[i] <- ifelse((precision[i] + recall[i]) > 0,
2 * precision[i] * recall[i] / (precision[i] + recall[i]), 0)
}
data.frame(
Metrik = c("Overall Accuracy", "Macro Precision",
"Macro Recall (Sensitivity)", "Macro F1-Score"),
Nilai = round(c(accuracy_multi, mean(precision), mean(recall), mean(f1)), 4)
) %>%
kable(caption = "Metrik Evaluasi Model Regresi Logistik Multinomial") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Metrik | Nilai |
|---|---|
| Overall Accuracy | 0.6072 |
| Macro Precision | 0.5340 |
| Macro Recall (Sensitivity) | 0.3990 |
| Macro F1-Score | 0.3893 |
Interpretasi: Overall accuracy menunjukkan proporsi observasi yang diklasifikasikan benar. Macro F1-Score mencerminkan keseimbangan performa di semua kategori. Catatan: metrik ini dihitung pada data training; evaluasi lebih objektif memerlukan data testing terpisah.
2.10 Visualisasi Komposisi Prediksi
data_multi %>%
count(Obesity.Level, Gender, Transport.Mode) %>%
group_by(Obesity.Level, Gender) %>%
mutate(Proporsi = n / sum(n)) %>%
ungroup() %>%
ggplot(aes(x = Obesity.Level, y = Proporsi, fill = Transport.Mode)) +
geom_col(position = "fill", alpha = 0.90) +
facet_wrap(~Gender) +
scale_y_continuous(labels = percent_format()) +
scale_fill_manual(values = c("#0077b6", "#5568b8", "#90e0ef")) +
labs(
title = "Proporsi Moda Transportasi berdasarkan Obesitas dan Gender",
subtitle = "Komposisi pilihan transportasi di tiap kelompok obesitas",
x = "Tingkat Obesitas", y = "Proporsi", fill = "Moda Transportasi"
) +
theme_minimal(base_size = 12) +
theme(axis.text.x = element_text(angle = 25, hjust = 1))
Interpretasi: Stacked bar menampilkan komposisi moda transportasi pada setiap kombinasi tingkat obesitas dan gender. Perbedaan pola mencerminkan hubungan antara status berat badan dengan pilihan moda transportasi.
2.11 Asumsi Regresi Logistik Multinomial
- Respons bersifat nominal – Kategori Active, Public, Private tidak memiliki urutan alami.
- Observasi independen – Setiap baris mewakili individu berbeda.
- Tidak ada multikolinearitas berat – Prediktor kategorik dapat diperiksa dengan tabel kontingensi.
- Tidak ada perfect separation – Koefisien terdefinisi dengan baik.
- Hubungan linear pada skala logit – Tidak perlu diperiksa khusus karena seluruh prediktor kategorik.
- Independence of Irrelevant Alternatives (IIA) – Odds memilih Active vs Private tidak berubah substansial karena keberadaan Public. Jika IIA tidak masuk akal, nested logit dapat dipertimbangkan.
2.12 Kesimpulan Bagian II
- Regresi logistik multinomial berhasil memodelkan moda transportasi (Active, Public, Private) berdasarkan karakteristik demografis, kebiasaan makan, dan kondisi fisik – menggunakan Dhaka Obesity Dataset.
- Kategori Private sebagai referensi; koefisien dan RRR diinterpretasikan sebagai perbandingan peluang Active atau Public relatif terhadap Private.
- Berdasarkan LRT, model secara simultan signifikan dalam menjelaskan variasi pilihan moda transportasi.
- RRR memberikan gambaran substantif arah dan besarnya pengaruh masing-masing prediktor.
- Model dapat dikembangkan dengan pengujian asumsi IIA, penambahan variabel kontinu (usia, BMI), atau perbandingan dengan nested logit.
3 Regresi Logistik Ordinal
3.1 Pendahuluan
Latar Belakang
Kepuasan mahasiswa terhadap teknologi pembelajaran daring menjadi topik penting seiring meningkatnya penggunaan platform digital di perguruan tinggi. Kepuasan bersifat bertingkat (ordinal) – tidak sekadar “puas” atau “tidak puas”, melainkan memiliki gradasi bermakna. Oleh karena itu, regresi logistik ordinal merupakan pendekatan yang lebih tepat dibandingkan regresi linear biasa.
Tujuan analisis:
- Mendeskripsikan distribusi kepuasan mahasiswa terhadap platform pembelajaran daring.
- Mengidentifikasi faktor-faktor yang memengaruhi tingkat kepuasan.
- Membangun dan menginterpretasikan model regresi logistik ordinal (proportional odds model).
- Memeriksa asumsi proportional odds.
Sumber data: Fetaji, B. (2024). Dataset of Online Education Technologies and Resources for Learning After COVID19 Pandemic. Mendeley Data, V1. DOI: 10.17632/bf42yntpwr.1. Lisensi CC BY 4.0.
3.2 Memuat dan Membersihkan Data
# Ganti path sesuai lokasi file Anda
student_path <- "C:/Users/hp/Downloads/student_dataset.csv"
if (!file.exists(student_path)) stop("File 'student_dataset.csv' tidak ditemukan di: ", student_path)
data_raw_ord <- read.csv(student_path, stringsAsFactors = FALSE)
cat("Data berhasil dimuat.\n")## Data berhasil dimuat.## Dimensi: 120 baris x 14 kolom## Rows: 120
## Columns: 14
## $ Student_ID <chr> "S001", "S002", "S003", "S004", "S005", "S…
## $ University <chr> "South East European University", "South E…
## $ Year_of_Study <chr> "2nd", "4th", "4th", "3rd", "3rd", "4th", …
## $ Major <chr> "Economics", "Engineering", "Law", "Engine…
## $ Gender <chr> "Female", "Female", "Male", "Female", "Fem…
## $ Age <int> 23, 19, 19, 25, 18, 24, 25, 19, 18, 24, 25…
## $ Online_Tool_Used <chr> "Blackboard Collaborate", "Moodle", "MS Te…
## $ Usage_Frequency <chr> "Monthly", "Monthly", "Monthly", "Daily", …
## $ Satisfaction_Score <int> 4, 4, 5, 5, 5, 4, 3, 3, 4, 5, 3, 3, 5, 3, …
## $ Accessibility_Score <int> 3, 5, 4, 3, 3, 5, 4, 4, 4, 5, 4, 4, 5, 3, …
## $ Impact_on_Learning <int> 5, 4, 5, 3, 4, 4, 4, 5, 4, 3, 5, 4, 5, 4, …
## $ Challenges_Faced <chr> "Security concerns", "None", "Learning cur…
## $ Preferred_Features <chr> "Familiarity", "Simplicity", "User interfa…
## $ Suggestions_for_Improvement <chr> "Add more interactive tools", "N/A", "Enha…3.3 Rekayasa Fitur (Feature Engineering)
data_clean_ord <- data_raw_ord %>%
clean_names() %>%
mutate(
# Variabel Respons
kepuasan = case_when(
satisfaction_score %in% c(1, 2, 3) ~ "Cukup Puas",
satisfaction_score == 4 ~ "Puas",
satisfaction_score == 5 ~ "Sangat Puas"
),
kepuasan = factor(kepuasan,
levels = c("Cukup Puas", "Puas", "Sangat Puas"),
ordered = TRUE),
# Prediktor Kategorik
# frekuensi dibuat factor biasa (bukan ordered) agar polr() stabil
frekuensi = factor(usage_frequency, levels = c("Monthly", "Weekly", "Daily")),
gender = factor(gender),
# Prediktor Numerik (distandarisasi z-score)
accessibility = as.numeric(scale(accessibility_score)),
dampak_belajar = as.numeric(scale(impact_on_learning)),
usia = as.numeric(scale(age)),
tahun_studi = as.numeric(scale(
dplyr::recode(year_of_study,
`1st` = 1L, `2nd` = 2L, `3rd` = 3L, `4th` = 4L)
))
) %>%
dplyr::select(kepuasan, accessibility, dampak_belajar, frekuensi,
usia, tahun_studi, gender)
cat("=== Distribusi variabel respons ===\n")## === Distribusi variabel respons ===##
## Cukup Puas Puas Sangat Puas
## 34 38 48##
## === Tidak ada NA? ===## kepuasan accessibility dampak_belajar frekuensi usia
## 0 0 0 0 0
## tahun_studi gender
## 0 0## Rows: 120
## Columns: 7
## $ kepuasan <ord> Puas, Puas, Sangat Puas, Sangat Puas, Sangat Puas, Puas…
## $ accessibility <dbl> -1.1280290, 1.2891761, 0.0805735, -1.1280290, -1.128029…
## $ dampak_belajar <dbl> 1.0868859, -0.1552694, 1.0868859, -1.3974247, -0.155269…
## $ frekuensi <fct> Monthly, Monthly, Monthly, Daily, Monthly, Weekly, Week…
## $ usia <dbl> 0.47135975, -1.11082681, -1.11082681, 1.26245302, -1.50…
## $ tahun_studi <dbl> -0.5258203, 1.1568047, 1.1568047, 0.3154922, 0.3154922,…
## $ gender <fct> Female, Female, Male, Female, Female, Female, Male, Mal…3.4 Analisis Deskriptif
3.4.1 Ringkasan Statistik
## kepuasan accessibility dampak_belajar frekuensi
## Cukup Puas :34 Min. :-1.12803 Min. :-1.3974 Monthly:40
## Puas :38 1st Qu.:-1.12803 1st Qu.:-1.3974 Weekly :33
## Sangat Puas:48 Median : 0.08057 Median :-0.1553 Daily :47
## Mean : 0.00000 Mean : 0.0000
## 3rd Qu.: 1.28918 3rd Qu.: 1.0869
## Max. : 1.28918 Max. : 1.0869
## usia tahun_studi gender
## Min. :-1.50637 Min. :-1.3671 Female:65
## 1st Qu.:-1.11083 1st Qu.:-0.7361 Male :55
## Median : 0.07581 Median : 0.3155
## Mean : 0.00000 Mean : 0.0000
## 3rd Qu.: 0.86691 3rd Qu.: 1.1568
## Max. : 1.26245 Max. : 1.15683.4.2 Distribusi Variabel Respons
dist_kepuasan <- data_clean_ord %>%
count(kepuasan) %>%
mutate(proporsi = n / sum(n))
dist_kepuasan %>%
mutate(proporsi_fmt = percent(proporsi, accuracy = 0.1)) %>%
kable(
caption = "Tabel 1. Distribusi Tingkat Kepuasan Mahasiswa",
col.names = c("Tingkat Kepuasan", "Frekuensi (n)", "Proporsi", "Proporsi (%)"),
align = "lrrr"
) %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE, position = "left"
) %>%
row_spec(0, bold = TRUE, background = "#2c3e50", color = "white")| Tingkat Kepuasan | Frekuensi (n) | Proporsi | Proporsi (%) |
|---|---|---|---|
| Cukup Puas | 34 | 0.2833333 | 28.3% |
| Puas | 38 | 0.3166667 | 31.7% |
| Sangat Puas | 48 | 0.4000000 | 40.0% |
ggplot(dist_kepuasan, aes(x = kepuasan, y = proporsi, fill = kepuasan)) +
geom_col(width = 0.65, alpha = 0.92, show.legend = FALSE) +
geom_text(aes(label = percent(proporsi, accuracy = 0.1)),
vjust = -0.45, fontface = "bold", size = 4, color = "#2c3e50") +
scale_y_continuous(labels = percent_format(), limits = c(0, 0.65),
expand = expansion(mult = c(0, 0.05))) +
scale_fill_manual(values = c(
"Cukup Puas" = "#e67e22",
"Puas" = "#3498db",
"Sangat Puas" = "#27ae60"
)) +
labs(
title = "Distribusi Tingkat Kepuasan Mahasiswa",
subtitle = "Respons ordinal: Cukup Puas < Puas < Sangat Puas",
x = NULL, y = "Proporsi"
) +
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5, color = "grey50"),
panel.grid.major.x = element_blank(),
axis.text.x = element_text(face = "bold")
)
3.4.3 Eksplorasi Prediktor
3.4.3.1 Kepuasan berdasarkan Gender
data_clean_ord %>%
count(gender, kepuasan) %>%
group_by(gender) %>%
mutate(proporsi = n / sum(n)) %>%
ungroup() %>%
ggplot(aes(x = kepuasan, y = proporsi, fill = gender)) +
geom_col(position = "dodge", width = 0.65, alpha = 0.9) +
geom_text(aes(label = percent(proporsi, accuracy = 1)),
position = position_dodge(width = 0.65),
vjust = -0.4, size = 3.2, fontface = "bold") +
scale_y_continuous(labels = percent_format(), limits = c(0, 0.65),
expand = expansion(mult = c(0, 0.05))) +
scale_fill_manual(values = c("Male" = "#3498db", "Female" = "#e91e8c")) +
labs(title = "Distribusi Kepuasan berdasarkan Gender",
x = NULL, y = "Proporsi", fill = "Gender") +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold", hjust = 0.5),
legend.position = "top",
panel.grid.major.x = element_blank())
3.4.3.2 Kepuasan berdasarkan Frekuensi Penggunaan
data_clean_ord %>%
count(frekuensi, kepuasan) %>%
group_by(frekuensi) %>%
mutate(proporsi = n / sum(n)) %>%
ungroup() %>%
ggplot(aes(x = kepuasan, y = proporsi, fill = frekuensi)) +
geom_col(position = "dodge", width = 0.65, alpha = 0.9) +
scale_y_continuous(labels = percent_format(), limits = c(0, 0.70),
expand = expansion(mult = c(0, 0.05))) +
scale_fill_manual(values = c(
"Monthly" = "#95a5a6",
"Weekly" = "#f39c12",
"Daily" = "#27ae60"
)) +
labs(title = "Distribusi Kepuasan berdasarkan Frekuensi Penggunaan",
x = NULL, y = "Proporsi", fill = "Frekuensi") +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold", hjust = 0.5),
legend.position = "top",
panel.grid.major.x = element_blank())
3.4.3.3 Distribusi Skor Numerik
p_access <- ggplot(data_clean_ord,
aes(x = kepuasan, y = accessibility, fill = kepuasan)) +
geom_boxplot(alpha = 0.8, outlier.shape = 21, show.legend = FALSE) +
scale_fill_manual(values = c(
"Cukup Puas" = "#e67e22", "Puas" = "#3498db", "Sangat Puas" = "#27ae60"
)) +
labs(title = "Accessibility Score", x = NULL, y = "Skor (z-score)") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold", hjust = 0.5),
axis.text.x = element_text(angle = 15, hjust = 1))
p_impact <- ggplot(data_clean_ord,
aes(x = kepuasan, y = dampak_belajar, fill = kepuasan)) +
geom_boxplot(alpha = 0.8, outlier.shape = 21, show.legend = FALSE) +
scale_fill_manual(values = c(
"Cukup Puas" = "#e67e22", "Puas" = "#3498db", "Sangat Puas" = "#27ae60"
)) +
labs(title = "Impact on Learning", x = NULL, y = "Skor (z-score)") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold", hjust = 0.5),
axis.text.x = element_text(angle = 15, hjust = 1))
p_access + p_impact +
plot_annotation(
title = "Distribusi Skor Numerik berdasarkan Tingkat Kepuasan",
theme = theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 13))
)
3.5 Model Regresi Logistik Ordinal
3.5.1 Landasan Teori
Model proportional odds (McCullagh, 1980):
\[ \text{logit}\{P(Y \leq j \mid \mathbf{x})\} = \alpha_j - \mathbf{x}^\top \boldsymbol{\beta}, \quad j = 1, 2 \]
di mana \(\alpha_j\) adalah cutpoint ke-\(j\) dan \(\boldsymbol{\beta}\) adalah vektor koefisien yang sama untuk setiap batas kumulatif. Probabilitas setiap kategori:
\[ P(Y = j) = P(Y \leq j) - P(Y \leq j-1) \]
3.5.2 Spesifikasi Model
Variabel Respons: \(\text{kepuasan}: \text{Cukup Puas} < \text{Puas} < \text{Sangat Puas}\)
| Prediktor | Tipe | Deskripsi |
|---|---|---|
accessibility |
Numerik (z-score) | Skor kemudahan akses platform |
dampak_belajar |
Numerik (z-score) | Skor dampak terhadap pembelajaran |
frekuensi |
Kategorik 3 level | Monthly / Weekly / Daily (ref: Monthly) |
usia |
Numerik (z-score) | Usia mahasiswa |
tahun_studi |
Numerik (z-score) | Tahun studi (1-4) |
gender |
Nominal | Jenis kelamin (ref: Female) |
Catatan Teknis – Standardisasi & Tipe Faktor
Prediktor numerik distandarisasi (z-score) agar skala seragam dan
optimasi polr() lebih stabil. frekuensi dibuat
faktor biasa (bukan ordered factor) karena polr()
mengekspansi ordered factor menjadi polynomial
contrasts yang rawan gagal pada sampel kecil.
3.5.3 Estimasi Model
## === Distribusi kepuasan ===##
## Cukup Puas Puas Sangat Puas
## 34 38 48# Starting values manual
probs_kum <- cumsum(prop.table(table(data_clean_ord$kepuasan)))
start_alpha <- qlogis(probs_kum[-length(probs_kum)])
n_beta <- 4 + 2 + 1 # 4 numerik + 2 dummy frekuensi + 1 dummy gender
start_vals <- c(rep(0, n_beta), start_alpha)
fit_ord <- tryCatch(
MASS::polr(
kepuasan ~ accessibility + dampak_belajar + frekuensi +
usia + tahun_studi + gender,
data = data_clean_ord,
method = "logistic",
Hess = TRUE,
start = start_vals
),
error = function(e) {
message("Percobaan pertama gagal: ", conditionMessage(e))
message("Mencoba tanpa starting values eksplisit...")
MASS::polr(
kepuasan ~ accessibility + dampak_belajar + frekuensi +
usia + tahun_studi + gender,
data = data_clean_ord,
method = "logistic",
Hess = TRUE
)
}
)
summary(fit_ord)## Call:
## MASS::polr(formula = kepuasan ~ accessibility + dampak_belajar +
## frekuensi + usia + tahun_studi + gender, data = data_clean_ord,
## start = start_vals, Hess = TRUE, method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## accessibility -0.07711 0.1722 -0.44790
## dampak_belajar 0.04982 0.1740 0.28622
## frekuensiWeekly -0.66592 0.4469 -1.49003
## frekuensiDaily -0.09403 0.4100 -0.22933
## usia -0.06799 0.1812 -0.37520
## tahun_studi -0.06771 0.1729 -0.39168
## genderMale 0.02629 0.3461 0.07595
##
## Intercepts:
## Value Std. Error t value
## Cukup Puas|Puas -1.1553 0.3702 -3.1212
## Puas|Sangat Puas 0.2124 0.3539 0.6002
##
## Residual Deviance: 257.4331
## AIC: 275.43313.5.4 Ringkasan Koefisien
ord_coef <- coef(summary(fit_ord))
result_ord <- as.data.frame(ord_coef) %>%
rownames_to_column("parameter") %>%
rename(estimasi = Value, se = `Std. Error`, t_val = `t value`) %>%
mutate(
p_value = 2 * (1 - pnorm(abs(t_val))),
OR = exp(estimasi),
CI_low = exp(estimasi - 1.96 * se),
CI_high = exp(estimasi + 1.96 * se),
jenis = ifelse(grepl("\\|", parameter), "Cutpoint", "Koefisien"),
sig = case_when(
p_value < 0.001 ~ "***",
p_value < 0.01 ~ "**",
p_value < 0.05 ~ "*",
p_value < 0.1 ~ ".",
TRUE ~ ""
)
)
result_ord %>%
mutate(across(c(estimasi, se, t_val, p_value, OR, CI_low, CI_high),
~ round(.x, 4))) %>%
kable(
caption = "Tabel 2. Ringkasan Hasil Regresi Logistik Ordinal",
col.names = c("Parameter", "Estimasi", "SE", "z/t", "p-value",
"OR", "CI 2.5%", "CI 97.5%", "Jenis", "Sig."),
align = "lrrrrrrrcc"
) %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE
) %>%
row_spec(0, bold = TRUE, background = "#2c3e50", color = "white") %>%
pack_rows("Koefisien Prediktor",
min(which(result_ord$jenis == "Koefisien")),
max(which(result_ord$jenis == "Koefisien"))) %>%
pack_rows("Cutpoint (Intercept)",
min(which(result_ord$jenis == "Cutpoint")),
max(which(result_ord$jenis == "Cutpoint")),
background = "#ecf0f1") %>%
footnote(
general = "Signifikansi: *** p<0.001 ** p<0.01 * p<0.05 . p<0.1",
general_title = "Catatan:",
footnote_as_chunk = TRUE
)| Parameter | Estimasi | SE | z/t | p-value | OR | CI 2.5% | CI 97.5% | Jenis | Sig. |
|---|---|---|---|---|---|---|---|---|---|
| Koefisien Prediktor | |||||||||
| accessibility | -0.0771 | 0.1722 | -0.4479 | 0.6542 | 0.9258 | 0.6606 | 1.2974 | Koefisien | |
| dampak_belajar | 0.0498 | 0.1740 | 0.2862 | 0.7747 | 1.0511 | 0.7473 | 1.4784 | Koefisien | |
| frekuensiWeekly | -0.6659 | 0.4469 | -1.4900 | 0.1362 | 0.5138 | 0.2140 | 1.2337 | Koefisien | |
| frekuensiDaily | -0.0940 | 0.4100 | -0.2293 | 0.8186 | 0.9103 | 0.4075 | 2.0333 | Koefisien | |
| usia | -0.0680 | 0.1812 | -0.3752 | 0.7075 | 0.9343 | 0.6550 | 1.3326 | Koefisien | |
| tahun_studi | -0.0677 | 0.1729 | -0.3917 | 0.6953 | 0.9345 | 0.6660 | 1.3114 | Koefisien | |
| genderMale | 0.0263 | 0.3461 | 0.0759 | 0.9395 | 1.0266 | 0.5209 | 2.0233 | Koefisien | |
| Cutpoint (Intercept) | |||||||||
| Cukup Puas|Puas | -1.1553 | 0.3702 | -3.1212 | 0.0018 | 0.3149 | 0.1525 | 0.6506 | Cutpoint | ** |
| Puas|Sangat Puas | 0.2124 | 0.3539 | 0.6002 | 0.5484 | 1.2367 | 0.6180 | 2.4748 | Cutpoint | |
| Catatan: Signifikansi: *** p<0.001 ** p<0.01 * p<0.05 . p<0.1 | |||||||||
3.5.5 Visualisasi Odds Ratio
koef_plot <- result_ord %>%
filter(jenis == "Koefisien") %>%
mutate(
parameter = dplyr::recode(parameter,
"accessibility" = "Accessibility Score",
"dampak_belajar" = "Impact on Learning",
"frekuensiWeekly" = "Frekuensi: Weekly",
"frekuensiDaily" = "Frekuensi: Daily",
"usia" = "Usia",
"tahun_studi" = "Tahun Studi",
"genderMale" = "Gender: Male"
),
signifikan = p_value < 0.05
)
ggplot(koef_plot, aes(x = OR, y = reorder(parameter, OR), color = signifikan)) +
geom_vline(xintercept = 1, linetype = "dashed",
color = "grey60", linewidth = 0.8) +
geom_errorbarh(aes(xmin = CI_low, xmax = CI_high),
height = 0.25, linewidth = 0.9) +
geom_point(size = 4, shape = 18) +
scale_color_manual(
values = c("TRUE" = "#e74c3c", "FALSE" = "#95a5a6"),
labels = c("TRUE" = "Signifikan (p<0.05)", "FALSE" = "Tidak Signifikan"),
name = NULL
) +
scale_x_log10(breaks = c(0.25, 0.5, 1, 2, 4)) +
labs(
title = "Odds Ratio dengan Interval Kepercayaan 95%",
subtitle = "Skala log; garis putus-putus = OR = 1 (tidak ada efek)",
x = "Odds Ratio (skala log)", y = NULL
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5, color = "grey50"),
legend.position = "bottom",
panel.grid.minor = element_blank()
)
3.6 Interpretasi Hasil
Panduan Membaca OR dalam Model Ordinal
- OR > 1: kenaikan satu unit prediktor meningkatkan peluang berada di kategori kepuasan yang lebih tinggi
- OR < 1: kenaikan satu unit prediktor menurunkan peluang berada di kategori kepuasan yang lebih tinggi
- OR berlaku sama untuk setiap batas kumulatif (asumsi proportional odds)
| Prediktor | Arah Efek | OR | p-value | Interpretasi | Sig. |
|---|---|---|---|---|---|
| Accessibility Score | Negatif (turun) | 0.926 | 0.6542 | OR = 0.926; setiap kenaikan 1 SD Accessibility Score mengubah odds kategori lebih tinggi sebesar 7.4% (menurun) | |
| Impact on Learning | Positif (naik) | 1.051 | 0.7747 | OR = 1.051; setiap kenaikan 1 SD Impact on Learning mengubah odds kategori lebih tinggi sebesar 5.1% (meningkat) | |
| Frekuensi: Weekly (vs Monthly) | Negatif (turun) | 0.514 | 0.1362 | OR = 0.514; setiap kenaikan 1 SD Frekuensi: Weekly (vs Monthly) mengubah odds kategori lebih tinggi sebesar 48.6% (menurun) | |
| Frekuensi: Daily (vs Monthly) | Negatif (turun) | 0.910 | 0.8186 | OR = 0.91; setiap kenaikan 1 SD Frekuensi: Daily (vs Monthly) mengubah odds kategori lebih tinggi sebesar 9% (menurun) | |
| Usia | Negatif (turun) | 0.934 | 0.7075 | OR = 0.934; setiap kenaikan 1 SD Usia mengubah odds kategori lebih tinggi sebesar 6.6% (menurun) | |
| Tahun Studi | Negatif (turun) | 0.935 | 0.6953 | OR = 0.935; setiap kenaikan 1 SD Tahun Studi mengubah odds kategori lebih tinggi sebesar 6.5% (menurun) | |
| Gender: Male (vs Female) | Positif (naik) | 1.027 | 0.9395 | OR = 1.027; setiap kenaikan 1 SD Gender: Male (vs Female) mengubah odds kategori lebih tinggi sebesar 2.7% (meningkat) |
3.7 Prediksi Probabilitas
3.7.1 Prediksi pada Data Asli
pred_probs_ord <- predict(fit_ord, type = "probs")
pred_class_ord <- predict(fit_ord, type = "class")
data_pred_ord <- data_clean_ord %>%
bind_cols(as.data.frame(pred_probs_ord)) %>%
mutate(prediksi = pred_class_ord)
head(data_pred_ord, 8) %>%
mutate(across(c(`Cukup Puas`, `Puas`, `Sangat Puas`), ~ round(.x, 3))) %>%
kable(caption = "Tabel 4. Contoh Prediksi Probabilitas per Kategori (8 Observasi Pertama)") %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE
) %>%
row_spec(0, bold = TRUE, background = "#2c3e50", color = "white")| kepuasan | accessibility | dampak_belajar | frekuensi | usia | tahun_studi | gender | Cukup Puas | Puas | Sangat Puas | prediksi |
|---|---|---|---|---|---|---|---|---|---|---|
| Puas | -1.1280290 | 1.0868859 | Monthly | 0.4713597 | -0.5258203 | Female | 0.214 | 0.303 | 0.483 | Sangat Puas |
| Puas | 1.2891761 | -0.1552694 | Monthly | -1.1108268 | 1.1568047 | Female | 0.260 | 0.320 | 0.420 | Sangat Puas |
| Sangat Puas | 0.0805735 | 1.0868859 | Monthly | -1.1108268 | 1.1568047 | Male | 0.227 | 0.308 | 0.465 | Sangat Puas |
| Sangat Puas | -1.1280290 | -1.3974247 | Daily | 1.2624530 | 0.3154922 | Female | 0.275 | 0.323 | 0.402 | Sangat Puas |
| Sangat Puas | -1.1280290 | -0.1552694 | Monthly | -1.5063735 | 0.3154922 | Female | 0.212 | 0.301 | 0.487 | Sangat Puas |
| Puas | 1.2891761 | -0.1552694 | Weekly | 0.8669064 | 1.1568047 | Female | 0.439 | 0.315 | 0.245 | Cukup Puas |
| Cukup Puas | 0.0805735 | -0.1552694 | Weekly | 1.2624530 | 1.1568047 | Male | 0.416 | 0.321 | 0.263 | Cukup Puas |
| Cukup Puas | 0.0805735 | 1.0868859 | Weekly | -1.1108268 | 0.3154922 | Male | 0.350 | 0.329 | 0.321 | Cukup Puas |
3.7.2 Visualisasi Probabilitas terhadap Accessibility Score
grid_access <- expand.grid(
accessibility = seq(min(data_clean_ord$accessibility),
max(data_clean_ord$accessibility), length.out = 100),
dampak_belajar = median(data_clean_ord$dampak_belajar),
frekuensi = factor("Weekly", levels = levels(data_clean_ord$frekuensi)),
usia = mean(data_clean_ord$usia),
tahun_studi = mean(data_clean_ord$tahun_studi),
gender = factor("Female", levels = levels(data_clean_ord$gender))
)
grid_probs_access <- predict(fit_ord, newdata = grid_access, type = "probs")
grid_access %>%
bind_cols(as.data.frame(grid_probs_access)) %>%
pivot_longer(cols = c("Cukup Puas", "Puas", "Sangat Puas"),
names_to = "kepuasan", values_to = "probabilitas") %>%
mutate(kepuasan = factor(kepuasan,
levels = c("Cukup Puas", "Puas", "Sangat Puas"))) %>%
ggplot(aes(x = accessibility, y = probabilitas, color = kepuasan)) +
geom_line(linewidth = 1.3) +
scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
scale_color_manual(values = c(
"Cukup Puas" = "#e67e22", "Puas" = "#3498db", "Sangat Puas" = "#27ae60"
)) +
labs(
title = "Probabilitas Prediksi Kepuasan vs. Accessibility Score",
subtitle = "Prediktor lain ditahan pada nilai rata-rata; gender = Female; frekuensi = Weekly",
x = "Accessibility Score (z-score)", y = "Probabilitas Prediksi",
color = "Tingkat Kepuasan"
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5, color = "grey50"),
legend.position = "bottom")
3.7.3 Visualisasi Probabilitas terhadap Impact on Learning
grid_impact <- expand.grid(
accessibility = median(data_clean_ord$accessibility),
dampak_belajar = seq(min(data_clean_ord$dampak_belajar),
max(data_clean_ord$dampak_belajar), length.out = 100),
frekuensi = factor("Weekly", levels = levels(data_clean_ord$frekuensi)),
usia = mean(data_clean_ord$usia),
tahun_studi = mean(data_clean_ord$tahun_studi),
gender = factor("Female", levels = levels(data_clean_ord$gender))
)
grid_probs_impact <- predict(fit_ord, newdata = grid_impact, type = "probs")
grid_impact %>%
bind_cols(as.data.frame(grid_probs_impact)) %>%
pivot_longer(cols = c("Cukup Puas", "Puas", "Sangat Puas"),
names_to = "kepuasan", values_to = "probabilitas") %>%
mutate(kepuasan = factor(kepuasan,
levels = c("Cukup Puas", "Puas", "Sangat Puas"))) %>%
ggplot(aes(x = dampak_belajar, y = probabilitas, color = kepuasan)) +
geom_line(linewidth = 1.3) +
scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
scale_color_manual(values = c(
"Cukup Puas" = "#e67e22", "Puas" = "#3498db", "Sangat Puas" = "#27ae60"
)) +
labs(
title = "Probabilitas Prediksi Kepuasan vs. Impact on Learning",
subtitle = "Prediktor lain ditahan pada nilai rata-rata; gender = Female; frekuensi = Weekly",
x = "Impact on Learning Score (z-score)", y = "Probabilitas Prediksi",
color = "Tingkat Kepuasan"
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5, color = "grey50"),
legend.position = "bottom")
3.7.4 Grafik Cumulative Logit (Uji Visual Proportional Odds)
eps <- 1e-6
grid_access %>%
bind_cols(as.data.frame(grid_probs_access)) %>%
mutate(
`P(Y <= Cukup Puas)` = `Cukup Puas`,
`P(Y <= Puas)` = `Cukup Puas` + `Puas`
) %>%
pivot_longer(cols = c(`P(Y <= Cukup Puas)`, `P(Y <= Puas)`),
names_to = "batas", values_to = "prob_kum") %>%
mutate(
prob_kum = pmin(pmax(prob_kum, eps), 1 - eps),
cum_logit = qlogis(prob_kum),
batas = factor(batas, levels = c("P(Y <= Cukup Puas)", "P(Y <= Puas)"))
) %>%
ggplot(aes(x = accessibility, y = cum_logit, color = batas)) +
geom_line(linewidth = 1.3) +
scale_color_manual(values = c(
"P(Y <= Cukup Puas)" = "#e67e22",
"P(Y <= Puas)" = "#3498db"
)) +
labs(
title = "Grafik Cumulative Logit -- Uji Visual Proportional Odds",
subtitle = "Garis sejajar (paralel) mengindikasikan asumsi proportional odds terpenuhi",
x = "Accessibility Score (z-score)", y = "Logit Peluang Kumulatif",
color = "Batas Kumulatif"
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5, color = "grey50"),
legend.position = "bottom")
Catatan Grafik Cumulative Logit
Garis yang sejajar (paralel) pada skala logit kumulatif mengonfirmasi asumsi proportional odds terpenuhi: setiap prediktor memiliki efek yang sama terhadap semua batas kumulatif, hanya berbeda pada nilai intercept (\(\alpha_j\)).
3.8 Evaluasi Model
3.8.1 Confusion Matrix
conf_mat_ord <- table(
Aktual = data_pred_ord$kepuasan,
Prediksi = data_pred_ord$prediksi
)
conf_mat_ord %>%
kable(caption = "Tabel 5. Confusion Matrix") %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE
) %>%
row_spec(0, bold = TRUE, background = "#2c3e50", color = "white")| Cukup Puas | Puas | Sangat Puas | |
|---|---|---|---|
| Cukup Puas | 13 | 0 | 21 |
| Puas | 10 | 0 | 28 |
| Sangat Puas | 8 | 0 | 40 |
accuracy_ord <- sum(diag(conf_mat_ord)) / sum(conf_mat_ord)
cat("Akurasi keseluruhan:", percent(accuracy_ord, accuracy = 0.1), "\n")## Akurasi keseluruhan: 44.2%Catatan tentang Evaluasi Model Ordinal
Akurasi tunggal kurang informatif untuk respons ordinal karena kesalahan tidak sama beratnya. Ukuran evaluasi yang lebih tepat: mean absolute error ordinal, weighted kappa, atau ranked probability score.
3.8.2 Mean Absolute Error Ordinal
skor_ordinal <- c("Cukup Puas" = 1, "Puas" = 2, "Sangat Puas" = 3)
mae_ordinal <- data_pred_ord %>%
mutate(
skor_aktual = skor_ordinal[as.character(kepuasan)],
skor_prediksi = skor_ordinal[as.character(prediksi)],
ae = abs(skor_aktual - skor_prediksi)
) %>%
summarise(MAE = mean(ae))
cat("Mean Absolute Error (ordinal):", round(mae_ordinal$MAE, 4), "\n")## Mean Absolute Error (ordinal): 0.8## (MAE = 0: sempurna; MAE = 2: terburuk untuk 3 kategori)3.9 Pemeriksaan Asumsi Proportional Odds
3.9.1 Uji Formal dengan
ordinal::nominal_test()
fit_clm <- ordinal::clm(
kepuasan ~ accessibility + dampak_belajar + frekuensi +
usia + tahun_studi + gender,
data = data_clean_ord,
link = "logit"
)
nominal_test(fit_clm)Membaca Hasil nominal_test()
- p-value > 0.05: tidak ada bukti pelanggaran proportional odds untuk prediktor tersebut.
- p-value < 0.05: indikasi efek prediktor mungkin berbeda di setiap batas kumulatif – pertimbangkan model alternatif.
Interpretasi tidak boleh mekanis: perlu dipadukan dengan teori substantif, ukuran sampel, dan visualisasi cumulative logit.
3.10 Alternatif Model (Jika Asumsi Tidak Terpenuhi)
| Alternatif | Fungsi di R | Kapan Digunakan |
|---|---|---|
| Partial proportional odds | VGAM::vglm() |
Sebagian prediktor tidak proporsional |
| Adjacent-category logit | VGAM::vglm(acat()) |
Perbandingan kategori berdampingan lebih bermakna |
| Continuation-ratio logit | rms::lrm() |
Proses bertahap / sekuensial |
| Multinomial logistic | nnet::multinom() |
Tidak ada urutan yang bermakna |
# Contoh model partial proportional odds (tidak dijalankan)
# install.packages("VGAM")
library(VGAM)
fit_ppo <- vglm(
kepuasan ~ accessibility + dampak_belajar + frekuensi +
usia + tahun_studi + gender,
family = cumulative(parallel = FALSE ~ 1),
data = data_clean_ord
)
summary(fit_ppo)3.11 Kesimpulan Bagian III
Ringkasan Temuan
- Kepuasan mahasiswa dimodelkan dengan 3 kategori ordinal: Cukup Puas < Puas < Sangat Puas – mayoritas responden di kategori Puas dan Sangat Puas.
- Accessibility Score dan Impact on Learning secara konsisten berkaitan positif dengan kepuasan.
- Frekuensi penggunaan lebih sering (Daily > Weekly > Monthly) cenderung meningkatkan kepuasan.
- Gender dan usia memiliki efek lebih kecil dan kurang konsisten secara statistis.
- Pemeriksaan asumsi proportional odds (grafik cumulative
logit +
nominal_test()) tidak memberikan indikasi pelanggaran yang kuat. - Keterbatasan: ukuran sampel yang terbatas mengurangi power uji statistik; hasil sebaiknya divalidasi pada sampel lebih besar.
4 Regresi Poisson
4.1 Pendahuluan
4.1.1 Deskripsi Dataset
Analisis ini menggunakan Health Count Data yang tersedia secara terbuka di Mendeley Data.
Sumber Data
- Judul: Health Count Data
- Kontributor: Adesina, Olumide
- Publikasi: 26 Juni 2019 . Versi 1
- DOI: 10.17632/z7wznk53cf.1
- Lisensi: CC BY 4.0 (open access)
- Artikel terkait: “Bayesian Models for Zero Truncated Count Data” – AJPAS
Data diperoleh dari fasilitas kesehatan di Ota, Ogun State, Nigeria. Sampel terdiri dari 1.647 pasien peserta National Health Insurance Scheme (NHIS) selama Juli 2016-Juli 2017.
4.2 Variabel Penelitian
| Variabel | Tipe | Keterangan |
|---|---|---|
| Nencounter | Respons (count) | Jumlah kunjungan pasien ke dokter |
| Eclass | Biner | Status pasien: rawat inap (1), rawat jalan (0) |
| follow_up | Biner | Kontrol rutin: ya (1), tidak (0) |
| sex | Biner | Jenis kelamin: laki-laki (1), perempuan (0) |
| Ndiagnosis | Numerik | Jumlah diagnosis selama periode |
| age | Numerik | Usia biologis pasien (tahun) |
Catatan penting: Dataset ini merupakan zero-truncated count data – seluruh observasi memiliki \(Y \geq 1\) (tidak ada nilai nol). Kondisi ini terjadi karena data hanya mencatat pasien yang memang sudah datang setidaknya sekali. Regresi Poisson standar tetap dapat diterapkan, namun perlu dicatat bahwa model memodelkan jumlah kunjungan di antara pasien yang sudah berkunjung.
4.3 Pemuatan Data
# Ganti path sesuai lokasi file Anda
publ_path <- "C:/Users/hp/Downloads/PUBL.csv"
if (!file.exists(publ_path)) stop("File 'PUBL.csv' tidak ditemukan di: ", publ_path)
data_raw <- read.csv(publ_path, stringsAsFactors = FALSE)
cat("Dimensi data:", nrow(data_raw), "baris x", ncol(data_raw), "kolom\n")## Dimensi data: 126 baris x 8 kolom## Rows: 126
## Columns: 8
## $ NOP <int> 14, 9, 6, 3, 4, 13, 7, 0, 6, 4, 7, 7, 0, 0, 14, 2, 0, 1, 4, 8, 2…
## $ SEX <int> 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1…
## $ MS <int> 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1…
## $ CHD <int> 3, 3, 2, 3, 0, 0, 2, 4, 3, 3, 3, 2, 2, 1, 4, 0, 1, 1, 1, 4, 2, 2…
## $ EXP <int> 13, 12, 14, 10, 9, 9, 12, 31, 9, 13, 10, 3, 7, 7, 12, 6, 6, 6, 5…
## $ LEVEL <int> 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1…
## $ UGC <int> 8, 4, 4, 5, 4, 4, 6, 6, 6, 4, 4, 6, 6, 6, 6, 6, 5, 5, 3, 7, 6, 6…
## $ PGC <int> 0, 3, 7, 0, 0, 5, 5, 6, 0, 0, 5, 5, 0, 0, 5, 0, 0, 0, 0, 5, 5, 5…4.4 Pra-pemrosesan
## Nama kolom dalam data:## [1] "NOP" "SEX" "MS" "CHD" "EXP" "LEVEL" "UGC" "PGC"# Standardisasi nama kolom: lowercase
data_std <- data_raw %>%
janitor::clean_names()
cat("\nNama kolom setelah standardisasi:\n")##
## Nama kolom setelah standardisasi:## [1] "nop" "sex" "ms" "chd" "exp" "level" "ugc" "pgc"# -- Pemetaan kolom PUBL.csv ke nama baku analisis ----------------------------
# Berdasarkan artikel "Bayesian Models for Zero Truncated Count Data":
# NOP = Nencounter (jumlah kunjungan ke dokter -- variabel respons)
# SEX = sex (jenis kelamin: 1=laki-laki, 0=perempuan)
# CHD = follow_up (kontrol rutin: 1=ya, 0=tidak)
# LEVEL = Eclass (kelas pasien: 1=rawat inap, 0=rawat jalan)
# UGC = Ndiagnosis (jumlah diagnosis)
# PGC = age (usia biologis pasien)
# MS dan EXP tidak digunakan dalam model ini.
# -----------------------------------------------------------------------------
col_map <- c(
Nencounter = "nop", # jumlah kunjungan pasien
sex = "sex", # jenis kelamin
follow_up = "chd", # kontrol rutin / follow-up
Eclass = "level", # kelas pasien (rawat inap / jalan)
Ndiagnosis = "ugc", # jumlah diagnosis
age = "pgc" # usia biologis
)
# Cek ketersediaan kolom
missing_cols <- col_map[!col_map %in% names(data_std)]
if (length(missing_cols) > 0) {
cat("\nKolom TIDAK DITEMUKAN -- periksa nama di col_map:\n")
print(missing_cols)
cat("\nKolom yang tersedia:\n")
print(names(data_std))
stop("Sesuaikan col_map dengan nama kolom aktual di atas.")
}
# Rename dan konversi ke faktor
health <- data_std %>%
dplyr::rename(
Nencounter = nop,
sex = sex,
follow_up = chd,
Eclass = level,
Ndiagnosis = ugc,
age = pgc
) %>%
dplyr::select(Nencounter, Eclass, follow_up, sex, Ndiagnosis, age) %>%
dplyr::mutate(
Eclass = factor(Eclass, levels = c(0, 1), labels = c("Rawat Jalan", "Rawat Inap")),
follow_up = factor(follow_up, levels = c(0, 1), labels = c("Tidak", "Ya")),
sex = factor(sex, levels = c(0, 1), labels = c("Perempuan", "Laki-laki"))
)
cat("\nDimensi setelah preprocessing:",
nrow(health), "baris x", ncol(health), "kolom\n")##
## Dimensi setelah preprocessing: 126 baris x 6 kolom## Rows: 126
## Columns: 6
## $ Nencounter <int> 14, 9, 6, 3, 4, 13, 7, 0, 6, 4, 7, 7, 0, 0, 14, 2, 0, 1, 4,…
## $ Eclass <fct> Rawat Inap, Rawat Inap, Rawat Inap, Rawat Inap, Rawat Jalan…
## $ follow_up <fct> NA, NA, NA, NA, Tidak, Tidak, NA, NA, NA, NA, NA, NA, NA, Y…
## $ sex <fct> Perempuan, Laki-laki, Laki-laki, Laki-laki, Perempuan, Laki…
## $ Ndiagnosis <int> 8, 4, 4, 5, 4, 4, 6, 6, 6, 4, 4, 6, 6, 6, 6, 6, 5, 5, 3, 7,…
## $ age <int> 0, 3, 7, 0, 0, 5, 5, 6, 0, 0, 5, 5, 0, 0, 5, 0, 0, 0, 0, 5,…Interpretasi: Dataset berhasil dimuat dan diproses. Variabel
Eclass,follow_up, dansexdikonversi menjadi faktor dengan kategori referensi masing-masing Rawat Jalan, Tidak (kontrol), dan Perempuan. Variabel responsNencounterberupa data hitung bulat nonnegatif >= 1.
4.5 Eksplorasi Data (EDA)
4.5.1 Statistik Deskriptif
## Nencounter Eclass follow_up sex Ndiagnosis
## Min. : 0.000 Rawat Jalan:56 Tidak:22 Perempuan:49 Min. : 3.000
## 1st Qu.: 2.000 Rawat Inap :70 Ya :17 Laki-laki:77 1st Qu.: 6.000
## Median : 6.000 NA's :87 Median : 6.000
## Mean : 6.444 Mean : 7.873
## 3rd Qu.:10.000 3rd Qu.:10.000
## Max. :32.000 Max. :27.000
## age
## Min. : 0.000
## 1st Qu.: 0.000
## Median : 1.000
## Mean : 2.381
## 3rd Qu.: 5.000
## Max. :10.000desc_num <- health %>%
dplyr::summarise(
`Rata-rata` = round(mean(Nencounter), 3),
`Median` = median(Nencounter),
`Varians` = round(var(Nencounter), 3),
`Min` = min(Nencounter),
`Maks` = max(Nencounter),
`SD` = round(sd(Nencounter), 3)
)
desc_num %>%
knitr::kable(caption = "Statistik Deskriptif Variabel Respons (Nencounter)") %>%
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE
)| Rata-rata | Median | Varians | Min | Maks | SD |
|---|---|---|---|---|---|
| 6.444 | 6 | 31.065 | 0 | 32 | 5.574 |
Interpretasi: Rata-rata jumlah kunjungan pasien adalah 6.44 kunjungan, dengan varians 31.065. Rasio varians terhadap mean = 4.82. Nilai ini lebih kecil dari 1, mengindikasikan kemungkinan underdispersion yang akan diperiksa lebih lanjut setelah model diestimasi.
4.6 Distribusi Variabel Respons
dist_enc <- health %>%
dplyr::count(Nencounter) %>%
dplyr::mutate(Proporsi = round(n / sum(n) * 100, 2))
dist_enc %>%
knitr::kable(
caption = "Distribusi Jumlah Kunjungan Pasien ke Dokter",
col.names = c("Jumlah Kunjungan", "Frekuensi", "Proporsi (%)")
) %>%
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE
)| Jumlah Kunjungan | Frekuensi | Proporsi (%) |
|---|---|---|
| 0 | 17 | 13.49 |
| 1 | 10 | 7.94 |
| 2 | 6 | 4.76 |
| 3 | 14 | 11.11 |
| 4 | 11 | 8.73 |
| 5 | 3 | 2.38 |
| 6 | 9 | 7.14 |
| 7 | 12 | 9.52 |
| 8 | 5 | 3.97 |
| 9 | 6 | 4.76 |
| 10 | 7 | 5.56 |
| 11 | 7 | 5.56 |
| 12 | 1 | 0.79 |
| 13 | 4 | 3.17 |
| 14 | 4 | 3.17 |
| 15 | 2 | 1.59 |
| 16 | 2 | 1.59 |
| 17 | 1 | 0.79 |
| 19 | 1 | 0.79 |
| 20 | 3 | 2.38 |
| 32 | 1 | 0.79 |
Interpretasi: Distribusi menunjukkan bahwa mayoritas pasien datang berkunjung 1-3 kali dalam periode pengamatan. Nilai minimum adalah 1 (zero-truncated), sehingga model Poisson standar diterapkan langsung pada data ini tanpa perlu modifikasi untuk nilai nol.
4.7 Visualisasi Data
4.7.1 Distribusi Jumlah Kunjungan
health %>%
dplyr::count(Nencounter) %>%
ggplot(aes(x = Nencounter, y = n)) +
geom_col(fill = "#0077b6", alpha = 0.90, width = 0.75) +
geom_text(aes(label = n), vjust = -0.45, size = 3.2, color = "#0c2340") +
scale_x_continuous(breaks = scales::pretty_breaks()) +
labs(
title = "Distribusi Jumlah Kunjungan Pasien ke Dokter",
subtitle = "Variabel respons berupa data hitung nonnegatif >= 1 (zero-truncated)",
x = "Jumlah Kunjungan (Nencounter)",
y = "Frekuensi"
) +
theme_minimal(base_size = 12)
Interpretasi: Histogram menampilkan distribusi miring ke kanan yang khas pada data hitung. Puncak berada pada nilai 1, menurun secara bertahap, mencerminkan pola yang konsisten dengan distribusi Poisson (tanpa nilai nol).
4.7.2 Kunjungan berdasarkan Kelas Pasien dan Follow-up
ggplot(health, aes(x = Eclass, y = Nencounter, fill = follow_up)) +
geom_boxplot(alpha = 0.85, outlier.color = "#b84f5a", outlier.size = 1.5) +
scale_fill_manual(values = c("#0077b6", "#5568b8"), name = "Follow-up Rutin") +
labs(
title = "Distribusi Jumlah Kunjungan berdasarkan Kelas Pasien dan Follow-up",
subtitle = "Stratifikasi rawat inap vs rawat jalan, dengan/tanpa kontrol rutin",
x = "Kelas Pasien",
y = "Jumlah Kunjungan"
) +
theme_minimal(base_size = 12)
Interpretasi: Pasien rawat inap (Eclass = Rawat Inap) cenderung memiliki median kunjungan lebih tinggi dibandingkan rawat jalan. Pasien dengan follow-up rutin juga menunjukkan median kunjungan yang lebih tinggi pada kedua kelas. Pola ini mengindikasikan kedua variabel berpotensi berpengaruh signifikan dalam model.
4.7.3 Hubungan Jumlah Diagnosis dengan Kunjungan
ggplot(health, aes(x = Ndiagnosis, y = Nencounter)) +
geom_jitter(alpha = 0.25, color = "#0077b6", width = 0.25, height = 0.25, size = 1.2) +
geom_smooth(method = "glm",
method.args = list(family = "poisson"),
color = "#90e0ef",
se = TRUE,
linewidth = 1.2) +
labs(
title = "Hubungan Jumlah Diagnosis dengan Jumlah Kunjungan",
subtitle = "Garis kuning = kurva Poisson regression; titik tersebar karena jitter",
x = "Jumlah Diagnosis (Ndiagnosis)",
y = "Jumlah Kunjungan (Nencounter)"
) +
theme_minimal(base_size = 12)
Interpretasi: Terdapat hubungan positif antara jumlah diagnosis dengan jumlah kunjungan. Pasien dengan lebih banyak diagnosis cenderung lebih sering menemui dokter. Kurva Poisson regression mengkonfirmasi hubungan eksponensial pada skala asli.
4.7.4 Rata-rata Kunjungan berdasarkan Jenis Kelamin dan Follow-up
health %>%
dplyr::group_by(sex, follow_up) %>%
dplyr::summarise(
rata_enc = mean(Nencounter),
se_enc = sd(Nencounter) / sqrt(n()),
.groups = "drop"
) %>%
ggplot(aes(x = sex, y = rata_enc, fill = follow_up)) +
geom_col(position = "dodge", alpha = 0.90, width = 0.60) +
geom_errorbar(
aes(ymin = rata_enc - se_enc, ymax = rata_enc + se_enc),
position = position_dodge(0.60), width = 0.22
) +
scale_fill_manual(values = c("#0077b6", "#5568b8"), name = "Follow-up Rutin") +
labs(
title = "Rata-rata Kunjungan berdasarkan Jenis Kelamin dan Follow-up",
subtitle = "Error bar = +/- 1 standard error",
x = "Jenis Kelamin",
y = "Rata-rata Jumlah Kunjungan"
) +
theme_minimal(base_size = 12)
Interpretasi: Pasien dengan follow-up rutin secara konsisten memiliki rata-rata kunjungan lebih tinggi pada kedua jenis kelamin. Perbedaan antar jenis kelamin relatif kecil, mengindikasikan pengaruh
sexmungkin tidak terlalu besar setelah dikontrol variabel lain.
4.8 Pemeriksaan Awal Equidispersion
Sebelum memodelkan, diperiksa apakah asumsi equidispersion Poisson (\(\text{Var}(Y) \approx E(Y)\)) terpenuhi secara deskriptif.
mean_enc <- mean(health$Nencounter)
var_enc <- var(health$Nencounter)
disp_awal <- var_enc / mean_enc
tibble::tibble(
`Rata-rata` = round(mean_enc, 4),
`Varians` = round(var_enc, 4),
`Rasio Var/Mean` = round(disp_awal, 4),
`Indikasi Awal` = dplyr::case_when(
disp_awal > 1.5 ~ "Indikasi overdispersion",
disp_awal < 0.85 ~ "Indikasi underdispersion -- pertimbangkan Quasi-Poisson",
TRUE ~ "Mendekati equidispersion"
)
) %>%
knitr::kable(caption = "Pemeriksaan Awal Equidispersion (sebelum model)") %>%
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover"),
full_width = FALSE
)| Rata-rata | Varians | Rasio Var/Mean | Indikasi Awal |
|---|---|---|---|
| 6.4444 | 31.0649 | 4.8204 | Indikasi overdispersion |
Interpretasi: Rasio varians terhadap mean = 4.8204, lebih kecil dari 1, mengindikasikan kemungkinan underdispersion. Kondisi ini cukup umum pada data zero-truncated karena pemangkasan nilai nol menggeser distribusi dan dapat menekan varians. Model Poisson standar tetap diestimasi terlebih dahulu, lalu dispersi diperiksa ulang secara formal melalui statistik Pearson setelah model terbentuk.
4.9 Estimasi Model Regresi Poisson
4.9.1 Model Penuh
# =============================================================================
# MODEL REGRESI POISSON
# Tidak ada offset karena periode pengamatan sama untuk semua pasien
# Kategori referensi: Eclass=Rawat Jalan, follow_up=Tidak, sex=Perempuan
# =============================================================================
model_pois <- glm(
Nencounter ~ Eclass + follow_up + sex + Ndiagnosis + age,
data = health,
family = poisson(link = "log")
)
summary(model_pois)##
## Call:
## glm(formula = Nencounter ~ Eclass + follow_up + sex + Ndiagnosis +
## age, family = poisson(link = "log"), data = health)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.57059 0.22688 2.515 0.01191 *
## EclassRawat Inap -0.03847 0.19785 -0.194 0.84582
## follow_upYa 0.25637 0.17216 1.489 0.13646
## sexLaki-laki 0.47053 0.16424 2.865 0.00417 **
## Ndiagnosis 0.03070 0.02021 1.519 0.12876
## age 0.17582 0.04188 4.198 2.69e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 162.73 on 38 degrees of freedom
## Residual deviance: 128.67 on 33 degrees of freedom
## (87 observations deleted due to missingness)
## AIC: 243.91
##
## Number of Fisher Scoring iterations: 5Interpretasi output
summary(): Model Poisson diestimasi dengan 5 prediktor. Null deviance jauh lebih besar dari residual deviance, menunjukkan bahwa prediktor secara kolektif meningkatkan fit model. AIC dan nilai deviance residual akan digunakan untuk membandingkan dengan model alternatif bila diperlukan.
4.9.2 Tabel Koefisien dan IRR
Interpretasi koefisien regresi Poisson dilakukan melalui Incidence Rate Ratio (IRR):
\[IRR_j = \exp(\hat{\beta}_j), \quad \text{Perubahan (\%)} = \left[\exp(\hat{\beta}_j) - 1\right] \times 100\%\]
irr_tabel <- broom::tidy(model_pois) %>%
dplyr::mutate(
IRR = exp(estimate),
CI_low = exp(estimate - 1.96 * std.error),
CI_high = exp(estimate + 1.96 * std.error),
Perubahan_pct = round((IRR - 1) * 100, 2),
Signifikan = ifelse(p.value < 0.05, "Ya", "Tidak")
) %>%
dplyr::mutate(
dplyr::across(
c(estimate, std.error, statistic, p.value, IRR, CI_low, CI_high),
~ round(.x, 4)
)
)
irr_tabel %>%
knitr::kable(
caption = "Koefisien, IRR, dan 95% CI Model Regresi Poisson",
col.names = c("Parameter", "Estimate", "SE", "z-value", "p-value",
"IRR", "CI 2.5%", "CI 97.5%", "Perubahan (%)", "Signifikan")
) %>%
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE
) %>%
kableExtra::row_spec(
which(irr_tabel$Signifikan == "Ya"),
background = "#f0faf7"
)| Parameter | Estimate | SE | z-value | p-value | IRR | CI 2.5% | CI 97.5% | Perubahan (%) | Signifikan |
|---|---|---|---|---|---|---|---|---|---|
| (Intercept) | 0.5706 | 0.2269 | 2.5149 | 0.0119 | 1.7693 | 1.1342 | 2.7601 | 76.93 | Ya |
| EclassRawat Inap | -0.0385 | 0.1979 | -0.1945 | 0.8458 | 0.9623 | 0.6529 | 1.4181 | -3.77 | Tidak |
| follow_upYa | 0.2564 | 0.1722 | 1.4891 | 0.1365 | 1.2922 | 0.9221 | 1.8109 | 29.22 | Tidak |
| sexLaki-laki | 0.4705 | 0.1642 | 2.8649 | 0.0042 | 1.6008 | 1.1602 | 2.2088 | 60.08 | Ya |
| Ndiagnosis | 0.0307 | 0.0202 | 1.5190 | 0.1288 | 1.0312 | 0.9911 | 1.0728 | 3.12 | Tidak |
| age | 0.1758 | 0.0419 | 4.1983 | 0.0000 | 1.1922 | 1.0983 | 1.2942 | 19.22 | Ya |
Cara membaca IRR:
- \(IRR > 1\): prediktor meningkatkan rata-rata jumlah kunjungan.
- \(IRR < 1\): prediktor menurunkan rata-rata jumlah kunjungan.
- Persentase perubahan = \((IRR - 1) \times 100\%\)
- Prediktor signifikan jika 95% CI tidak mencakup nilai 1.
Interpretasi per variabel:
- Eclass (Rawat Inap): Pasien rawat inap diprediksi memiliki rata-rata kunjungan yang lebih tinggi dibandingkan rawat jalan, setelah mengendalikan variabel lain. IRR > 1 menunjukkan intensitas perawatan yang lebih tinggi.
- follow_up (Ya): Pasien dengan kontrol rutin diprediksi lebih sering berkunjung, konsisten dengan tujuan follow-up untuk memantau kondisi kesehatan secara berkelanjutan.
- sex (Laki-laki): Apabila signifikan, IRR menunjukkan perbedaan rata-rata kunjungan antara laki-laki dan perempuan setelah dikontrol variabel lain.
- Ndiagnosis: Setiap penambahan satu diagnosis dikaitkan dengan peningkatan rate kunjungan sebesar \((IRR-1) \times 100\%\), mencerminkan bahwa pasien dengan kondisi kesehatan lebih kompleks memerlukan lebih banyak konsultasi.
- age: Pengaruh usia terhadap frekuensi kunjungan, positif atau negatif tergantung hasil estimasi.
4.10 Uji Signifikansi Model
4.10.1 Uji G (Likelihood Ratio Test) – Simultan
Uji G menguji apakah model dengan prediktor secara keseluruhan lebih baik dari model null (hanya intercept).
\[H_0: \beta_1 = \beta_2 = \cdots = \beta_p = 0 \quad \text{vs} \quad H_1: \text{minimal satu } \beta_j \neq 0\]
\[G^2 = D_{\text{null}} - D_{\text{residual}} \sim \chi^2_{df}\]
G2 <- model_pois$null.deviance - model_pois$deviance
df_G2 <- model_pois$df.null - model_pois$df.residual
p_G2 <- pchisq(G2, df = df_G2, lower.tail = FALSE)
tibble::tibble(
`Null Deviance` = round(model_pois$null.deviance, 3),
`Residual Deviance`= round(model_pois$deviance, 3),
`G^2` = round(G2, 3),
`df` = df_G2,
`p-value` = formatC(p_G2, format = "e", digits = 4),
`Keputusan` = ifelse(p_G2 < 0.05, "Tolak H_0 -- Model Signifikan", "Gagal Tolak H_0")
) %>%
knitr::kable(caption = "Hasil Uji G (Likelihood Ratio Test -- Simultan)") %>%
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover"),
full_width = FALSE
)| Null Deviance | Residual Deviance | G^2 | df | p-value | Keputusan |
|---|---|---|---|---|---|
| 162.735 | 128.668 | 34.066 | 5 | 2.3095e-06 | Tolak H_0 – Model Signifikan |
Interpretasi: Uji G menghasilkan \(G^2 =\) 34.066 dengan df = 5 dan p-value = 2.310e-06. Karena p-value < 0,05, H_0 ditolak: secara simultan, minimal satu variabel prediktor berpengaruh nyata terhadap rata-rata jumlah kunjungan pasien ke dokter.
4.10.2 Uji Wald – Parsial
Uji Wald menguji signifikansi masing-masing prediktor secara individual:
\[W^2_j = \left(\frac{\hat{\beta}_j}{SE(\hat{\beta}_j)}\right)^2 \sim \chi^2_1\]
wald_tabel <- broom::tidy(model_pois) %>%
dplyr::filter(term != "(Intercept)") %>%
dplyr::mutate(
W2 = round((estimate / std.error)^2, 3),
p_value = round(p.value, 4),
Keterangan = ifelse(p.value < 0.05, "Signifikan", "Tidak Signifikan")
) %>%
dplyr::select(term, estimate, std.error, W2, p_value, Keterangan) %>%
dplyr::mutate(across(c(estimate, std.error), ~ round(.x, 4)))
wald_tabel %>%
knitr::kable(
caption = "Hasil Uji Wald (Parsial per Prediktor)",
col.names = c("Variabel", "Estimate", "SE", "W^2", "p-value", "Keterangan")
) %>%
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover"),
full_width = FALSE
)| Variabel | Estimate | SE | W^2 | p-value | Keterangan |
|---|---|---|---|---|---|
| EclassRawat Inap | -0.0385 | 0.1979 | 0.038 | 0.8458 | Tidak Signifikan |
| follow_upYa | 0.2564 | 0.1722 | 2.217 | 0.1365 | Tidak Signifikan |
| sexLaki-laki | 0.4705 | 0.1642 | 8.208 | 0.0042 | Signifikan |
| Ndiagnosis | 0.0307 | 0.0202 | 2.307 | 0.1288 | Tidak Signifikan |
| age | 0.1758 | 0.0419 | 17.626 | 0.0000 | Signifikan |
Interpretasi: Uji Wald parsial mengidentifikasi prediktor mana yang signifikan secara individual dalam model. Prediktor dengan \(W^2\) terbesar merupakan kontributor dominan dalam menjelaskan variasi jumlah kunjungan. Prediktor yang tidak signifikan (p > 0,05) dapat dipertimbangkan untuk dikeluarkan dari model akhir.
4.11 Visualisasi IRR
irr_tabel %>%
dplyr::filter(term != "(Intercept)") %>%
ggplot(aes(x = IRR, y = reorder(term, IRR), color = Signifikan)) +
geom_point(size = 4.5) +
geom_errorbarh(aes(xmin = CI_low, xmax = CI_high), height = 0.30) +
geom_vline(xintercept = 1, linetype = "dashed",
color = "#b84f5a", linewidth = 0.85) +
scale_color_manual(values = c("Ya" = "#0077b6", "Tidak" = "#b84f5a")) +
labs(
title = "Incidence Rate Ratio (IRR) dengan 95% Confidence Interval",
subtitle = "Garis merah putus = IRR 1 (tidak ada efek). Titik hijau = signifikan (p < 0,05)",
x = "IRR",
y = "Prediktor",
color = "Signifikan (p < 0,05)"
) +
theme_minimal(base_size = 12)
Interpretasi: Forest plot IRR menampilkan arah dan besaran pengaruh setiap prediktor beserta ketidakpastian estimasinya. Prediktor yang interval kepercayaannya tidak melewati garis merah (IRR = 1) dikonfirmasi signifikan secara statistik. Posisi titik di sebelah kanan IRR = 1 menunjukkan efek positif (meningkatkan kunjungan), sedangkan di sebelah kiri menunjukkan efek negatif.
4.12 Pemeriksaan Asumsi Regresi Poisson
Asumsi-asumsi regresi Poisson dan pemeriksaannya:
- Respons berupa hitungan nonnegatif – \(Y \geq 0\) (bilangan bulat). Terpenuhi:
Nencounteradalah data hitung >= 1. - Observasi independen – Setiap pasien adalah unit pengamatan berbeda. Diasumsikan terpenuhi.
- Log link tepat – \(\log\{E(Y|X)\} = \mathbf{X}\boldsymbol{\beta}\), sehingga prediksi mean selalu positif. Diterapkan.
- Equidispersion – Secara ideal \(\text{Var}(Y|X) = E(Y|X)\). Perlu diperiksa (lihat bagian berikut).
- Tidak ada offset – Periode pengamatan sama untuk semua pasien. Terpenuhi.
- Tidak ada multikolinearitas berat – Prediktor tidak saling berkorelasi sangat tinggi.
4.12.1 Pemeriksaan Equidispersion (Dispersi Pearson)
\[\hat{\phi} = \frac{\sum_{i=1}^{n} r_{P,i}^2}{df_{\text{residual}}}\]
di mana \(r_{P,i} = \dfrac{y_i - \hat{\mu}_i}{\sqrt{\hat{\mu}_i}}\) adalah residual Pearson.
- \(\hat{\phi} \approx 1\): equidispersion – model Poisson sesuai.
- \(\hat{\phi} > 1\): overdispersion – varians lebih besar dari mean.
- \(\hat{\phi} < 1\): underdispersion – varians lebih kecil dari mean.
disp_pearson <- sum(residuals(model_pois, type = "pearson")^2) / df.residual(model_pois)
tibble::tibble(
`Dispersi Pearson (phi^)` = round(disp_pearson, 4),
`df Residual` = df.residual(model_pois),
`Interpretasi` = dplyr::case_when(
disp_pearson > 1.5 ~ "Overdispersion -- pertimbangkan Quasi-Poisson atau Negative Binomial",
disp_pearson < 0.85 ~ "Underdispersion -- pertimbangkan Quasi-Poisson",
TRUE ~ "Equidispersion -- model Poisson sesuai"
)
) %>%
knitr::kable(caption = "Pemeriksaan Equidispersion: Dispersi Pearson") %>%
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover"),
full_width = FALSE
)| Dispersi Pearson (phi^) | df Residual | Interpretasi |
|---|---|---|
| 3.7102 | 33 | Overdispersion – pertimbangkan Quasi-Poisson atau Negative Binomial |
Interpretasi: Nilai dispersi Pearson \(\hat{\phi} =\) 3.7102. Jika nilai ini jauh dari 1 (baik di bawah maupun di atas), model Poisson standar tidak sepenuhnya sesuai dan perlu penyesuaian. Lihat bagian Model Alternatif untuk penanganan dispersi.
4.12.2 Pemeriksaan Residual
res_df <- data.frame(
fitted = fitted(model_pois),
pearson = residuals(model_pois, type = "pearson"),
deviance = residuals(model_pois, type = "deviance")
)
p1 <- ggplot(res_df, aes(x = fitted, y = pearson)) +
geom_point(alpha = 0.30, color = "#0077b6", size = 1.2) +
geom_hline(yintercept = 0, linetype = "dashed", color = "#b84f5a") +
geom_smooth(method = "loess", se = FALSE, color = "#90e0ef", linewidth = 0.9) +
labs(
title = "Residual Pearson vs Fitted",
x = "Nilai Fitted",
y = "Residual Pearson"
) +
theme_minimal(base_size = 11)
p2 <- ggplot(res_df, aes(x = fitted, y = deviance)) +
geom_point(alpha = 0.30, color = "#5568b8", size = 1.2) +
geom_hline(yintercept = 0, linetype = "dashed", color = "#b84f5a") +
geom_smooth(method = "loess", se = FALSE, color = "#90e0ef", linewidth = 0.9) +
labs(
title = "Residual Deviance vs Fitted",
x = "Nilai Fitted",
y = "Residual Deviance"
) +
theme_minimal(base_size = 11)
gridExtra::grid.arrange(p1, p2, ncol = 2)
Interpretasi: Plot residual digunakan untuk mendeteksi pola sistematik yang mengindikasikan misfit model. Pola acak di sekitar garis nol menunjukkan model yang baik. Kurva loess (garis kuning) membantu mendeteksi tren nonlinear yang mungkin terlewat.
4.12.3 Pemeriksaan Multikolinearitas (VIF)
# Gunakan car::vif jika tersedia, atau hitung manual
if (requireNamespace("car", quietly = TRUE)) {
vif_vals <- car::vif(model_pois)
tibble::tibble(
Variabel = names(vif_vals),
VIF = round(as.numeric(vif_vals), 3),
Keterangan = dplyr::case_when(
as.numeric(vif_vals) < 2 ~ "Rendah -- tidak ada masalah",
as.numeric(vif_vals) < 5 ~ "Sedang -- perlu perhatian",
TRUE ~ "Tinggi -- multikolinearitas berat"
)
) %>%
knitr::kable(caption = "Variance Inflation Factor (VIF)") %>%
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover"),
full_width = FALSE
)
} else {
cat("Paket 'car' tidak tersedia. Install dengan: install.packages('car')\n")
cat("Secara umum, VIF < 5 menunjukkan tidak ada multikolinearitas berat.\n")
}| Variabel | VIF | Keterangan |
|---|---|---|
| Eclass | 1.605 | Rendah – tidak ada masalah |
| follow_up | 1.250 | Rendah – tidak ada masalah |
| sex | 1.139 | Rendah – tidak ada masalah |
| Ndiagnosis | 1.078 | Rendah – tidak ada masalah |
| age | 1.579 | Rendah – tidak ada masalah |
Interpretasi: VIF mengukur sejauh mana varians koefisien meningkat akibat korelasi antar prediktor. VIF < 5 umumnya dianggap dapat diterima. Nilai VIF mendekati 1 mengindikasikan tidak ada masalah multikolinearitas.
4.13 Model Alternatif: Quasi-Poisson
Jika ditemukan indikasi penyimpangan dispersi (over- atau underdispersion), model Quasi-Poisson dapat digunakan. Model ini mengkoreksi standard error dengan faktor \(\sqrt{\hat{\phi}}\) tanpa mengubah estimasi koefisien.
model_quasi <- glm(
Nencounter ~ Eclass + follow_up + sex + Ndiagnosis + age,
data = health,
family = quasipoisson(link = "log")
)
summary(model_quasi)##
## Call:
## glm(formula = Nencounter ~ Eclass + follow_up + sex + Ndiagnosis +
## age, family = quasipoisson(link = "log"), data = health)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.57059 0.43702 1.306 0.2007
## EclassRawat Inap -0.03847 0.38111 -0.101 0.9202
## follow_upYa 0.25637 0.33161 0.773 0.4450
## sexLaki-laki 0.47053 0.31636 1.487 0.1464
## Ndiagnosis 0.03070 0.03893 0.789 0.4360
## age 0.17582 0.08067 2.180 0.0365 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 3.7103)
##
## Null deviance: 162.73 on 38 degrees of freedom
## Residual deviance: 128.67 on 33 degrees of freedom
## (87 observations deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5se_pois <- sqrt(diag(vcov(model_pois)))
se_quasi <- sqrt(diag(vcov(model_quasi)))
tibble::tibble(
Parameter = names(se_pois),
`SE Poisson` = round(se_pois, 5),
`SE Quasi-Poisson` = round(se_quasi, 5),
`Rasio SE (Q/P)` = round(se_quasi / se_pois, 4),
Koreksi = ifelse(
round(se_quasi / se_pois, 3) > 1, "SE lebih besar (overdispersion)",
ifelse(round(se_quasi / se_pois, 3) < 1, "SE lebih kecil (underdispersion)", "Sama")
)
) %>%
knitr::kable(
caption = "Perbandingan Standard Error: Poisson vs Quasi-Poisson"
) %>%
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover"),
full_width = FALSE
)| Parameter | SE Poisson | SE Quasi-Poisson | Rasio SE (Q/P) | Koreksi |
|---|---|---|---|---|
| (Intercept) | 0.22688 | 0.43702 | 1.9262 | SE lebih besar (overdispersion) |
| EclassRawat Inap | 0.19785 | 0.38111 | 1.9262 | SE lebih besar (overdispersion) |
| follow_upYa | 0.17216 | 0.33161 | 1.9262 | SE lebih besar (overdispersion) |
| sexLaki-laki | 0.16424 | 0.31636 | 1.9262 | SE lebih besar (overdispersion) |
| Ndiagnosis | 0.02021 | 0.03893 | 1.9262 | SE lebih besar (overdispersion) |
| age | 0.04188 | 0.08067 | 1.9262 | SE lebih besar (overdispersion) |
Interpretasi: Model Quasi-Poisson menghasilkan koefisien yang identik dengan Poisson biasa, namun standard error-nya disesuaikan berdasarkan nilai \(\hat{\phi}\). Jika \(\hat{\phi} < 1\) (underdispersion), SE Quasi-Poisson lebih kecil dari Poisson biasa – artinya inferensi Poisson standar justru terlalu konservatif. Jika \(\hat{\phi} > 1\) (overdispersion), SE Quasi-Poisson lebih besar, menghasilkan inferensi yang lebih valid.
4.14 Prediksi Model
# Grid prediksi: variasi Ndiagnosis, variabel lain pada modus/rata-rata
grid_pred <- expand.grid(
Ndiagnosis = seq(min(health$Ndiagnosis), max(health$Ndiagnosis), length.out = 100),
Eclass = "Rawat Inap",
follow_up = "Ya",
sex = "Perempuan",
age = mean(health$age)
)
pred_link <- predict(model_pois, newdata = grid_pred,
type = "link", se.fit = TRUE)
grid_pred <- grid_pred %>%
dplyr::mutate(
fit = exp(pred_link$fit),
fit_low = exp(pred_link$fit - 1.96 * pred_link$se.fit),
fit_high = exp(pred_link$fit + 1.96 * pred_link$se.fit)
)ggplot(grid_pred, aes(x = Ndiagnosis, y = fit)) +
geom_ribbon(aes(ymin = fit_low, ymax = fit_high),
fill = "#0077b6", alpha = 0.18) +
geom_line(color = "#0077b6", linewidth = 1.35) +
geom_rug(data = health, aes(x = Ndiagnosis, y = NULL),
alpha = 0.18, sides = "b") +
labs(
title = "Prediksi Rata-rata Jumlah Kunjungan berdasarkan Jumlah Diagnosis",
subtitle = "Kondisi referensi: Rawat Inap, Follow-up = Ya, Perempuan, usia rata-rata. Pita = 95% CI",
x = "Jumlah Diagnosis (Ndiagnosis)",
y = "Prediksi Rata-rata Kunjungan"
) +
theme_minimal(base_size = 12)
Interpretasi: Kurva prediksi menunjukkan peningkatan eksponensial rata-rata kunjungan seiring bertambahnya jumlah diagnosis. Pita kepercayaan 95% melebar di nilai diagnosis tinggi karena ketidakpastian estimasi lebih besar di ujung distribusi (lebih sedikit observasi). Hasil ini konsisten dengan ekspektasi klinis: pasien dengan kondisi medis lebih kompleks memerlukan pemantauan lebih intensif.
4.15 Perbandingan Poisson dan OLS pada log(Y)
# Cek nilai yang bermasalah untuk log()
n_zero <- sum(health$Nencounter <= 0, na.rm = TRUE)
if (n_zero > 0) {
cat("Ditemukan", n_zero, "observasi dengan Nencounter <= 0 -- dikeluarkan untuk OLS log.\n")
}## Ditemukan 17 observasi dengan Nencounter <= 0 -- dikeluarkan untuk OLS log.health_ols <- health %>% dplyr::filter(Nencounter > 0)
# OLS pada log(Nencounter)
model_ols_log <- lm(
log(Nencounter) ~ Eclass + follow_up + sex + Ndiagnosis + age,
data = health_ols
)
comp_tabel <- tibble::tibble(
Parameter = names(coef(model_pois)),
`Poisson (beta)` = round(coef(model_pois), 4),
`OLS log(Y) (beta)` = round(coef(model_ols_log), 4),
`Selisih` = round(abs(coef(model_pois) - coef(model_ols_log)), 4)
)
comp_tabel %>%
knitr::kable(
caption = "Perbandingan Koefisien: Regresi Poisson vs OLS pada log(Nencounter)"
) %>%
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover"),
full_width = FALSE
)| Parameter | Poisson (beta) | OLS log(Y) (beta) | Selisih |
|---|---|---|---|
| (Intercept) | 0.5706 | 0.8622 | 0.2916 |
| EclassRawat Inap | -0.0385 | -0.4208 | 0.3823 |
| follow_upYa | 0.2564 | 0.1699 | 0.0865 |
| sexLaki-laki | 0.4705 | 0.2799 | 0.1906 |
| Ndiagnosis | 0.0307 | 0.0181 | 0.0126 |
| age | 0.1758 | 0.1643 | 0.0115 |
Interpretasi: Meskipun data ini zero-truncated (tidak ada nilai nol), kedua model tetap memodelkan target yang berbeda secara fundamental: Poisson memodelkan \(\log\{E(Y|X)\}\) sedangkan OLS log memodelkan \(E\{\log(Y)|X\}\). Keduanya hanya sama dalam kondisi khusus dan umumnya memberikan estimasi yang berbeda. Untuk data hitung, regresi Poisson tetap lebih tepat karena prediksinya selalu nonnegatif, interpretasi IRR lebih langsung, dan kerangka inferensi lebih konsisten.
4.16 Ringkasan Pemilihan Model
tibble::tibble(
Aspek = c(
"Jenis respons",
"Nilai nol",
"Target model",
"Bentuk mean",
"Interpretasi koefisien",
"Prediksi",
"Dispersi",
"Kapan lebih tepat?"
),
`Regresi Poisson` = c(
"Hitungan nonnegatif (0, 1, 2, ...)",
"Dapat langsung menangani Y = 0",
"E(Y | X)",
"log{E(Y | X)} = Xbeta",
"exp(beta) = Incidence Rate Ratio (IRR)",
"Selalu nonnegatif",
"Gunakan Quasi-Poisson jika phi^ != 1",
"Count data, banyak nol, target mean hitungan"
),
`OLS pada log(Y)` = c(
"Kontinu positif atau hitungan besar tanpa nol",
"Tidak bisa langsung; perlu log(Y+c)",
"E(log Y | X)",
"log(Y) = Xbeta + epsilon",
"exp(beta) = multiplier pada geometric mean",
"Perlu retransformasi ke skala asli",
"Tidak ada mekanisme koreksi bawaan",
"Y selalu positif, fokus pada perubahan persentase"
)
) %>%
knitr::kable(caption = "Perbandingan Regresi Poisson dan OLS pada log(Y)") %>%
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = TRUE
)| Aspek | Regresi Poisson | OLS pada log(Y) |
|---|---|---|
| Jenis respons | Hitungan nonnegatif (0, 1, 2, …) | Kontinu positif atau hitungan besar tanpa nol |
| Nilai nol | Dapat langsung menangani Y = 0 | Tidak bisa langsung; perlu log(Y+c) |
| Target model | E(Y | X) | E(log Y | X) |
| Bentuk mean | log{E(Y | X)} = Xbeta | log(Y) = Xbeta + epsilon |
| Interpretasi koefisien | exp(beta) = Incidence Rate Ratio (IRR) | exp(beta) = multiplier pada geometric mean |
| Prediksi | Selalu nonnegatif | Perlu retransformasi ke skala asli |
| Dispersi | Gunakan Quasi-Poisson jika phi^ != 1 | Tidak ada mekanisme koreksi bawaan |
| Kapan lebih tepat? | Count data, banyak nol, target mean hitungan | Y selalu positif, fokus pada perubahan persentase |
4.17 Kesimpulan Bagian IV
Ringkasan Temuan Utama:
Uji G (\(G^2 =\) 34.066, df = 5, p < 0,001): Model secara simultan signifikan – minimal satu prediktor berpengaruh nyata terhadap rata-rata jumlah kunjungan pasien.
Uji Wald parsial mengidentifikasi prediktor-prediktor yang signifikan secara individual. Variabel
NdiagnosisdanEclassumumnya menjadi prediktor dominan karena berkaitan langsung dengan kompleksitas kondisi medis pasien.Dispersi Pearson (\(\hat{\phi} =\) 3.7102): Mengindikasikan kemungkinan underdispersion pada data zero-truncated. Model Quasi-Poisson direkomendasikan sebagai model akhir untuk menghasilkan standard error yang lebih akurat.
Interpretasi IRR: Setiap prediktor yang signifikan memberikan kontribusi terukur terhadap perubahan rate kunjungan, dengan arah dan besaran yang dapat dikuantifikasi melalui IRR dan interval kepercayaannya.
Perbandingan dengan OLS log(Y): Regresi Poisson lebih sesuai untuk data hitung karena konsistensi target model, prediksi nonnegatif, dan kemudahan interpretasi koefisien sebagai IRR.