Tugas Analisis Data Kategori

Regresi Logistik Biner, Multinomial, Ordinal, dan Regresi Poisson

1 Regresi Logistik Biner

1.1 Pendahuluan

Regresi logistik biner digunakan ketika variabel dependen hanya memiliki dua kategori. Model ini tidak memprediksi nilai Y secara langsung, melainkan memperkirakan probabilitas suatu observasi masuk ke salah satu kategori, menggunakan fungsi glm() dengan argumen family = binomial(link = "logit").

Tujuan analisis ini adalah memodelkan tipe kepribadian seseorang, Introvert atau Extrovert, berdasarkan karakteristik perilaku sosialnya.

Variabel yang digunakan:

Tipe	Variabel	Keterangan
Dependen (Y)	Personality	Introvert / Extrovert
Independen Numerik	Time_spent_Alone	Waktu yang dihabiskan sendirian (jam/hari)
Independen Numerik	Social_event_attendance	Frekuensi kehadiran di acara sosial
Independen Numerik	Friends_circle_size	Ukuran lingkaran pertemanan
Independen Numerik	Post_frequency	Frekuensi posting di media sosial
Independen Kategorik	Drained_after_socializing	Merasa lelah setelah bersosialisasi (Yes/No)

1.2 Memuat Paket

library(readxl)
library(tidyverse)
library(scales)
library(knitr)
library(kableExtra)
library(gtsummary)
library(car)
library(lmtest)
library(ResourceSelection)
library(pROC)
library(caret)

Paket utama yang digunakan adalah glm() dari base R untuk estimasi model, gtsummary untuk tabel deskriptif, car untuk uji VIF, lmtest untuk Likelihood Ratio Test, ResourceSelection untuk uji Hosmer-Lemeshow, dan pROC untuk kurva ROC.

1.3 Memuat dan Memeriksa Data

df_b <- readxl::read_excel("Data Biner.xlsx")
df_b$Personality <- as.factor(df_b$Personality)

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

## Jumlah baris   : 2900

cat("Jumlah kolom   :", ncol(df_b), "\n")

## Jumlah kolom   : 8

head(df_b, 6)

str(df_b)

## tibble [2,900 × 8] (S3: tbl_df/tbl/data.frame)
##  $ Time_spent_Alone         : num [1:2900] 4 9 9 0 3 1 4 2 10 0 ...
##  $ Stage_fear               : chr [1:2900] "No" "Yes" "Yes" "No" ...
##  $ Social_event_attendance  : num [1:2900] 4 0 1 6 9 7 9 8 1 8 ...
##  $ Going_outside            : num [1:2900] 6 0 2 7 4 5 3 4 3 6 ...
##  $ Drained_after_socializing: chr [1:2900] "No" "Yes" "Yes" "No" ...
##  $ Friends_circle_size      : num [1:2900] 13 0 5 14 8 6 7 7 0 13 ...
##  $ Post_frequency           : num [1:2900] 5 3 2 8 5 6 7 8 3 8 ...
##  $ Personality              : Factor w/ 2 levels "Extrovert","Introvert": 1 2 2 1 1 1 1 1 2 1 ...

Variabel Personality memiliki 2 level yaitu Extrovert dan Introvert. Variabel numerik bertipe num, sementara Drained_after_socializing masih bertipe chr dan akan diubah pada tahap pra-pemrosesan.

1.4 Pra-Pemrosesan Data

Langkah pra-pemrosesan yang dilakukan meliputi memilih variabel yang relevan, mengatur kategori referensi pada tiap variabel kategorik, dan menghapus baris yang memiliki missing value menggunakan drop_na().

Kategori referensi ditetapkan sebagai berikut: Personality = “Extrovert” sehingga model memprediksi log-odds menjadi Introvert, dan Drained_after_socializing = “No”.

df_model_b <- df_b %>%
  select(Personality, Time_spent_Alone, Social_event_attendance,
         Friends_circle_size, Post_frequency, Drained_after_socializing) %>%
  mutate(
    Personality               = relevel(factor(Personality), ref = "Extrovert"),
    Drained_after_socializing = relevel(factor(Drained_after_socializing), ref = "No")
  ) %>%
  drop_na()

cat("\nJumlah observasi setelah drop_na():", nrow(df_model_b), "\n")

## 
## Jumlah observasi setelah drop_na(): 2900

cat("Distribusi Y (Personality):\n")

## Distribusi Y (Personality):

print(table(df_model_b$Personality))

## 
## Extrovert Introvert 
##      1491      1409

Setelah drop_na(), data bersih dari missing value. Semua koefisien model akan diinterpretasikan sebagai log-odds menjadi Introvert relatif terhadap Extrovert.

1.5 Analisis Deskriptif

1.5.1 Distribusi Variabel Dependen

df_model_b %>%
  count(Personality) %>%
  mutate(proporsi = n / sum(n)) %>%
  kable(digits = 3, caption = "Distribusi tipe kepribadian") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE)

Distribusi tipe kepribadian

Personality	n	proporsi
Extrovert	1491	0.514
Introvert	1409	0.486

df_model_b %>%
  count(Personality) %>%
  mutate(proporsi = n / sum(n)) %>%
  ggplot(aes(x = Personality, y = proporsi, fill = Personality)) +
  geom_col(width = 0.55, alpha = 0.93) +
  geom_text(aes(label = percent(proporsi, accuracy = 0.1)),
            vjust = -0.4, fontface = "bold", color = "#243142") +
  scale_fill_manual(values = c("Extrovert" = "#2f7f73", "Introvert" = "#5568b8")) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 0.65)) +
  labs(title    = "Distribusi Tipe Kepribadian",
       subtitle = paste0("Total observasi: ", nrow(df_model_b)),
       x = "Personality", y = "Proporsi") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none",
        plot.title = element_text(face = "bold"))

Distribusi dua kelas relatif seimbang. Keseimbangan ini menguntungkan karena model tidak akan condong ke satu kelas saja.

1.5.2 Distribusi Variabel Kategorik vs Kepribadian

for (v in c("Drained_after_socializing")) {
  p <- ggplot(df_model_b, aes(x = .data[[v]], fill = Personality)) +
    geom_bar(position = "fill", color = "white") +
    scale_fill_manual(values = c("Extrovert" = "#2f7f73", "Introvert" = "#5568b8")) +
    scale_y_continuous(labels = percent_format()) +
    labs(title = paste(v, "vs Personality"), x = v, y = "Proporsi") +
    theme_minimal(base_size = 11) +
    theme(plot.title = element_text(face = "bold"))
  print(p)
}

Jika komposisi warna berbeda nyata antar kategori, berarti variabel tersebut berhubungan dengan tipe kepribadian.

1.5.3 Distribusi Variabel Numerik vs Kepribadian

for (v in c("Time_spent_Alone", "Social_event_attendance",
            "Friends_circle_size", "Post_frequency")) {
  p <- ggplot(df_model_b, aes(x = Personality, y = .data[[v]], fill = Personality)) +
    geom_boxplot(outlier.alpha = 0.3, outlier.size = 1) +
    scale_fill_manual(values = c("Extrovert" = "#2f7f73", "Introvert" = "#5568b8")) +
    labs(title = paste(v, "vs Personality"), x = "Personality", y = v) +
    theme_minimal(base_size = 11) +
    theme(legend.position = "none", plot.title = element_text(face = "bold"))
  print(p)
}

Perbedaan median atau sebaran antar tipe kepribadian menunjukkan bahwa variabel numerik tersebut berperan dalam membedakan Introvert dan Extrovert.

1.6 Pemeriksaan Asumsi

1.6.1 Multikolinearitas (VIF)

model_vif_b <- lm(as.numeric(Personality) ~
                    Time_spent_Alone + Social_event_attendance +
                    Friends_circle_size + Post_frequency +
                    Drained_after_socializing,
                  data = df_model_b)

vif_vals_b <- vif(model_vif_b)

data.frame(
  Variabel = names(vif_vals_b),
  VIF      = round(as.numeric(vif_vals_b), 3),
  Status   = ifelse(as.numeric(vif_vals_b) < 5, "Aman (< 5)", "Perhatian (>= 5)")
) %>%
  kable(caption = "Variance Inflation Factor (VIF)") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Variance Inflation Factor (VIF)

Variabel	VIF	Status
Time_spent_Alone	3.326	Aman (< 5)
Social_event_attendance	3.369	Aman (< 5)
Friends_circle_size	2.940	Aman (< 5)
Post_frequency	3.199	Aman (< 5)
Drained_after_socializing	7.338	Perhatian (>= 5)

VIF < 5 untuk semua prediktor menunjukkan tidak ada multikolinearitas yang mengkhawatirkan.

Asumsi lain yang perlu dipenuhi dalam regresi logistik biner:

Asumsi	Status
Variabel respons bersifat biner	Terpenuhi: Introvert / Extrovert
Observasi independen	Diasumsikan terpenuhi
Tidak ada multikolinearitas berat (VIF < 5)	Terpenuhi (lihat tabel VIF)
Tidak ada perfect separation	Model berhasil diestimasi
Hubungan linear pada skala logit untuk prediktor kontinu	Perlu diperiksa lebih lanjut

1.7 Estimasi Model Regresi Logistik Biner

Model diestimasi menggunakan glm() dengan family = binomial(link = "logit"). Extrovert ditetapkan sebagai kategori referensi sehingga model menghasilkan satu set koefisien yang merepresentasikan log-odds menjadi Introvert.

fit_logit <- glm(
  Personality ~ Time_spent_Alone + Social_event_attendance +
    Friends_circle_size + Post_frequency + Drained_after_socializing,
  data   = df_model_b,
  family = binomial(link = "logit")
)

summary(fit_logit)

## 
## Call:
## glm(formula = Personality ~ Time_spent_Alone + Social_event_attendance + 
##     Friends_circle_size + Post_frequency + Drained_after_socializing, 
##     family = binomial(link = "logit"), data = df_model_b)
## 
## Coefficients:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  -0.978972   0.393859  -2.486 0.012934 *  
## Time_spent_Alone              0.019742   0.035282   0.560 0.575778    
## Social_event_attendance      -0.111795   0.044411  -2.517 0.011826 *  
## Friends_circle_size           0.001742   0.028058   0.062 0.950486    
## Post_frequency               -0.163337   0.044140  -3.700 0.000215 ***
## Drained_after_socializingYes  3.651520   0.305485  11.953  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 4017.9  on 2899  degrees of freedom
## Residual deviance: 1549.5  on 2894  degrees of freedom
## AIC: 1561.5
## 
## Number of Fisher Scoring iterations: 5

Output summary() memberikan koefisien log-odds, standard error, z-value, dan p-value untuk setiap prediktor. Null Deviance adalah devians model tanpa prediktor; Residual Deviance adalah devians model dengan prediktor. Selisih keduanya menunjukkan kontribusi prediktor dalam menjelaskan variabilitas Y.

1.8 Ringkasan Koefisien, SE, z-Value, dan p-Value

Fungsi glm() langsung menyediakan p-value melalui uji Wald: \[z = \frac{\hat{\beta}}{SE(\hat{\beta})}, \quad p = 2 \times \left(1 - \Phi(|z|)\right)\]

coef_sum_b <- as.data.frame(summary(fit_logit)$coefficients) %>%
  tibble::rownames_to_column("Variabel") %>%
  rename(Beta = Estimate, Std_Error = `Std. Error`,
         z_value = `z value`, p_value = `Pr(>|z|)`) %>%
  mutate(
    OR         = exp(Beta),
    CI_low     = exp(Beta - 1.96 * Std_Error),
    CI_high    = exp(Beta + 1.96 * Std_Error),
    Signifikan = case_when(
      p_value < 0.001 ~ "***", p_value < 0.01 ~ "**",
      p_value < 0.05  ~ "*",   p_value < 0.10 ~ ".",
      TRUE            ~ "ns"
    )
  ) %>%
  mutate(across(where(is.numeric), ~round(., 4)))

coef_sum_b %>%
  kable(caption = "Ringkasan koefisien model regresi logistik biner",
        col.names = c("Variabel", "β", "SE", "z", "p-value",
                      "OR", "CI Low", "CI High", "Sig.")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = TRUE) %>%
  footnote(general = "Kode sig.: *** p<0.001 | ** p<0.01 | * p<0.05 | . p<0.10 | ns tidak signifikan")

Ringkasan koefisien model regresi logistik biner

Variabel	β	SE	z	p-value	OR	CI Low	CI High	Sig.
(Intercept)	-0.9790	0.3939	-2.4856	0.0129	0.3757	0.1736	0.8130
Time_spent_Alone	0.0197	0.0353	0.5596	0.5758	1.0199	0.9518	1.0930	ns
Social_event_attendance	-0.1118	0.0444	-2.5173	0.0118	0.8942	0.8197	0.9756
Friends_circle_size	0.0017	0.0281	0.0621	0.9505	1.0017	0.9481	1.0584	ns
Post_frequency	-0.1633	0.0441	-3.7004	0.0002	0.8493	0.7789	0.9261	***
Drained_after_socializingYes	3.6515	0.3055	11.9532	0.0000	38.5332	21.1739	70.1244	***
Note:
Kode sig.: *** p<0.001 | ** p<0.01 | * p<0.05 | . p<0.10 | ns tidak signifikan

Kolom β adalah log-odds; tanda positif berarti meningkatkan peluang menjadi Introvert. Kolom OR adalah exp(β); nilai di atas 1 menunjukkan kenaikan prediktor meningkatkan odds menjadi Introvert, nilai di bawah 1 sebaliknya. Interval kepercayaan 95% yang tidak mencakup nilai 1 menandakan efek yang signifikan.

1.9 Interpretasi Koefisien

coef_sum_b %>%
  filter(Variabel != "(Intercept)", p_value < 0.05) %>%
  select(Variabel, Beta, OR, CI_low, CI_high, p_value, Signifikan) %>%
  kable(caption = "Variabel signifikan (p < 0.05)",
        col.names = c("Variabel", "β", "OR", "CI 2.5%", "CI 97.5%", "p-value", "Sig.")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = TRUE)

Variabel signifikan (p < 0.05)

Variabel	β	OR	CI 2.5%	CI 97.5%	p-value	Sig.
Social_event_attendance	-0.1118	0.8942	0.8197	0.9756	0.0118
Post_frequency	-0.1633	0.8493	0.7789	0.9261	0.0002	***
Drained_after_socializingYes	3.6515	38.5332	21.1739	70.1244	0.0000	***

coef_sum_b %>%
  filter(Variabel != "(Intercept)") %>%
  ggplot(aes(x = reorder(Variabel, OR), y = OR)) +
  geom_point(size = 3.5, color = "#243142") +
  geom_errorbar(aes(ymin = CI_low, ymax = CI_high),
                width = 0.25, color = "#243142", linewidth = 0.8) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "#b84f5a", linewidth = 1) +
  coord_flip() +
  scale_y_log10() +
  labs(title    = "Odds Ratio Prediktor terhadap Kepribadian Introvert",
       subtitle = "Garis merah putus-putus = OR 1 (tidak ada efek); skala log",
       x = NULL, y = "Odds Ratio") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"), panel.grid.minor = element_blank())

OR > 1 berarti prediktor meningkatkan odds seseorang menjadi Introvert relatif terhadap Extrovert. OR < 1 berarti sebaliknya. Kategori referensi model adalah Extrovert untuk Y dan “No” untuk Drained_after_socializing.

1.10 Uji Simultan

Uji ini membandingkan model lengkap dengan model null untuk menguji apakah minimal satu prediktor berpengaruh signifikan.

\[H_0: \beta_1 = \beta_2 = \cdots = \beta_p = 0\] \[H_1: \text{minimal satu } \beta_j \neq 0\]

model_null_b <- glm(Personality ~ 1, data = df_model_b, family = binomial(link = "logit"))
lrt_result_b <- lmtest::lrtest(model_null_b, fit_logit)
print(lrt_result_b)

## Likelihood ratio test
## 
## Model 1: Personality ~ 1
## Model 2: Personality ~ Time_spent_Alone + Social_event_attendance + Friends_circle_size + 
##     Post_frequency + Drained_after_socializing
##   #Df   LogLik Df  Chisq Pr(>Chisq)    
## 1   1 -2008.97                         
## 2   6  -774.77  5 2468.4  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

lrt_stat_b <- lrt_result_b$Chisq[2]
lrt_df_b   <- lrt_result_b$Df[2]
lrt_p_b    <- lrt_result_b$`Pr(>Chisq)`[2]

tibble(
  `Chi-square` = round(lrt_stat_b, 4),
  df           = lrt_df_b,
  `p-value`    = format(lrt_p_b, scientific = TRUE),
  Keputusan    = ifelse(lrt_p_b < 0.05,
                        "Tolak H0 - Model signifikan secara simultan",
                        "Gagal tolak H0")
) %>%
  kable(caption = "Likelihood Ratio Test: model penuh vs model null") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Likelihood Ratio Test: model penuh vs model null

Chi-square	df	p-value	Keputusan
2468.4	5	0e+00	Tolak H0 - Model signifikan secara simultan

Nilai p-value yang sangat kecil (< 0,05) berarti kita menolak H0. Secara simultan, setidaknya satu prediktor berpengaruh signifikan terhadap tipe kepribadian.

1.11 Uji Kelayakan Model

1.11.1 Uji Hosmer-Lemeshow

\[H_0: \text{Model fit}\] \[H_1: \text{Model tidak fit}\]

hl_b <- ResourceSelection::hoslem.test(
  as.numeric(df_model_b$Personality) - 1,
  fitted(fit_logit), g = 10
)

tibble(
  `Chi-square` = round(hl_b$statistic, 4),
  df           = hl_b$parameter,
  `p-value`    = round(hl_b$p.value, 4),
  Keputusan    = ifelse(hl_b$p.value > 0.05,
                        "Gagal tolak H0 - Model FIT",
                        "Tolak H0 - Model TIDAK FIT")
) %>%
  kable(caption = "Uji Hosmer-Lemeshow (g = 10 kelompok)") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Uji Hosmer-Lemeshow (g = 10 kelompok)

Chi-square	df	p-value	Keputusan
102.6157	8	0	Tolak H0 - Model TIDAK FIT

Nilai p > 0,05 menandakan model fit terhadap data, yaitu tidak ada perbedaan signifikan antara frekuensi yang diprediksi dan yang diamati.

1.11.2 Pseudo R² (Nagelkerke)

n_b          <- nrow(df_model_b)
ll_null_b    <- as.numeric(logLik(model_null_b))
ll_full_b    <- as.numeric(logLik(fit_logit))
cox_snell_b  <- 1 - exp((2 / n_b) * (ll_null_b - ll_full_b))
nagelkerke_b <- cox_snell_b / (1 - exp((2 / n_b) * ll_null_b))
mcfadden_b   <- 1 - (ll_full_b / ll_null_b)

tibble(
  Metode  = c("Cox & Snell R²", "Nagelkerke R²", "McFadden R²"),
  Nilai   = round(c(cox_snell_b, nagelkerke_b, mcfadden_b), 4),
  Catatan = c("Maksimum tidak mencapai 1",
              "Disesuaikan agar maksimum = 1 (direkomendasikan)",
              "Proporsi peningkatan log-likelihood")
) %>%
  kable(caption = "Pseudo R²") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Pseudo R²

Metode	Nilai	Catatan
Cox & Snell R²	0.5731	Maksimum tidak mencapai 1
Nagelkerke R²	0.7643	Disesuaikan agar maksimum = 1 (direkomendasikan)
McFadden R²	0.6143	Proporsi peningkatan log-likelihood

Nagelkerke R² sebesar 0.7643 menunjukkan bahwa sekitar 76.4% variabilitas tipe kepribadian dapat dijelaskan oleh kelima prediktor dalam model.

1.12 Prediksi Probabilitas

pred_prob_b <- predict(fit_logit, type = "response")

df_pred_b <- df_model_b %>%
  mutate(
    prob_introvert = pred_prob_b,
    prediksi       = factor(ifelse(pred_prob_b >= 0.5, "Introvert", "Extrovert"),
                            levels = c("Extrovert", "Introvert"))
  )

head(df_pred_b %>% select(Personality, prediksi, prob_introvert))

Kolom prob_introvert adalah probabilitas prediksi menjadi Introvert. Klasifikasi menggunakan threshold 0,5 yang dapat disesuaikan tergantung kebutuhan.

1.13 Confusion Matrix dan Akurasi Model

conf_biner <- table(Aktual = df_pred_b$Personality, Prediksi = df_pred_b$prediksi)
conf_biner

##            Prediksi
## Aktual      Extrovert Introvert
##   Extrovert      1380       111
##   Introvert       113      1296

accuracy_biner <- sum(diag(conf_biner)) / sum(conf_biner)
cat("\nAkurasi model:", round(accuracy_biner, 4), "\n")

## 
## Akurasi model: 0.9228

cm_b <- caret::confusionMatrix(df_pred_b$prediksi, df_pred_b$Personality, positive = "Introvert")
print(cm_b)

## Confusion Matrix and Statistics
## 
##            Reference
## Prediction  Extrovert Introvert
##   Extrovert      1380       113
##   Introvert       111      1296
##                                           
##                Accuracy : 0.9228          
##                  95% CI : (0.9124, 0.9322)
##     No Information Rate : 0.5141          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.8454          
##                                           
##  Mcnemar's Test P-Value : 0.9467          
##                                           
##             Sensitivity : 0.9198          
##             Specificity : 0.9256          
##          Pos Pred Value : 0.9211          
##          Neg Pred Value : 0.9243          
##              Prevalence : 0.4859          
##          Detection Rate : 0.4469          
##    Detection Prevalence : 0.4852          
##       Balanced Accuracy : 0.9227          
##                                           
##        'Positive' Class : Introvert       
## 

as.data.frame(conf_biner) %>%
  ggplot(aes(x = Aktual, y = Prediksi, fill = Freq)) +
  geom_tile(color = "white", linewidth = 0.8) +
  geom_text(aes(label = Freq), size = 8, fontface = "bold", color = "white") +
  scale_fill_gradient(low = "#AED6F1", high = "#1A5276") +
  labs(title    = "Confusion Matrix - Regresi Logistik Biner",
       subtitle = paste0("Akurasi: ", percent(accuracy_biner, accuracy = 0.1)),
       x = "Aktual", y = "Prediksi") +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold"), legend.position = "none")

Diagonal confusion matrix menunjukkan prediksi benar; nilai di luar diagonal adalah kesalahan klasifikasi. Perhatikan juga Sensitivity dan Specificity dari output confusionMatrix() di atas.

1.14 Kurva ROC dan AUC

roc_obj_b <- pROC::roc(df_model_b$Personality, pred_prob_b, levels = c("Extrovert", "Introvert"))
auc_val_b <- round(pROC::auc(roc_obj_b), 4)

plot(roc_obj_b, col = "#5568b8", lwd = 2.5,
     main = paste0("Kurva ROC - AUC = ", auc_val_b),
     xlab = "1 - Specificity (False Positive Rate)",
     ylab = "Sensitivity (True Positive Rate)",
     print.auc = TRUE, print.auc.y = 0.45)
abline(a = 0, b = 1, lty = 2, col = "gray60")

kat_auc_b <- dplyr::case_when(
  auc_val_b >= 0.90 ~ "Luar biasa (Excellent)", auc_val_b >= 0.80 ~ "Sangat baik (Very Good)",
  auc_val_b >= 0.70 ~ "Baik (Good)",            auc_val_b >= 0.60 ~ "Cukup (Fair)",
  TRUE              ~ "Buruk (Poor)"
)
tibble(AUC = auc_val_b, Kategori = kat_auc_b) %>%
  kable(caption = "Interpretasi nilai AUC") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Interpretasi nilai AUC

AUC	Kategori
0.9209	Luar biasa (Excellent)

AUC = 0.9209 menunjukkan kemampuan model membedakan Introvert dan Extrovert termasuk kategori Luar biasa (Excellent). Semakin mendekati 1 semakin baik; AUC = 0,5 berarti tidak lebih baik dari tebakan acak.

1.15 Visualisasi Probabilitas Prediksi

1.15.1 Probabilitas terhadap Waktu Sendiri

grid_alone <- data.frame(
  Time_spent_Alone          = seq(min(df_model_b$Time_spent_Alone),
                                   max(df_model_b$Time_spent_Alone), length.out = 100),
  Social_event_attendance   = mean(df_model_b$Social_event_attendance),
  Friends_circle_size       = mean(df_model_b$Friends_circle_size),
  Post_frequency            = mean(df_model_b$Post_frequency),
  Drained_after_socializing = factor("No", levels = c("No", "Yes"))
)
grid_alone$prob_introvert <- predict(fit_logit, newdata = grid_alone, type = "response")

ggplot(grid_alone, aes(x = Time_spent_Alone, y = prob_introvert)) +
  geom_line(color = "#5568b8", linewidth = 1.4) +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "#b84f5a", linewidth = 0.8) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
  labs(title    = "Prediksi P(Introvert) terhadap Waktu yang Dihabiskan Sendiri",
       subtitle = "Prediktor lain di-hold pada nilai rata-rata / referensi",
       x = "Time_spent_Alone (jam)", y = "P(Introvert)") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"), panel.grid.minor = element_blank())

Grafik ini menggunakan pendekatan ceteris paribus, yaitu prediktor lain ditahan pada nilai tetap. Jika garis naik ke kanan, semakin banyak waktu yang dihabiskan sendiri semakin tinggi probabilitas seseorang menjadi Introvert.

1.15.2 Probabilitas berdasarkan Drained After Socializing

grid_drained <- data.frame(
  Time_spent_Alone          = mean(df_model_b$Time_spent_Alone),
  Social_event_attendance   = mean(df_model_b$Social_event_attendance),
  Friends_circle_size       = mean(df_model_b$Friends_circle_size),
  Post_frequency            = mean(df_model_b$Post_frequency),
  Drained_after_socializing = factor(c("No", "Yes"), levels = c("No", "Yes"))
)
grid_drained$prob_introvert <- predict(fit_logit, newdata = grid_drained, type = "response")

ggplot(grid_drained, aes(x = Drained_after_socializing,
                          y = prob_introvert, fill = Drained_after_socializing)) +
  geom_col(width = 0.50, alpha = 0.93) +
  geom_text(aes(label = percent(prob_introvert, accuracy = 0.1)),
            vjust = -0.4, fontface = "bold", color = "#243142") +
  scale_fill_manual(values = c("No" = "#2f7f73", "Yes" = "#b84f5a")) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
  labs(title    = "Prediksi P(Introvert) berdasarkan Drained After Socializing",
       subtitle = "Prediktor lain di-hold pada nilai rata-rata",
       x = "Drained After Socializing", y = "P(Introvert)") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none", plot.title = element_text(face = "bold"))

Jika probabilitas Introvert pada kelompok “Yes” jauh lebih tinggi, ini sejalan dengan logika bahwa merasa lelah setelah bersosialisasi merupakan ciri khas kepribadian Introvert.

1.16 Ringkasan dan Kesimpulan

Ringkasan hasil analisis regresi logistik biner

Aspek	Nilai
Variabel Dependen	Personality (Extrovert / Introvert) - biner
Kategori Referensi	Extrovert
Jumlah Observasi	2900
Jumlah Prediktor	5 prediktor (4 numerik + 1 kategorik)
Null Deviance	4017.93
Residual Deviance	1549.53
AIC	1561.53
Uji Simultan (LRT)	Chi-square = 2468.4, p < 0.05, Signifikan
Hosmer-Lemeshow	p = 0, Tidak FIT
Nagelkerke R²	76.4% variabilitas Y dijelaskan model
AUC	0.9209 - Luar biasa (Excellent)
Akurasi Klasifikasi	92.3%

Beberapa poin utama dari hasil analisis ini:

1. Model signifikan secara simultan berdasarkan Likelihood Ratio Test, artinya kelima prediktor secara bersama-sama mampu membedakan kepribadian Introvert dan Extrovert.

2. Nilai AUC 0.9209 menunjukkan kemampuan diskriminasi model termasuk kategori Luar biasa (Excellent).

3. Semua koefisien diinterpretasikan relatif terhadap Extrovert, yaitu sebagai log-odds menjadi Introvert dengan semua prediktor lain konstan.

4. Asumsi terpenuhi: VIF < 5 untuk semua prediktor, distribusi Y bersifat biner, tidak ada missing value setelah drop_na().

---

2 Regresi Logistik Multinomial

2.1 Pendahuluan

Regresi logistik multinomial merupakan perluasan dari regresi logistik biner untuk kasus di mana variabel dependen bersifat nominal dengan tiga kategori atau lebih. Model ini membangun K-1 persamaan log-odds terhadap satu kategori referensi yang dipilih.

Tujuan analisis ini adalah memodelkan segmentasi pelanggan (A, B, C, D) berdasarkan karakteristik demografis dan perilaku mereka, menggunakan fungsi multinom() dari paket nnet.

Variabel yang digunakan:

Tipe	Variabel
Dependen (Y)	Segmentation - kategori A, B, C, D
Independen Kategorik	Gender, Ever_Married, Graduated, Spending_Score, Profession, Var_1
Independen Numerik	Age, Work_Experience, Family_Size

2.2 Memuat Paket

library(nnet)

Paket utama yang digunakan di bagian ini adalah nnet untuk estimasi model multinomial via multinom().

2.3 Memuat dan Memeriksa Data

df_m <- readxl::read_excel("Data Multinomial ADK.xlsx")
df_m$Segmentation <- as.factor(df_m$Segmentation)

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

## Jumlah baris   : 2154

cat("Jumlah kolom   :", ncol(df_m), "\n")

## Jumlah kolom   : 11

Dataset terdiri dari 2.154 observasi dengan 11 variabel.

head(df_m, 6)

str(df_m)

## tibble [2,154 × 11] (S3: tbl_df/tbl/data.frame)
##  $ ID             : num [1:2154] 459032 461899 467915 463923 467050 ...
##  $ Gender         : chr [1:2154] "Male" "Female" "Female" "Female" ...
##  $ Ever_Married   : chr [1:2154] "No" "No" "No" "No" ...
##  $ Age            : num [1:2154] 18 18 18 18 18 18 18 18 18 18 ...
##  $ Graduated      : chr [1:2154] "No" "No" "No" "No" ...
##  $ Profession     : chr [1:2154] "Doctor" "Doctor" "Doctor" "Engineer" ...
##  $ Work_Experience: num [1:2154] 0 0 1 0 1 0 0 0 0 0 ...
##  $ Spending_Score : chr [1:2154] "Low" "Low" "Low" "Low" ...
##  $ Family_Size    : num [1:2154] 3 6 4 4 5 4 4 4 2 3 ...
##  $ Var_1          : chr [1:2154] "Cat_6" "Cat_6" "Cat_6" "Cat_6" ...
##  $ Segmentation   : Factor w/ 4 levels "A","B","C","D": 4 1 2 4 4 1 2 4 3 2 ...

Variabel Segmentation memiliki 4 level (A, B, C, D), sedangkan Gender memiliki 2 level (Female, Male). Variabel numerik bertipe num, sementara beberapa variabel kategorik lain masih bertipe chr dan akan diubah pada tahap pra-pemrosesan.

2.4 Pra-Pemrosesan Data

Langkah pra-pemrosesan meliputi menghapus kolom ID, mengatur kategori referensi untuk setiap variabel kategorik, dan listwise deletion menggunakan drop_na().

df_model_m <- df_m %>%
  select(-ID) %>%
  mutate(
    Segmentation   = relevel(factor(Segmentation),   ref = "A"),
    Gender         = relevel(factor(Gender),         ref = "Male"),
    Ever_Married   = relevel(factor(Ever_Married),   ref = "No"),
    Graduated      = relevel(factor(Graduated),      ref = "No"),
    Spending_Score = relevel(factor(Spending_Score), ref = "Low"),
    Var_1          = factor(Var_1),
    Profession     = relevel(factor(Profession),     ref = "Artist")
  ) %>%
  drop_na()

cat("\nJumlah observasi setelah drop_na():", nrow(df_model_m), "\n")

## 
## Jumlah observasi setelah drop_na(): 2154

cat("Distribusi Y (Segmentation):\n")

## Distribusi Y (Segmentation):

print(table(df_model_m$Segmentation))

## 
##   A   B   C   D 
## 692 450 381 631

Setelah drop_na(), data tidak memiliki missing value. Distribusi awal menunjukkan Segmen A (692) adalah yang terbanyak, diikuti D (631), B (450), dan C (381). Kategori referensi ditetapkan sebagai Segmen A.

2.5 Analisis Deskriptif

2.5.1 Distribusi Variabel Dependen

df_model_m %>%
  count(Segmentation) %>%
  mutate(proporsi = n / sum(n)) %>%
  kable(digits = 3, caption = "Distribusi segmentasi pelanggan") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE)

Distribusi segmentasi pelanggan

Segmentation	n	proporsi
A	692	0.321
B	450	0.209
C	381	0.177
D	631	0.293

df_model_m %>%
  count(Segmentation) %>%
  mutate(proporsi = n / sum(n)) %>%
  ggplot(aes(x = Segmentation, y = proporsi, fill = Segmentation)) +
  geom_col(width = 0.60, alpha = 0.93) +
  geom_text(aes(label = percent(proporsi, accuracy = 0.1)),
            vjust = -0.4, fontface = "bold", color = "#243142") +
  scale_y_continuous(labels = percent_format(), limits = c(0, 0.45)) +
  labs(title    = "Distribusi Segmentasi Pelanggan",
       subtitle = paste0("Total observasi: ", nrow(df_model_m)),
       x = "Segmentation", y = "Proporsi") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none", plot.title = element_text(face = "bold"))

Distribusi kategori tidak seimbang. Segmen A mendominasi dengan proporsi sekitar 32,1%, diikuti Segmen D (29,3%), B (20,9%), dan C (17,7%). Ketidakseimbangan ini perlu diperhatikan saat mengevaluasi performa model.

2.5.2 Distribusi Variabel Kategorik vs Segmentasi

for (v in c("Gender", "Ever_Married", "Graduated", "Spending_Score", "Profession", "Var_1")) {
  p <- ggplot(df_model_m, aes(x = .data[[v]], fill = Segmentation)) +
    geom_bar(position = "fill", color = "white") +
    scale_y_continuous(labels = percent_format()) +
    labs(title = paste(v, "vs Segmentation"), x = v, y = "Proporsi") +
    theme_minimal(base_size = 11) +
    theme(plot.title = element_text(face = "bold"),
          axis.text.x = element_text(angle = 30, hjust = 1))
  print(p)
}

Jika komposisi warna berbeda nyata antar kategori, berarti variabel tersebut berhubungan dengan segmentasi.

2.5.3 Distribusi Variabel Numerik vs Segmentasi

for (v in c("Age", "Work_Experience", "Family_Size")) {
  p <- ggplot(df_model_m, aes(x = Segmentation, y = .data[[v]], fill = Segmentation)) +
    geom_boxplot(outlier.alpha = 0.3, outlier.size = 1) +
    labs(title = paste(v, "vs Segmentation"), x = "Segmentation", y = v) +
    theme_minimal(base_size = 11) +
    theme(legend.position = "none", plot.title = element_text(face = "bold"))
  print(p)
}

Perbedaan median atau sebaran antar segmen mengindikasikan bahwa variabel numerik tersebut berperan dalam membedakan segmen.

2.6 Pemeriksaan Asumsi

2.6.1 Multikolinearitas (VIF)

model_vif_m <- lm(as.numeric(Segmentation) ~
                    Gender + Ever_Married + Age + Graduated +
                    Profession + Work_Experience + Spending_Score +
                    Family_Size + Var_1,
                  data = df_model_m)

vif_vals_m <- vif(model_vif_m)

data.frame(
  Variabel = rownames(vif_vals_m),
  VIF      = round(vif_vals_m[, 1], 3),
  Status   = ifelse(vif_vals_m[, 1] < 5, "Aman (< 5)", "Perhatian (>= 5)")
) %>%
  kable(caption = "Variance Inflation Factor (VIF)") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Variance Inflation Factor (VIF)

	Variabel	VIF	Status
Gender	Gender	1.182	Aman (< 5)
Ever_Married	Ever_Married	2.441	Aman (< 5)
Age	Age	2.746	Aman (< 5)
Graduated	Graduated	1.336	Aman (< 5)
Profession	Profession	4.328	Aman (< 5)
Work_Experience	Work_Experience	1.114	Aman (< 5)
Spending_Score	Spending_Score	2.733	Aman (< 5)
Family_Size	Family_Size	1.338	Aman (< 5)
Var_1	Var_1	1.214	Aman (< 5)

VIF < 5 untuk semua prediktor menunjukkan tidak ada multikolinearitas yang mengkhawatirkan.

Asumsi lain yang perlu dipenuhi:

Asumsi	Status
Variabel respons bersifat nominal (tidak berurutan)	Terpenuhi: Segmen A/B/C/D tidak berurutan
Observasi independen	Diasumsikan terpenuhi
Tidak ada multikolinearitas berat (VIF < 5)	Terpenuhi (lihat tabel VIF)
Tidak ada perfect separation	Model berhasil diestimasi
Independence of Irrelevant Alternatives (IIA)	Asumsi kunci yang perlu diperhatikan

2.7 Estimasi Model Regresi Logistik Multinomial

Model diestimasi menggunakan multinom() dari paket nnet dengan Segmen A sebagai kategori referensi, menghasilkan K-1 = 3 set koefisien.

fit_multi <- nnet::multinom(
  Segmentation ~ Gender + Ever_Married + Age + Graduated +
    Profession + Work_Experience + Spending_Score + Family_Size + Var_1,
  data  = df_model_m,
  trace = FALSE
)

summary(fit_multi)

## Call:
## nnet::multinom(formula = Segmentation ~ Gender + Ever_Married + 
##     Age + Graduated + Profession + Work_Experience + Spending_Score + 
##     Family_Size + Var_1, data = df_model_m, trace = FALSE)
## 
## Coefficients:
##   (Intercept) GenderFemale Ever_MarriedYes         Age GraduatedYes
## B  -1.2701447  -0.19482277      0.09462304  0.01291404  -0.01897905
## C  -1.7858378   0.04315797     -0.06168062  0.02170057   0.13996868
## D   0.1880661   0.05485543     -0.06638361 -0.01298968  -0.06587659
##   ProfessionDoctor ProfessionEngineer ProfessionEntertainment
## B       0.07143787      -0.0065387842              0.13913858
## C      -0.22084773      -0.0274583124             -0.14923864
## D       0.03906531       0.0005041315              0.08673519
##   ProfessionExecutive ProfessionHealthcare ProfessionHomemaker ProfessionLawyer
## B         -0.02965350           0.02759224           0.4872480       -0.3385210
## C         -0.02729965          -0.43784561           0.3630431       -0.5237040
## D          0.65888639           0.57247297           0.4534136        0.6944297
##   ProfessionMarketing Work_Experience Spending_ScoreAverage Spending_ScoreHigh
## B           0.4490288     -0.02247512            0.14208647         -0.1293985
## C          -0.4948929     -0.02292035            0.35259281          0.2121841
## D           0.5945434     -0.02355472            0.09791452         -0.3455605
##   Family_Size Var_1Cat_2 Var_1Cat_3 Var_1Cat_4 Var_1Cat_5 Var_1Cat_6 Var_1Cat_7
## B -0.09661794  0.6887618 0.50748131  0.8775085 -0.2118283  0.5328975  0.7881988
## C -0.06378399  0.4974148 0.17314519  0.6233587 -0.9679091  0.3379020  0.1292760
## D -0.05012197  0.2602361 0.04940856  0.2951497  0.9292392  0.2839611  0.4278135
## 
## Std. Errors:
##   (Intercept) GenderFemale Ever_MarriedYes         Age GraduatedYes
## B   0.6720075    0.1338001       0.1899720 0.005967699    0.1452624
## C   0.6885084    0.1419684       0.2092670 0.006241037    0.1575556
## D   0.5590163    0.1216585       0.1765691 0.005838405    0.1323529
##   ProfessionDoctor ProfessionEngineer ProfessionEntertainment
## B        0.2236149          0.2420733               0.2054019
## C        0.2449923          0.2467863               0.2252111
## D        0.2126353          0.2275997               0.2010705
##   ProfessionExecutive ProfessionHealthcare ProfessionHomemaker ProfessionLawyer
## B           0.3044998            0.2400665           0.4184089        0.3066608
## C           0.3046114            0.2854205           0.4407164        0.3056672
## D           0.2841969            0.1992146           0.3839272        0.2988256
##   ProfessionMarketing Work_Experience Spending_ScoreAverage Spending_ScoreHigh
## B           0.3212151      0.01933018             0.1866082          0.2316774
## C           0.4320443      0.02100183             0.1988453          0.2351289
## D           0.2933815      0.01716602             0.1839680          0.2299233
##   Family_Size Var_1Cat_2 Var_1Cat_3 Var_1Cat_4 Var_1Cat_5 Var_1Cat_6 Var_1Cat_7
## B  0.04612512  0.6478517  0.6212026  0.6091552  0.9987976  0.5944269  0.7052194
## C  0.04968559  0.6594698  0.6336980  0.6157883  1.2281427  0.5982910  0.7663454
## D  0.03992669  0.5288130  0.5060175  0.4949472  0.6692369  0.4777975  0.5883043
## 
## Residual Deviance: 5642.048 
## AIC: 5780.048

Output summary() memberikan koefisien log-odds dan standard error untuk setiap prediktor pada masing-masing perbandingan (B vs A, C vs A, D vs A). Perhatikan bahwa multinom() tidak langsung menyediakan p-value sehingga dihitung secara manual pada bagian berikutnya.

2.8 Ringkasan Koefisien, SE, z-Value, dan p-Value

Karena multinom() tidak menyediakan p-value secara langsung, nilai tersebut dihitung menggunakan pendekatan Wald: \[z = \frac{\hat{\beta}}{SE(\hat{\beta})}, \quad p = 2 \times \left(1 - \Phi(|z|)\right)\]

multi_sum <- summary(fit_multi)

coef_long_m <- as.data.frame(multi_sum$coefficients) %>%
  tibble::rownames_to_column("kategori") %>%
  tidyr::pivot_longer(cols = -kategori, names_to = "variabel", values_to = "estimate")

se_long_m <- as.data.frame(multi_sum$standard.errors) %>%
  tibble::rownames_to_column("kategori") %>%
  tidyr::pivot_longer(cols = -kategori, names_to = "variabel", values_to = "std_error")

result_multi <- coef_long_m %>%
  left_join(se_long_m, by = c("kategori", "variabel")) %>%
  mutate(
    z_value    = estimate / std_error,
    p_value    = 2 * (1 - pnorm(abs(z_value))),
    OR         = exp(estimate),
    CI_low     = exp(estimate - 1.96 * std_error),
    CI_high    = exp(estimate + 1.96 * std_error),
    Signifikan = case_when(
      p_value < 0.001 ~ "***", p_value < 0.01 ~ "**",
      p_value < 0.05  ~ "*",   p_value < 0.10 ~ ".",
      TRUE            ~ "ns"
    )
  ) %>%
  mutate(across(where(is.numeric), ~round(., 4)))

result_multi %>%
  kable(caption = "Ringkasan koefisien model regresi logistik multinomial",
        col.names = c("Kategori vs A", "Variabel", "β", "SE", "z", "p-value",
                      "OR", "CI Low", "CI High", "Sig.")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = TRUE) %>%
  footnote(general = "Kode sig.: *** p<0.001 | ** p<0.01 | * p<0.05 | . p<0.10 | ns tidak signifikan")

Ringkasan koefisien model regresi logistik multinomial

Kategori vs A	Variabel	β	SE	z	p-value	OR	CI Low	CI High	Sig.
B	(Intercept)	-1.2701	0.6720	-1.8901	0.0587	0.2808	0.0752	1.0481	.
B	GenderFemale	-0.1948	0.1338	-1.4561	0.1454	0.8230	0.6331	1.0698	ns
B	Ever_MarriedYes	0.0946	0.1900	0.4981	0.6184	1.0992	0.7575	1.5952	ns
B	Age	0.0129	0.0060	2.1640	0.0305	1.0130	1.0012	1.0249
B	GraduatedYes	-0.0190	0.1453	-0.1307	0.8960	0.9812	0.7381	1.3044	ns
B	ProfessionDoctor	0.0714	0.2236	0.3195	0.7494	1.0741	0.6929	1.6648	ns
B	ProfessionEngineer	-0.0065	0.2421	-0.0270	0.9785	0.9935	0.6182	1.5967	ns
B	ProfessionEntertainment	0.1391	0.2054	0.6774	0.4982	1.1493	0.7684	1.7190	ns
B	ProfessionExecutive	-0.0297	0.3045	-0.0974	0.9224	0.9708	0.5345	1.7633	ns
B	ProfessionHealthcare	0.0276	0.2401	0.1149	0.9085	1.0280	0.6421	1.6456	ns
B	ProfessionHomemaker	0.4872	0.4184	1.1645	0.2442	1.6278	0.7169	3.6963	ns
B	ProfessionLawyer	-0.3385	0.3067	-1.1039	0.2696	0.7128	0.3908	1.3002	ns
B	ProfessionMarketing	0.4490	0.3212	1.3979	0.1621	1.5668	0.8348	2.9406	ns
B	Work_Experience	-0.0225	0.0193	-1.1627	0.2450	0.9778	0.9414	1.0155	ns
B	Spending_ScoreAverage	0.1421	0.1866	0.7614	0.4464	1.1527	0.7996	1.6617	ns
B	Spending_ScoreHigh	-0.1294	0.2317	-0.5585	0.5765	0.8786	0.5579	1.3836	ns
B	Family_Size	-0.0966	0.0461	-2.0947	0.0362	0.9079	0.8294	0.9938
B	Var_1Cat_2	0.6888	0.6479	1.0631	0.2877	1.9912	0.5593	7.0890	ns
B	Var_1Cat_3	0.5075	0.6212	0.8169	0.4140	1.6611	0.4916	5.6127	ns
B	Var_1Cat_4	0.8775	0.6092	1.4405	0.1497	2.4049	0.7287	7.9363	ns
B	Var_1Cat_5	-0.2118	0.9988	-0.2121	0.8320	0.8091	0.1142	5.7306	ns
B	Var_1Cat_6	0.5329	0.5944	0.8965	0.3700	1.7039	0.5314	5.4629	ns
B	Var_1Cat_7	0.7882	0.7052	1.1177	0.2637	2.1994	0.5521	8.7620	ns
C	(Intercept)	-1.7858	0.6885	-2.5938	0.0095	0.1677	0.0435	0.6464	**
C	GenderFemale	0.0432	0.1420	0.3040	0.7611	1.0441	0.7905	1.3791	ns
C	Ever_MarriedYes	-0.0617	0.2093	-0.2947	0.7682	0.9402	0.6239	1.4169	ns
C	Age	0.0217	0.0062	3.4771	0.0005	1.0219	1.0095	1.0345	***
C	GraduatedYes	0.1400	0.1576	0.8884	0.3743	1.1502	0.8446	1.5664	ns
C	ProfessionDoctor	-0.2208	0.2450	-0.9014	0.3674	0.8018	0.4961	1.2961	ns
C	ProfessionEngineer	-0.0275	0.2468	-0.1113	0.9114	0.9729	0.5998	1.5781	ns
C	ProfessionEntertainment	-0.1492	0.2252	-0.6627	0.5075	0.8614	0.5540	1.3393	ns
C	ProfessionExecutive	-0.0273	0.3046	-0.0896	0.9286	0.9731	0.5356	1.7678	ns
C	ProfessionHealthcare	-0.4378	0.2854	-1.5340	0.1250	0.6454	0.3689	1.1293	ns
C	ProfessionHomemaker	0.3630	0.4407	0.8238	0.4101	1.4377	0.6061	3.4105	ns
C	ProfessionLawyer	-0.5237	0.3057	-1.7133	0.0867	0.5923	0.3254	1.0783	.
C	ProfessionMarketing	-0.4949	0.4320	-1.1455	0.2520	0.6096	0.2614	1.4218	ns
C	Work_Experience	-0.0229	0.0210	-1.0914	0.2751	0.9773	0.9379	1.0184	ns
C	Spending_ScoreAverage	0.3526	0.1988	1.7732	0.0762	1.4228	0.9635	2.1008	.
C	Spending_ScoreHigh	0.2122	0.2351	0.9024	0.3668	1.2364	0.7798	1.9602	ns
C	Family_Size	-0.0638	0.0497	-1.2838	0.1992	0.9382	0.8511	1.0342	ns
C	Var_1Cat_2	0.4974	0.6595	0.7543	0.4507	1.6445	0.4515	5.9893	ns
C	Var_1Cat_3	0.1731	0.6337	0.2732	0.7847	1.1890	0.3434	4.1173	ns
C	Var_1Cat_4	0.6234	0.6158	1.0123	0.3114	1.8652	0.5579	6.2358	ns
C	Var_1Cat_5	-0.9679	1.2281	-0.7881	0.4306	0.3799	0.0342	4.2175	ns
C	Var_1Cat_6	0.3379	0.5983	0.5648	0.5722	1.4020	0.4340	4.5292	ns
C	Var_1Cat_7	0.1293	0.7663	0.1687	0.8660	1.1380	0.2534	5.1106	ns
D	(Intercept)	0.1881	0.5590	0.3364	0.7366	1.2069	0.4035	3.6101	ns
D	GenderFemale	0.0549	0.1217	0.4509	0.6521	1.0564	0.8323	1.3409	ns
D	Ever_MarriedYes	-0.0664	0.1766	-0.3760	0.7069	0.9358	0.6620	1.3227	ns
D	Age	-0.0130	0.0058	-2.2249	0.0261	0.9871	0.9759	0.9985
D	GraduatedYes	-0.0659	0.1324	-0.4977	0.6187	0.9362	0.7223	1.2135	ns
D	ProfessionDoctor	0.0391	0.2126	0.1837	0.8542	1.0398	0.6854	1.5775	ns
D	ProfessionEngineer	0.0005	0.2276	0.0022	0.9982	1.0005	0.6404	1.5630	ns
D	ProfessionEntertainment	0.0867	0.2011	0.4314	0.6662	1.0906	0.7354	1.6174	ns
D	ProfessionExecutive	0.6589	0.2842	2.3184	0.0204	1.9326	1.1072	3.3734
D	ProfessionHealthcare	0.5725	0.1992	2.8736	0.0041	1.7726	1.1996	2.6194	**
D	ProfessionHomemaker	0.4534	0.3839	1.1810	0.2376	1.5737	0.7415	3.3398	ns
D	ProfessionLawyer	0.6944	0.2988	2.3239	0.0201	2.0026	1.1149	3.5971
D	ProfessionMarketing	0.5945	0.2934	2.0265	0.0427	1.8122	1.0197	3.2206
D	Work_Experience	-0.0236	0.0172	-1.3722	0.1700	0.9767	0.9444	1.0101	ns
D	Spending_ScoreAverage	0.0979	0.1840	0.5322	0.5946	1.1029	0.7690	1.5817	ns
D	Spending_ScoreHigh	-0.3456	0.2299	-1.5029	0.1329	0.7078	0.4510	1.1108	ns
D	Family_Size	-0.0501	0.0399	-1.2553	0.2094	0.9511	0.8795	1.0285	ns
D	Var_1Cat_2	0.2602	0.5288	0.4921	0.6226	1.2972	0.4601	3.6572	ns
D	Var_1Cat_3	0.0494	0.5060	0.0976	0.9222	1.0506	0.3897	2.8326	ns
D	Var_1Cat_4	0.2951	0.4949	0.5963	0.5510	1.3433	0.5092	3.5440	ns
D	Var_1Cat_5	0.9292	0.6692	1.3885	0.1650	2.5326	0.6822	9.4022	ns
D	Var_1Cat_6	0.2840	0.4778	0.5943	0.5523	1.3284	0.5207	3.3887	ns
D	Var_1Cat_7	0.4278	0.5883	0.7272	0.4671	1.5339	0.4842	4.8593	ns
Note:
Kode sig.: *** p<0.001 | ** p<0.01 | * p<0.05 | . p<0.10 | ns tidak signifikan

Kolom β adalah log-odds; tanda positif berarti meningkatkan peluang masuk kategori tersebut relatif terhadap Segmen A. OR > 1 berarti kenaikan prediktor meningkatkan odds relatif; OR < 1 berarti sebaliknya. Interval kepercayaan 95% yang tidak mencakup nilai 1 menandakan efek yang signifikan.

2.9 Interpretasi Koefisien

result_multi %>%
  filter(variabel != "(Intercept)", p_value < 0.05) %>%
  select(kategori, variabel, estimate, OR, CI_low, CI_high, p_value, Signifikan) %>%
  kable(caption = "Variabel signifikan (p < 0.05) per kategori",
        col.names = c("Kategori vs A", "Variabel", "β", "OR", "CI 2.5%", "CI 97.5%", "p-value", "Sig.")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = TRUE)

Variabel signifikan (p < 0.05) per kategori

Kategori vs A	Variabel	β	OR	CI 2.5%	CI 97.5%	p-value	Sig.
B	Age	0.0129	1.0130	1.0012	1.0249	0.0305
B	Family_Size	-0.0966	0.9079	0.8294	0.9938	0.0362
C	Age	0.0217	1.0219	1.0095	1.0345	0.0005	***
D	Age	-0.0130	0.9871	0.9759	0.9985	0.0261
D	ProfessionExecutive	0.6589	1.9326	1.1072	3.3734	0.0204
D	ProfessionHealthcare	0.5725	1.7726	1.1996	2.6194	0.0041	**
D	ProfessionLawyer	0.6944	2.0026	1.1149	3.5971	0.0201
D	ProfessionMarketing	0.5945	1.8122	1.0197	3.2206	0.0427

OR > 1 berarti prediktor meningkatkan odds pelanggan masuk segmen tersebut relatif terhadap Segmen A. OR < 1 berarti sebaliknya. Kategori referensi: Segmen A (Y), Male (Gender), No (Ever_Married & Graduated), Low (Spending_Score), Artist (Profession).

2.10 Uji Simultan

\[H_0: \beta_1 = \beta_2 = \cdots = \beta_p = 0\] \[H_1: \text{minimal satu } \beta_j \neq 0\]

model_null_m <- nnet::multinom(Segmentation ~ 1, data = df_model_m, trace = FALSE)

lrt_stat_m <- 2 * (as.numeric(logLik(fit_multi)) - as.numeric(logLik(model_null_m)))
lrt_df_m   <- length(coef(fit_multi)) - length(coef(model_null_m))
lrt_p_m    <- 1 - pchisq(lrt_stat_m, lrt_df_m)

tibble(
  `Chi-square` = round(lrt_stat_m, 4),
  df           = lrt_df_m,
  `p-value`    = format(lrt_p_m, scientific = TRUE),
  Keputusan    = ifelse(lrt_p_m < 0.05,
                        "Tolak H0 - Model signifikan secara simultan",
                        "Gagal tolak H0")
) %>%
  kable(caption = "Likelihood Ratio Test: model penuh vs model null") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Likelihood Ratio Test: model penuh vs model null

Chi-square	df	p-value	Keputusan
208.1821	66	1.110223e-16	Tolak H0 - Model signifikan secara simultan

Nilai p-value yang sangat kecil berarti kita menolak H0. Secara simultan, setidaknya satu prediktor berpengaruh signifikan terhadap segmentasi pelanggan.

2.11 Prediksi Probabilitas

pred_prob_multi <- predict(fit_multi, type = "probs")
head(pred_prob_multi)

##           A         B          C         D
## 1 0.3508102 0.1702103 0.08069851 0.3982809
## 2 0.3954009 0.1181569 0.07842799 0.4080142
## 3 0.3719057 0.1318301 0.08190571 0.4143584
## 4 0.3695076 0.1239085 0.10102918 0.4055547
## 5 0.3791251 0.1318955 0.11148754 0.3774919
## 6 0.3441595 0.1560254 0.09771703 0.4020980

df_pred_m <- df_model_m %>%
  bind_cols(as.data.frame(pred_prob_multi)) %>%
  mutate(prediksi = predict(fit_multi, type = "class"))

head(df_pred_m %>% select(Segmentation, prediksi, A, B, C, D))

Setiap baris menunjukkan probabilitas prediksi untuk keempat segmen. Kolom prediksi adalah segmen dengan probabilitas tertinggi.

2.12 Confusion Matrix dan Akurasi Model

conf_multi <- table(Aktual = df_pred_m$Segmentation, Prediksi = df_pred_m$prediksi)
conf_multi

##       Prediksi
## Aktual   A   B   C   D
##      A 441  29  44 178
##      B 279  44  30  97
##      C 252  31  40  58
##      D 306  32  33 260

accuracy_multi <- sum(diag(conf_multi)) / sum(conf_multi)
cat("\nAkurasi model:", round(accuracy_multi, 4), "\n")

## 
## Akurasi model: 0.3644

as.data.frame(conf_multi) %>%
  ggplot(aes(x = Aktual, y = Prediksi, fill = Freq)) +
  geom_tile(color = "white", linewidth = 0.8) +
  geom_text(aes(label = Freq), size = 6, fontface = "bold", color = "white") +
  scale_fill_gradient(low = "#AED6F1", high = "#1A5276") +
  labs(title    = "Confusion Matrix - Regresi Logistik Multinomial",
       subtitle = paste0("Akurasi: ", percent(accuracy_multi, accuracy = 0.1)),
       x = "Aktual", y = "Prediksi") +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold"), legend.position = "none")

Diagonal confusion matrix menunjukkan prediksi benar. Nilai di luar diagonal adalah kesalahan klasifikasi. Model paling sering salah mengklasifikasikan observasi ke Segmen A atau D, yang merupakan segmen terbesar.

2.13 Visualisasi Probabilitas Prediksi

2.13.1 Probabilitas terhadap Usia

grid_multi <- expand.grid(
  Gender          = "Male", Ever_Married = "Yes",
  Age             = seq(min(df_model_m$Age), max(df_model_m$Age), length.out = 100),
  Graduated       = "Yes", Profession = "Artist",
  Work_Experience = mean(df_model_m$Work_Experience, na.rm = TRUE),
  Spending_Score  = "Low",
  Family_Size     = median(df_model_m$Family_Size, na.rm = TRUE),
  Var_1           = "Cat_6"
)

grid_prob_m <- predict(fit_multi, newdata = grid_multi, type = "probs")

grid_multi %>%
  bind_cols(as.data.frame(grid_prob_m)) %>%
  tidyr::pivot_longer(cols = c("A", "B", "C", "D"),
                      names_to = "Segmentation", values_to = "Probabilitas") %>%
  ggplot(aes(x = Age, y = Probabilitas, color = Segmentation)) +
  geom_line(linewidth = 1.2) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
  scale_color_manual(values = c("A" = "#2f7f73", "B" = "#5568b8",
                                "C" = "#d18b2f", "D" = "#b84f5a")) +
  labs(title    = "Prediksi Probabilitas Segmentasi terhadap Usia",
       subtitle = "Prediktor lain di-hold pada nilai referensi / rata-rata",
       x = "Age", y = "Probabilitas", color = "Segmentation") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"), panel.grid.minor = element_blank())

Grafik ini menggunakan pendekatan ceteris paribus. Tren garis menunjukkan arah pengaruh usia terhadap keanggotaan segmen. Jika garis Segmen D naik seiring bertambahnya usia sementara Segmen A menurun, pelanggan yang lebih tua cenderung berada di Segmen D.

2.13.2 Probabilitas berdasarkan Spending Score

grid_spend <- expand.grid(
  Gender = "Male", Ever_Married = "Yes",
  Age    = round(mean(df_model_m$Age)), Graduated = "Yes", Profession = "Artist",
  Work_Experience = mean(df_model_m$Work_Experience, na.rm = TRUE),
  Spending_Score  = c("Low", "Average", "High"),
  Family_Size     = median(df_model_m$Family_Size, na.rm = TRUE),
  Var_1           = "Cat_6"
)

grid_spend_prob <- predict(fit_multi, newdata = grid_spend, type = "probs")

grid_spend %>%
  bind_cols(as.data.frame(grid_spend_prob)) %>%
  tidyr::pivot_longer(cols = c("A", "B", "C", "D"),
                      names_to = "Segmentation", values_to = "Probabilitas") %>%
  mutate(Spending_Score = factor(Spending_Score, levels = c("Low", "Average", "High"))) %>%
  ggplot(aes(x = Spending_Score, y = Probabilitas, fill = Segmentation)) +
  geom_col(position = "dodge", width = 0.65, alpha = 0.92) +
  scale_y_continuous(labels = percent_format()) +
  scale_fill_manual(values = c("A" = "#2f7f73", "B" = "#5568b8",
                               "C" = "#d18b2f", "D" = "#b84f5a")) +
  labs(title    = "Prediksi Probabilitas Segmentasi berdasarkan Spending Score",
       subtitle = "Prediktor lain di-hold pada nilai referensi / rata-rata",
       x = "Spending Score", y = "Probabilitas", fill = "Segmentation") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))

Pergeseran distribusi batang antar level Spending Score mengindikasikan seberapa besar skor pengeluaran memengaruhi keanggotaan segmen.

2.14 Ringkasan dan Kesimpulan

Ringkasan hasil analisis regresi logistik multinomial

Aspek	Nilai
Variabel Dependen	Segmentation (A, B, C, D) - nominal
Kategori Referensi	Segmen A
Jumlah Observasi	2.154
Jumlah Prediktor	9 prediktor (6 kategorik + 3 numerik)
Residual Deviance	5642.05
AIC	5780.05
Uji Simultan (LRT)	Chi-square = 208.18, p < 0.05, Signifikan
Akurasi Klasifikasi	36.4% (relatif rendah)

Beberapa poin utama dari hasil analisis ini:

1. Model signifikan secara simultan berdasarkan Likelihood Ratio Test.

2. Akurasi sekitar 36,4% terbilang rendah, mengindikasikan tumpang tindih antar segmen yang cukup besar.

3. Semua koefisien diinterpretasikan relatif terhadap Segmen A.

4. Asumsi terpenuhi: VIF < 5 untuk semua prediktor, tidak ada missing value setelah drop_na().

---

3 Regresi Logistik Ordinal

3.1 Pendahuluan

Regresi logistik ordinal digunakan ketika variabel dependen bersifat ordinal, yaitu memiliki urutan alami yang bermakna antar kategorinya. Model ini menggunakan cumulative logit untuk memperkirakan peluang kumulatif P(Y ≤ j) pada setiap cutpoint j.

Tujuan analisis ini adalah memodelkan variant kendaraan listrik (Base, Standard, Long Range, Premium, Performance) berdasarkan spesifikasi teknis dan harga, menggunakan fungsi polr() dari paket MASS.

Variabel yang digunakan:

Tipe	Variabel
Dependen (Y)	variant - ordinal: Base < Standard < Long Range < Premium < Performance
Independen (X)	battery_capacity_kwh, range_miles, charging_speed_kw, price_usd (numerik)

3.2 Memuat Paket

library(MASS)
library(ordinal)

select  <- dplyr::select
filter  <- dplyr::filter
rename  <- dplyr::rename

Paket utama yang digunakan adalah MASS untuk estimasi model ordinal via polr() dan ordinal untuk uji asumsi proportional odds via nominal_test().

3.3 Memuat dan Memeriksa Data

df_o_raw <- read_excel("Data Ordinal.xls")

df_o <- df_o_raw %>%
  dplyr::select(variant, battery_capacity_kwh, range_miles,
                charging_speed_kw, price_usd) %>%
  drop_na()

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

## Jumlah baris   : 2000

cat("Jumlah kolom   :", ncol(df_o), "\n")

## Jumlah kolom   : 5

head(df_o, 6)

str(df_o)

## tibble [2,000 × 5] (S3: tbl_df/tbl/data.frame)
##  $ variant             : chr [1:2000] "Performance" "Premium" "Premium" "Long Range" ...
##  $ battery_capacity_kwh: num [1:2000] 118.7 58.8 58.2 102.5 93.9 ...
##  $ range_miles         : num [1:2000] 400 219 225 349 314 251 222 292 184 251 ...
##  $ charging_speed_kw   : num [1:2000] 234.5 148.1 104.9 66.5 298.5 ...
##  $ price_usd           : num [1:2000] 104881 48217 49651 38132 144080 ...

Variabel variant memiliki 5 kategori dengan urutan alami. Pada tahap berikutnya, variabel ini akan dikonversi menjadi ordered factor sesuai urutan yang bermakna.

3.4 Pra-Pemrosesan Data

Dua langkah utama pra-pemrosesan: mengubah Y menjadi ordered factor sesuai hierarki variant, dan standarisasi prediktor karena price_usd memiliki skala yang jauh berbeda sehingga dapat menyebabkan masalah numerik pada polr().

df_o <- df_o %>%
  mutate(
    variant = ordered(variant,
                      levels = c("Base", "Standard", "Long Range",
                                 "Premium", "Performance"))
  )

cat("Distribusi variant (ordered factor):\n")

## Distribusi variant (ordered factor):

print(table(df_o$variant))

## 
##        Base    Standard  Long Range     Premium Performance 
##         420         403         407         392         378

cat("\nUrutan level:\n")

## 
## Urutan level:

print(levels(df_o$variant))

## [1] "Base"        "Standard"    "Long Range"  "Premium"     "Performance"

df_o_scaled <- df_o %>%
  mutate(
    battery_capacity_kwh = as.numeric(scale(battery_capacity_kwh)),
    range_miles          = as.numeric(scale(range_miles)),
    charging_speed_kw    = as.numeric(scale(charging_speed_kw)),
    price_usd            = as.numeric(scale(price_usd))
  )

cat("Cek setelah scaling (harus tidak ada NA/Inf):\n")

## Cek setelah scaling (harus tidak ada NA/Inf):

print(summary(df_o_scaled[, c("battery_capacity_kwh", "range_miles",
                               "charging_speed_kw", "price_usd")]))

##  battery_capacity_kwh  range_miles       charging_speed_kw   price_usd      
##  Min.   :-1.72599     Min.   :-1.57438   Min.   :-1.2991   Min.   :-1.7701  
##  1st Qu.:-0.80605     1st Qu.:-0.81412   1st Qu.:-0.7932   1st Qu.:-0.7521  
##  Median : 0.02177     Median :-0.01162   Median :-0.2853   Median :-0.1660  
##  Mean   : 0.00000     Mean   : 0.00000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.64170     3rd Qu.: 0.67824   3rd Qu.: 0.7668   3rd Qu.: 0.5923  
##  Max.   : 2.25253     Max.   : 2.60705   Max.   : 2.4042   Max.   : 5.4085

Setelah standarisasi, setiap prediktor memiliki mean = 0 dan sd = 1. Koefisien yang dihasilkan merefleksikan perubahan log-odds per 1 SD kenaikan prediktor.

3.5 Analisis Deskriptif

3.5.1 Distribusi Variabel Dependen

df_o %>%
  count(variant) %>%
  mutate(proporsi = n / sum(n)) %>%
  kable(digits = 3, caption = "Distribusi variant kendaraan listrik") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE)

Distribusi variant kendaraan listrik

variant	n	proporsi
Base	420	0.210
Standard	403	0.202
Long Range	407	0.203
Premium	392	0.196
Performance	378	0.189

df_o %>%
  count(variant) %>%
  mutate(proporsi = n / sum(n)) %>%
  ggplot(aes(x = variant, y = proporsi, fill = variant)) +
  geom_col(width = 0.65, alpha = 0.94) +
  geom_text(aes(label = percent(proporsi, accuracy = 0.1)),
            vjust = -0.4, fontface = "bold", color = "#243142") +
  scale_y_continuous(labels = percent_format(), limits = c(0, 0.35)) +
  scale_fill_manual(values = c("Base" = "#2f7f73", "Standard" = "#5568b8",
                               "Long Range" = "#d18b2f", "Premium" = "#b84f5a",
                               "Performance" = "#243142")) +
  labs(title    = "Distribusi Variant Kendaraan Listrik",
       subtitle = paste0("Total observasi: ", nrow(df_o),
                         " | Respons ordinal: Base < Standard < Long Range < Premium < Performance"),
       x = "Variant", y = "Proporsi") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none", plot.title = element_text(face = "bold"))

Distribusi yang relatif merata antar kategori idealnya mendukung estimasi model yang stabil.

3.5.2 Statistik Deskriptif Prediktor per Variant

df_o %>%
  group_by(variant) %>%
  summarise(across(
    c(battery_capacity_kwh, range_miles, charging_speed_kw, price_usd),
    list(mean = ~round(mean(.x), 2), sd = ~round(sd(.x), 2))
  )) %>%
  kable(caption = "Rata-rata dan SD prediktor per kategori variant") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = TRUE)

Rata-rata dan SD prediktor per kategori variant

variant	battery_capacity_kwh_mean	battery_capacity_kwh_sd	range_miles_mean	range_miles_sd	charging_speed_kw_mean	charging_speed_kw_sd	price_usd_mean	price_usd_sd
Base	62.94	13.13	221.92	47.39	144.03	82.07	61567.79	26246.27
Standard	61.64	13.27	216.52	47.71	136.80	79.95	68175.93	29524.52
Long Range	92.51	9.96	323.65	37.99	137.61	77.36	80786.13	33190.99
Premium	62.20	12.59	218.01	45.26	147.60	85.55	92192.48	37288.88
Performance	95.28	14.90	333.33	54.14	214.08	45.46	93637.66	37675.69

Perbedaan rata-rata yang konsisten dan monoton antar variant mengindikasikan prediktor tersebut berkorelasi kuat dengan variabel dependen.

3.5.3 Boxplot Prediktor vs Variant

for (v in c("battery_capacity_kwh", "range_miles", "charging_speed_kw", "price_usd")) {
  p <- ggplot(df_o, aes(x = variant, y = .data[[v]], fill = variant)) +
    geom_boxplot(outlier.alpha = 0.3, outlier.size = 1) +
    scale_fill_manual(values = c("Base" = "#2f7f73", "Standard" = "#5568b8",
                                 "Long Range" = "#d18b2f", "Premium" = "#b84f5a",
                                 "Performance" = "#243142")) +
    labs(title = paste(v, "vs Variant"), x = "Variant", y = v) +
    theme_minimal(base_size = 11) +
    theme(legend.position = "none", plot.title = element_text(face = "bold"))
  print(p)
}

Pola yang diharapkan pada data ordinal adalah distribusi prediktor yang bergeser secara sistematis mengikuti urutan variant.

3.6 Pemeriksaan Asumsi

3.6.1 Multikolinearitas (VIF)

model_vif_o <- lm(as.numeric(variant) ~
                    battery_capacity_kwh + range_miles +
                    charging_speed_kw + price_usd,
                  data = df_o_scaled)

vif_vals_o <- vif(model_vif_o)

data.frame(
  Variabel = names(vif_vals_o),
  VIF      = round(vif_vals_o, 3),
  Status   = ifelse(vif_vals_o < 5, "Aman (< 5)", "Perhatian (>= 5)")
) %>%
  kable(caption = "Variance Inflation Factor (VIF)") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Variance Inflation Factor (VIF)

	Variabel	VIF	Status
battery_capacity_kwh	battery_capacity_kwh	22.923	Perhatian (>= 5)
range_miles	range_miles	22.845	Perhatian (>= 5)
charging_speed_kw	charging_speed_kw	1.115	Aman (< 5)
price_usd	price_usd	1.118	Aman (< 5)

VIF < 5 untuk semua prediktor menunjukkan tidak ada multikolinearitas yang mengkhawatirkan.

Asumsi lain:

Asumsi	Status
Variabel respons bersifat ordinal dengan urutan bermakna	Terpenuhi: Base < Standard < Long Range < Premium < Performance
Observasi independen	Diasumsikan terpenuhi
Tidak ada multikolinearitas berat (VIF < 5)	Terpenuhi (lihat tabel VIF)
Proportional Odds - efek prediktor sama di semua cutpoint	Diuji formal pada bagian selanjutnya

3.7 Estimasi Model Regresi Logistik Ordinal

Model diestimasi menggunakan polr() dari paket MASS. Perlu diperhatikan bahwa konvensi polr() menggunakan persamaan:

\[\text{logit}\{P(Y \leq j)\} = \alpha_j - \mathbf{x}'\boldsymbol{\beta}\]

sehingga tanda koefisien polr berlawanan dengan konvensi slide/Agresti, yaitu \(\beta_{\text{slide}} = -\beta_{\text{polr}}\).

fit_ord <- MASS::polr(
  variant ~ battery_capacity_kwh + range_miles + charging_speed_kw + price_usd,
  data   = df_o_scaled,
  method = "logistic",
  Hess   = TRUE
)

summary(fit_ord)

## Call:
## MASS::polr(formula = variant ~ battery_capacity_kwh + range_miles + 
##     charging_speed_kw + price_usd, data = df_o_scaled, Hess = TRUE, 
##     method = "logistic")
## 
## Coefficients:
##                        Value Std. Error t value
## battery_capacity_kwh  0.9527    0.19661  4.8458
## range_miles          -0.1141    0.19487 -0.5853
## charging_speed_kw     0.2593    0.04410  5.8797
## price_usd             0.5468    0.04564 11.9800
## 
## Intercepts:
##                     Value    Std. Error t value 
## Base|Standard        -1.6868   0.0620   -27.2019
## Standard|Long Range  -0.4434   0.0508    -8.7297
## Long Range|Premium    0.5793   0.0513    11.3000
## Premium|Performance   1.7596   0.0645    27.2710
## 
## Residual Deviance: 5740.981 
## AIC: 5756.981

Output terdiri dari dua bagian: Coefficients (slope untuk setiap prediktor) dan Intercepts (cutpoints yang memisahkan antar kategori ordinal). Perlu diingat bahwa polr() tidak langsung menyediakan p-value.

3.8 Ringkasan Koefisien, SE, z-Value, p-Value, dan Odds Ratio

ord_coef <- coef(summary(fit_ord))

result_ord <- as.data.frame(ord_coef) %>%
  tibble::rownames_to_column("parameter") %>%
  dplyr::rename(estimate = Value, std_error = `Std. Error`, t_value = `t value`) %>%
  mutate(
    p_value        = 2 * (1 - pnorm(abs(t_value))),
    jenis          = ifelse(grepl("\\|", parameter), "Cutpoint", "Koefisien"),
    estimate_slide = ifelse(jenis == "Koefisien", -estimate, estimate),
    OR_polr        = ifelse(jenis == "Koefisien", exp(estimate),  NA_real_),
    OR_slide       = ifelse(jenis == "Koefisien", exp(-estimate), NA_real_),
    CI_low_polr    = ifelse(jenis == "Koefisien", exp(estimate - 1.96 * std_error), NA_real_),
    CI_high_polr   = ifelse(jenis == "Koefisien", exp(estimate + 1.96 * std_error), NA_real_),
    Signifikan     = case_when(
      jenis == "Cutpoint" ~ "-",
      p_value < 0.001     ~ "***", p_value < 0.01 ~ "**",
      p_value < 0.05      ~ "*",   p_value < 0.10 ~ ".",
      TRUE                ~ "ns"
    )
  )

result_ord %>%
  mutate(across(c(estimate, estimate_slide, std_error, t_value, p_value,
                  OR_polr, OR_slide, CI_low_polr, CI_high_polr), ~round(.x, 4))) %>%
  kable(caption = "Ringkasan hasil regresi logistik ordinal",
        col.names = c("Parameter", "Estimate polr", "SE", "z/t-value", "p-value",
                      "Jenis", "Estimate slide", "OR polr", "OR slide",
                      "CI polr 2.5%", "CI polr 97.5%", "Sig.")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = TRUE) %>%
  footnote(general = "OR polr: exp(beta_polr). OR slide: exp(-beta_polr). Tanda koefisien polr berlawanan dengan slide/PDF.")

Ringkasan hasil regresi logistik ordinal

Parameter	Estimate polr	SE	z/t-value	p-value	Jenis	Estimate slide	OR polr	OR slide	CI polr 2.5%	CI polr 97.5%	Sig.
battery_capacity_kwh	0.9527	0.1966	4.8458	0.0000	Koefisien	-0.9527	2.5928	0.3857	1.7636	3.8117	***
range_miles	-0.1141	0.1949	-0.5853	0.5583	Koefisien	0.1141	0.8922	1.1208	0.6090	1.3072	ns
charging_speed_kw	0.2593	0.0441	5.8797	0.0000	Koefisien	-0.2593	1.2960	0.7716	1.1887	1.4130	***
price_usd	0.5468	0.0456	11.9800	0.0000	Koefisien	-0.5468	1.7277	0.5788	1.5799	1.8894	***
Base|Standard	-1.6868	0.0620	-27.2019	0.0000	Cutpoint	-1.6868	NA	NA	NA	NA
Standard|Long Range	-0.4434	0.0508	-8.7297	0.0000	Cutpoint	-0.4434	NA	NA	NA	NA
Long Range|Premium	0.5793	0.0513	11.3000	0.0000	Cutpoint	0.5793	NA	NA	NA	NA
Premium|Performance	1.7596	0.0645	27.2710	0.0000	Cutpoint	1.7596	NA	NA	NA	NA
Note:
OR polr: exp(beta_polr). OR slide: exp(-beta_polr). Tanda koefisien polr berlawanan dengan slide/PDF.

Estimate polr yang positif berarti kenaikan prediktor menurunkan peluang kumulatif P(Y <= j), sehingga respons cenderung bergeser ke kategori lebih tinggi. OR slide > 1 berarti prediktor meningkatkan odds kumulatif P(Y <= j)/P(Y > j), atau respons cenderung ke variant lebih rendah.

3.9 Koefisien Signifikan

result_ord %>%
  filter(jenis == "Koefisien", p_value < 0.05) %>%
  mutate(across(c(estimate, estimate_slide, std_error, t_value, p_value,
                  OR_polr, OR_slide, CI_low_polr, CI_high_polr), ~round(.x, 4))) %>%
  select(parameter, estimate, OR_polr, OR_slide, CI_low_polr, CI_high_polr, p_value, Signifikan) %>%
  kable(caption = "Koefisien signifikan (p < 0.05)",
        col.names = c("Variabel", "β (polr)", "OR polr", "OR slide",
                      "CI 2.5%", "CI 97.5%", "p-value", "Sig.")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Koefisien signifikan (p < 0.05)

Variabel	β (polr)	OR polr	OR slide	CI 2.5%	CI 97.5%	p-value	Sig.
battery_capacity_kwh	0.9527	2.5928	0.3857	1.7636	3.8117	0	***
charging_speed_kw	0.2593	1.2960	0.7716	1.1887	1.4130	0	***
price_usd	0.5468	1.7277	0.5788	1.5799	1.8894	0	***

OR slide < 1 berarti kenaikan prediktor menurunkan odds kumulatif, sehingga respons cenderung ke variant yang lebih tinggi. Contoh: jika OR slide untuk price_usd = 0,4, artinya setiap kenaikan 1 SD harga menurunkan odds kumulatif sebesar 60%, atau kendaraan yang lebih mahal cenderung berada di variant yang lebih tinggi.

3.10 Uji Simultan

\[H_0: \beta_1 = \beta_2 = \cdots = \beta_p = 0\] \[H_1: \text{minimal satu } \beta_j \neq 0\]

model_null_o <- MASS::polr(variant ~ 1, data = df_o_scaled, method = "logistic", Hess = TRUE)

lrt_stat_o <- 2 * (as.numeric(logLik(fit_ord)) - as.numeric(logLik(model_null_o)))
lrt_df_o   <- length(coef(fit_ord)) - length(coef(model_null_o))
lrt_p_o    <- 1 - pchisq(lrt_stat_o, lrt_df_o)

tibble(
  `Chi-square` = round(lrt_stat_o, 4),
  df           = lrt_df_o,
  `p-value`    = format(lrt_p_o, scientific = TRUE),
  Keputusan    = ifelse(lrt_p_o < 0.05,
                        "Tolak H0 - Model signifikan secara simultan",
                        "Gagal tolak H0")
) %>%
  kable(caption = "Likelihood Ratio Test: model penuh vs model null") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Likelihood Ratio Test: model penuh vs model null

Chi-square	df	p-value	Keputusan
694.2484	4	0e+00	Tolak H0 - Model signifikan secara simultan

Nilai p-value yang sangat kecil berarti kita menolak H0. Secara simultan, setidaknya satu prediktor berpengaruh signifikan terhadap variant kendaraan listrik.

3.11 Prediksi Probabilitas Ordinal

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

##         Base   Standard Long Range    Premium Performance
## 1 0.01461842 0.03430315 0.07621705 0.19256865  0.68229274
## 2 0.37628370 0.30027822 0.17674178 0.09653584  0.05016045
## 3 0.41310187 0.29624446 0.16222748 0.08512304  0.04330315
## 4 0.12431247 0.20554219 0.24798039 0.23887389  0.18329107
## 5 0.01826199 0.04232734 0.09148815 0.21655590  0.63136662
## 6 0.04386402 0.09337221 0.16944046 0.28347988  0.40984343

df_o_pred <- df_o_scaled %>%
  bind_cols(as.data.frame(pred_prob_ord)) %>%
  mutate(prediksi = predict(fit_ord, type = "class"))

head(df_o_pred %>% select(variant, prediksi, Base, Standard, `Long Range`, Premium, Performance))

Setiap baris menunjukkan probabilitas prediksi untuk kelima variant. Pada model ordinal, prediksi yang meleset ke kategori yang berdekatan secara konseptual lebih dapat diterima dibanding meleset jauh.

3.12 Confusion Matrix dan Akurasi Model

conf_ord <- table(Aktual = df_o_pred$variant, Prediksi = df_o_pred$prediksi)
conf_ord

##              Prediksi
## Aktual        Base Standard Long Range Premium Performance
##   Base         198      149         43      18          12
##   Standard     205      106         55      24          13
##   Long Range     1       53         85      95         173
##   Premium      126      109         67      48          42
##   Performance    0       11         44      73         250

accuracy_ord <- sum(diag(conf_ord)) / sum(conf_ord)
cat("\nAkurasi model:", round(accuracy_ord, 4), "\n")

## 
## Akurasi model: 0.3435

as.data.frame(conf_ord) %>%
  ggplot(aes(x = Aktual, y = Prediksi, fill = Freq)) +
  geom_tile(color = "white", linewidth = 0.8) +
  geom_text(aes(label = Freq), size = 5, fontface = "bold", color = "white") +
  scale_fill_gradient(low = "#AED6F1", high = "#1A5276") +
  labs(title    = "Confusion Matrix - Regresi Logistik Ordinal",
       subtitle = paste0("Akurasi: ", percent(accuracy_ord, accuracy = 0.1)),
       x = "Aktual", y = "Prediksi") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"), legend.position = "none",
        axis.text.x = element_text(angle = 20, hjust = 1))

Untuk model ordinal, perlu diperhatikan apakah kesalahan cenderung terjadi antar kategori yang berdekatan (kesalahan ringan) atau berjauhan (kesalahan berat).

3.13 Visualisasi Probabilitas Prediksi

3.13.1 Probabilitas terhadap Harga

price_mean_o <- mean(df_o$price_usd)
price_sd_o   <- sd(df_o$price_usd)

grid_price_o <- expand.grid(
  battery_capacity_kwh = 0, range_miles = 0, charging_speed_kw = 0,
  price_usd = seq(-2.5, 2.5, length.out = 150)
)

grid_prob_price_o <- predict(fit_ord, newdata = grid_price_o, type = "probs")

grid_price_o %>%
  bind_cols(as.data.frame(grid_prob_price_o)) %>%
  mutate(price_usd_orig = price_usd * price_sd_o + price_mean_o) %>%
  tidyr::pivot_longer(cols = c("Base", "Standard", "Long Range", "Premium", "Performance"),
                      names_to = "Variant", values_to = "Probabilitas") %>%
  mutate(Variant = ordered(Variant, levels = c("Base", "Standard", "Long Range",
                                               "Premium", "Performance"))) %>%
  ggplot(aes(x = price_usd_orig, y = Probabilitas, color = Variant)) +
  geom_line(linewidth = 1.2) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
  scale_x_continuous(labels = label_dollar(scale = 1e-3, suffix = "K")) +
  scale_color_manual(values = c("Base" = "#2f7f73", "Standard" = "#5568b8",
                                "Long Range" = "#d18b2f", "Premium" = "#b84f5a",
                                "Performance" = "#243142")) +
  labs(title    = "Prediksi Probabilitas Variant terhadap Harga (price_usd)",
       subtitle = "Prediktor lain di-hold pada nilai rata-rata",
       x = "Harga (USD)", y = "Probabilitas Prediksi", color = "Variant") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"), panel.grid.minor = element_blank())

Jika garis Performance naik seiring kenaikan harga sementara garis Base menurun, kendaraan yang lebih mahal secara khas tergolong variant yang lebih tinggi, yang sejalan dengan logika pasar EV.

3.13.2 Probabilitas terhadap Kapasitas Baterai

bat_mean_o <- mean(df_o$battery_capacity_kwh)
bat_sd_o   <- sd(df_o$battery_capacity_kwh)

grid_bat_o <- expand.grid(
  battery_capacity_kwh = seq(-2.5, 2.5, length.out = 150),
  range_miles = 0, charging_speed_kw = 0, price_usd = 0
)

grid_prob_bat_o <- predict(fit_ord, newdata = grid_bat_o, type = "probs")

grid_bat_o %>%
  bind_cols(as.data.frame(grid_prob_bat_o)) %>%
  mutate(bat_orig = battery_capacity_kwh * bat_sd_o + bat_mean_o) %>%
  tidyr::pivot_longer(cols = c("Base", "Standard", "Long Range", "Premium", "Performance"),
                      names_to = "Variant", values_to = "Probabilitas") %>%
  mutate(Variant = ordered(Variant, levels = c("Base", "Standard", "Long Range",
                                               "Premium", "Performance"))) %>%
  ggplot(aes(x = bat_orig, y = Probabilitas, color = Variant)) +
  geom_line(linewidth = 1.2) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
  scale_color_manual(values = c("Base" = "#2f7f73", "Standard" = "#5568b8",
                                "Long Range" = "#d18b2f", "Premium" = "#b84f5a",
                                "Performance" = "#243142")) +
  labs(title    = "Prediksi Probabilitas Variant terhadap Kapasitas Baterai",
       subtitle = "Prediktor lain di-hold pada nilai rata-rata",
       x = "Battery Capacity (kWh)", y = "Probabilitas Prediksi", color = "Variant") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"), panel.grid.minor = element_blank())

Kendaraan dengan baterai lebih besar umumnya merupakan variant yang lebih canggih, sehingga garis Performance diharapkan mendominasi pada nilai kWh yang tinggi.

3.14 Uji Asumsi Proportional Odds

3.14.1 Grafik Parallel Lines

Asumsi proportional odds mensyaratkan bahwa efek setiap prediktor sama di semua cutpoint, yang secara visual terlihat sebagai garis-garis cumulative logit yang paralel.

eps <- 1e-6

grid_bat_o %>%
  bind_cols(as.data.frame(grid_prob_bat_o)) %>%
  mutate(bat_orig = battery_capacity_kwh * bat_sd_o + bat_mean_o) %>%
  mutate(
    `P(Y <= Base)`       = Base,
    `P(Y <= Standard)`   = Base + Standard,
    `P(Y <= Long Range)` = Base + Standard + `Long Range`,
    `P(Y <= Premium)`    = Base + Standard + `Long Range` + Premium
  ) %>%
  tidyr::pivot_longer(cols = c(`P(Y <= Base)`, `P(Y <= Standard)`,
                               `P(Y <= Long Range)`, `P(Y <= Premium)`),
                      names_to = "batas_kumulatif", values_to = "prob_kumulatif") %>%
  mutate(
    prob_kumulatif   = pmin(pmax(prob_kumulatif, eps), 1 - eps),
    cumulative_logit = qlogis(prob_kumulatif)
  ) %>%
  ggplot(aes(x = bat_orig, y = cumulative_logit, color = batas_kumulatif)) +
  geom_line(linewidth = 1.2) +
  scale_color_manual(values = c("P(Y <= Base)" = "#2f7f73", "P(Y <= Standard)" = "#5568b8",
                                "P(Y <= Long Range)" = "#d18b2f", "P(Y <= Premium)" = "#b84f5a")) +
  labs(title    = "Parallel Lines - Model Ordinal (Cumulative Logit)",
       subtitle = "Jika garis sejajar (paralel), asumsi proportional odds terpenuhi.",
       x = "Battery Capacity (kWh)", y = "Logit peluang kumulatif",
       color = "Batas kumulatif") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"), panel.grid.minor = element_blank())

Jika keempat garis berjalan secara paralel, asumsi proportional odds terpenuhi secara visual. Garis yang menyilang atau memiliki kemiringan yang sangat berbeda menjadi indikasi pelanggaran asumsi.

3.14.2 Uji Formal Proportional Odds

\[H_0: \text{Efek prediktor sama di semua batas kumulatif}\] \[H_1: \text{Efek prediktor berbeda antar cutpoint}\]

fit_clm <- ordinal::clm(
  variant ~ battery_capacity_kwh + range_miles + charging_speed_kw + price_usd,
  data = df_o_scaled, link = "logit"
)

summary(fit_clm)

## formula: 
## variant ~ battery_capacity_kwh + range_miles + charging_speed_kw + price_usd
## data:    df_o_scaled
## 
##  link  threshold nobs logLik   AIC     niter max.grad cond.H 
##  logit flexible  2000 -2870.49 5756.98 4(0)  6.20e-07 9.2e+01
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## battery_capacity_kwh  0.95274    0.19661   4.846 1.26e-06 ***
## range_miles          -0.11408    0.19487  -0.585    0.558    
## charging_speed_kw     0.25927    0.04410   5.880 4.11e-09 ***
## price_usd             0.54680    0.04564  11.980  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Threshold coefficients:
##                     Estimate Std. Error z value
## Base|Standard       -1.68678    0.06201 -27.202
## Standard|Long Range -0.44339    0.05079  -8.729
## Long Range|Premium   0.57933    0.05127  11.300
## Premium|Performance  1.75963    0.06452  27.271

cat("\n--- Uji Nominal (Proportional Odds) ---\n")

## 
## --- Uji Nominal (Proportional Odds) ---

ordinal::nominal_test(fit_clm)

P-value yang besar (> 0,05) untuk suatu variabel menunjukkan tidak ada bukti pelanggaran asumsi proportional odds. Jika p-value < 0,05, efek prediktor tersebut berbeda antar cutpoint dan model perlu dieksplorasi lebih lanjut, misalnya dengan partial proportional odds model.

3.15 Ringkasan dan Kesimpulan

Ringkasan hasil analisis regresi logistik ordinal

Aspek	Nilai
Variabel Dependen	variant (5 kategori) - ordinal
Urutan Ordinal	Base < Standard < Long Range < Premium < Performance
Jumlah Observasi	2000
Jumlah Prediktor	4 prediktor numerik
Metode Estimasi	MASS::polr() - logistic
Standarisasi Prediktor	Ya: mean = 0, sd = 1 (wajib untuk stabilitas numerik)
Uji Simultan (LRT)	Chi-square = 694.25, p < 0.05, Signifikan
Akurasi Klasifikasi	34.4%
Asumsi Proportional Odds	Diperiksa via grafik parallel lines dan nominal_test()

Beberapa poin utama dari hasil analisis ini:

1. Model signifikan secara simultan berdasarkan Likelihood Ratio Test.

2. Koefisien polr() memiliki tanda yang berlawanan dengan konvensi Agresti/slide. Gunakan estimate_slide = -estimate_polr untuk interpretasi yang konsisten dengan materi kuliah.

3. Standarisasi prediktor wajib dilakukan sebelum polr() karena perbedaan skala ekstrem (terutama price_usd) dapat menyebabkan kegagalan komputasi numerik.

4. Asumsi proportional odds merupakan asumsi kunci model ordinal dan perlu diperiksa, baik secara visual maupun melalui uji formal.

---

4 Regresi Poisson

4.1 Pendahuluan

Regresi Poisson digunakan untuk variabel respons berupa data hitung, yaitu nilai bilangan bulat nonnegatif. Berbeda dari regresi linear yang mengasumsikan respons kontinu dan berdistribusi normal, regresi Poisson memodelkan rata-rata hitungan menggunakan fungsi penghubung log:

\[\log(\mu_i) = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip}\]

sehingga prediksi \(\hat{\mu}_i = \exp(\hat{\eta}_i)\) selalu bernilai positif.

Tujuan analisis ini adalah memodelkan jumlah penyewaan sepeda per jam di Seoul sebagai fungsi dari kondisi cuaca menggunakan regresi Poisson.

Variabel yang digunakan:

Tipe	Variabel
Dependen (Y)	Rented.Bike.Count - jumlah sepeda yang disewa per jam
Independen Numerik	Temperature - suhu udara (°C)
	Humidity - kelembapan udara (%)
	Wind.speed - kecepatan angin (m/s)
	Rainfall - curah hujan (mm)

4.2 Memuat Paket

Bagian ini menggunakan paket yang sudah dimuat sebelumnya (readxl, tidyverse, knitr, kableExtra, scales).

4.3 Memuat dan Memeriksa Data

df_p_raw <- readxl::read_excel("Data Poisson.xls")

df_p <- df_p_raw %>%
  rename(
    Rented.Bike.Count = `Rented Bike Count`,
    Temperature       = `Temperature(C)`,
    Humidity          = `Humidity(%)`,
    Wind.speed        = `Wind speed (m/s)`,
    Rainfall          = `Rainfall(mm)`
  ) %>%
  select(Rented.Bike.Count, Temperature, Humidity, Wind.speed, Rainfall)

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

## Jumlah baris   : 8760

cat("Jumlah kolom   :", ncol(df_p), "\n")

## Jumlah kolom   : 5

Dataset terdiri dari 8.760 observasi (satu per jam selama satu tahun) dengan 5 variabel yang digunakan.

head(df_p, 6)

str(df_p)

## tibble [8,760 × 5] (S3: tbl_df/tbl/data.frame)
##  $ Rented.Bike.Count: num [1:8760] 254 204 173 107 78 100 181 460 930 490 ...
##  $ Temperature      : num [1:8760] -5.2 -5.5 -6 -6.2 -6 -6.4 -6.6 -7.4 -7.6 -6.5 ...
##  $ Humidity         : num [1:8760] 37 38 39 40 36 37 35 38 37 27 ...
##  $ Wind.speed       : num [1:8760] 2.2 0.8 1 0.9 2.3 1.5 1.3 0.9 1.1 0.5 ...
##  $ Rainfall         : num [1:8760] 0 0 0 0 0 0 0 0 0 0 ...

Variabel respons Rented.Bike.Count bertipe integer (bilangan bulat nonnegatif), yang merupakan syarat utama untuk regresi Poisson. Semua prediktor bertipe numerik kontinu dan tidak ditemukan missing value.

4.4 Analisis Deskriptif

4.4.1 Statistik Ringkas

df_p %>%
  summarise(across(
    everything(),
    list(Min = ~min(.x, na.rm = TRUE), Mean = ~round(mean(.x, na.rm = TRUE), 2),
         Median = ~median(.x, na.rm = TRUE), Max = ~max(.x, na.rm = TRUE),
         SD = ~round(sd(.x, na.rm = TRUE), 2)),
    .names = "{.col}__{.fn}"
  )) %>%
  tidyr::pivot_longer(everything(), names_to = c("Variabel", "Statistik"), names_sep = "__") %>%
  tidyr::pivot_wider(names_from = Statistik, values_from = value) %>%
  kable(caption = "Statistik deskriptif variabel penelitian", digits = 2) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE)

Statistik deskriptif variabel penelitian

Variabel	Min	Mean	Median	Max	SD
Rented.Bike.Count	0.0	704.60	504.5	3556.0	645.00
Temperature	-17.8	12.88	13.7	39.4	11.94
Humidity	0.0	58.23	57.0	98.0	20.36
Wind.speed	0.0	1.72	1.5	7.4	1.04
Rainfall	0.0	0.15	0.0	35.0	1.13

Rata-rata jumlah penyewaan sepeda per jam sekitar 705 unit dengan standar deviasi yang cukup besar, mengindikasikan variasi yang tinggi. Curah hujan memiliki nilai rata-rata sangat rendah karena mayoritas jam tidak hujan sehingga distribusinya sangat miring ke kanan.

4.4.2 Distribusi Variabel Respons

df_p %>%
  ggplot(aes(x = Rented.Bike.Count)) +
  geom_histogram(bins = 50, fill = "#2f7f73", alpha = 0.85, color = "white") +
  scale_x_continuous(labels = comma_format()) +
  scale_y_continuous(labels = comma_format()) +
  labs(title    = "Distribusi Jumlah Penyewaan Sepeda per Jam",
       subtitle = "Respons berupa data hitung (count data) - bernilai bilangan bulat nonnegatif",
       x = "Jumlah sepeda yang disewa", y = "Frekuensi") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))

Distribusi respons sangat miring ke kanan. Sebagian besar jam memiliki jumlah penyewaan rendah, sementara sebagian kecil jam memiliki penyewaan sangat tinggi. Pola ini konsisten dengan karakteristik data hitung dan menunjukkan regresi Poisson lebih tepat dibanding regresi linear biasa.

4.4.3 Hubungan Prediktor dengan Respons

predictors_p  <- c("Temperature", "Humidity", "Wind.speed", "Rainfall")
pred_labels_p <- c("Suhu (°C)", "Kelembapan (%)", "Kecepatan Angin (m/s)", "Curah Hujan (mm)")

par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
for (i in seq_along(predictors_p)) {
  plot(df_p[[predictors_p[i]]], df_p$Rented.Bike.Count,
       xlab = pred_labels_p[i], ylab = "Rented Bike Count",
       main = paste("Rented Bike Count vs", pred_labels_p[i]),
       pch = 16, col = adjustcolor("#2f7f73", alpha.f = 0.15), cex = 0.6)
  abline(lm(df_p$Rented.Bike.Count ~ df_p[[predictors_p[i]]]),
         col = "#b84f5a", lwd = 2, lty = 2)
}

par(mfrow = c(1, 1))

Suhu memiliki hubungan positif dengan penyewaan. Kelembapan dan curah hujan cenderung negatif. Kecepatan angin menunjukkan pola lemah namun sedikit negatif.

4.5 Pemeriksaan Asumsi

Daftar asumsi regresi Poisson dan statusnya

Asumsi	Status
Respons berupa hitungan nonnegatif	Terpenuhi: Y adalah jumlah sepeda (0, 1, 2, …)
Observasi independen	Diasumsikan terpenuhi: setiap jam berdiri sendiri
Log-link: hubungan linear pada skala log	Akan diperiksa melalui residual plot
Equidispersion - E(Y) = Var(Y)	Perlu diperiksa: data penyewaan sering overdispersi
Tidak ada nilai negatif pada Y	Terpenuhi: tidak ada jumlah negatif

4.5.1 Indikasi Awal: Perbandingan Mean vs Varians

Pada distribusi Poisson yang ideal berlaku equidispersion, yaitu E(Y) = Var(Y), sehingga rasio varians terhadap mean seharusnya mendekati 1.

tibble(
  Statistik = c("Mean Y", "Varians Y", "Rasio Varians/Mean"),
  Nilai     = c(round(mean(df_p$Rented.Bike.Count), 2),
                round(var(df_p$Rented.Bike.Count), 2),
                round(var(df_p$Rented.Bike.Count) / mean(df_p$Rented.Bike.Count), 2))
) %>%
  kable(caption = "Perbandingan mean dan varians variabel respons") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Perbandingan mean dan varians variabel respons

Statistik	Nilai
Mean Y	704.60
Varians Y	416021.73
Rasio Varians/Mean	590.44

Varians jauh lebih besar dari mean (rasio >> 1), mengindikasikan overdispersion. Ini umum terjadi pada data penyewaan yang dipengaruhi banyak faktor tidak teramati. Model Poisson tetap dapat diestimasi, namun overdispersion akan dikonfirmasi kembali setelah model difit.

4.6 Estimasi Model Regresi Poisson

fit_pois <- glm(
  Rented.Bike.Count ~ Temperature + Humidity + Wind.speed + Rainfall,
  data   = df_p,
  family = poisson(link = "log")
)

summary(fit_pois)

## 
## Call:
## glm(formula = Rented.Bike.Count ~ Temperature + Humidity + Wind.speed + 
##     Rainfall, family = poisson(link = "log"), data = df_p)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  6.302e+00  1.822e-03  3458.1   <2e-16 ***
## Temperature  4.545e-02  3.636e-05  1249.9   <2e-16 ***
## Humidity    -9.125e-03  2.328e-05  -392.0   <2e-16 ***
## Wind.speed   4.230e-02  4.211e-04   100.4   <2e-16 ***
## Rainfall    -4.934e-01  2.129e-03  -231.8   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 4979261  on 8759  degrees of freedom
## Residual deviance: 2917946  on 8755  degrees of freedom
## AIC: 2985050
## 
## Number of Fisher Scoring iterations: 6

Tanda positif pada koefisien berarti prediktor meningkatkan log rata-rata penyewaan; tanda negatif berarti menguranginya. Selisih antara Null deviance dan Residual deviance menunjukkan seberapa besar variasi yang dapat dijelaskan model.

4.7 Ringkasan Koefisien, IRR, dan Interval Kepercayaan

Koefisien dieksponensialkan untuk memperoleh Incidence Rate Ratio (IRR), yaitu faktor perkalian rata-rata jumlah penyewaan untuk setiap kenaikan satu unit prediktor:

\[IRR_j = \exp(\hat{\beta}_j)\]

Persentase perubahan rata-rata: \(100 \times \{\exp(\hat{\beta}_j) - 1\}\%\)

pois_coef <- as.data.frame(coef(summary(fit_pois))) %>%
  tibble::rownames_to_column("Parameter") %>%
  dplyr::rename(Estimate = Estimate, Std.Error = `Std. Error`,
                z.value = `z value`, p.value = `Pr(>|z|)`) %>%
  mutate(
    IRR           = exp(Estimate),
    CI_low        = exp(Estimate - 1.96 * Std.Error),
    CI_high       = exp(Estimate + 1.96 * Std.Error),
    Perubahan_pct = round(100 * (IRR - 1), 4),
    Signifikan    = case_when(
      p.value < 0.001 ~ "***", p.value < 0.01 ~ "**",
      p.value < 0.05  ~ "*",   p.value < 0.10 ~ ".",
      TRUE            ~ "ns"
    )
  ) %>%
  mutate(across(c(Estimate, Std.Error, z.value, p.value, IRR, CI_low, CI_high), ~round(.x, 4)))

pois_coef %>%
  kable(caption = "Ringkasan hasil regresi Poisson - Koefisien, IRR, dan Interval Kepercayaan 95%",
        col.names = c("Parameter", "β", "SE", "z-value", "p-value",
                      "IRR", "CI 2,5%", "CI 97,5%", "Perubahan (%)", "Sig.")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = TRUE) %>%
  footnote(general = "Kode sig.: *** p<0.001 | ** p<0.01 | * p<0.05 | . p<0.10 | ns tidak signifikan")

Ringkasan hasil regresi Poisson - Koefisien, IRR, dan Interval Kepercayaan 95%

Parameter	β	SE	z-value	p-value	IRR	CI 2,5%	CI 97,5%	Perubahan (%)	Sig.
(Intercept)	6.3017	0.0018	3458.1326	0	545.4980	543.5531	547.4498	54449.7960	***
Temperature	0.0454	0.0000	1249.8976	0	1.0465	1.0464	1.0466	4.6496	***
Humidity	-0.0091	0.0000	-392.0145	0	0.9909	0.9909	0.9910	-0.9084	***
Wind.speed	0.0423	0.0004	100.4478	0	1.0432	1.0423	1.0441	4.3210	***
Rainfall	-0.4934	0.0021	-231.7566	0	0.6106	0.6080	0.6131	-38.9446	***
Note:
Kode sig.: *** p<0.001 | ** p<0.01 | * p<0.05 | . p<0.10 | ns tidak signifikan

4.8 Interpretasi Koefisien

4.8.1 Intercept

Ketika semua prediktor bernilai nol, log rata-rata jumlah penyewaan sepeda adalah sebesar \(\hat{\beta}_0\). Nilai ini kurang bermakna secara praktis karena kondisi tersebut tidak realistis, namun diperlukan sebagai konstanta dasar model.

4.8.2 Suhu (Temperature)

Koefisien Temperature bernilai positif dan signifikan (p < 0,001). Setiap kenaikan suhu sebesar 1°C, dengan prediktor lain konstan, dikaitkan dengan rata-rata jumlah penyewaan yang naik sekitar 4.65% (IRR = 1.0465). Cuaca yang lebih hangat mendorong lebih banyak orang untuk bersepeda.

4.8.3 Kelembapan (Humidity)

Koefisien Humidity bernilai negatif dan signifikan (p < 0,001). Setiap kenaikan kelembapan sebesar 1%, dengan prediktor lain konstan, dikaitkan dengan rata-rata jumlah penyewaan yang turun sekitar 0.91% (IRR = 0.9909). Udara yang lebih lembap terasa tidak nyaman untuk bersepeda.

4.8.4 Kecepatan Angin (Wind.speed)

Koefisien Wind.speed bernilai negatif dan signifikan (p < 0,001). Setiap kenaikan kecepatan angin sebesar 1 m/s, dengan prediktor lain konstan, dikaitkan dengan rata-rata jumlah penyewaan yang turun sekitar 4.32% (IRR = 1.0432). Angin kencang mempersulit dan mengurangi kenyamanan bersepeda.

4.8.5 Curah Hujan (Rainfall)

Koefisien Rainfall bernilai negatif dan signifikan (p < 0,001). Setiap kenaikan curah hujan sebesar 1 mm, dengan prediktor lain konstan, dikaitkan dengan rata-rata jumlah penyewaan yang turun sekitar 38.94% (IRR = 0.6106). Hujan merupakan hambatan terbesar bagi pengguna sepeda.

4.9 Visualisasi Prediksi Rate Poisson

4.9.1 Prediksi Berdasarkan Suhu

grid_temp_p <- data.frame(
  Temperature = seq(min(df_p$Temperature), max(df_p$Temperature), length.out = 200),
  Humidity    = mean(df_p$Humidity), Wind.speed = mean(df_p$Wind.speed), Rainfall = 0
)
pred_temp_p <- predict(fit_pois, newdata = grid_temp_p, type = "link", se.fit = TRUE)

grid_temp_p %>%
  mutate(fit_link = pred_temp_p$fit, se_link = pred_temp_p$se.fit,
         rate = exp(fit_link), rate_low = exp(fit_link - 1.96 * se_link),
         rate_high = exp(fit_link + 1.96 * se_link)) %>%
  ggplot(aes(x = Temperature, y = rate)) +
  geom_ribbon(aes(ymin = rate_low, ymax = rate_high), fill = "#2f7f73", alpha = 0.20) +
  geom_line(color = "#2f7f73", linewidth = 1.4) +
  scale_y_continuous(labels = comma_format()) +
  labs(title    = "Prediksi Rata-Rata Penyewaan Sepeda terhadap Suhu",
       subtitle = "Kelembapan dan kecepatan angin ditahan pada rata-rata; curah hujan = 0 mm",
       x = "Suhu (°C)", y = "Prediksi jumlah penyewaan per jam") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))

Kurva prediksi memperlihatkan kenaikan eksponensial rata-rata penyewaan seiring meningkatnya suhu. Pita adalah interval kepercayaan 95%. Pola ini mengkonfirmasi bahwa cuaca hangat secara signifikan mendorong permintaan penyewaan sepeda.

4.9.2 Prediksi Berdasarkan Kelembapan

grid_humid_p <- data.frame(
  Temperature = mean(df_p$Temperature),
  Humidity    = seq(min(df_p$Humidity), max(df_p$Humidity), length.out = 200),
  Wind.speed  = mean(df_p$Wind.speed), Rainfall = 0
)
pred_humid_p <- predict(fit_pois, newdata = grid_humid_p, type = "link", se.fit = TRUE)

grid_humid_p %>%
  mutate(fit_link = pred_humid_p$fit, se_link = pred_humid_p$se.fit,
         rate = exp(fit_link), rate_low = exp(fit_link - 1.96 * se_link),
         rate_high = exp(fit_link + 1.96 * se_link)) %>%
  ggplot(aes(x = Humidity, y = rate)) +
  geom_ribbon(aes(ymin = rate_low, ymax = rate_high), fill = "#5568b8", alpha = 0.20) +
  geom_line(color = "#5568b8", linewidth = 1.4) +
  scale_y_continuous(labels = comma_format()) +
  labs(title    = "Prediksi Rata-Rata Penyewaan Sepeda terhadap Kelembapan",
       subtitle = "Suhu dan kecepatan angin ditahan pada rata-rata; curah hujan = 0 mm",
       x = "Kelembapan (%)", y = "Prediksi jumlah penyewaan per jam") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))

Kurva menurun secara eksponensial dari kelembapan rendah ke tinggi, mengkonfirmasi efek negatif udara lembap terhadap permintaan penyewaan sepeda.

4.9.3 Prediksi Berdasarkan Curah Hujan

grid_rain_p <- data.frame(
  Temperature = mean(df_p$Temperature), Humidity = mean(df_p$Humidity),
  Wind.speed  = mean(df_p$Wind.speed), Rainfall = seq(0, 10, length.out = 200)
)
pred_rain_p <- predict(fit_pois, newdata = grid_rain_p, type = "link", se.fit = TRUE)

grid_rain_p %>%
  mutate(fit_link = pred_rain_p$fit, se_link = pred_rain_p$se.fit,
         rate = exp(fit_link), rate_low = exp(fit_link - 1.96 * se_link),
         rate_high = exp(fit_link + 1.96 * se_link)) %>%
  ggplot(aes(x = Rainfall, y = rate)) +
  geom_ribbon(aes(ymin = rate_low, ymax = rate_high), fill = "#d18b2f", alpha = 0.20) +
  geom_line(color = "#d18b2f", linewidth = 1.4) +
  scale_y_continuous(labels = comma_format()) +
  labs(title    = "Prediksi Rata-Rata Penyewaan Sepeda terhadap Curah Hujan",
       subtitle = "Suhu, kelembapan, dan kecepatan angin ditahan pada rata-rata",
       x = "Curah hujan (mm)", y = "Prediksi jumlah penyewaan per jam") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))

Penurunan paling dramatis terlihat di sini. Bahkan curah hujan kecil (1-2 mm) sudah mengurangi prediksi rata-rata penyewaan secara signifikan, konsisten dengan perilaku nyata pengguna sepeda yang sangat sensitif terhadap hujan.

4.10 Uji Simultan

\[H_0: \beta_1 = \beta_2 = \beta_3 = \beta_4 = 0\] \[H_1: \text{minimal satu } \beta_j \neq 0\]

model_null_p <- glm(Rented.Bike.Count ~ 1, data = df_p, family = poisson(link = "log"))

lrt_stat_p <- 2 * (as.numeric(logLik(fit_pois)) - as.numeric(logLik(model_null_p)))
lrt_df_p   <- length(coef(fit_pois)) - length(coef(model_null_p))
lrt_p_p    <- 1 - pchisq(lrt_stat_p, lrt_df_p)

tibble(
  `Chi-square` = round(lrt_stat_p, 4),
  df           = lrt_df_p,
  `p-value`    = format(lrt_p_p, scientific = TRUE),
  Keputusan    = ifelse(lrt_p_p < 0.05,
                        "Tolak H0 - Model signifikan secara simultan",
                        "Gagal tolak H0")
) %>%
  kable(caption = "Likelihood Ratio Test: model penuh vs model null") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Likelihood Ratio Test: model penuh vs model null

Chi-square	df	p-value	Keputusan
2061315	4	0e+00	Tolak H0 - Model signifikan secara simultan

Nilai G² yang sangat besar dengan p-value mendekati 0 berarti kita menolak H0. Secara simultan, minimal satu prediktor cuaca berpengaruh signifikan terhadap rata-rata jumlah penyewaan sepeda.

4.11 Pemeriksaan Overdispersion

Overdispersion terjadi ketika \(\hat{\phi} = \sum r_{P,i}^2 / df \gg 1\), di mana \(r_{P,i}\) adalah residual Pearson dan \(df\) adalah derajat bebas residual.

disp_pois <- sum(residuals(fit_pois, type = "pearson")^2) / df.residual(fit_pois)

tibble(
  `Dispersion Pearson` = round(disp_pois, 3),
  Interpretasi = dplyr::case_when(
    disp_pois < 1.5 ~ "Tidak ada indikasi overdispersion berat",
    disp_pois < 2.5 ~ "Ada indikasi overdispersion sedang",
    TRUE            ~ "Ada indikasi overdispersion kuat - pertimbangkan quasi-Poisson / Negative Binomial"
  )
) %>%
  kable(caption = "Pemeriksaan overdispersion pada model Poisson") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Pemeriksaan overdispersion pada model Poisson

Dispersion Pearson	Interpretasi
131030.2	Ada indikasi overdispersion kuat - pertimbangkan quasi-Poisson / Negative Binomial

Nilai Dispersion Pearson yang jauh di atas 1 mengkonfirmasi adanya overdispersion kuat. Ini berarti varians aktual data jauh melebihi asumsi distribusi Poisson. Akibatnya, standard error koefisien dari model Poisson standar cenderung terlalu kecil sehingga uji signifikansi individual perlu dicermati. Pertimbangkan penggunaan quasi-Poisson atau Negative Binomial regression sebagai langkah lanjutan.

4.12 Prediksi Model

df_pred_p <- df_p %>%
  mutate(
    fitted_log  = predict(fit_pois, type = "link"),
    fitted_mean = predict(fit_pois, type = "response")
  )

df_pred_p %>%
  select(Rented.Bike.Count, Temperature, Humidity, Wind.speed, Rainfall,
         fitted_log, fitted_mean) %>%
  head(10) %>%
  mutate(across(c(fitted_log, fitted_mean), ~round(.x, 2))) %>%
  kable(caption = "Sepuluh baris pertama: nilai aktual vs prediksi model Poisson",
        col.names = c("Y Aktual", "Suhu", "Kelembapan", "Angin", "Hujan",
                      "η̂ (log scale)", "μ̂ (skala asli)")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = TRUE)

Sepuluh baris pertama: nilai aktual vs prediksi model Poisson

Y Aktual	Suhu	Kelembapan	Angin	Hujan	η̂ (log scale	μ̂ (skala asl
254	-5.2	37	2.2	0	5.82	337.24
204	-5.5	38	0.8	0	5.74	310.70
173	-6.0	39	1.0	0	5.72	303.52
107	-6.2	40	0.9	0	5.69	296.78
78	-6.0	36	2.3	0	5.80	329.58
100	-6.4	37	1.5	0	5.74	310.03
181	-6.6	35	1.3	0	5.74	310.24
460	-7.4	38	0.9	0	5.66	286.20
930	-7.6	37	1.1	0	5.67	288.65
490	-6.5	27	0.5	0	5.78	324.10

Kolom η̂ adalah prediksi pada skala log-link, sedangkan μ̂ = exp(η̂) adalah prediksi rata-rata jumlah penyewaan pada skala asli.

4.12.1 Perbandingan Nilai Aktual vs Prediksi

df_pred_p %>%
  slice_sample(n = 1000) %>%
  ggplot(aes(x = fitted_mean, y = Rented.Bike.Count)) +
  geom_point(alpha = 0.25, color = "#2f7f73", size = 1.2) +
  geom_abline(slope = 1, intercept = 0, color = "#b84f5a", linewidth = 1.2, linetype = "dashed") +
  scale_x_continuous(labels = comma_format()) +
  scale_y_continuous(labels = comma_format()) +
  labs(title    = "Nilai Aktual vs Prediksi Model Poisson",
       subtitle = "Garis merah putus-putus = garis sempurna (aktual = prediksi); sampel 1.000 observasi",
       x = "Prediksi (μ̂)", y = "Aktual (Y)") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))

Jika model sempurna, semua titik akan berada tepat di garis merah. Sebaran titik yang lebar di sekitar garis menunjukkan variasi yang belum tertangkap model, konsisten dengan overdispersion yang terdeteksi sebelumnya.

4.13 Ringkasan dan Kesimpulan

Ringkasan hasil analisis regresi Poisson

Aspek	Nilai
Variabel Dependen	Rented Bike Count (jumlah penyewaan per jam)
Tipe Data Respons	Count data - bilangan bulat nonnegatif (0, 1, 2, …)
Distribusi Asumsi	Poisson
Fungsi Penghubung	log
Jumlah Observasi	8760
Jumlah Prediktor	4 prediktor numerik
Null Deviance	4979261.4
Residual Deviance	2917946.19
AIC	2985049.92
Uji Simultan (LRT)	G² = 2061315.21, p ≈ 0, Tolak H0, model signifikan
Overdispersion	Dispersion Pearson = 131030.22, Overdispersion kuat

Beberapa poin utama dari hasil analisis ini:

1. Model Poisson signifikan secara simultan. Secara bersama-sama, kondisi cuaca berpengaruh nyata terhadap rata-rata jumlah penyewaan sepeda per jam.

2. Suhu berpengaruh positif. Setiap kenaikan 1°C dikaitkan dengan kenaikan rata-rata penyewaan sekitar 4.65%. Cuaca hangat mendorong lebih banyak orang bersepeda.

3. Kelembapan, kecepatan angin, dan curah hujan berpengaruh negatif. Curah hujan memiliki efek paling dramatis di antara ketiganya.

4. Overdispersion terdeteksi kuat. Varians data jauh melebihi mean, melanggar asumsi equidispersion Poisson. Penggunaan quasi-Poisson atau Negative Binomial direkomendasikan untuk inferensi yang lebih tepat.

5. Koefisien diinterpretasikan sebagai IRR. Nilai exp(β̂_j) memberikan faktor perkalian rata-rata hitungan untuk setiap kenaikan satu unit prediktor, dengan prediktor lain konstan.

---