I have greatly enjoyed the conversation here, however I feel that the answers did not correctly address all the (very good) components to the question you put forth. The second half of the example page for polr
is all about profiling. A good technical reference here is Venerables and Ripley who discuss profiling and what it does. This is a critical technique when you step beyond the comfort zone of fitting exponential family models with full likelihood (regular GLMs).
The key departure here is the use of categorical thresholds. You will notice POLR doesn't estimate a usual intercept term. Rather, there are $k-1$ nuisance parameters: thresholds for which the fitted risk tends to fall in a certain cumulative of the $k$ possible categories. Since these thresholds are never jointly estimated, their covariance with model parameters is unknown. Unlike GLMs we cannot "perturb" a coefficient by an amount and be certain of how it might affect other estimates. We use profiling to do this accounting for the nuisance thresholds. Profiling is an immense subject, but basically the goal is robustly measuring the covariance of regression coefficients when the model is maximizing an irregular likelihood, like with lmer
, nls
, polr
, and glm.nb
.
The help page for ?profile.glm
should be of some use as polr
objects are essentially GLMs (plus the categorical thresholds). Lastly, you can actually peak at the source code, if it's of any use, using getS3method('profile', 'polr')
. I use this getS3method
function a lot because, while R seems to insist many methods should be hidden, one can learn surprisingly much about implementation and methods by reviewing code.
•What information does pr contain? The help page on profile is
generic, and gives no guidance for polr.
pr
is a profile.polr, profile
object (inherited class profile
). There's an entry for each covariate. The profiler loops over each covariate, recalculates the optimal model fit with that covariate fixed to some slightly different amount. The output shows the covariate's fixed value measured as a scaled "z-score" difference from its estimated value and the resulting fixed effects in other covariates. For instance, if you look at pr$InflMedium
, you'll note that, when "z" is 0, the other fixed effects are the same as are found in the original fit.
•What is plot(pr) showing? I see six graphs. Each has an X axis that
is numeric, although the label is an indicator variable (looks like an
input variable that is an indicator for an ordinal value). Then the Y
axis is "tau" which is completely unexplained.
Again, ?plot.profile
gives the description. The plot roughly shows how the regression coefficients covary. tau is the scaled difference, the z score before, so it's 0 value gives the optimal fit coefficients, depicted with a tick mark. You wouldn't tell for this fit is so well behaved, but those "lines" are actually splines. If the likelihood were very irregularly behaved at the optimal fit, you would observe strange and unpredictable behavior in the plot. This would behoove you to estimate the output using a more robust error estimation (bootstrap/jackknife), to calculate CIs using method='profile'
, to recode variables, or to perform other diagnostics.
•What is pairs(pr) showing? It looks like a plot for each pair of
input variables, but again I see no explanation of the X or Y axes.
The help file says: "The pairs method shows, for each pair of parameters x and y, two curves intersecting at the maximum likelihood estimate, which give the loci of the points at which the tangents to the contours of the bivariate profile likelihood become vertical and horizontal, respectively. In the case of an exactly bivariate normal profile likelihood, these two curves would be straight lines giving the conditional means of y|x and x|y, and the contours would be exactly elliptical." Basically, they again help you to visualize the confidence ellipses. Non-orthogonal axes indicate highly covariable measures, such as InfMedium and InfHigh which are intuitively very related. Again, irregular likelihoods would lead to images that are quite baffling here.
•How can one understand if the model gave a good fit?
summary(house.plr) shows Residual Deviance 3479.149 and AIC (Akaike
Information Criterion?) of 3495.149. Is that good? In the case those
are only useful as relative measures (i.e. to compare to another model
fit), what is a good absolute measure? Is the residual deviance
approximately chi-squared distributed? Can one use "% correctly
predicted" on the original data or some cross-validation? What is the
easiest way to do that?
One assumption that is good to assess is the proportional odds assumption. This is reflected somewhat in the global test (which assesses polr against a saturated loglinear model). A limitation here is that with large data, global tests always fail. As a result, using graphics and inspecting estimates (betas) and precision (SEs) for the loglinear model and polr fit is a good idea. If they massively disagree, something is perhaps wrong.
With ordered outcomes, it is hard to define percent agreement. How will you choose a classifier based on the model, and if you do how will you suss poor performance from a poor classifier. mode
is a bad choice. If I have 10 category logits and my prediction is always but one category off, perhaps that is not a bad thing. Further, my model may correctly predict a 40% chance of a 0 response, but also 20% chances of 8, 9, 10. So if I observe 9 is that good or bad? If you must measure agreement use a weighted kappa, or even MSE. The loglinear model will always produce the best agreement. That is not what the POLR does.
•How does one apply and interpret anova on this model? The docs say
"There are methods for the standard model-fitting functions, including
predict, summary, vcov, anova." However, running anova(house.plr)
results in anova is not implemented for a single "polr" object
You can test nested models with waldtest
and lrtest
in the lmtest
package in R. This is equivalent to ANOVA. The interpretation is exactly the same as with GLMs.
•How does one interpret the t values for each coefficient? Unlike some
model fits, there are no P values here.
Again, unlike linear models, the POLR model is capable of having issues with irregular likelihood so inference based on the Hessian can be very unstable. It is analogous to fitting mixed models, see for instance the helpfile on confint.merMod
for the lme4 package. Here, the assessments made with profiling show the covariance is well behaved. The programmers would have done this by default, except that profiling can be computationally very intensive, and thus they leave it to your hands. If you must see the Wald based inference, use coeftest(house.plr)
from the lrtest
package.
methods("profile")
will give you the (S3 in this case) methods associated to an Rprofile
object, then you will see than there is a dedicated method forpolr
results, that you can browse on-line by typinggetAnywhere("profile.polr")
at the R prompt. $\endgroup$ – chl Feb 28 '11 at 21:13