Perbandingan tingkat kerusakan buah pada tiga sistem penyimpanan dengan uji Kruskal-Wallis

# Contoh kode R: Nonparametric Tests (Mann-Whitney and Kruskal-Wallis)
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(scales)
})
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: skor menceng pada tiga kelompok + satu variabel kontinu ----
n <- 200
lv <- c("A", "B", "C")
d <- data.frame(storage_system = factor(sample(rep(lv, length.out = n)), levels = lv),
                storage_duration = round(pmax(rnorm(n, 30, 6), 0), 2))
z2 <- as.numeric(scale(d$storage_duration))
efek <- c(0, runif(1, 0.3, 0.6), runif(1, 0.6, 1.1))[sample(3)]      # pergeseran lokasi antarkelompok (dalam SD)
bentuk <- runif(1, 1.5, 3)
galat <- (rgamma(n, shape = bentuk) - bentuk) / sqrt(bentuk)          # galat menceng kanan
z <- efek[as.integer(d$storage_system)] - mean(efek) + runif(1, 0.25, 0.5) * z2 + galat * runif(1, 0.7, 0.9)
d$fruit_damage_score <- pmax(100 + 15 * z, 0)
d$fruit_damage_score <- if (15 >= 3) round(d$fruit_damage_score) else round(d$fruit_damage_score, 2)    # skor bulat bila skalanya lebar
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Pemeriksaan normalitas (alasan memakai uji nonparametrik) ----
sw <- tapply(d$fruit_damage_score, d$storage_system, function(v) shapiro.test(v)$p.value)
print(sw)
print(aggregate(fruit_damage_score ~ storage_system, data = d, FUN = function(v) c(median = median(v), IQR = IQR(v), n = length(v))))

# ---- Kruskal-Wallis + ukuran efek epsilon kuadrat ----
kw <- kruskal.test(fruit_damage_score ~ storage_system, data = d)
print(kw)
eps2 <- unname(kw$statistic) / (nrow(d) - 1)

# ---- Perbandingan berpasangan Wilcoxon (koreksi Holm) ----
pw <- pairwise.wilcox.test(d$fruit_damage_score, d$storage_system, p.adjust.method = "holm", exact = FALSE)
print(pw)

# ---- Mann-Whitney: dua kelompok dengan median paling jauh + ukuran efek r ----
med <- tapply(d$fruit_damage_score, d$storage_system, median)
ekstrem <- names(c(which.min(med), which.max(med)))
if (ekstrem[1] == ekstrem[2]) ekstrem <- lv[c(1, 3)]
sub <- droplevels(d[d$storage_system %in% ekstrem, ])
mw <- wilcox.test(fruit_damage_score ~ storage_system, data = sub, exact = FALSE)
print(mw)
z_mw <- qnorm(mw$p.value / 2)
r_mw <- abs(z_mw) / sqrt(nrow(sub))

# ---- Korelasi peringkat Spearman ----
sp <- cor.test(d$fruit_damage_score, d$storage_duration, method = "spearman", exact = FALSE)
print(sp)


# ---- Grafik ----
bungkus <- scales::label_wrap(12)
pendek <- function(s, w = 12) ifelse(nchar(s) > w, paste0(trimws(substr(s, 1, w - 1)), "."), s)
dmed <- data.frame(storage_system = factor(names(med), levels = lv), median = as.numeric(med))
p1 <- ggplot(d, aes(storage_system, fruit_damage_score)) +
  geom_violin(aes(fill = storage_system), alpha = 0.3, colour = NA, trim = TRUE) +
  geom_boxplot(width = 0.18, outlier.size = 0.8, colour = warna$redup, fill = warna$panel, alpha = 0.8) +
  geom_point(data = dmed, aes(y = median), shape = 21, size = 3, fill = pal[2], colour = warna$teks) +
  scale_fill_manual(values = pal[c(1, 3, 4)], guide = "none") +
  scale_x_discrete(labels = bungkus) +
  labs(title = "Distribution by group (dot = median)", x = NULL, y = "fruit_damage_score")
mp <- matrix(NA_real_, 3, 3, dimnames = list(lv, lv))
mp[rownames(pw$p.value), colnames(pw$p.value)] <- pw$p.value
mp[upper.tri(mp)] <- t(mp)[upper.tri(mp)]
dp <- expand.grid(a = 1:3, b = 1:3)
dp$p <- mp[cbind(dp$a, dp$b)]
dp$teks <- ifelse(is.na(dp$p), "-", ifelse(dp$p < 0.001, "<0.001", sprintf("%.3f", dp$p)))
dp$skor <- ifelse(is.na(dp$p), NA, pmin(-log10(dp$p), 6))
dp$a <- factor(pendek(lv)[dp$a], levels = pendek(lv)); dp$b <- factor(pendek(lv)[dp$b], levels = rev(pendek(lv)))
p2 <- ggplot(dp, aes(a, b, fill = skor)) +
  geom_tile(colour = warna$panel, linewidth = 1.2) +
  geom_label(aes(label = teks), fill = warna$panel, colour = warna$teks, size = 3.2, label.size = 0, alpha = 0.85) +
  scale_fill_gradient(low = warna$grid, high = pal[2], na.value = warna$panel, limits = c(0, 6), guide = "none") +
  scale_x_discrete(labels = bungkus) + scale_y_discrete(labels = bungkus) +
  labs(title = "Pairwise Wilcoxon p-values (Holm)", x = NULL, y = NULL) +
  theme(panel.grid = element_blank())
dr <- data.frame(ry = rank(d$fruit_damage_score), rx = rank(d$storage_duration), kelompok = d$storage_system)
p3 <- ggplot(dr, aes(rx, ry)) +
  geom_point(aes(colour = kelompok), alpha = 0.55, size = 1.4) +
  geom_smooth(method = "lm", formula = y ~ x, colour = warna$teks, fill = warna$redup, alpha = 0.2, linewidth = 1) +
  scale_colour_manual(values = pal[c(1, 3, 4)], name = NULL, labels = pendek) +
  labs(title = sprintf("Spearman rank scatter (rho = %.2f)", sp$estimate), x = "Rank of storage_duration", y = "Rank of fruit_damage_score") +
  theme(legend.position = "bottom")
p_hasil <- if (tata == 1) p1 | (p2 / p3) else (p2 / p3) | p1

print(p_hasil)