###################################################
#
# Edps/Psych/Stat 587
# Fall 2020
# C.J.Anderson
#
#  LSAT6 -- Rasch & why 2pl can't be fit by nlme
#
#  Need packages:
#      lme4
#      nlme
#
#  Models fit to data in individual level format
#
###################################################

lsat <- read.table(file="C:/Users/cja/Documents/Dropbox/edps587/lectures/10 logreg/lsat6_long_data.txt", header=TRUE)
head(lsat,25)


###############################################
# glmer and Rasch 
###############################################
# Default is LaPlace

rasch.laplace <- glmer(y ~ -1 + i1 + i2 + i3 +i4 + i5 + (1 | id),data=lsat, family=binomial)
summary(rasch.laplace)

# Guass quadrature -- MLE gold standard

rasch.quad <- glmer(y ~ -1 + i1 + i2 + i3 +i4 + i5 + (1 | id),data=lsat,family=binomial, nAGQ=10 )
summary(rasch.quad)

# What is modeled above is P(Y=0), to get P(Y=1) either
#  (a) flip the signs of the estimated parameters
#  (b) recode
#

lsat$ycorrect <- ifelse (lsat$y==0, 1, 0)

rasch2.quad <- glmer(ycorrect ~ -1 + i1 + i2 + i3 +i4 + i5 + (1 | id),data=lsat,family=binomial, nAGQ=10 )
summary(rasch2.quad)



# Look at curves:
x <- seq(from=0, to=8, by=0.1)
theta<- x-4
item1 <-  exp(theta +2.73001)/(1 +exp(theta +2.73001))
item2 <-  exp(theta +0.99862)/(1 +exp(theta +0.99862))
item3 <-  exp(theta +0.23986)/(1 +exp(theta +0.23986))
item4 <-  exp(theta +1.30646)/(1 +exp(theta +1.30646))
item5 <-  exp(theta +2.09940)/(1 +exp(theta +2.09940))


plot(theta,item1, type='l', 
     col='blue', ylim=c(0,1),
     ylab='Probability Correct', 
     main='Item Characteristic Curves: Rasch (lsat6)')
lines(theta,item2,col='red')
lines(theta,item3,col='green')
lines(theta,item4,col='black')
lines(theta,item5,color='cyan')
abline(h=.5)



############################################
# nlme and 1pl  (rasch) 
#
#
############################################

attach(lsat)

# model to fit
onePL <- function(b1,b2,b3,b4,b5,theta) {
           b=b1*i1 + b2*i2 + b3*i3 + b4*i4 + b5*i5
          exp(theta-b)/( 1 + exp(theta-b) )  
}

rasch3 <- nlme(y ~ onePL(b1,b2,b3,b4,b5,theta),
              data=lsat,
              fixed= b1+b2+b3+b4+b5 ~ 1,
              random = theta ~ 1 |id,
              start=c(b1=1, b2=1,b3=1, b4=1, b5=1) )
summary(rasch3)

#######################################################
# nlme 2pl 
#
#  Model is not identified
#     o Need to fix variance to 1 or someother constraint.
#     o This does not seem possible, which
#        (and some other tests) leads me to think
#         that...
#     o I do not seem to be able to fit any model
#         that multiplies the random effect
#
#  Use SAS NLMIXED or other program where can fix variance
#
#######################################################

twoPL <- function(b1,b2,b3,b4,b5,a1,a2,a3,a4,a5,theta) {
         a = a1*i1 + a2*i2 + a3*i3 + a4*i4 + a5*i5
         b = b1*i1 + b2*i2 + b3*i3 + b4*i4 + b5*i5
         exp(a*theta-b)/( 1 + exp(a*theta-b) )  
}


two.pl <- nlme(y ~ twoPL(b1,b2,b3,b4,b5,a1,a2,a3,a4,a5,theta),
          data=lsat,
          fixed= b1+b2+b3+b4+b5+ a1+a2+a3+a4+a5~ 1 ,
          random= theta ~ 1 | id,
          weight=varFixed( )
          start=c(b1=1, b2=1, b3=1, b4=1, b5=1, a1=1, a2=1, a3=1, a4=1, a5=1)
       )

###########################################################
# In R, use Bayesian estimation of use IRT software 
#    (e.g, mirt package)
#
# Switch to SAS
###########################################################