Simulation shows gain of clmm over ANOVA is small
[This article was first published on Wiekvoet, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
After last post’s setting up for a simulation, it is now time to look how the models compare. To my disappointment with my simple simulations of assessors behavior the gain is minimal. Unfortunately, the simulation took much more time than I expected, so I will not expand it.Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Background
I have been looking at the ordered logistic model in a number of postings. The reason is that in general I find people use ANOVA for analyzing data on a nine point scale, whereas you would think an ordered logistic model works better. Two posts (first and second) showed the current methods in R, and and JAGS are easy to use and with some tweaking provide suitable output to present. Subsequently I made the simulated data generator. Now it is time to make the final comparison.
Simulations
The core of the simulator is explained elsewhere so I won’t explain here again. I did however notice a small error, so the corrected code is given here. Some new parts are added, wrappers around the data generator and the analysis. And to my big disappointment I could not even build that as desired. The call anova(Res.clmm2,Res.clmm) with subsequent extraction has been replaced by the ugly pchisq(2*(Res.clmm$logLik-Res.clmm2$logLik),2,lower.tail=FALSE ). 2 represents the degrees of freedom for products in my simulations. Somehow that call to ANOVA did not run within a function, after trying too many variations I choose the short cut.
library(ordinal)
library(ggplot2)
num2scale <- function(x,scale) {
cut(x,c(-Inf,scale,Inf),1:(length(scale)+1))
}
pop.limits2ind.limits <- function(scale,sd=.5) {
newscale <- scale+rnorm(length(scale),sd=sd)
newscale[order(newscale)]
}
obs.score <- function(obs,pop.score,pop.limits,
sensitivity.sd=.1, precision.sd=1,
additive.sd=2,center.scale=5,
labels=LETTERS[1:length(pop.score)]) {
# individual sensitivity (multiplicative)
obs.errorfreeintensity <- center.scale +
(pop.score-center.scale)*rlnorm(1,sd=sensitivity.sd)
#individual (additive)
obs.errorfreeintensity <- obs.errorfreeintensity +
rnorm(1,mean=0,sd=additive.sd)
# individual observation error
obs.intensity <- obs.errorfreeintensity+
rnorm(length(pop.score),sd=precision.sd)
# individual cut offs between categories
obs.limits <- pop.limits2ind.limits(pop.limits)
obs.score <- num2scale(obs.intensity,obs.limits)
data.frame(obs=obs,
score = obs.score,
product=labels
)
}
panel.score <- function(n=100,pop.score,pop.limits,
sensitivity.sd=.1,precision.sd=1,
additive.sd=2,center.scale=5,
labels=LETTERS[1:length(pop.score)]) {
la <- lapply(1:n,function(x) {
obs.score(x,pop.score,pop.limits,sensitivity.sd=sensitivity.sd,
precision.sd=precision.sd,additive.sd=additive.sd,labels=labels)
})
dc <- do.call(rbind,la)
dc$obs <- factor(dc$obs)
dc
}
overallP <- function(scores) {
Res.aov <- aov( numresponse ~ obs + product , data=scores)
paov <- anova(Res.aov)['product','Pr(>F)’]
Res.clmm <- clmm(score ~ product + (1|obs),data=scores)
Res.clmm2 <- clmm(score ~ 1 + (1|obs),data=scores)
c(ANOVA=paov,
clmm = pchisq(2*(Res.clmm$logLik-Res.clmm2$logLik),2,lower.tail=FALSE ) #.1687
)
}
onesim <- function(prodmeans,pop.limits,center.scale,additive.sd=2) {
scores <- panel.score(40,prodmeans,pop.limits,
center.scale=center.scale,additive.sd=additive.sd)
scores$numresponse <- as.numeric(levels(scores$score))[scores$score]
overallP(scores)
}
Simulation I, 5 categories
The first simulation is with 5 categories. This represents the Just About Right (JAR) scale. The final plot shows the difference between ANOVA and clmm is minimal.
pop.limits <- c(1,2.5,4.5,6)
nsim <- 250
sim5cat1 <- lapply(seq(0,.6,.05),function(dif) {
prodmeans=rep(3.5,3)+c(-dif,0,dif)
sa <- sapply(1:nsim,function(x) onesim(prodmeans,pop.limits))
data.frame(dif=dif,nreject=as.numeric(rowSums(sa<0.05)),
method=rownames(sa))
}
)
sim5cat1tb <- do.call(rbind,sim5cat1)
ggplot(sim5cat1tb, aes(dif,nreject/nsim ,colour=method)) +
geom_line() + xlab(‘Difference between products’) +
ylab(‘Proportion significant (at 5% Test)’) +
theme(legend.position = “bottom”) + ylim(c(0,1)) +
guides(colour = guide_legend(‘Analysis Method’))
pop.limits <- c(1,3:8,10)
prodmeans <- c(7,7,7)
scores <- panel.score(40,prodmeans,pop.limits,additive.sd=1)
scores$numresponse <- as.numeric(levels(scores$score))[scores$score]
nsim <- 250
sim9cat <- lapply(seq(0.,.6,.05),function(dif) {
prodmeans=c(7,7,7)+c(-dif,0,dif)
sa <- sapply(1:nsim,function(x) onesim(prodmeans,pop.limits,
additive.sd=1))
data.frame(dif=dif,nsign=as.numeric(rowSums(sa>0.05)),
method=rownames(sa))
}
)
sim9cattb <- do.call(rbind,sim9cat)
ggplot(sim9cattb, aes(dif,(nsim-nsign)/nsim ,colour=method)) +
geom_line() + xlab(‘Difference between products’) +
ylab(‘Proportion significant (at 5% Test)’) +
theme(legend.position = “bottom”) + ylim(c(0,1)) +
guides(colour = guide_legend(‘Analysis Method’))
pop.limits <- c(1,2.5,4.5,6)
nsim <- 250
sim5cat1 <- lapply(seq(0,.6,.05),function(dif) {
prodmeans=rep(3.5,3)+c(-dif,0,dif)
sa <- sapply(1:nsim,function(x) onesim(prodmeans,pop.limits))
data.frame(dif=dif,nreject=as.numeric(rowSums(sa<0.05)),
method=rownames(sa))
}
)
sim5cat1tb <- do.call(rbind,sim5cat1)
ggplot(sim5cat1tb, aes(dif,nreject/nsim ,colour=method)) +
geom_line() + xlab(‘Difference between products’) +
ylab(‘Proportion significant (at 5% Test)’) +
theme(legend.position = “bottom”) + ylim(c(0,1)) +
guides(colour = guide_legend(‘Analysis Method’))
Simulation 2, 9 categories
This simulation represents the intensity and liking scales. Again the difference between ANOVA and clmm are minimal.pop.limits <- c(1,3:8,10)
prodmeans <- c(7,7,7)
scores <- panel.score(40,prodmeans,pop.limits,additive.sd=1)
scores$numresponse <- as.numeric(levels(scores$score))[scores$score]
nsim <- 250
sim9cat <- lapply(seq(0.,.6,.05),function(dif) {
prodmeans=c(7,7,7)+c(-dif,0,dif)
sa <- sapply(1:nsim,function(x) onesim(prodmeans,pop.limits,
additive.sd=1))
data.frame(dif=dif,nsign=as.numeric(rowSums(sa>0.05)),
method=rownames(sa))
}
)
sim9cattb <- do.call(rbind,sim9cat)
ggplot(sim9cattb, aes(dif,(nsim-nsign)/nsim ,colour=method)) +
geom_line() + xlab(‘Difference between products’) +
ylab(‘Proportion significant (at 5% Test)’) +
theme(legend.position = “bottom”) + ylim(c(0,1)) +
guides(colour = guide_legend(‘Analysis Method’))
Conclusion
The differences between clmm and ANOVA seem to be minimal. I had not expected large differences, but especially at 5 categories my expectation were to find real differences as the continuous data is more violated. Obviously, a load more simulations would be needed to draw final conclusions. This is beyond the scope of my blog.
To conclude, in theory clmm is much more suitable than ANOVA for ordinal data. There are no reasons in terms of presentation to prefer ANOVA over clmm. But in practice the difference may be minimal.
To leave a comment for the author, please follow the link and comment on their blog: Wiekvoet.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.