R version 3.2.2 (2015-08-14) -- "Fire Safety" Copyright (C) 2015 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > > # 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) race trauma happy "factor" "integer" "factor" > levels(foo$race) [1] "Black" "White" > levels(foo$happy) [1] "Not too happy" "Pretty happy" "Very happy" > > foo <- transform(foo, happy = ordered(happy)) > sapply(foo, class) $race [1] "factor" $trauma [1] "integer" $happy [1] "ordered" "factor" > levels(foo$happy) [1] "Not too happy" "Pretty happy" "Very happy" > > library(MASS) > pout <- polr(happy ~ race + trauma, data = foo) > summary(pout) Re-fitting to get Hessian Call: polr(formula = happy ~ race + trauma, data = foo) Coefficients: Value Std. Error t value raceWhite 2.0362 0.6859 2.968 trauma -0.4056 0.1830 -2.216 Intercepts: Value Std. Error t value Not too happy|Pretty happy -1.3644 0.6540 -2.0863 Pretty happy|Very happy 2.5542 0.7201 3.5473 Residual Deviance: 148.407 AIC: 156.407 > > # 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. > > proc.time() user system elapsed 0.263 0.037 1.036