模拟研究

模拟研究是教育测量研究中常见的组成部分，也是重要组成部分。其目的在于说明在模拟情况为真的情况下，模型的统计性能如何。我们可以根据我们的理论假设来设定或者调整数据生成的过程，由此探讨在不同条件下模型的统计性能如何。

此外，模拟研究还有以下优势

a) 可以操纵多个条件(e.g., 被试数量; 题目难度分布)，如果使用真实数据，其耗费将是巨大到不可承受的。

b) 模拟研究可以知道能力和难度等模型参数的真值，可以比较真值和估计值之间的差异。而实证研究是不知道的，不能比较。

c) 在模拟研究中，由于数据生成过程是我们自己设定的，是已知的，可以完全排除其他无关因素的干扰。而真实数据中，由于测验过程的复杂性，可能纳入了太多的无关变量。

因此无论是学术研究还是实践应用，模拟研究均被广泛采用。

在下面的部分中，将展示如何对1PL进行模拟研究。主要使用到的语言为R，使用到的包为JAGS。

模拟研究可分为3部分，数据生成，参数估计，模型评估。

数据生成负责根据真值生成观测数据，参数估计则根据观测数据估计真值(注意，参数估计程序不会用到真值)，模型评估则通常评估参数估计和真值的差异。

1.数据生成

数据生成的部分要根据理论假设去模仿学生真实的作答情况。在当前这个例子中，我们假设学生的数据生成过程符合1PL。因此我们根据1PL的公式去生成数据，如下。

\(P_{i}(\theta_n)= \frac{e^{\theta_n-b_{i}}}{1+e^{\theta_n-b_{i}}}\)

其中能力\(\theta_n\)为人层面的参数能力，\(b_{i}\)为题目层面的因素题目难度。学生能力高且题目难度低时，作答正确概率最高。

我们假设学生人数N=500，题目数量I=20 。将数据生成过程分成3步，人层面参数的生成，题目层面参数的生成，计算作答正确概率。

# 设定题目参数个数和人的参数个数
N <- 500
I <- 20

# 人层面参数的生成(N)
theta <- as.matrix(rnorm(N, mean = 0, sd = 1), ncol=1)
# 题目层面参数的生成(I)
b <- as.matrix(rnorm(I, mean = 0, sd = 1), ncol=1)

# 计算作答正确概率(N*I)和作答结果(N*I)
correct_p <- matrix(NA, N, I)
response <- matrix(NA, N, I)
for (n in 1:N){
  for (i in 1:I){
    # 作答正确概率
    correct_p[n,i] <- exp(theta[n]-b[i])/(1+exp(theta[n]-b[i]))
    # 作答结果
    response[n,i] <- rbinom(1,1,correct_p[n,i])
  }
}

2.参数估计

在这部分中，我们介绍如何使用JAGS实现参数估计。

JAGS代码

首先我们要用JAGS语法定义1PL并保存至1PL.bug。JAGS语法可见JAGS文档

https://people.stat.sc.edu/hansont/stat740/jags_user_manual.pdf

以下是JAGS的1PL代码：


model{
  # 计算过程
  for (n in 1:N){
    for (i in 1:I){
      p[n,i] <- ilogit(theta[n]-b[i])
      response[n,i] ~ dbern(p[n,i])
    }
  }

# 参数先验
for (n in 1:N){
  theta[n] ~ dnorm(0,1)
}
for (i in 1:I){
  b[i] ~ dnorm(0, 0.25)
}

}

在【计算过程】这一部分和数据生成很类似，仅仅是从R的语法替换成了JGAS的语法。其中ilogit(p) = exp(p)/(1+exp(p))。

在【参数先验】这一部分中，我们给定了参数的先验，dnorm为正态先验，第一个参数为均值，而第二个参数为密度(1/方差)。

R代码

随后我们利用R2jags这个包调用JAGS对我们的模型进行估计。

library(R2jags)
# 模型输入数据，这里包括作答response，学生数量N，以及题目数量I
jags.data <- list("response", "N", "I") 

# 配置并运行参数估计程序
# parameters.to.save:需要保存的变量|model：模型路径| n.chains：独立运行的链条数
# n.iter：总迭代次数| n.burnin：预热次数
sim<-jags(data=jags.data,inits=NULL,parameters.to.save=c("b","theta")
          ,model="1PL.bug", n.chains=2,n.iter=5000,n.burnin=2500,DIC=T)
 
# 获取EAP作为点估计
theta_hat <- sim$BUGSoutput$mean$theta
b_hat <-  sim$BUGSoutput$mean$b

3. 模型评估

通常我们采用所有参数Bias的均值Mean(Bias)和RMSE的均值Mean(RMSE)来评估模型的返真性能，也就是真值和我们估计值之间的差异

\[Mean(Bias) =\frac{\sum_{i=1}^{P}\sum_{i=1}^{K} \left(\hat{X_{pk}}-X_{pk}\right)}{P \times K} \]

其中\(\hat{X_i}\)为估计值，而\(X_{i}\)为真值，P为参数个数，K为重复次数。在我们的例子中，仅重复了一次，因此K=1 。

Bias计算的其实是所有参数和估计值差异的平均值，评估的是参数估计是否有偏，而非参数估计的精确性。

\[Mean(R M S E)=\sum_{i=1}^{P} \frac{{\sqrt{\frac{\sum_{i=1}^{K}\left(\hat{X_{pk}}-X_{pk}\right)^{2}}{K}}}}{P} \]

RMSE计算的是每个参数在不同重复次数的估计精度如何，评估的则是估计的精确程度。

在当前例子中K=1

代码实现如下：


Bias_theta <- mean(as.matrix(theta_hat, ncol=1) - theta)
RMSE_theta <- mean(((as.matrix(theta_hat, ncol=1) - theta)^2)^0.5)

整体代码：

#############################
########【数据生成】###########
#############################
# 设定题目参数个数和人的参数个数
N <- 500
I <- 20

# 人层面参数的生成(N)
theta <- as.matrix(rnorm(N, mean = 0, sd = 1), ncol=1)
# 题目层面参数的生成(I)
b <- as.matrix(rnorm(I, mean = 0, sd = 1), ncol=1)

# 计算作答正确概率(N*I)和作答结果(N*I)
correct_p <- matrix(NA, N, I)
response <- matrix(NA, N, I)
for (n in 1:N){
  for (i in 1:I){
    # 作答正确概率
    correct_p[n,i] <- exp(theta[n]-b[i])/(1+exp(theta[n]-b[i]))
    # 作答结果
    response[n,i] <- rbinom(1,1,correct_p[n,i])
  }
}

#############################
########【参数估计】###########
#############################
library(R2jags)
# 模型输入数据，这里包括作答response，学生数量N，以及题目数量I
jags.data <- list("response", "N", "I") 

# 配置并运行参数估计程序
# parameters.to.save:需要保存的变量|model：模型路径| n.chains：独立运行的链条数
# n.iter：总迭代次数| n.burnin：预热次数
sim<-jags(data=jags.data,inits=NULL,parameters.to.save=c("b","theta")
          ,model="1PL.bug", n.chains=2,n.iter=5000,n.burnin=2500,DIC=T)
 
# 获取EAP作为点估计
theta_hat <- sim$BUGSoutput$mean$theta
b_hat <-  sim$BUGSoutput$mean$b

#############################
########【模型评估】###########
#############################
Bias_theta <- mean(as.matrix(theta_hat, ncol=1) - theta)
RMSE_theta <- mean(((as.matrix(theta_hat, ncol=1) - theta)^2)^0.5)

526互联

1PL模拟研究代码片段

模拟研究

1.数据生成

2.参数估计

JAGS代码

R代码

3. 模型评估