# R / Ordinal Scripts

Update: Added visualizations produced with the scripts
Update 2: Updated plot.probabilities() so that response variable can have arbitrary levels

I have been using ordinal package to crunch data. Tutorial for mixed models is only for clmm2 and not for clmm. Here are code for visualizing predicted probabilities for clmm. All the code is  based on  clmm2 tutorial.

A function to convert the effect to probability:

```pred <- function(eta, theta, cat = 1:(length(theta) + 1), inv.link = plogis) {
Theta <- c(-1000, theta, 1000)

Then actual drawing function (based on the clmm2 tutorial).

EDIT: fixed the code to check the ordinal scale of coefficients instead of assuming 5-point scale.

```plot.probabilities<-function(grid, model, leg, draw.these=NULL, main="", xlab="", legend.pos="topleft", ylim=NULL, col=NULL, lty=NULL) {
co <- model\$coefficients[1:length(model\$y.levels)-1]
pre.mat <- pred(eta=rowSums(grid), theta=co)
n.total<-length(pre.mat)
n.rows<-length(pre.mat[1,])
n <- n.total/n.rows
ylim <- if(is.null(ylim)) c(0,1)
else ylim
draw.these <- if(is.null(draw.these)) 1:n
else draw.these

plot(model\$y.levels, pre.mat[draw.these,], lty=1, type="l", ylim=ylim, xlab=xlab, axes=FALSE, ylab="Probability", las=1, main=main)
axis(1)
axis(2)
i <- 1
for(k in draw.these) {
draw_color <- if(is.null(col)) "black"
else col[i]
curr_lty <- if(is.null(lty)) i
else lty[i]
lines(model\$y.levels, pre.mat[k,], lty=curr_lty, col=draw_color)
i <- i + 1
}
if(is.null(lty)) {
legend(legend.pos, leg, lty=1:i, bty="n")
}
else {
legend(legend.pos, leg, lty=lty, bty="n")
}
}```

To draw fixed effects and random effects (this will reproduce figure 3 from clmm2 tutorial):

```library("ordinal")
data(wine)
fm2 <- clmm(rating ~ temp + contact + (1|judge), data = wine,Hess = TRUE, nAGQ = 10)
mat_no_cold <- expand.grid(judge = qnorm(0.95) * c(-1, 0, 1) * fm2\$stDev, contact = c(0), temp=c(0))
mat_yes_cold <- expand.grid(judge = qnorm(0.95) * c(-1, 0, 1) * fm2\$stDev, contact = c(fm2\$beta), temp=c(0))
mat_no_warm <- expand.grid(judge = qnorm(0.95) * c(-1, 0, 1) * fm2\$stDev, contact = c(0), temp=c(fm2\$beta))
mat_yes_warm <- expand.grid(judge = qnorm(0.95) * c(-1, 0, 1) * fm2\$stDev, contact = c(fm2\$beta), temp=c(fm2\$beta))
par(mfrow = c(2,2))
plot.probabilities(mat_no_cold, fm2, c("5th %-tile judge","avg. judge","95th %-tile judge"), main="contact=no, temp=cold")
plot.probabilities(mat_yes_cold, fm2, c("5th %-tile judge%","avg. judge","95th %-tile judge"), main="contact=yes, temp=cold")
plot.probabilities(mat_no_warm, fm2, c("5th %-tile judge","avg. judge","95th %-tile judge"), main="contact=no, temp=warm")
plot.probabilities(mat_yes_warm, fm2, c("5th %-tile judge","avg. judge","95th %-tile judge"), main="contact=yes, temp=warm")```

And then drawing all fixed effect probabilities on a figure without random effects:

```mat <- expand.grid(contact = c(0, fm2\$beta), temp = c(0, fm2\$beta))
plot.probabilities(mat, fm2, c("contact=no,temp=cold","contact=yes,temp=cold","contact=no,temp=warm","contact=yes,temp=warm"), main="Rating probabilities for average judge")```

Petri Lankoski, D.Arts, is a Associate Professor in Game Studies at the school of Communication, Media and IT at the Södertörn University, Sweden. His research focuses on game design, game characters, role-playing, and playing experience. Petri has been concentrating on single-player video games but researched also (multi-player) pnp and live-action role-playing games. This blog focuses on his research on games and related things.

## 4 thoughts on “R / Ordinal Scripts”

1. lankoski says:

clmm fits cumulative link mixed models,
Tutz, G. and W. Hennevogl (1996). Random effects in ordinal regression models. Computational
Statistics & Data Analysis 22, 537–557.

2. Maggie Wisniewska says:

I am trying to use your code to graph clmm results for the following model
cb3 <- clmm(DIST ~ AGE + LOC + MOM + (1|ID), data = cb, Hess=TRUE, nAGQ=10)
I have 1/2, 2/3 treshold coefficients. I do not understand the following in you code:
library("ordinal")
data(wine)
My data, named cb, is already attached, so why I need to specify data()?
The error message that I get next to each plot.probability line is the following:
Error in xy.coords(x, y, xlabe l, ylabel, log) :
'x' and 'y' lengths differ

This is my code. Could you please let me know what is wrong with my code?
pred <- function(eta, theta, cat = 1:(length(theta) + 1), inv.link = plogis) {
Theta <- c(-1000, theta, 1000)

plot.probabilities<-function(grid, model, leg, draw.these=NULL, main="Calf's Distance from Mother", xlab="Distance from Mother", legend.pos="topleft",
ylim=NULL) {
co <- model\$coefficients[1:3]
pre.mat <- pred(eta=rowSums(grid), theta=co)
n.total<-length(pre.mat)
n.rows<-length(pre.mat[1,])
n <- n.total/n.rows
ylim <- if(is.null(ylim)) c(0,1)
else ylim
draw.these <- if(is.null(draw.these)) 1:n
else draw.these
plot(1:3, pre.mat[draw.these,], lty=1, type="l", ylim=ylim, xlab=xlab, axes=FALSE, ylab="Probability", las=1, main=main)
axis(1)
axis(2)
i <- 1
for(k in draw.these) {
lines(1:3, pre.mat[k,], lty=i)
i <- i + 1
}
legend(legend.pos, leg, lty=1:i, bty="n")
}

cb3 <- clmm(DIST ~ AGE + LOC + MOM +(1|ID), data = cb,Hess = TRUE, nAGQ = 10)
mat_M_MUS_UNW <- expand.grid(ID = qnorm(0.95) * c(-1, 0, 1) * cb3\$stDev, MOM=c(0),LOC=c(0),AGE=c(0))
mat_NM_MUS_UNW <- expand.grid(ID = qnorm(0.95) * c(-1, 0, 1) * cb3\$stDev, MOM= c(cb3\$beta), LOC=c(0),AGE=c(0))
mat_M_WAP_UNW <- expand.grid(ID = qnorm(0.95) * c(-1, 0, 1) * cb3\$stDev, MOM=c(0), LOC=c(cb3\$beta),AGE=c(0))
mat_NM_WAP_UNW <- expand.grid(ID = qnorm(0.95) * c(-1, 0, 1) * cb3\$stDev, MOM= c(cb3\$beta), LOC=c(cb3\$beta),AGE=c(0))
mat_M_WAP_WEA <- expand.grid(ID = qnorm(0.95) * c(-1, 0, 1) * cb3\$stDev, MOM=c(0), LOC=c(cb3\$beta),AGE=c(cb3\$beta))
mat_NM_WAP_WEA <- expand.grid(ID = qnorm(0.95) * c(-1, 0, 1) * cb3\$stDev, MOM= c(cb3\$beta), LOC=c(cb3\$beta),AGE=c(cb3\$beta))
mat_M_MUS_WEA <- expand.grid(ID = qnorm(0.95) * c(-1, 0, 1) * cb3\$stDev, MOM=c(0), LOC=c(0),AGE=c(cb3\$beta))
mat_NM_MUS_WEA <- expand.grid(ID = qnorm(0.95) * c(-1, 0, 1) * cb3\$stDev, MOM= c(cb3\$beta), LOC=c(0),AGE=c(cb3\$beta))
par(mfrow = c(4,4))
plot.probabilities(mat_M_MUS_UNW, cb3, c("5th %-tile ID","avg. ID","95th %-tile ID"), main="MOM=M, LOC=MUS, AGE=UNW")
plot.probabilities(mat_NM_MUS_UNW, cb3, c("5th %-tile ID","avg. ID","95th %-tile ID"), main="MOM=NM, LOC=MUS, AGE=UNW")
plot.probabilities(mat_M_WAP_UNW, cb3, c("5th %-tile ID","avg. ID","95th %-tile ID"), main="MOM=M, LOC=WAP, AGE=UNW")
plot.probabilities(mat_NM_WAP_UNW, cb3, c("5th %-tile ID","avg. ID","95th %-tile ID"), main="MOM=NM, LOC=WAP, AGE=UNW")
plot.probabilities(mat_M_WAP_WEA, cb3, c("5th %-tile ID","avg. ID","95th %-tile ID"), main="MOM=M, LOC=WAP, AGE=WEA")
plot.probabilities(mat_NM_WAP_WEA, cb3, c("5th %-tile ID","avg. ID","95th %-tile ID"), main="MOM=NM, LOC=WAP, AGE=WEA")
plot.probabilities(mat_M_MUS_WEA, cb3, c("5th %-tile ID","avg. ID","95th %-tile ID"), main="MOM=M, LOC=MUS, AGE=WEA")
plot.probabilities(mat_NM_MUS_WEA, cb3, c("5th %-tile ID","avg. ID","95th %-tile ID"), main="MOM=NM, LOC=MUS, AGE=WEA")

mat <- expand.grid(AGE = c(0, cb3\$beta), MOT = c(0, cb3\$beta), LOC=c(0, cb3\$beta))
plot.probabilities(mat, cb3, c("AGE=WEA,LOC=WAP,MOM=NM","AGE=WEA,LOC=WAP,MOM=M","AGE=UNW,LOC=WAP,MOM=M","AGE=UNW,LOC=WAP,MOM=NM","AGE=WEA,LOC=MUS,MOM=M",
"AGE=WEA,LOC=MUS,MOM=NM","AGE=UNW,LOC=MUS,MOM=M","AGE=UNW,LOC=MUS,MOM=NM"),
main="Distance from Mother Probabilities for average Calf")

1. lankoski says:

The example code uses wine data from ordinal packacge. The data is loaded to environment using data(“wine”) command. If you are using your own data, that line is not needed.

It is hard to say what is the problem with your code just by reading it and I have never tested the code with something else than 5-point Likert scale reponce . If you sent me data or a subset of data (it can be simulated or some columns can be changed) so that I can take a look at what is happening there and hopefully help you with the issue.

1. lankoski says:

The issue is in the line of plot-probabilites:
co <- model\$coefficients[1:3]

This should be
co <- model\$coefficients[1:2]

My plot.probabilites() above has been updated so that it can determine the number of response coefficients automatically.

This site uses Akismet to reduce spam. Learn how your comment data is processed.