I have tried to plot random effects after clmm
, but I got error message: "Error in sort.list(y):'x' must be atomic for `sort.list' Have you called 'sort' an a list?" The following codes are representative of my actual codes.
library(ordinal)
library(lattice)
###I am using the wine data in the ordinal package
d<-data.frame(wine)
result<-clmm(rating~1+temp+contact+(1+temp|judge), data=d)
###In my actual codes, I put "as.ordered(rating)" instead of "rating".
Then, I am trying to plot random effects of temp and judge:
dotplot(ranef(result, condVar=TRUE))
Then, the error message appears: "Error in sort.list(y):'x' must be atomic for `sort.list' Have you called 'sort' on a list?"
My best guess is that rating is implemented as ordered under clmm
, which seems for me to make sense given the error message. Yet, the guess is highly speculative, and I have no idea how to handle this situation. Specifically, what I would like to plot is the random effects of temp and judge (intercept) with their CIs. Please refer to the following plot that I generated using the codes below
result2<-lmer(as.numeric(rating)~1+temp+contact+(1+temp|judge), data=d)
dotplot(ranef(result2, condVar=TRUE))
If you could give any comment or suggestion on how to circumvent the seemingly-conflicting situation in which using clmm
prevents me from using dotplot
using the code above, it will be really appreciated.
I started out trying to hack the guts of lme4:::dotplot.ranef.mer
, but it turned out to be easier to do this with reshape2::melt
and ggplot2
. Sorry if you have a really strong preference for lattice
...
Note by the way that, technically, these aren't "confidence intervals" (since the predicted conditional modes are not "estimates" in the technical, frequentist sense) ...
Rearrange point estimates and SEs of conditional modes into a useful shape:
library(reshape2)
melt.ranef.clmm <- function(re,cv) {
reList <- lapply(re,
function(x) {
x$id <- reorder(factor(rownames(x)),x[,1])
return(melt(x,id.var="id"))
})
cvList <- lapply(cv,melt,id.var=NULL,value.name="se")
mm <- Map(cbind,reList,cvList)
return(mm)
}
Apply this to a model (m1
is a fitted model):
ss <- melt.ranef.clmm(ranef(m1),condVar(m1))
Plot:
library(ggplot2); theme_set(theme_bw())
ggplot(ss[[1]],aes(value,id))+
geom_errorbarh(aes(xmin=value-1.96*se,xmax=value+1.96*se),
height=0)+
ggtitle(names(ss)[1])+
geom_point(colour="blue")+
facet_wrap(~variable,scale="free_x")