diff --git a/RProject/Scripts/Analysis/aggregation.R b/RProject/Scripts/Analysis/aggregation.R index 52688e4..b808859 100644 --- a/RProject/Scripts/Analysis/aggregation.R +++ b/RProject/Scripts/Analysis/aggregation.R @@ -215,6 +215,41 @@ compile_simulation <- function(scenario) { } N <- ifelse(substr(scenario,start=nchar(scenario)-1,stop=nchar(scenario))=="50","50",substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))) zz <- ifelse(N=="50",substr(scenario,start=0,stop=nchar(scenario)-3),substr(scenario,start=0,stop=nchar(scenario)-4)) + zz.type <- substr(zz,start=nchar(zz),stop=nchar(zz)) + + if (unique(is.na(s$dif.size))) { + s$dif.size <- 0 + } + if(zz=="10B") { + s$dif.size <- 0.3 + } + if (substr(zz,1,1)%in%c("6","8")) { + s$nb.dif <- 1 + } + if (substr(zz,1,2)%in%c("10","12","14",'16')) { + s$nb.dif <- 2 + } + if (substr(zz,1,1)%in%c("18","20")) { + s$nb.dif <- 3 + } + + if (zz.type!="A") { + if (zz.type=="B") { + s$eff.size <- 0.2 + } else if (zz.type=="C" & unique(s$dif.size)==0) { + s$eff.size <- 0.4 + } else if (zz.type=="C" & unique(s$dif.size)!=0) { + s$eff.size <- 0.2 + } else if (zz.type=="D" & unique(s$dif.size)!=0) { + s$eff.size <- 0.4 + } else if (zz.type=="E" & unique(s$dif.size)!=0) { + s$eff.size <- 0.4 + } + } else { + s$eff.size <- 0 + s$dif.size <- -1*s$dif.size + } + b <- data.frame(scenario=zz, scenario.type=substr(zz,start=nchar(zz),stop=nchar(zz)), N=N, @@ -229,6 +264,7 @@ compile_simulation <- function(scenario) { se.analytical.beta=mean(s$se.beta), m.low.ci.beta=mean(s$low.ci.beta), m.high.ci.beta=mean(s$high.ci.beta), + bias=mean(s$beta-s$eff.size), true.value.in.ci.p=mean(s$true.value.in.ci), h0.rejected.p=mean(s$h0.rejected), beta.same.sign.truebeta.p=mean(s$beta.same.sign.truebeta,na.rm=T), @@ -247,18 +283,8 @@ for (x in results[seq(2,length(results))]) { res.dat <- bind_rows(res.dat,y) } -res.dat[res.dat$scenario.type=='A','dif.size'] <- -res.dat[res.dat$scenario.type=='A','dif.size'] -res.dat[is.na(res.dat$dif.size),'dif.size'] <- 0 -res.dat[res.dat$scenario=="10B",]$dif.size <- 0.3 -res.dat[substr(res.dat$scenario,1,1)%in%c("6","8"),'nb.dif'] <- 1 -res.dat[substr(res.dat$scenario,1,2)%in%seq(10,16,2),'nb.dif'] <- 2 -res.dat[substr(res.dat$scenario,1,2)%in%seq(18,20,2),'nb.dif'] <- 3 res.dat[res.dat$N==50,"dif.size"] <- res.dat[which(res.dat$N==50)+1,"dif.size"] -res.dat[res.dat$scenario.type=="B",]$eff.size <- 0.2 -res.dat[res.dat$scenario.type=="C" & res.dat$dif.size==0,]$eff.size <- 0.4 -res.dat[res.dat$scenario.type=="C" & res.dat$dif.size!=0,]$eff.size <- 0.2 -res.dat[res.dat$scenario.type=="D" & res.dat$dif.size!=0,]$eff.size <- 0.4 -res.dat[res.dat$scenario.type=="E" & res.dat$dif.size!=0,]$eff.size <- 0.4 + is.nan.data.frame <- function(x) { do.call(cbind, lapply(x, is.nan)) @@ -266,7 +292,6 @@ is.nan.data.frame <- function(x) { res.dat[is.nan(res.dat)] <- NA -res.dat$bias <- res.dat$eff.size-res.dat$m.beta ############################################################################## diff --git a/RProject/Scripts/Analysis/article.R b/RProject/Scripts/Analysis/article.R index cb82c1c..194e623 100644 --- a/RProject/Scripts/Analysis/article.R +++ b/RProject/Scripts/Analysis/article.R @@ -6,129 +6,6 @@ boxplot(as.numeric(res.dat.article$typeIerror)~as.numeric(res.dat.article$N), ylim=c(0,1),xlab="N",ylab="Type-I error",col="#ff7777",pch=3) abline(h=0.05,col="red",lty=2,lwd=2) -########################## -# 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"]) - - - - -########################## -# 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) - - - - - - - - - ########################## # ICC / CCC BASE ########################## @@ -148,7 +25,7 @@ plot.tam.2 <- function(x, items=1:x$nitems, type="expected", # device.Option <- getOption("device") time1 <- NULL - pall <- c('#200c23', '#62403d', '#a87b5e', '#e9bf98' + pall <- c('#f7c611', '#61b5b7', '#ca346d', '#204776' ) if ( fix.devices ){ old.opt.dev <- getOption("device") @@ -335,11 +212,11 @@ plot.tam.2 <- function(x, items=1:x$nitems, type="expected", dfr <- dat2 dfr1a <- dfr[ dfr$cat==kk, ] graphics::plot( dfr1a$Theta, dfr1a$P, ylim=c(-.1,1.1), - xlab=expression(theta), + xlab="", col=pall[kk], type="l", xpd=TRUE,axes=F, ... ) axis(1) - axis(2) + axis(2,las=2) grid(nx = NA, ny = NULL, lty = 3, col = "lightgray", lwd = 1) @@ -420,35 +297,35 @@ segments(x0=zaza$xsi[7,]$xsi,x1=zaza$xsi[7,]$xsi,y0=0,y1=1.1,lty=3) segments(x0=zaza$xsi[8,]$xsi,x1=zaza$xsi[8,]$xsi,y0=0,y1=1.1,lty=3) segments(x0=zaza$xsi[9,]$xsi,x1=zaza$xsi[9,]$xsi,y0=0,y1=1.1,lty=3) -text(x=-2.5,y=0.85,"0",col="#200c23") -text(x=-0.65,y=0.55,"1",col="#62403d") -text(x=0.95,y=0.55,"2",col="#a87b5e") -text(x=2.5,y=0.85,"3",col="#e9bf98") +text(x=-2.5,y=0.85,"0",col="#f7c611") +text(x=-0.65,y=0.55,"1",col="#61b5b7") +text(x=0.95,y=0.55,"2",col="#ca346d") +text(x=2.5,y=0.85,"3",col="#204776") par(xpd=T,mar=c(5.1,4.1,4.1,2.1)) -segments(x0=-3,x1=zaza$xsi[7,]$xsi,y0=1.1,y1=1.1,col="#200c23",lwd=2) -segments(x0=zaza$xsi[7,]$xsi,x1=zaza$xsi[8,]$xsi,y0=1.1,y1=1.1,col="#62403d",lwd=2) -segments(x0=zaza$xsi[8,]$xsi,x1=zaza$xsi[9,]$xsi,y0=1.1,y1=1.1,col="#a87b5e",lwd=2) -segments(x0=zaza$xsi[9,]$xsi,x1=3,y0=1.1,y1=1.1,col="#e9bf98",lwd=2) +segments(x0=-3,x1=zaza$xsi[7,]$xsi,y0=1.1,y1=1.1,col="#f7c611",lwd=2) +segments(x0=zaza$xsi[7,]$xsi,x1=zaza$xsi[8,]$xsi,y0=1.1,y1=1.1,col="#61b5b7",lwd=2) +segments(x0=zaza$xsi[8,]$xsi,x1=zaza$xsi[9,]$xsi,y0=1.1,y1=1.1,col="#ca346d",lwd=2) +segments(x0=zaza$xsi[9,]$xsi,x1=3,y0=1.1,y1=1.1,col="#204776",lwd=2) points(x = zaza$xsi[7,]$xsi, y=1.1,pch=9,cex=1) points(x = zaza$xsi[8,]$xsi, y=1.1,pch=9,cex=1) points(x = zaza$xsi[9,]$xsi, y=1.1,pch=9,cex=1) -text( x=mean(c(-3,zaza$xsi[7,]$xsi)), y=1.2,"0",col="#200c23" ) -text( x=mean(c(zaza$xsi[7,]$xsi,zaza$xsi[8,]$xsi)), y=1.2,"1",col="#62403d" ) -text( x=mean(c(zaza$xsi[8,]$xsi,zaza$xsi[9,]$xsi)), y=1.2,"2",col="#a87b5e" ) -text( x=mean(c(zaza$xsi[9,]$xsi,3)), y=1.2,"3",col="#e9bf98" ) -text( x=mean(c(-3,zaza$xsi[7,]$xsi)),cex=0.7, y=1.15,"Much less than usual",col="#200c23" ) -text( x=mean(c(zaza$xsi[7,]$xsi,zaza$xsi[8,]$xsi)),cex=0.7, y=1.15,"Less so than usual",col="#62403d" ) -text( x=mean(c(zaza$xsi[8,]$xsi,zaza$xsi[9,]$xsi)),cex=0.7, y=1.15,"As much as usual",col="#a87b5e" ) -text( x=mean(c(zaza$xsi[9,]$xsi,3)),cex=0.7, y=1.15,"More so than usual",col="#e9bf98" ) +text( x=mean(c(-3,zaza$xsi[7,]$xsi)), y=1.2,"0",col="#f7c611" ) +text( x=mean(c(zaza$xsi[7,]$xsi,zaza$xsi[8,]$xsi)), y=1.2,"1",col="#61b5b7" ) +text( x=mean(c(zaza$xsi[8,]$xsi,zaza$xsi[9,]$xsi)), y=1.2,"2",col="#ca346d" ) +text( x=mean(c(zaza$xsi[9,]$xsi,3)), y=1.2,"3",col="#204776" ) +text( x=mean(c(-3,zaza$xsi[7,]$xsi)),cex=0.7, y=1.15,"Much less than usual",col="#f7c611" ) +text( x=mean(c(zaza$xsi[7,]$xsi,zaza$xsi[8,]$xsi)),cex=0.7, y=1.15,"Less so than usual",col="#61b5b7" ) +text( x=mean(c(zaza$xsi[8,]$xsi,zaza$xsi[9,]$xsi)),cex=0.7, y=1.15,"As much as usual",col="#ca346d" ) +text( x=mean(c(zaza$xsi[9,]$xsi,3)),cex=0.7, y=1.15,"More so than usual",col="#204776" ) text(x=zaza$xsi[7,]$xsi,y=1.15,expression(delta["j,1"])) text(x=zaza$xsi[8,]$xsi,y=1.15,expression(delta["j,2"])) text(x=zaza$xsi[9,]$xsi,y=1.15,expression(delta["j,3"])) -text(x = 0,y=1.3,"Most probable response catgory") +text(x = 0,y=1.3,"Most probable response category") CurlyBraces(x0=-2.5, x1=2.5, y0=1.225, y1=1.225, pos = 2, direction = 1, depth=0.05) arrows(x0=-2.65,x1=-3,y0=-0.35,length = 0.15,lwd = 2) arrows(x0=2.65,x1=3,y0=-0.35,length = 0.15,lwd = 2) @@ -459,20 +336,27 @@ rect(xleft = 3,xright=5,ybottom = 0,ytop=1.1,col = "white",border = "white") lines(x=c(zaza$xsi[7,]$xsi,zaza$xsi[7,]$xsi),y=c(0,-0.15),lty=3) lines(x=c(zaza$xsi[8,]$xsi,zaza$xsi[8,]$xsi),y=c(0,-0.15),lty=3) lines(x=c(zaza$xsi[9,]$xsi,zaza$xsi[9,]$xsi),y=c(0,-0.15),lty=3) -title(main='Example item: \n "Have you been able to enjoy your normal daily activities?" ', font.main=2) +title(main='Item: \n "Have you been able to enjoy your normal daily activities?" ', font.main=2) +title(xlab=expression("Latent variable " * theta * "\n (mental health)")) par(xpd=F) dev.off() -########################## -# ICC / CCC DIF -########################## -plot.tam.dif <- function(x, items=1:x$nitems, type="expected", - low=-3, high=3, ngroups=6, groups_by_item=FALSE, - wle=NULL, export=TRUE, export.type="png", - export.args=list(), observed=TRUE, overlay=FALSE, - ask=FALSE, package="lattice", - fix.devices=TRUE, nnodes=100, ...) + +#### ICC BASE + +zaz <- read.csv("/home/corentin/Documents/These/Recherche/Simulations/Data/NoDIF/N300/scenario_2A_300.csv") +zaza <- tam.mml(resp = zaz[zaz$TT==0,paste0("item",1:4)]) +zaza2 <- tam.mml(resp = zaz[zaz$TT==1,paste0("item",1:4)]) +zaza3 <- tam.mml(resp = zaz[,paste0("item",1:4)]) + + +icc.tam.base <- function(x, items=1:x$nitems, type="expected",TT=F, + low=-3, high=3, ngroups=6, groups_by_item=FALSE, + wle=NULL, export=TRUE, export.type="png", + export.args=list(), observed=TRUE, overlay=FALSE, + ask=FALSE, package="lattice", + fix.devices=TRUE, nnodes=100, ...) { require_namespace_msg("grDevices") if ( package=="lattice"){ @@ -480,11 +364,10 @@ plot.tam.dif <- function(x, items=1:x$nitems, type="expected", } # device.Option <- getOption("device") - pall <- c('#200c23', '#62403d', '#a87b5e', '#e9bf98' - ) - low <- low-2 - high <- high+2 + time1 <- NULL + pall <- c('#204776', '#204776', '#204776', '#204776' + ) if ( fix.devices ){ old.opt.dev <- getOption("device") old.opt.err <- c( getOption("show.error.messages")) @@ -518,7 +401,7 @@ plot.tam.dif <- function(x, items=1:x$nitems, type="expected", if (ndim==1 ){ theta <- matrix(seq(low, high, length=nnodes), nrow=nnodes, ncol=ndim) } else { - nodes <- seq(low, high, length=nnodes) + nodes <- seq(-3, 3, length=nnodes) theta <- as.matrix( expand.grid( as.data.frame( matrix( rep(nodes, ndim), ncol=ndim ) ) ) ) nnodes <- nrow(theta) B <- tamobj$B @@ -609,9 +492,19 @@ plot.tam.dif <- function(x, items=1:x$nitems, type="expected", if ( type=="expected"){ if (i==1 || !overlay) { ylim2 <- c(0,max( tamobj$resp[,i], na.rm=TRUE ) ) - graphics::plot(theta, expScore[,i],,col=12, type="l", lwd=3, las=1, ylab="Score", xlab="Ability", - main=paste("Expected Scores Curve - Item ", colnames(tamobj$resp)[i] ) , - ylim=ylim2, ... ) + if (TT==0) { + graphics::plot(theta, expScore[,i], type="l", lwd=3, las=1, ylab="Expected item response", + main=NULL,xlab="", + col="#204776",ylim=ylim2,xlim= ... ) + } else { + if (TT==1) { + lines(theta,expScore[,i], type="l", lwd=3, las=1,col="#204776") + } else { + lines(theta,expScore[,i], type="l", lwd=2, las=1,col="#204776",lty=2) + } + } + + } else { graphics::lines(theta, expScore[,i],type="l", col=i, lwd=3, pch=1) } @@ -629,7 +522,7 @@ plot.tam.dif <- function(x, items=1:x$nitems, type="expected", aggr <- stats::aggregate(resp_i, list(groupnumber_i), mean, na.rm=TRUE ) obScore_i <- aggr[,2] } - graphics::lines(theta2_i, obScore_i, type="o", lwd=2, pch=1) + #graphics::lines(theta2_i, obScore_i, type="o", lwd=2, pch=1) } } #*********************************************************** @@ -669,21 +562,21 @@ plot.tam.dif <- function(x, items=1:x$nitems, type="expected", kk <- 1 dfr <- dat2 dfr1a <- dfr[ dfr$cat==kk, ] - graphics::plot( ifelse(dfr1a$Theta+0.5>-3,dfr1a$Theta+0.5,NA), dfr1a$P, ylim=c(-.1,1.1),xlim=c(-3,3), - xlab=expression(theta), + graphics::plot( dfr1a$Theta, dfr1a$P, ylim=c(-.1,1.1), + xlab="", col=pall[kk], type="l", xpd=TRUE,axes=F, ... ) axis(1) - axis(2) + axis(2,las=2) grid(nx = NA, ny = NULL, lty = 3, col = "lightgray", lwd = 1) - graphics::lines( ifelse(dfr1a$Theta+0.5>-3,dfr1a$Theta+0.5,NA), dfr1a$P,xlim=c(-3,3), - col=pall[kk], type="l", ... - ) + #graphics::lines( dfr1a$Theta, dfr1a$P, + #col=pall[kk], type="l", ... + #) for (kk in seq(2,K) ){ dfr1a <- dfr[ dfr$cat==kk, ] - graphics::lines( ifelse(dfr1a$Theta+0.5>-3,dfr1a$Theta+0.5,NA), dfr1a$P, col=pall[kk] ) + #graphics::lines( dfr1a$Theta, dfr1a$P, col=pall[kk] ) # graphics::points( dfr1a$Theta, dfr1a$P, pch=kk, col=kk+1 ) } @@ -700,6 +593,24 @@ plot.tam.dif <- function(x, items=1:x$nitems, type="expected", } +png(file = '/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/PNG/icc_base.png',width=800,height=800) +par(xpd=T,mar=c(6.1,5.1,7.6,2.1)) +icc.tam.base(zaza,type = "expected",export=F,package = "graphics",items = 3,axes=F,TT=0) +axis(1) +axis(2,las=2) +title(main='Item: \n "Have you been able to enjoy your normal daily activities?" ', font.main=2) +arrows(x0=-2.65,x1=-3,y0=-0.45,length = 0.15,lwd = 2) +arrows(x0=2.65,x1=3,y0=-0.45,length = 0.15,lwd = 2) +text(x=-2.55,y=-0.45,"Worse\nmental\nhealth",adj=0) +text(x=2.5,y=-0.45,"Better\nmental\nhealth",adj=1) +title(xlab=expression("Latent variable " * theta * "\n (mental health)")) +par(xpd=F) +dev.off() + + +########################## +# CCC DIF +########################## #### CCC @@ -715,132 +626,57 @@ zaza <- tam.mml(resp = zaz[,paste0("item",1:4)]) # base -pdf(file = '/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/PDF/ccc_dif_1.pdf') +pdf(file = '/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/PDF/ccc_dif_1_homo.pdf') par(xpd=F,mar=c(6.1,5.1,7.6,2.1)) -plot.tam.2(zaza,type = "items",export=F,ylab="Probability of response",main=NULL,package = "graphics",items = 3) +plot.tam.2(zaza,type = "items",export=F,ylab="Probability of response in control group (X=0)",main=NULL,package = "graphics",items = 3) segments(x0=zaza$xsi[7,]$xsi,x1=zaza$xsi[7,]$xsi,y0=0,y1=1.1,lty=3) segments(x0=zaza$xsi[8,]$xsi,x1=zaza$xsi[8,]$xsi,y0=0,y1=1.1,lty=3) segments(x0=zaza$xsi[9,]$xsi,x1=zaza$xsi[9,]$xsi,y0=0,y1=1.1,lty=3) -text(x=-2.5,y=0.85,"0",col="#200c23") -text(x=-0.65,y=0.55,"1",col="#62403d") -text(x=0.95,y=0.55,"2",col="#a87b5e") -text(x=2.5,y=0.85,"3",col="#e9bf98") +text(x=-2.5,y=0.85,"0",col="#f7c611") +text(x=-0.65,y=0.55,"1",col="#61b5b7") +text(x=0.95,y=0.55,"2",col="#ca346d") +text(x=2.5,y=0.85,"3",col="#204776") par(xpd=T,mar=c(5.1,4.1,4.1,2.1)) -segments(x0=-3,x1=zaza$xsi[7,]$xsi,y0=1.1,y1=1.1,col="#200c23",lwd=2) -segments(x0=zaza$xsi[7,]$xsi,x1=zaza$xsi[8,]$xsi,y0=1.1,y1=1.1,col="#62403d",lwd=2) -segments(x0=zaza$xsi[8,]$xsi,x1=zaza$xsi[9,]$xsi,y0=1.1,y1=1.1,col="#a87b5e",lwd=2) -segments(x0=zaza$xsi[9,]$xsi,x1=3,y0=1.1,y1=1.1,col="#e9bf98",lwd=2) +segments(x0=-3,x1=zaza$xsi[7,]$xsi,y0=1.1,y1=1.1,col="#f7c611",lwd=2) +segments(x0=zaza$xsi[7,]$xsi,x1=zaza$xsi[8,]$xsi,y0=1.1,y1=1.1,col="#61b5b7",lwd=2) +segments(x0=zaza$xsi[8,]$xsi,x1=zaza$xsi[9,]$xsi,y0=1.1,y1=1.1,col="#ca346d",lwd=2) +segments(x0=zaza$xsi[9,]$xsi,x1=3,y0=1.1,y1=1.1,col="#204776",lwd=2) points(x = zaza$xsi[7,]$xsi, y=1.1,pch=9,cex=1) points(x = zaza$xsi[8,]$xsi, y=1.1,pch=9,cex=1) points(x = zaza$xsi[9,]$xsi, y=1.1,pch=9,cex=1) -text( x=mean(c(-3,zaza$xsi[7,]$xsi)), y=1.2,"0",col="#200c23" ) -text( x=mean(c(zaza$xsi[7,]$xsi,zaza$xsi[8,]$xsi)), y=1.2,"1",col="#62403d" ) -text( x=mean(c(zaza$xsi[8,]$xsi,zaza$xsi[9,]$xsi)), y=1.2,"2",col="#a87b5e" ) -text( x=mean(c(zaza$xsi[9,]$xsi,3)), y=1.2,"3",col="#e9bf98" ) -text( x=mean(c(-3,zaza$xsi[7,]$xsi)),cex=0.7, y=1.15,"Much less than usual",col="#200c23" ) -text( x=mean(c(zaza$xsi[7,]$xsi,zaza$xsi[8,]$xsi)),cex=0.7, y=1.15,"Less so than usual",col="#62403d" ) -text( x=mean(c(zaza$xsi[8,]$xsi,zaza$xsi[9,]$xsi)),cex=0.7, y=1.15,"As much as usual",col="#a87b5e" ) -text( x=mean(c(zaza$xsi[9,]$xsi,3)),cex=0.7, y=1.15,"More so than usual",col="#e9bf98" ) +text( x=mean(c(-3,zaza$xsi[7,]$xsi)), y=1.15,"0",col="#f7c611" ) +text( x=mean(c(zaza$xsi[7,]$xsi,zaza$xsi[8,]$xsi)), y=1.15,"1",col="#61b5b7" ) +text( x=mean(c(zaza$xsi[8,]$xsi,zaza$xsi[9,]$xsi)), y=1.15,"2",col="#ca346d" ) +text( x=mean(c(zaza$xsi[9,]$xsi,3)), y=1.15,"3",col="#204776" ) text(x=zaza$xsi[7,]$xsi,y=1.15,expression(delta["j,1"])) text(x=zaza$xsi[8,]$xsi,y=1.15,expression(delta["j,2"])) text(x=zaza$xsi[9,]$xsi,y=1.15,expression(delta["j,3"])) -text(x = 0,y=1.3,"Most probable response catgory") -CurlyBraces(x0=-2.5, x1=2.5, y0=1.225, y1=1.225, pos = 2, direction = 1, depth=0.05) +text(x = 0,y=1.275,"Most probable response category") +CurlyBraces(x0=-2.5, x1=2.5, y0=1.2, y1=1.2, pos = 2, direction = 1, depth=0.05) arrows(x0=-2.65,x1=-3,y0=-0.35,length = 0.15,lwd = 2) arrows(x0=2.65,x1=3,y0=-0.35,length = 0.15,lwd = 2) text(x=-2.55,y=-0.35,"Worse\nmental\nhealth",adj=0) text(x=2.5,y=-0.35,"Better\nmental\nhealth",adj=1) rect(xleft = 3,xright=5,ybottom = 0,ytop=1.1,col = "white",border = "white") -lines(x=c(zaza$xsi[7,]$xsi,zaza$xsi[7,]$xsi),y=c(0,-0.6),lty=3) -lines(x=c(zaza$xsi[8,]$xsi,zaza$xsi[8,]$xsi),y=c(0,-0.6),lty=3) +lines(x=c(zaza$xsi[7,]$xsi,zaza$xsi[7,]$xsi),y=c(0,-0.275),lty=3) +lines(x=c(zaza$xsi[7,]$xsi,zaza$xsi[7,]$xsi),y=c(-0.4,-0.6),lty=3) +lines(x=c(zaza$xsi[8,]$xsi,zaza$xsi[8,]$xsi),y=c(0,-0.275),lty=3) +lines(x=c(zaza$xsi[8,]$xsi,zaza$xsi[8,]$xsi),y=c(-0.4,-0.6),lty=3) lines(x=c(zaza$xsi[9,]$xsi,zaza$xsi[9,]$xsi),y=c(0,-0.6),lty=3) -title(main='Example item: \n "Have you been able to enjoy your normal daily activities?" ', font.main=2) +title(xlab=expression("Latent variable " * theta * "\n (mental health)"),line=1.8) +title(main='A - Homogeneous DIF', font.main=2) par(xpd=F) dev.off() # DIF homogène - - -pdf(file = '/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/PDF/ccc_dif_2.pdf') -par(xpd=F,mar=c(12.6,5.1,1.1,2.1)) -plot.tam.dif(zaza,type = "items",export=F,ylab="Probability of response",main=NULL,package = "graphics",items = 3) - - -text(x=-2,y=0.85,"0",col="#200c23") -text(x=-.05,y=0.55,"1",col="#62403d") -text(x=1.45,y=0.55,"2",col="#a87b5e") -text(x=2.5,y=0.7,"3",col="#e9bf98") -par(xpd=T,mar=c(5.1,4.1,4.1,2.1)) -segments(x0=-3,x1=.5+zaza$xsi[7,]$xsi,y0=-0.5,col="#200c23",lwd=2) -segments(x0=.5+zaza$xsi[7,]$xsi,x1=.5+zaza$xsi[8,]$xsi,y0=-0.5,col="#62403d",lwd=2) -segments(x0=.5+zaza$xsi[8,]$xsi,x1=.5+zaza$xsi[9,]$xsi,y0=-0.5,col="#a87b5e",lwd=2) -segments(x0=.5+zaza$xsi[9,]$xsi,x1=3,y0=-0.5,col="#e9bf98",lwd=2) - -points(x =.5 + zaza$xsi[7,]$xsi, y=-0.5,pch=9,cex=1) -points(x =.5 + zaza$xsi[8,]$xsi, y=-0.5,pch=9,cex=1) -points(x =.5 + zaza$xsi[9,]$xsi, y=-0.5,pch=9,cex=1) - - -text( x=mean(c(-2.5,zaza$xsi[7,]$xsi)), y=-0.65,"0",col="#200c23" ) -text( x=0.5+mean(c(zaza$xsi[7,]$xsi,zaza$xsi[8,]$xsi)), y=-0.65,"1",col="#62403d" ) -text( x=0.5+mean(c(zaza$xsi[8,]$xsi,zaza$xsi[9,]$xsi)), y=-0.65,"2",col="#a87b5e" ) -text( x=0.5+mean(c(zaza$xsi[9,]$xsi,2.5)), y=-0.65,"3",col="#e9bf98" ) -text( x=mean(c(-2.5,zaza$xsi[7,]$xsi)),cex=0.7, y=-0.6,"Much less than usual",col="#200c23" ) -text( x=.5+mean(c(zaza$xsi[7,]$xsi,zaza$xsi[8,]$xsi)),cex=0.7, y=-0.6,"Less so than usual",col="#62403d" ) -text( x=.5+mean(c(zaza$xsi[8,]$xsi,zaza$xsi[9,]$xsi)),cex=0.7, y=-0.6,"As much as usual",col="#a87b5e" ) -text( x=.5+mean(c(zaza$xsi[9,]$xsi,2.5)),cex=0.7, y=-0.6,"More so than usual",col="#e9bf98" ) - -text(x=0.5+zaza$xsi[7,]$xsi,y=-0.55,expression(delta["j,1"])) -text(x=0.5+zaza$xsi[8,]$xsi,y=-0.55,expression(delta["j,2"])) -text(x=0.5+zaza$xsi[9,]$xsi,y=-0.55,expression(delta["j,3"])) - -arrows(x0=zaza$xsi[7,]$xsi+0.05,x1=zaza$xsi[7,]$xsi+0.5-0.05,y0=0.625,length = 0.1, lwd = 2) -arrows(x0=zaza$xsi[8,]$xsi+0.05,x1=zaza$xsi[8,]$xsi+0.5-0.05,y0=0.625,length = 0.1, lwd = 2) -arrows(x0=zaza$xsi[9,]$xsi+0.05,x1=zaza$xsi[9,]$xsi+0.5-0.05,y0=0.625,length = 0.1, lwd = 2) -text(x=0.25+zaza$xsi[7,]$xsi,y=0.675,expression(gamma["j,1"])) -text(x=0.25+zaza$xsi[8,]$xsi,y=0.675,expression(gamma["j,2"])) -text(x=0.25+zaza$xsi[9,]$xsi,y=0.675,expression(gamma["j,3"])) - -text(x = 0.25,y=-0.75,"Most probable response catgory") -CurlyBraces(x0=-2.25, x1=2.75, y0=-0.675, y1=-0.675, pos = 2, direction = 2, depth=0.05) -arrows(x0=-2.65,x1=-3,y0=-0.35,length = 0.15,lwd = 2) -arrows(x0=2.65,x1=3,y0=-0.35,length = 0.15,lwd = 2) -text(x=-2.55,y=-0.35,"Worse\nmental\nhealth",adj=0) -text(x=2.5,y=-0.35,"Better\nmental\nhealth",adj=1) - -segments(x0=zaza$xsi[7,]$xsi,x1=zaza$xsi[7,]$xsi,y0=0.6,y1=1.6,lty=3) -segments(x0=zaza$xsi[8,]$xsi,x1=zaza$xsi[8,]$xsi,y0=0.6,y1=1.6,lty=3) -segments(x0=zaza$xsi[9,]$xsi,x1=zaza$xsi[9,]$xsi,y0=0.6,y1=1.6,lty=3) -rect(xleft = 3,xright=5,ybottom = -.1,ytop=1.1,col = "white",border = "white") - - -lines(x=.5+c(zaza$xsi[7,]$xsi,zaza$xsi[7,]$xsi),y=c(0.65,-0.5),lty=3) -lines(x=.5+c(zaza$xsi[8,]$xsi,zaza$xsi[8,]$xsi),y=c(0.65,-0.5),lty=3) -lines(x=.5+c(zaza$xsi[9,]$xsi,zaza$xsi[9,]$xsi),y=c(0.65,-0.25),lty=3) -lines(x=.5+c(zaza$xsi[9,]$xsi,zaza$xsi[9,]$xsi),y=c(-0.45,-0.5),lty=3) -par(xpd=F) -dev.off() - - - - - - - - - - - -# CCC DIF HETEROGENE 1 - -plot.tam.difhet1 <- function(x, items=1:x$nitems, type="expected", +plot.tam.dif <- function(x, items=1:x$nitems, type="expected", low=-3, high=3, ngroups=6, groups_by_item=FALSE, wle=NULL, export=TRUE, export.type="png", export.args=list(), observed=TRUE, overlay=FALSE, @@ -851,11 +687,12 @@ plot.tam.difhet1 <- function(x, items=1:x$nitems, type="expected", if ( package=="lattice"){ require_namespace_msg("lattice") } - low <- low-2 - high <- high+2 + # device.Option <- getOption("device") - pall <- c('#200c23', '#62403d', '#a87b5e', '#e9bf98' + pall <- c('#f7c611', '#61b5b7', '#ca346d', '#204776' ) + low <- low-2 + high <- high+2 time1 <- NULL if ( fix.devices ){ old.opt.dev <- getOption("device") @@ -1041,21 +878,21 @@ plot.tam.difhet1 <- function(x, items=1:x$nitems, type="expected", kk <- 1 dfr <- dat2 dfr1a <- dfr[ dfr$cat==kk, ] - graphics::plot( ifelse(dfr1a$Theta+0.5>-3 & dfr1a$Theta+0.5<3, dfr1a$Theta+0.5, NA), dfr1a$P, ylim=c(-.1,1.1),xlim=c(-3,3), - xlab=expression(theta), + graphics::plot( ifelse(dfr1a$Theta+0.5>-3,dfr1a$Theta+0.5,NA), dfr1a$P, ylim=c(-.1,1.1),xlim=c(-3,3), + xlab="", col=pall[kk], type="l", xpd=TRUE,axes=F, ... ) axis(1) - axis(2) + axis(2,las=2) grid(nx = NA, ny = NULL, lty = 3, col = "lightgray", lwd = 1) - graphics::lines( ifelse(dfr1a$Theta+0.5>-3 & dfr1a$Theta+0.5<3, dfr1a$Theta+0.5, NA), dfr1a$P,xlim=c(-3,3), + graphics::lines( ifelse(dfr1a$Theta+0.5>-3,dfr1a$Theta+0.5,NA), dfr1a$P,xlim=c(-3,3), col=pall[kk], type="l", ... ) for (kk in seq(2,K) ){ dfr1a <- dfr[ dfr$cat==kk, ] - graphics::lines(ifelse(dfr1a$Theta+ifelse(kk==3,1,ifelse(kk==4,0.35,0.5))>-3 & dfr1a$Theta+ifelse(kk==3,1,ifelse(kk==4,0.35,0.5))<3, dfr1a$Theta, NA)+ifelse(kk==3,1,ifelse(kk==4,0.35,0.5)), dfr1a$P, col=pall[kk] ) + graphics::lines( ifelse(dfr1a$Theta+0.5>-3,dfr1a$Theta+0.5,NA), dfr1a$P, col=pall[kk] ) # graphics::points( dfr1a$Theta, dfr1a$P, pch=kk, col=kk+1 ) } @@ -1072,52 +909,43 @@ plot.tam.difhet1 <- function(x, items=1:x$nitems, type="expected", } - -# DIF heterogene convergent - - -pdf(file = '/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/PDF/ccc_dif_het1.pdf') +pdf(file = '/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/PDF/ccc_dif_2_homo.pdf') par(xpd=F,mar=c(12.6,5.1,1.1,2.1)) -plot.tam.difhet1(zaza,type = "items",export=F,ylab="Probability of response",main=NULL,package = "graphics",items = 3) +plot.tam.dif(zaza,type = "items",export=F,ylab="Probability of response in treatment group (X=1)",main=NULL,package = "graphics",items = 3) - -text(x=-2,y=0.85,"0",col="#200c23") -text(x=-0.15,y=0.55,"1",col="#62403d") -text(x=1.45+0.475,y=0.55,"2",col="#a87b5e") -text(x=2.5,y=0.7,"3",col="#e9bf98") +text(x=-2,y=0.85,"0",col="#f7c611") +text(x=-.05,y=0.55,"1",col="#61b5b7") +text(x=1.45,y=0.55,"2",col="#ca346d") +text(x=2.5,y=0.7,"3",col="#204776") par(xpd=T,mar=c(5.1,4.1,4.1,2.1)) -segments(x0=-3,x1=.5+zaza$xsi[7,]$xsi,y0=-0.5,col="#200c23",lwd=2) -segments(x0=.5+zaza$xsi[7,]$xsi,x1=.5+zaza$xsi[8,]$xsi+0.25,y0=-0.5,col="#62403d",lwd=2) -segments(x0=.5+zaza$xsi[8,]$xsi+0.25,x1=.5+zaza$xsi[9,]$xsi,y0=-0.5,col="#a87b5e",lwd=2) -segments(x0=.5+zaza$xsi[9,]$xsi,x1=3,y0=-0.5,col="#e9bf98",lwd=2) +segments(x0=-3,x1=.5+zaza$xsi[7,]$xsi,y0=-0.5,col="#f7c611",lwd=2) +segments(x0=.5+zaza$xsi[7,]$xsi,x1=.5+zaza$xsi[8,]$xsi,y0=-0.5,col="#61b5b7",lwd=2) +segments(x0=.5+zaza$xsi[8,]$xsi,x1=.5+zaza$xsi[9,]$xsi,y0=-0.5,col="#ca346d",lwd=2) +segments(x0=.5+zaza$xsi[9,]$xsi,x1=3,y0=-0.5,col="#204776",lwd=2) points(x =.5 + zaza$xsi[7,]$xsi, y=-0.5,pch=9,cex=1) -points(x =.5 + zaza$xsi[8,]$xsi+0.25, y=-0.5,pch=9,cex=1) +points(x =.5 + zaza$xsi[8,]$xsi, y=-0.5,pch=9,cex=1) points(x =.5 + zaza$xsi[9,]$xsi, y=-0.5,pch=9,cex=1) -rect(xleft = 3,xright=5,ybottom = 0,ytop=1.1,col = "white",border = "white") -text( x=mean(c(-2.5,zaza$xsi[7,]$xsi)), y=-0.65,"0",col="#200c23" ) -text( x=0.125+0.5+mean(c(zaza$xsi[7,]$xsi,zaza$xsi[8,]$xsi)), y=-0.65,"1",col="#62403d" ) -text( x=0.5+0.125+mean(c(zaza$xsi[8,]$xsi,zaza$xsi[9,]$xsi)), y=-0.65,"2",col="#a87b5e" ) -text( x=0.5+mean(c(zaza$xsi[9,]$xsi,2.5)), y=-0.65,"3",col="#e9bf98" ) -text( x=mean(c(-2.5,zaza$xsi[7,]$xsi)),cex=0.7, y=-0.6,"Much less than usual",col="#200c23" ) -text( x=0.125+.5+mean(c(zaza$xsi[7,]$xsi,zaza$xsi[8,]$xsi)),cex=0.7, y=-0.6,"Less so than usual",col="#62403d" ) -text( x=.125+.5+mean(c(zaza$xsi[8,]$xsi,zaza$xsi[9,]$xsi)),cex=0.7, y=-0.6,"As much as usual",col="#a87b5e" ) -text( x=.5+mean(c(zaza$xsi[9,]$xsi,2.5)),cex=0.7, y=-0.6,"More so than usual",col="#e9bf98" ) + +text( x=mean(c(-2.5,zaza$xsi[7,]$xsi)), y=-0.55,"0",col="#f7c611" ) +text( x=0.5+mean(c(zaza$xsi[7,]$xsi,zaza$xsi[8,]$xsi)), y=-0.55,"1",col="#61b5b7" ) +text( x=0.5+mean(c(zaza$xsi[8,]$xsi,zaza$xsi[9,]$xsi)), y=-0.55,"2",col="#ca346d" ) +text( x=0.5+mean(c(zaza$xsi[9,]$xsi,2.5)), y=-0.55,"3",col="#204776" ) text(x=0.5+zaza$xsi[7,]$xsi,y=-0.55,expression(delta["j,1"])) -text(x=0.5+zaza$xsi[8,]$xsi+0.25,y=-0.55,expression(delta["j,2"])) +text(x=0.5+zaza$xsi[8,]$xsi,y=-0.55,expression(delta["j,2"])) text(x=0.5+zaza$xsi[9,]$xsi,y=-0.55,expression(delta["j,3"])) arrows(x0=zaza$xsi[7,]$xsi+0.05,x1=zaza$xsi[7,]$xsi+0.5-0.05,y0=0.625,length = 0.1, lwd = 2) -arrows(x0=zaza$xsi[8,]$xsi+0.05,x1=zaza$xsi[8,]$xsi+0.75-0.05,y0=0.625,length = 0.1, lwd = 2,col="darkred") +arrows(x0=zaza$xsi[8,]$xsi+0.05,x1=zaza$xsi[8,]$xsi+0.5-0.05,y0=0.625,length = 0.1, lwd = 2) arrows(x0=zaza$xsi[9,]$xsi+0.05,x1=zaza$xsi[9,]$xsi+0.5-0.05,y0=0.625,length = 0.1, lwd = 2) -text(x=0.25+zaza$xsi[7,]$xsi,y=0.675,expression(gamma["j,1"])) -text(x=0.25+zaza$xsi[8,]$xsi+0.125,y=0.675,expression(gamma["j,2"]),col="darkred") -text(x=0.25+zaza$xsi[9,]$xsi,y=0.675,expression(gamma["j,3"])) +text(x=0.25+zaza$xsi[7,]$xsi,y=0.675,"DIF") +text(x=0.25+zaza$xsi[8,]$xsi,y=0.675,"DIF") +text(x=0.25+zaza$xsi[9,]$xsi,y=0.675,"DIF") -text(x = 0.25,y=-0.75,"Most probable response catgory") -CurlyBraces(x0=-2.25, x1=2.75, y0=-0.675, y1=-0.675, pos = 2, direction = 2, depth=0.05) +text(x = 0.25,y=-0.675,"Most probable response category") +CurlyBraces(x0=-2.25, x1=2.75, y0=-0.6, y1=-0.6, pos = 2, direction = 2, depth=0.05) arrows(x0=-2.65,x1=-3,y0=-0.35,length = 0.15,lwd = 2) arrows(x0=2.65,x1=3,y0=-0.35,length = 0.15,lwd = 2) text(x=-2.55,y=-0.35,"Worse\nmental\nhealth",adj=0) @@ -1126,12 +954,15 @@ text(x=2.5,y=-0.35,"Better\nmental\nhealth",adj=1) segments(x0=zaza$xsi[7,]$xsi,x1=zaza$xsi[7,]$xsi,y0=0.6,y1=1.6,lty=3) segments(x0=zaza$xsi[8,]$xsi,x1=zaza$xsi[8,]$xsi,y0=0.6,y1=1.6,lty=3) segments(x0=zaza$xsi[9,]$xsi,x1=zaza$xsi[9,]$xsi,y0=0.6,y1=1.6,lty=3) -rect(xleft = 3,xright=5,ybottom = -.1,ytop=1.1,col = "white",border = "white") +rect(xleft = 3,xright=5,ybottom = -.1,ytop=1.1,col = "white",border = "white") +title(xlab=expression("Latent variable " * theta * "\n (mental health)"),line=-4.3) -lines(x=.5+c(zaza$xsi[7,]$xsi,zaza$xsi[7,]$xsi),y=c(0.65,-0.5),lty=3) -lines(x=.5+c(zaza$xsi[8,]$xsi+0.25,zaza$xsi[8,]$xsi+0.25),y=c(0.65,-0.5),lty=3) +lines(x=.5+c(zaza$xsi[7,]$xsi,zaza$xsi[7,]$xsi),y=c(0.65,-0.25),lty=3) +lines(x=.5+c(zaza$xsi[8,]$xsi,zaza$xsi[8,]$xsi),y=c(0.65,-0.25),lty=3) lines(x=.5+c(zaza$xsi[9,]$xsi,zaza$xsi[9,]$xsi),y=c(0.65,-0.25),lty=3) +lines(x=.5+c(zaza$xsi[7,]$xsi,zaza$xsi[7,]$xsi),y=c(-0.45,-0.5),lty=3) +lines(x=.5+c(zaza$xsi[8,]$xsi,zaza$xsi[8,]$xsi),y=c(-0.45,-0.5),lty=3) lines(x=.5+c(zaza$xsi[9,]$xsi,zaza$xsi[9,]$xsi),y=c(-0.45,-0.5),lty=3) par(xpd=F) dev.off() @@ -1140,14 +971,20 @@ dev.off() -# CCC DIF HETEROGENE 2 -plot.tam.difhet2 <- function(x, items=1:x$nitems, type="expected", - low=-3, high=3, ngroups=6, groups_by_item=FALSE, - wle=NULL, export=TRUE, export.type="png", - export.args=list(), observed=TRUE, overlay=FALSE, - ask=FALSE, package="lattice", - fix.devices=TRUE, nnodes=100, ...) + + + + + +# CCC DIF HETEROGENE 1 + +plot.tam.difhet1 <- function(x, items=1:x$nitems, type="expected", + low=-3, high=3, ngroups=6, groups_by_item=FALSE, + wle=NULL, export=TRUE, export.type="png", + export.args=list(), observed=TRUE, overlay=FALSE, + ask=FALSE, package="lattice", + fix.devices=TRUE, nnodes=100, ...) { require_namespace_msg("grDevices") if ( package=="lattice"){ @@ -1156,7 +993,7 @@ plot.tam.difhet2 <- function(x, items=1:x$nitems, type="expected", low <- low-2 high <- high+2 # device.Option <- getOption("device") - pall <- c('#200c23', '#62403d', '#a87b5e', '#e9bf98' + pall <- c('#f7c611', '#61b5b7', '#ca346d', '#204776' ) time1 <- NULL if ( fix.devices ){ @@ -1343,21 +1180,21 @@ plot.tam.difhet2 <- function(x, items=1:x$nitems, type="expected", kk <- 1 dfr <- dat2 dfr1a <- dfr[ dfr$cat==kk, ] - graphics::plot(ifelse(dfr1a$Theta-1.125>-3,dfr1a$Theta-1.125,NA) , dfr1a$P, ylim=c(-.1,1.1),xlim=c(-3,3), - xlab=expression(theta), + graphics::plot( ifelse(dfr1a$Theta+0.5>-3 & dfr1a$Theta+0.5<3, dfr1a$Theta+0.5, NA), dfr1a$P, ylim=c(-.1,1.1),xlim=c(-3,3), + xlab="", col=pall[kk], type="l", xpd=TRUE,axes=F, ... ) axis(1) - axis(2) + axis(2,las=2) grid(nx = NA, ny = NULL, lty = 3, col = "lightgray", lwd = 1) - graphics::lines( ifelse(dfr1a$Theta-1.125>-3,dfr1a$Theta-1.125,NA), dfr1a$P,xlim=c(-3,3), - col=pall[kk], type="l", xlim=c(-3,3), ... + graphics::lines( ifelse(dfr1a$Theta+0.5>-3 & dfr1a$Theta+0.5<3, dfr1a$Theta+0.5, NA), dfr1a$P,xlim=c(-3,3), + col=pall[kk], type="l", ... ) for (kk in seq(2,K) ){ dfr1a <- dfr[ dfr$cat==kk, ] - graphics::lines( ifelse(dfr1a$Theta+ifelse(kk==3,1,ifelse(kk==4,0.35,0.5))>-3,dfr1a$Theta,NA)+ifelse(kk==3,1,ifelse(kk==4,0.35,0.5)), dfr1a$P, col=pall[kk] , xlim=c(-3,3)) + graphics::lines(ifelse(dfr1a$Theta+ifelse(kk==3,1.5,ifelse(kk==4,0.4,0.5))>-3 & dfr1a$Theta+ifelse(kk==3,1.5,ifelse(kk==4,0.4,0.5))<3, dfr1a$Theta, NA)+ifelse(kk==3,1.5,ifelse(kk==4,0.4,0.5)), dfr1a$P, col=pall[kk] ) # graphics::points( dfr1a$Theta, dfr1a$P, pch=kk, col=kk+1 ) } @@ -1375,67 +1212,65 @@ plot.tam.difhet2 <- function(x, items=1:x$nitems, type="expected", } -# DIF heterogene divergent +# DIF heterogene convergent -pdf(file = '/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/PDF/ccc_dif_het2.pdf') +pdf(file = '/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/PDF/ccc_dif_2_conv.pdf') par(xpd=F,mar=c(12.6,5.1,1.1,2.1)) -plot.tam.difhet2(zaza,type = "items",export=F,ylab="Probability of response",main=NULL,package = "graphics",items = 3) +plot.tam.difhet1(zaza,type = "items",export=F,ylab="Probability of response in treatment group (X=1)",main=NULL,package = "graphics",items = 3) -text(x=-2.5,y=0.6,"0",col="#200c23") -text(x=-0.15,y=0.55,"1",col="#62403d") -text(x=1.45+0.475,y=0.55,"2",col="#a87b5e") -text(x=2.5,y=0.7,"3",col="#e9bf98") +text(x=-2,y=0.85,"0",col="#f7c611") +text(x=-0.15,y=0.55,"1",col="#61b5b7") +text(x=1.45+1.,y=0.55,"2",col="#ca346d") +text(x=2.5,y=0.7,"3",col="#204776") par(xpd=T,mar=c(5.1,4.1,4.1,2.1)) -segments(x0=-3,x1=-.5+zaza$xsi[7,]$xsi,y0=-0.5,col="#200c23",lwd=2) -segments(x0=-.5+zaza$xsi[7,]$xsi,x1=.5+zaza$xsi[8,]$xsi+0.25,y0=-0.5,col="#62403d",lwd=2) -segments(x0=.5+zaza$xsi[8,]$xsi+0.25,x1=.5+zaza$xsi[9,]$xsi,y0=-0.5,col="#a87b5e",lwd=2) -segments(x0=.5+zaza$xsi[9,]$xsi,x1=3,y0=-0.5,col="#e9bf98",lwd=2) +segments(x0=-3,x1=.5+zaza$xsi[7,]$xsi,y0=-0.5,col="#f7c611",lwd=2) +segments(x0=.5+zaza$xsi[7,]$xsi,x1=.75+zaza$xsi[8,]$xsi+0.25,y0=-0.5,col="#61b5b7",lwd=2) +segments(x0=.75+zaza$xsi[8,]$xsi+0.25,x1=.5+zaza$xsi[9,]$xsi,y0=-0.5,col="#ca346d",lwd=2) +segments(x0=.5+zaza$xsi[9,]$xsi,x1=3,y0=-0.5,col="#204776",lwd=2) -points(x =.5 + zaza$xsi[7,]$xsi-1, y=-0.5,pch=9,cex=1) -points(x =.5 + zaza$xsi[8,]$xsi+0.25, y=-0.5,pch=9,cex=1) +points(x =.5 + zaza$xsi[7,]$xsi, y=-0.5,pch=9,cex=1) +points(x =.75 + zaza$xsi[8,]$xsi+0.25, y=-0.5,pch=9,cex=1) points(x =.5 + zaza$xsi[9,]$xsi, y=-0.5,pch=9,cex=1) rect(xleft = 3,xright=5,ybottom = 0,ytop=1.1,col = "white",border = "white") +text( x=mean(c(-2.5,zaza$xsi[7,]$xsi)), y=-0.55,"0",col="#f7c611" ) +text( x=0.75+mean(c(zaza$xsi[7,]$xsi,zaza$xsi[8,]$xsi)), y=-0.55,"1",col="#61b5b7" ) +text( x=0.75+mean(c(zaza$xsi[8,]$xsi,zaza$xsi[9,]$xsi)), y=-0.55,"2",col="#ca346d" ) +text( x=0.5+mean(c(zaza$xsi[9,]$xsi,2.5)), y=-0.55,"3",col="#204776" ) -text( x=-0.25+mean(c(-3,zaza$xsi[7,]$xsi)), y=-0.65,"0",col="#200c23" ) -text( x=-0.5+0.125+0.5+mean(c(zaza$xsi[7,]$xsi,zaza$xsi[8,]$xsi)), y=-0.65,"1",col="#62403d" ) -text( x=0.5+0.125+mean(c(zaza$xsi[8,]$xsi,zaza$xsi[9,]$xsi)), y=-0.65,"2",col="#a87b5e" ) -text( x=0.5+mean(c(zaza$xsi[9,]$xsi,2.5)), y=-0.65,"3",col="#e9bf98" ) -text( x=-0.25+mean(c(-3,zaza$xsi[7,]$xsi)),cex=0.7, y=-0.6,"Much less than usual",col="#200c23" ) -text( x=-0.5+0.125+.5+mean(c(zaza$xsi[7,]$xsi,zaza$xsi[8,]$xsi)),cex=0.7, y=-0.6,"Less so than usual",col="#62403d" ) -text( x=.125+.5+mean(c(zaza$xsi[8,]$xsi,zaza$xsi[9,]$xsi)),cex=0.7, y=-0.6,"As much as usual",col="#a87b5e" ) -text( x=.5+mean(c(zaza$xsi[9,]$xsi,2.5)),cex=0.7, y=-0.6,"More so than usual",col="#e9bf98" ) - -text(x=-0.5+zaza$xsi[7,]$xsi,y=-0.55,expression(delta["j,1"])) -text(x=0.5+zaza$xsi[8,]$xsi+0.25,y=-0.55,expression(delta["j,2"])) +text(x=0.5+zaza$xsi[7,]$xsi,y=-0.55,expression(delta["j,1"])) +text(x=0.75+zaza$xsi[8,]$xsi+0.25,y=-0.55,expression(delta["j,2"])) text(x=0.5+zaza$xsi[9,]$xsi,y=-0.55,expression(delta["j,3"])) -arrows(x0=zaza$xsi[7,]$xsi-0.05,x1=zaza$xsi[7,]$xsi-0.5+0.05,y0=0.625,length = 0.1, lwd = 2,col="darkred") -arrows(x0=zaza$xsi[8,]$xsi+0.05,x1=zaza$xsi[8,]$xsi+0.75-0.05,y0=0.625,length = 0.1, lwd = 2,col="darkred") +arrows(x0=zaza$xsi[7,]$xsi+0.05,x1=zaza$xsi[7,]$xsi+0.5-0.05,y0=0.625,length = 0.1, lwd = 2) +arrows(x0=zaza$xsi[8,]$xsi+0.05,x1=zaza$xsi[8,]$xsi+1-0.05,y0=0.625,length = 0.1, lwd = 2) arrows(x0=zaza$xsi[9,]$xsi+0.05,x1=zaza$xsi[9,]$xsi+0.5-0.05,y0=0.625,length = 0.1, lwd = 2) -text(x=-0.25+zaza$xsi[7,]$xsi,y=0.675,expression(gamma["j,1"]),col="darkred") -text(x=0.25+zaza$xsi[8,]$xsi+0.125,y=0.675,expression(gamma["j,2"]),col="darkred") -text(x=0.25+zaza$xsi[9,]$xsi,y=0.675,expression(gamma["j,3"])) +text(x=0.25+zaza$xsi[7,]$xsi,y=0.675,"DIF") +text(x=0.25+zaza$xsi[8,]$xsi+0.25,y=0.675,"DIF") +text(x=0.25+zaza$xsi[9,]$xsi,y=0.675,"DIF") + -text(x = 0.25,y=-0.75,"Most probable response catgory") -CurlyBraces(x0=-2.75, x1=2.75, y0=-0.675, y1=-0.675, pos = 2, direction = 2, depth=0.05) +text(x = 0.25,y=-0.675,"Most probable response category") +CurlyBraces(x0=-2.25, x1=2.75, y0=-0.6, y1=-0.6, pos = 2, direction = 2, depth=0.05) arrows(x0=-2.65,x1=-3,y0=-0.35,length = 0.15,lwd = 2) arrows(x0=2.65,x1=3,y0=-0.35,length = 0.15,lwd = 2) text(x=-2.55,y=-0.35,"Worse\nmental\nhealth",adj=0) text(x=2.5,y=-0.35,"Better\nmental\nhealth",adj=1) -rect(xleft = 3,xright=5,ybottom = -.1,ytop=1.1,col = "white",border = "white") segments(x0=zaza$xsi[7,]$xsi,x1=zaza$xsi[7,]$xsi,y0=0.6,y1=1.6,lty=3) segments(x0=zaza$xsi[8,]$xsi,x1=zaza$xsi[8,]$xsi,y0=0.6,y1=1.6,lty=3) segments(x0=zaza$xsi[9,]$xsi,x1=zaza$xsi[9,]$xsi,y0=0.6,y1=1.6,lty=3) +rect(xleft = 3,xright=5,ybottom = -.1,ytop=1.1,col = "white",border = "white") +title(xlab=expression("Latent variable " * theta * "\n (mental health)"),line=-4.3) -lines(x=-.25+c(zaza$xsi[7,]$xsi,zaza$xsi[7,]$xsi)-.25,y=c(0.65,-0.325),lty=3) -lines(x=-.25+c(zaza$xsi[7,]$xsi,zaza$xsi[7,]$xsi)-.25,y=c(-0.375,-0.5),lty=3) -lines(x=.5+c(zaza$xsi[8,]$xsi+0.25,zaza$xsi[8,]$xsi+0.25),y=c(0.65,-0.5),lty=3) +lines(x=.5+c(zaza$xsi[7,]$xsi,zaza$xsi[7,]$xsi),y=c(0.65,-0.25),lty=3) +lines(x=1+c(zaza$xsi[8,]$xsi,zaza$xsi[8,]$xsi),y=c(0.65,-0.25),lty=3) +lines(x=1+c(zaza$xsi[8,]$xsi,zaza$xsi[8,]$xsi),y=c(-0.45,-0.5),lty=3) lines(x=.5+c(zaza$xsi[9,]$xsi,zaza$xsi[9,]$xsi),y=c(0.65,-0.25),lty=3) +lines(x=.5+c(zaza$xsi[7,]$xsi,zaza$xsi[7,]$xsi),y=c(-0.45,-0.5),lty=3) lines(x=.5+c(zaza$xsi[9,]$xsi,zaza$xsi[9,]$xsi),y=c(-0.45,-0.5),lty=3) par(xpd=F) dev.off() @@ -1443,99 +1278,444 @@ dev.off() -########################## -# LM FACTEURS IGNORE DIF -########################## - - -####### SCENARIOS SANS TE - -res.dat.article.ignore.h0 <- res.dat.article.ignore[res.dat.article.ignore$true.beta==0,] -res.dat.article.ignore.h0$prop.dif <- res.dat.article.ignore.h0$nb.dif/res.dat.article.ignore.h0$J - -res.dat.article.ignore.h0.long <- reshape(res.dat.article.ignore.h0,idvar=c("J",'true.beta',"true.gamma","nb.dif","prop.dif"),v.names=c('betahat','bias','typeIerror','power',"coverage")) -rownames(res.dat.article.ignore.h0.long) <- NULL -colnames(res.dat.article.ignore.h0.long)[7:12] <- c("betahat","bias","typeIerror",'power','theoretical.power','coverage') -res.dat.article.ignore.h0.long$prop.dif <- as.numeric(res.dat.article.ignore.h0.long$prop.dif) -res.dat.article.ignore.h0.long$N <- as.numeric(res.dat.article.ignore.h0.long$N) -res.dat.article.ignore.h0.long$true.gamma <- as.numeric(res.dat.article.ignore.h0.long$true.gamma) -res.dat.article.ignore.h0.long$J <- as.numeric(res.dat.article.ignore.h0.long$J) - -res.dat.article.ignore.h0$true.gamma <- as.numeric(res.dat.article.ignore.h0$true.gamma) - -# bias -res.dat.article.ignore.h0.long$abs.bias <- abs(res.dat.article.ignore.h0.long$bias) -summary(lm(abs.bias~true.gamma+prop.dif+N+J,data = res.dat.article.ignore.h0.long)) -summary(lm(abs.bias~true.gamma+prop.dif,data = res.dat.article.ignore.h0.long)) -summary(res.dat.article.ignore.h1.long[res.dat.article.ignore.h1.long$abs.gamma==0.5,"abs.bias"]) -summary(res.dat.article.ignore.h1.long[res.dat.article.ignore.h1.long$abs.gamma==0.3 & res.dat.article.ignore.h1.long$prop.dif<0.3,"abs.bias"]) - -# type I -summary(lm(typeIerror~J+true.gamma+prop.dif+N,data = res.dat.article.ignore.h0.long)) -summary(lm(typeIerror~true.gamma+prop.dif+N,data = res.dat.article.ignore.h0.long)) -res.dat.article.ignore.h0.long$abs.gamma <- abs(res.dat.article.ignore.h0.long$true.gamma) -summary(lm(typeIerror~abs.gamma+prop.dif+N,data = res.dat.article.ignore.h0.long)) - -res.dat.article.ignore.h0[res.dat.article.ignore.h0$true.gamma==-0.5 & - res.dat.article.ignore.h0$prop.dif>0.3,]$typeIerror.300 -res.dat.article.ignore.h0[res.dat.article.ignore.h0$true.gamma==-0.3 & - res.dat.article.ignore.h0$prop.dif<0.3,]$typeIerror.50 - -# coverage -summary(lm(coverage~abs.gamma+prop.dif+N,data = res.dat.article.ignore.h0.long)) - - -####### SCENARIOS AVEC TE -res.dat.article.ignore.h1 <- res.dat.article.ignore[res.dat.article.ignore$true.beta!=0,] -res.dat.article.ignore.h1$prop.dif <- res.dat.article.ignore.h1$nb.dif/res.dat.article.ignore.h1$J +# CCC DIF HETEROGENE 2 -res.dat.article.ignore.h1.long <- reshape(res.dat.article.ignore.h1,idvar=c("J",'true.beta',"true.gamma","nb.dif","prop.dif"),v.names=c('betahat','bias','typeIerror','power',"coverage")) -rownames(res.dat.article.ignore.h1.long) <- NULL -colnames(res.dat.article.ignore.h1.long)[7:12] <- c("betahat","bias","typeIerror",'power','theoretical.power','coverage') -res.dat.article.ignore.h1.long$prop.dif <- as.numeric(res.dat.article.ignore.h1.long$prop.dif) -res.dat.article.ignore.h1.long$N <- as.numeric(res.dat.article.ignore.h1.long$N) -res.dat.article.ignore.h1.long$true.gamma <- as.numeric(res.dat.article.ignore.h1.long$true.gamma) -res.dat.article.ignore.h1.long$J <- as.numeric(res.dat.article.ignore.h1.long$J) +plot.tam.difhet2 <- function(x, items=1:x$nitems, type="expected", + low=-3, high=3, ngroups=6, groups_by_item=FALSE, + wle=NULL, export=TRUE, export.type="png", + export.args=list(), observed=TRUE, overlay=FALSE, + ask=FALSE, package="lattice", + fix.devices=TRUE, nnodes=100, ...) +{ + require_namespace_msg("grDevices") + if ( package=="lattice"){ + require_namespace_msg("lattice") + } + low <- low-2 + high <- high+2 + # device.Option <- getOption("device") + pall <- c('#f7c611', '#61b5b7', '#ca346d', '#204776' + ) + time1 <- NULL + if ( fix.devices ){ + old.opt.dev <- getOption("device") + old.opt.err <- c( getOption("show.error.messages")) + old.par.ask <- graphics::par("ask") + # remember new pars' values + old.par.xpd <- graphics::par("xpd") + old.par.mar <- graphics::par("mar") + on.exit( options("device"=old.opt.dev)) + on.exit( options("show.error.messages"=old.opt.err), add=TRUE) + on.exit( graphics::par("ask"=old.par.ask), add=TRUE) + # restore new pars' values + on.exit( graphics::par("xpd"=old.par.xpd), add=TRUE) + on.exit( graphics::par("mar"=old.par.mar), add=TRUE) + } + + tamobj <- x + ndim <- tamobj$ndim + tammodel <- "mml" + if(is.null(ndim)) { + ndim <- 1 + tammodel <- "jml" + } + if (ndim > 1 ) { + if ( type=="expected"){ + stop ("Expected scores curves are only available for uni-dimensional models") + } + } + + nitems <- tamobj$nitems + + if (ndim==1 ){ + theta <- matrix(seq(low, high, length=nnodes), nrow=nnodes, ncol=ndim) + } else { + nodes <- seq(low, high, length=nnodes) + theta <- as.matrix( expand.grid( as.data.frame( matrix( rep(nodes, ndim), ncol=ndim ) ) ) ) + nnodes <- nrow(theta) + B <- tamobj$B + } + + iIndex <- 1:nitems + A <- tamobj$A + B <- tamobj$B + if (tammodel=="mml") { + xsi <- tamobj$xsi$xsi + } else { + xsi <- tamobj$xsi + } + maxK <- tamobj$maxK + resp <- tamobj$resp + resp.ind <- tamobj$resp.ind + resp[resp.ind==0] <- NA + AXsi <- matrix(0,nrow=nitems,ncol=maxK ) + res <- tam_mml_calc_prob(iIndex=iIndex, A=A, AXsi=AXsi, B=B, xsi=xsi, theta=theta, + nnodes=nnodes, maxK=maxK, recalc=TRUE ) + rprobs <- res[["rprobs"]] + AXsi <- res[["AXsi"]] + cat <- 1:maxK - 1 + + #@@@ define initial empty objects + expScore <- obScore <- wle_intervals <- NULL + theta2 <- NULL + + #**** type='expected' + if ( type=="expected" ){ + expScore <- sapply(1:nitems, function(i) colSums(cat*rprobs[i,,], na.rm=TRUE)) + #-- compute WLE score groups + res <- plot_tam_grouped_wle( tamobj=tamobj, tammodel=tammodel, + wle=wle, ngroups=ngroups, resp=resp ) + wle <- res$wle + theta2 <- res$theta2 + d <- res$d + d1 <- res$d1 + d2 <- res$d2 + groupnumber <- res$groupnumber + ngroups <- res$ngroups + wle_intervals <- res$wle_intervals + #-- compute observed scores + obScore <- apply(d2,2, function(x){ + stats::aggregate(x, list(groupnumber), mean, na.rm=TRUE) + } ) + } + + #---------------------------------------------------- + # adds observed score for type="items" + if (type=="items") { + require_namespace_msg("plyr") + #-- compute WLE score groups + res <- plot_tam_grouped_wle( tamobj=tamobj, tammodel=tammodel, + wle=wle, ngroups=ngroups, resp=resp ) + wle <- res$wle + theta2 <- res$theta2 + d <- res$d + d1 <- res$d1 + d2 <- res$d2 + groupnumber <- res$groupnumber + ngroups <- res$ngroups + + obScore <- lapply(d2, function(item) { + comp_case=stats::complete.cases(item) + item=item[comp_case] + uniq_cats=sort(unique(item)) + plyr::ldply(split(item, groupnumber[comp_case]), .id="group", + function (group) { + ngroup=length(group) + cat_freq=list() + for (catt in uniq_cats) { + cat_freq[[paste0("cat_", catt)]]=sum(group==catt)/ngroup + } + data.frame(cat_freq) + }) + }) + } + + #************************************************* + # begin plot function + probs_plot <- as.list(1:nitems) + names(probs_plot) <- items + + for (i in (1:nitems)[items]) { + #*********************************************************** + #** expected item response curves + if ( type=="expected"){ + if (i==1 || !overlay) { + ylim2 <- c(0,max( tamobj$resp[,i], na.rm=TRUE ) ) + graphics::plot(theta, expScore[,i],,col=12, type="l", lwd=3, las=1, ylab="Score", xlab="Ability", + main=paste("Expected Scores Curve - Item ", colnames(tamobj$resp)[i] ) , + ylim=ylim2, ... ) + } else { + graphics::lines(theta, expScore[,i],type="l", col=i, lwd=3, pch=1) + } + if (observed){ + theta2_i <- theta2 + obScore_i <- obScore[[i]]$x + if (groups_by_item){ + ind_i <- ! is.na(resp[,i]) + resp_i <- resp[ind_i, i, drop=FALSE] + wle_i <- wle[ ind_i ] + res <- plot_tam_grouped_wle( tamobj=tamobj, tammodel=tammodel, + wle=wle_i, ngroups=ngroups, resp=resp_i ) + theta2_i <- res$theta2 + groupnumber_i <- res$groupnumber + aggr <- stats::aggregate(resp_i, list(groupnumber_i), mean, na.rm=TRUE ) + obScore_i <- aggr[,2] + } + graphics::lines(theta2_i, obScore_i, type="o", lwd=2, pch=1) + } + } + #*********************************************************** + + if ( ndim==1 ){ theta0 <- theta } + + if ( type=="items"){ + rprobs.ii <- rprobs[i,,] + rprobs.ii <- rprobs.ii[ rowMeans( is.na(rprobs.ii) ) < 1, ] + K <- nrow(rprobs.ii) + dat2 <- NULL + #************ + if ( ndim > 1 ){ + B.ii <- B[i,,] + ind.ii <- which( colSums( B.ii ) > 0 )[1] + rprobs0.ii <- rprobs.ii + rprobs0.ii <- stats::aggregate( t(rprobs0.ii), list( theta[,ind.ii] ), mean ) + theta0 <- rprobs0.ii[,1,drop=FALSE] + rprobs.ii <- t( rprobs0.ii[,-1] ) + } + probs_plot[[i]] <- rprobs.ii + #************** + for (kk in 1:K){ + dat2a <- data.frame( "Theta"=theta0[,1], "cat"=kk, "P"=rprobs.ii[kk,] ) + dat2 <- rbind(dat2, dat2a) + } + auto.key <- NULL + simple.key <- paste0("Cat", 1:K - 1) + auto.key <- simple.key + dat2$time <- dat2$cat + dat2$time1 <- paste0("Cat", dat2$time ) + + simple.key <- FALSE + Kpercol <- K + # package graphics + if ( package=="graphics" ){ + kk <- 1 + dfr <- dat2 + dfr1a <- dfr[ dfr$cat==kk, ] + graphics::plot(ifelse(dfr1a$Theta-1.125>-3,dfr1a$Theta-1.125,NA) , dfr1a$P, ylim=c(-.1,1.1),xlim=c(-3,3), + xlab="", + col=pall[kk], type="l", xpd=TRUE,axes=F, ... + ) + axis(1) + axis(2,las=2) + grid(nx = NA, + ny = NULL, + lty = 3, col = "lightgray", lwd = 1) + graphics::lines( ifelse(dfr1a$Theta-1.125>-3,dfr1a$Theta-1.125,NA), dfr1a$P,xlim=c(-3,3), + col=pall[kk], type="l", xlim=c(-3,3), ... + ) + for (kk in seq(2,K) ){ + dfr1a <- dfr[ dfr$cat==kk, ] + graphics::lines(ifelse(dfr1a$Theta+ifelse(kk==3,1.5,ifelse(kk==4,0.4,0.5))>-3 & dfr1a$Theta+ifelse(kk==3,1.5,ifelse(kk==4,0.4,0.5))<3, dfr1a$Theta, NA)+ifelse(kk==3,1.5,ifelse(kk==4,0.4,0.5)), dfr1a$P, col=pall[kk] ) + # graphics::points( dfr1a$Theta, dfr1a$P, pch=kk, col=kk+1 ) + + } + } + + + #*************************************** + + } + #*************** + graphics::par(ask=ask) + } # end item ii + #************************************************* + +} -res.dat.article.ignore.h1$true.gamma <- as.numeric(res.dat.article.ignore.h1$true.gamma) -# bias -res.dat.article.ignore.h1.long$abs.bias <- abs(res.dat.article.ignore.h1.long$bias) -res.dat.article.ignore.h1.long$abs.gamma <- abs(res.dat.article.ignore.h1.long$true.gamma) -res.dat.article.ignore.h1.long$sign.gamma <- sign(res.dat.article.ignore.h1.long$true.gamma) -summary(lm(abs.bias~abs.gamma+sign.gamma+prop.dif+true.beta+N,data = res.dat.article.ignore.h1.long)) -summary(lm(abs.bias~abs.gamma+prop.dif,data = res.dat.article.ignore.h1.long)) -summary(res.dat.article.ignore.h1.long$abs.bias) +# DIF heterogene divergent +pdf(file = '/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/PDF/ccc_dif_2_div.pdf') +par(xpd=F,mar=c(12.6,5.1,1.1,2.1)) +plot.tam.difhet2(zaza,type = "items",export=F,ylab="Probability of response in treatment group (X=1)",main=NULL,package = "graphics",items = 3) -# coverage -res.dat.article.ignore.h1.long$masks <- res.dat.article.ignore.h1.long$true.beta/res.dat.article.ignore.h1.long$true.gamma>0 -res.dat.article.ignore.h1.long$masks <- 1*res.dat.article.ignore.h1.long$masks -summary(lm(coverage~abs.gamma+prop.dif+true.beta+N+masks,data = res.dat.article.ignore.h1.long)) -summary(lm(coverage~abs.gamma+prop.dif+N+masks,data = res.dat.article.ignore.h1.long)) -summary(lm(coverage~abs.gamma+prop.dif+N,data = res.dat.article.ignore.h1.long)) -summary(res.dat.article.ignore.h1.long$coverage) +text(x=-2.5,y=0.6,"0",col="#f7c611") +text(x=-0.15,y=0.55,"1",col="#61b5b7") +text(x=1.45+0.475,y=0.55,"2",col="#ca346d") +text(x=2.5,y=0.7,"3",col="#204776") +par(xpd=T,mar=c(5.1,4.1,4.1,2.1)) +segments(x0=-3,x1=-.5+zaza$xsi[7,]$xsi,y0=-0.5,col="#f7c611",lwd=2) +segments(x0=-.5+zaza$xsi[7,]$xsi,x1=.75+zaza$xsi[8,]$xsi+0.25,y0=-0.5,col="#61b5b7",lwd=2) +segments(x0=.75+zaza$xsi[8,]$xsi+0.25,x1=.5+zaza$xsi[9,]$xsi,y0=-0.5,col="#ca346d",lwd=2) +segments(x0=.5+zaza$xsi[9,]$xsi,x1=3,y0=-0.5,col="#204776",lwd=2) +points(x =.5 + zaza$xsi[7,]$xsi-1, y=-0.5,pch=9,cex=1) +points(x =.75 + zaza$xsi[8,]$xsi+0.25, y=-0.5,pch=9,cex=1) +points(x =.5 + zaza$xsi[9,]$xsi, y=-0.5,pch=9,cex=1) +rect(xleft = 3,xright=5,ybottom = 0,ytop=1.1,col = "white",border = "white") -# power -res.dat.article.ignore.h1.long$powerdif <- as.numeric(res.dat.article.ignore.h1.long$power)-as.numeric(res.dat.article.ignore.h1.long$theoretical.power) -summary(lm(powerdif~masks*prop.dif+true.beta+N+masks*abs.gamma,data = res.dat.article.ignore.h1.long)) -summary(res.dat.article.ignore.h1.long[res.dat.article.ignore.h1.long$masks==0,]$powerdif) -summary(res.dat.article.ignore.h1.long[res.dat.article.ignore.h1.long$masks==1,]$powerdif) -# bias +- -summary(lm(bias~abs.gamma+prop.dif+true.beta+N+masks,data = res.dat.article.ignore.h1.long)) -summary(lm(bias~masks,data = res.dat.article.ignore.h1.long)) -summary(res.dat.article.ignore.h1.long[res.dat.article.ignore.h1.long$masks==1,]$bias) -summary(res.dat.article.ignore.h1.long[res.dat.article.ignore.h1.long$masks==0,]$bias) +text( x=-0.25+mean(c(-3,zaza$xsi[7,]$xsi)), y=-0.55,"0",col="#f7c611" ) +text( x=-0.25+0.125+0.5+mean(c(zaza$xsi[7,]$xsi,zaza$xsi[8,]$xsi)), y=-0.55,"1",col="#61b5b7" ) +text( x=0.5+0.25+mean(c(zaza$xsi[8,]$xsi,zaza$xsi[9,]$xsi)), y=-0.55,"2",col="#ca346d" ) +text( x=0.5+mean(c(zaza$xsi[9,]$xsi,2.5)), y=-0.55,"3",col="#204776" ) -########################## -# DESCRIPTION PCM-DIF -########################## -####### SCENARIOS SANS TE +text(x=-0.5+zaza$xsi[7,]$xsi,y=-0.55,expression(delta["j,1"])) +text(x=0.75+zaza$xsi[8,]$xsi+0.25,y=-0.55,expression(delta["j,2"])) +text(x=0.5+zaza$xsi[9,]$xsi,y=-0.55,expression(delta["j,3"])) + +arrows(x0=zaza$xsi[7,]$xsi-0.05,x1=zaza$xsi[7,]$xsi-0.5+0.05,y0=0.625,length = 0.1, lwd = 2) +arrows(x0=zaza$xsi[8,]$xsi+0.05,x1=zaza$xsi[8,]$xsi+1-0.05,y0=0.625,length = 0.1, lwd = 2) +arrows(x0=zaza$xsi[9,]$xsi+0.05,x1=zaza$xsi[9,]$xsi+0.5-0.05,y0=0.625,length = 0.1, lwd = 2) +text(x=-0.25+zaza$xsi[7,]$xsi,y=0.675,"DIF") +text(x=0.25+zaza$xsi[8,]$xsi+0.25,y=0.675,"DIF") +text(x=0.25+zaza$xsi[9,]$xsi,y=0.675,"DIF") + + +text(x = 0.25,y=-0.675,"Most probable response category") +CurlyBraces(x0=-2.25, x1=2.75, y0=-0.6, y1=-0.6, pos = 2, direction = 2, depth=0.05) +arrows(x0=-2.65,x1=-3,y0=-0.35,length = 0.15,lwd = 2) +arrows(x0=2.65,x1=3,y0=-0.35,length = 0.15,lwd = 2) +text(x=-2.55,y=-0.35,"Worse\nmental\nhealth",adj=0) +text(x=2.5,y=-0.35,"Better\nmental\nhealth",adj=1) + +segments(x0=zaza$xsi[7,]$xsi,x1=zaza$xsi[7,]$xsi,y0=0.6,y1=1.6,lty=3) +segments(x0=zaza$xsi[8,]$xsi,x1=zaza$xsi[8,]$xsi,y0=0.6,y1=1.6,lty=3) +segments(x0=zaza$xsi[9,]$xsi,x1=zaza$xsi[9,]$xsi,y0=0.6,y1=1.6,lty=3) + +rect(xleft = 3,xright=5,ybottom = -.1,ytop=1.1,col = "white",border = "white") +title(xlab=expression("Latent variable " * theta * "\n (mental health)"),line=-4.3) + +lines(x=-.5+c(zaza$xsi[7,]$xsi,zaza$xsi[7,]$xsi),y=c(0.65,-0.5),lty=3) +lines(x=1+c(zaza$xsi[8,]$xsi,zaza$xsi[8,]$xsi),y=c(0.65,-0.25),lty=3) +lines(x=1+c(zaza$xsi[8,]$xsi,zaza$xsi[8,]$xsi),y=c(-0.45,-0.5),lty=3) +lines(x=.5+c(zaza$xsi[9,]$xsi,zaza$xsi[9,]$xsi),y=c(0.65,-0.25),lty=3) +lines(x=.5+c(zaza$xsi[9,]$xsi,zaza$xsi[9,]$xsi),y=c(-0.45,-0.5),lty=3) +par(xpd=F) +dev.off() + + + + +########################## +# LM FACTEURS IGNORE DIF +########################## + + +####### SCENARIOS SANS TE + +res.dat.article.ignore.long <- reshape(res.dat.article.ignore,idvar=c("J",'true.beta',"true.gamma","nb.dif","prop.dif"),v.names=c('betahat','bias','typeIerror','power',"coverage")) +rownames(res.dat.article.ignore.long) <- NULL +colnames(res.dat.article.ignore.long)[6:11] <- c("betahat","bias","typeIerror",'power','theoretical.power','coverage') + +res.dat.article.ignore.h0 <- res.dat.article.ignore[res.dat.article.ignore$true.beta==0,] +res.dat.article.ignore.h0$prop.dif <- res.dat.article.ignore.h0$nb.dif/res.dat.article.ignore.h0$J + +res.dat.article.ignore.h0.long <- reshape(res.dat.article.ignore.h0,idvar=c("J",'true.beta',"true.gamma","nb.dif","prop.dif"),v.names=c('betahat','bias','typeIerror','power',"coverage")) +rownames(res.dat.article.ignore.h0.long) <- NULL +colnames(res.dat.article.ignore.h0.long)[7:12] <- c("betahat","bias","typeIerror",'power','theoretical.power','coverage') +res.dat.article.ignore.h0.long$prop.dif <- as.numeric(res.dat.article.ignore.h0.long$prop.dif) +res.dat.article.ignore.h0.long$N <- as.numeric(res.dat.article.ignore.h0.long$N) +res.dat.article.ignore.h0.long$true.gamma <- as.numeric(res.dat.article.ignore.h0.long$true.gamma) +res.dat.article.ignore.h0.long$J <- as.numeric(res.dat.article.ignore.h0.long$J) + +res.dat.article.ignore.h0$true.gamma <- as.numeric(res.dat.article.ignore.h0$true.gamma) + +res.dat.article.ignore.h0.long$bias <- abs(res.dat.article.ignore.h0.long$bias) +res.dat.article.ignore.h0.long$betahat <- abs(res.dat.article.ignore.h0.long$betahat) + +# bias +res.dat.article.ignore.h0.long$abs.bias <- abs(res.dat.article.ignore.h0.long$bias) +summary(lm(abs.bias~true.gamma+prop.dif+N+J,data = res.dat.article.ignore.h0.long)) +summary(lm(abs.bias~true.gamma+prop.dif,data = res.dat.article.ignore.h0.long)) +summary(res.dat.article.ignore.h1.long[res.dat.article.ignore.h1.long$abs.gamma==0.5,"abs.bias"]) +summary(res.dat.article.ignore.h1.long[res.dat.article.ignore.h1.long$abs.gamma==0.3 & res.dat.article.ignore.h1.long$prop.dif<0.3,"abs.bias"]) + +# type I +summary(lm(typeIerror~J+true.gamma+prop.dif+N,data = res.dat.article.ignore.h0.long)) +summary(lm(typeIerror~true.gamma+prop.dif+N,data = res.dat.article.ignore.h0.long)) +res.dat.article.ignore.h0.long$abs.gamma <- abs(res.dat.article.ignore.h0.long$true.gamma) +summary(lm(typeIerror~abs.gamma+prop.dif+N,data = res.dat.article.ignore.h0.long)) + +res.dat.article.ignore.h0[res.dat.article.ignore.h0$true.gamma==-0.5 & + res.dat.article.ignore.h0$prop.dif>0.3,]$typeIerror.300 +res.dat.article.ignore.h0[res.dat.article.ignore.h0$true.gamma==-0.3 & + res.dat.article.ignore.h0$prop.dif<0.3,]$typeIerror.50 + +# coverage +summary(lm(coverage~abs.gamma+prop.dif+N,data = res.dat.article.ignore.h0.long)) + + +####### SCENARIOS AVEC TE MASK + + + +res.dat.article.ignore.h1 <- res.dat.article.ignore[res.dat.article.ignore$true.beta!=0 & res.dat.article.ignore$true.gamma>0,] +res.dat.article.ignore.h1$prop.dif <- res.dat.article.ignore.h1$nb.dif/res.dat.article.ignore.h1$J + +res.dat.article.ignore.h1.long <- reshape(res.dat.article.ignore.h1,idvar=c("J",'true.beta',"true.gamma","nb.dif","prop.dif"),v.names=c('betahat','bias','typeIerror','power',"coverage")) +rownames(res.dat.article.ignore.h1.long) <- NULL +colnames(res.dat.article.ignore.h1.long)[7:12] <- c("betahat","bias","typeIerror",'power','theoretical.power','coverage') +res.dat.article.ignore.h1.long$prop.dif <- as.numeric(res.dat.article.ignore.h1.long$prop.dif) +res.dat.article.ignore.h1.long$N <- as.numeric(res.dat.article.ignore.h1.long$N) +res.dat.article.ignore.h1.long$true.gamma <- as.numeric(res.dat.article.ignore.h1.long$true.gamma) +res.dat.article.ignore.h1.long$J <- as.numeric(res.dat.article.ignore.h1.long$J) + +res.dat.article.ignore.h1$true.gamma <- as.numeric(res.dat.article.ignore.h1$true.gamma) + +# bias +res.dat.article.ignore.h1.long$abs.bias <- abs(res.dat.article.ignore.h1.long$bias) +res.dat.article.ignore.h1.long$abs.gamma <- abs(res.dat.article.ignore.h1.long$true.gamma) +res.dat.article.ignore.h1.long$sign.gamma <- sign(res.dat.article.ignore.h1.long$true.gamma) +summary(lm(abs.bias~true.gamma+prop.dif+true.beta+N+J,data = res.dat.article.ignore.h1.long)) + + +# coverage +res.dat.article.ignore.h1.long$masks <- res.dat.article.ignore.h1.long$true.beta/res.dat.article.ignore.h1.long$true.gamma>0 +res.dat.article.ignore.h1.long$masks <- 1*res.dat.article.ignore.h1.long$masks +summary(lm(coverage~true.gamma+prop.dif+true.beta+N+J,data = res.dat.article.ignore.h1.long)) + +# power +res.dat.article.ignore.h1.long$powerdif <- as.numeric(res.dat.article.ignore.h1.long$power)-as.numeric(res.dat.article.ignore.h1.long$theoretical.power) +summary(lm(powerdif~true.gamma+prop.dif+true.beta+N+J,data = res.dat.article.ignore.h1.long)) + +# bias +- +summary(lm(bias~abs.gamma+prop.dif+true.beta+N+masks,data = res.dat.article.ignore.h1.long)) +summary(lm(bias~masks,data = res.dat.article.ignore.h1.long)) +summary(res.dat.article.ignore.h1.long[res.dat.article.ignore.h1.long$masks==1,]$bias) +summary(res.dat.article.ignore.h1.long[res.dat.article.ignore.h1.long$masks==0,]$bias) + + + +####### SCENARIOS AVEC TE AMPLIFY + + + +res.dat.article.ignore.h1 <- res.dat.article.ignore[res.dat.article.ignore$true.beta!=0 & res.dat.article.ignore$true.gamma<0,] +res.dat.article.ignore.h1$prop.dif <- res.dat.article.ignore.h1$nb.dif/res.dat.article.ignore.h1$J + +res.dat.article.ignore.h1.long <- reshape(res.dat.article.ignore.h1,idvar=c("J",'true.beta',"true.gamma","nb.dif","prop.dif"),v.names=c('betahat','bias','typeIerror','power',"coverage")) +rownames(res.dat.article.ignore.h1.long) <- NULL +colnames(res.dat.article.ignore.h1.long)[7:12] <- c("betahat","bias","typeIerror",'power','theoretical.power','coverage') +res.dat.article.ignore.h1.long$prop.dif <- as.numeric(res.dat.article.ignore.h1.long$prop.dif) +res.dat.article.ignore.h1.long$N <- as.numeric(res.dat.article.ignore.h1.long$N) +res.dat.article.ignore.h1.long$true.gamma <- as.numeric(res.dat.article.ignore.h1.long$true.gamma) +res.dat.article.ignore.h1.long$J <- as.numeric(res.dat.article.ignore.h1.long$J) + +res.dat.article.ignore.h1$true.gamma <- as.numeric(res.dat.article.ignore.h1$true.gamma) +res.dat.article.ignore.h1$true.beta <- as.numeric(res.dat.article.ignore.h1$true.beta) + +# bias +res.dat.article.ignore.h1.long$abs.bias <- abs(res.dat.article.ignore.h1.long$bias) +res.dat.article.ignore.h1.long$abs.gamma <- abs(res.dat.article.ignore.h1.long$true.gamma) +res.dat.article.ignore.h1.long$sign.gamma <- sign(res.dat.article.ignore.h1.long$true.gamma) +summary(lm(abs.bias~abs.gamma+prop.dif+true.beta+N+J,data = res.dat.article.ignore.h1.long)) + + +# coverage +res.dat.article.ignore.h1.long$masks <- res.dat.article.ignore.h1.long$true.beta/res.dat.article.ignore.h1.long$true.gamma>0 +res.dat.article.ignore.h1.long$masks <- 1*res.dat.article.ignore.h1.long$masks +summary(lm(coverage~abs.gamma+prop.dif+true.beta+N+J,data = res.dat.article.ignore.h1.long)) + +# power +res.dat.article.ignore.h1.long$powerdif <- as.numeric(res.dat.article.ignore.h1.long$power)-as.numeric(res.dat.article.ignore.h1.long$theoretical.power) +summary(lm(powerdif~abs.gamma+prop.dif+true.beta+N+J,data = res.dat.article.ignore.h1.long)) + +# bias +- +summary(lm(bias~abs.gamma+prop.dif+true.beta+N+masks,data = res.dat.article.ignore.h1.long)) +summary(lm(bias~masks,data = res.dat.article.ignore.h1.long)) +summary(res.dat.article.ignore.h1.long[res.dat.article.ignore.h1.long$masks==1,]$bias) +summary(res.dat.article.ignore.h1.long[res.dat.article.ignore.h1.long$masks==0,]$bias) + + +########################## +# DESCRIPTION PCM-DIF +########################## + +####### SCENARIOS SANS TE res.dat.article.dif.h0 <- res.dat.article.dif[res.dat.article.dif$true.beta==0,] res.dat.article.dif.h0$prop.dif <- res.dat.article.dif.h0$nb.dif/res.dat.article.dif.h0$J @@ -1598,15 +1778,13 @@ res.dat.article.rosali.dif.h0.long$prop.dif <- as.numeric(res.dat.article.rosali res.dat.article.rosali.dif.h0.long$N <- as.numeric(res.dat.article.rosali.dif.h0.long$N) res.dat.article.rosali.dif.h0.long$true.gamma <- as.numeric(res.dat.article.rosali.dif.h0.long$true.gamma) res.dat.article.rosali.dif.h0.long$J <- as.numeric(res.dat.article.rosali.dif.h0.long$J) +res.dat.article.rosali.dif.h0.long$abs.gamma <- abs(res.dat.article.rosali.dif.h0.long$true.gamma) res.dat.article.rosali.dif.h0$true.gamma <- as.numeric(res.dat.article.rosali.dif.h0$true.gamma) # bias -res.dat.article.rosali.dif.h0.long$abs.gamma <- abs(res.dat.article.rosali.dif.h0.long$true.gamma) res.dat.article.rosali.dif.h0.long$abs.bias <- abs(res.dat.article.rosali.dif.h0.long$bias) -summary(lm(abs.bias~abs.gamma+prop.dif+N+J,data = res.dat.article.rosali.dif.h0.long)) -summary(lm(abs.bias~abs.gamma+prop.dif+N,data = res.dat.article.rosali.dif.h0.long)) -summary(res.dat.article.rosali.dif.h0.long$abs.bias) +summary(lm(abs.bias~true.gamma+prop.dif+true.beta+N+J,data = res.dat.article.rosali.dif.h0.long)) # type I summary(lm(typeIerror~abs.gamma+prop.dif+N+J,data = res.dat.article.rosali.dif.h0.long)) @@ -1618,9 +1796,9 @@ summary(as.numeric(res.dat.article.rosali.dif.h0.long$typeIerror)) summary(lm(coverage~abs.gamma+prop.dif+N,data = res.dat.article.ignore.h0.long)) -####### SCENARIOS AVEC TE +####### SCENARIOS AVEC TE MASK -res.dat.article.rosali.dif.h1 <- res.dat.article.rosali.dif[res.dat.article.rosali.dif$true.beta!=0,] +res.dat.article.rosali.dif.h1 <- res.dat.article.rosali.dif[res.dat.article.rosali.dif$true.beta!=0 & res.dat.article.rosali.dif$true.gamma>0,] res.dat.article.rosali.dif.h1$prop.dif <- res.dat.article.rosali.dif.h1$nb.dif/res.dat.article.rosali.dif.h1$J res.dat.article.rosali.dif.h1.long <- reshape(res.dat.article.rosali.dif.h1,idvar=c("J",'true.beta',"true.gamma","nb.dif","prop.dif"),v.names=c('betahat','bias','typeIerror','power',"coverage")) @@ -1637,27 +1815,17 @@ res.dat.article.rosali.dif.h1$true.gamma <- as.numeric(res.dat.article.rosali.di res.dat.article.rosali.dif.h1.long$abs.bias <- abs(res.dat.article.rosali.dif.h1.long$bias) res.dat.article.rosali.dif.h1.long$abs.gamma <- abs(res.dat.article.rosali.dif.h1.long$true.gamma) res.dat.article.rosali.dif.h1.long$sign.gamma <- sign(res.dat.article.rosali.dif.h1.long$true.gamma) -summary(lm(abs.bias~abs.gamma+sign.gamma+prop.dif+true.beta+N+J,data = res.dat.article.rosali.dif.h1.long)) -summary(lm(abs.bias~abs.gamma+prop.dif+N,data = res.dat.article.rosali.dif.h1.long)) -summary(res.dat.article.rosali.dif.h1.long$abs.bias) - +summary(lm(abs.bias~true.gamma+prop.dif+true.beta+N+J,data = res.dat.article.rosali.dif.h1.long)) # coverage -res.dat.article.rosali.dif.h1.long$masks <- res.dat.article.rosali.dif.h1.long$true.beta/res.dat.article.rosali.dif.h1.long$true.gamma>0 -res.dat.article.rosali.dif.h1.long$masks <- 1*res.dat.article.rosali.dif.h1.long$masks -summary(lm(coverage~abs.gamma+prop.dif+true.beta+N+masks+J,data = res.dat.article.rosali.dif.h1.long)) -summary(lm(coverage~abs.gamma+prop.dif+N,data = res.dat.article.rosali.dif.h1.long)) +summary(lm(coverage~true.gamma+prop.dif+true.beta+N+J,data = res.dat.article.rosali.dif.h1.long)) summary(res.dat.article.rosali.dif.h1.long$coverage) # power res.dat.article.rosali.dif.h1.long$powerdif <- as.numeric(res.dat.article.rosali.dif.h1.long$power)-as.numeric(res.dat.article.rosali.dif.h1.long$theoretical.power) -summary(lm(powerdif~masks+true.beta+N,data = res.dat.article.rosali.dif.h1.long)) -summary(lm(abs(powerdif)~N+prop.dif,data = res.dat.article.rosali.dif.h1.long)) - -summary(res.dat.article.rosali.dif.h1.long[res.dat.article.rosali.dif.h1.long$masks==0,]$powerdif) -summary(res.dat.article.rosali.dif.h1.long[res.dat.article.rosali.dif.h1.long$masks==1,]$powerdif) +summary(lm(powerdif~true.gamma+prop.dif+true.beta+N+J,data = res.dat.article.rosali.dif.h1.long)) # bias +- summary(lm(bias~abs.gamma+prop.dif+true.beta+N+masks,data = res.dat.article.rosali.dif.h1.long)) @@ -1666,6 +1834,40 @@ summary(res.dat.article.rosali.dif.h1.long[res.dat.article.rosali.dif.h1.long$ma summary(res.dat.article.rosali.dif.h1.long[res.dat.article.rosali.dif.h1.long$masks==0,]$bias) + + +####### SCENARIOS AVEC TE AMPLIFY + +res.dat.article.rosali.dif.h1 <- res.dat.article.rosali.dif[res.dat.article.rosali.dif$true.beta!=0 & res.dat.article.rosali.dif$true.gamma<0,] +res.dat.article.rosali.dif.h1$prop.dif <- res.dat.article.rosali.dif.h1$nb.dif/res.dat.article.rosali.dif.h1$J + +res.dat.article.rosali.dif.h1.long <- reshape(res.dat.article.rosali.dif.h1,idvar=c("J",'true.beta',"true.gamma","nb.dif","prop.dif"),v.names=c('betahat','bias','typeIerror','power',"coverage")) +rownames(res.dat.article.rosali.dif.h1.long) <- NULL +colnames(res.dat.article.rosali.dif.h1.long)[7:12] <- c("betahat","bias","typeIerror",'power','theoretical.power','coverage') +res.dat.article.rosali.dif.h1.long$prop.dif <- as.numeric(res.dat.article.rosali.dif.h1.long$prop.dif) +res.dat.article.rosali.dif.h1.long$N <- as.numeric(res.dat.article.rosali.dif.h1.long$N) +res.dat.article.rosali.dif.h1.long$true.gamma <- as.numeric(res.dat.article.rosali.dif.h1.long$true.gamma) +res.dat.article.rosali.dif.h1.long$J <- as.numeric(res.dat.article.rosali.dif.h1.long$J) + +res.dat.article.rosali.dif.h1$true.gamma <- as.numeric(res.dat.article.rosali.dif.h1$true.gamma) + +# bias +res.dat.article.rosali.dif.h1.long$abs.bias <- abs(res.dat.article.rosali.dif.h1.long$bias) +res.dat.article.rosali.dif.h1.long$abs.gamma <- abs(res.dat.article.rosali.dif.h1.long$true.gamma) +res.dat.article.rosali.dif.h1.long$sign.gamma <- sign(res.dat.article.rosali.dif.h1.long$true.gamma) +summary(lm(abs.bias~abs.gamma+prop.dif+true.beta+N+J,data = res.dat.article.rosali.dif.h1.long)) + +# coverage +summary(lm(coverage~abs.gamma+prop.dif+true.beta+N+J,data = res.dat.article.rosali.dif.h1.long)) + +summary(res.dat.article.rosali.dif.h1.long$coverage) + + +# power +res.dat.article.rosali.dif.h1.long$powerdif <- as.numeric(res.dat.article.rosali.dif.h1.long$power)-as.numeric(res.dat.article.rosali.dif.h1.long$theoretical.power) +summary(lm(powerdif~abs.gamma+prop.dif+true.beta+N+J,data = res.dat.article.rosali.dif.h1.long)) + + ########################## # Plots ROSALI vs ignore ########################## @@ -1744,6 +1946,7 @@ res.dat.article.residif.dif.h0.long$prop.dif <- as.numeric(res.dat.article.resid res.dat.article.residif.dif.h0.long$N <- as.numeric(res.dat.article.residif.dif.h0.long$N) res.dat.article.residif.dif.h0.long$true.gamma <- as.numeric(res.dat.article.residif.dif.h0.long$true.gamma) res.dat.article.residif.dif.h0.long$J <- as.numeric(res.dat.article.residif.dif.h0.long$J) +res.dat.article.residif.dif.h0.long$abs.gamma <- as.numeric(res.dat.article.residif.dif.h0.long$true.gamma) res.dat.article.residif.dif.h0$true.gamma <- as.numeric(res.dat.article.residif.dif.h0$true.gamma) @@ -1764,9 +1967,9 @@ summary(as.numeric(res.dat.article.residif.dif.h0.long$typeIerror)) summary(lm(coverage~abs.gamma+prop.dif+N,data = res.dat.article.ignore.h0.long)) -####### SCENARIOS AVEC TE +####### SCENARIOS AVEC TE MASK -res.dat.article.residif.dif.h1 <- res.dat.article.residif.dif[res.dat.article.residif.dif$true.beta!=0,] +res.dat.article.residif.dif.h1 <- res.dat.article.residif.dif[res.dat.article.residif.dif$true.beta!=0 & res.dat.article.residif.dif$true.gamma>0,] res.dat.article.residif.dif.h1$prop.dif <- res.dat.article.residif.dif.h1$nb.dif/res.dat.article.residif.dif.h1$J res.dat.article.residif.dif.h1.long <- reshape(res.dat.article.residif.dif.h1,idvar=c("J",'true.beta',"true.gamma","nb.dif","prop.dif"),v.names=c('betahat','bias','typeIerror','power',"coverage")) @@ -1782,45 +1985,50 @@ res.dat.article.residif.dif.h1$true.gamma <- as.numeric(res.dat.article.residif. # bias res.dat.article.residif.dif.h1.long$abs.bias <- abs(res.dat.article.residif.dif.h1.long$bias) res.dat.article.residif.dif.h1.long$abs.gamma <- abs(res.dat.article.residif.dif.h1.long$true.gamma) -res.dat.article.residif.dif.h1.long$sign.gamma <- sign(res.dat.article.residif.dif.h1.long$true.gamma) -res.dat.article.residif.dif.h1.long$masks <- res.dat.article.residif.dif.h1.long$true.beta/res.dat.article.residif.dif.h1.long$true.gamma>0 -res.dat.article.residif.dif.h1.long$masks <- 1*res.dat.article.residif.dif.h1.long$masks -summary(lm(abs.bias~abs.gamma+masks+prop.dif+true.beta+N+J,data = res.dat.article.residif.dif.h1.long)) -summary(lm(abs.bias~abs.gamma+masks+prop.dif+N,data = res.dat.article.residif.dif.h1.long)) -summary(res.dat.article.residif.dif.h1.long$abs.bias) +summary(lm(abs.bias~true.gamma+prop.dif+true.beta+N+J,data = res.dat.article.residif.dif.h1.long)) -summary(res.dat.article.residif.dif.h1.long[res.dat.article.residif.dif.h1.long$masks==1,]$abs.bias) -summary(res.dat.article.residif.dif.h1.long[res.dat.article.residif.dif.h1.long$masks==0,]$abs.bias) +# coverage -res.dat.article.rosali.dif.h1.long$masks <- 1*(res.dat.article.rosali.dif.h1.long$true.gamma>0) -res.dat.article.rosali.dif.h1.long$abs.bias <- abs(as.numeric(res.dat.article.rosali.dif.h1.long$bias)) -summary(res.dat.article.rosali.dif.h1.long[res.dat.article.rosali.dif.h1.long$masks==1,]$abs.bias) -summary(res.dat.article.rosali.dif.h1.long[res.dat.article.rosali.dif.h1.long$masks==0,]$abs.bias) +summary(lm(coverage~true.gamma+prop.dif+true.beta+N+J,data = res.dat.article.residif.dif.h1.long)) -# coverage +# power +res.dat.article.residif.dif.h1.long$powerdif <- as.numeric(res.dat.article.residif.dif.h1.long$power)-as.numeric(res.dat.article.residif.dif.h1.long$theoretical.power) +summary(lm(powerdif~true.gamma+prop.dif+true.beta+N+J,data = res.dat.article.residif.dif.h1.long)) -summary(lm(coverage~abs.gamma+prop.dif+true.beta+N+masks+J,data = res.dat.article.residif.dif.h1.long)) -summary(lm(coverage~abs.gamma+prop.dif+N+masks,data = res.dat.article.residif.dif.h1.long)) -summary(res.dat.article.residif.dif.h1.long$coverage) -# power -res.dat.article.residif.dif.h1.long$powerdif <- as.numeric(res.dat.article.residif.dif.h1.long$power)-as.numeric(res.dat.article.residif.dif.h1.long$theoretical.power) -summary(lm(powerdif~masks+true.beta+N+abs.gamma+prop.dif+J,data = res.dat.article.residif.dif.h1.long)) -summary(lm(powerdif~masks+true.beta+N,data = res.dat.article.residif.dif.h1.long)) -summary(res.dat.article.residif.dif.h1.long[res.dat.article.residif.dif.h1.long$masks==0,]$powerdif) -summary(res.dat.article.residif.dif.h1.long[res.dat.article.residif.dif.h1.long$masks==1,]$powerdif) +####### SCENARIOS AVEC TE AMPLIFY -# bias +- -summary(lm(bias~abs.gamma+prop.dif+true.beta+N+masks,data = res.dat.article.residif.dif.h1.long)) -summary(lm(bias~masks,data = res.dat.article.residif.dif.h1.long)) -summary(res.dat.article.residif.dif.h1.long[res.dat.article.residif.dif.h1.long$masks==1,]$bias) -summary(res.dat.article.residif.dif.h1.long[res.dat.article.residif.dif.h1.long$masks==0,]$bias) +res.dat.article.residif.dif.h1 <- res.dat.article.residif.dif[res.dat.article.residif.dif$true.beta!=0 & res.dat.article.residif.dif$true.gamma<0,] +res.dat.article.residif.dif.h1$prop.dif <- res.dat.article.residif.dif.h1$nb.dif/res.dat.article.residif.dif.h1$J + +res.dat.article.residif.dif.h1.long <- reshape(res.dat.article.residif.dif.h1,idvar=c("J",'true.beta',"true.gamma","nb.dif","prop.dif"),v.names=c('betahat','bias','typeIerror','power',"coverage")) +rownames(res.dat.article.residif.dif.h1.long) <- NULL +colnames(res.dat.article.residif.dif.h1.long)[7:12] <- c("betahat","bias","typeIerror",'power','theoretical.power','coverage') +res.dat.article.residif.dif.h1.long$prop.dif <- as.numeric(res.dat.article.residif.dif.h1.long$prop.dif) +res.dat.article.residif.dif.h1.long$N <- as.numeric(res.dat.article.residif.dif.h1.long$N) +res.dat.article.residif.dif.h1.long$true.gamma <- as.numeric(res.dat.article.residif.dif.h1.long$true.gamma) +res.dat.article.residif.dif.h1.long$J <- as.numeric(res.dat.article.residif.dif.h1.long$J) + +res.dat.article.residif.dif.h1$true.gamma <- as.numeric(res.dat.article.residif.dif.h1$true.gamma) + +# bias +res.dat.article.residif.dif.h1.long$abs.bias <- abs(res.dat.article.residif.dif.h1.long$bias) +res.dat.article.residif.dif.h1.long$abs.gamma <- abs(res.dat.article.residif.dif.h1.long$true.gamma) +summary(lm(bias~abs.gamma+prop.dif+true.beta+N+J,data = res.dat.article.residif.dif.h1.long)) + +# coverage + +summary(lm(coverage~abs.gamma+prop.dif+true.beta+N+J,data = res.dat.article.residif.dif.h1.long)) +# power +res.dat.article.residif.dif.h1.long$powerdif <- as.numeric(res.dat.article.residif.dif.h1.long$power)-as.numeric(res.dat.article.residif.dif.h1.long$theoretical.power) +summary(lm(powerdif~abs.gamma+prop.dif+true.beta+N+J,data = res.dat.article.residif.dif.h1.long)) + ########################## # Plots RESIDIF vs ignore ########################## @@ -3703,73 +3911,1659 @@ bp.dat.power.ignore.magnif <- as.numeric(res.dat.article.2[res.dat.article.2$N== bp.dat.power.nodif <- as.numeric(res.dat.article.nodif.2[res.dat.article.nodif.2$true.beta>0,"power"])-as.numeric(res.dat.article.nodif.2[res.dat.article.nodif.2$true.beta>0,"theoretical.power"]) -bp.dat.power.rosali.mask <- as.numeric(res.dat.article.rosali.2[res.dat.article.rosali.2$N==300 & res.dat.article.rosali.2$abs.gamma==0.5 & res.dat.article.rosali.2$prop.dif<0.3 & res.dat.article.rosali.2$true.beta>0 & res.dat.article.rosali.2$true.gamma>0 & res.dat.article.rosali.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.rosali.2[res.dat.article.rosali.2$N==300 & res.dat.article.rosali.2$abs.gamma==0.5 & res.dat.article.rosali.2$prop.dif<0.3 & res.dat.article.rosali.2$true.beta>0 & res.dat.article.rosali.2$true.gamma>0 & res.dat.article.rosali.2$nb.dif!=0,"theoretical.power"]) -bp.dat.power.rosali.magnif <- as.numeric(res.dat.article.rosali.2[res.dat.article.rosali.2$N==300 & res.dat.article.rosali.2$abs.gamma==0.5 & res.dat.article.rosali.2$prop.dif<0.3 & res.dat.article.rosali.2$true.beta>0 & res.dat.article.rosali.2$true.gamma<0 & res.dat.article.rosali.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.rosali.2[res.dat.article.rosali.2$N==300 & res.dat.article.rosali.2$abs.gamma==0.5 & res.dat.article.rosali.2$prop.dif<0.3 & res.dat.article.rosali.2$true.beta>0 & res.dat.article.rosali.2$true.gamma<0 & res.dat.article.rosali.2$nb.dif!=0,"theoretical.power"]) +bp.dat.power.rosali.mask <- as.numeric(res.dat.article.rosali.2[res.dat.article.rosali.2$N==300 & res.dat.article.rosali.2$abs.gamma==0.5 & res.dat.article.rosali.2$prop.dif<0.3 & res.dat.article.rosali.2$true.beta>0 & res.dat.article.rosali.2$true.gamma>0 & res.dat.article.rosali.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.rosali.2[res.dat.article.rosali.2$N==300 & res.dat.article.rosali.2$abs.gamma==0.5 & res.dat.article.rosali.2$prop.dif<0.3 & res.dat.article.rosali.2$true.beta>0 & res.dat.article.rosali.2$true.gamma>0 & res.dat.article.rosali.2$nb.dif!=0,"theoretical.power"]) +bp.dat.power.rosali.magnif <- as.numeric(res.dat.article.rosali.2[res.dat.article.rosali.2$N==300 & res.dat.article.rosali.2$abs.gamma==0.5 & res.dat.article.rosali.2$prop.dif<0.3 & res.dat.article.rosali.2$true.beta>0 & res.dat.article.rosali.2$true.gamma<0 & res.dat.article.rosali.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.rosali.2[res.dat.article.rosali.2$N==300 & res.dat.article.rosali.2$abs.gamma==0.5 & res.dat.article.rosali.2$prop.dif<0.3 & res.dat.article.rosali.2$true.beta>0 & res.dat.article.rosali.2$true.gamma<0 & res.dat.article.rosali.2$nb.dif!=0,"theoretical.power"]) + +bp.dat.power.residif.mask <- as.numeric(res.dat.article.residif.2[res.dat.article.residif.2$N==300 & res.dat.article.residif.2$abs.gamma==0.5 & res.dat.article.residif.2$prop.dif<0.3 & res.dat.article.residif.2$true.beta>0 & res.dat.article.residif.2$true.gamma>0 & res.dat.article.residif.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.residif.2[res.dat.article.residif.2$N==300 & res.dat.article.residif.2$abs.gamma==0.5 & res.dat.article.residif.2$prop.dif<0.3 & res.dat.article.residif.2$true.beta>0 & res.dat.article.residif.2$true.gamma>0 & res.dat.article.residif.2$nb.dif!=0,"theoretical.power"]) +bp.dat.power.residif.magnif <- as.numeric(res.dat.article.residif.2[res.dat.article.residif.2$N==300 & res.dat.article.residif.2$abs.gamma==0.5 & res.dat.article.residif.2$prop.dif<0.3 & res.dat.article.residif.2$true.beta>0 & res.dat.article.residif.2$true.gamma<0 & res.dat.article.residif.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.residif.2[res.dat.article.residif.2$N==300 & res.dat.article.residif.2$abs.gamma==0.5 & res.dat.article.residif.2$prop.dif<0.3 & res.dat.article.residif.2$true.beta>0 & res.dat.article.residif.2$true.gamma<0 & res.dat.article.residif.2$nb.dif!=0,"theoretical.power"]) + +bp.dat.power.dif.mask <- as.numeric(res.dat.article.dif.2[res.dat.article.dif.2$N==300 & res.dat.article.dif.2$abs.gamma==0.5 & res.dat.article.dif.2$prop.dif<0.3 & res.dat.article.dif.2$true.beta>0 & res.dat.article.dif.2$true.gamma>0 & res.dat.article.dif.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.dif.2[res.dat.article.dif.2$N==300 & res.dat.article.dif.2$abs.gamma==0.5 & res.dat.article.dif.2$prop.dif<0.3 & res.dat.article.dif.2$true.beta>0 & res.dat.article.dif.2$true.gamma>0 & res.dat.article.dif.2$nb.dif!=0,"theoretical.power"]) +bp.dat.power.dif.magnif <- as.numeric(res.dat.article.dif.2[res.dat.article.dif.2$N==300 & res.dat.article.dif.2$abs.gamma==0.5 & res.dat.article.dif.2$prop.dif<0.3 & res.dat.article.dif.2$true.beta>0 & res.dat.article.dif.2$true.gamma<0 & res.dat.article.dif.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.dif.2[res.dat.article.dif.2$N==300 & res.dat.article.dif.2$abs.gamma==0.5 & res.dat.article.dif.2$prop.dif<0.3 & res.dat.article.dif.2$true.beta>0 & res.dat.article.dif.2$true.gamma<0 & res.dat.article.dif.2$nb.dif!=0,"theoretical.power"]) + +bp.dat.power <- data.frame(power=c(bp.dat.power.nodif,bp.dat.power.ignore.mask,bp.dat.power.ignore.magnif,bp.dat.power.rosali.mask,bp.dat.power.rosali.magnif,bp.dat.power.residif.mask,bp.dat.power.residif.magnif, + bp.dat.power.dif.mask,bp.dat.power.dif.magnif), + method=c(rep("NO DIF",12), + rep("MASK1",4),rep("AMPLIFY1",4), + rep("MASK2",4),rep("AMPLIFY2",4), + rep("MASK3",4),rep("AMPLIFY3",4), + rep("MASK4",4),rep("AMPLIFY4",4) )) +bp.dat.power$method <- factor(bp.dat.power$method,levels=c("NO DIF","MASK1","AMPLIFY1","MASK2","AMPLIFY2","MASK3","AMPLIFY3","MASK4","AMPLIFY4")) + + +boxplot(bp.dat.power$power~bp.dat.power$method,xlab="",pch=3,main="", + ylab="",ylim=c(-1,1),yaxt="n",xaxt="n", + cex.lab=1.45,cex.main=1.5,cex.axis=1.15, + col=c("#e69875","#798A5D","#D4E8B5","#798A5D","#D4E8B5","#798A5D","#D4E8B5","#798A5D","#D4E8B5") + ,border=c("#CD5E35","#697850","#697850","#697850","#697850","#697850","#697850","#697850","#697850"), + width=c(0.8,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4), + at=c(1,2,2.5,3.25,3.75,4.5,5,5.75,6.25)) +axis(2,seq(-1,1,0.25),cex.axis=1.45) +axis(1,c(1,2.25,3.5,4.75,6),labels=c("NO DIF","IGNORE-DIF","ROSALI","RESIDIF","PCM-DIF"),cex.axis=1.45) +abline(h=0,lty=2,col='#595959',lwd=2) + + +## DIF 05 50% + +bp.dat.power.ignore.mask <- as.numeric(res.dat.article.2[res.dat.article.2$N==300 & res.dat.article.2$abs.gamma==0.5 & res.dat.article.2$prop.dif<0.3 & res.dat.article.2$true.gamma>0 & res.dat.article.2$true.beta>0 & res.dat.article.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.2[res.dat.article.2$N==300 & res.dat.article.2$abs.gamma==0.5 & res.dat.article.2$prop.dif<0.3 & res.dat.article.2$true.gamma>0 & res.dat.article.2$true.beta>0 & res.dat.article.2$nb.dif!=0,"theoretical.power"]) +bp.dat.power.ignore.magnif <- as.numeric(res.dat.article.2[res.dat.article.2$N==300 & res.dat.article.2$abs.gamma==0.5 & res.dat.article.2$prop.dif<0.3 & res.dat.article.2$true.gamma<0 & res.dat.article.2$true.beta>0 & res.dat.article.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.2[res.dat.article.2$N==300 & res.dat.article.2$abs.gamma==0.5 & res.dat.article.2$prop.dif<0.3 & res.dat.article.2$true.gamma<0 & res.dat.article.2$true.beta>0 & res.dat.article.2$nb.dif!=0,"theoretical.power"]) + +bp.dat.power.nodif <- as.numeric(res.dat.article.nodif.2[res.dat.article.nodif.2$true.beta>0,"power"])-as.numeric(res.dat.article.nodif.2[res.dat.article.nodif.2$true.beta>0,"theoretical.power"]) + +bp.dat.power.rosali.mask <- as.numeric(res.dat.article.rosali.2[res.dat.article.rosali.2$N==300 & res.dat.article.rosali.2$abs.gamma==0.5 & res.dat.article.rosali.2$prop.dif<0.3 & res.dat.article.rosali.2$true.beta>0 & res.dat.article.rosali.2$true.gamma>0 & res.dat.article.rosali.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.rosali.2[res.dat.article.rosali.2$N==300 & res.dat.article.rosali.2$abs.gamma==0.5 & res.dat.article.rosali.2$prop.dif<0.3 & res.dat.article.rosali.2$true.beta>0 & res.dat.article.rosali.2$true.gamma>0 & res.dat.article.rosali.2$nb.dif!=0,"theoretical.power"]) +bp.dat.power.rosali.magnif <- as.numeric(res.dat.article.rosali.2[res.dat.article.rosali.2$N==300 & res.dat.article.rosali.2$abs.gamma==0.5 & res.dat.article.rosali.2$prop.dif<0.3 & res.dat.article.rosali.2$true.beta>0 & res.dat.article.rosali.2$true.gamma<0 & res.dat.article.rosali.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.rosali.2[res.dat.article.rosali.2$N==300 & res.dat.article.rosali.2$abs.gamma==0.5 & res.dat.article.rosali.2$prop.dif<0.3 & res.dat.article.rosali.2$true.beta>0 & res.dat.article.rosali.2$true.gamma<0 & res.dat.article.rosali.2$nb.dif!=0,"theoretical.power"]) + +bp.dat.power.residif.mask <- as.numeric(res.dat.article.residif.2[res.dat.article.residif.2$N==300 & res.dat.article.residif.2$abs.gamma==0.5 & res.dat.article.residif.2$prop.dif<0.3 & res.dat.article.residif.2$true.beta>0 & res.dat.article.residif.2$true.gamma>0 & res.dat.article.residif.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.residif.2[res.dat.article.residif.2$N==300 & res.dat.article.residif.2$abs.gamma==0.5 & res.dat.article.residif.2$prop.dif<0.3 & res.dat.article.residif.2$true.beta>0 & res.dat.article.residif.2$true.gamma>0 & res.dat.article.residif.2$nb.dif!=0,"theoretical.power"]) +bp.dat.power.residif.magnif <- as.numeric(res.dat.article.residif.2[res.dat.article.residif.2$N==300 & res.dat.article.residif.2$abs.gamma==0.5 & res.dat.article.residif.2$prop.dif<0.3 & res.dat.article.residif.2$true.beta>0 & res.dat.article.residif.2$true.gamma<0 & res.dat.article.residif.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.residif.2[res.dat.article.residif.2$N==300 & res.dat.article.residif.2$abs.gamma==0.5 & res.dat.article.residif.2$prop.dif<0.3 & res.dat.article.residif.2$true.beta>0 & res.dat.article.residif.2$true.gamma<0 & res.dat.article.residif.2$nb.dif!=0,"theoretical.power"]) + +bp.dat.power.dif.mask <- as.numeric(res.dat.article.dif.2[res.dat.article.dif.2$N==300 & res.dat.article.dif.2$abs.gamma==0.5 & res.dat.article.dif.2$prop.dif<0.3 & res.dat.article.dif.2$true.beta>0 & res.dat.article.dif.2$true.gamma>0 & res.dat.article.dif.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.dif.2[res.dat.article.dif.2$N==300 & res.dat.article.dif.2$abs.gamma==0.5 & res.dat.article.dif.2$prop.dif<0.3 & res.dat.article.dif.2$true.beta>0 & res.dat.article.dif.2$true.gamma>0 & res.dat.article.dif.2$nb.dif!=0,"theoretical.power"]) +bp.dat.power.dif.magnif <- as.numeric(res.dat.article.dif.2[res.dat.article.dif.2$N==300 & res.dat.article.dif.2$abs.gamma==0.5 & res.dat.article.dif.2$prop.dif<0.3 & res.dat.article.dif.2$true.beta>0 & res.dat.article.dif.2$true.gamma<0 & res.dat.article.dif.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.dif.2[res.dat.article.dif.2$N==300 & res.dat.article.dif.2$abs.gamma==0.5 & res.dat.article.dif.2$prop.dif<0.3 & res.dat.article.dif.2$true.beta>0 & res.dat.article.dif.2$true.gamma<0 & res.dat.article.dif.2$nb.dif!=0,"theoretical.power"]) + +bp.dat.power <- data.frame(power=c(bp.dat.power.nodif,bp.dat.power.ignore.mask,bp.dat.power.ignore.magnif,bp.dat.power.rosali.mask,bp.dat.power.rosali.magnif,bp.dat.power.residif.mask,bp.dat.power.residif.magnif, + bp.dat.power.dif.mask,bp.dat.power.dif.magnif), + method=c(rep("NO DIF",12), + rep("MASK1",4),rep("AMPLIFY1",4), + rep("MASK2",4),rep("AMPLIFY2",4), + rep("MASK3",4),rep("AMPLIFY3",4), + rep("MASK4",4),rep("AMPLIFY4",4) )) +bp.dat.power$method <- factor(bp.dat.power$method,levels=c("NO DIF","MASK1","AMPLIFY1","MASK2","AMPLIFY2","MASK3","AMPLIFY3","MASK4","AMPLIFY4")) + + +boxplot(bp.dat.power$power~bp.dat.power$method,xlab="",pch=3,main="", + ylab="",ylim=c(-1,1),yaxt="n",xaxt="n", + cex.lab=1.45,cex.main=1.5,cex.axis=1.15, + col=c("#e69875","#798A5D","#D4E8B5","#798A5D","#D4E8B5","#798A5D","#D4E8B5","#798A5D","#D4E8B5") + ,border=c("#CD5E35","#697850","#697850","#697850","#697850","#697850","#697850","#697850","#697850"), + width=c(0.8,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4), + at=c(1,2,2.5,3.25,3.75,4.5,5,5.75,6.25)) +axis(2,seq(-1,1,0.25),cex.axis=1.45) +axis(1,c(1,2.25,3.5,4.75,6),labels=c("NO DIF","IGNORE-DIF","ROSALI","RESIDIF","PCM-DIF"),cex.axis=1.45) +abline(h=0,lty=2,col='#595959',lwd=2) + + +par(mfrow=c(1,1)) + + +########################## +# ICC DIF +########################## + + +# DIF HOMOGENE + +zaz <- read.csv("/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N300/scenario_XA_300.csv") +zaza0 <- tam.mml(resp = zaz[zaz$TT==0,paste0("item",1:4)]) +zaza1 <- tam.mml(resp = zaz[zaz$TT==1,paste0("item",1:4)]) + +icc.tam.difh <- function(x, items=1:x$nitems, type="expected",TT=F, + low=-4, high=4, ngroups=6, groups_by_item=FALSE, + wle=NULL, export=TRUE, export.type="png", + export.args=list(), observed=TRUE, overlay=FALSE, + ask=FALSE, package="lattice", + fix.devices=TRUE, nnodes=100, ...) +{ + require_namespace_msg("grDevices") + if ( package=="lattice"){ + require_namespace_msg("lattice") + } + + # device.Option <- getOption("device") + + time1 <- NULL + pall <- c('#f7c611', '#61b5b7', '#ca346d', '#204776' + ) + if ( fix.devices ){ + old.opt.dev <- getOption("device") + old.opt.err <- c( getOption("show.error.messages")) + old.par.ask <- graphics::par("ask") + # remember new pars' values + old.par.xpd <- graphics::par("xpd") + old.par.mar <- graphics::par("mar") + on.exit( options("device"=old.opt.dev)) + on.exit( options("show.error.messages"=old.opt.err), add=TRUE) + on.exit( graphics::par("ask"=old.par.ask), add=TRUE) + # restore new pars' values + on.exit( graphics::par("xpd"=old.par.xpd), add=TRUE) + on.exit( graphics::par("mar"=old.par.mar), add=TRUE) + } + + tamobj <- x + ndim <- tamobj$ndim + tammodel <- "mml" + if(is.null(ndim)) { + ndim <- 1 + tammodel <- "jml" + } + if (ndim > 1 ) { + if ( type=="expected"){ + stop ("Expected scores curves are only available for uni-dimensional models") + } + } + + nitems <- tamobj$nitems + + if (ndim==1 ){ + theta <- matrix(seq(low, high, length=nnodes), nrow=nnodes, ncol=ndim) + } else { + nodes <- seq(-3, 3, length=nnodes) + theta <- as.matrix( expand.grid( as.data.frame( matrix( rep(nodes, ndim), ncol=ndim ) ) ) ) + nnodes <- nrow(theta) + B <- tamobj$B + } + + iIndex <- 1:nitems + A <- tamobj$A + B <- tamobj$B + if (tammodel=="mml") { + xsi <- tamobj$xsi$xsi + } else { + xsi <- tamobj$xsi + } + maxK <- tamobj$maxK + resp <- tamobj$resp + resp.ind <- tamobj$resp.ind + resp[resp.ind==0] <- NA + AXsi <- matrix(0,nrow=nitems,ncol=maxK ) + res <- tam_mml_calc_prob(iIndex=iIndex, A=A, AXsi=AXsi, B=B, xsi=xsi, theta=theta, + nnodes=nnodes, maxK=maxK, recalc=TRUE ) + rprobs <- res[["rprobs"]] + AXsi <- res[["AXsi"]] + cat <- 1:maxK - 1 + + #@@@ define initial empty objects + expScore <- obScore <- wle_intervals <- NULL + theta2 <- NULL + + #**** type='expected' + if ( type=="expected" ){ + expScore <- sapply(1:nitems, function(i) colSums(cat*rprobs[i,,], na.rm=TRUE)) + #-- compute WLE score groups + res <- plot_tam_grouped_wle( tamobj=tamobj, tammodel=tammodel, + wle=wle, ngroups=ngroups, resp=resp ) + wle <- res$wle + theta2 <- res$theta2 + d <- res$d + d1 <- res$d1 + d2 <- res$d2 + groupnumber <- res$groupnumber + ngroups <- res$ngroups + wle_intervals <- res$wle_intervals + #-- compute observed scores + obScore <- apply(d2,2, function(x){ + stats::aggregate(x, list(groupnumber), mean, na.rm=TRUE) + } ) + } + + #---------------------------------------------------- + # adds observed score for type="items" + if (type=="items") { + require_namespace_msg("plyr") + #-- compute WLE score groups + res <- plot_tam_grouped_wle( tamobj=tamobj, tammodel=tammodel, + wle=wle, ngroups=ngroups, resp=resp ) + wle <- res$wle + theta2 <- res$theta2 + d <- res$d + d1 <- res$d1 + d2 <- res$d2 + groupnumber <- res$groupnumber + ngroups <- res$ngroups + + obScore <- lapply(d2, function(item) { + comp_case=stats::complete.cases(item) + item=item[comp_case] + uniq_cats=sort(unique(item)) + plyr::ldply(split(item, groupnumber[comp_case]), .id="group", + function (group) { + ngroup=length(group) + cat_freq=list() + for (catt in uniq_cats) { + cat_freq[[paste0("cat_", catt)]]=sum(group==catt)/ngroup + } + data.frame(cat_freq) + }) + }) + } + + #************************************************* + # begin plot function + probs_plot <- as.list(1:nitems) + names(probs_plot) <- items + + for (i in (1:nitems)[items]) { + #*********************************************************** + #** expected item response curves + if ( type=="expected"){ + if (i==1 || !overlay) { + ylim2 <- c(0,max( tamobj$resp[,i], na.rm=TRUE ) ) + if (!TT) { + graphics::plot(theta, expScore[,i], type="l", lwd=3, las=1, ylab="Expected item response", + main=NULL, + col='#204776',ylim=ylim2,xlim=c(-4,4),xlab="", ... ) + } else { + lines(theta,expScore[,i], type="l", lwd=3, las=1,col='#ca346d') + } + + + } else { + graphics::lines(theta, expScore[,i],type="l", col=i, lwd=3, pch=1) + } + if (observed){ + theta2_i <- theta2 + obScore_i <- obScore[[i]]$x + if (groups_by_item){ + ind_i <- ! is.na(resp[,i]) + resp_i <- resp[ind_i, i, drop=FALSE] + wle_i <- wle[ ind_i ] + res <- plot_tam_grouped_wle( tamobj=tamobj, tammodel=tammodel, + wle=wle_i, ngroups=ngroups, resp=resp_i ) + theta2_i <- res$theta2 + groupnumber_i <- res$groupnumber + aggr <- stats::aggregate(resp_i, list(groupnumber_i), mean, na.rm=TRUE ) + obScore_i <- aggr[,2] + } + #graphics::lines(theta2_i, obScore_i, type="o", lwd=2, pch=1) + } + } + #*********************************************************** + + if ( ndim==1 ){ theta0 <- theta } + + if ( type=="items"){ + rprobs.ii <- rprobs[i,,] + rprobs.ii <- rprobs.ii[ rowMeans( is.na(rprobs.ii) ) < 1, ] + K <- nrow(rprobs.ii) + dat2 <- NULL + #************ + if ( ndim > 1 ){ + B.ii <- B[i,,] + ind.ii <- which( colSums( B.ii ) > 0 )[1] + rprobs0.ii <- rprobs.ii + rprobs0.ii <- stats::aggregate( t(rprobs0.ii), list( theta[,ind.ii] ), mean ) + theta0 <- rprobs0.ii[,1,drop=FALSE] + rprobs.ii <- t( rprobs0.ii[,-1] ) + } + probs_plot[[i]] <- rprobs.ii + #************** + for (kk in 1:K){ + dat2a <- data.frame( "Theta"=theta0[,1], "cat"=kk, "P"=rprobs.ii[kk,] ) + dat2 <- rbind(dat2, dat2a) + } + auto.key <- NULL + simple.key <- paste0("Cat", 1:K - 1) + auto.key <- simple.key + dat2$time <- dat2$cat + dat2$time1 <- paste0("Cat", dat2$time ) + + simple.key <- FALSE + Kpercol <- K + # package graphics + if ( package=="graphics" ){ + kk <- 1 + dfr <- dat2 + dfr1a <- dfr[ dfr$cat==kk, ] + graphics::plot( dfr1a$Theta, dfr1a$P, ylim=c(-.1,1.1), + xlab=expression(theta), + col=pall[kk], type="l", xpd=TRUE,axes=F, ... + ) + axis(1) + axis(2) + grid(nx = NA, + ny = NULL, + lty = 3, col = "lightgray", lwd = 1) + graphics::lines( dfr1a$Theta, dfr1a$P, + col=pall[kk], type="l", ... + ) + for (kk in seq(2,K) ){ + dfr1a <- dfr[ dfr$cat==kk, ] + graphics::lines( dfr1a$Theta, dfr1a$P, col=pall[kk] ) + # graphics::points( dfr1a$Theta, dfr1a$P, pch=kk, col=kk+1 ) + + } + } + + + #*************************************** + + } + #*************** + graphics::par(ask=ask) + } # end item ii + #************************************************* + +} + +pdf(file = '/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/PDF/icc_difh.pdf') +par(xpd=T,mar=c(6.1,5.1,7.6,2.1)) +icc.tam.difh(zaza0,type = "expected",export=F,package = "graphics",items = 3,axes=F,TT=F) +icc.tam.difh(zaza1,type = "expected",export=F,package = "graphics",items = 3,axes=F,TT=T) +axis(1) +axis(2,las=2) +title(main='A \n \n ', font.main=2) +arrows(x0=-3.65,x1=-4,y0=-0.55,length = 0.15,lwd = 2) +arrows(x0=3.65,x1=4,y0=-0.55,length = 0.15,lwd = 2) +text(x=-3.55,y=-0.55,"Worse\nmental\nhealth",adj=0) +text(x=3.5,y=-0.55,"Better\nmental\nhealth",adj=1) +arrows(x1=-.65,x0=.15,y0=1.5,length = 0.15,lwd = 2) +arrows(x1=-1.7,x0=-0.9,y0=.75,length = 0.15,lwd = 2) +arrows(x1=0.4,x0=1.2,y0=2.25,length = 0.15,lwd = 2) +legend(y=0.5,x=1.5,lty=c(1,1),lwd=c(2,2),col=c('#204776', '#ca346d'),c("Control group","Treatment group")) +title(xlab=expression("Latent variable " * theta * "\n (mental health)")) + +par(xpd=F) +dev.off() + + +# DIF CONVERGENT + +zaz <- read.csv("/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N300/scenario_YA_300.csv") +zaza0 <- tam.mml(resp = zaz[zaz$TT==0,paste0("item",1:4)]) +zaza1 <- tam.mml(resp = zaz[zaz$TT==1,paste0("item",1:4)]) + + +pdf(file = '/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/PDF/icc_difnh1.pdf') +par(xpd=T,mar=c(6.1,5.1,7.6,2.1)) +icc.tam.difh(zaza0,type = "expected",export=F,package = "graphics",items = 3,axes=F,TT=F) +icc.tam.difh(zaza1,type = "expected",export=F,package = "graphics",items = 3,axes=F,TT=T) +axis(1) +axis(2,las=2) +title(main='B \n \n Item: \n "Have you been able to enjoy your normal daily activities?" ', font.main=2) +arrows(x0=-3.65,x1=-4,y0=-0.55,length = 0.15,lwd = 2) +arrows(x0=3.65,x1=4,y0=-0.55,length = 0.15,lwd = 2) +text(x=-3.55,y=-0.55,"Worse\nmental\nhealth",adj=0) +text(x=3.5,y=-0.55,"Better\nmental\nhealth",adj=1) +arrows(x1=-0.1,x0=1.2,y0=2.25,length = 0.15,lwd = 2) +arrows(x1=-.9,x0=.15,y0=1.5,length = 0.15,lwd = 2) +arrows(x1=-1.65,x0=-0.9,y0=.75,length = 0.15,lwd = 2) +legend(y=0.5,x=1.5,lty=c(1,1),lwd=c(2,2),col=c('#204776', '#ca346d'),c("Control group","Treatment group")) +title(xlab=expression("Latent variable " * theta * "\n (mental health)")) + +par(xpd=F) +dev.off() + +# DIF DIVERGENT + +zaz <- read.csv("/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N300/scenario_ZA_300.csv") +zaza0 <- tam.mml(resp = zaz[zaz$TT==0,paste0("item",1:4)]) +zaza1 <- tam.mml(resp = zaz[zaz$TT==1,paste0("item",1:4)]) + + +pdf(file = '/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/PDF/icc_difnh2.pdf') +par(xpd=T,mar=c(6.1,5.1,7.6,2.1)) +icc.tam.difh(zaza0,type = "expected",export=F,package = "graphics",items = 3,axes=F,TT=F) +icc.tam.difh(zaza1,type = "expected",export=F,package = "graphics",items = 3,axes=F,TT=T) +axis(1) +axis(2,las=2) +title(main='C \n \n ', font.main=2) +arrows(x0=-3.65,x1=-4,y0=-0.55,length = 0.15,lwd = 2) +arrows(x0=3.65,x1=4,y0=-0.55,length = 0.15,lwd = 2) +text(x=-3.55,y=-0.55,"Worse\nmental\nhealth",adj=0) +text(x=3.5,y=-0.55,"Better\nmental\nhealth",adj=1) +legend(y=0.5,x=1.5,lty=c(1,1),lwd=c(2,2),col=c('#204776', '#ca346d'),c("Control group","Treatment group")) +arrows(x1=1.9,x0=1.4,y0=2.25,length = 0.15,lwd = 2) +arrows(x1=-1.4,x0=-.9,y0=0.75,length = 0.15,lwd = 2) +title(xlab=expression("Latent variable " * theta * "\n (mental health)")) + +par(xpd=F) +dev.off() + + + +#### ICC DIF BASE + +zaz <- read.csv("/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N300/scenario_XB_300.csv") +zaza <- tam.mml(resp = zaz[zaz$TT==0,paste0("item",1:4)]) +zaza2 <- tam.mml(resp = zaz[zaz$TT==1,paste0("item",1:4)]) +zaza3 <- tam.mml(resp = zaz[,paste0("item",1:4)]) + + +icc.tam <- function(x, items=1:x$nitems, type="expected",TT=F, + low=-3, high=3, ngroups=6, groups_by_item=FALSE, + wle=NULL, export=TRUE, export.type="png", + export.args=list(), observed=TRUE, overlay=FALSE, + ask=FALSE, package="lattice", + fix.devices=TRUE, nnodes=100, ...) +{ + require_namespace_msg("grDevices") + if ( package=="lattice"){ + require_namespace_msg("lattice") + } + + # device.Option <- getOption("device") + + time1 <- NULL + pall <- c('#200c23', '#62403d', '#a87b5e', '#e9bf98' + ) + if ( fix.devices ){ + old.opt.dev <- getOption("device") + old.opt.err <- c( getOption("show.error.messages")) + old.par.ask <- graphics::par("ask") + # remember new pars' values + old.par.xpd <- graphics::par("xpd") + old.par.mar <- graphics::par("mar") + on.exit( options("device"=old.opt.dev)) + on.exit( options("show.error.messages"=old.opt.err), add=TRUE) + on.exit( graphics::par("ask"=old.par.ask), add=TRUE) + # restore new pars' values + on.exit( graphics::par("xpd"=old.par.xpd), add=TRUE) + on.exit( graphics::par("mar"=old.par.mar), add=TRUE) + } + + tamobj <- x + ndim <- tamobj$ndim + tammodel <- "mml" + if(is.null(ndim)) { + ndim <- 1 + tammodel <- "jml" + } + if (ndim > 1 ) { + if ( type=="expected"){ + stop ("Expected scores curves are only available for uni-dimensional models") + } + } + + nitems <- tamobj$nitems + + if (ndim==1 ){ + theta <- matrix(seq(low, high, length=nnodes), nrow=nnodes, ncol=ndim) + } else { + nodes <- seq(-3, 3, length=nnodes) + theta <- as.matrix( expand.grid( as.data.frame( matrix( rep(nodes, ndim), ncol=ndim ) ) ) ) + nnodes <- nrow(theta) + B <- tamobj$B + } + + iIndex <- 1:nitems + A <- tamobj$A + B <- tamobj$B + if (tammodel=="mml") { + xsi <- tamobj$xsi$xsi + } else { + xsi <- tamobj$xsi + } + maxK <- tamobj$maxK + resp <- tamobj$resp + resp.ind <- tamobj$resp.ind + resp[resp.ind==0] <- NA + AXsi <- matrix(0,nrow=nitems,ncol=maxK ) + res <- tam_mml_calc_prob(iIndex=iIndex, A=A, AXsi=AXsi, B=B, xsi=xsi, theta=theta, + nnodes=nnodes, maxK=maxK, recalc=TRUE ) + rprobs <- res[["rprobs"]] + AXsi <- res[["AXsi"]] + cat <- 1:maxK - 1 + + #@@@ define initial empty objects + expScore <- obScore <- wle_intervals <- NULL + theta2 <- NULL + + #**** type='expected' + if ( type=="expected" ){ + expScore <- sapply(1:nitems, function(i) colSums(cat*rprobs[i,,], na.rm=TRUE)) + #-- compute WLE score groups + res <- plot_tam_grouped_wle( tamobj=tamobj, tammodel=tammodel, + wle=wle, ngroups=ngroups, resp=resp ) + wle <- res$wle + theta2 <- res$theta2 + d <- res$d + d1 <- res$d1 + d2 <- res$d2 + groupnumber <- res$groupnumber + ngroups <- res$ngroups + wle_intervals <- res$wle_intervals + #-- compute observed scores + obScore <- apply(d2,2, function(x){ + stats::aggregate(x, list(groupnumber), mean, na.rm=TRUE) + } ) + } + + #---------------------------------------------------- + # adds observed score for type="items" + if (type=="items") { + require_namespace_msg("plyr") + #-- compute WLE score groups + res <- plot_tam_grouped_wle( tamobj=tamobj, tammodel=tammodel, + wle=wle, ngroups=ngroups, resp=resp ) + wle <- res$wle + theta2 <- res$theta2 + d <- res$d + d1 <- res$d1 + d2 <- res$d2 + groupnumber <- res$groupnumber + ngroups <- res$ngroups + + obScore <- lapply(d2, function(item) { + comp_case=stats::complete.cases(item) + item=item[comp_case] + uniq_cats=sort(unique(item)) + plyr::ldply(split(item, groupnumber[comp_case]), .id="group", + function (group) { + ngroup=length(group) + cat_freq=list() + for (catt in uniq_cats) { + cat_freq[[paste0("cat_", catt)]]=sum(group==catt)/ngroup + } + data.frame(cat_freq) + }) + }) + } + + #************************************************* + # begin plot function + probs_plot <- as.list(1:nitems) + names(probs_plot) <- items + + for (i in (1:nitems)[items]) { + #*********************************************************** + #** expected item response curves + if ( type=="expected"){ + if (i==1 || !overlay) { + ylim2 <- c(0,max( tamobj$resp[,i], na.rm=TRUE ) ) + if (TT==0) { + graphics::plot(theta, expScore[,i], type="l", lwd=3, las=1, ylab="Expected item response", + main=NULL,xlab="", + col="#204776",ylim=ylim2,xlim= ... ) + } else { + if (TT==1) { + lines(theta,expScore[,i], type="l", lwd=3, las=1,col="#ca346d") + } else { + lines(theta,expScore[,i], type="l", lwd=2, las=1,col="#999999",lty=2) + } + } + + + } else { + graphics::lines(theta, expScore[,i],type="l", col=i, lwd=3, pch=1) + } + if (observed){ + theta2_i <- theta2 + obScore_i <- obScore[[i]]$x + if (groups_by_item){ + ind_i <- ! is.na(resp[,i]) + resp_i <- resp[ind_i, i, drop=FALSE] + wle_i <- wle[ ind_i ] + res <- plot_tam_grouped_wle( tamobj=tamobj, tammodel=tammodel, + wle=wle_i, ngroups=ngroups, resp=resp_i ) + theta2_i <- res$theta2 + groupnumber_i <- res$groupnumber + aggr <- stats::aggregate(resp_i, list(groupnumber_i), mean, na.rm=TRUE ) + obScore_i <- aggr[,2] + } + #graphics::lines(theta2_i, obScore_i, type="o", lwd=2, pch=1) + } + } + #*********************************************************** + + if ( ndim==1 ){ theta0 <- theta } + + if ( type=="items"){ + rprobs.ii <- rprobs[i,,] + rprobs.ii <- rprobs.ii[ rowMeans( is.na(rprobs.ii) ) < 1, ] + K <- nrow(rprobs.ii) + dat2 <- NULL + #************ + if ( ndim > 1 ){ + B.ii <- B[i,,] + ind.ii <- which( colSums( B.ii ) > 0 )[1] + rprobs0.ii <- rprobs.ii + rprobs0.ii <- stats::aggregate( t(rprobs0.ii), list( theta[,ind.ii] ), mean ) + theta0 <- rprobs0.ii[,1,drop=FALSE] + rprobs.ii <- t( rprobs0.ii[,-1] ) + } + probs_plot[[i]] <- rprobs.ii + #************** + for (kk in 1:K){ + dat2a <- data.frame( "Theta"=theta0[,1], "cat"=kk, "P"=rprobs.ii[kk,] ) + dat2 <- rbind(dat2, dat2a) + } + auto.key <- NULL + simple.key <- paste0("Cat", 1:K - 1) + auto.key <- simple.key + dat2$time <- dat2$cat + dat2$time1 <- paste0("Cat", dat2$time ) + + simple.key <- FALSE + Kpercol <- K + # package graphics + if ( package=="graphics" ){ + kk <- 1 + dfr <- dat2 + dfr1a <- dfr[ dfr$cat==kk, ] + graphics::plot( dfr1a$Theta, dfr1a$P, ylim=c(-.1,1.1), + xlab=expression(theta), + col=pall[kk], type="l", xpd=TRUE,axes=F, ... + ) + axis(1) + axis(2) + grid(nx = NA, + ny = NULL, + lty = 3, col = "lightgray", lwd = 1) + graphics::lines( dfr1a$Theta, dfr1a$P, + col=pall[kk], type="l", ... + ) + for (kk in seq(2,K) ){ + dfr1a <- dfr[ dfr$cat==kk, ] + graphics::lines( dfr1a$Theta, dfr1a$P, col=pall[kk] ) + # graphics::points( dfr1a$Theta, dfr1a$P, pch=kk, col=kk+1 ) + + } + } + + + #*************************************** + + } + #*************** + graphics::par(ask=ask) + } # end item ii + #************************************************* + +} + +png(file = '/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/PNG/icc_beta.png',width=800,height=800) +par(xpd=T,mar=c(6.1,5.1,7.6,2.1)) +icc.tam(zaza,type = "expected",export=F,package = "graphics",items = 3,axes=F,TT=0) +icc.tam(zaza2,type = "expected",export=F,package = "graphics",items = 3,axes=F,TT=1) +icc.tam(zaza3,type = "expected",export=F,package = "graphics",items = 3,axes=F,TT=2) +axis(1) +axis(2,las=2) +title(main='A - Item Characteristic Curve \n Item: \n "Have you been able to enjoy your normal daily activities?" ', font.main=2) +arrows(x0=-2.65,x1=-3,y0=-0.45,length = 0.15,lwd = 2) +arrows(x0=2.65,x1=3,y0=-0.45,length = 0.15,lwd = 2) +text(x=-2.55,y=-0.45,"Worse\nmental\nhealth",adj=0) +text(x=2.5,y=-0.45,"Better\nmental\nhealth",adj=1) +arrows(x0=0.175,x1=-0.175,y0=1.5,length = 0.15,lwd = 2,col="black") +text(x=0.05,y=1.575,"DIF",col="black") +legend(y=3,x=-3,lty=c(1,1,2),lwd=c(2,2,2),col=c('#204776', '#ca346d',"#999999"),c("Control group","Treatment group","True ICC")) +title(xlab=expression("Latent variable " * theta * "\n (mental health)")) + +par(xpd=F) +dev.off() + +icc.tam.2 <- function(x, items=1:x$nitems, type="expected",TT=F, + low=-3, high=3, ngroups=6, groups_by_item=FALSE, + wle=NULL, export=TRUE, export.type="png", + export.args=list(), observed=TRUE, overlay=FALSE, + ask=FALSE, package="lattice", + fix.devices=TRUE, nnodes=100, ...) +{ + require_namespace_msg("grDevices") + if ( package=="lattice"){ + require_namespace_msg("lattice") + } + + # device.Option <- getOption("device") + + time1 <- NULL + pall <- c('#204776', '#ca346d' + ) + if ( fix.devices ){ + old.opt.dev <- getOption("device") + old.opt.err <- c( getOption("show.error.messages")) + old.par.ask <- graphics::par("ask") + # remember new pars' values + old.par.xpd <- graphics::par("xpd") + old.par.mar <- graphics::par("mar") + on.exit( options("device"=old.opt.dev)) + on.exit( options("show.error.messages"=old.opt.err), add=TRUE) + on.exit( graphics::par("ask"=old.par.ask), add=TRUE) + # restore new pars' values + on.exit( graphics::par("xpd"=old.par.xpd), add=TRUE) + on.exit( graphics::par("mar"=old.par.mar), add=TRUE) + } + + tamobj <- x + ndim <- tamobj$ndim + tammodel <- "mml" + if(is.null(ndim)) { + ndim <- 1 + tammodel <- "jml" + } + if (ndim > 1 ) { + if ( type=="expected"){ + stop ("Expected scores curves are only available for uni-dimensional models") + } + } + + nitems <- tamobj$nitems + + if (ndim==1 ){ + theta <- matrix(seq(low, high, length=nnodes), nrow=nnodes, ncol=ndim) + } else { + nodes <- seq(-3, 3, length=nnodes) + theta <- as.matrix( expand.grid( as.data.frame( matrix( rep(nodes, ndim), ncol=ndim ) ) ) ) + nnodes <- nrow(theta) + B <- tamobj$B + } + + iIndex <- 1:nitems + A <- tamobj$A + B <- tamobj$B + if (tammodel=="mml") { + xsi <- tamobj$xsi$xsi + } else { + xsi <- tamobj$xsi + } + maxK <- tamobj$maxK + resp <- tamobj$resp + resp.ind <- tamobj$resp.ind + resp[resp.ind==0] <- NA + AXsi <- matrix(0,nrow=nitems,ncol=maxK ) + res <- tam_mml_calc_prob(iIndex=iIndex, A=A, AXsi=AXsi, B=B, xsi=xsi, theta=theta, + nnodes=nnodes, maxK=maxK, recalc=TRUE ) + rprobs <- res[["rprobs"]] + AXsi <- res[["AXsi"]] + cat <- 1:maxK - 1 + + #@@@ define initial empty objects + expScore <- obScore <- wle_intervals <- NULL + theta2 <- NULL + + #**** type='expected' + if ( type=="expected" ){ + expScore <- sapply(1:nitems, function(i) colSums(cat*rprobs[i,,], na.rm=TRUE)) + #-- compute WLE score groups + res <- plot_tam_grouped_wle( tamobj=tamobj, tammodel=tammodel, + wle=wle, ngroups=ngroups, resp=resp ) + wle <- res$wle + theta2 <- res$theta2 + d <- res$d + d1 <- res$d1 + d2 <- res$d2 + groupnumber <- res$groupnumber + ngroups <- res$ngroups + wle_intervals <- res$wle_intervals + #-- compute observed scores + obScore <- apply(d2,2, function(x){ + stats::aggregate(x, list(groupnumber), mean, na.rm=TRUE) + } ) + } + + #---------------------------------------------------- + # adds observed score for type="items" + if (type=="items") { + require_namespace_msg("plyr") + #-- compute WLE score groups + res <- plot_tam_grouped_wle( tamobj=tamobj, tammodel=tammodel, + wle=wle, ngroups=ngroups, resp=resp ) + wle <- res$wle + theta2 <- res$theta2 + d <- res$d + d1 <- res$d1 + d2 <- res$d2 + groupnumber <- res$groupnumber + ngroups <- res$ngroups + + obScore <- lapply(d2, function(item) { + comp_case=stats::complete.cases(item) + item=item[comp_case] + uniq_cats=sort(unique(item)) + plyr::ldply(split(item, groupnumber[comp_case]), .id="group", + function (group) { + ngroup=length(group) + cat_freq=list() + for (catt in uniq_cats) { + cat_freq[[paste0("cat_", catt)]]=sum(group==catt)/ngroup + } + data.frame(cat_freq) + }) + }) + } + + #************************************************* + # begin plot function + probs_plot <- as.list(1:nitems) + names(probs_plot) <- items + + for (i in (1:nitems)[items]) { + #*********************************************************** + #** expected item response curves + if ( type=="expected"){ + if (i==1 || !overlay) { + ylim2 <- c(0,max( tamobj$resp[,i], na.rm=TRUE ) ) + if (TT==0) { + graphics::plot(theta, expScore[,i], type="l", lwd=3, las=1, ylab="Expected item response", + main=NULL,xlab="", + col="#ca346d",ylim=ylim2,xlim= ... ) + } else { + if (TT==1) { + lines(theta,expScore[,i], type="l", lwd=3, las=1,col="#204776") + } else { + lines(theta,expScore[,i], type="l", lwd=2, las=1,col="#999999",lty=2) + } + } + + + } else { + graphics::lines(theta, expScore[,i],type="l", col=i, lwd=3, pch=1) + } + if (observed){ + theta2_i <- theta2 + obScore_i <- obScore[[i]]$x + if (groups_by_item){ + ind_i <- ! is.na(resp[,i]) + resp_i <- resp[ind_i, i, drop=FALSE] + wle_i <- wle[ ind_i ] + res <- plot_tam_grouped_wle( tamobj=tamobj, tammodel=tammodel, + wle=wle_i, ngroups=ngroups, resp=resp_i ) + theta2_i <- res$theta2 + groupnumber_i <- res$groupnumber + aggr <- stats::aggregate(resp_i, list(groupnumber_i), mean, na.rm=TRUE ) + obScore_i <- aggr[,2] + } + #graphics::lines(theta2_i, obScore_i, type="o", lwd=2, pch=1) + } + } + #*********************************************************** + + if ( ndim==1 ){ theta0 <- theta } + + if ( type=="items"){ + rprobs.ii <- rprobs[i,,] + rprobs.ii <- rprobs.ii[ rowMeans( is.na(rprobs.ii) ) < 1, ] + K <- nrow(rprobs.ii) + dat2 <- NULL + #************ + if ( ndim > 1 ){ + B.ii <- B[i,,] + ind.ii <- which( colSums( B.ii ) > 0 )[1] + rprobs0.ii <- rprobs.ii + rprobs0.ii <- stats::aggregate( t(rprobs0.ii), list( theta[,ind.ii] ), mean ) + theta0 <- rprobs0.ii[,1,drop=FALSE] + rprobs.ii <- t( rprobs0.ii[,-1] ) + } + probs_plot[[i]] <- rprobs.ii + #************** + for (kk in 1:K){ + dat2a <- data.frame( "Theta"=theta0[,1], "cat"=kk, "P"=rprobs.ii[kk,] ) + dat2 <- rbind(dat2, dat2a) + } + auto.key <- NULL + simple.key <- paste0("Cat", 1:K - 1) + auto.key <- simple.key + dat2$time <- dat2$cat + dat2$time1 <- paste0("Cat", dat2$time ) + + simple.key <- FALSE + Kpercol <- K + # package graphics + if ( package=="graphics" ){ + kk <- 1 + dfr <- dat2 + dfr1a <- dfr[ dfr$cat==kk, ] + graphics::plot( dfr1a$Theta, dfr1a$P, ylim=c(-.1,1.1), + xlab=expression(theta), + col=pall[kk], type="l", xpd=TRUE,axes=F, ... + ) + axis(1) + axis(2) + grid(nx = NA, + ny = NULL, + lty = 3, col = "lightgray", lwd = 1) + graphics::lines( dfr1a$Theta, dfr1a$P, + col=pall[kk], type="l", ... + ) + for (kk in seq(2,K) ){ + dfr1a <- dfr[ dfr$cat==kk, ] + graphics::lines( dfr1a$Theta, dfr1a$P, col=pall[kk] ) + # graphics::points( dfr1a$Theta, dfr1a$P, pch=kk, col=kk+1 ) + + } + } + + + #*************************************** + + } + #*************** + graphics::par(ask=ask) + } # end item ii + #************************************************* + +} + +png(file = '/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/PNG/icc_beta2.png',width=800,height=800) +par(xpd=T,mar=c(6.1,5.1,7.6,2.1)) +icc.tam.2(zaza,type = "expected",export=F,package = "graphics",items = 3,axes=F,TT=0) +icc.tam.2(zaza2,type = "expected",export=F,package = "graphics",items = 3,axes=F,TT=1) +icc.tam.2(zaza3,type = "expected",export=F,package = "graphics",items = 3,axes=F,TT=2) +axis(1) +axis(2,las=2) +title(main='A - Item Characteristic Curve \n Item: \n "Have you been able to enjoy your normal daily activities?" ', font.main=2) +arrows(x0=-2.65,x1=-3,y0=-0.45,length = 0.15,lwd = 2) +arrows(x0=2.65,x1=3,y0=-0.45,length = 0.15,lwd = 2) +text(x=-2.55,y=-0.45,"Worse\nmental\nhealth",adj=0) +text(x=2.5,y=-0.45,"Better\nmental\nhealth",adj=1) +arrows(x1=0.175,x0=-0.175,y0=1.5,length = 0.15,lwd = 2,col="black") +text(x=0.05,y=1.575,"DIF",col="black") +legend(y=3,x=-3,lty=c(1,1,2),lwd=c(2,2,2),col=c('#204776', '#ca346d',"#999999"),c("Control group","Treatment group","True ICC")) +title(xlab=expression("Latent variable " * theta * "\n (mental health)")) +title(xlab=expression("Latent variable " * theta * "\n (mental health)")) + +par(xpd=F) +dev.off() + + +#### Densité theta beta0 + + +png(file = '/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/PNG/densite_theta_null.png',width=800,height=800) +par(xpd=T,mar=c(6.1,5.1,7.6,2.1)) +plot(seq(-3,3,0.001),dnorm(seq(-3,3,0.001),mean = 0,sd=1),type="l",lwd=2,col="#204776",axes=F, + xlab="",ylab="Probability density") +axis(1) +axis(2,las=2) +title(main='B - true vs. observed latent variable density \n Item: \n "Have you been able to enjoy your normal daily activities?" ', font.main=2,line=2.8) +lines(seq(-3,3,0.001),dnorm(seq(-3,3,0.001),mean = 0,sd=1),type="l",lwd=0.5,col="#ca346d",lty=1) +text(x=0,y=0.4225,expression(beta==0),adj=0.5) +legend(y=0.38,x=-3,lty=c(1,1,2,2),lwd=c(2,2,2,2),cex=1,col=c('#204776', '#ca346d','#204776', '#ca346d'),c("Control group (true)","Treatment group (true)","Control group (observed)","Treatment group (observed)")) +lines(seq(-3,3,0.001),dnorm(seq(-3,3,0.001),mean = 1,sd=1),type="l",lwd=2,col="#ca346d",lty=2) +arrows(x0=0,x1=1,y0=.4125,length = 0.1,lwd = 2,col="black") +text(x=0.5,y=0.425,expression(widehat(beta)>0),adj=0.5,col="black") +segments(x0=0,x1=0,y0=0.4,y1=0.4125,lty=3,col="black") +segments(x0=1,x1=1,y0=0.4,y1=0.4125,lty=3,col="black") +arrows(x0=-2.65,x1=-3,y0=-0.05,length = 0.15,lwd = 2) +arrows(x0=2.65,x1=3,y0=-0.05,length = 0.15,lwd = 2) +text(x=-2.55,y=-0.05,"Worse\nmental\nhealth",adj=0) +text(x=2.5,y=-0.05,"Better\nmental\nhealth",adj=1) +title(xlab=expression("Latent variable " * theta * "\n (mental health)")) +par(xpd=F) +dev.off() + + +#### Densité theta mask + + +png(file = '/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/PNG/densite_theta_masks.png',width=800,height=800) +par(xpd=T,mar=c(6.1,5.1,7.6,2.1)) +plot(seq(-3,3,0.001),dnorm(seq(-3,3,0.001),mean = 0,sd=1),type="l",lwd=2,col="#204776",axes=F, + xlab="",ylab="Probability density") +axis(1) +axis(2,las=2) +title(main='B - true vs. observed latent variable density \n Item: \n "Have you been able to enjoy your normal daily activities?" ', font.main=2,line = 3.4) +lines(seq(-3,3,0.001),dnorm(seq(-3,3,0.001),mean = 1,sd=1),type="l",lwd=2,col="#ca346d",lty=1) +lines(seq(-3,3,0.001),dnorm(seq(-3,3,0.001),mean = 0.5,sd=1),type="l",lwd=2,col="#ca346d",lty=2) +arrows(x0=0,x1=1,y0=.4325,length = 0.1,lwd = 2,col="black") +text(x=0.5,y=0.4425,expression(beta),adj=0.5) +segments(x0=0,x1=0,y0=0.4,y1=0.4325,lty=3) +segments(x0=1,x1=1,y0=0.4,y1=0.4325,lty=3) + +arrows(x0=1,x1=0.51,y0=.415,length = 0.1,lwd = 2,col="black") +text(x=0.75,y=0.4225,"DIF",adj=0.5,col="black") +segments(x0=0.5,x1=0.5,y0=0.4,y1=.415,lty=3,col="black") + +arrows(x0=0,x1=0.49,y0=.405,length = 0.1,lwd = 2,col="black") +text(x=0.25,y=0.4175,expression(widehat(beta)),adj=0.5,col="black") +segments(x0=0,x1=0,y0=0.4,y1=.405,lty=3,col="black") +segments(x0=0.5,x1=0.5,y0=0.4,y1=.405,lty=3,col="black") +title(xlab=expression("Latent variable " * theta * "\n (mental health)")) + + + +legend(y=0.38,x=-3,lty=c(1,1,2,2),lwd=c(2,2,2,2),cex=1,col=c('#204776', '#ca346d','#204776', '#ca346d'),c("Control group (true)","Treatment group (true)","Control group (observed)","Treatment group (observed)")) +arrows(x0=-2.65,x1=-3,y0=-0.05,length = 0.15,lwd = 2) +arrows(x0=2.65,x1=3,y0=-0.05,length = 0.15,lwd = 2) +text(x=-2.55,y=-0.05,"Worse\nmental\nhealth",adj=0) +text(x=2.5,y=-0.05,"Better\nmental\nhealth",adj=1) +par(xpd=F) +dev.off() + + +#### Densité theta amplify + + +png(file = '/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/PNG/densite_theta_amplify.png',width=800,height=800) +par(xpd=T,mar=c(6.1,5.1,7.6,2.1)) +plot(seq(-3,3,0.001),dnorm(seq(-3,3,0.001),mean = 0,sd=1),type="l",lwd=2,col="#204776",axes=F, + xlab="",ylab="Probability density") +title(xlab=expression("Latent variable " * theta * "\n (mental health)")) +axis(1) +axis(2,las=2) +title(main='B - true vs. observed latent variable density \n Item: \n "Have you been able to enjoy your normal daily activities?" ', font.main=2,line = 3.4) +lines(seq(-3,3,0.001),dnorm(seq(-3,3,0.001),mean = 1,sd=1),type="l",lwd=2,col="#ca346d",lty=1) +lines(seq(-3,3,0.001),dnorm(seq(-3,3,0.001),mean = 1.5,sd=1),type="l",lwd=2,col="#ca346d",lty=2) +arrows(x0=0,x1=0.98,y0=.4325,length = 0.1,lwd = 2,col="black") +text(x=0.5,y=0.4425,expression(beta),adj=0.5) +segments(x0=0,x1=0,y0=0.4,y1=0.4325,lty=3) +segments(x0=1,x1=1,y0=0.4,y1=0.4325,lty=3) + +arrows(x0=1.02,x1=1.5,y0=.4325,length = 0.1,lwd = 2,col="black") +text(x=1.25,y=0.4425,"DIF",adj=0.5,col="black") +segments(x0=1.5,x1=1.5,y0=0.4,y1=.4325,lty=3,col="black") + +arrows(x0=0,x1=1.5,y0=.41,length = 0.1,lwd = 2,col="black") +text(x=0.75,y=0.4215,expression(widehat(beta)),adj=0.5,col="black") +segments(x0=0,x1=0,y0=0.4,y1=.41,lty=3,col="black") +segments(x0=1.5,x1=1.5,y0=0.4,y1=.41,lty=3,col="black") + +legend(y=0.38,x=-3,lty=c(1,1,2,2),lwd=c(2,2,2,2),cex=1,col=c('#204776', '#ca346d','#204776', '#ca346d'),c("Control group (true)","Treatment group (true)","Control group (observed)","Treatment group (observed)")) +arrows(x0=-2.65,x1=-3,y0=-0.05,length = 0.15,lwd = 2) +arrows(x0=2.65,x1=3,y0=-0.05,length = 0.15,lwd = 2) +text(x=-2.55,y=-0.05,"Worse\nmental\nhealth",adj=0) +text(x=2.5,y=-0.05,"Better\nmental\nhealth",adj=1) +par(xpd=F) +dev.off() + +#### ICC DIF SUPPL RESIDIF 1 + +zaz <- read.csv("/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N300/scenario_XB_300.csv") +zaza <- tam.mml(resp = zaz[zaz$TT==0,paste0("item",1:4)]) +zaza2 <- tam.mml(resp = zaz[zaz$TT==1,paste0("item",1:4)]) +zaza3 <- tam.mml(resp = zaz[,paste0("item",1:4)]) + + +icc.tam.sup <- function(x, items=1:x$nitems, type="expected",TT=F, + low=-3, high=3, ngroups=6, groups_by_item=FALSE, + wle=NULL, export=TRUE, export.type="png", + export.args=list(), observed=TRUE, overlay=FALSE, + ask=FALSE, package="lattice", + fix.devices=TRUE, nnodes=100, ...) +{ + require_namespace_msg("grDevices") + if ( package=="lattice"){ + require_namespace_msg("lattice") + } + + # device.Option <- getOption("device") + + time1 <- NULL + pall <- c('#200c23', '#62403d', '#a87b5e', '#e9bf98' + ) + if ( fix.devices ){ + old.opt.dev <- getOption("device") + old.opt.err <- c( getOption("show.error.messages")) + old.par.ask <- graphics::par("ask") + # remember new pars' values + old.par.xpd <- graphics::par("xpd") + old.par.mar <- graphics::par("mar") + on.exit( options("device"=old.opt.dev)) + on.exit( options("show.error.messages"=old.opt.err), add=TRUE) + on.exit( graphics::par("ask"=old.par.ask), add=TRUE) + # restore new pars' values + on.exit( graphics::par("xpd"=old.par.xpd), add=TRUE) + on.exit( graphics::par("mar"=old.par.mar), add=TRUE) + } + + tamobj <- x + ndim <- tamobj$ndim + tammodel <- "mml" + if(is.null(ndim)) { + ndim <- 1 + tammodel <- "jml" + } + if (ndim > 1 ) { + if ( type=="expected"){ + stop ("Expected scores curves are only available for uni-dimensional models") + } + } + + nitems <- tamobj$nitems + + if (ndim==1 ){ + theta <- matrix(seq(low, high, length=nnodes), nrow=nnodes, ncol=ndim) + } else { + nodes <- seq(-3, 3, length=nnodes) + theta <- as.matrix( expand.grid( as.data.frame( matrix( rep(nodes, ndim), ncol=ndim ) ) ) ) + nnodes <- nrow(theta) + B <- tamobj$B + } + + iIndex <- 1:nitems + A <- tamobj$A + B <- tamobj$B + if (tammodel=="mml") { + xsi <- tamobj$xsi$xsi + } else { + xsi <- tamobj$xsi + } + maxK <- tamobj$maxK + resp <- tamobj$resp + resp.ind <- tamobj$resp.ind + resp[resp.ind==0] <- NA + AXsi <- matrix(0,nrow=nitems,ncol=maxK ) + res <- tam_mml_calc_prob(iIndex=iIndex, A=A, AXsi=AXsi, B=B, xsi=xsi, theta=theta, + nnodes=nnodes, maxK=maxK, recalc=TRUE ) + rprobs <- res[["rprobs"]] + AXsi <- res[["AXsi"]] + cat <- 1:maxK - 1 + + #@@@ define initial empty objects + expScore <- obScore <- wle_intervals <- NULL + theta2 <- NULL + + #**** type='expected' + if ( type=="expected" ){ + expScore <- sapply(1:nitems, function(i) colSums(cat*rprobs[i,,], na.rm=TRUE)) + #-- compute WLE score groups + res <- plot_tam_grouped_wle( tamobj=tamobj, tammodel=tammodel, + wle=wle, ngroups=ngroups, resp=resp ) + wle <- res$wle + theta2 <- res$theta2 + d <- res$d + d1 <- res$d1 + d2 <- res$d2 + groupnumber <- res$groupnumber + ngroups <- res$ngroups + wle_intervals <- res$wle_intervals + #-- compute observed scores + obScore <- apply(d2,2, function(x){ + stats::aggregate(x, list(groupnumber), mean, na.rm=TRUE) + } ) + } + + #---------------------------------------------------- + # adds observed score for type="items" + if (type=="items") { + require_namespace_msg("plyr") + #-- compute WLE score groups + res <- plot_tam_grouped_wle( tamobj=tamobj, tammodel=tammodel, + wle=wle, ngroups=ngroups, resp=resp ) + wle <- res$wle + theta2 <- res$theta2 + d <- res$d + d1 <- res$d1 + d2 <- res$d2 + groupnumber <- res$groupnumber + ngroups <- res$ngroups + + obScore <- lapply(d2, function(item) { + comp_case=stats::complete.cases(item) + item=item[comp_case] + uniq_cats=sort(unique(item)) + plyr::ldply(split(item, groupnumber[comp_case]), .id="group", + function (group) { + ngroup=length(group) + cat_freq=list() + for (catt in uniq_cats) { + cat_freq[[paste0("cat_", catt)]]=sum(group==catt)/ngroup + } + data.frame(cat_freq) + }) + }) + } + + #************************************************* + # begin plot function + probs_plot <- as.list(1:nitems) + names(probs_plot) <- items + + for (i in (1:nitems)[items]) { + #*********************************************************** + #** expected item response curves + if ( type=="expected"){ + if (i==1 || !overlay) { + ylim2 <- c(0,max( tamobj$resp[,i], na.rm=TRUE ) ) + if (TT==0) { + graphics::plot(theta, expScore[,i]-0.2, type="l", lwd=3, las=1, ylab="Expected item response", + main=NULL,xlab="", + col="#204776",ylim=ylim2,xlim= ... ) + } else { + if (TT==1) { + lines(theta,expScore[,i]+0.2, type="l", lwd=3, las=1,col="#ca346d") + } else { + lines(theta,expScore[,i], type="l", lwd=2, las=1,col="#999999",lty=2) + } + } + + + } else { + graphics::lines(theta, expScore[,i],type="l", col=i, lwd=3, pch=1) + } + if (observed){ + theta2_i <- theta2 + obScore_i <- obScore[[i]]$x + if (groups_by_item){ + ind_i <- ! is.na(resp[,i]) + resp_i <- resp[ind_i, i, drop=FALSE] + wle_i <- wle[ ind_i ] + res <- plot_tam_grouped_wle( tamobj=tamobj, tammodel=tammodel, + wle=wle_i, ngroups=ngroups, resp=resp_i ) + theta2_i <- res$theta2 + groupnumber_i <- res$groupnumber + aggr <- stats::aggregate(resp_i, list(groupnumber_i), mean, na.rm=TRUE ) + obScore_i <- aggr[,2] + } + #graphics::lines(theta2_i, obScore_i, type="o", lwd=2, pch=1) + } + } + #*********************************************************** + + if ( ndim==1 ){ theta0 <- theta } + + if ( type=="items"){ + rprobs.ii <- rprobs[i,,] + rprobs.ii <- rprobs.ii[ rowMeans( is.na(rprobs.ii) ) < 1, ] + K <- nrow(rprobs.ii) + dat2 <- NULL + #************ + if ( ndim > 1 ){ + B.ii <- B[i,,] + ind.ii <- which( colSums( B.ii ) > 0 )[1] + rprobs0.ii <- rprobs.ii + rprobs0.ii <- stats::aggregate( t(rprobs0.ii), list( theta[,ind.ii] ), mean ) + theta0 <- rprobs0.ii[,1,drop=FALSE] + rprobs.ii <- t( rprobs0.ii[,-1] ) + } + probs_plot[[i]] <- rprobs.ii + #************** + for (kk in 1:K){ + dat2a <- data.frame( "Theta"=theta0[,1], "cat"=kk, "P"=rprobs.ii[kk,] ) + dat2 <- rbind(dat2, dat2a) + } + auto.key <- NULL + simple.key <- paste0("Cat", 1:K - 1) + auto.key <- simple.key + dat2$time <- dat2$cat + dat2$time1 <- paste0("Cat", dat2$time ) + + simple.key <- FALSE + Kpercol <- K + # package graphics + if ( package=="graphics" ){ + kk <- 1 + dfr <- dat2 + dfr1a <- dfr[ dfr$cat==kk, ] + graphics::plot( dfr1a$Theta, dfr1a$P, ylim=c(-.1,1.1), + xlab=expression(theta), + col=pall[kk], type="l", xpd=TRUE,axes=F, ... + ) + axis(1) + axis(2) + grid(nx = NA, + ny = NULL, + lty = 3, col = "lightgray", lwd = 1) + graphics::lines( dfr1a$Theta, dfr1a$P, + col=pall[kk], type="l", ... + ) + for (kk in seq(2,K) ){ + dfr1a <- dfr[ dfr$cat==kk, ] + graphics::lines( dfr1a$Theta, dfr1a$P, col=pall[kk] ) + # graphics::points( dfr1a$Theta, dfr1a$P, pch=kk, col=kk+1 ) + + } + } + + + #*************************************** + + } + #*************** + graphics::par(ask=ask) + } # end item ii + #************************************************* + +} + +png(file = '/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/PNG/icc_residif_1.png',width=800,height=800) +par(xpd=T,mar=c(6.1,5.1,7.6,2.1)) +icc.tam.sup(zaza,type = "expected",export=F,package = "graphics",items = 3,axes=F,TT=0) +icc.tam.sup(zaza2,type = "expected",export=F,package = "graphics",items = 3,axes=F,TT=1) +icc.tam.sup(zaza3,type = "expected",export=F,package = "graphics",items = 3,axes=F,TT=2) +axis(1) +axis(2,las=2) +title(main='Item: \n "Have you been able to enjoy your normal daily activities?" ', font.main=2) +arrows(x0=-2.65,x1=-3,y0=-0.45,length = 0.15,lwd = 2) +arrows(x0=2.65,x1=3,y0=-0.45,length = 0.15,lwd = 2) +arrows(x0=-2,x1=-2,y1=0.375,y0=0.575,length = 0.15,lwd = 2,col="#ca346d") +arrows(x0=0,x1=0,y1=1.6,y0=1.8,length = 0.15,lwd = 2,col="#ca346d") +arrows(x0=2,x1=2,y1=2.7,y0=2.9,length = 0.15,lwd = 2,col="#ca346d") + +text(x=-2,y=0.725,col="#ca346d","-",cex=1.5) +text(x=0,y=2.0,col="#ca346d","-",cex=1.5) +text(x=2,y=3,col="#ca346d","-",cex=1.5) +text(x=-2,y=0.0,col="#204776","+",cex=1.5) +text(x=0,y=0.975,col="#204776","+",cex=1.5) +text(x=2,y=2.3,col="#204776","+",cex=1.5) + +arrows(x0=-2,x1=-2,y1=0.305,y0=0.105,length = 0.15,lwd = 2,col="#204776") +arrows(x0=0,x1=0,y1=1.4,y0=1.2,length = 0.15,lwd = 2,col="#204776") +arrows(x0=2,x1=2,y1=2.625,y0=2.425,length = 0.15,lwd = 2,col="#204776") +text(x=-2.55,y=-0.45,"Worse\nmental\nhealth",adj=0) +text(x=2.5,y=-0.45,"Better\nmental\nhealth",adj=1) + +legend(cex=0.9,y=3,x=-3,lty=c(1,1,2,NA,NA),lwd=c(2,2,2,NA,NA),pch = c(NA,NA,NA,NA,NA),col=c('#204776', '#ca346d',"#999999",NA,NA),c("Control group ICC","Treatment group ICC","Overall ICC","Person-item residuals - Control","Person-item residuals - Treatment")) +par(font=5) +legend(cex=0.9,y=3,x=-3,lty=c(NA,NA,NA,NA,NA),lwd=c(NA,NA,NA,NA,NA),pt.cex=c(NA,NA,NA,1.5,1.5),pch = c(NA,NA,NA,174,174),col=c(NA,NA,NA,'#204776', '#ca346d'),c(NA,NA,NA,NA,NA),bty="n") +par(font=1) +title(xlab=expression("Latent variable " * theta * "\n (mental health)")) +par(xpd=F) +dev.off() + + + + +#### ICC DIF SUPPL RESIDIF 2 + +zaz <- read.csv("/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N300/scenario_YA_300.csv") +zaza <- tam.mml(resp = zaz[zaz$TT==0,paste0("item",1:4)]) +zaza2 <- tam.mml(resp = zaz[zaz$TT==1,paste0("item",1:4)]) +zaza3 <- tam.mml(resp = zaz[,paste0("item",1:4)]) + + +icc.tam.sup <- function(x, items=1:x$nitems, type="expected",TT=F, + low=-3, high=3, ngroups=6, groups_by_item=FALSE, + wle=NULL, export=TRUE, export.type="png", + export.args=list(), observed=TRUE, overlay=FALSE, + ask=FALSE, package="lattice", + fix.devices=TRUE, nnodes=100, ...) +{ + require_namespace_msg("grDevices") + if ( package=="lattice"){ + require_namespace_msg("lattice") + } + + # device.Option <- getOption("device") + + time1 <- NULL + pall <- c('#200c23', '#62403d', '#a87b5e', '#e9bf98' + ) + if ( fix.devices ){ + old.opt.dev <- getOption("device") + old.opt.err <- c( getOption("show.error.messages")) + old.par.ask <- graphics::par("ask") + # remember new pars' values + old.par.xpd <- graphics::par("xpd") + old.par.mar <- graphics::par("mar") + on.exit( options("device"=old.opt.dev)) + on.exit( options("show.error.messages"=old.opt.err), add=TRUE) + on.exit( graphics::par("ask"=old.par.ask), add=TRUE) + # restore new pars' values + on.exit( graphics::par("xpd"=old.par.xpd), add=TRUE) + on.exit( graphics::par("mar"=old.par.mar), add=TRUE) + } + + tamobj <- x + ndim <- tamobj$ndim + tammodel <- "mml" + if(is.null(ndim)) { + ndim <- 1 + tammodel <- "jml" + } + if (ndim > 1 ) { + if ( type=="expected"){ + stop ("Expected scores curves are only available for uni-dimensional models") + } + } + + nitems <- tamobj$nitems + + if (ndim==1 ){ + theta <- matrix(seq(low, high, length=nnodes), nrow=nnodes, ncol=ndim) + } else { + nodes <- seq(-3, 3, length=nnodes) + theta <- as.matrix( expand.grid( as.data.frame( matrix( rep(nodes, ndim), ncol=ndim ) ) ) ) + nnodes <- nrow(theta) + B <- tamobj$B + } + + iIndex <- 1:nitems + A <- tamobj$A + B <- tamobj$B + if (tammodel=="mml") { + xsi <- tamobj$xsi$xsi + } else { + xsi <- tamobj$xsi + } + maxK <- tamobj$maxK + resp <- tamobj$resp + resp.ind <- tamobj$resp.ind + resp[resp.ind==0] <- NA + AXsi <- matrix(0,nrow=nitems,ncol=maxK ) + res <- tam_mml_calc_prob(iIndex=iIndex, A=A, AXsi=AXsi, B=B, xsi=xsi, theta=theta, + nnodes=nnodes, maxK=maxK, recalc=TRUE ) + rprobs <- res[["rprobs"]] + AXsi <- res[["AXsi"]] + cat <- 1:maxK - 1 + + #@@@ define initial empty objects + expScore <- obScore <- wle_intervals <- NULL + theta2 <- NULL + + #**** type='expected' + if ( type=="expected" ){ + expScore <- sapply(1:nitems, function(i) colSums(cat*rprobs[i,,], na.rm=TRUE)) + #-- compute WLE score groups + res <- plot_tam_grouped_wle( tamobj=tamobj, tammodel=tammodel, + wle=wle, ngroups=ngroups, resp=resp ) + wle <- res$wle + theta2 <- res$theta2 + d <- res$d + d1 <- res$d1 + d2 <- res$d2 + groupnumber <- res$groupnumber + ngroups <- res$ngroups + wle_intervals <- res$wle_intervals + #-- compute observed scores + obScore <- apply(d2,2, function(x){ + stats::aggregate(x, list(groupnumber), mean, na.rm=TRUE) + } ) + } + + #---------------------------------------------------- + # adds observed score for type="items" + if (type=="items") { + require_namespace_msg("plyr") + #-- compute WLE score groups + res <- plot_tam_grouped_wle( tamobj=tamobj, tammodel=tammodel, + wle=wle, ngroups=ngroups, resp=resp ) + wle <- res$wle + theta2 <- res$theta2 + d <- res$d + d1 <- res$d1 + d2 <- res$d2 + groupnumber <- res$groupnumber + ngroups <- res$ngroups + + obScore <- lapply(d2, function(item) { + comp_case=stats::complete.cases(item) + item=item[comp_case] + uniq_cats=sort(unique(item)) + plyr::ldply(split(item, groupnumber[comp_case]), .id="group", + function (group) { + ngroup=length(group) + cat_freq=list() + for (catt in uniq_cats) { + cat_freq[[paste0("cat_", catt)]]=sum(group==catt)/ngroup + } + data.frame(cat_freq) + }) + }) + } + + #************************************************* + # begin plot function + probs_plot <- as.list(1:nitems) + names(probs_plot) <- items + + for (i in (1:nitems)[items]) { + #*********************************************************** + #** expected item response curves + if ( type=="expected"){ + if (i==1 || !overlay) { + ylim2 <- c(0,max( tamobj$resp[,i], na.rm=TRUE ) ) + if (TT==0) { + graphics::plot(theta, expScore[,i]-0.2, type="l", lwd=3, las=1, ylab="Expected item response", + main=NULL,xlab="", + col="#204776",ylim=ylim2,xlim= ... ) + } else { + if (TT==1) { + lines(theta,expScore[,i]+0.2, type="l", lwd=3, las=1,col="#ca346d") + } else { + lines(theta,expScore[,i], type="l", lwd=2, las=1,col="#999999",lty=2) + } + } + + + } else { + graphics::lines(theta, expScore[,i],type="l", col=i, lwd=3, pch=1) + } + if (observed){ + theta2_i <- theta2 + obScore_i <- obScore[[i]]$x + if (groups_by_item){ + ind_i <- ! is.na(resp[,i]) + resp_i <- resp[ind_i, i, drop=FALSE] + wle_i <- wle[ ind_i ] + res <- plot_tam_grouped_wle( tamobj=tamobj, tammodel=tammodel, + wle=wle_i, ngroups=ngroups, resp=resp_i ) + theta2_i <- res$theta2 + groupnumber_i <- res$groupnumber + aggr <- stats::aggregate(resp_i, list(groupnumber_i), mean, na.rm=TRUE ) + obScore_i <- aggr[,2] + } + #graphics::lines(theta2_i, obScore_i, type="o", lwd=2, pch=1) + } + } + #*********************************************************** + + if ( ndim==1 ){ theta0 <- theta } + + if ( type=="items"){ + rprobs.ii <- rprobs[i,,] + rprobs.ii <- rprobs.ii[ rowMeans( is.na(rprobs.ii) ) < 1, ] + K <- nrow(rprobs.ii) + dat2 <- NULL + #************ + if ( ndim > 1 ){ + B.ii <- B[i,,] + ind.ii <- which( colSums( B.ii ) > 0 )[1] + rprobs0.ii <- rprobs.ii + rprobs0.ii <- stats::aggregate( t(rprobs0.ii), list( theta[,ind.ii] ), mean ) + theta0 <- rprobs0.ii[,1,drop=FALSE] + rprobs.ii <- t( rprobs0.ii[,-1] ) + } + probs_plot[[i]] <- rprobs.ii + #************** + for (kk in 1:K){ + dat2a <- data.frame( "Theta"=theta0[,1], "cat"=kk, "P"=rprobs.ii[kk,] ) + dat2 <- rbind(dat2, dat2a) + } + auto.key <- NULL + simple.key <- paste0("Cat", 1:K - 1) + auto.key <- simple.key + dat2$time <- dat2$cat + dat2$time1 <- paste0("Cat", dat2$time ) + + simple.key <- FALSE + Kpercol <- K + # package graphics + if ( package=="graphics" ){ + kk <- 1 + dfr <- dat2 + dfr1a <- dfr[ dfr$cat==kk, ] + graphics::plot( dfr1a$Theta, dfr1a$P, ylim=c(-.1,1.1), + xlab=expression(theta), + col=pall[kk], type="l", xpd=TRUE,axes=F, ... + ) + axis(1) + axis(2) + grid(nx = NA, + ny = NULL, + lty = 3, col = "lightgray", lwd = 1) + graphics::lines( dfr1a$Theta, dfr1a$P, + col=pall[kk], type="l", ... + ) + for (kk in seq(2,K) ){ + dfr1a <- dfr[ dfr$cat==kk, ] + graphics::lines( dfr1a$Theta, dfr1a$P, col=pall[kk] ) + # graphics::points( dfr1a$Theta, dfr1a$P, pch=kk, col=kk+1 ) + + } + } + + + #*************************************** + + } + #*************** + graphics::par(ask=ask) + } # end item ii + #************************************************* + +} + +png(file = '/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Figures/PNG/icc_residif_2.png',width=800,height=800) +par(xpd=T,mar=c(6.1,5.1,7.6,2.1)) +icc.tam.sup(zaza,type = "expected",export=F,package = "graphics",items = 3,axes=F,TT=0) +icc.tam.sup(zaza2,type = "expected",export=F,package = "graphics",items = 3,axes=F,TT=1) +icc.tam.sup(zaza3,type = "expected",export=F,package = "graphics",items = 3,axes=F,TT=2) +axis(1) +axis(2,las=2) +title(main='Item: \n "Have you been able to enjoy your normal daily activities?" ', font.main=2) +arrows(x0=-2.65,x1=-3,y0=-0.45,length = 0.15,lwd = 2) +arrows(x0=2.65,x1=3,y0=-0.45,length = 0.15,lwd = 2) +arrows(x0=-2,x1=-2,y1=0.375,y0=0.675,length = 0.15,lwd = 2,col="#ca346d") +arrows(x0=0,x1=0,y1=1.9,y0=2.45,length = 0.15,lwd = 2,col="#ca346d") +arrows(x0=2,x1=2,y1=2.875,y0=3.075,length = 0.15,lwd = 2,col="#ca346d") + +text(x=-2,y=0.85,col="#ca346d","-",cex=1.5) +text(x=0,y=2.65,col="#ca346d","-",cex=1.5) +text(x=2,y=3.15,col="#ca346d","-",cex=1.5) +text(x=-2,y=0.0,col="#204776","+",cex=1.5) +text(x=0,y=0.975,col="#204776","+",cex=1.5) +text(x=2,y=2.3,col="#204776","+",cex=1.5) + +arrows(x0=-2,x1=-2,y1=0.305,y0=0.105,length = 0.15,lwd = 2,col="#204776") +arrows(x0=0,x1=0,y1=1.75,y0=1.15,length = 0.15,lwd = 2,col="#204776") +arrows(x0=2,x1=2,y1=2.8,y0=2.425,length = 0.15,lwd = 2,col="#204776") +text(x=-2.55,y=-0.45,"Worse\nmental\nhealth",adj=0) +text(x=2.5,y=-0.45,"Better\nmental\nhealth",adj=1) + +legend(cex=0.8,y=3,x=-3,lty=c(1,1,2,NA,NA),lwd=c(2,2,2,NA,NA),pch = c(NA,NA,NA,NA,NA),col=c('#204776', '#ca346d',"#999999",NA,NA),c("Control group ICC","Treatment group ICC","Overall ICC","Person-item residuals - Control","Person-item residuals - Treatment")) +par(font=5) +legend(cex=0.8,y=3,x=-3,lty=c(NA,NA,NA,NA,NA),lwd=c(NA,NA,NA,NA,NA),pt.cex=c(NA,NA,NA,1.5,1.5),pch = c(NA,NA,NA,174,174),col=c(NA,NA,NA,'#204776', '#ca346d'),c(NA,NA,NA,NA,NA),bty="n") +par(font=1) +title(xlab=expression("Latent variable " * theta * "\n (mental health)")) +par(xpd=F) +dev.off() -bp.dat.power.residif.mask <- as.numeric(res.dat.article.residif.2[res.dat.article.residif.2$N==300 & res.dat.article.residif.2$abs.gamma==0.5 & res.dat.article.residif.2$prop.dif<0.3 & res.dat.article.residif.2$true.beta>0 & res.dat.article.residif.2$true.gamma>0 & res.dat.article.residif.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.residif.2[res.dat.article.residif.2$N==300 & res.dat.article.residif.2$abs.gamma==0.5 & res.dat.article.residif.2$prop.dif<0.3 & res.dat.article.residif.2$true.beta>0 & res.dat.article.residif.2$true.gamma>0 & res.dat.article.residif.2$nb.dif!=0,"theoretical.power"]) -bp.dat.power.residif.magnif <- as.numeric(res.dat.article.residif.2[res.dat.article.residif.2$N==300 & res.dat.article.residif.2$abs.gamma==0.5 & res.dat.article.residif.2$prop.dif<0.3 & res.dat.article.residif.2$true.beta>0 & res.dat.article.residif.2$true.gamma<0 & res.dat.article.residif.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.residif.2[res.dat.article.residif.2$N==300 & res.dat.article.residif.2$abs.gamma==0.5 & res.dat.article.residif.2$prop.dif<0.3 & res.dat.article.residif.2$true.beta>0 & res.dat.article.residif.2$true.gamma<0 & res.dat.article.residif.2$nb.dif!=0,"theoretical.power"]) -bp.dat.power.dif.mask <- as.numeric(res.dat.article.dif.2[res.dat.article.dif.2$N==300 & res.dat.article.dif.2$abs.gamma==0.5 & res.dat.article.dif.2$prop.dif<0.3 & res.dat.article.dif.2$true.beta>0 & res.dat.article.dif.2$true.gamma>0 & res.dat.article.dif.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.dif.2[res.dat.article.dif.2$N==300 & res.dat.article.dif.2$abs.gamma==0.5 & res.dat.article.dif.2$prop.dif<0.3 & res.dat.article.dif.2$true.beta>0 & res.dat.article.dif.2$true.gamma>0 & res.dat.article.dif.2$nb.dif!=0,"theoretical.power"]) -bp.dat.power.dif.magnif <- as.numeric(res.dat.article.dif.2[res.dat.article.dif.2$N==300 & res.dat.article.dif.2$abs.gamma==0.5 & res.dat.article.dif.2$prop.dif<0.3 & res.dat.article.dif.2$true.beta>0 & res.dat.article.dif.2$true.gamma<0 & res.dat.article.dif.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.dif.2[res.dat.article.dif.2$N==300 & res.dat.article.dif.2$abs.gamma==0.5 & res.dat.article.dif.2$prop.dif<0.3 & res.dat.article.dif.2$true.beta>0 & res.dat.article.dif.2$true.gamma<0 & res.dat.article.dif.2$nb.dif!=0,"theoretical.power"]) -bp.dat.power <- data.frame(power=c(bp.dat.power.nodif,bp.dat.power.ignore.mask,bp.dat.power.ignore.magnif,bp.dat.power.rosali.mask,bp.dat.power.rosali.magnif,bp.dat.power.residif.mask,bp.dat.power.residif.magnif, - bp.dat.power.dif.mask,bp.dat.power.dif.magnif), - method=c(rep("NO DIF",12), - rep("MASK1",4),rep("AMPLIFY1",4), - rep("MASK2",4),rep("AMPLIFY2",4), - rep("MASK3",4),rep("AMPLIFY3",4), - rep("MASK4",4),rep("AMPLIFY4",4) )) -bp.dat.power$method <- factor(bp.dat.power$method,levels=c("NO DIF","MASK1","AMPLIFY1","MASK2","AMPLIFY2","MASK3","AMPLIFY3","MASK4","AMPLIFY4")) -boxplot(bp.dat.power$power~bp.dat.power$method,xlab="",pch=3,main="", - ylab="",ylim=c(-1,1),yaxt="n",xaxt="n", - cex.lab=1.45,cex.main=1.5,cex.axis=1.15, - col=c("#e69875","#798A5D","#D4E8B5","#798A5D","#D4E8B5","#798A5D","#D4E8B5","#798A5D","#D4E8B5") - ,border=c("#CD5E35","#697850","#697850","#697850","#697850","#697850","#697850","#697850","#697850"), - width=c(0.8,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4), - at=c(1,2,2.5,3.25,3.75,4.5,5,5.75,6.25)) -axis(2,seq(-1,1,0.25),cex.axis=1.45) -axis(1,c(1,2.25,3.5,4.75,6),labels=c("NO DIF","IGNORE-DIF","ROSALI","RESIDIF","PCM-DIF"),cex.axis=1.45) -abline(h=0,lty=2,col='#595959',lwd=2) +########################## +# DETECTION PERF NODIF +########################## -## DIF 05 50% +#### res.dat ROSALI -bp.dat.power.ignore.mask <- as.numeric(res.dat.article.2[res.dat.article.2$N==300 & res.dat.article.2$abs.gamma==0.5 & res.dat.article.2$prop.dif<0.3 & res.dat.article.2$true.gamma>0 & res.dat.article.2$true.beta>0 & res.dat.article.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.2[res.dat.article.2$N==300 & res.dat.article.2$abs.gamma==0.5 & res.dat.article.2$prop.dif<0.3 & res.dat.article.2$true.gamma>0 & res.dat.article.2$true.beta>0 & res.dat.article.2$nb.dif!=0,"theoretical.power"]) -bp.dat.power.ignore.magnif <- as.numeric(res.dat.article.2[res.dat.article.2$N==300 & res.dat.article.2$abs.gamma==0.5 & res.dat.article.2$prop.dif<0.3 & res.dat.article.2$true.gamma<0 & res.dat.article.2$true.beta>0 & res.dat.article.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.2[res.dat.article.2$N==300 & res.dat.article.2$abs.gamma==0.5 & res.dat.article.2$prop.dif<0.3 & res.dat.article.2$true.gamma<0 & res.dat.article.2$true.beta>0 & res.dat.article.2$nb.dif!=0,"theoretical.power"]) +res.dat.sup.rosali <- res.dat.dif.rosali[,c("N","J","eff.size","nb.dif","dif.size", + "prop.perfect","flexible.detect","moreflexible.detect","dif.detected")] +colnames(res.dat.sup.rosali)[3] <- "true.beta" +colnames(res.dat.sup.rosali)[5] <- "true.gamma" +colnames(res.dat.sup.rosali)[6] <- "perfect" +colnames(res.dat.sup.rosali)[7] <- "flexible" +colnames(res.dat.sup.rosali)[8] <- "mostflexible" +res.dat.sup.rosali[res.dat.sup.rosali$nb.dif==0,"true.gamma"] <- NA +res.dat.sup.rosali[is.na(res.dat.sup.rosali)] <- " " +res.dat.sup.rosali.w <- reshape(res.dat.sup.rosali, + direction = "wide", idvar = c("J","true.beta","nb.dif",'true.gamma'),timevar = "N" ) -bp.dat.power.nodif <- as.numeric(res.dat.article.nodif.2[res.dat.article.nodif.2$true.beta>0,"power"])-as.numeric(res.dat.article.nodif.2[res.dat.article.nodif.2$true.beta>0,"theoretical.power"]) +res.dat.sup.rosali.nodif <- res.dat.sup.rosali.w[res.dat.sup.rosali.w$nb.dif==0,] +res.dat.sup.rosali.dif <- res.dat.sup.rosali.w[res.dat.sup.rosali.w$nb.dif!=0,] -bp.dat.power.rosali.mask <- as.numeric(res.dat.article.rosali.2[res.dat.article.rosali.2$N==300 & res.dat.article.rosali.2$abs.gamma==0.5 & res.dat.article.rosali.2$prop.dif<0.3 & res.dat.article.rosali.2$true.beta>0 & res.dat.article.rosali.2$true.gamma>0 & res.dat.article.rosali.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.rosali.2[res.dat.article.rosali.2$N==300 & res.dat.article.rosali.2$abs.gamma==0.5 & res.dat.article.rosali.2$prop.dif<0.3 & res.dat.article.rosali.2$true.beta>0 & res.dat.article.rosali.2$true.gamma>0 & res.dat.article.rosali.2$nb.dif!=0,"theoretical.power"]) -bp.dat.power.rosali.magnif <- as.numeric(res.dat.article.rosali.2[res.dat.article.rosali.2$N==300 & res.dat.article.rosali.2$abs.gamma==0.5 & res.dat.article.rosali.2$prop.dif<0.3 & res.dat.article.rosali.2$true.beta>0 & res.dat.article.rosali.2$true.gamma<0 & res.dat.article.rosali.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.rosali.2[res.dat.article.rosali.2$N==300 & res.dat.article.rosali.2$abs.gamma==0.5 & res.dat.article.rosali.2$prop.dif<0.3 & res.dat.article.rosali.2$true.beta>0 & res.dat.article.rosali.2$true.gamma<0 & res.dat.article.rosali.2$nb.dif!=0,"theoretical.power"]) +#### res.dat RESIDIF -bp.dat.power.residif.mask <- as.numeric(res.dat.article.residif.2[res.dat.article.residif.2$N==300 & res.dat.article.residif.2$abs.gamma==0.5 & res.dat.article.residif.2$prop.dif<0.3 & res.dat.article.residif.2$true.beta>0 & res.dat.article.residif.2$true.gamma>0 & res.dat.article.residif.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.residif.2[res.dat.article.residif.2$N==300 & res.dat.article.residif.2$abs.gamma==0.5 & res.dat.article.residif.2$prop.dif<0.3 & res.dat.article.residif.2$true.beta>0 & res.dat.article.residif.2$true.gamma>0 & res.dat.article.residif.2$nb.dif!=0,"theoretical.power"]) -bp.dat.power.residif.magnif <- as.numeric(res.dat.article.residif.2[res.dat.article.residif.2$N==300 & res.dat.article.residif.2$abs.gamma==0.5 & res.dat.article.residif.2$prop.dif<0.3 & res.dat.article.residif.2$true.beta>0 & res.dat.article.residif.2$true.gamma<0 & res.dat.article.residif.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.residif.2[res.dat.article.residif.2$N==300 & res.dat.article.residif.2$abs.gamma==0.5 & res.dat.article.residif.2$prop.dif<0.3 & res.dat.article.residif.2$true.beta>0 & res.dat.article.residif.2$true.gamma<0 & res.dat.article.residif.2$nb.dif!=0,"theoretical.power"]) +res.dat.sup.residif <- res.dat.dif.resali[,c("N","J","eff.size","nb.dif","dif.size", + "prop.perfect","flexible.detect","moreflexible.detect","dif.detected")] +colnames(res.dat.sup.residif)[3] <- "true.beta" +colnames(res.dat.sup.residif)[5] <- "true.gamma" +colnames(res.dat.sup.residif)[6] <- "perfect" +colnames(res.dat.sup.residif)[7] <- "flexible" +colnames(res.dat.sup.residif)[8] <- "mostflexible" +res.dat.sup.residif[res.dat.sup.residif$nb.dif==0,"true.gamma"] <- NA +res.dat.sup.residif[is.na(res.dat.sup.residif)] <- " " +res.dat.sup.residif.w <- reshape(res.dat.sup.residif, + direction = "wide", idvar = c("J","true.beta","nb.dif",'true.gamma'),timevar = "N" ) -bp.dat.power.dif.mask <- as.numeric(res.dat.article.dif.2[res.dat.article.dif.2$N==300 & res.dat.article.dif.2$abs.gamma==0.5 & res.dat.article.dif.2$prop.dif<0.3 & res.dat.article.dif.2$true.beta>0 & res.dat.article.dif.2$true.gamma>0 & res.dat.article.dif.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.dif.2[res.dat.article.dif.2$N==300 & res.dat.article.dif.2$abs.gamma==0.5 & res.dat.article.dif.2$prop.dif<0.3 & res.dat.article.dif.2$true.beta>0 & res.dat.article.dif.2$true.gamma>0 & res.dat.article.dif.2$nb.dif!=0,"theoretical.power"]) -bp.dat.power.dif.magnif <- as.numeric(res.dat.article.dif.2[res.dat.article.dif.2$N==300 & res.dat.article.dif.2$abs.gamma==0.5 & res.dat.article.dif.2$prop.dif<0.3 & res.dat.article.dif.2$true.beta>0 & res.dat.article.dif.2$true.gamma<0 & res.dat.article.dif.2$nb.dif!=0,"power"])-as.numeric(res.dat.article.dif.2[res.dat.article.dif.2$N==300 & res.dat.article.dif.2$abs.gamma==0.5 & res.dat.article.dif.2$prop.dif<0.3 & res.dat.article.dif.2$true.beta>0 & res.dat.article.dif.2$true.gamma<0 & res.dat.article.dif.2$nb.dif!=0,"theoretical.power"]) +res.dat.sup.residif.nodif <- res.dat.sup.residif.w[res.dat.sup.residif.w$nb.dif==0,] +res.dat.sup.residif.dif <- res.dat.sup.residif.w[res.dat.sup.residif.w$nb.dif!=0,] -bp.dat.power <- data.frame(power=c(bp.dat.power.nodif,bp.dat.power.ignore.mask,bp.dat.power.ignore.magnif,bp.dat.power.rosali.mask,bp.dat.power.rosali.magnif,bp.dat.power.residif.mask,bp.dat.power.residif.magnif, - bp.dat.power.dif.mask,bp.dat.power.dif.magnif), - method=c(rep("NO DIF",12), - rep("MASK1",4),rep("AMPLIFY1",4), - rep("MASK2",4),rep("AMPLIFY2",4), - rep("MASK3",4),rep("AMPLIFY3",4), - rep("MASK4",4),rep("AMPLIFY4",4) )) -bp.dat.power$method <- factor(bp.dat.power$method,levels=c("NO DIF","MASK1","AMPLIFY1","MASK2","AMPLIFY2","MASK3","AMPLIFY3","MASK4","AMPLIFY4")) +#### Sortie table nodif +write.csv(cbind(res.dat.sup.rosali.nodif[,c(1:3,8,12,16)],res.dat.sup.residif.nodif[,c(8,12,16)]),"/home/corentin/Documents/These/Valorisation/Articles/Simulations 1/Tables/nodif_detect.csv") -boxplot(bp.dat.power$power~bp.dat.power$method,xlab="",pch=3,main="", - ylab="",ylim=c(-1,1),yaxt="n",xaxt="n", - cex.lab=1.45,cex.main=1.5,cex.axis=1.15, - col=c("#e69875","#798A5D","#D4E8B5","#798A5D","#D4E8B5","#798A5D","#D4E8B5","#798A5D","#D4E8B5") - ,border=c("#CD5E35","#697850","#697850","#697850","#697850","#697850","#697850","#697850","#697850"), - width=c(0.8,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4), - at=c(1,2,2.5,3.25,3.75,4.5,5,5.75,6.25)) -axis(2,seq(-1,1,0.25),cex.axis=1.45) -axis(1,c(1,2.25,3.5,4.75,6),labels=c("NO DIF","IGNORE-DIF","ROSALI","RESIDIF","PCM-DIF"),cex.axis=1.45) -abline(h=0,lty=2,col='#595959',lwd=2) -par(mfrow=c(1,1)) +########################## +# DETECTION PERF DIF +########################## \ No newline at end of file diff --git a/RProject/Scripts/Analysis/functions/resali.R b/RProject/Scripts/Analysis/functions/resali.R index 7fefac8..02ad6f1 100644 --- a/RProject/Scripts/Analysis/functions/resali.R +++ b/RProject/Scripts/Analysis/functions/resali.R @@ -30,6 +30,7 @@ resali <- function(df=NULL,items=NULL,group=NULL,verbose=T) { pval[i] <- res.anova[[i]][1,"Pr(>F)"] fval[i] <- res.anova[[i]][1,'F value'] } + print(res.anova) if (verbose) { cat('DONE\n') cat('-----------------------------------------------------------\n') diff --git a/RProject/Scripts/Analysis/resali_analysis.R b/RProject/Scripts/Analysis/resali_analysis.R index 5e32ac3..b4f774f 100644 --- a/RProject/Scripts/Analysis/resali_analysis.R +++ b/RProject/Scripts/Analysis/resali_analysis.R @@ -4,7 +4,7 @@ #----------------------------------------------------------------------------# ############################################################################## -source(paste0(getwd(),"/functions/resali.R")) +source("/home/corentin/Documents/These/Recherche/Simulations/RProject/Scripts/Analysis/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))) diff --git a/RProject/Scripts/Analysis/residif_analysis.R b/RProject/Scripts/Analysis/residif_analysis.R new file mode 100644 index 0000000..71d9833 --- /dev/null +++ b/RProject/Scripts/Analysis/residif_analysis.R @@ -0,0 +1,279 @@ +############################################################################## +#----------------------------------------------------------------------------# +################################### RESIDIF ################################## +#----------------------------------------------------------------------------# +############################################################################## + +generate_residif <- 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 <- residif(df=dat[dat$replication==1,],items = paste0("item",1:7),grp="TT",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 <- residif(df=dat[dat$replication==k,],items = paste0("item",1:7),grp="TT",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 <- residif(df=dat[dat$replication==1,],items = paste0("item",1:4),grp="TT",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 <- residif(df=dat[dat$replication==k,],items = paste0("item",1:4),grp="TT",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(c(2,4),function(x) paste0(x,c('A','B','C'))),sapply(c(6,8),function(x) paste0(x,c('A','B','C','D','E')))) + +results2 <- c(sapply(seq(10,20,2),function(x) paste0(x,c('A','B','C','D','E')))) + +results <- c(sapply(c(50,100,300),function(x) paste0(results,'_',x))) + +results2 <- c(sapply(c(50,100,300),function(x) paste0(results2,'_',x))) + +results <- sort(results) + +results2 <- sort(results2) + +results <- c(results,results2) + +for (r in c(results[73:138])) { + cat(paste0(r,"\n")) + cat(paste0("-------------------------------------------","\n")) + write.csv(generate_residif(r,"TT"),paste0("/home/corentin/Documents/These/Recherche/residif/detection/",r,".csv")) + cat(paste0("-------------------------------------------","\n")) +} + +############################################################################## +#----------------------------------------------------------------------------# +################################### NEWDATA ################################## +#----------------------------------------------------------------------------# +############################################################################## + +## Liste des scenarios + +results <- c(sapply(c(2,4),function(x) paste0(x,c('A','B','C'))),sapply(c(6,8),function(x) paste0(x,c('A','B','C','D','E')))) + +results2 <- c(sapply(seq(10,20,2),function(x) paste0(x,c('A','B','C','D','E')))) + +results <- c(sapply(c(50,100,300),function(x) paste0(results,'_',x))) + +results2 <- c(sapply(c(50,100,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/residif/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/residif/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/residif/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/residif/detection_data/",r,".csv")) +} + diff --git a/RProject/Scripts/Tools/generate_item_difficulties.R b/RProject/Scripts/Tools/generate_item_difficulties.R index ca53b8c..5303cab 100644 --- a/RProject/Scripts/Tools/generate_item_difficulties.R +++ b/RProject/Scripts/Tools/generate_item_difficulties.R @@ -15,7 +15,7 @@ generate_diff_irt <- function(J=7,M=4) { } difficulties = matrix(c(0), J,M-1) rownames(difficulties)=paste("item",1:J) - colnames(difficulties)=paste("Moda", 1:(M-1)) + colnames(difficulties)=paste("delta", 1:(M-1)) for (j in 1:J){ difficulties[j,1] = qnorm(j/(J+1)) }