# =============================================================================
# Title:
# - Illustrating False Discovery Rate (FDR) Control with Simulated p-values. 
#
# Authors:
# - Oliver Gutiérrez-Hernández, Department of Geography, University of Málaga (UMA).
# - Luis V. García, Institute of Natural Resources and Agrobiology of Seville (IRNAS), 
#   Spanish National Research Council (CSIC).
#
# Objective:
# - To demonstrate how the False Discovery Rate can be controlled using:
#   1) the classical FDR-BH procedure (Benjamini & Hochberg, 1995), and
#   2) the adaptive FDR-BKY procedure (Benjamini, Krieger & Yekutieli, 2006).
# - To assess the reliability of the m₀/m ratio estimation through Monte Carlo simulation
#   (10,000 repetitions) under repeated random sampling with a known true m₀ = 60.
# - The example uses simulated ordered p-values, showing:
#   • how BH and BKY procedures differ in rejection thresholds and discoveries, and
#   • how stable the m₀/m ratio remains under repeated sampling.
#
# References:
# - Benjamini, Y., & Hochberg, Y. (1995). "Controlling the False Discovery Rate:
#   A Practical and Powerful Approach to Multiple Testing."
#   Journal of the Royal Statistical Society: Series B (Methodological), 57, 289–300.
#   https://doi.org/10.1111/j.2517-6161.1995.tb02031.x
#
# - Benjamini, Y., Krieger, A. M., & Yekutieli, D. (2006). 
#   "Adaptive Linear Step-Up Procedures that Control the False Discovery Rate."
#   Biometrika, 93(3), 491–507. https://doi.org/10.1093/biomet/93.3.491
#
# Practical Applications:
# - Provides a clear illustration of how adaptive FDR (BKY) improves upon classical BH:
#     • Increases discoveries when many true alternatives exist.
#     • Maintains rigorous control of false positives.
# - Includes a quantitative evaluation of the stability of the m₀/m ratio.
# - Relevant for any field requiring multiple hypothesis testing:
#     • Environmental sciences (trend detection in gridded data).
#     • Remote sensing (pixel-wise analysis of NDVI, EVI, or land cover).
#     • Genomics and proteomics (large-scale testing of markers or variants).
#     • Social sciences or economics (high-dimensional hypothesis testing).
#
# Input:
# - Simulated p-values sorted in ascending order.
#
# Outputs:
# 1. A plot of simulated p-values with BH and BKY thresholds.
# 2. Classification of discoveries (shared, additional under BKY, non-significant).
# 3. Monte Carlo summary statistics for the stability of m₀/m estimation:
#    mean, bias, SD, RMSE, and 95% empirical interval.
# =============================================================================


# --- Parameters
m <- 100
q <- 0.05
m0_hat <- 60  # adaptive estimation

# --- Simulate ordered p-values
set.seed(123)
pvals <- sort(c(runif(40, 0, 0.01), runif(60, 0, 0.2)))
i <- 1:m

# --- Thresholds
t_bh  <- (i / m) * q
q_star <- (q * m) / m0_hat
t_bky <- (i / m) * q_star

# --- Classification
discoveries_bh <- pvals <= t_bh
discoveries_bky <- pvals <= t_bky
shared <- discoveries_bh & discoveries_bky
only_bky <- discoveries_bky & !discoveries_bh
nonsignificant <- !discoveries_bky

# --- Plot
plot(i, pvals, type = "n",
     xlab = expression("Rank of "*italic(p)*"-value"),
     ylab = expression(italic(p)*"-value"),
     ylim = c(0, max(pvals, t_bky)))

# --- Title inside the plot with q in italics and standard size
text(x = par("usr")[1] + 2, 
     y = par("usr")[4] - 0.005,
     labels = expression("False Discovery Rate (FDR) control ("*italic(q)*" = 0.05)"),
     font = 2, adj = c(0, 1))

# --- Threshold lines
lines(i, t_bky, col = "skyblue", lwd = 2)
lines(i, t_bh, col = "skyblue", lwd = 2, lty = 2)

# --- Points (blue ones made finer)
points(i[shared], pvals[shared], pch = 19, col = "blue", cex = 0.6)
points(i[only_bky], pvals[only_bky], pch = 17, col = "blue", cex = 0.7)
points(i[nonsignificant], pvals[nonsignificant], pch = 1, col = "black")

# --- Final spacious and complete legend
legend("topleft",
       inset = c(0.01, 0.05),
       legend = c(expression("FDR-BKY threshold ("*italic(m)[0]*" = 60)"),
                  "FDR-BH threshold",
                  "Discoveries under BH and BKY",
                  "Additional discoveries under BKY",
                  "Not discoveries"),
       col = c("skyblue", "skyblue", "blue", "blue", "black"),
       lty = c(1, 2, NA, NA, NA),
       pch = c(NA, NA, 19, 17, 1),
       lwd = c(2, 2, NA, NA, NA),
       pt.cex = 1.2,
       bty = "n",
       x.intersp = 1.1,
       y.intersp = 1.4)


# =============================================================================
# --- Reliability assessment of m0/m ratio
# =============================================================================
# Objective:
# - To evaluate the stability of the ratio m0/m (= 0.6) under repeated random sampling
#   with the same true m0 = 60 and m = 100.

set.seed(123)
n_iter <- 10000
m0_true <- 60
m_true <- 100
pi0_true <- m0_true / m_true

pi0_hat <- numeric(n_iter)

for (k in 1:n_iter) {
  # 60 nulls (uniform 0,1) + 40 alternatives (uniform 0,0.01)
  pvals <- c(runif(40, 0, 0.01), runif(60, 0, 1))
  i <- 1:m_true
  t_bh <- (i / m_true) * q
  R <- max(which(sort(pvals) <= t_bh), 0)
  m0_hat_iter <- m_true - R * (1 - q)
  pi0_hat[k] <- m0_hat_iter / m_true
}

# --- Summary statistics
mean_hat <- mean(pi0_hat)
bias <- mean_hat - pi0_true
sd_hat <- sd(pi0_hat)
rmse <- sqrt(mean((pi0_hat - pi0_true)^2))
ci <- quantile(pi0_hat, c(0.025, 0.975))

cat("\n--- Reliability of m0/m ratio estimation (Monte Carlo, n = 10000) ---\n")
cat(sprintf("True m0/m: %.2f\n", pi0_true))
cat(sprintf("Mean estimated ratio: %.3f\n", mean_hat))
cat(sprintf("Bias: %.3f, SD: %.3f, RMSE: %.3f\n", bias, sd_hat, rmse))
cat(sprintf("95%% empirical interval: [%.3f, %.3f]\n\n", ci[1], ci[2]))
# =============================================================================
