#############################################################################
# eirm_c2_summary.R
#   adaptation of IRT analyses for the book "Explanatory Item
#   Response Models" (De Boeck and Wilson, Springer 2004)
#
#   !!! This file has to be called after eirm_c2.R !!!
#
# -------------------------------------------
# Christophe Lalanne
# Centre international d'études pédagogiques
# lalanne@ciep.fr
# source code avalaible on: www.aliquote.org
# -------------------------------------------
#
# Time-stamp: Time-stamp: <2008-02-26 16:37:52 chl>
#############################################################################


####
#### Summary of the results
####

# The four models and their estimates (beta + theta) are summarized in the
# following figures, with the help of item-person maps.

get.mid.interval <- function(x) {
        #
        # *input*
        # a vector computed using R cut() function
        # *output*
        # return a vector containing the corresponding mid positions 
        #
        # Note: not really clean coding...
        #
        out <- NULL
        for (i in 1:length(x)) {
                comma.pos <- regexpr(",",as.character(x[i]))
                end <- regexpr("]",as.character(x[i]))
                xmin <- as.numeric(substr(as.character(x[i]),2,comma.pos[1]-1))
                xmax <- as.numeric(substr(as.character(x[i]),comma.pos[1]+1,end[1]-1))
                out <- append(out,(xmin+xmax)/2)
        }
        return(out)
}

ipm <- function(x, y, k=8, n=10, ...) {
        #
        # *input*
        # x = beta estimates (difficulty)
        # y = theta estimates (ability)
        # k = number of classes to consider
        # n = number of person per subclass
        # *output*
        # graphical side-effect
        #
        dx <- .1
        op <- par(mar=c(1,1,3,1),lab=c(1,1,2))
        plot(c(-.5,.5),c(-5,5),type="n",xlab="",ylab="",axes=FALSE, ...)
        axis(4,at=-5:5,labels=FALSE,pos=0)
        axis(4,at=-4.5:4.5,labels=FALSE,tcl=-.25,pos=0)
        text(rep(-dx/3,11),-5:5,-5:5,srt=90,cex=.8)
        text(c(-.25,.25),c(5,5),c("Person","Item"),font=2)
        text(c(-.25,.25),c(4,4),c(paste("(N=",length(y),")",sep=""),paste("(k=",length(x),")",sep="")))
        abline(v=c(-dx,dx),col="lightgrey")
        h.y <- hist(y,breaks=k,prob=TRUE,plot=FALSE)
        x.breaks <- h.y$breaks
        if (min(x) < h.y$breaks[1]) x.breaks <- c(min(x),x.breaks)
        if (max(x) > h.y$breaks[length(h.y$breaks)]) x.breaks <- c(x.breaks,max(x))     
        h.x <- hist(x,breaks=x.breaks,prob=TRUE,plot=FALSE)
        ## xd <- density(x); lines(dx+xd$y,xd$x)
        yd <- density(y); lines(-dx-yd$y,yd$x)
        x.idx <- cut(x,breaks=x.breaks)
        y.idx <- cut(y,breaks=h.y$breaks)
        x.tab <- table(get.mid.interval(x.idx))
        y.tab <- ceiling(table(get.mid.interval(y.idx))/n)
        dz <- max(h.y$density)/max(y.tab)
        j <- 0
        for (i in 1:length(y.tab)) {        # person 
                if (y.tab[i]>1) {
                        ## xx <- -dx+seq(0,-h.y$density[i],length=y.tab[i])
                        xx <- -dx+seq(0,-h.y$density[i],by=-dz)
                        yy <- rep(as.numeric(names(y.tab[i])),y.tab[i])
                  text(xx[-length(xx)],yy,"x",cex=.8)
                }
                else
                text(-dx,as.numeric(names(y.tab[i])),"x",cex=.8)                  
        }
        j <- 0
        for (i in 1:length(table(x.idx))) { # item
          if (table(x.idx)[i]!=0) {
                        j <- j+1
                        #text(seq(dx,h.x$density[i]+dx,length=x.tab[j]),rep(as.numeric(names(x.tab[j])),x.tab[j]),names(xx[which(as.numeric(x.idx)==i)]),cex=.8,pos=4)
                        if (x.tab[j]>1) {
                                xx <- dx+seq(0,x.tab[j]*.5/10,by=.5/10)
                                yy <- rep(as.numeric(names(x.tab[j])),x.tab[j])
                          text(xx[-length(xx)],yy,names(x[which(as.numeric(x.idx)==i)]),cex=.8,pos=4)
                        }
                        else
                        text(dx,as.numeric(names(x.tab[j])),names(x[which(as.numeric(x.idx)==i)]),cex=.8,pos=4)                 
                }
        }
        # # item statistics
        # segments(.1,-4.5,.55,-4.5)
        # text(.1,-4.7,paste(names(summary(x)),collapse=" "),cex=.8,pos=4)
        # text(.1,-5,paste(round(summary(x),2),collapse=" "),cex=.8,pos=4)
        # # person statistics
        # segments(-.1,-4.5,-.55,-4.5)
        # text(-.55,-4.7,paste(names(summary(x)),collapse=" "),cex=.8,pos=4)
        # text(-.55,-5,paste(round(summary(y),2),collapse=" "),cex=.8,pos=4)
        # par(op)
}

# Rasch (ltm)
xx1 <- coef(a.1.rasch.1)[,1]
names(xx1) <- 1:24
yy1 <- factor.scores(a.1.rasch.1)$score.dat[,"z1"]
# Rasch (lmer)
xx2 <- b.lme.diff
names(xx2) <- 1:24
yy2 <- b.lme.abil
# LRRM
xx3 <- c.lme.diff[4:26]
names(xx3) <- 2:24
yy3 <- c.lme.abil
# LLTM
xx4 <- a.1.lltm.diff
names(xx4) <- 1:24
yy4 <- a.1.lltm.abil

op <- par(mfrow=c(2,2))
ipm(xx1,yy1,main="Rasch (ltm)")
ipm(xx2,yy2,main="Rasch (lmer)")
ipm(xx3,yy3,main="LRRM (lmer)")
ipm(xx4,yy4,main="LLTM (eRm)")
par(op)