ordinal (version 2011.05-10)

clm: Cumulative link models

Description

Fits cumulative link models with an additive model for the location and a multiplicative model for the scale. The function allows for structured thresholds. A popular special case of a CLM is the proportional odds model. In addition to the standard link functions, two flexible link functions, "Arandar-Ordaz" and "log-gamma" are available, where an extra link function parameter provides additional flexibility. A subset of the predictors can be allowed to have nominal rather than ordinal effects. This has been termed "partial proportional odds" when the link is the logistic.

Usage

clm(location, scale, nominal, data, weights, start, subset,
    na.action, contrasts, Hess = TRUE, model,
    link = c("logistic", "probit", "cloglog", "loglog",
    "cauchit", "Aranda-Ordaz", "log-gamma"), lambda,
    doFit = TRUE, control,
    threshold = c("flexible", "symmetric", "equidistant"), ...)

Arguments

location
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.
scale
a optional formula expression as for the location part, of the form ~ predictors, i.e. with an empty left hand side. If left unspecified, the model assumes a constant scale for models with factor response and estimates the residual s
nominal
an optional formula of the form ~ predictors, i.e. with an empty left hand side. The effects of the predictors in this formula are assumed to nominal.
data
an optional data frame in which to interpret the variables occurring in the formulas.
weights
optional case weights in fitting. Defaults to 1.
start
initial values for the parameters in the format c(alpha, beta, log(zeta), lambda).
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. Applies to terms in all three formulae.
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 inverse of the observed information matrix) should be computed. Use Hess = TRUE if you intend to call summary or vcov on the fit and Hess = FALSE in all o
model
logical for whether the model frames should be part of the returned object.
link
link function, i.e. the type of location-scale distribution assumed for the latent distribution. The Aranda-Ordaz and log-gamma links add additional flexibility with a link function parameter, lambda. The
lambda
numerical scalar: the link function parameter. Used in combination with link Aranda-Ordaz or log-gamma and otherwise ignored. If lambda is specified, the model is estimated with lambda fixed at this value and otherwise l
doFit
logical for whether the model should be fit or the model environment should be returned.
control
a call to clm.control.
threshold
specifies a potential structure for the thresholds (cut-points). "flexible" provides the standard unstructured thresholds, "symmetric" restricts the distance between the thresholds to be symmetric around the central one
...
additional arguments are passed on to clm.control and possibly further on to the optimizer, which can lead to surprising error or warning messages when mistyping arguments etc.

Value

  • If doFit = FALSE the result is an environment representing the model ready to be optimized. If doFit = TRUE the result is an object of class "clm" with the following components:
  • betathe parameter estimates of the location part.
  • zetathe parameter estimates of the scale part on the log scale; the scale parameter estimates on the original scale are given by exp(zeta).
  • Alphavector or matrix of the threshold parameters.
  • Thetavector or matrix of the thresholds.
  • xivector of threshold parameters, which, given a threshold function (e.g. "equidistant"), and possible nominal effects define the class boundaries, Theta.
  • lambdathe value of lambda if lambda is supplied or estimated, otherwise missing.
  • coefficientsthe coefficients of the intercepts (theta), the location (beta), the scale (zeta), and the link function parameter (lambda).
  • df.residualthe number of residual degrees of freedoms, calculated using the weights.
  • fitted.valuesvector of fitted values for each observation. An observation here is each of the scalar elements of the multinomial table and not a multinomial vector.
  • convergenceTRUE if the gradient based convergence criterion is met and FALSE otherwise.
  • gradientvector of gradients for all the parameters at termination of the optimizer.
  • optReslist with results from the optimizer. The contents of the list depends on the choice of optimizer.
  • logLikthe log likelihood of the model at optimizer termination.
  • Hessianif the model was fitted with Hess = TRUE, this is the Hessian matrix of the parameters at the optimum.
  • scalemodel.frame for the scale model.
  • locationmodel.frame for the location model.
  • nominalmodel.frame for the nominal model.
  • edfthe (effective) number of degrees of freedom used by the model.
  • startthe starting values.
  • convTolconvergence tolerance for the maximum absolute gradient of the parameters at termination of the optimizer.
  • methodcharacter, the optimizer.
  • ythe response variable.
  • levthe names of the levels of the response variable.
  • nobsthe (effective) number of observations, calculated as the sum of the weights.
  • thresholdcharacter, the threshold function used in the model.
  • estimLambda1 if lambda is estimated in one of the flexible link functions and 0 otherwise.
  • linkcharacter, the link function used in the model.
  • callthe matched call.
  • contrastscontrasts applied to terms in location and scale models.
  • na.actionthe function used to filter missing data.

Details

There are methods for the standard model-fitting functions, including summary, vcov, predict, anova, logLik, profile, plot.profile, confint, update, dropterm, addterm, and an extractAIC method. The design of the implementation is inspired by an idea proposed by Douglas Bates in the talk "Exploiting sparsity in model matrices" presented at the DSC conference in Copenhagen, July 14 2009. Basically an environment is set up with all the information needed to optimize the likelihood function. Extractor functions are then used to get the value of likelihood at current or given parameter values and to extract current values of the parameters. All computations are performed inside the environment and relevant variables are updated during the fitting process. After optimizer termination relevant variables are extracted from the environment and the remaining are discarded. Some aspects of clm, for instance, how starting values are obtained, and of the associated methods are inspired by polr from package MASS.

References

Agresti, A. (2002) Categorical Data Analysis. Second edition. Wiley. Aranda-Ordaz, F. J. (1983) An Extension of the Proportional-Hazards Model for Grouped Data. Biometrics, 39, 109-117. Genter, F. C. and Farewell, V. T. (1985) Goodness-of-link testing in ordinal regression models. The Canadian Journal of Statistics, 13(1), 37-44. Christensen, R. H. B., Brockhoff, P. B. and Cleaver, G. (2010) Statistical and Thurstonian models for the A-not A protocol with and without sureness. Manuscript for Food Quality and Preference.

Examples

Run this code
options(contrasts = c("contr.treatment", "contr.poly"))
data(soup)

## More manageable data set:
(tab26 <- with(soup, table("Product" = PROD, "Response" = SURENESS)))
dimnames(tab26)[[2]] <- c("Sure", "Not Sure", "Guess", "Guess", "Not Sure", "Sure")
dat26 <- expand.grid(sureness = as.factor(1:6), prod = c("Ref", "Test"))
dat26$wghts <- c(t(tab26))

m1 <- clm(sureness ~ prod, scale = ~prod, data = dat26,
          weights = wghts, link = "logistic")

## print, summary, vcov, logLik, AIC:
m1
summary(m1)
vcov(m1)
logLik(m1)
AIC(m1)
coef(m1)
coef(summary(m1))

## link functions:
m2 <- update(m1, link = "probit")
m3 <- update(m1, link = "cloglog")
m4 <- update(m1, link = "loglog")
m5 <- update(m1, link = "cauchit", start = coef(m1))
m6 <- update(m1, link = "Aranda-Ordaz", lambda = 1)
m7 <- update(m1, link = "Aranda-Ordaz")
m8 <- update(m1, link = "log-gamma", lambda = 1)
m9 <- update(m1, link = "log-gamma")

## nominal effects:
mN1 <-  clm(sureness ~ 1, nominal = ~ prod, data = dat26,
            weights = wghts, link = "logistic")
anova(m1, mN1)

## optimizer / method:
update(m1, scale = ~ 1, method = "Newton")
update(m1, scale = ~ 1, method = "nlminb")
update(m1, scale = ~ 1, method = "optim")
update(m1, scale = ~ 1, method = "model.frame")
update(m1, location = ~.-prod, scale = ~ 1,
       nominal = ~ prod, method = "model.frame")

## threshold functions
mT1 <- update(m1, threshold = "symmetric")
mT2 <- update(m1, threshold = "equidistant")
anova(m1, mT1, mT2)

## Extend example from polr in package MASS:
## Fit model from polr example:
data(housing, package = "MASS")
fm1 <- clm(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
fm1
summary(fm1)
## With probit link:
summary(update(fm1, link = "probit"))
## Allow scale to depend on Cont-variable
summary(fm2 <- update(fm1, scale =~ Cont))
anova(fm1, fm2)
## which seems to improve the fit

#################################
## It is possible to fit multinomial models (i.e. with nominal
## effects) as the following example shows:
library(nnet)
(hous1.mu <- multinom(Sat ~ 1, weights = Freq, data = housing))
(hous1.clm <- clm(Sat ~ 1, weights = Freq, data = housing))

## It is the same likelihood:
all.equal(logLik(hous1.mu), logLik(hous1.clm))

## and the same fitted values:
fitHous.mu <-
    t(fitted(hous1.mu))[t(col(fitted(hous1.mu)) == unclass(housing$Sat))]
all.equal(fitted(hous1.clm), fitHous.mu)

## The coefficients of multinom can be retrieved from the clm-object
## by:
Pi <- diff(c(0, plogis(hous1.clm$xi), 1))
log(Pi[2:3]/Pi[1])

## A larger model with explanatory variables:
(hous.mu <- multinom(Sat ~ Infl + Type + Cont, weights = Freq, data = housing))
(hous.clm <- clm(Sat ~ 1, nominal = ~ Infl + Type + Cont, weights = Freq,
                 data = housing))

## Almost the same likelihood:
all.equal(logLik(hous.mu), logLik(hous.clm))

## And almost the same fitted values:
fitHous.mu <-
    t(fitted(hous.mu))[t(col(fitted(hous.mu)) == unclass(housing$Sat))]
all.equal(fitted(hous.clm), fitHous.mu)
all.equal(round(fitted(hous.clm), 5), round(fitHous.mu), 5)

Run the code above in your browser using DataCamp Workspace