From de01f2d837501f8dabe9eedef9bfb58d30ec5b94 Mon Sep 17 00:00:00 2001
From: corentinchoisy <corentin.choisy@proton.me>
Date: Tue, 3 Jun 2025 16:25:15 +0200
Subject: [PATCH] corrected theta and residuals computation whith dif in pcm

---
 R/pcm.R | 30 ++++++++++++++++++++++++------
 1 file changed, 24 insertions(+), 6 deletions(-)

diff --git a/R/pcm.R b/R/pcm.R
index 9a0680d..88674cc 100644
--- a/R/pcm.R
+++ b/R/pcm.R
@@ -45,11 +45,11 @@ pcm <- function(df=NULL,items=NULL,grp=NULL,dif.items=NULL,type.dif=NULL,weights
     stop('ERROR: no column named id provided')
   }
   suppressWarnings(
-  if (!is.null(grp) & any(apply(df[df[,grp]==0,items],2,function(k) max(k,na.rm = T))<max(df[,items],na.rm = T)) | any(apply(df[df[,grp]==1,items],2,function(k) max(k,na.rm = T))<max(df[,items],na.rm=T)) ) {
-    if (fit=="ucminf") {
-      fit <- "optim"
+    if (!is.null(grp) & any(apply(df[df[,grp]==0,items],2,function(k) max(k,na.rm = T))<max(df[,items],na.rm = T)) | any(apply(df[df[,grp]==1,items],2,function(k) max(k,na.rm = T))<max(df[,items],na.rm=T)) ) {
+      if (fit=="ucminf") {
+        fit <- "optim"
+      }
     }
-  }
   )
   ##### Analysis
   restab.diftype <- NULL
@@ -271,6 +271,8 @@ pcm <- function(df=NULL,items=NULL,grp=NULL,dif.items=NULL,type.dif=NULL,weights
     }
 
   }
+
+  if (is.null(dif.items)) {
     if (method.theta=="eap") {
       theta <- c(-1*ranef(mod,norm=F)+ ifelse(is.null(grpo),0, ifelse(grp==1,beta,0) ) )
     } else if (method.theta=="wle") {
@@ -278,9 +280,25 @@ pcm <- function(df=NULL,items=NULL,grp=NULL,dif.items=NULL,type.dif=NULL,weights
     } else if (method.theta=="mle") {
       theta <- PP::PP_gpcm(as.matrix(df[,items]),t(restab),rep(1,length(items)),type="mle")$resPP$resPP[,1]
     }
+    resid <- apply(matrix(1:nbitems,ncol=length(nbitems)),1, function(k) sapply(1:nrow(df), function(j) resi(theta[j],restab[k,],df[j,items[k]],beta=0)))
+  } else {
+    restaba <- restab
+    restaba[rownames(restab)[dif.items],] <- restaba[rownames(restab)[dif.items],] + restab.dif[paste0("dif.",rownames(restab)[dif.items]),]
+    theta <- rep(NA,nrow(df))
+    resid <- matrix(NA,nrow=nrow(df),ncol=nbitems)
+    if (method.theta=="eap") {
+      theta <- c(-1*ranef(mod,norm=F)+ ifelse(is.null(grpo),0, ifelse(grp==1,beta,0) ) )
+    } else if (method.theta=="wle") {
+      theta[grp==0] <- PP::PP_gpcm(as.matrix(df[grp==0,items]),t(restab),rep(1,length(items)))$resPP$resPP[,1]
+      theta[grp==1] <- PP::PP_gpcm(as.matrix(df[grp==1,items]),t(restaba),rep(1,length(items)))$resPP$resPP[,1]
+    } else if (method.theta=="mle") {
+      theta[grp==0] <- PP::PP_gpcm(as.matrix(df[grp==0,items]),t(restab),rep(1,length(items)),type="mle")$resPP$resPP[,1]
+      theta[grp==1] <- PP::PP_gpcm(as.matrix(df[grp==1,items]),t(restaba),rep(1,length(items)),type="mle")$resPP$resPP[,1]
+    }
+    resid[grp==0,] <- apply(matrix(1:nbitems,ncol=length(nbitems)),1, function(k) sapply(which(grp==0), function(j) resi(theta[j],restab[k,],df[j,items[k]],beta=0)))
+    resid[grp==1,] <- apply(matrix(1:nbitems,ncol=length(nbitems)),1, function(k) sapply(which(grp==1), function(j) resi(theta[j],restaba[k,],df[j,items[k]],beta=0)))
+  }
 
-
-  resid <- apply(matrix(1:nbitems,ncol=length(nbitems)),1, function(k) sapply(1:nrow(df), function(j) resi(theta[j],restab[k,],df[j,items[k]],beta=0)))
   colnames(resid) <- items_o
 
   ##### Output