polr.plot <- function( model , x , category.labels=TRUE , main=names(model$model)[1] , ylab="Probability" , xlab=x , lty=2 , ... ) {
    xvalues <- model$model[,x]
    yvalues <- model$model[,1]
    b <- coef(model)[x]
    z <- model$zeta
    nz <- length(z)
    
    # add up terms for all x variables, using mean value of each
    xmeans <- rep(0,ncol(model$model)-1)
    for ( i in 2:ncol(model$model) ) {
        xmeans[i-1] <- mean( as.numeric(model$model[,i]) )
    }
    ncoef <- length(coef(model))
    if ( ncoef > length(xmeans) ) {
        ncoef <- length(xmeans)-1
    }
    modelsum <- sum( -coef(model)[1:ncoef] * xmeans[1:ncoef] )
    # subtract x variable we will use as x-axis
    modelsum <- modelsum - (-coef(model)[x])*mean(model$model[,x])
    
    xseq <- seq( from=min(xvalues) , to=max(xvalues) , by=(max(xvalues)-min(xvalues))/10 )
    
    p1 <- exp( z[1] - b * xseq + modelsum ) / ( 1 + exp( z[1] - b * xseq + modelsum ) )
    
    plot( xseq , p1 , type="l" , ylim=c(0,1) , lty=lty , main=main , xlab=xlab , ylab=ylab , ... )
    
    if ( category.labels==TRUE ) {
        # find widest vertical distance, and place category label there
        dif <- p1 - 0
        l.x <- xseq[dif==max(dif)]
        # y coord is half way between the two lines
        l.y <- p1[dif==max(dif)] - max(dif)/2
        # draw label
        text( l.x , l.y , levels(model$model[,1])[1] )
    }
    
    pprev <- p1
    if ( nz > 1 ) {
        for ( i in 2:nz ) {
            p <- exp( z[i] - b * xseq + modelsum ) / ( 1 + exp( z[i] - b * xseq + modelsum ) )
            lines( xseq , p , lty=lty )
            if ( category.labels==TRUE ) {
                # find widest vertical distance, and place category label there
                dif <- p - pprev
                l.x <- xseq[dif==max(dif)]
                # y coord is half way between the two lines
                l.y <- pprev[dif==max(dif)] + max(dif)/2
                # draw label
                text( l.x , l.y , levels(model$model[,1])[i] )
            }
            pprev <- p
        }
    }
    
    if ( category.labels==TRUE ) {
        # find widest vertical distance, and place category label there
        dif <- 1 - pprev
        l.x <- xseq[dif==max(dif)]
        # y coord is half way between the two lines
        l.y <- pprev[dif==max(dif)] + max(dif)/2
        # draw label
        text( l.x , l.y , levels(model$model[,1])[nz+1] )
    }
    lines( c(min(xvalues),max(xvalues)), c(0,0) , lty=lty )
    lines( c(min(xvalues),max(xvalues)), c(1,1) , lty=lty )
}

### example
#x <- rnorm(50,mean=0)
#x2 <- rnorm(50,mean=0)
#x3 <- rnorm(50,mean=5)
#x4 <- rnorm(50,mean=-5)
#y <- rbinom(50,prob=1/(1+exp(0 -1/2*x+1/4*x2)),size=5)
#y <- ifelse( y>0 , y , 1 )

#library(MASS)
#m <- polr( as.ordered(y) ~ x + x2 + x3 + x4 )

#par(mfcol=c(1,4))
#polr.plot( m , x="x" , category.labels=TRUE )
#polr.plot( m , x="x2" )
#polr.plot( m , x="x3" )
#polr.plot( m , x="x4" )