Added na.action option to pcm

This commit is contained in:
2025-06-03 13:40:05 +02:00
parent 0ce205a6d4
commit 873020b813
2 changed files with 12 additions and 10 deletions

19
R/pcm.R
View File

@ -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") {