########## generate random group transformation
#' generate random group transformation
#' @return ranGT a random group transformation
#' @export
genGT <- function(n) {
ranGT = rbinom(n, 1, 0.5) * 2 - 1
return(ranGT)
}
#' generate k sample permutation
#' @return a list of k component
#' @export
kSamplePermutation <- function(k, Ns,ran =sample.int(sum(Ns))) {
ans <- 0
out <- list()
for (i in 1:k) {
name <- paste0('S', i)
out[[name]] <- ran[(ans + 1):(ans + Ns[i])]
ans <- ans + Ns[i]
}
return(out)
}
#' randomization test
#' @param initGT the identity of group transformation
#' @param genGT a function to generate random group transformation
#' @param GTStat a function of random group trasformation. It calculate the randomization test statistic
#' @param alpha test level
#' @param B randomization times
#' @export
randomizationTest <- function(GTStat,
genGT,
theTestStat,
B = 100,
alpha = 0.05,
needPvalue = FALSE) {
theSum <- 0
if (needPvalue) {
for (i in 1:B) {
if (GTStat(genGT()) >= theTestStat) {
theSum <- theSum + 1
}
}
theSum <- theSum + 1
pvalue <- theSum / (B + 1)
return(pvalue)
}
else{
for (i in 1:B) {
if (GTStat(genGT()) >= theTestStat) {
theSum <- theSum + 1
if (theSum > (B + 1) * alpha - 1) {
return(1)
}
}
else{
if (i - theSum >= B - (B + 1) * alpha + 1) {
return(0)
}
}
}
if (theSum > (B + 1) * alpha - 1) {
return(1)
}
if (i - theSum >= B - (B + 1) * alpha + 1) {
return(0)
}
}
}
#' simulation of randomization test, bootstrap test and CQ test
#' @param n n
#' @param p p
#' @param M simulate M times test procedure to compute empirical power
#' @param mu mean value
#' @param B resample B times to get randomization test
#' @param SNR the signal noise ratio we need, mean value will scale to fit this SNR
#' @export
simRanCQ <- function(n,
p,
M,
mu,
sparse,
B,
SNR,
innov = list("type" = "MA", "innov" = "normal"),
maParam = NULL,
factorParam = NULL,
compoundParam = NULL,
needPvalue = FALSE) {
if (innov$type == "MA") {
popVar <- varMovingAverage(p, maParam$k, maParam$rho)
}
else if(innov$type == "factor"){
popVar <- varFanFactorModel(p, factorParam)
}
else if(innov$type == "compound"){
popVar <- varCompoundModel(p,compoundParam)
}
trSquaredPopVar <- sum(popVar ^ 2)
if (sparse) {
mu[-sample(p, floor(p * 0.05))] <- 0
}
theSNR <-
sqrt(n * (n - 1)) * sum(mu ^ 2) / sqrt(2 * trSquaredPopVar)
# scale the mean
mu <- mu / sqrt(theSNR) * sqrt(SNR)
# the SNR
theNewSNR <-
sqrt(n * (n - 1)) * sum(mu ^ 2) / sqrt(2 * trSquaredPopVar)
thePvalue <- rep(0, M)
CQPvalue <- rep(0, M)
bootstrapPvalue <- rep(0, M)
pb <- txtProgressBar(style=3)
for (i in 1:M) {
if (innov$type == "MA") {
theData <- movingAverageGen(n, p, maParam, mu, innov = innov$innov)
}
else if (innov$type == "factor") {
theData <- fanFactorModelGen(n, p, theParam = factorParam, mu = mu)
}
else if (innov$type == "compound"){
theData <- compoundModelGen(n,p, compoundParam, mu= mu)
}
# randomization p value
tmpGTStat <- chenTempStatistic(theData)
thePvalue[i] <- randomizationTest(
B = B,
GTStat = tmpGTStat,
genGT = function() {
ranGT <- rbinom(n, 1, 0.5) * 2 - 1
return(ranGT)
},
theTestStat = tmpGTStat(rep(1,n)),
alpha = 0.05,
needPvalue = needPvalue
)
# Bootstrap p value
bootstrapPvalue[i] <- randomizationTest(
B = B,
GTStat = bootstrapTempStatistic(theData),
genGT = function() {
ranGT <- sample(n,replace=TRUE)
return(ranGT)
},
theTestStat = tmpGTStat(rep(1,n)),
alpha = 0.05,
needPvalue = needPvalue
)
# asymptotic p value
CQPvalue[i] <- CQTest(theData, trSquaredPopVar)
setTxtProgressBar(pb,i/M)
}
close(pb)
theResult <- list()
theResult$n <- n
theResult$p <- p
theResult$k <- k
theResult$SNR <- theNewSNR
#theResult$theoreticalPower <- pnorm(qnorm(0.05) + theNewSNR)
theResult$randomizationPower <- mean(thePvalue <= 0.05)
theResult$CQPower <- mean(CQPvalue <= 0.05)
theResult$bootstrapPower <- mean(bootstrapPvalue <= 0.05)
unlist(theResult)
}
#' simulation of two sample permutation test for spiked model
#' @param n n
#' @param p p
#' @param M simulate M times test procedure to compute empirical power
#' @param B resample B times to get randomization test
#' @param SNR the signal noise ratio we need, mean value will scale to fit this SNR
#' @param modelParam the parameters of spiked model, including theSigmaSq, V, D, r
#' @export
simTwoSample <- function(n1,
n2,
p,
M,
#mu,
sparse,
B,
SNR,
modelParam,
needPvalue = FALSE) {
tau <- 1 / n1 + 1 / n2
if (sparse) {
mu <- runif(p, 2, 3)
mu[-sample(p, floor(p * 0.05))] <- 0
}
else{
mu <- runif(p, 2, 3)
}
theSNR <-
(sum(mu ^ 2) - sum((t(V) %*% mu) ^ 2)) / modelParam$theSigmaSq / sqrt(2 * tau ^
2 * p)
# scale the mean
mu <- mu / sqrt(theSNR) * sqrt(SNR)
# the SNR
theNewSNR <-
(sum(mu ^ 2) - sum((t(V) %*% mu) ^ 2)) / modelParam$theSigmaSq / sqrt(2 * tau ^
2 * p)
# permutation function
perFun <- function(){
kSamplePermutation(2,c(n1,n2))
}
iniPer <- kSamplePermutation(2,c(n1,n2),seq(n1+n2))
thePvalue <- rep(0, M)
CQPvalue <- rep(0, M)
for (i in 1:M) {
theData1 <- genSpikedModel(
n = n1,
p = p,
mu = 0,
modelParam = modelParam
)
theData2 <- genSpikedModel(
n = n2,
p = p,
mu = mu,
modelParam = modelParam
)
theData <- rbind(theData1,theData2)
thePvalue[i] <- randomizationTest(
B = B,
GTStat = spikedTempStatistic(n1,n2,r=modelParam$r,theData),
genGT = perFun,
initGT = iniPer,
alpha = 0.05,
needPvalue = needPvalue
)
CQPvalue[i] <- randomizationTest(
B = B,
GTStat = twoSampleTempStatistic(n1,n2,theData),
genGT = perFun,
initGT = iniPer,
alpha = 0.05,
needPvalue = needPvalue
)
}
theResult <- list()
theResult$n1 <- n1
theResult$n2 <- n2
theResult$p <- p
theResult$theoreticalPower <- pnorm(qnorm(0.05) + theNewSNR)
theResult$randomizationPower <- mean(thePvalue <= 0.05)
theResult$CQPower <- mean(CQPvalue <= 0.05)
theResult$SNR <- theNewSNR
unlist(theResult)
}
#' simulation of the level of A_1,...,A_4
#' @param n n
#' @param p p
#' @param M simulate M times test procedure
#' @param k Moving average parameter
#' @param rho Moving average parameter
#' @param mu mean value
#' @param B resample B times to get randomization test
#' @export
simLevelA14 <- function(n,
p,
M,
B,
innov = list("type" = "MA", "innov" = "normal"),
maParam = NULL,
factorParam = NULL,
needPvalue = TRUE) {
A1 <- rep(0, M)
A2 <- rep(0, M)
A3 <- rep(0, M)
A4 <- rep(0, M)
pb <- txtProgressBar(style=3)
for (i in 1:M) {
if (innov$type == "MA") {
theData <- movingAverageGen(n, p, maParam, mu=rep(0,p), innov = innov$innov)
}
else if (innov$type == "factor") {
theData <- fanFactorModelGen(n, p, theParam = factorParam, mu = rep(0,p))
}
S <- t(theData)%*% theData
A1[i] <- randomizationTest(
B = B,
GTStat = A1TempStatistic(theData,S),
genGT = function() {
ranGT = rbinom(n, 1, 0.5) * 2 - 1
return(ranGT)
},
initGT = rep(1, n),
alpha = 0.05,
needPvalue = needPvalue
)
A2[i] <- randomizationTest(
B = B,
GTStat = A2TempStatistic(theData,S),
genGT = function() {
ranGT = rbinom(n, 1, 0.5) * 2 - 1
return(ranGT)
},
initGT = rep(1, n),
alpha = 0.05,
needPvalue = needPvalue
)
A3[i] <- randomizationTest(
B = B,
GTStat = A3TempStatistic(theData,S),
genGT = function() {
ranGT = rbinom(n, 1, 0.5) * 2 - 1
return(ranGT)
},
initGT = rep(1, n),
alpha = 0.05,
needPvalue = needPvalue
)
A4[i] <- randomizationTest(
B = B,
GTStat = A4TempStatistic(theData,S),
genGT = function() {
ranGT = rbinom(n, 1, 0.5) * 2 - 1
return(ranGT)
},
initGT = rep(1, n),
alpha = 0.05,
needPvalue = needPvalue
)
setTxtProgressBar(pb,i/M)
}
close(pb)
theResult <- data.frame(A1,A2,A3,A4)
}
#' simulation of the level of the proposed randomization test
#' @param n n
#' @param p p
#' @param M simulate M times test procedure
#' @param k Moving average parameter
#' @param rho Moving average parameter
#' @param mu mean value
#' @param B resample B times to get randomization test
#' @export
simLevelProposed <- function(n,
p,
M,
B,
innov = list("type" = "MA", "innov" = "normal"),
maParam = NULL,
factorParam = NULL,
compoundParam = NULL,
needPvalue = TRUE) {
thePvalue <- rep(0, M)
CQPvalue <- rep(0, M)
bootstrapPvalue <- rep(0, M)
pb <- txtProgressBar(style=3)
for (i in 1:M) {
if (innov$type == "MA") {
theData <- movingAverageGen(n, p, maParam, mu=rep(0,p), innov = innov$innov)
}
else if (innov$type == "factor") {
theData <- fanFactorModelGen(n, p, theParam = factorParam, mu = rep(0,p))
}
else if (innov$type == "compound"){
theData <- compoundModelGen(n,p, compoundParam, mu= rep(0,p))
}
tmpGTStat <- chenTempStatistic(theData)
thePvalue[i] <- randomizationTest(
B = B,
GTStat = tmpGTStat,
genGT = function() {
ranGT <- rbinom(n, 1, 0.5) * 2 - 1
return(ranGT)
},
theTestStat = tmpGTStat(rep(1,n)),
alpha = 0.05,
needPvalue = needPvalue
)
# Bootstrap
bootstrapPvalue[i] <- randomizationTest(
B = B,
GTStat = bootstrapTempStatistic(theData),
genGT = function() {
ranGT <- sample(n,replace=TRUE)
return(ranGT)
},
theTestStat = tmpGTStat(rep(1,n)),
alpha = 0.05,
needPvalue = needPvalue
)
CQPvalue[i] <- CQTest(theData)
setTxtProgressBar(pb,i/M)
}
close(pb)
theResult <- data.frame(thePvalue,CQPvalue,bootstrapPvalue)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.