clmm2_tutorial.R
In ordinal: Regression Models for Ordinal Data

### R code from vignette source 'clmm2_tutorial.Rnw'

###################################################
### code chunk number 1: Initialize
###################################################

## 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=" ")



###################################################
### code chunk number 2: clmm2_tutorial.Rnw:152-155
###################################################
data(wine)
head(wine)
str(wine)


###################################################
### code chunk number 3: clmm2_tutorial.Rnw:176-190
###################################################
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)


###################################################
### code chunk number 4: clmm2_tutorial.Rnw:217-219
###################################################
fm1 <- clmm2(rating ~ temp + contact, random=judge, data=wine)
fm1


###################################################
### code chunk number 5: clmm2_tutorial.Rnw:226-229
###################################################
fm2 <- clmm2(rating ~ temp + contact, random=judge, data=wine,
            Hess=TRUE, nAGQ=10)
summary(fm2)


###################################################
### code chunk number 6: clmm2_tutorial.Rnw:265-266
###################################################
exp(coef(fm2)[5])


###################################################
### code chunk number 7: clmm2_tutorial.Rnw:274-276
###################################################
fm3 <- clmm2(rating ~ temp, random=judge, data=wine, nAGQ=10)
anova(fm3, fm2)


###################################################
### code chunk number 8: clmm2_tutorial.Rnw:282-284
###################################################
fm4 <- clm2(rating ~ temp + contact, data=wine)
anova(fm4, fm2)


###################################################
### code chunk number 9: clmm2_tutorial.Rnw:295-297
###################################################
pr2 <- profile(fm2, range=c(.1, 4), nSteps=30, trace=0)
confint(pr2)


###################################################
### code chunk number 10: profilePlot
###################################################
getOption("SweaveHooks")[["fig"]]()
plot(pr2)


###################################################
### code chunk number 11: profileFig
###################################################
getOption("SweaveHooks")[["fig"]]()
plot(pr2)


###################################################
### code chunk number 12: ranefPlot
###################################################
getOption("SweaveHooks")[["fig"]]()
ci <- fm2$ranef + qnorm(0.975) * sqrt(fm2$condVar) %o% c(-1, 1)
ord.re <- order(fm2$ranef)
ci <- ci[order(fm2$ranef),]
plot(1:9, fm2$ranef[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)


###################################################
### code chunk number 13: clmm2_tutorial.Rnw:348-349
###################################################
getOption("SweaveHooks")[["fig"]]()
ci <- fm2$ranef + qnorm(0.975) * sqrt(fm2$condVar) %o% c(-1, 1)
ord.re <- order(fm2$ranef)
ci <- ci[order(fm2$ranef),]
plot(1:9, fm2$ranef[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)


###################################################
### code chunk number 14: clmm2_tutorial.Rnw:361-362
###################################################
head(cbind(wine, fitted(fm2)))


###################################################
### code chunk number 15: clmm2_tutorial.Rnw:367-368
###################################################
head(cbind(wine, pred=predict(fm2, newdata=wine)))


###################################################
### code chunk number 16: clmm2_tutorial.Rnw:386-388
###################################################
plogis(fm2$Theta[3] - fm2$beta[2]) -
  plogis(fm2$Theta[2] - fm2$beta[2])


###################################################
### code chunk number 17: clmm2_tutorial.Rnw:396-397
###################################################
qnorm(0.95) * c(-1, 1) * fm2$stDev


###################################################
### code chunk number 18: clmm2_tutorial.Rnw:402-410
###################################################
pred <-
  function(eta, theta, cat = 1:(length(theta)+1), inv.link = plogis)
{
  Theta <- c(-1e3, theta, 1e3)
  sapply(cat, function(j)
         inv.link(Theta[j+1] - eta) - inv.link(Theta[j] - eta) )
}
pred(qnorm(0.05) * fm2$stDev, fm2$Theta)


###################################################
### code chunk number 19: clmm2_tutorial.Rnw:416-434
###################################################
mat <- expand.grid(judge = qnorm(0.95) * c(-1, 0, 1) * fm2$stDev,
                   contact = c(0, fm2$beta[2]),
                   temp = c(0, fm2$beta[1]))
pred.mat <- pred(eta=rowSums(mat), theta=fm2$Theta)
lab <- paste("contact=", rep(levels(wine$contact), 2), ", ",
             "temp=", rep(levels(wine$temp), each=2), sep="")
par(mfrow=c(2, 2))
for(k in c(1, 4, 7, 10)) {
  plot(1:5, pred.mat[k,], lty=2, type = "l", ylim=c(0,1),
       xlab="Bitterness rating scale", axes=FALSE,
       ylab="Probability", main=lab[ceiling(k/3)], las=1)
  axis(1); axis(2)
  lines(1:5, pred.mat[k+1, ], lty=1)
  lines(1:5, pred.mat[k+2, ], lty=3)
  legend("topright",
         c("avg. judge", "5th %-tile judge", "95th %-tile judge"),
         lty=1:3, bty="n")
}


###################################################
### code chunk number 20: clmm2_tutorial.Rnw:439-449
###################################################
getOption("SweaveHooks")[["fig"]]()
k <- 1
plot(1:5, pred.mat[k,], lty=2, type = "l", ylim=c(0,1),
     xlab="Bitterness rating scale", axes=FALSE,
     ylab="Probability", main=lab[ceiling(k/3)], las=1)
axis(1); axis(2)
lines(1:5, pred.mat[k+1, ], lty=1)
lines(1:5, pred.mat[k+2, ], lty=3)
legend("topright",
       c("avg. judge", "5th %-tile judge", "95th %-tile judge"),
       lty=1:3, bty="n")


###################################################
### code chunk number 21: clmm2_tutorial.Rnw:451-461
###################################################
getOption("SweaveHooks")[["fig"]]()
k <- 4
plot(1:5, pred.mat[k,], lty=2, type = "l", ylim=c(0,1),
     xlab="Bitterness rating scale", axes=FALSE,
     ylab="Probability", main=lab[ceiling(k/3)], las=1)
axis(1); axis(2)
lines(1:5, pred.mat[k+1, ], lty=1)
lines(1:5, pred.mat[k+2, ], lty=3)
legend("topright",
       c("avg. judge", "5th %-tile judge", "95th %-tile judge"),
       lty=1:3, bty="n")


###################################################
### code chunk number 22: clmm2_tutorial.Rnw:463-473
###################################################
getOption("SweaveHooks")[["fig"]]()
k <- 7
plot(1:5, pred.mat[k,], lty=2, type = "l", ylim=c(0,1),
     xlab="Bitterness rating scale", axes=FALSE,
     ylab="Probability", main=lab[ceiling(k/3)], las=1)
axis(1); axis(2)
lines(1:5, pred.mat[k+1, ], lty=1)
lines(1:5, pred.mat[k+2, ], lty=3)
legend("topright",
       c("avg. judge", "5th %-tile judge", "95th %-tile judge"),
       lty=1:3, bty="n")


###################################################
### code chunk number 23: clmm2_tutorial.Rnw:475-485
###################################################
getOption("SweaveHooks")[["fig"]]()
k <- 10
plot(1:5, pred.mat[k,], lty=2, type = "l", ylim=c(0,1),
     xlab="Bitterness rating scale", axes=FALSE,
     ylab="Probability", main=lab[ceiling(k/3)], las=1)
axis(1); axis(2)
lines(1:5, pred.mat[k+1, ], lty=1)
lines(1:5, pred.mat[k+2, ], lty=3)
legend("topright",
       c("avg. judge", "5th %-tile judge", "95th %-tile judge"),
       lty=1:3, bty="n")


###################################################
### code chunk number 24: clmm2_tutorial.Rnw:495-496
###################################################
exp(2*qnorm(0.95) * fm2$stDev)


###################################################
### code chunk number 25: clmm2_tutorial.Rnw:502-503
###################################################
exp(2*qnorm(0.75) * fm2$stDev)


###################################################
### code chunk number 26: misc (eval = FALSE)
###################################################
##