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

```