# we do example 9.2.4 in Agresti

foo <- read.table("http://www.stat.umn.edu/geyer/5421/data/table-8.5.txt",
    header = TRUE)
sapply(foo, class)
levels(foo$race)
levels(foo$happy)

foo <- transform(foo, happy = ordered(happy))
sapply(foo, class)
levels(foo$happy)

library(MASS)
pout <- polr(happy ~ race + trauma, data = foo)
summary(pout)

# see that the coefficient for trauma agrees with that in the SAS output
# shown in Agresti Table 8.6 and the coefficient for race agrees except
# for sign (which can be changed by convention).  So we must be fitting
# the same model.  But the "intercepts" seem to have nothing to do with
# each other, so the definitions of the models must be slightly different.
#
# The polr web page says
#
#     This model is what Agresti (2002) calls a _cumulative link_ model.
#     The basic interpretation is as a _coarsened_ version of a latent
#     variable Y_i which has a logistic or normal or extreme-value or
#     Cauchy distribution with scale parameter one and a linear model
#     for the mean.  The ordered factor which is observed is which bin
#     Y_i falls into with breakpoints
#
#                 zeta_0 = -Inf < zeta_1 < ... < zeta_K = Inf            
#     
#     This leads to the model
#
#                     logit P(Y <= k | x) = zeta_k - eta                 
#     
#     with _logit_ replaced by _probit_ for a normal latent variable,
#     and eta being the linear predictor, a linear function of the
#     explanatory variables (with no intercept).  Note that it is quite
#     common for other software to use the opposite sign for eta (and
#     hence the coefficients ‘beta’).
#
# and the second displayed equation quoted here certainly looks like
# the first displayed equation on the page in Agresti where Table 8.6 appears.