# Ordered Logistic or Probit Regression Fits a logistic or probit regression model to an ordered factor response. The default logistic case is **proportional odds logistic regression**, after which the function is named. ```r polr(formula, data, weights, start, ..., subset, na.action, contrasts = NULL, Hess = FALSE, model = TRUE, method = c("logistic", "probit", "loglog", "cloglog", "cauchit")) ``` ## Arguments - `formula`: a formula expression as for regression models, of the form `response ~ predictors`. The response should be a factor (preferably an ordered factor), which will be interpreted as an ordinal response, with levels ordered as in the factor. The model must have an intercept: attempts to remove one will lead to a warning and be ignored. An offset may be used. See the documentation of `formula` for other details. - `data`: an optional data frame, list or environment in which to interpret the variables occurring in `formula`. - `weights`: optional case weights in fitting. Default to 1. - `start`: initial values for the parameters. This is in the format `c(coefficients, zeta)`: see the Values section. - ``...``: additional arguments to be passed to `optim`, most often a `control` argument. - `subset`: expression saying which subset of the rows of the data should be used in the fit. All observations are included by default. - `na.action`: a function to filter missing data. - `contrasts`: a list of contrasts to be used for some or all of the factors appearing as variables in the model formula. - `Hess`: logical for whether the Hessian (the observed information matrix) should be returned. Use this if you intend to call `summary` or `vcov` on the fit. - `model`: logical for whether the model matrix should be returned. - `method`: logistic or probit or (complementary) log-log or cauchit (corresponding to a Cauchy latent variable). ## Details 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 = -\infty < \zeta_1 < \cdots < \zeta_K = \inftyzeta_0 = -Inf < zeta_1 < \dots < zeta_K = Inf $$ This leads to the model $$ \mbox{logit} P(Y \le k | x) = \zeta_k - \etalogit 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`). In the logistic case, the left-hand side of the last display is the log odds of category $k$ or less, and since these are log odds which differ only by a constant for different $k$, the odds are proportional. Hence the term **proportional odds logistic regression**. The log-log and complementary log-log links are the increasing functions $F^-1(p) = -log(-log(p))$ and $F^-1(p) = log(-log(1-p))$; some call the first the negative log-log link. These correspond to a latent variable with the extreme-value distribution for the maximum and minimum respectively. A **proportional hazards** model for grouped survival times can be obtained by using the complementary log-log link with grouping ordered by increasing times. `predict`, `summary`, `vcov`, `anova`, `model.frame` and an `extractAIC` method for use with `stepAIC` (and `step`). There are also `profile` and `confint` methods. ## Returns A object of class `"polr"`. This has components - **coefficients**: the coefficients of the linear predictor, which has no intercept. - **zeta**: the intercepts for the class boundaries. - **deviance**: the residual deviance. - **fitted.values**: a matrix, with a column for each level of the response. - **lev**: the names of the response levels. - **terms**: the `terms` structure describing the model. - **df.residual**: the number of residual degrees of freedoms, calculated using the weights. - **edf**: the (effective) number of degrees of freedom used by the model - **n, nobs**: the (effective) number of observations, calculated using the weights. (`nobs` is for use by `stepAIC`. - **call**: the matched call. - **method**: the matched method used. - **convergence**: the convergence code returned by `optim`. - **niter**: the number of function and gradient evaluations used by `optim`. - **lp**: the linear predictor (including any offset). - **Hessian**: (if `Hess` is true). Note that this is a numerical approximation derived from the optimization proces. - **model**: (if `model` is true). ## Note The `vcov` method uses the approximate Hessian: for reliable results the model matrix should be sensibly scaled with all columns having range the order of one. Prior to version 7.3-32, `method = "cloglog"` confusingly gave the log-log link, implicitly assuming the first response level was the best . ## References Agresti, A. (2002) **Categorical Data.** Second edition. Wiley. Venables, W. N. and Ripley, B. D. (2002) **Modern Applied Statistics with S.** Fourth edition. Springer. ## See Also `optim`, `glm`, `multinom`. ## Examples ```r options(contrasts = c("contr.treatment", "contr.poly")) house.plr <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing) house.plr summary(house.plr, digits = 3) ## slightly worse fit from summary(update(house.plr, method = "probit", Hess = TRUE), digits = 3) ## although it is not really appropriate, can fit summary(update(house.plr, method = "loglog", Hess = TRUE), digits = 3) summary(update(house.plr, method = "cloglog", Hess = TRUE), digits = 3) predict(house.plr, housing, type = "p") addterm(house.plr, ~.^2, test = "Chisq") house.plr2 <- stepAIC(house.plr, ~.^2) house.plr2$anova anova(house.plr, house.plr2) house.plr <- update(house.plr, Hess=TRUE) pr <- profile(house.plr) confint(pr) plot(pr) pairs(pr) ```

polr function