R/irt_item_fit.R

In Pflegermeister/wizirt: A United API for IRT

#' A function used to get additional item-level information, such as ICC probabilities, and item information.
#'
#' @param wizirt_fit An object coming from the fit_wizirt function.
#' @param stats A character or character string identifying item-level fit measures. Must be at least one of c('Zh', 'X2', 'G2', 'infit'). More are coming very soon.  Default is 'X2'.
#' @return A list with item_stats (item-level statistics), item_information (data for item and test information), and item_probabilities (data for ICC curves).
#' @examples
#' ifa <- wizirt2:::irt_item_fit(my_model)
#' @export
irt_item_fit <- function(wizirt_fit, stats = c("X2")) {
  if (!wizirt_fit$fit$model$item_type %in% c("Rasch", "1PL", "2PL")) {
    rlang::abort(glue::glue("Cannot calculate item-fit statistics for models of type: {wizirt_fit$fit$model$item_type}"))
  }

  out <- list(
    item_stats = NULL,
    item_information = NULL,
    item_probabilities = NULL,
    spec = list(stats = stats)
  )

  # Item_stats ---------------------------------------------------------------------------------------
  # If anyone tries to use groups this won't work. xxxx

  mod <- to_mirt(wizirt_fit)

  item_stats <- mirt::itemfit(mod,
    fit_stats = stats,
    Theta = matrix(wizirt_fit$fit$parameters$persons$ability)
  )

  item_stats <- dplyr::left_join(wizirt_fit$fit$parameters$coefficients,
    item_stats,
    by = "item"
  )

  # item_information ---------------------------------------------------------------------------------


  item_information <- data.frame(item = names(wizirt_fit$fit$data), data.frame(t(sapply(1:ncol(wizirt_fit$fit$data), function(i) {
    item <- mirt::extract.item(mod, i)
    mirt::iteminfo(item, Theta = seq(-6, 6, length.out = 100))
  }))))
  item_information <- item_information %>% `colnames<-`(c("item", seq(-6, 6, length.out = 100)))
  item_information <- item_information %>%
    tidyr::pivot_longer(cols = -1, names_to = "theta", values_to = "info") %>%
    dplyr::mutate(theta = as.numeric(theta))

  # item_probabilities ---------------------------------------------------------------------------------
  plt_list <- list()

  for (i in 1:ncol(wizirt_fit$fit$data)) {
    plt <- mirt::itemplot(mod, i)
    pltdata <- data.frame(lapply(plt$panel.args, function(x) do.call(cbind, x))[[1]])
    y <- pltdata$y
    names(y) <- pltdata$x
    plt_list[[i]] <- y
  }

  plt_data <- tibble::as_tibble(do.call(rbind, plt_list))
  plt_data$item <- colnames(wizirt_fit$fit$data)
  item_probabilities <- plt_data %>%
    tidyr::pivot_longer(cols = -item, names_to = "x", values_to = "y") %>%
    dplyr::mutate(x = as.numeric(x))
  # joining ------------------------------------------------------------------------------------------

  out$item_stats <- item_stats
  out$item_information <- item_information
  out$item_probabilities <- item_probabilities

  class(out) <- c("wizirt_ifa", class(out))
  out
}

to_slope <- function(b, a, g = 0, u = 1) {
  d <- -a * b
  tibble::tibble(a1 = a, d = d, g = g, u = u)
}

to_mirt <- function(wizirt_fit) {
  fd <- `colnames<-`(
    cbind(
      matrix(0,
        ncol = ncol(wizirt_fit$fit$data),
        nrow = nrow(wizirt_fit$fit$data)
      ),
      matrix(1,
        ncol = ncol(wizirt_fit$fit$data),
        nrow = nrow(wizirt_fit$fit$data)
      )
    ),
    c(
      paste(colnames(wizirt_fit$fit$data), 1, sep = "_"),
      paste(colnames(wizirt_fit$fit$data), 2, sep = "_")
    )
  )

  fd <- fd[, as.vector(rbind(1:ncol(wizirt_fit$fit$data), (ncol(wizirt_fit$fit$data) + 1):ncol(fd)))]
  fd <- `colnames<-`(matrix(fd, ncol = ncol(fd)), colnames(fd))

  Call <- wizirt_fit$fit$model$engine$call %>% `class<-`("call")

  Data <- list(
    data = data.frame(wizirt_fit$fit$data),
    group = factor(rep("all", nrow(wizirt_fit$fit$data)), levels = "all", labels = 1),
    groupNames = "all",
    ngroups = 1,
    mins = rep(0, ncol(wizirt_fit$fit$data)),
    fulldata = list(fd),
    K = rep(2L, ncol(wizirt_fit$fit$data)),
    nitems = ncol(wizirt_fit$fit$data),
    model = list(tabdata = matrix(wizirt_fit$fit$data))
  )

  Model <- list(
    itemloc = seq(from = 1, to = (ncol(wizirt_fit$fit$data) * 2 + 1), by = 2),
    nfact = 1,
    Theta = wizirt_fit$fit$parameters$persons$ability
  )

  mirt_coefs <- to_slope(
    b = wizirt_fit$fit$parameters$coefficients$difficulty,
    a = wizirt_fit$fit$parameters$coefficients$discrimination,
    g = wizirt_fit$fit$parameters$coefficients$guessing
  )

  a1 <- mirt_coefs$a1
  d <- mirt_coefs$d
  g <- mirt_coefs$g
  u <- mirt_coefs$u

  est <- `names<-`(dplyr::case_when(
    wizirt_fit$fit$model$item_type == "Rasch" ~ c(F, T, F, F),
    T ~ c(T, T, F, F)
  ), c("a1", "d", "g", "u"))

  pars_den <- function(obj, Theta) {
    gpars <- ExtractGroupPars(obj)
    mu <- gpars$gmeans
    sigma <- gpars$gcov
    d <- mirt_dmvnorm(Theta, mean = mu, sigma = sigma)
    d <- ifelse(d < 1e-300, 1e-300, d)
    d
  }
  pars <- list()
  for (i in 1:ncol(wizirt_fit$fit$data)) {
    pars[[i]] <- new("dich",
      par = c(a1[[i]], d[[i]], g[[i]], u[[i]]),
      est = est,
      nfact = 1L,
      itemclass = 1L,
      ncat = length(unique(wizirt_fit$fit$data[[i]]))
    )
  }
  pars[[i + 1]] <- new("GroupPars",
    density = 0,
    dentype = "Gaussian",
    den = pars_den,
    nfact = 1L,
    itemclass = 0L
  )

  # I need to put this somewhere for future review.
  # Here I have really hacked up the mirt classes to make them fit.
  # If I can make the translation more elegant than I can go back and forth easily.
  # However, this is not elegant at all.

  ParObjects <- list(pars = pars)

  Internals <- list(CUSTIM.IND = 0)

  Fit <- list(F = matrix(0, ncol = wizirt_fit$fit$model$n_factors, nrow = nrow(wizirt_fit$fit$parameters$persons)))

  Options <- list(exploratory = FALSE, quadpts = mirt:::select_quadpts(nfact = wizirt_fit$fit$model$n_factors))

  mod <- new("SingleGroupClass",
    Call = Call,
    Data = Data,
    Model = Model,
    Fit = Fit,
    ParObjects = ParObjects,
    Options = Options,
    Internals = Internals
  )
  return(mod)
}