# this function generates pValues
myGen <- function(n, n0) {
list(procInput=list(pValues = c(runif(n-n0, 0, 0.01), runif(n0))),
groundTruth = c(rep(FALSE, times=n-n0), rep(TRUE, times=n0)))
}
sim <- simulation(replications = 10, list(funName="myGen", fun=myGen, n=200, n0=c(50,100)),
list(list(funName="BH",
fun=function(pValues, alpha) BH(pValues, alpha, silent=TRUE), alpha=c(0.25, 0.5)),
list(funName="holm", fun=holm, alpha=c(0.25, 0.5),silent=TRUE)))
# the following happend:
# Call myGen(200,50) and append result to sim$data
# Apply bonferroni and holm each with alpha 0.25 and 0.5 to this data set
# Append the results to sim$restults
# Repeat this 10 times.
# Call myGen(200, 100) and append restult to sim$data
# Apply bonferroni and holm each with alpha 0.25 and 0.5 to this data set
# Append the results to sim$restults
# Repeat this 10 times.
length(sim$data)
length(sim$results)
# we now reproduce the 6th item in results
print(sim$results[[6]]$data.set.number)
print(sim$results[[6]]$parameters)
all(BH(sim$data[[2]]$procInput$pValues, 0.5, silent=TRUE)$adjPValues == sim$results[[6]]$adjPValues)
#
# Just calculating some statistics and making some plots
NumberOfType1Error <- function(data, result) sum(data$groundTruth * result$rejected)
result.all <- gatherStatistics(sim, list(NumOfType1Err = NumberOfType1Error))
result <- gatherStatistics(sim, list(NumOfType1Err = NumberOfType1Error),
list(median=median, mean=mean, sd=sd))
print(result)
require(lattice)
histogram(~NumOfType1Err | method*alpha, data = result.all$statisticDF)
barchart(NumOfType1Err.median + NumOfType1Err.mean ~ method | alpha, data = result$statisticDF)
Run the code above in your browser using DataLab