first commit

R/bpcm.R

@@ -0,0 +1,416 @@
+## File Name: bpcm.R
+## File version: 1.0
+
+
+#' Compute Bayesian Partial Credit Model (BPCM) for polytomous and dichotomous items
+#'
+#' This function computes a bayesian PCM, potentially accounting for DIF on specified items
+#'
+#' @param df data.frame containing the response data
+#' @param grp vector containing the column where an optional group membership variable is stored in df
+#' @param is.dif indicator vector containing 1 at the inded of each DIF item in df and 0 otherwise
+#' @param is.unif indicator vector containing 1 at the inded of each uniform DIF item in df and 0 otherwise
+#' @param prior prior function to be used
+#' @param diag.plots boolean indicating whether the JAGS diagnosis plots should be displayed
+#' @return A data.frame containing various model outputs
+#' @import rjags
+#' @export
+
+bpcm <- function(df=NULL,grp=NULL,is.dif=NULL,is.unif=NULL,priors=NULL,param=list(),verbose=T,diag.plots=F) {
+
+  ##### Detecting errors
+  if ( (sum(is.dif)!=0 | sum(is.unif)!=0) & any(is.null(grp))) {
+    stop('ERROR: no group variable provided, but is.dif or is.unif are not NULL')
+  }
+  if ( sum(is.dif)==0 & sum(is.unif)!=0) {
+    warning('WARNING: no DIF item specified, but is.unif is not NULL. Ignoring is.unif')
+  }
+  if ( sum(is.dif)!=0 ) {
+    if (length(is.dif)!=ncol(df)) {
+      stop('ERROR: is.dif must have one element per item in df')
+    }
+    if (length(is.unif)!=ncol(df)) {
+      stop('ERROR: is.unif must have one element per item in df')
+    }
+    if ( !("m.gamma"%in%names(priors)) ) {
+      stop('ERROR: DIF requested. Please provide mean of prior gamma distribution as priors$m.gamma')
+    }
+    if ( !("s.gamma"%in%names(priors)) ) {
+      stop('ERROR: DIF requested. Please provide SD of prior gamma distribution as priors$s.gamma')
+    }
+  }
+  if ( !any(is.null(grp)) ) {
+    if ( !("m.beta"%in%names(priors)) ) {
+      stop('ERROR: Group effect requested. Please provide mean of prior beta distribution as priors$m.beta')
+    }
+    if ( !("s.gamma"%in%names(priors)) ) {
+      stop('ERROR: Group effect requested. Please provide SD of prior beta distribution as priors$s.beta')
+    }
+    if (nrow(df)!=length(grp)) {
+      stop('ERROR: grp must be of length nrow(df)')
+    }
+  }
+  if ( !("m.delta"%in%names(priors)) ) {
+    stop('ERROR: Please provide mean of prior delta distribution as priors$m.delta')
+  }
+  if ( !("s.delta"%in%names(priors)) ) {
+    stop('ERROR: Please provide SD of prior delta distribution as priors$s.delta')
+  }
+
+
+  if (verbose) {
+    cat('\n')
+    cat("#################################################################################################\n")
+    cat("######################################### FITTING MODEL #########################################\n")
+    cat("#################################################################################################\n")
+  }
+
+  ###############################################
+  ############### NO GROUP EFFECT ###############
+  ###############################################
+  if (is.null(grp)) {
+    nam <- colnames(df)
+    nam_o <- nam
+    Y <- matrix(unlist(df),nrow=nrow(df))+1
+    n <- nrow(Y)
+    p <- ncol(Y)
+    K <- apply(Y,2,max,na.rm=TRUE)
+    m.delta <- priors$m.delta
+    s.delta <- priors$s.delta
+
+    namlist <- c(sapply(1:p, function(x) paste0("delta[",x,',',paste0(1:K[p],"]"))))
+    if (sum(is.dif)!=0) {
+      namlist <- c(namlist,sapply(1:p, function(x) paste0("gamma[",x,',',paste0(1:K[p],"]"))),"beta")
+    } else {
+      if (!is.null(grp)) {
+        namlist <- c(namlist,"beta")
+      }
+    }
+
+    data <- list(Y=Y,n=n,p=p,K=K,m.delta=m.delta,s.delta=s.delta)
+    params <- c("delta")
+
+
+    if ("n.burn"%in%names(param)) {
+      n.burn <- param$n.burn
+    } else {
+      n.burn <- 4000
+    }
+    if ("n.thin"%in%names(param)) {
+      n.thin <- param$n.thin
+    } else {
+      n.thin <- 20
+    }
+    if ("n.chains"%in%names(param)) {
+      n.chains <- param$n.chains
+    } else {
+      n.chains <- 2
+    }
+    if ("n.sim"%in%names(param)) {
+      n.sim <- param$n.sim
+    } else {
+      n.sim <- 20000
+    }
+    if ("nb.cores"%in%names(param)) {
+      nb.cores <- param$nb.cores
+    } else {
+      nb.cores <- 6
+    }
+    cl <- parallel::makeForkCluster(nb.cores)
+
+    if (verbose) {
+      cat('Initializing Chains...')
+    }
+    dclone::parJagsModel(cl,name="mod.jags",data = data,inits = NULL,file = "/home/corentin/Documents/These/Recherche/SPT/src/bpcm",n.chains=n.chains)
+    if (verbose) {
+      cat('DONE \n')
+      cat('Applying burn-in iterations...')
+    }
+    dclone::parUpdate(cl,object="mod.jags",n.burn=n.burn)
+    if (verbose) {
+      cat('DONE \n')
+      cat('Running Markov Chains...')
+    }
+    mod.samples <- dclone::parCodaSamples(cl,model="mod.jags",variable.names=params,n.iter = n.sim,thin = n.thin)
+    if (verbose) {
+      cat('DONE \n')
+    }
+    parallel::stopCluster(cl)
+    if (diag.plots) {
+      cat("Displaying traceplots...")
+      traceplot(mod.samples)
+      readline(prompt="Use the arrows to navigate between traceplots. Press [enter] to continue")
+      cat("DONE\n")
+      cat("Displaying autocorrelation plots...")
+      autocorr.plot(mod.samples,ask=F)
+      readline(prompt="Use the arrows to navigate between autocorr plots. Press [enter] to continue")
+      cat("DONE\n ")
+    }
+    res <- mod.samples[[1]]
+    res <- as.data.frame(res)[,namlist]
+    namlist2 <- unlist(c(sapply(nam,function(x) paste0(x,"_",1:K[which(nam_o==x)]-1))))
+    if (sum(is.dif)!=0) {
+      namlist2 <- c(namlist2,unlist(c(sapply(nam,function(x) paste0(x,"_",1:K[which(nam_o==x)]-1,":grp")))),"beta")
+    } else {
+      if (!is.null(grp)) {
+        namlist2 <- c(namlist2,"beta")
+      }
+    }
+    colnames(res) <- namlist2
+    res <- res[,namlist2]
+    res <- res[,apply(res,2,function(x) all(x==0))==0]
+    xsi <- apply(res,2,function(x) c(mean(x),sd(x),quantile(x,0.05),quantile(x,0.95)) )
+    rownames(xsi) <- c("post.mean","post.sd","post.90.cred.low","post.90.cred.high")
+    xsi <- round(t(xsi),4)
+    if ("beta" %in% rownames(xsi)) {
+      beta <- xsi[rownames(xsi)=="beta",]
+      xsi <- xsi[rownames(xsi)!="beta",]
+    }
+    out <- list(mcmc.res=res,
+                dif.items=nam_o[which(is.dif==1)],
+                beta=beta,
+                thresholds=xsi)
+
+
+  #####################################################
+  ############### GROUP EFFECT / NO DIF ###############
+  #####################################################
+  } else if (is.null(is.dif) | sum(is.dif)==0) {
+    nam <- colnames(df)
+    nam_o <- nam
+    Y <- matrix(unlist(df),nrow=nrow(df))+1
+    Z <- matrix(unlist(grp),nrow=length(grp))
+    n <- nrow(Y)
+    p <- ncol(Y)
+    K <- apply(Y,2,max,na.rm=TRUE)
+    m.delta <- priors$m.delta
+    s.delta <- priors$s.delta
+    m.beta <- priors$m.beta
+    s.beta <- priors$s.beta
+
+    namlist <- c(sapply(1:p, function(x) paste0("delta[",x,',',paste0(1:K[p],"]"))))
+    if (sum(is.dif)!=0) {
+      namlist <- c(namlist,sapply(1:p, function(x) paste0("gamma[",x,',',paste0(1:K[p],"]"))),"beta")
+    } else {
+      if (!is.null(grp)) {
+        namlist <- c(namlist,"beta")
+      }
+    }
+
+    data <- list(Y=Y,Z=Z,n=n,p=p,K=K,m.beta=m.beta,s.beta=s.beta,m.delta=m.delta,s.delta=s.delta)
+    params <- c("delta","beta")
+
+
+    if ("n.burn"%in%names(param)) {
+      n.burn <- param$n.burn
+    } else {
+      n.burn <- 4000
+    }
+    if ("n.thin"%in%names(param)) {
+      n.thin <- param$n.thin
+    } else {
+      n.thin <- 20
+    }
+    if ("n.chains"%in%names(param)) {
+      n.chains <- param$n.chains
+    } else {
+      n.chains <- 2
+    }
+    if ("n.sim"%in%names(param)) {
+      n.sim <- param$n.sim
+    } else {
+      n.sim <- 20000
+    }
+    if ("nb.cores"%in%names(param)) {
+      nb.cores <- param$nb.cores
+    } else {
+      nb.cores <- 6
+    }
+    cl <- parallel::makeForkCluster(nb.cores)
+
+    if (verbose) {
+      cat('Initializing Chains...')
+    }
+    dclone::parJagsModel(cl,name="mod.jags",data = data,inits = NULL,file = "/home/corentin/Documents/These/Recherche/SPT/src/bpcm_beta",n.chains=n.chains)
+    if (verbose) {
+      cat('DONE \n')
+      cat('Applying burn-in iterations...')
+    }
+    dclone::parUpdate(cl,object="mod.jags",n.burn=n.burn)
+    if (verbose) {
+      cat('DONE \n')
+      cat('Running Markov Chains...')
+    }
+    mod.samples <- dclone::parCodaSamples(cl,model="mod.jags",variable.names=params,n.iter = n.sim,thin = n.thin)
+    if (verbose) {
+      cat('DONE \n')
+    }
+    parallel::stopCluster(cl)
+    if (diag.plots) {
+      cat("Displaying traceplots...")
+      traceplot(mod.samples)
+      readline(prompt="Use the arrows to navigate between traceplots. Press [enter] to continue")
+      cat("DONE\n")
+      cat("Displaying autocorrelation plots...")
+      autocorr.plot(mod.samples,ask=F)
+      readline(prompt="Use the arrows to navigate between autocorr plots. Press [enter] to continue")
+      cat("DONE\n ")
+    }
+    res <- mod.samples[[1]]
+    res <- as.data.frame(res)[,namlist]
+    namlist2 <- unlist(c(sapply(nam,function(x) paste0(x,"_",1:K[which(nam_o==x)]-1))))
+    if (sum(is.dif)!=0) {
+      namlist2 <- c(namlist2,unlist(c(sapply(nam,function(x) paste0(x,"_",1:K[which(nam_o==x)]-1,":grp")))),"beta")
+    } else {
+      if (!is.null(grp)) {
+        namlist2 <- c(namlist2,"beta")
+      }
+    }
+    colnames(res) <- namlist2
+    res <- res[,namlist2]
+    res <- res[,apply(res,2,function(x) all(x==0))==0]
+    xsi <- apply(res,2,function(x) c(mean(x),sd(x),quantile(x,0.05),quantile(x,0.95)) )
+    rownames(xsi) <- c("post.mean","post.sd","post.90.cred.low","post.90.cred.high")
+    xsi <- round(t(xsi),4)
+    if ("beta" %in% rownames(xsi)) {
+      beta <- xsi[rownames(xsi)=="beta",]
+      xsi <- xsi[rownames(xsi)!="beta",]
+    }
+    out <- list(mcmc.res=res,
+                dif.items=nam_o[which(is.dif==1)],
+                beta=beta,
+                thresholds=xsi)
+
+  #####################################################
+  ############### GROUP EFFECT / NO DIF ###############
+  #####################################################
+  } else {
+
+
+    nam <- colnames(df)
+    nam_o <- nam
+    Y <- matrix(unlist(df),nrow=nrow(df))+1
+    Z <- matrix(unlist(grp),nrow=length(grp))
+    n <- nrow(Y)
+    p <- ncol(Y)
+
+    pnodif <- p-sum(is.dif)
+    pnodif1 <- p-sum(is.dif)+1
+    pdif <- sum(is.dif)
+    pnounif <- pnodif+pdif-sum(is.unif)
+    pnounif1 <- pnodif+pdif-sum(is.unif)+1
+
+    K <- apply(Y,2,max,na.rm=TRUE)
+    m.delta <- priors$m.delta
+    s.delta <- priors$s.delta
+    m.beta <- priors$m.beta
+    s.beta <- priors$s.beta
+
+    m.gamma <- priors$m.gamma
+    s.gamma <- priors$s.gamma
+
+    namlist <- c(sapply(1:p, function(x) paste0("delta[",x,',',paste0(1:K[p],"]"))))
+    if (sum(is.dif)!=0) {
+      namlist <- c(namlist,sapply(1:p, function(x) paste0("gamma[",x,',',paste0(1:K[p],"]"))),"beta")
+    } else {
+      if (!is.null(grp)) {
+        namlist <- c(namlist,"beta")
+      }
+    }
+
+    Y <- Y[,c(which(is.dif+is.unif==0),which(is.dif+is.unif==1),which(is.dif+is.unif==2))]
+    nam <- nam[c(which(is.dif+is.unif==0),which(is.dif+is.unif==1),which(is.dif+is.unif==2))]
+
+    data <- list(Y=Y,Z=Z,n=n,p=p,pnounif=pnounif,pnounif1=pnounif1,pdif=pdif,pnodif1=pnodif1,pnodif=pnodif,K=K,m.beta=m.beta,s.beta=s.beta,m.gamma=m.gamma,s.gamma=s.gamma,m.delta=m.delta,s.delta=s.delta,difff=as.factor(is.dif),unif=as.factor(is.unif))
+    params <- c("delta","gamma","beta")
+
+
+    if ("n.burn"%in%names(param)) {
+      n.burn <- param$n.burn
+    } else {
+      n.burn <- 4000
+    }
+    if ("n.thin"%in%names(param)) {
+      n.thin <- param$n.thin
+    } else {
+      n.thin <- 20
+    }
+    if ("n.chains"%in%names(param)) {
+      n.chains <- param$n.chains
+    } else {
+      n.chains <- 2
+    }
+    if ("n.sim"%in%names(param)) {
+      n.sim <- param$n.sim
+    } else {
+      n.sim <- 20000
+    }
+    if ("nb.cores"%in%names(param)) {
+      nb.cores <- param$nb.cores
+    } else {
+      nb.cores <- 6
+    }
+    cl <- parallel::makeForkCluster(nb.cores)
+
+
+    if (verbose) {
+      cat('Initializing Chains...')
+    }
+    dclone::parJagsModel(cl,name="mod.jags",data = data,inits = NULL,file = "/home/corentin/Documents/These/Recherche/SPT/src/bpcm_dif",n.chains=n.chains)
+    if (verbose) {
+      cat('DONE \n')
+      cat('Applying burn-in iterations...')
+    }
+    dclone::parUpdate(cl,object="mod.jags",n.burn=n.burn)
+    if (verbose) {
+      cat('DONE \n')
+      cat('Running Markov Chains...')
+    }
+    mod.samples <- dclone::parCodaSamples(cl,model="mod.jags",variable.names=params,n.iter = n.sim,thin = n.thin)
+    if (verbose) {
+      cat('DONE \n')
+    }
+    parallel::stopCluster(cl)
+    if (diag.plots) {
+      cat("Displaying traceplots...")
+      traceplot(mod.samples)
+      readline(prompt="Use the arrows to navigate between traceplots. Press [enter] to continue")
+      cat("DONE\n")
+      cat("Displaying autocorrelation plots...")
+      autocorr.plot(mod.samples,ask=F)
+      readline(prompt="Use the arrows to navigate between autocorr plots. Press [enter] to continue")
+      cat("DONE\n ")
+    }
+    res <- mod.samples[[1]]
+    res <- as.data.frame(res)[,namlist]
+    namlist2 <- unlist(c(sapply(nam,function(x) paste0(x,"_",1:K[which(nam_o==x)]-1))))
+    if (sum(is.dif)!=0) {
+      namlist2 <- c(namlist2,unlist(c(sapply(nam,function(x) paste0(x,"_",1:K[which(nam_o==x)]-1,":grp")))),"beta")
+    } else {
+      if (!is.null(grp)) {
+        namlist2 <- c(namlist2,"beta")
+      }
+    }
+    colnames(res) <- namlist2
+    res <- res[,namlist2]
+    res <- res[,apply(res,2,function(x) all(x==0))==0]
+    xsi <- apply(res,2,function(x) c(mean(x),sd(x),quantile(x,0.05),quantile(x,0.95)) )
+    rownames(xsi) <- c("post.mean","post.sd","post.90.cred.low","post.90.cred.high")
+    xsi <- round(t(xsi),4)
+    if ("beta" %in% rownames(xsi)) {
+      beta <- xsi[rownames(xsi)=="beta",]
+      if (is.null(grp)) {
+        beta <- NA
+      }
+      xsi <- xsi[rownames(xsi)!="beta",]
+    }
+    out <- list(mcmc.res=res,
+                dif.items=nam_o[which(is.dif==1)],
+                beta=beta,
+                thresholds=xsi)
+  }
+  if (is.null(dim(out$beta))) {
+    out$beta <- NA
+    }
+  return(out)
+}

R/pcm.R

@@ -0,0 +1,236 @@
+## File Name: pcm.R
+## File version: 1.0
+
+#' Compute Partial Credit Model (PCM) for polytomous and dichotomous items
+#'
+#' This function computes a frequentist PCM, potentially accounting for DIF on specified items
+#'
+#' @param df data.frame containing the data
+#' @param items vector containing the names of columns where item responses are stored in df
+#' @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
+#' @return A data.frame containing various model outputs
+#' @import vcrpart
+#' @export
+
+pcm <- function(df=NULL,items=NULL,grp=NULL,dif.items=NULL,type.dif=NULL,verbose=T,fit="ucminf") {
+  ##### Detecting errors
+
+  if (any(!(items %in% colnames(df)))) {
+    stop("ERROR: provided item name does not exist in df")
+  }
+  if (any(!(grp %in% colnames(df)))) {
+    stop("ERROR: provided group variable name does not exist in df")
+  }
+  if (any(!is.null(grp))) {
+    if (any(!(grp%in%colnames(df)))) {
+      stop("ERROR: group name does not exist in df")
+    }
+  }
+  if (!is.null(dif.items) & length(dif.items)!=length(type.dif)) {
+    stop('ERROR: type.dif is not the same length as dif.items')
+  }
+  if (!is.null(dif.items) & is.null(type.dif)) {
+    warning("WARNING: no type.dif provided, assuming non-homogeneous DIF on all items")
+  }
+  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 (fit=="ucminf") {
+      fit <- "optim"
+    }
+  }
+  ##### Analysis
+  restab.diftype <- NULL
+  se.beta <- NULL
+  beta.ci <- NULL
+  beta.p <- NULL
+  nbitems <- length(items)
+  items_o <- items
+  colnames(df)[which(colnames(df)%in%items_o)] <- paste0("item",1:nbitems)
+  items <- paste0("item",1:nbitems)
+  # If no group
+  if (is.null(grp)) {
+    if (verbose) {
+      cat('\n')
+      cat("#################################################################################################\n")
+      cat("######################################### FITTING MODEL #########################################\n")
+      cat("#################################################################################################\n")
+    }
+      grp <- NULL
+      # prepare data
+      df <- df[,c('id',items)]
+      print(df)
+      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))
+      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
+      mod <- olmm(resp ~ 0 + ce(item) + re(0|id),data=df.long,family = adjacent(link = "logit"))
+      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 {
+    grp <- df[,grp]
+    df$grp <- grp
+
+    # If group and DIF
+    if (!is.null(dif.items)) {
+      if (verbose) {
+        cat('\n')
+        cat("#################################################################################################\n")
+        cat("######################################### FITTING MODEL #########################################\n")
+        cat("#################################################################################################\n")
+      }
+      # prepare data
+      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))
+      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)
+
+      # Create 1 dif column per dif item
+      for (i in 1:length(dif.items)) {
+        df.long[,paste0("dif",i)] <- ifelse(df.long$item==dif.items[i],1,0)
+      }
+      difvar <- sapply(1:length(dif.items),function(x) paste0("dif",x))
+      difvar.unif <- difvar[type.dif==1]
+      difvar.nonunif <- difvar[type.dif==0]
+      # fit pcm
+      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)
+      mod <- olmm(formudif,data=df.long,family = adjacent(link = "logit"),control=olmm_control(fit=fit))
+      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))]))
+      difcoef.unif <- NULL
+      if (length(difvar.unif)>0) {
+        difcoef.unif <- comod[(length(comod)-length(difvar.unif)):(length(comod)-1)]
+        if (length(difvar.unif)!=1) {
+          difcoef.unif <- as.matrix(difcoef.unif)
+        } else {
+          difcoef.unif <- t(as.matrix(difcoef.unif))
+        }
+        rname <- paste0("item",dif.items[type.dif==1])
+        rownames(difcoef.unif) <- paste0("dif.",items_o[which(items%in%rname)])
+        colnames(difcoef.unif) <- "gamma"
+        difcoef.unif <- as.data.frame(difcoef.unif)
+        for (k in 1:maxmod) {
+          difcoef.unif[,paste0("gamma_",k)] <- difcoef.unif[,"gamma"]
+        }
+        difcoef.unif <- as.matrix(difcoef.unif[,2:ncol(difcoef.unif)])
+      }
+      difcoef.nonunif <- NULL
+      if (length(difvar.nonunif)>0) {
+        difcoef.nonunif <- restab[nbitems+c(1:length(difvar.nonunif)),]
+        if (length(difvar.nonunif)==1) {
+          difcoef.nonunif <- t(as.matrix(difcoef.nonunif))
+        } else {
+          difcoef.nonunif <- as.matrix(difcoef.nonunif)
+        }
+        rname <- paste0("item",dif.items[type.dif==0])
+        rownames(difcoef.nonunif) <- paste0("dif.",items_o[which(items%in%rname)])
+        colnames(difcoef.nonunif) <- paste0("gamma_",1:maxmod)
+      }
+      restab <- restab[1:nbitems,]
+      rownames(restab) <- items_o
+      colnames(restab) <- paste0("delta_",1:maxmod)
+      restab.dif <- rbind(difcoef.nonunif,difcoef.unif)
+      restab.diftype <- matrix(ifelse(type.dif==1,"HOMOGENEOUS","NON-HOMOGENEOUS"))
+      restab.diftype <- noquote(restab.diftype)
+      rownames(restab.diftype) <- rownames(restab.dif)
+      colnames(restab.diftype) <- "dif.type"
+      beta <- comod["grp"]
+      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 {
+      # If group no DIF
+      if (verbose) {
+        cat('\n')
+        cat("#################################################################################################\n")
+        cat("######################################### FITTING MODEL #########################################\n")
+        cat("#################################################################################################\n")
+      }
+        # prepare data
+        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))
+        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
+        mod <- olmm(resp ~ 0 + ge(grp) + ce(item) + re(0|id),data=df.long,family = adjacent(link = "logit"),control=olmm_control(fit=fit))
+        comod <- coef(mod)
+        # output results
+        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])
+    }
+
+  }
+  theta <- -1*ranef(mod,norm=F)+ifelse(grp==1,beta,0)
+  resid <- apply(matrix(1:nbitems,ncol=length(nbitems)),1, function(k) sapply(1:nrow(df), function(j) res_ij(theta[j],restab[k,],df[j,items[k]],beta=0)))
+  colnames(resid) <- items_o
+
+  ##### Output
+  if (verbose) {
+    cat(paste0('Number of individuals: ',nrow(df),"\n"))
+    cat(paste0('Number of items: ',length(items),"\n"))
+    cat(paste0('Item Thresholds and DIF parameters: ',"\n"))
+  }
+
+
+  out <- list(
+    beta=beta,
+    beta.se=se.beta,
+    beta.ci=beta.ci,
+    beta.p=beta.p,
+    dif.items=dif.items,
+    dif.type=restab.diftype,
+    thresholds=restab,
+    dif.param=restab.dif,
+    theta=theta,
+    residuals=resid
+    )
+  return(out)
+}
+
+
+

R/res_ij.R

@@ -0,0 +1,37 @@
+## File Name: res_ij.R
+## File version: 1.0
+
+#' Compute rasch person-item residuals
+#'
+#' This function computes person-item residuals for rasch family adjacent-categories models
+#' like the Rasch model or the PCM
+#'
+#' @param thetahat_i Person parameter estimate for person i
+#' @param delta_hat_j Item threshold parameter estimates for item j
+#' @param res Response of person i to item j
+#' @param beta Estimated latent mean difference between groups, can be ignored
+#' @return A single scalar equal to the residual for person i on item j
+#' @import vcrpart
+#' @export
+
+res_ij <- function(thetahat_i,delta_hat_j,res,beta=0){
+  var=0
+  denomin=1
+  for (i in 1:length(delta_hat_j)) {
+
+    denomin=denomin+exp(i*(thetahat_i+beta)-sum(delta_hat_j[1:i]))
+  }
+  posterior=c(1)
+  for (i in 1:length(delta_hat_j)) {
+    posterior=c(posterior,exp(i*(thetahat_i+beta)-sum(delta_hat_j[1:i])))
+  }
+
+  posterior=posterior/denomin
+
+  for(i in 0:length(delta_hat_j)){
+
+    var=var+((i-sum(posterior*0:length(delta_hat_j)))**2)*posterior[i+1]
+  }
+  return((res-sum(posterior*0:length(delta_hat_j)))/sqrt(var))
+
+}

R/residif.R

@@ -0,0 +1,117 @@
+## File Name: residif.R
+## File version: 1.0
+
+#' RESIDIF procedure for DIF detection as per Andrich and Hagquist (2015)
+#'
+#' This function detects DIF on PCM items using ANOVA of person-item residuals
+#'
+#' @param df data.frame containing the data
+#' @param items vector containing the names of columns where item responses are stored in df
+#' @param grp vector containing the name of the column where an optional group membership variable is stored in df
+#' @return A data.frame containing a column listing the detected DIF item and another listing detected DIF forms
+#' @import vcrpart
+#' @export
+
+residif <- function(df=NULL,items=NULL,grp=NULL,verbose=T) {
+  if (any(!(items %in% colnames(df)))) {
+    stop("ERROR: provided item name does not exist in df")
+  }
+  if (any(!(grp %in% colnames(df)))) {
+    stop("ERROR: provided grp variable name does not exist in df")
+  }
+  if (any(is.null(grp))) {
+    stop("ERROR: no grp variable provided")
+  }
+  if (any(is.null(items))) {
+    stop("ERROR: no items provided")
+  }
+
+  maxcat <- max(df[,items])
+  nbitems <- length(items)
+  nbitems_o <- nbitems
+
+  if (verbose) {
+    cat('\n')
+    cat("#################################################################################################\n")
+    cat("##################################### COMPUTING INITIAL PCM #####################################\n")
+    cat("#################################################################################################\n")
+  }
+  startt <- Sys.time()
+  pcm_initial <- pcm(df = df,items = items,grp = grp,verbose=F)
+  dat <- df
+  dat$score <- rowSums(dat[,items])
+  nqt <- ifelse(length(unique(quantile(dat$score,seq(0,1,0.2))))==6,5,length(unique(quantile(dat$score,seq(0,1,0.2))))-1)
+  while (length(unique(quantile(dat$score,seq(0,1,1/nqt))))!=nqt+1) {
+    nqt <- nqt-1
+  }
+  dat$score_q5 <- cut(dat$score,unique(quantile(dat$score,seq(0,1,1/nqt))),labels=1:nqt,include.lowest=T)
+  res.anova <- rep(NA,nbitems)
+  pval <- rep(NA,nbitems)
+  fval <- rep(NA,nbitems)
+  for (i in 1:nbitems) {
+    dat[,paste0('res_',i)] <- pcm_initial$residuals[,i]
+    res.anova[i] <- summary(aov(dat[,paste0('res_',i)]~get(grp)*score_q5,data=dat))
+    pval[c(i,i+nbitems)] <- c(res.anova[[i]][1,"Pr(>F)"],res.anova[[i]][3,"Pr(>F)"])
+    fval[c(i,i+nbitems)] <- c(res.anova[[i]][1,'F value'],res.anova[[i]][3,"F value"])
+  }
+  if (verbose) {
+    cat('DONE\n')
+    cat('#################################################################################################\n')
+  }
+  res.items <- c()
+  res.uniform <- c()
+  resp <- df[,items]
+  k <- 1
+  while(any(pval<0.05/(nbitems_o*3))) {
+    k <- k+1
+    if (verbose) {
+      cat(paste("######################################## COMPUTING STEP",k,"#######################################\n"))
+      cat("#################################################################################################\n")
+    }
+    numitem <- ifelse(which.max(fval)%%(length(fval)/2)!=0,which.max(fval)%%(length(fval)/2),length(fval)/2)
+    res.item <- gsub("[a-z]", "",colnames(resp)[numitem])
+    res.items <- c(res.items,res.item)
+    res.uni <- res.anova[[numitem]][3,"Pr(>F)"]>0.05
+    res.uniform <- c(res.uniform,res.uni)
+    pcm_while <- pcm(df = df,items = items,grp = grp,dif.items = res.items,type.dif = res.uniform,verbose=F)
+    res.anova <- rep(NA,nbitems)
+    pval <- rep(NA,nbitems_o)
+    fval <- rep(NA,nbitems_o)
+    numitems <- 1:nbitems_o
+    numitems <- numitems[-which(numitems%in%res.items)]
+    for (i in numitems) {
+      dat[,paste0('res_',i)] <- pcm_while$residuals[,i]
+      res.anova[i] <- summary(aov(dat[,paste0('res_',i)]~dat[,grp]*score_q5,data=dat))
+      pval[c(i,i+nbitems)] <- c(res.anova[[i]][1,"Pr(>F)"],res.anova[[i]][3,"Pr(>F)"])
+      fval[c(i,i+nbitems)] <- c(res.anova[[i]][1,'F value'],res.anova[[i]][3,"F value"])
+    }
+    for (i in 1:nbitems_o) {
+      pval[i] <- ifelse(is.na(pval[i]),999,pval[i])
+      fval[i] <- ifelse(is.na(fval[i]),-999,fval[i])
+    }
+    if (verbose) {
+      cat('DONE\n')
+      if (any(pval<0.05/(nbitems_o*3))) {
+        cat('#################################################################################################\n')
+      }
+    }
+  }
+  endt <- Sys.time()
+  cat(paste(c('Algorithm ran for',round(endt-startt,4),"seconds\n")))
+  if (verbose) {
+    cat('#################################################################################################\n')
+    cat("###################################### DETECTED DIF ITEMS #######################################\n")
+    cat("#################################################################################################\n")
+  }
+  if (length(res.items>0)) {
+    results <- data.frame(dif.items=res.items,
+                          uniform=ifelse(res.uniform==1,TRUE,FALSE))
+    return(results)
+  }
+  else {
+    if (verbose) {
+      cat("No DIF was detected\n")
+    }
+    return(NULL)
+  }
+}