####################################################################################
#                         L7 - IRT MODELS FOR BINARY DATA                          #
#                                    NMST570                                       #
#                             Last update: 16/11/2018                              #
####################################################################################

# Slides available at http://www.cs.cas.cz/martinkova/NMST570.html
# HW assignment available at http://www.cs.cas.cz/drabinova/documents/NMST570_HW7.pdf

####################################################################################
# Loading data #####################################################################
####################################################################################

data <- read.csv("epi_escore.csv")
summary(data)

####################################################################################
# Loading packages #################################################################
####################################################################################

library(mirt)
library(ltm)

####################################################################################
# Ex. 2.1. #########################################################################
####################################################################################

#-----------------------------------------------------------------------------------
# fitting Rasch model with ltm package

### open help for function rasch()
?rasch

### by default, rasch() function fit 1PL model, you need to use argument constraints
### to fit Rasch model
### replace ???
fitR.ltm <- rasch(data = ???, constraint = cbind(ncol(data) + 1, 1))

### print estimated parameters, you can use coef() function on model

#-----------------------------------------------------------------------------------
# fitting Rasch model with mirt package

### open help for function mirt()
?mirt

### fit Rasch model with mirt() function. You need to specify data, model and itemtype
### In our case model = 1 means that we assumed that latent trait is unidimensional
### itemtype = "Rasch" means that we assumed Rasch model for all items
### replace ???
fitR.mirt <- mirt(data = ???, model = ???, itemtype = ???)

### print estimated parameters, you can use coef() function on model
### use simplify = T argument
### for more options open help for coef method
??coef.SingleGroupClass

####################################################################################
# Ex. 2.2. #########################################################################
####################################################################################

#-----------------------------------------------------------------------------------
# fitting 2PL model with ltm package

### open help for function ltm()
?ltm

### you have to introduce formula for latent trait
fit2PL.ltm <- ltm(data ~ z1)

### print estimated parameters, use coef() function on model


#-----------------------------------------------------------------------------------
# fitting 2PL model with mirt package

### fitting 2PL IRT model model with mirt() function is the same as fitting Rasch 
### model, only difference is in itemtype = "2PL"
### replace ???
fit2PL.mirt <- mirt(data = ???, model = ???, itemtype = ???)

### print estimated parameters, use coef() function on model
### use simplify = T argument

####################################################################################
# Ex. 2.3. #########################################################################
####################################################################################

#-----------------------------------------------------------------------------------
# 1. reversing some items

### this is function to reverse some items of dataset
reverse.items <- function(data, which.items){
  rev.data <- sapply(data, function(x) as.numeric(x == 0))
  new.data <- data
  new.data[, which.items] <- rev.data[, which.items]
  return(new.data)
}

### specify which items do you want to reverse (their colnames)
### replace ???
colnames(data)
which.items <- c(???)

### now create new data with reverse.items() function
### replace ???
newdata <- reverse.items(data = ???, which.items = ???)

#-----------------------------------------------------------------------------------
# 2. fitting 2PL model with mirt package on data with reversed items

### replace ??? 
fit2PL.mirt.rev <- mirt(data = ???, model = ???, itemtype = ???)

#-----------------------------------------------------------------------------------
# 3. print estimated parameters

### use coef() function on model
### use simplify = T argument

#-----------------------------------------------------------------------------------
# 4. extract factor scores

### open help for fscores() function
?fscores

### use fscores() to extract factor scores
###### 2PL IRT model on original data, replace ???
fs2PL <- as.vector(fscores(???))

###### 2PL IRT model on data with reversed items, replace ???
fs2PL.rev <- as.vector(fscores(???))

### compare factor scores, you can for example print summary on both, plot them using 
### function plot() or calculate their correlation with corr() function


#-----------------------------------------------------------------------------------
# 5. unidimensionality of latent trait on data with reversed items

### we have actually already fitted unidimensional 2PL IRT model, this was specified 
### with model = 1 argument in mirt() function
fit2PL.mirt.rev
summary(fit2PL.mirt.rev)
plot(fit2PL.mirt.rev)

### we need to fit 2PL IRT model with 2 factors, i.e., with model = 2
fit2PL.mirt.rev2 <- mirt(data = ???, model = ???, itemtype = ???)
summary(fit2PL.mirt.rev2)
plot(fit2PL.mirt.rev2)

### help for plot
??plot.SingleGroupClass

### help for anova
??anova.SingleGroupClass

### now compare both models with anova() function
anova(fit2PL.mirt.rev, fit2PL.mirt.rev2)

#-----------------------------------------------------------------------------------
# 6. item-fit of models on data with reversed items

### open help for itemfit() function
?itemfit

### item-fit for unidimensional model, replace ???
itemfit(???)

### item-fit for model with 2 factors, replace ???
itemfit(???)

#-----------------------------------------------------------------------------------
# 7. estimation of parameters in mirt()

### open help for mirt() function
?mirt

### fit 2PL IRT models using various estimation algorithms:
###### EM, replace ???
fit2PL.mirt.EM <- mirt(data = ???, model = ???, itemtype = ???, method = ???)

###### QMCEM, replace ???
fit2PL.mirt.EM <- mirt(data = ???, model = ???, itemtype = ???, method = ???)

### print estimated parameters, use coef() function on model
### use simplify = T argument


#-----------------------------------------------------------------------------------
# 8. different distributions of latent trait

### fit 2PL IRT models using different distributions of the latent trait
###### Gaussian, replace ???
fit2PL.mirt.Gauss <- mirt(data = ???, model = ???, itemtype = ???, dentype = ???)

###### empirical histogram, replace ???
fit2PL.mirt.EH <- mirt(data = ???, model = ???, itemtype = ???, dentype = ???)

### print estimated parameters, use coef() function on model
### use simplify = T argument


### use fscores() to extract factor scores


### draw histograms of factor scores