ordinal (version 2010.03-04)

clmm: Cumulative link mixed models

Description

Fits cumulative link mixed models, i.e. cumulative link models with random effects via the Laplace approximation or the adaptive Gauss-Hermite quadrature approximation. The functionality in clm is also implemented here. Currently only a single random term is allowed in the location-part of the model.

Usage

clmm(location, scale, nominal, random, data, weights, start, subset,
     na.action, contrasts = NULL, Hess = FALSE, model = TRUE,
     method = c("ucminf", "nlminb"),
     link = c("logistic", "probit", "cloglog", "loglog",
     "cauchit", "Aranda-Ordaz", "log-gamma"), lambda = NULL,
     doFit = TRUE,  control = list(), nAGQ = 1,
     threshold = c("flexible", "symmetric", "equidistant"), ...)

Arguments

location

as in clm.

scale

as in clm.

nominal

as in clm.

random

a factor for the random effects in the location-part of the model.

data

as in clm.

weights

as in clm.

start

initial values for the parameters in the format c(alpha, beta, log(zeta), lambda, log(stDev)) where stDev is the standard deviation of the random effects.

subset

as in clm.

na.action

as in clm.

contrasts

as in clm.

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 othe

model

as in clm.

method

the optimizer used to maximize the likelihood function. Note that "Newton" and "optim" are not available for clmms. model.frame simply returns a list of model frames with the location, scale and

link

as in clm.

lambda

as in clm.

doFit

as in clm.

control

as in clm.

threshold

as in clm.

nAGQ

the number of quadrature points to be used in the adaptive Gauss-Hermite quadrature approximation to the marginal likelihood. Defaults to 1 which leads to the Laplace approximation. An odd number of quadrature points is encoura

...

additional arguments---currently not used.

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 "clmm" with the following components:
• stDev
  the standard deviation of the random effects.
• Niter
  the total number of iterations in the Newton updates of the conditional modes of the random effects.
• grFac
  the grouping factor defining the random effects.
• nAGQ
  the number of quadrature points used in the adaptive Gauss-Hermite Quadrature approximation to the marginal likelihood.
• condVar
  the conditional variances of the random effects at their conditional modes.
• beta
  the parameter estimates of the location part.
• zeta
  the parameter estimates of the scale part on the log scale; the scale parameter estimates on the original scale are given by exp(zeta).
• Alpha
  vector or matrix of the threshold parameters.
• Theta
  vector or matrix of the thresholds.
• xi
  vector of threshold parameters, which, given a threshold function (e.g. "equidistant"), and possible nominal effects define the class boundaries, Theta.
• lambda
  the value of lambda if lambda is supplied or estimated, otherwise missing.
• coefficients
  the coefficients of the intercepts (theta), the location (beta), the scale (zeta), and the link function parameter (lambda).
• df.residual
  the number of residual degrees of freedoms, calculated using the weights.
• fitted.values
  vector of fitted values for each observation. An observation here is each of the scalar elements of the multinomial table and not a multinomial vector.
• convergence
  TRUE if the optimizer terminates wihtout error and FALSE otherwise.
• gradient
  vector of gradients for the unit-variance random effects at their conditional modes.
• optRes
  list with results from the optimizer. The contents of the list depends on the choice of optimizer.
• logLik
  the log likelihood of the model at optimizer termination.
• Hessian
  if the model was fitted with Hess = TRUE, this is the Hessian matrix of the parameters at the optimum.
• scale
  model.frame for the scale model.
• location
  model.frame for the location model.
• nominal
  model.frame for the nominal model.
• edf
  the (effective) number of degrees of freedom used by the model.
• start
  the starting values.
• method
  character, the optimizer.
• y
  the response variable.
• lev
  the names of the levels of the response variable.
• nobs
  the (effective) number of observations, calculated as the sum of the weights.
• threshold
  character, the threshold function used in the model.
• estimLambda
  1 if lambda is estimated in one of the flexible link functions and 0 otherwise.
• link
  character, the link function used in the model.
• call
  the matched call.
• contrasts
  contrasts applied to terms in location and scale models.
• na.action
  the function used to filter missing data.

Details

The function calls clm to set op the computing environment and to get starting values. The function is implemented in pure R-code, so it will not be as fast as glmer, which have of its code in C, for binomial fits. A Newton scheme is used to obtain the conditional modes of the random effects and a non-linear optimization is performed over the fixed parameter set to get the maximum likelihood estimates. The Newton scheme uses the observed Hessian rather than the expected as is done in e.g. glmer, so results from the Laplace approximation for binomial fits may differ slightly. Control parameters can be supplied to manage the optimization process. The supplied set of control parameters should match those relevant for the chosen optimizer; see ucminf, and nlminb. There are methods for the standard model-fitting functions, including summary, vcov, anova, logLik, and an extractAIC method.

References

Agresti, A. (2002) Categorical Data. Second edition. Wiley.

Examples

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

## More manageable data set:
dat <- subset(soup, as.numeric(as.character(RESP)) <=  24)
dat$RESP <- dat$RESP[drop=TRUE]

m1 <- clmm(SURENESS ~ PROD, random = RESP, data = dat, link="probit",
           Hess = TRUE, method="ucminf", threshold = "symmetric")

m1
summary(m1)
logLik(m1)
vcov(m1)
extractAIC(m1)
anova(m1, update(m1, location = SURENESS ~ 1, Hess = FALSE))
anova(m1, update(m1, random = NULL))

update(mms1, Hess = FALSE, nAGQ = 3)

## Binomial example with data from the lme4-package:
data(cbpp, package = "lme4")
cbpp2 <- rbind(cbpp[,-(2:3)], cbpp[,-(2:3)])
cbpp2 <- within(cbpp2, {
    incidence <- as.factor(rep(0:1, each=nrow(cbpp)))
    freq <- with(cbpp, c(incidence, size - incidence))
})

## Fit with Laplace approximation:
fm1 <- clmm(incidence ~ period, random = herd, weights = freq,
            data = cbpp2, Hess = 1)
summary(fm1)

## Fit with the adaptive Gauss-Hermite quadrature approximation:
fm2 <- clmm(incidence ~ period, random = herd, weights = freq,
            data = cbpp2, Hess = 1, nAGQ = 7)
summary(fm2)