`plotPWmap` <-
function(object, pmap=FALSE, imap=TRUE, item.subset="all", person.subset="all",
                 mainitem="Item Map", mainperson="Person Map",
                 mainboth="Item/Person Map", latdim="Latent Dimension",
                 tlab="Infit t statistic", pp=NULL, cex.gen=0.6, cex.pch=1,
                 person.pch=1, item.pch=16, personCI=NULL, itemCI=NULL, horiz=FALSE)
{
  def.par <- par(no.readonly = TRUE) ## save default, for resetting...

  ## Pathway map currently only for RM, PCM and RSM

  ## The next part of the code finds locations and standard errors for
  ## the item thresholds
  if ((object$model == "LLTM") || (object$model == "LRSM") || (object$model == "LPCM"))
    stop("Pathway Map can only be computed for RM, RSM, and PCM!")

  if (!pmap && !imap)
    stop("Pathway Map requires you to request at least one map (item or person)!")

  ## compute threshtable (from betapars for dichotomous models) and item names
  if (object$model == "RM" || max(object$X, na.rm=TRUE) < 2 ) { # dichotomous model
    dRm <- TRUE

    ## betapars are easiness parameters; only the pars need negating
    threshtable<-cbind(object$betapar * -1, object$se.beta)
    rownames(threshtable) <- colnames(object$X)

    ## shorter synonym
    tt<-threshtable
  } else { ## polytomous model
    dRm <- FALSE

    thresh <- thresholds(object)
    threshtable <- cbind(thresh$threshpar, thresh$se.thresh)
    tlevels<-apply(thresh$threshtable[[1]], 1,
                   function(x) length(na.exclude(x))) - 1
    if (!(sum(tlevels)==nrow(threshtable)))
      stop("Threshtable rows not equal to number of thresholds - oops!")

    ttl<-NULL ## threshtable labels
    for (i in rownames(as.matrix(tlevels)))
      if (tlevels[i]==1)
        ttl<-c(ttl,i)
      else
        ttl<-c(ttl,paste(i,1:tlevels[i],sep=":"))
    rownames(threshtable)<-ttl

    ## shorter synonyms
    tt<-threshtable
    tl<-tlevels
  }

  if (is.null(pp))
    suppressWarnings(pp<-person.parameter(object))
  else if (class(pp) != "ppar" || !identical(pp$X,object$X))
    stop("pp is not a person.parameter object which matches the main Rasch data object!")


  ## We will be plotting the infit data versus the parameters for
  ## both items and persons
  iloc<-tt[,1]
  ise<-tt[,2]
  ifit <- itemfit(pp)
  ifitZ <- ifit$i.infitZ

  ploc <- as.matrix(pp$theta.table['Person Parameter'])[,1]
  pse <- unlist(pp$se.theta, recursive=FALSE)
  names(pse) <- sub("^NAgroup[0-9]*\\.","",names(pse))
  pse <- pse[names(ploc)]
  pfit <- personfit(pp)
  pfitZ <- pfit$p.infitZ

  ## We can now do item and person subsetting; the item subsetting is
  ## pretty ugly as there are multiple cases.  (We dare not do it earlier
  ## as we have to take items from all of iloc, ise and ifitZ.)
  if (imap && is.character(item.subset)) {
    ## Case 1: item subsetting by item names
    if (dRm) {
      if (length(item.subset)>1 && all(item.subset %in% rownames(tt))) {
        iloc  <- iloc[item.subset]
        ise   <- ise[item.subset]
        ifitZ <- ifitZ[item.subset]
        tt    <- tt[item.subset,]
      }
      else if(length(item.subset)!=1 || !(item.subset=="all"))
        stop("item.subset misspecified. Use 'all' or vector of at least two valid item indices/names.")
    } else {
      if (length(item.subset)>1 && all(item.subset %in% rownames(as.matrix(tl)))) {
        tl    <- tl[item.subset]
        keep.subset <- c()
        for (i in rownames(as.matrix(tl)))
          if (tl[i]==1)
            keep.subset<-c(keep.subset,i)
          else
            keep.subset<-c(keep.subset,paste(i,1:tl[i],sep=":"))

        iloc  <- iloc[keep.subset]
        ise   <- ise[keep.subset]
        ifitZ <- ifitZ[item.subset]
        tt<-tt[keep.subset,]
      }
      else if(length(item.subset)!=1 || !(item.subset=="all"))
        stop("item.subset misspecified. Use 'all' or vector of at least two valid item indices/names.")
    }
  } else if (imap) {
    ## Case 2: item subsetting by item numbers
    if (dRm) {
      if (length(item.subset)>1 && all(item.subset %in% 1:nrow(tt))) {
        iloc  <- iloc[item.subset]
        ise   <- ise[item.subset]
        ifitZ <- ifitZ[item.subset]
        tt    <- tt[item.subset,]
      }
      else
        stop("item.subset misspecified. Use 'all' or vector of at least two valid item indices/names.")
    }
    else {
      if (length(item.subset)>1 && all(item.subset %in% 1:length(tl))) {
        tl    <- tl[item.subset]
        keep.subset <- c()
        for (i in rownames(as.matrix(tl)))
          if (tl[i]==1)
            keep.subset<-c(keep.subset,i)
          else
            keep.subset<-c(keep.subset,paste(i,1:tl[i],sep=":"))

        iloc  <- iloc[keep.subset]
        ise   <- ise[keep.subset]
        ifitZ <- ifitZ[item.subset]
        tt<-tt[keep.subset,]
      }
      else
        stop("item.subset misspecified. Use 'all' or vector of at least two valid item indices/names.")
    }
  }

  ## We can now do person subsetting; this is significantly easier than
  ## item subsetting, as there is no dRM/eRm distinction.
  if (pmap && is.character(person.subset)) {
    ## Case 1: person subsetting by person names
    if (length(person.subset)>1 && all(person.subset %in% names(ploc))) {
      ploc  <- ploc[person.subset]
      pse   <- pse[person.subset]
      pfitZ <- pfitZ[person.subset]
    }
    else if(length(person.subset)!=1 || !(person.subset=="all"))
      stop("person.subset misspecified. Use 'all' or vector of at least two valid person indices/names.")
  } else if (pmap) {
    ## Case 2: person subsetting by person numbers
    if (length(person.subset)>1 && all(person.subset %in% 1:length(ploc))) {
      ploc  <- ploc[person.subset]
      pse   <- pse[person.subset]
      pfitZ <- pfitZ[person.subset]
    }
    else
      stop("person.subset misspecified. Use 'all' or vector of at least two valid person indices/names.")
  }

  ## Confidence intervals for persons and items
  ##
  ## Need defaults for multiple of standard error for purpose of range
  ## calculation; these are zero as default is not to draw confidence
  ## intervals
  pci=0
  ici=0

  ## Our calculation is simplistic; we use the normal distribution to
  ## estimate our confidence interval from our standard error.  However,
  ## since this is likely to only be approximate and indicative anyway, we
  ## are not concerned by this.
  if(pmap && !is.null(personCI)) {
    if(is.null(personCI$clevel)) personCI$clevel <- 0.95
    if(is.null(personCI$col))    personCI$col    <- "orange"
    if(is.null(personCI$lty))    personCI$lty    <- "dotted"
    pci <- qnorm((1-personCI$clevel)/2, lower.tail=FALSE)
  }
  if(imap && !is.null(itemCI)) {
    if(is.null(itemCI$clevel)) itemCI$clevel <- 0.95
    if(is.null(itemCI$col))    itemCI$col    <- "red"
    if(is.null(itemCI$lty))    itemCI$lty    <- "dotted"
    ici <- qnorm((1-itemCI$clevel)/2, lower.tail=FALSE)
  }

  ## Now we can plot the Pathway Map

  if (pmap) { ## person map
    xrange.pmap <- range(pfitZ,finite=TRUE)
    xrange.pmap[1] <- min(-2.5,xrange.pmap[1])
    xrange.pmap[2] <- max(2.5,xrange.pmap[2]+1) ## need space for labels
    yrange.pmap<-range(ploc,finite=TRUE)
    yrange.pmap[1]<-yrange.pmap[1]-pci*max(pse)
    yrange.pmap[2]<-yrange.pmap[2]+pci*max(pse)
  }
  if (imap) { ## item map
    xrange.imap <- range(ifitZ,finite=TRUE)
    xrange.imap[1] <- min(-2.5,xrange.imap[1])
    xrange.imap[2] <- max(2.5,xrange.imap[2]+1) ## need space for labels
    yrange.imap<-range(iloc,finite=TRUE)
    yrange.imap[1]<-yrange.imap[1]-ici*max(ise)
    yrange.imap[2]<-yrange.imap[2]+ici*max(ise)
  }

  if (pmap && !imap) {
    xrange <- xrange.pmap
    yrange <- yrange.pmap
    maintitle <- mainperson
  } else if (!pmap && imap) {
    xrange <- xrange.imap
    yrange <- yrange.imap
    maintitle <- mainitem
  } else {
    xrange <- numeric(2)
    yrange <- numeric(2)
    xrange[1] <- min(xrange.pmap[1], xrange.imap[1])
    xrange[2] <- max(xrange.pmap[2], xrange.imap[2])
    yrange[1] <- min(yrange.pmap[1], yrange.imap[1])
    yrange[2] <- max(yrange.pmap[2], yrange.imap[2])
    maintitle <- mainboth
  }


  par(mar=c(5,4,4,2))

  if (!horiz){  # rh 2010-12-09
    plot(xrange,yrange, xlim=xrange, ylim=yrange, main=maintitle,
         ylab=latdim, xlab=tlab, type="n")
    abline(v=c(-2,2),col="lightgreen")
  } else {
    plot(yrange,xrange, xlim=yrange, ylim=xrange, main=maintitle,
         ylab=tlab, xlab=latdim, type="n")
    abline(h=c(-2,2),col="lightgreen")
  }

  if (pmap) { ## person map
    zt <- pfitZ
    if (!horiz){
      if (pci>0) ## draw confidence intervals
        arrows(zt,ploc+pci*pse, zt,ploc-pci*pse, angle=90, code=3, length=0.04,
               col=personCI$col, lty=personCI$lty)
      points(zt,ploc,pch=person.pch,cex=cex.pch)
      text(zt,ploc,names(ploc),cex=cex.gen,pos=4)
    } else {
      if (pci>0) ## draw confidence intervals
        arrows(ploc+pci*pse, zt,ploc-pci*pse, zt, angle=90, code=3, length=0.04,
               col=personCI$col, lty=personCI$lty)
      points(ploc, zt, pch=person.pch,cex=cex.pch)
      text(ploc, zt, names(ploc),cex=cex.gen,pos=4)
    }
  }



  if (imap) { ## item map
    if (dRm)
      zt <- ifitZ
    else
      zt <- rep(ifitZ,times=tl)

    if (!horiz){
      if (ici>0) ## draw confidence intervals
        arrows(zt,iloc+ici*ise, zt,iloc-ici*ise, angle=90, code=3, length=0.04,
               col=itemCI$col, lty=itemCI$lty)
      points(zt,iloc,pch=item.pch,cex=cex.pch)
      text(zt,iloc,rownames(tt),cex=cex.gen,pos=4)
    } else {
      if (ici>0) ## draw confidence intervals
        arrows(iloc+ici*ise, zt,iloc-ici*ise,zt, angle=90, code=3, length=0.04,
               col=itemCI$col, lty=itemCI$lty)
      points(iloc, zt,pch=item.pch,cex=cex.pch)
      text(iloc,zt, rownames(tt),cex=cex.gen,pos=4)
    }
  }

  par(def.par)

  invisible(NULL)
}