#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#
# This document contains WinBugs codes for three IRT guessing models
#
# The program was developed by Jing Cao and Lynne Stokes.
#
# For full details of the methods implemented in this program, see the paper:
#
# Cao, J. and Stokes, S.L. Bayesian IRT guessing models for partial guessing behaviors
#
# A copy of this paper is available from the website: http://www.smu.edu/statistics/TechReports/tech-rpts.asp
#
# Copyright Jing Cao and Lynne Stokes 2007
#
# You are free to use this program, for non-commercial purposes only, under two conditions:
#
# (1) This note is not to be removed;
# (2) Publications using this program should reference the manuscript mentioned above
#
# Please email all comments or questions to jcao@smu.edu
#
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 





################################################################################################################
#
# The IRT threshold guessing model (IRT-TG)
#
################################################################################################################


###############################
#
# IRT-TG model specification
#
###############################

model
{ for (i in 1:N){
      for (j in 1:J){
           logit(p[i,j]) <- gamma[j]*(theta[i]-delta[j]) -step(j-alpha[i]+.5)*(gamma[j]*(theta[i]-delta[j]) +1.1)
           x[i,j] ~ dbern(p[i,j])
            } 
      alpha[i] ~ dcat(PI[])
      }

  for(j in 1:J){
     delta[j] ~ dnorm(0,1)
     gamma[j] ~ dgamma(1,1)
     }
  for(i in 1:N){
     theta[i] ~ dnorm(0,1)
     }
  PI[J] ~ dbeta(1,1)
  temp1 <- ( 1.0 - PI[J] ) / pow((J-1), b)  
  PI[1] <- temp1
  for(j in 2:(J-1)){
     PI[j] <- (pow(j, b)-pow((j-1), b))*temp1
     }
  b ~ dgamma(2,1)
 }



###############################
#
# IRT-TG prior assignment
#
###############################
 

list(delta=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),  b=2,
     gamma=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1),
     theta=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
alpha=c(40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
              40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
             40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
              40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
            40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
              40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
           40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
              40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40),
PI=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
        NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,0.5))


###############################
#
# data 
#
###############################


list(N=200, J=40,  
x = structure(.Data = 
c(  0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0,
 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1,
 ......

 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0,
 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1), 
 .Dim = c(200,40)))





################################################################################################################
#
# The IRT difficulty-based guessing model (IRT-DG)
#
################################################################################################################


###############################
#
# IRT-DG model specification
#
###############################


model
{ for (i in 1:N){
      for (j in 1:J){
           logit(p[i,j]) <- gamma[j]*(theta[i]-delta[j]) -(alpha[i]-1)*step(delta[j]-theta[i]-b)*(gamma[j]*(theta[i]-delta[j]) +1.1)
           x[i,j] ~ dbern(p[i,j])
            } 
      alpha[i] ~ dcat(PI[])
      }

  for(j in 1:J){
     delta[j] ~ dnorm(0,1)
     gamma[j] ~ dgamma(1,1)
     }
  for(i in 1:N){
     theta[i] ~ dnorm(0,1)
     }
  PI[1] ~ dbeta(1,1)
  PI[2] <- 1.0 - PI[1]
  b ~ dnorm(0,taub)
  taub ~ dgamma(0.001, 0.001)
 }


###############################
#
# IRT-DG prior assignment
#
###############################
 

list(delta=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), PI=c(0.5, NA), b=0, taub=0.1, 
     gamma=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1),
     theta=c(1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
alpha=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1))





################################################################################################################
#
# The IRT continuous guessing model (IRT-CG)
#
################################################################################################################


###############################
#
# IRT-CG model specification
#
###############################


model
{ for (i in 1:N){
      for (j in 1:J){
           logit(p[i,j]) <-  gamma[j]*(theta[i]-delta[j]) -(alpha[i]-1)*beta[j]
           x[i,j] ~ dbern(p[i,j])
            } 
      alpha[i] ~ dcat(PI[])
      }

  for(j in 1:J){
     delta[j] ~ dnorm(0,1)
     gamma[j] ~ dgamma(1,1)
     }
  for(i in 1:N){
     theta[i] ~ dnorm(0,1)
     }
  PI[1] ~ dbeta(1,1)
  PI[2] <- 1.0 - PI[1]
  beta[1] <- 0
  beta[2] <- 0

  taub6 <- taub*6 
  temp[3] <- -(-4*beta[4]+beta[5])/6 
  temp[4] <- -(-4*beta[3]-4*beta[5]+beta[6])/6 
  beta[3] ~ dnorm( temp[3], taub6)
  beta[4] ~ dnorm( temp[4], taub6) 
   for(j in 5:(J-2)){
     temp[j] <- -(beta[j-2]-4*beta[j-1]-4*beta[j+1]+beta[j+2])/6 
     beta[j] ~ dnorm( temp[j], taub6)
     }
  temp[J-1] <- -(beta[J-3]-4*beta[J-2]-2*beta[J])/5 
  temp[J] <- -(beta[J-2]-2*beta[J-1])
  taub5 <- taub*5
  beta[J-1] ~ dnorm( temp[J-1], taub5)
  beta[J] ~ dnorm( temp[J], taub) 

   taub ~ dgamma(0.001, 0.001)
 }


###############################
#
# IRT-CG prior assignment
#
###############################
 


list(delta=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), PI=c(0.5, NA),  taub=1.0,
     gamma=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1),  
     theta=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
alpha=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1),
beta=c(NA,NA,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0))