Bayesian estimation of soybean yield response to organic fertiliser and soil moisture

# Contoh kode R: Bayesian Regression (rstanarm)
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(rstanarm)
})
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(organic_fertiliser_dose = rnorm(n, 50, 10),
                soil_moisture = rnorm(n, 30, 6))
z1 <- as.numeric(scale(d$organic_fertiliser_dose))
z2 <- as.numeric(scale(d$soil_moisture))
b1 <- runif(1, 0.4, 0.8) * sample(c(-1, 1), 1, prob = c(0.25, 0.75))
b2 <- runif(1, 0.12, 0.45) * sample(c(-1, 1), 1)
d$soybean_yield <- 100 + 15 * (b1 * z1 + b2 * z2 + rnorm(n, 0, runif(1, 0.5, 0.85)))
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Analisis: regresi linear Bayesian (rstanarm, MCMC) ----
# Prior normal lemah (autoscale menyesuaikan skala data), 2 rantai x 1500 iterasi (750 warmup)
fit <- stan_glm(soybean_yield ~ organic_fertiliser_dose + soil_moisture, data = d, family = gaussian(),
                prior = normal(0, 2.5, autoscale = TRUE),
                prior_intercept = normal(0, 5, autoscale = TRUE),
                prior_aux = exponential(1, autoscale = TRUE),
                chains = 2, iter = 1500, refresh = 0, cores = 1, seed = sample.int(1e6, 1))
print(fit, digits = 3)
sm <- summary(fit, probs = c(0.025, 0.5, 0.975), digits = 3)
print(sm)
ci <- posterior_interval(fit, prob = 0.95)        # interval kredibel 95%
print(ci)
post <- as.matrix(fit)                           # draw posterior (baris = draw)
med <- apply(post, 2, median)
rhat_maks <- max(sm[, "Rhat"], na.rm = TRUE)       # konvergensi: idealnya < 1.01
ess_min <- min(sm[, "n_eff"], na.rm = TRUE)        # ukuran sampel efektif
r2_draw <- bayes_R2(fit)                         # Bayesian R2 per draw
p_pos1 <- mean(post[, "organic_fertiliser_dose"] > 0)             # peluang posterior efek positif
p_pos2 <- mean(post[, "soil_moisture"] > 0)
cat("Median Bayesian R2:", round(median(r2_draw), 3), "\n")
cat("P(koefisien organic_fertiliser_dose > 0):", round(p_pos1, 3), " P(koefisien soil_moisture > 0):", round(p_pos2, 3), "\n")


# ---- Grafik ----
# 1) densitas posterior koefisien dengan interval kredibel 95%
koef <- c("organic_fertiliser_dose", "soil_moisture")
lab_koef <- make.unique(c("organic_fertiliser_dose", "soil_moisture"), sep = " ")
dp <- do.call(rbind, lapply(1:2, function(j) {
  dz <- density(post[, koef[j]], n = 256)
  data.frame(koef = lab_koef[j], x = dz$x, y = dz$y,
             dalam = dz$x >= quantile(post[, koef[j]], 0.025) & dz$x <= quantile(post[, koef[j]], 0.975))
}))
dp$koef <- factor(dp$koef, levels = lab_koef)
dm <- data.frame(koef = factor(lab_koef, levels = lab_koef), m = apply(post[, koef], 2, median))
p1 <- ggplot(dp, aes(x, y)) +
  geom_area(data = subset(dp, dalam), fill = pal[1], alpha = 0.35) +
  geom_line(colour = pal[1], linewidth = 0.9) +
  geom_vline(data = dm, aes(xintercept = m), colour = pal[2], linewidth = 0.9) +
  geom_vline(xintercept = 0, linetype = "dashed", colour = warna$redup) +
  facet_wrap(~koef, scales = "free", ncol = 1) +
  labs(title = "Posterior of coefficients (95% CrI shaded)", x = "Coefficient value", y = "Posterior density")
# 2) posterior predictive check: densitas y amatan vs 40 replikasi
yrep <- posterior_predict(fit, draws = 40)
dr <- do.call(rbind, lapply(seq_len(nrow(yrep)), function(i) {
  dz <- density(yrep[i, ], n = 128)
  data.frame(rep = i, x = dz$x, y = dz$y)
}))
dy <- density(d$soybean_yield, n = 256)
p2 <- ggplot() +
  geom_line(data = dr, aes(x, y, group = rep, colour = "Replicated"), alpha = 0.3, linewidth = 0.4) +
  geom_line(data = data.frame(x = dy$x, y = dy$y), aes(x, y, colour = "Observed"), linewidth = 1.1) +
  scale_colour_manual(values = c(Observed = warna$teks, Replicated = pal[3]), name = NULL) +
  labs(title = "Posterior predictive check", x = "soybean_yield", y = "Density")
# 3) trace plot koefisien organic_fertiliser_dose per rantai
arr <- as.array(fit)
dt <- data.frame(iter = rep(seq_len(dim(arr)[1]), 2),
                 rantai = factor(rep(1:2, each = dim(arr)[1])),
                 nilai = c(arr[, 1, "organic_fertiliser_dose"], arr[, 2, "organic_fertiliser_dose"]))
p3 <- ggplot(dt, aes(iter, nilai, colour = rantai)) +
  geom_line(alpha = 0.7, linewidth = 0.35) +
  scale_colour_manual(values = pal[c(1, 4)], name = "Chain") +
  labs(title = "Trace of organic_fertiliser_dose coefficient", x = "Post-warmup iteration", y = "Coefficient")
p_hasil <- if (tata == 1) p1 | (p2 / p3) else (p1 | p2) / p3 + plot_layout(heights = c(1.4, 1))

print(p_hasil)