R/comparePercICCsMirt.R

In joakimwallmark/PolyOptimalIRT:

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

# mirtThetas and optThetas need to correspond to the same if transformed
# binvalues is matrix the number of categories as columns
comparePercICCsMirt <- function(item, mirtModel, mirtThetas, thetasToShow=seq(0, 100, 1),
                            bincenters=NULL, binvalues = NULL,
                            surprisal=F, maxy=6, title="ICC curves",
                            grayscale=F, legend=T) {
    # Compute the probability trace lines
    extr <- extract.item(mirtModel, item)
    categories <- extr@ncat
    # Evalue ICC prob/surp. corresponding to each percentile theta value
    mirtItemICC <- probtrace(extr, mirtThetas)

    if (surprisal) {
        yText <- expression(S(theta))
        mirtItemICC <- -log(mirtItemICC, base=categories)
        maxy <- min(maxy, max(mirtItemICC))
    }
    else {
        maxy <- 1
        yText <- expression(P(theta))
    }

    itemForPlot <- data.frame(Theta = thetasToShow)
    itemForPlot1 <- cbind(itemForPlot, rep("Parametric", 101), mirtItemICC)
    colnames(itemForPlot1) <- c(colnames(itemForPlot)[1], "type", 0:(categories-1))
    # gather IRT probs
    longer.format <- tidyr::gather(itemForPlot1, category, probability,
                                    3:(3+categories-1))

    p <- ggplot(NULL) +
        geom_line(data = longer.format, aes(Theta, probability,
                                            color=category)) +
        #ggtitle(title) +
        xlab(expression("percentage of" ~ 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.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
}