# Simulation in the profiling model

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

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

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 |

### Low level parameters

#### Products mean scores

#### Panelist mean scores

#### Panelist scale use (multiplicative factor)

#### Panelist standard error

## Discussion

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