# NOT RUN {
## this is the simulation example in Wang et al. (2015).
n <- 100
p <- 1000
r <- 2
set.seed(1)
data <- gen.sim.data(n = n, p = p, r = r,
alpha = rep(1 / sqrt(r), r),
beta.strength = 3 * sqrt(1 + 1) / sqrt(n),
Gamma.strength = c(seq(3, 1, length = r)) * sqrt(p),
Sigma = 1 / rgamma(p, 3, rate = 2),
beta.nonzero.frac = 0.05)
X.data <- data.frame(X1 = data$X1)
sva.results <- sva.wrapper(~ X1 | 1, X.data, data$Y,
r = r, sva.method = "irw")
ruv.results <- ruv.wrapper(~ X1 | 1, X.data, data$Y, r = r,
nc = sample(data$beta.zero.pos, 30), ruv.method = "RUV4")
leapp.results <- leapp.wrapper(~ X1 | 1, X.data, data$Y, r = r)
cate.results <- cate(~ X1 | 1, X.data, data$Y, r = r)
## p-values after adjustment
par(mfrow = c(2, 2))
hist(sva.results$beta.p.value)
hist(ruv.results$beta.p.value)
hist(leapp.results$beta.p.value)
hist(cate.results$beta.p.value)
## type I error
mean(sva.results$beta.p.value[data$beta.zero.pos] < 0.05)
## power
mean(sva.results$beta.p.value[data$beta.nonzero.pos] < 0.05)
## false discovery proportion for sva
discoveries.sva <- which(p.adjust(sva.results$beta.p.value, "BH") < 0.2)
fdp.sva <- length(setdiff(discoveries.sva, data$beta.nonzero.pos)) / max(length(discoveries.sva), 1)
fdp.sva
# }
Run the code above in your browser using DataLab