Pemilihan jenis akomodasi oleh wisatawan milenial di Ubud dengan regresi multinomial

# Contoh kode R: Multinomial Logistic Regression
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(nnet)
})
set.seed(42)
# ---- Gaya grafik ----
pal <- c("#1f6feb", "#e8590c", "#2f9e44", "#ae3ec9", "#f59f00")
warna <- list(bg = "#f4f6fa", panel = "#ffffff", teks = "#1f2933", redup = "#5f6b78", grid = "#e3e8ee", bar = "#1d3557")
tata <- 1
theme_set(theme_minimal(base_size = 11) + theme(legend.position = "bottom"))

# ---- Data sintetis ----
n <- 200
d <- data.frame(travel_budget = rnorm(n, 50, 10),
                number_of_nights = rnorm(n, 30, 6))
z1 <- as.numeric(scale(d$travel_budget))
z2 <- as.numeric(scale(d$number_of_nights))
kat <- c("A", "B", "C")                     # kategori nominal, kategori pertama = referensi
# Utilitas tiap kategori relatif terhadap kategori referensi (utilitas 0)
u2 <- runif(1, -0.3, 0.3) + runif(1, 0.9, 1.6) * z1 + runif(1, -0.5, 0.2) * z2
u3 <- runif(1, -0.3, 0.3) + runif(1, -0.2, 0.6) * z1 + runif(1, 0.9, 1.6) * z2
pm <- cbind(1, exp(u2), exp(u3)); pm <- pm / rowSums(pm)
d$chosen_accommodation_type <- factor(kat[apply(pm, 1, function(p) sample(1:3, 1, prob = p))], levels = kat)
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Analisis: regresi logistik multinomial (nnet::multinom) ----
d$chosen_accommodation_type <- relevel(factor(d$chosen_accommodation_type, levels = kat), ref = kat[1])   # tetapkan kategori referensi
fit <- nnet::multinom(chosen_accommodation_type ~ travel_budget + number_of_nights, data = d, trace = FALSE, maxit = 500)
s <- summary(fit)
print(s)
zval <- s$coefficients / s$standard.errors
pval <- 2 * pnorm(abs(zval), lower.tail = FALSE)            # uji Wald tiap koefisien
print(round(pval, 4))
print(round(exp(coef(fit)), 4))                              # relative risk ratio (RRR) per satuan

# RRR per 1 simpangan baku untuk tiap kategori vs referensi
sdx <- c(sd(d$travel_budget), sd(d$number_of_nights))
rr <- do.call(rbind, lapply(1:2, function(j) {
  b <- s$coefficients[j, 2:3]; se <- s$standard.errors[j, 2:3]
  data.frame(kategori = kat[j + 1], var = c("travel_budget", "number_of_nights"), rrr = exp(b * sdx),
             lo = exp((b - 1.96 * se) * sdx), hi = exp((b + 1.96 * se) * sdx), p = pval[j, 2:3])
}))
print(rr)

# ---- Uji model ----
fit0 <- nnet::multinom(chosen_accommodation_type ~ 1, data = d, trace = FALSE)
lr <- anova(fit0, fit)                                        # uji likelihood ratio
print(lr)
p_lr <- lr[2, "Pr(Chi)"]
r2_mcf <- 1 - as.numeric(logLik(fit)) / as.numeric(logLik(fit0))
d$prediksi <- factor(predict(fit, type = "class"), levels = kat)
tab <- table(Aktual = d$chosen_accommodation_type, Prediksi = d$prediksi)
print(tab)
akurasi <- mean(d$prediksi == d$chosen_accommodation_type)
nir <- max(table(d$chosen_accommodation_type)) / n                               # no-information rate (kelas mayoritas)


# ---- Grafik ----
warna_kat <- setNames(pal[1:3], kat)
lab_pendek <- function(s) substr(s, 1, 16)
lab_lipat <- function(s) vapply(s, function(x) paste(strwrap(x, 11), collapse = "\n"), "")
# Peluang prediksi tiap kategori menurut X1 (X2 pada rata-rata)
g <- data.frame(travel_budget = seq(min(d$travel_budget), max(d$travel_budget), length.out = 120), number_of_nights = mean(d$number_of_nights))
pg <- predict(fit, newdata = g, type = "probs")
dg <- data.frame(x = rep(g$travel_budget, 3), kategori = factor(rep(kat, each = nrow(g)), levels = kat), p = as.numeric(pg))
p1 <- ggplot(dg, aes(x, p, colour = kategori)) +
  geom_line(linewidth = 1.2) +
  scale_colour_manual(values = warna_kat, name = NULL, labels = lab_pendek) +
  scale_y_continuous(labels = scales::percent, limits = c(0, 1)) +
  labs(title = "Predicted probability by category", x = "travel_budget", y = "Probability")
# Wilayah prediksi pada bidang X1-X2
gr <- expand.grid(travel_budget = seq(min(d$travel_budget), max(d$travel_budget), length.out = 90),
                  number_of_nights = seq(min(d$number_of_nights), max(d$number_of_nights), length.out = 90))
gr$kelas <- factor(predict(fit, newdata = gr, type = "class"), levels = kat)
p2 <- ggplot() +
  geom_raster(data = gr, aes(travel_budget, number_of_nights, fill = kelas), alpha = 0.22, show.legend = FALSE) +
  geom_point(data = d, aes(travel_budget, number_of_nights, colour = chosen_accommodation_type), size = 1.5, alpha = 0.8) +
  scale_fill_manual(values = warna_kat) +
  scale_colour_manual(values = warna_kat, name = NULL, labels = lab_pendek) +
  labs(title = "Observed classes and predicted regions", x = "travel_budget", y = "number_of_nights") +
  theme(legend.position = "none")
# Matriks konfusi
dt <- as.data.frame(tab)
p3 <- ggplot(dt, aes(Prediksi, Aktual)) +
  geom_tile(aes(alpha = Freq), fill = pal[1], colour = warna$panel, linewidth = 1) +
  geom_label(aes(label = Freq), fill = warna$panel, colour = warna$teks, label.size = 0, size = 3.3) +
  scale_alpha(range = c(0.08, 0.75), guide = "none") +
  scale_x_discrete(labels = lab_lipat) +
  scale_y_discrete(labels = lab_lipat, limits = rev(kat)) +
  labs(title = sprintf("Confusion matrix (accuracy %.1f%%)", 100 * mean(d$prediksi == d$chosen_accommodation_type)),
       x = "Predicted", y = "Observed")
p_hasil <- if (tata == 1) p1 / (p2 | p3) else p1 | (p2 / p3)

print(p_hasil)