(This article was first published on

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.**Wiekvoet**, and kindly contributed to R-bloggers)### 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 his blog:**Wiekvoet**.R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...