From d071500c446e6d63dcc4ef22e275bbee6561c8e9 Mon Sep 17 00:00:00 2001 From: corentinchoisy Date: Mon, 26 May 2025 16:00:29 +0200 Subject: [PATCH] Correted pcm weights --- R/pcm.R | 168 +++++++++++++++++++++++++++++------------------------ man/pcm.Rd | 2 +- 2 files changed, 94 insertions(+), 76 deletions(-) diff --git a/R/pcm.R b/R/pcm.R index 2f184a6..07b286c 100644 --- a/R/pcm.R +++ b/R/pcm.R @@ -1,5 +1,5 @@ ## File Name: pcm.R -## File version: 1.0 +## File version: 1.1 #' Compute Partial Credit Model (PCM) for polytomous and dichotomous items #' @@ -10,7 +10,7 @@ #' @param grp string containing the name of the column where an optional group membership variable is stored in df #' @param dif.items vector containing the list of indexes in "items" corresponding to dif items #' @param type.dif vector containing DIF form for each item specified in dif.items. 1 is homogeneous DIF, 0 is heterogeneous DIF -#' @param weights string containing the name of the column where optional weights are stored in df +#' @param weights string containing the name of the column where an optional variable containing weights is stored in df #' @param verbose set to TRUE to print a detailed output, FALSE otherwise #' @param fit string determining the optimization algorithm. Values "ucminf" or "nlminb" ar recommended #' @param method.theta string determining the estimation method for individual latent variable values. Either "eap", "mle" or "wle" @@ -19,6 +19,7 @@ #' @import PP #' @export + pcm <- function(df=NULL,items=NULL,grp=NULL,dif.items=NULL,type.dif=NULL,weights=NULL,verbose=T,fit="ucminf",method.theta="eap") { ##### Detecting errors @@ -64,39 +65,42 @@ pcm <- function(df=NULL,items=NULL,grp=NULL,dif.items=NULL,type.dif=NULL,weights cat("######################################### FITTING MODEL #########################################\n") cat("#################################################################################################\n") } - grp <- NULL - # prepare data - if (is.null(weights)) { - df <- df[,c('id',items)] - } else { - df <- df[,c('id',items,weights)] - } - print(df) + grp <- NULL + # prepare data + if (is.null(weights)) { + df <- df[,c('id',items)] colnames(df)[2:(length(colnames(df)))] <- paste0("item",seq(1,length(colnames(df))-1)) - df.long <- reshape(df,v.names=c("item"),direction="long",varying=c(items)) - if (is.null(weights)) { - colnames(df.long) <- c("id","item","resp") - } else { - colnames(df.long) <- c("id","item","resp","weights") - } + } else { + df <- df[,c('id',items,weights)] + colnames(df)[2:(length(colnames(df)-1))] <- paste0("item",seq(1,length(colnames(df))-1)) + } + df.long <- reshape(df,v.names=c("item"),direction="long",varying=c(items)) + if (is.null(weights)) { + colnames(df.long) <- c("id","item","resp") nbitems <- length(2:(length(colnames(df)))) maxmod <- max(df[,2:(length(colnames(df)))]) df.long$item <- factor(df.long$item,levels=seq(1,length(colnames(df))-1),ordered = F) - df.long$resp <- factor(df.long$resp,0:maxmod,ordered=T) - 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")) - } else { - mod <- olmm(resp ~ 0 + ce(item) + re(0|id),data=df.long,family = adjacent(link = "logit"),weights = df.long$weights) - } - comod <- coef(mod) - # output results - restab <- t(sapply(1:nbitems,function(x) comod[seq(x,length(comod)-1,nbitems)])) - rownames(restab) <- paste0("item",1:nbitems) - colnames(restab) <- paste0("delta_",1:maxmod) - restab.dif <- NULL - beta <- NULL + } else { + colnames(df.long) <- c("id","weights","item","resp") + nbitems <- length(2:(length(colnames(df))-1)) + maxmod <- max(df[,2:(length(colnames(df))-1)]) + df.long$item <- factor(df.long$item,levels=seq(1,length(colnames(df))-2),ordered = F) + } + df.long$resp <- factor(df.long$resp,0:maxmod,ordered=T) + 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")) + } else { + mod <- olmm(resp ~ 0 + ce(item) + re(0|id),data=df.long,family = adjacent(link = "logit"),weights = df.long$weights) + } + comod <- coef(mod) + # output results + restab <- t(sapply(1:nbitems,function(x) comod[seq(x,length(comod)-1,nbitems)])) + rownames(restab) <- paste0("item",1:nbitems) + colnames(restab) <- paste0("delta_",1:maxmod) + restab.dif <- NULL + beta <- NULL } # If group else { @@ -114,19 +118,23 @@ pcm <- function(df=NULL,items=NULL,grp=NULL,dif.items=NULL,type.dif=NULL,weights # prepare data if (is.null(weights)) { df <- df[,c('id',items,"grp")] + colnames(df)[2:(length(colnames(df))-1)] <- paste0("item",seq(1,length(colnames(df))-2)) } else { df <- df[,c('id',items,"grp",weights)] + colnames(df)[2:(length(colnames(df))-2)] <- paste0("item",seq(1,length(colnames(df))-3)) } - colnames(df)[2:(length(colnames(df))-1)] <- paste0("item",seq(1,length(colnames(df))-2)) df.long <- reshape(df,v.names=c("item"),direction="long",varying=c(items)) if (is.null(weights)) { colnames(df.long) <- c("id","grp","item","resp") + nbitems <- length(2:(length(colnames(df))-1)) + maxmod <- max(df[,2:(length(colnames(df))-1)]) + df.long$item <- factor(df.long$item,levels=seq(1,length(colnames(df))-2),ordered = F) } else { - colnames(df.long) <- c("id","grp","item","resp","weights") + colnames(df.long) <- c("id","grp","weights","item","resp") + nbitems <- length(2:(length(colnames(df))-2)) + maxmod <- max(df[,2:(length(colnames(df))-2)]) + df.long$item <- factor(df.long$item,levels=seq(1,length(colnames(df))-3),ordered = F) } - nbitems <- length(2:(length(colnames(df))-1)) - maxmod <- max(df[,2:(length(colnames(df))-1)]) - df.long$item <- factor(df.long$item,levels=seq(1,length(colnames(df))-2),ordered = F) df.long$resp <- factor(df.long$resp,0:maxmod,ordered=T) df.long$id <- factor(df.long$id) @@ -148,7 +156,11 @@ pcm <- function(df=NULL,items=NULL,grp=NULL,dif.items=NULL,type.dif=NULL,weights comod <- coef(mod) # output results nbcoef <- nbitems+length(difvar.nonunif) - restab <- t(sapply(1:nbcoef,function(x) comod[seq(x,length(comod)-2-length(difvar.unif),nbitems+length(difvar.nonunif))])) + if (is.null(weights)) { + restab <- t(sapply(1:nbcoef,function(x) comod[seq(x,length(comod)-2-length(difvar.unif),nbitems+length(difvar.nonunif))])) + } else { + restab <- t(sapply(1:nbcoef,function(x) comod[seq(x,length(comod)-2-length(difvar.unif),nbitems+length(difvar.nonunif))])) + } difcoef.unif <- NULL if (length(difvar.unif)>0) { difcoef.unif <- comod[(length(comod)-length(difvar.unif)):(length(comod)-1)] @@ -204,45 +216,54 @@ pcm <- function(df=NULL,items=NULL,grp=NULL,dif.items=NULL,type.dif=NULL,weights cat("######################################### FITTING MODEL #########################################\n") cat("#################################################################################################\n") } - # prepare data - if (is.null(weights)) { - df <- df[,c('id',items,"grp")] - } else { - df <- df[,c('id',items,"grp",weights)] - } + # prepare data + if (is.null(weights)) { + df <- df[,c('id',items,"grp")] colnames(df)[2:(length(colnames(df))-1)] <- paste0("item",seq(1,length(colnames(df))-2)) - df.long <- reshape(df,v.names=c("item"),direction="long",varying=c(items)) - if (is.null(weights)) { - colnames(df.long) <- c("id","grp","item","resp") - } else { - colnames(df.long) <- c("id","grp","item","resp","weights") - } + + } else { + df <- df[,c('id',items,"grp",weights)] + colnames(df)[2:(length(colnames(df))-2)] <- paste0("item",seq(1,length(colnames(df))-3)) + } + df.long <- reshape(df,v.names=c("item"),direction="long",varying=c(items)) + if (is.null(weights)) { + colnames(df.long) <- c("id","grp","item","resp") nbitems <- length(2:(length(colnames(df))-1)) maxmod <- max(df[,2:(length(colnames(df))-1)]) df.long$item <- factor(df.long$item,levels=seq(1,length(colnames(df))-2),ordered = F) - df.long$resp <- factor(df.long$resp,0:maxmod,ordered=T) - 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)) - } 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) - } - comod <- coef(mod) - # output results + } else { + colnames(df.long) <- c("id","grp","weights","item","resp") + nbitems <- length(2:(length(colnames(df))-2)) + maxmod <- max(df[,2:(length(colnames(df))-2)]) + df.long$item <- factor(df.long$item,levels=seq(1,length(colnames(df))-3),ordered = F) + } + df.long$resp <- factor(df.long$resp,0:maxmod,ordered=T) + 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)) + } 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) + } + comod <- coef(mod) + # output results + if (is.null(weights)) { restab <- t(sapply(1:nbitems,function(x) comod[seq(x,length(comod)-2,nbitems)])) - rownames(restab) <- items_o - colnames(restab) <- paste0("delta_",1:maxmod) - restab.dif <- NULL - beta <- comod[length(comod)-1] - se.beta <- (confint(mod)["grp",2]-beta)/1.96 - beta.ci <- confint(mod)["grp",] - beta.p <- 2*pnorm(-abs(beta/se.beta)) - beta <- as.numeric(beta) - se.beta <- as.numeric(se.beta) - beta.p <- as.numeric(beta.p) - beta <- -1*beta - beta.ci <- -1*c(beta.ci[2],beta.ci[1]) + } else { + restab <- t(sapply(1:nbitems,function(x) comod[seq(x,length(comod)-2,nbitems)])) + } + rownames(restab) <- items_o + colnames(restab) <- paste0("delta_",1:maxmod) + restab.dif <- NULL + beta <- comod[length(comod)-1] + se.beta <- (confint(mod)["grp",2]-beta)/1.96 + beta.ci <- confint(mod)["grp",] + beta.p <- 2*pnorm(-abs(beta/se.beta)) + beta <- as.numeric(beta) + se.beta <- as.numeric(se.beta) + beta.p <- as.numeric(beta.p) + beta <- -1*beta + beta.ci <- -1*c(beta.ci[2],beta.ci[1]) } } @@ -275,9 +296,6 @@ pcm <- function(df=NULL,items=NULL,grp=NULL,dif.items=NULL,type.dif=NULL,weights dif.param=restab.dif, theta=theta, residuals=resid - ) + ) return(out) } - - - diff --git a/man/pcm.Rd b/man/pcm.Rd index 46b892b..4086cf7 100644 --- a/man/pcm.Rd +++ b/man/pcm.Rd @@ -27,7 +27,7 @@ pcm( \item{type.dif}{vector containing DIF form for each item specified in dif.items. 1 is homogeneous DIF, 0 is heterogeneous DIF} -\item{weights}{string containing the name of the column where optional weights are stored in df} +\item{weights}{string containing the name of the column where an optional variable containing weights is stored in df} \item{verbose}{set to TRUE to print a detailed output, FALSE otherwise}