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

library(MASS)

fibrosis$stage3vs012 <- ifelse(fibrosis$stage == 3, 1, 0)
fibrosis$stage23vs01 <- ifelse(fibrosis$stage %in% 2:3, 1, 0)
fibrosis$stage123vs0 <- ifelse(fibrosis$stage %in% 1:3, 1, 0)

## throwing away 3 missing observations
fibrosis.compl <- na.omit(fibrosis)

propodds <- polr(stage ~ log2(ha) + log2(p3np) + log2(ykl40),
                 data = fibrosis.compl,
                 Hess = T ## not to recalculate Hessian later
                 )

round(exp(cbind(Estimate = propodds$coef, confint.default(propodds))),
      digits = 2)


## SD of standardized covariates.
## pi/sqrt(3) is the SD of the Bernoulli distribution

SD <- apply(subset(fibrosis.compl, select = c(ha, p3np, ykl40)), 2,
            function(x)sd(log2(x)/(pi/sqrt(3))))

exp(propodds$coef*SD)

## LR tests for effects of single covariates
dropterm(propodds, test = "Chisq")

## Wald tests can be carried out by comparing the squared t-values
## to a chisq distribution with 1 df.

## The centered model mentioned in the book is
propoddsCent <- polr(stage ~ log2(ha/100) + log2(p3np/10) + log2(ykl40/500),
                     data = fibrosis.compl)