Which emmeans to choose between full and main effect mixed effect models with heteroscedasticity corrected

by Christophe Nguyen   Last Updated May 24, 2019 22:19 PM

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"), 
                                random=~1|num_bloc/E,data=dat,method='ML',
                                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"), 
                                random=~1|num_bloc/E,data=dat,method='ML',
                                contrasts = list(F='contr.treatment'),na.action="na.exclude")
main <- update(main0,weight=varIdent(form = ~ 1 | E*F)) 

anova(full,main)
     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.



Related Questions


Updated June 30, 2015 21:09 PM

Updated May 28, 2015 00:08 AM

Updated July 27, 2015 13:08 PM

Updated September 28, 2018 14:19 PM

Updated June 21, 2019 12:19 PM