####################################################################################################
#    http://www.cs.cas.cz/martinkova/documents/MartinkovaZvara_2007kybernetikaRcode_selected.R     #
#                                                                                                  #
# Selected code for paper:                                                                         #
# Martinkova P & Zvara K (2007): Reliability in the Rasch model, Kybernetika, 43(3),               #
#  pp. 315-326. URL https://dml.cz/bitstream/handle/10338.dmlcz/135776/Kybernetika_43-2007-3_4.pdf #
#                                                                                                  #
# and for IMPS 2013 poster:                                                                        #
# Martinkova P & Turcicova M: Generalized reliability used for the comparison of reliability       #
#  estimates in tests with binary items. IMPS 2013, Arnhem, Netherlands.                           # 
#  URL www.cs.cas.cz/martinkova/documents/posterIMPS2013.pdf                                       #
#                                                                                                  #
# Last edit: July 18, 2018 (Patricia Martinkova, martinkova@cs.cas.cz)                             #
####################################################################################################

###########################################################################################
#                        Calculation of Logistic alpha                                    #
#   as proposed in Eq.10 in the Kybernetika paper by Martinkova & Zvara (2007)            # 
###########################################################################################

logistic.alpha=function(data){
  n=dim(data)[1]	# number of students
  m=dim(data)[2]	# number of items
  success = c(data)	# creates long vector
  item = as.factor(rep(1:m,each=n))
  student = as.factor(rep(1:n,m))
  model1 = glm(success ~ item + student, family = "binomial")
  model2 = glm(success ~ item, family = "binomial")
  Tab = anova(model2, model1, test="Chisq")
  a_Log = 1 - Tab[2,3] / Tab[2,4]
  return(a_Log)
}

####################################################################################
## Definition and calculation of true reliability for Rasch model                 ##
## as proposed in Eq. 13 in the Kybernetika paper by Martinkova and Zvara (2007)  ##
## (newly generalized for 3PL IRT model here                                      ## 
####################################################################################
Rm3=function(b,c,d,sigma){
  J=length(b)
  B=rep(NA,J)
  D=rep(NA,J)
  C=matrix(NA,J,J)
  for (j in 1:J){
    # calculation of Bj terms
    integrand.Bj = function(A) {(exp(-c[j]*(A-b[j])-A^2/(2*sigma^2))+d[j]*exp(-2*c[j]*(A-b[j])-A^2/(2*sigma^2)))*(1-d[j])/((1+exp(-c[j]*(A-b[j])))^2)*(1/sqrt(2*pi*sigma^2))}
    B[j]=integrate(integrand.Bj,-Inf,Inf)$val
    # calculation of Dj terms
    integrand.Dj = function(A) {(exp(-A^2/(2*sigma^2))+d[j]*exp(-c[j]*(A-b[j])-A^2/(2*sigma^2)))/((1+exp(-c[j]*(A-b[j])))*sqrt(2*pi*sigma^2))}
    D[j]=integrate(integrand.Dj,-Inf,Inf)$val
    # calculation of Cjt terms
    for (t in 1:J){
      integrand.Cjt = function(A) {(d[j]+(1-d[j])/(1+exp(-c[j]*(A-b[j]))))*(d[t]+(1-d[t])/((1+exp(-c[t]*(A-b[t])))))*(1/sqrt(2*pi*sigma^2))*exp(-A^2/(2*sigma^2))}
      C[j,t]=integrate(integrand.Cjt,-Inf,Inf)$val
    }
  }
  numerator = sum(C - D%*%t(D))
  denominator = numerator + sum(B)
  Rm = numerator/denominator
  return(Rm)
}
#####################################################################################
############################   Simulation study #####################################
#####################################################################################

library(ltm) # only used for calculation of Cronbach's alpha 
# (can also be done easily without ltm package)

s = c(seq(0.001,1,length=21), seq(1,1.5,length=5), seq(1.5,3,length=6)) #variances
K = length(s)  #number of variance terms
n = 20 # number of students
m = 50 # number of items
nsim = 10#00 # number of simulation runs
b = seq(-0.1,0.1,length=m) # item dificulty parameters
c = rep(1,m) # item discrimination parameters
d = rep(0,m) # probability of correct answer
o = c("Cronbach","logAlpha") # compared estimators of reliability
a = array(dim = c(K, nsim+1, 2))

set.seed(1234)
date()
for (k in 1:K){      			
  for (i in 1:nsim){
    A = rnorm(n,mean=0,sd=s[k])
    p1 = matrix(0,n,m)
    for(l in 1:n){
      for(z in 1:m){
        p1[l,z] = d[z]+(1-d[z])/(1+exp(-c[z]*(A[l]-b[z])))
      }
    }
    Y = matrix(rbinom(n*m,1,prob=p1),ncol=m,nrow=n)
    a[k,i+1,1] = cronbach.alpha(Y)$alpha # Cronbach's alpha calculation
    a[k,i+1,2] = logistic.alpha(Y) # calculation of Logistic alpha
    
  }   
  a[k,1,]=Rm3(b,c,d,s[k]) # True reliability
}
date()

# Simulation for nsim=10 may take few minutes
# Increase nsim to 500 or 1000 to get more precise results

###########################################################################################
###############################  PLOT  ####################################################
###########################################################################################
mean.o1 = apply(a[,2:(nsim+1),1],1,mean)
mean.o2 = apply(a[,2:(nsim+1),2],1,mean)

sd.o1 = apply(a[,2:(nsim+1),1],1,sd)
sd.o2 = apply(a[,2:(nsim+1),2],1,sd)

bias.o1 = mean.o1-a[,1,1]  				
bias.o2 = mean.o2-a[,1,1]

krit = 1.96/sqrt(nsim-1) # confidence interval, as presented in IMPS poster
eps = 0.0015		# horizontal shift for the second estimate and interval
colors = c("black","red") 
# limist for y axis
fromB = min(bias.o1-sd.o1*krit, bias.o2-sd.o2*krit) - 0.01;
toB = max(bias.o1+sd.o1*krit, bias.o2+sd.o2*krit) + 0.01

par(mfrow = c(1,1), mar = c(5.1,4.1,4.1,2.1))
plot(a[,1,1]-eps, bias.o1, axes=F, col=colors[1], pch=16, xlim=c(0,1), ylim=c(fromB,toB), xlab="", ylab="")	
axis(1, tck=-0.01)
axis(2, tck=-0.01)
mtext("reliability", 1, line=3)
mtext("bias", 2, line=2.5)
box()
points(a[,1,1] + eps, bias.o2, col=colors[2], pch=16)
arrows(x0 = a[,1,1]-eps, y0=bias.o1-sd.o1*krit, x1=a[,1,1]-eps, y1=bias.o1+sd.o1*krit, code=3, length=0.04, col=colors[1], angle=90)
arrows(x0 = a[,1,1]+eps, y0=bias.o2-sd.o2*krit, x1=a[,1,1]+eps, y1=bias.o2+sd.o2*krit, code=3, length=0.04, col=colors[2], angle=90)
abline(h=0, col="grey64")	
legend("bottomright", legend=o, col=colors, pch=rep(16,length(colors)))