## ----setup, include = FALSE------------------------------------------------------------- knitr::opts_chunk$set(collapse = TRUE, comment = "#>") old_options <- options(width = 90) ## ----data, message = FALSE-------------------------------------------------------------- library(tppcm) # also attaches TAM data(data.Students, package = "CDM") sc <- data.Students[, c("sc1", "sc2", "sc3", "sc4")] sapply(sc, table, useNA = "ifany") ## ----fit-tppcm-------------------------------------------------------------------------- m_tppcm <- tam.mml.3pl(sc, E = tppcm(sc), est.variance = FALSE, verbose = FALSE) matrix(m_tppcm$gammaslope, nrow = 4, byrow = TRUE) # rows = items, columns = steps ## ----fit-nrm---------------------------------------------------------------------------- m_nrm <- tam.mml.2pl(sc, irtmodel = "2PL", verbose = FALSE) c(TPPCM = logLik(m_tppcm), NRM = logLik(m_nrm)) ## ----bfi, eval = FALSE------------------------------------------------------------------ # bfi <- psych::bfi[, paste0("A", 1:5)] - 1 # bfi is coded 1-6; tppcm needs 0-5 # bfi$A1 <- 5 - bfi$A1 # reverse-keyed item # m_bad <- tam.mml.3pl(bfi, E = tppcm(bfi), est.variance = FALSE, verbose = FALSE, # control = list(maxiter = 5000)) # m_bad$iter # 5000: not converged # range(m_bad$gammaslope) # -9.99 and 9.99: at the bound ## ----bfi-fix, eval = FALSE-------------------------------------------------------------- # m_ok <- tam.mml.3pl(bfi, E = tppcm(bfi), est.variance = FALSE, verbose = FALSE, # control = list(maxiter = 3000, max.increment = 0.1)) # m_nrm_bfi <- tam.mml.2pl(bfi, irtmodel = "2PL", verbose = FALSE) # same model, converges as is # c(TPPCM = logLik(m_ok), NRM = logLik(m_nrm_bfi)) ## ----fit-models------------------------------------------------------------------------- m_gpcmE <- tam.mml.3pl(sc, E = tppcm(sc, design = ~ item), est.variance = FALSE, verbose = FALSE) m_step <- tam.mml.3pl(sc, E = tppcm(sc, design = ~ step), est.variance = FALSE, verbose = FALSE) m_pcmE <- tam.mml.3pl(sc, E = tppcm(sc, design = ~ 1), est.variance = FALSE, verbose = FALSE) attr(tppcm(sc, design = ~ step), "index") # one parameter per column (step) m_step$gammaslope # gamma_1, gamma_2, gamma_3 ## ----pcm-equiv-------------------------------------------------------------------------- m_pcm <- tam.mml(sc, verbose = FALSE) c(common_a = unname(m_pcmE$gammaslope), tam_sd = sqrt(m_pcm$variance[1])) c(tppcm = logLik(m_pcmE), tam = logLik(m_pcm)) ## ----unequal---------------------------------------------------------------------------- data(data.timssAusTwn.scored, package = "TAM") timss <- data.timssAusTwn.scored[, 1:11] attr(tppcm(timss), "index") m_timss <- tam.mml.3pl(timss, E = tppcm(timss), est.variance = FALSE, verbose = FALSE) irt_pars(m_timss, fit = FALSE, se = FALSE) ## ----fit-mix---------------------------------------------------------------------------- idx <- rbind(1:3, 4:6, rep(7, 3), rep(8, 3)) m_mix <- tam.mml.3pl(sc, E = tppcm(sc, index = idx), est.variance = FALSE, verbose = FALSE) ## ----fit-fix---------------------------------------------------------------------------- m_fix <- tam.mml.3pl(sc, E = tppcm(sc), est.variance = FALSE, verbose = FALSE, gammaslope.fixed = cbind(2, 3)) ## ----fit-pcm-gpcm----------------------------------------------------------------------- m_gpcm <- tam.mml.2pl(sc, irtmodel = "GPCM", verbose = FALSE) lr <- IRT.compareModels(m_pcmE, m_gpcm, m_step, m_mix, m_tppcm)$LRtest nested <- !paste(lr$Model1, lr$Model2) %in% c("m_gpcm m_step", "m_step m_gpcm", "m_step m_mix", "m_mix m_step") lr[nested, ] ## ----irt-pars--------------------------------------------------------------------------- irt_pars(m_tppcm) ## ----irt-pars-si------------------------------------------------------------------------ irt_pars(m_gpcm, IRTpars = FALSE) ## ----irt-pars-pcm----------------------------------------------------------------------- irt_pars(m_pcm) ## ----steps-ci--------------------------------------------------------------------------- head(irt_pars(m_tppcm, long = TRUE, level = 0.90)) ## ----m2--------------------------------------------------------------------------------- m2(m_gpcm) # full M2: all univariate and bivariate margins m2(m_tppcm) m2(m_gpcm, "C2") # univariate margins plus bivariate cross products ## ----step-info, fig.width = 7, fig.height = 3.5----------------------------------------- inf <- step_info(m_tppcm) inf op <- par(mfrow = c(1, 2)); plot(inf); plot(inf, item = "sc1"); par(op) ## ----wald------------------------------------------------------------------------------- wald_test(m_tppcm) ## ----wald-mix--------------------------------------------------------------------------- wald_test(m_mix) # items 3 and 4 are GPCM items in m_mix ## ----score-tppcm------------------------------------------------------------------------ score_test(m_gpcm) ## ----score-rank1------------------------------------------------------------------------ score_test(m_gpcm, against = ~ item + step) ## ----item-tests------------------------------------------------------------------------- item_test(m_gpcm) ## ----mi--------------------------------------------------------------------------------- mi(m_gpcm) ## ----rank1------------------------------------------------------------------------------ m_rank1 <- xxirt_tppcm(as.matrix(sc), design = ~ item + step) irt_pars(m_rank1, restricted = TRUE) # alpha, gamma and their SEs irt_pars(m_rank1) ## ----rank1-score------------------------------------------------------------------------ score_test(m_rank1) ## ----rank1-lr--------------------------------------------------------------------------- suppressWarnings(IRT.compareModels(m_gpcm, m_rank1, m_tppcm))$LRtest ## ----rank1-vcov------------------------------------------------------------------------- round(sqrt(diag(xxirt_vcov(m_rank1)))[1:4], 3) ## ----dif-test--------------------------------------------------------------------------- cov <- data.timssAusTwn.scored[, c("IDCNTRY", "ITSEX")] cov$country <- factor(cov$IDCNTRY, labels = c("AUS", "TWN")) cov$sex <- factor(cov$ITSEX, labels = c("girl", "boy")) sp_timss <- get_parts(m_timss) # compute once, reuse below dif_test(sp_timss, cov$sex, by = "item") ## ----dif-impact------------------------------------------------------------------------- dif_test(sp_timss, cov$country, parm = "diff", by = "item") ## ----dif-mg----------------------------------------------------------------------------- m_mg <- tam.mml.3pl(timss, E = tppcm(timss), group = cov$country, est.variance = FALSE, verbose = FALSE) irt_pars(m_mg, fit = FALSE)$group ## ----dsf-screen------------------------------------------------------------------------- st <- dif_test(m_mg, cov$country, by = "step") st flagged <- st$unit[st$p_adj < 0.05] flagged ## ----dsf-leak--------------------------------------------------------------------------- dif_test(m_mg, cov$country, parm = c("M032166_disc1", "M032166_diff1")) ## ----tree-va---------------------------------------------------------------------------- data("VerbalAggression", package = "psychotools") tr_va <- tppcmtree(resp ~ gender + anger, data = VerbalAggression, design = ~ 1, parm = "diff", minsize = 50) partykit::character_split(tr_va$node$split, tr_va$data)$name ## ----get-parts-------------------------------------------------------------------------- get_parts(m_gpcm, "est") # slope parameters of the fitted model (item slopes for the GPCM) get_parts(m_gpcm, c("loglik", "npar")) dim(get_parts(m_gpcm, "information")) # observed information of the saturated model ## ----get-parts-all---------------------------------------------------------------------- sp <- get_parts(m_gpcm) sp logLik(sp) score_test(sp); mi(sp)[1:3, ] ## ----reset-options, include = FALSE------------------------------------------- options(old_options)