source("http://192.38.117.59/~linearpredictors/datafiles/readFibrosis.R")

## proportional odds logistic regression model
library(MASS)
prop.odds <- polr(stage ~ log2(ha) + log2(p3np) + log2(ykl40),
                   data = fibrosis)

used.data <- na.omit(fibrosis)

plot(0,0, type = "n",
     xlim = c(-0.25, 3.25),
     xaxp = c(0, 3, 3),
     ylim = c(0, 0.5),
     xlab = "Observed Fibrosis degree",
     ylab = "Average probabilities")

## fitted probabilities for each stage split according to observed stage
fitted.p <<- split(data.frame(prop.odds$fitted), used.data$stage)

## means of the fitted probabilities. columns = observed stage
mean.fitted.p <<- sapply(fitted.p, function(d) apply(d, 2, mean))

for (s in 0:3){
  points(0:3, mean.fitted.p[s+1,], type = "b")
  points(s, mean.fitted.p[s+1,s+1], pch = 19)
  
  ## a curve signifies the average probability for hitting a level
  text(3.2, mean.fitted.p[s+1, 4], bquote(p[.(s)]))
}