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.

Updated June 30, 2015 21:09 PM

Updated July 27, 2015 13:08 PM

Updated September 28, 2018 14:19 PM

- Serverfault Query
- Superuser Query
- Ubuntu Query
- Webapps Query
- Webmasters Query
- Programmers Query
- Dba Query
- Drupal Query
- Wordpress Query
- Magento Query
- Joomla Query
- Android Query
- Apple Query
- Game Query
- Gaming Query
- Blender Query
- Ux Query
- Cooking Query
- Photo Query
- Stats Query
- Math Query
- Diy Query
- Gis Query
- Tex Query
- Meta Query
- Electronics Query
- Stackoverflow Query
- Bitcoin Query
- Ethereum Query