Experiment: I collected data from N participants, each was shown 50 photos and asked to provide sharing likelihood (dependent variable). I also measured personal traits (e.g., affiliative score which where then grouped), and each photo had an associated valence score (also grouped). There were three experimental conditions. I want to measure the impact of condition, valence, and personal traits on the sharing likelihood.

I created a lmer model:

lmer1.model=lmer(share ~ gender+ age_group + overall_photo_share_freq_group+
     (1|photo), data = picshare_df, REML = TRUE)


Should the random part be participant_id instead of photo?

I want to use it for post-hoc tests. But the mean values I am getting from the model are very different than the mean values I directly calculate. In the model, the dependent variable is share and I want to test hypotheses for different levels of affiliative_score_group, which has three levels.


aggregate(picshare_df$share, by=list(picshare_df$affiliative_score_group),FUN=mean)

I get

Low -0.8699349          
Medium  -0.9134223          
High    -0.8120141  

But from the model:

 affiliative_score_group lsmean     SE  df asymp.LCL asymp.UCL
 Low                     -0.707 0.0568 Inf    -0.819       -0.596
 Medium                  -1.038 0.0565 Inf    -1.149    -0.927
 High                    -1.088 0.0563 Inf    -1.198    -0.977

These means are different than before, even the order is also different, previously Medium had the lowest value and High had the highest.

The test results using lsmeans

 $`differences from control of affiliative_score_group`
 contrast     estimate     SE  df z.ratio p.value
 Medium - Low   -0.331 0.0262 Inf -12.639 <.0001 
 High - Low     -0.380 0.0285 Inf -13.343 <.0001 

If I do the tests directly:

pairwise.t.test( g=picshare_df$affiliative_score_group , x=picshare_df$share ,p.adjust.method="bonferroni" ,pool.sd=TRUE )

       Low     Medium 
Medium 0.28793 -      
High   0.05501 0.00015

Now only High-Medium contrast is significant.

Both ways of testing shows results consistent with the respective means (from lsmeans and arithmetic average, respectively), but they give contradictory results for the same output. Also I am very confused about the mean values.


What I am doing wrong or not understanding?


With a simpler model, I get this:

lmer1.model=lmer(share ~ gender+ age_group + overall_photo_share_freq_group+
     (1|photo), data = picshare_df, REML = FALSE)

lsmeans(lmer1.model, list(trt.vs.ctrl1~affiliative_score_group))
$`lsmeans of affiliative_score_group`
 affiliative_score_group lsmean     SE  df asymp.LCL asymp.UCL
 Low                     -0.892 0.0550 Inf    -1.000    -0.785
 Medium                  -0.995 0.0554 Inf    -1.103    -0.886
 High                    -0.889 0.0547 Inf    -0.997    -0.782

 Results are averaged over the levels of: gender, age_group,     overall_photo_share_freq_group, valence_group, condition 
Degrees-of-freedom method: asymptotic 
Confidence level used: 0.95 

$`differences from control of affiliative_score_group`
 contrast     estimate     SE  df z.ratio p.value
 Medium - Low   -0.102 0.0250 Inf -4.093  0.0001 
 High - Low      0.003 0.0235 Inf  0.128  0.9817 

These are consistent with the observed means and results from pairwise.t.test. Is it valid to run 4 such simpler models than the complex one? I do not think so. Also, to me, intuitively the means from lsmean from the complex model make more sense (it has the order Lowaffiliative_score_group predictor) than the observed mean.

More questions

If I use the complex model, should I report the estimates from lsmeans when I talk about the average effect, or plot them (e.g., bar plot or interaction plot)?

  • 1
    $\begingroup$ What is your model? (i.e., can you edit your question to include the formula?) $\endgroup$
    – Ben Bolker
    Feb 7, 2020 at 15:37
  • $\begingroup$ It seems that the model presented to lsmeans was a mixed model, while your other calculations only represent the fixed effects. Please describe the mixed model, too. Providing more details about the nature of the share variable would also help determine whether treating it as a continuous variable is the best way to proceed. $\endgroup$
    – EdM
    Feb 7, 2020 at 16:08
  • $\begingroup$ @BenBolker, I have included the model formula. @EdM, share is an integer from a 7-point Likert item (this is not ideal but in CS often such data is used as continuous variable). $\endgroup$
    – Rakib
    Feb 7, 2020 at 16:17
  • $\begingroup$ The model you have fitted seems rather complicated. Have you tried comparing the sample means with the fitted models from a mixed model that only contains the affiliative_score_group variable? $\endgroup$ Feb 7, 2020 at 20:19
  • $\begingroup$ @DimitrisRizopoulos, I have updated my question $\endgroup$
    – Rakib
    Feb 7, 2020 at 20:41

1 Answer 1


For your second question:

What I am doing wrong or not understanding?

The mean values that you directly calculate don't take into account any of the other predictors included in your models, such as gender, age_group, condition, valence, or any group other than affiliative_score_group. They don't involve the further correction that you attempt with your inclusion of photo as a random effect.

So you shouldn't be surprised that your directly calculated marginal means disagree from the estimated marginal means provided by lsmeans. As a vignette for the emmeans package (a newer package that now provides the underlying code for lsmeans) puts it:

Estimated marginal means are based on a model – not directly on data.

The results from lsmeans thus do take into account all the predictors and interactions included in your model, in a way that your directly calculated means don't.

The answer to your first question:

Should the random part be participant_id instead of photo?

is to include participant_id as a random effect in addition to photo. The models that you showed don't have any way of accounting for the correlations of responses to the photos within each of the N participants. The standard p-value and confidence-interval calculations assume that all observations are independent, which isn't the case if there are multiple responses by the same individuals. Inclusion of participant_id as a random effect is a standard way to take such within-individual correlations into account, just as you tried to take within-photo correlations into account by including photo as a random effect.

Finally, I fear that your first model will be too complex unless you have an extremely large number of individuals tested. There is a very large number of predictors, as each interaction term adds a number of predictors equal to the product of the number of levels, minus 1, of each categorical variable. For example, with 3 levels of condition, 3 levels of valence_group, and 3 levels of affiliative_score_group the one term


would add 2*2*2=8 predictors to your model. My guess is, including interaction terms, you have an the order of 30 to 50 fixed-effect predictors in your model beyond the random effects. Unless you have on the order of 300-500 participants, and have a reasonable number of participants in each of the various groups, I would be seriously worried about overfitting. In that case, consider cutting back on the 3-way interaction terms in the model and only include interactions that you expect to be important based on your knowledge of the subject matter.

You certainly do not want to do separate models for each of the 4 groups of interest, as each of them would ignore information from the other groups and would lead to bias in estimates if the group memberships are at all correlated. Use all the groups in a single model and cut back complexity by cutting back on the interaction terms.

  • $\begingroup$ Thanks @EdM. Your answer is very illuminating, but I still have some confusion about the first question. If I had only one photo and 400 participants, there would be no random effect, since there is no repeated measure. Now I have 50 images, each rated by all the participants, so I have repeated measure within them. Hence I should add participant_id , but why also photo_no? Note: each photo was shown to each participant, but only once. $\endgroup$
    – Rakib
    Feb 18, 2020 at 4:47
  • $\begingroup$ @Rakib this choice depends on your knowledge of the subject matter. A random effect for photo_no could allow for photos with the same valence_score to have different mean values for share when the predictors associated with the participants are taken into account, potentially improving the precision of your estimates. To correct for differences among 50 photos, modeling a random effect would be more efficient than treating them as 50 fixed effects. If you think that all that matters about a photo is its valence_score then a random or fixed effect for photo_no would not be necessary. $\endgroup$
    – EdM
    Feb 20, 2020 at 2:08
  • $\begingroup$ Could you explain the last sentence and/or point to some tutorials where I can learn more? I want to study how valence_score affects the outcome. Secondly, it seems to me to include that random effects, valence_score should be added instead of photo_no? Finally, to be even more precise, there have to be random slope terms in addition to random intercepts? I tried that for participant_id, but not enough data, as you pointed out. Could you plz also refer to some resources on how to estimate the required number of participants for mixed models before conducting an experiment? $\endgroup$
    – Rakib
    Feb 20, 2020 at 18:54
  • 1
    $\begingroup$ @Rakib a random effect for photo would allow for different photos with the same valence_score to have different mean values for share after the other factors are also taken into account. So adding a random effect for photo_id would give an overall fixed effect for valence_score, but you could gain power by allowing for differences around that shared fixed effect for individual photos having the same score. I would suggest adding that random effect. The simr package in R provides power analysis for mixed models. $\endgroup$
    – EdM
    Feb 20, 2020 at 19:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.