######################################################################
##  Copyright 2022 Sébastien Béland & Carl F. Falk
##
##  This program is free software: you can redistribute it and/or
##    modify it under the terms of the GNU General Public License as
##    published by the Free Software Foundation, either version 3 of
##    the License, or (at your option) any later version.
##
##    This program is distributed in the hope that it will be useful,
##    but WITHOUT ANY WARRANTY; without even the implied warranty of
##    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##    GNU General Public License for more details.
##    <http://www.gnu.org/licenses/>
######################################################################
##
## This file simulates a dichotomous item dataset based on the 2PLM
## then computes its true reliability and several estimates

source("ReliabilityFunctions.R")

require(cacIRT) # for generating fake data
require(lavaan) # for fitting linear factor analysis models and Green & Yang (2009)
require(ltm) # used for calibration for Dimitrov's model-based reliability approach
require(mirt) # alternative estimation package for Dimitrov's model-based reliability
require(sirt) # contains functions for Green & Yang (2009) approach
require(psych) # glb, alpha
require(Rcsdp) # necessary for glb


##---------------------------------------------------------------------------------------
## Simulate data & compute true reliability

set.seed(9887)

# number of items and subjects
nItems <-45
nSubjs <- 1000

# item parameters
a     <- rlnorm(nItems,-.75,.25)
b     <- rnorm(nItems,0,.75)
c     <- rep(0,nItems) # guessing

abil  <- rnorm(nSubjs)
pa    <- cbind(a,b,c)
Xfull <- cacIRT:::sim(pa, abil)

# True reliability
true.rel <- DQ(a,b, nItems)


##---------------------------------------------------------------------------------------
## Estimate reliability using all approaches in paper

# Dimitrov Analytic approach, using true item parameters
# Just to know how much bias it has even if true parameters are known
DATruepar <- DA(a,b,nItems)

# Estimation of 2PL parameters from Xfull
irt.mtf <- ltm(Xfull ~ z1, IRT.param = TRUE)
Param   <-coef(irt.mtf)
a2PL <- Param[,2]/1.7 # scaling constant necessary here
b2PL <- Param[,1]

# Dimitrov Quadrature approach, using estimated parameters
DQest <- DQ(a2PL, b2PL)

# Dimitrov Analytic approach, using estimated parameters
DAest <- DA(a2PL, b2PL, nItems)

# alpha
alphaest <- alpha(Xfull)$total$raw_alpha

# glb
glbest <- glb.algebraic(Xfull)$glb

# Linear factor analysis (CFA) approach using lavaan
CFAest <- modelrel(Xfull)

# Green & Yang (2009) approach using both lavaan and sirt
GYest <- greenyang2009(Xfull,estimator="WLSMV")

# All estimates, including true reliability as first entry
data.frame(true.rel,DATruepar,DQest,DAest,alphaest,glbest,CFAest,GYest)