From 6086d74551bb039b317fe5e01a65961e2a0bbd66 Mon Sep 17 00:00:00 2001
From: corentinchoisy <corentin.choisy@proton.me>
Date: Wed, 29 May 2024 14:34:03 +0200
Subject: [PATCH] Final code before cleaning
---
RProject/Scripts/aggregation.R | 2 +-
RProject/Scripts/article.R | 283 +++++++++++++++++++++++++++++----
2 files changed, 250 insertions(+), 35 deletions(-)
diff --git a/RProject/Scripts/aggregation.R b/RProject/Scripts/aggregation.R
index 5aae908..2339800 100644
--- a/RProject/Scripts/aggregation.R
+++ b/RProject/Scripts/aggregation.R
@@ -14,7 +14,7 @@ library(dplyr)
library(readxl)
lastChar <- function(str){
- substr(str, nchar(str)-2, nchar(str))
+ substr(str, nchar(str), nchar(str))
}
##############################################################################
diff --git a/RProject/Scripts/article.R b/RProject/Scripts/article.R
index 769872f..b396463 100644
--- a/RProject/Scripts/article.R
+++ b/RProject/Scripts/article.R
@@ -110,31 +110,92 @@ 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
+####################################################
+
+
+##########################
+# TABLES NO DIF RECOVERY
+##########################
+res.dat$dif.dir <- sign(res.dat$dif.size)
+res.dat.dif$dif.dir <- sign(res.dat.dif$dif.size)
+
+tabs1 <- res.dat[res.dat$dif.size==0,
+ c("scenario","N","J","M","eff.size",
+ "h0.rejected.p","theoretical.power","true.value.in.ci.p","beta.same.sign.truebeta.p","beta.same.sign.truebeta.signif.p","bias"
+ )]
+tabs1 <-reshape(data = tabs1,direction = "wide", idvar = c("scenario","J","M","eff.size"),timevar = "N")
+write.csv(tabs1,"/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/tabs1.csv")
+
+
+tabs2 <- res.dat[res.dat$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"
+ )]
+tabs2 <-reshape(data = tabs2,direction = "wide", idvar = c("scenario","J","M","eff.size","dif.size","dif.dir","nb.dif"),timevar = "N")
+tabs2$dif.size <- abs(tabs2$dif.size)
+write.csv(tabs2,"/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/tabs2.csv")
+
+##########################
+# TABLES PERF DIF RECOVERY
+##########################
+tabs3 <- res.dat.dif[res.dat.dif$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"
+ )]
+tabs3 <-reshape(data = tabs3,direction = "wide", idvar = c("scenario","J","M","eff.size","dif.size","dif.dir","nb.dif"),timevar = "N")
+tabs3$dif.size <- abs(tabs3$dif.size)
+write.csv(tabs3,"/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/tabs3.csv")
+
+
+
+
+
+
+
+
##########################
# TABLES DETECT
##########################
+
+# false dif detection
+
+tab3.resali <- res.dat.dif.resali[res.dat.dif.resali$dif.size==0,
+ c("scenario","N",
+ "dif.detected"
+ )]
+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.detected"
+ )]
+tab3.rosali <-reshape(data = tab3.rosali,direction = "wide", idvar = c("scenario","J","M","eff.size"),timevar = "N")
+tab3 <- merge(tab3.rosali,tab3.resali,by="scenario",suffixes = c(".rosali",".residuals"))
+write.csv(tab3,"/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/tab3.csv")
+
+
+# dif detection
+
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,
+tabs4.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,
+tabs4.resali <-reshape(data = tabs4.resali,direction = "wide", idvar = c("scenario"),timevar = "N")
+tabs4.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")
-
-
-
+tabs4.rosali <-reshape(data = tabs4.rosali,direction = "wide", idvar = c("scenario","J","M","eff.size","dif.size","dif.dir","nb.dif"),timevar = "N")
+tabs4 <- merge(tabs4.rosali,tabs4.resali,by="scenario",suffixes = c(".rosali",".residuals"))
+tabs4 <- rbind(tabs4[78:112,],tabs4[1:77,])
+tabs4$dif.size <- abs(tabs4$dif.size)
+write.csv(tabs4,"/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/tabs4.csv")
##########################
# TABLES CAUSAL
@@ -147,37 +208,22 @@ res.dat.dif.rosali$diff.power <- ifelse(res.dat.dif.rosali$scenario.type!="A",re
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,
+tabs5.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,
+tabs5.resali <-reshape(data = tabs5.resali,direction = "wide", idvar = c("scenario"),timevar = "N")
+tabs5.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")
+tabs5.rosali <-reshape(data = tabs5.rosali,direction = "wide", idvar = c("scenario","J","M","eff.size","dif.size","dif.dir","nb.dif"),timevar = "N")
+tabs5 <- merge(tabs5.rosali,tabs5.resali,by="scenario",suffixes = c(".rosali",".residuals"))
+tabs5 <- rbind(tabs5[78:112,],tabs5[1:77,])
+tabs5$dif.size <- abs(tabs5$dif.size)
+write.csv(tabs5,"/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/tabs5.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")
@@ -222,3 +268,172 @@ 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)
+
+summary(tab3[tab3$dif.size==0.3,]$prop.perfect.300.rosali)
+summary(tab3[tab3$dif.size==0.3,]$prop.perfect.300.residuals)
+
+summary(tab3[tab3$dif.size==0.5,]$prop.perfect.300.rosali)
+summary(tab3[tab3$dif.size==0.5,]$prop.perfect.300.residuals)
+
+summary(tab3[tab3$dif.dir==sign(tab3$eff.size),]$prop.perfect.300.rosali)
+summary(tab3[tab3$dif.dir!=sign(tab3$eff.size),]$prop.perfect.300.rosali)
+
+summary(tab3[tab3$dif.dir==sign(tab3$eff.size),]$prop.perfect.300.residuals)
+summary(tab3[tab3$dif.dir==-sign(tab3$eff.size),]$prop.perfect.300.residuals)
+
+summary(tab3$moreflexible.detect.300.rosali-tab3$flexible.detect.300.rosali)
+summary(tab3$moreflexible.detect.300.residuals-tab3$flexible.detect.300.residuals)
+
+
+res.dat[res.dat$N=="300" & res.dat$scenario.type=="A" & abs(res.dat$dif.size)==0.5 &
+ res.dat$nb.dif==2 & res.dat$J==4,]
+summary(tab3[tab3$eff.size==0,]$prop.perfect.300.residuals)
+summary(tab3[tab3$eff.size==0,]$prop.perfect.300.rosali)
+
+
+
+length(res.dat.dif.rosali[res.dat.dif.rosali$nb.dif>0 &
+ res.dat.dif.rosali$prop.perfect>0.3,]$prop.perfect)/nrow(res.dat.dif.rosali[res.dat.dif.rosali$nb.dif>0,])
+
+summary(res.dat.dif.rosali[res.dat.dif.rosali$nb.dif>0 &
+ res.dat.dif.rosali$prop.perfect>0.3,]$prop.perfect)
+
+length(res.dat.dif.resali[res.dat.dif.resali$nb.dif>0 &
+ res.dat.dif.resali$prop.perfect>0.3,]$prop.perfect)/nrow(res.dat.dif.resali[res.dat.dif.resali$nb.dif>0,])
+
+summary(res.dat.dif.resali[res.dat.dif.resali$nb.dif>0 &
+ res.dat.dif.resali$prop.perfect>0.3,]$prop.perfect)
+
+
+
+
+
+
+
+
+##########################
+# PLOTS CAUSAL
+##########################
+
+
+
+###
+
+# Plot bias vs perf detect
+
+plot(seq(0,0.85,0.001),sapply(seq(0,0.85,0.001),
+ function(x) mean(abs(res.dat.dif.rosali[res.dat.dif.rosali$prop.perfect>=x-0.05 & res.dat.dif.rosali$prop.perfect<=x+0.05,]$bias),na.rm = T)),
+ type="l",ylim=c(0,0.2),lwd=2,col="red",xaxs = "i",yaxs="i",xlab = "Perfect detection rate",ylab="Average absolute value bias")
+lines(seq(0,0.85,0.001),sapply(seq(0,0.85,0.001),
+ function(x) mean(abs(res.dat.dif.resali[res.dat.dif.resali$prop.perfect>=x-0.05 & res.dat.dif.resali$prop.perfect<=x+0.05,]$bias),na.rm = T)),
+ type="l",lwd=2,col="blue",lty=2,xlab = "Perfect detection rate",ylab="Average absolute value bias")
+#title("Average absolute value bias in scenarios at given perfect detection rate")
+
+# Plot true bias vs perf detect
+
+lines(seq(0,0.85,0.001),sapply(seq(0,0.85,0.001),
+ function(x) mean(abs(res.dat[res.dat.dif.rosali$prop.perfect>=x-0.05 & res.dat.dif.rosali$prop.perfect<=x+0.05,]$bias),na.rm = T)),
+ type="l",ylim=c(0,0.2),lwd=2,col="pink",xlab = "Perfect detection rate",ylab="Average absolute value bias")
+lines(seq(0,0.85,0.001),sapply(seq(0,0.85,0.001),
+ function(x) mean(abs(res.dat[res.dat.dif.resali$prop.perfect>=x-0.05 & res.dat.dif.resali$prop.perfect<=x+0.05,]$bias),na.rm = T)),
+ type="l",lwd=2,col="#8193f1",lty=2,xlab = "Perfect detection rate",ylab="Average absolute value bias")
+legend(x=0.535,y=0.195,legend=c('ROSALI - accounting for DIF','Residuals - accounting for DIF',
+ 'ROSALI - not accounting for DIF','Residuals - not accounting for DIF'),
+ col=c("red","blue","pink","#8193f1"),lty=c(1,2,1,2),lwd=c(2,2,2,2))
+###
+
+
+###
+
+# Plot alpha vs. perfect
+
+plot(seq(0,0.85,0.001),sapply(seq(0,0.85,0.001),
+ function(x) mean(abs(res.dat.dif.rosali[res.dat.dif.rosali$scenario.type=="A" & res.dat.dif.rosali$prop.perfect>=x-0.1 & res.dat.dif.rosali$prop.perfect<=x+0.1,]$h0.rejected.p),na.rm = T)),
+ type="l",ylim=c(0,1),lwd=2,col="red",xaxs = "i",yaxs="i",xlab = "Perfect detection rate",ylab="Type-I error rate")
+lines(seq(0,0.85,0.001),sapply(seq(0,0.85,0.001),
+ function(x) mean(abs(res.dat.dif.resali[res.dat.dif.resali$scenario.type=="A" & res.dat.dif.resali$prop.perfect>=x-0.1 & res.dat.dif.resali$prop.perfect<=x+0.1,]$h0.rejected.p),na.rm = T)),
+ type="l",lwd=2,col="blue",lty=2,xlab = "Perfect detection rate",ylab="Average absolute value bias")
+lines(seq(0,0.85,0.001),sapply(seq(0,0.85,0.001),
+ function(x) mean(abs(res.dat[res.dat.dif.rosali$scenario.type=="A" & res.dat.dif.rosali$prop.perfect>=x-0.1 & res.dat.dif.rosali$prop.perfect<=x+0.1,]$h0.rejected.p),na.rm = T)),
+ type="l",ylim=c(0,1),lwd=2,col="pink",xlab = "Perfect detection rate",ylab="Type-I error rate")
+lines(seq(0,0.85,0.001),sapply(seq(0,0.85,0.001),
+ function(x) mean(abs(res.dat[res.dat.dif.resali$scenario.type=="A" & res.dat.dif.resali$prop.perfect>=x-0.1 & res.dat.dif.resali$prop.perfect<=x+0.1,]$h0.rejected.p),na.rm = T)),
+ type="l",lwd=2,col="#8193f1",lty=2,xlab = "Perfect detection rate",ylab="Average absolute value bias")
+#title("Average type-I error rate in scenarios at given perfect detection rate")
+legend(x=0.535,y=0.98,legend=c('ROSALI - accounting for DIF','Residuals - accounting for DIF',
+ 'ROSALI - not accounting for DIF','Residuals - not accounting for DIF'),
+ col=c("red","blue","pink","#8193f1"),lty=c(1,2,1,2),lwd=c(2,2,2,2))
+abline(h=0.05,lty=3)
+###
+
+
+
+###
+
+# Plot truevalueinci vs. perfect
+
+plot(seq(0,0.85,0.001),sapply(seq(0,0.85,0.001),
+ function(x) mean(abs(res.dat.dif.rosali[res.dat.dif.rosali$prop.perfect>=x-0.05 & res.dat.dif.rosali$prop.perfect<=x+0.05,]$true.value.in.ci.p),na.rm = T)),
+ type="l",ylim=c(0.5,1),lwd=2,col="red",xaxs = "i",yaxs="i",xlab = "Perfect detection rate",ylab="Average proportion of true effect in estimate CI")
+lines(seq(0,0.85,0.001),sapply(seq(0,0.85,0.001),
+ function(x) mean(abs(res.dat.dif.resali[res.dat.dif.resali$prop.perfect>=x-0.05 & res.dat.dif.resali$prop.perfect<=x+0.05,]$true.value.in.ci.p),na.rm = T)),
+ type="l",lwd=2,col="blue",lty=2,xlab = "Perfect detection rate",ylab="Average absolute value bias")
+lines(seq(0,0.85,0.001),sapply(seq(0,0.85,0.001),
+ function(x) mean(abs(res.dat[res.dat.dif.rosali$prop.perfect>=x-0.05 & res.dat.dif.rosali$prop.perfect<=x+0.05,]$true.value.in.ci.p),na.rm = T)),
+ type="l",ylim=c(0,0.2),lwd=2,col="pink",xlab = "Perfect detection rate",ylab="Average absolute value bias")
+lines(seq(0,0.85,0.001),sapply(seq(0,0.85,0.001),
+ function(x) mean(abs(res.dat[res.dat.dif.resali$prop.perfect>=x-0.05 & res.dat.dif.resali$prop.perfect<=x+0.05,]$true.value.in.ci.p),na.rm = T)),
+ type="l",lwd=2,col="#8193f1",lty=2,xlab = "Perfect detection rate",ylab="Average absolute value bias")
+
+#title("Average proportion of true effect in estimate CI in scenarios at given perfect detection rate")
+legend(x=0.535,y=0.6,legend=c('ROSALI - accounting for DIF','Residuals - accounting for DIF',
+ 'ROSALI - not accounting for DIF','Residuals - not accounting for DIF'),
+ col=c("red","blue","pink","#8193f1"),lty=c(1,2,1,2),lwd=c(2,2,2,2))
+###
+
+
+###
+
+# Plot powerdif vs. perfect
+
+plot(seq(0,0.85,0.001),sapply(seq(0,0.85,0.001),
+ function(x) mean(res.dat.dif.rosali[res.dat.dif.rosali$prop.perfect>=x-0.05 & res.dat.dif.rosali$prop.perfect<=x+0.05,]$diff.power,na.rm = T)),
+ type="l",col="red",xaxs = "i",yaxs="i",lwd=2,ylim=c(-0.5,0.5),xlab = "Perfect detection rate",ylab="Observed power - theoretical power (average)")
+lines(seq(0,0.85,0.001),sapply(seq(0,0.85,0.001),
+ function(x) mean(res.dat.dif.resali[res.dat.dif.resali$prop.perfect>=x-0.05 & res.dat.dif.resali$prop.perfect<=x+0.05,]$diff.power,na.rm = T)),
+ type="l",col="blue",lwd=2,lty=2,ylim=c(-0.2,0),xlab = "Perfect detection rate",ylab="Average absolute value bias")
+lines(seq(0,0.85,0.001),sapply(seq(0,0.85,0.001),
+ function(x) mean(res.dat[res.dat.dif.rosali$prop.perfect>=x-0.05 & res.dat.dif.rosali$prop.perfect<=x+0.05,]$dif.power,na.rm = T)),
+ type="l",col="pink",lwd=2,ylim=c(-0.2,0.1),xlab = "Perfect detection rate",ylab="Observed power - theoretical power (average)")
+lines(seq(0,0.85,0.001),sapply(seq(0,0.85,0.001),
+ function(x) mean(res.dat[res.dat.dif.resali$prop.perfect>=x-0.05 & res.dat.dif.resali$prop.perfect<=x+0.05,]$dif.power,na.rm = T)),
+ type="l",col="#8193f1",lwd=2,lty=2,ylim=c(-0.2,0),xlab = "Perfect detection rate",ylab="Average absolute value bias")
+#title("Average difference with theoretical power in scenarios at given perfect detection rate")
+legend(x=0.54,y=0.48,legend=c('ROSALI - accounting for DIF','Residuals - accounting for DIF',
+ 'ROSALI - not accounting for DIF','Residuals - not accounting for DIF'),
+ col=c("red","blue","pink","#8193f1"),lty=c(1,2,1,2),lwd=c(2,2,2,2))
+###
+
+
+###
+
+# Plot betasame vs. perfect
+
+plot(seq(0,0.85,0.001),sapply(seq(0,0.85,0.001),
+ function(x) mean(abs(res.dat.dif.rosali[res.dat.dif.rosali$prop.perfect>=x-0.05 & res.dat.dif.rosali$prop.perfect<=x+0.05,]$beta.same.sign.truebeta.p),na.rm = T)),
+ type="l",ylim=c(0.5,1),lwd=2,col="red",xaxs = "i",yaxs="i",xlab = "Perfect detection rate",ylab="Average proportion of true effect of same sign as estimate")
+lines(seq(0,0.85,0.001),sapply(seq(0,0.85,0.001),
+ function(x) mean(abs(res.dat.dif.resali[res.dat.dif.resali$prop.perfect>=x-0.05 & res.dat.dif.resali$prop.perfect<=x+0.05,]$beta.same.sign.truebeta.p),na.rm = T)),
+ type="l",lwd=2,col="blue",lty=2,xlab = "Perfect detection rate",ylab="Average absolute value bias")
+lines(seq(0,0.85,0.001),sapply(seq(0,0.85,0.001),
+ function(x) mean(abs(res.dat[res.dat.dif.rosali$prop.perfect>=x-0.05 & res.dat.dif.rosali$prop.perfect<=x+0.05,]$beta.same.sign.truebeta.p),na.rm = T)),
+ type="l",ylim=c(0,0.2),lwd=2,col="pink",xlab = "Perfect detection rate",ylab="Average absolute value bias")
+lines(seq(0,0.85,0.001),sapply(seq(0,0.85,0.001),
+ function(x) mean(abs(res.dat[res.dat.dif.resali$prop.perfect>=x-0.05 & res.dat.dif.resali$prop.perfect<=x+0.05,]$beta.same.sign.truebeta.p),na.rm = T)),
+ type="l",lwd=2,col="#8193f1",lty=2,xlab = "Perfect detection rate",ylab="Average absolute value bias")
+
+#title("Average proportion of true effect in estimate CI in scenarios at given perfect detection rate")
+legend(x=0.535,y=0.6,legend=c('ROSALI - accounting for DIF','Residuals - accounting for DIF',
+ 'ROSALI - not accounting for DIF','Residuals - not accounting for DIF'),
+ col=c("red","blue","pink","#8193f1"),lty=c(1,2,1,2),lwd=c(2,2,2,2))
+###