########################## # 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 #################################################### ########################## # 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) tabs4.resali <- res.dat.dif.resali[res.dat.dif.resali$dif.size!=0, c("scenario","N", "dif.detected","prop.perfect","flexible.detect","moreflexible.detect" )] 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" )] 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 ########################## 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) 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" )] 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" )] 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") ########################## # 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) 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)) ###