From a87a080993d974ed1c7b4f6be4d9382dde91566b Mon Sep 17 00:00:00 2001 From: corentinchoisy Date: Thu, 21 Mar 2024 13:18:00 +0100 Subject: [PATCH] Analysis scripts --- Rproject/Rproject.Rproj | 13 ++++ Rproject/analyse_rosali.R | 143 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 156 insertions(+) create mode 100644 Rproject/Rproject.Rproj create mode 100644 Rproject/analyse_rosali.R diff --git a/Rproject/Rproject.Rproj b/Rproject/Rproject.Rproj new file mode 100644 index 0000000..5183317 --- /dev/null +++ b/Rproject/Rproject.Rproj @@ -0,0 +1,13 @@ +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: knitr +LaTeX: XeLaTeX diff --git a/Rproject/analyse_rosali.R b/Rproject/analyse_rosali.R new file mode 100644 index 0000000..564eb8c --- /dev/null +++ b/Rproject/analyse_rosali.R @@ -0,0 +1,143 @@ +############################################################################## +#----------------------------------------------------------------------------# +################################## LIBRARIES ################################# +#----------------------------------------------------------------------------# +############################################################################## + +library(TAM) +library(doMC) +library(parallel) +library(pbmcapply) +library(funprog) +library(dplyr) +library(readxl) + +lastChar <- function(str){ + substr(str, nchar(str), nchar(str)) +} + + + +############################################################################## +#----------------------------------------------------------------------------# +########################### AGGREGATION DIF MATRICES ######################### +#----------------------------------------------------------------------------# +############################################################################## + +#### Create data.frame + +results <- c(sapply(c(2,4),function(x) paste0(x,c('A','B','C'))),sapply(8,function(x) paste0(x,c('A','B','C')))) + +results2 <- c(sapply(16,function(x) paste0(x,c('A','B','C')))) + +results <- c(sapply(c(50,300),function(x) paste0(results,'_',x))) + +results2 <- c(sapply(c(50,300),function(x) paste0(results2,'_',x))) + +results <- sort(results) + +results2 <- sort(results2) + +results <- c(results,results2) + +results <- c(sapply(c('noBF','noLRT','noLRT_noBF','original'),function(x) paste0(results,'_',x))) + +#### Compiler function + +compile_simulation2 <- function(scenario) { + name <- as.numeric(gsub("[^0-9.-]", "", substr(scenario,start=0,stop=2))) + if (substr(scenario,start=nchar(scenario)-5,stop=nchar(scenario))=="0_noBF") { + s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/ROSALI-SIM/Analysis/noBF/',scenario,'.xls')) + s$version <- "noBF" + } + if (substr(scenario,start=nchar(scenario)-4,stop=nchar(scenario))=="noLRT") { + s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/ROSALI-SIM/Analysis/noLRT/',scenario,'.xls')) + s$version <- "noLRT" + } + if (substr(scenario,start=nchar(scenario)-9,stop=nchar(scenario))=="noLRT_noBF") { + s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/ROSALI-SIM/Analysis/noLRTnoBF/',scenario,'.xls')) + s$version <- "noLRT & noBF" + } + if (substr(scenario,start=nchar(scenario)-7,stop=nchar(scenario))=="original") { + s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/ROSALI-SIM/Analysis/original/',scenario,'.xls')) + s$version <- "original" + } + name <- gsub("_", "", substr(scenario,start=0,stop=3)) + if (as.numeric(gsub("[A,B,C]","",name))==2){ + colnames(s) <- c("dif.detect.1","dif.detect.2","dif.detect.3","dif.detect.4", + "N","J","nb.dif","version") + s$dif.detected <- sapply(1:1000,function(x) any(!is.na(s[x,1:4]))) + s$prop.perfect <- NA + s$prop.flexible <- NA + s$prop.moreflexible <- NA + } + if (as.numeric(gsub("[A,B,C]","",name))==4){ + colnames(s) <- c("dif.detect.1","dif.detect.2","dif.detect.3","dif.detect.4", + "dif.detect.5","dif.detect.6","dif.detect.7","N","J","nb.dif","version") + s$dif.detected <- sapply(1:1000,function(x) any(!is.na(s[x,1:7]))) + s$prop.perfect <- NA + s$prop.flexible <- NA + s$prop.moreflexible <- NA + } + if (as.numeric(gsub("[A,B,C]","",name))==8){ + colnames(s) <- c("dif.detect.1","dif.detect.2","dif.detect.3","dif.detect.4","N","J","nb.dif","true.dif","version") + s$dif.detected <- sapply(1:1000,function(x) any(!is.na(s[x,1:4]))) + s$perfect.detection <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,1:4]))==1,s[x,which(sapply(1:4,function(y) !is.na(s[x,y])))]==s[x,"true.dif"],0) ) + s$prop.perfect <- mean(s$perfect.detection) + s$flexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,1:4]))!=0,s[x,"true.dif"]%in%s[x,which(sapply(1:4,function(y) !is.na(s[x,y])))],0) ) + s$prop.flexible <- mean(s$flexible.detect) + s$prop.moreflexible <- s$prop.flexible + } + if (as.numeric(gsub("[A,B,C]","",name))==16){ + colnames(s) <- c("dif.detect.1","dif.detect.2","dif.detect.3","dif.detect.4", + "dif.detect.5","dif.detect.6","dif.detect.7","N","J","nb.dif","true.dif.1","true.dif.2","version") + s$dif.detected <- sapply(1:1000,function(x) any(!is.na(s[x,1:7]))) + perfect.detection <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,1:7]))==2,s[x,which(sapply(1:7,function(y) !is.na(s[x,y])))[1]]%in%c(s[x,c("true.dif.1","true.dif.2")]) & s[x,which(sapply(1:7,function(y) !is.na(s[x,y])))[2]]%in%c(s[x,c("true.dif.1","true.dif.2")]) + ,0) ) + s$prop.perfect <- mean(perfect.detection) + s$flexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,1:7]))!=0,s[x,"true.dif.1"]%in%c(s[x,which(sapply(1:7,function(y) !is.na(s[x,y])))]) & + s[x,"true.dif.2"]%in%c(s[x,which(sapply(1:7,function(y) !is.na(s[x,y])))]),0) ) + s$prop.flexible <- mean(s$flexible.detect) + s$moreflexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,1:7]))!=0,s[x,"true.dif.1"]%in%c(s[x,which(sapply(1:7,function(y) !is.na(s[x,y])))]) | + s[x,"true.dif.2"]%in%c(s[x,which(sapply(1:7,function(y) !is.na(s[x,y])))]),0) ) + s$prop.moreflexible <- mean(s$moreflexible.detect) + } + d <- data.frame(scenario=name, + version=unique(s$version), + N=unique(s$N), + J=unique(s$J), + M=4, + nb.dif=unique(s$nb.dif), + prop.dif.detect=mean(s$dif.detected), + prop.perfect=ifelse(substr(name,start=nchar(name),stop=nchar(name))!="D",unique(s$prop.perfect),NA), + prop.flexible=ifelse(substr(name,start=nchar(name),stop=nchar(name))!="D",unique(s$prop.flexible),NA), + prop.moreflexible=ifelse(substr(name,start=nchar(name),stop=nchar(name))!="D",unique(s$prop.moreflexible),NA) + ) + return(d) +} + + +#### Compiled results + +results <- results[-c(21,23,93,95)] + +res.dat.dif <- compile_simulation2("2A_300_noBF") + +for (x in results[seq(2,length(results))]) { + y <- compile_simulation2(x) + res.dat.dif <- bind_rows(res.dat.dif,y) +} + +res.dat.dif + + +############################################################################## +#----------------------------------------------------------------------------# +############################### PLOTS OF RESULTS ############################# +#----------------------------------------------------------------------------# +############################################################################## + +# False positives + +plot.dat <- res.dat.dif[res.dat.dif$J==4 & res.dat.dif$M==4 & res.dat.dif$nb.dif==0 & res.dat.dif$version=="original",] +boxplot(plot.dat$prop.dif.detect~plot.dat$N,ylim=c(0,0.05))