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

Experimented some MDC methods

Browse Source
This commit is contained in:
Corentin Choisy
2024-05-07 10:15:41 +02:00
parent c52800cf53
commit 49d465573f
3 changed files with 160 additions and 18 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored
View File

@ -3,3 +3,4 @@
.Rproj.user
*.RData
*.Rhistory
*.log

88
RProject/Scripts/aggregation.R
View File

@ -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,])
##############################################################################

89
RProject/Scripts/mdc.R Normal file
View File

@ -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 ######################
#----------------------------------------------------------------------------#
##############################################################################