1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
|
% file MASS/man/polr.Rd
% copyright (C) 1994-2014 W. N. Venables and B. D. Ripley
%
\name{polr}
\alias{polr}
\title{
Ordered Logistic or Probit Regression
}
\description{
Fits a logistic or probit regression model to an ordered factor
response. The default logistic case is \emph{proportional odds
logistic regression}, after which the function is named.
}
\usage{
polr(formula, data, weights, start, \dots, subset, na.action,
contrasts = NULL, Hess = FALSE, model = TRUE,
method = c("logistic", "probit", "loglog", "cloglog", "cauchit"))
}
\arguments{
\item{formula}{
a formula expression as for regression models, of the form
\code{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 \code{\link{formula}} for other details.
}
\item{data}{
an optional data frame in which to interpret the variables occurring
in \code{formula}.
}
\item{weights}{
optional case weights in fitting. Default to 1.
}
\item{start}{
initial values for the parameters. This is in the format
\code{c(coefficients, zeta)}: see the Values section.
}
\item{\dots}{
additional arguments to be passed to \code{\link{optim}}, most often a
\code{control} argument.
}
\item{subset}{
expression saying which subset of the rows of the data should be used
in the fit. All observations are included by default.
}
\item{na.action}{
a function to filter missing data.
}
\item{contrasts}{
a list of contrasts to be used for some or all of
the factors appearing as variables in the model formula.
}
\item{Hess}{
logical for whether the Hessian (the observed information matrix)
should be returned. Use this if you intend to call \code{summary} or
\code{vcov} on the fit.
}
\item{model}{
logical for whether the model matrix should be returned.
}
\item{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 \emph{cumulative link}
model. The basic interpretation is as a \emph{coarsened} version of a
latent variable \eqn{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 \eqn{Y_i} falls into with breakpoints
\deqn{\zeta_0 = -\infty < \zeta_1 < \cdots < \zeta_K = \infty}{zeta_0 = -Inf < zeta_1 < \dots < zeta_K = Inf}
This leads to the model
\deqn{\mbox{logit} P(Y \le k | x) = \zeta_k - \eta}{logit P(Y <= k | x) = zeta_k - eta}
with \emph{logit} replaced by \emph{probit} for a normal latent
variable, and \eqn{\eta}{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 \eqn{\eta}{eta} (and hence the coefficients \code{beta}).
In the logistic case, the left-hand side of the last display is the
log odds of category \eqn{k} or less, and since these are log odds
which differ only by a constant for different \eqn{k}, the odds are
proportional. Hence the term \emph{proportional odds logistic
regression}.
The log-log and complementary log-log links are the increasing functions
\eqn{F^{-1}(p) = -log(-log(p))}{F^-1(p) = -log(-log(p))} and
\eqn{F^{-1}(p) = log(-log(1-p))}{F^-1(p) = log(-log(1-p))};
some call the first the \sQuote{negative log-log} link. These
correspond to a latent variable with the extreme-value distribution for
the maximum and minimum respectively.
A \emph{proportional hazards} model for grouped survival times can be
obtained by using the complementary log-log link with grouping ordered
by increasing times.
\code{\link{predict}}, \code{\link{summary}}, \code{\link{vcov}},
\code{\link{anova}}, \code{\link{model.frame}} and an
\code{extractAIC} method for use with \code{\link{stepAIC}} (and
\code{\link{step}}). There are also \code{\link{profile}} and
\code{\link{confint}} methods.
}
\value{
A object of class \code{"polr"}. This has components
\item{coefficients}{the coefficients of the linear predictor, which has no
intercept.}
\item{zeta}{the intercepts for the class boundaries.}
\item{deviance}{the residual deviance.}
\item{fitted.values}{a matrix, with a column for each level of the response.}
\item{lev}{the names of the response levels.}
\item{terms}{the \code{terms} structure describing the model.}
\item{df.residual}{the number of residual degrees of freedoms,
calculated using the weights.}
\item{edf}{the (effective) number of degrees of freedom used by the model}
\item{n, nobs}{the (effective) number of observations, calculated using the
weights. (\code{nobs} is for use by \code{\link{stepAIC}}.}
\item{call}{the matched call.}
\item{method}{the matched method used.}
\item{convergence}{the convergence code returned by \code{optim}.}
\item{niter}{the number of function and gradient evaluations used by
\code{optim}.}
\item{lp}{the linear predictor (including any offset).}
\item{Hessian}{(if \code{Hess} is true). Note that this is a
numerical approximation derived from the optimization proces.}
\item{model}{(if \code{model} is true).}
}
\note{
The \code{\link{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, \code{method = "cloglog"} confusingly gave
the log-log link, implicitly assuming the first response level was the
\sQuote{best}.
}
\references{
Agresti, A. (2002) \emph{Categorical Data.} Second edition. Wiley.
Venables, W. N. and Ripley, B. D. (2002)
\emph{Modern Applied Statistics with S.} Fourth edition. Springer.
}
\seealso{
\code{\link{optim}}, \code{\link{glm}}, \code{\link[nnet]{multinom}}.
}
\examples{
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)
}
\keyword{models}
|
