用mirt包估计2PL参数 难度参数不正常恶化

用mirt包跑2PL的异常数据 设置了10%的被试会在%30的题目上随机作答，也就是P=0.5，在总异常作答只有3%的情况下，估出来的参数应该和无异常相差无几，但是难度参数却极大的恶化了，θ和区分度只有轻微的变化，不知道哪里出了差错

library(mirt)

set.seed(1)
N <- 1000  # number of participants
J <- 30  # number of items
D <- 1.7
U <- matrix(1, N, J)  # response m atrix
p <- matrix(NA, N, J)  # probability of answering all item correctly

# Simulate response -------------------------------------------------------
theta <- rnorm(N)  # ability
b <- rnorm(J)  # difficulty
theta[theta > 3] <- 3
theta[theta < -3] <- -3
b[b > 3] <- 3
b[b < -3] <- -3
b <- matrix(b,1,J)
theta <- matrix(theta,N,1)
a <- runif(30,0,2)

# method 1, conventional way

  
for(i in 1:N){
  for(j in 1 :J){
    p[i, j] <- 1 / (1 + exp(-a[j]* (theta[i] - b[j])))
    rand_num <- runif(1)
    if (p[i, j] < rand_num) {
      U[i, j] <- 0
    }
  }
}
rn <- sample(N,N*0.1)
for(i in 1:N){
  if(i %in% rn)
    p[i, sample(J,J*0.3)] <-(0.5)
}
for(i in 1:N){
  for(j in 1 :J){
    rand_num <- runif(1)
    if (p[i, j] < rand_num) {
      U[i, j] <- 0
    }
  }
}
  # 建立作答矩阵 ------------------------------------------------------------------
  
colnames(U) <- c("i1","i2","i3","i4","i5","i6","i7","i8","i9","i10","i11","i12",
                   "i13","i14","i15","i16","i17","i18","i19","20","i21","i22","i23","i24",
                   "i25","i26","i27","i28","i29","i30")

#................拟合............
nihe <- function(data){
  mod <- mirt(data,1)
  theta1<-round(fscores(mod),2)
  
  #..........能力，题目参数...............
  t <- coef(mod,IRTpars=T)
  vv<-matrix(NA,30,4)
  for (z in 1:J) {
    vv[z,] <- t[[z]]
  }
  a1 <- vv[,1]
  b1 <- vv[,2]
  
  #................结果输出.................
  bias <- mean(theta - theta1)
  t_mse<-sum((theta - theta1)^2)/N
  t_Rmse<-t_mse^0.5
  
  a_bias <- mean(a - a1)
  a_mse<-sum((a - a1)^2)/J
  a_Rmse<-a_mse^0.5
  
  b_bias <- mean(b - b1)
  b_mse<-sum((b - b1)^2)/J
  b_Rmse<-b_mse^0.5
  shuju <- c(bias,t_Rmse,a_bias,a_Rmse,b_bias,b_Rmse)
  return(shuju)
}


z_0_0 <- matrix(NA,100,6)
for(i in 1:100){
  z_0_0[i,] <- nihe(U)
}
mz_0_0 <- apply(z_0_0,2,mean)

无异常和异常的结果如图 a的RMSE无太大变化但是bias确恶化严重，在1000个被试只有3%异常作答的情况下真的很奇怪。

IRT模型的参数估计方法（EM算法和MCMC算法）
https://blog.csdn.net/weixin_41524411/article/details/90702670