# set.seed(123)
# 
# ## for the dataset ECPE
# library(CDM)
# Q <- data.ecpe$q.matrix
# Y <- data.ecpe$data[, -1]
# 
# # library(openxlsx)
# # wb <- createWorkbook()
# # addWorksheet(wb, "Q_matrix")
# # writeData(wb, sheet = "Q_matrix", x = Q)
# # addWorksheet(wb, "Response_Data")
# # writeData(wb, sheet = "Response_Data", x = Y)
# # saveWorkbook(wb, file = "ECPE_data.xlsx", overwrite = TRUE)
# 
# print(Q)
# 
# ######## CDM ########
# library(Qval)
# 
# is.Qident(Q)
# 
# CDM.obj <- CDM(Y, Q, model = "GDINA", mono.constraint = TRUE)
# summary(CDM.obj)
# 
# plot(CDM.obj)
# 
# 
# ######## validation with CDM ########
# Qval.obj <- validation(Y, Q, CDM.obj=CDM.obj,
#                        method="GDI", iter.level="no")
# 
# 
# ######## GDI ########
# Qval.GDI.obj <- validation(Y, Q, model="GDINA",
#                            mono.constraint = TRUE,
#                            method="GDI", search.method="SSA",
#                            iter.level="test.att", maxitr=150,
#                            verbose=1)
# summary(Qval.GDI.obj)
# 
# 
# ######## Hull ########
# Qval.Hull.obj <- validation(Y, Q, model="GDINA",
#                             mono.constraint = TRUE,
#                             method="Hull", criter="PVAF",
#                             search.method="ESA",
#                             iter.level="item", maxitr=150,
#                             verbose=1)
# 
# Q.Hull.sug <- extract(Qval.Hull.obj, "Q.sug")
# head(Q.Hull.sug)
# 
# 
# ######## Hull Plot ########
# plot(Qval.Hull.obj, i=1)
# 
# 
# ######## Wald ########
# Qval.Wald.obj <- validation(Y, Q, model="GDINA",
#                             mono.constraint = TRUE,
#                             method="Wald", alpha.level=0.05,
#                             search.method="stepwise",
#                             iter.level="no", verbose=1)
# Q.Wald.sug <- extract(Qval.Wald.obj, "Q.sug")
# head(Q.Wald.sug)
# 
# 
# ######## fit ########
# fit.obj1 <- fit(Y, Q, model="GDINA")
# summary(fit.obj1)
# 
# Qval.MLR.obj <- validation(Y, Q, model="GDINA",
#                            mono.constraint = TRUE,
#                            method="MLR-B",
#                            search.method="PAA",
#                            iter.level="test",
#                            maxitr=20, verbose=1)
# Q.MLR.sug <- extract(Qval.MLR.obj, "Q.sug")
# 
# fit.obj2 <- fit(Y, Q.MLR.sug, model="GDINA")
# summary(fit.obj2)
# 
# 
# ######## beta ########
# Qval.beta.obj <- validation(Y, Q, model="GDINA",
#                             mono.constraint = TRUE,
#                             method="beta",
#                             search.method="PAA",
#                             criter="AIC",
#                             iter.level="no", verbose=1)
# Q.beta.sug <- extract(Qval.beta.obj, "Q.sug")
# head(Q.beta.sug)
# 
# ############################## different fit between Q-matrices ##############################
# fit.original.obj <- fit(Y, Q, model="GDINA")
# fit.GDI.obj <- fit(Y, Qval.GDI.obj$Q.sug, model="GDINA")
# fit.Hull.obj <- fit(Y, Q.Hull.sug, model="GDINA")
# fit.Wald.obj <- fit(Y, Q.Wald.sug, model="GDINA")
# fit.MLR.obj <- fit(Y, Q.MLR.sug, model="GDINA")
# fit.beta.obj <- fit(Y, Q.beta.sug, model="GDINA")
# 
# res <- rbind(unlist(fit.original.obj[-c(9, 10)]), unlist(fit.GDI.obj[-c(9, 10)]), 
#              unlist(fit.Hull.obj[-c(9, 10)]), unlist(fit.Wald.obj[-c(9, 10)]), 
#              unlist(fit.MLR.obj[-c(9, 10)]), unlist(fit.beta.obj[-c(9, 10)]))
# print(res)
# 
# 
# ############################# Hull and experts ################################
# summary(Qval.Hull.obj)
# 
# is.Qident(Q.Hull.sug)
# 
# ######## PVAF ########
# PVAF <- get.PVAF(Y, Q.Hull.sug, model="GDINA",
#                  mono.constraint = TRUE)
# PVAF <- round(PVAF, 3)
# head(PVAF)
# 
# ######## R2 ########
# R2 <- get.R2(Y, Q.Hull.sug, model="GDINA", mono.constraint = TRUE)
# R2 <- round(R2, 3)
# R2[c(4, 6), ]
# 
# ######## beta ########
# beta <- get.beta(Y, Q.Hull.sug, model="GDINA", mono.constraint = TRUE)
# beta <- round(beta, 3)
# beta[c(4, 6), ]
# 
# ######## priority ########
# priority <- get.priority(Y, Q.Hull.sug, model="GDINA", mono.constraint = TRUE)
# priority <- round(priority, 3)
# priority[c(4, 6), ]
# 
# ######## Wald test ########
# CDM.obj <- CDM(Y, Q.Hull.sug, model="GDINA", mono.constraint = TRUE)
# q1 <- c(0, 0, 1)
# q2 <- c(1, 0, 1)
# Wald.test(CDM.obj, q1, q2, i=4)
# Wald.test(CDM.obj, q1, q2, i=6)
# 
# ######## att.hierarchy ########
# hiera.obj <- att.hierarchy(Y, Q, model="GDINA", 
#                            mono.constraint = TRUE)
# summary(hiera.obj)
# 
# hiera.Hull.obj <- att.hierarchy(Y, Q.Hull.sug, 
#                                 model="GDINA", 
#                                 mono.constraint = TRUE)
# summary(hiera.Hull.obj)