## ----setup, include = FALSE-------------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
old_options <- options(width = 90)

## ----load, message = FALSE--------------------------------------------------------------
library(tppcm)   # also attaches TAM
data(data.Students, package = "CDM")
sc <- data.Students[, c("sc1", "sc2", "sc3", "sc4")]

## ----a1-raw, error = TRUE---------------------------------------------------------------
try({
data(Science, package = "mirt")
head(Science, 3)
tppcm(Science)                                   # scored 1-4
sci <- Science - 1
tppcm(cbind(id = seq_len(nrow(sci)), sci))       # an ID column slipped in
s_fac <- sci; s_fac$Comfort <- factor(s_fac$Comfort)
tppcm(s_fac)                                     # a factor column
s_one <- sci; s_one$Work <- 0
tppcm(s_one)                                     # an item with one category
s_gap <- sci; s_gap$Work[s_gap$Work == 2] <- 3
tppcm(s_gap)                                     # nobody in category 2 of Work
})

## ----a1-clean---------------------------------------------------------------------------
s_gap$Work[s_gap$Work == 3] <- 2
attr(tppcm(s_gap), "index")                      # Work now has 2 steps

## ----a2---------------------------------------------------------------------------------
fit_sat <- tam.mml.3pl(sc, E = tppcm(sc), est.variance = FALSE, verbose = FALSE)
fit_sat$gammaslope[1:6]               # the step discriminations

## ----a2-broken--------------------------------------------------------------------------
fit_bad <- tam.mml.3pl(sc, E = tppcm(sc), verbose = FALSE)
p_bad <- irt_pars(fit_bad, fit = FALSE)

## ----a3---------------------------------------------------------------------------------
fit_gpcm <- tam.mml.3pl(sc, E = tppcm(sc, design = ~ item), est.variance = FALSE, verbose = FALSE)
fit_step <- tam.mml.3pl(sc, E = tppcm(sc, design = ~ step), est.variance = FALSE, verbose = FALSE)
fit_pcm  <- tam.mml.3pl(sc, E = tppcm(sc, design = ~ 1),  est.variance = FALSE, verbose = FALSE)
dimnames(tppcm(sc, design = ~ step))[[4]]                  # parameter names
idx <- rbind(1:3, 4:6, rep(7, 3), rep(8, 3))      # only sc1 and sc2 saturated
fit_mix <- tam.mml.3pl(sc, E = tppcm(sc, index = idx), est.variance = FALSE, verbose = FALSE)
fit_fix <- tam.mml.3pl(sc, E = tppcm(sc), est.variance = FALSE, verbose = FALSE,
                       gammaslope.fixed = cbind(2, 3))                  # sc1 step 2 fixed at 3

## ----a3-broken, error = TRUE------------------------------------------------------------
try({
tppcm(sc, design = ~ item + step)
})

## ----a4---------------------------------------------------------------------------------
data(data.timssAusTwn.scored, package = "TAM")
timss <- data.timssAusTwn.scored[, 1:11]
fit_timss <- tam.mml.3pl(timss, E = tppcm(timss), est.variance = FALSE, verbose = FALSE)
irt_pars(fit_timss, se = FALSE)

## ----a5---------------------------------------------------------------------------------
fit_r1 <- xxirt_tppcm(as.matrix(sc), design = ~ item + step)
irt_pars(fit_r1, restricted = TRUE)              # alpha_i and gamma_l
dim(xxirt_vcov(fit_r1))                          # covariance matrix of the xxirt fit

## ----b1---------------------------------------------------------------------------------
irt_pars(fit_sat)

## ----b1-forms---------------------------------------------------------------------------
irt_pars(fit_sat, IRTpars = FALSE, fit = FALSE, se = FALSE)   # slope-intercept: slope * theta + int
head(irt_pars(fit_sat, long = TRUE, level = 0.90))            # long format with intervals
irt_pars(fit_gpcm, restricted = TRUE)                         # parameters of the fitted model itself

## ----b1-kable---------------------------------------------------------------------------
p <- irt_pars(fit_sat, fit = FALSE)
knitr::kable(cbind(p$items[, 1:3], p$se[, 1:3]), digits = 2,
             col.names = c(paste0("a", 1:3), paste0("SE(a", 1:3, ")")))

## ----b2---------------------------------------------------------------------------------
st <- irt_pars(fit_sat, long = TRUE, level = 0.90)
head(st[st$par == "disc", ])

## ----b3---------------------------------------------------------------------------------
get_parts(fit_gpcm, "est")
get_parts(fit_gpcm, c("loglik", "npar", "nobs"))
dim(get_parts(fit_gpcm, "information"))
head(get_parts(fit_gpcm, "estfun")[, 1:4], 3)

## ----b3-reuse---------------------------------------------------------------------------
sp <- get_parts(fit_gpcm)
sp
logLik(sp); nobs(sp)

## ----b3-broken, error = TRUE------------------------------------------------------------
try({
get_parts(fit_gpcm, "hessian")      # no such component
})

## ----b4, fig.width = 7, fig.height = 3.5------------------------------------------------
inf <- step_info(fit_sat)
inf
op <- par(mfrow = c(1, 2)); plot(inf); plot(inf, item = "sc1"); par(op)

## ----c1---------------------------------------------------------------------------------
wald_test(fit_sat)

## ----c1-broken--------------------------------------------------------------------------
wald_test(fit_gpcm)          # in a GPCM the steps are equal by construction

## ----c2---------------------------------------------------------------------------------
score_test(fit_gpcm)                    # GPCM against the saturated TPPCM
score_test(fit_gpcm, against = ~ item + step)           # GPCM against the product form
item_test(fit_gpcm)                     # which items
head(mi(fit_gpcm), 5)                   # which steps (modification indices and EPCs)

## ----c2-broken, error = TRUE------------------------------------------------------------
try({
score_test(fit_sat)                     # already saturated
score_test(fit_timss, against = ~ item + step)          # TIMSS items have 1 or 2 steps
})

## ----c3---------------------------------------------------------------------------------
lr <- IRT.compareModels(fit_pcm, fit_gpcm, fit_step, fit_sat)$LRtest
lr

## ----c3-nested--------------------------------------------------------------------------
lr[-4, ]

## ----c3-mixed---------------------------------------------------------------------------
suppressWarnings(IRT.compareModels(fit_gpcm, fit_r1, fit_sat))$LRtest

## ----c4---------------------------------------------------------------------------------
m2(fit_gpcm)                 # full M2
m2(fit_sat)
m2(fit_gpcm, "C2")           # C2, margins reduced to moments
m2(fit_timss)                # missing data and unequal categories are fine

## ----c4-broken, error = TRUE------------------------------------------------------------
try({
w <- runif(nrow(sc), 0.5, 2)
fit_w <- tam.mml.2pl(sc, irtmodel = "GPCM", pweights = w, verbose = FALSE)
m2(fit_w)                    # not available with weights
})

## ----d1---------------------------------------------------------------------------------
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(fit_timss)                 # reused below
d_disc <- dif_test(sp_timss, cov$sex, parm = "disc")   # every step discrimination
d_disc[d_disc$p < 0.10, ]
dif_test(sp_timss, cov$sex, parm = "diff", by = "item")

## ----d2-impact--------------------------------------------------------------------------
dif_test(sp_timss, cov$country, parm = "diff", by = "item")

## ----d2-mg------------------------------------------------------------------------------
fit_mg <- tam.mml.3pl(timss, E = tppcm(timss), group = cov$country,
                      est.variance = FALSE, verbose = FALSE)
irt_pars(fit_mg, fit = FALSE)$group
sp_mg <- get_parts(fit_mg)
dif_test(sp_mg, cov$country, by = "item")             # all parameters of an item jointly

## ----d3---------------------------------------------------------------------------------
data("VerbalAggression", package = "psychotools")
va <- VerbalAggression
fit_va <- tam.mml.3pl(va$resp, E = tppcm(va$resp, design = ~ item), est.variance = FALSE, verbose = FALSE)
sp_va <- get_parts(fit_va)
d_anger <- dif_test(sp_va, va$anger, parm = "diff", by = "item")        # numeric: Cramer-von Mises
head(d_anger[order(d_anger$p), ], 5)
anger3 <- cut(va$anger, quantile(va$anger, 0:3 / 3), include.lowest = TRUE, ordered_result = TRUE)
dif_test(sp_va, anger3, parm = "disc", by = "block")                     # ordered: weighted max LM
dif_test(sp_va, va$anger, parm = "disc", by = "block", functional = "maxLM")

## ----d3-broken, error = TRUE------------------------------------------------------------
try({
dif_test(sp_va, replace(va$anger, 1, NA))                                # covariate with NA
dif_test(sp_va, va$anger[-1])                                            # wrong length
})

## ----d4---------------------------------------------------------------------------------
tr_diff <- tppcmtree(resp ~ gender + anger, data = va, design = ~ 1, parm = "diff", minsize = 50)
partykit::width(tr_diff)
partykit::character_split(tr_diff$node$split, tr_diff$data)$name
tr_disc <- tppcmtree(resp ~ gender + anger, data = va, design = ~ item, parm = "disc", minsize = 50)
partykit::width(tr_disc)

## ----d4-broken, error = TRUE------------------------------------------------------------
try({
tppcmtree(resp ~ gender + anger, data = va, design = ~ item, parm = "diff", minsize = 50)
})

## ----d5-broken, error = TRUE------------------------------------------------------------
try({
dif_test(fit_mg, cov$country, parm = "disc", by = "step")   # a step is tested jointly
dsf(fit_timss)                                              # needs a multigroup fit
})

## ----e1---------------------------------------------------------------------------------
disc <- rbind(c(1, 1.5, 0.7), c(1, 1.5, 0.7), c(1, 1, NA))
diff <- rbind(c(-1, 0, 1),    c(-1, 0, 1),    c(-0.5, 0.5, NA))
head(sim_tppcm(5, disc, diff))

## ----e1-power---------------------------------------------------------------------------
set.seed(3)
A <- outer(c(1, 1, 1.2, 0.8, 1), c(1, 1.5, 0.7)); B <- matrix(c(-1, 0, 1), 5, 3, byrow = TRUE)
p <- replicate(10, {
  d <- sim_tppcm(300, A, B)
  score_test(tam.mml.2pl(d, irtmodel = "GPCM", verbose = FALSE))$p.value
})
mean(p < 0.05)

## ----e1-broken, error = TRUE------------------------------------------------------------
try({
sim_tppcm(10, disc, diff[, 1:2])                 # shapes differ
})

## ----e3---------------------------------------------------------------------------------
irt_pars(fit_w, fit = FALSE)$se[1:2, ]                                            # weights enter the SEs
fit_w1 <- tam.mml.2pl(sc, irtmodel = "GPCM", pweights = w / mean(w), verbose = FALSE)
irt_pars(fit_w1, fit = FALSE)$se[1:2, ]                                           # after normalizing
irt_pars(fit_mg, fit = FALSE)$group                                               # group means and variances
fit_lr <- tam.mml.3pl(sc, E = tppcm(sc), Y = data.Students[, c("female", "urban")],
                      est.variance = FALSE, verbose = FALSE)
p_lr <- irt_pars(fit_lr, fit = FALSE)                                             # warns

## ----reset-options, include = FALSE-------------------------------------------
options(old_options)