Classifying cocoa land suitability with Random Forest from elevation and soil pH

# Contoh kode R: Random Forest
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(randomForest)
  library(pROC)
})
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: kelas biner dari X1, X2 dan fitur tambahan ----
n <- max(200, 300)
d <- data.frame(elevation = round(rnorm(n, 50, 10), 3),
                soil_ph = round(rnorm(n, 30, 6), 3),
                fitur_3 = round(rnorm(n, 50, 10), 2),            # fitur pendukung (berpengaruh sedang)
                fitur_4 = round(runif(n, 0, 100), 2),            # fitur pendukung (berpengaruh lemah)
                fitur_5 = round(rnorm(n, 20, 5), 2),             # fitur pengganggu (tanpa pengaruh)
                kategori = factor(sample(c("A", "B", "C"), n, replace = TRUE)))
z1 <- as.numeric(scale(d$elevation)); z2 <- as.numeric(scale(d$soil_ph))
z3 <- as.numeric(scale(d$fitur_3)); z4 <- as.numeric(scale(d$fitur_4))
eta <- runif(1, -0.2, 0.4) + runif(1, 1.8, 2.6) * z1 + runif(1, 1.1, 1.7) * sample(c(-1, 1), 1) * z2 +
  runif(1, 0.6, 1) * z3 + 0.3 * z4 - runif(1, 0.4, 0.7) * z1^2 + 0.6 * (d$kategori == "C") * z2
d$cocoa_suitable_land <- factor(ifelse(runif(n) < plogis(eta), "ya", "tidak"), levels = c("tidak", "ya"))
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Bagi data latih 70% dan data uji 30% (bertingkat per kelas) ----
id_latih <- unlist(lapply(split(seq_len(n), d$cocoa_suitable_land), function(i) sample(i, round(0.7 * length(i)))))
latih <- d[id_latih, ]
uji <- d[-id_latih, ]
print(table(latih$cocoa_suitable_land)); print(table(uji$cocoa_suitable_land))

# ---- Random Forest (500 pohon) ----
rf <- randomForest::randomForest(cocoa_suitable_land ~ ., data = latih, ntree = 500, importance = TRUE)
print(rf)
oob <- unname(rf$err.rate[rf$ntree, "OOB"])

# ---- Evaluasi pada data uji ----
prob <- predict(rf, uji, type = "prob")[, "ya"]
pred <- factor(ifelse(prob >= 0.5, "ya", "tidak"), levels = c("tidak", "ya"))
cm <- table(Aktual = uji$cocoa_suitable_land, Prediksi = pred)
print(cm)
tp <- cm["ya", "ya"]; fp <- cm["tidak", "ya"]; fn <- cm["ya", "tidak"]; tn <- cm["tidak", "tidak"]
akurasi <- (tp + tn) / sum(cm)
presisi <- if (tp + fp > 0) tp / (tp + fp) else 0
sensitivitas <- if (tp + fn > 0) tp / (tp + fn) else 0
spesifisitas <- if (tn + fp > 0) tn / (tn + fp) else 0
f1 <- if (presisi + sensitivitas > 0) 2 * presisi * sensitivitas / (presisi + sensitivitas) else 0
roc_uji <- pROC::roc(uji$cocoa_suitable_land, prob, levels = c("tidak", "ya"), direction = "<", quiet = TRUE)
auc <- as.numeric(pROC::auc(roc_uji))
print(roc_uji)

# ---- Kepentingan variabel ----
imp <- randomForest::importance(rf)
print(round(imp, 3))


# ---- Grafik ----
nama_tampil <- c("elevation" = "elevation", "soil_ph" = "soil_ph", fitur_3 = "Feature 3", fitur_4 = "Feature 4",
                 fitur_5 = "Feature 5", kategori = "Category")
nama_tampil[] <- make.unique(nama_tampil, sep = " ")              # cegah label kembar
di <- data.frame(fitur = nama_tampil[rownames(imp)], mda = imp[, "MeanDecreaseAccuracy"], gini = imp[, "MeanDecreaseGini"])
di$fitur <- factor(di$fitur, levels = di$fitur[order(di$mda)])
p1 <- ggplot(di, aes(mda, fitur)) +
  geom_segment(aes(x = 0, xend = mda, yend = fitur), colour = warna$redup, linewidth = 0.8) +
  geom_point(aes(size = gini), colour = pal[1]) +
  scale_size_continuous(range = c(2, 6), name = "Gini decrease") +
  labs(title = "Variable importance (mean decrease accuracy)", x = "Mean decrease in accuracy", y = NULL) +
  theme(legend.position = "bottom")
dr <- data.frame(fpr = 1 - roc_uji$specificities, tpr = roc_uji$sensitivities)
dr <- dr[order(dr$fpr, dr$tpr), ]
p2 <- ggplot(dr, aes(fpr, tpr)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = warna$redup) +
  geom_step(colour = pal[2], linewidth = 1.1) +
  annotate("label", x = 0.62, y = 0.18, label = sprintf("AUC = %.3f", auc), colour = warna$teks, fill = warna$panel,
           size = 3.5, label.size = 0) +
  coord_equal() +
  labs(title = "ROC curve (test set)", x = "False positive rate", y = "True positive rate")
pd <- randomForest::partialPlot(rf, pred.data = latih, x.var = "elevation", which.class = "ya", plot = FALSE, n.pt = 40)
p3 <- ggplot(data.frame(x = pd$x, y = pd$y), aes(x, y)) +
  geom_rug(data = latih, aes(x = elevation), inherit.aes = FALSE, colour = warna$redup, alpha = 0.35) +
  geom_line(colour = pal[3], linewidth = 1.2) +
  labs(title = "Partial dependence on elevation", x = "elevation", y = "Partial dependence (logit)")
de <- data.frame(pohon = seq_len(rf$ntree), oob = rf$err.rate[, "OOB"])
p4 <- ggplot(de, aes(pohon, oob)) +
  geom_line(colour = pal[4], linewidth = 0.9) +
  geom_hline(yintercept = oob, linetype = "dashed", colour = warna$redup) +
  labs(title = "OOB error vs number of trees", x = "Number of trees", y = "OOB error rate")
p_hasil <- if (tata == 1) (p1 | p2) / (p3 | p4) else (p2 | p1) / (p4 | p3)

print(p_hasil)