diff --git a/RProject/Scripts/aggregation.R b/RProject/Scripts/aggregation.R
index f1fb458..5aae908 100644
--- a/RProject/Scripts/aggregation.R
+++ b/RProject/Scripts/aggregation.R
@@ -1190,6 +1190,12 @@ res.dat.dif.resali$theoretical.power <- res.dat$theoretical.power
 #----------------------------------------------------------------------------#
 ##############################################################################
 
+# Items dichotomiques
+
+res.dat$method <- "NONE"
+res.dat.dif$method <- "PERFECT"
+res.dat.dif.rosali$method <- "ROSALI"
+res.dat.dif.resali$method <- "RESIDUS"
 
 
 # Correction of N=50 scenarios
@@ -1201,20 +1207,19 @@ res.dat.dif.resali[res.dat.dif.resali$N==50,]$dif.size <- sapply(which(res.dat.d
 res.dat[res.dat$dif.size!=0 & res.dat$nb.dif==0,]$nb.dif <- 1
 res.dat.dif <- res.dat.dif %>%
   relocate(method, .after = theoretical.power)
-
-#aa <- res.dat.dif[res.dat.dif$scenario=="10B",]
-#bb <- res.dat[res.dat$scenario=="10B",]
-#res.dat.dif[res.dat.dif$scenario=="10B",] <- bb
-#res.dat[res.dat$scenario=="10B",] <- aa
-#res.dat[res.dat$method=="PERFECT",]$method <- "NONE"
-#res.dat.dif[res.dat.dif$method=="NONE",]$method <- "PERFECT"
-
-# Items dichotomiques
-
-res.dat$method <- "NONE"
-res.dat.dif$method <- "PERFECT"
-res.dat.dif.rosali$method <- "ROSALI"
-res.dat.dif.resali$method <- "RESIDUS"
+res.dat[res.dat$scenario=="10B",]$dif.size <- 0.3
+res.dat.dif[res.dat.dif$scenario=="10B",]$dif.size <- 0.3
+res.dat.dif.rosali[res.dat.dif.rosali$scenario=="10B",]$dif.size <- 0.3
+res.dat.dif.resali[res.dat.dif.resali$scenario=="10B",]$dif.size <- 0.3
+res.dat.dif[res.dat.dif$N=="50" & res.dat.dif$nb.dif>0,]$eff.size <- rep(c(0,0.2,0.2,0.4,0.4,-0.2,-0.4),16)
+res.dat.dif[res.dat.dif$N=="50" & res.dat.dif$nb.dif>0 & res.dat.dif$scenario.type=="C",]$bias <- res.dat.dif[res.dat.dif$N=="50" & res.dat.dif$nb.dif>0 & res.dat.dif$scenario.type=="C",]$bias -0.2
+res.dat.dif[res.dat.dif$N=="50" & res.dat.dif$nb.dif>0 & res.dat.dif$scenario.type=="D",]$bias <- res.dat.dif[res.dat.dif$N=="50" & res.dat.dif$nb.dif>0 & res.dat.dif$scenario.type=="D",]$bias +0.6
+res.dat.dif[res.dat.dif$N=="50" & res.dat.dif$nb.dif>0 & res.dat.dif$scenario.type=="E",]$bias <- res.dat.dif[res.dat.dif$N=="50" & res.dat.dif$nb.dif>0 & res.dat.dif$scenario.type=="E",]$bias +0.8
+
+res.dat[res.dat$N=="50" & res.dat$nb.dif>0,]$eff.size <- rep(c(0,0.2,0.2,0.4,0.4,-0.2,-0.4),16)
+res.dat[res.dat$N=="50" & res.dat$nb.dif>0 & res.dat$scenario.type=="C",]$bias <- res.dat[res.dat$N=="50" & res.dat$nb.dif>0 & res.dat$scenario.type=="C",]$bias -0.2
+res.dat[res.dat$N=="50" & res.dat$nb.dif>0 & res.dat$scenario.type=="D",]$bias <- res.dat[res.dat$N=="50" & res.dat$nb.dif>0 & res.dat$scenario.type=="D",]$bias +0.6
+res.dat[res.dat$N=="50" & res.dat$nb.dif>0 & res.dat$scenario.type=="E",]$bias <- res.dat[res.dat$N=="50" & res.dat$nb.dif>0 & res.dat$scenario.type=="E",]$bias +0.8
 
 res.dat.dicho <- res.dat[res.dat$M==2,]
 res.dat.dicho <- rbind(res.dat.dicho,res.dat.dif[res.dat.dif$M==2,])
diff --git a/RProject/Scripts/article.R b/RProject/Scripts/article.R
new file mode 100644
index 0000000..769872f
--- /dev/null
+++ b/RProject/Scripts/article.R
@@ -0,0 +1,224 @@
+
+
+
+
+
+
+
+##########################
+# IGNORING DIF
+##########################
+
+########### Power
+
+# Prepare
+res.dat$dif.agrees.tt <- ifelse(res.dat$eff.size!=0 & res.dat$dif.size!=0, res.dat$dif.size/res.dat$eff.size<0,NA)
+res.dat[res.dat$scenario.type!="A" & res.dat$nb.dif>0 
+        & res.dat$dif.agrees.tt,"dif.power"] <- res.dat[res.dat$scenario.type!="A" 
+                & res.dat$nb.dif>0 & res.dat$dif.agrees.tt,]$h0.rejected.p-res.dat[res.dat$scenario.type!="A" 
+                              & res.dat$nb.dif>0 & res.dat$dif.agrees.tt,]$theoretical.power
+res.dat[res.dat$scenario.type!="A" & res.dat$nb.dif>0 
+        & !res.dat$dif.agrees.tt,"dif.power"] <- res.dat[res.dat$scenario.type!="A" 
+                                                        & res.dat$nb.dif>0 & !res.dat$dif.agrees.tt,]$h0.rejected.p-res.dat[res.dat$scenario.type!="A" 
+                                                                                                                           & res.dat$nb.dif>0 & !res.dat$dif.agrees.tt,]$theoretical.power
+
+
+
+# Histo coloré par typo
+par(mfrow=c(2,1))
+hist(res.dat[abs(res.dat$dif.size)==0.3 & !res.dat$dif.agrees.tt,]$dif.power,breaks = seq(-0.7,0.6,0.05),freq=F,xlim = c(-0.7,0.7),ylim=c(0,4),col=rgb(1,0,0,1/4),
+           main="real power - theoretical power in scenarios with DIF size 0.3",xlab="Real power - theoretical power (raw % difference)")
+hist(res.dat[abs(res.dat$dif.size)==0.3 & res.dat$dif.agrees.tt,]$dif.power,breaks = seq(-0.7,0.6,0.05),freq=F,xlim = c(-0.7,0.7),ylim=c(0,4),col=rgb(0,0,1,1/4),add=T)
+abline(v=0,lty=2,col="black",lwd=2)
+
+hist(res.dat[abs(res.dat$dif.size)==0.5 & !res.dat$dif.agrees.tt,]$dif.power,breaks = seq(-0.7,0.6,0.05),freq=F,xlim = c(-0.7,0.7),ylim=c(0,4),col=rgb(1,0,0,1/4),
+     main="real power - theoretical power in scenarios with DIF size 0.5",xlab="Real power - theoretical power (raw % difference)")
+hist(res.dat[abs(res.dat$dif.size)==0.5 & res.dat$dif.agrees.tt,]$dif.power,breaks = seq(-0.7,0.6,0.05),freq=F,xlim = c(-0.7,0.7),ylim=c(0,4),col=rgb(0,0,1,1/4),add=T)
+abline(v=0,lty=2,col="black",lwd=2)
+
+par(xpd=NA)
+legend(x = -0.825,y=6.25,fill = c(rgb(1,0,0,1/4),rgb(0,0,1,1/4)),c('DIF effect contradicts treatment effect',"DIF effect concurs with treatment effect"),ncol=2)
+
+par(mfrow=c(1,1))
+
+# DIF and treatment opposite signs
+summary(res.dat[res.dat$scenario.type!="A" & res.dat$nb.dif>0 & res.dat$dif.agrees.tt,c("dif.power")])
+
+# DIF and treatment same signs
+summary(res.dat[res.dat$scenario.type!="A" & res.dat$nb.dif>0 & 1-res.dat$dif.agrees.tt,c("dif.power")])
+
+# N=50 vs 300
+summary(res.dat[res.dat$scenario.type!="A" & res.dat$N=="50" & res.dat$nb.dif>0 & res.dat$dif.agrees.tt,c("dif.power")])
+summary(res.dat[res.dat$scenario.type!="A" & res.dat$N==100 & res.dat$nb.dif>0 & res.dat$dif.agrees.tt,c("dif.power")])
+summary(res.dat[res.dat$scenario.type!="A" & res.dat$N==200 & res.dat$nb.dif>0 & res.dat$dif.agrees.tt,c("dif.power")])
+summary(res.dat[res.dat$scenario.type!="A" & res.dat$N==300 & res.dat$nb.dif>0 & res.dat$dif.agrees.tt,c("dif.power")])
+
+########### Treatment effect estimation sign
+
+# Overall
+summary(res.dat[res.dat$scenario.type!="A" & res.dat$nb.dif>0,c("beta.same.sign.truebeta.p")])
+
+# Worst case scenario
+summary(res.dat[res.dat$scenario.type!="A" & res.dat$nb.dif>0 & res.dat$dif.agrees.tt==FALSE & abs(res.dat$dif.size)>0.3 & abs(res.dat$eff.size)==0.2,c("beta.same.sign.truebeta.signif.p")])
+
+########### Bias
+
+summary(res.dat[res.dat$nb.dif>0,c("bias")])
+
+########### true value in CI
+
+summary(abs(res.dat[res.dat$nb.dif>0,c("true.value.in.ci.p")]))
+summary(abs(res.dat[res.dat$N=="50" & res.dat$nb.dif>0,c("true.value.in.ci.p")]))
+
+
+
+
+
+##########################
+# DETECTION
+##########################
+
+# Which performed better
+summary(res.dat.dif.rosali$prop.perfect-res.dat.dif.resali$prop.perfect)
+
+# ROSALI better more than 10% ?
+res.dat.dif.rosali$better <- res.dat.dif.rosali$prop.perfect-res.dat.dif.resali$prop.perfect>0.1
+table(res.dat.dif.rosali[res.dat.dif.rosali$better & res.dat.dif.rosali$nb.dif!=0,"scenario.type"])
+
+# ROSALI worse more than 10% ?
+res.dat.dif.rosali$worse <- res.dat.dif.rosali$prop.perfect-res.dat.dif.resali$prop.perfect< -0.1
+res.dat.dif.rosali[res.dat.dif.rosali$worse & res.dat.dif.rosali$nb.dif!=0,]
+
+# ROSALI perf per subsc
+summary(res.dat.dif.rosali[ res.dat.dif.rosali$N==300 & res.dat.dif.rosali$nb.dif>0 & res.dat.dif.rosali$scenario.type%in%c("C","E"),]$prop.perfect)
+summary(res.dat.dif.rosali[ res.dat.dif.rosali$N==300 & res.dat.dif.rosali$nb.dif>0 & res.dat.dif.rosali$scenario.type%in%c("B","D","F","G"),]$prop.perfect)
+
+# AHRM perf per subsc
+summary(res.dat.dif.resali[ res.dat.dif.resali$N==300 & res.dat.dif.resali$nb.dif>0 & res.dat.dif.resali$scenario.type%in%c("C","E"),]$prop.perfect)
+summary(res.dat.dif.resali[ res.dat.dif.resali$N==300 & res.dat.dif.resali$nb.dif>0 & res.dat.dif.resali$scenario.type%in%c("B","D","F","G"),]$prop.perfect)
+
+# False DIF detect
+summary(res.dat.dif.rosali[res.dat.dif.rosali$nb.dif==0,"dif.detected"])
+summary(res.dat.dif.resali[res.dat.dif.resali$nb.dif==0,"dif.detected"])
+
+# Causal inference NO DIF
+summary(res.dat.dif.resali[res.dat.dif.resali$nb.dif==0,"bias"])
+summary(res.dat.dif.rosali[res.dat.dif.rosali$nb.dif==0,"bias"])
+summary(res.dat.dif.rosali[res.dat.dif.rosali$nb.dif==0 & res.dat.dif.rosali$eff.size==0,"h0.rejected.p"])
+summary(res.dat.dif.resali[res.dat.dif.resali$nb.dif==0 & res.dat.dif.resali$eff.size==0,"h0.rejected.p"])
+summary(res.dat.dif.resali[res.dat.dif.resali$nb.dif==0,"true.value.in.ci.p"])
+summary(res.dat.dif.rosali[res.dat.dif.rosali$nb.dif==0,"true.value.in.ci.p"])
+
+
+##########################
+# TABLES DETECT
+##########################
+res.dat.dif.rosali$dif.agrees.tt <- ifelse(res.dat.dif.rosali$eff.size!=0 & res.dat.dif.rosali$dif.size!=0, res.dat.dif.rosali$dif.size/res.dat.dif.rosali$eff.size<0,NA)
+res.dat.dif.resali$dif.agrees.tt <- ifelse(res.dat.dif.resali$eff.size!=0 & res.dat.dif.resali$dif.size!=0, res.dat.dif.resali$dif.size/res.dat.dif.resali$eff.size<0,NA)
+res.dat.dif.rosali$dif.dir <- sign(res.dat.dif.rosali$dif.size)
+res.dat.dif.resali$dif.dir <- sign(res.dat.dif.resali$dif.size)
+
+tab3.resali <- res.dat.dif.resali[res.dat.dif.resali$dif.size!=0,
+                                  c("scenario","N",
+                                    "dif.detected","prop.perfect","flexible.detect","moreflexible.detect"
+                                    )]
+tab3.resali <-reshape(data = tab3.resali,direction = "wide", idvar = c("scenario"),timevar = "N")
+tab3.rosali <- res.dat.dif.rosali[res.dat.dif.rosali$dif.size!=0,
+                                  c("scenario","N","J","M","eff.size","dif.size","dif.dir","nb.dif",
+                                    "dif.detected","prop.perfect","flexible.detect","moreflexible.detect"
+                                  )]
+tab3.rosali <-reshape(data = tab3.rosali,direction = "wide", idvar = c("scenario","J","M","eff.size","dif.size","dif.dir","nb.dif"),timevar = "N")
+tab3 <- merge(tab3.rosali,tab3.resali,by="scenario",suffixes = c(".rosali",".residuals"))
+tab3 <- rbind(tab3[78:112,],tab3[1:77,])
+tab3$dif.size <- abs(tab3$dif.size)
+write.csv(tab3,"/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/tab3.csv")
+
+
+
+
+##########################
+# TABLES CAUSAL
+##########################
+res.dat.dif.rosali$dif.power <- res.dat.dif.rosali$h0.rejected.p-res.dat.dif.rosali$theoretical.power
+res.dat.dif.resali$dif.power <- res.dat.dif.resali$h0.rejected.p-res.dat.dif.resali$theoretical.power
+
+res.dat.dif.rosali$typeI.error <- ifelse(res.dat.dif.rosali$scenario.type=="A",res.dat.dif.rosali$h0.rejected.p,NA)
+res.dat.dif.rosali$diff.power <- ifelse(res.dat.dif.rosali$scenario.type!="A",res.dat.dif.rosali$dif.power,NA)
+res.dat.dif.resali$typeI.error <- ifelse(res.dat.dif.resali$scenario.type=="A",res.dat.dif.resali$h0.rejected.p,NA)
+res.dat.dif.resali$diff.power <- ifelse(res.dat.dif.resali$scenario.type!="A",res.dat.dif.resali$dif.power,NA)
+
+tab4.resali <- res.dat.dif.resali[res.dat.dif.resali$dif.size!=0,
+                                  c("scenario","N",
+                                    "h0.rejected.p","theoretical.power","true.value.in.ci.p","beta.same.sign.truebeta.p","beta.same.sign.truebeta.signif.p","bias"
+                                  )]
+tab4.resali <-reshape(data = tab4.resali,direction = "wide", idvar = c("scenario"),timevar = "N")
+tab4.rosali <- res.dat.dif.rosali[res.dat.dif.rosali$dif.size!=0,
+                                  c("scenario","N","J","M","eff.size","dif.size","dif.dir","nb.dif",
+                                    "h0.rejected.p","theoretical.power","true.value.in.ci.p","beta.same.sign.truebeta.p","beta.same.sign.truebeta.signif.p","bias"
+                                  )]
+tab4.rosali <-reshape(data = tab4.rosali,direction = "wide", idvar = c("scenario","J","M","eff.size","dif.size","dif.dir","nb.dif"),timevar = "N")
+tab4 <- merge(tab4.rosali,tab4.resali,by="scenario",suffixes = c(".rosali",".residuals"))
+tab4 <- rbind(tab4[78:112,],tab4[1:77,])
+tab4$dif.size <- abs(tab4$dif.size)
+write.csv(tab4,"/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/tab4.csv")
+
+
+##########################
+# TABLES NODIF
+##########################
+tabs2.resali <- res.dat.dif.resali[res.dat.dif.resali$dif.size==0,
+                                  c("scenario","N",
+                                    "dif.detected","h0.rejected.p","theoretical.power","true.value.in.ci.p","beta.same.sign.truebeta.p","beta.same.sign.truebeta.signif.p","bias"
+                                  )]
+tabs2.resali <-reshape(data = tabs2.resali,direction = "wide", idvar = c("scenario"),timevar = "N")
+tabs2.rosali <- res.dat.dif.rosali[res.dat.dif.rosali$dif.size==0,
+                                  c("scenario","N","J","M","eff.size",
+                                    "dif.detected","h0.rejected.p","theoretical.power","true.value.in.ci.p","beta.same.sign.truebeta.p","beta.same.sign.truebeta.signif.p","bias"
+                                  )]
+tabs2.rosali <-reshape(data = tabs2.rosali,direction = "wide", idvar = c("scenario","J","M","eff.size"),timevar = "N")
+tabs2 <- merge(tabs2.rosali,tabs2.resali,by="scenario",suffixes = c(".rosali",".residuals"))
+write.csv(tabs2,"/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/tabs2.csv")
+
+
+
+##########################
+# STATS DIF DETECTION
+##########################
+
+# sample size
+summary(tab3$moreflexible.detect.50.rosali)
+summary(tab3$moreflexible.detect.50.residuals)
+tab3[which.max(tab3$prop.perfect.50.residuals),]
+
+summary(tab3$moreflexible.detect.100.rosali)
+summary(tab3$moreflexible.detect.100.residuals)
+summary(tab3[tab3$nb.dif==1,]$moreflexible.detect.100.residuals)
+summary(tab3[tab3$nb.dif==1,]$moreflexible.detect.100.rosali)
+
+summary(tab3$moreflexible.detect.200.rosali)
+summary(tab3$moreflexible.detect.200.residuals)
+
+summary(tab3$moreflexible.detect.300.rosali)
+summary(tab3$moreflexible.detect.300.residuals)
+summary(tab3[tab3$nb.dif==1,]$moreflexible.detect.300.rosali)
+summary(tab3[tab3$nb.dif==1,]$moreflexible.detect.300.residuals)
+
+summary(tab3[tab3$nb.dif==3 & tab3$J==7,]$flexible.detect.300.rosali)
+summary(tab3[tab3$nb.dif==3 & tab3$J==7,]$flexible.detect.300.residuals)
+
+summary(tab3[tab3$nb.dif==2 & tab3$J==4,]$prop.perfect.300.rosali)
+summary(tab3[tab3$nb.dif==2 & tab3$J==4,]$prop.perfect.300.residuals)
+summary(tab3[tab3$nb.dif==2 & tab3$J==7,]$prop.perfect.300.rosali)
+summary(tab3[tab3$nb.dif==2 & tab3$J==7,]$prop.perfect.300.residuals)
+nrow(tab3[tab3$nb.dif==2 & tab3$prop.perfect.300.residuals>0.5,])
+nrow(tab3[tab3$nb.dif==2 & tab3$prop.perfect.300.rosali>0.5,])
+
+summary(tab3[tab3$M==2,]$prop.perfect.300.rosali)
+summary(tab3[tab3$M==2,]$prop.perfect.300.residuals)
+summary(tab3[tab3$M==4,]$prop.perfect.300.rosali)
+summary(tab3[tab3$M==4,]$prop.perfect.300.residuals)
+
+summary(tab3[tab3$M==2,]$prop.perfect.200.rosali)
+summary(tab3[tab3$M==2,]$prop.perfect.200.residuals)
+summary(tab3[tab3$M==4,]$prop.perfect.200.rosali)
+summary(tab3[tab3$M==4,]$prop.perfect.200.residuals)