#############################################################################
# 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)