## ----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)