Updated RESALI with new bonferroni correction

RProject/resali.R

@@ -14,7 +14,186 @@ resali <- function(df=NULL,items=NULL,group=NULL,verbose=T) {
   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)
+  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)] <- 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)
@@ -31,6 +210,7 @@ resali <- function(df=NULL,items=NULL,group=NULL,verbose=T) {
   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) {
@@ -72,7 +252,16 @@ resali <- function(df=NULL,items=NULL,group=NULL,verbose=T) {
       return(results)
     }
     else {
-    cat("No DIF was detected")
+      if (verbose) {
+        cat("No DIF was detected\n")
+      }
+      return(NULL)
+    }
+  }
+  else {
+    if (verbose) {
+      cat("No DIF was detected\n")
+    }
     return(NULL)
   }
 }