You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
513 lines
35 KiB
R
513 lines
35 KiB
R
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
|
|
points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='#053305',pch=3)
|
|
# J=7 / 2 DIF
|
|
res.dat.temp <- res.dat[res.dat$J==7 & res.dat$nb.dif==2,]
|
|
boxplot(true.value.in.ci.p~dif.size,data=res.dat.temp,col=c(2,3,3,2),xlab='DIF size',
|
|
ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
|
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.5,]
|
|
points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
|
|
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.5,]
|
|
points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
|
|
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.3,]
|
|
points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
|
|
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
|
|
points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='#053305',pch=3)
|
|
# J=7 / 3 DIF
|
|
res.dat.temp <- res.dat[res.dat$J==7 & res.dat$nb.dif==3,]
|
|
boxplot(true.value.in.ci.p~dif.size,data=res.dat.temp,col=c(2,3,3,2),xlab='DIF size',
|
|
ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
|
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.5,]
|
|
points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
|
|
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.5,]
|
|
points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
|
|
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.3,]
|
|
points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
|
|
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
|
|
points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='#053305',pch=3)
|
|
##############################################################################
|
|
#----------------------------------------------------------------------------#
|
|
########################## BETA SIGN CHANGE BOXPLOTS #########################
|
|
#----------------------------------------------------------------------------#
|
|
##############################################################################
|
|
# Overall
|
|
boxplot(beta.same.sign.truebeta.p~dif.size,data=res.dat,col=c(2,3),xlab='DIF size',
|
|
ylab='Proportion of estimates with same sign as true value in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
|
res.null3 <- res.dat[res.dat$dif.size==-0.5,]
|
|
points(y=res.null3$beta.same.sign.truebeta.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
|
|
res.null3 <- res.dat[res.dat$dif.size==0.5,]
|
|
points(y=res.null3$beta.same.sign.truebeta.p,x=rep(5,nrow(res.null3)),col='#590b0c',pch=3)
|
|
res.null3 <- res.dat[res.dat$dif.size==-0.3,]
|
|
points(y=res.null3$beta.same.sign.truebeta.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
|
|
res.null3 <- res.dat[res.dat$dif.size==0.3,]
|
|
points(y=res.null3$beta.same.sign.truebeta.p,x=rep(4,nrow(res.null3)),col='#053305',pch=3)
|
|
res.null3 <- res.dat[res.dat$dif.size==0,]
|
|
points(y=res.null3$beta.same.sign.truebeta.p,x=rep(3,nrow(res.null3)),col='gray',pch=3)
|
|
# Overall // H0 rejected
|
|
boxplot(beta.same.sign.truebeta.signif.p~dif.size,data=res.dat,col=c(2,3),xlab='DIF size',
|
|
ylab='Proportion of estimates with same sign as true value in target scenario',main='When H0 rejected',ylim=c(0,1))
|
|
res.null3 <- res.dat[res.dat$dif.size==-0.5,]
|
|
points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
|
|
res.null3 <- res.dat[res.dat$dif.size==0.5,]
|
|
points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(5,nrow(res.null3)),col='#590b0c',pch=3)
|
|
res.null3 <- res.dat[res.dat$dif.size==-0.3,]
|
|
points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
|
|
res.null3 <- res.dat[res.dat$dif.size==0.3,]
|
|
points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(4,nrow(res.null3)),col='#053305',pch=3)
|
|
res.null3 <- res.dat[res.dat$dif.size==0,]
|
|
points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(3,nrow(res.null3)),col='gray',pch=3)
|
|
## 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')
|
|
unique(plot.dat[plot.dat$scenario.type=="A",]$N)
|
|
View(res.dat)
|
|
##############################################################################
|
|
#----------------------------------------------------------------------------#
|
|
################################# 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),4),
|
|
xlab='N',ylab='Theoretical power',
|
|
xaxt='n',yaxt='n',type='l',xlim=c(50,300),ylim=c(0,1),lty=4,col='#03a18a')
|
|
axis(1,c(50,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)
|
|
##############################################################################
|
|
#----------------------------------------------------------------------------#
|
|
################################# 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),4),
|
|
xlab='N',ylab='Theoretical power',
|
|
xaxt='n',yaxt='n',type='l',xlim=c(50,300),ylim=c(0,1),lty=4,col='#03a18a')
|
|
axis(1,c(50,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)
|
|
## 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),4),
|
|
xlab='N',ylab='Theoretical power',
|
|
xaxt='n',yaxt='n',type='l',xlim=c(50,300),ylim=c(0,1),lty=4,col='#03a18a')
|
|
axis(1,c(50,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)
|
|
##############################################################################
|
|
#----------------------------------------------------------------------------#
|
|
################################# 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),4),
|
|
xlab='N',ylab='Theoretical power',
|
|
xaxt='n',yaxt='n',type='l',xlim=c(50,300),ylim=c(0,1),lty=4,col='#03a18a')
|
|
axis(1,c(50,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)
|
|
## 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),4),
|
|
xlab='N',ylab='Theoretical power',
|
|
xaxt='n',yaxt='n',type='l',xlim=c(50,300),ylim=c(0,1),lty=4,col='#03a18a')
|
|
axis(1,c(50,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),4),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(50,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(50,300),ylim=c(0,1),lty=4,col='#03a18a')
|
|
axis(1,c(50,100,200,300))
|
|
axis(2,seq(0,1,0.1))
|
|
mean.A <- sapply(c(50,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(50,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(50,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(50,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(50,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(50,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(50,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(50,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(50,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(50,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(50,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(50,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(50,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(50,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(50,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(50,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(50,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)
|
|
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(50,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(50,300),ylim=c(0,1),lty=4,col='#03a18a')
|
|
axis(1,c(50,100,200,300))
|
|
axis(2,seq(0,1,0.1))
|
|
unique(plot.dat[plot.dat$scenario.type=="A",]$N)
|
|
mean.A
|
|
mean.A[c(4,1:3)]
|
|
# 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(50,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[c(4,1:3)],
|
|
xlab='N',ylab='Proportion of null hypothesis rejection',
|
|
xaxt='n',yaxt='n',type='l',xlim=c(50,300),ylim=c(0,1),lty=4,col='#03a18a')
|
|
axis(1,c(50,100,200,300))
|
|
axis(2,seq(0,1,0.1))
|
|
mean.A <- sapply(c(50,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(50,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(50,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(50,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(50,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(50,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(50,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(50,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(50,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[c(4,1:3)],pch=17,col='#03a18a')
|
|
mean.A <- sapply(c(50,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(50,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(50,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(50,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(50,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(50,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(50,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(50,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)
|
|
# 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(50,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[c(4,1:3)],
|
|
xlab='N',ylab='Proportion of null hypothesis rejection',
|
|
xaxt='n',yaxt='n',type='l',xlim=c(50,300),ylim=c(0,1),lty=4,col='#03a18a')
|
|
axis(1,c(50,100,200,300))
|
|
axis(2,seq(0,1,0.1))
|
|
unique(plot.dat[plot.dat$scenario.type=="A",]$N)
|
|
mean.A[c(4,1:3)]
|
|
mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1B","1D") & plot.dat$N==x,]$h0.rejected.p))
|
|
unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N)
|
|
##############################################################################
|
|
#----------------------------------------------------------------------------#
|
|
################################# 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),4),
|
|
xlab='N',ylab='Theoretical power',
|
|
xaxt='n',yaxt='n',type='l',xlim=c(50,300),ylim=c(0,1),lty=4,col='#03a18a')
|
|
axis(1,c(50,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),4),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(50,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[c(4,1:3)],
|
|
xlab='N',ylab='Proportion of null hypothesis rejection',
|
|
xaxt='n',yaxt='n',type='l',xlim=c(50,300),ylim=c(0,1),lty=4,col='#03a18a')
|
|
axis(1,c(50,100,200,300))
|
|
axis(2,seq(0,1,0.1))
|
|
mean.A <- sapply(c(50,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[c(4,1:3)],col='#c6d18d',lty=4)
|
|
mean.A <- sapply(c(50,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[c(4,1:3)],col='#c0c23b',lty=4)
|
|
mean.A <- sapply(c(50,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[c(4,1:3)],col='#da77c7',lty=4)
|
|
mean.A <- sapply(c(50,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[c(4,1:3)],col='#a12471',lty=4)
|
|
mean.A <- sapply(c(50,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[c(4,1:3)],col='#b5a180',lty=4)
|
|
mean.A <- sapply(c(50,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[c(4,1:3)],col='#9b6541',lty=4)
|
|
mean.A <- sapply(c(50,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[c(4,1:3)],col='#30a466',lty=4)
|
|
mean.A <- sapply(c(50,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[c(4,1:3)],col='#1a342b',lty=4)
|
|
mean.A <- sapply(c(50,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[c(4,1:3)],pch=17,col='#03a18a')
|
|
mean.A <- sapply(c(50,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[c(4,1:3)],col='#c6d18d',pch=17)
|
|
mean.A <- sapply(c(50,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[c(4,1:3)],col='#c0c23b',pch=17)
|
|
mean.A <- sapply(c(50,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[c(4,1:3)],col='#da77c7',pch=17)
|
|
mean.A <- sapply(c(50,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[c(4,1:3)],col='#a12471',pch=17)
|
|
mean.A <- sapply(c(50,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[c(4,1:3)],col='#b5a180',pch=17)
|
|
mean.A <- sapply(c(50,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[c(4,1:3)],col='#9b6541',pch=17)
|
|
mean.A <- sapply(c(50,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[c(4,1:3)],col='#30a466',pch=17)
|
|
mean.A <- sapply(c(50,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[c(4,1:3)],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)
|
|
View(res.dat.dif[res.dat.dif$scenario=="20A",])
|