Confusing results from lsmeans and pairwise tests

Question

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+
     affiliative_score_group*valence_group*condition+
     self_enhancing_score_group*valence_group*condition+
     self_defeating_score_group*valence_group*condition+
     aggressive_score_group*valence_group*condition+ 
     (1|photo), data = picshare_df, REML = TRUE)

Question

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.

From

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:

lsmeans(lmer.model,list(trt.vs.ctrl1~affiliative_score_group))
 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.

Question

What I am doing wrong or not understanding?

Update

With a simpler model, I get this:

lmer1.model=lmer(share ~ gender+ age_group + overall_photo_share_freq_group+
     affiliative_score_group*valence_group*condition+
     (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)?

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

affiliative_score_group*valence_group*condition

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.