Package {tirt}

Binary (Dichotomous) Item Response Theory Estimation Using Likelihood or Bayesian

Description

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

Usage

binary_irt(data, model = "2PL", method = "EM", control = list())

Arguments

Value

A list containing:

• item_params: A data frame of estimated item parameters (discrimination, difficulty, guessing) and their standard errors.

• 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), the Bayesian Information Criterion (BIC), and Log-Likelihood.

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

Examples

  # # Simulate data
  set.seed(123)
  N <- 500; J <- 10
  true_theta <- rnorm(N)
  true_b <- seq(-2, 2, length.out=J)
  true_a <- runif(J, 0.8, 1.2)
  data_mat <- matrix(NA, N, J)
  for(i in 1:N) {
    p <- 1 / (1 + exp(-true_a * (true_theta[i] - true_b)))
    data_mat[i,] <- rbinom(J, 1, p)
  }
  df <- as.data.frame(data_mat)
  names(df) <- paste0("Q", 1:J)

  # # Run Function without prior
  res <- binary_irt(df, model="2PL", method="EM")


  # # Run Function with prior (function-based)
  res_prior <- binary_irt(df, model="2PL", method="EM",
                          control=list(prior=list(
                            a = function(x) dlnorm(x, 0, 0.5, log=TRUE),
                            b = function(x) dnorm(x, 0, 2, log=TRUE)
                          )))

  # # Run Function with fixed value prior
  res_fixed <- binary_irt(df, model="2PL", method="EM",
                          control=list(prior=list(
                            a = 1,  # Fix all discrimination prior to 1
                            b = function(x) dnorm(x, 0, 2, log=TRUE)
                          )))

  # View Results
  head(res$item_params)
  head(res$person_params)
  print(res$model_fit)
  # --- Example 2: With Package Data ---
  data("ela1", package = "tirt")
  # Subset the first 30 columns (must use the object name 'data_binary')
  df <- ela1[, 1:30]
  # Run Function on package data
  real_res <- binary_irt(df, model="2PL", method="EM")
  head(real_res$item_params)
  

Differential Item Functioning (Mantel-Haenszel and Logistic Regression)

Description

Detects Differential Item Functioning (DIF) for dichotomous items, that is, items that behave differently for two groups of examinees (for example, a reference and a focal group) after matching on ability. Two widely used approaches are provided: the Mantel-Haenszel (1959) procedure, with the ETS delta effect size and A/B/C flagging scheme (Holland & Thayer, 1988), and the logistic regression approach (Swaminathan & Rogers, 1990), which separates uniform and non-uniform DIF.

Usage

dif(data, group, focal = NULL, match = NULL, purify = FALSE)

Arguments

Value

A data frame with one row per item and the columns:

• item: the item name.

• MH_chisq, MH_p: the Mantel-Haenszel chi-square (with continuity correction) and its p-value.

• MH_OR: the Mantel-Haenszel common odds ratio.

• MH_delta: the ETS delta effect size, -2.35 \ln(OR).

• ETS_class: the ETS DIF classification, "A" (negligible), "B" (moderate), or "C" (large).

• LR_chisq, LR_p: the 2-degree-of-freedom logistic regression test for combined (uniform + non-uniform) DIF.

• LR_uniform_p: p-value for uniform DIF (group main effect).

• LR_nonuniform_p: p-value for non-uniform DIF (group-by-match interaction).

References

Holland, P. W., & Thayer, D. T. (1988). Differential item performance and the Mantel-Haenszel procedure. In H. Wainer & H. I. Braun (Eds.), Test validity (pp. 129-145). Erlbaum.

Swaminathan, H., & Rogers, H. J. (1990). Detecting differential item functioning using logistic regression procedures. Journal of Educational Measurement, 27(4), 361-370.

Examples

  set.seed(123)
  sim <- sim_irt(n_people = 600,
                 item_structure = list(list(model = "2PL", n_items = 10)))
  resp <- sim$resp

  # Create two groups and plant DIF in item 3 (harder for group "B")
  grp <- rep(c("A", "B"), each = 300)
  flip <- grp == "B" & resp[[3]] == 1
  resp[[3]][flip] <- rbinom(sum(flip), 1, 0.6)

  dif_res <- dif(resp, group = grp)
  dif_res

Mixed-Format English Language Arts (ELA) Assessment Data (Form 1)

Description

A dataset containing binary and polytomous responses for demonstration.

Usage

ela1

Format

A data frame with 52417 rows and 47 columns:

• ITEM1 - ITEM30: Binary responses (0 = Incorrect, 1 = Correct).

• ITEM31 - ITEM45: Polytomous responses (scored 0-5).

• THETA: Latent ability estimates.

• COVARIATE: Person-level background variable.

Source

Tang, C., Xiong, J., & Engelhard, G. (2025). Identification of writing strategies in educational assessments with an unsupervised learning measurement framework. Education Sciences, 15(7), 912. doi:10.3390/educsci15070912

Examples

data(ela1)
head(ela1)

Mixed-Format English Language Arts (ELA) Assessment Data (Form 2)

Description

A smaller dataset containing item responses.

Usage

ela2

Format

A data frame with columns representing item responses.

• ITEM1 - ITEM7: Binary responses (0 = Incorrect, 1 = Correct).

• ITEM8: Polytomous response (scored 0-2).

• ITEM9: Polytomous response (scored 0-5).

• ITEM10: Polytomous response (scored 0-5).

Source

Tang, C., Xiong, J., & Engelhard, G. (2025). Identification of writing strategies in educational assessments with an unsupervised learning measurement framework. Education Sciences, 15(7), 912. doi:10.3390/educsci15070912

Examples

data(ela2)
head(ela2)

Large-Scale Mixed-Format English Language Arts (ELA) Assessment Data (Form 3)

Description

A long format dataset containing binary and polytomous responses for demonstration.

Usage

ela3

Format

A data frame with 2,434,185 observations and 5 variables:

• STUDENTID: Unique student identifier.

• FORMID: Unique test form identifier.

• ITEMID: Unique item identifier, where _T1 indicates trait 1 and _T2 indicates trait 2.

• SEQ: Position of the item in a form.

• SCORE: Dichotomously and polytomously scored item responses (0,1,2...).

Source

Tang, C., Xiong, J., & Engelhard, G. (2025). Identification of writing strategies in educational assessments with an unsupervised learning measurement framework. Education Sciences, 15(7), 912. doi:10.3390/educsci15070912

Examples

data(ela3)
head(ela3)

Large-Scale Mixed-Format English Language Arts (ELA) Assessment Data Testmap (Form 3 Testmap)

Description

Test map for the long format ela3 response data

Usage

ela3_testmap

Format

A data frame with 328 rows and 5 variables:

• FORMID: Unique test form identifier.

• SEQ: Position of the item in a form.

• ITEMID: Unique item identifier, where _T1 indicates trait 1 and _T2 indicates trait 2.

• TYPE: Type of an item, where OP indicates operational items and FT indicates field-test items.

• MAX_SCORE: Max score of that item.

Source

Tang, C., Xiong, J., & Engelhard, G. (2025). Identification of writing strategies in educational assessments with an unsupervised learning measurement framework. Education Sciences, 15(7), 912. doi:10.3390/educsci15070912

Examples

data(ela3_testmap)
head(ela3_testmap)

Item Response Theory Equating / Linking

Description

Conducts item response theory scale linking using Mean-Mean, Mean-Sigma, and Stocking-Lord methods. Supports mixed formats of both dichotomous and polytomous models. Automatically detects anchor items and validates model consistency.

Usage

equate_irt(base_params, new_params, person_params = NULL, methods = NULL)

Arguments

Value

A list containing three data frames:

Examples

  # ===========================================================================
  # Example: Equating Form Y (New) to Form X (Base)
  # ===========================================================================
  set.seed(123)

  # 1. Generate "True" Base Parameters (Form X)
  # ---------------------------------------------------------------------------
  # 10 Common Items (Anchors) + 10 Unique Items
  # 2PL and GRM mixed

  gen_item_params <- function(n, type="2PL") {
    if(type=="2PL") {
      data.frame(
        item = paste0("Item_", 1:n),
        model = "2PL",
        a = round(runif(n, 0.8, 1.5), 2),
        b = round(rnorm(n, 0, 1), 2),
        stringsAsFactors = FALSE
      )
    } else {
      # GRM with 3 thresholds
      d <- t(apply(matrix(rnorm(n*3, 0, 0.5), n, 3), 1, sort))
      df <- data.frame(
        item = paste0("Poly_", 1:n),
        model = "GRM",
        a = round(runif(n, 0.8, 1.5), 2),
        stringsAsFactors = FALSE
      )
      df <- cbind(df, setNames(as.data.frame(d), paste0("step_", 1:3)))
      df
    }
  }

  # Anchors
  anchor_2pl <- gen_item_params(5, "2PL")
  anchor_grm <- gen_item_params(3, "GRM")
  # Unique Form X
  unique_x <- gen_item_params(5, "2PL")
  unique_x$item <- paste0("X_", unique_x$item)

  base_params <- dplyr::bind_rows(anchor_2pl, anchor_grm, unique_x)

  # 2. Generate "New" Form Y Parameters (with Scale Shift)
  # ---------------------------------------------------------------------------
  # Scale Transformation: Theta_base = 1.2 * Theta_new + 0.5
  # True Constants: A = 1.2, B = 0.5
  TRUE_A <- 1.2
  TRUE_B <- 0.5

  # Transform Anchor Parameters to "New" scale (Inverse Logic)
  # a_new = a_base * A
  # b_new = (b_base - B) / A

  anchor_2pl_new <- anchor_2pl
  anchor_2pl_new$a <- anchor_2pl$a * TRUE_A
  anchor_2pl_new$b <- (anchor_2pl$b - TRUE_B) / TRUE_A

  anchor_grm_new <- anchor_grm
  anchor_grm_new$a <- anchor_grm$a * TRUE_A
  step_cols <- grep("step_", names(anchor_grm_new))
  anchor_grm_new[, step_cols] <- (anchor_grm[, step_cols] - TRUE_B) / TRUE_A

  # Unique Form Y
  unique_y <- gen_item_params(5, "2PL")
  unique_y$item <- paste0("Y_", unique_y$item)

  new_params <- dplyr::bind_rows(anchor_2pl_new, anchor_grm_new, unique_y)

  # 3. Create Dummy Person Parameters for Form Y
  # ---------------------------------------------------------------------------
  person_params <- data.frame(
    id = paste0("P", 1:50),
    theta = rnorm(50, 0, 1),
    theta_se = runif(50, 0.2, 0.5)
  )

  # 4. Perform Equating
  # ---------------------------------------------------------------------------
  # We expect to recover A approx 1.2 and B approx 0.5
  results <- equate_irt(
    base_params = base_params,
    new_params = new_params,
    person_params = person_params,
    methods = c("Mean-Mean", "Stocking-Lord")
  )

  # 5. Inspect Results
  # ---------------------------------------------------------------------------
  # Linking Constants
  print(results$linking_constants)

  # Transformed Items (Form Y items on Form X scale)
  head(results$transformed_item_params)

  # Transformed Persons
  head(results$transformed_person_params)

Fixed Person Calibration with or without Covariate

Description

Estimates item parameters (difficulty, discrimination) given fixed person parameters (theta), with an optional person-level covariate. Supports Rasch and 2-Parameter Logistic models.

Usage

fix_person(df, theta, model = c("Rasch", "2PL"), covariate = NULL)

Arguments

Value

A data frame containing:

• Item statistics (difficulty, standard errors, z-values, p-values).

• Discrimination parameters (for 2PL model).

• Global covariate effect (if covariate is provided).

• Classical item statistics (p-value, count, point-biserial correlation).

• Mean theta per item (average ability of persons answering the item).

• Infit and Outfit statistics (for Rasch model only).

Examples

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

  # Subset data for a manageable example
  # Select the first 500 examinees and 30 item responses
  df_real <- ela1[1:500, 1:30]

  # Extract pre-estimated latent traits and covariates
  fixed_theta <- ela1$THETA[1:500]
  fixed_cov <- ela1$COVARIATE[1:500]

  # Estimate item parameters given fixed ability levels
  # fitting a 2-parameter logistic (2PL) model
  real_res <- fix_person(df = df_real,
                                   theta = fixed_theta,
                                   model = "2PL",
                                   covariate = fixed_cov)
  head(real_res)

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

  # Select Item Responses (Cols 1-30)
  df_real <- ela1[, 1:30]

  fixed_theta <- ela1$THETA
  fixed_cov <- ela1$COVARIATE

  real_res <- fix_person(df = df_real,
                                   theta = fixed_theta,
                                   model = "2PL",
                                   covariate = fixed_cov)
  head(real_res)

Fixed Item Calibration

Description

Estimates unknown item parameters using Marginal Maximum Likelihood via Expectation-Maximization Algorithm. Uses a custom Bounded Newton-Raphson solver. Supports mixed-format data containing dichotomous and polytomous responses

Usage

fixed_item(response_df, item_params_df, control = list())

Arguments

Value

A list containing:

• item_params: Estimated parameters for unknown items.

• person_params: Estimated person parameters.

• model_fit: A data frame containing number off estimated parameters and fit statistics such as Akaike’s Information Criterion (AIC), the Bayesian Information Criterion (BIC), and Log-Likelihood.

Examples

# 1. TOY EXAMPLE
# ===========================================================================
set.seed(123)
# Create a very small dataset (N=50, J=4)
N_toy <- 50
df_toy <- data.frame(
  I1 = rbinom(N_toy, 1, 0.5), I2 = rbinom(N_toy, 1, 0.6), # Known items
  U1 = rbinom(N_toy, 1, 0.5), U2 = rbinom(N_toy, 1, 0.4)  # Unknown items
)

# Define the "Known" parameters for I1 and I2
known_params <- data.frame(
  item = c("I1", "I2"),
  model = c("2PL", "2PL"),
  a = c(1.0, 1.2),
  b = c(-0.5, 0.5)
)

# Run Fixed Item Calibration with very low iterations
fit_toy <- fixed_item(df_toy, known_params, control=list(max_iter=2, verbose=FALSE))
print(head(fit_toy$item_params))
  
  # --- Example 2: Simulation ---
  set.seed(123)
  N <- 500
  true_theta <- rnorm(N, 0, 1)

  # 1. Simulation Helpers
  sim_2pl <- function(theta, a, b) {
    p <- 1 / (1 + exp(-1.7 * a * (theta - b)))
    rbinom(N, 1, p)
  }
  sim_poly <- function(theta, a, steps) {
    n_cat <- length(steps) + 1
    probs <- matrix(0, length(theta), n_cat)
    for(k in 1:n_cat) {
      score <- k - 1
      if(score == 0) num <- 0
      else num <- a * (score * theta - sum(steps[1:score]))
      probs[, k] <- exp(num)
    }
    probs <- probs / rowSums(probs)
    apply(probs, 1, function(x) sample(0:(n_cat-1), 1, prob=x))
  }

  # 2. Generate Data (Mixed Known/Unknown Items)
  # Items 1-5: Known Binary (2PL)
  # Items 6-10: Unknown Binary (2PL)
  # Items 11-12: Known Poly (GPCM)
  # Items 13-15: Unknown Poly (GPCM)

  resp_mat <- matrix(NA, N, 15)
  colnames(resp_mat) <- paste0("Item_", 1:15)

  # Known Binary Parameters
  a_bin <- c(1.0, 1.2, 0.9, 1.1, 0.8)
  b_bin <- c(-1, -0.5, 0, 0.5, 1)

  for(i in 1:5) resp_mat[,i] <- sim_2pl(true_theta, a_bin[i], b_bin[i])
  for(i in 6:10) resp_mat[,i] <- sim_2pl(true_theta, runif(1,0.8,1.2), rnorm(1))

  # Known Poly Parameters
  a_poly <- c(1.0, 0.8)
  d_poly <- list(c(-1, 1), c(-0.5, 0.5))

  resp_mat[,11] <- sim_poly(true_theta, a_poly[1], d_poly[[1]])
  resp_mat[,12] <- sim_poly(true_theta, a_poly[2], d_poly[[2]])
  for(i in 13:15) resp_mat[,i] <- sim_poly(true_theta, 1.0, c(-0.5, 0.5))

  df_resp <- as.data.frame(resp_mat)

  # 3. Create 'Known Parameters' Dataframe
  # This tells the function: "Fix these, Estimate the rest"
  known_df <- data.frame(
    item = c(paste0("Item_", 1:5), "Item_11", "Item_12"),
    model = c(rep("2PL", 5), rep("GPCM", 2)),
    a = c(a_bin, a_poly),
    b = c(b_bin, NA, NA),     # Binary difficulty
    step_1 = c(rep(NA, 5), -1, -0.5), # Poly steps
    step_2 = c(rep(NA, 5),  1,  0.5),
    stringsAsFactors = FALSE
  )

  # 4. Run Estimation
  res <- fixed_item(df_resp, known_df, control=list(max_iter=20))

  # View Results
  # Notice Items 1-5 and 11-12 have Status "Fixed"
  head(res$item_params, 12)

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

  # Let's treat the first 5 items as "Known" with arbitrary parameters
  # just to demonstrate syntax.
  df_real <- ela1[, 1:20]

  known_real <- data.frame(
    item = paste0("Q", 1:5),
    model = "2PL",
    a = 1.0,
    b = seq(-1, 1, length.out=5)
  )

  # Ideally, column names in df_real should match 'item' column in known_real
  colnames(df_real)[1:5] <- paste0("Q", 1:5)

  real_res <- fixed_item(df_real, known_real, control=list(max_iter=10))
  head(real_res$item_params)
  

Joint Item Response Theory and Testlet Response Theory Estimation (Dichotomous & Polytomous)

Description

Provides a unified marginal maximum likelihood estimation framework for a broad class of item response theory and testlet response theory models. The function automatically detects data structures to apply appropriate models, along with their testlet-effect extensions (Bradlow et al., 1999).

Usage

irt_trt(data, item_spec, method = "EM", control = list())

Arguments

Details

The estimation utilizes a robust Newton-Raphson update within the M-step. For testlet models, dimension reduction is achieved through the integration of the nuisance testlet effect (Li et al., 2006). The function automatically corrects model specifications if the data levels (binary vs. polytomous) do not align with the requested model string.

Value

A list containing three components:

References

Bradlow, E. T., Wainer, H., & Wang, X. (1999). A testlet response model for multidimensionality in item response theory. Psychometrika, 64(2), 147-168.

Li, Y., Bolt, D. M., & Fu, J. (2006). A comparison of methods for estimating secondary dimensions in testlet-based data. Applied Psychological Measurement, 30(3), 203-223.

Examples

  # --- Quick Example (small data) ---
  set.seed(1)
  n <- 50
  resp <- data.frame(
    I1 = rbinom(n, 1, 0.6),
    I2 = rbinom(n, 1, 0.5),
    I3 = rbinom(n, 1, 0.4),
    I4 = sample(0:2, n, replace = TRUE),
    I5 = rbinom(n, 1, 0.5),
    I6 = rbinom(n, 1, 0.6)
  )
  spec <- data.frame(
    item = colnames(resp),
    model = c("2PL", "2PL", "2PL", "GRM", "2PLT", "2PLT"),
    testlet = c(NA, NA, NA, NA, "T1", "T1"),
    stringsAsFactors = FALSE
  )
  res <- irt_trt(resp, spec, method = "EM",
                 control = list(max_iter = 5, verbose = FALSE))
  head(res$item_params)

  
  # --- Full Example: Simulation (Binary + Poly + Testlets) ---
  set.seed(2025)
  N <- 100; J <- 20

  # 1. Generate Parameters
  theta <- rnorm(N, 0, 1)
  gamma_1 <- rnorm(N, 0, 0.5) # Testlet 1 effect
  gamma_2 <- rnorm(N, 0, 0.6) # Testlet 2 effect

  a_true <- runif(J, 0.8, 1.5)
  b_true <- seq(-1.5, 1.5, length.out = J)

  resp_matrix <- matrix(NA, N, J)
  colnames(resp_matrix) <- paste0("Item_", 1:J)

  # 2. Simulate Responses
  # Items 1-10: Binary Independent (Model: 2PL)
  for(j in 1:10) {
    p <- 1 / (1 + exp(-a_true[j] * (theta - b_true[j])))
    resp_matrix[,j] <- rbinom(N, 1, p)
  }

  # Items 11-15: Poly Independent (Model: GRM)
  for(j in 11:15) {
    thresh <- sort(c(b_true[j] - 0.7, b_true[j] + 0.7))
    p1 <- 1 / (1 + exp(-a_true[j] * (theta - thresh[1])))
    p2 <- 1 / (1 + exp(-a_true[j] * (theta - thresh[2])))
    probs <- cbind(1-p1, p1-p2, p2)
    resp_matrix[,j] <- apply(probs, 1, function(p) sample(0:2, 1, prob=p))
  }

  # Items 16-17: Binary Testlet 1 (Model: 2PLT)
  for(j in 16:17) {
    eff_theta <- theta + gamma_1
    p <- 1 / (1 + exp(-a_true[j] * (eff_theta - b_true[j])))
    resp_matrix[,j] <- rbinom(N, 1, p)
  }

  # Items 18-20: Poly Testlet 2 (Model: GRT)
  for(j in 18:20) {
    eff_theta <- theta + gamma_2
    thresh <- sort(c(b_true[j] - 0.5, b_true[j] + 0.5))
    p1 <- 1 / (1 + exp(-a_true[j] * (eff_theta - thresh[1])))
    p2 <- 1 / (1 + exp(-a_true[j] * (eff_theta - thresh[2])))
    probs <- cbind(1-p1, p1-p2, p2)
    resp_matrix[,j] <- apply(probs, 1, function(p) sample(0:2, 1, prob=p))
  }

  df_sim <- as.data.frame(resp_matrix)

  # 3. Create Item Specification
  # STRICT naming: Independent=2PL/GRM, Testlet=2PLT/GRT
  spec <- data.frame(
    item = colnames(df_sim),
    model = c(rep("2PL", 10), rep("GRM", 5), rep("2PLT", 2), rep("GRT", 3)),
    testlet = c(rep(NA, 15), rep("T1", 2), rep("T2", 3)),
    stringsAsFactors = FALSE
  )

  # 4. Run Estimation
  res <- irt_trt(df_sim, spec, method = "EM",
                 control = list(max_iter = 20, verbose = FALSE))

  head(res$item_params)
  head(res$person_params)
  

Item Fit Statistics (Infit and Outfit)

Description

Computes item-level fit statistics based on standardized response residuals: the information-weighted mean square (infit) and the unweighted mean square (outfit), together with their standardized (ZSTD) transformations. These statistics flag items whose observed responses depart from the fitted model, either because of noise in the tails of the ability distribution (outfit) or near the item location (infit). The statistics are computed for dichotomous and polytomous items alike.

Usage

item_fit(
  data,
  item_params,
  theta = NULL,
  theta_range = c(-4, 4),
  n_points = 41,
  model = NULL,
  D = 1
)

Arguments

Value

A data frame with one row per item and the columns:

• item: item name.

• n: number of responses used.

• outfit: outfit mean square (unweighted).

• outfit_z: standardized outfit (ZSTD).

• infit: infit mean square (information-weighted).

• infit_z: standardized infit (ZSTD).

Mean-square values near 1 indicate good fit; values above about 1.3-1.5 suggest underfit (noise), and values below about 0.7 suggest overfit (dependency or redundancy).

References

Wright, B. D., & Masters, G. N. (1982). Rating scale analysis. MESA Press.

Examples

  set.seed(123)
  sim <- sim_irt(n_people = 400,
                 item_structure = list(list(model = "2PL", n_items = 8)))
  fit <- binary_irt(sim$resp, model = "2PL", method = "EM",
                    control = list(max_iter = 15, verbose = FALSE))

  item_fit(sim$resp, fit$item_params)

Item Information Function

Description

Computes the Fisher item information function for every item at a set of ability (theta) points. The item information function shows how precisely an item measures ability across the latent trait continuum: high information at a given theta means the item discriminates well among examinees located there. It is the building block of the test information function (test_info) and of the conditional standard error of measurement.

Usage

item_info(item_params, theta = seq(-4, 4, by = 0.1), model = NULL, D = 1)

Arguments

Details

For dichotomous models the information is I(\theta) = (Da)^2 P(\theta)(1 - P(\theta)) for the Rasch and 2PL models, and the standard Birnbaum form I(\theta) = (Da)^2 (1 - c)(1 - P^*)P^{*2} / P for the 3PL model, where P^* is the 2PL probability and P the full 3PL probability (this reduces to the 2PL expression when the guessing parameter c = 0). For polytomous models the general expected-information formula I(\theta) = \sum_k (\partial P_k / \partial \theta)^2 / P_k is used, which for the GPCM and PCM equals (Da)^2 \mathrm{Var}(X \mid \theta).

Value

A numeric matrix with one row per item and one column per theta point. Row names are the item names and column names are the theta values.

See Also

test_info, binary_irt, polytomous_irt

Examples

  # --- Example 1: dichotomous items ---
  set.seed(123)
  sim <- sim_irt(n_people = 400,
                 item_structure = list(list(model = "2PL", n_items = 6)))
  fit <- binary_irt(sim$resp, model = "2PL", method = "EM",
                    control = list(max_iter = 15, verbose = FALSE))

  # Information for each item at a grid of abilities
  info <- item_info(fit$item_params, theta = seq(-3, 3, by = 0.5))
  round(info, 3)

  
  # --- Example 2: polytomous (GRM) items ---
  simp <- sim_irt(n_people = 400,
                  item_structure = list(list(model = "GRM", n_items = 5,
                                             categories = 4)))
  fitp <- polytomous_irt(simp$resp, model = "GRM", method = "EM",
                         control = list(max_iter = 15, verbose = FALSE))
  item_info(fitp$item_params, theta = c(-2, -1, 0, 1, 2))
  

Local Dependence Statistics (Yen's Q3)

Description

Computes Yen's (1984) Q_3 statistic for every item pair. Q_3 is the correlation between the model residuals of two items after the latent trait has been partialled out. Under the local-independence assumption of unidimensional item response theory these residual correlations should be small and negative (around -1/(J-1)); large positive values signal that two items share something beyond the common trait, such as a passage or a scenario. Q_3 is therefore used both to decide whether testlet modeling is needed and to check whether a testlet model has adequately absorbed the residual dependence.

Usage

ld_stats(
  data,
  item_params,
  theta = NULL,
  theta_range = c(-4, 4),
  n_points = 41,
  model = NULL,
  D = 1
)

Arguments

Value

A symmetric numeric matrix of Q_3 statistics with one row and one column per item (diagonal set to 1). The off-diagonal entries are the residual correlations for each item pair. The average off-diagonal Q_3 is attached as the attribute "mean_q3" and the largest absolute value as "max_abs_q3".

References

Yen, W. M. (1984). Effects of local item dependence on the fit and equating performance of the three-parameter logistic model. Applied Psychological Measurement, 8(2), 125-145.

Examples

  set.seed(123)
  sim <- sim_irt(n_people = 400,
                 item_structure = list(list(model = "2PL", n_items = 8)))
  fit <- binary_irt(sim$resp, model = "2PL", method = "EM",
                    control = list(max_iter = 15, verbose = FALSE))

  q3 <- ld_stats(sim$resp, fit$item_params)
  round(q3, 3)
  attr(q3, "max_abs_q3")

Multidimensional Binary Item Response Theory Estimation

Description

Estimates item and person parameters for multidimensional binary item response models (M-Rasch, M-2PL, M-3PL). The function is fully self-contained, using custom Newton-Raphson optimization without relying on external optimizers.

Usage

mirt_binary(
  data,
  model = "2PL",
  dimension = 2,
  method = "MML",
  control = list()
)

Arguments

Value

A list containing:

• item_params: A data frame of estimated item parameters (a_1...a_D, d, g) and standard errors.

• person_params: A data frame of estimated multidimensional person abilities (if ability=TRUE).

• model_fit: Fit statistics (Log-Likelihood, AIC, BIC).

• settings: Control parameters used in the run.

Examples

  # --- Simulation for Multidimensional Data ---
  set.seed(202)
  N <- 800
  J <- 20
  D <- 2

  # Simulate Abilities (2 Dimensions, correlated)
  Sigma <- matrix(c(1, 0.4, 0.4, 1), 2, 2)
  Z <- matrix(rnorm(N * D), N, D)
  theta <- Z %*% chol(Sigma)

  # Define Q-matrix (Items 1-10 on Dim1, Items 11-20 on Dim 2)
  Q <- matrix(0, J, D)
  Q[1:10, 1] <- 1
  Q[11:20, 2] <- 1

  # Simulate Item Parameters (Multidimensional 2PL)
  a_true <- matrix(runif(J * D, 0.8, 1.8), J, D) * Q
  d_true <- seq(-2, 2, length.out = J)

  # Generate Responses
  data_mat <- matrix(NA, N, J)
  for(i in 1:N) {
    # Matrix multiply JxD slopes by Dx1 person abilities, add Jx1 intercepts
    z <- as.vector(a_true %*% theta[i, ]) + d_true
    p <- 1 / (1 + exp(-z))
    data_mat[i, ] <- rbinom(J, 1, p)
  }
  df <- as.data.frame(data_mat)
  names(df) <- paste0("Item", 1:J)

  # --- Run Custom Multidimensional Function ---
  res <- mirt_binary(df, model="2PL", dimension=2, method="MML",
                     control=list(Q_matrix=Q, ability=TRUE, max_iter=10))

  print(head(res$item_params))
  print(head(res$person_params))
  print(res$model_fit)
  
  # --- Validation with 'mirt' (For comparison, if installed) ---
  # library(mirt)
  # mirt_model <- paste0("F1 = 1-10\nF2 = 11-20\nCOV = F1*F2")
  # mirt_res <- mirt(df, mirt.model(mirt_model), itemtype='2PL', method='EM')
  # coef(mirt_res, IRTpars=FALSE, simplify=TRUE)$items
  

Mixed Item Response Model Estimation (Dichotomous & Polytomous) with Prior Support

Description

Provides a estimation framework for a broad class of different item response theory models. This function can model different combinations of item categories. Now supports flexible prior distributions for Bayesian estimation (MAP estimation).

Usage

mixed_irt(data, model = "2PL", method = "EM", control = list())

Arguments

Value

A list containing:

• item_params: Data frame of item parameters (discrimination, difficulty/thresholds, guessing).

• 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.

Examples

  # --- Example 1: Simulation (Mixed 2PL + GPCM) ---
  set.seed(2025)
  N <- 500
  n_bin <- 5
  n_poly <- 2
  J <- n_bin + n_poly

  # 1. Generate Theta (Wide range to match user request)
  true_theta <- rnorm(N, mean = 0, sd = 3)

  # 2. Simulation Helper: GPCM
  sim_gpcm <- function(theta, a, steps) {
    n_cat <- length(steps) + 1
    probs <- matrix(0, length(theta), n_cat)
    for(k in 1:n_cat) {
      score <- k - 1
      if(score == 0) numer <- rep(0, length(theta))
      else numer <- a * (score * theta - sum(steps[1:score]))
      probs[, k] <- exp(numer)
    }
    probs <- probs / rowSums(probs)
    apply(probs, 1, function(p) sample(0:(n_cat-1), 1, prob=p))
  }

  # 3. Create Data
  data_sim <- data.frame(matrix(NA, nrow = N, ncol = J))
  colnames(data_sim) <- paste0("Item_", 1:J)

  # Binary Items (2PL)
  a_bin <- runif(n_bin, 0.8, 1.5)
  b_bin <- seq(-3, 3, length.out = n_bin)
  for(j in 1:n_bin) {
    prob <- 1 / (1 + exp(-(a_bin[j] * (true_theta - b_bin[j]))))
    data_sim[, j] <- rbinom(N, 1, prob)
  }

  # Polytomous Items (GPCM)
  # Item 6: 2 steps (-2, 2)
  data_sim[, 6] <- sim_gpcm(true_theta, a=1.0, steps=c(-2, 2))
  # Item 7: 5 steps
  data_sim[, 7] <- sim_gpcm(true_theta, a=1.2, steps=c(-5, -2.5, 0, 2.5, 5))

  # 4. Run Estimation without prior
  # Note: Wide theta_range needed due to SD=3 in simulation
  my_models <- c(rep("2PL", n_bin), rep("GPCM", n_poly))

  res <- mixed_irt(data = data_sim, model = my_models, method = "EM",
                   control = list(max_iter = 20, theta_range = c(-6, 6)))

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

  
  # 5. Run Estimation with prior (MAP)
  res_prior <- mixed_irt(data = data_sim, model = my_models, method = "EM",
                         control = list(max_iter = 20, theta_range = c(-6, 6),
                                       prior = list(
                                         "2PL" = list(
                                           a = function(x) dlnorm(x, 0, 0.5, log=TRUE),
                                           b = function(x) dnorm(x, 0, 2, log=TRUE)
                                         ),
                                         "GPCM" = list(
                                           a = function(x) dlnorm(x, 0, 0.5, log=TRUE),
                                           d = function(x) dnorm(x, 0, 2, log=TRUE)
                                         )
                                       )))
  head(res_prior$item_params)
  print(res_prior$model_fit)
  # --- Example 2: With Package Data ---
  data("ela2", package = "tirt")

  # Define Models (7 Binary, 3 Poly)
  real_models <- c(rep("2PL", 7), rep("GRM", 3))

  # Run Estimation
  real_res <- mixed_irt(ela2, model = real_models, method = "EM",
                        control = list(max_iter = 10))

  head(real_res$item_params)
  print(real_res$model_fit)
  

Mixture Item Response Theory Model (Latent-Class Rasch / 2PL)

Description

Fits a mixture item response theory model (Rost, 1990) in which the examinee population is assumed to consist of a small number of unobserved latent classes, each with its own set of item parameters. Mixture IRT models are used to detect qualitatively different response strategies, unmodeled subpopulations, or classes for which item difficulty ordering differs. The model is estimated by the Expectation-Maximization algorithm with a fixed standard-normal ability distribution within each class.

Usage

mixture_irt(data, n_class = 2, model = "Rasch", control = list())

Arguments

Value

A list containing:

• item_params: A data frame of class-specific item difficulties (and discriminations for the 2PL model), with one column per class.

• class_params: A data frame of estimated mixing proportions for each latent class.

• person_params: A data frame of posterior class-membership probabilities and the modal (most likely) class for each person.

• model_fit: A data frame with the log-likelihood, AIC, BIC, number of classes, and classification entropy.

References

Rost, J. (1990). Rasch models in latent classes: An integration of two approaches to item analysis. Applied Psychological Measurement, 14(3), 271-282.

Examples

  # Two classes with reversed difficulty ordering
  set.seed(2025)
  N <- 300; J <- 8
  b1 <- seq(-1.5, 1.5, length.out = J)
  b2 <- rev(b1)
  theta <- rnorm(N)
  cls <- rep(1:2, each = N / 2)
  resp <- matrix(0, N, J)
  for (i in 1:N) {
    b <- if (cls[i] == 1) b1 else b2
    resp[i, ] <- rbinom(J, 1, 1 / (1 + exp(-(theta[i] - b))))
  }
  df <- as.data.frame(resp); names(df) <- paste0("I", 1:J)

  fit <- mixture_irt(df, n_class = 2, model = "Rasch",
                     control = list(max_iter = 30, verbose = FALSE))
  fit$class_params
  head(fit$item_params)

Person Fit Statistics (Standardized Log-Likelihood, lz)

Description

Computes person-level fit statistics to detect examinees whose response patterns are unlikely under the fitted item response model (for example, careless responding, cheating, or unusual guessing). The function returns the l_z standardized log-likelihood index of Drasgow, Levine, and Williams (1985), which under model fit is approximately standard normal, and flags persons whose l_z falls below a critical value.

Usage

person_fit(data, item_params, theta, critical = -1.96, model = NULL, D = 1)

Arguments

Details

For a person with ability \theta, let l_0 be the log-likelihood of the observed responses. The index is l_z = (l_0 - E[l_0]) / \sqrt{\mathrm{Var}(l_0)}, where the expectation and variance are taken over the model-implied category probabilities at \theta. Large negative values indicate response patterns that are less likely than the model predicts.

Value

A data frame with one row per person and the columns:

• person: person index (row number).

• n_items: number of items answered.

• theta: the ability estimate used.

• loglik: the observed log-likelihood of the response pattern.

• lz: the standardized log-likelihood person-fit index.

• flag: logical; TRUE if lz < critical.

References

Drasgow, F., Levine, M. V., & Williams, E. A. (1985). Appropriateness measurement with polychotomous item response models and standardized indices. British Journal of Mathematical and Statistical Psychology, 38(1), 67-86.

Examples

  set.seed(123)
  sim <- sim_irt(n_people = 300,
                 item_structure = list(list(model = "2PL", n_items = 15)))
  fit <- binary_irt(sim$resp, model = "2PL", method = "EM",
                    control = list(max_iter = 15, verbose = FALSE))

  pf <- person_fit(sim$resp, fit$item_params, fit$person_params)
  head(pf)

  # How many examinees are flagged as misfitting?
  sum(pf$flag, na.rm = TRUE)

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())

Arguments

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.

Examples

  # --- Quick Example (small data) ---
  set.seed(1)
  n <- 50
  resp <- data.frame(
    I1 = sample(0:2, n, replace = TRUE),
    I2 = sample(0:2, n, replace = TRUE),
    I3 = sample(0:3, n, replace = TRUE)
  )
  res <- polytomous_irt(resp, model = "GPCM", method = "EM",
                        control = list(max_iter = 5, verbose = FALSE))
  head(res$item_params)
 
  # --- Full 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)

  # 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)
  

Test Reliability Indices

Description

Computes commonly reported reliability coefficients for a test. Depending on the information supplied it returns the empirical (marginal) reliability from latent-trait estimates and their standard errors, the model-based marginal reliability from the test information function, and the classical Cronbach's alpha from the raw responses. Reliability summarizes how consistently the test orders examinees and is a standard entry in every technical report.

Usage

reliability(
  person_params = NULL,
  data = NULL,
  item_params = NULL,
  prior_sd = 1,
  theta_range = c(-4, 4),
  n_points = 61,
  model = NULL,
  D = 1
)

Arguments

Details

The empirical reliability is \mathrm{Var}(\hat\theta) / (\mathrm{Var}(\hat\theta) + \overline{SE^2}). The model-based marginal reliability is \sigma^2 / (\sigma^2 + \overline{1/I(\theta)}), where the average error variance is taken over the population ability distribution. Cronbach's alpha is \frac{J}{J-1}\left(1 - \frac{\sum_j s_j^2}{s_T^2}\right).

Value

A data frame with the columns Index and Value, reporting whichever reliability coefficients could be computed from the supplied arguments:

• empirical_reliability: from person estimates and their SEs.

• marginal_reliability: from the test information function.

• cronbach_alpha: classical internal-consistency reliability.

• n_items, n_persons: sample descriptors when available.

Examples

  set.seed(123)
  sim <- sim_irt(n_people = 500,
                 item_structure = list(list(model = "2PL", n_items = 12)))
  fit <- binary_irt(sim$resp, model = "2PL", method = "EM",
                    control = list(max_iter = 15, verbose = FALSE))

  reliability(person_params = fit$person_params,
              data = sim$resp,
              item_params = fit$item_params)

Summed-Score to Theta Conversion Table

Description

Generates a raw-score (summed-score) to scale-score conversion table, mapping every possible integer summed score to its corresponding ability estimate and standard error. Such tables are standard in operational testing programs because they let each examinee's number-correct (or number-of-points) score be translated directly into a theta estimate without needing the full response pattern. Three scoring methods are supported: expected a posteriori ("EAP"), weighted likelihood ("WLE"), and maximum likelihood ("MLE").

Usage

score_table(
  item_params,
  method = c("EAP", "WLE", "MLE"),
  prior_mean = 0,
  prior_sd = 1,
  theta_range = c(-4, 4),
  n_points = 81,
  model = NULL,
  D = 1
)

Arguments

Details

The distribution of the summed score at each ability grid point is obtained with the Lord-Wingersky (1984) recursion, which works for any mixture of dichotomous and polytomous items. The "EAP" estimate is the mean of the posterior of theta given the summed score, with se equal to the posterior standard deviation. The "MLE" estimate maximizes the summed-score likelihood; the "WLE" estimate (Warm, 1989) maximizes the summed-score likelihood weighted by the square root of the test information, which reduces the bias of the MLE. At the two extreme summed scores the MLE is not finite and is reported at the boundary of theta_range with a missing standard error; the WLE remains finite there and is estimated normally.

Value

A data frame with one row per possible summed score and the columns:

• summed_score: the integer summed score (0 to the maximum possible score).

• theta: the ability estimate for that summed score.

• se: the standard error of the ability estimate.

References

Lord, F. M., & Wingersky, M. S. (1984). Comparison of IRT true-score and equipercentile observed-score equatings. Applied Psychological Measurement, 8(4), 453-461.

Warm, T. A. (1989). Weighted likelihood estimation of ability in item response theory. Psychometrika, 54(3), 427-450.

Examples

  set.seed(123)
  sim <- sim_irt(n_people = 400,
                 item_structure = list(list(model = "2PL", n_items = 10)))
  fit <- binary_irt(sim$resp, model = "2PL", method = "EM",
                    control = list(max_iter = 15, verbose = FALSE))

  # EAP conversion table (0 to 10 correct)
  score_table(fit$item_params, method = "EAP")

  
  # Maximum-likelihood conversion table
  score_table(fit$item_params, method = "MLE")
  

Simulate Item Response Theory Data

Description

Simulate item responses data. Support both dichotomous and polytomous responses. Provide an easy implementation with a few default settings.

Usage

sim_irt(
  n_people = 1000,
  item_structure = list(),
  theta = NULL,
  theta_mean = 0,
  theta_sd = 1
)

Arguments

Value

A list containing:

Examples

  # 1. Define the Test Blueprint
  # We want:
  # - 10 items using 2PL (medium difficulty)
  # - 5 items using 3PL (difficult, with guessing)
  # - 5 items using GPCM (4-point Likert scale)
  # - 5 items using GRM (5-point Likert scale)

  my_test_structure <- list(
    # Block 1: 2PL
    list(model = "2PL", n_items = 10, a = c(0.8, 1.2), b = c(-1, 1)),

    # Block 2: 3PL (Harder items, b from 1 to 2.5, fixing guessing at 0.2)
    list(model = "3PL", n_items = 5, a = c(1.0, 1.5), b = c(1.0, 2.5), c = 0.2),

    # Block 3: GPCM (Polytomous, 4 categories 0-3)
    list(model = "GPCM", n_items = 5, categories = 4, a = c(0.7, 1.3), b = c(-1, 1)),

    # Block 4: GRM (Polytomous, 5 categories 0-4)
    list(model = "GRM", n_items = 5, categories = 5, a = c(1.0, 2.0))
  )

  # 2. Run the Simulation
  # Define N and a specific Theta vector
  N <- 2000
  theta_vec <- rnorm(N, 0, 2)

  sim_data <- sim_irt(
    n_people = N,
    theta = theta_vec,
    item_structure = my_test_structure
  )

  # 3. Inspect the Output
  # The Response Matrix
  head(sim_data$resp)

  # The True Parameters (Useful for recovery studies)
  # Note how it aligns a, b, and threshold parameters (step_1, step_2...)
  head(sim_data$true_params)

Simulate Multidimensional Item Response Theory Data

Description

Simulates responses from compensatory multidimensional item response models, the multidimensional counterpart of sim_irt. It supports the multidimensional binary models estimated by mirt_binary (M-Rasch, M-2PL, M-3PL) as well as multidimensional polytomous models (multidimensional GPCM and GRM). Each block of items may load on any subset of the latent dimensions, so both simple-structure and complex-structure designs are easy to generate.

Usage

sim_mirt(
  n_people = 1000,
  item_structure = list(),
  dimension = 2,
  theta = NULL,
  theta_mean = 0,
  Sigma = NULL
)

Arguments

Details

Each element of item_structure is a list describing one block:

• model: one of "M2PL", "M3PL", "MRasch" (dichotomous) or "MGPCM", "MGRM" (polytomous). The non-prefixed names ("2PL", "GRM", ...) are also accepted.

• n_items: number of items in the block.

• dims: integer vector of the dimensions the block loads on (default: all dimensions). Slopes on the remaining dimensions are 0.

• a: discrimination/slope, given as a range c(lo, hi) to sample from, or a single fixed value (default c(0.8, 1.8)).

• d: intercept, given as a range or a fixed value (default c(-1.5, 1.5)). Note the compensatory model uses the intercept metric z = a'\theta + d, matching mirt_binary().

• c: lower asymptote for "M3PL" (default 0.2).

• categories: number of categories for polytomous blocks (default 3).

Value

A list containing:

See Also

mirt_binary, sim_irt

Examples

  # Two correlated dimensions, simple structure
  set.seed(2025)
  Sigma <- matrix(c(1, 0.4, 0.4, 1), 2, 2)

  design <- list(
    list(model = "M2PL", n_items = 8, dims = 1),          # loads on Dim 1
    list(model = "M2PL", n_items = 8, dims = 2),          # loads on Dim 2
    list(model = "M3PL", n_items = 4, dims = c(1, 2), c = 0.2)  # both dims
  )

  sim <- sim_mirt(n_people = 500, item_structure = design,
                  dimension = 2, Sigma = Sigma)

  head(sim$resp)
  head(sim$true_params)

  
  # Recover with the multidimensional estimator
  Q <- as.matrix(sim$true_params[, c("a_Dim1", "a_Dim2")] != 0) * 1
  fit <- mirt_binary(sim$resp, model = "2PL", dimension = 2,
                     control = list(Q_matrix = Q, max_iter = 10, verbose = FALSE))
  head(fit$item_params)
  

Simulate Mixed IRT and Testlet (TRT) Data

Description

Simulates a data set that mixes ordinary independent item response theory items with testlet (locally dependent) items in a single test form, the simulation counterpart of irt_trt. Independent blocks depend only on the primary trait (theta); testlet blocks additionally depend on a testlet-specific nuisance effect (gamma). Both dichotomous and polytomous formats are supported for the independent and the testlet parts, so the whole family of models handled by irt_trt can be generated at once.

Usage

sim_tirt(
  n_people = 1000,
  item_structure = list(),
  theta = NULL,
  theta_mean = 0,
  theta_sd = 1
)

Arguments

Details

Each element of item_structure is a list describing one block. The model name determines whether the block is independent or a testlet:

• Independent (theta only): "Rasch", "2PL", "3PL", "PCM", "GPCM", "GRM".

• Testlet (theta + gamma): "RaschT", "2PLT", "3PLT", "BiFT", "PCT", "GPCT", "GRT".

Other recognized keys are n_items, categories (polytomous), a, b, c (given as a fixed value or a c(lo, hi) range to sample from), s (testlet loading for "BiFT"), testlet_id (the testlet label, required for testlet blocks), testlet_var (variance of the gamma effect), and gamma_vector (a user-supplied gamma effect of length n_people).

Value

A list containing:

See Also

irt_trt, sim_trt, sim_irt

Examples

  # A form with independent items and two testlets
  set.seed(2025)
  design <- list(
    list(model = "2PL",  n_items = 8),                       # independent
    list(model = "GRM",  n_items = 4, categories = 3),       # independent poly
    list(model = "2PLT", n_items = 4, testlet_id = "P1",
         testlet_var = 0.6),                                 # testlet
    list(model = "GPCT", n_items = 3, categories = 3,
         testlet_id = "P2", testlet_var = 0.5)               # testlet poly
  )

  sim <- sim_tirt(n_people = 600, item_structure = design)

  head(sim$resp)
  sim$true_item_params[, c("item_id", "model", "testlet")]
  head(sim$true_person_params)

  
  # Feed the true structure straight into the joint estimator
  spec <- data.frame(
    item = sim$true_item_params$item_id,
    model = sim$true_item_params$model,
    testlet = sim$true_item_params$testlet,
    stringsAsFactors = FALSE
  )
  fit <- irt_trt(sim$resp, spec, method = "EM",
                 control = list(max_iter = 15, verbose = FALSE))
  head(fit$item_params)
  

Simulate Testlet Response Theory Data (Vector Supported Version)

Description

Simulate testlet responses data. Support both dichotomous and polytomous responses. Provide an easy implementation with a few default settings.

Usage

sim_trt(
  n_people = 1000,
  item_structure = list(),
  theta = NULL,
  theta_mean = 0,
  theta_sd = 1
)

Arguments

Value

A list containing:

Examples

  # =========================================================================
  # Example 1: Complex Testlet Design
  # =========================================================================
  # Define the Testlet Blueprint
  trt_design <- list(
    # Testlet 1: Rasch Testlet Model (High dependence: var=0.8)
    list(model = "RaschT", n_items = 5, testlet_id = "Read_A",
         testlet_var = 0.8, b = c(-1, 1)),

    # Testlet 2: 2PL Testlet Model (Default dependence: var=0.5)
    list(model = "2PLT", n_items = 5, testlet_id = "Read_B",
         a = c(0.7, 1.3)),

    # Testlet 3: Graded Response Testlet (Polytomous, 4 categories)
    list(model = "GRT", n_items = 4, testlet_id = "Survey",
         categories = 4, testlet_var = 0.2)
  )

  # Run Simulation
  trt_data <- sim_trt(n_people = 500, item_structure = trt_design)

  # Inspect Results
  # 1. Responses
  head(trt_data$resp)

  # 2. Item Parameters
  # (Notice 'testlet_loading' equals 'discrimination' for standard models)
  head(trt_data$true_item_params)

  # 3. Person Parameters (Ability + Gamma for each testlet)
  head(trt_data$true_person_params)

  # =========================================================================
  # Example 2: Manual Control (Theta, Gamma, and Parameters)
  # =========================================================================

  # 1. Manual Theta (e.g., everyone has high ability)
  manual_theta <- rep(2.0, 100)

  # 2. Manual Gamma (e.g., zero effect for T1)
  manual_gamma <- rep(0, 100)

  # 3. Item Parameters: Exact Match vs Range Sampling
  custom_structure <- list(
    # Case A: Manual Gamma Vector
    list(model = "2PLT", n_items = 5, testlet_id = "T1",
         gamma_vector = manual_gamma),

    # Case B: Exact Parameter Match (Length of 'a' equals n_items)
    list(model = "2PLT", n_items = 2, testlet_id = "T2",
         a = c(0.5, 2.5)),

    # Case C: Range Sampling (Length of 'a' is 2, but n_items != 2)
    list(model = "2PLT", n_items = 5, testlet_id = "T3",
         a = c(0.5, 2.5))
  )

  res_custom <- sim_trt(n_people = 100, theta = manual_theta,
                        item_structure = custom_structure)

  # Verify Manual Theta
  print(mean(res_custom$true_person_params$ability)) # Should be 2.0

  # Verify Manual Gamma (T1 should be 0)
  print(head(res_custom$true_person_params$testlet_T1))

  # Verify Exact Match (T2 discrimination should be 0.5 and 2.5)
  print(res_custom$true_item_params[res_custom$true_item_params$testlet_id == "T2",
                                    "discrimination"])

Test and Item Characteristic Curves (Expected Scores)

Description

Computes item characteristic curves (the expected score of each item as a function of ability) and the test characteristic curve (the expected total score as a function of ability). The test characteristic curve maps the latent ability scale onto the number-correct (true-score) scale and underlies true-score equating and score reporting; the item curves show how each item's expected score rises across the trait.

Usage

tcc(item_params, theta = seq(-4, 4, by = 0.1), model = NULL, D = 1)

Arguments

Value

A list with two elements:

• test_curve: a data frame with columns theta and expected_score (the test characteristic curve).

• item_curves: a numeric matrix with one row per item and one column per theta value, giving the expected score of each item.

See Also

item_info, score_table

Examples

  set.seed(123)
  sim <- sim_irt(n_people = 400,
                 item_structure = list(list(model = "2PL", n_items = 10)))
  fit <- binary_irt(sim$resp, model = "2PL", method = "EM",
                    control = list(max_iter = 15, verbose = FALSE))

  curves <- tcc(fit$item_params, theta = seq(-3, 3, by = 0.5))
  curves$test_curve

  # Expected score on the first item across ability
  curves$item_curves[1, ]

Test Information Function and Conditional Standard Error

Description

Computes the test information function (TIF) as the sum of all item information values at each ability point, together with the conditional standard error of measurement (SEM). The test information function summarizes where along the ability scale the test measures most precisely, and the conditional SEM (1 / \sqrt{TIF}) translates that precision back onto the theta metric. This is a routine part of test evaluation and form assembly in operational testing programs.

Usage

test_info(item_params, theta = seq(-4, 4, by = 0.1), model = NULL, D = 1)

Arguments

Value

A data frame with one row per theta point and the columns:

• theta: the ability value.

• test_info: the test information (sum of item information).

• sem: the conditional standard error of measurement, 1 / \sqrt{test\_info}.

• reliability: the marginal-style conditional reliability TIF / (TIF + 1) for a standard-normal ability scale.

See Also

item_info

Examples

  set.seed(123)
  sim <- sim_irt(n_people = 400,
                 item_structure = list(list(model = "2PL", n_items = 10)))
  fit <- binary_irt(sim$resp, model = "2PL", method = "EM",
                    control = list(max_iter = 15, verbose = FALSE))

  tif <- test_info(fit$item_params, theta = seq(-3, 3, by = 0.5))
  print(tif)

  # The ability where the test is most informative
  tif$theta[which.max(tif$test_info)]

Unidimensional Binary (Dichotomous) Testlet Response Theory Estimation

Description

Estimates item and person parameters for Unidimensional Binary (Dichotomous) Testlet response models using Penalized Expectation-Maximization or Joint Maximum Likelihood Estimation with stabilization.

Usage

trt_binary(
  data,
  group,
  model = c("RaschT", "2PLT", "3PLT", "BiFT"),
  method = c("EM", "MLE"),
  control = list()
)

Arguments

Value

A list containing:

• item_params: Estimated item parameters.

• person_params: Estimated person abilities and testlet effects.

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

Examples

  # --- Quick Example (small data) ---
  set.seed(1)
  n <- 50
  J <- 4
  theta <- rnorm(n)
  gamma <- rnorm(n, 0, 0.5)
  resp <- matrix(NA, n, J)
  colnames(resp) <- paste0("Item_", 1:J)
  b <- c(-0.5, 0, 0.5, 1)
  for(j in 1:J) {
    p <- 1 / (1 + exp(-(theta + gamma - b[j])))
    resp[, j] <- rbinom(n, 1, p)
  }
  res <- trt_binary(
    data = as.data.frame(resp),
    group = list(1:2, 3:4),
    model = "2PLT",
    method = "EM",
    control = list(max_iter = 5, verbose = FALSE)
  )
  head(res$item_params)

 
  # --- Full Example: Simulation (2PLT) ---
  set.seed(2025)
  n_persons <- 500
  n_testlets <- 3
  items_per_testlet <- 3
  n_items <- n_testlets * items_per_testlet

  # 1. Generate Parameters
  # Discrimination (a): Varying -> 2PLT
  a_true <- runif(n_items, 0.8, 1.5)
  # Difficulty (b)
  b_true <- seq(-1, 1, length.out = n_items)
  # Testlet Variances (Sigma)
  sigma_true <- c(1.0, 1.5, 2.0)

  # 2. Generate Person Params
  theta_true <- rnorm(n_persons, 0, 1)
  gamma_matrix <- matrix(0, nrow = n_persons, ncol = n_testlets)
  for(d in 1:n_testlets) {
    gamma_matrix[, d] <- rnorm(n_persons, 0, sigma_true[d])
  }

  # 3. Generate Responses
  resp_matrix <- matrix(0, nrow = n_persons, ncol = n_items)
  colnames(resp_matrix) <- paste0("Item_", 1:n_items)
  group_list <- list()

  idx_counter <- 1
  for(d in 1:n_testlets) {
    indices <- idx_counter:(idx_counter + items_per_testlet - 1)
    group_list[[d]] <- indices

    for(i in indices) {
      # 2PLT Model: a * (theta + gamma - b)
      lin <- a_true[i] * (theta_true + gamma_matrix[, d] - b_true[i])
      prob <- 1 / (1 + exp(-lin))
      resp_matrix[, i] <- rbinom(n_persons, 1, prob)
    }
    idx_counter <- idx_counter + items_per_testlet
  }
  df_sim <- as.data.frame(resp_matrix)

  # 4. Run Estimation
  # We use "2PLT" because data was generated with varying 'a'
  res <- trt_binary(
    data = df_sim,
    group = group_list,
    model = "2PLT",
    method = "EM",
    control = list(max_iter = 20, verbose = FALSE)
  )

  head(res$item_params)
  head(res$person_params)
  

Unidimensional Polytomous Testlet Response Theory Estimation

Description

Estimates item and person parameters for Polytomous Testlet models using Robust Newton-Raphson optimization.

Usage

trt_poly(
  data,
  group,
  model = c("GRT", "PCMT", "BiFT"),
  method = c("MLE", "EM"),
  control = list()
)

Arguments

Value

A list containing:

• item_params: A data frame of estimated item parameters.

• person_params: A data frame of estimated person abilities and testlet effects .

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

Examples

  # --- Example: Simulation (Mixed Categories GRT) ---
  set.seed(42)
  N <- 500; J <- 16

  # Define Groups (4 Testlets)
  groups <- list(c(1:4), c(5:8), c(9:12), c(13:16))

  # Define Categories (Binary, 3-cat, 4-cat, Mixed)
  # Items 1-4: 2 cats; 5-8: 3 cats; 9-12: 4 cats; 13-16: mixed
  cats <- c(rep(2, 4), rep(3, 4), rep(4, 4), 3, 5, 3, 5)

  # 1. Generate Parameters
  theta <- rnorm(N)
  # Gamma for 4 testlets (SD = 0.8)
  gamma <- matrix(rnorm(N * 4, 0, 0.8), N, 4)

  a <- rlnorm(J, 0, 0.2)
  b_list <- vector("list", J)

  # Generate Thresholds based on category count
  for(j in 1:J) {
    n_thresh <- cats[j] - 1
    if(n_thresh == 1) {
      b_list[[j]] <- rnorm(1)
    } else {
      # Spread thresholds
      b_list[[j]] <- sort(rnorm(1) + seq(-1, 1, length.out=n_thresh))
    }
  }

  # 2. Generate Responses (GRT Logic)
  resp <- matrix(NA, N, J)
  colnames(resp) <- paste0("Item_", 1:J)

  for(i in 1:N) {
    for(j in 1:J) {
      # Identify Testlet ID
      tid <- which(sapply(groups, function(x) j %in% x))
      eff <- theta[i] + gamma[i, tid]

      # Calculate Probabilities (Graded Response)
      K <- cats[j]
      probs <- numeric(K)
      P_prev <- 1
      for(k in 1:(K-1)) {
        term <- a[j] * (eff - b_list[[j]][k])
        P_star <- 1 / (1 + exp(-term))
        probs[k] <- P_prev - P_star
        P_prev <- P_star
      }
      probs[K] <- P_prev

      # Sample Response
      resp[i, j] <- sample(0:(K-1), 1, prob = probs)
    }
  }
  df_sim <- as.data.frame(resp)

  # 3. Run Estimation
  fit <- trt_poly(
    data = df_sim,
    group = groups,
    model = "GRT",
    method = "EM",
    control = list(max_iter = 20, verbose = FALSE)
  )

  head(fit$item_params)
  head(fit$person_params)