pch=17,col='#9b6541')
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="9A" & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='#1a342b')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
pch=17,col='#1a342b')
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="11A" & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='#c0c23b')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
pch=17,col='#c0c23b')
## Plot 1 - baseline scenarios vs theoretical power
par(mfrow=c(1,2))
# theoretical
plot.dat <- res.dat[res.dat$scenario %in% sapply(c('A','B','C','D','E'),function(x) paste0(1:4,x)),]
plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),rep(unique(plot.dat[plot.dat$scenario.type=="A",]$theoretical.power),3),
xlab='N',ylab='Theoretical power',
xaxt='n',yaxt='n',type='l',xlim=c(100,300),ylim=c(0,1),lty=4,col='#03a18a')
axis(1,c(100,200,300))
axis(2,seq(0,1,0.1))
lines(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$theoretical.power),col='#c0c23b',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$theoretical.power),col='#c0c23b',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$theoretical.power),col='#a12471',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$theoretical.power),col='#a12471',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$theoretical.power),col='#9b6541',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$theoretical.power),col='#9b6541',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$theoretical.power),col='#1a342b',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$theoretical.power),col='#1a342b',lty=4)
points(unique(plot.dat[plot.dat$scenario %in% c("1A",'1A'),]$N),
rep(unique(plot.dat[plot.dat$scenario %in% c("1A",'1A'),]$theoretical.power),3),col='#03a18a',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$theoretical.power),col='#c0c23b',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$theoretical.power),col='#c0c23b',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$theoretical.power),col='#a12471',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$theoretical.power),col='#a12471',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$theoretical.power),col='#9b6541',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$theoretical.power),col='#9b6541',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$theoretical.power),col='#1a342b',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$theoretical.power),col='#1a342b',pch=17)
legend("topleft",lty=rep(4,5),col=c("#03a18a","#c0c23b","#a12471","#9b6541","#1a342b"),
pch=rep(17,5),
legend=c("Scenario A",
"Scenario 1-2 / B-D",
"Scenario 3-4 / B-D",
"Scenario 1-2 / C-E",
"Scenario 3-4 / C-E"),cex=0.7)
# real
plot.dat <- res.dat[res.dat$scenario %in% sapply(c('A','B','C','D','E'),function(x) paste0(1:4,x)),]
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario.type=="A" & plot.dat$N==x,]$h0.rejected.p))
plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
xlab='N',ylab='Proportion of null hypothesis rejection',
xaxt='n',yaxt='n',type='l',xlim=c(100,300),ylim=c(0,1),lty=4,col='#03a18a')
axis(1,c(100,200,300))
axis(2,seq(0,1,0.1))
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1B","1D") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
mean.A,col='#c6d18d',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2B","2D") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
mean.A,col='#c0c23b',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3B","3D") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
mean.A,col='#da77c7',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4B","4D") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
mean.A,col='#a12471',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1C","1E") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
mean.A,col='#b5a180',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2C","2E") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
mean.A,col='#9b6541',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3C","3E") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
mean.A,col='#30a466',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4C","4E") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
mean.A,col='#1a342b',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario.type=="A" & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,pch=17,col='#03a18a')
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1B","1D") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
mean.A,col='#c6d18d',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2B","2D") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
mean.A,col='#c0c23b',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3B","3D") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
mean.A,col='#da77c7',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4B","4D") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
mean.A,col='#a12471',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1C","1E") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
mean.A,col='#b5a180',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2C","2E") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
mean.A,col='#9b6541',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3C","3E") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
mean.A,col='#30a466',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4C","4E") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
mean.A,col='#1a342b',pch=17)
legend("topleft",lty=rep(4,5),col=c("#03a18a","#c0c23b","#a12471","#9b6541","#1a342b"),
pch=rep(17,5),
legend=c("Scenario A",
"Scenario 1-2 / B-D",
"Scenario 3-4 / B-D",
"Scenario 1-2 / C-E",
"Scenario 3-4 / C-E"),cex=0.7)
library(tinytex)
tlmgr_install('soul')
res.dat[res.dat$scenario=='2B',]$theoretical.power
res.dat[res.dat$scenario=='2B',]$h0.rejected.p
res.dat[res.dat$scenario=='1B',]$h0.rejected.p
res.dat[res.dat$scenario=='1B',]$theoretical.power
library(TAM)
library(doMC)
library(parallel)
library(pbmcapply)
library(funprog)
library(dplyr)
library(readxl)
res.dat[res.dat$scenario=="2D",]$h0.rejected.p
##############################################################################
#----------------------------------------------------------------------------#
################################# POWER PLOTS ################################
#----------------------------------------------------------------------------#
##############################################################################
## Plot 1 - baseline scenarios vs theoretical power
par(mfrow=c(1,2))
# theoretical
plot.dat <- res.dat[res.dat$scenario %in% sapply(c('A','B','C','D','E'),function(x) paste0(1:4,x)),]
plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),rep(unique(plot.dat[plot.dat$scenario.type=="A",]$theoretical.power),3),
xlab='N',ylab='Theoretical power',
xaxt='n',yaxt='n',type='l',xlim=c(100,300),ylim=c(0,1),lty=4,col='#03a18a')
axis(1,c(100,200,300))
axis(2,seq(0,1,0.1))
lines(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$theoretical.power),col='#c0c23b',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$theoretical.power),col='#c0c23b',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$theoretical.power),col='#a12471',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$theoretical.power),col='#a12471',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$theoretical.power),col='#9b6541',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$theoretical.power),col='#9b6541',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$theoretical.power),col='#1a342b',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$theoretical.power),col='#1a342b',lty=4)
points(unique(plot.dat[plot.dat$scenario %in% c("1A",'1A'),]$N),
rep(unique(plot.dat[plot.dat$scenario %in% c("1A",'1A'),]$theoretical.power),3),col='#03a18a',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$theoretical.power),col='#c0c23b',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$theoretical.power),col='#c0c23b',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$theoretical.power),col='#a12471',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$theoretical.power),col='#a12471',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$theoretical.power),col='#9b6541',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$theoretical.power),col='#9b6541',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$theoretical.power),col='#1a342b',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$theoretical.power),col='#1a342b',pch=17)
legend("topleft",lty=rep(4,5),col=c("#03a18a","#c0c23b","#a12471","#9b6541","#1a342b"),
pch=rep(17,5),
legend=c("Scenario A",
"Scenario 1-2 / B-D",
"Scenario 3-4 / B-D",
"Scenario 1-2 / C-E",
"Scenario 3-4 / C-E"),cex=0.7)
# real
plot.dat <- res.dat[res.dat$scenario %in% sapply(c('A','B','C','D','E'),function(x) paste0(1:4,x)),]
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario.type=="A" & plot.dat$N==x,]$h0.rejected.p))
plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
xlab='N',ylab='Proportion of null hypothesis rejection',
xaxt='n',yaxt='n',type='l',xlim=c(100,300),ylim=c(0,1),lty=4,col='#03a18a')
axis(1,c(100,200,300))
axis(2,seq(0,1,0.1))
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1B","1D") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
mean.A,col='#c6d18d',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2B","2D") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
mean.A,col='#c0c23b',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3B","3D") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
mean.A,col='#da77c7',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4B","4D") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
mean.A,col='#a12471',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1C","1E") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
mean.A,col='#b5a180',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2C","2E") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
mean.A,col='#9b6541',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3C","3E") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
mean.A,col='#30a466',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4C","4E") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
mean.A,col='#1a342b',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario.type=="A" & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,pch=17,col='#03a18a')
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1B","1D") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
mean.A,col='#c6d18d',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2B","2D") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
mean.A,col='#c0c23b',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3B","3D") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
mean.A,col='#da77c7',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4B","4D") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
mean.A,col='#a12471',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1C","1E") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
mean.A,col='#b5a180',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2C","2E") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
mean.A,col='#9b6541',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3C","3E") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
mean.A,col='#30a466',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4C","4E") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
mean.A,col='#1a342b',pch=17)
legend("topleft",lty=rep(4,5),col=c("#03a18a","#c0c23b","#a12471","#9b6541","#1a342b"),
pch=rep(17,5),
legend=c("Scenario A",
"Scenario 1-2 / B-D",
"Scenario 3-4 / B-D",
"Scenario 1-2 / C-E",
"Scenario 3-4 / C-E"),cex=0.7)
###### Puissance théorique
res.dat$theoretical.power <- 0
### Scénarios N=100
## Scénarios J=4 / M=2
res.dat[res.dat$scenario %in% paste0(c(1,5,7,9,11),'A') & res.dat$N==100,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(1,5,7,9,11),'B') & res.dat$N==100,]$theoretical.power <- 0.1543
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'C') & res.dat$N==100,]$theoretical.power <- 0.1543
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'D') & res.dat$N==100,]$theoretical.power <- 0.4627
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'E') & res.dat$N==100,]$theoretical.power <- 0.4627
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'F') & res.dat$N==100,]$theoretical.power <- 0.1543
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'G') & res.dat$N==100,]$theoretical.power <- 0.4627
res.dat[res.dat$scenario %in% paste0(1,'C') & res.dat$N==100,]$theoretical.power <- 0.4627
res.dat[res.dat$scenario %in% paste0(1,'D') & res.dat$N==100,]$theoretical.power <- 0.1543
res.dat[res.dat$scenario %in% paste0(1,'E') & res.dat$N==100,]$theoretical.power <- 0.4627
## Scénarios J=4 / M=4
res.dat[res.dat$scenario %in% paste0(c(2,6,8,10,12),'A') & res.dat$N==100,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(2,6,8,10,12),'B') & res.dat$N==100,]$theoretical.power <- 0.2177
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'C') & res.dat$N==100,]$theoretical.power <- 0.2177
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'D') & res.dat$N==100,]$theoretical.power <- 0.6586
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'E') & res.dat$N==100,]$theoretical.power <- 0.6586
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'F') & res.dat$N==100,]$theoretical.power <- 0.2177
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'G') & res.dat$N==100,]$theoretical.power <- 0.6586
res.dat[res.dat$scenario %in% paste0(2,'C') & res.dat$N==100,]$theoretical.power <- 0.6586
res.dat[res.dat$scenario %in% paste0(2,'D') & res.dat$N==100,]$theoretical.power <- 0.2177
res.dat[res.dat$scenario %in% paste0(2,'E') & res.dat$N==100,]$theoretical.power <- 0.6586
## Scénarios J=7 / M=2
res.dat[res.dat$scenario %in% paste0(c(3,13,15,17,19),'A') & res.dat$N==100,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(3,13,15,17,19),'B') & res.dat$N==100,]$theoretical.power <- 0.1870
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'C') & res.dat$N==100,]$theoretical.power <- 0.1870
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'D') & res.dat$N==100,]$theoretical.power <- 0.5666
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'E') & res.dat$N==100,]$theoretical.power <- 0.5666
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'F') & res.dat$N==100,]$theoretical.power <- 0.1870
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'G') & res.dat$N==100,]$theoretical.power <- 0.5666
res.dat[res.dat$scenario %in% paste0(3,'C') & res.dat$N==100,]$theoretical.power <- 0.5666
res.dat[res.dat$scenario %in% paste0(3,'D') & res.dat$N==100,]$theoretical.power <- 0.1870
res.dat[res.dat$scenario %in% paste0(3,'E') & res.dat$N==100,]$theoretical.power <- 0.5666
## Scénarios J=7 / M=4
res.dat[res.dat$scenario %in% paste0(c(4,14,16,18,20),'A') & res.dat$N==100,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(4,14,16,18,20),'B') & res.dat$N==100,]$theoretical.power <- 0.2450
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'C') & res.dat$N==100,]$theoretical.power <- 0.2450
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'D') & res.dat$N==100,]$theoretical.power <- 0.7136
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'E') & res.dat$N==100,]$theoretical.power <- 0.7136
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'F') & res.dat$N==100,]$theoretical.power <- 0.2450
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'G') & res.dat$N==100,]$theoretical.power <- 0.7136
res.dat[res.dat$scenario %in% paste0(4,'C') & res.dat$N==100,]$theoretical.power <- 0.7136
res.dat[res.dat$scenario %in% paste0(4,'D') & res.dat$N==100,]$theoretical.power <- 0.2450
res.dat[res.dat$scenario %in% paste0(4,'E') & res.dat$N==100,]$theoretical.power <- 0.7136
### Scénarios N=200
## Scénarios J=4 / M=2
res.dat[res.dat$scenario %in% paste0(c(1,5,7,9,11),'A') & res.dat$N==200,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(1,5,7,9,11),'B') & res.dat$N==200,]$theoretical.power <- 0.2618
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'C') & res.dat$N==200,]$theoretical.power <- 0.2618
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'D') & res.dat$N==200,]$theoretical.power <- 0.7507
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'E') & res.dat$N==200,]$theoretical.power <- 0.7507
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'F') & res.dat$N==200,]$theoretical.power <- 0.2618
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'G') & res.dat$N==200,]$theoretical.power <- 0.7507
res.dat[res.dat$scenario %in% paste0(1,'C') & res.dat$N==200,]$theoretical.power <- 0.7507
res.dat[res.dat$scenario %in% paste0(1,'D') & res.dat$N==200,]$theoretical.power <- 0.2618
res.dat[res.dat$scenario %in% paste0(1,'E') & res.dat$N==200,]$theoretical.power <- 0.7507
## Scénarios J=4 / M=4
res.dat[res.dat$scenario %in% paste0(c(2,6,8,10,12),'A') & res.dat$N==200,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(2,6,8,10,12),'B') & res.dat$N==200,]$theoretical.power <- 0.3875
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'C') & res.dat$N==200,]$theoretical.power <- 0.3875
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'D') & res.dat$N==200,]$theoretical.power <- 0.9161
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'E') & res.dat$N==200,]$theoretical.power <- 0.9161
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'F') & res.dat$N==200,]$theoretical.power <- 0.3875
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'G') & res.dat$N==200,]$theoretical.power <- 0.9161
res.dat[res.dat$scenario %in% paste0(2,'C') & res.dat$N==200,]$theoretical.power <- 0.9161
res.dat[res.dat$scenario %in% paste0(2,'D') & res.dat$N==200,]$theoretical.power <- 0.3875
res.dat[res.dat$scenario %in% paste0(2,'E') & res.dat$N==200,]$theoretical.power <- 0.9161
## Scénarios J=7 / M=2
res.dat[res.dat$scenario %in% paste0(c(3,13,15,17,19),'A') & res.dat$N==200,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(3,13,15,17,19),'B') & res.dat$N==200,]$theoretical.power <- 0.3258
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'C') & res.dat$N==200,]$theoretical.power <- 0.3258
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'D') & res.dat$N==200,]$theoretical.power <- 0.8538
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'E') & res.dat$N==200,]$theoretical.power <- 0.8538
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'F') & res.dat$N==200,]$theoretical.power <- 0.3258
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'G') & res.dat$N==200,]$theoretical.power <- 0.8538
res.dat[res.dat$scenario %in% paste0(3,'C') & res.dat$N==200,]$theoretical.power <- 0.8538
res.dat[res.dat$scenario %in% paste0(3,'D') & res.dat$N==200,]$theoretical.power <- 0.3258
res.dat[res.dat$scenario %in% paste0(3,'E') & res.dat$N==200,]$theoretical.power <- 0.8538
## Scénarios J=7 / M=4
res.dat[res.dat$scenario %in% paste0(c(4,14,16,18,20),'A') & res.dat$N==200,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(4,14,16,18,20),'B') & res.dat$N==200,]$theoretical.power <- 0.4321
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'C') & res.dat$N==200,]$theoretical.power <- 0.4321
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'D') & res.dat$N==200,]$theoretical.power <- 0.9471
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'E') & res.dat$N==200,]$theoretical.power <- 0.9471
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'F') & res.dat$N==200,]$theoretical.power <- 0.4321
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'G') & res.dat$N==200,]$theoretical.power <- 0.9471
res.dat[res.dat$scenario %in% paste0(4,'C') & res.dat$N==200,]$theoretical.power <- 0.9471
res.dat[res.dat$scenario %in% paste0(4,'D') & res.dat$N==200,]$theoretical.power <- 0.4321
res.dat[res.dat$scenario %in% paste0(4,'E') & res.dat$N==200,]$theoretical.power <- 0.9471
### Scénarios N=300
## Scénarios J=4 / M=2
res.dat[res.dat$scenario %in% paste0(c(1,5,7,9,11),'A') & res.dat$N==300,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(1,5,7,9,11),'B') & res.dat$N==300,]$theoretical.power <- 0.3660
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'C') & res.dat$N==300,]$theoretical.power <- 0.3660
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'D') & res.dat$N==300,]$theoretical.power <- 0.8981
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'E') & res.dat$N==300,]$theoretical.power <- 0.8981
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'F') & res.dat$N==300,]$theoretical.power <- 0.3660
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'G') & res.dat$N==300,]$theoretical.power <- 0.8981
res.dat[res.dat$scenario %in% paste0(1,'C') & res.dat$N==300,]$theoretical.power <- 0.8981
res.dat[res.dat$scenario %in% paste0(1,'D') & res.dat$N==300,]$theoretical.power <- 0.3660
res.dat[res.dat$scenario %in% paste0(1,'E') & res.dat$N==300,]$theoretical.power <- 0.8981
## Scénarios J=4 / M=4
res.dat[res.dat$scenario %in% paste0(c(2,6,8,10,12),'A') & res.dat$N==300,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(2,6,8,10,12),'B') & res.dat$N==300,]$theoretical.power <- 0.5373
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'C') & res.dat$N==300,]$theoretical.power <- 0.5373
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'D') & res.dat$N==300,]$theoretical.power <- 0.9834
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'E') & res.dat$N==300,]$theoretical.power <- 0.9834
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'F') & res.dat$N==300,]$theoretical.power <- 0.5373
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'G') & res.dat$N==300,]$theoretical.power <- 0.9834
res.dat[res.dat$scenario %in% paste0(2,'C') & res.dat$N==300,]$theoretical.power <- 0.9834
res.dat[res.dat$scenario %in% paste0(2,'D') & res.dat$N==300,]$theoretical.power <- 0.5373
res.dat[res.dat$scenario %in% paste0(2,'E') & res.dat$N==300,]$theoretical.power <- 0.9834
## Scénarios J=7 / M=2
res.dat[res.dat$scenario %in% paste0(c(3,13,15,17,19),'A') & res.dat$N==300,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(3,13,15,17,19),'B') & res.dat$N==300,]$theoretical.power <- 0.4550
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'C') & res.dat$N==300,]$theoretical.power <- 0.4550
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'D') & res.dat$N==300,]$theoretical.power <- 0.9584
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'E') & res.dat$N==300,]$theoretical.power <- 0.9584
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'F') & res.dat$N==300,]$theoretical.power <- 0.4550
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'G') & res.dat$N==300,]$theoretical.power <- 0.9584
res.dat[res.dat$scenario %in% paste0(3,'C') & res.dat$N==300,]$theoretical.power <- 0.9584
res.dat[res.dat$scenario %in% paste0(3,'D') & res.dat$N==300,]$theoretical.power <- 0.4550
res.dat[res.dat$scenario %in% paste0(3,'E') & res.dat$N==300,]$theoretical.power <- 0.9584
## Scénarios J=7 / M=4
res.dat[res.dat$scenario %in% paste0(c(4,14,16,18,20),'A') & res.dat$N==300,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(4,14,16,18,20),'B') & res.dat$N==300,]$theoretical.power <- 0.5907
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'C') & res.dat$N==300,]$theoretical.power <- 0.5907
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'D') & res.dat$N==300,]$theoretical.power <- 0.9919
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'E') & res.dat$N==300,]$theoretical.power <- 0.9919
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'F') & res.dat$N==300,]$theoretical.power <- 0.5907
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'G') & res.dat$N==300,]$theoretical.power <- 0.9919
res.dat[res.dat$scenario %in% paste0(4,'C') & res.dat$N==300,]$theoretical.power <- 0.9919
res.dat[res.dat$scenario %in% paste0(4,'D') & res.dat$N==300,]$theoretical.power <- 0.5907
res.dat[res.dat$scenario %in% paste0(4,'E') & res.dat$N==300,]$theoretical.power <- 0.9919
### DIF scenarios
res.dat.dif$theoretical.power <- res.dat[61:nrow(res.dat),]$theoretical.power
## Plot 1 - baseline scenarios vs theoretical power
par(mfrow=c(1,2))
# theoretical
plot.dat <- res.dat[res.dat$scenario %in% sapply(c('A','B','C','D','E'),function(x) paste0(1:4,x)),]
plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),rep(unique(plot.dat[plot.dat$scenario.type=="A",]$theoretical.power),3),
xlab='N',ylab='Theoretical power',
xaxt='n',yaxt='n',type='l',xlim=c(100,300),ylim=c(0,1),lty=4,col='#03a18a')
axis(1,c(100,200,300))
axis(2,seq(0,1,0.1))
lines(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$theoretical.power),col='#c6d18d',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$theoretical.power),col='#c0c23b',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$theoretical.power),col='#da77c7',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$theoretical.power),col='#a12471',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$theoretical.power),col='#b5a180',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$theoretical.power),col='#9b6541',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$theoretical.power),col='#30a466',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$theoretical.power),col='#1a342b',lty=4)
points(unique(plot.dat[plot.dat$scenario %in% c("1A",'1A'),]$N),
rep(unique(plot.dat[plot.dat$scenario %in% c("1A",'1A'),]$theoretical.power),3),col='#03a18a',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$theoretical.power),col='#c6d18d',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$theoretical.power),col='#c0c23b',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$theoretical.power),col='#da77c7',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$theoretical.power),col='#a12471',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$theoretical.power),col='#b5a180',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$theoretical.power),col='#9b6541',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$theoretical.power),col='#30a466',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$theoretical.power),col='#1a342b',pch=17)
legend("topleft",lty=rep(4,5),col=c("#03a18a","#c0c23b","#a12471","#9b6541","#1a342b"),
pch=rep(17,5),
legend=c("Scenario A",
"Scenario 1-2 / B-D",
"Scenario 3-4 / B-D",
"Scenario 1-2 / C-E",
"Scenario 3-4 / C-E"),cex=0.7)
# real
plot.dat <- res.dat[res.dat$scenario %in% sapply(c('A','B','C','D','E'),function(x) paste0(1:4,x)),]
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario.type=="A" & plot.dat$N==x,]$h0.rejected.p))
plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
xlab='N',ylab='Proportion of null hypothesis rejection',
xaxt='n',yaxt='n',type='l',xlim=c(100,300),ylim=c(0,1),lty=4,col='#03a18a')
axis(1,c(100,200,300))
axis(2,seq(0,1,0.1))
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1B","1D") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
mean.A,col='#c6d18d',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2B","2D") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
mean.A,col='#c0c23b',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3B","3D") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
mean.A,col='#da77c7',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4B","4D") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
mean.A,col='#a12471',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1C","1E") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
mean.A,col='#b5a180',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2C","2E") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
mean.A,col='#9b6541',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3C","3E") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
mean.A,col='#30a466',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4C","4E") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
mean.A,col='#1a342b',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario.type=="A" & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,pch=17,col='#03a18a')
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1B","1D") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
mean.A,col='#c6d18d',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2B","2D") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
mean.A,col='#c0c23b',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3B","3D") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
mean.A,col='#da77c7',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4B","4D") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
mean.A,col='#a12471',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1C","1E") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
mean.A,col='#b5a180',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2C","2E") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
mean.A,col='#9b6541',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3C","3E") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
mean.A,col='#30a466',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4C","4E") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
mean.A,col='#1a342b',pch=17)
legend("topleft",lty=rep(4,5),col=c("#03a18a","#c0c23b","#a12471","#9b6541","#1a342b"),
pch=rep(17,5),
legend=c("Scenario A",
"Scenario 1-2 / B-D",
"Scenario 3-4 / B-D",
"Scenario 1-2 / C-E",
"Scenario 3-4 / C-E"),cex=0.7)