# posledni uprava: 1. dubna 2014

#########################################################################################
############### Vypocet reliability znalostnich a psychologickych testu #################
#########################################################################################


#########################################################################################
########################## Vice latentnich promennych ###################################
#########################################################################################
library(psych)
data(Tucker)
omega(Tucker,nfactors=2)

#########################################################################################
################# Korelovane chyby polozkovych mereni####################################
#########################################################################################
library(sem)
data(Thurstone)

#zjednoduseni modelu z http://cran.r-project.org/web/packages/sem/sem.pdf
model.thur <- specifyModel()
F1 -> Sentences, lam11
F1 -> Vocabulary, lam21
F1 -> Sent.Completion, lam31
F2 -> First.Letters, lam42
F2 -> 4.Letter.Words, lam52
F2 -> Suffixes, lam62
F3 -> Letter.Series, lam7348 sem
F3 -> Pedigrees, lam83
F3 -> Letter.Group, lam9
F1 <-> F1, NA, 1
F2 <-> F2, NA, 1
F3 <-> F3, NA, 1
F1 <-> F2, cov12
F1 <-> F3, cov13
F2 <-> F3, cov23

##### stiskni Enter

sem.thur <- sem(model.thur, Thurstone, 213)
summary(sem.thur)

# vypocet odhadu korelacni matice dle modelu sem.thur
sigma=matrix(0,9,9)
correlationF1=c(rep(1,3),rep(sem.thur$coeff["cov12"],3),rep(sem.thur$coeff["cov13"],3))
for(i in 1:3){
  for(j in i:9){
    sigma[i,j]=sem.thur$coeff[i]*sem.thur$coeff[j]*correlationF1[j]
  }
}
correlationF2=c(rep(1,3),rep(sem.thur$coeff["cov23"],3))
for(i in 4:6){
  for(j in (i-3):6){
    sigma[i,3+j]=sem.thur$coeff[i]*sem.thur$coeff[3+j]*correlationF2[j]
  }
}
for(i in 7:9){
  for(j in i:9){
    sigma[i,j]=sem.thur$coeff[i]*sem.thur$coeff[j]
  }
}
for (i in 1:9){
  for(j in 1:9){
    if (sigma[i,j]==0){sigma[i,j]=sigma[j,i]}
  }
}
sigmaT=sum(sigma)

sigmaX=sum(Thurstone)
reliability=sigmaT/sigmaX
reliability
#########################################################################################
################# Generalizability theory: vice zdroju chyb #############################
#########################################################################################
library(lme4)
#library(sirt)  # knihovna s daty
#data(data.ratings1)
# data ke stazeni take zde:
data.ratings1<-read.csv2("http://www.cs.cas.cz/martinkova/data.ratings1.csv")
m=5
l=length(levels(data.ratings1$rater))
data.long=reshape(data.ratings1,direction="long",varying=list(3:7),v.names = "ratings", timevar = "items")
#jestli mate pamet mensi nez 6000
#memory.limit(size = NA) 
#pouzijte prosim prikaz
#memory.limit(size = 6000)

(l1<-lmer(ratings~(1|items)+(1|rater)+(1|idstud)+(1|items:rater)+(1|idstud:rater)+(1|idstud:items),data=data.long))

summary(l1)
matrixSD<-as.matrix(lme4::VarCorr(l1)) 
Var<-data.frame(c(matrixSD,attr(matrixSD,"sc")^2)) #rozptyly z lmer
names(Var)=c("idstud:items","idstud:rater","idstud","items:rater",
             "rater","items","Residual")

delta_ik=m*Var["items"]+m^2*Var["rater"]+m*Var["idstud:items"]+m^2*Var["idstud:rater"]+
  +m*Var["items:rater"]+m*Var["Residual"]
(rel=m^2*Var["idstud"]/(m^2*Var["idstud"]+delta_ik))

l=2
delta_i=m*Var["items"]+m^2/l*Var["rater"]+m*Var["idstud:items"]+m^2/l*Var["idstud:rater"]+
  +m/l*Var["items:rater"]+m/l*Var["Residual"]
(rel2=m^2*Var["idstud"]/(m^2*Var["idstud"]+delta_i)) 

l=3
delta_i=m*Var["items"]+m^2/l*Var["rater"]+m*Var["idstud:items"]+m^2/l*Var["idstud:rater"]+
  +m/l*Var["items:rater"]+m/l*Var["Residual"]
(rel2=m^2*Var["idstud"]/(m^2*Var["idstud"]+delta_i)) 

###########################################################################################
##################################Logistic alpha###########################################
###########################################################################################
logistic.alpha=function(data){
  n=dim(data)[1]	# pocet studentu
  m=dim(data)[2]	# pocet polozek
  uspech=c(data)	# vytvori dlouhy vektor
  polozka=as.factor(rep(1:m,each=n))
  student=as.factor(rep(1:n,m))
  model1=glm(uspech~polozka+student,family = "binomial")
  model2=glm(uspech~polozka,family = "binomial")
  Tab=anova(model2,model1,test="Chisq")
  a_Log=1-Tab[2,3]/Tab[2,4]
  return(a_Log)
}
###########################################################################################
################ Vypocet skutecne reliability v 3-parametrickem IRT modelu ################
###########################################################################################
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){
    # vypocet Bj
    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
    # vypocet Dj
    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
    # vypocet Cjt
    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
    }
  }
  citatel=sum(C-D%*%t(D))
  jmenovatel=citatel+sum(B)
  Rm=citatel/jmenovatel
  return(Rm)
}
#####################################################################################
############################   Simulace   ###########################################
#####################################################################################
date()
library(ltm)
s=c(seq(0.001,1,length=21),seq(1,1.5,length=5),seq(1.5,3,length=6))#rozptyly
K=length(s)  #pocet rozptylu
n=20# pocet respondentu
m=50# pocet polozek
nsim=10 #pocet simulaci
b=seq(-0.1,0.1,length=m)#obtiznosti polozek
c=rep(1,m)#diskriminace
d=rep(0,m)#pravdepodobnost uhodnuti spravne odpovedi
o=c("Cronbach","logAlpha") # nazvy odhadu, ktere maji byt porovnany      
a=array(dim = c( K,nsim+1, 2))
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 #vypocet cronbachova alfa
    a[k,i+1,2]=logistic.alpha(Y)#vypocet logistickeho alfa
    
  }   
  a[k,1,]=Rm3(b,c,d,s[k])# skutecna reliabilita
}
date()

# simulace trva priblizne 2 min, pro presnejsi vysledky nutno zvysit nsim

###########################################################################################
###############################  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)#interval
eps = 0.0015		# bocni odstup prumeru pro stejnou hodnotu na vodorovné ose
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)))