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])
  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)
  fval <- rep(NA,nbitems)
  for (i in items) {
    dat[,paste0('res_',i)] <- IRT.residuals(pcm_initial)$stand_residuals[,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*3))) {
    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)] <- IRT.residuals(pcm_while)$stand_residuals[,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 {
    if (verbose) {
      cat("No DIF was detected\n")
    }
    return(NULL)
  }
}





resali_BF2 <- 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])
  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)
  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)
  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+1))) {
    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 {
    if (verbose) {
      cat("No DIF was detected\n")
    }
    return(NULL)
  }
}








resali_LRT <- 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]
  items_n_alt <- paste0(items_n,c("_TT","_noTT"))
  for (i in items_n) {
    df[df$TT==0,paste0(i,"_noTT")] <- df[df$TT==0,i]
    df[df$TT==1,paste0(i,"_TT")] <- df[df$TT==1,i]
  }
  resp_alt <- df[,items_n_alt]
  pcm_initial <- TAM::tam.mml(resp=resp,Y=grp,irtmodel = "PCM",est.variance = T,verbose=F)
  pcm_alt <- TAM::tam.mml(resp=resp_alt,Y=grp,irtmodel = "PCM",est.variance = T,verbose=F)
  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)
  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)
  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
  if (anova(pcm_initial,pcm_alt)$p[1]<0.05) {
    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 {
      if (verbose) {
        cat("No DIF was detected\n")
      }
      return(NULL)
    }
  }
  else {
    if (verbose) {
      cat("No DIF was detected\n")
    }
    return(NULL)
  }
}