##############################################################
# Edps/Psych/Stat 587
# Fall 2020
# c.j.anderson
#
# Using GSS 2004 vocabulary items show how to fit Rasch model 
#   with explanatory variables for theta
#
# Note 2pl is a better model for these data, but this is to
#  show how to fit rasch model with expanatory varibles
#
# Packages
#   lme4
#   nlme
############################################################

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

attach(vocb)

# Rasch via laplace
rasch1 <- glmer(Y ~ -1 + A + B + C + D + E + F + G + H + I + J + (1 | id),data=vocab, family=binomial)
summary(rasch1)


# Rasch via guass newton quadrature
rasch2 <- glmer(Y ~ -1 + A + B + C + D + E + F + G + H + I + J + (1 | id),data=vocab, family=binomial,nAGQ=10)
summary(rasch2)

# Rasch via nlme

# Basic
onePL <- function(b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,theta) {
          b=b1*A + b2*B + b3*C + b4*D + b5*E + b6*F + b7*G + b8*H + b9*I + b10*J 
          exp(theta-b)/( 1 + exp(theta-b) )  
}

rasch3 <- nlme(Y ~ onePL(b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,theta),
              data=vocab,
              fixed= b1+b2+b3+b4+b5+b6+b7+b8+b9+b10 ~ 1,
              random = theta ~ 1 |id,
              start=c(b1=1, b2=1,b3=1, b4=1, b5=1,b6=1, b7=1, b8=1, b9=1, b10=1) )
summary(rasch3)


# Rasch with explanatory variables for theta
onePL <- function(b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,v1,v2,v3,epsilon) {
          b=b1*A + b2*B + b3*C + b4*D + b5*E + b6*F + b7*G + b8*H + b9*I + b10*J 
          theta = v1*age + v2*elementary + v3*hsplus + epsilon
          exp(theta-b)/( 1 + exp(theta-b) )  
}

rasch4 <- nlme(Y ~ onePL(b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,v1,v2,v3,epsilon),
              data=vocab,
              fixed= b1 +b2 +b3 +b4 +b5 +b6 +b7 +b8+b9 +b10 +v1 +v2 +v3 ~ 1,
              random = epsilon ~ 1 |id,
              start=c(b1=1, b2=1,b3=1, b4=1, b5=1,b6=1, b7=1, b8=1, b9=1, b10=1, v1=.1, v2=.1, v3=.1 ) )
summary(rasch4)