OneWayAnova.java

/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.apache.commons.statistics.inference;

import java.util.Collection;
import java.util.Iterator;
import org.apache.commons.numbers.core.Sum;
import org.apache.commons.statistics.distribution.FDistribution;

/**
 * Implements one-way ANOVA (analysis of variance) statistics.
 *
 * <p>Tests for differences between two or more categories of univariate data
 * (for example, the body mass index of accountants, lawyers, doctors and
 * computer programmers). When two categories are given, this is equivalent to
 * the {@link TTest}.
 *
 * <p>This implementation computes the F statistic using the definitional formula:
 *
 * <p>\[ F = \frac{\text{between-group variability}}{\text{within-group variability}} \]
 *
 * @see <a href="https://en.wikipedia.org/wiki/Analysis_of_variance">Analysis of variance (Wikipedia)</a>
 * @see <a href="https://en.wikipedia.org/wiki/F-test#Multiple-comparison_ANOVA_problems">
 * Multiple-comparison ANOVA problems (Wikipedia)</a>
 * @see <a href="https://www.biostathandbook.com/onewayanova.html">
 * McDonald, J.H. 2014. Handbook of Biological Statistics (3rd ed.). Sparky House Publishing, Baltimore, Maryland.
 * One-way anova. pp 145-156.</a>
 * @since 1.1
 */
public final class OneWayAnova {
    /** Default instance. */
    private static final OneWayAnova DEFAULT = new OneWayAnova();

    /**
     * Result for the one-way ANOVA.
     *
     * <p>This class is immutable.
     *
     * @since 1.1
     */
    public static final class Result extends BaseSignificanceResult {
        /** Degrees of freedom in numerator (between groups). */
        private final int dfbg;
        /** Degrees of freedom in denominator (within groups). */
        private final long dfwg;
        /** Mean square between groups. */
        private final double msbg;
        /** Mean square within groups. */
        private final double mswg;
        /** nO value used to partition the variance. */
        private final double nO;

        /**
         * @param dfbg Degrees of freedom in numerator (between groups).
         * @param dfwg Degrees of freedom in denominator (within groups).
         * @param msbg Mean square between groups.
         * @param mswg Mean square within groups.
         * @param nO Factor for partitioning the variance.
         * @param f F statistic
         * @param p P-value.
         */
        Result(int dfbg, long dfwg, double msbg, double mswg, double nO, double f, double p) {
            super(f, p);
            this.dfbg = dfbg;
            this.dfwg = dfwg;
            this.msbg = msbg;
            this.mswg = mswg;
            this.nO = nO;
        }

        /**
         * Gets the degrees of freedom in the numerator (between groups).
         *
         * @return degrees of freedom between groups
         */
        int getDFBG() {
            return dfbg;
        }

        /**
         * Gets the degrees of freedom in the denominator (within groups).
         *
         * @return degrees of freedom within groups
         */
        long getDFWG() {
            return dfwg;
        }

        /**
         * Gets the mean square between groups.
         *
         * @return mean square between groups
         */
        public double getMSBG() {
            return msbg;
        }

        /**
         * Gets the mean square within groups.
         *
         * @return mean square within groups
         */
        public double getMSWG() {
            return mswg;
        }

        /**
         * Gets the variance component between groups.
         *
         * <p>The value is a partitioning of the variance.
         * It is the complement of {@link #getVCWG()}.
         *
         * <p>Partitioning the variance applies only to a model II
         * (random effects) one-way anova. This applies when the
         * groups are random samples from a larger set of groups;
         * partitioning the variance allows comparison of the
         * variation between groups to the variation within groups.
         *
         * <p>If the {@linkplain #getMSBG() MSBG} is less than the
         * {@linkplain #getMSWG() MSWG} this returns 0. Otherwise this
         * creates an estimate of the added variance component
         * between groups as:
         *
         * <p>\[ \text{between-group variance} = A = (\text{MS}_{\text{bg}} - \text{MS}_{\text{wg}}) / n_o \]
         *
         * <p>where \( n_o \) is a number close to, but usually less than,
         * the arithmetic mean of the sample size \(n_i\) of each
         * of the \( a \) groups:
         *
         * <p>\[ n_o = \frac{1}{a-1} \left( \sum_i{n_i} - \frac{\sum_i{n_i^2}}{\sum_i{n_i}} \right) \]
         *
         * <p>The added variance component among groups \( A \) is expressed
         * as a fraction of the total variance components \( A + B \) where
         * \( B \) is the {@linkplain #getMSWG() MSWG}.
         *
         * @return variance component between groups (in [0, 1]).
         */
        public double getVCBG() {
            if (msbg <= mswg) {
                return 0;
            }
            // a is an estimate of the between-group variance
            final double a = (msbg - mswg) / nO;
            final double b = mswg;
            return a / (a + b);
        }

        /**
         * Gets the variance component within groups.
         *
         * <p>The value is a partitioning of the variance.
         * It is the complement of {@link #getVCBG()}. See
         * that method for details.
         *
         * @return variance component within groups (in [0, 1]).
         */
        public double getVCWG() {
            if (msbg <= mswg) {
                return 1;
            }
            final double a = (msbg - mswg) / nO;
            final double b = mswg;
            return b / (a + b);
        }
    }

    /** Private constructor. */
    private OneWayAnova() {
        // Do nothing
    }

    /**
     * Return an instance using the default options.
     *
     * @return default instance
     */
    public static OneWayAnova withDefaults() {
        return DEFAULT;
    }

    /**
     * Computes the F statistic for an ANOVA test for a collection of category data,
     * evaluating the null hypothesis that there is no difference among the means of
     * the data categories.
     *
     * <p>Special cases:
     * <ul>
     * <li>If the value in each category is the same (no variance within groups) but different
     * between groups, the f-value is {@linkplain Double#POSITIVE_INFINITY infinity}.
     * <li>If the value in every group is the same (no variance within or between groups),
     * the f-value is {@link Double#NaN NaN}.
     * </ul>
     *
     * @param data Category summary data.
     * @return F statistic
     * @throws IllegalArgumentException if the number of categories is less than
     * two; a contained category does not have at least one value; or all
     * categories have only one value (zero degrees of freedom within groups)
     */
    public double statistic(Collection<double[]> data) {
        final double[] f = new double[1];
        aov(data, f);
        return f[0];
    }

    /**
     * Performs an ANOVA test for a collection of category data,
     * evaluating the null hypothesis that there is no difference among the means of
     * the data categories.
     *
     * <p>Special cases:
     * <ul>
     * <li>If the value in each category is the same (no variance within groups) but different
     * between groups, the f-value is {@linkplain Double#POSITIVE_INFINITY infinity} and the p-value is zero.
     * <li>If the value in every group is the same (no variance within or between groups),
     * the f-value and p-value are {@link Double#NaN NaN}.
     * </ul>
     *
     * @param data Category summary data.
     * @return test result
     * @throws IllegalArgumentException if the number of categories is less than
     * two; a contained category does not have at least one value; or all
     * categories have only one value (zero degrees of freedom within groups)
     */
    public Result test(Collection<double[]> data) {
        return aov(data, null);
    }

    /**
     * Performs an ANOVA test for a collection of category data, evaluating the null
     * hypothesis that there is no difference among the means of the data categories.
     *
     * <p>This is a utility method to allow computation of the F statistic without
     * the p-value or partitioning of the variance. If the {@code statistic} is not null
     * the method will record the F statistic in the array and return null.
     *
     * @param data Category summary data.
     * @param statistic Result for the F statistic (or null).
     * @return test result (or null)
     * @throws IllegalArgumentException if the number of categories is less than two; a
     * contained category does not have at least one value; or all categories have only
     * one value (zero degrees of freedom within groups)
     */
    private static Result aov(Collection<double[]> data, double[] statistic) {
        Arguments.checkCategoriesRequiredSize(data.size(), 2);
        long n = 0;
        for (final double[] array : data) {
            n += array.length;
            Arguments.checkValuesRequiredSize(array.length, 1);
        }
        final long dfwg = n - data.size();
        if (dfwg == 0) {
            throw new InferenceException(InferenceException.ZERO, "Degrees of freedom within groups");
        }

        // wg = within group
        // bg = between group

        // F = Var(bg) / Var(wg)
        // Var = SS / df
        // SStotal = sum((x - u)^2) = sum(x^2) - sum(x)^2/n
        //         = SSwg + SSbg
        // Some cancellation of terms reduces the computation to 3 sums:
        // SSwg = [ sum(x^2) - sum(x)^2/n ] - [ sum_g { sum(sum(x^2) - sum(x)^2/n) } ]
        // SSbg = SStotal - SSwg
        //      = sum_g { sum(x)^2/n) } - sum(x)^2/n
        // SSwg = SStotal - SSbg
        //      = sum(x^2) - sum_g { sum(x)^2/n) }

        // Stabilize the computation by shifting all to a common mean of zero.
        // This minimise the magnitude of x^2 terms.
        // The terms sum(x)^2/n -> 0. Included them to capture the round-off.
        final double m = StatisticUtils.mean(data);
        final Sum sxx = Sum.create();
        final Sum sx = Sum.create();
        final Sum sg = Sum.create();
        // Track if each group had the same value
        boolean eachSame = true;
        for (final double[] array : data) {
            eachSame = eachSame && allMatch(array[0], array);
            final Sum s = Sum.create();
            for (final double v : array) {
                final double x = v - m;
                s.add(x);
                // sum(x)
                sx.add(x);
                // sum(x^2)
                sxx.add(x * x);
            }
            // Create the negative sum so we can subtract it via 'add'
            // -sum_g { sum(x)^2/n) }
            sg.add(-pow2(s.getAsDouble()) / array.length);
        }

        // Note: SS terms should not be negative given:
        // SS = sum((x - u)^2)
        // This can happen due to floating-point error in sum(x^2) - sum(x)^2/n
        final double sswg = Math.max(0, sxx.add(sg).getAsDouble());
        // Flip the sign back
        final double ssbg = Math.max(0, -sg.add(pow2(sx.getAsDouble()) / n).getAsDouble());
        final int dfbg = data.size() - 1;
        // Handle edge-cases:
        // Note: 0 / 0 -> NaN : x / 0 -> inf
        // These are documented results and should output p=NaN or 0.
        // This result will occur naturally.
        // However the SS totals may not be 0.0 so we correct these cases.
        final boolean allSame = eachSame && allMatch(data);
        final double msbg = allSame ? 0 : ssbg / dfbg;
        final double mswg = eachSame ? 0 : sswg / dfwg;
        final double f = msbg / mswg;
        if (statistic != null) {
            statistic[0] = f;
            return null;
        }
        final double p = FDistribution.of(dfbg, dfwg).survivalProbability(f);

        // Support partitioning the variance
        // ni = size of each of the groups
        // nO=(1/(a−1))*(sum(ni)−(sum(ni^2)/sum(ni))
        final double nO = (n - data.stream()
                .mapToDouble(x -> pow2(x.length)).sum() / n) / dfbg;

        return new Result(dfbg, dfwg, msbg, mswg, nO, f, p);
    }

    /**
     * Return true if all values in the array match the specified value.
     *
     * @param v Value.
     * @param a Array.
     * @return true if all match
     */
    private static boolean allMatch(double v, double[] a) {
        for (final double w : a) {
            if (v != w) {
                return false;
            }
        }
        return true;
    }

    /**
     * Return true if all values in the arrays match.
     *
     * <p>Assumes that there are at least two arrays and that each array has the same
     * value throughout. Thus only the first element in each array is checked.
     *
     * @param data Arrays.
     * @return true if all match
     */
    private static boolean allMatch(Collection<double[]> data) {
        final Iterator<double[]> iter = data.iterator();
        final double v = iter.next()[0];
        while (iter.hasNext()) {
            if (iter.next()[0] != v) {
                return false;
            }
        }
        return true;
    }

    /**
     * Compute {@code x^2}.
     *
     * @param x Value.
     * @return {@code x^2}
     */
    private static double pow2(double x) {
        return x * x;
    }
}