From 7454fa129531a3028afdd874402ceb68968f5b58 Mon Sep 17 00:00:00 2001
From: corentinchoisy <corentin.choisy@proton.me>
Date: Thu, 14 Mar 2024 17:02:10 +0100
Subject: [PATCH] final update

---
 RProject/.Rhistory                   | 1022 +++++++++++++-------------
 RProject/Scripts/desc_analysis_dif.R |  204 ++---
 2 files changed, 613 insertions(+), 613 deletions(-)

diff --git a/RProject/.Rhistory b/RProject/.Rhistory
index d810d6f..e2fd56d 100644
--- a/RProject/.Rhistory
+++ b/RProject/.Rhistory
@@ -1,512 +1,512 @@
-splits_evtl[[count+1]][[item]][[variable]][(knoten+1),splits_evtl[[count+1]][[item]][[variable]][(knoten+1),]<=split] <- NA
-# any split?
-anysplit <- !all(is.na(unlist(splits_evtl[[count+1]])))
-# passe vars_evtl an
-vars_evtl[[count+1]]                             <- vars_evtl[[count]]
-vars_evtl[[count+1]][[item]]                     <- rep(0,n_knots)
-vars_evtl[[count+1]][[item]][c(knoten,knoten+1)] <- rep(vars_evtl[[count]][[item]][knoten],2)
-vars_evtl[[count+1]][[item]][-c(knoten,knoten+1)]<- vars_evtl[[count]][[item]][-knoten]
-if(length(which(!is.na(splits_evtl[[count+1]][[item]][[variable]][knoten,])))==0){
-vars_evtl[[count+1]][[item]][knoten] <- vars_evtl[[count+1]][[item]][knoten]-1
-}
-if(length(which(!is.na(splits_evtl[[count+1]][[item]][[variable]][knoten+1,])))==0){
-vars_evtl[[count+1]][[item]][knoten+1] <- vars_evtl[[count+1]][[item]][knoten+1]-1
-}
-# passe which_obs an
-which_obs[[count+1]]                               <- which_obs[[count]]
-which_obs[[count+1]][[item]]                       <- matrix(0,nrow=n_knots,ncol=npersons)
-which_obs[[count+1]][[item]][c(knoten,knoten+1),]  <- matrix(rep(which_obs[[count]][[item]][knoten,],2),nrow=2,byrow=T)
-which_obs[[count+1]][[item]][-c(knoten,knoten+1),] <- which_obs[[count]][[item]][-knoten,]
-thresh <- ordered_values[[variable]][1:n_s[variable]][split]
-which_obs[[count+1]][[item]][knoten,DM_kov[,variable]>thresh] <- NA
-which_obs[[count+1]][[item]][(knoten+1),DM_kov[,variable]<=thresh] <- NA
-# passe numbers an
-numbers[[count+1]]                              <- numbers[[count]]
-numbers[[count+1]][[item]]                      <- numeric(length=n_knots)
-numbers[[count+1]][[item]][c(knoten,knoten+1)]  <- c(left,right)
-numbers[[count+1]][[item]][-c(knoten,knoten+1)] <- numbers[[count]][[item]][-knoten]
-# trace
-if(trace){
-cat(paste0("\n Split"," ",count,";"," ","Item"," ",item,"\n"))
-}
-# erhoehe counter
-count <- count+1
-} else{
-sig <- FALSE
-}
-}
-###################################################################################
-# prettify results
-mod_opt     <- mod_potential[[count]]
-ip_opt      <- names(coef(mod_opt))[-c(1:(npersons-1))]
-theta_hat   <- c(coef(mod_opt)[1:(npersons-1)],0)
-delta_hat   <- coef(mod_opt)[npersons:length(coef(mod_opt))]
-if(count>1){
-dif_items   <- unique(splits[,2])
-nodif_items <- c(1:nitems)[-dif_items]
-delta_hat_nodif <- sapply(nodif_items,function(j) delta_hat[grep(paste0("delta",j,":"),ip_opt)])
-rownames(delta_hat_nodif) <- 1:nrow(delta_hat_nodif)
-colnames(delta_hat_nodif) <- paste0("delta", nodif_items)
-delta_hat_dif   <- lapply(dif_items, function(j)  delta_hat[grep(paste0("delta",j,":"),ip_opt)])
-names(delta_hat_dif) <- dif_items
-help9 <- cumsum(c(0,(n_levels-1)))
-colnames(splits) <- c("var","item","split","level","node","number","left","right")
-splits <- data.frame(cbind(splits[,1:5,drop=FALSE],"variable"=rep(NA,nrow(splits)),"threshold"=rep(NA,nrow(splits)),splits[,6:8,drop=FALSE]))
-for(i in 1:nrow(splits)){
-if(!is.null(colnames(DM_kov))){
-splits[i,6] <- colnames(DM_kov)[splits[i,1]]
-} else{
-splits[i,6] <- splits[i,1]
-}
-v2 <- lapply(1:nvar,function(j) ordered_values[[j]][-length(ordered_values[[j]])])
-splits[i,7] <- v2[[splits[i,1]]][splits[i,3]]
-}
-splits <- splits[,-1]
-nthres <- length(unique(y))-1
-for(i in dif_items){
-info <- splits[splits[,"item"]==i,]
-endnodes <- get_endnodes(info)
-names(delta_hat_dif[[paste(i)]]) <- paste(rep(endnodes,each=nthres),rep(1:nthres,length(endnodes)),sep=":")
-}
-} else{
-delta_hat_nodif <- sapply(1:nitems,function(j) delta_hat[grep(paste0("delta",j,":"),ip_opt)])
-delta_hat_dif   <- c()
-}
-to_return <- list("splits"=splits,
-"thetas"=theta_hat,
-"deltas_nodif"=delta_hat_nodif,
-"deltas_dif"=delta_hat_dif,
-"pvalues"=pvalues,
-"devs"=devs,
-"crits"=crits)
-return(to_return)
-}
-DIFtree(Y=dat1[,c('item1','item2','item3','item4','item5','item6','item7')],X=as.data.frame(dat1[,'TT']),model="PCM")
-DIFtree
-View(DIFtree)
-DIFtree <- function (Y, X, model = c("Rasch", "Logistic", "PCM"), type = c("udif",
-"dif", "nudif"), alpha = 0.05, nperm = 1000, trace = FALSE,
-penalize = FALSE, ...)
-{
-UseMethod("DIFtree")
-}
-DIFtree <- function (Y, X, model = c("Rasch", "Logistic", "PCM"), type = c("udif",
-"dif", "nudif"), alpha = 0.05, nperm = 1000, trace = FALSE,
-penalize = FALSE, ...)
-{
-browser()
-UseMethod("DIFtree")
-}
-DIFtree(Y=dat1[,c('item1','item2','item3','item4','item5','item6','item7')],X=as.data.frame(dat1[,'TT']),model="PCM")
-library(DIFtree)
-dat1 <- read.csv(file = '/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N300/scenario_13A_300.csv')
-dat1 <- dat1[dat1$replication==1]
-dat1 <- dat1[dat1$replication==1,]
-DIFtree(Y=dat1[,c('item1','item2','item3','item4','item5','item6','item7')],X=as.data.frame(dat1[,'TT']),model="PCM")
-function (Y, X, model = c("Rasch", "Logistic", "PCM"), type = c("udif",
-"dif", "nudif"), alpha = 0.05, nperm = 1000, trace = FALSE,
-penalize = FALSE, ...)
-{
-UseMethod("DIFtree")
-}
-DIFtree <- function (Y, X, model = c("Rasch", "Logistic", "PCM"), type = c("udif",
-"dif", "nudif"), alpha = 0.05, nperm = 1000, trace = FALSE,
-penalize = FALSE, ...)
-{
-UseMethod("DIFtree")
-}
-DIFtree(Y=dat1[,c('item1','item2','item3','item4','item5','item6','item7')],X=as.data.frame(dat1[,'TT']),model="PCM")
-library(DIFtree)
-dat1 <- read.csv(file = '/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N300/scenario_13A_300.csv')
-dat1 <- dat1[dat1$replication==1,]
-DIFtree(Y=dat1[,c('item1','item2','item3','item4','item5','item6','item7')],X=as.data.frame(dat1[,'TT']),model="PCM")
-DIFtree(Y=dat1[,c('item1','item2','item3','item4','item5','item6','item7')],X=as.data.frame(dat1[,'TT']),model="PCM",trace=T)
-dat1 <- read.csv(file = '/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N300/scenario_14A_300.csv')
-dat1 <- dat1[dat1$replication==1,]
-DIFtree(Y=dat1[,c('item1','item2','item3','item4','item5','item6','item7')],X=as.data.frame(dat1[,'TT']),model="PCM",trace=T)
-dat1$item1 <- as.factor(dat1$item1)
-dat1$item2 <- as.factor(dat1$item2)
-dat1$item3 <- as.factor(dat1$item3)
-dat1$item4 <- as.factor(dat1$item4)
-dat1$item5 <- as.factor(dat1$item5)
-dat1$item6 <- as.factor(dat1$item6)
-dat1$item7 <- as.factor(dat1$item7)
-DIFtree(Y=dat1[,c('item1','item2','item3','item4','item5','item6','item7')],X=as.data.frame(dat1[,'TT']),model="PCM",trace=T)
-data("data_sim_PCM")
-data_sim_PCM
-dat1 <- read.csv(file = '/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N300/scenario_14A_300.csv')
-dat1 <- dat1[dat1$replication==1,]
-dat1
-DIFtree(Y=dat1[,c('item1','item2','item3','item4','item5','item6','item7')],X=as.data.frame(dat1[,'TT']),model="PCM",trace=T)
-dat1[,c('item1','item2','item3','item4','item5','item6','item7')]
-dat1[,c('item1','item2','item3','item4','item5','item6','item7')] <- dat1[,c('item1','item2','item3','item4','item5','item6','item7')]+1
-dat1[,c('item1','item2','item3','item4','item5','item6','item7')]
-DIFtree(Y=dat1[,c('item1','item2','item3','item4','item5','item6','item7')],X=as.data.frame(dat1[,'TT']),model="PCM",trace=T)
-DIFtree(Y=dat1[,c('item1')],X=as.data.frame(dat1[,'TT']),model="PCM",trace=T)
-dat1[,c('item1')]
-dat1[,c('TT')]
-DIFtree(Y=dat1[,c('item1')],X=as.data.frame(dat1[,'TT']),model="PCM",trace=T)
-as.data.frame(dat1[,'TT'])
-DIFtree(Y=data.matrix(dat1[,c('item1','item2','item3','item4','item5','item6','item7')]),X=as.matrix(dat1[,'TT']),model="PCM",trace=T)
-DIFtree(Y=data.matrix(dat1[,c('item1','item2','item3','item4','item5','item6','item7')]),X=data.frame(x1=dat1[,'TT']),model="PCM",trace=T)
-data.matrix(dat1[,c('item1','item2','item3','item4','item5','item6','item7')])
-Y2 <- data_sim_PCM[,1]
-X2 <- data_sim_PCM[,-1]
-Y2
-X2
-DIFtree(Y=as.matrix(dat1[,c('item1','item2','item3','item4','item5','item6','item7')]),X=data.frame(x1=dat1[,'TT']),model="PCM",trace=T)
-mod2 <- DIFtree(Y=Y2,X=X2,model="PCM",alpha=0.05,nperm=100,trace=TRUE)
-library(DIFtree)
-data(data_sim_Rasch)
-data(data_sim_PCM)
-Y1 <- data_sim_Rasch[,1]
-X1 <- data_sim_Rasch[,-1]
-Y2 <- data_sim_PCM[,1]
-X2 <- data_sim_PCM[,-1]
-## Not run:
-mod1 <- DIFtree(Y=Y1,X=X1,model="Logistic",type="udif",alpha=0.05,nperm=1000,trace=TRUE)
-mod2 <- DIFtree(Y=Y2,X=X2,model="PCM",alpha=0.05,nperm=100,trace=TRUE)
-remove.packages("DIFtree")
-library(devtools)
-install_version("DIFtree","3.1.4")
-library(DIFtree)
-data(data_sim_Rasch)
-data(data_sim_PCM)
-Y1 <- data_sim_Rasch[,1]
-X1 <- data_sim_Rasch[,-1]
-Y2 <- data_sim_PCM[,1]
-X2 <- data_sim_PCM[,-1]
-## Not run:
-mod1 <- DIFtree(Y=Y1,X=X1,model="Logistic",type="udif",alpha=0.05,nperm=1000,trace=TRUE)
-mod2 <- DIFtree(Y=Y2,X=X2,model="PCM",alpha=0.05,nperm=100,trace=TRUE)
-remove.packages("DIFtree")
-install_version("DIFtree","2.0.1")
-library(DIFtree)
-data(data_sim_Rasch)
-data(data_sim_PCM)
-Y1 <- data_sim_Rasch[,1]
-X1 <- data_sim_Rasch[,-1]
-Y2 <- data_sim_PCM[,1]
-X2 <- data_sim_PCM[,-1]
-## Not run:
-mod1 <- DIFtree(Y=Y1,X=X1,model="Logistic",type="udif",alpha=0.05,nperm=1000,trace=TRUE)
-mod2 <- DIFtree(Y=Y2,X=X2,model="PCM",alpha=0.05,nperm=100,trace=TRUE)
-tree_PCM <-
-function(y,
-DM_kov,
-npersons,
-nitems,
-nvar,
-ordered_values,
-n_levels,
-n_s,
-alpha,
-nperm,
-trace
-){
-# design of PCM
-pp_design <- diag(npersons)                   # persons, person P reference
-pp_design <- pp_design[rep(1:nrow(pp_design),each=nitems),]
-pp_design <- pp_design[,-npersons]
-ip_design <- -1*diag(nitems)                  # item parameter
-ip_design <- ip_design[rep(1:nrow(ip_design),times=npersons),]
-dm_pcm <- cbind(pp_design,ip_design)
-names_pcm    <- c(paste("theta",1:(npersons-1),sep=""),paste("delta",1:nitems,sep=""))
-colnames(dm_pcm) <- names_pcm
-# functions to build design
-thresholds <- lapply(1:nvar, function(j) ordered_values[[j]][-length(ordered_values[[j]])])
-v <- lapply(1:nvar,function(j) 1:(n_levels[j]-1))
-w <- lapply(1:nvar, function(j) rep(paste0("s",j),n_s[j]))
-design_one  <- function(x,threshold,upper){
-if(upper){
-ret <- ifelse(x > threshold,1,0)
-} else{
-ret <- ifelse(x > threshold,0,1)
-}
-return(ret)
-}
-design <- function(x,thresholds,upper){
-ret <- sapply(thresholds, function(j) design_one(x,j,upper))
-return(ret)
-}
-whole_design <- function(X,var,item,thresholds,upper=TRUE){
-design_tree <- matrix(0,nrow=nitems*npersons,ncol=length(thresholds[[var]]))
-rows        <- seq(item,(nitems*npersons),by=nitems)
-design_tree[rows,] <- design(X[,var],thresholds[[var]],upper)
-z <- rep(paste0(ifelse(upper,"_u","_l"),item),length(thresholds[[var]]))
-colnames(design_tree) <- paste0(w[[var]],v[[var]],z)
-return(design_tree)
-}
-designlists <- function(X,thresholds,upper=TRUE){
-ret <- lapply(1:nitems, function(j){
-lapply(1:nvar, function(var){
-whole_design(X,var,j,thresholds,upper)
-})
-})
-return(ret)
-}
-#########################################################################################
-mod_potential <- list()
-devs          <- c()
-crits          <- c()
-splits        <- c()
-pvalues       <- c()
-ip            <- list()
-vars_evtl     <- list()
-splits_evtl   <- list()
-which_obs     <- list()
-numbers       <- list()
-count <- 1
-numbers[[1]]       <- lapply(1:nitems,function(j) 1)
-which_obs[[1]]     <- lapply(1:nitems,function(j) matrix(1:npersons,nrow=1))
-splits_evtl[[1]]   <- lapply(1:nitems,function(j) lapply(1:nvar, function(var) matrix(1:n_s[var],nrow=1)))
-vars_evtl[[1]]     <- lapply(1:nitems,function(j) nvar)
-ip[[1]]            <- lapply(1:nitems,function(j) paste0("delta",j))
-### PCM ###
-pp     <- paste("theta",1:(npersons-1),sep="")
-help_p <- paste0(pp,collapse="+")
-help01           <- formula(paste("y~",help_p,"+",paste0(unlist(ip[[1]]),collapse="+"),"-1"))
-help02           <- formula(paste0("FALSE~",paste0(unlist(ip[[1]]),collapse="+")))
-dat0             <- data.frame(y,dm_pcm)
-mod0             <- tryCatch(vglm(help01,
-family=acat(parallel=help02, reverse=FALSE),
-data=dat0,
-na.action=na.omit,
-checkwz=FALSE),
-error = function(e) stop("PCM not identified!", call. =FALSE))
-start              <- VGAM::predict(mod0)
-mod_potential[[1]] <- mod0
-design_upper <- designlists(DM_kov,thresholds)
-design_lower <- designlists(DM_kov,thresholds,upper=FALSE)
-sig   <- TRUE
-anysplit <- TRUE
-# function to compute all models in one knot
-allmodels <- function(i,var,kn,design_lower,design_upper){
-deviances <- rep(0,n_s[var])
-help_kn <- ip[[count]][[i]][kn]
-help1   <- paste0(unlist(ip[[count]])[-which(unlist(ip[[count]])==help_kn)],collapse="+")
-splits_aktuell <- splits_evtl[[count]][[i]][[var]][kn,]
-splits_aktuell <- splits_aktuell[!is.na(splits_aktuell)]
-obs0 <- which(!is.na(which_obs[[count]][[i]][kn,]))
-if(length(splits_aktuell)>0){
-for(j in splits_aktuell){
-n_lower <- sum(DM_kov[obs0,var]<=ordered_values[[var]][j])
-n_upper <- sum(DM_kov[obs0,var]>ordered_values[[var]][j])
-if(n_lower>=30 & n_upper>=30){
-dat   <- data.frame(dat0,design_lower[[i]][[var]][,j,drop=FALSE],design_upper[[i]][[var]][,j,drop=FALSE])
-help2 <- paste(ip[[count]][[i]][kn],c(colnames(design_lower[[i]][[var]])[j],colnames(design_upper[[i]][[var]])[j]),sep=":")
-help3 <- paste(help2,collapse="+")
-help41 <- formula(paste("y~",help1,"+",help3,"-1"))
-help42 <- formula(paste0("FALSE~",help1,"+",help3))
-suppressWarnings(
-mod <- try(vglm(help41,
-family=acat(parallel=help42, reverse=FALSE),
-data=dat,
-checkwz=FALSE,
-na.action=na.omit,
-offset=start))
-)
-if(class(mod)!="try-error"){
-deviances[j] <-  deviance(mod0)-deviance(mod)
-}
-}
-}
-}
-return(deviances)
-}
-# estimate tree
-while(sig & anysplit){
-# compute all models
-dv <- lapply(1:nvar,function(var) {
-lapply(1:nitems,function(i) {
-n_knots   <- length(ip[[count]][[i]])
-deviances <- matrix(rep(0,n_s[var]*n_knots),ncol=n_knots)
-for(kn in 1:n_knots){
-deviances[,kn] <- allmodels(i,var,kn,design_lower,design_upper)
-}
-return(deviances)
-})
-})
-# select optimum
-variable <- which.max(lapply(1:nvar,function(j) max(unlist(dv[[j]]))))
-item     <- which.max(lapply(1:nitems, function(j) max(dv[[variable]][[j]])))
-split    <- as.numeric(which(dv[[variable]][[item]]==max(dv[[variable]][[item]]),arr.ind=TRUE)[,1])
-knoten   <- as.numeric(which(dv[[variable]][[item]]==max(dv[[variable]][[item]]),arr.ind=TRUE)[,2])
-if(length(split)>1){
-split  <- split[1]
-knoten <- knoten[1]
-warning(paste("Maximum in iteration ",count," not uniquely defined"))
-}
-ip_old   <- ip[[count]][[item]][knoten]
-level    <- length(strsplit(ip_old,":")[[1]])
-number   <- numbers[[count]][[item]][knoten]
-left     <- max(numbers[[count]][[item]])+1
-right    <- max(numbers[[count]][[item]])+2
-# compute permutation test
-dev <- rep(NA,nperm)
-for(perm in 1:nperm){
-dv_perm <- rep(0,n_s[variable])
-obs_aktuell <- which_obs[[count]][[item]][knoten,]
-obs_aktuell <- obs_aktuell[!is.na(obs_aktuell)]
-DM_kov_perm <- DM_kov
-DM_kov_perm[obs_aktuell,variable] <- sample(DM_kov_perm[obs_aktuell,variable],length(obs_aktuell))
-design_upper_perm      <- design_upper
-design_upper_perm[[item]][[variable]] <- whole_design(DM_kov_perm,variable,item,thresholds)
-design_lower_perm      <- design_lower
-design_lower_perm[[item]][[variable]] <- whole_design(DM_kov_perm,variable,item,thresholds,upper=FALSE)
-dv_perm <- allmodels(item,variable,knoten,design_lower_perm,design_upper_perm)
-dev[perm] <- max(dv_perm)
-if(trace){
-cat(".")
-}
-}
-# test decision
-crit_val <- quantile(dev,1-(alpha/vars_evtl[[count]][[item]][knoten]))
-proof    <- max(dv[[variable]][[item]]) > crit_val
-devs[count]   <- max(dv[[variable]][[item]])
-crits[count]   <- crit_val
-pvalues[count] <- length(which(dev>max(dv[[variable]][[item]])))/nperm
-if(proof){
-# get new formula
-help_kn2 <- ip[[count]][[item]][knoten]
-help5    <- paste0(unlist(ip[[count]])[-which(unlist(ip[[count]])==help_kn2)],collapse="+")
-help6    <- paste(ip[[count]][[item]][knoten],c(colnames(design_lower[[item]][[variable]])[split],colnames(design_upper[[item]][[variable]])[split]),sep=":")
-help7    <- paste(help6,collapse="+")
-help81   <- formula(paste("y~",help_p,"+",help5,"+",help7,"-1"))
-help82   <- formula(paste0("FALSE~",help5,"+",help7))
-######################
-if(level>1){
-help_kn4 <- lu(c(),1,level-1,c())
-help_kn5 <- unlist(strsplit(help_kn2,""))
-help_kn6 <- paste0(help_kn5[which(help_kn5=="_")+1],collapse="")
-knoten2  <- which(help_kn4==help_kn6)
-} else{
-knoten2 <- knoten
-}
-######################
-splits <- rbind(splits,c(variable,item,split,level,knoten2,number,left,right))
-# fit new model
-dat <- dat0 <- data.frame(dat0,design_lower[[item]][[variable]][,split,drop=FALSE],design_upper[[item]][[variable]][,split,drop=FALSE])
-suppressWarnings(
-mod0  <- mod_potential[[count+1]] <-  tryCatch(vglm(help81,
-family=acat(parallel=help82, reverse=FALSE),
-data=dat,
-na.action=na.omit,
-checkwz=FALSE,
-etastart=start),
-error = function(e) stop("IFT_PCM not identified!", call. =FALSE))
-)
-start <- VGAM::predict(mod0)
-# generiere neue itemparameter
-ip[[count+1]]                             <- ip[[count]]
-ip[[count+1]][[item]]                     <- rep("",length(ip[[count]][[item]])+1)
-ip[[count+1]][[item]][c(knoten,knoten+1)] <- help6
-ip[[count+1]][[item]][-c(knoten,knoten+1)]<- ip[[count]][[item]][-knoten]
-# passe splits_evtl an
-n_knots                                                       <- length(ip[[count+1]][[item]])
-splits_evtl[[count+1]]                                        <- splits_evtl[[count]]
-for(var in 1:nvar){
-splits_evtl[[count+1]][[item]][[var]]                       <- matrix(0,nrow=n_knots,ncol=n_s[var])
-splits_evtl[[count+1]][[item]][[var]][c(knoten,knoten+1),]  <- matrix(rep(splits_evtl[[count]][[item]][[var]][knoten,],2),nrow=2,byrow=T)
-splits_evtl[[count+1]][[item]][[var]][-c(knoten,knoten+1),] <- splits_evtl[[count]][[item]][[var]][-knoten,]
-}
-splits_evtl[[count+1]][[item]][[variable]][knoten,splits_evtl[[count+1]][[item]][[variable]][knoten,]>=split] <- NA
-splits_evtl[[count+1]][[item]][[variable]][(knoten+1),splits_evtl[[count+1]][[item]][[variable]][(knoten+1),]<=split] <- NA
-# any split?
-anysplit <- !all(is.na(unlist(splits_evtl[[count+1]])))
-# passe vars_evtl an
-vars_evtl[[count+1]]                             <- vars_evtl[[count]]
-vars_evtl[[count+1]][[item]]                     <- rep(0,n_knots)
-vars_evtl[[count+1]][[item]][c(knoten,knoten+1)] <- rep(vars_evtl[[count]][[item]][knoten],2)
-vars_evtl[[count+1]][[item]][-c(knoten,knoten+1)]<- vars_evtl[[count]][[item]][-knoten]
-if(length(which(!is.na(splits_evtl[[count+1]][[item]][[variable]][knoten,])))==0){
-vars_evtl[[count+1]][[item]][knoten] <- vars_evtl[[count+1]][[item]][knoten]-1
-}
-if(length(which(!is.na(splits_evtl[[count+1]][[item]][[variable]][knoten+1,])))==0){
-vars_evtl[[count+1]][[item]][knoten+1] <- vars_evtl[[count+1]][[item]][knoten+1]-1
-}
-# passe which_obs an
-which_obs[[count+1]]                               <- which_obs[[count]]
-which_obs[[count+1]][[item]]                       <- matrix(0,nrow=n_knots,ncol=npersons)
-which_obs[[count+1]][[item]][c(knoten,knoten+1),]  <- matrix(rep(which_obs[[count]][[item]][knoten,],2),nrow=2,byrow=T)
-which_obs[[count+1]][[item]][-c(knoten,knoten+1),] <- which_obs[[count]][[item]][-knoten,]
-thresh <- ordered_values[[variable]][1:n_s[variable]][split]
-which_obs[[count+1]][[item]][knoten,DM_kov[,variable]>thresh] <- NA
-which_obs[[count+1]][[item]][(knoten+1),DM_kov[,variable]<=thresh] <- NA
-# passe numbers an
-numbers[[count+1]]                              <- numbers[[count]]
-numbers[[count+1]][[item]]                      <- numeric(length=n_knots)
-numbers[[count+1]][[item]][c(knoten,knoten+1)]  <- c(left,right)
-numbers[[count+1]][[item]][-c(knoten,knoten+1)] <- numbers[[count]][[item]][-knoten]
-# trace
-if(trace){
-cat(paste0("\n Split"," ",count,";"," ","Item"," ",item,"\n"))
-}
-# erhoehe counter
-count <- count+1
-} else{
-sig <- FALSE
-}
-}
-###################################################################################
-# prettify results
-mod_opt     <- mod_potential[[count]]
-ip_opt      <- names(coef(mod_opt))[-c(1:(npersons-1))]
-theta_hat   <- c(coef(mod_opt)[1:(npersons-1)],0)
-delta_hat   <- coef(mod_opt)[npersons:length(coef(mod_opt))]
-if(count>1){
-dif_items   <- unique(splits[,2])
-nodif_items <- c(1:nitems)[-dif_items]
-delta_hat_nodif <- sapply(nodif_items,function(j) delta_hat[grep(paste0("delta",j,":"),ip_opt)])
-rownames(delta_hat_nodif) <- 1:nrow(delta_hat_nodif)
-colnames(delta_hat_nodif) <- paste0("delta", nodif_items)
-delta_hat_dif   <- lapply(dif_items, function(j)  delta_hat[grep(paste0("delta",j,":"),ip_opt)])
-names(delta_hat_dif) <- dif_items
-help9 <- cumsum(c(0,(n_levels-1)))
-colnames(splits) <- c("var","item","split","level","node","number","left","right")
-splits <- data.frame(cbind(splits[,1:5,drop=FALSE],"variable"=rep(NA,nrow(splits)),"threshold"=rep(NA,nrow(splits)),splits[,6:8,drop=FALSE]))
-for(i in 1:nrow(splits)){
-if(!is.null(colnames(DM_kov))){
-splits[i,6] <- colnames(DM_kov)[splits[i,1]]
-} else{
-splits[i,6] <- splits[i,1]
-}
-v2 <- lapply(1:nvar,function(j) ordered_values[[j]][-length(ordered_values[[j]])])
-splits[i,7] <- v2[[splits[i,1]]][splits[i,3]]
-}
-splits <- splits[,-1]
-nthres <- length(unique(y))-1
-for(i in dif_items){
-info <- splits[splits[,"item"]==i,]
-endnodes <- get_endnodes(info)
-names(delta_hat_dif[[paste(i)]]) <- paste(rep(endnodes,each=nthres),rep(1:nthres,length(endnodes)),sep=":")
-}
-} else{
-delta_hat_nodif <- sapply(1:nitems,function(j) delta_hat[grep(paste0("delta",j,":"),ip_opt)])
-delta_hat_dif   <- c()
-}
-to_return <- list("splits"=splits,
-"thetas"=theta_hat,
-"deltas_nodif"=delta_hat_nodif,
-"deltas_dif"=delta_hat_dif,
-"pvalues"=pvalues,
-"devs"=devs,
-"crits"=crits)
-return(to_return)
-}
-library(VGAM)
-mod2 <- DIFtree(Y=Y2,X=X2,model="PCM",alpha=0.05,nperm=100,trace=TRUE)
-colnames(res.dat.dif)
-max(res.dat[res.dat$scenario.type=="A",]$h0.rejected.p)
-max(1-res.dat[res.dat$scenario.type!="A",]$beta.same.sign.truebeta.signif.p)
-max(res.dat[res.dat$scenario.type!="A",]$beta.same.sign.truebeta.signif.p)
-res.dat[res.dat$scenario.type!="A",]$beta.same.sign.truebeta.signif.p
+points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
+points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='#053305',pch=3)
+# J=7 / 2 DIF
+res.dat.temp <- res.dat[res.dat$J==7 & res.dat$nb.dif==2,]
+boxplot(true.value.in.ci.p~dif.size,data=res.dat.temp,col=c(2,3,3,2),xlab='DIF size',
+ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.5,]
+points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.5,]
+points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.3,]
+points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
+points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='#053305',pch=3)
+# J=7 / 3 DIF
+res.dat.temp <- res.dat[res.dat$J==7 & res.dat$nb.dif==3,]
+boxplot(true.value.in.ci.p~dif.size,data=res.dat.temp,col=c(2,3,3,2),xlab='DIF size',
+ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.5,]
+points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.5,]
+points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.3,]
+points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
+points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='#053305',pch=3)
+##############################################################################
+#----------------------------------------------------------------------------#
+########################## BETA SIGN CHANGE BOXPLOTS #########################
+#----------------------------------------------------------------------------#
+##############################################################################
+# Overall
+boxplot(beta.same.sign.truebeta.p~dif.size,data=res.dat,col=c(2,3),xlab='DIF size',
+ylab='Proportion of estimates with same sign as true value in target scenario',main='DIF on 3 items',ylim=c(0,1))
+res.null3 <- res.dat[res.dat$dif.size==-0.5,]
+points(y=res.null3$beta.same.sign.truebeta.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
+res.null3 <- res.dat[res.dat$dif.size==0.5,]
+points(y=res.null3$beta.same.sign.truebeta.p,x=rep(5,nrow(res.null3)),col='#590b0c',pch=3)
+res.null3 <- res.dat[res.dat$dif.size==-0.3,]
+points(y=res.null3$beta.same.sign.truebeta.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
+res.null3 <- res.dat[res.dat$dif.size==0.3,]
+points(y=res.null3$beta.same.sign.truebeta.p,x=rep(4,nrow(res.null3)),col='#053305',pch=3)
+res.null3 <- res.dat[res.dat$dif.size==0,]
+points(y=res.null3$beta.same.sign.truebeta.p,x=rep(3,nrow(res.null3)),col='gray',pch=3)
+# Overall // H0 rejected
+boxplot(beta.same.sign.truebeta.signif.p~dif.size,data=res.dat,col=c(2,3),xlab='DIF size',
+ylab='Proportion of estimates with same sign as true value in target scenario',main='When H0 rejected',ylim=c(0,1))
+res.null3 <- res.dat[res.dat$dif.size==-0.5,]
+points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
+res.null3 <- res.dat[res.dat$dif.size==0.5,]
+points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(5,nrow(res.null3)),col='#590b0c',pch=3)
+res.null3 <- res.dat[res.dat$dif.size==-0.3,]
+points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
+res.null3 <- res.dat[res.dat$dif.size==0.3,]
+points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(4,nrow(res.null3)),col='#053305',pch=3)
+res.null3 <- res.dat[res.dat$dif.size==0,]
+points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(3,nrow(res.null3)),col='gray',pch=3)
+## Plot 1 - baseline scenarios vs theoretical power
+par(mfrow=c(1,2))
+# theoretical
+plot.dat <- res.dat[res.dat$scenario %in% sapply(c('A','B','C','D','E'),function(x) paste0(1:4,x)),]
+plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),rep(unique(plot.dat[plot.dat$scenario.type=="A",]$theoretical.power),3),
+xlab='N',ylab='Theoretical power',
+xaxt='n',yaxt='n',type='l',xlim=c(100,300),ylim=c(0,1),lty=4,col='#03a18a')
+unique(plot.dat[plot.dat$scenario.type=="A",]$N)
 View(res.dat)
-res.dat[res.dat$beta.same.sign.truebeta.p==NaN & res.dat$eff.size!=0,]
-res.dat[res.dat$beta.same.sign.truebeta.p==NaN & res.dat$eff.size!=0,]$N
-res.dat[is.nan(res.dat$beta.same.sign.truebeta.p) & res.dat$eff.size!=0,]$N
-res.dat[res.dat$N==50,]
-res.dat[res.dat$N==50,]$eff.size
-max(res.dat[res.dat$scenario.type!="A" & res.dat$N!=50,]$beta.same.sign.truebeta.signif.p)
-max(1-res.dat[res.dat$scenario.type!="A" & res.dat$N!=50,]$beta.same.sign.truebeta.signif.p)
-max(1-res.dat[res.dat$scenario.type!="A" & res.dat$N!=50,]$beta.same.sign.truebeta.p)
+##############################################################################
+#----------------------------------------------------------------------------#
+################################# POWER PLOTS ################################
+#----------------------------------------------------------------------------#
+##############################################################################
+## Plot 1 - baseline scenarios vs theoretical power
+par(mfrow=c(1,2))
+# theoretical
+plot.dat <- res.dat[res.dat$scenario %in% sapply(c('A','B','C','D','E'),function(x) paste0(1:4,x)),]
+plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),rep(unique(plot.dat[plot.dat$scenario.type=="A",]$theoretical.power),4),
+xlab='N',ylab='Theoretical power',
+xaxt='n',yaxt='n',type='l',xlim=c(50,300),ylim=c(0,1),lty=4,col='#03a18a')
+axis(1,c(50,100,200,300))
+axis(2,seq(0,1,0.1))
+lines(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$theoretical.power),col='#c6d18d',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$theoretical.power),col='#c0c23b',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$theoretical.power),col='#da77c7',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$theoretical.power),col='#a12471',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$theoretical.power),col='#b5a180',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$theoretical.power),col='#9b6541',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$theoretical.power),col='#30a466',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$theoretical.power),col='#1a342b',lty=4)
+points(unique(plot.dat[plot.dat$scenario %in% c("1A",'1A'),]$N),
+rep(unique(plot.dat[plot.dat$scenario %in% c("1A",'1A'),]$theoretical.power),3),col='#03a18a',pch=17)
+##############################################################################
+#----------------------------------------------------------------------------#
+################################# POWER PLOTS ################################
+#----------------------------------------------------------------------------#
+##############################################################################
+## Plot 1 - baseline scenarios vs theoretical power
+par(mfrow=c(1,2))
+# theoretical
+plot.dat <- res.dat[res.dat$scenario %in% sapply(c('A','B','C','D','E'),function(x) paste0(1:4,x)),]
+plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),rep(unique(plot.dat[plot.dat$scenario.type=="A",]$theoretical.power),4),
+xlab='N',ylab='Theoretical power',
+xaxt='n',yaxt='n',type='l',xlim=c(50,300),ylim=c(0,1),lty=4,col='#03a18a')
+axis(1,c(50,100,200,300))
+axis(2,seq(0,1,0.1))
+lines(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$theoretical.power),col='#c6d18d',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$theoretical.power),col='#c0c23b',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$theoretical.power),col='#da77c7',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$theoretical.power),col='#a12471',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$theoretical.power),col='#b5a180',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$theoretical.power),col='#9b6541',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$theoretical.power),col='#30a466',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$theoretical.power),col='#1a342b',lty=4)
+points(unique(plot.dat[plot.dat$scenario %in% c("1A",'1A'),]$N),
+rep(unique(plot.dat[plot.dat$scenario %in% c("1A",'1A'),]$theoretical.power),3),col='#03a18a',pch=17)
+## Plot 1 - baseline scenarios vs theoretical power
+par(mfrow=c(1,2))
+# theoretical
+plot.dat <- res.dat[res.dat$scenario %in% sapply(c('A','B','C','D','E'),function(x) paste0(1:4,x)),]
+plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),rep(unique(plot.dat[plot.dat$scenario.type=="A",]$theoretical.power),4),
+xlab='N',ylab='Theoretical power',
+xaxt='n',yaxt='n',type='l',xlim=c(50,300),ylim=c(0,1),lty=4,col='#03a18a')
+axis(1,c(50,100,200,300))
+axis(2,seq(0,1,0.1))
+lines(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$theoretical.power),col='#c6d18d',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$theoretical.power),col='#c0c23b',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$theoretical.power),col='#da77c7',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$theoretical.power),col='#a12471',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$theoretical.power),col='#b5a180',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$theoretical.power),col='#9b6541',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$theoretical.power),col='#30a466',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$theoretical.power),col='#1a342b',lty=4)
+points(unique(plot.dat[plot.dat$scenario %in% c("1A",'1A'),]$N),
+rep(unique(plot.dat[plot.dat$scenario %in% c("1A",'1A'),]$theoretical.power),3),col='#03a18a',pch=17)
+##############################################################################
+#----------------------------------------------------------------------------#
+################################# POWER PLOTS ################################
+#----------------------------------------------------------------------------#
+##############################################################################
+## Plot 1 - baseline scenarios vs theoretical power
+par(mfrow=c(1,2))
+# theoretical
+plot.dat <- res.dat[res.dat$scenario %in% sapply(c('A','B','C','D','E'),function(x) paste0(1:4,x)),]
+plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),rep(unique(plot.dat[plot.dat$scenario.type=="A",]$theoretical.power),4),
+xlab='N',ylab='Theoretical power',
+xaxt='n',yaxt='n',type='l',xlim=c(50,300),ylim=c(0,1),lty=4,col='#03a18a')
+axis(1,c(50,100,200,300))
+axis(2,seq(0,1,0.1))
+lines(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$theoretical.power),col='#c6d18d',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$theoretical.power),col='#c0c23b',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$theoretical.power),col='#da77c7',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$theoretical.power),col='#a12471',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$theoretical.power),col='#b5a180',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$theoretical.power),col='#9b6541',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$theoretical.power),col='#30a466',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$theoretical.power),col='#1a342b',lty=4)
+points(unique(plot.dat[plot.dat$scenario %in% c("1A",'1A'),]$N),
+rep(unique(plot.dat[plot.dat$scenario %in% c("1A",'1A'),]$theoretical.power),3),col='#03a18a',pch=17)
+## Plot 1 - baseline scenarios vs theoretical power
+par(mfrow=c(1,2))
+# theoretical
+plot.dat <- res.dat[res.dat$scenario %in% sapply(c('A','B','C','D','E'),function(x) paste0(1:4,x)),]
+plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),rep(unique(plot.dat[plot.dat$scenario.type=="A",]$theoretical.power),4),
+xlab='N',ylab='Theoretical power',
+xaxt='n',yaxt='n',type='l',xlim=c(50,300),ylim=c(0,1),lty=4,col='#03a18a')
+axis(1,c(50,100,200,300))
+axis(2,seq(0,1,0.1))
+lines(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$theoretical.power),col='#c6d18d',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$theoretical.power),col='#c0c23b',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$theoretical.power),col='#da77c7',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$theoretical.power),col='#a12471',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$theoretical.power),col='#b5a180',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$theoretical.power),col='#9b6541',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$theoretical.power),col='#30a466',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$theoretical.power),col='#1a342b',lty=4)
+points(unique(plot.dat[plot.dat$scenario %in% c("1A",'1A'),]$N),
+rep(unique(plot.dat[plot.dat$scenario %in% c("1A",'1A'),]$theoretical.power),4),col='#03a18a',pch=17)
+points(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$theoretical.power),col='#c6d18d',pch=17)
+points(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$theoretical.power),col='#c0c23b',pch=17)
+points(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$theoretical.power),col='#da77c7',pch=17)
+points(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$theoretical.power),col='#a12471',pch=17)
+points(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$theoretical.power),col='#b5a180',pch=17)
+points(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$theoretical.power),col='#9b6541',pch=17)
+points(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$theoretical.power),col='#30a466',pch=17)
+points(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$theoretical.power),col='#1a342b',pch=17)
+legend("topleft",lty=rep(4,5),col=c("#03a18a","#c0c23b","#a12471","#9b6541","#1a342b"),
+pch=rep(17,5),
+legend=c("Scenario A",
+"Scenario 1-2 / B-D",
+"Scenario 3-4 / B-D",
+"Scenario 1-2 / C-E",
+"Scenario 3-4 / C-E"),cex=0.7)
+# real
+plot.dat <- res.dat[res.dat$scenario %in% sapply(c('A','B','C','D','E'),function(x) paste0(1:4,x)),]
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario.type=="A" & plot.dat$N==x,]$h0.rejected.p))
+plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
+xlab='N',ylab='Proportion of null hypothesis rejection',
+xaxt='n',yaxt='n',type='l',xlim=c(50,300),ylim=c(0,1),lty=4,col='#03a18a')
+axis(1,c(50,100,200,300))
+axis(2,seq(0,1,0.1))
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1B","1D") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
+mean.A,col='#c6d18d',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2B","2D") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
+mean.A,col='#c0c23b',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3B","3D") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
+mean.A,col='#da77c7',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4B","4D") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
+mean.A,col='#a12471',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1C","1E") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
+mean.A,col='#b5a180',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2C","2E") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
+mean.A,col='#9b6541',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3C","3E") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
+mean.A,col='#30a466',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4C","4E") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
+mean.A,col='#1a342b',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario.type=="A" & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,pch=17,col='#03a18a')
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1B","1D") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
+mean.A,col='#c6d18d',pch=17)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2B","2D") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
+mean.A,col='#c0c23b',pch=17)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3B","3D") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
+mean.A,col='#da77c7',pch=17)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4B","4D") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
+mean.A,col='#a12471',pch=17)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1C","1E") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
+mean.A,col='#b5a180',pch=17)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2C","2E") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
+mean.A,col='#9b6541',pch=17)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3C","3E") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
+mean.A,col='#30a466',pch=17)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4C","4E") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
+mean.A,col='#1a342b',pch=17)
+legend("topleft",lty=rep(4,5),col=c("#03a18a","#c0c23b","#a12471","#9b6541","#1a342b"),
+pch=rep(17,5),
+legend=c("Scenario A",
+"Scenario 1-2 / B-D",
+"Scenario 3-4 / B-D",
+"Scenario 1-2 / C-E",
+"Scenario 3-4 / C-E"),cex=0.7)
+plot.dat <- res.dat[res.dat$scenario %in% sapply(c('A','B','C','D','E'),function(x) paste0(1:4,x)),]
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario.type=="A" & plot.dat$N==x,]$h0.rejected.p))
+plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
+xlab='N',ylab='Proportion of null hypothesis rejection',
+xaxt='n',yaxt='n',type='l',xlim=c(50,300),ylim=c(0,1),lty=4,col='#03a18a')
+axis(1,c(50,100,200,300))
+axis(2,seq(0,1,0.1))
+unique(plot.dat[plot.dat$scenario.type=="A",]$N)
+mean.A
+mean.A[c(4,1:3)]
+# real
+plot.dat <- res.dat[res.dat$scenario %in% sapply(c('A','B','C','D','E'),function(x) paste0(1:4,x)),]
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario.type=="A" & plot.dat$N==x,]$h0.rejected.p))
+plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A[c(4,1:3)],
+xlab='N',ylab='Proportion of null hypothesis rejection',
+xaxt='n',yaxt='n',type='l',xlim=c(50,300),ylim=c(0,1),lty=4,col='#03a18a')
+axis(1,c(50,100,200,300))
+axis(2,seq(0,1,0.1))
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1B","1D") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
+mean.A,col='#c6d18d',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2B","2D") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
+mean.A,col='#c0c23b',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3B","3D") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
+mean.A,col='#da77c7',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4B","4D") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
+mean.A,col='#a12471',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1C","1E") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
+mean.A,col='#b5a180',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2C","2E") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
+mean.A,col='#9b6541',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3C","3E") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
+mean.A,col='#30a466',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4C","4E") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
+mean.A,col='#1a342b',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario.type=="A" & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A[c(4,1:3)],pch=17,col='#03a18a')
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1B","1D") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
+mean.A,col='#c6d18d',pch=17)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2B","2D") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
+mean.A,col='#c0c23b',pch=17)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3B","3D") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
+mean.A,col='#da77c7',pch=17)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4B","4D") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
+mean.A,col='#a12471',pch=17)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1C","1E") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
+mean.A,col='#b5a180',pch=17)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2C","2E") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
+mean.A,col='#9b6541',pch=17)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3C","3E") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
+mean.A,col='#30a466',pch=17)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4C","4E") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
+mean.A,col='#1a342b',pch=17)
+legend("topleft",lty=rep(4,5),col=c("#03a18a","#c0c23b","#a12471","#9b6541","#1a342b"),
+pch=rep(17,5),
+legend=c("Scenario A",
+"Scenario 1-2 / B-D",
+"Scenario 3-4 / B-D",
+"Scenario 1-2 / C-E",
+"Scenario 3-4 / C-E"),cex=0.7)
+# real
+plot.dat <- res.dat[res.dat$scenario %in% sapply(c('A','B','C','D','E'),function(x) paste0(1:4,x)),]
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario.type=="A" & plot.dat$N==x,]$h0.rejected.p))
+plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A[c(4,1:3)],
+xlab='N',ylab='Proportion of null hypothesis rejection',
+xaxt='n',yaxt='n',type='l',xlim=c(50,300),ylim=c(0,1),lty=4,col='#03a18a')
+axis(1,c(50,100,200,300))
+axis(2,seq(0,1,0.1))
+unique(plot.dat[plot.dat$scenario.type=="A",]$N)
+mean.A[c(4,1:3)]
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1B","1D") & plot.dat$N==x,]$h0.rejected.p))
+unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N)
+##############################################################################
+#----------------------------------------------------------------------------#
+################################# POWER PLOTS ################################
+#----------------------------------------------------------------------------#
+##############################################################################
+## Plot 1 - baseline scenarios vs theoretical power
+par(mfrow=c(1,2))
+# theoretical
+plot.dat <- res.dat[res.dat$scenario %in% sapply(c('A','B','C','D','E'),function(x) paste0(1:4,x)),]
+plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),rep(unique(plot.dat[plot.dat$scenario.type=="A",]$theoretical.power),4),
+xlab='N',ylab='Theoretical power',
+xaxt='n',yaxt='n',type='l',xlim=c(50,300),ylim=c(0,1),lty=4,col='#03a18a')
+axis(1,c(50,100,200,300))
+axis(2,seq(0,1,0.1))
+lines(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$theoretical.power),col='#c6d18d',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$theoretical.power),col='#c0c23b',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$theoretical.power),col='#da77c7',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$theoretical.power),col='#a12471',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$theoretical.power),col='#b5a180',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$theoretical.power),col='#9b6541',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$theoretical.power),col='#30a466',lty=4)
+lines(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$theoretical.power),col='#1a342b',lty=4)
+points(unique(plot.dat[plot.dat$scenario %in% c("1A",'1A'),]$N),
+rep(unique(plot.dat[plot.dat$scenario %in% c("1A",'1A'),]$theoretical.power),4),col='#03a18a',pch=17)
+points(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$theoretical.power),col='#c6d18d',pch=17)
+points(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$theoretical.power),col='#c0c23b',pch=17)
+points(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$theoretical.power),col='#da77c7',pch=17)
+points(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$theoretical.power),col='#a12471',pch=17)
+points(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$theoretical.power),col='#b5a180',pch=17)
+points(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$theoretical.power),col='#9b6541',pch=17)
+points(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$theoretical.power),col='#30a466',pch=17)
+points(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
+unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$theoretical.power),col='#1a342b',pch=17)
+legend("topleft",lty=rep(4,5),col=c("#03a18a","#c0c23b","#a12471","#9b6541","#1a342b"),
+pch=rep(17,5),
+legend=c("Scenario A",
+"Scenario 1-2 / B-D",
+"Scenario 3-4 / B-D",
+"Scenario 1-2 / C-E",
+"Scenario 3-4 / C-E"),cex=0.7)
+# real
+plot.dat <- res.dat[res.dat$scenario %in% sapply(c('A','B','C','D','E'),function(x) paste0(1:4,x)),]
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario.type=="A" & plot.dat$N==x,]$h0.rejected.p))
+plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A[c(4,1:3)],
+xlab='N',ylab='Proportion of null hypothesis rejection',
+xaxt='n',yaxt='n',type='l',xlim=c(50,300),ylim=c(0,1),lty=4,col='#03a18a')
+axis(1,c(50,100,200,300))
+axis(2,seq(0,1,0.1))
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1B","1D") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
+mean.A[c(4,1:3)],col='#c6d18d',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2B","2D") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
+mean.A[c(4,1:3)],col='#c0c23b',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3B","3D") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
+mean.A[c(4,1:3)],col='#da77c7',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4B","4D") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
+mean.A[c(4,1:3)],col='#a12471',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1C","1E") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
+mean.A[c(4,1:3)],col='#b5a180',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2C","2E") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
+mean.A[c(4,1:3)],col='#9b6541',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3C","3E") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
+mean.A[c(4,1:3)],col='#30a466',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4C","4E") & plot.dat$N==x,]$h0.rejected.p))
+lines(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
+mean.A[c(4,1:3)],col='#1a342b',lty=4)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario.type=="A" & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A[c(4,1:3)],pch=17,col='#03a18a')
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1B","1D") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
+mean.A[c(4,1:3)],col='#c6d18d',pch=17)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2B","2D") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
+mean.A[c(4,1:3)],col='#c0c23b',pch=17)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3B","3D") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
+mean.A[c(4,1:3)],col='#da77c7',pch=17)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4B","4D") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
+mean.A[c(4,1:3)],col='#a12471',pch=17)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1C","1E") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
+mean.A[c(4,1:3)],col='#b5a180',pch=17)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2C","2E") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
+mean.A[c(4,1:3)],col='#9b6541',pch=17)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3C","3E") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
+mean.A[c(4,1:3)],col='#30a466',pch=17)
+mean.A <- sapply(c(50,100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4C","4E") & plot.dat$N==x,]$h0.rejected.p))
+points(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
+mean.A[c(4,1:3)],col='#1a342b',pch=17)
+legend("topleft",lty=rep(4,5),col=c("#03a18a","#c0c23b","#a12471","#9b6541","#1a342b"),
+pch=rep(17,5),
+legend=c("Scenario A",
+"Scenario 1-2 / B-D",
+"Scenario 3-4 / B-D",
+"Scenario 1-2 / C-E",
+"Scenario 3-4 / C-E"),cex=0.7)
diff --git a/RProject/Scripts/desc_analysis_dif.R b/RProject/Scripts/desc_analysis_dif.R
index c07bb54..e7c8480 100644
--- a/RProject/Scripts/desc_analysis_dif.R
+++ b/RProject/Scripts/desc_analysis_dif.R
@@ -6,13 +6,13 @@
 
 ## Proportion of rejected h0 per dif value in h0 scenarios (A) by DIF size
 
-res.null <- res.dat.dif[res.dat.dif$eff.size==0,]
+res.null <- res.dat[res.dat$eff.size==0,]
 boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2,1),xlab='DIF size',ylab='H0 rejection proportion in target scenario')
-res.null0 <- res.dat.dif[res.dat.dif$eff.size==0 & res.dat.dif$dif.size==0,]
+res.null0 <- res.dat[res.dat$eff.size==0 & res.dat$dif.size==0,]
 points(y=res.null0$h0.rejected.p,x=rep(3,nrow(res.null0)),col='gray',pch=3)
-res.null3 <- res.dat.dif[res.dat.dif$eff.size==0 & res.dat.dif$dif.size==-0.3,]
+res.null3 <- res.dat[res.dat$eff.size==0 & res.dat$dif.size==-0.3,]
 points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
-res.null5 <- res.dat.dif[res.dat.dif$eff.size==0 & res.dat.dif$dif.size==-0.5,]
+res.null5 <- res.dat[res.dat$eff.size==0 & res.dat$dif.size==-0.5,]
 points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
 
 ## Proportion of rejected h0 per dif value in h0 scenarios (A) by DIF size
@@ -20,7 +20,7 @@ points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
 par(mfrow=c(2,2))
 
 # 1 item
-res.null <- res.dat.dif[res.dat.dif$eff.size==0 & res.dat.dif$nb.dif==1,]
+res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==1,]
 boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
         ylab='H0 rejection proportion in target scenario',main='DIF on 1 item',ylim=c(0,1))
 res.null3 <- res.null[res.null$dif.size==-0.3,]
@@ -29,7 +29,7 @@ res.null5 <- res.null[res.null$dif.size==-0.5,]
 points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
 
 # 2 items
-res.null <- res.dat.dif[res.dat.dif$eff.size==0 & res.dat.dif$nb.dif==2,]
+res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==2,]
 boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
         ylab='H0 rejection proportion in target scenario',main='DIF on 2 items',ylim=c(0,1))
 res.null3 <- res.null[res.null$dif.size==-0.3,]
@@ -38,7 +38,7 @@ res.null5 <- res.null[res.null$dif.size==-0.5,]
 points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
 
 # 3 items
-res.null <- res.dat.dif[res.dat.dif$eff.size==0 & res.dat.dif$nb.dif==3,]
+res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==3,]
 boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
         ylab='H0 rejection proportion in target scenario',main='DIF on 3 items',ylim=c(0,1))
 res.null3 <- res.null[res.null$dif.size==-0.3,]
@@ -53,7 +53,7 @@ par(mfrow=c(1,1))
 par(mfrow=c(2,2))
 
 # 1 item
-res.null <- res.dat.dif[res.dat.dif$eff.size==0 & res.dat.dif$nb.dif==1 & res.dat.dif$J==4,]
+res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==1 & res.dat$J==4,]
 boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
         ylab='H0 rejection proportion in target scenario',main='DIF on 1 item',ylim=c(0,1))
 res.null3 <- res.null[res.null$dif.size==-0.3,]
@@ -62,7 +62,7 @@ res.null5 <- res.null[res.null$dif.size==-0.5,]
 points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
 
 # 2 items
-res.null <- res.dat.dif[res.dat.dif$eff.size==0 & res.dat.dif$nb.dif==2 & res.dat.dif$J==4,]
+res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==2 & res.dat$J==4,]
 boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
         ylab='H0 rejection proportion in target scenario',main='DIF on 2 items',ylim=c(0,1))
 res.null3 <- res.null[res.null$dif.size==-0.3,]
@@ -78,7 +78,7 @@ par(mfrow=c(1,1))
 par(mfrow=c(2,2))
 
 # 2 items
-res.null <- res.dat.dif[res.dat.dif$eff.size==0 & res.dat.dif$nb.dif==2 & res.dat.dif$J==7,]
+res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==2 & res.dat$J==7,]
 boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
         ylab='H0 rejection proportion in target scenario',main='DIF on 2 items',ylim=c(0,1))
 res.null3 <- res.null[res.null$dif.size==-0.3,]
@@ -87,7 +87,7 @@ res.null5 <- res.null[res.null$dif.size==-0.5,]
 points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
 
 # 3 items
-res.null <- res.dat.dif[res.dat.dif$eff.size==0 & res.dat.dif$nb.dif==3 & res.dat.dif$J==7,]
+res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==3 & res.dat$J==7,]
 boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
         ylab='H0 rejection proportion in target scenario',main='DIF on 3 items',ylim=c(0,1))
 res.null3 <- res.null[res.null$dif.size==-0.3,]
@@ -105,7 +105,7 @@ par(mfrow=c(1,1))
 par(mfrow=c(2,2))
 
 # 1 item
-res.null <- res.dat.dif[res.dat.dif$eff.size==0 & res.dat.dif$nb.dif==1 & res.dat.dif$N==100,]
+res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==1 & res.dat$N==100,]
 boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
         ylab='H0 rejection proportion in target scenario',main='DIF on 1 item',ylim=c(0,1))
 res.null3 <- res.null[res.null$dif.size==-0.3,]
@@ -114,7 +114,7 @@ res.null5 <- res.null[res.null$dif.size==-0.5,]
 points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
 
 # 2 items
-res.null <- res.dat.dif[res.dat.dif$eff.size==0 & res.dat.dif$nb.dif==2 & res.dat.dif$N==100,]
+res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==2 & res.dat$N==100,]
 boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
         ylab='H0 rejection proportion in target scenario',main='DIF on 2 items',ylim=c(0,1))
 res.null3 <- res.null[res.null$dif.size==-0.3,]
@@ -123,7 +123,7 @@ res.null5 <- res.null[res.null$dif.size==-0.5,]
 points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
 
 # 3 items
-res.null <- res.dat.dif[res.dat.dif$eff.size==0 & res.dat.dif$nb.dif==3 & res.dat.dif$N==100,]
+res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==3 & res.dat$N==100,]
 boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
         ylab='H0 rejection proportion in target scenario',main='DIF on 3 items',ylim=c(0,1))
 res.null3 <- res.null[res.null$dif.size==-0.3,]
@@ -140,7 +140,7 @@ par(mfrow=c(1,1))
 par(mfrow=c(2,2))
 
 # 1 item
-res.null <- res.dat.dif[res.dat.dif$eff.size==0 & res.dat.dif$nb.dif==1 & res.dat.dif$N==200,]
+res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==1 & res.dat$N==200,]
 boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
         ylab='H0 rejection proportion in target scenario',main='DIF on 1 item',ylim=c(0,1))
 res.null3 <- res.null[res.null$dif.size==-0.3,]
@@ -149,7 +149,7 @@ res.null5 <- res.null[res.null$dif.size==-0.5,]
 points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
 
 # 2 items
-res.null <- res.dat.dif[res.dat.dif$eff.size==0 & res.dat.dif$nb.dif==2 & res.dat.dif$N==200,]
+res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==2 & res.dat$N==200,]
 boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
         ylab='H0 rejection proportion in target scenario',main='DIF on 2 items',ylim=c(0,1))
 res.null3 <- res.null[res.null$dif.size==-0.3,]
@@ -158,7 +158,7 @@ res.null5 <- res.null[res.null$dif.size==-0.5,]
 points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
 
 # 3 items
-res.null <- res.dat.dif[res.dat.dif$eff.size==0 & res.dat.dif$nb.dif==3 & res.dat.dif$N==200,]
+res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==3 & res.dat$N==200,]
 boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
         ylab='H0 rejection proportion in target scenario',main='DIF on 3 items',ylim=c(0,1))
 res.null3 <- res.null[res.null$dif.size==-0.3,]
@@ -175,7 +175,7 @@ par(mfrow=c(1,1))
 par(mfrow=c(2,2))
 
 # 1 item
-res.null <- res.dat.dif[res.dat.dif$eff.size==0 & res.dat.dif$nb.dif==1 & res.dat.dif$N==300,]
+res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==1 & res.dat$N==300,]
 boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
         ylab='H0 rejection proportion in target scenario',main='DIF on 1 item',ylim=c(0,1))
 res.null3 <- res.null[res.null$dif.size==-0.3,]
@@ -184,7 +184,7 @@ res.null5 <- res.null[res.null$dif.size==-0.5,]
 points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
 
 # 2 items
-res.null <- res.dat.dif[res.dat.dif$eff.size==0 & res.dat.dif$nb.dif==2 & res.dat.dif$N==300,]
+res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==2 & res.dat$N==300,]
 boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
         ylab='H0 rejection proportion in target scenario',main='DIF on 2 items',ylim=c(0,1))
 res.null3 <- res.null[res.null$dif.size==-0.3,]
@@ -193,7 +193,7 @@ res.null5 <- res.null[res.null$dif.size==-0.5,]
 points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
 
 # 3 items
-res.null <- res.dat.dif[res.dat.dif$eff.size==0 & res.dat.dif$nb.dif==3 & res.dat.dif$N==300,]
+res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==3 & res.dat$N==300,]
 boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
         ylab='H0 rejection proportion in target scenario',main='DIF on 3 items',ylim=c(0,1))
 res.null3 <- res.null[res.null$dif.size==-0.3,]
@@ -213,17 +213,17 @@ par(mfrow=c(1,1))
 
 ## Proportion of rejected h0 per dif value in h1 scenarios by DIF size // eff.size positive
 
-res.null <- res.dat.dif[res.dat.dif$eff.size>0,]
+res.null <- res.dat[res.dat$eff.size>0,]
 boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2,4,2,3),xlab='DIF size',ylab='H0 rejection proportion in target scenario')
-res.null0 <- res.dat.dif[res.dat.dif$eff.size>0 & res.dat.dif$dif.size==0,]
+res.null0 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0,]
 points(y=res.null0$h0.rejected.p,x=rep(3,nrow(res.null0)),col='darkblue',pch=3)
-res.null3 <- res.dat.dif[res.dat.dif$eff.size>0 & res.dat.dif$dif.size==-0.3,]
+res.null3 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==-0.3,]
 points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
-res.null5 <- res.dat.dif[res.dat.dif$eff.size>0 & res.dat.dif$dif.size==-0.5,]
+res.null5 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==-0.5,]
 points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
-res.null3 <- res.dat.dif[res.dat.dif$eff.size>0 & res.dat.dif$dif.size==0.3,]
+res.null3 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0.3,]
 points(y=res.null3$h0.rejected.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
-res.null5 <- res.dat.dif[res.dat.dif$eff.size>0 & res.dat.dif$dif.size==0.5,]
+res.null5 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0.5,]
 points(y=res.null5$h0.rejected.p,x=rep(5,nrow(res.null5)),col='#053305',pch=3)
 
 
@@ -231,32 +231,32 @@ points(y=res.null5$h0.rejected.p,x=rep(5,nrow(res.null5)),col='#053305',pch=3)
 
 ####### N=100
 
-res.null <- res.dat.dif[res.dat.dif$eff.size>0 & res.dat.dif$N==100,]
+res.null <- res.dat[res.dat$eff.size>0 & res.dat$N==100,]
 boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2,4,2,3),xlab='DIF size',ylab='H0 rejection proportion in target scenario')
-res.null0 <- res.dat.dif[res.dat.dif$eff.size>0 & res.dat.dif$dif.size==0 & res.dat.dif$N==100,]
+res.null0 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0 & res.dat$N==100,]
 points(y=res.null0$h0.rejected.p,x=rep(3,nrow(res.null0)),col='darkblue',pch=3)
-res.null3 <- res.dat.dif[res.dat.dif$eff.size>0 & res.dat.dif$dif.size==-0.3 & res.dat.dif$N==100,]
+res.null3 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==-0.3 & res.dat$N==100,]
 points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
-res.null5 <- res.dat.dif[res.dat.dif$eff.size>0 & res.dat.dif$dif.size==-0.5 & res.dat.dif$N==100,]
+res.null5 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==-0.5 & res.dat$N==100,]
 points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
-res.null3 <- res.dat.dif[res.dat.dif$eff.size>0 & res.dat.dif$dif.size==0.3 & res.dat.dif$N==100,]
+res.null3 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0.3 & res.dat$N==100,]
 points(y=res.null3$h0.rejected.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
-res.null5 <- res.dat.dif[res.dat.dif$eff.size>0 & res.dat.dif$dif.size==0.5 & res.dat.dif$N==100,]
+res.null5 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0.5 & res.dat$N==100,]
 points(y=res.null5$h0.rejected.p,x=rep(5,nrow(res.null5)),col='#053305',pch=3)
 
 ####### N=300 // DIF à 0.5 - QUELS SONT LES SCENARIOS EN HAUT
 
-res.null <- res.dat.dif[res.dat.dif$eff.size>0 & res.dat.dif$N==300,]
+res.null <- res.dat[res.dat$eff.size>0 & res.dat$N==300,]
 boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2,4,2,3),xlab='DIF size',ylab='H0 rejection proportion in target scenario')
-res.null0 <- res.dat.dif[res.dat.dif$eff.size>0 & res.dat.dif$dif.size==0 & res.dat.dif$N==300,]
+res.null0 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0 & res.dat$N==300,]
 points(y=res.null0$h0.rejected.p,x=rep(3,nrow(res.null0)),col='darkblue',pch=3)
-res.null3 <- res.dat.dif[res.dat.dif$eff.size>0 & res.dat.dif$dif.size==-0.3 & res.dat.dif$N==300,]
+res.null3 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==-0.3 & res.dat$N==300,]
 points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
-res.null5 <- res.dat.dif[res.dat.dif$eff.size>0 & res.dat.dif$dif.size==-0.5 & res.dat.dif$N==300,]
+res.null5 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==-0.5 & res.dat$N==300,]
 points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
-res.null3 <- res.dat.dif[res.dat.dif$eff.size>0 & res.dat.dif$dif.size==0.3 & res.dat.dif$N==300,]
+res.null3 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0.3 & res.dat$N==300,]
 points(y=res.null3$h0.rejected.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
-res.null5 <- res.dat.dif[res.dat.dif$eff.size>0 & res.dat.dif$dif.size==0.5 & res.dat.dif$N==300,]
+res.null5 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0.5 & res.dat$N==300,]
 points(y=res.null5$h0.rejected.p,x=rep(5,nrow(res.null5)),col='#053305',pch=3)
 
 
@@ -264,25 +264,25 @@ points(y=res.null5$h0.rejected.p,x=rep(5,nrow(res.null5)),col='#053305',pch=3)
 
 ####### N=100
 
-res.null <- res.dat.dif[res.dat.dif$eff.size<0 & res.dat.dif$N==100,]
+res.null <- res.dat[res.dat$eff.size<0 & res.dat$N==100,]
 boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2,4,2,3),xlab='DIF size',ylab='H0 rejection proportion in target scenario')
-res.null0 <- res.dat.dif[res.dat.dif$eff.size<0 & res.dat.dif$dif.size==0 & res.dat.dif$N==100,]
+res.null0 <- res.dat[res.dat$eff.size<0 & res.dat$dif.size==0 & res.dat$N==100,]
 points(y=res.null0$h0.rejected.p,x=rep(3,nrow(res.null0)),col='darkblue',pch=3)
-res.null3 <- res.dat.dif[res.dat.dif$eff.size<0 & res.dat.dif$dif.size==-0.3 & res.dat.dif$N==100,]
+res.null3 <- res.dat[res.dat$eff.size<0 & res.dat$dif.size==-0.3 & res.dat$N==100,]
 points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
-res.null5 <- res.dat.dif[res.dat.dif$eff.size<0 & res.dat.dif$dif.size==-0.5 & res.dat.dif$N==100,]
+res.null5 <- res.dat[res.dat$eff.size<0 & res.dat$dif.size==-0.5 & res.dat$N==100,]
 points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
 
 
 ####### N=300 // DIF à 0.5 - QUELS SONT LES SCENARIOS EN HAUT
 
-res.null <- res.dat.dif[res.dat.dif$eff.size<0 & res.dat.dif$N==300,]
+res.null <- res.dat[res.dat$eff.size<0 & res.dat$N==300,]
 boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2,4,2,3),xlab='DIF size',ylab='H0 rejection proportion in target scenario')
-res.null0 <- res.dat.dif[res.dat.dif$eff.size<0 & res.dat.dif$dif.size==0 & res.dat.dif$N==300,]
+res.null0 <- res.dat[res.dat$eff.size<0 & res.dat$dif.size==0 & res.dat$N==300,]
 points(y=res.null0$h0.rejected.p,x=rep(3,nrow(res.null0)),col='darkblue',pch=3)
-res.null3 <- res.dat.dif[res.dat.dif$eff.size<0 & res.dat.dif$dif.size==-0.3 & res.dat.dif$N==300,]
+res.null3 <- res.dat[res.dat$eff.size<0 & res.dat$dif.size==-0.3 & res.dat$N==300,]
 points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
-res.null5 <- res.dat.dif[res.dat.dif$eff.size<0 & res.dat.dif$dif.size==-0.5 & res.dat.dif$N==300,]
+res.null5 <- res.dat[res.dat$eff.size<0 & res.dat$dif.size==-0.5 & res.dat$N==300,]
 points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
 
 
@@ -294,110 +294,110 @@ points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
 
 # Overall
 
-boxplot(true.value.in.ci.p~dif.size,data=res.dat.dif,col=c(2,3),xlab='DIF size',
+boxplot(true.value.in.ci.p~dif.size,data=res.dat,col=c(2,3),xlab='DIF size',
         ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
-res.null3 <- res.dat.dif[res.dat.dif$dif.size==-0.5,]
+res.null3 <- res.dat[res.dat$dif.size==-0.5,]
 points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
-res.null3 <- res.dat.dif[res.dat.dif$dif.size==0.5,]
+res.null3 <- res.dat[res.dat$dif.size==0.5,]
 points(y=res.null3$true.value.in.ci.p,x=rep(5,nrow(res.null3)),col='#590b0c',pch=3)
-res.null3 <- res.dat.dif[res.dat.dif$dif.size==-0.3,]
+res.null3 <- res.dat[res.dat$dif.size==-0.3,]
 points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
-res.null3 <- res.dat.dif[res.dat.dif$dif.size==0.3,]
+res.null3 <- res.dat[res.dat$dif.size==0.3,]
 points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#053305',pch=3)
-res.null3 <- res.dat.dif[res.dat.dif$dif.size==0,]
+res.null3 <- res.dat[res.dat$dif.size==0,]
 points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='gray',pch=3)
 
 # J=4
-res.dat.dif.temp <- res.dat.dif[res.dat.dif$J==4,]
+res.dat.temp <- res.dat[res.dat$J==4,]
 
-boxplot(true.value.in.ci.p~dif.size,data=res.dat.dif.temp,col=c(2,3),xlab='DIF size',
+boxplot(true.value.in.ci.p~dif.size,data=res.dat.temp,col=c(2,3),xlab='DIF size',
         ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==-0.5,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.5,]
 points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==0.5,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.5,]
 points(y=res.null3$true.value.in.ci.p,x=rep(5,nrow(res.null3)),col='#590b0c',pch=3)
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==-0.3,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.3,]
 points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==0.3,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
 points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#053305',pch=3)
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==0,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==0,]
 points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='gray',pch=3)
 
 
 # J=7
-res.dat.dif.temp <- res.dat.dif[res.dat.dif$J==7,]
+res.dat.temp <- res.dat[res.dat$J==7,]
 
-boxplot(true.value.in.ci.p~dif.size,data=res.dat.dif.temp,col=c(2,3),xlab='DIF size',
+boxplot(true.value.in.ci.p~dif.size,data=res.dat.temp,col=c(2,3),xlab='DIF size',
         ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==-0.5,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.5,]
 points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==0.5,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.5,]
 points(y=res.null3$true.value.in.ci.p,x=rep(5,nrow(res.null3)),col='#590b0c',pch=3)
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==-0.3,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.3,]
 points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==0.3,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
 points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#053305',pch=3)
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==0,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==0,]
 points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='gray',pch=3)
 
 # J=4 / 1 DIF
-res.dat.dif.temp <- res.dat.dif[res.dat.dif$J==4 & res.dat.dif$nb.dif==1,]
+res.dat.temp <- res.dat[res.dat$J==4 & res.dat$nb.dif==1,]
 
-boxplot(true.value.in.ci.p~dif.size,data=res.dat.dif.temp,col=c(2,3,3,2),xlab='DIF size',
+boxplot(true.value.in.ci.p~dif.size,data=res.dat.temp,col=c(2,3,3,2),xlab='DIF size',
         ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==-0.5,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.5,]
 points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==0.5,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.5,]
 points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==-0.3,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.3,]
 points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==0.3,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
 points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='#053305',pch=3)
 
 # J=4 / 2 DIF
-res.dat.dif.temp <- res.dat.dif[res.dat.dif$J==4 & res.dat.dif$nb.dif==2,]
+res.dat.temp <- res.dat[res.dat$J==4 & res.dat$nb.dif==2,]
 
-boxplot(true.value.in.ci.p~dif.size,data=res.dat.dif.temp,col=c(2,3,3,2),xlab='DIF size',
+boxplot(true.value.in.ci.p~dif.size,data=res.dat.temp,col=c(2,3,3,2),xlab='DIF size',
         ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==-0.5,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.5,]
 points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==0.5,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.5,]
 points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==-0.3,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.3,]
 points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==0.3,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
 points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='#053305',pch=3)
 
 
 # J=7 / 2 DIF
 
-res.dat.dif.temp <- res.dat.dif[res.dat.dif$J==7 & res.dat.dif$nb.dif==2,]
+res.dat.temp <- res.dat[res.dat$J==7 & res.dat$nb.dif==2,]
 
-boxplot(true.value.in.ci.p~dif.size,data=res.dat.dif.temp,col=c(2,3,3,2),xlab='DIF size',
+boxplot(true.value.in.ci.p~dif.size,data=res.dat.temp,col=c(2,3,3,2),xlab='DIF size',
         ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==-0.5,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.5,]
 points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==0.5,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.5,]
 points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==-0.3,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.3,]
 points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==0.3,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
 points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='#053305',pch=3)
 
 
 # J=7 / 3 DIF
 
-res.dat.dif.temp <- res.dat.dif[res.dat.dif$J==7 & res.dat.dif$nb.dif==3,]
+res.dat.temp <- res.dat[res.dat$J==7 & res.dat$nb.dif==3,]
 
-boxplot(true.value.in.ci.p~dif.size,data=res.dat.dif.temp,col=c(2,3,3,2),xlab='DIF size',
+boxplot(true.value.in.ci.p~dif.size,data=res.dat.temp,col=c(2,3,3,2),xlab='DIF size',
         ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==-0.5,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.5,]
 points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==0.5,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.5,]
 points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==-0.3,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.3,]
 points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
-res.null3 <- res.dat.dif.temp[res.dat.dif.temp$dif.size==0.3,]
+res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
 points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='#053305',pch=3)
 
 
@@ -409,30 +409,30 @@ points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='#053305',pch
 
 # Overall
 
-boxplot(beta.same.sign.truebeta.p~dif.size,data=res.dat.dif,col=c(2,3),xlab='DIF size',
+boxplot(beta.same.sign.truebeta.p~dif.size,data=res.dat,col=c(2,3),xlab='DIF size',
         ylab='Proportion of estimates with same sign as true value in target scenario',main='DIF on 3 items',ylim=c(0,1))
-res.null3 <- res.dat.dif[res.dat.dif$dif.size==-0.5,]
+res.null3 <- res.dat[res.dat$dif.size==-0.5,]
 points(y=res.null3$beta.same.sign.truebeta.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
-res.null3 <- res.dat.dif[res.dat.dif$dif.size==0.5,]
+res.null3 <- res.dat[res.dat$dif.size==0.5,]
 points(y=res.null3$beta.same.sign.truebeta.p,x=rep(5,nrow(res.null3)),col='#590b0c',pch=3)
-res.null3 <- res.dat.dif[res.dat.dif$dif.size==-0.3,]
+res.null3 <- res.dat[res.dat$dif.size==-0.3,]
 points(y=res.null3$beta.same.sign.truebeta.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
-res.null3 <- res.dat.dif[res.dat.dif$dif.size==0.3,]
+res.null3 <- res.dat[res.dat$dif.size==0.3,]
 points(y=res.null3$beta.same.sign.truebeta.p,x=rep(4,nrow(res.null3)),col='#053305',pch=3)
-res.null3 <- res.dat.dif[res.dat.dif$dif.size==0,]
+res.null3 <- res.dat[res.dat$dif.size==0,]
 points(y=res.null3$beta.same.sign.truebeta.p,x=rep(3,nrow(res.null3)),col='gray',pch=3)
 
 # Overall // H0 rejected 
 
-boxplot(beta.same.sign.truebeta.signif.p~dif.size,data=res.dat.dif,col=c(2,3),xlab='DIF size',
+boxplot(beta.same.sign.truebeta.signif.p~dif.size,data=res.dat,col=c(2,3),xlab='DIF size',
         ylab='Proportion of estimates with same sign as true value in target scenario',main='When H0 rejected',ylim=c(0,1))
-res.null3 <- res.dat.dif[res.dat.dif$dif.size==-0.5,]
+res.null3 <- res.dat[res.dat$dif.size==-0.5,]
 points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
-res.null3 <- res.dat.dif[res.dat.dif$dif.size==0.5,]
+res.null3 <- res.dat[res.dat$dif.size==0.5,]
 points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(5,nrow(res.null3)),col='#590b0c',pch=3)
-res.null3 <- res.dat.dif[res.dat.dif$dif.size==-0.3,]
+res.null3 <- res.dat[res.dat$dif.size==-0.3,]
 points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
-res.null3 <- res.dat.dif[res.dat.dif$dif.size==0.3,]
+res.null3 <- res.dat[res.dat$dif.size==0.3,]
 points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(4,nrow(res.null3)),col='#053305',pch=3)
-res.null3 <- res.dat.dif[res.dat.dif$dif.size==0,]
+res.null3 <- res.dat[res.dat$dif.size==0,]
 points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(3,nrow(res.null3)),col='gray',pch=3)
\ No newline at end of file