added rosali R script

RProject/Scripts/rosali.R

@@ -0,0 +1,31 @@
+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)
+}