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)
     sapply(cat, function(j) inv.link(Theta[j + 1] - eta) - inv.link(Theta[j] - eta))}

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[1],], 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[2]), 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[1]))
mat_yes_warm <- expand.grid(judge = qnorm(0.95) * c(-1, 0, 1) * fm2$stDev, contact = c(fm2$beta[2]), temp=c(fm2$beta[1]))
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")

Rating probabilities for different conditions
Rating probabilities for different conditions

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

mat <- expand.grid(contact = c(0, fm2$beta[2]), temp = c(0, fm2$beta[1]))
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")

Rating probabilities for an average judge
Rating probabilities for an average judge
Advertisements

4 thoughts on “R / Ordinal Scripts

  1. 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. 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)
    sapply(cat, function(j) inv.link(Theta[j + 1] – eta) – inv.link(Theta[j] – eta))}

    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[1],], 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[3]), 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[2]),AGE=c(0))
    mat_NM_WAP_UNW <- expand.grid(ID = qnorm(0.95) * c(-1, 0, 1) * cb3$stDev, MOM= c(cb3$beta[3]), LOC=c(cb3$beta[2]),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[2]),AGE=c(cb3$beta[1]))
    mat_NM_WAP_WEA <- expand.grid(ID = qnorm(0.95) * c(-1, 0, 1) * cb3$stDev, MOM= c(cb3$beta[3]), LOC=c(cb3$beta[2]),AGE=c(cb3$beta[1]))
    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[1]))
    mat_NM_MUS_WEA <- expand.grid(ID = qnorm(0.95) * c(-1, 0, 1) * cb3$stDev, MOM= c(cb3$beta[3]), LOC=c(0),AGE=c(cb3$beta[1]))
    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[1]), MOT = c(0, cb3$beta[3]), LOC=c(0, cb3$beta[2]))
    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")

    • 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.

      • 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.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s