# Multivariate linear mixed model using lmer

I want to detect differences between 4 treatments of a food product based on the results from a sensory evaluation, where 9 panelists assessed the 4 products (i.e. treatments) following a number of sensory attributes. The output variables are scores (on a scale from 0 to 100) reflecting the intensity of each attribute. For each attribute I want to know how the 4 products differ (or not) from each other. To achieve this, I use a multivariate linear mixed model using the lmer() function. Here’s how the data looks like:

panelist treatment repetition Odor1 Odor2 Odor3 Flavor1 Flavor2 Flavor3 Flavor4 Texture1
1     CTRL           1    35    25     0      64      83      21      20       84
1     T1             1    29    19    48      29      38      53      53       58
1     T2             1    28    22    38      36      44      51      51       39
1     T3             1    28    28    67      27      66      46      45       72
1     CTRL           2    19    32     0      42      43      22      21       42
1     T1             2     6    32    11      65      57      79      79       71
1     T2             2    46    32    78      38      63      63      63       63
1     T3             2    31    31    78      85      31      61      60       48
1     CTRL           3    64    51    50      50      67      67      67       76
1     T1             3    42    37    59      82      58      35      34       67
1     T2             3     0    26    58      57      70      21      21       54
1     T3             3     0    41    69      49      70      24      24       57


In order to perform the multivariate analysis with lmer(), I reshaped the dataset to a long-format table as explained here (and there) so that the scores are grouped in one column and indexed to the "attribute" variable (listing all sensory attributes):

panelist treatment repetition attribute value
1    CTRL          1        Odor1     35
1    T1            1        Odor1     29
1    T2            1        Odor1     28
1    T3            1        Odor1     28
1    CTRL          2        Odor1     19
1    T1            2        Odor1     6
1    T2            2        Odor1     46
1    T3            2        Odor1     31
1    CTRL          3        Odor1     64
...


"Panelist" is treated as a random factor to adjust for different uses of the scale among the assessors. I then used the following model as a starting point:

mod <- lmer(score ~ treatment + (1|panelist) + (1|attribute), data)


I feel that treating “attribute” as a random effect is the only way to go and the results obtained make sense with my data as opposed to treating it as fixed effect. But it only returns the significance between the 4 groups based on the 12 attributes all together. So my questions are:

• How can I figure out the significance (or not) among treatment groups for each attribute?

• Is this model the right choice in this case or should I use a generalized linear mixed model instead?

• Thanks for your answer. I am fairly new to this method but this sounds contradictory since (1|attribute) defines the random slope shared by scores from the same attribute while we define that the effect of attribute is constant across score. The results from lmer(score ~ treatment * attribute + (1|panelist) + (1|attribute), data) returns rather similar significance values for the fixed effect of treatment (from Satterthwaite's method) but not surprisingly all p-values for the fixed effect of attribute equals 1.0. – Drosof Aug 30 '18 at 22:43
• The truth is I am not sure how to handle the index variable attribute in this case. Including it as a fixed factor in the model as you suggested could be an easy way out as I would get the answer to my first question directly from the summary output. However, based on the results from lmer(score ~ treatment * attribute + (1|panelist), data) this model do not seem appropriate to my data due to highly significant results for attribute in which scores are very similar across treatment. – Drosof Aug 30 '18 at 22:58
• Should the fixed effect of attribute be nested into treatment? And how would I define this in the formula? – Drosof Sep 3 '18 at 20:58
I actually found the solution using the SensMixed package (cran.r-project.org/web/packages/SensMixed/SensMixed.pdf) with no need to reshape the dataset.