I have a split splot full factorial design : 3 blocks, each block contains 4 plots for factor E and each plot is divided into 3 subplots for factor F. I use a mixed effect model with random effects on both blocks and plot nested in blocks.

> dat
   num_bloc F      E      y
1     bloc1 T      B 0.0695
2     bloc2 T      B 0.1540
3     bloc3 T      B 0.0634
4     bloc1 C      B 0.0770
5     bloc2 C      B 0.0746
6     bloc3 C      B 0.1020
7     bloc1 P      B 0.0825
8     bloc2 P      B 0.0559
9     bloc3 P      B 0.0832
10    bloc1 T   B.Br 0.0891
11    bloc2 T   B.Br 0.1050
12    bloc3 T   B.Br 0.1150
13    bloc1 C   B.Br 0.1580
14    bloc2 C   B.Br 0.0989
15    bloc3 C   B.Br 0.1510
16    bloc1 P   B.Br 0.1020
17    bloc2 P   B.Br 0.0751
18    bloc3 P   B.Br 0.0655
19    bloc1 T    B.S 0.1020
20    bloc2 T    B.S 0.0755
21    bloc3 T    B.S 0.0631
22    bloc1 C    B.S 0.0705
23    bloc2 C    B.S 0.0782
24    bloc3 C    B.S 0.0751
25    bloc1 P    B.S 0.0826
26    bloc2 P    B.S 0.0800
27    bloc3 P    B.S 0.0996
28    bloc1 T B.S.Br 0.0850
29    bloc2 T B.S.Br 0.0688
30    bloc3 T B.S.Br 0.0727
31    bloc1 C B.S.Br 0.0762
32    bloc2 C B.S.Br 0.0880
33    bloc3 C B.S.Br 0.0751
34    bloc1 P B.S.Br 0.0694
35    bloc2 P B.S.Br 0.0619
36    bloc3 P B.S.Br 0.0627

full0 <- lme(y~F*E,control=lmeControl(opt = "optim"), 
                                contrasts = list(F='contr.treatment'),na.action="na.exclude")

Because I found differences in residual variances between groups I updated the model as follows:

full <- update(full0,weight=varIdent(form = ~ 1 | E*F))

I also test the main effect model since interaction is not significant

main0 <- lme(y~F+E,control=lmeControl(opt = "optim"), 
                                contrasts = list(F='contr.treatment'),na.action="na.exclude")
main <- update(main0,weight=varIdent(form = ~ 1 | E*F)) 

     Model df       AIC       BIC   logLik   Test  L.Ratio p-value
full     1 26 -174.1895 -133.0180 113.0948                        
main     2 20 -175.5688 -143.8985 107.7844 1 vs 2 10.62064  0.1008

I am interested in the contrasts comparing treatments versus control for factor F. I can calculate them for the 4 models:

> contrast(emmeans(full0,~F),method="trt.vs.ctrl")
NOTE: Results may be misleading due to involvement in interactions
 contrast estimate      SE df t.ratio p.value
 C - T     0.00513 0.00856 16  0.599  0.7690 
 P - T    -0.01189 0.00856 16 -1.390  0.3123 

Results are averaged over the levels of: E 
P value adjustment: dunnettx method for 2 tests 

> contrast(emmeans(full,~F),method="trt.vs.ctrl")
NOTE: Results may be misleading due to involvement in interactions
 contrast estimate      SE df t.ratio p.value
 C - T     0.00513 0.00973 16  0.527  0.8110 
 P - T    -0.01189 0.00905 16 -1.314  0.3482 

Results are averaged over the levels of: E 
P value adjustment: dunnettx method for 2 tests 

> contrast(emmeans(main0,~F),method="trt.vs.ctrl")
 contrast estimate      SE df t.ratio p.value
 C - T     0.00513 0.00902 22  0.568  0.7857 
 P - T    -0.01189 0.00902 22 -1.318  0.3392 

Results are averaged over the levels of: E 
P value adjustment: dunnettx method for 2 tests 

> contrast(emmeans(main,~F),method="trt.vs.ctrl")
 contrast  estimate      SE df t.ratio p.value
 C - T     0.000686 0.00469 22  0.146  0.9773 
 P - T    -0.011254 0.00222 22 -5.072  0.0001 

Results are averaged over the levels of: E 
P value adjustment: dunnettx method for 2 tests 

The design is balanced with no missing data, so I do not understand why emmeans differ between the full and main model and even between main and main0 models. It seems to come from the weights given to treatments by the varIndent function. My questions is can anyone help me to understand this? And second question is which model is to be used to correctly calculate the contrasts? My feeling would be to use the main model because it drops the non significant interaction and it corrects the heteroscedasticity. However, I am confused by the difference in both the emmeans and p.values depending on the models (P-T can shift from ns to highly significant). Furthermore, for some other variables, the calculated emmeans of the main model differ much more from the experimental data. I am aware that emmeans are modelled values and not experimental data but it is not comfortable to argue for a 15% difference between two treatments (p<0.001) whereas the boxplots of experimental data do not show that!! Thanks for any help.

  • 1
    $\begingroup$ I’d say that a p value if 0.1 for an interaction isn’t necessarily justification for ignoring it. What you show is notably devoid of graphics. Look at emmip(full, F ~ E) to get an idea of what is going on, and keep in mind that the ENMs for this model are equally-weighted averages of the points you see. $\endgroup$
    – Russ Lenth
    May 25, 2019 at 2:51
  • $\begingroup$ And I was wrong. The interaction has a p value of .043 via anova(full). $\endgroup$
    – Russ Lenth
    May 25, 2019 at 3:00
  • $\begingroup$ Actually, the significance of the interaction depends on how we assess it : anova(full) indeed gives a p=0.0428, full versus main model comparison using the Lratio test gives p=0.1 and finally Anova(full,main, type=3) gives p=0.00024 from the Chisquare test. Which option to use?? $\endgroup$
    – Chris
    May 27, 2019 at 19:00
  • $\begingroup$ Additionally, the anova(full) or anova(full,main) or Anova(full,3) test for all interactions. I aim planned comparisons and I am not interested in comparing all treatment combinations. Therefore I consider significant interaction if one of the p-values for the specific contrasts P-T|E is p<0.05. Does it make sense? $\endgroup$
    – Chris
    May 27, 2019 at 19:19
  • $\begingroup$ If you've found the answer helpful, please don't to forget to upvote and accept it - it seems @rvl put a notable effort in it. $\endgroup$ May 29, 2019 at 12:32

1 Answer 1


First, I note that you used method = "ML". That is fine for model comparisons, but when it comes time to do post hoc comparisons of means, I recommend re-fitting with method = REML" because this reduces biases in the estimates of the fixed effects.

Second, a graph of the full model is really useful for understanding the nature of the interaction and whether it can be set aside.

fullR <- update(full, method = "REML")
emmip(fullR, F ~ E)

Interaction plot of predictions

Note that the C level of F comes out as the highest, intermediate, and lowest prediction, depending on the level of E. It seems like not a good idea then to average over E to obtain estimated marginal means for F, and comparisons thereof. Instead, we should compare them separately for each E:

emm <- emmeans(fullR, ~ F | E)
emm            ## look at the estimates
pairs(emm)     ## compare them

(Results not shown to save space)

  • $\begingroup$ OK to refit the model with REML, I agree. I also agree that looking at plots help to decide if an interaction makes sense. Let's assume we believe the interaction is NOT significant. My strategy is to choose the main model because it more parsimonious (has less df) and avoid a possible overparameterization. Is it OK? $\endgroup$
    – Chris
    May 27, 2019 at 19:09
  • $\begingroup$ I can’t put myself in your shoes and be willing to ignore a subjectively strong interaction. Yes, interactions make things more complicated, but if that’s what the situation demands, that’s what you do. But it’s your research; you decide. $\endgroup$
    – Russ Lenth
    May 27, 2019 at 21:16

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.