R/pfa_mlm.R

In Pflegermeister/wizirt2: A United API for IRT in R

#' A function to explore reasons for person misfit.
#' @param pfa_object An object of class wizirt_pfa (i.e. from irt_person_fit)
#' @param wizirt_fit An object of class wizirt_irt (i.e. from fit_wizirt)
#' @param bins For method == 'Conijn' only. Either a number indicating how many bins to break up items into,
#' or a vector of proportions for as many groups as are desired.  If a vector, must sum to 1.
#' @param predictors Either a matrix or dataframe with a row per person and a column per predictor.
#' @param method Either 'Reise' (the default for Rasch) or 'Conijn' (2PL default)
#' @param rownames A vector of unique IDs to identify the rows (i.i person IDs).
#' @export
pfa_mlm <- function(pfa_object, wizirt_fit, rownames = NULL, predictors = NULL, method = NULL, bins = 5){
  if (is.null(method)) {
    message(glue::glue('Defaulting to Conijn method with {bins} bins.'))
    method <- 'Conijn'
  }
  if (!method %in% c('Reise', 'Conijn')){
    rlang::abort(glue::glue("Uknown method '{method}'."))
  }

  if(is.null(rownames)){
  rownames <- 1:nrow(wizirt_fit$fit$data)
  }


  mlm_data <- data.frame(cbind(ids = rownames, wizirt_fit$fit$parameters$persons, wizirt_fit$fit$data)) %>%
    tidyr::pivot_longer(cols = -1:-3, names_to = 'item') %>%
    dplyr::left_join(wizirt_fit$fit$parameters$coefficients, by = 'item')
  if(!is.null(predictors)){
    mlm_data <- dplyr::left_join(mlm_data,
                                 data.frame(cbind(ids = rownames,
                                                  predictors)), by = 'ids')
  }
  mlm_data <- mlm_data %>% tibble::as_tibble() %>% dplyr::mutate(prob = exp(discrimination*(ability-difficulty))/
                                 (1 + exp(discrimination*(ability-difficulty))) )

  if(method == 'Reise'){

    mod1 <- mlm_data %>%
      blme::bglmer(value ~
                     (1 + difficulty|ids) +
                     difficulty,
                   data = ., family = binomial)

    mod2 <- mlm_data %>%
      blme::bglmer(value ~
                     (1 + difficulty|ids) +
                     ability + difficulty,
                   data = ., family = binomial)

    mod3 <- mlm_data %>%
      blme::bglmer(value ~
                     (0 + difficulty + ability|ids) +
                     ability + difficulty,
                   data = ., family = binomial)

    out <- list(models = list(level_1 = mod1,
                              level_2 = mod2,
                              level_3 = mod3)) %>%
      `class<-`(c('list', 'pfa_mlm'))

    return(out)
  } else if (method == 'Conijn') {
    mlm_data <- mlm_data %>% dplyr::left_join(get_bins(wizirt_fit = wizirt_fit, bins = bins), by = 'item')
    out <- lapply(1:bins, pfa_fit_subset, data = mlm_data, stats = pfa_object$spec$stats)
    mlm_data <- mlm_data %>% dplyr::left_join(do.call(rbind.data.frame, out), by = c('ids', 'item_bin' = 'bin')) %>% dplyr::select(-std_err:-prob) %>% dplyr::distinct()
    mod_list <- list()
    for (i in stats){
      mod_list[[i]] <- eval(parse(text = glue::glue('blme::blmer(',
                                                    '{i} ~',
                                                    '(1|ids),',
                                                    'data = mlm_data)')))
    }


    lapply(mod_list, plot)
    out <- list(models = mod_list, icc = sapply(mod_list, performance::icc)) %>% `class<-`(c('list', 'pfa_mlm'))
    return(out)
  }
}


pfa_fit_subset <- function(bin, data = mlm_data, stats = pfa_object$spec$stats){

  mlm_sub <- data %>%
    dplyr::filter(item_bin == bin) %>%
    tidyr::pivot_wider(id_cols = ids, names_from = item, values_from = value) %>%
    dplyr::select(-ids)

  mlm_coefs <- data %>%
    dplyr::filter(item %in% colnames(mlm_sub)) %>%
    dplyr::distinct(item, .keep_all = T) %>%
    dplyr::select(difficulty, discrimination)

  stats_list <- list()
  for (i in stats) {
    fit <- eval(parse(text = glue::glue('PerFit::{i}(mlm_sub,',
                                        'IP = cbind(mlm_coefs, guessing = 0),',
                                        'Ability = mlm_data %>% dplyr::distinct(ids, ability) %>% dplyr::pull(ability)',
                                        ')')))
    stats_list[[i]] <- fit$PFscores$PFscores
  }

  tibble::as_tibble(stats_list) %>%
    dplyr::mutate(ids = rownames, bin = bin)

}

get_bins <- function(wizirt_fit, bins){
  r <- wizirt_fit$fit$parameters$coefficients %>% dplyr::arrange(difficulty) %>% nrow()
  remainder <- rep(1:0, times = c((r%%bins), bins-(r%%bins)))
  data.frame(item = wizirt_fit$fit$parameters$coefficients$item,
             item_bin = rep(1:bins, times = (rep(r%/%bins, times = bins) + remainder)))
}
# I am going to make a new function for the multilevel of Reise and Conijn
# mlm = "Conijn",
# mlm_covariates = NULL
# Reise will be the default when Rasch is used, otherwise Conijn

# At level 1, the person’s response to an item is a function
# of a person intercept and a person slope (Walker, Engelhard, 2015).

#
#
# mlm_data %>% dplyr::mutate(prob = exp(discrimination*(theta-difficulty))/
#                              (1 + exp(discrimination*(theta-difficulty))) ) %>%
#   ggplot2::ggplot(aes(x = difficulty, y = prob)) +
#   geom_point() +
#   facet_wrap(~unique_id)
#
# mod1 <- mlm_data %>%
#   blme::bglmer(response ~
#                 (1 + difficulty|unique_id) +
#                 difficulty,
#               data = ., family = binomial)
#
#
# mod1
# coef(mod1)
# performance::icc(mod1)
#
# mlm_coefs <- coef(mod1)$unique_id %>%
#   as.data.frame() %>%
#   tibble::rownames_to_column()
#
# mod2 <- mlm_data %>%
#   blme::bglmer(response ~
#                 (1 + difficulty|unique_id) +
#                 theta + difficulty,
#               data = ., family = binomial)
# mod2
# coef(mod2)
#
#
# # conijn, to implement at a later time
#
#
# # step 1: break the test up into different subsets
# # step 2: calculate a fit stat for each.
# # step 3: predict the fit per person using person
#
#