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

### Background

### Simulations

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

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

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