From 873020b813ce31db291c812b126f87c80591e367 Mon Sep 17 00:00:00 2001
From: corentinchoisy <corentin.choisy@proton.me>
Date: Tue, 3 Jun 2025 13:40:05 +0200
Subject: [PATCH] Added na.action option to pcm

---
 R/pcm.R    | 19 ++++++++++---------
 man/pcm.Rd |  3 ++-
 2 files changed, 12 insertions(+), 10 deletions(-)

diff --git a/R/pcm.R b/R/pcm.R
index b6c7ea2..c65e62b 100644
--- a/R/pcm.R
+++ b/R/pcm.R
@@ -20,7 +20,7 @@
 #' @export
 
 
-pcm <- function(df=NULL,items=NULL,grp=NULL,dif.items=NULL,type.dif=NULL,weights=NULL,verbose=T,fit="ucminf",method.theta="eap") {
+pcm <- function(df=NULL,items=NULL,grp=NULL,dif.items=NULL,type.dif=NULL,weights=NULL,verbose=T,fit="ucminf",method.theta="eap",na.action=na.omit) {
   ##### Detecting errors
 
   if (any(!(items %in% colnames(df)))) {
@@ -43,7 +43,7 @@ pcm <- function(df=NULL,items=NULL,grp=NULL,dif.items=NULL,type.dif=NULL,weights
   if (!("id"%in%colnames(df))) {
     stop('ERROR: no column named id provided')
   }
-  if ( any(apply(df[df[,grp]==0,items],2,max)<max(df[,items])) | any(apply(df[df[,grp]==1,items],2,max)<max(df[,items])) ) {
+  if (!is.null(grp) & any(apply(df[df[,grp]==0,items],2,max)<max(df[,items])) | any(apply(df[df[,grp]==1,items],2,max)<max(df[,items])) ) {
     if (fit=="ucminf") {
       fit <- "optim"
     }
@@ -57,6 +57,7 @@ pcm <- function(df=NULL,items=NULL,grp=NULL,dif.items=NULL,type.dif=NULL,weights
   items_o <- items
   colnames(df)[which(colnames(df)%in%items_o)] <- paste0("item",1:nbitems)
   items <- paste0("item",1:nbitems)
+  grpo <- grp
   # If no group
   if (is.null(grp)) {
     if (verbose) {
@@ -90,9 +91,9 @@ pcm <- function(df=NULL,items=NULL,grp=NULL,dif.items=NULL,type.dif=NULL,weights
     df.long$id <- factor(df.long$id)
     # fit pcm
     if (is.null(weights)) {
-      mod <- olmm(resp ~ 0 + ce(item) + re(0|id),data=df.long,family = adjacent(link = "logit"))
+      mod <- olmm(resp ~ 0 + ce(item) + re(0|id),data=df.long,family = adjacent(link = "logit"),na.action = na.action)
     } else {
-      mod <- olmm(resp ~ 0 + ce(item) + re(0|id),data=df.long,family = adjacent(link = "logit"),weights = df.long$weights)
+      mod <- olmm(resp ~ 0 + ce(item) + re(0|id),data=df.long,family = adjacent(link = "logit"),weights = df.long$weights,na.action = na.action)
     }
     comod <- coef(mod)
     # output results
@@ -149,9 +150,9 @@ pcm <- function(df=NULL,items=NULL,grp=NULL,dif.items=NULL,type.dif=NULL,weights
       formudif <- paste0("resp ~ 0 + ge(grp",ifelse(length(difvar.unif>0),"+",""),ifelse(length(difvar.unif>0),paste0(difvar.unif,":grp",collapse="+"),""),")+ce(item",ifelse(length(difvar.nonunif>0),"+",""),ifelse(length(difvar.nonunif)>0,paste0(difvar.nonunif,":grp",collapse="+"),""),")+re(0|id)")
       formudif <- as.formula(formudif)
       if (is.null(weights)) {
-        mod <- olmm(formudif,data=df.long,family = adjacent(link = "logit"),control=olmm_control(fit=fit))
+        mod <- olmm(formudif,data=df.long,family = adjacent(link = "logit"),control=olmm_control(fit=fit),na.action = na.action)
       } else {
-        mod <- olmm(formudif,data=df.long,family = adjacent(link = "logit"),control=olmm_control(fit=fit),weights = df.long$weights)
+        mod <- olmm(formudif,data=df.long,family = adjacent(link = "logit"),control=olmm_control(fit=fit),weights = df.long$weights,na.action = na.action)
       }
       comod <- coef(mod)
       # output results
@@ -241,9 +242,9 @@ pcm <- function(df=NULL,items=NULL,grp=NULL,dif.items=NULL,type.dif=NULL,weights
       df.long$id <- factor(df.long$id)
       # fit pcm
       if (is.null(weights)) {
-        mod <- olmm(resp ~ 0 + ge(grp) + ce(item) + re(0|id),data=df.long,family = adjacent(link = "logit"),control=olmm_control(fit=fit))
+        mod <- olmm(resp ~ 0 + ge(grp) + ce(item) + re(0|id),data=df.long,family = adjacent(link = "logit"),control=olmm_control(fit=fit),na.action = na.action)
       } else {
-        mod <- olmm(resp ~ 0 + ge(grp) + ce(item) + re(0|id),data=df.long,family = adjacent(link = "logit"),control=olmm_control(fit=fit),weights=df.long$weights)
+        mod <- olmm(resp ~ 0 + ge(grp) + ce(item) + re(0|id),data=df.long,family = adjacent(link = "logit"),control=olmm_control(fit=fit),weights=df.long$weights,na.action = na.action)
       }
       comod <- coef(mod)
       # output results
@@ -268,7 +269,7 @@ pcm <- function(df=NULL,items=NULL,grp=NULL,dif.items=NULL,type.dif=NULL,weights
 
   }
   if (method.theta=="eap") {
-    theta <- c(-1*ranef(mod,norm=F)+ifelse(grp==1,beta,0))
+    theta <- c(-1*ranef(mod,norm=F)+ ifelse(is.null(grpo),0, ifelse(grp==1,beta,0) ) )
   } else if (method.theta=="wle") {
     theta <- PP::PP_gpcm(as.matrix(df[,items]),t(restab),rep(1,length(items)))$resPP$resPP[,1]
   } else if (method.theta=="mle") {
diff --git a/man/pcm.Rd b/man/pcm.Rd
index 4086cf7..1aa1a0e 100644
--- a/man/pcm.Rd
+++ b/man/pcm.Rd
@@ -13,7 +13,8 @@ pcm(
   weights = NULL,
   verbose = T,
   fit = "ucminf",
-  method.theta = "eap"
+  method.theta = "eap",
+  na.action = na.omit
 )
 }
 \arguments{