Example 1

A toy example for estimating the confirmatory IRT model using the proposed USP method.

## load the package
library(lvmcomp2)

## simulate data from the M2PL model under a confirmatory setting
## set parameters 
N <- 1000
J <- 20
K <- 2
set.seed(1)
# A_skelaton <- t(matrix(as.integer(intToBits(1:2^K)),32))[1:(2^K-1),1:K]
# Q <- A_skelaton[sample(2^K-1,J,replace = T),]
Q <- matrix(0, 20, 2)
for(k in 1:2){
  Q[(5*(k-1)+1):(5*k),k] <- 1
  Q[11:20,k] <- 1
}
A <- Q
A[which(Q>0)] <- runif(sum(Q),0.5,1.5)
d <- rnorm(J, 0, 1)

sigma <- matrix(0.4,K,K)
diag(sigma) <- 1
B <- chol(sigma)

## generate latent values and response data
theta <- matrix(rnorm(N*K,0,1), N)
theta <- theta %*% chol(sigma)
temp = cbind(rep(1,N), theta) %*% t(cbind(d, A))
prob_ori = 1/(1+exp(-temp))
response = matrix(rbinom(N*J, 1, prob = prob_ori), N, J)

# Initial values
Q <- A > 0
A0 <- (A>0)*1
d0 <- rep(0, J)
B0 <- diag(K)
theta0 <- matrix(rnorm(N*K,0,1), N)

## perform estimation
result_sa_conf <- sa_mirt_conf(response, A0, Q, d0, B0, theta, alpha = 0.51, step = 1, ave_from=500, max_steps = 1000)

## estimated loading matrix vs. the true loadings
data.frame("est_A" = result_sa_conf$A_hat, "true_A" = A)
## estimated intercept paramters vs. the true values
data.frame("est_d" = result_sa_conf$d_hat, "true_d" = d)
## estimated correlations vs. the true values
data.frame("est_sigma" = result_sa_conf$sigma_hat, "true_sigma"=sigma)

## standard errors for estimated parameters can be obtained
se_all <- sqrt(diag(solve(result_sa_conf$oakes)))