# Equivalent of Rs `aov(Y ~ B + Error(A/B), data=d)`

Is there a julia equivalent of Rโs way to compute classical anovas with random effects?

For example, is there an easy way to obtain the following result without calling R?

``````julia> using DataFrames, RCall

julia> df = DataFrame(group = [fill("A", 24); fill("B", 24); fill("C", 24)],
order = repeat([0, 0, 1, 1], 18),
id = vcat([repeat(repeat(4i+1:4(i+1), inner = 2), 3) for i in 0:2]...),
trial = repeat(vcat([fill(i, 8) for i in 1:3]...), 3),
treatment = repeat(["X", "Y"], 36),
data = rand(72)
)
72ร6 DataFrame
Row โ group   order  id     trial  treatment  data
โ String  Int64  Int64  Int64  String     Float64
โโโโโโผโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ
1 โ A           0      1      1  X          0.207635
2 โ A           0      1      1  Y          0.376889
3 โ A           1      2      1  X          0.172655
4 โ A           1      2      1  Y          0.711968
5 โ A           0      3      1  X          0.451885
6 โ A           0      3      1  Y          0.344559
7 โ A           1      4      1  X          0.909206
8 โ A           1      4      1  Y          0.0705808
9 โ A           0      1      2  X          0.924089
10 โ A           0      1      2  Y          0.207873
11 โ A           1      2      2  X          0.312672
โฎ  โ   โฎ       โฎ      โฎ      โฎ        โฎ          โฎ
63 โ C           1     12      2  X          0.0685433
64 โ C           1     12      2  Y          0.732906
65 โ C           0      9      3  X          0.492607
66 โ C           0      9      3  Y          0.851092
67 โ C           1     10      3  X          0.162958
68 โ C           1     10      3  Y          0.190779
69 โ C           0     11      3  X          0.588648
70 โ C           0     11      3  Y          0.911885
71 โ C           1     12      3  X          0.169196
72 โ C           1     12      3  Y          0.87193
51 rows omitted

julia> categorical!(df, [:group, :order, :id, :trial, :treatment])
72ร6 DataFrame
Row โ group  order  id    trial  treatment  data
โ Catโฆ   Catโฆ   Catโฆ  Catโฆ   Catโฆ       Float64
โโโโโโผโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ
1 โ A      0      1     1      X          0.207635
2 โ A      0      1     1      Y          0.376889
3 โ A      1      2     1      X          0.172655
4 โ A      1      2     1      Y          0.711968
5 โ A      0      3     1      X          0.451885
6 โ A      0      3     1      Y          0.344559
7 โ A      1      4     1      X          0.909206
8 โ A      1      4     1      Y          0.0705808
9 โ A      0      1     2      X          0.924089
10 โ A      0      1     2      Y          0.207873
11 โ A      1      2     2      X          0.312672
โฎ  โ   โฎ      โฎ     โฎ      โฎ        โฎ          โฎ
63 โ C      1      12    2      X          0.0685433
64 โ C      1      12    2      Y          0.732906
65 โ C      0      9     3      X          0.492607
66 โ C      0      9     3      Y          0.851092
67 โ C      1      10    3      X          0.162958
68 โ C      1      10    3      Y          0.190779
69 โ C      0      11    3      X          0.588648
70 โ C      0      11    3      Y          0.911885
71 โ C      1      12    3      X          0.169196
72 โ C      1      12    3      Y          0.87193
51 rows omitted

julia> @rput df
72ร6 DataFrame
Row โ group  order  id    trial  treatment  data
โ Catโฆ   Catโฆ   Catโฆ  Catโฆ   Catโฆ       Float64
โโโโโโผโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ
1 โ A      0      1     1      X          0.207635
2 โ A      0      1     1      Y          0.376889
3 โ A      1      2     1      X          0.172655
4 โ A      1      2     1      Y          0.711968
5 โ A      0      3     1      X          0.451885
6 โ A      0      3     1      Y          0.344559
7 โ A      1      4     1      X          0.909206
8 โ A      1      4     1      Y          0.0705808
9 โ A      0      1     2      X          0.924089
10 โ A      0      1     2      Y          0.207873
11 โ A      1      2     2      X          0.312672
โฎ  โ   โฎ      โฎ     โฎ      โฎ        โฎ          โฎ
63 โ C      1      12    2      X          0.0685433
64 โ C      1      12    2      Y          0.732906
65 โ C      0      9     3      X          0.492607
66 โ C      0      9     3      Y          0.851092
67 โ C      1      10    3      X          0.162958
68 โ C      1      10    3      Y          0.190779
69 โ C      0      11    3      X          0.588648
70 โ C      0      11    3      Y          0.911885
71 โ C      1      12    3      X          0.169196
72 โ C      1      12    3      Y          0.87193
51 rows omitted

julia> R"""
summary(aov(data ~ group*trial*treatment + Error(id/(trial*treatment)), df))
"""
RObject{VecSxp}

Error: id
Df Sum Sq Mean Sq F value Pr(>F)
group      2 0.0273 0.01363   0.312   0.74
Residuals  9 0.3932 0.04368

Error: id:trial
Df Sum Sq Mean Sq F value Pr(>F)
trial        2 0.0944 0.04718   0.373  0.694
group:trial  4 0.0733 0.01833   0.145  0.963
Residuals   18 2.2776 0.12653

Error: id:treatment
Df Sum Sq Mean Sq F value Pr(>F)
treatment        1 0.0744 0.07444   1.776 0.2154
group:treatment  2 0.4652 0.23259   5.549 0.0269 *
Residuals        9 0.3773 0.04192
---
Signif. codes:  0 โ***โ 0.001 โ**โ 0.01 โ*โ 0.05 โ.โ 0.1 โ โ 1

Error: id:trial:treatment
Df Sum Sq Mean Sq F value Pr(>F)
trial:treatment        2 0.0986 0.04932   0.555  0.584
group:trial:treatment  4 0.7584 0.18960   2.134  0.118
Residuals             18 1.5995 0.08886
``````

Itโs a bit out of date but you could have a look as to which of the packages mentioned in the summary towards the end of the thread are still actively maintained & fit your use case.

1 Like