Adapted from the sirt package examples (sim.qm.ramsay).
library(mirt)
library(sirt)
set.seed(657)
# simulate data according to the Ramsay model
N <- 1000
# persons
I <- 11
# items
theta <- exp( rnorm( N ) ) # person ability
b <- exp( seq(-2,2,len=I)) # item difficulty
K <- rep( 3 , I )
# K parameter (=> guessing)
# apply simulation function
dat <- sim.qm.ramsay( theta , b , K )
mmliter <- 50
# maximum number of iterations
I <- ncol(dat)
fixed.K <- rep( 3 , I )
# Ramsay QM with fixed K parameter (K=3 in fixed.K specification)
mod1 <- rasch.mml2( dat , mmliter = mmliter , irtmodel = "ramsay.qm",progress = FALSE,
fixed.K = fixed.K )
summary(mod1)
## ------------------------------------------------------------
## sirt 3.12-66 (2022-05-16 12:27:54)
##
## Date of Analysis: 2023-07-05 16:12:11
## Time difference of 0.1008546 secs
## Computation time: 0.1008546
##
## Call:
## rasch.mml2(dat = dat, mmliter = mmliter, progress = FALSE, fixed.K = fixed.K,
## irtmodel = "ramsay.qm")
##
## Semiparametric Marginal Maximum Likelihood Estimation
## Function 'rasch.mml2'
##
## irtmodel= ramsay.qm
## Quotient Model (Ramsay, 1989)
## ------------------------------------------------------------
## Number of iterations = 34
## Deviance = 11746.2 | Log Likelihood = -5873.1
## Number of persons = 1000
## Number of estimated parameters = 12
## AIC = 11770.2 | penalty= 24 | AIC=-2*LL + 2*p
## AICc = 11770.52 | penalty= 24.32 | AICc=-2*LL + 2*p + 2*p*(p+1)/(n-p-1) (bias corrected AIC)
## BIC = 11829.09 | penalty= 82.89 | BIC=-2*LL + log(n)*p
## CAIC = 11841.09 | penalty= 94.89 | CAIC=-2*LL + [log(n)+1]*p (consistent AIC)
##
## Trait Distribution ( 21 Knots )
## Mean = 0
## SD = 0.991
## Skewness = 0 Note: log theta distribution is parametrized!
## Item Difficulty Distribution ( 11 Items )
## Mean = -0.064 SD = 1.317
## Distribution of Items Administered ( 11 Items )
## Mean = 11 SD = 0
##
## EAP Reliability: 0.713
## ------------------------------------------------------------
## Item Parameter
## item N p K est.K b log_b est.b guess.K emp.discrim
## 1 I01 1000 0.944 3 0 0.116 -2.158 1 0.25 0.453
## 2 I02 1000 0.889 3 0 0.185 -1.689 2 0.25 0.543
## 3 I03 1000 0.816 3 0 0.288 -1.246 3 0.25 0.569
## 4 I04 1000 0.730 3 0 0.433 -0.837 4 0.25 0.594
## 5 I05 1000 0.632 3 0 0.672 -0.397 5 0.25 0.551
## 6 I06 1000 0.540 3 0 1.012 0.012 6 0.25 0.497
## 7 I07 1000 0.469 3 0 1.536 0.429 7 0.25 0.395
## 8 I08 1000 0.402 3 0 2.191 0.784 8 0.25 0.361
## 9 I09 1000 0.378 3 0 2.954 1.083 9 0.25 0.237
## 10 I10 1000 0.315 3 0 5.038 1.617 10 0.25 0.160
## 11 I11 1000 0.297 3 0 5.462 1.698 11 0.25 0.182
#*** estimate Model 1 with user defined function in mirt package
# create user defined function for Ramsays quotient model
name <- "ramsayqm"
par <- c("K" = 3 , "b" = 1 )
est <- c(TRUE, TRUE)
P.ramsay <- function(par,Theta,ncat){
eps <- .01
K <- par[1]
b <- par[2]
num <- exp( exp( Theta[,1] ) / b )
denom <- K + num
P1 <- num / denom
P1 <- eps + ( 1 - 2*eps ) * P1
cbind(1-P1, P1)
}
# create item response function
ramsayqm <- mirt::createItem(name, par=par, est=est, P=P.ramsay)
# define parameters to be estimated
mod1m.pars <- mirt::mirt(dat, 1, rep( "ramsayqm",I) ,
customItems=list("ramsayqm"=ramsayqm), pars = "values")
mod1m.pars[ mod1m.pars$name == "K" , "est" ] <- FALSE
# define Theta design matrix
Theta <- matrix( seq(-3,3,len=10) , ncol=1)
# estimate model
mod1m <- mirt::mirt(dat, 1, rep( "ramsayqm",I) , customItems=list("ramsayqm"=ramsayqm ),
pars = mod1m.pars , verbose=FALSE ,
technical = list( customTheta=Theta , NCYCLES=50)
)
print(mod1m)
##
## Call:
## mirt::mirt(data = dat, model = 1, itemtype = rep("ramsayqm",
## I), pars = mod1m.pars, verbose = FALSE, technical = list(customTheta = Theta,
## NCYCLES = 50), customItems = list(ramsayqm = ramsayqm))
##
## Full-information item factor analysis with 1 factor(s).
## Converged within 1e-04 tolerance after 15 EM iterations.
## mirt version: 1.39.4
## M-step optimizer: BFGS
## EM acceleration: Ramsay
## Number of rectangular quadrature: 10
## Latent density type: Gaussian
##
## Log-likelihood = -5880.067
## Estimated parameters: 11
## AIC = 11782.13
## BIC = 11836.12; SABIC = 11801.18
## G2 (2036) = 925.39, p = 1
## RMSEA = 0, CFI = NaN, TLI = NaN
cmod1m <- coef( mod1m , simplify=TRUE)$items
# compare simulated and estimated values
dfr <- cbind( b , cmod1m[,'b'] , exp(mod1$item$b ) )
colnames(dfr) <- c("simulated" , "mirt" , "sirt_rasch.mml2")
round( dfr , 2 )
## simulated mirt sirt_rasch.mml2
## I01 0.14 0.11 0.12
## I02 0.20 0.17 0.18
## I03 0.30 0.27 0.29
## I04 0.45 0.42 0.43
## I05 0.67 0.66 0.67
## I06 1.00 1.00 1.01
## I07 1.49 1.54 1.54
## I08 2.23 2.22 2.19
## I09 3.32 3.02 2.95
## I10 4.95 5.24 5.04
## I11 7.39 5.65 5.46
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