\documentclass[a4paper]{article}
\usepackage{amsmath}%the AMS math extension of LaTeX.
\usepackage{amssymb}%the extended AMS math symbols.
%% \usepackage{amsthm}
\usepackage{bm}%Use 'bm.sty' to get `bold math' symbols
\usepackage{natbib}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage{Sweave}
\usepackage{url}
\usepackage{float}%Use `float.sty'
\usepackage[left=3.5cm,right=3.5cm]{geometry}
\usepackage{algorithmic}
\usepackage[amsmath,thmmarks,standard,thref]{ntheorem}

%%\VignetteIndexEntry{clm tutorial}
%%\VignetteDepends{ordinal, xtable}
\title{A Tutorial on fitting Cumulative Link Models with
  the \textsf{ordinal} Package}
\author{Rune Haubo B Christensen}

%% \numberwithin{equation}{section}
\setlength{\parskip}{2mm}%.8\baselineskip}
\setlength{\parindent}{0in}

%%  \DefineVerbatimEnvironment{Sinput}{Verbatim}%{}
%%  {fontshape=sl, xleftmargin=1em}
%%  \DefineVerbatimEnvironment{Soutput}{Verbatim}%{}
%%  {xleftmargin=1em}
%%  \DefineVerbatimEnvironment{Scode}{Verbatim}%{}
%%  {fontshape=sl, xleftmargin=1em}
\fvset{listparameters={\setlength{\topsep}{0pt}}}
%% \fvset{listparameters={\setlength{\botsep}{0pt}}}
\renewenvironment{Schunk}{\vspace{-1mm}}{\vspace{-1mm}}

%RE-DEFINE marginpar
\setlength{\marginparwidth}{1in}
\let\oldmarginpar\marginpar
\renewcommand\marginpar[1]{\oldmarginpar[\-\raggedleft\tiny #1]%
{\tiny #1}}
%uncomment to _HIDE_MARGINPAR_:
%\renewcommand\marginpar[1]{}

\newcommand{\var}{\textup{var}}
\newcommand{\I}{\mathcal{I}}
\newcommand{\bta}{\bm \theta}
\newcommand{\ta}{\theta}
\newcommand{\tah}{\hat \theta}
\newcommand{\di}{~\textup{d}}
\newcommand{\td}{\textup{d}}
\newcommand{\Si}{\Sigma}
\newcommand{\si}{\sigma}
\newcommand{\bpi}{\bm \pi}
\newcommand{\bmeta}{\bm \eta}
\newcommand{\tdots}{\hspace{10mm} \texttt{....}}
\newcommand{\FL}[1]{\fvset{firstline= #1}}
\newcommand{\LL}[1]{\fvset{lastline= #1}}
\newcommand{\s}{\square}
\newcommand{\bs}{\blacksquare}

% figurer bagerst i artikel
%% \usepackage[tablesfirst, nolists]{endfloat}
%% \renewcommand{\efloatseparator}{\vspace{.5cm}}

\theoremstyle{plain} %% {break}
\theoremseparator{:}
\theoremsymbol{{\tiny $\square$}}
%%\theoremstyle{plain}
\theorembodyfont{\small}
\theoremindent5mm
\renewtheorem{example}{Example}

%% \newtheoremstyle{example}{\topsep}{\topsep}%
%% {}%         Body font
%% {}%         Indent amount (empty = no indent, \parindent = para indent)
%% {\bfseries}% Thm head font
%% {}%        Punctuation after thm head
%% {\newline}%     Space after thm head (\newline = linebreak)
%% {\thmname{#1}\thmnumber{ #2}\thmnote{ #3}}%         Thm head spec
%%
%% \theoremstyle{example}
%% %% \newtheorem{example}{Example}[subsection]
%% \newtheorem{example}{Example}[section]

\usepackage{lineno}
% \linenumbers
\newcommand*\patchAmsMathEnvironmentForLineno[1]{%
\expandafter\let\csname old#1\expandafter\endcsname\csname #1\endcsname
\expandafter\let\csname oldend#1\expandafter\endcsname\csname end#1\endcsname
\renewenvironment{#1}%
{\linenomath\csname old#1\endcsname}%
{\csname oldend#1\endcsname\endlinenomath}}%
\newcommand*\patchBothAmsMathEnvironmentsForLineno[1]{%
 \patchAmsMathEnvironmentForLineno{#1}%
 \patchAmsMathEnvironmentForLineno{#1*}}%
\AtBeginDocument{%
\patchBothAmsMathEnvironmentsForLineno{equation}%
\patchBothAmsMathEnvironmentsForLineno{align}%
\patchBothAmsMathEnvironmentsForLineno{flalign}%
\patchBothAmsMathEnvironmentsForLineno{alignat}%
\patchBothAmsMathEnvironmentsForLineno{gather}%
\patchBothAmsMathEnvironmentsForLineno{multline}%
}

\begin{document}
\bibliographystyle{chicago}
\maketitle

\begin{abstract}

  It is shown by example how a cumulative link mixed model is fitted
  with the \texttt{clm} function in package \textsf{ordinal}. Model
  interpretation and inference is briefly discussed.

\end{abstract}

%% \newpage
\tableofcontents
\newpage

\SweaveOpts{echo=TRUE, results=verb, width=4.5, height=4.5}
\SweaveOpts{prefix.string=figs}
\fvset{listparameters={\setlength{\topsep}{0pt}}, gobble=0, fontsize=\small}
%% \fvset{gobble=0, fontsize=\small}
\setkeys{Gin}{width=.49\textwidth}

<<Initialize, echo=FALSE, results=hide>>=

## Load common packages, functions and set settings:
library(ordinal)
library(xtable)
##
RUN <- FALSE    #redo computations and write .RData files
## Change options:
op <- options() ## To be able to reset settings
options("digits" = 7)
options(help_type = "html")
## options("width" = 75)
options("SweaveHooks" = list(fig=function()
        par(mar=c(4,4,.5,0)+.5)))
options(continue=" ")

@

\section{Introduction}
\label{sec:introduction}

We will consider the data on the bitterness of wine from
\citet{randall89} presented in Table~\ref{tab:winedata} and available
as the object \texttt{wine} in package \textsf{ordinal}. The data were
also analyzed with mixed effects models by \citet{tutz96}.
The following gives an impression of the wine data object:
<<>>=
data(wine)
head(wine)
str(wine)
@
The data represent a factorial experiment on factors determining the
bitterness of wine with 1 = ``least bitter'' and 5 = ``most bitter''.
Two treatment factors (temperature and contact)
each have two levels. Temperature and contact between juice and
skins can be controlled when crushing grapes during wine
production. Nine judges each assessed wine from two bottles from
each of the four treatment conditions, hence there are 72
observations in all. For more information see the manual entry for the
wine data: \texttt{help(wine)}.

\begin{table}
  \centering
  \caption{Ratings of the bitterness of some white wines. Data are
    adopted from \citet{randall89}.}
  \label{tab:winedata}
  \begin{tabular}{lllrrrrrrrrr}
    \hline
    & & & \multicolumn{9}{c}{Judge} \\
    \cline{4-12}
<<echo=FALSE, results=tex>>=
data(wine)
temp.contact.bottle <- with(wine, temp:contact:bottle)[drop=TRUE]
tab <- xtabs(as.numeric(rating) ~ temp.contact.bottle + judge,
             data=wine)
class(tab) <- "matrix"
attr(tab, "call") <- NULL
mat <- cbind(rep(c("cold", "warm"), each = 4),
             rep(rep(c("no", "yes"), each=2), 2),
             1:8, tab)
colnames(mat) <-
  c("Temperature", "Contact", "Bottle", 1:9)
xtab <- xtable(mat)
print(xtab, only.contents=TRUE, include.rownames=FALSE,
      sanitize.text.function = function(x) x)
@
\end{tabular}
\end{table}

The intention with this tutorial is not to explain cumulative link
models in details, rather the main aim is to briefly cover the main
functions and methods in the \textsf{ordinal} package to analyze
ordinal data. A more thorough introduction to cumulative link models
and the \textsf{ordinal} package is given in
\citet{christensen11}; a book length treatment of ordinal data
analysis is given by \citet{agresti10} although not related to the
\textsf{ordinal} package.

\section{Fitting Cumulative Link Models}
\label{sec:fitt-cumul-link}

We will fit the following cumulative link model to the wine
data:
\begin{equation}
  \label{eq:CLM}
  \begin{array}{c}
    \textup{logit}(P(Y_i \leq j)) = \theta_j - \beta_1 (\mathtt{temp}_i)
    - \beta_2(\mathtt{contact}_i) \\
    i = 1,\ldots, n, \quad j = 1, \ldots, J-1
  \end{array}
\end{equation}
This is a model for the cumulative probability of the $i$th rating
falling in the $j$th category or below, where $i$ index
all observations ($n=72$) and $j = 1, \ldots, J$ index the response
categories ($J = 5$). The $theta_j$ parameter is the intercept
for the $j$th cumulative logit $\textup{logit}(P(Y_i \leq j))$; they
are known as threshold parameters, intercepts or cut-points.

This model is also known as the \emph{proportional odds model}, a
\emph{cumulative logit model}, and an \emph{ordered logit model}.

We fit this cumulative link model by maximum likelihood with the
\texttt{clm} function in package \textsf{ordinal}. Here we save the
fitted \texttt{clm} model
in the object \texttt{fm1} (short for \texttt{f}itted \texttt{m}odel
\texttt{1}) and \texttt{print} the model by simply typing its name:
<<>>=
fm1 <- clm(rating ~ temp + contact, data=wine)
fm1
@
Additional information is provided with the \texttt{summary} method:
<<>>=
summary(fm1)
@
The primary result is the coefficient table with parameter estimates,
standard errors and Wald (or normal) based $p$-values for tests of the
parameters being zero.
The maximum likelihood estimates of the parameters are:
\begin{equation}
  \label{eq:parameters}
  \hat\beta_1 = 2.50, ~~\hat\beta_2 = 1.53,
  ~~\{\hat\theta_j\} = \{-1.34,~ 1.25,~ 3.47,~ 5.01\}.
\end{equation}

The number of Newton-Raphson iterations is given below \texttt{niter}
with the number of step-halvings in parenthesis.
\texttt{max.grad} is the maximum absolute gradient of the
log-likelihood function with respect to the parameters. A small
absolute gradient is a necessary condition for convergence of the
model. The iterative procedure will declare convergence whenever the
maximum absolute gradient is below
<<>>=
clm.control()$gradTol
@
which may be altered --- see \texttt{help(clm.control)}.

The condition number of the Hessian (\texttt{cond.H}) measures the empirical
identifiability of the model. High numbers, say larger than $10^4$ or
$10^6$ indicate that the model is ill defined. This could indicate that
the model can be simplified, that possibly some parameters are not
identifiable, and that optimization of the model can be difficult. In
this case the condition number of the Hessian does not indicate a
problem with the model.

The coefficients for \texttt{temp} and \texttt{contact} are positive
indicating that higher temperature and more contact increase the
bitterness of wine, i.e., rating in higher categories is more likely.
The odds ratio of the event $Y \geq j$ is
$\exp(\beta_{\textup{treatment}})$, thus the odds ratio of bitterness
being rated in
category $j$ or above at warm relative to cold temperatures is
<<>>=
exp(coef(fm1)[5])
@

The $p$-values for the location coefficients provided by the
\texttt{summary} method are based on the so-called Wald
statistic. More accurate tests are provided by likelihood ratio
tests. These can be obtained with the \texttt{anova} method, for
example, the likelihood ratio test of \texttt{contact} is
<<>>=
fm2 <- clm(rating ~ temp, data=wine)
anova(fm2, fm1)
@
which in this case produces a slightly lower $p$-value.
Equivalently we can use \texttt{drop1} to obtain likelihood ratio
tests of the explanatory variables while \emph{controlling} for the
remaining variables:
<<>>=
drop1(fm1, test = "Chi")
@
Likelihood ratio tests of the explanatory variables while
\emph{ignoring} the remaining variables are provided by the
\texttt{add1} method:
<<>>=
fm0 <- clm(rating ~ 1, data=wine)
add1(fm0, scope = ~ temp + contact, test = "Chi")
@
In this case these latter tests are not as strong as the tests
controlling for the other variable.

Confidence intervals are provided by the \texttt{confint} method:
<<>>=
confint(fm1)
@
These are based on the profile likelihood function and generally
fairly accurate. Less accurate, but simple and symmetric confidence
intervals based on the standard errors of the parameters (so-called
Wald confidence intervals) can be obtained with
<<>>=
confint(fm1, type="Wald")
@

In addition to the logit link, which is the default, the probit,
log-log, complementary log-log and cauchit links are also
available. For instance, a proportional hazards model for grouped
survival times is fitted using the complementary log-log link:
<<>>=
fm.cll <- clm(rating ~ contact + temp, data=wine, link="cloglog")
@

The cumulative link model in \eqref{eq:CLM} assume that the
thresholds, $\{\theta_j\}$ are constant for all values of the
remaining explanatory variables, here \texttt{temp} and
\texttt{contact}. This is generally referred to as the
\emph{proportional odds assumption} or \emph{equal slopes
  assumption}.
We can relax that
assumption in two general ways: with nominal effects and scale effects
which we will now discuss in turn.

\section{Nominal Effects}
\label{sec:nominal-effects}

The CLM in \eqref{eq:CLM} specifies a structure in which the
regression parameters, $\bm\beta$ are not allowed to vary with $j$.
Nominal effects relax this assumption by allowing one or more
regression parameters to vary with $j$. In the following model we
allow the regression parameter for \texttt{contact} to vary with $j$:
\begin{equation}
  \label{eq:CLM_nominal}
  \begin{array}{c}
    \textup{logit}(P(Y_i \leq j)) = \theta_j - \beta_1 (\mathtt{temp}_i)
    - \beta_{2j} (\mathtt{contact}_i) \\
    i = 1,\ldots, n, \quad j = 1, \ldots, J-1
  \end{array}
\end{equation}
This means that there is one estimate of $\beta_2$ for each $j$.
This model is specified as follows with \texttt{clm}:
<<>>=
fm.nom <- clm(rating ~ temp, nominal=~contact, data=wine)
summary(fm.nom)
@
As can be seen from the output of \texttt{summary} there is no
regression coefficient estimated for \texttt{contact}, but there are
two sets of threshold parameters estimated.

The first five threshold
parameters have \texttt{.(Intercept)} appended their names indicating
that these are the estimates of $\theta_j$.
The following five threshold parameters
have \texttt{.contactyes} appended their name indicating that these
parameters are differences between the threshold parameters at the two
levels of contact. This interpretation corresponds to the default
treatment contrasts; if other types of contrasts are specified, the
interpretation is a little different. As can be seen from the output,
the effect of \texttt{contact} is almost constant across thresholds
and around $1.5$ corresponding to the estimate from \texttt{fm1} on
page \pageref{eq:parameters}, so probably there is not much evidence
that the effect of \texttt{contact} varies with $j$

We can perform a likelihood ratio test of the equal slopes or
proportional odds assumption for \texttt{contact} by comparing the
likelihoods of models \eqref{eq:CLM} and \eqref{eq:CLM_nominal} as
follows:
<<>>=
anova(fm1, fm.nom)
@
There is only little difference in the log-likelihoods of the two
models, so the test is insignificant. There is therefore no evidence
that the proportional odds assumption is violated for
\texttt{contact}.

It is not possible to estimate both $\beta_2$ and $\beta_{2j}$ in the
same model. Consequently variables that appear in \texttt{nominal}
cannot enter elsewhere as well. For instance not all parameters are
identifiable in the following model:
<<>>=
fm.nom2 <- clm(rating ~ temp + contact, nominal=~contact, data=wine)
@
We are made aware of this when summarizing or printing the model:
<<>>=
summary(fm.nom2)
@

\section{Scale Effects}
\label{sec:scale-effects}

Scale effects are usually motivated from the latent variable
interpretation of a CLM. Assume the following model for a latent
variable:
\begin{equation}
  \label{eq:latent}
  S_i = \alpha^* + \bm x_i^T \bm\beta^* + \varepsilon_i, \quad
  \varepsilon_i \sim N(0, \sigma^{*^2})
\end{equation}
If the ordinal variable $Y_i$ is observed such that $Y_i = j$ is
recorded if $\theta_{j-1}^* < S_i \leq \theta_j^*$, where
\begin{equation}
  \label{eq:thresholds}
  -\infty \equiv \theta_0^* < \theta_1^* < \ldots < \theta^*_{J-1} <
  \theta_{J}^* \equiv \infty
\end{equation}
then we have the cumulative link model for $Y_i$:
\begin{equation*}
  P(Y_i \leq j) = \Phi ( \theta_j - \bm x_i^T \bm\beta )
\end{equation*}
where we have used $\theta_j = (\theta_j^* + \alpha^*) / \sigma^*$ and
$\bm\beta = \bm\beta^* / \sigma^*$ (parameters with a ``$^*$'' exist on
the latent scale, while those without are identifiable),
and $\Phi$ is the inverse probit
link and denotes the standard normal CDF.
Other assumptions on the distribution of the latent variable, $S_i$
lead to other link functions, e.g., a logistic distribution for the
latent variable corresponds to a logit link.

If the scale (or dispersion) of the latent distribution is described
by a log-linear model such that $\log( \sigma_i) = \bm z_i^T \bm
\zeta$ (equivalently $\sigma_i = \exp(\bm z_i^T \bm \zeta)$; also note
that $\bm z_i^T \bm \zeta$ is a linear model just like $\bm x_i^T \bm
\beta$), then the resulting CLM reads (for more details, see e.g.,
\citet{christensen11} or \citet{agresti10}):
\begin{equation}
  \label{eq:CLM_scale}
  P(Y_i \leq j) = \Phi \left( \frac{\theta_j - \bm x_i^T \bm\beta}
    {\sigma_i} \right)
\end{equation}
The conventional link functions in cumulative link models correspond
to distributions for the latent variable that are members of the
location-scale family of distributions
(cf. \url{http://en.wikipedia.org/wiki/Location-scale_family}).
They have the common form $F(\mu,
\sigma)$, where $\mu$ denotes the location of the distribution and
$\sigma$ denotes the scale of the distribution. For instance in the
normal distribution (probit link) $\mu$ is the mean, and $\sigma$ is
the spread, while in the logistic distribution (logit link), $\mu$ is
the mean and $\sigma\pi / \sqrt 3$ is the spread
(cf. \url{http://en.wikipedia.org/wiki/Logistic_distribution}).

Thus allowing for scale effects corresponds to modelling not only the
location of the latent distribution, but also the scale.
Just as the absolute location ($\alpha^*$) is not identifiable, the
absolute scale ($\sigma^*$) is not identifiable either in the CLM,
however \emph{differences} in location modelled with $\bm x_i^T
\bm\beta$ and \emph{ratios} of scale modelled with $\exp(\bm z_i^T
\bm\zeta)$ are identifiable.

Now we turn to our running example and fit a model where we allow the
scale of the latent distribution to depend on
temperature:
\begin{equation}
  \label{eq:CLM_scale_wine}
  \textup{logit}(P(Y_i \leq j)) = \frac{\theta_j - \beta_1 (\mathtt{temp}_i)
    - \beta_{2} (\mathtt{contact}_i)} {\exp( \zeta_1 (\mathtt{temp}_i))}
\end{equation}
\begin{equation*}
  i = 1,\ldots, n, \quad j = 1, \ldots, J-1
\end{equation*}
We can estimate this model with
<<>>=
fm.sca <- clm(rating ~ temp + contact, scale=~temp, data=wine)
summary(fm.sca)
@
Notice that both location and scale effects of \texttt{temp} are
identifiable. Also notice that the scale coefficient for \texttt{temp}
is given on the
log-scale, where the Wald test is more appropriate.
The absolute scale of the latent distribution is not estimable, but we
can estimate the scale at warm conditions relative to cold
conditions. Therefore the estimate of $\kappa$ in the relation
$\sigma_{warm} = \kappa \sigma_{cold}$ is given by
<<>>=
exp(fm.sca$zeta)
@
However, the scale difference is not significant in this case as
judged by the $p$-value in the summary output. \texttt{confint} and
\texttt{anova} apply with no change to models with scale, but
\texttt{drop1}, \texttt{add1} and \texttt{step} methods will only drop
or add terms to the (location) \texttt{formula} and not to
\texttt{scale}.

\section{Structured Thresholds}
\label{sec:struct-thresh}

In section~\ref{sec:nominal-effects} we relaxed the assumption that
regression parameters have the same effect across all thresholds. In
this section we will instead impose additional restrictions on the
thresholds. In the following model we require that the thresholds,
$\theta_j$ are equidistant or equally spaced ($\theta_j - \theta_{j-1}
= \textup{constant}$ for $j = 2, ..., J-1$):
<<>>=
fm.equi <- clm(rating ~ temp + contact, data=wine,
               threshold="equidistant")
summary(fm.equi)
@
The parameters determining the thresholds are now the first threshold
and the spacing among consecutive thresholds. The mapping to this
parameterization is stored in the transpose of the Jacobian
(\texttt{tJac}) component of the model fit. This makes it possible to
get the thresholds imposed by the equidistance structure with
<<>>=
c(with(fm.equi, tJac %*% alpha))
@
The following shows that the distance between consecutive thresholds
in \texttt{fm1} is very close to the \texttt{spacing} parameter from
\texttt{fm.equi}:
<<>>=
mean(diff(fm1$alpha))
@

The gain in imposing additional restrictions on the thresholds is the
use of fewer parameters. Whether the restrictions are warranted by the
data can be tested in a likelihood ratio test:
<<>>=
anova(fm1, fm.equi)
@
In this case the test is non-significant, so there is no considerable
loss of fit at the gain of saving two parameters, hence we may retain
the model with equally spaced thresholds.


\section{Predictions}
\label{sec:predictions}

Fitted values are extracted with e.g., \texttt{fitted(fm1)} and
produce fitted probabilities, i.e., the $i$th fitted probability is the
probability that the $i$th observation falls in the response category
that it did. The predictions of which response class the observations
would be most likely to fall in given the model are obtained with:
<<>>=
predict(fm1, type = "class")
@

Say we just wanted the predictions for the four combinations of
\texttt{temp} and \texttt{contact}. The probability that an
observation falls in each of the five response categories based on the
fitted model is given by:
<<>>=
newData <- expand.grid(temp=levels(wine$temp),
                       contact=levels(wine$contact))
cbind(newData, predict(fm1, newdata=newData)$fit)
@

Standard errors and confidence intervals of predictions are also
available, for example, the predictions, standard errors and 95\%
confidence intervals are given by (illustrated here using
\texttt{do.call} for the first six observations):
<<>>=
head(do.call("cbind", predict(fm1, se.fit=TRUE, interval=TRUE)))
@
The confidence level can be set with the \texttt{level} argument and
other types of predictions are available with the \texttt{type}
argument.

\section{Infinite Parameter Estimates}
\label{sec:infin-param-estim}

If we attempt to test the proportional odds assumption for
\texttt{temp}, some peculiarities show up:
<<>>=
fm.nom2 <- clm(rating ~ contact, nominal=~temp, data=wine)
summary(fm.nom2)
@
Several of the threshold coefficients are extremely large with huge
standard errors. Also the condition number of the Hessian is very
large and a larger number of iterations was used all indicating that
something is not well-behaved. The problem is that the the ML estimates of
some of the threshold parameters are at (plus/minus)
infinity. \texttt{clm} is not able to detect this and just stops
iterating when the likelihood function is flat enough for the
termination criterion to be satisfied, i.e., when the maximum absolute
gradient is small enough.

Even though some parameter estimates are not at (plus/minus) infinity
while they should have been, the remaining parameters are accurately
determined and the value of the log-likelihood is also accurately
determined. This means that likelihood ratio tests are still
available, for example, it is still possible to test the proportional
odds assumption for \texttt{temp}:
<<>>=
anova(fm1, fm.nom2)
@

\section{Unidentified parameters}
\label{sec:unid-param}

In the following example (using another data set) one parameter is not
identifiable:
<<>>=
data(soup)
fm.soup <- clm(SURENESS ~ PRODID * DAY, data=soup)
summary(fm.soup)
@
The undefined parameter shows up as \texttt{NA} in the coefficient
table. The
source of the singularity is revealed in the following table:
<<>>=
with(soup, table(DAY, PRODID))
@
which shows that the third \texttt{PRODID} was not presented at the
second day at all. The design matrix will in this case be column rank
deficient (also referred to as singular). This is detected by
\texttt{clm} using the \texttt{drop.coef} function from
\textsf{ordinal}. The following illustrates that the column rank of
the design matrix is less than its number of columns:
<<>>=
mm <- model.matrix( ~ PRODID * DAY, data=soup)
ncol(mm)
qr(mm, LAPACK = FALSE)$rank
@

A similar type of rank deficiency occurs when variables in
\texttt{nominal} are also present in \texttt{formula} or
\texttt{scale} as illustrated in section \ref{sec:nominal-effects}.

\section{Assessment of Model Convergence}
\label{sec:assessm-model-conv}

The maximum likelihood estimates of the parameters in cumulative link
models do not have closed form expressions, so iterative methods have
to be applied to fit the models. Such iterative methods can fail to
converge simply because an optimum cannot be found or because the
parameter estimates are not determined accurately enough.

An optimum has been found if the maximum absolute gradient is small
and if the condition number of the Hessian (defined here as the ratio
of the largest to the smallest eigenvalues of the Hessian evaluated at
convergence) is positive and not very large, say smaller than $10^4$ or
$10^6$. The maximum absolute gradient (\texttt{max.grad}) and the
condition number of the Hessian (\texttt{cond.H}) are reported by the
summary method, for an example see page \pageref{eq:parameters}.
A large condition number of the Hessian is not necessarily a
problem, but it can be. A small condition number of the Hessian on the
other hand is a good insurance that an optimum has been reached.
Thus the maximum absolute gradient and the condition number
of the Hessian can be used to check if the optimization has reached a
well-defined optimum.

To determine the accuracy of the parameter estimates we use the
\texttt{convergence} method:
<<>>=
convergence(fm1)
@
The most important information is the number of correct decimals
(\texttt{Cor.Dec}) and the number of significant digits
(\texttt{Sig.Dig}) with which the parameters are determined. In this
case all parameters are very accurately determined, so there is no
reason to lower the convergence tolerance. The \texttt{logLik.error}
shows that the error in the reported value of the log-likelihood is
below $10^{-10}$, which is by far small enough that likelihood
ratio tests based on this model are accurate.

The convergence properties of the fitted model may be illustrated by
plotting slices of the log-likelihood function for the parameters. The
following code produce the slices in Figure~\ref{fig:slice1}.
<<>>=
slice.fm1 <- slice(fm1, lambda = 5)
par(mfrow = c(2, 3))
plot(slice.fm1)
@
The slices illustrates the log-likelihood function plotted as a
function each parameter in turn while the remaining parameters are
fixed at the ML estimates.
The \texttt{lambda} argument controls how far from the ML estimates
the slices should be computed; it can be interpreted as a multiplier
in curvature units, where a curvature unit is similar to a standard
error.

\setkeys{Gin}{width=.32\textwidth}
\begin{figure}
  \centering
<<slice11, echo=FALSE, results=hide, fig=TRUE, include=TRUE, width=3, height=3>>=
plot(slice.fm1, parm = 1)
@
<<slice12, echo=FALSE, results=hide, fig=TRUE, include=TRUE, width=3, height=3>>=
plot(slice.fm1, parm = 2)
@
<<slice13, echo=FALSE, results=hide, fig=TRUE, include=TRUE, width=3, height=3>>=
plot(slice.fm1, parm = 3)
@
<<slice14, echo=FALSE, results=hide, fig=TRUE, include=TRUE, width=3, height=3>>=
plot(slice.fm1, parm = 4)
@
<<slice15, echo=FALSE, results=hide, fig=TRUE, include=TRUE, width=3, height=3>>=
plot(slice.fm1, parm = 5)
@
<<slice16, echo=FALSE, results=hide, fig=TRUE, include=TRUE, width=3, height=3>>=
plot(slice.fm1, parm = 6)
@
\caption{Slices of the (negative) log-likelihood function for
  parameters in a
  model for the bitterness-of-wine data. Dashed lines indicate
  quadratic approximations to the log-likelihood function and
  vertical bars indicate maximum likelihood estimates.}
\label{fig:slice1}
\end{figure}

For an inspection of the log-likelihood function closer to the optimum
we can use a smaller \texttt{lambda}:
<<slice2, echo=TRUE, results=hide>>=
slice2.fm1 <- slice(fm1, lambda = 1e-5)
par(mfrow = c(2, 3))
plot(slice2.fm1)
@
The resulting figure is shown in Fig.~\ref{fig:slice2}.

\setkeys{Gin}{width=.32\textwidth}
\begin{figure}
  \centering
<<slice21, echo=FALSE, results=hide, fig=TRUE, include=TRUE, width=3, height=3>>=
plot(slice2.fm1, parm = 1)
@
<<slice22, echo=FALSE, results=hide, fig=TRUE, include=TRUE, width=3, height=3>>=
plot(slice2.fm1, parm = 2)
@
<<slice23, echo=FALSE, results=hide, fig=TRUE, include=TRUE, width=3, height=3>>=
plot(slice2.fm1, parm = 3)
@
<<slice24, echo=FALSE, results=hide, fig=TRUE, include=TRUE, width=3, height=3>>=
plot(slice2.fm1, parm = 4)
@
<<slice25, echo=FALSE, results=hide, fig=TRUE, include=TRUE, width=3, height=3>>=
plot(slice2.fm1, parm = 5)
@
<<slice26, echo=FALSE, results=hide, fig=TRUE, include=TRUE, width=3, height=3>>=
plot(slice2.fm1, parm = 6)
@
  \caption{Slices of the log-likelihood function for
    parameters in a model for the bitterness-of-wine data very close
    to the MLEs. Dashed lines indicate
    quadratic approximations to the log-likelihood function and
    vertical bars the indicate maximum likelihood estimates.}
\label{fig:slice2}
\end{figure}

\section{Profile Likelihood}
\label{sec:profile-likelihood}

The profile likelihood can be used for several things. Two of the most
important objectives are to provide accurate likelihood confidence
intervals and to illustrate effects of parameters in the fitted
model.

Confidence intervals based on the profile likelihood were already
obtained in section~\ref{sec:fitt-cumul-link} and will not be treated
any further here.

The effects of \texttt{contact} and \texttt{temp} can be illustrated
with
<<profileLikelihood, echo=TRUE>>=
pr1 <- profile(fm1, alpha=1e-4)
plot(pr1)
@
and provided in Figure~\ref{fig:ProfileLikelihood}.
The \texttt{alpha} argument is the significance level controling how
far from the maximum likelihood estimate the likelihood function
should be profiled. Learn more about the arguments to \texttt{profile}
with \texttt{help(profile.clm)}.
From the relative
profile likelihood for \texttt{tempwarm} we see that parameter values
between 1 and 4 are reasonably well supported by the data, and values
outside this range has little likelihood. Values between 2 and 3 are
very well supported by the data and all have high likelihood.

\setkeys{Gin}{width=.32\textwidth}
\begin{figure}
  \centering
<<prof1, echo=FALSE, results=hide, fig=TRUE, include=TRUE, width=3, height=3>>=
plot(pr1, which.par=1)
@
<<prof2, echo=FALSE, results=hide, fig=TRUE, include=TRUE, width=3, height=3>>=
plot(pr1, which.par=2)
@
  \caption{Relative profile likelihoods for the regression parameters
    in the Wine study. Horizontal lines indicate 95\% and 99\%
    confidence bounds.}
  \label{fig:ProfileLikelihood}
\end{figure}



\newpage
\bibliography{ordinal}
%% \newpage



\end{document}

<<misc, eval=FALSE, echo=FALSE>>=
@


