CASH Simulations
Lei Sun
2018-02-10
Last updated: 2018-05-18
workflowr checks: (Click a bullet for more information)-
✔ R Markdown file: up-to-date
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
-
✔ Environment: empty
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
-
✔ Seed:
set.seed(12345)The command
set.seed(12345)was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible. -
✔ Session information: recorded
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
-
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility. The version displayed above was the version of the Git repository at the time these results were generated.✔ Repository version: 7c1e2f8
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can usewflow_publishorwflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.Ignored files: Ignored: .DS_Store Ignored: .Rhistory Ignored: .Rproj.user/ Ignored: analysis/.DS_Store Ignored: analysis/BH_robustness_cache/ Ignored: analysis/FDR_Null_cache/ Ignored: analysis/FDR_null_betahat_cache/ Ignored: analysis/Rmosek_cache/ Ignored: analysis/StepDown_cache/ Ignored: analysis/alternative2_cache/ Ignored: analysis/alternative_cache/ Ignored: analysis/ash_gd_cache/ Ignored: analysis/average_cor_gtex_2_cache/ Ignored: analysis/average_cor_gtex_cache/ Ignored: analysis/brca_cache/ Ignored: analysis/cash_deconv_cache/ Ignored: analysis/cash_fdr_1_cache/ Ignored: analysis/cash_fdr_2_cache/ Ignored: analysis/cash_fdr_3_cache/ Ignored: analysis/cash_fdr_4_cache/ Ignored: analysis/cash_fdr_5_cache/ Ignored: analysis/cash_fdr_6_cache/ Ignored: analysis/cash_plots_2_cache/ Ignored: analysis/cash_plots_cache/ Ignored: analysis/cash_sim_1_cache/ Ignored: analysis/cash_sim_2_cache/ Ignored: analysis/cash_sim_3_cache/ Ignored: analysis/cash_sim_4_cache/ Ignored: analysis/cash_sim_5_cache/ Ignored: analysis/cash_sim_6_cache/ Ignored: analysis/cash_sim_7_cache/ Ignored: analysis/correlated_z_2_cache/ Ignored: analysis/correlated_z_3_cache/ Ignored: analysis/correlated_z_cache/ Ignored: analysis/create_null_cache/ Ignored: analysis/cutoff_null_cache/ Ignored: analysis/design_matrix_2_cache/ Ignored: analysis/design_matrix_cache/ Ignored: analysis/diagnostic_ash_cache/ Ignored: analysis/diagnostic_correlated_z_2_cache/ Ignored: analysis/diagnostic_correlated_z_3_cache/ Ignored: analysis/diagnostic_correlated_z_cache/ Ignored: analysis/diagnostic_plot_2_cache/ Ignored: analysis/diagnostic_plot_cache/ Ignored: analysis/efron_leukemia_cache/ Ignored: analysis/fitting_normal_cache/ Ignored: analysis/gaussian_derivatives_2_cache/ Ignored: analysis/gaussian_derivatives_3_cache/ Ignored: analysis/gaussian_derivatives_4_cache/ Ignored: analysis/gaussian_derivatives_5_cache/ Ignored: analysis/gaussian_derivatives_cache/ Ignored: analysis/gd-ash_cache/ Ignored: analysis/gd_delta_cache/ Ignored: analysis/gd_lik_2_cache/ Ignored: analysis/gd_lik_cache/ Ignored: analysis/gd_w_cache/ Ignored: analysis/knockoff_10_cache/ Ignored: analysis/knockoff_2_cache/ Ignored: analysis/knockoff_3_cache/ Ignored: analysis/knockoff_4_cache/ Ignored: analysis/knockoff_5_cache/ Ignored: analysis/knockoff_6_cache/ Ignored: analysis/knockoff_7_cache/ Ignored: analysis/knockoff_8_cache/ Ignored: analysis/knockoff_9_cache/ Ignored: analysis/knockoff_cache/ Ignored: analysis/knockoff_var_cache/ Ignored: analysis/marginal_z_alternative_cache/ Ignored: analysis/marginal_z_cache/ Ignored: analysis/mosek_reg_2_cache/ Ignored: analysis/mosek_reg_4_cache/ Ignored: analysis/mosek_reg_5_cache/ Ignored: analysis/mosek_reg_6_cache/ Ignored: analysis/mosek_reg_cache/ Ignored: analysis/pihat0_null_cache/ Ignored: analysis/plot_diagnostic_cache/ Ignored: analysis/poster_obayes17_cache/ Ignored: analysis/real_data_simulation_2_cache/ Ignored: analysis/real_data_simulation_3_cache/ Ignored: analysis/real_data_simulation_4_cache/ Ignored: analysis/real_data_simulation_5_cache/ Ignored: analysis/real_data_simulation_cache/ Ignored: analysis/rmosek_primal_dual_2_cache/ Ignored: analysis/rmosek_primal_dual_cache/ Ignored: analysis/seqgendiff_cache/ Ignored: analysis/simulated_correlated_null_2_cache/ Ignored: analysis/simulated_correlated_null_3_cache/ Ignored: analysis/simulated_correlated_null_cache/ Ignored: analysis/simulation_real_se_2_cache/ Ignored: analysis/simulation_real_se_cache/ Ignored: analysis/smemo_2_cache/ Ignored: data/LSI/ Ignored: docs/.DS_Store Ignored: docs/figure/.DS_Store Ignored: output/fig/
Expand here to see past versions:
source("../code/gdfit.R")
source("../code/gdash_lik.R")
source("../code/count_to_summary.R")
library(ashr)
library(locfdr)
library(qvalue)
library(reshape2)
library(ggplot2)
library(grid)
library(gridExtra)
library(RColorBrewer)
library(scales)
library(cowplot)
library(ggpubr)mean_sdp <- function (x) {
m <- mean(x)
ymax <- m + sd(x)
return(c(y = m, ymax = ymax, ymin = m))
}
mad.mean <- function (x) {
return(mean(abs(x - median(x))))
}
FDP <- function (FDR, qvalue, beta) {
return(sum(qvalue <= FDR & beta == 0) / max(sum(qvalue <= FDR), 1))
}
pFDP <- function (FDR, qvalue, beta) {
return(sum(qvalue <= FDR & beta == 0) / sum(qvalue <= FDR))
}
power <- function (FDR, qvalue, beta) {
return(sum(qvalue <= FDR & beta != 0) / sum(beta != 0))
}
FSP <- function (FSR, svalue, beta, betahat) {
return(sum(sign(betahat[svalue <= FSR]) != sign(beta[svalue <= FSR])) / max(sum(svalue <= FSR), 1))
}r <- readRDS("../data/liver.rds")top_genes_index = function (g, X) {
return(order(rowSums(X), decreasing = TRUE)[1 : g])
}
lcpm = function (r) {
R = colSums(r)
t(log2(((t(r) + 0.5) / (R + 1)) * 10^6))
}nsamp <- 5
ngene <- 1e4
pi0.vec <- c(0.5, 0.9, 0.99)Y = lcpm(r)
subset = top_genes_index(ngene, Y)
r = r[subset,]Bimodal: \(g_2 = 0.5 N\left(-2, 1\right) + 0.5 N\left(2, 1\right)\)
q.vec <- seq(0.001, 0.20, by = 0.001)
method.name <- c("BHq", "qvalue", "locfdr", "ASH.Q", "ASH.S", "CASH.Q", "CASH.S")
method.name.S <- c("ASH", "CASH")FDP.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)))
power.summary <- array(0, dim = c(5, length(q.vec), length(method.name), length(pi0.vec)))
FSP.array <- array(0, dim = c(nsim, length(q.vec), 2, length(pi0.vec)))
FSP.summary <- array(0, dim = c(7, length(q.vec), 2, 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]])
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)
)
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))) {
for (k in seq(2)) {
for (i in seq(nsim)) {
FSP.array[i, , k, j] <- sapply(q.vec, FSP, svalue = svalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]], betahat = betahat.pi0.list[[j]][[i]])
}
FSP.summary[, , k, j] <- rbind(
avg <- colMeans(FSP.array[, , k, j], na.rm = TRUE),
sd <- apply(FSP.array[, , k, j], 2, sd, na.rm = TRUE),
n <- colSums(!is.na(FSP.array[, , k, j])),
q975 <- apply(FSP.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
q025 <- apply(FSP.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE),
q750 <- apply(FSP.array[, , k, j], 2, quantile, probs = 0.75, na.rm = TRUE),
q250 <- apply(FSP.array[, , k, j], 2, quantile, probs = 0.25, na.rm = TRUE)
)
}
}q <- 0.1
method.col <- scales::hue_pal()(length(method.name))
method.col.S <- method.col[c(5, 7)]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))
FSP.q <- FSP.array[, which(round(q.vec, 4) == q), , j]
FSP.q.noise <- rbind.data.frame(cbind.data.frame(Noise, FSP.q), cbind.data.frame(Noise = rep("All", length(Noise)), FSP.q))
FSR.plot <- ggplot(data = melt(FSP.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.S) +
facet_wrap(~Noise, nrow = 1, ncol = 4) +
geom_hline(yintercept = q, col = "black", linetype = "dashed", size = 1) +
scale_x_discrete(labels = method.name.S) +
labs(x = "", y = "FSP") +
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 + rremove("x.text"),
FSR.plot,
align = "v", ncol = 1, nrow = 5,
heights = c(1.5, 1, 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)
}
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 <- 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)))
power.summary <- array(0, dim = c(5, length(q.vec), length(method.name), length(pi0.vec)))
FSP.array <- array(0, dim = c(nsim, length(q.vec), 2, length(pi0.vec)))
FSP.summary <- array(0, dim = c(7, length(q.vec), 2, 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]])
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)
)
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))) {
for (k in seq(2)) {
for (i in seq(nsim)) {
FSP.array[i, , k, j] <- sapply(q.vec, FSP, svalue = svalue.pi0.list[[j]][[i]][, k], beta = beta.pi0.list[[j]][[i]], betahat = betahat.pi0.list[[j]][[i]])
}
FSP.summary[, , k, j] <- rbind(
avg <- colMeans(FSP.array[, , k, j], na.rm = TRUE),
sd <- apply(FSP.array[, , k, j], 2, sd, na.rm = TRUE),
n <- colSums(!is.na(FSP.array[, , k, j])),
q975 <- apply(FSP.array[, , k, j], 2, quantile, probs = 0.975, na.rm = TRUE),
q025 <- apply(FSP.array[, , k, j], 2, quantile, probs = 0.025, na.rm = TRUE),
q750 <- apply(FSP.array[, , k, j], 2, quantile, probs = 0.75, na.rm = TRUE),
q250 <- apply(FSP.array[, , k, j], 2, quantile, probs = 0.25, 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))
FSP.q <- FSP.array[, which(round(q.vec, 4) == q), , j]
FSP.q.noise <- rbind.data.frame(cbind.data.frame(Noise, FSP.q), cbind.data.frame(Noise = rep("All", length(Noise)), FSP.q))
FSR.plot <- ggplot(data = melt(FSP.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.S) +
facet_wrap(~Noise, nrow = 1, ncol = 4) +
geom_hline(yintercept = q, col = "black", linetype = "dashed", size = 1) +
scale_x_discrete(labels = method.name.S) +
labs(x = "", y = "FSP") +
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 + rremove("x.text"),
FSR.plot,
align = "v", ncol = 1, nrow = 5,
heights = c(1.5, 1, 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)
}
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-7 Rmosek_8.0.69
[13] CVXR_0.95 REBayes_1.3 Matrix_1.2-14
[16] SQUAREM_2017.10-1 EQL_1.0-0 ttutils_1.0-1
[19] PolynomF_1.0-2
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.19
[7] gmp_0.5-13.1 rlang_0.2.0 R.oo_1.22.0
[10] pillar_1.2.2 glue_1.2.0 Rmpfr_0.7-0
[13] R.utils_2.6.0 bit64_0.9-7 bindrcpp_0.2.2
[16] bindr_0.1.1 scs_1.1-1 foreach_1.4.4
[19] plyr_1.8.4 stringr_1.3.1 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-8 bit_1.1-13 digest_0.6.15
[37] stringi_1.2.2 dplyr_0.7.4 rprojroot_1.3-2
[40] ECOSolveR_0.4 tools_3.4.3 lazyeval_0.2.1
[43] tibble_1.4.2 pkgconfig_2.0.1 whisker_0.3-2
[46] MASS_7.3-50 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 This reproducible R Markdown analysis was created with workflowr 1.0.1