Skip to content

Example 1

Siliang Zhang edited this page May 2, 2022 · 4 revisions

Simulation for an M2PL model under the confirmatory setting

using the "Unified stochastic proximal algorithm" USP()

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(NK,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(NJ, 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(NK,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)))

Clone this wiki locally