Adapted from the sirt package examples (data.read).
library(mirt)
library(sirt)
library(mvtnorm)
data( data.read )
dat <- data.read
Q <- matrix( 0 , nrow=12,ncol=3)
Q[1:4,1] <- 1
Q[5:8,2] <- 1
Q[9:12,3] <- 1
# define discrete theta distribution with 3 dimensions
Theta <- matrix(c(0,0,0, 1,0,0, 0,1,0, 0,0,1, 1,1,0, 1,0,1, 0,1,1, 1,1,1), ncol=3, byrow=TRUE)
Theta
## [,1] [,2] [,3]
## [1,] 0 0 0
## [2,] 1 0 0
## [3,] 0 1 0
## [4,] 0 0 1
## [5,] 1 1 0
## [6,] 1 0 1
## [7,] 0 1 1
## [8,] 1 1 1
#-- Model 13a: din (in CDM)
mod13a <- CDM::din( dat , q.matrix=Q , rule="DINA", progress = FALSE)
summary(mod13a)
## 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:12:07
## Time difference of 0.03082705 secs
## Computation Time: 0.03082705
##
##
## Deviance = 3847.473 | Log-Likelihood= -1923.737
##
## Number of iterations: 40
##
## Number of item parameters: 24
## Number of skill class parameters: 7
##
## Information criteria:
## AIC = 3909.473
## BIC = 4027.056
##
## Mean of RMSEA item fit: 0.053
##
## Item parameters
## item guess slip IDI rmsea
## 1 A1 0.691 0.000 0.309 0.024
## 2 A2 0.491 0.031 0.478 0.040
## 3 A3 0.302 0.184 0.514 0.016
## 4 A4 0.244 0.337 0.419 0.049
## 5 B1 0.568 0.163 0.269 0.042
## 6 B2 0.329 0.344 0.327 0.071
## 7 B3 0.817 0.014 0.170 0.091
## 8 B4 0.431 0.104 0.465 0.044
## 9 C1 0.188 0.013 0.799 0.023
## 10 C2 0.050 0.239 0.711 0.140
## 11 C3 0.000 0.065 0.935 0.028
## 12 C4 0.000 0.212 0.788 0.070
##
## Marginal skill probabilities:
## skill.prob
## Skill1 0.5159
## Skill2 0.5421
## Skill3 0.9326
##
## Tetrachoric correlations among skill dimensions
## Skill1 Skill2 Skill3
## Skill1 1.0000 0.9727 0.3044
## Skill2 0.9727 1.0000 0.3129
## Skill3 0.3044 0.3129 1.0000
##
## Skill Pattern Probabilities
##
## 000 100 010 001 110 101 011 111
## 0.03879 0.00832 0.00948 0.39361 0.01080 0.01720 0.04220 0.47959
# compare used Theta distributions
cbind( Theta , mod13a$attribute.patt.splitted)
## Skill1 Skill2 Skill3
## [1,] 0 0 0 0 0 0
## [2,] 1 0 0 1 0 0
## [3,] 0 1 0 0 1 0
## [4,] 0 0 1 0 0 1
## [5,] 1 1 0 1 1 0
## [6,] 1 0 1 1 0 1
## [7,] 0 1 1 0 1 1
## [8,] 1 1 1 1 1 1
#-- Model 13b: gdm (in CDM)
mod13b <- CDM::gdm( dat , Qmatrix=Q , theta.k=Theta , skillspace="full", progress=FALSE)
summary(mod13b)
## -----------------------------------------------------------------------------
## CDM 8.2-6 (2022-08-25 15:43:23)
##
## Date of Analysis: 2023-07-05 16:12:10
## Time difference of 3.131896 secs
## Computation Time: 3.131896
##
## General Diagnostic Model
##
## Call:
## CDM::gdm(data = dat, theta.k = Theta, Qmatrix = Q, skillspace = "full",
## progress = FALSE)
##
## 328 Cases, 12 Items, 1 Group(s) , 3 Dimension(s)
## Saturated skill space
##
## -----------------------------------------------------------------------------
## Number of iterations= 1000
## Maximum number of iterations was reached.
##
## Deviance = 3865.63 | Log Likelihood = -1932.82
## Number of persons = 328
## Number of estimated parameters = 31
## Number of estimated item parameters = 24
## 12 Intercepts and 12 Slopes
## 0 centered intercepts and 0 centered slopes
## Number of estimated distribution parameters = 7
##
## AIC = 3928 | penalty = 62 | AIC = -2*LL + 2*p
## AICc = 3934 | penalty = 68.7 | AICc = -2*LL + 2*p + 2*p*(p+1)/(n-p-1) (bias corrected AIC)
## BIC = 4045 | penalty = 179.58 | BIC = -2*LL + log(n)*p
## CAIC = 4076 | penalty = 210.58 | CAIC = -2*LL + [log(n)+1]*p (consistent AIC)
##
## -----------------------------------------------------------------------------
## Trait Distribution
##
## M Trait:
## F1 F2 F3
## Group1 0.524 0.507 0.816
##
## SD Trait:
## F1 F2 F3
## Group1 0.499 0.5 0.387
##
## Skewness Trait:
## F1 F2 F3
## Group1 -0.098 -0.028 -1.633
##
## Correlations Trait:
## Group 1
## F1 F2 F3
## F1 1.000 0.861 0.368
## F2 0.861 1.000 0.349
## F3 0.368 0.349 1.000
##
## EAP Reliability:
## F1 F2 F3
## [1,] 0.709 0.651 0.657
## -----------------------------------------------------------------------------
## Item Parameters
## item N M b.Cat1 a.F1 a.F2 a.F3 itemfit.rmsea
## 1 A1 328 0.851 0.781 33.673 0.000 0.000 0.039
## 2 A2 328 0.738 -0.062 3.456 0.000 0.000 0.020
## 3 A3 328 0.567 -0.843 2.283 0.000 0.000 0.031
## 4 A4 328 0.460 -1.160 1.830 0.000 0.000 0.018
## 5 B1 328 0.713 0.273 0.000 1.500 0.000 0.045
## 6 B2 328 0.506 -0.657 0.000 1.348 0.000 0.054
## 7 B3 328 0.909 1.551 0.000 3.022 0.000 0.041
## 8 B4 328 0.683 -0.169 0.000 2.388 0.000 0.044
## 9 C1 328 0.933 0.887 0.000 0.000 33.995 0.066
## 10 C2 328 0.713 -0.685 0.000 0.000 2.079 0.103
## 11 C3 328 0.872 -0.104 0.000 0.000 3.323 0.049
## 12 C4 328 0.735 -1.040 0.000 0.000 2.709 0.033
##
## Mean of RMSEA item fit: 0.045
#-- Model 13c: mirt (in mirt)
# define mirt model
I <- ncol(dat) # I = 12
mirtmodel <- mirt::mirt.model("
F1 = 1-4
F2 = 5-8
F3 = 9-12
")
# 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 MEAN_3 75 0.000000 -Inf Inf FALSE none
## 76 all GROUP GroupPars COV_11 76 1.000000 1e-04 Inf FALSE none
## 77 all GROUP GroupPars COV_21 77 0.000000 -Inf Inf FALSE none
## 78 all GROUP GroupPars COV_31 78 0.000000 -Inf Inf FALSE none
## 79 all GROUP GroupPars COV_22 79 1.000000 1e-04 Inf FALSE none
## 80 all GROUP GroupPars COV_32 80 0.000000 -Inf Inf FALSE none
## 81 all GROUP GroupPars COV_33 81 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
## 78 NaN NaN
## 79 NaN NaN
## 80 NaN NaN
## 81 NaN NaN
#* define prior for latent class analysis
lca_prior <- function(Theta,Etable){
TP <- nrow(Theta)
if ( is.null(Etable) ){
prior <- rep( 1/TP , TP )
}
if ( ! is.null(Etable) ){
prior <- ( rowSums( Etable[ , seq(1,2*I,2)] ) + rowSums( Etable[ , seq(2,2*I,2)]))}
prior <- prior / sum(prior)
return(prior)
}
#* estimate model
mod13c <- 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
coef13c <- coef(mod13c, simplify=TRUE)$items
#* inspect estimated distribution
cbind(Theta, extract.mirt(mod13c, "Prior")[[1]])
## [,1] [,2] [,3] [,4]
## [1,] 0 0 0 0.040517629
## [2,] 1 0 0 0.008688677
## [3,] 0 1 0 0.007497199
## [4,] 0 0 1 0.415904645
## [5,] 1 1 0 0.010802342
## [6,] 1 0 1 0.041193529
## [7,] 0 1 1 0.009551117
## [8,] 1 1 1 0.465844863
#-* comparisons of estimated parameters
# extract guessing and slipping parameters from din
dfr <- as.data.frame(matrix(coef(mod13a)[1:24], ncol=2))
colnames(dfr) <- paste0("din.",c("guess","slip") )
# estimated parameters from gdm
dfr$gdm.guess <- plogis(mod13b$item$b)
dfr$gdm.slip <- 1 - plogis( rowSums(mod13b$item[,c("b.Cat1","a.F1","a.F2","a.F3")]))# estimated parameters from mirt
dfr$mirt.guess <- plogis( coef13c[,'d'] )
dfr$mirt.slip <- 1 - plogis( rowSums(coef13c[,c("d","a1","a2","a3")]) )
# comparison
round(dfr[, c(1,3,5,2,4,6)],3)
## din.guess gdm.guess mirt.guess din.slip gdm.slip mirt.slip
## 1 0.691 0.686 0.685 0.817 0.000 0.000
## 2 0.000 0.485 0.488 0.014 0.032 0.038
## 3 0.491 0.301 0.300 0.431 0.192 0.193
## 4 0.031 0.239 0.239 0.104 0.338 0.340
## 5 0.302 0.568 0.579 0.188 0.145 0.148
## 6 0.184 0.341 0.343 0.013 0.334 0.327
## 7 0.244 0.825 0.827 0.050 0.010 0.008
## 8 0.337 0.458 0.461 0.239 0.098 0.090
## 9 0.568 0.708 0.189 0.000 0.000 0.013
## 10 0.163 0.335 0.050 0.065 0.199 0.239
## 11 0.329 0.474 0.001 0.000 0.038 0.065
## 12 0.344 0.261 0.000 0.212 0.159 0.212
# estimated class sizes
dfr <- data.frame( "Theta" = Theta , "din"=mod13a$attribute.patt$class.prob ,
"gdm"=mod13b$pi.k , "mirt" = extract.mirt(mod13c, "Prior")[[1]])
# comparison
round(dfr,3)
## Theta.1 Theta.2 Theta.3 din gdm mirt
## 1 0 0 0 0.039 0.147 0.041
## 2 1 0 0 0.008 0.011 0.009
## 3 0 1 0 0.009 0.011 0.007
## 4 0 0 1 0.394 0.302 0.416
## 5 1 1 0 0.011 0.014 0.011
## 6 1 0 1 0.017 0.033 0.041
## 7 0 1 1 0.042 0.015 0.010
## 8 1 1 1 0.480 0.467 0.466
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] CDM_8.2-6 mvtnorm_1.2-2 sirt_3.12-66 ltm_1.2-0
## [5] polycor_0.8-1 msm_1.7 MASS_7.3-58.2 mirt_1.39.4
## [9] lattice_0.20-45 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 GPArotation_2023.3-1
## [28] splines_4.2.2 mime_0.12 permute_0.9-7
## [31] Deriv_4.1.3 markdown_1.7 vegan_2.6-4