corentinchoisy/simul_these
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Cleaned up code

Browse Source
This commit is contained in:
Corentin Choisy
2024-05-30 16:50:38 +02:00
parent 6086d74551
commit fd0d2de43f
50 changed files with 370 additions and 8057 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1236
RProject/Scripts/Analysis/aggregation.R Normal file
View File

File diff suppressed because it is too large Load Diff

430
RProject/Scripts/Analysis/article.R Normal file
View File

@ -0,0 +1,430 @@
source(paste0(getwd(),"/functions/resali.R"))
##########################
# 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 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))
###

166
RProject/Scripts/Analysis/functions/resali.R Normal file
View File

@ -0,0 +1,166 @@
library(TAM)
resali <- function(df=NULL,items=NULL,group=NULL,verbose=T) {
if (verbose) {
cat('-----------------------------------------------------------\n')
cat('COMPUTING INITIAL PCM\n')
cat('-----------------------------------------------------------\n')
}
nbitems <- length(items)
nbitems_o <- nbitems
items_n <- paste0('item',items)
resp <- df[,items_n]
grp <- df[,group]
pcm_initial <- TAM::tam.mml(resp=resp,Y=grp,irtmodel = "PCM",est.variance = T,verbose=F)
dat <- df
dat$score <- rowSums(dat[,items_n])
nqt <- ifelse(length(unique(quantile(dat$score,seq(0,1,0.2))))==6,5,length(unique(quantile(dat$score,seq(0,1,0.2))))-1)
while (length(unique(quantile(dat$score,seq(0,1,1/nqt))))!=nqt+1) {
nqt <- nqt-1
}
# ITEM POLYTOMIQUE
if (max(resp)>1) {
dat$score_q5 <- cut(dat$score,unique(quantile(dat$score,seq(0,1,1/nqt))),labels=1:nqt,include.lowest=T)
res.anova <- rep(NA,nbitems)
pval <- rep(NA,nbitems)
fval <- rep(NA,nbitems)
for (i in items) {
dat[,paste0('res_',i)] <- IRT.residuals(pcm_initial)$stand_residuals[,i]
res.anova[i] <- summary(aov(dat[,paste0('res_',i)]~TT*score_q5,data=dat))
pval[i] <- res.anova[[i]][1,"Pr(>F)"]
fval[i] <- res.anova[[i]][1,'F value']
}
if (verbose) {
cat('DONE\n')
cat('-----------------------------------------------------------\n')
}
res.items <- c()
res.uniform <- c()
k <- 1
while(any(pval<0.05/(nbitems_o*3))) {
k <- k+1
if (verbose) {
cat(paste('COMPUTING STEP',k,'\n'))
cat('-----------------------------------------------------------\n')
}
res.item <- gsub("[a-z]", "",colnames(resp)[which.max(fval)])
res.items <- c(res.items,res.item)
res.uni <- res.anova[[which.max(fval)]][3,"Pr(>F)"]>0.05
res.uniform <- c(res.uniform,res.uni)
items_n <- c(items_n[items_n!=paste0('item',res.item)],paste0("item",res.item,c("noTT","TT")))
dat[dat$TT==1,paste0("item",res.item,'TT')] <- dat[dat$TT==1,paste0('item',res.item)]
dat[dat$TT==0,paste0("item",res.item,'noTT')] <- dat[dat$TT==0,paste0('item',res.item)]
resp <- dat[,items_n]
grp <- dat[,group]
pcm_while <- TAM::tam.mml(resp=resp,Y=grp,irtmodel = "PCM",est.variance = T,verbose=F)
nbitems <- length(items_n)
res.anova <- rep(NA,nbitems)
pval <- rep(NA,nbitems)
fval <- rep(NA,nbitems)
for (i in 1:nbitems) {
dat[,paste0('res_',i)] <- IRT.residuals(pcm_while)$stand_residuals[,i]
res.anova[i] <- summary(aov(dat[,paste0('res_',i)]~TT*score_q5,data=dat))
pval[i] <- res.anova[[i]][1,"Pr(>F)"]
fval[i] <- res.anova[[i]][1,'F value']
}
zz <- 0
for (name_i in items_n) {
zz <- zz+1
if (grepl("TT",name_i)) {
pval[zz] <- 1
fval[zz] <- 0
}
}
if (verbose) {
cat('DONE\n')
cat('-----------------------------------------------------------\n')
}
}
if (verbose) {
cat("DETECTED DIF ITEMS\n")
cat('-----------------------------------------------------------\n')
}
if (length(res.items>0)) {
results <- data.frame(dif.items=res.items,
uniform=1*res.uniform)
return(results)
}
else {
if (verbose) {
cat("No DIF was detected\n")
}
return(NULL)
}
# ITEM DICHOTOMIQUE
} else {
res.anova <- rep(NA,nbitems)
pval <- rep(NA,nbitems)
fval <- rep(NA,nbitems)
for (i in items) {
dat[,paste0('res_',i)] <- IRT.residuals(pcm_initial)$stand_residuals[,i]
res.anova[i] <- summary(aov(dat[,paste0('res_',i)]~TT,data=dat))
pval[i] <- res.anova[[i]][1,"Pr(>F)"]
fval[i] <- res.anova[[i]][1,'F value']
}
if (verbose) {
cat('DONE\n')
cat('-----------------------------------------------------------\n')
}
res.items <- c()
k <- 1
while(any(pval<0.05/(nbitems_o))) {
k <- k+1
if (verbose) {
cat(paste('COMPUTING STEP',k,'\n'))
cat('-----------------------------------------------------------\n')
}
res.item <- gsub("[a-z]", "",colnames(resp)[which.max(fval)])
res.items <- c(res.items,res.item)
items_n <- c(items_n[items_n!=paste0('item',res.item)],paste0("item",res.item,c("noTT","TT")))
dat[dat$TT==1,paste0("item",res.item,'TT')] <- dat[dat$TT==1,paste0('item',res.item)]
dat[dat$TT==0,paste0("item",res.item,'noTT')] <- dat[dat$TT==0,paste0('item',res.item)]
resp <- dat[,items_n]
grp <- dat[,group]
pcm_while <- TAM::tam.mml(resp=resp,Y=grp,irtmodel = "PCM",est.variance = T,verbose=F)
nbitems <- length(items_n)
res.anova <- rep(NA,nbitems)
pval <- rep(NA,nbitems)
fval <- rep(NA,nbitems)
for (i in 1:nbitems) {
dat[,paste0('res_',i)] <- IRT.residuals(pcm_while)$stand_residuals[,i]
res.anova[i] <- summary(aov(dat[,paste0('res_',i)]~TT,data=dat))
pval[i] <- res.anova[[i]][1,"Pr(>F)"]
fval[i] <- res.anova[[i]][1,'F value']
}
zz <- 0
for (name_i in items_n) {
zz <- zz+1
if (grepl("TT",name_i)) {
pval[zz] <- 1
fval[zz] <- 0
}
}
if (verbose) {
cat('DONE\n')
cat('-----------------------------------------------------------\n')
}
}
if (verbose) {
cat("DETECTED DIF ITEMS\n")
cat('-----------------------------------------------------------\n')
}
if (length(res.items>0)) {
results <- data.frame(dif.items=res.items,
uniform=rep(1,length(res.items)))
return(results)
}
else {
if (verbose) {
cat("No DIF was detected\n")
}
return(NULL)
}
}
}

2381
RProject/Scripts/Analysis/pcm_nodif.R Normal file
View File

File diff suppressed because it is too large Load Diff

281
RProject/Scripts/Analysis/resali_analysis.R Normal file
View File

@ -0,0 +1,281 @@
##############################################################################
#----------------------------------------------------------------------------#
################################### RESALI ###################################
#----------------------------------------------------------------------------#
##############################################################################
source(paste0(getwd(),"/functions/resali.R"))
generate_resali <- function(scenario=NULL,grp=NULL) {
scen <- as.numeric(gsub("[A,B,C,D,E,F,G,_]","",substr(scenario,0,3)))
if (substr(scenario,start=nchar(scenario)-1,stop=nchar(scenario))=="50") {
N <- 50
}
if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="100") {
N <- 100
}
if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="200") {
N <- 200
}
if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="300") {
N <- 300
}
if (scen<5) {
dat <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Data/NoDIF/N',N,'/scenario_',scenario,'.csv'))
}
if (scen>=5) {
dat <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N',N,'/scenario_',scenario,'.csv'))
}
if (scen%in%c(3,4,13:20)) {
res <- resali(df=dat[dat$replication==1,],items = seq(1,7),group=grp,verbose=FALSE)
df_res <- data.frame(dif.detect.1=ifelse(length(res$dif.items)>=1,res$dif.items[1],NA),
dif.detect.2=ifelse(length(res$dif.items)>=2,res$dif.items[2],NA),
dif.detect.3=ifelse(length(res$dif.items)>=3,res$dif.items[3],NA),
dif.detect.4=ifelse(length(res$dif.items)>=4,res$dif.items[4],NA),
dif.detect.5=ifelse(length(res$dif.items)>=5,res$dif.items[5],NA),
dif.detect.6=ifelse(length(res$dif.items)>=6,res$dif.items[6],NA),
dif.detect.7=ifelse(length(res$dif.items)>=7,res$dif.items[7],NA),
dif.detect.unif.1=ifelse(length(res$uniform)>=1,res$uniform[1],NA),
dif.detect.unif.2=ifelse(length(res$uniform)>=2,res$uniform[2],NA),
dif.detect.unif.3=ifelse(length(res$uniform)>=3,res$uniform[3],NA),
dif.detect.unif.4=ifelse(length(res$uniform)>=4,res$uniform[4],NA),
dif.detect.unif.5=ifelse(length(res$uniform)>=5,res$uniform[5],NA),
dif.detect.unif.6=ifelse(length(res$uniform)>=6,res$uniform[6],NA),
dif.detect.unif.7=ifelse(length(res$uniform)>=7,res$uniform[7],NA),
N=N,
nbdif=ifelse(scen<=4,0,ifelse(scen<=16,2,3)),
true.dif.1=ifelse(scen<=4,NA,unique(dat[dat$replication==1,]$dif1)),
true.dif.2=ifelse(scen<=4,NA,unique(dat[dat$replication==1,]$dif2)),
true.dif.3=ifelse(scen<=16,NA,unique(dat[dat$replication==1,]$dif3))
)
for (k in 2:1000) {
if (k%%100==0) {
cat(paste0('N=',k,'/1000\n'))
}
res <- resali(df=dat[dat$replication==k,],items = seq(1,7),group=grp,verbose=FALSE)
df_res2 <- data.frame(dif.detect.1=ifelse(length(res$dif.items)>=1,res$dif.items[1],NA),
dif.detect.2=ifelse(length(res$dif.items)>=2,res$dif.items[2],NA),
dif.detect.3=ifelse(length(res$dif.items)>=3,res$dif.items[3],NA),
dif.detect.4=ifelse(length(res$dif.items)>=4,res$dif.items[4],NA),
dif.detect.5=ifelse(length(res$dif.items)>=5,res$dif.items[5],NA),
dif.detect.6=ifelse(length(res$dif.items)>=6,res$dif.items[6],NA),
dif.detect.7=ifelse(length(res$dif.items)>=7,res$dif.items[7],NA),
dif.detect.unif.1=ifelse(length(res$uniform)>=1,res$uniform[1],NA),
dif.detect.unif.2=ifelse(length(res$uniform)>=2,res$uniform[2],NA),
dif.detect.unif.3=ifelse(length(res$uniform)>=3,res$uniform[3],NA),
dif.detect.unif.4=ifelse(length(res$uniform)>=4,res$uniform[4],NA),
dif.detect.unif.5=ifelse(length(res$uniform)>=5,res$uniform[5],NA),
dif.detect.unif.6=ifelse(length(res$uniform)>=6,res$uniform[6],NA),
dif.detect.unif.7=ifelse(length(res$uniform)>=7,res$uniform[7],NA),
N=N,
nbdif=ifelse(scen<=4,0,ifelse(scen<=16,2,3)),
true.dif.1=ifelse(scen<=4,NA,unique(dat[dat$replication==k,]$dif1)),
true.dif.2=ifelse(scen<=4,NA,unique(dat[dat$replication==k,]$dif2)),
true.dif.3=ifelse(scen<=16,NA,unique(dat[dat$replication==k,]$dif3)))
df_res <- rbind(df_res,df_res2)
}
}
else if (scen%in%c(1,2,5:12)) {
res <- resali(df=dat[dat$replication==1,],items = seq(1,4),group=grp,verbose=FALSE)
df_res <- data.frame(dif.detect.1=ifelse(length(res$dif.items)>=1,res$dif.items[1],NA),
dif.detect.2=ifelse(length(res$dif.items)>=2,res$dif.items[2],NA),
dif.detect.3=ifelse(length(res$dif.items)>=3,res$dif.items[3],NA),
dif.detect.4=ifelse(length(res$dif.items)>=4,res$dif.items[4],NA),
dif.detect.unif.1=ifelse(length(res$uniform)>=1,res$uniform[1],NA),
dif.detect.unif.2=ifelse(length(res$uniform)>=2,res$uniform[2],NA),
dif.detect.unif.3=ifelse(length(res$uniform)>=3,res$uniform[3],NA),
dif.detect.unif.4=ifelse(length(res$uniform)>=4,res$uniform[4],NA),
N=N,
nbdif=ifelse(scen<=4,0,ifelse(scen<=8,1,2)),
true.dif.1=ifelse(scen<=4,NA,unique(dat[dat$replication==1,]$dif1)),
true.dif.2=ifelse(scen<=8,NA,unique(dat[dat$replication==1,]$dif2))
)
for (k in 2:1000) {
if (k%%100==0) {
cat(paste0('N=',k,'/1000\n'))
}
res <- resali(df=dat[dat$replication==k,],items = seq(1,4),group=grp,verbose=FALSE)
df_res2 <- data.frame(dif.detect.1=ifelse(length(res$dif.items)>=1,res$dif.items[1],NA),
dif.detect.2=ifelse(length(res$dif.items)>=2,res$dif.items[2],NA),
dif.detect.3=ifelse(length(res$dif.items)>=3,res$dif.items[3],NA),
dif.detect.4=ifelse(length(res$dif.items)>=4,res$dif.items[4],NA),
dif.detect.unif.1=ifelse(length(res$uniform)>=1,res$uniform[1],NA),
dif.detect.unif.2=ifelse(length(res$uniform)>=2,res$uniform[2],NA),
dif.detect.unif.3=ifelse(length(res$uniform)>=3,res$uniform[3],NA),
dif.detect.unif.4=ifelse(length(res$uniform)>=4,res$uniform[4],NA),
N=N,
nbdif=ifelse(scen<=4,0,ifelse(scen<=8,1,2)),
true.dif.1=ifelse(scen<=4,NA,unique(dat[dat$replication==k,]$dif1)),
true.dif.2=ifelse(scen<=8,NA,unique(dat[dat$replication==k,]$dif2)))
df_res <- rbind(df_res,df_res2)
}
}
return(df_res)
}
results <- c(sapply(1:4,function(x) paste0(x,c('A','B','C','D','E'))),sapply(5:9,function(x) paste0(x,c('A','B','C','D','E','F','G'))))
results2 <- c(sapply(10:20,function(x) paste0(x,c('A','B','C','D','E','F','G'))))
results <- c(sapply(c(50,100,200,300),function(x) paste0(results,'_',x)))
results2 <- c(sapply(c(50,100,200,300),function(x) paste0(results2,'_',x)))
results <- sort(results)
results2 <- sort(results2)
results <- c(results,results2)
for (r in results) {
cat(paste0(r,"\n"))
cat(paste0("-------------------------------------------","\n"))
write.csv(generate_resali(r,"TT"),paste0("/home/corentin/Documents/These/Recherche/Simulations/Analysis/RESALI/Detection/",r,".csv"))
cat(paste0("-------------------------------------------","\n"))
}
##############################################################################
#----------------------------------------------------------------------------#
################################### NEWDATA ##################################
#----------------------------------------------------------------------------#
##############################################################################
## Liste des scenarios
results <- c(sapply(1:4,function(x) paste0(x,c('A','B','C','D','E'))),sapply(5:9,function(x) paste0(x,c('A','B','C','D','E','F','G'))))
results2 <- c(sapply(10:20,function(x) paste0(x,c('A','B','C','D','E','F','G'))))
results <- c(sapply(c(50,100,200,300),function(x) paste0(results,'_',x)))
results2 <- c(sapply(c(50,100,200,300),function(x) paste0(results2,'_',x)))
results <- sort(results)
results2 <- sort(results2)
results <- c(results,results2)
## Importer l'analyse resali pour chaque scenario
for (r in results) {
cat('--------------------------------------------------------------------------\n')
cat(paste0(r,"\n"))
cat('--------------------------------------------------------------------------\n')
#### Importer les datas
scen <- as.numeric(gsub("[A,B,C,D,E,F,G,_]","",substr(r,0,3)))
if (substr(r,start=nchar(r)-1,stop=nchar(r))=="50") {
N <- 50
}
if (substr(r,start=nchar(r)-2,stop=nchar(r))=="100") {
N <- 100
}
if (substr(r,start=nchar(r)-2,stop=nchar(r))=="200") {
N <- 200
}
if (substr(r,start=nchar(r)-2,stop=nchar(r))=="300") {
N <- 300
}
if (scen<5) {
datt <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Data/NoDIF/N',N,'/scenario_',r,'.csv'))
}
if (scen>=5) {
datt <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N',N,'/scenario_',r,'.csv'))
}
#### Importer l'analyse
analyse <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/RESALI/Detection/',r,".csv"))
#### Pour chaque replication
for (k in 1:1000) {
if (k%%100==0) {
cat(paste0("N = ",k," / 1000\n"))
}
datt[datt$replication==k,"dif.detect.1"] <- analyse[analyse$X==k,"dif.detect.1"]
datt[datt$replication==k,"dif.detect.2"] <- analyse[analyse$X==k,"dif.detect.2"]
datt[datt$replication==k,"dif.detect.3"] <- analyse[analyse$X==k,"dif.detect.3"]
datt[datt$replication==k,"dif.detect.4"] <- analyse[analyse$X==k,"dif.detect.4"]
datt[datt$replication==k,"dif.detect.unif.1"] <- analyse[analyse$X==k,"dif.detect.unif.1"]
datt[datt$replication==k,"dif.detect.unif.2"] <- analyse[analyse$X==k,"dif.detect.unif.2"]
datt[datt$replication==k,"dif.detect.unif.3"] <- analyse[analyse$X==k,"dif.detect.unif.3"]
datt[datt$replication==k,"dif.detect.unif.4"] <- analyse[analyse$X==k,"dif.detect.unif.4"]
if (scen==3 | scen==4 | scen>=13) {
datt[datt$replication==k,"dif.detect.5"] <- analyse[analyse$X==k,"dif.detect.5"]
datt[datt$replication==k,"dif.detect.6"] <- analyse[analyse$X==k,"dif.detect.6"]
datt[datt$replication==k,"dif.detect.7"] <- analyse[analyse$X==k,"dif.detect.7"]
datt[datt$replication==k,"dif.detect.unif.5"] <- analyse[analyse$X==k,"dif.detect.unif.5"]
datt[datt$replication==k,"dif.detect.unif.6"] <- analyse[analyse$X==k,"dif.detect.unif.6"]
datt[datt$replication==k,"dif.detect.unif.7"] <- analyse[analyse$X==k,"dif.detect.unif.7"]
}
}
datt[is.na(datt$dif.detect.1),"dif.detect.1"] <- ""
datt[is.na(datt$dif.detect.2),"dif.detect.2"] <- ""
datt[is.na(datt$dif.detect.3),"dif.detect.3"] <- ""
datt[is.na(datt$dif.detect.4),"dif.detect.4"] <- ""
datt[is.na(datt$dif.detect.unif.1),"dif.detect.unif.1"] <- ""
datt[is.na(datt$dif.detect.unif.2),"dif.detect.unif.2"] <- ""
datt[is.na(datt$dif.detect.unif.3),"dif.detect.unif.3"] <- ""
datt[is.na(datt$dif.detect.unif.4),"dif.detect.unif.4"] <- ""
if (scen==3 | scen==4 | scen>=13) {
datt[is.na(datt$dif.detect.5),"dif.detect.5"] <- ""
datt[is.na(datt$dif.detect.6),"dif.detect.6"] <- ""
datt[is.na(datt$dif.detect.7),"dif.detect.7"] <- ""
datt[is.na(datt$dif.detect.unif.5),"dif.detect.unif.5"] <- ""
datt[is.na(datt$dif.detect.unif.6),"dif.detect.unif.6"] <- ""
datt[is.na(datt$dif.detect.unif.7),"dif.detect.unif.7"] <- ""
}
write.csv(datt,paste0("/home/corentin/Documents/These/Recherche/Simulations/Analysis/RESALI/Detection_data/",r,".csv"))
}
for (r in results) {
cat('--------------------------------------------------------------------------\n')
cat(paste0(r,"\n"))
cat('--------------------------------------------------------------------------\n')
#### Importer les datas
scen <- as.numeric(gsub("[A,B,C,D,E,F,G,_]","",substr(r,0,3)))
if (substr(r,start=nchar(r)-1,stop=nchar(r))=="50") {
N <- 50
}
if (substr(r,start=nchar(r)-2,stop=nchar(r))=="100") {
N <- 100
}
if (substr(r,start=nchar(r)-2,stop=nchar(r))=="200") {
N <- 200
}
if (substr(r,start=nchar(r)-2,stop=nchar(r))=="300") {
N <- 300
}
#### Importer l'analyse
analyse <- read.csv(paste0("/home/corentin/Documents/These/Recherche/Simulations/Analysis/RESALI/Detection_data/",r,".csv"))
analyse[is.na(analyse)] <- ""
names(analyse)[names(analyse)=="dif.detect.1"] <- "dif_detect_1"
names(analyse)[names(analyse)=="dif.detect.2"] <- "dif_detect_2"
names(analyse)[names(analyse)=="dif.detect.3"] <- "dif_detect_3"
names(analyse)[names(analyse)=="dif.detect.4"] <- "dif_detect_4"
names(analyse)[names(analyse)=="dif.detect.unif.1"] <- "dif_detect_unif_1"
names(analyse)[names(analyse)=="dif.detect.unif.2"] <- "dif_detect_unif_2"
names(analyse)[names(analyse)=="dif.detect.unif.3"] <- "dif_detect_unif_3"
names(analyse)[names(analyse)=="dif.detect.unif.4"] <- "dif_detect_unif_4"
if (scen==3 | scen==4 | scen>=13) {
names(analyse)[names(analyse)=="dif.detect.5"] <- "dif_detect_5"
names(analyse)[names(analyse)=="dif.detect.6"] <- "dif_detect_6"
names(analyse)[names(analyse)=="dif.detect.7"] <- "dif_detect_7"
names(analyse)[names(analyse)=="dif.detect.unif.5"] <- "dif_detect_unif_5"
names(analyse)[names(analyse)=="dif.detect.unif.6"] <- "dif_detect_unif_6"
names(analyse)[names(analyse)=="dif.detect.unif.7"] <- "dif_detect_unif_7"
}
analyse <- analyse[,!names(analyse) %in% c("X","X.1","X.2")]
write.csv(analyse,paste0("/home/corentin/Documents/These/Recherche/Simulations/Analysis/RESALI/Detection_data/",r,".csv"))
}