6
$\begingroup$

I am trying to run a Friedman Test in R with a repeated measure, however, my data do not qualify as an unreplicated complete block design. I am wondering what alternative test I can run as a 3-way repeated measures ANOVA and Friedman test are not appropriate.

I am interested in knowing the interaction that soil_type, treatment, and days have on the cl_conc of my subjects. There are 48 subjects that were tested over a period of 3 days (4, 11, and 18).

Description of my data:

  • cl_conc (response variable)
  • soil_type (explanatory variable 1; 4 levels)
  • treatment (explanatory variable 2; 4 levels)
  • days (repeated measure; 3 levels)
  • core_id (subject variable; 48 levels)

Data for reproducibility:

leach2 <-
  structure(
    list(
      core = c(
        "MS",
        "MS",
        "MS",
        "ML",
        "ML",
        "ML",
        "MK",
        "MK",
        "MK",
        "MC",
        "MC",
        "MC",
        "FS",
        "FS",
        "FS",
        "FL",
        "FL",
        "FL",
        "FK",
        "FK",
        "FK",
        "FC",
        "FC",
        "FC",
        "MS",
        "MS",
        "MS",
        "ML",
        "ML",
        "ML",
        "MK",
        "MK",
        "MK",
        "MC",
        "MC",
        "MC",
        "FS",
        "FS",
        "FS",
        "FL",
        "FL",
        "FL",
        "FK",
        "FK",
        "FK",
        "FC",
        "FC",
        "FC",
        "MS",
        "MS",
        "MS",
        "ML",
        "ML",
        "ML",
        "MK",
        "MK",
        "MK",
        "MC",
        "MC",
        "MC",
        "FS",
        "FS",
        "FS",
        "FL",
        "FL",
        "FL",
        "FK",
        "FK",
        "FK",
        "FC",
        "FC",
        "FC",
        "CS",
        "CL",
        "CK",
        "CC",
        "PS",
        "PL",
        "PK",
        "PC",
        "CS",
        "CL",
        "CK",
        "CC",
        "PS",
        "PL",
        "PK",
        "PC",
        "CS",
        "CL",
        "CK",
        "CC",
        "PS",
        "PL",
        "PK",
        "PC",
        "CS",
        "CL",
        "CK",
        "CC",
        "PS",
        "PL",
        "PK",
        "PC",
        "CS",
        "CL",
        "CK",
        "CC",
        "PS",
        "PL",
        "PK",
        "PC",
        "CS",
        "CL",
        "CK",
        "CC",
        "PS",
        "PL",
        "PK",
        "PC",
        "CS",
        "CK",
        "CL",
        "CC",
        "PS",
        "PL",
        "PK",
        "PC",
        "CS",
        "CL",
        "CK",
        "CC",
        "PS",
        "PL",
        "PK",
        "PC",
        "CS",
        "CL",
        "CK",
        "CC",
        "PS",
        "PL",
        "PK",
        "PC"
      ),
      core_id = c(
        "MS1",
        "MS1",
        "MS1",
        "ML1",
        "ML1",
        "ML1",
        "MK1",
        "MK1",
        "MK1",
        "MC1",
        "MC1",
        "MC1",
        "FS1",
        "FS1",
        "FS1",
        "FL1",
        "FL1",
        "FL1",
        "FK1",
        "FK1",
        "FK1",
        "FC1",
        "FC1",
        "FC1",
        "MS2",
        "MS2",
        "MS2",
        "ML2",
        "ML2",
        "ML2",
        "MK2",
        "MK2",
        "MK2",
        "MC2",
        "MC2",
        "MC2",
        "FS2",
        "FS2",
        "FS2",
        "FL2",
        "FL2",
        "FL2",
        "FK2",
        "FK2",
        "FK2",
        "FC2",
        "FC2",
        "FC2",
        "MS3",
        "MS3",
        "MS3",
        "ML3",
        "ML3",
        "ML3",
        "MK3",
        "MK3",
        "MK3",
        "MC3",
        "MC3",
        "MC3",
        "FS3",
        "FS3",
        "FS3",
        "FL3",
        "FL3",
        "FL3",
        "FK3",
        "FK3",
        "FK3",
        "FC3",
        "FC3",
        "FC3",
        "CS1",
        "CL1",
        "CK1",
        "CC1",
        "PS1",
        "PL1",
        "PK1",
        "PC1",
        "CS2",
        "CL2",
        "CK2",
        "CC2",
        "PS2",
        "PL2",
        "PK2",
        "PC2",
        "CS3",
        "CL3",
        "CK3",
        "CC3",
        "PS3",
        "PL3",
        "PK3",
        "PC3",
        "CS1",
        "CL1",
        "CK1",
        "CC1",
        "PS1",
        "PL1",
        "PK1",
        "PC1",
        "CS2",
        "CL2",
        "CK2",
        "CC2",
        "PS2",
        "PL2",
        "PK2",
        "PC2",
        "CS3",
        "CL3",
        "CK3",
        "CC3",
        "PS3",
        "PL3",
        "PK3",
        "PC3",
        "CS1",
        "CK1",
        "CL1",
        "CC1",
        "PS1",
        "PL1",
        "PK1",
        "PC1",
        "CS2",
        "CL2",
        "CK2",
        "CC2",
        "PS2",
        "PL2",
        "PK2",
        "PC2",
        "CS3",
        "CL3",
        "CK3",
        "CC3",
        "PS3",
        "PL3",
        "PK3",
        "PC3"
      ),
      soil_type = c(
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "WSL",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "NCL ",
        "ESL",
        "ESL",
        "ESL",
        "ESL",
        "FAH",
        "FAH",
        "FAH",
        "FAH",
        "ESL",
        "ESL",
        "ESL",
        "ESL",
        "FAH",
        "FAH",
        "FAH",
        "FAH",
        "ESL",
        "ESL",
        "ESL",
        "ESL",
        "FAH",
        "FAH",
        "FAH",
        "FAH",
        "ESL",
        "ESL",
        "ESL",
        "ESL",
        "FAH",
        "FAH",
        "FAH",
        "FAH",
        "ESL",
        "ESL",
        "ESL",
        "ESL",
        "FAH",
        "FAH",
        "FAH",
        "FAH",
        "ESL",
        "ESL",
        "ESL",
        "ESL",
        "FAH",
        "FAH",
        "FAH",
        "FAH",
        "ESL",
        "ESL",
        "ESL",
        "ESL",
        "FAH",
        "FAH",
        "FAH",
        "FAH",
        "ESL",
        "ESL",
        "ESL",
        "ESL",
        "FAH",
        "FAH",
        "FAH",
        "FAH",
        "ESL",
        "ESL",
        "ESL",
        "ESL",
        "FAH",
        "FAH",
        "FAH",
        "FAH"
      ),
      treatment = c(
        "Tl",
        "Tl",
        "Tl",
        "SM",
        "SM",
        "SM",
        "KCl",
        "KCl",
        "KCl",
        "Control",
        "Control",
        "Control",
        "Tl",
        "Tl",
        "Tl",
        "SM",
        "SM",
        "SM",
        "KCl",
        "KCl",
        "KCl",
        "Control",
        "Control",
        "Control",
        "Tl",
        "Tl",
        "Tl",
        "SM",
        "SM",
        "SM",
        "KCl",
        "KCl",
        "KCl",
        "Control",
        "Control",
        "Control",
        "Tl",
        "Tl",
        "Tl",
        "SM",
        "SM",
        "SM",
        "KCl",
        "KCl",
        "KCl",
        "Control",
        "Control",
        "Control",
        "Tl",
        "Tl",
        "Tl",
        "SM",
        "SM",
        "SM",
        "KCl",
        "KCl",
        "KCl",
        "Control",
        "Control",
        "Control",
        "Tl",
        "Tl",
        "Tl",
        "SM",
        "SM",
        "SM",
        "KCl",
        "KCl",
        "KCl",
        "Control",
        "Control",
        "Control",
        "Tl",
        "SM",
        "KCl",
        "Control",
        "Tl",
        "SM",
        "KCl",
        "Control",
        "Tl",
        "SM",
        "KCl",
        "Control",
        "Tl",
        "SM",
        "KCl",
        "Control",
        "Tl",
        "SM",
        "KCl",
        "Control",
        "Tl",
        "SM",
        "KCl",
        "Control",
        "Tl",
        "SM",
        "KCl",
        "Control",
        "Tl",
        "SM",
        "KCl",
        "Control",
        "Tl",
        "SM",
        "KCl",
        "Control",
        "Tl",
        "SM",
        "KCl",
        "Control",
        "Tl",
        "SM",
        "KCl",
        "Control",
        "Tl",
        "SM",
        "KCl",
        "Control",
        "Tl",
        "KCl",
        "SM",
        "Control",
        "Tl",
        "SM",
        "KCl",
        "Control",
        "Tl",
        "SM",
        "KCl",
        "Control",
        "Tl",
        "SM",
        "KCl",
        "Control",
        "Tl",
        "SM",
        "KCl",
        "Control",
        "Tl",
        "SM",
        "KCl",
        "Control"
      ),
      days = c(
        4L,
        11L,
        18L,
        4L,
        11L,
        18L,
        4L,
        11L,
        18L,
        4L,
        11L,
        18L,
        4L,
        11L,
        18L,
        4L,
        11L,
        18L,
        4L,
        11L,
        18L,
        4L,
        11L,
        18L,
        4L,
        11L,
        18L,
        4L,
        11L,
        18L,
        4L,
        11L,
        18L,
        4L,
        11L,
        18L,
        4L,
        11L,
        18L,
        4L,
        11L,
        18L,
        4L,
        11L,
        18L,
        4L,
        11L,
        18L,
        4L,
        11L,
        18L,
        4L,
        11L,
        18L,
        4L,
        11L,
        18L,
        4L,
        11L,
        18L,
        4L,
        11L,
        18L,
        4L,
        11L,
        18L,
        4L,
        11L,
        18L,
        4L,
        11L,
        18L,
        4L,
        4L,
        4L,
        4L,
        4L,
        4L,
        4L,
        4L,
        4L,
        4L,
        4L,
        4L,
        4L,
        4L,
        4L,
        4L,
        4L,
        4L,
        4L,
        4L,
        4L,
        4L,
        4L,
        4L,
        11L,
        11L,
        11L,
        11L,
        11L,
        11L,
        11L,
        11L,
        11L,
        11L,
        11L,
        11L,
        11L,
        11L,
        11L,
        11L,
        11L,
        11L,
        11L,
        11L,
        11L,
        11L,
        11L,
        11L,
        18L,
        18L,
        18L,
        18L,
        18L,
        18L,
        18L,
        18L,
        18L,
        18L,
        18L,
        18L,
        18L,
        18L,
        18L,
        18L,
        18L,
        18L,
        18L,
        18L,
        18L,
        18L,
        18L,
        18L
      ),
      cl_conc = c(
        18.1,
        18.1,
        17.4,
        77.1,
        81.4,
        66.8,
        19.4,
        22.3,
        36.9,
        1.9,
        1.2,
        0.6,
        27.8,
        28.3,
        28.3,
        107.8,
        150.3,
        94.6,
        84.8,
        53.4,
        51.9,
        9.1,
        4.25,
        1.9,
        19.8,
        20.7,
        20.5,
        102,
        56.7,
        47.4,
        33.4,
        15.3,
        19.9,
        2,
        1.2,
        0.8,
        37.1,
        39.8,
        34.8,
        81.9,
        67.5,
        56,
        41.1,
        38.3,
        30.9,
        12.4,
        6,
        3.1,
        27.8,
        27.8,
        24.9,
        79.7,
        65.5,
        55.2,
        13.5,
        20.4,
        24.7,
        1.6,
        1.2,
        0.7,
        42.7,
        40.5,
        30.1,
        121.2,
        73.6,
        38,
        53,
        38.5,
        22.3,
        4.7,
        1.9,
        0.85,
        46.5,
        52.6,
        32.9,
        2.8,
        45.1,
        1.3,
        51.2,
        2.6,
        47.59251129,
        68.3,
        38.8,
        5.4,
        34.1,
        66.7,
        23.51266468,
        0.6,
        34.2,
        55.7,
        23.8,
        5,
        42.1,
        47.9,
        44.3,
        0.8,
        56.23151874,
        81.2,
        36.1,
        1.6,
        36.3,
        48.2,
        35.6,
        1.5,
        44.8,
        80.9,
        34.66600908,
        3.1,
        33.3,
        81.5,
        20.2,
        0.4,
        40.1,
        66.8,
        24.5,
        3.6,
        39,
        68.2,
        36,
        0.303367677,
        31.1,
        23.2,
        75.7,
        0.6,
        26.2,
        45.3,
        21.3,
        0.6,
        33.76030379,
        47.5,
        20.5,
        1.1,
        28.6,
        65.9,
        18.9,
        0.2,
        30.2,
        65.5,
        23.3,
        2.7,
        23.9,
        64,
        24.7,
        0.1
      ),
      cl_load = c(
        0.058825,
        0.0543,
        0.0609,
        0.26985,
        0.26455,
        0.2171,
        0.0582,
        0.0669,
        0.119925,
        0.0057,
        0.0036,
        0.00195,
        0.09035,
        0.0849,
        0.0849,
        0.3773,
        0.4509,
        0.3311,
        0.2756,
        0.1602,
        0.1557,
        0.03185,
        0.010625,
        0.006175,
        0.06435,
        0.07245,
        0.07175,
        0.3315,
        0.19845,
        0.1659,
        0.10855,
        0.05355,
        0.064675,
        0.0065,
        0.0039,
        0.0026,
        0.1113,
        0.1393,
        0.1044,
        0.266175,
        0.185625,
        0.168,
        0.113025,
        0.105325,
        0.0927,
        0.0372,
        0.018,
        0.010075,
        0.09035,
        0.09035,
        0.0747,
        0.2391,
        0.22925,
        0.1656,
        0.0405,
        0.0663,
        0.06175,
        0.0048,
        0.0036,
        0.0021,
        0.1281,
        0.131625,
        0.097825,
        0.3939,
        0.2392,
        0.114,
        0.17225,
        0.125125,
        0.0669,
        0.015275,
        0.006175,
        0.0023375,
        0.11625,
        0.14465,
        0.0987,
        0.0084,
        0.0902,
        0.0039,
        0.1152,
        0.00715,
        0.118981278,
        0.2049,
        0.1164,
        0.0135,
        0.093775,
        0.2001,
        0.064659828,
        0.0018,
        0.0855,
        0.1671,
        0.06545,
        0.0125,
        0.094725,
        0.131725,
        0.099675,
        0.0024,
        0.177129284,
        0.2436,
        0.1083,
        0.0048,
        0.1089,
        0.15665,
        0.1068,
        0.0045,
        0.1344,
        0.2427,
        0.11266453,
        0.010075,
        0.0999,
        0.2445,
        0.0606,
        0.0012,
        0.1203,
        0.2004,
        0.0735,
        0.0126,
        0.117,
        0.2046,
        0.117,
        0.000985945,
        0.0933,
        0.0696,
        0.2268,
        0.0018,
        0.0917,
        0.15855,
        0.0639,
        0.00195,
        0.101280911,
        0.1425,
        0.0615,
        0.003025,
        0.0858,
        0.214175,
        0.0567,
        0.00065,
        0.1057,
        0.212875,
        0.075725,
        0.0081,
        0.0717,
        0.192,
        0.0741,
        3e-04
      )
    ),
    row.names = c(
      2L,
      3L,
      4L,
      6L,
      7L,
      8L,
      10L,
      11L,
      12L,
      14L,
      15L,
      16L,
      18L,
      19L,
      20L,
      22L,
      23L,
      24L,
      26L,
      27L,
      28L,
      30L,
      31L,
      32L,
      34L,
      35L,
      36L,
      38L,
      39L,
      40L,
      42L,
      43L,
      44L,
      46L,
      47L,
      48L,
      50L,
      51L,
      52L,
      54L,
      55L,
      56L,
      58L,
      59L,
      60L,
      62L,
      63L,
      64L,
      66L,
      67L,
      68L,
      70L,
      71L,
      72L,
      74L,
      75L,
      76L,
      78L,
      79L,
      80L,
      82L,
      83L,
      84L,
      86L,
      87L,
      88L,
      90L,
      91L,
      92L,
      94L,
      95L,
      96L,
      121L,
      122L,
      123L,
      124L,
      125L,
      126L,
      127L,
      128L,
      129L,
      130L,
      131L,
      132L,
      133L,
      134L,
      135L,
      136L,
      137L,
      138L,
      139L,
      140L,
      141L,
      142L,
      143L,
      144L,
      145L,
      146L,
      147L,
      148L,
      149L,
      150L,
      151L,
      152L,
      153L,
      154L,
      155L,
      156L,
      157L,
      158L,
      159L,
      160L,
      161L,
      162L,
      163L,
      164L,
      165L,
      166L,
      167L,
      168L,
      169L,
      170L,
      171L,
      172L,
      173L,
      174L,
      175L,
      176L,
      177L,
      178L,
      179L,
      180L,
      181L,
      182L,
      183L,
      184L,
      185L,
      186L,
      187L,
      188L,
      189L,
      190L,
      191L,
      192L
    ),
    class = "data.frame"
  )
Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ please, explain which differences are you interested in. Interaction, marginal effect, etc. ? $\endgroup$
    – utobi
    Commented Mar 27, 2023 at 13:35
  • 1
    $\begingroup$ Of course, I am interested in the interaction (if any) that soil_type, treatment, and days [since application] have on cl_conc in my experiment. I don't believe I am interested in marginal effects. Hope that helps, but I can try to elaborate more if need be $\endgroup$
    – mbelanger081
    Commented Mar 27, 2023 at 13:38
  • 1
    $\begingroup$ See also: stats.stackexchange.com/questions/577873/is-there-a-non-parametric-form-of-a-3-way-anova/610996 $\endgroup$
    – Sal Mangiafico
    Commented Mar 28, 2023 at 12:58

1 Answer 1

Reset to default
10
$\begingroup$

The tests you cited are not appropriate due to the presence of repeated measures. The common way to deal with repeated measures is via mixed-effects linear models.

I'm considering here the most general model, borrowing from one of your earlier posts.

> leach_lme <- lme(fixed = cl_conc ~ soil_type*treatment*days,
+                  random =~1|core_id, data = leach2, 
+                  method = "ML")
> anova(leach_lme)
                         numDF denDF  F-value p-value
(Intercept)                  1    80 677.4590  <.0001
soil_type                    3    32   6.1510  0.0020
treatment                    3    32 109.4603  <.0001
days                         1    80  17.3933  0.0001
soil_type:treatment          9    32   2.2588  0.0436
soil_type:days               3    80   3.1330  0.0301
treatment:days               3    80   1.0310  0.3834
soil_type:treatment:days     9    80   3.9676  0.0003

As you can see, the three-way interaction is significant, and so is the two-way interaction soil_type:days, etc. Linear mixed-effects is a parametric model, so as usual in the context of a linear model, one needs to check the residuals are well-behaving.

As per request in the comments, here is a quick residual check.

plot(leach_lme)

enter image description here

The message here is that residuals may be heteroscedastic. Now let's log-transforming the response

leach_lme2 <- lme(fixed = log(cl_conc) ~ soil_type*treatment*days,
                 random =~1|core_id, data = leach2, 
                 method = "ML")
plot(leach_lme2)

enter image description here

A part of a single observation which appears to be far from the bulk of the data, the residuals look fine to me, i.e. homoscedastic. The QQ-plot as well (not shown here but you can plot it yourself via qqnorm(leach_lme2)) doesn't seem that bad.

P.S. In this answer, I treated days as a numerical variable. To treat it as a factor, as you seem to be interested in (thanks Sal Mangiacifo for pointing it out), use

leach2$days <- factor(leach2$days)

and redo the analyses. The output will be slightly different compared to the one shown above; there will be two additional parameters to be estimated.

Improve this answer
$\endgroup$
8
  • $\begingroup$ Thanks for the advice. It has been very helpful. I ran the code, and the residuals were not behaving correctly and do not get better with data transformation. It's my understanding I would need a nonparametric test, but the Friedman Test does not work on my data. I am unsure what is appropriate for this data as I have been having trouble finding a test that works for a while. Any advice for me to try going forward? $\endgroup$
    – mbelanger081
    Commented Mar 27, 2023 at 14:05
  • 1
    $\begingroup$ From what I can see, the residuals here are a bit heteroscedastic. here are some quick fixes that come to my mind right now: (1) Log-transform the response; I've checked and seem to fix the issue. (2) include further covariates if you have any. (3) work with the original model but compute the p-value via bootstrap; this would solve the heteroscedasticity issue. $\endgroup$
    – utobi
    Commented Mar 27, 2023 at 14:12
  • 1
    $\begingroup$ @mbelanger081 see my updated answer. As you'll see, log-transformation does fix the heteroscedasticity issue. There is, however, a problem with a single observation having an exceptionally high residual. If this outlier bothers you, try fitting lme via a robust method; check the robustlmm package in R. $\endgroup$
    – utobi
    Commented Mar 27, 2023 at 15:09
  • 1
    $\begingroup$ I would use normality tests for checking residuals with great care; see here $\endgroup$
    – utobi
    Commented Mar 27, 2023 at 15:12
  • 1
    $\begingroup$ +1, although it sounds like OP wants to treat Days as a factor variable. (It has only 3 levels). $\endgroup$
    – Sal Mangiafico
    Commented Mar 28, 2023 at 12:15

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