1
$\begingroup$

My experiment has the following factors:

  • "year" (fixed) with levels "1" and "2"
  • "age" (fixed) with levels "old" and "new"
  • "treatment" (fixed) with levels "control" and "treated"
  • "block" (random) with levels ("1_1", "2_1", "3_1", "1_2", "2_2", "3_2")
  • "variety" (random) with levels ("A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L")

"Year", "age", and "treatment" are crossed, "id" is nested within "age", and "block" is nested within "year":

$$year \times block(year) \times age \times treatment \times id(age)$$

A visualization of my data:

enter image description here

I am analyzing this experiment with generalized linear mixed effects models using the lmerTest library in R. In the following, I show this for an example of my data:

dat <- structure(list(value = c(579, 575, 740, 546, 692, 490, 535, 436, 
962, 474, 867, 824, 549, 954, 511, 876, 556, 462, 443, 620, 573, 
626, 664, 975, 545, 846, 868, 973, 744, 597, 685, 768, 535, 460, 
585, 589, 1085, 755, 1004, 651, 971, 612, 682, 711, 652, 502, 
603, 609, 677, 604, 604, 488, 441, 742, 753, 532, 551, 432, 614, 
556, 615, 639, 718, 550, 642, 845, 426, 791, 499, 603, 606, 570, 
518, 638, 594, 598, 401, 479, 503, 738, 929, 480, 891, 606, 705, 
679, 941, 758, 623, 446, 619, 1077, 873, 826, 1259, 892, 890, 
470, 793, 617, 611, 501, 729, 549, 851, 460, 1032, 802, 526, 
928, 582, 790, 680, 759, 558, 367, 617, 589, 394, 529, 825, 536, 
748, 757, 655, 593, 767, 639, 675, 872, 665, 787, 660, 924, 1031, 
590, 813, 783, 736, 472, 629, 504, 898, 420, 607, 781, 687, 761, 
1083, 770, 508, 598, 923, 542, 548, 591, 526, 537, 526, 570, 
632, 919, 797, 637, 855, 811, 613), year = structure(c(2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 
1L, 1L, 1L, 1L, 1L), levels = c("1", "2"), class = "factor"), 
    age = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L), levels = c("old", "new"), class = "factor"), 
    treatment = c("control", "treated", "control", "treated", 
    "control", "treated", "treated", "control", "control", "treated", 
    "treated", "control", "treated", "control", "treated", "control", 
    "control", "treated", "treated", "control", "treated", "control", 
    "treated", "control", "treated", "control", "control", "treated", 
    "control", "treated", "control", "treated", "control", "treated", 
    "treated", "control", "control", "treated", "treated", "control", 
    "control", "treated", "treated", "control", "control", "treated", 
    "treated", "control", "control", "treated", "control", "treated", 
    "treated", "control", "control", "treated", "control", "treated", 
    "treated", "control", "treated", "control", "control", "treated", 
    "treated", "control", "treated", "control", "control", "treated", 
    "treated", "control", "treated", "control", "treated", "control", 
    "treated", "control", "treated", "control", "control", "treated", 
    "treated", "treated", "control", "treated", "control", "control", 
    "treated", "treated", "control", "control", "treated", "treated", 
    "control", "control", "treated", "treated", "control", "treated", 
    "control", "treated", "control", "treated", "control", "treated", 
    "control", "control", "treated", "control", "treated", "control", 
    "treated", "control", "treated", "treated", "control", "control", 
    "treated", "treated", "control", "treated", "control", "control", 
    "treated", "treated", "control", "treated", "control", "treated", 
    "control", "treated", "control", "control", "treated", "treated", 
    "control", "treated", "control", "treated", "control", "treated", 
    "control", "treated", "control", "control", "treated", "treated", 
    "control", "control", "treated", "treated", "control", "treated", 
    "control", "control", "treated", "treated", "control", "treated", 
    "control", "control", "treated", "control", "treated", "treated", 
    "control"), block = c("1_2", "1_2", "2_2", "2_2", "3_2", 
    "3_2", "4_2", "4_2", "1_1_ru", "1_1_ru", "2_1_ru", "2_1_ru", 
    "3_1_ru", "3_1_ru", "1_2", "1_2", "2_2", "2_2", "3_2", "3_2", 
    "4_2", "4_2", "1_1_ru", "1_1_ru", "3_1_ru", "2_1_ru", "3_1_ru", 
    "3_1_ru", "1_2", "1_2", "2_2", "2_2", "3_2", "3_2", "4_2", 
    "4_2", "1_1_ru", "1_1_ru", "2_1_ru", "2_1_ru", "3_1_ru", 
    "3_1_ru", "1_2", "1_2", "2_2", "2_2", "3_2", "3_2", "4_2", 
    "4_2", "1_1_ru", "1_1_ru", "2_1_ru", "2_1_ru", "3_1_ru", 
    "3_1_ru", "1_2", "1_2", "2_2", "2_2", "3_2", "3_2", "4_2", 
    "4_2", "1_1_ru", "1_1_ru", "2_1_ru", "2_1_ru", "3_1_ru", 
    "3_1_ru", "1_2", "1_2", "2_2", "2_2", "3_2", "3_2", "4_2", 
    "4_2", "1_1_ru", "1_1_ru", "2_1_ru", "2_1_ru", "3_1_ru", 
    "1_2", "1_2", "2_2", "2_2", "3_2", "3_2", "4_2", "4_2", "1_1_ru", 
    "1_1_ru", "2_1_ru", "2_1_ru", "3_1_ru", "3_1_ru", "1_2", 
    "1_2", "2_2", "2_2", "3_2", "3_2", "4_2", "4_2", "1_1_ru", 
    "1_1_ru", "2_1_ru", "2_1_ru", "3_1_ru", "3_1_ru", "1_2", 
    "1_2", "2_2", "2_2", "3_2", "3_2", "4_2", "4_2", "1_1_ru", 
    "1_1_ru", "2_1_ru", "2_1_ru", "3_1_ru", "3_1_ru", "1_2", 
    "1_2", "2_2", "2_2", "3_2", "3_2", "4_2", "4_2", "1_1_ru", 
    "1_1_ru", "2_1_ru", "2_1_ru", "3_1_ru", "3_1_ru", "1_2", 
    "1_2", "2_2", "2_2", "3_2", "3_2", "4_2", "4_2", "1_1_ru", 
    "1_1_ru", "2_1_ru", "2_1_ru", "3_1_ru", "3_1_ru", "1_2", 
    "1_2", "2_2", "2_2", "3_2", "3_2", "4_2", "4_2", "1_1_ru", 
    "1_1_ru", "2_1_ru", "2_1_ru", "3_1_ru", "3_1_ru"), id = c("A", 
    "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", 
    "A", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", 
    "B", "B", "B", "C", "C", "C", "C", "C", "C", "C", "C", "C", 
    "C", "C", "C", "C", "C", "D", "D", "D", "D", "D", "D", "D", 
    "D", "D", "D", "D", "D", "D", "D", "E", "E", "E", "E", "E", 
    "E", "E", "E", "E", "E", "E", "E", "E", "E", "F", "F", "F", 
    "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "G", "G", 
    "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", 
    "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", 
    "H", "H", "I", "I", "I", "I", "I", "I", "I", "I", "I", "I", 
    "I", "I", "I", "I", "J", "J", "J", "J", "J", "J", "J", "J", 
    "J", "J", "J", "J", "J", "J", "K", "K", "K", "K", "K", "K", 
    "K", "K", "K", "K", "K", "K", "K", "K", "L", "L", "L", "L", 
    "L", "L", "L", "L", "L", "L", "L", "L", "L", "L")), row.names = c(NA, 
167L), class = "data.frame")

I fit the full design as:

model_full <- lmerTest::lmer(value ~ treatment*age*year*(1|block)*(1|id),
                             data = dat,
                             REML = TRUE)

summary(model_full)

Is this the correct way to write the design formula? Can I indicate that the random factor "id" is nested in the fixed factor "age"?

My model returns:

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: value ~ treatment * age * year * (1 | block) * (1 | id)
   Data: dat

REML criterion at convergence: 2012.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.42972 -0.63655 -0.08531  0.65485  3.06605 

Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept)  1881.3   43.37  
 block    (Intercept)   126.7   11.26  
 Residual             14691.6  121.21  
Number of obs: 167, groups:  id, 12; block, 7

Fixed effects:
                              Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)                     821.72      34.23   36.46  24.003  < 2e-16 ***
treatmenttreated               -240.38      40.41  144.50  -5.948 1.95e-08 ***
agenew                           56.86      48.06   48.99   1.183 0.242457    
year2                          -136.51      38.76   49.24  -3.522 0.000934 ***
treatmenttreated:agenew         116.85      57.58  144.14   2.029 0.044248 *  
treatmenttreated:year2           88.92      53.45  144.27   1.664 0.098378 .  
agenew:year2                   -101.57      53.92  144.20  -1.884 0.061597 .  
treatmenttreated:agenew:year2   -18.60      75.92  144.06  -0.245 0.806794    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) trtmnt agenew year2  trtmn: trtm:2 agnw:2
tretmnttrtd -0.590                                          
agenew      -0.687  0.420                                   
year2       -0.647  0.521  0.438                            
trtmnttrtd:  0.414 -0.702 -0.608 -0.366                     
trtmnttrt:2  0.446 -0.756 -0.318 -0.689  0.530              
agenew:yer2  0.442 -0.374 -0.649 -0.684  0.542  0.495       
trtmnttr::2 -0.314  0.532  0.461  0.485 -0.758 -0.704 -0.710

I interpret the significant (P = 0.044) interaction between "treatment" and "age" as that the change of "value" between "control" and "treated" conditions depends on whether a sample belongs to the "old" or "new" level. How would I compare this change in response to treatment between the "old" and "new" levels? I am interested in the direction of the change and whether it is smaller/larger. Could I generate subsets of my data for the "old" and "new" conditions and use generalized linear mixed-effects models on these to test the effect of "treatment"?

dat_old <- subset(dat, age == "old")
dat_new <- subset(dat, age == "new")

model_old <- lmerTest::lmer(value ~ treatment*year*(1|block)*(1|id),
                            data = dat_old,
                            REML = TRUE)


model_new <- lmerTest::lmer(value ~ treatment*year*(1|block)*(1|id),
                            data = dat_new,
                            REML = TRUE)

Also, is there a good way to add the results of these generalized linear mixed-effects models to a plot like the one above?

$\endgroup$
0

1 Answer 1

4
$\begingroup$

I have some comments before looking at the treatment effect comparison between "old" and "new".

This is a linear mixed model (a linear regression with random effects), not a generalized linear mixed model:

  • linear model: $\operatorname{E}(Y | X) = X\beta$
  • generalized linear model: $g(\operatorname{E}(Y | X)) = X\beta$ where $g$ is the (non-linear) link function

Furthermore, {lme4} ignores the product terms between the fixed effects and the random effects, and between the two random effects. This is clear in the summary output; see the "Random effects" and the "Fixed effects" parts. The actual LMM formula is:

value ~ treatment * age * year + (1 | block) + (1 | id)
  • The mean structure, $\operatorname{E}(Y | X)$, is treatment * age * year.
  • There are two crossed random effects, block and id.

Nesting describes a structuring of random effects (actually, of the variability attributed to these random effects). A random factor (id) cannot be nested within a fixed factor (age). You can let old and new id's have different variances by replacing (1 | id) with (0 + age || id). This doesn't improve the model. Similarly, block is not nested in year.

Most importantly, you don't need to split the data by age and fit two separate models to answer questions about the treatment effect for new vs old individuals. I use the {marginaleffects} package to illustrate how to do this analysis. You'll need to decide how to deal with the the year covariate as it interacts with treatment and age. Do you want to look at the treatment effect for years 1 and 2 separately, average over the two years, or compare the treatment effect between years?

I interpret the significant (P = 0.044) interaction between "treatment" and "age" as that the change of "value" between "control" and "treated" conditions depends on whether a sample belongs to the "old" or "new" level.

Not quite. The interaction term treatmenttreated:agenew corresponds to this comparison:

$$ \begin{aligned} \left(\operatorname{E}Y_{\text{treated, new, year=1}} - \operatorname{E}Y_{\text{control, new, year=1}}\right) - \left(\operatorname{E}Y_{\text{treated, old, year=1}} - \operatorname{E}Y_{\text{control, old, year=1}}\right) \end{aligned} $$

The comparison is complex due to the three-way interaction. First, it concerns the expected response in year 1 only. The interaction itself is a difference in differences: it tells you whether there is a difference in the expected improvement due to treatment when we compare new and old individuals.

Let's do this calculation explicitly with the help of {marginaleffects}. See also the tutorial on Comparisons in the Marginal Effects Zoo (the package documentation).

# Calculate the treatment effect, treated - control,
# for each age in year 1.
comparisons(
  fit,
  variables = "treatment",
  newdata = datagrid(year = "1", age = c("old", "new")),
  re.form = NA
)
#>       Term          Contrast year age Estimate Std. Error     z Pr(>|z|)
#>  treatment treated - control    1 old     -240       40.4 -5.95   <0.001
#>  treatment treated - control    1 new     -124       41.0 -3.01   0.0026

# Compare the treatment effect between old (b1) and new (b2) in year 1.
comparisons(
  fit,
  hypothesis = "b2 - b1 = 0",
  variables = "treatment",
  newdata = datagrid(year = "1", age = c("new", "old")),
  re.form = NA
)
#>     Term Estimate Std. Error     z Pr(>|z|)   S 2.5 % 97.5 %
#>  b2-b1=0     -117       57.6 -2.03   0.0424 4.6  -230     -4

Same as the treatmenttreated:agenew row in the summary output.

Let's also look at how the change in response to treatment between old and new changes itself from year 1 to year 2.

# Compare the treatment and the control (ie. estimate the difference `treated - control`
# for each combination of year and age.
comparisons(
  fit,
  variables = "treatment",
  newdata = datagrid(year = c("1", "2"), age = c("old", "new")),
  re.form = NA
)
#>       Term          Contrast year age Estimate Std. Error     z Pr(>|z|)
#>  treatment treated - control    1 old   -240.4       40.4 -5.95   <0.001
#>  treatment treated - control    1 new   -123.5       41.0 -3.01   0.0026
#>  treatment treated - control    2 old   -151.5       35.0 -4.33   <0.001
#>  treatment treated - control    2 new    -53.2       35.0 -1.52   0.1283

# Suppose the four rows/comparisons in the table above are labeled b1 to b4, and
# formulate a more complex contrast as a difference in differences.
# Of course, each b term itself is a difference between treated and controls.
comparisons(
  fit,
  hypothesis = "(b4 - b3) - (b2 - b1) = 0",
  variables = "treatment",
  newdata = datagrid(year = c("1", "2"), age = c("old", "new")),
  re.form = NA
)
#>               Term Estimate Std. Error      z Pr(>|z|)   S 2.5 % 97.5 %
#>  (b4-b3)-(b2-b1)=0    -18.6       75.9 -0.245    0.806 0.3  -167    130

Hopefully, you are not surprised that this computation gives the same result as the treatmenttreated:agenew:year2 term in the summary table. If you are surprised and want to understand this better, you'll have to read up on design matrices and the different ways to parametrize a (generalized) linear model. We're using a treatment contrast matrix, which is the default in R.

However, there is no need to do this reading! Instead, learn to express the comparisons you are interested in as a combination of treatment, age and year. For example, you might be interested in the treatment effect for "old" individuals averaged over years 1 and 2. (This comparison doesn't correspond exactly to one of the coefficients in the summary table.)

Now, with this great power comes great responsibility. Unless the comparisons you make are specified in advance (ie. the analysis is pre-registered), if you look at multiple contrasts in a kind of exploratory analysis to discover "significant" ones, you will have to apply a correction for multiple hypothesis testing.

$\endgroup$

Your Answer

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

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