# ANOVA Tests in Julia?

**URL:** <https://discourse.julialang.org/t/anova-tests-in-julia/1421>\
**Category:** Statistics\
**Created:** [January 11, 2017, 5:45pm UTC](https://discourse.julialang.org/t/anova-tests-in-julia/1421 "2017-01-11T17:45:15Z")\
**Posts on this page:** 17\
**Page:** 4

<div class="post-metadata">

**Author:** ![BioTurboNick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bioturbonick/32/6380_2.png) [@BioTurboNick](https://discourse.julialang.org/u/BioTurboNick)\
**Post date:** [April 12, 2021, 4:06pm UTC](https://discourse.julialang.org/t/anova-tests-in-julia/1421/61 "2021-04-12T16:06:06Z")

</div>

Since people keep discovering this post, I think the only two packages worth talking about here are now:

[marcpabst](https://github.com/marcpabst)/ **[ANOVA.jl](https://github.com/marcpabst/ANOVA.jl)** requires DataFrames, external setup and fitting of linear model (GLM.jl), and I _think_ is only fixed factors, with “type” argument referring to variant ways to calculate the sums of squares. I don’t know if there’s any way to do nested factors with linear models. I’m assuming robust to unbalanced data due to reliance on regression?, updated 8 months ago

and mine:

[BioTurboNick / **SimpleANOVA.jl**](https://github.com/BioTurboNick/SimpleANOVA.jl) which does a non-GLM-requiring calculation, does random effects, nesting and repeated measures, can use DataFrames or vectors or multidimensional arrays, but requires balanced data and only does one type of sum of squares. Also can do contrasts and omega-square effect size.

---

<div class="post-metadata">

**Author:** ![babaq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/babaq/32/891_2.png) [@babaq](https://discourse.julialang.org/u/babaq)\
**Post date:** [December 31, 2021, 3:43am UTC](https://discourse.julialang.org/t/anova-tests-in-julia/1421/62 "2021-12-31T03:43:09Z")

</div>

I’ve used the [ANOVA](https://github.com/marcpabst/ANOVA.jl) successfully, but the package seems not maintained anymore, and incompatible with some of the packages.

There is another package [MixedAnova](https://github.com/yufongpeng/MixedAnova.jl), so i am curious that if anyone tried it, and the overall experience of using it.

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [July 23, 2022, 9:43am UTC](https://discourse.julialang.org/t/anova-tests-in-julia/1421/63 "2022-07-23T09:43:34Z")

</div>

Is really no any actual package for ANOVA?

---

<div class="post-metadata">

**Author:** ![maxkapur](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maxkapur/32/21208_2.png) [@maxkapur](https://discourse.julialang.org/u/maxkapur)\
**Post date:** [July 23, 2022, 1:59pm UTC](https://discourse.julialang.org/t/anova-tests-in-julia/1421/64 "2022-07-23T13:59:04Z")

</div>

HypothesisTests.jl has one-way ANOVA now:

[https://juliastats.org/HypothesisTests.jl/stable/parametric/#One-way-ANOVA-Test](https://juliastats.org/HypothesisTests.jl/stable/parametric/#One-way-ANOVA-Test)

And GLM.jl’s `ftest` model claims to be usable for ANOVA:

[https://juliastats.org/GLM.jl/stable/api/#GLM.ftest](https://juliastats.org/GLM.jl/stable/api/#GLM.ftest)

I think these are both basic cases, however, and I will leave it to those who know more about statistics to speak to the availability of more general tools.

---

<div class="post-metadata">

**Author:** ![EvoArt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evoart/32/25357_2.png) [@EvoArt](https://discourse.julialang.org/u/EvoArt)\
**Post date:** [July 23, 2022, 9:54pm UTC](https://discourse.julialang.org/t/anova-tests-in-julia/1421/65 "2022-07-23T21:54:49Z")

</div>

> **[GitHub - yufongpeng/AnovaBase.jl: A simple framework for ANOVA on various...](https://github.com/yufongpeng/AnovaBase.jl)**
>
> A simple framework for ANOVA on various types of Julia statistical models - GitHub - yufongpeng/AnovaBase.jl: A simple framework for ANOVA on various types of Julia statistical models

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [July 24, 2022, 8:30pm UTC](https://discourse.julialang.org/t/anova-tests-in-julia/1421/66 "2022-07-24T20:30:50Z")

</div>

It seems work. But MixedAnova.jl should be intalled. I think there is ssue with package installation.

---

<div class="post-metadata">

**Author:** ![BioTurboNick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bioturbonick/32/6380_2.png) [@BioTurboNick](https://discourse.julialang.org/u/BioTurboNick)\
**Post date:** [August 10, 2022, 1:26pm UTC](https://discourse.julialang.org/t/anova-tests-in-julia/1421/67 "2022-08-10T13:26:56Z")

</div>

Not sure how you define “actual package” but mine does many different variants of ANOVA: [https://github.com/BioTurboNick/SimpleANOVA.jl](https://github.com/BioTurboNick/SimpleANOVA.jl)

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [August 10, 2022, 2:23pm UTC](https://discourse.julialang.org/t/anova-tests-in-julia/1421/68 "2022-08-10T14:23:39Z")

</div>

@BioTurboNick having just stumbled on this thread because it popped up due to recent comments, I see you implemented this but your background is in biology and not mathematics etc. My biggest concern is that calculations of sums of squares are numerically unstable and it’s easily possible to get wrong answers due to numerical computation issues. How much were you aware of numerical stability issues when coding your package?

---

<div class="post-metadata">

**Author:** ![BioTurboNick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bioturbonick/32/6380_2.png) [@BioTurboNick](https://discourse.julialang.org/u/BioTurboNick)\
**Post date:** [August 10, 2022, 2:33pm UTC](https://discourse.julialang.org/t/anova-tests-in-julia/1421/69 "2022-08-10T14:33:03Z")

</div>

At the time I wrote it, nearly zero awareness of numerical stability I must admit.

I have recently had a headlong crash-course in dealing with numerical stability in another context, so it would be good for me to revisit the package with that in mind.

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [August 10, 2022, 9:42pm UTC](https://discourse.julialang.org/t/anova-tests-in-julia/1421/70 "2022-08-10T21:42:11Z")

</div>

Hi! Good point! But what about `var` from StatsBase? This is really just sum of squares… Do we need adjust StatsBase? What kind of algorithm could be prefered? (Welford? Weighted incremental? )

```julia
##### General central moment
function _moment2(v::RealArray, m::Real; corrected=false)
    n = length(v)
    s = 0.0
    for i = 1:n
        @inbounds z = v[i] - m
        s += z * z
    end
    varcorrection(n, corrected) * s
end

```

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [August 11, 2022, 4:45am UTC](https://discourse.julialang.org/t/anova-tests-in-julia/1421/72 "2022-08-11T04:45:31Z")

</div>

> [@PharmCat](#):
>
> Welford

Seems like a reasonable choice.

---

<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:** [August 11, 2022, 6:50am UTC](https://discourse.julialang.org/t/anova-tests-in-julia/1421/73 "2022-08-11T06:50:25Z")

</div>

> [@BioTurboNick](#):
>
> so it would be good for me to revisit the package with that in mind.

Let me know, if I can contribute in some way, when you revisit your package. [My implementation of R’s aov](https://discourse.julialang.org/t/repeated-measures-anova/60765/12) shouldn’t suffer from numerical instability.

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [August 11, 2022, 10:29am UTC](https://discourse.julialang.org/t/anova-tests-in-julia/1421/74 "2022-08-11T10:29:55Z")

</div>

Do we really have problems with this?

Seems `var` work… Can you kindly provide more robust test?

```julia
using StatsBase

v = rand(1000000)
v2 = v .+ 1000_000_000
v3 = append!(copy(v2), v2)
tv = rand(8)

function var2(v)
    M = v[1]
    S = 0.0
    N = length(v)
    for k = 2:N
    x = v[k]
    Mn = M
    M += (x - M) / k
    S += (x - M) * (x - Mn)
    end
  return S / (N - 1)
end

function var3(v)
    M = v[1]
    S = 0.0
    N = length(v)
    for k = 2:N
        x = v[k]
        Mn = M
        M = Mn + (x - Mn) / k
        S = S + ((x - Mn) * (x - M) - S) / k
    end
  return S * N / (N - 1)
end

println("check: ", isapprox(var(tv), var2(tv)), " ", isapprox(var(tv), var3(tv)))
println("var: ", var(v) - var(v2), " ", var(v) - var(v3))
println("var2: ", var(v) - var2(v2), " ", var(v) - var2(v3))
println("var3: ", var(v) - var3(v2), " ", var(v) - var3(v3))

```

```julia
check: true true
var: 2.2700591406632498e-12 4.16894628540998e-8
var2: 1.3641951304710354e-8 6.139953312445101e-8
var3: 1.3641966167821096e-8 6.1399576617438e-8

```

---

<div class="post-metadata">

**Author:** ![BioTurboNick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bioturbonick/32/6380_2.png) [@BioTurboNick](https://discourse.julialang.org/u/BioTurboNick)\
**Post date:** [August 11, 2022, 2:02pm UTC](https://discourse.julialang.org/t/anova-tests-in-julia/1421/75 "2022-08-11T14:02:26Z")

</div>

I appreciate the offer! I did some checking and made PR.

[https://github.com/BioTurboNick/SimpleANOVA.jl/pull/23](https://github.com/BioTurboNick/SimpleANOVA.jl/pull/23)

It would be great if you could see if what I did makes any sense, or if I’m missing other sources of instability.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [August 11, 2022, 2:11pm UTC](https://discourse.julialang.org/t/anova-tests-in-julia/1421/76 "2022-08-11T14:11:12Z")

</div>

> [@PharmCat](#):
>
> provide more robust test?

Just some thoughts… Haven’t tested yet.

Try sampling uniform on 0 to max float value. Try sampling from a t dist with 2 degrees of freedom. Try sampling from a t dist with 2 dof and center at half max float. Try sampling from a distribution where variance is something like 4x machine epsilon for center of 1/16 max float…

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [August 11, 2022, 5:58pm UTC](https://discourse.julialang.org/t/anova-tests-in-julia/1421/77 "2022-08-11T17:58:11Z")

</div>

I made some experiments, and in the worst case when variable is `exp.(rand(Uniform(1,log(prevfloat(typemax(Float64)))), 10000000))` we have `Inf` for `var` and `var2` ; NaN for `var3`.

> Try sampling from a distribution where variance is something like 4x machine epsilon for center of 1/16 max float

For this, I have `Inf` for `var` and `0.0` for `var2` and `var3`

So for `var2` and `var3` - these are appropriate results, but `Inf` for `var` is discussible.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [August 11, 2022, 7:23pm UTC](https://discourse.julialang.org/t/anova-tests-in-julia/1421/78 "2022-08-11T19:23:02Z")

</div>

> [@PharmCat](#):
>
> `exp.(rand(Uniform(1,log(prevfloat(typemax(Float64)))), 10000000))` we have `Inf` for `var` and `var2` ; NaN for `var3`.

Since var is giving things of scale x^2, try something like:

```julia
var(exp.(rand(Uniform(1,log(prevfloat(typemax(Float64)))/2.1), 10000000)))

```

[Previous page](https://discourse.julialang.org/t/anova-tests-in-julia/1421.md?page=3)