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

**URL:** https://discourse.julialang.org/t/equivalent-of-rs-aov-y-b-error-a-b-data-d/56629
**Category:** General Usage
**Tags:** question, statistics
**Created:** [March 6, 2021, 8:36pm UTC](https://discourse.julialang.org/t/equivalent-of-rs-aov-y-b-error-a-b-data-d/56629 "2021-03-06T20:36:28Z")
**Posts on this page:** 2
**Page:** 1

<div class="post-metadata">

### Author: ![jbrea](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jbrea/32/3879_2.png) [@jbrea](https://discourse.julialang.org/u/jbrea)
#### Post date: [March 6, 2021, 8:36pm UTC](https://discourse.julialang.org/t/equivalent-of-rs-aov-y-b-error-a-b-data-d/56629/1 "2021-03-06T20:36:28Z")

</div>

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
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  

```

---

<div class="post-metadata">

### Author: ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)
#### Post date: [March 7, 2021, 10:22am UTC](https://discourse.julialang.org/t/equivalent-of-rs-aov-y-b-error-a-b-data-d/56629/2 "2021-03-07T10:22:19Z")

</div>

See: [ANOVA Tests in Julia? - #58 by olangsrud](https://discourse.julialang.org/t/anova-tests-in-julia/1421/58)

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.