From 49d465573f0070723598fb6ae9bda8193c389c7f Mon Sep 17 00:00:00 2001 From: corentinchoisy Date: Tue, 7 May 2024 10:15:41 +0200 Subject: [PATCH] Experimented some MDC methods --- .gitignore | 1 + RProject/Scripts/aggregation.R | 88 ++++++++++++++++++++++++++------- RProject/Scripts/mdc.R | 89 ++++++++++++++++++++++++++++++++++ 3 files changed, 160 insertions(+), 18 deletions(-) create mode 100644 RProject/Scripts/mdc.R diff --git a/.gitignore b/.gitignore index a52f3bc..10c3631 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,4 @@ .Rproj.user *.RData *.Rhistory +*.log diff --git a/RProject/Scripts/aggregation.R b/RProject/Scripts/aggregation.R index ffa4915..fcd9141 100644 --- a/RProject/Scripts/aggregation.R +++ b/RProject/Scripts/aggregation.R @@ -9,6 +9,7 @@ library(doMC) library(parallel) library(pbmcapply) library(funprog) +library(plyr) library(dplyr) library(readxl) @@ -480,6 +481,7 @@ compile_simulation2_rosali <- function(scenario) { prop.perfect <- NA flexible.detect <- NA moreflexible.detect <- NA + any.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")]) @@ -490,6 +492,7 @@ compile_simulation2_rosali <- function(scenario) { 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) + any.detect <- 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")]) @@ -501,6 +504,9 @@ compile_simulation2_rosali <- function(scenario) { 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) + any.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) ) + any.detect <- mean(any.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")]) @@ -512,6 +518,9 @@ compile_simulation2_rosali <- function(scenario) { 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) + any.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) ) + any.detect <- mean(any.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")]) @@ -523,6 +532,9 @@ compile_simulation2_rosali <- function(scenario) { 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) + any.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) ) + any.detect <- mean(any.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")]) @@ -531,8 +543,7 @@ compile_simulation2_rosali <- function(scenario) { 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) + any.detect <- moreflexible.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")]) @@ -542,8 +553,9 @@ compile_simulation2_rosali <- function(scenario) { 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) + any.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) ) + any.detect <- mean(any.detect) } if (nb.dif.true==2 & unique(b$J)==7 & unique(b$M)==2) { @@ -554,8 +566,9 @@ compile_simulation2_rosali <- function(scenario) { 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) + any.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) ) + any.detect <- mean(any.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")]) @@ -565,8 +578,9 @@ compile_simulation2_rosali <- function(scenario) { 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) + any.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) ) + any.detect <- mean(any.detect) } z <- data.frame(m.beta=mean(s$beta), se.empirical.beta=sd(s$beta), @@ -581,7 +595,7 @@ compile_simulation2_rosali <- function(scenario) { prop.perfect=prop.perfect, flexible.detect=flexible.detect, moreflexible.detect=moreflexible.detect, - percent.detect=ifelse(name%%2==0,NA,percent.detect) + any.detect=any.detect ) d <- cbind(b,a,z) d$prop. @@ -706,6 +720,7 @@ compile_simulation2_resali <- function(scenario) { prop.perfect <- NA flexible.detect <- NA moreflexible.detect <- NA + any.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")]) @@ -716,6 +731,7 @@ compile_simulation2_resali <- function(scenario) { 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) + any.detect <- 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")]) @@ -727,6 +743,9 @@ compile_simulation2_resali <- function(scenario) { 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) + any.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) ) + any.detect <- mean(any.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")]) @@ -738,6 +757,9 @@ compile_simulation2_resali <- function(scenario) { 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) + any.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) ) + any.detect <- mean(any.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")]) @@ -749,6 +771,9 @@ compile_simulation2_resali <- function(scenario) { 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) + any.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) ) + any.detect <- mean(any.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")]) @@ -757,8 +782,7 @@ compile_simulation2_resali <- function(scenario) { 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) + any.detect <- moreflexible.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")]) @@ -768,8 +792,9 @@ compile_simulation2_resali <- function(scenario) { 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) + any.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) ) + any.detect <- mean(any.detect) } if (nb.dif.true==2 & unique(b$J)==7 & unique(b$M)==2) { @@ -780,8 +805,9 @@ compile_simulation2_resali <- function(scenario) { 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) + any.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) ) + any.detect <- mean(any.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")]) @@ -791,8 +817,9 @@ compile_simulation2_resali <- function(scenario) { 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) + any.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) ) + any.detect <- mean(any.detect) } z <- data.frame(m.beta=mean(s$beta), se.empirical.beta=sd(s$beta), @@ -807,7 +834,7 @@ compile_simulation2_resali <- function(scenario) { prop.perfect=prop.perfect, flexible.detect=flexible.detect, moreflexible.detect=moreflexible.detect, - percent.detect=ifelse(name%%2==0,NA,percent.detect) + any.detect=any.detect ) d <- cbind(b,a,z) d$prop. @@ -827,6 +854,31 @@ for (x in results[seq(2,length(results))]) { res.dat.dif.resali$bias <- res.dat.dif.resali$eff.size-res.dat.dif.resali$m.beta +############################################################################## +#----------------------------------------------------------------------------# +######################### AGGREGATION OF ALL METHODS ######################### +#----------------------------------------------------------------------------# +############################################################################## + +# Items dichotomiques + +res.dat$method <- "NONE" +res.dat.dif$method <- "PERFECT" +res.dat.dif.rosali$method <- "ROSALI" +res.dat.dif.resali$method <- "RESIDUS" + +res.dat.dicho <- res.dat[res.dat$M==2,] +res.dat.dicho <- rbind(res.dat.dicho,res.dat.dif[res.dat.dif$M==2,]) +res.dat.dicho <- rbind.fill(res.dat.dicho,res.dat.dif.rosali[res.dat.dif.rosali$M==2,]) +res.dat.dicho <- rbind(res.dat.dicho,res.dat.dif.rosali[res.dat.dif.resali$M==2,]) + +# Items polytomiques + +res.dat.poly <- res.dat[res.dat$M==4,] +res.dat.poly <- rbind(res.dat.poly,res.dat.dif[res.dat.dif$M==4,]) +res.dat.poly <- rbind.fill(res.dat.poly,res.dat.dif.rosali[res.dat.dif.rosali$M==4,]) +res.dat.poly <- rbind(res.dat.poly,res.dat.dif.rosali[res.dat.dif.resali$M==4,]) + ############################################################################## diff --git a/RProject/Scripts/mdc.R b/RProject/Scripts/mdc.R new file mode 100644 index 0000000..4d8e0e5 --- /dev/null +++ b/RProject/Scripts/mdc.R @@ -0,0 +1,89 @@ +############################################################################## +#----------------------------------------------------------------------------# +############################# DATA TRANSFORMATION ############################ +#----------------------------------------------------------------------------# +############################################################################## + +# Import ROSALI and RESALI + +ros_mdc <- read_excel("/home/corentin/Documents/These/Recherche/Simulations/Analysis/ROSALI-DIF/N300/6A_300_original.xls") +res_mdc <- read_excel("/home/corentin/Documents/These/Recherche/Simulations/Analysis/RESALI/Results/N300/6A_300_original.xls") + + +# Perform MH + +library(difR) + +dat_mh <- read.csv('/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N300/scenario_6A_300.csv')[,c("item1","item2","item3","item4",'replication',"TT")] + +det_mh <- c() +for (k in 1:1000) { + if (k%%1000==0) { + cat(paste0(k,'/1000\n')) + } + dat_mh_temp <- dat_mh[dat_mh$replication==k,c("item1",'item2',"item3","item4",'TT')] + aa <- difMH(Data=dat_mh_temp,group = "TT",focal.name = 0,exact=F) + det_mh <- c(det_mh,1:4 %in% aa$DIFitems) +} + +# Create 1 line per item per replication in df +library(tidyr) + +da <- as.data.frame(sapply(1:4, function(k) sapply(1:1000,function(x) k%in%ros_mdc[x,paste0("dif_detect_",1:4)]))) +db <- as.data.frame(sapply(1:4, function(k) sapply(1:1000,function(x) k%in%res_mdc[x,paste0("dif_detect_",1:4)]))) +dc <- as.data.frame(sapply(1:4, function(k) sapply(1:1000,function(x) k%in%res_mdc[x,paste0("real_dif_",1)]))) + +data_mdca <- data.frame(rosali=da) +data_mdca <- pivot_longer(data_mdca,cols=1:4) +data_mdcb <- data.frame(resali=db) +data_mdcb <- pivot_longer(data_mdcb,cols=1:4) +data_mdcc <- data.frame(real=dc) +data_mdcc <- pivot_longer(data_mdcc,cols=1:4) + +data_mdc <- cbind(data_mdca,data_mdcb,data_mdcc)[,c(2,4,6)] +colnames(data_mdc) <- c("rosali","resali","real") + +make_repl <- function(kk) { + b <- c() + for (k in kk) { + a <- rep(k,4) + b <- c(b,a) + } + return(b) +} + +data_mdc$mh <- det_mh + +data_mdc$replication <- make_repl(1:1000) + + +############################################################################## +#----------------------------------------------------------------------------# +########################### FIT DIF DETECTION MODEL ########################## +#----------------------------------------------------------------------------# +############################################################################## + +# Fit TAN model + + + +# Fit logistic model, stratified on replication + +mod_glm <- glm(formula = real~rosali+resali,data = data_mdc[1:2000,],family = binomial()) +data_valid <- data_mdc[2000:4000,] +data_valid$predict <- predict(mod_glm,newdata = data_valid) +roc_c <- pROC::roc(response=data_valid$real,predictor=data_valid$predict) + +data_mdc$logit_pred <- predict(mod_glm,newdata = data_mdc)>=-0.6275167 + +perf_moreflex <- c() +for (k in 1:1000) { + dattt <- data_mdc[4*(k-1)+1:4,] + perf_moreflex <- c(perf_moreflex,all(rownames(dattt[dattt$real==TRUE,])%in%rownames(dattt[dattt$logit_pred==TRUE,]))) +} + +############################################################################## +#----------------------------------------------------------------------------# +######################## FIT UNIFORMITY DETECTION MODEL ###################### +#----------------------------------------------------------------------------# +############################################################################## \ No newline at end of file