diff --git a/RProject/.RData b/RProject/.RData
index 127fc83..2bf4797 100644
Binary files a/RProject/.RData and b/RProject/.RData differ
diff --git a/RProject/.Rhistory b/RProject/.Rhistory
index 6e2f07c..5f66ad1 100644
--- a/RProject/.Rhistory
+++ b/RProject/.Rhistory
@@ -1,512 +1,512 @@
-if (is.null(dif.items)) {
-tam1 <-  TAM::tam.mml(resp=resp,Y=grp,irtmodel = model,est.variance = T,verbose=F)
+stop("Models RSM and GRSM cannot be used for unequal numbers of response categories!")
+}
+design_list <- design_GPCMlasso(formula = formula, data = data,
+Y = Y, RSM = RSM, GPCM = GPCM, DSF = DSF, all.dummies = control$all.dummies,
+main.effects = main.effects)
+if (design_list$m == 0) {
+control$lambda <- 0
+control$adaptive <- FALSE
+cv <- FALSE
+}
+if (length(control$lambda) == 1) {
+cv <- FALSE
+}
+loglik_fun <- loglikPCMlasso
+score_fun <- scorePCMlasso
+log_score_fun <- loglikscorePCMlasso2
+if (all(k == 2)) {
+loglik_fun <- loglikDIFlasso
+score_fun <- scoreDIFlasso
+log_score_fun <- loglikscoreDIFlasso2
+}
+fit <- fit_GPCMlasso(model = model, loglik_fun = loglik_fun, acoefs2 = acoefs2,
+score_fun = score_fun, log_score_fun, design_list = design_list,
+control = control, start = NULL, scale_fac = 1, main.effects = main.effects)
+if (is.null(control$lambda)) {
+control$lambda <- fit$lambda
+}
+coefficients <- fit$coefficients
+coef.rescal <- fit$coef.rescal
+logLik <- fit$logLik
+df <- fit$df
+BIC <- -2 * logLik + log(design_list$n) * df
+AIC <- -2 * logLik + 2 * df
+cAIC <- -2 * logLik + 2 * df + 2 * df * (df - 1)/(design_list$n -
+df - 1)
+if (cv) {
+cv_error <- fit_cv_GPCMlasso(model, design_list, control,
+score_fun, loglik_fun, log_score_fun, Y)
 }
 else {
-nodif.items <- setdiff(items,dif.items)
-nodif.items.list <- paste0("item",nodif.items)
-dif.items.list <- c(paste0("item",dif.items,"_noTT"),paste0("item",dif.items,"_TT"))
-if (!is.null(dif.items)) {
-if (length(nodif.items)!=0) {
-resp <- df[,c(nodif.items.list,dif.items.list)]
+cv_error <- NULL
+}
+ret.list <- list(coefficients = coefficients, logLik = logLik,
+cv_error = cv_error, call = match.call(), model = model,
+data = data, control = control, DSF = DSF, formula = formula,
+item.names = item.names, Y = Y, design_list = design_list,
+AIC = AIC, BIC = BIC, cAIC = cAIC, df = df, coef.rescal = coef.rescal,
+main.effects = main.effects)
+class(ret.list) <- "GPCMlasso"
+return(ret.list)
+}
+GPCMlasso2 <- function (formula, data, DSF = FALSE, model = c("PCM", "RSM",
+"GPCM", "GRSM", "RM", "2PL"), control = ctrl_GPCMlasso(),
+cv = FALSE, main.effects = TRUE,acoefs2=NULL)
+{
+if (!is.data.frame(data))
+stop("data has to be a data.frame!")
+Y <- model.response(model.frame(formula, data = data))
+for (i in 1:ncol(Y)) {
+Y[, i] <- as.factor(Y[, i])
+}
+item.names <- colnames(Y)
+model <- match.arg(model)
+RSM <- GPCM <- FALSE
+if (model %in% c("GRSM", "RSM")) {
+RSM <- TRUE
+}
+if (model %in% c("GRSM", "GPCM", "2PL")) {
+GPCM <- TRUE
+}
+k <- apply(Y, 2, function(x) {
+length(levels(as.factor(x)))
+})
+if (all(k == 2) | model %in% c("RM", "2PL")) {
+DSF <- FALSE
+control$cv.crit <- "deviance"
+}
+if (length(unique(k)) > 1 & RSM) {
+stop("Models RSM and GRSM cannot be used for unequal numbers of response categories!")
+}
+design_list <- design_GPCMlasso(formula = formula, data = data,
+Y = Y, RSM = RSM, GPCM = GPCM, DSF = DSF, all.dummies = control$all.dummies,
+main.effects = main.effects)
+if (design_list$m == 0) {
+control$lambda <- 0
+control$adaptive <- FALSE
+cv <- FALSE
+}
+if (length(control$lambda) == 1) {
+cv <- FALSE
+}
+loglik_fun <- loglikPCMlasso
+score_fun <- scorePCMlasso
+log_score_fun <- loglikscorePCMlasso2
+if (all(k == 2)) {
+loglik_fun <- loglikDIFlasso
+score_fun <- scoreDIFlasso
+log_score_fun <- loglikscoreDIFlasso2
+}
+fit <- fit_GPCMlasso(model = model, loglik_fun = loglik_fun, acoefs2 = acoefs2,
+score_fun = score_fun, log_score_fun, design_list = design_list,
+control = control, start = NULL, scale_fac = 1, main.effects = main.effects)
+if (is.null(control$lambda)) {
+control$lambda <- fit$lambda
+}
+coefficients <- fit$coefficients
+coef.rescal <- fit$coef.rescal
+logLik <- fit$logLik
+df <- fit$df
+BIC <- -2 * logLik + log(design_list$n) * df
+AIC <- -2 * logLik + 2 * df
+cAIC <- -2 * logLik + 2 * df + 2 * df * (df - 1)/(design_list$n -
+df - 1)
+if (cv) {
+cv_error <- fit_cv_GPCMlasso(model, design_list, control,
+score_fun, loglik_fun, log_score_fun, Y)
 }
 else {
-resp <- df[,c(dif.items.list)]
-}
-}
-else {
-resp <- df[,c(nodif.items.list)]
-}
-tam1 <- TAM::tam.mml(resp=resp,Y=grp,irtmodel = model,est.variance = T,verbose=F)
-}
-}
-return(summary(tam1))
-}
-rosali <- function(dat=NULL,items=NULL,group=NULL) {
-nbitems <- length(items)
-items2 <- items
-# Un seul groupe
-if (length(group)!=1) {
-stop('Only one variable can be used for the group option')
-}
-# Recoder groupe en 0/1
-dat[,group] <- as.factor(dat[,group])
-if (!(all(names(levels(dat[,group]))==c("0","1")))) {
-levels(dat[,group]) <- c("0","1")
-}
-dat[,group] <- as.numeric(dat[,group])-1
-model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
-model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
-return(c(model_a))
-}
-rosali(dat=dat,items=1:7,group="TT")
-rosali <- function(dat=NULL,items=NULL,group=NULL) {
-nbitems <- length(items)
-items2 <- items
-# Un seul groupe
-if (length(group)!=1) {
-stop('Only one variable can be used for the group option')
-}
-# Recoder groupe en 0/1
-dat[,group] <- as.factor(dat[,group])
-if (!(all(names(levels(dat[,group]))==c("0","1")))) {
-levels(dat[,group]) <- c("0","1")
-}
-dat[,group] <- as.numeric(dat[,group])-1
-model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
-model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
-return(c(model_b))
-}
-rosali(dat=dat,items=1:7,group="TT")
-rosali <- function(dat=NULL,items=NULL,group=NULL) {
-nbitems <- length(items)
-items2 <- items
-# Un seul groupe
-if (length(group)!=1) {
-stop('Only one variable can be used for the group option')
-}
-# Recoder groupe en 0/1
-dat[,group] <- as.factor(dat[,group])
-if (!(all(names(levels(dat[,group]))==c("0","1")))) {
-levels(dat[,group]) <- c("0","1")
-}
-dat[,group] <- as.numeric(dat[,group])-1
-model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
-model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
-}
-rosali(dat=dat,items=1:7,group="TT")
-model_b <- pcm(df=dat,items=1:7,group = "TT",dif.items = 1:7)
-model_a <- pcm(df=dat,items=1:7,group = "TT",dif.items = NULL)
-lrt
-lrt''
-??lrt
-loglik(model_a)
-logLik(model_a)
-## File Name: pcm.R
-## File version: 1.0
-#' @import TAM
-#' @export
-pcm <- function(df=NULL,items=NULL,group=NULL,model="PCM2",method="MML",dif.items=NULL) {
-# Prepare analysis
-if (is.null(items)) {
-nbitems <- sum(sapply(1:100,function(x) paste0('item',x)) %in% colnames(df))
-items <- paste0('item',seq(1,nbitems))
-resp <- df[,items]
-}
-else {
-nbitems <- length(items)
-resp <- df[,paste0("item",items)]
-}
-if (is.null(group)) {
-grp <- NULL
-}
-else {
-grp <- df[,group]
-df$grp <- grp
-if (!is.null(dif.items)) {
-for (i in dif.items) {
-df[,paste0('item',i,'_noTT')] <- NA
-df[,paste0('item',i,'_TT')] <- NA
-df[df$grp==0,paste0('item',i,'_noTT')] <- df[df$grp==0,paste0('item',i)]
-df[df$grp==1,paste0('item',i,'_TT')] <- df[df$grp==1,paste0('item',i)]
-}
-}
-}
-# Analyze
-if (method=='MML') {
-if (is.null(dif.items)) {
-tam1 <-  TAM::tam.mml(resp=resp,Y=grp,irtmodel = model,est.variance = T,verbose=F)
-}
-else {
-nodif.items <- setdiff(items,dif.items)
-nodif.items.list <- paste0("item",nodif.items)
-dif.items.list <- c(paste0("item",dif.items,"_noTT"),paste0("item",dif.items,"_TT"))
-if (!is.null(dif.items)) {
-if (length(nodif.items)!=0) {
-resp <- df[,c(nodif.items.list,dif.items.list)]
-}
-else {
-resp <- df[,c(dif.items.list)]
-}
-}
-else {
-resp <- df[,c(nodif.items.list)]
-}
-tam1 <- TAM::tam.mml(resp=resp,Y=grp,irtmodel = model,est.variance = T,verbose=F)
-}
-}
-return(tam1)
-}
-rosali <- function(dat=NULL,items=NULL,group=NULL) {
-nbitems <- length(items)
-items2 <- items
-# Un seul groupe
-if (length(group)!=1) {
-stop('Only one variable can be used for the group option')
-}
-# Recoder groupe en 0/1
-dat[,group] <- as.factor(dat[,group])
-if (!(all(names(levels(dat[,group]))==c("0","1")))) {
-levels(dat[,group]) <- c("0","1")
-}
-dat[,group] <- as.numeric(dat[,group])-1
-model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
-model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
-}
-model_b <- pcm(df=dat,items=1:7,group = "TT",dif.items = 1:7)
-model_a <- pcm(df=dat,items=1:7,group = "TT",dif.items = NULL)
-logLik(model_a)
-logLik(model_b)
-library(lmtest)
-lrtest(model_a,model_b)
-lrtest(model_a,model_b)$p
-ss <- lrtest(model_a,model_b)
-ss$`Pr(>Chisq)`
-ss$`Pr(>Chisq)`[2]
-a <- c()
-a
-c(a,1)
-which(1:4,3)
-?which
-which.min(1:4)
-which.min(2:4)
-rosali <- function(dat=NULL,items=NULL,group=NULL) {
-nbitems <- length(items)
-items2 <- items
-# Un seul groupe
-if (length(group)!=1) {
-stop('Only one variable can be used for the group option')
-}
-# Recoder groupe en 0/1
-dat[,group] <- as.factor(dat[,group])
-if (!(all(names(levels(dat[,group]))==c("0","1")))) {
-levels(dat[,group]) <- c("0","1")
-}
-dat[,group] <- as.numeric(dat[,group])-1
-model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
-model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
-ltest <- lmtest::lrtest(model_a,model_b)
-if (ltest$`Pr(>Chisq)`[2]<1) {
-model_c <- model_b
-difit <- c()
-a <- sapply(1:nbitems,function(x)
-lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
-while (min(a)<0.4) {
-difit <- c(difit,which.min(a))
-model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
-a <- sapply(1:nbitems,function(x)
-lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
-}
-}
-}
-rosali(dat=dat,items=1:7,group="TT")
-rosali <- function(dat=NULL,items=NULL,group=NULL) {
-nbitems <- length(items)
-items2 <- items
-# Un seul groupe
-if (length(group)!=1) {
-stop('Only one variable can be used for the group option')
-}
-# Recoder groupe en 0/1
-dat[,group] <- as.factor(dat[,group])
-if (!(all(names(levels(dat[,group]))==c("0","1")))) {
-levels(dat[,group]) <- c("0","1")
-}
-dat[,group] <- as.numeric(dat[,group])-1
-model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
-model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
-ltest <- lmtest::lrtest(model_a,model_b)
-if (ltest$`Pr(>Chisq)`[2]<1) {
-model_c <- model_b
-difit <- c()
-a <- sapply(1:nbitems,function(x)
-lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
-while (min(a)<0.4) {
-difit <- c(difit,which.min(a))
-model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
-print(difit)
-a <- sapply(1:nbitems,function(x)
-lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
-}
-}
-}
-rosali(dat=dat,items=1:7,group="TT")
-rosali <- function(dat=NULL,items=NULL,group=NULL) {
-nbitems <- length(items)
-items2 <- items
-# Un seul groupe
-if (length(group)!=1) {
-stop('Only one variable can be used for the group option')
-}
-# Recoder groupe en 0/1
-dat[,group] <- as.factor(dat[,group])
-if (!(all(names(levels(dat[,group]))==c("0","1")))) {
-levels(dat[,group]) <- c("0","1")
-}
-dat[,group] <- as.numeric(dat[,group])-1
-model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
-model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
-ltest <- lmtest::lrtest(model_a,model_b)
-if (ltest$`Pr(>Chisq)`[2]<1) {
-model_c <- model_b
-difit <- c()
-a <- sapply(1:nbitems,function(x)
-lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
-while (min(a)<0.4) {
-difit <- c(difit,which.min(a))
-model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
-print(difit)
-a <- sapply(1:nbitems,function(x)
-lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
-print(a)
-}
-}
-}
-rosali(dat=dat,items=1:7,group="TT")
-rosali <- function(dat=NULL,items=NULL,group=NULL) {
-nbitems <- length(items)
-items2 <- items
-# Un seul groupe
-if (length(group)!=1) {
-stop('Only one variable can be used for the group option')
-}
-# Recoder groupe en 0/1
-dat[,group] <- as.factor(dat[,group])
-if (!(all(names(levels(dat[,group]))==c("0","1")))) {
-levels(dat[,group]) <- c("0","1")
-}
-dat[,group] <- as.numeric(dat[,group])-1
-model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
-model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
-ltest <- lmtest::lrtest(model_a,model_b)
-if (ltest$`Pr(>Chisq)`[2]<1) {
-model_c <- model_b
-difit <- c()
-a <- sapply(1:nbitems,function(x)
-lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
-while (min(a)<0.4) {
-difit <- c(difit,which.min(a))
-model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
-print(difit)
-a <- sapply(1:nbitems,function(x)
-ifelse(x%in%difit,lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2],1))
-print(a)
-}
-}
-}
-rosali(dat=dat,items=1:7,group="TT")
-rosali <- function(dat=NULL,items=NULL,group=NULL) {
-nbitems <- length(items)
-items2 <- items
-# Un seul groupe
-if (length(group)!=1) {
-stop('Only one variable can be used for the group option')
-}
-# Recoder groupe en 0/1
-dat[,group] <- as.factor(dat[,group])
-if (!(all(names(levels(dat[,group]))==c("0","1")))) {
-levels(dat[,group]) <- c("0","1")
-}
-dat[,group] <- as.numeric(dat[,group])-1
-model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
-model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
-ltest <- lmtest::lrtest(model_a,model_b)
-if (ltest$`Pr(>Chisq)`[2]<1) {
-model_c <- model_b
-difit <- c()
-a <- sapply(1:nbitems,function(x)
-lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
-while (min(a)<0.4) {
-difit <- c(difit,which.min(a))
-model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
-print(difit)
-a <- sapply(1:nbitems,function(x)
-ifelse(!(x%in%difit),lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2],1))
-print(a)
-}
-}
-}
-rosali(dat=dat,items=1:7,group="TT")
-rosali <- function(dat=NULL,items=NULL,group=NULL) {
-nbitems <- length(items)
-items2 <- items
-# Un seul groupe
-if (length(group)!=1) {
-stop('Only one variable can be used for the group option')
-}
-# Recoder groupe en 0/1
-dat[,group] <- as.factor(dat[,group])
-if (!(all(names(levels(dat[,group]))==c("0","1")))) {
-levels(dat[,group]) <- c("0","1")
-}
-dat[,group] <- as.numeric(dat[,group])-1
-model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
-model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
-ltest <- lmtest::lrtest(model_a,model_b)
-difit <- c()
-if (ltest$`Pr(>Chisq)`[2]<0.05) {
-model_c <- model_b
-a <- sapply(1:nbitems,function(x)
-lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
-while (min(a)<0.05) {
-difit <- c(difit,which.min(a))
-model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
-print(difit)
-a <- sapply(1:nbitems,function(x)
-ifelse(!(x%in%difit),lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2],1))
-}
-}
-return(difit)
-}
-rosali(dat=dat,items=1:7,group="TT")
-dat <- read.csv("/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N200/scenario_20A_200.csv")
-dat <- dat[dat$replication==2,]
-rosali(dat=dat,items=1:7,group="TT")
-rosali <- function(dat=NULL,items=NULL,group=NULL) {
-nbitems <- length(items)
-items2 <- items
-# Un seul groupe
-if (length(group)!=1) {
-stop('Only one variable can be used for the group option')
-}
-# Recoder groupe en 0/1
-dat[,group] <- as.factor(dat[,group])
-if (!(all(names(levels(dat[,group]))==c("0","1")))) {
-levels(dat[,group]) <- c("0","1")
-}
-dat[,group] <- as.numeric(dat[,group])-1
-model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
-model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
-ltest <- lmtest::lrtest(model_a,model_b)
-difit <- c()
-if (ltest$`Pr(>Chisq)`[2]<0.05) {
-model_c <- model_b
-a <- sapply(1:nbitems,function(x)
-lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
-while (min(a)<0.05) {
-difit <- c(difit,which.min(a))
-model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
-a <- sapply(1:nbitems,function(x)
-ifelse(!(x%in%difit),lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2],1))
-}
-}
-return(difit)
-}
-rosali(dat=dat,items=1:7,group="TT")
-dat <- read.csv("/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N200/scenario_20A_200.csv")
-dat <- dat[dat$replication==3,]
-rosali(dat=dat,items=1:7,group="TT")
-rosali <- function(dat=NULL,items=NULL,group=NULL) {
-nbitems <- length(items)
-items2 <- items
-# Un seul groupe
-if (length(group)!=1) {
-stop('Only one variable can be used for the group option')
-}
-# Recoder groupe en 0/1
-dat[,group] <- as.factor(dat[,group])
-if (!(all(names(levels(dat[,group]))==c("0","1")))) {
-levels(dat[,group]) <- c("0","1")
-}
-dat[,group] <- as.numeric(dat[,group])-1
-model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
-model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
-ltest <- lmtest::lrtest(model_a,model_b)
-difit <- c()
-if (ltest$`Pr(>Chisq)`[2]<0.05) {
-model_c <- model_b
-a <- sapply(1:nbitems,function(x)
-lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
-k <- 0
-while (min(a)<0.05/(nbitems-k)) {
-difit <- c(difit,which.min(a))
-model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
-a <- sapply(1:nbitems,function(x)
-ifelse(!(x%in%difit),lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2],1))
-}
-}
-return(difit)
-}
-rosali(dat=dat,items=1:7,group="TT")
-pcm(dat,items=1:7,group="TT",dif.items=1:2)
-summary(pcm(dat,items=1:7,group="TT",dif.items=1:2))
+cv_error <- NULL
+}
+ret.list <- list(coefficients = coefficients, logLik = logLik,
+cv_error = cv_error, call = match.call(), model = model,
+data = data, control = control, DSF = DSF, formula = formula,
+item.names = item.names, Y = Y, design_list = design_list,
+AIC = AIC, BIC = BIC, cAIC = cAIC, df = df, coef.rescal = coef.rescal,
+main.effects = main.effects)
+class(ret.list) <- "GPCMlasso"
+return(ret.list)
+}
+GPCMlasso2(formula = cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control=ctrl_GPCMlasso(lambda=0.0000000000001),acoefs2 = u)
+u
+## create responses for acat from ordinal values
+createResponse <- function(Y){
+Y <- as.factor(Y)
+model.matrix(~0+Y)[,-length(levels(Y))]
+}
+design_GPCMlasso <- function(formula = formula, Y=Y, data = data, RSM = RSM,
+GPCM = GPCM, DSF = DSF, all.dummies = TRUE, main.effects = main.effects){
+## extract covariates
+if(all.dummies){
+term.labels <- attr(terms(formula), "term.labels")
+X <- matrix(rep(1,nrow(data)))
+for(ij in 1:length(term.labels)){
+x.now <- data[,term.labels[ij], drop = FALSE]
+if(is.factor(x.now[,1])){
+if(nlevels(x.now[,1])==2){
+X <- cbind(X,model.matrix(~.,data=x.now)[,-1,drop=FALSE])
+}else{
+X <- cbind(X,model.matrix(~0+.,data=x.now))
+}
+}else{
+X <- cbind(X,model.matrix(~.,data=x.now)[,-1,drop=FALSE])
+}
+}
+X <- X[,-1,drop = FALSE]
+}else{
+X <- model.matrix(formula, data = data)
+if(ncol(X)>=1){
+if(colnames(X)[1]=="(Intercept)"){
+X <- X[,-1,drop = FALSE]
+}
+}
+}
+x.names <- colnames(X)
+## factorize response vector
+## initialize basic parameters
+k <- apply(Y, 2, function(x){length(levels(as.factor(x)))})
+q <- k-1
+n <- nrow(Y)
+I <- ncol(Y)
+m <- ncol(X)
+n_sigma <- 1
+if(GPCM){
+n_sigma <- I
+}
+## get final response vector
+if(all(q==1)){
+response <- as.numeric(as.factor(c(t(Y))))-1
+}else{
+response <- lapply(as.data.frame(Y),createResponse)
+response <- c(t(do.call("cbind",response)))
+}
+## total number of parameters to be optimized
+px <- sum(q)+n_sigma
+if(RSM){
+px <- q[1]+I+n_sigma-1
+}
+## scale X and save standard deviations for re-scaling
+# X <- scale(X, center = FALSE)
+X <- scale(X, center = FALSE, scale = apply(X, 2, sd, na.rm = TRUE))
+sd.vec <- attributes(X)$scale
+sd.vec <- create.sd.vec(sd.vec, DSF, px, n_sigma, I, q, main.effects)
+## create design matrix for basic item parameters
+if(!RSM){
+design <- -diag(px-n_sigma)
+}else{
+design <- -cbind(matrix(rep(diag(I),each=q[1]),ncol=I),
+matrix(rep(diag(q[1]),I),ncol=q[1],byrow = TRUE))[,-(I+1)]
+}
+## create design matrix for covariate part
+if(m>=1){
+designX <- -get_designX(X, DSF, m, I, q, n)
+## create penalization matrix
+acoefs <- get_acoefs(RSM, DSF, m, I, q, n_sigma)
+## update number of parameters to be optimized
+px <- px + ncol(designX)
+}else{
+## dummy matrix in case no covariates are specified
+designX <- matrix(0,0,0)
+acoefs <- matrix(0,nrow=px,ncol=1)
+}
+if(main.effects){
+design.main <- matrix(rep(X, each = sum(q)),ncol = ncol(X))
+designX <- cbind(design.main, designX)
+acoefs <- rbind(matrix(0, ncol = ncol(acoefs), nrow = ncol(design.main)), acoefs)
+px <- px + ncol(X)
+}
+ret.list <- list(q = q, I = I, m = m, px = px, n = n, response = response,
+design = design, designX = designX, sd.vec = sd.vec,
+acoefs = acoefs, n_sigma = n_sigma, x.names = x.names,
+RSM = RSM, GPCM = GPCM, Y = Y, DSF = DSF)
+return(ret.list)
+}
+create.sd.vec <- function(sd.vec, DSF, px, n_sigma, I, q, main.effects){
+if(DSF){
+new_vec <- rep(sd.vec, sum(q))
+}else{
+new_vec <- rep(sd.vec, I)
+}
+if(main.effects){
+new_vec <- c(sd.vec, new_vec)
+}
+new_vec <- c(rep(1,px-n_sigma),new_vec,rep(1,n_sigma))
+new_vec
+}
+GPCMlasso2(formula = cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control=ctrl_GPCMlasso(lambda=0.0000000000001),acoefs2 = u)
+get_designX <- function(X, DSF, m, I, q, n){
+if(!DSF){
+designX <- matrix(0, ncol = m * I, nrow = n * sum(q))
+pos_u <- 1
+for (u in 1:n) {
+for(uu in 1:I){
+designX[pos_u:(pos_u+q[uu]-1), ((uu - 1) * m + 1):(uu * m)] <-
+matrix(rep(X[u,], q[uu]), byrow = TRUE, ncol = m, nrow = q[uu])
+pos_u <- pos_u+q[uu]
+}
+}
+}else{
+designX <- matrix(0, ncol = m * sum(q), nrow = n * sum(q))
+pos_u <- 1
+for (u in 1:n) {
+pos_uu <- 1
+for(uu in 1:I){
+for(uuu in 1:q[uu]){
+designX[pos_u, pos_uu:(pos_uu+m-1)] <-
+X[u,]
+pos_u <- pos_u+1
+pos_uu <- pos_uu+m
+}
+}
+}
+}
+return(designX)
+}
+get_acoefs_old <- function(RSM, DSF, m, I, q, n_sigma){
+if(!DSF){
+pen1 <- diag(m*I)
+}else{
+pen1 <- matrix(0,nrow=m*sum(q),ncol=m*sum(choose(q,2)))
+pos1 <-1
+pos_pos <- 1
+for(u in 1:I){
+n_comb <- choose(q[u],2)
+if(n_comb>0){
+combis <- combn(q[u],2)-1
+for(uuu in 1:m){
+for(uu in 1:n_comb){
+pen1[combis[1,uu]*m+pos_pos,pos1] <- 1
+pen1[combis[2,uu]*m+pos_pos,pos1] <- -1
+pos1 <- pos1+1
+}
+pos_pos <- pos_pos+1
+}
+pos_pos <- pos_pos+(q[u]-1)*m
+}
+}
+pen1 <- cbind(diag(m*sum(q)),pen1)
+}
+if(RSM){
+acoefs <- rbind(matrix(0,nrow=q[1]+I-1,ncol=ncol(pen1)),pen1,
+matrix(0, ncol = ncol(pen1), nrow = n_sigma))
+}else{
+acoefs <- rbind(matrix(0,nrow=sum(q),ncol=ncol(pen1)),pen1,
+matrix(0, ncol = ncol(pen1), nrow = n_sigma))
+}
+}
+get_acoefs <- function(RSM, DSF, m, I, q, n_sigma){
+if(!DSF){
+pen1 <- diag(m*I)
+}else{
+pen1 <- matrix(0,nrow=m*sum(q),ncol=m*sum(q-1))
+pos1 <-1
+pos_pos <- 1
+for(u in 1:I){
+n_comb <- q[u] - 1
+if(n_comb>0){
+combis <- rep(1:q[u], each = 2)
+combis <- matrix(combis[-c(1,length(combis))]-1, nrow = 2)
+for(uuu in 1:m){
+for(uu in 1:n_comb){
+pen1[combis[1,uu]*m+pos_pos,pos1] <- 1
+pen1[combis[2,uu]*m+pos_pos,pos1] <- -1
+pos1 <- pos1+1
+}
+pos_pos <- pos_pos+1
+}
+pos_pos <- pos_pos+(q[u]-1)*m
+}
+}
+pen1 <- cbind(diag(m*sum(q)),pen1)
+}
+if(RSM){
+acoefs <- rbind(matrix(0,nrow=q[1]+I-1,ncol=ncol(pen1)),pen1,
+matrix(0, ncol = ncol(pen1), nrow = n_sigma))
+}else{
+acoefs <- rbind(matrix(0,nrow=sum(q),ncol=ncol(pen1)),pen1,
+matrix(0, ncol = ncol(pen1), nrow = n_sigma))
+}
+}
+GPCMlasso2(formula = cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control=ctrl_GPCMlasso(lambda=0.0000000000001),acoefs2 = u)
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+loglikPCMlasso <- function(alpha, Y, X, Z, Q, q, n, I, px, GHweights, GHnodes, acoefs, lambda, lambda2, cvalue, cores, weight, n_sigma, scale_fac) {
+.Call(`_GPCMlasso_loglikPCMlasso`, alpha, Y, X, Z, Q, q, n, I, px, GHweights, GHnodes, acoefs, lambda, lambda2, cvalue, cores, weight, n_sigma, scale_fac)
+}
+scorePCMlasso <- function(alpha, Y, X, Z, Q, q, n, I, px, GHweights, GHnodes, acoefs, lambda, lambda2, cvalue, cores, weight, n_sigma, scale_fac) {
+.Call(`_GPCMlasso_scorePCMlasso`, alpha, Y, X, Z, Q, q, n, I, px, GHweights, GHnodes, acoefs, lambda, lambda2, cvalue, cores, weight, n_sigma, scale_fac)
+}
+loglikDIFlasso <- function(alpha, Y, X, Z, Q, q, n, I, px, GHweights, GHnodes, acoefs, lambda, lambda2, cvalue, cores, weight, n_sigma, scale_fac) {
+.Call(`_GPCMlasso_loglikDIFlasso`, alpha, Y, X, Z, Q, q, n, I, px, GHweights, GHnodes, acoefs, lambda, lambda2, cvalue, cores, weight, n_sigma, scale_fac)
+}
+scoreDIFlasso <- function(alpha, Y, X, Z, Q, q, n, I, px, GHweights, GHnodes, acoefs, lambda, lambda2, cvalue, cores, weight, n_sigma, scale_fac) {
+.Call(`_GPCMlasso_scoreDIFlasso`, alpha, Y, X, Z, Q, q, n, I, px, GHweights, GHnodes, acoefs, lambda, lambda2, cvalue, cores, weight, n_sigma, scale_fac)
+}
+loglikscorePCMlasso <- function(alpha, Y, X, Z, Q, q, n, I, px, GHweights, GHnodes, acoefs, lambda, lambda2, cvalue, cores, weight, n_sigma, scale_fac) {
+.Call(`_GPCMlasso_loglikscorePCMlasso`, alpha, Y, X, Z, Q, q, n, I, px, GHweights, GHnodes, acoefs, lambda, lambda2, cvalue, cores, weight, n_sigma, scale_fac)
+}
+loglikscoreDIFlasso <- function(alpha, Y, X, Z, Q, q, n, I, px, GHweights, GHnodes, acoefs, lambda, lambda2, cvalue, cores, weight, n_sigma, scale_fac) {
+.Call(`_GPCMlasso_loglikscoreDIFlasso`, alpha, Y, X, Z, Q, q, n, I, px, GHweights, GHnodes, acoefs, lambda, lambda2, cvalue, cores, weight, n_sigma, scale_fac)
+}
+GPCMlasso2(formula = cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control=ctrl_GPCMlasso(lambda=0.0000000000001),acoefs2 = u)
+library(GPCMlasso)
+trace(GPCMlasso,edit=T)
+library(GPCMlasso)
+View(GPCMlasso)
+u
+aaa
+aaa <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Data/NoDIF/N',100,'/scenario_',"5A_100",'.csv'))
+aaa <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N',100,'/scenario_',"5A_100",'.csv'))
+aaaa <- aaa[aaa$replication==1,]
+GPCMlasso(formula=cbind('item1','item2',"item3","item4")~TT,data=aaaa,model="PCM")
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM")
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.00001))
+acc <- matrix(nrow = 10,ncol=4,0)
+acc
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.00001),acoefs2 = acc)
+acc[6,1]1
+acc[6,1]<-1
+acc
+acc[7,2]<-1
+acc[8,3]<-1
+acc[9,4]<-1
+acc
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.00001),acoefs2 = acc)
 library(TAM)
-resali <- function(df=NULL,items=NULL,group=NULL,verbose=T) {
-if (verbose) {
-cat('-----------------------------------------------------------\n')
-cat('COMPUTING INITIAL PCM\n')
-cat('-----------------------------------------------------------\n')
-}
-nbitems <- length(items)
-nbitems_o <- nbitems
-items_n <- paste0('item',items)
-resp <- df[,items_n]
-grp <- df[,group]
-pcm_initial <- TAM::tam.mml(resp=resp,Y=grp,irtmodel = "PCM",est.variance = T,verbose=F)
-dat <- df
-dat$score <- rowSums(dat[,items_n])
-dat$score_q5 <- cut(dat$score,quantile(dat$score,seq(0,1,0.2)),labels=1:5,include.lowest=T)
-res.anova <- rep(NA,nbitems)
-pval <- rep(NA,nbitems)
-fval <- rep(NA,nbitems)
-for (i in items) {
-dat[,paste0('res_',i)] <- rowMeans(predict(pcm_initial)$stand.resid[,,i])
-res.anova[i] <- summary(aov(dat[,paste0('res_',i)]~TT*score_q5,data=dat))
-pval[i] <- res.anova[[i]][1,"Pr(>F)"]
-fval[i] <- res.anova[[i]][1,'F value']
-}
-if (verbose) {
-cat('DONE\n')
-cat('-----------------------------------------------------------\n')
-}
-res.items <- c()
-res.uniform <- c()
-k <- 1
-while(any(pval<0.05/nbitems_o)) {
-k <- k+1
-if (verbose) {
-cat(paste('COMPUTING STEP',k,'\n'))
-cat('-----------------------------------------------------------\n')
-}
-res.item <- gsub("[a-z]", "",colnames(resp)[which.max(fval)])
-res.items <- c(res.items,res.item)
-res.uni <- res.anova[[which.max(fval)]][3,"Pr(>F)"]>0.05
-res.uniform <- c(res.uniform,res.uni)
-items_n <- c(items_n[items_n!=paste0('item',res.item)],paste0("item",res.item,c("noTT","TT")))
-dat[,paste0("item",res.item,'TT')] <- dat[dat$TT==1,paste0('item',res.item)]
-dat[,paste0("item",res.item,'noTT')] <- dat[dat$TT==0,paste0('item',res.item)]
-resp <- dat[,items_n]
-grp <- dat[,group]
-pcm_while <- TAM::tam.mml(resp=resp,Y=grp,irtmodel = "PCM",est.variance = T,verbose=F)
-nbitems <- length(items_n)
-res.anova <- rep(NA,nbitems)
-pval <- rep(NA,nbitems)
-fval <- rep(NA,nbitems)
-for (i in 1:nbitems) {
-dat[,paste0('res_',i)] <- rowMeans(predict(pcm_while)$stand.resid[,,i])
-res.anova[i] <- summary(aov(dat[,paste0('res_',i)]~TT*score_q5,data=dat))
-pval[i] <- res.anova[[i]][1,"Pr(>F)"]
-fval[i] <- res.anova[[i]][1,'F value']
-}
-if (verbose) {
-cat('DONE\n')
-cat('-----------------------------------------------------------\n')
-}
-}
-if (verbose) {
-cat("DETECTED DIF ITEMS\n")
-cat('-----------------------------------------------------------\n')
-}
-if (length(res.items>0)) {
-results <- data.frame(dif.items=res.items,
-uniform=1*res.uniform)
-return(results)
-}
-else {
-cat("No DIF was detected")
-return(NULL)
-}
-}
-resali(df=dat,items=1:7,group="TT")
-dat$dif1
-dat$dif2
-dat$dif3
+tam.mml(resp=aaaa[,c('item1',"item2","item3","item4")],group=aaaa[,'TT'])
+zzz <- tam.mml(resp=aaaa[,c('item1',"item2","item3","item4")],group=aaaa[,'TT'])
+summary(zzz)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.00001),acoefs2 = acc)
+acc <- matrix(nrow = 10,ncol=4,0)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.00001),acoefs2 = acc)
+summary(zzz)
+remove.packages(GPCMlasso)
+remove.packages("GPCMlasso")
+install.packages("GPCMlasso")
+install.packages("GPCMlasso")
+library(GPCMlasso)
+View(GPCMlasso)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.00001))
+tam.mml(resp=aaaa[,c('item1',"item2","item3","item4")],group=aaaa[,'TT'])
+library(tam.mml)
+library(TAM)
+tam.mml(resp=aaaa[,c('item1',"item2","item3","item4")],group=aaaa[,'TT'])
+tam.mml(resp=aaaa[,c('item1',"item2","item3","item4")],group=aaaa[,'TT'])
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.00001))
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.000000000000000000001))
+remove.packages("GPCMlasso")
+install.packages('/home/corentin/Documents/These/Recherche/GPCMlasso_0.1-7.tar.gz')
+library(GPCMlasso)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.000000000000000000001))
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0000000000000000000000000000000001))
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0000000000000000000000000000000001),acoefs=acc)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0000000000000000000000000000000001),acoefs2=acc)
+View(GPCMlasso)
+library(GPCMlasso)
+View(GPCMlasso)
+aaa <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Data/NoDIF/N',100,'/scenario_',"5A_100",'.csv'))
+aaa <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N',100,'/scenario_',"5A_100",'.csv'))
+aaaa <- aaa[aaa$replication==1,]
+acc <- matrix(nrow = 10,ncol=4,0)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0000000000000000000000000000000001),acoefs2=acc)
+acc[6,1] <- 1;acc[7,2] <- 1;acc[8,3] <- 1;acc[9,4] <- 1
+acc
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0000000000000000000000000000000001),acoefs2=acc)
+remove.packages("GPCMlasso")
+install.packages('/home/corentin/Documents/These/Recherche/GPCMlasso_0.1-7.tar.gz')
+library(GPCMlasso)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0000000000000000000000000000000001))
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 100))
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 1000))
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0000000000000000000000000000000001))
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0000000000000000000000000000000001))$coefficients
+install.packages("geoR")
+library(geoR)
+nlmP
+nlmP()
+.nlmP()
+geoR::.nlmP
+remove.packages("GPCMlasso")
+install.packages('/home/corentin/Documents/These/Recherche/GPCMlasso_0.1-7.tar.gz')
+install.packages('/home/corentin/Documents/These/Recherche/GPCMlasso_0.1-7.tar.gz')
+library(GPCMlasso)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0000000000000000000000000000000001))
+remove.packages("GPCMlasso")
+install.packages('/home/corentin/Documents/These/Recherche/GPCMlasso_0.1-7.tar.gz')
+library(GPCMlasso)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0000000000000000000000000000000001))
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0000000000000000000000000000000001))$design_list
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0000000000000000000000000000000001),DSF = T)
+summary(GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0000000000000000000000000000000001),DSF = T))
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0000000000000000000000000000000001),DSF = T)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0000000000000000000000000000000001),DSF = F)
+library(TAM)
+tam.mml(resp = aaaa[,c("item1","item2","item3","item4")],group=aaaa[,"TT"],irtmodel="PCM2")
+summary(tam.mml(resp = aaaa[,c("item1","item2","item3","item4")],group=aaaa[,"TT"],irtmodel="PCM2"))
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0000000000000000000000000000000001),DSF = F)
+summary(tam.mml(resp = aaaa[,c("item1","item2","item3","item4")],group=aaaa[,"TT"],irtmodel="PCM2"))
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0000000000000000000000000000000001),DSF = F)
+summary(tam.mml(resp = aaaa[,c("item1","item2","item3","item4")],group=aaaa[,"TT"],irtmodel="PCM2"))
+summary(tam.jml(resp = aaaa[,c("item1","item2","item3","item4")],group=aaaa[,"TT"],irtmodel="PCM2"))
+summary(tam.jml(resp = aaaa[,c("item1","item2","item3","item4")],group=aaaa[,"TT"]))
+summary(tam(resp = aaaa[,c("item1","item2","item3","item4")],group=aaaa[,"TT"]))
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0000000000000000000000000000000001),DSF = F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0000000000000000000000000000000001),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0000000000000000000001),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 10),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 1),DSF = F,cv=F)
+summary(tam(resp = aaaa[,c("item1","item2","item3","item4")],group=aaaa[,"TT"]))
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 1),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.5),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.1),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.01),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.001),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0001),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.00000001),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0000001),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.01),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.00001),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.00000001),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0000001),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0.0000005),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0),DSF = F,cv=F)$design_list
+install.packages('/home/corentin/Documents/These/Recherche/GPCMlasso_0.1-7.tar.gz')
+library(GPCMlasso)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0),DSF = F,cv=F)
+remove.packages('GPCMlasso')
+install.packages('/home/corentin/Documents/These/Recherche/GPCMlasso_0.1-7.tar.gz')
+library(GPCMlasso)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0),DSF = F,cv=F)$design_list
+remove.packages("PCMlasso")
+remove.packages("GPCMlasso")
+install.packages('/home/corentin/Documents/These/Recherche/GPCMlasso_0.1-7.tar.gz')
+library(GPCMlasso)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0),DSF = F,cv=F)$design_list
+library(GPCMlasso)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0),DSF = F,cv=F)$design_list
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0),DSF = F,cv=F)
+remove.packages("GPCMlasso")
+install.packages('/home/corentin/Documents/These/Recherche/GPCMlasso_0.1-7.tar.gz')
+library(GPCMlasso)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0),DSF = F,cv=F)$design_list
+remove.packages("GPCMlasso")
+install.packages('/home/corentin/Documents/These/Recherche/GPCMlasso_0.1-7.tar.gz')
+library(GPCMlasso)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT:item1,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~I(TT:item1*0),data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT:item1*0,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0),DSF = F,cv=F)
+aaaa$item1tt <- 1-aaaa$item1
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT:item1tt,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0),DSF = F,cv=F)
+aaaa$item2tt <- 1-aaaa$item2
+aaaa$item3tt <- 1-aaaa$item3
+aaaa$item1
+aaaa$item4tt <- 1-aaaa$item4
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT:item1tt+TT:item2TT+TT:item3tt+TT:item4tt,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT:item1tt+TT:item2tt+TT:item3tt+TT:item4tt,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT*item1tt,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0),DSF = F,cv=F)
+GPCMlasso(formula=cbind(item1,item2,item3,item4)~TT*item1tt+TT*item2tt+TT*item3tt+TT*item4tt,data=aaaa,model="PCM",control = ctrl_GPCMlasso(lambda = 0),DSF = F,cv=F)
diff --git a/RProject/Scripts/pcm.R b/RProject/Scripts/pcm.R
index 4d83804..d15e24e 100644
--- a/RProject/Scripts/pcm.R
+++ b/RProject/Scripts/pcm.R
@@ -2050,7 +2050,7 @@ write.csv(res_nodif[[3]],'/home/corentin/Documents/These/Recherche/Simulations/A
 
 ##############################################################################
 #----------------------------------------------------------------------------#
-########################## SCENARIO ANALYSIS N=50 ##########################
+########################### SCENARIO ANALYSIS N=50 ###########################
 #----------------------------------------------------------------------------#
 ##############################################################################
 
@@ -2706,6 +2706,374 @@ for (x in results[seq(2,length(results))]) {
 res.dat$bias <- res.dat$eff.size-res.dat$m.beta
 res.dat.dif$bias <- res.dat.dif$eff.size-res.dat.dif$m.beta
 
+##############################################################################
+#----------------------------------------------------------------------------#
+####################### AGGREGATION DIF MATRICES ROSALI ######################
+#----------------------------------------------------------------------------#
+##############################################################################
+
+#### Create data.frame
+
+results <- c(sapply(1:4,function(x) paste0(x,c('A','B','C','D','E'))),sapply(5:9,function(x) paste0(x,c('A','B','C','D','E','F','G'))))
+
+results2 <- c(sapply(10:20,function(x) paste0(x,c('A','B','C','D','E','F','G'))))
+
+results <- c(sapply(c(50,100,200,300),function(x) paste0(results,'_',x)))
+
+results2 <- c(sapply(c(50,100,200,300),function(x) paste0(results2,'_',x)))
+
+results <- sort(results)
+
+results2 <- sort(results2)
+
+results <- c(results,results2)
+
+
+#### Compiler function
+
+compile_simulation2_rosali <- function(scenario) {
+  name <- as.numeric(gsub("[^0-9.-]", "", substr(scenario,start=0,stop=2)))
+  if (substr(scenario,start=nchar(scenario)-1,stop=nchar(scenario))=="50" & name>0) {
+    s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/ROSALI-DIF/N50/',scenario,'_original.xls'))
+  }
+  if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="100" & name>0) {
+    s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/ROSALI-DIF/N100/',scenario,'_original.xls'))
+  }
+  if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="200" & name>0) {
+    s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/ROSALI-DIF/N200/',scenario,'_original.xls'))
+  }
+  if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="300" & name>0) {
+    s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/ROSALI-DIF/N300/',scenario,'_original.xls'))
+  }
+  J <- max(which(sapply(1:7,function(x) paste0('item',x) %in% colnames(s) | paste0('item',x,'_1') %in% colnames(s))))
+  M <- 1+sum(sapply(1:3,function(x) paste0('item1_',x) %in% colnames(s) ))
+  if (M==1) {M <- 2}
+  nb.dif <- max(which(sapply(1:3,function(x) paste0('dif',x) %in% colnames(s) | paste0('dif',x,'_1') %in% colnames(s))))
+  if (J==4) {
+    if (M==2) {
+      a <- data.frame(m.item1=mean(s$item1),m.item2=mean(s$item2),m.item3=mean(s$item3),m.item4=mean(s$item4))
+    } else {
+      a <- data.frame(m.item1_1=mean(s$item1_1),m.item1_2=mean(s$item1_2),m.item1_3=mean(s$item1_3),
+                      m.item2_1=mean(s$item2_1),m.item2_2=mean(s$item2_2),m.item2_3=mean(s$item2_3),
+                      m.item3_1=mean(s$item3_1),m.item3_2=mean(s$item3_2),m.item3_3=mean(s$item3_3),
+                      m.item4_1=mean(s$item4_1),m.item4_2=mean(s$item4_2),m.item4_3=mean(s$item4_3)
+      )
+    }
+  } else {
+    if (M==2) {
+      a <- data.frame(m.item1=mean(s$item1),m.item2=mean(s$item2),m.item3=mean(s$item3),m.item4=mean(s$item4),
+                      m.item5=mean(s$item5),m.item6=mean(s$item6),m.item7=mean(s$item7))
+    } else {
+      a <- data.frame(m.item1_1=mean(s$item1_1),m.item1_2=mean(s$item1_2),m.item1_3=mean(s$item1_3),
+                      m.item2_1=mean(s$item2_1),m.item2_2=mean(s$item2_2),m.item2_3=mean(s$item2_3),
+                      m.item3_1=mean(s$item3_1),m.item3_2=mean(s$item3_2),m.item3_3=mean(s$item3_3),
+                      m.item4_1=mean(s$item4_1),m.item4_2=mean(s$item4_2),m.item4_3=mean(s$item4_3),
+                      m.item5_1=mean(s$item5_1),m.item5_2=mean(s$item5_2),m.item5_3=mean(s$item5_3),
+                      m.item6_1=mean(s$item6_1),m.item6_2=mean(s$item6_2),m.item6_3=mean(s$item6_3),
+                      m.item7_1=mean(s$item7_1),m.item7_2=mean(s$item7_2),m.item7_3=mean(s$item7_3)
+      )
+    }
+  }
+  N <- ifelse(substr(scenario,start=nchar(scenario)-1,stop=nchar(scenario))=="50","50",substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario)))
+  zz <- ifelse(N=="50",substr(scenario,start=0,stop=nchar(scenario)-3),substr(scenario,start=0,stop=nchar(scenario)-4))
+  eff.size <- unique(res.dat[res.dat$scenario==zz & res.dat$N==N,'eff.size'])
+  dif.size <- unique(res.dat[res.dat$scenario==zz & res.dat$N==N,'dif.size'])
+  b <- data.frame(scenario=zz,
+                  scenario.type=substr(zz,start=nchar(zz),stop=nchar(zz)),
+                  N=N,
+                  J=J,
+                  M=M,
+                  eff.size=eff.size,
+                  nb.dif=nb.dif,
+                  dif.size=dif.size
+  )
+  true.value.in.ci <- eff.size <= s$beta+1.96*s$se_beta & eff.size >= s$beta-1.96*s$se_beta
+  beta.same.sign.truebeta.p <- ifelse(rep(eff.size,nrow(s))==0,NA,(rep(eff.size,nrow(s))/s$beta)>0)
+  num.reject <- which((s$beta-1.96*s$se_beta)>0 | (s$beta+1.96*s$se_beta)<0)
+  z <- data.frame(m.beta=mean(s$beta),
+                  se.empirical.beta=sd(s$beta),
+                  se.analytical.beta=mean(s$se_beta),
+                  m.low.ci.beta=mean(s$beta-1.96*s$se_beta),
+                  m.high.ci.beta=mean(s$beta+1.96*s$se_beta),
+                  true.value.in.ci.p=mean(true.value.in.ci),
+                  h0.rejected.p=mean( (s$beta-1.96*s$se_beta)>0 | (s$beta+1.96*s$se_beta)<0 ),
+                  beta.same.sign.truebeta.p=mean(beta.same.sign.truebeta.p),
+                  beta.same.sign.truebeta.signif.p=mean(beta.same.sign.truebeta.p[num.reject])
+  )
+  d <- cbind(b,a,z)
+  d$prop.
+  return(d)
+}
+
+
+#### Compiled results
+
+res.dat.dif.rosali <- compile_simulation2_rosali('1A_100')
+
+for (x in results[seq(2,length(results))]) {
+  y <- compile_simulation2_rosali(x)
+  res.dat.dif.rosali <- bind_rows(res.dat.dif.rosali,y)
+}
+
+res.dat.dif.rosali$bias <- res.dat.dif.rosali$eff.size-res.dat.dif.rosali$m.beta
+
+
+##############################################################################
+#----------------------------------------------------------------------------#
+################################### RESALI ###################################
+#----------------------------------------------------------------------------#
+##############################################################################
+
+generate_resali <- function(scenario=NULL,grp=NULL) {
+  scen <- as.numeric(gsub("[A,B,C,D,E,F,G,_]","",substr(scenario,0,3)))
+  if (substr(scenario,start=nchar(scenario)-1,stop=nchar(scenario))=="50") {
+    N <- 50
+  }
+  if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="100") {
+    N <- 100
+  }
+  if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="200") {
+    N <- 200
+  }
+  if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="300") {
+    N <- 300
+  }
+  if (scen<5) {
+    dat <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Data/NoDIF/N',N,'/scenario_',scenario,'.csv'))
+  }
+  if (scen>=5) {
+    dat <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N',N,'/scenario_',scenario,'.csv'))
+  }
+  if (scen%in%c(3,4,13:20)) {
+    res <- resali(df=dat[dat$replication==1,],items = seq(1,7),group=grp,verbose=FALSE)
+    df_res <- data.frame(dif.detect.1=ifelse(length(res$dif.items)>=1,res$dif.items[1],NA),
+                         dif.detect.2=ifelse(length(res$dif.items)>=2,res$dif.items[2],NA),
+                         dif.detect.3=ifelse(length(res$dif.items)>=3,res$dif.items[3],NA),
+                         dif.detect.4=ifelse(length(res$dif.items)>=4,res$dif.items[4],NA),
+                         dif.detect.5=ifelse(length(res$dif.items)>=5,res$dif.items[5],NA),
+                         dif.detect.6=ifelse(length(res$dif.items)>=6,res$dif.items[6],NA),
+                         dif.detect.7=ifelse(length(res$dif.items)>=7,res$dif.items[7],NA),
+                         dif.detect.unif.1=ifelse(length(res$uniform)>=1,res$uniform[1],NA),
+                         dif.detect.unif.2=ifelse(length(res$uniform)>=2,res$uniform[2],NA),
+                         dif.detect.unif.3=ifelse(length(res$uniform)>=3,res$uniform[3],NA),
+                         dif.detect.unif.4=ifelse(length(res$uniform)>=4,res$uniform[4],NA),
+                         dif.detect.unif.5=ifelse(length(res$uniform)>=5,res$uniform[5],NA),
+                         dif.detect.unif.6=ifelse(length(res$uniform)>=6,res$uniform[6],NA),
+                         dif.detect.unif.7=ifelse(length(res$uniform)>=7,res$uniform[7],NA),
+                         N=N,
+                         nbdif=ifelse(scen<=4,0,ifelse(scen<=16,2,3)),
+                         true.dif.1=ifelse(scen<=4,NA,unique(dat[dat$replication==1,]$dif1)),
+                         true.dif.2=ifelse(scen<=4,NA,unique(dat[dat$replication==1,]$dif2)),
+                         true.dif.3=ifelse(scen<=16,NA,unique(dat[dat$replication==1,]$dif3))
+                         )
+    for (k in 2:1000) {
+      if (k%%100==0) {
+        cat(paste0('N=',k,'/1000\n'))
+      }
+      res <- resali(df=dat[dat$replication==k,],items = seq(1,7),group=grp,verbose=FALSE)
+      df_res2 <- data.frame(dif.detect.1=ifelse(length(res$dif.items)>=1,res$dif.items[1],NA),
+                            dif.detect.2=ifelse(length(res$dif.items)>=2,res$dif.items[2],NA),
+                            dif.detect.3=ifelse(length(res$dif.items)>=3,res$dif.items[3],NA),
+                            dif.detect.4=ifelse(length(res$dif.items)>=4,res$dif.items[4],NA),
+                            dif.detect.5=ifelse(length(res$dif.items)>=5,res$dif.items[5],NA),
+                            dif.detect.6=ifelse(length(res$dif.items)>=6,res$dif.items[6],NA),
+                            dif.detect.7=ifelse(length(res$dif.items)>=7,res$dif.items[7],NA),
+                            dif.detect.unif.1=ifelse(length(res$uniform)>=1,res$uniform[1],NA),
+                            dif.detect.unif.2=ifelse(length(res$uniform)>=2,res$uniform[2],NA),
+                            dif.detect.unif.3=ifelse(length(res$uniform)>=3,res$uniform[3],NA),
+                            dif.detect.unif.4=ifelse(length(res$uniform)>=4,res$uniform[4],NA),
+                            dif.detect.unif.5=ifelse(length(res$uniform)>=5,res$uniform[5],NA),
+                            dif.detect.unif.6=ifelse(length(res$uniform)>=6,res$uniform[6],NA),
+                            dif.detect.unif.7=ifelse(length(res$uniform)>=7,res$uniform[7],NA),
+                            N=N,
+                            nbdif=ifelse(scen<=4,0,ifelse(scen<=16,2,3)),
+                            true.dif.1=ifelse(scen<=4,NA,unique(dat[dat$replication==k,]$dif1)),
+                            true.dif.2=ifelse(scen<=4,NA,unique(dat[dat$replication==k,]$dif2)),
+                            true.dif.3=ifelse(scen<=16,NA,unique(dat[dat$replication==k,]$dif3)))
+      df_res <- rbind(df_res,df_res2)
+    }
+  }
+  else if (scen%in%c(1,2,5:12)) {
+    res <- resali(df=dat[dat$replication==1,],items = seq(1,4),group=grp,verbose=FALSE)
+    df_res <- data.frame(dif.detect.1=ifelse(length(res$dif.items)>=1,res$dif.items[1],NA),
+                         dif.detect.2=ifelse(length(res$dif.items)>=2,res$dif.items[2],NA),
+                         dif.detect.3=ifelse(length(res$dif.items)>=3,res$dif.items[3],NA),
+                         dif.detect.4=ifelse(length(res$dif.items)>=4,res$dif.items[4],NA),
+                         dif.detect.unif.1=ifelse(length(res$uniform)>=1,res$uniform[1],NA),
+                         dif.detect.unif.2=ifelse(length(res$uniform)>=2,res$uniform[2],NA),
+                         dif.detect.unif.3=ifelse(length(res$uniform)>=3,res$uniform[3],NA),
+                         dif.detect.unif.4=ifelse(length(res$uniform)>=4,res$uniform[4],NA),
+                         N=N,
+                         nbdif=ifelse(scen<=4,0,ifelse(scen<=8,1,2)),
+                         true.dif.1=ifelse(scen<=4,NA,unique(dat[dat$replication==1,]$dif1)),
+                         true.dif.2=ifelse(scen<=8,NA,unique(dat[dat$replication==1,]$dif2))
+    )
+    for (k in 2:1000) {
+      if (k%%100==0) {
+        cat(paste0('N=',k,'/1000\n'))
+      }
+      res <- resali(df=dat[dat$replication==k,],items = seq(1,4),group=grp,verbose=FALSE)
+      df_res2 <- data.frame(dif.detect.1=ifelse(length(res$dif.items)>=1,res$dif.items[1],NA),
+                            dif.detect.2=ifelse(length(res$dif.items)>=2,res$dif.items[2],NA),
+                            dif.detect.3=ifelse(length(res$dif.items)>=3,res$dif.items[3],NA),
+                            dif.detect.4=ifelse(length(res$dif.items)>=4,res$dif.items[4],NA),
+                            dif.detect.unif.1=ifelse(length(res$uniform)>=1,res$uniform[1],NA),
+                            dif.detect.unif.2=ifelse(length(res$uniform)>=2,res$uniform[2],NA),
+                            dif.detect.unif.3=ifelse(length(res$uniform)>=3,res$uniform[3],NA),
+                            dif.detect.unif.4=ifelse(length(res$uniform)>=4,res$uniform[4],NA),
+                            N=N,
+                            nbdif=ifelse(scen<=4,0,ifelse(scen<=8,1,2)),
+                            true.dif.1=ifelse(scen<=4,NA,unique(dat[dat$replication==k,]$dif1)),
+                            true.dif.2=ifelse(scen<=8,NA,unique(dat[dat$replication==k,]$dif2)))
+      df_res <- rbind(df_res,df_res2)
+    }
+  }
+  return(df_res)
+}
+
+
+
+
+
+
+results <- c(sapply(1:4,function(x) paste0(x,c('A','B','C','D','E'))),sapply(5:9,function(x) paste0(x,c('A','B','C','D','E','F','G'))))
+
+results2 <- c(sapply(10:20,function(x) paste0(x,c('A','B','C','D','E','F','G'))))
+
+results <- c(sapply(c(50,100,200,300),function(x) paste0(results,'_',x)))
+
+results2 <- c(sapply(c(50,100,200,300),function(x) paste0(results2,'_',x)))
+
+results <- sort(results)
+
+results2 <- sort(results2)
+
+results <- c(results,results2)
+
+for (r in results) {
+  cat(paste0(r,"\n"))
+  cat(paste0("-------------------------------------------","\n"))
+  write.csv(generate_resali(r,"TT"),paste0("/home/corentin/Documents/These/Recherche/Simulations/Analysis/RESALI/Detection/",r,".csv"))
+  cat(paste0("-------------------------------------------","\n"))
+}
+
+
+
+
+##############################################################################
+#----------------------------------------------------------------------------#
+####################### AGGREGATION DIF MATRICES RESALI ######################
+#----------------------------------------------------------------------------#
+##############################################################################
+
+#### Create data.frame
+
+results <- c(sapply(1:4,function(x) paste0(x,c('A','B','C','D','E'))),sapply(5:9,function(x) paste0(x,c('A','B','C','D','E','F','G'))))
+
+results2 <- c(sapply(10:20,function(x) paste0(x,c('A','B','C','D','E','F','G'))))
+
+results <- c(sapply(c(50,100,200,300),function(x) paste0(results,'_',x)))
+
+results2 <- c(sapply(c(50,100,200,300),function(x) paste0(results2,'_',x)))
+
+results <- sort(results)
+
+results2 <- sort(results2)
+
+results <- c(results,results2)
+
+
+#### Compiler function
+
+compile_simulation2_resali <- function(scenario) {
+  name <- as.numeric(gsub("[^0-9.-]", "", substr(scenario,start=0,stop=2)))
+  if (substr(scenario,start=nchar(scenario)-1,stop=nchar(scenario))=="50" & name>0) {
+    s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/resali-DIF/N50/',scenario,'_original.xls'))
+  }
+  if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="100" & name>0) {
+    s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/resali-DIF/N100/',scenario,'_original.xls'))
+  }
+  if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="200" & name>0) {
+    s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/resali-DIF/N200/',scenario,'_original.xls'))
+  }
+  if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="300" & name>0) {
+    s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/resali-DIF/N300/',scenario,'_original.xls'))
+  }
+  J <- max(which(sapply(1:7,function(x) paste0('item',x) %in% colnames(s) | paste0('item',x,'_1') %in% colnames(s))))
+  M <- 1+sum(sapply(1:3,function(x) paste0('item1_',x) %in% colnames(s) ))
+  if (M==1) {M <- 2}
+  nb.dif <- max(which(sapply(1:3,function(x) paste0('dif',x) %in% colnames(s) | paste0('dif',x,'_1') %in% colnames(s))))
+  if (J==4) {
+    if (M==2) {
+      a <- data.frame(m.item1=mean(s$item1),m.item2=mean(s$item2),m.item3=mean(s$item3),m.item4=mean(s$item4))
+    } else {
+      a <- data.frame(m.item1_1=mean(s$item1_1),m.item1_2=mean(s$item1_2),m.item1_3=mean(s$item1_3),
+                      m.item2_1=mean(s$item2_1),m.item2_2=mean(s$item2_2),m.item2_3=mean(s$item2_3),
+                      m.item3_1=mean(s$item3_1),m.item3_2=mean(s$item3_2),m.item3_3=mean(s$item3_3),
+                      m.item4_1=mean(s$item4_1),m.item4_2=mean(s$item4_2),m.item4_3=mean(s$item4_3)
+      )
+    }
+  } else {
+    if (M==2) {
+      a <- data.frame(m.item1=mean(s$item1),m.item2=mean(s$item2),m.item3=mean(s$item3),m.item4=mean(s$item4),
+                      m.item5=mean(s$item5),m.item6=mean(s$item6),m.item7=mean(s$item7))
+    } else {
+      a <- data.frame(m.item1_1=mean(s$item1_1),m.item1_2=mean(s$item1_2),m.item1_3=mean(s$item1_3),
+                      m.item2_1=mean(s$item2_1),m.item2_2=mean(s$item2_2),m.item2_3=mean(s$item2_3),
+                      m.item3_1=mean(s$item3_1),m.item3_2=mean(s$item3_2),m.item3_3=mean(s$item3_3),
+                      m.item4_1=mean(s$item4_1),m.item4_2=mean(s$item4_2),m.item4_3=mean(s$item4_3),
+                      m.item5_1=mean(s$item5_1),m.item5_2=mean(s$item5_2),m.item5_3=mean(s$item5_3),
+                      m.item6_1=mean(s$item6_1),m.item6_2=mean(s$item6_2),m.item6_3=mean(s$item6_3),
+                      m.item7_1=mean(s$item7_1),m.item7_2=mean(s$item7_2),m.item7_3=mean(s$item7_3)
+      )
+    }
+  }
+  N <- ifelse(substr(scenario,start=nchar(scenario)-1,stop=nchar(scenario))=="50","50",substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario)))
+  zz <- ifelse(N=="50",substr(scenario,start=0,stop=nchar(scenario)-3),substr(scenario,start=0,stop=nchar(scenario)-4))
+  eff.size <- unique(res.dat[res.dat$scenario==zz & res.dat$N==N,'eff.size'])
+  dif.size <- unique(res.dat[res.dat$scenario==zz & res.dat$N==N,'dif.size'])
+  b <- data.frame(scenario=zz,
+                  scenario.type=substr(zz,start=nchar(zz),stop=nchar(zz)),
+                  N=N,
+                  J=J,
+                  M=M,
+                  eff.size=eff.size,
+                  nb.dif=nb.dif,
+                  dif.size=dif.size
+  )
+  true.value.in.ci <- eff.size <= s$beta+1.96*s$se_beta & eff.size >= s$beta-1.96*s$se_beta
+  beta.same.sign.truebeta.p <- ifelse(rep(eff.size,nrow(s))==0,NA,(rep(eff.size,nrow(s))/s$beta)>0)
+  num.reject <- which((s$beta-1.96*s$se_beta)>0 | (s$beta+1.96*s$se_beta)<0)
+  z <- data.frame(m.beta=mean(s$beta),
+                  se.empirical.beta=sd(s$beta),
+                  se.analytical.beta=mean(s$se_beta),
+                  m.low.ci.beta=mean(s$beta-1.96*s$se_beta),
+                  m.high.ci.beta=mean(s$beta+1.96*s$se_beta),
+                  true.value.in.ci.p=mean(true.value.in.ci),
+                  h0.rejected.p=mean( (s$beta-1.96*s$se_beta)>0 | (s$beta+1.96*s$se_beta)<0 ),
+                  beta.same.sign.truebeta.p=mean(beta.same.sign.truebeta.p),
+                  beta.same.sign.truebeta.signif.p=mean(beta.same.sign.truebeta.p[num.reject])
+  )
+  d <- cbind(b,a,z)
+  d$prop.
+  return(d)
+}
+
+
+#### Compiled results
+
+res.dat.dif.resali <- compile_simulation2_resali('1A_100')
+
+for (x in results[seq(2,length(results))]) {
+  y <- compile_simulation2_resali(x)
+  res.dat.dif.resali <- bind_rows(res.dat.dif.resali,y)
+}
+
+res.dat.dif.resali$bias <- res.dat.dif.resali$eff.size-res.dat.dif.resali$m.beta
+
+
+
+
 ##############################################################################
 #----------------------------------------------------------------------------#
 ################################## RASCHPOWER ################################
diff --git a/RProject/resali.R b/RProject/resali.R
index 968a53c..1cf0a42 100644
--- a/RProject/resali.R
+++ b/RProject/resali.R
@@ -15,6 +15,9 @@ resali <- function(df=NULL,items=NULL,group=NULL,verbose=T) {
   dat <- df
   dat$score <- rowSums(dat[,items_n])
   nqt <- ifelse(length(unique(quantile(dat$score,seq(0,1,0.2))))==6,5,length(unique(quantile(dat$score,seq(0,1,0.2))))-1)
+  while (length(unique(quantile(dat$score,seq(0,1,1/nqt))))!=nqt+1) {
+    nqt <- nqt-1
+  }
   dat$score_q5 <- cut(dat$score,unique(quantile(dat$score,seq(0,1,1/nqt))),labels=1:nqt,include.lowest=T)
   res.anova <- rep(NA,nbitems)
   pval <- rep(NA,nbitems)