library(mvtnorm); library(quantreg)
p <- 200; n <- 100; h <- 3 # the number of variables, samples and factors
berlii <- rbinom(p, 1, 0.2) # 1 means the variable is non-null and 0 means it is null.
index0 <- which(berlii == 0); index1 <- which(berlii == 1)
mu <- matrix(rep(0, 1*p), nrow=p)
mu[index1] <- runif(length(index1), min=0.4, max=0.7) # expectation of data
B <- matrix(runif(h*p, min=-1, max=1), nrow=p) # factor loading matrix
t_error <- t(rmvt(n, sigma = diag(p), df = 10)) # error term followed t-distribution
f <- t(rmvt(n, diag(h), df = 4))/sqrt(4/(4-2)) # factor followed t-distribution
Y <- mu %*% matrix(rep(1, n*1), nrow=1) + B %*% f + t_error # data
res <- ACE(Z = Y, H0_indicator = berlii, gama = 0.05)
res$FDP # true FDP
res$Power # power
Run the code above in your browser using DataLab