R/plotICCcombo.R

In joakimwallmark/PolyOptimalIRT:

library(ggplot2)
library(tidyr)
library(mirt)

# thetas needs to be in optimal score scale
# binvalues is matrix the number of categories as columns
plotICCcombo <- function(item, mirtModel, WfdList, thetas,
                         bincenters=NULL, binvalues = NULL,
                         surprisal=F, infscale=F, maxy=6, title="ICC curves",
                         grayscale=F, legend=T) {
    # Compute the probability trace lines and information for all items
    # Put into lists
    categories <- ncol(WfdList[[item]]$Wmatfine)
    thetainf <- scaleThetaToInf(min(thetas), max(thetas), theta = thetas)
    extr <- extract.item(mirtModel, item)
    mirtItemICC <- probtrace(extr, thetainf)

    if (surprisal) {
        yText <- expression(S(theta))
        optimalItemICC <- WfdList[[item]]$Wmatfine
        maxy <- min(maxy, max(WfdList[[item]]$Wmatfine))
        mirtItemICC <- -log(mirtItemICC, base=categories)
    }
    else {
        maxy <- 1
        yText <- expression(P(theta))
        optimalItemICC <- WfdList[[item]]$Pmatfine
    }

    if (infscale) {
        title <- expression(paste("ICC curves ", theta, ":(-", infinity,
                         ", ", infinity, ")"))
        bincenters <- scaleThetaToInf(0, max(thetas), theta = bincenters)
        thetas <- thetainf
    }
    else {
        thetas <- thetas
    }

    itemForPlot <- data.frame(Theta = thetas)
    itemForPlot1 <- cbind(itemForPlot, rep("Parametric", 101), mirtItemICC)
    itemForPlot2 <- cbind(itemForPlot, rep("Optimal", 101), optimalItemICC)
    colnames(itemForPlot1) <- colnames(itemForPlot2) <- c(colnames(itemForPlot)[1],
                                                          "type", 0:(categories-1))
    # gather IRT probs
    longer.format1 <- tidyr::gather(itemForPlot1, category, probability,
                                    3:(3+categories-1))
    # gather optimal score probs
    longer.format2 <- tidyr::gather(itemForPlot2, category, probability,
                                    3:(3+categories-1))
    longer.format <- rbind(longer.format1, longer.format2)


    p <- ggplot(NULL) +
        geom_line(data = longer.format, aes(Theta, probability,
                                     color=category,
                                     linetype=type)) +
        #ggtitle(title) +
        xlab(expression(theta)) +
        ylab(yText) +
        scale_x_continuous(expand = c(0, 0)) +
        scale_y_continuous(limits = c(NA, maxy), expand = c(0, 0)) +
        #geom_hline(aes(yintercept = yinterc)) +
        theme_bw() +
        theme(plot.margin = unit(c(0, 0, 1, 0), "cm"))+
        theme(plot.margin = unit(c(0.2, 0.4, 0.2, 0.2), "cm"),
              text = element_text(size=16),
              axis.text.x=element_text(colour="black"),
              axis.text.y=element_text(colour="black"),
              legend.title=element_blank())
    if (grayscale)
        p <- p + scale_colour_grey()
    if (!legend)
        p <- p + theme(legend.position = "none")
    if (!is.null(bincenters)&!is.null(binvalues)) {
        binDF <- as.data.frame(cbind(bincenters, binvalues))
        colnames(binDF)[2:(categories+1)] <- seq(0, categories-1)
        binDF <- tidyr::gather(binDF, category, binval, 2:(categories+1))
        p <- p + geom_point(data = binDF,
                            aes(x = bincenters, y = binval, color=category))
    }
    p
}