# load the package and the data
library("equateMultiple")
data("mathTest", package = "equateMultiple")


# The first 5 rows the the first data frame in object mathTest
head(mathTest[[1]])


# estimation of the item parameters with the mirt package
library("mirt")
mods_mirt <- list()
for (i in 1:5)
  mods_mirt[[i]] <- mirt(mathTest[[i]], 1, itemtype = "2PL", SE = TRUE)


# conversion the item parameters and their covariance matrix to the parameterization in Equation (1)
mods_extract <- modIRT(mods_mirt, display = FALSE)


# number of items in common between pairs of forms
linkage_plan <- linkp(mods_extract)
linkage_plan


# a representation of the linkage plan
library("sna")
set.seed(6)
par(mar = c(0, 0, 0, 0))
gplot(linkage_plan, displaylabels = TRUE, vertex.sides = 4, 
      vertex.cex = 3, vertex.rot = 45, usearrows = FALSE, 
      label.pos = 5, label.cex = 1, vertex.col = 0, edge.lwd = 0.2)


# estimation of the equating coefficients with the mean-mean method
eq_mm <- multiec(mods_extract)
summary(eq_mm)


# estimation of the equating coefficients with the multiple item response function method
eq_irf <- multiec(mods_extract, method = "irf", base = 5)
summary(eq_irf)


# estimation of the equating coefficients with the multiple test response function method
eq_lik <- multiec(mods_extract, method = "lik", start = eq_mm, obsinf = FALSE)
summary(eq_lik)


# item parameters expressed on a common metric and their standard errors
items.com <- item.common(eq_lik)
items.com[1:4, ]
items.com[42:45,]


# equated scores obtained with true score equating
sc.eq.tse <- score(eq_lik)
round(sc.eq.tse, 3)[1:6,]


# equated scores obtained with observed score equating
sc.eq.ose <- score(eq_lik, method = "OSE")
round(sc.eq.ose, 3)[1:7,]



# -----------------------------------------------------------------------------
# simulated data

# true equating coefficients
A <- seq(1, 2, length = 5)
B <- seq(0, 2, length = 5)
A
B


# generation of ability values
set.seed(1)
n <- 100000 # sample size
theta <- list()
for (i in 1:5) theta[[i]] <- rnorm(n, B[i], A[i])


# generation of true item parameters
set.seed(1)
as <- runif(20, 0.7, 1.3)
bs <- sort(rnorm(20, 1, 1))


# function for generation of item responses
gen.resp <- function(theta, a, b) {
  lp <- a * (theta - b)
  pr <- plogis(lp)
  rn <- runif(length(theta))
  (pr > rn) * 1
}

# item labels
itms <- paste("I", formatC(1:20, width = 2, format = "d", flag = "0"), 
              sep = "")

# generation of item responses
set.seed(1)
resp <- list()
for (i in 1:5) {
  resp_i <- matrix(NA, n, 20)
  for (j in 1:20) resp_i[, j] <- gen.resp(theta[[i]], as[j], bs[j])
  colnames(resp_i) <- itms
  resp[[i]] <- resp_i
}


# histogrammes of row scores
row.scores <- lapply(resp, rowSums)
par(mfrow = c(1, 5), mar = c(2, 2, 2, 1))
for (i in 1:5) hist(row.scores[[i]], main = paste("Form", i), 
                    xlab = "", col = 5)


# estimation of item parameters
mods_mirt_sim <- list()
for (i in 1:5)
  mods_mirt_sim[[i]] <- mirt(resp[[i]], 1, itemtype = "2PL", SE = TRUE)


# conversion of item parameters
mods_extract_sim <- modIRT(mods_mirt_sim, display = FALSE)


# estimation of the equating coefficients
eq_irf_sim <- multiec(mods_extract_sim, method = "irf", se = TRUE)
summary(eq_irf_sim)


# plot of item parameter estimates before and after conversion
par(mfrow=c(2, 2))
plot(eq_irf_sim, form = "T5")


# equated scores
sc.eq.tse.sim <- score(eq_irf_sim)
round(sc.eq.tse.sim, 3)[1:6,]


# deletion of half of the items in Forms 1 and 5
resp[[1]] <- resp[[1]][, 1:10]
resp[[5]] <- resp[[5]][, 11:20]


# estimation of the item parameters
for (i in c(1, 5))
  mods_mirt_sim[[i]] <- mirt(resp[[i]], 1, itemtype = "2PL", SE = TRUE)
mods_extract_sim <- modIRT(mods_mirt_sim, display = FALSE)


# estimation of the equating coefficients
eq_irf_sim <- multiec(mods_extract_sim, method = "irf", se = TRUE)
summary(eq_irf_sim)


# equated scores
sc.eq.tse.sim <- score(eq_irf_sim)
round(sc.eq.tse.sim, 3)