diff --git a/RProject/.RData b/RProject/.RData
index 57c198b..cf9da9e 100644
Binary files a/RProject/.RData and b/RProject/.RData differ
diff --git a/RProject/.Rhistory b/RProject/.Rhistory
index 8f38e75..688f60e 100644
--- a/RProject/.Rhistory
+++ b/RProject/.Rhistory
@@ -1,512 +1,512 @@
-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",])
+resp <- df[,c(dif.items.list)]
+}
+}
+else {
+resp <- df[,c(nodif.items.list)]
+}
+print(grp)
+tam1 <- TAM::tam.mml(resp=resp,Y=grp,irtmodel = model,est.variance = T,verbose=F)
+}
+}
+return(summary(tam1))
+}
+rosali(dat=dat,items=1:7,group="TT")
+rosali <- function(dat=NULL,items=NULL,group=NULL) {
+nbitems <- length(items)
+items2 <- items
+# Un seul groupe
+if (length(group)!=1) {
+stop('Only one variable can be used for the group option')
+}
+# Recoder groupe en 0/1
+dat[,group] <- as.factor(dat[,group])
+if (!(all(names(levels(dat[,group]))==c("0","1")))) {
+levels(dat[,group]) <- c("0","1")
+}
+dat[,group] <- as.numeric(dat[,group])
+model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
+model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
+return(c(model_a,model_b))
+}
+rosali(dat=dat,items=1:7,group="TT")
+rosali <- function(dat=NULL,items=NULL,group=NULL) {
+nbitems <- length(items)
+items2 <- items
+# Un seul groupe
+if (length(group)!=1) {
+stop('Only one variable can be used for the group option')
+}
+# Recoder groupe en 0/1
+dat[,group] <- as.factor(dat[,group])
+if (!(all(names(levels(dat[,group]))==c("0","1")))) {
+levels(dat[,group]) <- c("0","1")
+}
+dat[,group] <- as.numeric(dat[,group])-1
+model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
+model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
+return(c(model_a,model_b))
+}
+rosali(dat=dat,items=1:7,group="TT")
+## File Name: pcm.R
+## File version: 1.0
+#' @import TAM
+#' @export
+pcm <- function(df=NULL,items=NULL,group=NULL,model="PCM2",method="MML",dif.items=NULL) {
+# Prepare analysis
+if (is.null(items)) {
+nbitems <- sum(sapply(1:100,function(x) paste0('item',x)) %in% colnames(df))
+items <- paste0('item',seq(1,nbitems))
+resp <- df[,items]
+}
+else {
+nbitems <- length(items)
+resp <- df[,paste0("item",items)]
+}
+if (is.null(group)) {
+grp <- NULL
+}
+else {
+grp <- df[,group]
+df$grp <- grp
+if (!is.null(dif.items)) {
+for (i in dif.items) {
+df[,paste0('item',i,'_noTT')] <- NA
+df[,paste0('item',i,'_TT')] <- NA
+df[df$grp==0,paste0('item',i,'_noTT')] <- df[df$grp==0,paste0('item',i)]
+df[df$grp==1,paste0('item',i,'_TT')] <- df[df$grp==1,paste0('item',i)]
+}
+}
+}
+# Analyze
+if (method=='MML') {
+if (is.null(dif.items)) {
+tam1 <-  TAM::tam.mml(resp=resp,Y=grp,irtmodel = model,est.variance = T,verbose=F)
+}
+else {
+nodif.items <- setdiff(items,dif.items)
+nodif.items.list <- paste0("item",nodif.items)
+dif.items.list <- c(paste0("item",dif.items,"_noTT"),paste0("item",dif.items,"_TT"))
+if (!is.null(dif.items)) {
+if (length(nodif.items)!=0) {
+resp <- df[,c(nodif.items.list,dif.items.list)]
+}
+else {
+resp <- df[,c(dif.items.list)]
+}
+}
+else {
+resp <- df[,c(nodif.items.list)]
+}
+tam1 <- TAM::tam.mml(resp=resp,Y=grp,irtmodel = model,est.variance = T,verbose=F)
+}
+}
+return(summary(tam1))
+}
+rosali <- function(dat=NULL,items=NULL,group=NULL) {
+nbitems <- length(items)
+items2 <- items
+# Un seul groupe
+if (length(group)!=1) {
+stop('Only one variable can be used for the group option')
+}
+# Recoder groupe en 0/1
+dat[,group] <- as.factor(dat[,group])
+if (!(all(names(levels(dat[,group]))==c("0","1")))) {
+levels(dat[,group]) <- c("0","1")
+}
+dat[,group] <- as.numeric(dat[,group])-1
+model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
+model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
+return(c(model_a))
+}
+rosali(dat=dat,items=1:7,group="TT")
+rosali <- function(dat=NULL,items=NULL,group=NULL) {
+nbitems <- length(items)
+items2 <- items
+# Un seul groupe
+if (length(group)!=1) {
+stop('Only one variable can be used for the group option')
+}
+# Recoder groupe en 0/1
+dat[,group] <- as.factor(dat[,group])
+if (!(all(names(levels(dat[,group]))==c("0","1")))) {
+levels(dat[,group]) <- c("0","1")
+}
+dat[,group] <- as.numeric(dat[,group])-1
+model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
+model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
+return(c(model_b))
+}
+rosali(dat=dat,items=1:7,group="TT")
+rosali <- function(dat=NULL,items=NULL,group=NULL) {
+nbitems <- length(items)
+items2 <- items
+# Un seul groupe
+if (length(group)!=1) {
+stop('Only one variable can be used for the group option')
+}
+# Recoder groupe en 0/1
+dat[,group] <- as.factor(dat[,group])
+if (!(all(names(levels(dat[,group]))==c("0","1")))) {
+levels(dat[,group]) <- c("0","1")
+}
+dat[,group] <- as.numeric(dat[,group])-1
+model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
+model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
+}
+rosali(dat=dat,items=1:7,group="TT")
+model_b <- pcm(df=dat,items=1:7,group = "TT",dif.items = 1:7)
+model_a <- pcm(df=dat,items=1:7,group = "TT",dif.items = NULL)
+lrt
+lrt''
+??lrt
+loglik(model_a)
+logLik(model_a)
+## File Name: pcm.R
+## File version: 1.0
+#' @import TAM
+#' @export
+pcm <- function(df=NULL,items=NULL,group=NULL,model="PCM2",method="MML",dif.items=NULL) {
+# Prepare analysis
+if (is.null(items)) {
+nbitems <- sum(sapply(1:100,function(x) paste0('item',x)) %in% colnames(df))
+items <- paste0('item',seq(1,nbitems))
+resp <- df[,items]
+}
+else {
+nbitems <- length(items)
+resp <- df[,paste0("item",items)]
+}
+if (is.null(group)) {
+grp <- NULL
+}
+else {
+grp <- df[,group]
+df$grp <- grp
+if (!is.null(dif.items)) {
+for (i in dif.items) {
+df[,paste0('item',i,'_noTT')] <- NA
+df[,paste0('item',i,'_TT')] <- NA
+df[df$grp==0,paste0('item',i,'_noTT')] <- df[df$grp==0,paste0('item',i)]
+df[df$grp==1,paste0('item',i,'_TT')] <- df[df$grp==1,paste0('item',i)]
+}
+}
+}
+# Analyze
+if (method=='MML') {
+if (is.null(dif.items)) {
+tam1 <-  TAM::tam.mml(resp=resp,Y=grp,irtmodel = model,est.variance = T,verbose=F)
+}
+else {
+nodif.items <- setdiff(items,dif.items)
+nodif.items.list <- paste0("item",nodif.items)
+dif.items.list <- c(paste0("item",dif.items,"_noTT"),paste0("item",dif.items,"_TT"))
+if (!is.null(dif.items)) {
+if (length(nodif.items)!=0) {
+resp <- df[,c(nodif.items.list,dif.items.list)]
+}
+else {
+resp <- df[,c(dif.items.list)]
+}
+}
+else {
+resp <- df[,c(nodif.items.list)]
+}
+tam1 <- TAM::tam.mml(resp=resp,Y=grp,irtmodel = model,est.variance = T,verbose=F)
+}
+}
+return(tam1)
+}
+rosali <- function(dat=NULL,items=NULL,group=NULL) {
+nbitems <- length(items)
+items2 <- items
+# Un seul groupe
+if (length(group)!=1) {
+stop('Only one variable can be used for the group option')
+}
+# Recoder groupe en 0/1
+dat[,group] <- as.factor(dat[,group])
+if (!(all(names(levels(dat[,group]))==c("0","1")))) {
+levels(dat[,group]) <- c("0","1")
+}
+dat[,group] <- as.numeric(dat[,group])-1
+model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
+model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
+}
+model_b <- pcm(df=dat,items=1:7,group = "TT",dif.items = 1:7)
+model_a <- pcm(df=dat,items=1:7,group = "TT",dif.items = NULL)
+logLik(model_a)
+logLik(model_b)
+library(lmtest)
+lrtest(model_a,model_b)
+lrtest(model_a,model_b)$p
+ss <- lrtest(model_a,model_b)
+ss$`Pr(>Chisq)`
+ss$`Pr(>Chisq)`[2]
+a <- c()
+a
+c(a,1)
+which(1:4,3)
+?which
+which.min(1:4)
+which.min(2:4)
+rosali <- function(dat=NULL,items=NULL,group=NULL) {
+nbitems <- length(items)
+items2 <- items
+# Un seul groupe
+if (length(group)!=1) {
+stop('Only one variable can be used for the group option')
+}
+# Recoder groupe en 0/1
+dat[,group] <- as.factor(dat[,group])
+if (!(all(names(levels(dat[,group]))==c("0","1")))) {
+levels(dat[,group]) <- c("0","1")
+}
+dat[,group] <- as.numeric(dat[,group])-1
+model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
+model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
+ltest <- lmtest::lrtest(model_a,model_b)
+if (ltest$`Pr(>Chisq)`[2]<1) {
+model_c <- model_b
+difit <- c()
+a <- sapply(1:nbitems,function(x)
+lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
+while (min(a)<0.4) {
+difit <- c(difit,which.min(a))
+model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
+a <- sapply(1:nbitems,function(x)
+lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
+}
+}
+}
+rosali(dat=dat,items=1:7,group="TT")
+rosali <- function(dat=NULL,items=NULL,group=NULL) {
+nbitems <- length(items)
+items2 <- items
+# Un seul groupe
+if (length(group)!=1) {
+stop('Only one variable can be used for the group option')
+}
+# Recoder groupe en 0/1
+dat[,group] <- as.factor(dat[,group])
+if (!(all(names(levels(dat[,group]))==c("0","1")))) {
+levels(dat[,group]) <- c("0","1")
+}
+dat[,group] <- as.numeric(dat[,group])-1
+model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
+model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
+ltest <- lmtest::lrtest(model_a,model_b)
+if (ltest$`Pr(>Chisq)`[2]<1) {
+model_c <- model_b
+difit <- c()
+a <- sapply(1:nbitems,function(x)
+lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
+while (min(a)<0.4) {
+difit <- c(difit,which.min(a))
+model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
+print(difit)
+a <- sapply(1:nbitems,function(x)
+lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
+}
+}
+}
+rosali(dat=dat,items=1:7,group="TT")
+rosali <- function(dat=NULL,items=NULL,group=NULL) {
+nbitems <- length(items)
+items2 <- items
+# Un seul groupe
+if (length(group)!=1) {
+stop('Only one variable can be used for the group option')
+}
+# Recoder groupe en 0/1
+dat[,group] <- as.factor(dat[,group])
+if (!(all(names(levels(dat[,group]))==c("0","1")))) {
+levels(dat[,group]) <- c("0","1")
+}
+dat[,group] <- as.numeric(dat[,group])-1
+model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
+model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
+ltest <- lmtest::lrtest(model_a,model_b)
+if (ltest$`Pr(>Chisq)`[2]<1) {
+model_c <- model_b
+difit <- c()
+a <- sapply(1:nbitems,function(x)
+lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
+while (min(a)<0.4) {
+difit <- c(difit,which.min(a))
+model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
+print(difit)
+a <- sapply(1:nbitems,function(x)
+lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
+print(a)
+}
+}
+}
+rosali(dat=dat,items=1:7,group="TT")
+rosali <- function(dat=NULL,items=NULL,group=NULL) {
+nbitems <- length(items)
+items2 <- items
+# Un seul groupe
+if (length(group)!=1) {
+stop('Only one variable can be used for the group option')
+}
+# Recoder groupe en 0/1
+dat[,group] <- as.factor(dat[,group])
+if (!(all(names(levels(dat[,group]))==c("0","1")))) {
+levels(dat[,group]) <- c("0","1")
+}
+dat[,group] <- as.numeric(dat[,group])-1
+model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
+model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
+ltest <- lmtest::lrtest(model_a,model_b)
+if (ltest$`Pr(>Chisq)`[2]<1) {
+model_c <- model_b
+difit <- c()
+a <- sapply(1:nbitems,function(x)
+lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
+while (min(a)<0.4) {
+difit <- c(difit,which.min(a))
+model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
+print(difit)
+a <- sapply(1:nbitems,function(x)
+ifelse(x%in%difit,lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2],1))
+print(a)
+}
+}
+}
+rosali(dat=dat,items=1:7,group="TT")
+rosali <- function(dat=NULL,items=NULL,group=NULL) {
+nbitems <- length(items)
+items2 <- items
+# Un seul groupe
+if (length(group)!=1) {
+stop('Only one variable can be used for the group option')
+}
+# Recoder groupe en 0/1
+dat[,group] <- as.factor(dat[,group])
+if (!(all(names(levels(dat[,group]))==c("0","1")))) {
+levels(dat[,group]) <- c("0","1")
+}
+dat[,group] <- as.numeric(dat[,group])-1
+model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
+model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
+ltest <- lmtest::lrtest(model_a,model_b)
+if (ltest$`Pr(>Chisq)`[2]<1) {
+model_c <- model_b
+difit <- c()
+a <- sapply(1:nbitems,function(x)
+lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
+while (min(a)<0.4) {
+difit <- c(difit,which.min(a))
+model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
+print(difit)
+a <- sapply(1:nbitems,function(x)
+ifelse(!(x%in%difit),lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2],1))
+print(a)
+}
+}
+}
+rosali(dat=dat,items=1:7,group="TT")
+rosali <- function(dat=NULL,items=NULL,group=NULL) {
+nbitems <- length(items)
+items2 <- items
+# Un seul groupe
+if (length(group)!=1) {
+stop('Only one variable can be used for the group option')
+}
+# Recoder groupe en 0/1
+dat[,group] <- as.factor(dat[,group])
+if (!(all(names(levels(dat[,group]))==c("0","1")))) {
+levels(dat[,group]) <- c("0","1")
+}
+dat[,group] <- as.numeric(dat[,group])-1
+model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
+model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
+ltest <- lmtest::lrtest(model_a,model_b)
+difit <- c()
+if (ltest$`Pr(>Chisq)`[2]<0.05) {
+model_c <- model_b
+a <- sapply(1:nbitems,function(x)
+lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
+while (min(a)<0.05) {
+difit <- c(difit,which.min(a))
+model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
+print(difit)
+a <- sapply(1:nbitems,function(x)
+ifelse(!(x%in%difit),lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2],1))
+}
+}
+return(difit)
+}
+rosali(dat=dat,items=1:7,group="TT")
+dat <- read.csv("/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N200/scenario_20A_200.csv")
+dat <- dat[dat$replication==2,]
+rosali(dat=dat,items=1:7,group="TT")
+rosali <- function(dat=NULL,items=NULL,group=NULL) {
+nbitems <- length(items)
+items2 <- items
+# Un seul groupe
+if (length(group)!=1) {
+stop('Only one variable can be used for the group option')
+}
+# Recoder groupe en 0/1
+dat[,group] <- as.factor(dat[,group])
+if (!(all(names(levels(dat[,group]))==c("0","1")))) {
+levels(dat[,group]) <- c("0","1")
+}
+dat[,group] <- as.numeric(dat[,group])-1
+model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
+model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
+ltest <- lmtest::lrtest(model_a,model_b)
+difit <- c()
+if (ltest$`Pr(>Chisq)`[2]<0.05) {
+model_c <- model_b
+a <- sapply(1:nbitems,function(x)
+lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
+while (min(a)<0.05) {
+difit <- c(difit,which.min(a))
+model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
+a <- sapply(1:nbitems,function(x)
+ifelse(!(x%in%difit),lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2],1))
+}
+}
+return(difit)
+}
+rosali(dat=dat,items=1:7,group="TT")
+dat <- read.csv("/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N200/scenario_20A_200.csv")
+dat <- dat[dat$replication==3,]
+rosali(dat=dat,items=1:7,group="TT")
+rosali <- function(dat=NULL,items=NULL,group=NULL) {
+nbitems <- length(items)
+items2 <- items
+# Un seul groupe
+if (length(group)!=1) {
+stop('Only one variable can be used for the group option')
+}
+# Recoder groupe en 0/1
+dat[,group] <- as.factor(dat[,group])
+if (!(all(names(levels(dat[,group]))==c("0","1")))) {
+levels(dat[,group]) <- c("0","1")
+}
+dat[,group] <- as.numeric(dat[,group])-1
+model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
+model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
+ltest <- lmtest::lrtest(model_a,model_b)
+difit <- c()
+if (ltest$`Pr(>Chisq)`[2]<0.05) {
+model_c <- model_b
+a <- sapply(1:nbitems,function(x)
+lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
+k <- 0
+while (min(a)<0.05/(nbitems-k)) {
+difit <- c(difit,which.min(a))
+model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
+a <- sapply(1:nbitems,function(x)
+ifelse(!(x%in%difit),lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2],1))
+}
+}
+return(difit)
+}
+rosali(dat=dat,items=1:7,group="TT")
+pcm(dat,items=1:7,group="TT",dif.items=1:2)
+summary(pcm(dat,items=1:7,group="TT",dif.items=1:2))
diff --git a/RProject/Scripts/rosali.R b/RProject/Scripts/rosali.R
new file mode 100644
index 0000000..0283576
--- /dev/null
+++ b/RProject/Scripts/rosali.R
@@ -0,0 +1,31 @@
+rosali <- function(dat=NULL,items=NULL,group=NULL) {
+  nbitems <- length(items)
+  items2 <- items
+  # Un seul groupe
+  if (length(group)!=1) {
+    stop('Only one variable can be used for the group option')
+  }
+  # Recoder groupe en 0/1
+  dat[,group] <- as.factor(dat[,group])
+  if (!(all(names(levels(dat[,group]))==c("0","1")))) {
+    levels(dat[,group]) <- c("0","1") 
+  }
+  dat[,group] <- as.numeric(dat[,group])-1
+  model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
+  model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
+  ltest <- lmtest::lrtest(model_a,model_b)
+  difit <- c()
+  if (ltest$`Pr(>Chisq)`[2]<0.05) {
+    model_c <- model_b
+    a <- sapply(1:nbitems,function(x) 
+      lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
+    k <- 0
+    while (min(a)<0.05/(nbitems-k)) {
+      difit <- c(difit,which.min(a))
+      model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
+      a <- sapply(1:nbitems,function(x) 
+        ifelse(!(x%in%difit),lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2],1))
+    }
+  }
+  return(difit)
+}
\ No newline at end of file