CASH Simulations

\(g_1 = N\left(0, 2^2\right)\)

q.vec <- seq(0.001, 0.20, by = 0.001)
method.name <- c("BHq", "qvalue", "locfdr", "ASH", "CASH")

FDP.array <- pFDP.array <- power.array <- array(0, dim = c(nsim, length(q.vec), length(method.name), length(pi0.vec)))
FDP.summary <- array(0, dim = c(7, length(q.vec), length(method.name), length(pi0.vec)))
pFDP.summary <- power.summary <- array(0, dim = c(5, length(q.vec), length(method.name), length(pi0.vec)))
for (j in seq(length(pi0.vec))) {
  for (k in seq(length(method.name))) {
    for (i in seq(nsim)) {
      FDP.array[i, , k, j] <- sapply(q.vec, FDP, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
      pFDP.array[i, , k, j] <- sapply(q.vec, pFDP, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
      power.array[i, , k, j] <- sapply(q.vec, power, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
    }
    FDP.summary[, , k, j] <- rbind(
      avg <- colMeans(FDP.array[, , k, j], na.rm = TRUE),
      sd <- apply(FDP.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(FDP.array[, , k, j])),
      q975 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE),
      q750 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.75, na.rm = TRUE),
      q250 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.25, na.rm = TRUE)
    )
    pFDP.summary[, , k, j] <- rbind(
      avg <- colMeans(pFDP.array[, , k, j], na.rm = TRUE),
      sd <- apply(pFDP.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(pFDP.array[, , k, j])),
      q975 <- apply(pFDP.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(pFDP.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE)
    )
    power.summary[, , k, j] <- rbind(
      avg <- colMeans(power.array[, , k, j], na.rm = TRUE),
      sd <- apply(power.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(power.array[, , k, j])),
      q975 <- apply(power.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(power.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE)
    )
  }
}

q <- 0.1
z.over <- 1.05
z.under <- 0.95
method.col <- scales::hue_pal()(5)
# method.col <- c("#377eb8", "#984ea3", "#4daf4a", "#ff7f00", "#e41a1c")

for (j in seq(length(pi0.vec))) {
  sd.z <- sapply(z.pi0.list[[j]], sd)
  Noise <- cut(sd.z, breaks = c(0, quantile(sd.z, probs = 1 : 2 / 3), Inf), labels = c("Deflated Noise", "In-between", "Inflated Noise"))
# Noise <- cut(sd.z, breaks = c(0, z.under, z.over, Inf), labels = c("Deflated Noise", "In-between", "Inflated Noise"))

  pi0.pi0 <- matrix(unlist(pi0.pi0.list[[j]]), byrow = TRUE, length(pi0.pi0.list[[j]]))
  pi0.pi0.noise <- rbind.data.frame(cbind.data.frame(Noise, pi0.pi0), cbind.data.frame(Noise = rep("All", length(Noise)), pi0.pi0))
  
  pi0.plot <- ggplot(data = melt(pi0.pi0.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col[-1]) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_hline(yintercept = pi0.vec[j], col = "black", linetype = "dashed", size = 1) +
  scale_x_discrete(labels = method.name[-1]) +
  labs(x = "", y = expression(hat(pi)[0])) +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))

  FDP.summary.pi0 <- aperm(FDP.summary[, , , j], c(2, 1, 3))
  FDP.summary.pi0.method <- FDP.summary.pi0[, , 1]
  for (kk in 2 : length(method.name)) {
    FDP.summary.pi0.method <- rbind.data.frame(FDP.summary.pi0.method, FDP.summary.pi0[, , kk])
  }
  FDP.summary.pi0.method <- cbind.data.frame(
    rep(factor(seq(method.name)), each = dim(FDP.summary.pi0)[1]),
    rep(q.vec, length(method.name)),
    FDP.summary.pi0.method
  )
  colnames(FDP.summary.pi0.method) <- c(
    "Method", "FDR", "FDP", "sd", "n", "q975", "q025", "q750", "q250"
  )
  
  FDP.array.pi0 <- aperm(FDP.array[, , , j], c(2, 1, 3))
  FDP.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, mean, na.rm = TRUE), c(2, 1, 3)))
  sd.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, sd, na.rm = TRUE), c(2, 1, 3)))
  n.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, function(x){sum(!is.na(x))}), c(2, 1, 3)))
  q975.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.975, na.rm = TRUE), c(2, 1, 3)))
  q025.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.025, na.rm = TRUE), c(2, 1, 3)))
  q750.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.75, na.rm = TRUE), c(2, 1, 3)))
  q250.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.25, na.rm = TRUE), c(2, 1, 3)))
  FDP.summary.pi0.method.noise <- cbind.data.frame(
    rep(rep(levels(Noise), each = length(q.vec)), length(method.name)),
    rep(factor(seq(method.name)), each = length(levels(Noise)) * length(q.vec)),
    rep(q.vec, length(levels(Noise)) * length(method.name)),
    FDP.pi0.noise,
    sd.pi0.noise,
    n.pi0.noise,
    q975.pi0.noise,
    q025.pi0.noise,
    q750.pi0.noise,
    q250.pi0.noise
  )
  colnames(FDP.summary.pi0.method.noise) <- c(
    "Noise", "Method", "FDR", "FDP", "sd", "n", "q975", "q025", "q750", "q250"
  )
  FDP.summary.pi0.method.noise <- rbind.data.frame(
    FDP.summary.pi0.method.noise,
    cbind.data.frame(Noise = rep("All", dim(FDP.summary.pi0.method)[1]), FDP.summary.pi0.method)
  )
  
  FDR.calib.plot <- ggplot(data = FDP.summary.pi0.method.noise, aes(x = FDR, y = FDP, group = Method, col = Method)) +
  geom_line() +
  geom_ribbon(aes(ymin = q025, ymax = q975, fill = Method), alpha = 0.35, linetype = "blank") +
  scale_color_manual(labels = method.name, values = method.col) +
  scale_fill_manual(labels = method.name, values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", size = 1, col = "black") +
  labs(x = "Nominal FDR", y = "FDP") +
  theme(legend.position = "top", legend.text = element_text(size = 15), plot.title = element_text(hjust = 0.5, size = 15), axis.title.x = element_text(size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(angle = 45, size = 15), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  FDP.q <- FDP.array[, which(round(q.vec, 4) == q), , j]
  FDP.q.noise <- rbind.data.frame(cbind.data.frame(Noise, FDP.q), cbind.data.frame(Noise = rep("All", length(Noise)), FDP.q))

  FDR.plot <- ggplot(data = melt(FDP.q.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_hline(yintercept = q, col = "black", linetype = "dashed", size = 1) +
  scale_x_discrete(labels = method.name) +
  labs(x = "", y = "FDP") +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  TDP.q <- power.array[, which(round(q.vec, 4) == q), , j]
  TDP.q.noise <- rbind.data.frame(cbind.data.frame(Noise, TDP.q), cbind.data.frame(Noise = rep("All", length(Noise)), TDP.q))

  power.plot <- ggplot(data = melt(TDP.q.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  scale_x_discrete(labels = method.name) +
  labs(x = "", y = "TPP") +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  joint <- ggarrange(FDR.calib.plot,
            pi0.plot + rremove("x.text"),
            FDR.plot + rremove("x.text"),
            power.plot,
            align = "v", ncol = 1, nrow = 4,
            heights = c(1.5, 1, 1, 1.2)
  )
  joint <- annotate_figure(joint,
    top = text_grob(bquote(pi[0] == .(pi0.vec[j])), size = 15)
  )
  print(joint)
  ggsave(paste0("../output/fig/g1_pi0_", pi0.vec[j], ".pdf"), joint, height = 10, width = 8)
}

Bimodal: \(g_2 = 0.5 N\left(-2, 1\right) + 0.5 N\left(2, 1\right)\)

FDP.array <- pFDP.array <- power.array <- array(0, dim = c(nsim, length(q.vec), length(method.name), length(pi0.vec)))
FDP.summary <- array(0, dim = c(7, length(q.vec), length(method.name), length(pi0.vec)))
pFDP.summary <- power.summary <- array(0, dim = c(5, length(q.vec), length(method.name), length(pi0.vec)))
for (j in seq(length(pi0.vec))) {
  for (k in seq(length(method.name))) {
    for (i in seq(nsim)) {
      FDP.array[i, , k, j] <- sapply(q.vec, FDP, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
      pFDP.array[i, , k, j] <- sapply(q.vec, pFDP, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
      power.array[i, , k, j] <- sapply(q.vec, power, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
    }
    FDP.summary[, , k, j] <- rbind(
      avg <- colMeans(FDP.array[, , k, j], na.rm = TRUE),
      sd <- apply(FDP.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(FDP.array[, , k, j])),
      q975 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE),
      q750 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.75, na.rm = TRUE),
      q250 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.25, na.rm = TRUE)
    )
    pFDP.summary[, , k, j] <- rbind(
      avg <- colMeans(pFDP.array[, , k, j], na.rm = TRUE),
      sd <- apply(pFDP.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(pFDP.array[, , k, j])),
      q975 <- apply(pFDP.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(pFDP.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE)
    )
    power.summary[, , k, j] <- rbind(
      avg <- colMeans(power.array[, , k, j], na.rm = TRUE),
      sd <- apply(power.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(power.array[, , k, j])),
      q975 <- apply(power.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(power.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE)
    )
  }
}

for (j in seq(length(pi0.vec))) {
  sd.z <- sapply(z.pi0.list[[j]], sd)
  Noise <- cut(sd.z, breaks = c(0, quantile(sd.z, probs = 1 : 2 / 3), Inf), labels = c("Deflated Noise", "In-between", "Inflated Noise"))
# Noise <- cut(sd.z, breaks = c(0, z.under, z.over, Inf), labels = c("Deflated Noise", "In-between", "Inflated Noise"))

  pi0.pi0 <- matrix(unlist(pi0.pi0.list[[j]]), byrow = TRUE, length(pi0.pi0.list[[j]]))
  pi0.pi0.noise <- rbind.data.frame(cbind.data.frame(Noise, pi0.pi0), cbind.data.frame(Noise = rep("All", length(Noise)), pi0.pi0))
  
  pi0.plot <- ggplot(data = melt(pi0.pi0.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col[-1]) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_hline(yintercept = pi0.vec[j], col = "black", linetype = "dashed", size = 1) +
  scale_x_discrete(labels = method.name[-1]) +
  labs(x = "", y = expression(hat(pi)[0])) +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))

  FDP.summary.pi0 <- aperm(FDP.summary[, , , j], c(2, 1, 3))
  FDP.summary.pi0.method <- FDP.summary.pi0[, , 1]
  for (kk in 2 : length(method.name)) {
    FDP.summary.pi0.method <- rbind.data.frame(FDP.summary.pi0.method, FDP.summary.pi0[, , kk])
  }
  FDP.summary.pi0.method <- cbind.data.frame(
    rep(factor(seq(method.name)), each = dim(FDP.summary.pi0)[1]),
    rep(q.vec, length(method.name)),
    FDP.summary.pi0.method
  )
  colnames(FDP.summary.pi0.method) <- c(
    "Method", "FDR", "FDP", "sd", "n", "q975", "q025", "q750", "q250"
  )
  
  FDP.array.pi0 <- aperm(FDP.array[, , , j], c(2, 1, 3))
  FDP.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, mean, na.rm = TRUE), c(2, 1, 3)))
  sd.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, sd, na.rm = TRUE), c(2, 1, 3)))
  n.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, function(x){sum(!is.na(x))}), c(2, 1, 3)))
  q975.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.975, na.rm = TRUE), c(2, 1, 3)))
  q025.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.025, na.rm = TRUE), c(2, 1, 3)))
  q750.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.75, na.rm = TRUE), c(2, 1, 3)))
  q250.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.25, na.rm = TRUE), c(2, 1, 3)))
  FDP.summary.pi0.method.noise <- cbind.data.frame(
    rep(rep(levels(Noise), each = length(q.vec)), length(method.name)),
    rep(factor(seq(method.name)), each = length(levels(Noise)) * length(q.vec)),
    rep(q.vec, length(levels(Noise)) * length(method.name)),
    FDP.pi0.noise,
    sd.pi0.noise,
    n.pi0.noise,
    q975.pi0.noise,
    q025.pi0.noise,
    q750.pi0.noise,
    q250.pi0.noise
  )
  colnames(FDP.summary.pi0.method.noise) <- c(
    "Noise", "Method", "FDR", "FDP", "sd", "n", "q975", "q025", "q750", "q250"
  )
  FDP.summary.pi0.method.noise <- rbind.data.frame(
    FDP.summary.pi0.method.noise,
    cbind.data.frame(Noise = rep("All", dim(FDP.summary.pi0.method)[1]), FDP.summary.pi0.method)
  )
  
  FDR.calib.plot <- ggplot(data = FDP.summary.pi0.method.noise, aes(x = FDR, y = FDP, group = Method, col = Method)) +
  geom_line() +
  geom_ribbon(aes(ymin = q025, ymax = q975, fill = Method), alpha = 0.35, linetype = "blank") +
  scale_color_manual(labels = method.name, values = method.col) +
  scale_fill_manual(labels = method.name, values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", size = 1, col = "black") +
  labs(x = "Nominal FDR", y = "FDP") +
  theme(legend.position = "top", legend.text = element_text(size = 15), plot.title = element_text(hjust = 0.5, size = 15), axis.title.x = element_text(size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(angle = 45, size = 15), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  FDP.q <- FDP.array[, which(round(q.vec, 4) == q), , j]
  FDP.q.noise <- rbind.data.frame(cbind.data.frame(Noise, FDP.q), cbind.data.frame(Noise = rep("All", length(Noise)), FDP.q))

  FDR.plot <- ggplot(data = melt(FDP.q.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_hline(yintercept = q, col = "black", linetype = "dashed", size = 1) +
  scale_x_discrete(labels = method.name) +
  labs(x = "", y = "FDP") +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  TDP.q <- power.array[, which(round(q.vec, 4) == q), , j]
  TDP.q.noise <- rbind.data.frame(cbind.data.frame(Noise, TDP.q), cbind.data.frame(Noise = rep("All", length(Noise)), TDP.q))

  power.plot <- ggplot(data = melt(TDP.q.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  scale_x_discrete(labels = method.name) +
  labs(x = "", y = "TPP") +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  joint <- ggarrange(FDR.calib.plot,
            pi0.plot + rremove("x.text"),
            FDR.plot + rremove("x.text"),
            power.plot,
            align = "v", ncol = 1, nrow = 4,
            heights = c(1.5, 1, 1, 1.2)
  )
  joint <- annotate_figure(joint,
    top = text_grob(bquote(pi[0] == .(pi0.vec[j])), size = 15)
  )
  print(joint)
  ggsave(paste0("../output/fig/g2_pi0_", pi0.vec[j], ".pdf"), joint, height = 10, width = 8)
}

Spiky: \(g_3 = 0.4 N\left(0, 0.5^2\right) + 0.2 N\left(0, 1^2\right) + 0.2 N\left(0, 2^2\right) + 0.2 N\left(0, 3^2\right)\)

FDP.array <- pFDP.array <- power.array <- array(0, dim = c(nsim, length(q.vec), length(method.name), length(pi0.vec)))
FDP.summary <- array(0, dim = c(7, length(q.vec), length(method.name), length(pi0.vec)))
pFDP.summary <- power.summary <- array(0, dim = c(5, length(q.vec), length(method.name), length(pi0.vec)))
for (j in seq(length(pi0.vec))) {
  for (k in seq(length(method.name))) {
    for (i in seq(nsim)) {
      FDP.array[i, , k, j] <- sapply(q.vec, FDP, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
      pFDP.array[i, , k, j] <- sapply(q.vec, pFDP, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
      power.array[i, , k, j] <- sapply(q.vec, power, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
    }
    FDP.summary[, , k, j] <- rbind(
      avg <- colMeans(FDP.array[, , k, j], na.rm = TRUE),
      sd <- apply(FDP.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(FDP.array[, , k, j])),
      q975 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE),
      q750 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.75, na.rm = TRUE),
      q250 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.25, na.rm = TRUE)
    )
    pFDP.summary[, , k, j] <- rbind(
      avg <- colMeans(pFDP.array[, , k, j], na.rm = TRUE),
      sd <- apply(pFDP.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(pFDP.array[, , k, j])),
      q975 <- apply(pFDP.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(pFDP.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE)
    )
    power.summary[, , k, j] <- rbind(
      avg <- colMeans(power.array[, , k, j], na.rm = TRUE),
      sd <- apply(power.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(power.array[, , k, j])),
      q975 <- apply(power.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(power.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE)
    )
  }
}

for (j in seq(length(pi0.vec))) {
  sd.z <- sapply(z.pi0.list[[j]], sd)
  Noise <- cut(sd.z, breaks = c(0, quantile(sd.z, probs = 1 : 2 / 3), Inf), labels = c("Deflated Noise", "In-between", "Inflated Noise"))
# Noise <- cut(sd.z, breaks = c(0, z.under, z.over, Inf), labels = c("Deflated Noise", "In-between", "Inflated Noise"))

  pi0.pi0 <- matrix(unlist(pi0.pi0.list[[j]]), byrow = TRUE, length(pi0.pi0.list[[j]]))
  pi0.pi0.noise <- rbind.data.frame(cbind.data.frame(Noise, pi0.pi0), cbind.data.frame(Noise = rep("All", length(Noise)), pi0.pi0))
  
  pi0.plot <- ggplot(data = melt(pi0.pi0.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col[-1]) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_hline(yintercept = pi0.vec[j], col = "black", linetype = "dashed", size = 1) +
  scale_x_discrete(labels = method.name[-1]) +
  labs(x = "", y = expression(hat(pi)[0])) +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))

  FDP.summary.pi0 <- aperm(FDP.summary[, , , j], c(2, 1, 3))
  FDP.summary.pi0.method <- FDP.summary.pi0[, , 1]
  for (kk in 2 : length(method.name)) {
    FDP.summary.pi0.method <- rbind.data.frame(FDP.summary.pi0.method, FDP.summary.pi0[, , kk])
  }
  FDP.summary.pi0.method <- cbind.data.frame(
    rep(factor(seq(method.name)), each = dim(FDP.summary.pi0)[1]),
    rep(q.vec, length(method.name)),
    FDP.summary.pi0.method
  )
  colnames(FDP.summary.pi0.method) <- c(
    "Method", "FDR", "FDP", "sd", "n", "q975", "q025", "q750", "q250"
  )
  
  FDP.array.pi0 <- aperm(FDP.array[, , , j], c(2, 1, 3))
  FDP.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, mean, na.rm = TRUE), c(2, 1, 3)))
  sd.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, sd, na.rm = TRUE), c(2, 1, 3)))
  n.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, function(x){sum(!is.na(x))}), c(2, 1, 3)))
  q975.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.975, na.rm = TRUE), c(2, 1, 3)))
  q025.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.025, na.rm = TRUE), c(2, 1, 3)))
  q750.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.75, na.rm = TRUE), c(2, 1, 3)))
  q250.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.25, na.rm = TRUE), c(2, 1, 3)))
  FDP.summary.pi0.method.noise <- cbind.data.frame(
    rep(rep(levels(Noise), each = length(q.vec)), length(method.name)),
    rep(factor(seq(method.name)), each = length(levels(Noise)) * length(q.vec)),
    rep(q.vec, length(levels(Noise)) * length(method.name)),
    FDP.pi0.noise,
    sd.pi0.noise,
    n.pi0.noise,
    q975.pi0.noise,
    q025.pi0.noise,
    q750.pi0.noise,
    q250.pi0.noise
  )
  colnames(FDP.summary.pi0.method.noise) <- c(
    "Noise", "Method", "FDR", "FDP", "sd", "n", "q975", "q025", "q750", "q250"
  )
  FDP.summary.pi0.method.noise <- rbind.data.frame(
    FDP.summary.pi0.method.noise,
    cbind.data.frame(Noise = rep("All", dim(FDP.summary.pi0.method)[1]), FDP.summary.pi0.method)
  )
  
  FDR.calib.plot <- ggplot(data = FDP.summary.pi0.method.noise, aes(x = FDR, y = FDP, group = Method, col = Method)) +
  geom_line() +
  geom_ribbon(aes(ymin = q025, ymax = q975, fill = Method), alpha = 0.35, linetype = "blank") +
  scale_color_manual(labels = method.name, values = method.col) +
  scale_fill_manual(labels = method.name, values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", size = 1, col = "black") +
  labs(x = "Nominal FDR", y = "FDP") +
  theme(legend.position = "top", legend.text = element_text(size = 15), plot.title = element_text(hjust = 0.5, size = 15), axis.title.x = element_text(size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(angle = 45, size = 15), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  FDP.q <- FDP.array[, which(round(q.vec, 4) == q), , j]
  FDP.q.noise <- rbind.data.frame(cbind.data.frame(Noise, FDP.q), cbind.data.frame(Noise = rep("All", length(Noise)), FDP.q))

  FDR.plot <- ggplot(data = melt(FDP.q.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_hline(yintercept = q, col = "black", linetype = "dashed", size = 1) +
  scale_x_discrete(labels = method.name) +
  labs(x = "", y = "FDP") +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  TDP.q <- power.array[, which(round(q.vec, 4) == q), , j]
  TDP.q.noise <- rbind.data.frame(cbind.data.frame(Noise, TDP.q), cbind.data.frame(Noise = rep("All", length(Noise)), TDP.q))

  power.plot <- ggplot(data = melt(TDP.q.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  scale_x_discrete(labels = method.name) +
  labs(x = "", y = "TPP") +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  joint <- ggarrange(FDR.calib.plot,
            pi0.plot + rremove("x.text"),
            FDR.plot + rremove("x.text"),
            power.plot,
            align = "v", ncol = 1, nrow = 4,
            heights = c(1.5, 1, 1, 1.2)
  )
  joint <- annotate_figure(joint,
    top = text_grob(bquote(pi[0] == .(pi0.vec[j])), size = 15)
  )
  print(joint)
  ggsave(paste0("../output/fig/g3_pi0_", pi0.vec[j], ".pdf"), joint, height = 10, width = 8)
}

Flattop: \(g_4 = 1 / 13 \left( N\left(-3, 0.5^2\right) + N\left(-2.5, 0.5^2\right) + N\left(-2, 0.5^2\right) + N\left(-1.5, 0.5^2\right) + N\left(-1, 0.5^2\right) + N\left(-0.5, 0.5^2\right) + N\left(0, 0.5^2\right) + N\left(0.5, 0.5^2\right) + N\left(1, 0.5^2\right) + N\left(1.5, 0.5^2\right) + N\left(2, 0.5^2\right) + N\left(2.5, 0.5^2\right) + N\left(3, 0.5^2\right) \right)\)

FDP.array <- pFDP.array <- power.array <- array(0, dim = c(nsim, length(q.vec), length(method.name), length(pi0.vec)))
FDP.summary <- array(0, dim = c(7, length(q.vec), length(method.name), length(pi0.vec)))
pFDP.summary <- power.summary <- array(0, dim = c(5, length(q.vec), length(method.name), length(pi0.vec)))
for (j in seq(length(pi0.vec))) {
  for (k in seq(length(method.name))) {
    for (i in seq(nsim)) {
      FDP.array[i, , k, j] <- sapply(q.vec, FDP, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
      pFDP.array[i, , k, j] <- sapply(q.vec, pFDP, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
      power.array[i, , k, j] <- sapply(q.vec, power, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
    }
    FDP.summary[, , k, j] <- rbind(
      avg <- colMeans(FDP.array[, , k, j], na.rm = TRUE),
      sd <- apply(FDP.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(FDP.array[, , k, j])),
      q975 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE),
      q750 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.75, na.rm = TRUE),
      q250 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.25, na.rm = TRUE)
    )
    pFDP.summary[, , k, j] <- rbind(
      avg <- colMeans(pFDP.array[, , k, j], na.rm = TRUE),
      sd <- apply(pFDP.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(pFDP.array[, , k, j])),
      q975 <- apply(pFDP.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(pFDP.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE)
    )
    power.summary[, , k, j] <- rbind(
      avg <- colMeans(power.array[, , k, j], na.rm = TRUE),
      sd <- apply(power.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(power.array[, , k, j])),
      q975 <- apply(power.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(power.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE)
    )
  }
}

for (j in seq(length(pi0.vec))) {
  sd.z <- sapply(z.pi0.list[[j]], sd)
  Noise <- cut(sd.z, breaks = c(0, quantile(sd.z, probs = 1 : 2 / 3), Inf), labels = c("Deflated Noise", "In-between", "Inflated Noise"))
# Noise <- cut(sd.z, breaks = c(0, z.under, z.over, Inf), labels = c("Deflated Noise", "In-between", "Inflated Noise"))

  pi0.pi0 <- matrix(unlist(pi0.pi0.list[[j]]), byrow = TRUE, length(pi0.pi0.list[[j]]))
  pi0.pi0.noise <- rbind.data.frame(cbind.data.frame(Noise, pi0.pi0), cbind.data.frame(Noise = rep("All", length(Noise)), pi0.pi0))
  
  pi0.plot <- ggplot(data = melt(pi0.pi0.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col[-1]) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_hline(yintercept = pi0.vec[j], col = "black", linetype = "dashed", size = 1) +
  scale_x_discrete(labels = method.name[-1]) +
  labs(x = "", y = expression(hat(pi)[0])) +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))

  FDP.summary.pi0 <- aperm(FDP.summary[, , , j], c(2, 1, 3))
  FDP.summary.pi0.method <- FDP.summary.pi0[, , 1]
  for (kk in 2 : length(method.name)) {
    FDP.summary.pi0.method <- rbind.data.frame(FDP.summary.pi0.method, FDP.summary.pi0[, , kk])
  }
  FDP.summary.pi0.method <- cbind.data.frame(
    rep(factor(seq(method.name)), each = dim(FDP.summary.pi0)[1]),
    rep(q.vec, length(method.name)),
    FDP.summary.pi0.method
  )
  colnames(FDP.summary.pi0.method) <- c(
    "Method", "FDR", "FDP", "sd", "n", "q975", "q025", "q750", "q250"
  )
  
  FDP.array.pi0 <- aperm(FDP.array[, , , j], c(2, 1, 3))
  FDP.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, mean, na.rm = TRUE), c(2, 1, 3)))
  sd.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, sd, na.rm = TRUE), c(2, 1, 3)))
  n.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, function(x){sum(!is.na(x))}), c(2, 1, 3)))
  q975.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.975, na.rm = TRUE), c(2, 1, 3)))
  q025.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.025, na.rm = TRUE), c(2, 1, 3)))
  q750.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.75, na.rm = TRUE), c(2, 1, 3)))
  q250.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.25, na.rm = TRUE), c(2, 1, 3)))
  FDP.summary.pi0.method.noise <- cbind.data.frame(
    rep(rep(levels(Noise), each = length(q.vec)), length(method.name)),
    rep(factor(seq(method.name)), each = length(levels(Noise)) * length(q.vec)),
    rep(q.vec, length(levels(Noise)) * length(method.name)),
    FDP.pi0.noise,
    sd.pi0.noise,
    n.pi0.noise,
    q975.pi0.noise,
    q025.pi0.noise,
    q750.pi0.noise,
    q250.pi0.noise
  )
  colnames(FDP.summary.pi0.method.noise) <- c(
    "Noise", "Method", "FDR", "FDP", "sd", "n", "q975", "q025", "q750", "q250"
  )
  FDP.summary.pi0.method.noise <- rbind.data.frame(
    FDP.summary.pi0.method.noise,
    cbind.data.frame(Noise = rep("All", dim(FDP.summary.pi0.method)[1]), FDP.summary.pi0.method)
  )
  
  FDR.calib.plot <- ggplot(data = FDP.summary.pi0.method.noise, aes(x = FDR, y = FDP, group = Method, col = Method)) +
  geom_line() +
  geom_ribbon(aes(ymin = q025, ymax = q975, fill = Method), alpha = 0.35, linetype = "blank") +
  scale_color_manual(labels = method.name, values = method.col) +
  scale_fill_manual(labels = method.name, values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", size = 1, col = "black") +
  labs(x = "Nominal FDR", y = "FDP") +
  theme(legend.position = "top", legend.text = element_text(size = 15), plot.title = element_text(hjust = 0.5, size = 15), axis.title.x = element_text(size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(angle = 45, size = 15), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  FDP.q <- FDP.array[, which(round(q.vec, 4) == q), , j]
  FDP.q.noise <- rbind.data.frame(cbind.data.frame(Noise, FDP.q), cbind.data.frame(Noise = rep("All", length(Noise)), FDP.q))

  FDR.plot <- ggplot(data = melt(FDP.q.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_hline(yintercept = q, col = "black", linetype = "dashed", size = 1) +
  scale_x_discrete(labels = method.name) +
  labs(x = "", y = "FDP") +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  TDP.q <- power.array[, which(round(q.vec, 4) == q), , j]
  TDP.q.noise <- rbind.data.frame(cbind.data.frame(Noise, TDP.q), cbind.data.frame(Noise = rep("All", length(Noise)), TDP.q))

  power.plot <- ggplot(data = melt(TDP.q.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  scale_x_discrete(labels = method.name) +
  labs(x = "", y = "TPP") +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  joint <- ggarrange(FDR.calib.plot,
            pi0.plot + rremove("x.text"),
            FDR.plot + rremove("x.text"),
            power.plot,
            align = "v", ncol = 1, nrow = 4,
            heights = c(1.5, 1, 1, 1.2)
  )
  joint <- annotate_figure(joint,
    top = text_grob(bquote(pi[0] == .(pi0.vec[j])), size = 15)
  )
  print(joint)
  ggsave(paste0("../output/fig/g4_pi0_", pi0.vec[j], ".pdf"), joint, height = 10, width = 8)
}

Near normal: \(g_5 = 0.6 N\left(0, 1^2\right) + 0.4 N\left(0, 3^2\right)\)

FDP.array <- pFDP.array <- power.array <- array(0, dim = c(nsim, length(q.vec), length(method.name), length(pi0.vec)))
FDP.summary <- array(0, dim = c(7, length(q.vec), length(method.name), length(pi0.vec)))
pFDP.summary <- power.summary <- array(0, dim = c(5, length(q.vec), length(method.name), length(pi0.vec)))
for (j in seq(length(pi0.vec))) {
  for (k in seq(length(method.name))) {
    for (i in seq(nsim)) {
      FDP.array[i, , k, j] <- sapply(q.vec, FDP, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
      pFDP.array[i, , k, j] <- sapply(q.vec, pFDP, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
      power.array[i, , k, j] <- sapply(q.vec, power, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
    }
    FDP.summary[, , k, j] <- rbind(
      avg <- colMeans(FDP.array[, , k, j], na.rm = TRUE),
      sd <- apply(FDP.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(FDP.array[, , k, j])),
      q975 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE),
      q750 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.75, na.rm = TRUE),
      q250 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.25, na.rm = TRUE)
    )
    pFDP.summary[, , k, j] <- rbind(
      avg <- colMeans(pFDP.array[, , k, j], na.rm = TRUE),
      sd <- apply(pFDP.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(pFDP.array[, , k, j])),
      q975 <- apply(pFDP.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(pFDP.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE)
    )
    power.summary[, , k, j] <- rbind(
      avg <- colMeans(power.array[, , k, j], na.rm = TRUE),
      sd <- apply(power.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(power.array[, , k, j])),
      q975 <- apply(power.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(power.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE)
    )
  }
}

for (j in seq(length(pi0.vec))) {
  sd.z <- sapply(z.pi0.list[[j]], sd)
  Noise <- cut(sd.z, breaks = c(0, quantile(sd.z, probs = 1 : 2 / 3), Inf), labels = c("Deflated Noise", "In-between", "Inflated Noise"))
# Noise <- cut(sd.z, breaks = c(0, z.under, z.over, Inf), labels = c("Deflated Noise", "In-between", "Inflated Noise"))

  pi0.pi0 <- matrix(unlist(pi0.pi0.list[[j]]), byrow = TRUE, length(pi0.pi0.list[[j]]))
  pi0.pi0.noise <- rbind.data.frame(cbind.data.frame(Noise, pi0.pi0), cbind.data.frame(Noise = rep("All", length(Noise)), pi0.pi0))
  
  pi0.plot <- ggplot(data = melt(pi0.pi0.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col[-1]) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_hline(yintercept = pi0.vec[j], col = "black", linetype = "dashed", size = 1) +
  scale_x_discrete(labels = method.name[-1]) +
  labs(x = "", y = expression(hat(pi)[0])) +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))

  FDP.summary.pi0 <- aperm(FDP.summary[, , , j], c(2, 1, 3))
  FDP.summary.pi0.method <- FDP.summary.pi0[, , 1]
  for (kk in 2 : length(method.name)) {
    FDP.summary.pi0.method <- rbind.data.frame(FDP.summary.pi0.method, FDP.summary.pi0[, , kk])
  }
  FDP.summary.pi0.method <- cbind.data.frame(
    rep(factor(seq(method.name)), each = dim(FDP.summary.pi0)[1]),
    rep(q.vec, length(method.name)),
    FDP.summary.pi0.method
  )
  colnames(FDP.summary.pi0.method) <- c(
    "Method", "FDR", "FDP", "sd", "n", "q975", "q025", "q750", "q250"
  )
  
  FDP.array.pi0 <- aperm(FDP.array[, , , j], c(2, 1, 3))
  FDP.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, mean, na.rm = TRUE), c(2, 1, 3)))
  sd.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, sd, na.rm = TRUE), c(2, 1, 3)))
  n.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, function(x){sum(!is.na(x))}), c(2, 1, 3)))
  q975.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.975, na.rm = TRUE), c(2, 1, 3)))
  q025.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.025, na.rm = TRUE), c(2, 1, 3)))
  q750.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.75, na.rm = TRUE), c(2, 1, 3)))
  q250.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.25, na.rm = TRUE), c(2, 1, 3)))
  FDP.summary.pi0.method.noise <- cbind.data.frame(
    rep(rep(levels(Noise), each = length(q.vec)), length(method.name)),
    rep(factor(seq(method.name)), each = length(levels(Noise)) * length(q.vec)),
    rep(q.vec, length(levels(Noise)) * length(method.name)),
    FDP.pi0.noise,
    sd.pi0.noise,
    n.pi0.noise,
    q975.pi0.noise,
    q025.pi0.noise,
    q750.pi0.noise,
    q250.pi0.noise
  )
  colnames(FDP.summary.pi0.method.noise) <- c(
    "Noise", "Method", "FDR", "FDP", "sd", "n", "q975", "q025", "q750", "q250"
  )
  FDP.summary.pi0.method.noise <- rbind.data.frame(
    FDP.summary.pi0.method.noise,
    cbind.data.frame(Noise = rep("All", dim(FDP.summary.pi0.method)[1]), FDP.summary.pi0.method)
  )
  
  FDR.calib.plot <- ggplot(data = FDP.summary.pi0.method.noise, aes(x = FDR, y = FDP, group = Method, col = Method)) +
  geom_line() +
  geom_ribbon(aes(ymin = q025, ymax = q975, fill = Method), alpha = 0.35, linetype = "blank") +
  scale_color_manual(labels = method.name, values = method.col) +
  scale_fill_manual(labels = method.name, values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", size = 1, col = "black") +
  labs(x = "Nominal FDR", y = "FDP") +
  theme(legend.position = "top", legend.text = element_text(size = 15), plot.title = element_text(hjust = 0.5, size = 15), axis.title.x = element_text(size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(angle = 45, size = 15), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  FDP.q <- FDP.array[, which(round(q.vec, 4) == q), , j]
  FDP.q.noise <- rbind.data.frame(cbind.data.frame(Noise, FDP.q), cbind.data.frame(Noise = rep("All", length(Noise)), FDP.q))

  FDR.plot <- ggplot(data = melt(FDP.q.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_hline(yintercept = q, col = "black", linetype = "dashed", size = 1) +
  scale_x_discrete(labels = method.name) +
  labs(x = "", y = "FDP") +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  TDP.q <- power.array[, which(round(q.vec, 4) == q), , j]
  TDP.q.noise <- rbind.data.frame(cbind.data.frame(Noise, TDP.q), cbind.data.frame(Noise = rep("All", length(Noise)), TDP.q))

  power.plot <- ggplot(data = melt(TDP.q.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  scale_x_discrete(labels = method.name) +
  labs(x = "", y = "TPP") +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  joint <- ggarrange(FDR.calib.plot,
            pi0.plot + rremove("x.text"),
            FDR.plot + rremove("x.text"),
            power.plot,
            align = "v", ncol = 1, nrow = 4,
            heights = c(1.5, 1, 1, 1.2)
  )
  joint <- annotate_figure(joint,
    top = text_grob(bquote(pi[0] == .(pi0.vec[j])), size = 15)
  )
  print(joint)
  ggsave(paste0("../output/fig/g5_pi0_", pi0.vec[j], ".pdf"), joint, height = 10, width = 8)
}

Big normal: \(g_6 = N\left(0, 5^2\right)\)

FDP.array <- pFDP.array <- power.array <- array(0, dim = c(nsim, length(q.vec), length(method.name), length(pi0.vec)))
FDP.summary <- array(0, dim = c(7, length(q.vec), length(method.name), length(pi0.vec)))
pFDP.summary <- power.summary <- array(0, dim = c(5, length(q.vec), length(method.name), length(pi0.vec)))
for (j in seq(length(pi0.vec))) {
  for (k in seq(length(method.name))) {
    for (i in seq(nsim)) {
      FDP.array[i, , k, j] <- sapply(q.vec, FDP, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
      pFDP.array[i, , k, j] <- sapply(q.vec, pFDP, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
      power.array[i, , k, j] <- sapply(q.vec, power, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
    }
    FDP.summary[, , k, j] <- rbind(
      avg <- colMeans(FDP.array[, , k, j], na.rm = TRUE),
      sd <- apply(FDP.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(FDP.array[, , k, j])),
      q975 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE),
      q750 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.75, na.rm = TRUE),
      q250 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.25, na.rm = TRUE)
    )
    pFDP.summary[, , k, j] <- rbind(
      avg <- colMeans(pFDP.array[, , k, j], na.rm = TRUE),
      sd <- apply(pFDP.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(pFDP.array[, , k, j])),
      q975 <- apply(pFDP.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(pFDP.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE)
    )
    power.summary[, , k, j] <- rbind(
      avg <- colMeans(power.array[, , k, j], na.rm = TRUE),
      sd <- apply(power.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(power.array[, , k, j])),
      q975 <- apply(power.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(power.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE)
    )
  }
}

for (j in seq(length(pi0.vec))) {
  sd.z <- sapply(z.pi0.list[[j]], sd)
  Noise <- cut(sd.z, breaks = c(0, quantile(sd.z, probs = 1 : 2 / 3), Inf), labels = c("Deflated Noise", "In-between", "Inflated Noise"))
# Noise <- cut(sd.z, breaks = c(0, z.under, z.over, Inf), labels = c("Deflated Noise", "In-between", "Inflated Noise"))

  pi0.pi0 <- matrix(unlist(pi0.pi0.list[[j]]), byrow = TRUE, length(pi0.pi0.list[[j]]))
  pi0.pi0.noise <- rbind.data.frame(cbind.data.frame(Noise, pi0.pi0), cbind.data.frame(Noise = rep("All", length(Noise)), pi0.pi0))
  
  pi0.plot <- ggplot(data = melt(pi0.pi0.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col[-1]) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_hline(yintercept = pi0.vec[j], col = "black", linetype = "dashed", size = 1) +
  scale_x_discrete(labels = method.name[-1]) +
  labs(x = "", y = expression(hat(pi)[0])) +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))

  FDP.summary.pi0 <- aperm(FDP.summary[, , , j], c(2, 1, 3))
  FDP.summary.pi0.method <- FDP.summary.pi0[, , 1]
  for (kk in 2 : length(method.name)) {
    FDP.summary.pi0.method <- rbind.data.frame(FDP.summary.pi0.method, FDP.summary.pi0[, , kk])
  }
  FDP.summary.pi0.method <- cbind.data.frame(
    rep(factor(seq(method.name)), each = dim(FDP.summary.pi0)[1]),
    rep(q.vec, length(method.name)),
    FDP.summary.pi0.method
  )
  colnames(FDP.summary.pi0.method) <- c(
    "Method", "FDR", "FDP", "sd", "n", "q975", "q025", "q750", "q250"
  )
  
  FDP.array.pi0 <- aperm(FDP.array[, , , j], c(2, 1, 3))
  FDP.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, mean, na.rm = TRUE), c(2, 1, 3)))
  sd.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, sd, na.rm = TRUE), c(2, 1, 3)))
  n.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, function(x){sum(!is.na(x))}), c(2, 1, 3)))
  q975.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.975, na.rm = TRUE), c(2, 1, 3)))
  q025.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.025, na.rm = TRUE), c(2, 1, 3)))
  q750.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.75, na.rm = TRUE), c(2, 1, 3)))
  q250.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.25, na.rm = TRUE), c(2, 1, 3)))
  FDP.summary.pi0.method.noise <- cbind.data.frame(
    rep(rep(levels(Noise), each = length(q.vec)), length(method.name)),
    rep(factor(seq(method.name)), each = length(levels(Noise)) * length(q.vec)),
    rep(q.vec, length(levels(Noise)) * length(method.name)),
    FDP.pi0.noise,
    sd.pi0.noise,
    n.pi0.noise,
    q975.pi0.noise,
    q025.pi0.noise,
    q750.pi0.noise,
    q250.pi0.noise
  )
  colnames(FDP.summary.pi0.method.noise) <- c(
    "Noise", "Method", "FDR", "FDP", "sd", "n", "q975", "q025", "q750", "q250"
  )
  FDP.summary.pi0.method.noise <- rbind.data.frame(
    FDP.summary.pi0.method.noise,
    cbind.data.frame(Noise = rep("All", dim(FDP.summary.pi0.method)[1]), FDP.summary.pi0.method)
  )
  
  FDR.calib.plot <- ggplot(data = FDP.summary.pi0.method.noise, aes(x = FDR, y = FDP, group = Method, col = Method)) +
  geom_line() +
  geom_ribbon(aes(ymin = q025, ymax = q975, fill = Method), alpha = 0.35, linetype = "blank") +
  scale_color_manual(labels = method.name, values = method.col) +
  scale_fill_manual(labels = method.name, values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", size = 1, col = "black") +
  labs(x = "Nominal FDR", y = "FDP") +
  theme(legend.position = "top", legend.text = element_text(size = 15), plot.title = element_text(hjust = 0.5, size = 15), axis.title.x = element_text(size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(angle = 45, size = 15), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  FDP.q <- FDP.array[, which(round(q.vec, 4) == q), , j]
  FDP.q.noise <- rbind.data.frame(cbind.data.frame(Noise, FDP.q), cbind.data.frame(Noise = rep("All", length(Noise)), FDP.q))

  FDR.plot <- ggplot(data = melt(FDP.q.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_hline(yintercept = q, col = "black", linetype = "dashed", size = 1) +
  scale_x_discrete(labels = method.name) +
  labs(x = "", y = "FDP") +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  TDP.q <- power.array[, which(round(q.vec, 4) == q), , j]
  TDP.q.noise <- rbind.data.frame(cbind.data.frame(Noise, TDP.q), cbind.data.frame(Noise = rep("All", length(Noise)), TDP.q))

  power.plot <- ggplot(data = melt(TDP.q.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  scale_x_discrete(labels = method.name) +
  labs(x = "", y = "TPP") +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  joint <- ggarrange(FDR.calib.plot,
            pi0.plot + rremove("x.text"),
            FDR.plot + rremove("x.text"),
            power.plot,
            align = "v", ncol = 1, nrow = 4,
            heights = c(1.5, 1, 1, 1.2)
  )
  joint <- annotate_figure(joint,
    top = text_grob(bquote(pi[0] == .(pi0.vec[j])), size = 15)
  )
  print(joint)
  ggsave(paste0("../output/fig/g6_pi0_", pi0.vec[j], ".pdf"), joint, height = 10, width = 8)
}

Skew: \(g_7 = 1/4 N\left(-2, 2^2\right) + 1/4 N\left(-1, 1.5^2\right) + 1/3 N\left(0, 1^2\right) + 1 / 6 N\left(1, 1^2\right)\)

FDP.array <- pFDP.array <- power.array <- array(0, dim = c(nsim, length(q.vec), length(method.name), length(pi0.vec)))
FDP.summary <- array(0, dim = c(7, length(q.vec), length(method.name), length(pi0.vec)))
pFDP.summary <- power.summary <- array(0, dim = c(5, length(q.vec), length(method.name), length(pi0.vec)))
for (j in seq(length(pi0.vec))) {
  for (k in seq(length(method.name))) {
    for (i in seq(nsim)) {
      FDP.array[i, , k, j] <- sapply(q.vec, FDP, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
      pFDP.array[i, , k, j] <- sapply(q.vec, pFDP, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
      power.array[i, , k, j] <- sapply(q.vec, power, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
    }
    FDP.summary[, , k, j] <- rbind(
      avg <- colMeans(FDP.array[, , k, j], na.rm = TRUE),
      sd <- apply(FDP.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(FDP.array[, , k, j])),
      q975 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE),
      q750 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.75, na.rm = TRUE),
      q250 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.25, na.rm = TRUE)
    )
    pFDP.summary[, , k, j] <- rbind(
      avg <- colMeans(pFDP.array[, , k, j], na.rm = TRUE),
      sd <- apply(pFDP.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(pFDP.array[, , k, j])),
      q975 <- apply(pFDP.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(pFDP.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE)
    )
    power.summary[, , k, j] <- rbind(
      avg <- colMeans(power.array[, , k, j], na.rm = TRUE),
      sd <- apply(power.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(power.array[, , k, j])),
      q975 <- apply(power.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(power.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE)
    )
  }
}

for (j in seq(length(pi0.vec))) {
  sd.z <- sapply(z.pi0.list[[j]], sd)
  Noise <- cut(sd.z, breaks = c(0, quantile(sd.z, probs = 1 : 2 / 3), Inf), labels = c("Deflated Noise", "In-between", "Inflated Noise"))
# Noise <- cut(sd.z, breaks = c(0, z.under, z.over, Inf), labels = c("Deflated Noise", "In-between", "Inflated Noise"))

  pi0.pi0 <- matrix(unlist(pi0.pi0.list[[j]]), byrow = TRUE, length(pi0.pi0.list[[j]]))
  pi0.pi0.noise <- rbind.data.frame(cbind.data.frame(Noise, pi0.pi0), cbind.data.frame(Noise = rep("All", length(Noise)), pi0.pi0))
  
  pi0.plot <- ggplot(data = melt(pi0.pi0.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col[-1]) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_hline(yintercept = pi0.vec[j], col = "black", linetype = "dashed", size = 1) +
  scale_x_discrete(labels = method.name[-1]) +
  labs(x = "", y = expression(hat(pi)[0])) +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))

  FDP.summary.pi0 <- aperm(FDP.summary[, , , j], c(2, 1, 3))
  FDP.summary.pi0.method <- FDP.summary.pi0[, , 1]
  for (kk in 2 : length(method.name)) {
    FDP.summary.pi0.method <- rbind.data.frame(FDP.summary.pi0.method, FDP.summary.pi0[, , kk])
  }
  FDP.summary.pi0.method <- cbind.data.frame(
    rep(factor(seq(method.name)), each = dim(FDP.summary.pi0)[1]),
    rep(q.vec, length(method.name)),
    FDP.summary.pi0.method
  )
  colnames(FDP.summary.pi0.method) <- c(
    "Method", "FDR", "FDP", "sd", "n", "q975", "q025", "q750", "q250"
  )
  
  FDP.array.pi0 <- aperm(FDP.array[, , , j], c(2, 1, 3))
  FDP.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, mean, na.rm = TRUE), c(2, 1, 3)))
  sd.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, sd, na.rm = TRUE), c(2, 1, 3)))
  n.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, function(x){sum(!is.na(x))}), c(2, 1, 3)))
  q975.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.975, na.rm = TRUE), c(2, 1, 3)))
  q025.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.025, na.rm = TRUE), c(2, 1, 3)))
  q750.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.75, na.rm = TRUE), c(2, 1, 3)))
  q250.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.25, na.rm = TRUE), c(2, 1, 3)))
  FDP.summary.pi0.method.noise <- cbind.data.frame(
    rep(rep(levels(Noise), each = length(q.vec)), length(method.name)),
    rep(factor(seq(method.name)), each = length(levels(Noise)) * length(q.vec)),
    rep(q.vec, length(levels(Noise)) * length(method.name)),
    FDP.pi0.noise,
    sd.pi0.noise,
    n.pi0.noise,
    q975.pi0.noise,
    q025.pi0.noise,
    q750.pi0.noise,
    q250.pi0.noise
  )
  colnames(FDP.summary.pi0.method.noise) <- c(
    "Noise", "Method", "FDR", "FDP", "sd", "n", "q975", "q025", "q750", "q250"
  )
  FDP.summary.pi0.method.noise <- rbind.data.frame(
    FDP.summary.pi0.method.noise,
    cbind.data.frame(Noise = rep("All", dim(FDP.summary.pi0.method)[1]), FDP.summary.pi0.method)
  )
  
  FDR.calib.plot <- ggplot(data = FDP.summary.pi0.method.noise, aes(x = FDR, y = FDP, group = Method, col = Method)) +
  geom_line() +
  geom_ribbon(aes(ymin = q025, ymax = q975, fill = Method), alpha = 0.35, linetype = "blank") +
  scale_color_manual(labels = method.name, values = method.col) +
  scale_fill_manual(labels = method.name, values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", size = 1, col = "black") +
  labs(x = "Nominal FDR", y = "FDP") +
  theme(legend.position = "top", legend.text = element_text(size = 15), plot.title = element_text(hjust = 0.5, size = 15), axis.title.x = element_text(size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(angle = 45, size = 15), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  FDP.q <- FDP.array[, which(round(q.vec, 4) == q), , j]
  FDP.q.noise <- rbind.data.frame(cbind.data.frame(Noise, FDP.q), cbind.data.frame(Noise = rep("All", length(Noise)), FDP.q))

  FDR.plot <- ggplot(data = melt(FDP.q.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_hline(yintercept = q, col = "black", linetype = "dashed", size = 1) +
  scale_x_discrete(labels = method.name) +
  labs(x = "", y = "FDP") +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  TDP.q <- power.array[, which(round(q.vec, 4) == q), , j]
  TDP.q.noise <- rbind.data.frame(cbind.data.frame(Noise, TDP.q), cbind.data.frame(Noise = rep("All", length(Noise)), TDP.q))

  power.plot <- ggplot(data = melt(TDP.q.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  scale_x_discrete(labels = method.name) +
  labs(x = "", y = "TPP") +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  joint <- ggarrange(FDR.calib.plot,
            pi0.plot + rremove("x.text"),
            FDR.plot + rremove("x.text"),
            power.plot,
            align = "v", ncol = 1, nrow = 4,
            heights = c(1.5, 1, 1, 1.2)
  )
  joint <- annotate_figure(joint,
    top = text_grob(bquote(pi[0] == .(pi0.vec[j])), size = 15)
  )
  print(joint)
  ggsave(paste0("../output/fig/g7_pi0_", pi0.vec[j], ".pdf"), joint, height = 10, width = 8)
}

Small normal: \(g_8 = N\left(0, 1^2\right)\)

FDP.array <- pFDP.array <- power.array <- array(0, dim = c(nsim, length(q.vec), length(method.name), length(pi0.vec)))
FDP.summary <- array(0, dim = c(7, length(q.vec), length(method.name), length(pi0.vec)))
pFDP.summary <- power.summary <- array(0, dim = c(5, length(q.vec), length(method.name), length(pi0.vec)))
for (j in seq(length(pi0.vec))) {
  for (k in seq(length(method.name))) {
    for (i in seq(nsim)) {
      FDP.array[i, , k, j] <- sapply(q.vec, FDP, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
      pFDP.array[i, , k, j] <- sapply(q.vec, pFDP, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
      power.array[i, , k, j] <- sapply(q.vec, power, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
    }
    FDP.summary[, , k, j] <- rbind(
      avg <- colMeans(FDP.array[, , k, j], na.rm = TRUE),
      sd <- apply(FDP.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(FDP.array[, , k, j])),
      q975 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE),
      q750 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.75, na.rm = TRUE),
      q250 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.25, na.rm = TRUE)
    )
    pFDP.summary[, , k, j] <- rbind(
      avg <- colMeans(pFDP.array[, , k, j], na.rm = TRUE),
      sd <- apply(pFDP.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(pFDP.array[, , k, j])),
      q975 <- apply(pFDP.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(pFDP.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE)
    )
    power.summary[, , k, j] <- rbind(
      avg <- colMeans(power.array[, , k, j], na.rm = TRUE),
      sd <- apply(power.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(power.array[, , k, j])),
      q975 <- apply(power.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(power.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE)
    )
  }
}

for (j in seq(length(pi0.vec))) {
  sd.z <- sapply(z.pi0.list[[j]], sd)
  Noise <- cut(sd.z, breaks = c(0, quantile(sd.z, probs = 1 : 2 / 3), Inf), labels = c("Deflated Noise", "In-between", "Inflated Noise"))
# Noise <- cut(sd.z, breaks = c(0, z.under, z.over, Inf), labels = c("Deflated Noise", "In-between", "Inflated Noise"))

  pi0.pi0 <- matrix(unlist(pi0.pi0.list[[j]]), byrow = TRUE, length(pi0.pi0.list[[j]]))
  pi0.pi0.noise <- rbind.data.frame(cbind.data.frame(Noise, pi0.pi0), cbind.data.frame(Noise = rep("All", length(Noise)), pi0.pi0))
  
  pi0.plot <- ggplot(data = melt(pi0.pi0.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col[-1]) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_hline(yintercept = pi0.vec[j], col = "black", linetype = "dashed", size = 1) +
  scale_x_discrete(labels = method.name[-1]) +
  labs(x = "", y = expression(hat(pi)[0])) +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))

  FDP.summary.pi0 <- aperm(FDP.summary[, , , j], c(2, 1, 3))
  FDP.summary.pi0.method <- FDP.summary.pi0[, , 1]
  for (kk in 2 : length(method.name)) {
    FDP.summary.pi0.method <- rbind.data.frame(FDP.summary.pi0.method, FDP.summary.pi0[, , kk])
  }
  FDP.summary.pi0.method <- cbind.data.frame(
    rep(factor(seq(method.name)), each = dim(FDP.summary.pi0)[1]),
    rep(q.vec, length(method.name)),
    FDP.summary.pi0.method
  )
  colnames(FDP.summary.pi0.method) <- c(
    "Method", "FDR", "FDP", "sd", "n", "q975", "q025", "q750", "q250"
  )
  
  FDP.array.pi0 <- aperm(FDP.array[, , , j], c(2, 1, 3))
  FDP.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, mean, na.rm = TRUE), c(2, 1, 3)))
  sd.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, sd, na.rm = TRUE), c(2, 1, 3)))
  n.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, function(x){sum(!is.na(x))}), c(2, 1, 3)))
  q975.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.975, na.rm = TRUE), c(2, 1, 3)))
  q025.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.025, na.rm = TRUE), c(2, 1, 3)))
  q750.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.75, na.rm = TRUE), c(2, 1, 3)))
  q250.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.25, na.rm = TRUE), c(2, 1, 3)))
  FDP.summary.pi0.method.noise <- cbind.data.frame(
    rep(rep(levels(Noise), each = length(q.vec)), length(method.name)),
    rep(factor(seq(method.name)), each = length(levels(Noise)) * length(q.vec)),
    rep(q.vec, length(levels(Noise)) * length(method.name)),
    FDP.pi0.noise,
    sd.pi0.noise,
    n.pi0.noise,
    q975.pi0.noise,
    q025.pi0.noise,
    q750.pi0.noise,
    q250.pi0.noise
  )
  colnames(FDP.summary.pi0.method.noise) <- c(
    "Noise", "Method", "FDR", "FDP", "sd", "n", "q975", "q025", "q750", "q250"
  )
  FDP.summary.pi0.method.noise <- rbind.data.frame(
    FDP.summary.pi0.method.noise,
    cbind.data.frame(Noise = rep("All", dim(FDP.summary.pi0.method)[1]), FDP.summary.pi0.method)
  )
  
  FDR.calib.plot <- ggplot(data = FDP.summary.pi0.method.noise, aes(x = FDR, y = FDP, group = Method, col = Method)) +
  geom_line() +
  geom_ribbon(aes(ymin = q025, ymax = q975, fill = Method), alpha = 0.35, linetype = "blank") +
  scale_color_manual(labels = method.name, values = method.col) +
  scale_fill_manual(labels = method.name, values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", size = 1, col = "black") +
  labs(x = "Nominal FDR", y = "FDP") +
  theme(legend.position = "top", legend.text = element_text(size = 15), plot.title = element_text(hjust = 0.5, size = 15), axis.title.x = element_text(size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(angle = 45, size = 15), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  FDP.q <- FDP.array[, which(round(q.vec, 4) == q), , j]
  FDP.q.noise <- rbind.data.frame(cbind.data.frame(Noise, FDP.q), cbind.data.frame(Noise = rep("All", length(Noise)), FDP.q))

  FDR.plot <- ggplot(data = melt(FDP.q.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_hline(yintercept = q, col = "black", linetype = "dashed", size = 1) +
  scale_x_discrete(labels = method.name) +
  labs(x = "", y = "FDP") +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  TDP.q <- power.array[, which(round(q.vec, 4) == q), , j]
  TDP.q.noise <- rbind.data.frame(cbind.data.frame(Noise, TDP.q), cbind.data.frame(Noise = rep("All", length(Noise)), TDP.q))

  power.plot <- ggplot(data = melt(TDP.q.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  scale_x_discrete(labels = method.name) +
  labs(x = "", y = "TPP") +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  joint <- ggarrange(FDR.calib.plot,
            pi0.plot + rremove("x.text"),
            FDR.plot + rremove("x.text"),
            power.plot,
            align = "v", ncol = 1, nrow = 4,
            heights = c(1.5, 1, 1, 1.2)
  )
  joint <- annotate_figure(joint,
    top = text_grob(bquote(pi[0] == .(pi0.vec[j])), size = 15)
  )
  print(joint)
  ggsave(paste0("../output/fig/g8_pi0_", pi0.vec[j], ".pdf"), joint, height = 10, width = 8)
}

Small normal: \(g_9 = N\left(0, 1.5^2\right)\)

FDP.array <- pFDP.array <- power.array <- array(0, dim = c(nsim, length(q.vec), length(method.name), length(pi0.vec)))
FDP.summary <- array(0, dim = c(7, length(q.vec), length(method.name), length(pi0.vec)))
pFDP.summary <- power.summary <- array(0, dim = c(5, length(q.vec), length(method.name), length(pi0.vec)))
for (j in seq(length(pi0.vec))) {
  for (k in seq(length(method.name))) {
    for (i in seq(nsim)) {
      FDP.array[i, , k, j] <- sapply(q.vec, FDP, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
      pFDP.array[i, , k, j] <- sapply(q.vec, pFDP, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
      power.array[i, , k, j] <- sapply(q.vec, power, qvalue = qvalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]])
    }
    FDP.summary[, , k, j] <- rbind(
      avg <- colMeans(FDP.array[, , k, j], na.rm = TRUE),
      sd <- apply(FDP.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(FDP.array[, , k, j])),
      q975 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE),
      q750 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.75, na.rm = TRUE),
      q250 <- apply(FDP.array[, , k, j], 2, quantile, probs = 0.25, na.rm = TRUE)
    )
    pFDP.summary[, , k, j] <- rbind(
      avg <- colMeans(pFDP.array[, , k, j], na.rm = TRUE),
      sd <- apply(pFDP.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(pFDP.array[, , k, j])),
      q975 <- apply(pFDP.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(pFDP.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE)
    )
    power.summary[, , k, j] <- rbind(
      avg <- colMeans(power.array[, , k, j], na.rm = TRUE),
      sd <- apply(power.array[, , k, j], 2, sd, na.rm = TRUE),
      n <- colSums(!is.na(power.array[, , k, j])),
      q975 <- apply(power.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
      q025 <- apply(power.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE)
    )
  }
}

for (j in seq(length(pi0.vec))) {
  sd.z <- sapply(z.pi0.list[[j]], sd)
  Noise <- cut(sd.z, breaks = c(0, quantile(sd.z, probs = 1 : 2 / 3), Inf), labels = c("Deflated Noise", "In-between", "Inflated Noise"))
# Noise <- cut(sd.z, breaks = c(0, z.under, z.over, Inf), labels = c("Deflated Noise", "In-between", "Inflated Noise"))

  pi0.pi0 <- matrix(unlist(pi0.pi0.list[[j]]), byrow = TRUE, length(pi0.pi0.list[[j]]))
  pi0.pi0.noise <- rbind.data.frame(cbind.data.frame(Noise, pi0.pi0), cbind.data.frame(Noise = rep("All", length(Noise)), pi0.pi0))
  
  pi0.plot <- ggplot(data = melt(pi0.pi0.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col[-1]) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_hline(yintercept = pi0.vec[j], col = "black", linetype = "dashed", size = 1) +
  scale_x_discrete(labels = method.name[-1]) +
  labs(x = "", y = expression(hat(pi)[0])) +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))

  FDP.summary.pi0 <- aperm(FDP.summary[, , , j], c(2, 1, 3))
  FDP.summary.pi0.method <- FDP.summary.pi0[, , 1]
  for (kk in 2 : length(method.name)) {
    FDP.summary.pi0.method <- rbind.data.frame(FDP.summary.pi0.method, FDP.summary.pi0[, , kk])
  }
  FDP.summary.pi0.method <- cbind.data.frame(
    rep(factor(seq(method.name)), each = dim(FDP.summary.pi0)[1]),
    rep(q.vec, length(method.name)),
    FDP.summary.pi0.method
  )
  colnames(FDP.summary.pi0.method) <- c(
    "Method", "FDR", "FDP", "sd", "n", "q975", "q025", "q750", "q250"
  )
  
  FDP.array.pi0 <- aperm(FDP.array[, , , j], c(2, 1, 3))
  FDP.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, mean, na.rm = TRUE), c(2, 1, 3)))
  sd.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, sd, na.rm = TRUE), c(2, 1, 3)))
  n.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, function(x){sum(!is.na(x))}), c(2, 1, 3)))
  q975.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.975, na.rm = TRUE), c(2, 1, 3)))
  q025.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.025, na.rm = TRUE), c(2, 1, 3)))
  q750.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.75, na.rm = TRUE), c(2, 1, 3)))
  q250.pi0.noise <- as.vector(aperm(apply(FDP.array.pi0, c(1, 3), tapply, Noise, quantile, probs = 0.25, na.rm = TRUE), c(2, 1, 3)))
  FDP.summary.pi0.method.noise <- cbind.data.frame(
    rep(rep(levels(Noise), each = length(q.vec)), length(method.name)),
    rep(factor(seq(method.name)), each = length(levels(Noise)) * length(q.vec)),
    rep(q.vec, length(levels(Noise)) * length(method.name)),
    FDP.pi0.noise,
    sd.pi0.noise,
    n.pi0.noise,
    q975.pi0.noise,
    q025.pi0.noise,
    q750.pi0.noise,
    q250.pi0.noise
  )
  colnames(FDP.summary.pi0.method.noise) <- c(
    "Noise", "Method", "FDR", "FDP", "sd", "n", "q975", "q025", "q750", "q250"
  )
  FDP.summary.pi0.method.noise <- rbind.data.frame(
    FDP.summary.pi0.method.noise,
    cbind.data.frame(Noise = rep("All", dim(FDP.summary.pi0.method)[1]), FDP.summary.pi0.method)
  )
  
  FDR.calib.plot <- ggplot(data = FDP.summary.pi0.method.noise, aes(x = FDR, y = FDP, group = Method, col = Method)) +
  geom_line() +
  geom_ribbon(aes(ymin = q025, ymax = q975, fill = Method), alpha = 0.35, linetype = "blank") +
  scale_color_manual(labels = method.name, values = method.col) +
  scale_fill_manual(labels = method.name, values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", size = 1, col = "black") +
  labs(x = "Nominal FDR", y = "FDP") +
  theme(legend.position = "top", legend.text = element_text(size = 15), plot.title = element_text(hjust = 0.5, size = 15), axis.title.x = element_text(size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(angle = 45, size = 15), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  FDP.q <- FDP.array[, which(round(q.vec, 4) == q), , j]
  FDP.q.noise <- rbind.data.frame(cbind.data.frame(Noise, FDP.q), cbind.data.frame(Noise = rep("All", length(Noise)), FDP.q))

  FDR.plot <- ggplot(data = melt(FDP.q.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  geom_hline(yintercept = q, col = "black", linetype = "dashed", size = 1) +
  scale_x_discrete(labels = method.name) +
  labs(x = "", y = "FDP") +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  TDP.q <- power.array[, which(round(q.vec, 4) == q), , j]
  TDP.q.noise <- rbind.data.frame(cbind.data.frame(Noise, TDP.q), cbind.data.frame(Noise = rep("All", length(Noise)), TDP.q))

  power.plot <- ggplot(data = melt(TDP.q.noise, id.vars = "Noise"), aes(x = variable, y = value, col = variable)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 3) +
  scale_color_manual(values = method.col) +
  facet_wrap(~Noise, nrow = 1, ncol = 4) +
  scale_x_discrete(labels = method.name) +
  labs(x = "", y = "TPP") +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5, size = 15), axis.title.y = element_text(size = 15), axis.text.x = element_text(size = 15, angle = 45, hjust = 1), axis.text.y = element_text(size = 15), strip.text = element_text(size = 15))
  
  joint <- ggarrange(FDR.calib.plot,
            pi0.plot + rremove("x.text"),
            FDR.plot + rremove("x.text"),
            power.plot,
            align = "v", ncol = 1, nrow = 4,
            heights = c(1.5, 1, 1, 1.2)
  )
  joint <- annotate_figure(joint,
    top = text_grob(bquote(pi[0] == .(pi0.vec[j])), size = 15)
  )
  print(joint)
  ggsave(paste0("../output/fig/g9_pi0_", pi0.vec[j], ".pdf"), joint, height = 10, width = 8)
}

Session information

sessionInfo()

R version 3.4.3 (2017-11-30)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.4

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] ggpubr_0.1.6       magrittr_1.5       cowplot_0.9.2     
 [4] scales_0.5.0       RColorBrewer_1.1-2 gridExtra_2.3     
 [7] ggplot2_2.2.1      reshape2_1.4.3     qvalue_2.10.0     
[10] locfdr_1.1-8       ashr_2.2-2         Rmosek_8.0.69     
[13] CVXR_0.95          REBayes_1.2        Matrix_1.2-12     
[16] SQUAREM_2017.10-1  EQL_1.0-0          ttutils_1.0-1     
[19] PolynomF_1.0-1    

loaded via a namespace (and not attached):
 [1] purrr_0.2.4       splines_3.4.3     lattice_0.20-35  
 [4] colorspace_1.3-2  htmltools_0.3.6   yaml_2.1.18      
 [7] gmp_0.5-13.1      rlang_0.1.6       R.oo_1.21.0      
[10] pillar_1.0.1      glue_1.2.0        Rmpfr_0.6-1      
[13] R.utils_2.6.0     bit64_0.9-7       bindrcpp_0.2     
[16] bindr_0.1         scs_1.1-1         foreach_1.4.4    
[19] plyr_1.8.4        stringr_1.3.0     munsell_0.4.3    
[22] gtable_0.2.0      workflowr_1.0.1   R.methodsS3_1.7.1
[25] codetools_0.2-15  evaluate_0.10.1   labeling_0.3     
[28] knitr_1.20        doParallel_1.0.11 pscl_1.5.2       
[31] parallel_3.4.3    Rcpp_0.12.16      backports_1.1.2  
[34] truncnorm_1.0-7   bit_1.1-12        digest_0.6.15    
[37] stringi_1.1.6     dplyr_0.7.4       rprojroot_1.3-2  
[40] ECOSolveR_0.4     tools_3.4.3       lazyeval_0.2.1   
[43] tibble_1.4.1      pkgconfig_2.0.1   whisker_0.3-2    
[46] MASS_7.3-47       assertthat_0.2.0  rmarkdown_1.9    
[49] iterators_1.0.9   R6_2.2.2          git2r_0.21.0     
[52] compiler_3.4.3