###################################################
### preliminaries
###################################################
library("strucchange")
library("psychotree")
library("psychotools")
library("psychomix")
library("multcomp")
set.seed(1090)


###################################################
### data
###################################################
library("psychotools")
data("MathExam14W", package = "psychotools")
mex <- subset(MathExam14W, nsolved > 0 & nsolved < 13)
head(mex$solved)
plot(mex$solved)


###################################################
### raschmodel
###################################################
mr <- raschmodel(mex$solved)
plot(mr, type = "profile")
plot(mr, type = "piplot")


###################################################
### raschmodel-group
###################################################
mr1 <- raschmodel(subset(mex, group == 1)$solved)
mr2 <- raschmodel(subset(mex, group == 2)$solved)
-2 * as.numeric(logLik(mr) - (logLik(mr1) + logLik(mr2)))


###################################################
### raschmodel-classicaltests
###################################################
mrc <- c(
  "LR" = -2 * as.numeric(logLik(mr) - (logLik(mr1) + logLik(mr2))),
  "Wald" = summary(anchortest(solved ~ group, data = mex)$final_tests, test = Chisqtest())$test$SSH,
  "Score/LM" = unname(sctest(mr, order.by = mex$group, vcov = "info", functional = "LMuo")$statistic)
)
round(mrc, digits = 1)


###################################################
### profileplot-1-fig
###################################################
plot(mr1, parg = list(ref = 1), ref = FALSE, ylim = c(-2.6, 2.6), col = 2)
plot(mr2, parg = list(ref = 1), ref = FALSE, add = TRUE, col = 4)
legend("topleft", paste("Group", 1:2), pch = 19, col = c(2, 4), bty = "n")


###################################################
### profileplot-10-fig
###################################################
plot(mr1, parg = list(ref = 10), ref = FALSE, ylim = c(-2.6, 2.6), col = 2)
plot(mr2, parg = list(ref = 10), ref = FALSE, add = TRUE, col = 4)
legend("topleft", paste("Group", 1:2), pch = 19, col = c(2, 4), bty = "n")


###################################################
### anchortest
###################################################
ma <- anchortest(solved ~ group, data = mex, adjust = "single-step")
plot(ma)


###################################################
### anchortest-fig
###################################################
par(mar = c(5, 8, 4, 1))
plot(ma, main = paste("Anchor item:", paste(ma$anchor_items, collapse = ", ")))


###################################################
### strucchange
###################################################
library("strucchange")
mex1 <- subset(mex, group == 1)
set.seed(1)
sctest(mr1, order.by = mex1$tests, vcov = "info", functional = "maxLMo", plot = TRUE)


###################################################
### raschtree-fig
###################################################
mex <- transform(mex,
  tests    = ordered(tests),
  nsolved  = ordered(nsolved),
  attempt  = ordered(attempt),
  semester = ordered(semester)
)
set.seed(1)
mrt <- raschtree(solved ~ group + tests + nsolved + gender + attempt + study + semester,
  data = mex, vcov = "info", minsize = 50, ordinal = "L2", nrep = 1e5)
plot(mrt)


###################################################
### raschmix-fig
###################################################
set.seed(1)
mrm <- raschmix(mex1$solved, k = 2, scores = "meanvar")
plot(mrm)
print(mrm)