# Unbalanced One-way ANOVA

**URL:** <https://discourse.julialang.org/t/unbalanced-one-way-anova/69489>\
**Category:** Statistics\
**Tags:** question\
**Created:** [October 9, 2021, 10:25pm UTC](https://discourse.julialang.org/t/unbalanced-one-way-anova/69489 "2021-10-09T22:25:40Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![AerospaceMark](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aerospacemark/32/29823_2.png) [@AerospaceMark](https://discourse.julialang.org/u/AerospaceMark)\
**Post date:** [October 9, 2021, 10:25pm UTC](https://discourse.julialang.org/t/unbalanced-one-way-anova/69489/1 "2021-10-09T22:25:40Z")

</div>

I’m a graduate student seeking to learn how to do statistics in Julia. From my internet searching and working with the SimpleANOVA package, I’ve been able to perform ANOVA tests on data where each group has the same number of samples. Now I’m trying to learn how to do ANOVA tests when each group has a different number of samples. This does not appear to be supported yet in the SimpleANOVA package, is there a good way to do this type of test in another package?

---

<div class="post-metadata">

**Author:** ![viraltux](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/viraltux/32/15236_2.png) [@viraltux](https://discourse.julialang.org/u/viraltux)\
**Post date:** [October 9, 2021, 11:14pm UTC](https://discourse.julialang.org/t/unbalanced-one-way-anova/69489/2 "2021-10-09T23:14:29Z")

</div>

Welcome @AerospaceMark!

Not sure about this specific case but I think overall it’d be good for you to have a look to this book and resources:

> **[Statistics with Julia Book](https://statisticswithjulia.org/index.html)**
>
> Statistics with Julia: Fundamentals for Data Science, Machine Learning and Artificial Intelligence

> **[GitHub - h-Klok/StatsWithJuliaBook](https://github.com/h-Klok/StatsWithJuliaBook)**
>
> Contribute to h-Klok/StatsWithJuliaBook development by creating an account on GitHub.

---

<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:** [October 11, 2021, 8:25am UTC](https://discourse.julialang.org/t/unbalanced-one-way-anova/69489/3 "2021-10-11T08:25:49Z")

</div>

There is also [ANOVA](https://github.com/marcpabst/ANOVA.jl) and here is some (messy but tested) code that implements R’s `aov` function in julia (I haven’t put it in a package yet). Use it with `aov(formula, data)`, where `formula` is a `StatsModel` formula like e.g. `formula = @formula(response ~ group * treatment)` and `data` is a `DataFrame`.

```julia
using StatsModels, DataFrames, LinearAlgebra, Distributions
toterm(s1::Term, s2::Term) = InteractionTerm((s1, s2))
toterm(s1::InteractionTerm, s2::Term) = InteractionTerm((s1.terms..., s2))
toterm(s1::Term, s2::InteractionTerm) = InteractionTerm((s1, s2.terms...))
toterm(x::Tuple, s2::Term) = map(t -> toterm(t, s2), x)
toterm(s1::Term, x::Tuple) = map(t -> toterm(s1, t), x)
function predictor_product(predictors)
    init, rest = Iterators.peel(predictors)
    reduce((x, y) -> x + y + toterm(x, y), term.(rest), init = term(init))
end
function errorterm(ef, nested)
    terms = term(1) + term(ef)
    if length(nested) > 0
        nestedterms = predictor_product(nested)
        terms += toterm(term(ef), nestedterms)
    end
    terms
end

name(t::AbstractTerm) = "$(t.sym)"
name(t::InteractionTerm) = join(name.(t.terms), " & ")
name(t::InterceptTerm) = "Intercept"
pivoted_asgn(assign, qrp) = assign[qrp.p[1:rank(qrp.R)]]
function mqr(f::FormulaTerm, data; contrasts = Dict{Symbol, Any}())
    mf = ModelFrame(f, data; contrasts)
    mm = ModelMatrix(mf)
    qrp = qr(mm.m, Val(true))
    (qrp = qrp, mf = mf, mm = mm)
end
function aov(f::FormulaTerm, data; return_locals = false)
    qrp, mf, mm = mqr(f, data)
    aov(qrp, pivoted_asgn(mm.assign, qrp), mf.f.rhs.terms,
        getproperty(data, f.lhs.sym); return_locals)
end
aov(l, y) = aov(l..., y)
function aov(qrp, asgn, terms, y; return_locals = false)
    effects = qrp.Q' * y
    r = rank(qrp.R)
    SS = Float64[]
    names = String[]
    DF = Int[]
    for i in union(asgn)
        n = name(terms[i])
        n == "Intercept" && continue
        idxs = findall(asgn .== i)
        push!(SS, sum(abs2, @view(effects[idxs])))
        push!(names, n)
        push!(DF, length(idxs))
    end
    push!(SS, sum(abs2, @view(effects[r+1:end])))
    push!(names, "Residuals")
    push!(DF, length(effects) - r)
    MSS = SS ./ DF
    F = MSS ./ MSS[end]
    F[end] = NaN
    DFres = DF[end]
    result = DataFrame(names = names, DF = DF, SS = SS, MSS = MSS, F = F,
                       p = @.(ccdf(FDist(DF, DFres), F)))
    if return_locals
        result, (qrp = qrp, asgn = asgn, terms = terms)
    else
        result
    end
end
function aoverror(f::FormulaTerm, e, nested::T, data; return_locals = false) where T
    nested = T <: Symbol ? [nested] : nested
    ef = errorterm(e, nested)
    efit = mqr(f.lhs ~ ef, data,
        contrasts = Dict(e => HelmertCoding(),
                         [t => HelmertCoding() for t in nested]...))
    efit.mm.assign = efit.mm.assign[efit.qrp.p[1:rank(efit.qrp.R)]]
    mf = ModelFrame(f, data)
    mm = ModelMatrix(mf)
    ys = efit.qrp.Q' * getproperty(data, f.lhs.sym)
    append!(efit.mm.assign,
            fill(maximum(efit.mm.assign) + 1, length(ys) - length(efit.mm.assign)))
    Xs = efit.qrp.Q' * mm.m
    rows = [findall(efit.mm.assign .== i)
            for i in Iterators.drop(sort!(union(efit.mm.assign)), 1)]
    result = [inner(Xs, ys, row, mm.assign, mf.f.rhs.terms; return_locals)
              for row in rows]
    if return_locals
        first.(result), (Q = efit.qrp.Q,
                         locals_t = [(rest..., row)
                                     for (row, rest) in zip(rows, last.(result))])
    else
        result
    end
end
function aoverror(l, y)
    ys = l.Q' * y
    [aov(qrp, asgn, terms, ys[rows]) for (qrp, asgn, terms, rows) in l.locals_t]
end
function inner(X, y, rows, asgn, terms; return_locals = false)
    T = X[rows, :]
    cols = findall(reshape(sum(abs2, T, dims = 1), :) .> 1e-5)
    qrp = qr(T[:, cols], Val(true))
    aov(qrp, pivoted_asgn(asgn[cols], qrp), terms, y[rows]; return_locals)
end

```

---

<div class="post-metadata">

**Author:** ![AerospaceMark](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aerospacemark/32/29823_2.png) [@AerospaceMark](https://discourse.julialang.org/u/AerospaceMark)\
**Post date:** [October 12, 2021, 12:46pm UTC](https://discourse.julialang.org/t/unbalanced-one-way-anova/69489/4 "2021-10-12T12:46:25Z")

</div>

Awesome, thanks!

---

<div class="post-metadata">

**Author:** ![AerospaceMark](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aerospacemark/32/29823_2.png) [@AerospaceMark](https://discourse.julialang.org/u/AerospaceMark)\
**Post date:** [October 12, 2021, 12:46pm UTC](https://discourse.julialang.org/t/unbalanced-one-way-anova/69489/5 "2021-10-12T12:46:55Z")

</div>

This looks like great code, thanks for sharing!