options(help_type="html")
require(ordinal)

## In old versions of ordinal you might need:
data(wine)

## Look at the wine data:
head(wine)
str(wine)
## Look at the manual page for the wine data:
## ?wine

#################################
## Fit a basic CLM:
fm1 <- clm(rating ~ temp + contact, data=wine, link="probit")
summary(fm1)
confint(fm1)
## Effect of the variables:
drop1(fm1, test="Chi")

# (some of the) ANODE table info: (uses logit link) 

fmlogit0 <- clm(rating ~ temp*contact, data=wine)
fmlogit1 <- clm(rating ~ temp + contact, data=wine)
fmlogit1a <- clm(rating ~ temp , data=wine)
fmlogit1b <- clm(rating ~ contact , data=wine)
fmlogit2 <- clm(rating ~ 1 , data=wine)

# Test for overall treatment effect:
anova(fmlogit0, fmlogit2)

# Test for TxC interaction:
anova(fmlogit0, fmlogit1)

# Test for main effect of C:
anova(fmlogit1, fmlogit1a)

# Test for main effect of T:
anova(fmlogit1, fmlogit1b)


#################################
## Equidistant thresholds:
fm2 <- clm(rating ~ temp + contact, data=wine, link="probit",
           threshold="equi")
summary(fm2)
## Compare to fm1 - can we reject equi-distant thresholds ?:
anova(fm2, fm1) # No

## Symmetric thresholds:
fm3 <- clm(rating ~ temp + contact, data=wine, link="probit",
           threshold="symmetric")
summary(fm3)

## Compare fm1 and fm3  - can we reject symmetric thresholds ?
anova(fm1, fm3) # no, but close

#################################
## Scale effects:

fm4 <- clm(rating ~ temp + contact, scale= ~ temp, data=wine,
           link="probit")
summary(fm4)

## Is the scale effect needed in our model?
anova(fm1, fm4) ## no

## Automatic test of scale effects for all variables:
scale_test(fm1)

#################################
## Mixed effects models:

## Random judge effects:
fm5 <- clmm(rating ~ temp + contact + (1|judge), data=wine,
            link="probit")
summary(fm5)

## Do we need judges in the model?
anova(fm1, fm5)

## More precise evaluation of the likelihood function:
fm6 <- clmm(rating ~ temp + contact + (1|judge), data=wine,
            link="probit", nAGQ=10)
summary(fm6)
## which does not make any difference here....

## Adding random bottle effects:
fm6 <- clmm(rating ~ temp + contact + (1|judge) + (1|bottle),
            data=wine, link="probit")
summary(fm6)

#################################
## Accessing random effects:

## Judge-specific effects:
(re <- ranef(fm5))
## their (conditional) variances (squared "standard-errors"):
condVar(fm5)
## Standard errors:
(se.re <- condVar(fm5)$judge[, 1]^.5)

## Plotting the judge performances:
ci <- unlist(ranef(fm5)) + qnorm(0.975) * se.re %o% c(-1, 1)
ord.re <- order(unlist(re))
ci <- ci[ord.re,]
plot(1:9, unlist(re)[ord.re], axes=FALSE, ylim=range(ci),
     xlab="Judge", ylab="Judge effect")
axis(1, at=1:9, labels = ord.re)
axis(2)
for(i in 1:9) segments(i, ci[i,1], i, ci[i, 2])
abline(h = 0, lty=2)
#################################

## Want to know more?

## - 3 Tutorials at http://cran.r-project.org/package=ordinal explains
##   more.
## - Look at the applications in the papers referenced in the slides.
## - Look at the source-code (R and C code) at
##   http://r-forge.r-project.org/R/?group_id=800