##############################################################################
#----------------------------------------------------------------------------#
################################# 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='#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)



##############################################################################
#----------------------------------------------------------------------------#
################################# ALPHA PLOTS ################################
#----------------------------------------------------------------------------#
##############################################################################

par(mfrow=c(1,1))

#### Scenarios with J=4 / M=2

# A / H=0

plot.dat <- res.dat

mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="1A" & 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',main='Scenarios where H0 is TRUE')
axis(1,c(100,200,300))
axis(2,seq(0,1,0.1))
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=="5A" & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='#a12471')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
       pch=17,col='#a12471')


mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="7A" & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='#9b6541')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
       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')






##############################################################################
#----------------------------------------------------------------------------#
################################# BIAS PLOTS #################################
#----------------------------------------------------------------------------#
##############################################################################

par(mfrow=c(1,1))

#### Scenarios with J=4 / M=2

# A / H=0

plot.dat <- res.dat

mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="1A" & plot.dat$N==x,]$bias))
plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
     xlab='N',ylab='Bias',
     xaxt='n',yaxt='n',type='l',xlim=c(100,300),ylim=c(-1,1),lty=4,col='#03a18a',main='Scenarios where H0 is TRUE')
axis(1,c(100,200,300))
axis(2,seq(0,1,0.1))
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=="5A" & plot.dat$N==x,]$bias))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='#a12471')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
       pch=17,col='#a12471')


mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="7A" & plot.dat$N==x,]$bias))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='#9b6541')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
       pch=17,col='#9b6541')

mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="9A" & plot.dat$N==x,]$bias))
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,]$bias))
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')


par(mfrow=c(1,1))

#### Scenarios with J=4 / M=2

# B/ H=0

plot.dat <- res.dat

mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="1B" & plot.dat$N==x,]$bias))
plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
     xlab='N',ylab='Bias',
     xaxt='n',yaxt='n',type='l',xlim=c(100,300),ylim=c(-1,1),lty=4,col='#03a18a',main='Scenarios where H0 is TRUE')
axis(1,c(100,200,300))
axis(2,seq(0,1,0.1))
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=="5B" & plot.dat$N==x,]$bias))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='#a12471')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
       pch=17,col='#a12471')


mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="7B" & plot.dat$N==x,]$bias))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='#9b6541')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
       pch=17,col='#9b6541')

mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="9B" & plot.dat$N==x,]$bias))
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=="11B" & plot.dat$N==x,]$bias))
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')


mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="5C" & plot.dat$N==x,]$bias))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='red')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
       pch=17,col='red')


mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="7C" & plot.dat$N==x,]$bias))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='blue')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
       pch=17,col='blue')

mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="9C" & plot.dat$N==x,]$bias))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='purple')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
       pch=17,col='purple')


mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="11C" & plot.dat$N==x,]$bias))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='orange')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
       pch=17,col='orange')