diff --git a/RProject/.Rhistory b/RProject/.Rhistory
index b38fa25..76b9b2d 100644
--- a/RProject/.Rhistory
+++ b/RProject/.Rhistory
@@ -1,512 +1,512 @@
-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)
-library(TAM)
-library(doMC)
-library(parallel)
-library(pbmcapply)
-library(funprog)
-library(dplyr)
-library(readxl)
-aaaa <- read_excel("/home/corentin/Documents/These/Recherche/Simulations/Analysis/DIF/N50/out/8A_50.xls")
-aaaa
-mean(aaaa$beta)
-mean(aaaa$dif1_1)
-mean(aaaa$dif1_2)
-mean(aaaa$dif1_3)
-mean(aaaa$dif1_3,na.rm = T)
-mean(aaaa$item1_1,na.rm = T)
-aaaa <- read_excel("/home/corentin/Documents/These/Recherche/Simulations/Analysis/DIF/N50/out/10A_50.xls")
-mean(aaaa$beta)
-mean(aaaa$dif1_3,na.rm = T)
-mean(aaaa$dif2_3,na.rm = T)
-mean(aaaa$dif1_3,na.rm = T)
-mean(aaaa$dif1_1,na.rm = T)
-aaaa
-aaaa <- read_excel("/home/corentin/Documents/These/Recherche/Simulations/Analysis/DIF/N50/out/10A_50.xls")
-library(TAM)
-library(doMC)
-library(parallel)
-library(pbmcapply)
-library(funprog)
-library(dplyr)
-library(readxl)
-aaaa <- read_excel("/home/corentin/Documents/These/Recherche/Simulations/Analysis/DIF/N50/out/10A_50.xls")
-mean(aaaa$item1_3,na.rm = T)
-mean(aaaa$item2_3,na.rm = T)
-mean(aaaa$item3_3,na.rm = T)
-mean(aaaa$item4_3,na.rm = T)
-mean(aaaa$item4_2,na.rm = T)
-mean(aaaa$item1_2,na.rm = T)
-mean(aaaa$item1_1,na.rm = T)
-mean(aaaa$item1_2,na.rm = T)-mean(aaaa$item1_1,na.rm = T)
-mean(aaaa$item1_3,na.rm = T)-mean(aaaa$item1_2,na.rm = T)
-mean(aaaa$item4_3,na.rm = T)-mean(aaaa$item4_2,na.rm = T)
+m.item5=mean(s$item5),m.item6=mean(s$item6),m.item7=mean(s$item7))
+} else {
+a <- data.frame(m.item1_1=mean(s$item1_1),m.item1_2=mean(s$item1_2),m.item1_3=mean(s$item1_3),
+m.item2_1=mean(s$item2_1),m.item2_2=mean(s$item2_2),m.item2_3=mean(s$item2_3),
+m.item3_1=mean(s$item3_1),m.item3_2=mean(s$item3_2),m.item3_3=mean(s$item3_3),
+m.item4_1=mean(s$item4_1),m.item4_2=mean(s$item4_2),m.item4_3=mean(s$item4_3),
+m.item5_1=mean(s$item5_1),m.item5_2=mean(s$item5_2),m.item5_3=mean(s$item5_3),
+m.item6_1=mean(s$item6_1),m.item6_2=mean(s$item6_2),m.item6_3=mean(s$item6_3),
+m.item7_1=mean(s$item7_1),m.item7_2=mean(s$item7_2),m.item7_3=mean(s$item7_3)
+)
+}
+}
+zz <- substr(scenario,start=0,stop=nchar(scenario)-4)
+b <- data.frame(scenario=zz,
+scenario.type=substr(zz,start=nchar(zz),stop=nchar(zz)),
+N=substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario)),
+J=unique(s$J),
+M=unique(s$M),
+eff.size=unique(s$eff.size),
+nb.dif=unique(s$nb.dif),
+dif.size=unique(s$dif.size)
+)
+z <- data.frame(m.beta=mean(s$beta),
+se.empirical.beta=sd(s$beta),
+se.analytical.beta=mean(s$se.beta),
+m.low.ci.beta=mean(s$low.ci.beta),
+m.high.ci.beta=mean(s$high.ci.beta),
+true.value.in.ci.p=mean(s$true.value.in.ci),
+h0.rejected.p=mean(s$h0.rejected),
+beta.same.sign.truebeta.p=mean(s$beta.same.sign.truebeta,na.rm=T),
+beta.same.sign.truebeta.signif.p=mean(s[s$h0.rejected==1,]$beta.same.sign.truebeta,na.rm=T))
+d <- cbind(b,a,z)
+d$prop.
+return(d)
+}
+#### Compiled results
+res.dat <- compile_simulation('1A_100')
+for (x in results[seq(2,length(results))]) {
+y <- compile_simulation(x)
+res.dat <- bind_rows(res.dat,y)
+}
+#### Create data.frame
+results <- c(sapply(1:4,function(x) paste0(x,c('A','B','C','D','E'))),sapply(5:9,function(x) paste0(x,c('A','B','C','D','E','F','G'))))
+results2 <- c(sapply(10:20,function(x) paste0(x,c('A','B','C','D','E','F','G'))))
+results <- c(sapply(c(50,100,200,300),function(x) paste0(results,'_',x)))
+results2 <- c(sapply(c(50,100,200,300),function(x) paste0(results2,'_',x)))
+results <- sort(results)
+results2 <- sort(results2)
+results <- c(results,results2)
+#### Compiler function
+compile_simulation <- function(scenario) {
+name <- as.numeric(gsub("[^0-9.-]", "", substr(scenario,start=0,stop=2)))
+if (substr(scenario,start=nchar(scenario)-1,stop=nchar(scenario))=="50" & name<=4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N50/scenario_',scenario,'.csv'))
+}
+if (substr(scenario,start=nchar(scenario)-1,stop=nchar(scenario))=="50" & name>4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N50/scenario_',scenario,'.csv'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="100" & name<=4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N100/scenario_',scenario,'.csv'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="100" & name>4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N100/scenario_',scenario,'_nodif.csv'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="200" & name<=4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N200/scenario_',scenario,'.csv'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="200" & name>4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N200/scenario_',scenario,'_nodif.csv'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="300" & name<=4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N300/scenario_',scenario,'.csv'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="300" & name>4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N300/scenario_',scenario,'_nodif.csv'))
+}
+if (unique(s$J)==4) {
+if (unique(s$M)==2) {
+a <- data.frame(m.item1=mean(s$item1),m.item2=mean(s$item2),m.item3=mean(s$item3),m.item4=mean(s$item4))
+} else {
+a <- data.frame(m.item1_1=mean(s$item1_1),m.item1_2=mean(s$item1_2),m.item1_3=mean(s$item1_3),
+m.item2_1=mean(s$item2_1),m.item2_2=mean(s$item2_2),m.item2_3=mean(s$item2_3),
+m.item3_1=mean(s$item3_1),m.item3_2=mean(s$item3_2),m.item3_3=mean(s$item3_3),
+m.item4_1=mean(s$item4_1),m.item4_2=mean(s$item4_2),m.item4_3=mean(s$item4_3)
+)
+}
+} else {
+if (unique(s$M)==2) {
+a <- data.frame(m.item1=mean(s$item1),m.item2=mean(s$item2),m.item3=mean(s$item3),m.item4=mean(s$item4),
+m.item5=mean(s$item5),m.item6=mean(s$item6),m.item7=mean(s$item7))
+} else {
+a <- data.frame(m.item1_1=mean(s$item1_1),m.item1_2=mean(s$item1_2),m.item1_3=mean(s$item1_3),
+m.item2_1=mean(s$item2_1),m.item2_2=mean(s$item2_2),m.item2_3=mean(s$item2_3),
+m.item3_1=mean(s$item3_1),m.item3_2=mean(s$item3_2),m.item3_3=mean(s$item3_3),
+m.item4_1=mean(s$item4_1),m.item4_2=mean(s$item4_2),m.item4_3=mean(s$item4_3),
+m.item5_1=mean(s$item5_1),m.item5_2=mean(s$item5_2),m.item5_3=mean(s$item5_3),
+m.item6_1=mean(s$item6_1),m.item6_2=mean(s$item6_2),m.item6_3=mean(s$item6_3),
+m.item7_1=mean(s$item7_1),m.item7_2=mean(s$item7_2),m.item7_3=mean(s$item7_3)
+)
+}
+}
+zz <- substr(scenario,start=0,stop=nchar(scenario)-4)
+b <- data.frame(scenario=zz,
+scenario.type=substr(zz,start=nchar(zz),stop=nchar(zz)),
+N=substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario)),
+J=unique(s$J),
+M=unique(s$M),
+eff.size=unique(s$eff.size),
+nb.dif=unique(s$nb.dif),
+dif.size=unique(s$dif.size)
+)
+z <- data.frame(m.beta=mean(s$beta),
+se.empirical.beta=sd(s$beta),
+se.analytical.beta=mean(s$se.beta),
+m.low.ci.beta=mean(s$low.ci.beta),
+m.high.ci.beta=mean(s$high.ci.beta),
+true.value.in.ci.p=mean(s$true.value.in.ci),
+h0.rejected.p=mean(s$h0.rejected),
+beta.same.sign.truebeta.p=mean(s$beta.same.sign.truebeta,na.rm=T),
+beta.same.sign.truebeta.signif.p=mean(s[s$h0.rejected==1,]$beta.same.sign.truebeta,na.rm=T))
+d <- cbind(b,a,z)
+d$prop.
+return(d)
+}
+#### Compiled results
+res.dat <- compile_simulation('1A_100')
+for (x in results[seq(2,length(results))]) {
+y <- compile_simulation(x)
+res.dat <- bind_rows(res.dat,y)
+}
+res.dat[res.dat$scenario.type=='A','dif.size'] <- -res.dat[res.dat$scenario.type=='A','dif.size']
+res.dat[is.na(res.dat$dif.size),'dif.size'] <- 0
+res.dat[193:417,'nb.dif'] <- 2
+res.dat[417:528,'nb.dif'] <- 3
+View(res.dat)
+res.dat.simple <- res.dat[,c(1:8,13,16:18)]
+res.dat.simple$m.beta <- round(res.dat.simple$m.beta,3)
+res.dat.simple
+#### Create data.frame
+results <- c(sapply(1:4,function(x) paste0(x,c('A','B','C','D','E'))),sapply(5:9,function(x) paste0(x,c('A','B','C','D','E','F','G'))))
+results2 <- c(sapply(10:20,function(x) paste0(x,c('A','B','C','D','E','F','G'))))
+results <- c(sapply(c(50,100,200,300),function(x) paste0(results,'_',x)))
+results2 <- c(sapply(c(50,100,200,300),function(x) paste0(results2,'_',x)))
+results <- sort(results)
+results2 <- sort(results2)
+results <- c(results,results2)
+#### Compiler function
+compile_simulation <- function(scenario) {
+name <- as.numeric(gsub("[^0-9.-]", "", substr(scenario,start=0,stop=2)))
+if (substr(scenario,start=nchar(scenario)-1,stop=nchar(scenario))=="50" & name<=4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N50/scenario_',scenario,'.csv'))
+}
+if (substr(scenario,start=nchar(scenario)-1,stop=nchar(scenario))=="50" & name>4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N50/scenario_',scenario,'.csv'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="100" & name<=4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N100/scenario_',scenario,'.csv'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="100" & name>4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N100/scenario_',scenario,'_nodif.csv'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="200" & name<=4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N200/scenario_',scenario,'.csv'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="200" & name>4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N200/scenario_',scenario,'_nodif.csv'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="300" & name<=4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N300/scenario_',scenario,'.csv'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="300" & name>4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N300/scenario_',scenario,'_nodif.csv'))
+}
+if (unique(s$J)==4) {
+if (unique(s$M)==2) {
+a <- data.frame(m.item1=mean(s$item1),m.item2=mean(s$item2),m.item3=mean(s$item3),m.item4=mean(s$item4))
+} else {
+a <- data.frame(m.item1_1=mean(s$item1_1),m.item1_2=mean(s$item1_2),m.item1_3=mean(s$item1_3),
+m.item2_1=mean(s$item2_1),m.item2_2=mean(s$item2_2),m.item2_3=mean(s$item2_3),
+m.item3_1=mean(s$item3_1),m.item3_2=mean(s$item3_2),m.item3_3=mean(s$item3_3),
+m.item4_1=mean(s$item4_1),m.item4_2=mean(s$item4_2),m.item4_3=mean(s$item4_3)
+)
+}
+} else {
+if (unique(s$M)==2) {
+a <- data.frame(m.item1=mean(s$item1),m.item2=mean(s$item2),m.item3=mean(s$item3),m.item4=mean(s$item4),
+m.item5=mean(s$item5),m.item6=mean(s$item6),m.item7=mean(s$item7))
+} else {
+a <- data.frame(m.item1_1=mean(s$item1_1),m.item1_2=mean(s$item1_2),m.item1_3=mean(s$item1_3),
+m.item2_1=mean(s$item2_1),m.item2_2=mean(s$item2_2),m.item2_3=mean(s$item2_3),
+m.item3_1=mean(s$item3_1),m.item3_2=mean(s$item3_2),m.item3_3=mean(s$item3_3),
+m.item4_1=mean(s$item4_1),m.item4_2=mean(s$item4_2),m.item4_3=mean(s$item4_3),
+m.item5_1=mean(s$item5_1),m.item5_2=mean(s$item5_2),m.item5_3=mean(s$item5_3),
+m.item6_1=mean(s$item6_1),m.item6_2=mean(s$item6_2),m.item6_3=mean(s$item6_3),
+m.item7_1=mean(s$item7_1),m.item7_2=mean(s$item7_2),m.item7_3=mean(s$item7_3)
+)
+}
+}
+N <- ifelse(substr(scenario,start=nchar(scenario)-1,stop=nchar(scenario))=="50",50,substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario)))
+zz <- ifelse(N=="50",substr(scenario,start=0,stop=nchar(scenario)-3),substr(scenario,start=0,stop=nchar(scenario)-4))
+b <- data.frame(scenario=zz,
+scenario.type=substr(zz,start=nchar(zz),stop=nchar(zz)),
+N=N,
+J=unique(s$J),
+M=unique(s$M),
+eff.size=unique(s$eff.size),
+nb.dif=unique(s$nb.dif),
+dif.size=unique(s$dif.size)
+)
+z <- data.frame(m.beta=mean(s$beta),
+se.empirical.beta=sd(s$beta),
+se.analytical.beta=mean(s$se.beta),
+m.low.ci.beta=mean(s$low.ci.beta),
+m.high.ci.beta=mean(s$high.ci.beta),
+true.value.in.ci.p=mean(s$true.value.in.ci),
+h0.rejected.p=mean(s$h0.rejected),
+beta.same.sign.truebeta.p=mean(s$beta.same.sign.truebeta,na.rm=T),
+beta.same.sign.truebeta.signif.p=mean(s[s$h0.rejected==1,]$beta.same.sign.truebeta,na.rm=T))
+d <- cbind(b,a,z)
+d$prop.
+return(d)
+}
+#### Compiled results
+res.dat <- compile_simulation('1A_100')
+for (x in results[seq(2,length(results))]) {
+y <- compile_simulation(x)
+res.dat <- bind_rows(res.dat,y)
+}
+#### Create data.frame
+results <- c(sapply(1:4,function(x) paste0(x,c('A','B','C','D','E'))),sapply(5:9,function(x) paste0(x,c('A','B','C','D','E','F','G'))))
+results2 <- c(sapply(10:20,function(x) paste0(x,c('A','B','C','D','E','F','G'))))
+results <- c(sapply(c(50,100,200,300),function(x) paste0(results,'_',x)))
+results2 <- c(sapply(c(50,100,200,300),function(x) paste0(results2,'_',x)))
+results <- sort(results)
+results2 <- sort(results2)
+results <- c(results,results2)
+#### Compiler function
+compile_simulation <- function(scenario) {
+name <- as.numeric(gsub("[^0-9.-]", "", substr(scenario,start=0,stop=2)))
+if (substr(scenario,start=nchar(scenario)-1,stop=nchar(scenario))=="50" & name<=4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N50/scenario_',scenario,'.csv'))
+}
+if (substr(scenario,start=nchar(scenario)-1,stop=nchar(scenario))=="50" & name>4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N50/scenario_',scenario,'.csv'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="100" & name<=4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N100/scenario_',scenario,'.csv'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="100" & name>4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N100/scenario_',scenario,'_nodif.csv'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="200" & name<=4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N200/scenario_',scenario,'.csv'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="200" & name>4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N200/scenario_',scenario,'_nodif.csv'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="300" & name<=4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N300/scenario_',scenario,'.csv'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="300" & name>4) {
+s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N300/scenario_',scenario,'_nodif.csv'))
+}
+if (unique(s$J)==4) {
+if (unique(s$M)==2) {
+a <- data.frame(m.item1=mean(s$item1),m.item2=mean(s$item2),m.item3=mean(s$item3),m.item4=mean(s$item4))
+} else {
+a <- data.frame(m.item1_1=mean(s$item1_1),m.item1_2=mean(s$item1_2),m.item1_3=mean(s$item1_3),
+m.item2_1=mean(s$item2_1),m.item2_2=mean(s$item2_2),m.item2_3=mean(s$item2_3),
+m.item3_1=mean(s$item3_1),m.item3_2=mean(s$item3_2),m.item3_3=mean(s$item3_3),
+m.item4_1=mean(s$item4_1),m.item4_2=mean(s$item4_2),m.item4_3=mean(s$item4_3)
+)
+}
+} else {
+if (unique(s$M)==2) {
+a <- data.frame(m.item1=mean(s$item1),m.item2=mean(s$item2),m.item3=mean(s$item3),m.item4=mean(s$item4),
+m.item5=mean(s$item5),m.item6=mean(s$item6),m.item7=mean(s$item7))
+} else {
+a <- data.frame(m.item1_1=mean(s$item1_1),m.item1_2=mean(s$item1_2),m.item1_3=mean(s$item1_3),
+m.item2_1=mean(s$item2_1),m.item2_2=mean(s$item2_2),m.item2_3=mean(s$item2_3),
+m.item3_1=mean(s$item3_1),m.item3_2=mean(s$item3_2),m.item3_3=mean(s$item3_3),
+m.item4_1=mean(s$item4_1),m.item4_2=mean(s$item4_2),m.item4_3=mean(s$item4_3),
+m.item5_1=mean(s$item5_1),m.item5_2=mean(s$item5_2),m.item5_3=mean(s$item5_3),
+m.item6_1=mean(s$item6_1),m.item6_2=mean(s$item6_2),m.item6_3=mean(s$item6_3),
+m.item7_1=mean(s$item7_1),m.item7_2=mean(s$item7_2),m.item7_3=mean(s$item7_3)
+)
+}
+}
+N <- ifelse(substr(scenario,start=nchar(scenario)-1,stop=nchar(scenario))=="50","50",substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario)))
+zz <- ifelse(N=="50",substr(scenario,start=0,stop=nchar(scenario)-3),substr(scenario,start=0,stop=nchar(scenario)-4))
+b <- data.frame(scenario=zz,
+scenario.type=substr(zz,start=nchar(zz),stop=nchar(zz)),
+N=N,
+J=unique(s$J),
+M=unique(s$M),
+eff.size=unique(s$eff.size),
+nb.dif=unique(s$nb.dif),
+dif.size=unique(s$dif.size)
+)
+z <- data.frame(m.beta=mean(s$beta),
+se.empirical.beta=sd(s$beta),
+se.analytical.beta=mean(s$se.beta),
+m.low.ci.beta=mean(s$low.ci.beta),
+m.high.ci.beta=mean(s$high.ci.beta),
+true.value.in.ci.p=mean(s$true.value.in.ci),
+h0.rejected.p=mean(s$h0.rejected),
+beta.same.sign.truebeta.p=mean(s$beta.same.sign.truebeta,na.rm=T),
+beta.same.sign.truebeta.signif.p=mean(s[s$h0.rejected==1,]$beta.same.sign.truebeta,na.rm=T))
+d <- cbind(b,a,z)
+d$prop.
+return(d)
+}
+#### Compiled results
+res.dat <- compile_simulation('1A_100')
+for (x in results[seq(2,length(results))]) {
+y <- compile_simulation(x)
+res.dat <- bind_rows(res.dat,y)
+}
+res.dat[res.dat$scenario.type=='A','dif.size'] <- -res.dat[res.dat$scenario.type=='A','dif.size']
+res.dat[is.na(res.dat$dif.size),'dif.size'] <- 0
+res.dat[193:417,'nb.dif'] <- 2
+res.dat[417:528,'nb.dif'] <- 3
+View(res.dat)
+res.dat.simple <- res.dat[,c(1:8,13,16:18)]
+res.dat.simple$m.beta <- round(res.dat.simple$m.beta,3)
+res.dat.simple
+#### Create data.frame
+results <- c(sapply(1:4,function(x) paste0(x,c('A','B','C','D','E'))),sapply(5:9,function(x) paste0(x,c('A','B','C','D','E','F','G'))))
+results2 <- c(sapply(10:20,function(x) paste0(x,c('A','B','C','D','E','F','G'))))
+results <- c(sapply(c(50,100,200,300),function(x) paste0(results,'_',x)))
+results2 <- c(sapply(c(50,100,200,300),function(x) paste0(results2,'_',x)))
+results <- sort(results)
+results2 <- sort(results2)
+results <- c(results,results2)[81:528]
+#### Compiler function
+compile_simulation2 <- function(scenario) {
+name <- as.numeric(gsub("[^0-9.-]", "", substr(scenario,start=0,stop=2)))
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="50" & name>4) {
+s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/DIF/N50/',scenario,'.xls'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="100" & name>4) {
+s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/DIF/N100/',scenario,'.xls'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="200" & name>4) {
+s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/DIF/N200/',scenario,'.xls'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="300" & name>4) {
+s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/DIF/N300/',scenario,'.xls'))
+}
+J <- max(which(sapply(1:7,function(x) paste0('item',x) %in% colnames(s) | paste0('item',x,'_1') %in% colnames(s))))
+M <- 1+sum(sapply(1:3,function(x) paste0('item1_',x) %in% colnames(s) ))
+if (M==1) {M <- 2}
+nb.dif <- max(which(sapply(1:3,function(x) paste0('dif',x) %in% colnames(s) | paste0('dif',x,'_1') %in% colnames(s))))
+if (J==4) {
+if (M==2) {
+a <- data.frame(m.item1=mean(s$item1),m.item2=mean(s$item2),m.item3=mean(s$item3),m.item4=mean(s$item4))
+} else {
+a <- data.frame(m.item1_1=mean(s$item1_1),m.item1_2=mean(s$item1_2),m.item1_3=mean(s$item1_3),
+m.item2_1=mean(s$item2_1),m.item2_2=mean(s$item2_2),m.item2_3=mean(s$item2_3),
+m.item3_1=mean(s$item3_1),m.item3_2=mean(s$item3_2),m.item3_3=mean(s$item3_3),
+m.item4_1=mean(s$item4_1),m.item4_2=mean(s$item4_2),m.item4_3=mean(s$item4_3)
+)
+}
+} else {
+if (M==2) {
+a <- data.frame(m.item1=mean(s$item1),m.item2=mean(s$item2),m.item3=mean(s$item3),m.item4=mean(s$item4),
+m.item5=mean(s$item5),m.item6=mean(s$item6),m.item7=mean(s$item7))
+} else {
+a <- data.frame(m.item1_1=mean(s$item1_1),m.item1_2=mean(s$item1_2),m.item1_3=mean(s$item1_3),
+m.item2_1=mean(s$item2_1),m.item2_2=mean(s$item2_2),m.item2_3=mean(s$item2_3),
+m.item3_1=mean(s$item3_1),m.item3_2=mean(s$item3_2),m.item3_3=mean(s$item3_3),
+m.item4_1=mean(s$item4_1),m.item4_2=mean(s$item4_2),m.item4_3=mean(s$item4_3),
+m.item5_1=mean(s$item5_1),m.item5_2=mean(s$item5_2),m.item5_3=mean(s$item5_3),
+m.item6_1=mean(s$item6_1),m.item6_2=mean(s$item6_2),m.item6_3=mean(s$item6_3),
+m.item7_1=mean(s$item7_1),m.item7_2=mean(s$item7_2),m.item7_3=mean(s$item7_3)
+)
+}
+}
+N <- ifelse(substr(scenario,start=nchar(scenario)-1,stop=nchar(scenario))=="50","50",substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario)))
+zz <- ifelse(N=="50",substr(scenario,start=0,stop=nchar(scenario)-3),substr(scenario,start=0,stop=nchar(scenario)-4))
+eff.size <- unique(res.dat[res.dat$scenario==zz & res.dat$N==N,'eff.size'])
+dif.size <- unique(res.dat[res.dat$scenario==zz & res.dat$N==N,'dif.size'])
+b <- data.frame(scenario=zz,
+scenario.type=substr(zz,start=nchar(zz),stop=nchar(zz)),
+N=N,
+J=J,
+M=M,
+eff.size=eff.size,
+nb.dif=nb.dif,
+dif.size=dif.size
+)
+true.value.in.ci <- eff.size <= s$beta+1.96*s$se_beta & eff.size >= s$beta-1.96*s$se_beta
+beta.same.sign.truebeta.p <- ifelse(rep(eff.size,nrow(s))==0,NA,(rep(eff.size,nrow(s))/s$beta)>0)
+num.reject <- which((s$beta-1.96*s$se_beta)>0 | (s$beta+1.96*s$se_beta)<0)
+z <- data.frame(m.beta=mean(s$beta),
+se.empirical.beta=sd(s$beta),
+se.analytical.beta=mean(s$se_beta),
+m.low.ci.beta=mean(s$beta-1.96*s$se_beta),
+m.high.ci.beta=mean(s$beta+1.96*s$se_beta),
+true.value.in.ci.p=mean(true.value.in.ci),
+h0.rejected.p=mean( (s$beta-1.96*s$se_beta)>0 | (s$beta+1.96*s$se_beta)<0 ),
+beta.same.sign.truebeta.p=mean(beta.same.sign.truebeta.p),
+beta.same.sign.truebeta.signif.p=mean(beta.same.sign.truebeta.p[num.reject])
+)
+d <- cbind(b,a,z)
+d$prop.
+return(d)
+}
+#### Compiled results
+res.dat.dif <- compile_simulation2('5A_100')
+for (x in results[seq(2,length(results))]) {
+y <- compile_simulation2(x)
+res.dat.dif <- bind_rows(res.dat.dif,y)
+}
+#### Create data.frame
+results <- c(sapply(1:4,function(x) paste0(x,c('A','B','C','D','E'))),sapply(5:9,function(x) paste0(x,c('A','B','C','D','E','F','G'))))
+results2 <- c(sapply(10:20,function(x) paste0(x,c('A','B','C','D','E','F','G'))))
+results <- c(sapply(c(50,100,200,300),function(x) paste0(results,'_',x)))
+results2 <- c(sapply(c(50,100,200,300),function(x) paste0(results2,'_',x)))
+results <- sort(results)
+results2 <- sort(results2)
+results <- c(results,results2)[81:528]
+#### Compiler function
+compile_simulation2 <- function(scenario) {
+name <- as.numeric(gsub("[^0-9.-]", "", substr(scenario,start=0,stop=2)))
+if (substr(scenario,start=nchar(scenario)-1,stop=nchar(scenario))=="50" & name>4) {
+s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/DIF/N50/',scenario,'.xls'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="100" & name>4) {
+s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/DIF/N100/',scenario,'.xls'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="200" & name>4) {
+s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/DIF/N200/',scenario,'.xls'))
+}
+if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="300" & name>4) {
+s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/DIF/N300/',scenario,'.xls'))
+}
+J <- max(which(sapply(1:7,function(x) paste0('item',x) %in% colnames(s) | paste0('item',x,'_1') %in% colnames(s))))
+M <- 1+sum(sapply(1:3,function(x) paste0('item1_',x) %in% colnames(s) ))
+if (M==1) {M <- 2}
+nb.dif <- max(which(sapply(1:3,function(x) paste0('dif',x) %in% colnames(s) | paste0('dif',x,'_1') %in% colnames(s))))
+if (J==4) {
+if (M==2) {
+a <- data.frame(m.item1=mean(s$item1),m.item2=mean(s$item2),m.item3=mean(s$item3),m.item4=mean(s$item4))
+} else {
+a <- data.frame(m.item1_1=mean(s$item1_1),m.item1_2=mean(s$item1_2),m.item1_3=mean(s$item1_3),
+m.item2_1=mean(s$item2_1),m.item2_2=mean(s$item2_2),m.item2_3=mean(s$item2_3),
+m.item3_1=mean(s$item3_1),m.item3_2=mean(s$item3_2),m.item3_3=mean(s$item3_3),
+m.item4_1=mean(s$item4_1),m.item4_2=mean(s$item4_2),m.item4_3=mean(s$item4_3)
+)
+}
+} else {
+if (M==2) {
+a <- data.frame(m.item1=mean(s$item1),m.item2=mean(s$item2),m.item3=mean(s$item3),m.item4=mean(s$item4),
+m.item5=mean(s$item5),m.item6=mean(s$item6),m.item7=mean(s$item7))
+} else {
+a <- data.frame(m.item1_1=mean(s$item1_1),m.item1_2=mean(s$item1_2),m.item1_3=mean(s$item1_3),
+m.item2_1=mean(s$item2_1),m.item2_2=mean(s$item2_2),m.item2_3=mean(s$item2_3),
+m.item3_1=mean(s$item3_1),m.item3_2=mean(s$item3_2),m.item3_3=mean(s$item3_3),
+m.item4_1=mean(s$item4_1),m.item4_2=mean(s$item4_2),m.item4_3=mean(s$item4_3),
+m.item5_1=mean(s$item5_1),m.item5_2=mean(s$item5_2),m.item5_3=mean(s$item5_3),
+m.item6_1=mean(s$item6_1),m.item6_2=mean(s$item6_2),m.item6_3=mean(s$item6_3),
+m.item7_1=mean(s$item7_1),m.item7_2=mean(s$item7_2),m.item7_3=mean(s$item7_3)
+)
+}
+}
+N <- ifelse(substr(scenario,start=nchar(scenario)-1,stop=nchar(scenario))=="50","50",substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario)))
+zz <- ifelse(N=="50",substr(scenario,start=0,stop=nchar(scenario)-3),substr(scenario,start=0,stop=nchar(scenario)-4))
+eff.size <- unique(res.dat[res.dat$scenario==zz & res.dat$N==N,'eff.size'])
+dif.size <- unique(res.dat[res.dat$scenario==zz & res.dat$N==N,'dif.size'])
+b <- data.frame(scenario=zz,
+scenario.type=substr(zz,start=nchar(zz),stop=nchar(zz)),
+N=N,
+J=J,
+M=M,
+eff.size=eff.size,
+nb.dif=nb.dif,
+dif.size=dif.size
+)
+true.value.in.ci <- eff.size <= s$beta+1.96*s$se_beta & eff.size >= s$beta-1.96*s$se_beta
+beta.same.sign.truebeta.p <- ifelse(rep(eff.size,nrow(s))==0,NA,(rep(eff.size,nrow(s))/s$beta)>0)
+num.reject <- which((s$beta-1.96*s$se_beta)>0 | (s$beta+1.96*s$se_beta)<0)
+z <- data.frame(m.beta=mean(s$beta),
+se.empirical.beta=sd(s$beta),
+se.analytical.beta=mean(s$se_beta),
+m.low.ci.beta=mean(s$beta-1.96*s$se_beta),
+m.high.ci.beta=mean(s$beta+1.96*s$se_beta),
+true.value.in.ci.p=mean(true.value.in.ci),
+h0.rejected.p=mean( (s$beta-1.96*s$se_beta)>0 | (s$beta+1.96*s$se_beta)<0 ),
+beta.same.sign.truebeta.p=mean(beta.same.sign.truebeta.p),
+beta.same.sign.truebeta.signif.p=mean(beta.same.sign.truebeta.p[num.reject])
+)
+d <- cbind(b,a,z)
+d$prop.
+return(d)
+}
+#### Compiled results
+res.dat.dif <- compile_simulation2('5A_100')
+for (x in results[seq(2,length(results))]) {
+y <- compile_simulation2(x)
+res.dat.dif <- bind_rows(res.dat.dif,y)
+}
+res.dat$bias <- res.dat$eff.size-res.dat$m.beta
+res.dat.dif$bias <- res.dat.dif$eff.size-res.dat.dif$m.beta
+res.dat.dif
+View(res.dat.dif)
+mean(res.dat.dif[res.dat.dif$N==50,]$bias)
+mean(res.dat.dif[res.dat.dif$N==50,]$bias,na.rm=T)
+mean(res.dat.dif$bias,na.rm=T)
+mean(res.dat.dif[res.dat.dif$N!=50,]$bias,na.rm=T)
+mean(res.dat[res.dat$N!=50,]$bias,na.rm=T)
+mean(res.dat[res.dat$N!=50,]$bias)
diff --git a/Scripts/Analysis/DIF/pcm_dif_50_J7.do b/Scripts/Analysis/DIF/pcm_dif_50_J7.do
index bd73ae7..da4bdab 100644
--- a/Scripts/Analysis/DIF/pcm_dif_50_J7.do
+++ b/Scripts/Analysis/DIF/pcm_dif_50_J7.do
@@ -1123,6 +1123,7 @@ forvalues k=1/1000 {
     if (`nummoda'==2) {
       qui recode item`j' (`checker'=0) (`checker2'=1)
     }
+    di `nummoda'
     if (`nummoda'==3) {
     if (`checker'!=0) {
       qui recode item`j' (`checker'=0) (`checker2'=1) (`checker3'=2)
@@ -1200,7 +1201,7 @@ forvalues i=1/`nbitems' {
     if (`checker2'==3) {
       local mod = "`mod'"+"(1.item`i'<-THETA@1 tt)(2.item`i'<-THETA@2 tt)(3.item`i'<-THETA@3 tt)"
   }
-    else if (`checker2'==3) {
+    else if (`checker2'==2) {
       local mod = "`mod'"+"(1.item`i'<-THETA@1 tt)(2.item`i'<-THETA@2 tt)"
     }
     else {
@@ -1580,6 +1581,7 @@ putexcel A1=matrix(outmat), colnames
 
 
 
+
 ************************************************** Scenario 20: J = 7 items / M = 4 modalities / DIF size 0.3 x3
 local N = "50"
 foreach Nnn in `N' {
@@ -1731,6 +1733,7 @@ forvalues k=1/1000 {
     if (`nummoda'==2) {
       qui recode item`j' (`checker'=0) (`checker2'=1)
     }
+    di `nummoda'
     if (`nummoda'==3) {
     if (`checker'!=0) {
       qui recode item`j' (`checker'=0) (`checker2'=1) (`checker3'=2)
@@ -1808,7 +1811,7 @@ forvalues i=1/`nbitems' {
     if (`checker2'==3) {
       local mod = "`mod'"+"(1.item`i'<-THETA@1 tt)(2.item`i'<-THETA@2 tt)(3.item`i'<-THETA@3 tt)"
   }
-    else if (`checker2'==3) {
+    else if (`checker2'==2) {
       local mod = "`mod'"+"(1.item`i'<-THETA@1 tt)(2.item`i'<-THETA@2 tt)"
     }
     else {