tirt (version 0.2.0)

polytomous_irt: Polytomous Item Response Theory Estimation Using Likelihood or Bayesian

Description

Estimates item and person parameters for polytomous item response theory models using either Marginal Maximum Likelihood or Joint Maximum Likelihood. Now supports flexible prior distributions for Bayesian estimation (MAP estimation).

Usage

polytomous_irt(data, model = "GPCM", method = "EM", control = list())

Value

A list containing:

• item_params: Data frame of estimated parameters (a, thresholds).

• person_params: A data frame of estimated person abilities (theta) and standard errors.

• model_fit: A data frame containing fit statistics such as Akaike’s Information Criterion (AIC) and the Bayesian Information Criterion (BIC).

• settings: A list of control parameters used in the estimation.

Arguments

data

A N x J data.frame of polytomous responses (0, 1, 2...). Missing values should be NA. Categories must be continuous integers.

model

String. "GPCM" (Generalized Partial Credit Model), "PCM" (Partial Credit Model), or "GRM" (Graded Response Model).

method

String. "EM" (Marginal Maximum Likelihood via Expectation-Maximization) or "MLE" (Joint Maximum Likelihood). However, using Bayesian will override the likelihood estimation.

control

A list of control parameters for the estimation algorithm:

• max_iter: Maximum number of EM iterations (default = 100).

• converge_tol: Convergence criterion for parameter change (default = 1e-4).

• theta_range: Numeric vector of length 2 specifying the integration grid bounds (default = c(-4, 4)).

• quad_points: Number of quadrature points (default = 21).

• verbose: Logical; if TRUE, prints progress to console.

• prior: A list specifying prior distributions for item parameters. Default is NULL (no priors). For GRM: list(a = function(x) dlnorm(x, 0, 0.5, log=TRUE), d = function(x) dnorm(x, 0, 2, log=TRUE)). For GPCM: list(a = function(x) dlnorm(x, 0, 0.5, log=TRUE), d = function(x) dnorm(x, 0, 2, log=TRUE)). For PCM: list(d = function(x) dnorm(x, 0, 2, log=TRUE)). Each prior is a function returning log-density and applies to ALL items. Fixed value priors are NOT supported. Item-specific priors are NOT supported.

Examples

  # --- Example 1: Simulation (GPCM) ---
  set.seed(2026)
  N <- 500; J <- 5
  n_cats <- c(3, 4, 3, 5, 4)

  true_theta <- rnorm(N)
  true_a <- runif(J, 0.8, 1.2)
  true_d <- list()

  # Generate Thresholds
  for(j in 1:J) {
    steps <- sort(rnorm(n_cats[j]-1, mean = 0, sd = 1.0))
    true_d[[j]] <- c(0, cumsum(steps))
  }

  # Simulation Helper (GPCM Logic)
  generate_resp <- function(theta, a, d_vec, n_cat) {
    probs <- matrix(0, length(theta), n_cat)
    for(k in 1:n_cat) {
      z <- a * (k-1) * theta - d_vec[k]
      probs[,k] <- exp(z)
    }
    probs <- probs / rowSums(probs)
    apply(probs, 1, function(p) sample(0:(n_cat-1), 1, prob=p))
  }

  # Create Data
  sim_data <- matrix(NA, nrow = N, ncol = J)
  for(j in 1:J) {
    sim_data[,j] <- generate_resp(true_theta, true_a[j], true_d[[j]], n_cats[j])
  }
  df_sim <- as.data.frame(sim_data)

  # Run Estimation (GPCM to match simulation logic without prior)
  res <- polytomous_irt(df_sim, model="GPCM", method="EM",
                        control=list(max_iter=20, verbose=TRUE))

  head(res$item_params)
  print(res$model_fit)

  # \donttest{
  # Run Estimation with prior (MAP)
  res_prior <- polytomous_irt(df_sim, model="PCM", method="EM",
                              control=list(max_iter=20, verbose=FALSE,
                                          prior=list(
                                            d = function(x) dnorm(x, 0, 2, log=TRUE)
                                          )))
  head(res$item_params)
  print(res$model_fit)

  # --- Example 2: With Package Data (GRM) ---
  data("ela1", package = "tirt")

  # Subset polytomous items (columns 31 to 45)
  df_poly <- ela1[, 31:45]

  # Run Estimation using GRM
  real_res <- polytomous_irt(df_poly, model="GRM", method="EM",
                             control = list(max_iter = 1000))

  head(real_res$item_params)
  head(real_res$person_params)
  print(real_res$model_fit)
  # Run Estimation using GRM with prior
  real_res2 <- polytomous_irt(df_poly, model="GRM", method="EM",
                             control = list(max_iter = 1000,
                                           prior = list(
                                             a = function(x) dlnorm(x, 0, 0.5, log=TRUE),
                                             d = function(x) dnorm(x, 0, 2, log=TRUE)
                                           )))

  head(real_res2$item_params)
  head(real_res2$person_params)
  print(real_res2$model_fit)
  # }