## Not run:
# # list of needed packages for the following examples
# packages <- scan(what="character")
# sirt TAM eRm CDM mirt ltm mokken psychotools psychomix
# psych
#
# # load packages. make an installation if necessary
# miceadds::library_install(packages)
#
# #############################################################################
# # EXAMPLE 1: Unidimensional models openness scale
# #############################################################################
#
# data(data.big5)
# # extract first 10 openness items
# items <- which( substring( colnames(data.big5) , 1 , 1 ) == "O" )[1:10]
# dat <- data.big5[ , items ]
# I <- ncol(dat)
# summary(dat)
# ## > colnames(dat)
# ## [1] "O3" "O8" "O13" "O18" "O23" "O28" "O33" "O38" "O43" "O48"
# # descriptive statistics
# psych::describe(dat)
#
# #****************
# # Model 1: Partial credit model
# #****************
#
# #-- M1a: rm.facets (in sirt)
# m1a <- sirt::rm.facets( dat )
# summary(m1a)
#
# #-- M1b: tam.mml (in TAM)
# m1b <- TAM::tam.mml( resp=dat )
# summary(m1b)
#
# #-- M1c: gdm (in CDM)
# theta.k <- seq(-6,6,len=21)
# m1c <- CDM::gdm( dat , irtmodel="1PL" ,theta.k=theta.k , skillspace="normal")
# summary(m1c)
# # compare results with loglinear skillspace
# m1c2 <- CDM::gdm( dat , irtmodel="1PL" ,theta.k=theta.k , skillspace="loglinear")
# summary(m1c2)
#
# #-- M1d: PCM (in eRm)
# m1d <- eRm::PCM( dat )
# summary(m1d)
#
# #-- M1e: gpcm (in ltm)
# m1e <- ltm::gpcm( dat , constraint = "1PL" , control=list(verbose=TRUE))
# summary(m1e)
#
# #-- M1f: mirt (in mirt)
# m1f <- mirt::mirt( dat , model=1 , itemtype="1PL" , verbose=TRUE)
# summary(m1f)
# coef(m1f)
#
# #-- M1g: PCModel.fit (in psychotools)
# mod1g <- psychotools::PCModel.fit(dat)
# summary(mod1g)
# plot(mod1g)
#
# #****************
# # Model 2: Generalized partial credit model
# #****************
#
# #-- M2a: rm.facets (in sirt)
# m2a <- sirt::rm.facets( dat , est.a.item=TRUE)
# summary(m2a)
# # Note that in rm.facets the mean of item discriminations is fixed to 1
#
# #-- M2b: tam.mml.2pl (in TAM)
# m2b <- TAM::tam.mml.2pl( resp=dat , irtmodel="GPCM")
# summary(m2b)
#
# #-- M2c: gdm (in CDM)
# m2c <- CDM::gdm( dat , irtmodel="2PL" ,theta.k=seq(-6,6,len=21) ,
# skillspace="normal" , standardized.latent=TRUE)
# summary(m2c)
#
# #-- M2d: gpcm (in ltm)
# m2d <- ltm::gpcm( dat , control=list(verbose=TRUE))
# summary(m2d)
#
# #-- M2e: mirt (in mirt)
# m2e <- mirt::mirt( dat , model=1 , itemtype="GPCM" , verbose=TRUE)
# summary(m2e)
# coef(m2e)
#
# #****************
# # Model 3: Nonparametric item response model
# #****************
#
# #-- M3a: ISOP and ADISOP model - isop.poly (in sirt)
# m3a <- sirt::isop.poly( dat )
# summary(m3a)
# plot(m3a)
#
# #-- M3b: Mokken scale analysis (in mokken)
# # Scalability coefficients
# mokken::coefH(dat)
# # Assumption of monotonicity
# monotonicity.list <- mokken::check.monotonicity(dat)
# summary(monotonicity.list)
# plot(monotonicity.list)
# # Assumption of non-intersecting ISRFs using method restscore
# restscore.list <- mokken::check.restscore(dat)
# summary(restscore.list)
# plot(restscore.list)
#
# #****************
# # Model 4: Graded response model
# #****************
#
# #-- M4a: mirt (in mirt)
# m4a <- mirt::mirt( dat , model=1 , itemtype="graded" , verbose=TRUE)
# print(m4a)
# mirt.wrapper.coef(m4a)
#
# #---- M4b: WLSMV estimation with cfa (in lavaan)
# lavmodel <- "F =~ O3__O48
# F ~~ 1*F
# "
# # transform lavaan syntax with lavaanify.IRT
# lavmodel <- TAM::lavaanify.IRT( lavmodel , items=colnames(dat) )$lavaan.syntax
# mod4b <- lavaan::cfa( data= as.data.frame(dat) , model=lavmodel, std.lv = TRUE,
# ordered=colnames(dat) , parameterization="theta")
# summary(mod4b , standardized=TRUE , fit.measures=TRUE , rsquare=TRUE)
# coef(mod4b)
#
# #****************
# # Model 5: Normally distributed residuals
# #****************
#
# #---- M5a: cfa (in lavaan)
# lavmodel <- "F =~ O3__O48
# F ~~ 1*F
# F ~ 0*1
# O3__O48 ~ 1
# "
# lavmodel <- TAM::lavaanify.IRT( lavmodel , items=colnames(dat) )$lavaan.syntax
# mod5a <- lavaan::cfa( data= as.data.frame(dat) , model=lavmodel, std.lv = TRUE ,
# estimator="MLR" )
# summary(mod5a , standardized=TRUE , fit.measures=TRUE , rsquare=TRUE)
#
# #---- M5b: mirt (in mirt)
#
# # create user defined function
# name <- 'normal'
# par <- c("d" = 1 , "a1" = 0.8 , "vy" = 1)
# est <- c(TRUE, TRUE,FALSE)
# P.normal <- function(par,Theta,ncat){
# d <- par[1]
# a1 <- par[2]
# vy <- par[3]
# psi <- vy - a1^2
# # expected values given Theta
# mui <- a1*Theta[,1] + d
# TP <- nrow(Theta)
# probs <- matrix( NA , nrow=TP, ncol= ncat )
# eps <- .01
# for (cc in 1:ncat){
# probs[,cc] <- stats::dnorm( cc , mean = mui , sd = sqrt( abs( psi + eps) ) )
# }
# psum <- matrix( rep(rowSums( probs ),each=ncat) , nrow=TP , ncol=ncat , byrow=TRUE)
# probs <- probs / psum
# return(probs)
# }
#
# # create item response function
# normal <- mirt::createItem(name, par=par, est=est, P=P.normal)
# customItems <- list("normal"=normal)
# itemtype <- rep( "normal",I)
# # define parameters to be estimated
# mod5b.pars <- mirt::mirt(dat, 1, itemtype=itemtype ,
# customItems=customItems , pars = "values")
# ind <- which( mod5b.pars$name == "vy")
# vy <- apply( dat , 2 , var , na.rm=TRUE )
# mod5b.pars[ ind , "value" ] <- vy
# ind <- which( mod5b.pars$name == "a1")
# mod5b.pars[ ind , "value" ] <- .5* sqrt(vy)
# ind <- which( mod5b.pars$name == "d")
# mod5b.pars[ ind , "value" ] <- colMeans( dat , na.rm=TRUE )
#
# # estimate model
# mod5b <- mirt::mirt(dat, 1, itemtype=itemtype , customItems=customItems ,
# pars = mod5b.pars , verbose=TRUE )
# sirt::mirt.wrapper.coef(mod5b)$coef
#
# # some item plots
# par(ask=TRUE)
# plot(mod5b, type = 'trace', layout = c(1,1))
# par(ask=FALSE)
# # Alternatively:
# sirt::mirt.wrapper.itemplot(mod5b)
# ## End(Not run)
Run the code above in your browser using DataLab