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