How to specify a repeated measure in a Mixed Effect ANOVA?

Question

I am trying to run a Repeated Measures Mixed Effect ANOVA in R. However, the syntax of the code is really tricky and it seems like there is not a lot of clarity on how to run it as I do my research.

My Fixed effects:

1. soil_type = 2 levels
2. treatment = 4 levels

Repeated measure:

1. days = 3 levels

Random effect:

1. rep = 3 levels (I know that this is probably too low of a replication for a random effect, unsure how big of an issue it will be)

The code I tried first was with the lmer package:

cl.mod <- lmer(cl_conc ~ soil_type * treatment + (days|rep), data = leach.conc, REML = FALSE)

and I received this output:

Type III Analysis of Variance Table with Satterthwaite's method
                    Sum Sq Mean Sq NumDF  DenDF  F value    Pr(>F)    
soil_type             3681  3680.8     1 67.858  20.9813 2.041e-05 ***
treatment            54038 18012.6     3 67.858 102.6751 < 2.2e-16 ***
soil_type:treatment    925   308.2     3 67.858   1.7567    0.1637    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

It looks right to me, but I do not think that I have the "days" (repeated measure) variable in the right location in the code. I am unsure where exactly this should go. Should the code reflect this instead?

cl.mod <- lmer(cl_conc ~ soil_type * treatment *days + (1|rep), data = leach.conc, REML = FALSE)

Alternatively, I have used the rstatix package running this code:

   cl.conc.aov <- anova_test(
  data = leach.conc, dv = cl_conc_trans, wid = core_id,
  between = c(soil_type, treatment), within = c(days))

but then I am unsure where the random variable of "rep" goes.

Any advice is appreciated! Just trying to learn some stats.

Below is the dataset that I am using. Each core_id is a combination of soil_type, treatment and replicate which was measured over 3 time periods. I intend to see how soil_type and treatment affects cl_conc and using time as a repeated measure.

 soil_type treatment days rep core_id cl_conc
1     Medium     Solid    4   1     MS1   18.10
2     Medium     Solid   11   1     MS1   18.10
3     Medium     Solid   18   1     MS1   17.40
4     Medium    Liquid    4   1     ML1   77.10
5     Medium    Liquid   11   1     ML1   81.40
6     Medium    Liquid   18   1     ML1   66.80
7     Medium       KCl    4   1     MK1   19.40
8     Medium       KCl   11   1     MK1   22.30
9     Medium       KCl   18   1     MK1   36.90
10    Medium   Control    4   1     MC1    1.90
11    Medium   Control   11   1     MC1    1.20
12    Medium   Control   18   1     MC1    0.60
13     Fine      Solid    4   1     FS1   27.80
14     Fine      Solid   11   1     FS1   28.30
15     Fine      Solid   18   1     FS1   28.30
16     Fine     Liquid    4   1     FL1  107.80
17     Fine     Liquid   11   1     FL1  150.30
18     Fine     Liquid   18   1     FL1   94.60
19     Fine        KCl    4   1     FK1   84.80
20     Fine        KCl   11   1     FK1   53.40
21     Fine        KCl   18   1     FK1   51.90
22     Fine    Control    4   1     FC1    9.10
23     Fine    Control   11   1     FC1    4.25
24     Fine    Control   18   1     FC1    1.90
25    Medium     Solid    4   2     MS2   19.80
26    Medium     Solid   11   2     MS2   20.70
27    Medium     Solid   18   2     MS2   20.50
28    Medium    Liquid    4   2     ML2  102.00
29    Medium    Liquid   11   2     ML2   56.70
30    Medium    Liquid   18   2     ML2   47.40
31    Medium       KCl    4   2     MK2   33.40
32    Medium       KCl   11   2     MK2   15.30
33    Medium       KCl   18   2     MK2   19.90
34    Medium   Control    4   2     MC2    2.00
35    Medium   Control   11   2     MC2    1.20
36    Medium   Control   18   2     MC2    0.80
37     Fine      Solid    4   2     FS2   37.10
38     Fine      Solid   11   2     FS2   39.80
39     Fine      Solid   18   2     FS2   34.80
40     Fine     Liquid    4   2     FL2   81.90
41     Fine     Liquid   11   2     FL2   67.50
42     Fine     Liquid   18   2     FL2   56.00
43     Fine        KCl    4   2     FK2   41.10
44     Fine        KCl   11   2     FK2   38.30
45     Fine        KCl   18   2     FK2   30.90
46     Fine    Control    4   2     FC2   12.40
47     Fine    Control   11   2     FC2    6.00
48     Fine    Control   18   2     FC2    3.10
49    Medium     Solid    4   3     MS3   27.80
50    Medium     Solid   11   3     MS3   27.80
51    Medium     Solid   18   3     MS3   24.90
52    Medium    Liquid    4   3     ML3   79.70
53    Medium    Liquid   11   3     ML3   65.50
54    Medium    Liquid   18   3     ML3   55.20
55    Medium       KCl    4   3     MK3   13.50
56    Medium       KCl   11   3     MK3   20.40
57    Medium       KCl   18   3     MK3   24.70
58    Medium   Control    4   3     MC3    1.60
59    Medium   Control   11   3     MC3    1.20
60    Medium   Control   18   3     MC3    0.70
61     Fine      Solid    4   3     FS3   42.70
62     Fine      Solid   11   3     FS3   40.50
63     Fine      Solid   18   3     FS3   30.10
64     Fine     Liquid    4   3     FL3  121.20
65     Fine     Liquid   11   3     FL3   73.60
66     Fine     Liquid   18   3     FL3   38.00
67     Fine        KCl    4   3     FK3   53.00
68     Fine        KCl   11   3     FK3   38.50
69     Fine        KCl   18   3     FK3   22.30
70     Fine    Control    4   3     FC3    4.70
71     Fine    Control   11   3     FC3    1.90
72     Fine    Control   18   3     FC3    0.85

Boxplot represents a visualization of how cl_conc changes through time grouped together by soil type and treatment. Each boxplot represents the 3 reps of each unique core_id in that given time period.

Answer 1

The core_id should be your random effect in this analysis. Those are the units over which there are multiple observations. By my count that's 24 "clusters" of observations, with 3 observations (on days 4, 11, and 18) within each cluster. That's OK.

Unless there's something systematic in the three different rep values (which are implicitly included in the core_id values), there's no need to include rep as a predictor at all. If there is something systematically different among the rep values, treat it as a fixed effect.

Your model did not have any fixed effect term for days. It needs that, in some form. I don't think that a random slope for days, as you attempted in the first model, will be reasonable with this data set.

You don't have enough data to evaluate all of the interaction terms that you would like to. You have 72 observations, which typically will allow you to fit about 4 to 7 coefficients without risk of overfitting. Your 4-level treatment predictor requires fitting 3 coefficients, your 2-level soil_type adds another, and your days predictor will require either 1 (if you treat is as linear) or 2 coefficients (if you treat it as a factor, which might make more sense). That already brings you up to 5 or 6 coefficients to estimate. Plus, you have to estimate the variance of the random (intercept) effect.

For example, a full 3-way interaction among days and treatment and soil_type, as in your second model, would require adding the following 2-way coefficients to the above, even if you treated days as linear: 3 more soil_type:treatment coefficients, 1 soil_type:days coefficient, and 3 days:treatment coefficients. You would also need to estimate 3 further 3-way coefficients. You probably don't have enough data to do that properly.

The overall structure thus probably needs to be limited to additive main effects of treatment, soil_type and days, with a random intercept for core_id.

Finally, your continuous-outcome variable doesn't seem to fit the linear-regression assumption of homogeneous variance. The within-condition variance increases massively at higher Cl outcomes. This type of outcome data might best be analyzed in a log scale, or with a generalized linear model under a log link. The following seems to work pretty well on your data:

cl.mod.additive <- lmer(log10(cl_conc) ~ soil_type +
    treatment + factor(days) +   
    (1|core_id), data=leach.conc)