Creating Williams designs with even number of products

A Williams design is a special Latin square with the additional property of first order carry over (each product is followed equally often by each other product). In R the package crossdes can be used to create them.

> williams(4)
     [,1] [,2] [,3] [,4]
[1,]    1    2    4    3
[2,]    2    3    1    4
[3,]    3    4    2    1
[4,]    4    1    3    2

As a consequence of the carry over restriction, the design has the property that a row in this design is also reversed present in the design. Example, row 3 is row 1 reversed. For small designs this property can be used to generate the designs by brute force. Example; a four by four design has without loss of generality the first row and column designated 1 to 4 (using . as unknown).

1 2 3 4

2 . . .

3 . . .

4 . . .

Adding the reversal of the first row gives:

1 2 3 4

2 . . .

3 . . .

4 3 2 1

In practice, for an even number of products the last column can be created as reversed first column.

1 2 3 4

2 . . 3

3 . . 2

4 3 2 1

It is rather obvious that only one solution remains for the 2,2 location and the rest is simple filling in the blanks. The resulting design is the same as the solution of williams(4) with '3' and '4' permuted.

1 2 3 4

2 4 . 3

3 . . 2

4 3 2 1

It is possible to use the same approach for small designs with an even number of products. For an odd number of products the number of rows has to be doubled so this will follow in a later post. To create the design with 6 products a program it is more convenient than manual work. It appears, unknown to many, that using this approach there are two solutions for 6 products. These two solutions are not permutations of each other!

[[1]]

     [,1] [,2] [,3] [,4] [,5] [,6]

[1,]    1    2    3    4    5    6

[2,]    2    4    1    6    3    5

[3,]    3    1    5    2    6    4

[4,]    4    6    2    5    1    3

[5,]    5    3    6    1    4    2

[6,]    6    5    4    3    2    1

[[2]]

     [,1] [,2] [,3] [,4] [,5] [,6]

[1,]    1    2    3    4    5    6

[2,]    2    4    6    1    3    5

[3,]    3    6    2    5    1    4

[4,]    4    1    5    2    6    3

[5,]    5    3    1    6    4    2

[6,]    6    5    4    3    2    1

There are 8 solutions for 8 products and 192 solutions for 10 products. Above that, calculation time is prohibitive. Even to get to 10, the program had to be streamlined considerably.

R code

gendesign <- function(n=6) {

nr <- as.integer(n)

nc <- nr

desmat <- matrix(NA,nrow=nr,ncol=nc)

desmat[1,] <- 1L:nc

desmat[,1] <- 1L:nr

desmat[nr,] <- nc:1L

desmat[,nc] <- nr:1L

carover <- matrix(0L,nrow=nr,ncol=nc)

for (i in 1L:(nc-1L)) carover[i+1,i] <- carover[i,i+1] <- 1L

desobject <- list(desmat=desmat,carover = carover)

desresult <- list()

addpoint(desobject,desresult)

}

nextpos <- function(desmat) which(is.na(desmat),arr.ind=TRUE)

checkdes <- function(desobject,row,col) {

# test for only once carry over 

all(desobject$carover<=1) & 

# each product only once in each row and each column

!any(desobject$desmat[row,-col]== desobject$desmat[row,col],na.rm=TRUE )  &

!any(desobject$desmat[-row,col]== desobject$desmat[row,col],na.rm=TRUE )  

}

addpoint <- function(desobject,desresult) {

todo <- nextpos(desobject$desmat)

if (length(todo)==0) {

l <- length(desresult)

desresult[[l+1]] <- desobject$desmat

return(desresult)

} 

row <- todo[1L,1L]

col <- todo[1L,2L]

nc <- ncol(desobject$desmat)

dob <- desobject

for (i in 1L:nc) {

desobject$desmat[row,col] <- i

desobject$desmat[nc-row+1,nc-col+1] <- i

desobject$carover[desobject$desmat[row,col-1L],i] <- desobject$carover[desobject$desmat[row,col-1],i] + 1L

desobject$carover[i,desobject$desmat[row,col-1L]] <- desobject$carover[i,desobject$desmat[row,col-1L]] + 1L

other <- desobject$desmat[row,col+1L]

if (!is.na(other)) {

desobject$carover[other,i] <- desobject$carover[other,i] + 1

desobject$carover[i,other] <- desobject$carover[i,other] + 1

}

if (checkdes(desobject,row,col)) desresult <- addpoint(desobject,desresult)

desobject <- dob

}

desresult

}

gendesign(6)

Simulation in the profiling model

In this post I try to make a small simulation of the sensory (flavour) profiling data, and examine if the parameters of simulated data can be retrieved by the Bayesian model build in the previous posts.
The conclusion is that it is difficult, the amount of uncertainty is too large for parts of the data. Especially, the multiplicative effect is difficult to retrieve,

Setting up the data

The data setup is the same as the data input: I will take the same design as previously used which is the chocolate data of sensominer. I will use the same kind of parameters as obtained from the chocolate data.
data(chocolates)
design <- sensochoc[,1:4]
# hyper parameters
sdPanelmean <- 0.978
sdRound <- 0.383
sdSessionPanelist <- 0.191
mu.lsdP    <- 0.387
sigma.lsdP <- 0.323
sPanelistmean <- 1.06
sPanelistsd <- 0.22
Using equivalent formulas as in the data analysis, the responses are simulated. For explanations of the terms involved I refer to the posts created in May while building the model.  I won't repeat the model here either.
nPanelist <- nlevels(sensochoc$Panelist)
nSession <-  length(unique(sensochoc$Session))
nRound = length(unique(sensochoc$Rank))
mPanelist <- rnorm(nPanelist,sd=sdPanelmean)
mround <- rnorm(nRound,sd=sdRound)
mSessionPanelist <- rnorm(nPanelist*nSession,sd=sdSessionPanelist)
meanProduct <- c(6.98,6.60,4.84,6.36,6.66,6.48)
sPanelist <- rnorm(nPanelist,mean=1,sd=0.22)
grandmean <- mean(meanProduct)
mProduct <- meanProduct-grandmean
sdPanelist <- rlnorm(nPanelist,mean=mu.lsdP,sd=sigma.lsdP)

design$score <- grandmean + mPanelist[design$Panelist] +  
 mSessionPanelist[interaction(design$Session,design$Panelist)] + 
 sPanelist[design$Panelist]*
                  (mProduct[design$Product] + mround[design$Rank])
design$score <- design$score + rnorm(nrow(design),sd=sdPanelist[design$Panelist]) 


jagsfit <- jags(data=data_list,inits=inits,model.file=model.file,parameters.to.save=parameters,n.chains=4,DIC=FALSE,n.iter=20000)

Results  

Hyperparameters

For the hyperparameters, the table shows the input data plus the resulting mean and sd of the sampler. 

	input	estimated mean	estimated sd
sdPanelmean	0.978	0.883	0.206
sdRound	0.383	0.262	0.142
sdSessionPanelist	0.191	0.198	0.138
mu.lsdP	0.387	0.332	0.08
sigma.lsdP	0.323	0.358	0.076

It would seem from these data, that the model does reasonable, except for the two parameters involved in lsdP. To remind, these parameters are hyperparameters for the panelists individual accuracy. 

Low level parameters

The low level parameters can be extracted and plotted with a structure such as this one:

ProductMean <- fitsummary$quantiles[ 

   grep('meanProduct',rownames(fitsummary$quantiles)),]

colnames(ProductMean) <-c('x2.5','x25','x50','x75','x97.5')

ProductMean <- as.data.frame(ProductMean)

ProductMean <- cbind(meanProduct,ProductMean)

limits <- aes(ymax = ProductMean$x97.5, ymin=ProductMean$x2.5) 

p <- ggplot(ProductMean, aes(y=x50, x=meanProduct)) 

p + geom_point() + geom_errorbar(limits, width=0.2) + 

scale_y_continuous('Product mean') + 

                scale_x_continuous("Input mean") +

geom_abline(slope=1)

I feel sure people will forgive me for not repeating this scrips 4 times with different parameters.

Products mean scores

It would seem that the product means are retrieved reasonable well. 

Panelist mean scores

Panelist mean scores are retrieved reasonable well. It appears the panelists mean is estimated so accurately, that some of the parameters are different from 0. The same happened in the chocolate data, which was used to provide hyper parameters for the simulation, which seems a good thing. 

Panelist scale use (multiplicative factor)

The panelist scale is a difficult thing to estimate. On top of that, these parameters seem closer to each other than the parameters obtained from the chocolate data. The parameters also have a large error. When one thinks about this, it is not strange. Look at the product means. One product is quite different form the remainder. This means the score for this product is mainly determining a panelists scale. In regression terms, it is a leverage point. Even though every product is scored twice by each panelist, this is not good news for determining panelists scale usage. Maybe this is only possible when the products have quite different scores. But, in that case any model will do the trick of finding product differences. In case of small differences between products the model may not be performing better than a more simple model after all. 

Panelist standard error

The panelists standard error can be retrieved reasonably well. The persons with the highest values obtained do indeed have high values. It could be better, but the persons who give the highest values are indeed the worst performers. It would seem that having this parameter is useful in the model, even though the hyper parameters are difficult to estimate.

Discussion

The model performs reasonably well. However, a case can be made of removing the panelists scale use parameters. The information does not seem to be there to estimate this for the relevant data. A parameter of panelists standard error is valuable. This is a bit disappointing. I wanted a structural component, which says this panelist has the systematic behaviour, but that seems to be too difficult. What I get is a model which says a bit about quality of panelists without ability to split this into its components. 

A complete Bayesian model for sensory profiling data

In this post I will try to add an important parts in the sensory profiling model I have been building. This concerns the question: 'Are all panelists equally reproducible?'. Obviously the answer is no, some are better than others. From this observation stems the approach in which under performing panelists are removed prior to the final analysis. The beauty of the Bayesian model is that this removal is not needed. The model can weigh down the panelists based on performance.
Finally, it could be chosen to introduce the concept that the data is not normal. An error distribution with fatter tails may give more robust results. For the current data, any look at the data will reveal that it is not normal. Panelists only use the numbers 0 to 10, suggesting a ordered logistic model, which is for another time. A continuous line scale with values 0 to 100 is used most often in sensory profiling. For a line scale using a t distribution rather than normal distribution can be appropriate. The beauty is, in the Bayesian model the data can be used to determine the degrees of freedom for the t distribution. To allow general application of the model I will leave the error normal though.

Modelling panelists' individual error

The model statements needed for this model part are taken from Data analysis using regression and multilevel/hierarchical models' by Gelman and Hill (2007).
First the error must be made panelist dependent:

y[i] ~ dnorm(fit[i],tau)

becomes

y[i] ~ dnorm(fit[i],tauPanelist[Panelist[i]])

Obviously tauPanelist needs its prior distribution, for which we need non informative hyperpriors. 

for (i in 1:nPanelist) {

         tauPanelist[i] <-  pow(sdPanelist[i],-2)

sdPanelist[i] ~ dlnorm(mu.lsdP,tau.lsdP)

        }

mu.lsdP ~ dnorm(0,.0001)

tau.lsdP <- pow(sigma.lsdP,-2)

sigma.lsdP ~ dunif(0,100)

Our interest is in sdPanelist[ ], which will be examined in the result phase.
A consequence of adding this model part is the need for more precise variable names. There was a variable sdPanelist before, which concerned the standard deviation between panelists' means. This has been renamed to sdPanel.

Result

Results consist of three parts, statements about the model, the panel and statements about the product. 

Model results

This section should contain the output of the fit. However, by now it is 128 lines of output and a stack of figures. Who wants that in a blog? For briefness a few remarks. I was very afraid that the model would misbehave and give correlation between panelists error and panelists scale usage values. This did not happen in any obvious way, for which I am very happy. 

What did happen was a lack of mixing for the term sdSessionPanelist. If this happens more often action is needed. As it is, this is a parameter of minor interest, so I let it be. However, the temptation is there to either attempt to reformulate the model or increase the number of samples.

Panel results

To start with the panel, we can draw three figures, location, scale and error.

Panelists' means reveal that panelist 10 scores a bit higher than the other panelists. Panelist 25 also has a bit high value and 21 a bit low value. 

Panelists' scales reveal that panelists 13 and 23 use a bit wider scale than the other panelists.  

Finally we have the figure for panelists' error. Panelist 29, 7 and 20 are less well performing than the other panelists. 

All things considered, I would focus on the most important parts of the panel performance. The few panelists with high standard error, such as panelist 29. After that I would look at location. Scale usage seems to be much more in order.

However, there are 14 descriptors in the data, so similar plots can be made for the other descriptors. So, based on those results I would reassess my statements about training. Hence we need all the analysis. While fitting all this in a loop is not much of a problem, this does not seem satisfactory. Panels run as a matter of standard on daily or weekly basis. How about a report for each product set, in full, within a day? This way panel feedback can be quick and efficient. To get there, it would be nice to make the analysis in a standard report. This seems to beg for usage of sweave, odfweave or knitr. Especially knitr seems to be the product to use, considering the recent postings on www.r-bloggers.com.

Product results

The added model terms have increased the number of product differences and decreased standard deviations. This is not so strange, the less performing panelists are weighed down. In this respect, it seems the model is performing admirably.
   v1 v2

2   1  3
3   2  3
4   1  4
6   3  4
9   3  5
11  1  6
13  3  6



          Mean        SD    Naive SE Time-series SE
choc1 6.987792 0.1833657 0.002899266    0.002886758
choc2 6.597215 0.1796334 0.002840253    0.003067686
choc3 4.842514 0.1975550 0.003123619    0.003345383
choc4 6.363283 0.1787572 0.002826399    0.002985240
choc5 6.665101 0.1833072 0.002898341    0.003356131
choc6 6.486780 0.1795052 0.002838226    0.002685428

R code

library(SensoMineR)

library(coda)

library(R2jags)

FullContrast <- function(n) {

UseMethod('FullContrast',n)

}

FullContrast.default <- function(n) stop('FullContrast only takes integer and factors')

FullContrast.integer <- function(n) {

mygrid <- expand.grid(v1=1:n,v2=1:n)

mygrid <- mygrid[mygrid$v1<mygrid$v2,]

rownames(mygrid) <- 1:nrow(mygrid)

as.matrix(mygrid)

}

FullContrast.factor <- function(n) {

FullContrast(nlevels(n))

}

data(chocolates)

data_list <- with(sensochoc,list(Panelist = as.numeric(Panelist),

nPanelist = nlevels(Panelist),

Product = as.numeric(Product),

nProduct = nlevels(Product),

Round = Rank,

nRound = length(unique(Rank)),

Session = Session,

nSession = length(unique(Session)),

y=CocoaA,

N=nrow(sensochoc)))

data_list$Productcontr <-  FullContrast(sensochoc$Product) 

data_list$nProductcontr <- nrow(data_list$Productcontr)

model.file <- file.path(tempdir(),'mijnmodel.txt')

mymodel <- function() {

# core of the model  

for (i in 1:N) {

fit[i] <-  grandmean +  mPanelist[Panelist[i]] +  mSessionPanelist[Session[i],Panelist[i]] + 

(mProduct[Product[i]] + mRound[Round[i]]) * sPanelist[Panelist[i]]

# 

y[i] ~ dnorm(fit[i],tauPanelist[Panelist[i]])

}

 

for (i in 1:nPanelist) {

    tauPanelist[i] <-  pow(sdPanelist[i],-2)

sdPanelist[i] ~ dlnorm(mu.lsdP,tau.lsdP)

    }

mu.lsdP ~ dnorm(0,.0001)

tau.lsdP <- pow(sigma.lsdP,-2)

sigma.lsdP ~ dunif(0,100)

grandmean ~ dnorm(0,.001)

# variable Panelist distribution  

for (i in 1:nPanelist) {

mPanelist[i] ~ dnorm(0,tauPanelMean) 

}

tauPanelMean ~ dgamma(0.001,0.001)

sdPanelMean <- sqrt(1/tauPanelMean)

# Product distribution 

for (i in 1:nProduct) {

mProduct[i] ~ dnorm(0,tauProduct)

}

tauProduct ~ dgamma(0.001,0.001)

sdProduct <- sqrt( 1/tauProduct)

for (i in 1:nRound) {

mRound[i] ~ dnorm(0,tauRound) 

}

tauRound <- 1/(sdRound^2)

sdRound ~ dnorm(0,4) %_% T(0,)

for (i in 1:nSession) {

for (j in 1:nPanelist) {

  mSessionPanelist[i,j] ~ dnorm(0,tauSessionPanelist)

   }

}

tauSessionPanelist <- 1/(sdSessionPanelist^2)

sdSessionPanelist ~ dnorm(0,4) %_% T(0,)

#

# distribution of the multiplicative effect

for (i in 1:nPanelist) {

sPanelist[i] <- exp(EsPanelist[i])

EsPanelist[i] ~ dnorm(0,9)  

}

# getting the interesting data

# true means for Panelist

for (i in 1:nPanelist) {

meanPanelist[i] <-  grandmean + mPanelist[i] +  

mean(mSessionPanelist[1:nSession,i]) + 

      ( mean(mProduct[1:nProduct]) + mean(mRound[1:nRound]))*sPanelist[i]

}

# true means for Product

for (i in 1:nProduct) {

meanProduct[i] <- grandmean + 

mean(mPanelist[1:nPanelist]) + 

mean(mSessionPanelist[1:nSession,1:nPanelist]) +

(mProduct[i] + mean(mRound[1:nRound]) )*exp(mean(EsPanelist[1:nPanelist]))  

}

for (i in 1:nProductcontr) {

Productdiff[i] <- meanProduct[Productcontr[i,1]]-meanProduct[Productcontr[i,2]]

}

}

write.model(mymodel,con=model.file)

inits <- function() list(

grandmean = rnorm(1,3,1),

mPanelist = c(0,rnorm(data_list$nPanelist-1)) ,

mProduct = c(0,rnorm(data_list$nProduct-1)) ,

EsPanelist = rnorm(data_list$nPanelist) ,

tau = runif(1,1,2),

tauPanel = runif(1,1,3),

tauProduct = runif(1,1,3),

mRound = rnorm(data_list$nRound),

sdRound = abs(rnorm(1)),

mSessionPanelist= matrix(rnorm(data_list$nPanelist*data_list$nSession),nrow=data_list$nSession),

sdSessionPanelist = abs(rnorm(1)),

sdPanelist=runif(data_list$nPanelist,.5,1.5),

mu.lsdP=.1,

sigma.lsdP=.1

)

#parameters <- c('sdPanelist','sdProduct','gsd','meanPanelist','meanProduct','Productdiff','sPanelist','sdRound','mRound')

parameters <- c('sdPanelMean','sdProduct','meanProduct','Productdiff','mPanelist','sPanelist','sdRound','mRound','sdSessionPanelist','sdPanelist','mu.lsdP','sigma.lsdP')

jagsfit <- jags(data=data_list,inits=inits,model.file=model.file,parameters.to.save=parameters,n.chains=4,DIC=FALSE,n.iter=20000)

# jagsfit # not shown

# plot(jagsfit) # not shown

jagsfit.mc <- as.mcmc(jagsfit)

# plot(jagsfit.mc) # not shown

fitsummary <- summary(jagsfit.mc)

# extract the scale effects and plot them

sPanelist <- fitsummary$quantiles[ grep('sPanelist',rownames(fitsummary$quantiles)),]

colnames(sPanelist) <-c('x2.5','x25','x50','x75','x97.5')

sPanelist <- as.data.frame(sPanelist)

sPanelist$pnum <- 1:nrow(sPanelist)

library(ggplot2)

limits <- aes(ymax = sPanelist$x97.5, ymin=sPanelist$x2.5) 

p <- ggplot(sPanelist, aes(y=x50, x=pnum)) 

p + geom_point() + geom_errorbar(limits, width=0.2) + scale_y_log10('Panelist scale') + scale_x_continuous("Panelist number") 

sdPanelist <- fitsummary$quantiles[ grep('sdPanelist',rownames(fitsummary$quantiles)),]

colnames(sdPanelist) <-c('x2.5','x25','x50','x75','x97.5')

sdPanelist <- as.data.frame(sdPanelist)

sdPanelist$pnum <- 1:nrow(sdPanelist)

limits <- aes(ymax = sdPanelist$x97.5, ymin=sdPanelist$x2.5) 

p <- ggplot(sdPanelist, aes(y=x50, x=pnum)) 

p + geom_point() + geom_errorbar(limits, width=0.2) + scale_y_log10('Panelist standard error') + scale_x_continuous("Panelist number") 

mPanelist <- fitsummary$quantiles[ grep('mPanelist',rownames(fitsummary$quantiles)),]

colnames(mPanelist) <-c('x2.5','x25','x50','x75','x97.5')

mPanelist <- as.data.frame(mPanelist)

mPanelist$pnum <- 1:nrow(mPanelist)

limits <- aes(ymax = mPanelist$x97.5, ymin=mPanelist$x2.5) 

p <- ggplot(mPanelist, aes(y=x50, x=pnum)) 

png('meanpanelist.png')

p + geom_point() + geom_errorbar(limits, width=0.2) + scale_y_continuous('Panelist mean') + scale_x_continuous("Panelist number") 

dev.off()

# extract product differences

Productdiff <- fitsummary$quantiles[ grep('Productdiff',rownames(fitsummary$quantiles)),]

# extract differences different from 0

data_list$Productcontr[Productdiff[,1]>0 | Productdiff[,5]<0,]

# get the product means

ProductMean <- fitsummary$statistics[ grep('meanProduct',rownames(fitsummary$quantiles)),]

rownames(ProductMean) <- levels(sensochoc$Product)

ProductMean

Extending the sensory profiling data model

In this post I extend the multiplicative Bayesian sensory profiling model with effects for rounds and sessions. Is is not a difficult extension, but it brings the need for informative priors into the model. I do believe round and session effects exist, but, they are small. The Bayesian paradigm allows to employ small directly in the model.

Rounds

This model envisions round effects as a way to model effects due to carry over and contrast. To explain these, envision a panelist tasting. He/she rinses to clean the mouth, a product is served and subsequently tasted. After scoring all the descriptors for the product a few minutes break is planned. In this break there is time for rinsing, maybe crackers or apples or another convenient product to clean the palate. After these few minutes the next product is served.

In the ideal situation, any residuals from the first product are removed when the second product is tasted. In the less ideal reality, this is not always so. Hence a product is influenced by the previous product, this is carry over. It is also possible, if the first product is very strong in taste and the second weak, that the panelist may feel to enhance the difference due to the contrast with the first product. Again, a property from the first product is influencing the second product. This is a way to get carry over effect. There are some variations there, but the core of the message is this. A product which is evaluated as second product, is evaluated differently than the first product. There are protocols such as rinsing to minimize this difference, there are experimental designs (Williams designs) to ensure these effects do not result is bias in product effects. In the end, however, they may still be there. 

In the current model I have chosen to keep the model simple, and make the round effect a simple round effect, independent on products tasted and panelists which taste. This is obviously a simplification, but the model is complex enough as it is.

Implementation of a round effect

Given the way I explained the round effect, it follows that I see round effect as a perception effect. The consequence is that the round effect is subjected to similar scaling effects as the product effect. This leads to the following effect

  fit[i] <- ...+ (mProduct[Product[i]] + mRound[Round[i]])*sPanelist[Panelist[i]]
The second part of the round effect is the hyperdistribution for mRound. Rather than making this uninformative, I was this to be informative, in order to suggest these effects are small. mRound is normal distributed with precision tauRound
  mRound[i] ~ dnorm(0,tauRound) 
The prior is on the distribution of tauRound. The associated standard deviation is from a half normal distribution with standard deviation 1/2 (precision 4, sd = 1/sqrt(precision)). Hence sdRound is larger than 0 and probably less than 1 (two standard deviations).
  tauRound <- 1/(sdRound^2)
  sdRound ~ dnorm(0,4) %_% T(0,)
The construction %_% T(0,) needs explanation. Basically T(0,) adapts the distribution to be larger than 0 lower than infinite. JAGS is happily using the format dnorm(0,4) T(0,). However since the model is written in R and stored as R model, this is an illegal statement: an operator is needed. Hence the dummy instruction %_%. R is happy and thinks %_% will be resolved at runtime. Subsequently write.model finds this operator and removes it before it writes the model to the temp file, so JAGS does not know it was there at some point. Everybody happy.

Interpretation of the round effect

In terms of interpretation, there are two aspects. First of all, what is the size of the specific rounds. If the first round has a clear effect, this may mean that an improved palate cleanser can improve the data. The second aspect is the size of the round effect as such. If the effect is large, then it is possible that the tasting protocol needs to be adapted. Both of these are a matter of judgement rather than science. 

Session effect

The session effect is a slightly different kind of beast as the round effect. In pure ANOVA terms, it states that scores of one day are different, lower or higher, than another day. Unless one thinks that a product changes between days, this is actually an improbable effect. It is much more probable that for each panelist one day is different from the second day. In terms of mixed models, panelist * day is a random effect. Given the way this is introduced, the effect is additive on top of the panelist effect. Hence it is not dependent on panelists scale usage. 
fit[i] <-  ... +  mSessionPanelist[Session[i],Panelist[i]] + ...
The prior for mSessionPanelist[ ]  is the same as for mRound[ ].
Interpretation of the SessionPanelist effect is more simple than for rounds. There is not really much point in looking at each individual result, unless there is reason to suspect high variation for a specific panelist. Hence only the standard deviation is extracted from JAGS. A big standard deviation means that the panel is unsure and needs more training on the associated descriptor.

Results

Jags output shows a low value for n.eff and Rhat a bit larger than 1 for the variable sdSessionPanelist. It appears this parameter has some difficulty mixing and updating. This is also the reason why the number of iterations has been increased to 20000.
It seems round is a bit more important than PanelistSession. The round effect is showing a large value for the first round, hence perhaps a palate cleanser might be in order. The second round has a similar negative effect, perhaps a contrast effect to the big effect of the first round.


Inference for Bugs model at "/tmp/RtmpJUFGM7/mijnmodel.txt", fit using jags,
 4 chains, each with 20000 iterations (first 10000 discarded), n.thin = 10
 n.sims = 4000 iterations saved
                  mu.vect sd.vect   2.5%    25%    50%    75%  97.5%  Rhat n.eff
Productdiff[1]      0.502   0.274 -0.022  0.318  0.501  0.684  1.052 1.001  4000
Productdiff[2]      2.109   0.288  1.542  1.924  2.110  2.294  2.681 1.001  4000
Productdiff[3]      1.607   0.282  1.059  1.419  1.610  1.796  2.159 1.001  4000
Productdiff[4]      0.657   0.272  0.132  0.478  0.651  0.842  1.202 1.001  4000
Productdiff[5]      0.155   0.268 -0.373 -0.026  0.153  0.337  0.674 1.001  4000
Productdiff[6]     -1.452   0.273 -1.989 -1.639 -1.454 -1.268 -0.920 1.001  4000
Productdiff[7]      0.281   0.271 -0.252  0.100  0.279  0.465  0.806 1.001  4000
Productdiff[8]     -0.222   0.270 -0.751 -0.400 -0.221 -0.037  0.313 1.001  4000
Productdiff[9]     -1.829   0.283 -2.391 -2.018 -1.826 -1.643 -1.268 1.001  4000
Productdiff[10]    -0.377   0.265 -0.897 -0.553 -0.381 -0.200  0.145 1.001  4000
Productdiff[11]     0.541   0.270  0.032  0.354  0.535  0.725  1.078 1.001  4000
Productdiff[12]     0.038   0.270 -0.482 -0.141  0.038  0.221  0.562 1.001  4000
Productdiff[13]    -1.569   0.276 -2.114 -1.753 -1.573 -1.383 -1.021 1.001  4000
Productdiff[14]    -0.117   0.266 -0.638 -0.300 -0.114  0.060  0.402 1.001  4000
Productdiff[15]     0.260   0.264 -0.255  0.080  0.257  0.436  0.795 1.001  4000
gsd                 1.611   0.068  1.482  1.564  1.608  1.655  1.747 1.001  4000
mRound[1]           0.424   0.258 -0.063  0.251  0.409  0.585  0.956 1.001  4000
mRound[2]          -0.417   0.264 -0.963 -0.582 -0.401 -0.241  0.060 1.001  3300
mRound[3]           0.272   0.250 -0.198  0.112  0.261  0.425  0.807 1.001  3300
mRound[4]          -0.191   0.255 -0.737 -0.346 -0.179 -0.024  0.296 1.001  4000
mRound[5]          -0.275   0.245 -0.796 -0.430 -0.262 -0.111  0.178 1.001  4000
mRound[6]           0.090   0.244 -0.406 -0.058  0.089  0.240  0.574 1.001  4000
meanProduct[1]      6.973   0.202  6.587  6.837  6.971  7.114  7.360 1.001  4000
meanProduct[2]      6.471   0.196  6.089  6.343  6.470  6.603  6.851 1.001  3500
meanProduct[3]      4.864   0.207  4.463  4.725  4.864  5.000  5.275 1.001  4000
meanProduct[4]      6.316   0.193  5.941  6.188  6.311  6.445  6.702 1.001  4000
meanProduct[5]      6.693   0.195  6.312  6.564  6.694  6.822  7.080 1.001  4000
meanProduct[6]      6.433   0.192  6.039  6.304  6.438  6.564  6.800 1.001  4000
sdPanelist          0.959   0.176  0.646  0.837  0.943  1.064  1.350 1.001  4000
sdProduct           0.896   0.393  0.427  0.638  0.804  1.049  1.903 1.001  4000
sdRound             0.425   0.170  0.167  0.307  0.401  0.512  0.830 1.002  3000
sdSessionPanelist   0.237   0.147  0.020  0.119  0.217  0.336  0.560 1.024   300


For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

Product differences are same as before. However, when the the round effect had a smaller precision and hence larger standard deviation, it was noted that product difference between 1 and 6 was not observed. It is a bit disturbing the effect is so clearly dependent on the prior employed.

Product means are slightly closer to each other than with the previous model, while the standard deviations are almost the same. 

          Mean        SD    Naive SE Time-series SE

choc1 6.973338 0.2018855 0.003192090    0.003270188

choc2 6.470882 0.1958250 0.003096265    0.003057404

choc3 4.863912 0.2069963 0.003272899    0.003379317

choc4 6.316028 0.1925660 0.003044736    0.003290025

choc5 6.692635 0.1949687 0.003082725    0.003631917

choc6 6.432695 0.1921969 0.003038899    0.003378925

Discussion

The beauty of building models in JAGS is revealed here. A few effects are added, but the structure is not overly influenced. 

However, also the problems are revealed. The prior in round effects has influence on our opinion of product differences. There are two aspects to this. One of these is the undesired effect of trying to do hypothesis testing. A small change in data or model can push a parameter from one side of the brink to the other. Reality is that there is some indication an effect exists and the black and white line should not be drawn. The second aspect relates to model building itself. On the one hand it is undesirable that a change in prior has influence on the interpretation. On the other hand, this is what we do anyway. If a classical ANOVA is used, adding or removing a factor has similar effects on significance's. Getting it from a prior is a novelty, but should not make much of a difference in looking at the data analysis.

A final discussion point is how to combine the round and session effects with the scale effect previously used.   The current approach is to make the choice based on philosophical base rather than on data. While this is not ideal, I have the feeling this is not much of a problem. The round and session effects are small. Given the amount of data, model complexity and size of effects, I don't even think it is possible to discriminate between model variations data based. So, I do not think this is a large problem.

R code

library(SensoMineR)
library(coda)
library(R2jags)


FullContrast <- function(n) {
UseMethod('FullContrast',n)
}
FullContrast.default <- function(n) stop('FullContrast only takes integer and factors')
FullContrast.integer <- function(n) {
mygrid <- expand.grid(v1=1:n,v2=1:n)
mygrid <- mygrid[mygrid$v1<mygrid$v2,]
rownames(mygrid) <- 1:nrow(mygrid)
as.matrix(mygrid)
}
FullContrast.factor <- function(n) {
FullContrast(nlevels(n))
}


data(chocolates)


data_list <- with(sensochoc,list(Panelist = as.numeric(Panelist),
nPanelist = nlevels(Panelist),
Product = as.numeric(Product),
nProduct = nlevels(Product),
Round = Rank,
nRound = length(unique(Rank)),
Session = Session,
nSession = length(unique(Session)),
y=CocoaA,
N=nrow(sensochoc)))
data_list$Productcontr <-  FullContrast(sensochoc$Product) 
data_list$nProductcontr <- nrow(data_list$Productcontr)




model.file <- file.path(tempdir(),'mijnmodel.txt')


mymodel <- function() {
# core of the model  
for (i in 1:N) {
fit[i] <-  grandmean +  mPanelist[Panelist[i]] +  mSessionPanelist[Session[i],Panelist[i]] + 
(mProduct[Product[i]] + mRound[Round[i]]) * sPanelist[Panelist[i]]
# 
y[i] ~ dnorm(fit[i],tau)
}
# grand mean and residual 
tau ~ dgamma(0.001,0.001)
gsd <-  sqrt(1/tau)
grandmean ~ dnorm(0,.001)
# variable Panelist distribution  
for (i in 1:nPanelist) {
mPanelist[i] ~ dnorm(0,tauPanelist) 
}
tauPanelist ~ dgamma(0.001,0.001)
sdPanelist <- sqrt(1/tauPanelist)
# Product distribution 
for (i in 1:nProduct) {
mProduct[i] ~ dnorm(0,tauProduct)
}
tauProduct ~ dgamma(0.001,0.001)
sdProduct <- sqrt( 1/tauProduct)
for (i in 1:nRound) {
mRound[i] ~ dnorm(0,tauRound) 
}
tauRound <- 1/(sdRound^2)
sdRound ~ dnorm(0,4) %_% T(0,)
for (i in 1:nSession) {
for (j in 1:nPanelist) {
  mSessionPanelist[i,j] ~ dnorm(0,tauSessionPanelist)
   }
}
tauSessionPanelist <- 1/(sdSessionPanelist^2)
sdSessionPanelist ~ dnorm(0,4) %_% T(0,)
#
# distribution of the multiplicative effect
for (i in 1:nPanelist) {
sPanelist[i] <- exp(EsPanelist[i])
EsPanelist[i] ~ dnorm(0,9)  
}
# getting the interesting data
# true means for Panelist
for (i in 1:nPanelist) {
meanPanelist[i] <-  grandmean + mPanelist[i] +  
mean(mSessionPanelist[1:nSession,i]) + 
      ( mean(mProduct[1:nProduct]) + mean(mRound[1:nRound]))*sPanelist[i]
}
# true means for Product
for (i in 1:nProduct) {
meanProduct[i] <- grandmean + 
mean(mPanelist[1:nPanelist]) + 
mean(mSessionPanelist[1:nSession,1:nPanelist]) +
(mProduct[i] + mean(mRound[1:nRound]) )*exp(mean(EsPanelist[1:nPanelist]))  
}
for (i in 1:nProductcontr) {
Productdiff[i] <- meanProduct[Productcontr[i,1]]-meanProduct[Productcontr[i,2]]
}

}


write.model(mymodel,con=model.file)


inits <- function() list(
grandmean = rnorm(1,3,1),
mPanelist = c(0,rnorm(data_list$nPanelist-1)) ,
mProduct = c(0,rnorm(data_list$nProduct-1)) ,
EsPanelist = rnorm(data_list$nPanelist) ,
tau = runif(1,1,2),
tauPanelist = runif(1,1,3),
tauProduct = runif(1,1,3),
mRound = rnorm(data_list$nRound),
sdRound = abs(rnorm(1)),
mSessionPanelist= matrix(rnorm(data_list$nPanelist*data_list$nSession),nrow=data_list$nSession),
sdSessionPanelist = abs(rnorm(1))
)


#parameters <- c('sdPanelist','sdProduct','gsd','meanPanelist','meanProduct','Productdiff','sPanelist','sdRound','mRound')
parameters <- c('sdPanelist','sdProduct','gsd','meanProduct','Productdiff','sdRound','mRound','sdSessionPanelist')


jagsfit <- jags(data=data_list,inits=inits,model.file=model.file,parameters.to.save=parameters,n.chains=4,DIC=FALSE,n.iter=20000)


jagsfit
#plot(jagsfit)


jagsfit.mc <- as.mcmc(jagsfit)
#plot(jagsfit.mc)
fitsummary <- summary(jagsfit.mc)


# extract differences
Productdiff <- fitsummary$quantiles[ grep('Productdiff',rownames(fitsummary$quantiles)),]
# extract differences different from 0
data_list$Productcontr[Productdiff[,1]>0 | Productdiff[,5]<0,]
# get the product means
ProductMean <- fitsummary$statistics[ grep('meanProduct',rownames(fitsummary$quantiles)),]
rownames(ProductMean) <- levels(sensochoc$Product)
ProductMean