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