library("future")
library("listenv")
##computeNodes <- c("cauchy", "leibniz", "bolzano", "shannon", "euler", "hamming", "bernoulli")
plan(remote, workers = rep("bernoulli", 2)) ## retrieve results (remote!)
plan(eager) ## retrieve results (local!)
test %<-% getwd()
dat %<-% {
rpath <- file.path("~/Documents/Packages/sanssouci", path)
fls <- list.files(rpath, full.names=TRUE)
id <- gsub(".rds$", "", basename(fls))
names(fls) <- id
plyr::ldply(fls, readRDS, .id="id")
}
head(dat)
kc <- as.character(dat$kMax)
kc[which(dat$kMax==m)] <- "m"
kc[which(dat$kMax==m/2)] <- "m/2"
kc[which(dat$kMax==(1-dat$pi0)*m)] <- "(1-pi0)*m"
dat$kMaxC <- kc
##dat$estPowM1 <- dat$estPow/(1-dat$pi0)
powerz <- list("Power1"="P(S(R,H1)>1", "Power"="P(S(R,H)>1", "s1"="E(S(R,H1))/m1", "s1"="E(S(R,H1))/m1")
## res <- rbind(JR=rej0>0, detPow1=rej1>0, detPow=rej01>0, v0, estPow1=s1, estPow=s01)
powerz <- list("detPow1"="P(S(R,H1)>1", "detPow"="P(S(R,H)>1",
"estPow1"="E(S(R,H1))/m1", "estPow"="E(S(R,H))/m1",
"powBH"="half of BH(0.05)", "pow0"="half of {p <= 0.05}")
## some reshaping
#datC <- subset(dat, kMax %in% as.character(c(10, 100, m, m/2)) & flavor != "unadjusted")
datC <- subset(dat, flavor != "unadjusted")
mm <- grep("Oracle|alpha", datC$flavor)
datC <- datC[-mm, ]
levels(datC$family) <- list("Balanced"="kFWER", "Linear"="Simes")
datC$alpha <- as.numeric(datC$alpha)
alphas <- unique(datC$alpha)
datC$ff <- factor(paste(datC$family, datC$flavor))
datC$ff <- factor(paste(datC$family, datC$kMaxC))
if (FALSE) {
## drop outliers?
ww <- which(with(datC, dep==0 & pi0==0.8 & alpha==0.025 & Power < 0.5))
datC <- datC[-ww, ]
}
figName <- sname0
figPath <- file.path("fig", sprintf("BNR,%s", ptag))
figPath <- R.utils::Arguments$getWritablePath(figPath)
pname2 <- gsub("0\\.", "", pname) ## to avoid '.' in LaTeX file names
## plot Power vs JFWER
library("ggplot2")
x <- c("JR", "alpha")[2]
SNRs <- unique(dat$SNR)
#SNRs <- 2
confs <- expand.grid(SNR=SNRs, power=names(powerz), x=x, stringsAsFactors=FALSE)
for (ii in 1:nrow(confs)) {
snr <- confs[ii, "SNR"]
pp <- confs[ii, "power"]
xx <- confs[ii, "x"]
pow <- powerz[[pp]]
ftag <- sprintf("SNR=%s,power=%s", snr, gsub("/", "--", pow))
if (xx=="JR") {
filename <- sprintf("PowerVsJFWER,%s,%s,%s.pdf", figName, pname2, ftag)
} else if (xx=="alpha") {
filename <- sprintf("PowerVsNominalJFWER,%s,%s,%s.pdf", figName, pname2, ftag)
} else {
stop("I don't know what to do when x=", xx)
}
pathname <- file.path(figPath, filename)
datI <- subset(datC, SNR==snr)
pdf(pathname)
p <- ggplot(datI, aes_string(x=xx, y=pp, group="ff", color="ff"))
##p <- p + geom_line() + geom_point(aes(shape=factor(alpha))) + scale_shape_identity()
#p <- p + geom_line(aes(linetype=flavor))
p <- p + geom_line()
##p <- p + facet_grid(pi0 ~ rho, scales="free_y")
p <- p + facet_grid(pi0 ~ dep,
scales="free_y",
labeller=label_bquote(
cols= rho==.(dep),
rows= pi[0]==.(pi0)))
vdat <- data.frame(x=unique(datC$alpha))
#p <- p + geom_vline(aes(xintercept=x), data=vdat, linetype="dashed", color="gray")
p <- p + scale_x_continuous(breaks=round(alphas, 2), minor_breaks=NULL, limits = range(alphas))
p <- p + scale_y_continuous(minor_breaks=NULL)
p <- p + theme(axis.text.x=element_text(angle=90))
p <- p + labs(color="Family",
linetype=expression(lambda-adjustment))
if (xx=="JR") { ## not useful otherwise
p <- p + geom_point(aes(shape=factor(alpha))) + scale_shape_manual(values=1:9)
p <- p + labs(shape="Target JR level")
} else {
p <- p + geom_point()
}
p <- p + labs(y=pow)
p <- p + scale_color_brewer(type="div")
print(p)
dev.off()
}
## sanity check wrt power definitions
pow1 <- subset(dat, family=="kFWER" & flavor =="Single Step")
pow1sd <- subset(dat, family=="kFWER" & flavor =="Step down")
identical(pow1$Power, pow1sd$Power)
identical(pow1$Power1, pow1sd$Power1)
table(sign(pow1sd$Power1-pow1$Power1))
identical(pow1$s1, pow1sd$s1)
table(sign(pow1sd$s1-pow1$s1))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.