#PRIMEIROS PASSOS NA TEORIA DA RESPOSTA AO ITEM
###############################################
Prof. Jorge Luis Bazán



#Gerando dados de um modelo logístico de 2 parâmetros
#####################################################

#função para gerar o modelo logístico de 2 parâmetros

gera2L<-function(thetasim,itemsim)
{
n<-length(thetasim)      #leitura do número de pessoas 
k<-nrow(itemsim)         #leitura do número de itens 
set.seed(245)            #com isto os dados gerados serão os mesmos  
thetas<-thetasim  # parâmetros das habilidades simuladas
items<-itemsim    # parâmetros dos itens simulados
names(items)<-c("b","a") 

# Cálculo das probabilidades de acerto e matriz de respostas 
f2L<-function(x) {plogis(t(x[2]%*%( t(thetas-x[1]) ) ) ) } 
p2L<-apply(items,1,f2L) 
u<-matrix(runif(n*k),n)
data2L<-ifelse(u<p2L,1,0)
return(data2L)
}


#Especificação dos valores dos dados

n<-1000  # Número de pessoas
k<-20    # Número de itens

#Especificação dos parâmetros verdadeiros a serem simulados

itemsim<-data.frame(cbind(runif(k,-2,2),runif(k,0.4,2))) #items
thetasim<-matrix(rnorm(n),ncol=1)  # habilidades 



#Gerando os dados usando a função para os valores simulados

data2L<-gera2L(thetasim,itemsim)
#visualizar os dados gerados
#data2L
#edit(data2L)

#salvando os dados e os parâmetros verdadeiros em arquivos
write.table(data2L,"C:/Users/Jorge Luis/Dropbox/MinicursoTRI/MinicursoVerão2017/Dia2/data2L.csv",sep=",") 
write.table(thetasim,"C:/Users/Jorge Luis/Dropbox/MinicursoTRI/MinicursoVerão2017/Dia2/Dia4pessoasim2L.csv",sep=",") 
write.table(itemsim,"C:/Users/Jorge Luis/Dropbox/MinicursoTRI/MinicursoVerão2017/Dia2/itensim2L.csv",sep=",") 

#definindo os dados como data frame
data2L=data.frame(data2L)      



#Ajustando um modelo 1L (Rasch), um modelo 2L e um modelo 3L (Calibração)
#usando estimação de máxima verossimilhança marginal 
#########################################################################  

require(ltm)
#Ajusta 2L 
fit2L<-ltm(data2L ~ z1,IRT.param = TRUE,control = list(GHk = 50) )
fit2L
fit2L$Coefficients


#Ajusta 1L(Rasch) 
fitrasch <- rasch(data2L, constraint = cbind(ncol(data2L) + 1, 1),control = list(GHk = 50) )
fitrasch
anova(fitrasch,fit2L)



#ML Estimation of Logistic IRT unidimensional models using ltm: Calibration of items
#with usual parameterization, models 1L, 2L and 3L
#################################################################################
#Using irtoys


require(irtoys)

#Estimation One Logistic IRT Model 1L
fit1La<- est(data2L, model="1PL", engine="ltm")

colnames(fit1La$est) <-c("a","b","c")
fit1La$est

#Estimation Two Logistic IRT Model 2L
fit2La<- est(data2L, model="2PL", engine="ltm")
colnames(fit2La$est) <-c("a","b","c")
fit2La$est

#Estimation Two Logistic IRT Model 2L
fit3La<- est(data2L, model="3PL", engine="ltm")
colnames(fit3La$est) <-c("a","b","c")
fit3La$est

#Comparison of estimation of difficults parameters under differents IRT models and Classical difficulty
##############################################################################
require(psychometric)

#classical difficulty (proportion of correct response)
score=apply(data2L,1,sum)
itemanalisis<-item.exam(data2L, y=score, discrim=TRUE)
p<-itemanalisis$Difficulty

difficulties=matrix(c(
1-p,
log((1-p)/p),
fit1La$est[,2],
fit2La$est[,2],
fit3La$est[,2]),
nrow=ncol(data2L))
colnames(difficulties) <- c("prop","prop tran","b 1L","b 2L","b 3L")
difficulties

#difuldades na escala ENEM
difficultiesest=500+10*difficulties[,-1]
difficultiesest

cor(difficultiesest)
par(mfrow=c(1,1))
matplot(difficultiesest)
#plot.matrix(difficultiesest)

#Estimation of abilities
#########################
detach("package:psych", unload=TRUE)

#2L using ltm
EAP=factor.scores(fit2L, method = "EAP")$score.dat
write.table(EAP,"C:/Users/Jorge Luis/Dropbox/MinicursoTRI/MinicursoVerão2017/Dia2/EAPmathpear.csv",sep=",") 

#2L using irtoys considering different methods

#Maximum Likelihood estimation
th.mle2L <- mlebme(resp=data2L, ip=fit2La$est, method="ML")
th.mle2L

#Expected a Posteriori
th.eap2L <- eap(resp=data2L, ip=fit2La$est, qu=normal.qu())
th.eap2L 

#Mode a Posteriori
th.map2L <- mlebme(resp=data2L, ip=fit2La$est, method="BM")
th.map2L

#Weight Likelihood estimation
th.wle2L <- wle(resp=data2L, ip=fit2La$est)
th.wle2L

#Comparison of estimation of abilities under differents methods and 
# Classical Score
###########################################################################


thetas=matrix(c(
score,
scale(score),
th.mle2L[,1],
th.eap2L[,1],
th.map2L[,1],
th.wle2L[,1]),
nrow=nrow(data2L))
colnames(thetas) <- c("Score","Z score","MLE","EAP","MAP","WLE")
thetas

#habilidades na escala ENEM

thetast=500+10*thetas[,-1]
thetast
cor(thetast)
matplot(thetast)

plot(thetast[,1],thetast[,3],xlim=c(470,525),ylim=c(470,525),xlab="Z score",ylab="EAP")
lines(thetast[,1],thetast[,1])


# Note. Reviewed your results with carefully  Compare your results with the obtained 
# using other software. For example you can use IRTPRO
# http://www.ssicentral.com/irt/student.html



#COMPLEMENTAR: Analise Classico dos Testes
##########################################

require(psych)
alpha(data2L)


# Avaliação da dimensionalidade
###############################
#Analise fatorial de variáveis dicotomicas

#calculando correlação tetracórica
require(psych)
c=tetrachoric(data2L)

fap1 <- fa(r=c$rho)  #factor analysis using tetrachoric correlation matrix
fap1
fap2 <- fa(r=c$rho,nfactors=2,rotate="varimax")
fap2

par(mfrow=c(1,2))
fa.diagram(fap2,cut=0.4)
factor.plot(fap2, cut=0.4,ylim=c(0,1),xlim=c(0,1))
abline(h=0.4,col=4,lty=2)
abline(v=0.4,col=4,lty=2)