diff --git a/RProject/Scripts/aggregation.R b/RProject/Scripts/aggregation.R index 41b1c20..ffa4915 100644 --- a/RProject/Scripts/aggregation.R +++ b/RProject/Scripts/aggregation.R @@ -629,24 +629,39 @@ results <- c(results,results2) compile_simulation2_resali <- function(scenario) { name <- as.numeric(gsub("[^0-9.-]", "", substr(scenario,start=0,stop=2))) if (substr(scenario,start=nchar(scenario)-1,stop=nchar(scenario))=="50" & name>0) { - s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/resali-DIF/N50/',scenario,'_original.xls')) + s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/RESALI/Results/N50/',scenario,'_original.xls')) } if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="100" & name>0) { - s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/resali-DIF/N100/',scenario,'_original.xls')) + s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/RESALI/Results/N100/',scenario,'_original.xls')) } if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="200" & name>0) { - s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/resali-DIF/N200/',scenario,'_original.xls')) + s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/RESALI/Results/N200/',scenario,'_original.xls')) } if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="300" & name>0) { - s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/resali-DIF/N300/',scenario,'_original.xls')) + s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/RESALI/Results/N300/',scenario,'_original.xls')) } J <- max(which(sapply(1:7,function(x) paste0('item',x) %in% colnames(s) | paste0('item',x,'_1') %in% colnames(s)))) M <- 1+sum(sapply(1:3,function(x) paste0('item1_',x) %in% colnames(s) )) if (M==1) {M <- 2} - nb.dif <- max(which(sapply(1:3,function(x) paste0('dif',x) %in% colnames(s) | paste0('dif',x,'_1') %in% colnames(s)))) + nb.dif.true <- ifelse(name<=4,0,ifelse(name<=8,1,ifelse(name<=16,2,3))) + if (name %in% c(3,4,13:20)) { + m.nb.dif.detect <- mean(ifelse(is.na(s$dif_1_1),0, + ifelse(is.na(s$dif_2_1),1, + ifelse(is.na(s$dif_3_1),2, + ifelse(is.na(s$dif_4_1),3, + ifelse(is.na(s$dif_5_1),4, + ifelse(is.na(s$dif_6_1),5, + ifelse(is.na(s$dif_7_1),6,7)))))))) + } + if (!(name %in% c(3,4,13:20))) { + m.nb.dif.detect <- mean(ifelse(is.na(s$dif_1_1),0, + ifelse(is.na(s$dif_2_1),1, + ifelse(is.na(s$dif_3_1),2, + ifelse(is.na(s$dif_4_1),3,4))))) + } if (J==4) { if (M==2) { - a <- data.frame(m.item1=mean(s$item1),m.item2=mean(s$item2),m.item3=mean(s$item3),m.item4=mean(s$item4)) + a <- data.frame(m.item1=mean(s$item1_1),m.item2=mean(s$item2_1),m.item3=mean(s$item3_1),m.item4=mean(s$item4_1)) } else { a <- data.frame(m.item1_1=mean(s$item1_1),m.item1_2=mean(s$item1_2),m.item1_3=mean(s$item1_3), m.item2_1=mean(s$item2_1),m.item2_2=mean(s$item2_2),m.item2_3=mean(s$item2_3), @@ -656,8 +671,8 @@ compile_simulation2_resali <- function(scenario) { } } else { if (M==2) { - a <- data.frame(m.item1=mean(s$item1),m.item2=mean(s$item2),m.item3=mean(s$item3),m.item4=mean(s$item4), - m.item5=mean(s$item5),m.item6=mean(s$item6),m.item7=mean(s$item7)) + a <- data.frame(m.item1=mean(s$item1_1),m.item2=mean(s$item2_1),m.item3=mean(s$item3_1),m.item4=mean(s$item4_1), + m.item5=mean(s$item5_1),m.item6=mean(s$item6_1),m.item7=mean(s$item7_1)) } else { a <- data.frame(m.item1_1=mean(s$item1_1),m.item1_2=mean(s$item1_2),m.item1_3=mean(s$item1_3), m.item2_1=mean(s$item2_1),m.item2_2=mean(s$item2_2),m.item2_3=mean(s$item2_3), @@ -679,12 +694,106 @@ compile_simulation2_resali <- function(scenario) { J=J, M=M, eff.size=eff.size, - nb.dif=nb.dif, + nb.dif=nb.dif.true, + m.nb.dif.detect=m.nb.dif.detect, dif.size=dif.size ) true.value.in.ci <- eff.size <= s$beta+1.96*s$se_beta & eff.size >= s$beta-1.96*s$se_beta beta.same.sign.truebeta.p <- ifelse(rep(eff.size,nrow(s))==0,NA,(rep(eff.size,nrow(s))/s$beta)>0) num.reject <- which((s$beta-1.96*s$se_beta)>0 | (s$beta+1.96*s$se_beta)<0) + dif.d <- mean(sapply(1:1000,function(x) any(!is.na(s[x,paste0("dif_",1:unique(b$J),"_1")])))) + if (nb.dif.true==0) { + prop.perfect <- NA + flexible.detect <- NA + moreflexible.detect <- NA + } + if (nb.dif.true==1 & unique(b$J)==4 & unique(b$M)==4) { + perfect.detection <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:4)]))==1,s[x,"dif_detect_unif_1"]==1 & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:4),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1")]) + ,0) ) + prop.perfect <- mean(perfect.detection) + flexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:4)]))==1,s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:4),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1")]) + ,0) ) + flexible.detect <- mean(flexible.detect) + moreflexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:4)]))!=0,s[x,"real_dif_1"]%in%c(s[x,paste0("dif_detect_",1:4)]),0) ) + moreflexible.detect <- mean(moreflexible.detect) + } + if (nb.dif.true==2 & unique(b$J)==4 & unique(b$M)==4) { + perfect.detection <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:4)]))==2,s[x,"dif_detect_unif_1"]==1 & s[x,"dif_detect_unif_2"]==1 & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:4),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1","real_dif_2")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:4),function(y) !is.na(s[x,y])))[2])]%in%c(s[x,c("real_dif_1","real_dif_2")]) + ,0) ) + prop.perfect <- mean(perfect.detection) + flexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:4)]))==2,s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:4),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1","real_dif_2")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:4),function(y) !is.na(s[x,y])))[2])]%in%c(s[x,c("real_dif_1","real_dif_2")]) + ,0) ) + flexible.detect <- mean(flexible.detect) + moreflexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:4)]))!=0,s[x,"real_dif_1"]%in%c(s[x,paste0("dif_detect_",1:4)]) & + s[x,"real_dif_2"]%in%c(s[x,paste0("dif_detect_",1:4)]),0) ) + moreflexible.detect <- mean(moreflexible.detect) + } + if (nb.dif.true==2 & unique(b$J)==7 & unique(b$M)==4) { + perfect.detection <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:7)]))==2,s[x,"dif_detect_unif_1"]==1 & s[x,"dif_detect_unif_2"]==1 & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1","real_dif_2")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[2])]%in%c(s[x,c("real_dif_1","real_dif_2")]) + ,0) ) + prop.perfect <- mean(perfect.detection) + flexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:7)]))==2,s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1","real_dif_2")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[2])]%in%c(s[x,c("real_dif_1","real_dif_2")]) + ,0) ) + flexible.detect <- mean(flexible.detect) + moreflexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:7)]))!=0,s[x,"real_dif_1"]%in%c(s[x,paste0("dif_detect_",1:7)]) & + s[x,"real_dif_2"]%in%c(s[x,paste0("dif_detect_",1:7)]),0) ) + moreflexible.detect <- mean(moreflexible.detect) + } + if (nb.dif.true==3 & unique(b$J)==7 & unique(b$M)==4) { + perfect.detection <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:7)]))==3,s[x,"dif_detect_unif_1"]==1 & s[x,"dif_detect_unif_2"]==1 & s[x,"dif_detect_unif_3"]==1 & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1","real_dif_2","real_dif_3")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[2])]%in%c(s[x,c("real_dif_1","real_dif_2","real_dif_3")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[3])]%in%c(s[x,c("real_dif_1","real_dif_2","real_dif_3")]) + ,0) ) + prop.perfect <- mean(perfect.detection) + flexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:7)]))==3,s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1","real_dif_2","real_dif_3")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[2])]%in%c(s[x,c("real_dif_1","real_dif_2","real_dif_3")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[3])]%in%c(s[x,c("real_dif_1","real_dif_2","real_dif_3")]) + ,0) ) + flexible.detect <- mean(flexible.detect) + moreflexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:7)]))!=0,s[x,"real_dif_1"]%in%c(s[x,paste0("dif_detect_",1:7)]) & + s[x,"real_dif_2"]%in%c(s[x,paste0("dif_detect_",1:7)]) & s[x,"real_dif_3"]%in%c(s[x,paste0("dif_detect_",1:7)]),0) ) + moreflexible.detect <- mean(moreflexible.detect) + } + if (nb.dif.true==1 & unique(b$J)==4 & unique(b$M)==2) { + perfect.detection <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:4)]))==1, s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:4),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1")]) + ,0) ) + prop.perfect <- mean(perfect.detection) + flexible.detect <- prop.perfect + moreflexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:4)]))!=0,s[x,"real_dif_1"]%in%c(s[x,paste0("dif_detect_",1:4)]),0) ) + moreflexible.detect <- mean(moreflexible.detect) + percent.detect <- sapply(1:1000,function(x) sum(s[x,paste0("dif_detect_",1:4)]%in%c(s[x,c("real_dif_1")]))/1) + percent.detect <- mean(percent.detect) + } + if (nb.dif.true==2 & unique(b$J)==4 & unique(b$M)==2) { + perfect.detection <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:4)]))==2,s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:4),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1","real_dif_2")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:4),function(y) !is.na(s[x,y])))[2])]%in%c(s[x,c("real_dif_1","real_dif_2")]) + ,0) ) + prop.perfect <- mean(perfect.detection) + flexible.detect <- prop.perfect + moreflexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:4)]))!=0,s[x,"real_dif_1"]%in%c(s[x,paste0("dif_detect_",1:4)]) & + s[x,"real_dif_2"]%in%c(s[x,paste0("dif_detect_",1:4)]),0) ) + moreflexible.detect <- mean(moreflexible.detect) + percent.detect <- sapply(1:1000,function(x) sum(s[x,paste0("dif_detect_",1:4)]%in%c(s[x,c("real_dif_1","real_dif_2")]))/2) + percent.detect <- mean(percent.detect) + + } + if (nb.dif.true==2 & unique(b$J)==7 & unique(b$M)==2) { + perfect.detection <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:7)]))==2,s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1","real_dif_2")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[2])]%in%c(s[x,c("real_dif_1","real_dif_2")]) + ,0) ) + prop.perfect <- mean(perfect.detection) + flexible.detect <- prop.perfect + moreflexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:7)]))!=0,s[x,"real_dif_1"]%in%c(s[x,paste0("dif_detect_",1:7)]) & + s[x,"real_dif_2"]%in%c(s[x,paste0("dif_detect_",1:7)]),0) ) + moreflexible.detect <- mean(moreflexible.detect) + percent.detect <- sapply(1:1000,function(x) sum(s[x,paste0("dif_detect_",1:7)]%in%c(s[x,c("real_dif_1","real_dif_2")]))/2) + percent.detect <- mean(percent.detect) + } + if (nb.dif.true==3 & unique(b$J)==7 & unique(b$M)==2) { + perfect.detection <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:7)]))==3,s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1","real_dif_2","real_dif_3")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[2])]%in%c(s[x,c("real_dif_1","real_dif_2","real_dif_3")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[3])]%in%c(s[x,c("real_dif_1","real_dif_2","real_dif_3")]) + ,0) ) + prop.perfect <- mean(perfect.detection) + flexible.detect <- prop.perfect + moreflexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:7)]))!=0,s[x,"real_dif_1"]%in%c(s[x,paste0("dif_detect_",1:7)]) & + s[x,"real_dif_2"]%in%c(s[x,paste0("dif_detect_",1:7)]) & s[x,"real_dif_3"]%in%c(s[x,paste0("dif_detect_",1:7)]),0) ) + moreflexible.detect <- mean(moreflexible.detect) + percent.detect <- sapply(1:1000,function(x) sum(s[x,paste0("dif_detect_",1:7)]%in%c(s[x,c("real_dif_1","real_dif_2","real_dif_3")]))/3) + percent.detect <- mean(percent.detect) + } z <- data.frame(m.beta=mean(s$beta), se.empirical.beta=sd(s$beta), se.analytical.beta=mean(s$se_beta), @@ -693,7 +802,12 @@ compile_simulation2_resali <- function(scenario) { true.value.in.ci.p=mean(true.value.in.ci), h0.rejected.p=mean( (s$beta-1.96*s$se_beta)>0 | (s$beta+1.96*s$se_beta)<0 ), beta.same.sign.truebeta.p=mean(beta.same.sign.truebeta.p), - beta.same.sign.truebeta.signif.p=mean(beta.same.sign.truebeta.p[num.reject]) + beta.same.sign.truebeta.signif.p=mean(beta.same.sign.truebeta.p[num.reject]), + dif.detected=dif.d, + prop.perfect=prop.perfect, + flexible.detect=flexible.detect, + moreflexible.detect=moreflexible.detect, + percent.detect=ifelse(name%%2==0,NA,percent.detect) ) d <- cbind(b,a,z) d$prop. diff --git a/RProject/Scripts/functions/resali.R b/RProject/Scripts/functions/resali.R index 36ba0f3..7fefac8 100644 --- a/RProject/Scripts/functions/resali.R +++ b/RProject/Scripts/functions/resali.R @@ -48,8 +48,8 @@ resali <- function(df=NULL,items=NULL,group=NULL,verbose=T) { res.uni <- res.anova[[which.max(fval)]][3,"Pr(>F)"]>0.05 res.uniform <- c(res.uniform,res.uni) items_n <- c(items_n[items_n!=paste0('item',res.item)],paste0("item",res.item,c("noTT","TT"))) - dat[,paste0("item",res.item,'TT')] <- dat[dat$TT==1,paste0('item',res.item)] - dat[,paste0("item",res.item,'noTT')] <- dat[dat$TT==0,paste0('item',res.item)] + dat[dat$TT==1,paste0("item",res.item,'TT')] <- dat[dat$TT==1,paste0('item',res.item)] + dat[dat$TT==0,paste0("item",res.item,'noTT')] <- dat[dat$TT==0,paste0('item',res.item)] resp <- dat[,items_n] grp <- dat[,group] pcm_while <- TAM::tam.mml(resp=resp,Y=grp,irtmodel = "PCM",est.variance = T,verbose=F) @@ -118,8 +118,8 @@ resali <- function(df=NULL,items=NULL,group=NULL,verbose=T) { res.item <- gsub("[a-z]", "",colnames(resp)[which.max(fval)]) res.items <- c(res.items,res.item) items_n <- c(items_n[items_n!=paste0('item',res.item)],paste0("item",res.item,c("noTT","TT"))) - dat[,paste0("item",res.item,'TT')] <- dat[dat$TT==1,paste0('item',res.item)] - dat[,paste0("item",res.item,'noTT')] <- dat[dat$TT==0,paste0('item',res.item)] + dat[dat$TT==1,paste0("item",res.item,'TT')] <- dat[dat$TT==1,paste0('item',res.item)] + dat[dat$TT==0,paste0("item",res.item,'noTT')] <- dat[dat$TT==0,paste0('item',res.item)] resp <- dat[,items_n] grp <- dat[,group] pcm_while <- TAM::tam.mml(resp=resp,Y=grp,irtmodel = "PCM",est.variance = T,verbose=F)