3
$\begingroup$

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
$\endgroup$
3
  • 2
    $\begingroup$ 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 <- lmer(y ~ maturity * type + block + (1|run) + (1|run:block)) where the first three factors are your fixed ones and the second two the random ones, since block has to be a fixed effect you can't really say it is nested within run, but you can add a possible interaction between both to your model like shown. $\endgroup$ Oct 27, 2020 at 7:10
  • 1
    $\begingroup$ @ThomasBaumgartner Of note, using standard R's formula, a/b expands to a + a:b. There's also an ongoing discussion at stats.stackexchange.com/q/286280/930. $\endgroup$
    – chl
    Oct 27, 2020 at 11:07
  • $\begingroup$ 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. $\endgroup$ Oct 27, 2020 at 12:09

1 Answer 1

5
$\begingroup$

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")
$\endgroup$
1
  • 1
    $\begingroup$ Does this answer your question ? If so, please consider marking it as accepted. If not, please let us know why. Thanks ! $\endgroup$ Nov 3, 2020 at 13:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.