Adapted from the sirt package examples (data.read).
library(mirt)
library(sirt)
## - sirt 3.12-66 (2022-05-16 12:27:54)
library(mvtnorm)
data(data.read)
dat <- data.read
I <- ncol(dat)
# define some simple Q-Matrix (does not really make in this application)
Q <- matrix( 0 , nrow=12,ncol=2)
Q[1:4,1] <- 1
Q[5:8,2] <- 1
Q[9:12,1:2] <- 1
# define discrete theta distribution with 3 dimensions
Theta <- matrix(c(0,0, 1,0, 0,1, 1,1), byrow=TRUE, ncol=2)
Theta
## [,1] [,2]
## [1,] 0 0
## [2,] 1 0
## [3,] 0 1
## [4,] 1 1
#-- Model 14a: din (in CDM)
mod14a <- CDM::din( dat , q.matrix=Q , rule="DINA", progress = FALSE)
summary(mod14a)
## CDM 8.2-6 (Built 2022-08-25 15:43:23)
## Call:
## CDM::din(data = dat, q.matrix = Q, rule = "DINA", progress = FALSE)
##
## Date of Analysis: 2023-07-05 16:11:51
## Time difference of 0.02211118 secs
## Computation Time: 0.02211118
##
##
## Deviance = 3915.31 | Log-Likelihood= -1957.655
##
## Number of iterations: 41
##
## Number of item parameters: 24
## Number of skill class parameters: 3
##
## Information criteria:
## AIC = 3969.31
## BIC = 4071.722
##
## Mean of RMSEA item fit: 0.028
##
## Item parameters
## item guess slip IDI rmsea
## 1 A1 0.674 0.030 0.295 0.011
## 2 A2 0.423 0.049 0.527 0.017
## 3 A3 0.258 0.224 0.519 0.012
## 4 A4 0.245 0.394 0.361 0.018
## 5 B1 0.534 0.166 0.300 0.026
## 6 B2 0.338 0.382 0.280 0.025
## 7 B3 0.796 0.016 0.188 0.019
## 8 B4 0.421 0.142 0.436 0.033
## 9 C1 0.850 0.000 0.150 0.043
## 10 C2 0.480 0.097 0.423 0.037
## 11 C3 0.746 0.026 0.228 0.054
## 12 C4 0.575 0.136 0.289 0.038
##
## Marginal skill probabilities:
## skill.prob
## Skill1 0.5964
## Skill2 0.5996
##
## Tetrachoric correlations among skill dimensions
## Skill1 Skill2
## Skill1 1.0000 0.9572
## Skill2 0.9572 1.0000
##
## Skill Pattern Probabilities
##
## 00 10 01 11
## 0.35663 0.04377 0.04693 0.55267
# compare used Theta distributions
cbind( Theta , mod14a$attribute.patt.splitted)
## Skill1 Skill2
## [1,] 0 0 0 0
## [2,] 1 0 1 0
## [3,] 0 1 0 1
## [4,] 1 1 1 1
#-- Model 14b: mirt (in mirt)
# define mirt model
I <- ncol(dat) # I = 12
mirtmodel <- mirt::mirt.model("
F1 = 1-4
F2 = 5-8
(F1*F2) = 9-12
")
#-> constructions like (F1*F2*F3) are also allowed in mirt.model
# get parameters
mod.pars <- mirt(dat, model=mirtmodel , pars = "values")
# starting values d parameters (transformed guessing parameters)
ind <- which( mod.pars$name == "d" )
mod.pars[ind,"value"] <- qlogis(.2)
# starting values transformed slipping parameters
ind <- which( ( mod.pars$name %in% paste0("a",1:3) ) & ( mod.pars$est ) )
mod.pars[ind,"value"] <- qlogis(.8) - qlogis(.2)
mod.pars
## group item class name parnum value lbound ubound est prior.type
## 1 all A1 dich a1 1 2.772589 -Inf Inf TRUE none
## 2 all A1 dich a2 2 0.000000 -Inf Inf FALSE none
## 3 all A1 dich a3 3 0.000000 -Inf Inf FALSE none
## 4 all A1 dich d 4 -1.386294 -Inf Inf TRUE none
## 5 all A1 dich g 5 0.000000 0e+00 1 FALSE none
## 6 all A1 dich u 6 1.000000 0e+00 1 FALSE none
## 7 all A2 dich a1 7 2.772589 -Inf Inf TRUE none
## 8 all A2 dich a2 8 0.000000 -Inf Inf FALSE none
## 9 all A2 dich a3 9 0.000000 -Inf Inf FALSE none
## 10 all A2 dich d 10 -1.386294 -Inf Inf TRUE none
## 11 all A2 dich g 11 0.000000 0e+00 1 FALSE none
## 12 all A2 dich u 12 1.000000 0e+00 1 FALSE none
## 13 all A3 dich a1 13 2.772589 -Inf Inf TRUE none
## 14 all A3 dich a2 14 0.000000 -Inf Inf FALSE none
## 15 all A3 dich a3 15 0.000000 -Inf Inf FALSE none
## 16 all A3 dich d 16 -1.386294 -Inf Inf TRUE none
## 17 all A3 dich g 17 0.000000 0e+00 1 FALSE none
## 18 all A3 dich u 18 1.000000 0e+00 1 FALSE none
## 19 all A4 dich a1 19 2.772589 -Inf Inf TRUE none
## 20 all A4 dich a2 20 0.000000 -Inf Inf FALSE none
## 21 all A4 dich a3 21 0.000000 -Inf Inf FALSE none
## 22 all A4 dich d 22 -1.386294 -Inf Inf TRUE none
## 23 all A4 dich g 23 0.000000 0e+00 1 FALSE none
## 24 all A4 dich u 24 1.000000 0e+00 1 FALSE none
## 25 all B1 dich a1 25 0.000000 -Inf Inf FALSE none
## 26 all B1 dich a2 26 2.772589 -Inf Inf TRUE none
## 27 all B1 dich a3 27 0.000000 -Inf Inf FALSE none
## 28 all B1 dich d 28 -1.386294 -Inf Inf TRUE none
## 29 all B1 dich g 29 0.000000 0e+00 1 FALSE none
## 30 all B1 dich u 30 1.000000 0e+00 1 FALSE none
## 31 all B2 dich a1 31 0.000000 -Inf Inf FALSE none
## 32 all B2 dich a2 32 2.772589 -Inf Inf TRUE none
## 33 all B2 dich a3 33 0.000000 -Inf Inf FALSE none
## 34 all B2 dich d 34 -1.386294 -Inf Inf TRUE none
## 35 all B2 dich g 35 0.000000 0e+00 1 FALSE none
## 36 all B2 dich u 36 1.000000 0e+00 1 FALSE none
## 37 all B3 dich a1 37 0.000000 -Inf Inf FALSE none
## 38 all B3 dich a2 38 2.772589 -Inf Inf TRUE none
## 39 all B3 dich a3 39 0.000000 -Inf Inf FALSE none
## 40 all B3 dich d 40 -1.386294 -Inf Inf TRUE none
## 41 all B3 dich g 41 0.000000 0e+00 1 FALSE none
## 42 all B3 dich u 42 1.000000 0e+00 1 FALSE none
## 43 all B4 dich a1 43 0.000000 -Inf Inf FALSE none
## 44 all B4 dich a2 44 2.772589 -Inf Inf TRUE none
## 45 all B4 dich a3 45 0.000000 -Inf Inf FALSE none
## 46 all B4 dich d 46 -1.386294 -Inf Inf TRUE none
## 47 all B4 dich g 47 0.000000 0e+00 1 FALSE none
## 48 all B4 dich u 48 1.000000 0e+00 1 FALSE none
## 49 all C1 dich a1 49 0.000000 -Inf Inf FALSE none
## 50 all C1 dich a2 50 0.000000 -Inf Inf FALSE none
## 51 all C1 dich a3 51 2.772589 -Inf Inf TRUE none
## 52 all C1 dich d 52 -1.386294 -Inf Inf TRUE none
## 53 all C1 dich g 53 0.000000 0e+00 1 FALSE none
## 54 all C1 dich u 54 1.000000 0e+00 1 FALSE none
## 55 all C2 dich a1 55 0.000000 -Inf Inf FALSE none
## 56 all C2 dich a2 56 0.000000 -Inf Inf FALSE none
## 57 all C2 dich a3 57 2.772589 -Inf Inf TRUE none
## 58 all C2 dich d 58 -1.386294 -Inf Inf TRUE none
## 59 all C2 dich g 59 0.000000 0e+00 1 FALSE none
## 60 all C2 dich u 60 1.000000 0e+00 1 FALSE none
## 61 all C3 dich a1 61 0.000000 -Inf Inf FALSE none
## 62 all C3 dich a2 62 0.000000 -Inf Inf FALSE none
## 63 all C3 dich a3 63 2.772589 -Inf Inf TRUE none
## 64 all C3 dich d 64 -1.386294 -Inf Inf TRUE none
## 65 all C3 dich g 65 0.000000 0e+00 1 FALSE none
## 66 all C3 dich u 66 1.000000 0e+00 1 FALSE none
## 67 all C4 dich a1 67 0.000000 -Inf Inf FALSE none
## 68 all C4 dich a2 68 0.000000 -Inf Inf FALSE none
## 69 all C4 dich a3 69 2.772589 -Inf Inf TRUE none
## 70 all C4 dich d 70 -1.386294 -Inf Inf TRUE none
## 71 all C4 dich g 71 0.000000 0e+00 1 FALSE none
## 72 all C4 dich u 72 1.000000 0e+00 1 FALSE none
## 73 all GROUP GroupPars MEAN_1 73 0.000000 -Inf Inf FALSE none
## 74 all GROUP GroupPars MEAN_2 74 0.000000 -Inf Inf FALSE none
## 75 all GROUP GroupPars COV_11 75 1.000000 1e-04 Inf FALSE none
## 76 all GROUP GroupPars COV_21 76 0.000000 -Inf Inf FALSE none
## 77 all GROUP GroupPars COV_22 77 1.000000 1e-04 Inf FALSE none
## prior_1 prior_2
## 1 NaN NaN
## 2 NaN NaN
## 3 NaN NaN
## 4 NaN NaN
## 5 NaN NaN
## 6 NaN NaN
## 7 NaN NaN
## 8 NaN NaN
## 9 NaN NaN
## 10 NaN NaN
## 11 NaN NaN
## 12 NaN NaN
## 13 NaN NaN
## 14 NaN NaN
## 15 NaN NaN
## 16 NaN NaN
## 17 NaN NaN
## 18 NaN NaN
## 19 NaN NaN
## 20 NaN NaN
## 21 NaN NaN
## 22 NaN NaN
## 23 NaN NaN
## 24 NaN NaN
## 25 NaN NaN
## 26 NaN NaN
## 27 NaN NaN
## 28 NaN NaN
## 29 NaN NaN
## 30 NaN NaN
## 31 NaN NaN
## 32 NaN NaN
## 33 NaN NaN
## 34 NaN NaN
## 35 NaN NaN
## 36 NaN NaN
## 37 NaN NaN
## 38 NaN NaN
## 39 NaN NaN
## 40 NaN NaN
## 41 NaN NaN
## 42 NaN NaN
## 43 NaN NaN
## 44 NaN NaN
## 45 NaN NaN
## 46 NaN NaN
## 47 NaN NaN
## 48 NaN NaN
## 49 NaN NaN
## 50 NaN NaN
## 51 NaN NaN
## 52 NaN NaN
## 53 NaN NaN
## 54 NaN NaN
## 55 NaN NaN
## 56 NaN NaN
## 57 NaN NaN
## 58 NaN NaN
## 59 NaN NaN
## 60 NaN NaN
## 61 NaN NaN
## 62 NaN NaN
## 63 NaN NaN
## 64 NaN NaN
## 65 NaN NaN
## 66 NaN NaN
## 67 NaN NaN
## 68 NaN NaN
## 69 NaN NaN
## 70 NaN NaN
## 71 NaN NaN
## 72 NaN NaN
## 73 NaN NaN
## 74 NaN NaN
## 75 NaN NaN
## 76 NaN NaN
## 77 NaN NaN
lca_prior <- function(Theta,Etable){
# number of latent Theta classes
TP <- nrow(Theta)
# prior in initial iteration
if ( is.null(Etable) ){
prior <- rep( 1/TP , TP )
}
# process Etable (this is correct for datasets without missing data)
if ( ! is.null(Etable) ){
# sum over correct and incorrect expected responses
prior <- ( rowSums( Etable[ , seq(1,2*I,2)] ) + rowSums( Etable[ , seq(2,2*I,2)]))}
prior <- prior / sum(prior)
return(prior)
}
#* estimate model
mod14b <- mirt(dat, mirtmodel , technical = list(
customTheta=Theta , customPriorFun = lca_prior) ,
pars = mod.pars , verbose=FALSE )
# estimated parameters from the user customized prior distribution.
nest <- as.integer(sum(mod.pars$est) + 2)
#* extract item parameters
coef14b <- coef(mod14b, simplify=TRUE)$items
#-* comparisons of estimated parameters
# extract guessing and slipping parameters from din
dfr <- data.frame(matrix(coef(mod14a)[1:24], ncol=2, byrow=T))
colnames(dfr) <- paste0("din.",c("guess","slip") )
# estimated parameters from mirt
dfr$mirt.guess <- plogis( coef14b[,'d'] )
dfr$mirt.slip <- 1 - plogis( rowSums(coef14b[,c("d","a1","a2","a3")]) )
# comparison
round(dfr[, c(1,3,2,4)],3)
## din.guess mirt.guess din.slip mirt.slip
## 1 0.674 0.671 0.030 0.030
## 2 0.423 0.420 0.049 0.050
## 3 0.258 0.254 0.224 0.225
## 4 0.245 0.243 0.394 0.395
## 5 0.534 0.543 0.166 0.164
## 6 0.338 0.348 0.382 0.380
## 7 0.796 0.803 0.016 0.015
## 8 0.421 0.437 0.142 0.140
## 9 0.850 0.851 0.000 0.000
## 10 0.480 0.480 0.097 0.097
## 11 0.746 0.746 0.026 0.026
## 12 0.575 0.577 0.136 0.137
dfr <- data.frame( "Theta" = Theta , "din"=mod14a$attribute.patt$class.prob ,
"mirt" = extract.mirt(mod14b, "Prior")[[1]])
# comparison
round(dfr,3)
## Theta.1 Theta.2 din mirt
## 1 0 0 0.357 0.370
## 2 1 0 0.044 0.049
## 3 0 1 0.047 0.030
## 4 1 1 0.553 0.551
sessionInfo()
## R version 4.2.2 Patched (2022-11-10 r83330)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.6 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] mvtnorm_1.2-2 sirt_3.12-66 ltm_1.2-0 polycor_0.8-1
## [5] msm_1.7 MASS_7.3-58.2 mirt_1.39.4 lattice_0.20-45
## [9] knitr_1.43
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.10 highr_0.10 compiler_4.2.2
## [4] dcurver_0.9.2 tools_4.2.2 evaluate_0.21
## [7] lifecycle_1.0.3 gtable_0.3.3 nlme_3.1-162
## [10] mgcv_1.8-41 rlang_1.1.1 Matrix_1.5-3
## [13] cli_3.6.1 rstudioapi_0.14 commonmark_1.9.0
## [16] parallel_4.2.2 TAM_4.1-4 expm_0.999-7
## [19] xfun_0.39 gridExtra_2.3 cluster_2.1.4
## [22] admisc_0.32 grid_4.2.2 glue_1.6.2
## [25] pbapply_1.7-0 survival_3.4-0 CDM_8.2-6
## [28] GPArotation_2023.3-1 splines_4.2.2 mime_0.12
## [31] permute_0.9-7 Deriv_4.1.3 markdown_1.7
## [34] vegan_2.6-4