请求协助将自定义R代码应用于目标数据文件

Alright, let's get your custom IRT ability estimation function connected to your data file properly. First, I noticed your code cuts off mid-calculation, so I'll start by completing that function with standard Newton-Raphson logic for theta estimation—this is the standard approach for IRT ability scoring. Then we'll walk through data prep, applying the function to your dataset, and checking the results.

Step 1: Complete the Custom Ability Function #

Here's the full, functional version of your code with the missing convergence logic added:

library(plyr); library(dplyr)

ability <- function(mdl, u, b, a, c){
  J <- length(b)
  
  # Set default parameters based on model type
  if (mdl == 1 | mdl == 2 | missing(c)) {
    c <- rep(0, J)
  }
  if (mdl == 1 | missing(a)) {
    a <- rep(1, J)
  }
  
  x <- sum(u)
  # Handle extreme response patterns (all correct or all incorrect)
  if (x == 0) {
    th <- -log(2*J)
  }
  if (x == J) {
    th <- log(2*J)
  }
  
  # For non-extreme patterns, use Newton-Raphson to estimate theta
  if (x > 0 & x < J) {
    th <- 0  # Initial guess for theta
    tol <- 1e-6  # Convergence threshold
    max_iter <- 100  # Max iterations to prevent infinite loops
    iter <- 0
    
    repeat {
      iter <- iter + 1
      # Calculate item response probabilities and their derivatives
      pstar <- 1 / (1 + exp(-a * (th - b)))
      phat <- c + (1 - c) * pstar
      
      # First derivative of log-likelihood (score function)
      d1 <- sum(a * (u - phat) * (1 - c) * exp(-a*(th - b)) / (1 + exp(-a*(th - b)))^2)
      # Second derivative (negative expected Fisher information)
      d2 <- -sum(a^2 * (1 - c)^2 * exp(-a*(th - b)) / (1 + exp(-a*(th - b)))^2 * phat * (1 - phat/(1 - c)))
      
      # Update theta estimate
      th_new <- th - d1/d2
      
      # Check for convergence
      if (abs(th_new - th) < tol || iter >= max_iter) {
        th <- th_new
        break
      }
      th <- th_new
    }
  }
  
  return(th)
}

Step 2: Prepare Your Data #

Your dataset needs to follow this structure:

• One row per survey respondent/test taker
• One column per item (binary responses: 1 = correct, 0 = incorrect)
• You'll also need pre-defined item parameters:
  • b: Vector of item difficulty scores (length = number of items)
  • a: Vector of item discrimination scores (only required for 2PL/3PL models; defaults to 1 for 1PL)
  • c: Vector of item guessing scores (only required for 3PL models; defaults to 0 for 1PL/2PL)

Load your data (example using a CSV file):

# Load your data file (replace with your file path)
irt_data <- read.csv("your_data_file.csv")

# Extract the response matrix (adjust column selection to match your item columns)
# Example: if items are named item_1 to item_10
response_matrix <- irt_data %>% select(starts_with("item_"))

# Define your item parameters (replace with your actual values)
item_b <- c(-1.2, -0.5, 0.0, 0.3, 0.8, 1.1, -0.7, 0.2, 0.6, 1.3)  # 10 items
item_a <- rep(1, 10)  # Discrimination (1PL model)
item_c <- rep(0, 10)  # Guessing (no guessing for 1PL)

Step 3: Apply the Function to All Respondents #

Use either dplyr (tidyverse style) or plyr to compute ability scores for every row in your dataset:

Option 1: Dplyr (Tidy Workflow)

# Add theta (ability) scores to your original dataset
irt_data_with_ability <- irt_data %>%
  rowwise() %>%
  mutate(theta = ability(
    mdl = 1,  # 1=1PL, 2=2PL, 3=3PL; adjust based on your model
    u = c_across(starts_with("item_")),  # Select all item response columns
    b = item_b,
    a = item_a,
    c = item_c
  )) %>%
  ungroup()

Option 2: Plyr (Legacy Workflow)

# Calculate theta scores for each respondent
theta_scores <- adply(response_matrix, 1, function(row) {
  ability(mdl = 1, u = unlist(row), b = item_b, a = item_a, c = item_c)
})

# Merge scores back with original data
irt_data_with_ability <- cbind(irt_data, theta = theta_scores$V1)

Step 4: Validate the Results #

Check that your ability estimates make sense with basic summaries and visualizations:

# Summary statistics of ability scores
summary(irt_data_with_ability$theta)

# Plot the distribution of theta estimates
hist(irt_data_with_ability$theta, 
     main = "Distribution of Respondent Ability Estimates",
     xlab = "Theta (Ability Score)",
     col = "lightblue")

Key Tips #

• Model Adjustments: Change the mdl argument to match your IRT model:
  • mdl=1: 1PL (Rasch) model (uses default a=1, c=0)
  • mdl=2: 2PL model (uses default c=0)
  • mdl=3: 3PL model (requires all three parameters)
• Extreme Responses: The function handles all-correct or all-incorrect responses with a conservative initial theta to avoid numerical issues.
• Convergence: If you need tighter or looser convergence, adjust the tol (tolerance) or max_iter values inside the ability function.

内容的提问来源于stack exchange，提问作者Omer