Block effect in mixed model

I have a problem i really can't wrap my head around.

I have a dataset from a randomised block field trial conducted over two years. DV is Yield and IV's are Accession, SD and Year. The random effects are the nested blocks within years, when both years are included in the model, and Block, when only one year is included.

I have no problem running a mixed model anova with the combined data from the two years, but when i try to run a model for each year they behave differently though my data looks very similar.

2016 dataframe:

> stat16
Block Accession SD   Yield
1      1        H5 S1  909.20
2      2        H5 S1  513.60
3      3        H5 S1 2565.07
4      1        H5 S2 2691.13
5      2        H5 S2 2665.33
6      3        H5 S2 2490.13
7      1        H5 S3 2808.27
8      2        H5 S3 2131.73
9      3        H5 S3 2499.53
10     1      Lisa S1  815.67
11     2      Lisa S1  439.73
12     3      Lisa S1 1369.67
13     1      Lisa S2 2738.80
14     2      Lisa S2 2942.27
15     3      Lisa S2 1811.60
16     1      Lisa S3 2961.00
17     2      Lisa S3 2197.33
18     3      Lisa S3 2710.73
19     1        KR S1  809.13
20     2        KR S1 1054.60
21     3        KR S1 3145.13
22     1        KR S2 3829.73
23     2        KR S2 3946.13
24     3        KR S2 2814.87
25     1        KR S3 3617.27
26     2        KR S3 2919.87
27     3        KR S3 2798.07

> str(stat16)
'data.frame':   27 obs. of  4 variables:
$$Block : Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1 ...$$ Accession: Factor w/ 3 levels "H5","KR","Lisa": 1 1 1 1 1 1 1 1 1 3 ...
$$SD : Factor w/ 3 levels "S1","S2","S3": 1 1 1 2 2 2 3 3 3 1 ...$$ Yield    : num  909 514 2565 2691 2665 ...

2017 dataframe:

> stat17
Block Accession SD   Yield
28     1        H5 S1 3632.47
29     2        H5 S1 1393.47
30     3        H5 S1  839.87
31     1        H5 S2 2595.67
32     2        H5 S2 2333.00
33     3        H5 S2 1713.27
34     1        H5 S3 2838.40
35     2        H5 S3 2419.33
36     3        H5 S3 2227.40
37     1      Lisa S1 2435.47
38     2      Lisa S1  716.60
39     3      Lisa S1  647.73
40     1      Lisa S2 2946.80
41     2      Lisa S2 2207.47
42     3      Lisa S2 1861.80
43     1      Lisa S3 2791.47
44     2      Lisa S3 1125.47
45     3      Lisa S3 2297.53
46     1        KR S1 1526.40
47     2        KR S1  925.80
48     3        KR S1      NA
49     1        KR S2 2397.33
50     2        KR S2 2346.67
51     3        KR S2 1714.13
52     1        KR S3 2473.33
53     2        KR S3 1861.13
54     3        KR S3 2252.93

> str(stat17)
'data.frame':   27 obs. of  4 variables:
$$Block : Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1 ...$$ Accession: Factor w/ 3 levels "H5","KR","Lisa": 1 1 1 1 1 1 1 1 1 3 ...
$$SD : Factor w/ 3 levels "S1","S2","S3": 1 1 1 2 2 2 3 3 3 1 ...$$ Yield    : num  3632 1393 840 2596 2333 ...

The following two models returns quite different results as one the p-value of the block effect in 2016 seems to be 1, which i don't believe.

> yield2016 <- lmer(Yield ~ Accession*SD + (1|Block), data=stat16)
> ranova(yield2016)
ANOVA-like table for random-effects: Single term deletions

Model:
Yield ~ Accession + SD + (1 | Block) + Accession:SD
npar  logLik    AIC        LRT Df Pr(>Chisq)
<none>        11 -148.12 318.25
(1 | Block)   10 -148.12 316.25 5.6843e-14  1          1

> yield2017 <- lmer(Yield ~ Accession*SD + (1|Block), data=stat17)
> ranova(yield2017)
ANOVA-like table for random-effects: Single term deletions

Model:
Yield ~ Accession + SD + (1 | Block) + Accession:SD
npar  logLik    AIC    LRT Df Pr(>Chisq)
<none>        11 -137.58 297.16
(1 | Block)   10 -141.59 303.18 8.0126  1   0.004645 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

I use the lmer_Test package.

How is only one of the mixed models returning a p-value? The only difference I can see in the two datasets besides the values of course is a missing value in the 2017 dataset, which i tried to introduce into the 2016 dataset just to test if that was contributing to the issue, without any changes in performance.

I hope someone can shed some light over the issue :) Thanks a lot.

Dennis

The issue arise because your models are estimating a variance term (the variance of the random intercepts) from just 3 cluster (Blocks), and because the differences across blocks are smaller in 2016, as indicated by the standard deviations

> sd(tapply(stat16$$Yield, stat16$$Block, mean))
[1] 193.4393
> sd(tapply(stat17$$Yield, stat17$$Block, mean, na.rm=T))
[1] 535.5668

As a result when considering 2016 alone, the model does not have enough evidence to believe that different Blocks have different average values of Yield, and the estimated standard deviation turns out to be zero:

> VarCorr(yield2016)
Groups   Name        Std.Dev.
Block    (Intercept)   0.00
Residual             641.15

This is different from what you get in 2017, where the estimated standard deviation is very close to the standard deviation computed above

> VarCorr(yield2017)
Groups   Name        Std.Dev.
Block    (Intercept) 536.05
Residual             505.32

(As a test, if you artificially increased the between-block variance, the model would estimate a non-zero variance also for 2016)

> stat16$$Yield[stat16$$Block==1] <- stat16$$Yield[stat16$$Block==1] + 800
> yield2016 <- lmer(Yield ~ Accession+SD + (1|Block), data=stat16)
> VarCorr(yield2016)
Groups   Name        Std.Dev.
Block    (Intercept) 494.07
Residual             646.92

If you do not believe that the variance across block is zero in 2016, I would suggest to use a Bayesian approach and incorporate your prior knowledge in the model (you could do that with rstanarm package and the stan_glmer() function)

• Hi matteo. Thanks for you reply. Very well explained! Would this mean that it is actually not possible to model the years separately? – den nis Nov 1 '18 at 12:32
• No, I think you can model them separately, but with only three Blocks the estimates random intercept variance will not be very precise, and in some cases may result equal to zero. – matteo Nov 1 '18 at 14:03