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")
|
|
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")
|
|
Petri Lankoski 
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.
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.