Mixed model with nesting

Question

I am new to working with mixed models in R and am hoping someone can tell me if this syntax is correct.

m <- lmer(y ~ maturity * type + (1|run/block))

Where block is a fixed effect nested in the random effect of run
Also interested in the main effects and interactions between maturity and type.

   Sample             maturity        type        rundate       run           
Length:120         Min.   :29.06     1  :42   5/13/2020:30   Length:120        
Class :character   1st Qu.:31.82     2  :36   6/16/2020:30   Class:character  
Mode  :character   Median :34.50     3  :42   6/2/2020 :30   Mode  :character  
                   Mean   :34.05              6/9/2020 :30                     
                   3rd Qu.:35.88                                                
                   Max.   :38.43                                                
 block      rep              y      
 A:30   Length:120         Min.   : 8.87  
 B:30   Class :character   1st Qu.:15.44  
 C:30   Mode  :character   Median :18.67  
 D:30                      Mean   :19.33  
                           3rd Qu.:23.89  
                           Max.   :32.26

    Sample  rep maturity  type  rundate   run  block  y
1   35795   1   38.38521    1   5/13/2020   1   B   27.54
2   35795   2   38.38521    1   5/13/2020   1   B   24.58
3   35795   3   38.38521    1   5/13/2020   1   B   29.51
4   28888   1   37.72189    3   5/13/2020   1   B   16.58
5   28888   2   37.72189    3   5/13/2020   1   B   17.56
6   28888   3   37.72189    3   5/13/2020   1   B   16.29
7   42582   1   36.17173    2   5/13/2020   1   B   13.22
8   42582   2   36.17173    2   5/13/2020   1   B   12.7
9   42582   3   36.17173    2   5/13/2020   1   B   15.08
10  34759   1   31.20564    1   5/13/2020   1   B   17.77
11  34759   2   31.20564    1   5/13/2020   1   B   18.17
12  34759   3   31.20564    1   5/13/2020   1   B   18.7
13  30623   1   29.16932    2   5/13/2020   1   B   15.5
14  30623   2   29.16932    2   5/13/2020   1   B   17.02
15  30623   3   29.16932    2   5/13/2020   1   B   18.66
16  31413   1   32.928      1   5/13/2020   1   D   27.11
17  31413   2   32.928      1   5/13/2020   1   D   26.62
18  31413   3   32.928      1   5/13/2020   1   D   23.88
19  31414   1   35.48105    3   5/13/2020   1   D   18.94
20  31414   2   35.48105    3   5/13/2020   1   D   16.37

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: ust8 ~ maturity * type + block + (1 | run) + (1 | run:block)
   Data: wpc

REML criterion at convergence: 583.7

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.92238 -0.76620 -0.05373  0.69237  2.30801 

Random effects:
 Groups    Name        Variance Std.Dev.
 run:block (Intercept) 0.000    0.000   
 run       (Intercept) 3.678    1.918   
 Residual              7.303    2.702   
Number of obs: 120, groups:  run:block, 8; run, 4

Fixed effects:
                 Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)      -6.11204    3.19194  92.82683  -1.915   0.0586 .  
maturity          0.78035    0.08827 108.20906   8.840 1.97e-14 ***
type1            -3.55955    4.83520 108.20906  -0.736   0.4632    
type2            11.14114    4.62248 108.20906   2.410   0.0176 *  
blockB           -4.69970    0.98479 110.98838  -4.772 5.59e-06 ***
blockC            0.87225    0.94700 110.93534   0.921   0.3590    
blockD           -0.51150    1.07615 109.27831  -0.475   0.6355    
maturity:type1    0.22890    0.14187 108.20906   1.613   0.1096    
maturity:type2   -0.33706    0.13346 108.20906  -2.525   0.0130 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
           (Intr) mat     type1 type2  blockB blockC blockD mat:ty1
mat        -0.928                                                 
type1      -0.004  0.098                                          
type2       0.060 -0.020 -0.387                                   
blockB     -0.096 -0.082 -0.280 -0.290                            
blockC     -0.204  0.035 -0.359 -0.090  0.472                     
blockD     -0.071 -0.121 -0.416  0.005  0.625  0.638              
mat:type1   0.002 -0.098 -0.997  0.375  0.286  0.371  0.419       
mat:type2  -0.048  0.008  0.383 -0.997  0.283  0.081 -0.003 -0.374

Please add more Information about your experiment, about its design, some of your data (like: ```summary(data); head(data,20)``` and a summary of your model ```summary(m)```, apart from that if block is a fixed factor your model should look like this: ```m — Thomas Baumgartner, Oct 27 '20 at 07:10
@ThomasBaumgartner Of note, using standard R's formula, `a/b` expands to `a + a:b`. There's also an ongoing discussion at https://stats.stackexchange.com/q/286280/930. — chl, Oct 27 '20 at 11:07
Thank you! That makes sense. I've added the information. I hope it's sufficient - it's a little tough because it's for work and I can't have too much floating out there.And thank you for the link that is a helpful discussion. — yaynikkiprograms, Oct 27 '20 at 12:09

score 5 · Accepted Answer · answered Oct 27 '20 at 12:58

Changing your initial model from:
m1 <- lmer(y ~ maturity * type + (1|run) + (1|run:block))

as correctly mentioned by @chl (1|run/block) expands to (1|run) + (1|run:block) you can read more about that here: Crossed vs nested random effects: how do they differ and how are they specified correctly in lme4?

to your new model shown in your summary:
m2 <- lmer(y ~ maturity * type + block + (1|run) + (1|run:block))
adds block as fixed effect while keeping the nesting/interaction term of block within run, as you wanted to

however if you look at the Variance and Std.Dev in the random effects part of your model summary:

Random effects:
 Groups    Name        Variance Std.Dev.
 run:block (Intercept) 0.000    0.000

the interaction does not explain any of your residual variance and thus should be removed from your model, you probably also got a boundary (singular) fit: see ?isSingular warning there, more about this topic is explained here: Dealing with singular fit in mixed models

So your final model could look like this:
m2 <- lmer(y ~ maturity * type + block + (1|run))

Additional steps I can recommend from personal experience are

to compare your models with anova(m1,m2) and take a look at the AIC or BIC values, which indicate model fit
and performing a model critique respectively a visual inspection of your
residual plot plot(m1)
and QQ plot qqnorm(resid(m1)); qqline(resid(m1), col = "red")

Does this answer your question ? If so, please consider marking it as accepted. If not, please let us know why. Thanks ! — Thomas Baumgartner, Nov 03 '20 at 13:49

Mixed model with nesting

1 Answers1