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