demo/toolbox2.R

In psychotools: Psychometric Modeling Infrastructure

library(strucchange)

set.seed(2906)

data("VerbalAggression", package = "psychotools")

refmodel <- nplmodel(VerbalAggression$resp2[, 1:6])

dgp <- function(model, N = 1000, G = 1, impact = FALSE, cotype = "random") {

  mu <- if(impact) rep_len(c(-0.5, 0.5), N) else rep_len(0, N)
  theta <- rnorm(N, mu, 1)

  covariate <- switch(as.character(cotype),
    "random" = rnorm(N, 0, 1),
    "linear" = theta + rnorm(N, 0, 1),
    "quadratic" = theta ^ 2 + rnorm(N, 0, 1)
  )

  d <- data.frame(theta = theta, covariate = covariate)
  if(G > 1) {
    d$impact <- cut(covariate, labels = 1:G, include.lowest = TRUE,
      breaks = quantile(covariate, probs = 0:G / G))
  }

  itempar_dif <- itempar(model)
  itempar_dif[1] <- itempar_dif[1] + sd(itempar_dif)

  dif_id <- covariate > median(covariate)

  d$resp <- matrix(NA, N, length(model$items))
  d$resp[!dif_id, ] <-
    rpl(theta[!dif_id],
      a = discrpar(model),
      b = itempar(model),
      g = guesspar(model),
      u = upperpar(model),
      return_setting = FALSE)
  d$resp[dif_id, ] <-
    rpl(theta[dif_id],
      a = discrpar(model),
      b = itempar_dif,
      g = guesspar(model),
      u = upperpar(model),
      return_setting = FALSE)

  return(d)
}

hitrate <- function(model, M = 1000, alpha = 0.05, parm = NULL, N = 1000,
  G = 1, impact = FALSE, cotype = "random") {

  pval <- replicate(M, {
    d <- dgp(model, N = N, G = G, impact = impact, cotype = cotype)
    m <- nplmodel(d$resp, type = "2PL", impact = d$impact,
      maxit = 5000, reltol = 1e-4, vcov = FALSE)
    sctest(m, order.by = d$covariate, functional = "DM", parm = parm)$p.value
  })
  mean(pval < alpha)
}

sim <- function(model, M = 1000, alpha = 0.05, parm = NULL, N = 1000,
  G = c(1, 2, 5, 25), impact = c(FALSE, TRUE),
  cotype = c("random", "linear", "quadratic")) {

  d <- expand.grid(G = G, impact = impact, cotype = cotype)
  d$hitrate <- NA
  for(i in seq_len(NROW(d))) {
    d$hitrate[i] <- hitrate(model, M = M, alpha = alpha, parm = parm, N = N,
    G = d$G[i], impact = d$impact[i], cotype = d$cotype[i])
  }
  return(d)
}

simulation2 <- sim(refmodel)