Однофакторный дисперсионный анализ (One-way ANOVA) на Python

Задача дисперсионного анализа — одна из самых распространенных в практике обработки данных.  Ее типичная формулировка состоит в том, что требуется статистически оценить факт различия средних значений для нескольких групп.

Формулировка задачи дисперсионного анализа

Распространенность дисперсион­ного анализа связана с тем, что исследование воздействия фактора на явление формулируется как задача обнаружения изменений в средних значениях некоторых параметров под действием различных градаций фактора. Например, рассматривается эксперимент по установлению влияния какого-либо удобрения в трех различных концентрациях на рост посевов. Проводится опыт, в котором одна группа образцов не подвергается никакому воздействию, т.е. является контрольной группой; а остальные три — обрабатываются удобрением различной концентрации. Через определенное время измеряются абсолютные размеры растений во всех группах. Если влияние имело место быть, то средние значения размеров растений для групп будут различны. Таким образом, задача исследования воздействия удобрения сводится к  сравнению нескольких средних значений. Ввиду различия всхожести семян и многих других факторов, крайне маловероятно, что средние значения будут одинаковыми даже в отсутствии реального воздействия со стороны удобрения. Решить вопрос о том, являются ли наблюдаемые различия результатом действия фактора (различных концентраций удобрения) или обусловлены природой семян (естественной изменчивостью их всхожести) призвана схема однофакторного дисперсионного анализа.

Yet another anova...:)

Базовым для статистической процедуры однофакторного дисперсионного анализа является ответ на вопрос: какова вероятность наблюдать имеющиеся данные (например, выполненные измерения средних для четырех различных групп объектов), если бы  группы формировались из распределений с одинаковыми (теоретическими) средними значениями.

В привычной формулировке задача дисперсионного анализа выглядит следующим образом: в условиях нормальности и равенства дисперсий для данных в каждой из групп, необходимо оценить вероятность реализации наблюдаемой (полученной по выборочным данным) картины в предположении, что средние значения (математические ожидания) для всех групп одинаковы. Если оценка такой вероятности оказывается очень малой, т.е. вероятность наблюдать такие выборочные данные при условии равенства средних мала, то естественно гипотеза о равенстве средних значений должна быть отвергнута; в противном случае — делают заключение, что наблюдаемые данные не противоречат гипотезе о равенстве средних. Оценку такой вероятности как правило обозначаются p-value.

Требование нормальности данных и равенства дисперсий для сравниваемых групп накладывают определенные ограничения на применение схемы дисперсионного анализа. И хотя  классический дисперсионный анализ проявляет устойчивость при уклонении внутригрупповых распределений от условий нормальности и даже равенства дисперсий, оценки p-value, получаемые в условиях таких уклонений, могут быть искажены.

В случае невыполнения условий нормальности существует непараметрическая форма однофакторного дисперсионного анализа, предложенная Крускалом и Уоллисом в 1952 г.

Однако, прежде чем переходить к непараметрическим критериям можно попробовать использовать различные преобразования данных (логарифмирование, взятие арксинуса и т.п.), способные привести их к распределениям близким к нормальному. В биологических исследованиях часто встречаются данные, распределенные по лог-нормальному закону (логарифмирование которых приводит к нормально распределенным данным).

Резюмируем. Строго, классический однофакторный дисперсионный анализ применим, когда:

  1. В отношении выборочных данных не отвергается гипотеза об их нормальности (для каждой выборки);
  2. Выборки имеют равные теоретические дисперсии;

На практике такие условия будут выполняться далеко не всегда, и в общем случае возможны следующие ситуации: 1) условия 1-2 выполняются; 2) условие 1 выполняется, условие 2 не выполняется; 3) условие 1 не выполняется, условие 2 выполняется; 4) оба условия не выполняются.

Таким образом, создавая "универсальную" реализацию схемы дисперсионного анализа необходимо учесть все эти ситуации.

В случае 2) можно воспользоваться критерием Уэлча (scipy.stats.ttest_ind c опцией equal_var=False); или прибегнуть к непараметрическому тесту Крускала-Уоллиса.

В случае 3), в силу устойчивости классической схемы при уконениях от нормальности можно использовать scipy.stats.f_oneway; Также альтернативой здесь является переход к непараметрическому тесту, как и случае 2).

В случае 4) — выбор падает на непараметрический тест Крускала-Уоллиса (scipy.stats.kruskal).

Для проверки нормальности данных можно использовать весьма чувствительный тест Шапиро-Уилка (scipy.stats.shapiro), а для статистического подтверждения равенства дисперсий тест Левена (scipy.stats.levene).

Итак, все необходимые тесты уже реализованы в пакете SciPy (можно еще посмотреть в сторону пакета statsmodels, содержит много полезных статистических алгоритмов), осталось только составить из них композицию.

Что должен уметь алгоритм комплексного однофакторного анализа: 1) проверять данные на нормальность при заданном уровне значимости; 2) проверять равенство дисперсий у выборок; 3) подбирать подходящее преобразование данных, при котором обеспечивается их нормальность и\или равенство дисперсий (включим в него простые преобразования по умолчанию — извлечение корня, логарифмирование, взятие арксинуса); 4) в зависимости от ситуации применять подходящий тест: 1) классический тест Фишера, 2) непараметрический тест Крускалла-Уоллиса; 3) попарное сравнение по тесту Уэлча.

Реализуя составную схему дисперсионного анализа, удовлетворяющую 1)-4), получим что-то подобное:

from __future__ import print_function
import scipy.stats as st
import numpy as np
from itertools import combinations


def anova_ext(*args, **kwargs):
    '''Extended anova algorithm.
    
    Parameters
    ----------
        arg1, arg2, ..., argn -- samples to be compared;
        transormations -- array of tuples ('trans_name', trans_fun);
                          used to choose appropriate data transformation
        alpha -- level of significance, used when data is checked for normality

    Returns
    -------
        (value of statistic selected,
        p-value,
        statistic's name,
        nonnormal or normal (result of Shapiro-Wilk's test for normality),
        equal or unequal (result of checking for equalness of variances)
        name of transformation used
        )
    '''
    
    
    def _arcsine_trans(*args):
        maxel = max(map(lambda x: np.max(np.abs(x)), args))
        return map(lambda x: np.arcsin(np.asarray(x)/maxel), args)
    
    def _log_trans(*args):
        minel = min(map(lambda x: np.min(x), args))
        return map(lambda x: np.log(np.asarray(x) - minel + 1.0), args)
        
    def _sqrt_trans(*args):
        minel = min(map(lambda x: np.min(x), args))
        return map(lambda x: np.sqrt(np.asarray(x) - minel + 1.0), args)
    
    def _identity_trans(*args): return args
        
    if 'alpha' not in kwargs:
        alpha = 0.05
    else: 
        alpha = kwargs['alpha']
    if 'transformations' not in kwargs:
        transformations = [('identity', _identity_trans),
                         ('arcsine', _arcsine_trans),
                         ('log', _log_trans),
                         ('sqrt', _sqrt_trans)]
    else:
        if 'identity' not in map(lambda x:x[1], transformations):
            transformations = [('identity', lambda x: x)] + kwargs['transformations']
        else:
            transformations = kwargs['transformations']
    
    def _check_normality(*args):
        for arg in args:
            if st.shapiro(arg)[1] < alpha:
                return False
        return True

    def _check_variance_equality_normal(*args):
        return st.bartlett(*args)[1] > alpha
    
    def _check_variance_equality_nonnormal(*args):
        return st.levene(*args)[1] > alpha

    dtuple = []
    priorities = []
    for tname, tfun in transformations:
        transformed = tfun(*args)
        if _check_normality(*transformed):
            if _check_variance_equality_normal(*transformed):
                dtuple += [(tname, transformed, 'normal', 'equal')]
                priorities.append(1)
                break
            else:
                if 1 in priorities: break
                dtuple += [('identity', args, 'normal', 'unequal')]
                priorities.append(2)
        else:
            if _check_variance_equality_nonnormal(*transformed):
                if 2 in priorities: break
                dtuple += [(tname, transformed, 'nonnormal', 'equal')]
                priorities.append(3)
            else:
                if 3 in priorities: break
                dtuple += [('identity', args, 'nonnormal', 'unequal')]
                priorities.append(4)
    if 1 in priorities:
        data = dtuple[priorities.index(1)]
        fstat, pval = st.f_oneway(*data[1])
        return (fstat, pval, 'fisher', 'normal', 'equal', data[0])
    elif 2 in priorities:
        data = dtuple[priorities.index(2)]
        _pval = 1.0
        _fstat = 0.0
        for a, b in combinations(data[1], 2):
            fstat, pval = st.ttest_ind(a, b, equal_var=False)
            if _pval > pval:
                _pval = pval
                _fstat = fstat
        return (_fstat, _pval, 'welch-paired', 'normal', 'unequal', 'identity')
    elif 3 in priorities:
        data = dtuple[priorities.index(3)]
        fstat, pval = st.f_oneway(*data[1])
        return (fstat, pval, 'fisher', 'nonnormal', 'equal', data[0])
    elif 4 in priorities:
        data = dtuple[priorities.index(4)]
        fstat, pval = st.kruskal(*data[1])
        return (fstat, pval, 'kruskal', 'nonnormal', 'unequal', data[0])


if __name__ == '__main__':
    
    # three normally distributed samples
    a, b, c = np.random.randn(100), np.random.randn(100), np.random.randn(100)
    print('Applying extended anova scheme to 3 n.d. samples:', anova_ext(a,b,c))

    # We still use Fisher test... dut to its stability
    a, b, c = np.random.rand(200), np.random.rand(100), np.random.rand(300)
    print('Applying extended anova scheme to 3 u.d. samples:', anova_ext(a,b,c))
    
    # Crazy samples...
    a, b, c = [23,3,1,2,3,2,1], range(1000), np.linspace(20, 30, 100)
    print('Applying extended anova scheme to 3 "crazy" samples:', anova_ext(a,b,c))
 

Тестовые результаты работы скрипта следующие:

Applying extended anova scheme to 3 n.d. samples: (0.2087309265768498, 0.81173254734007605, 'fisher', 'normal', 'equal', 'identity')
Applying extended anova scheme to 3 u.d. samples: (0.31199078361214472, 0.73210754860488425, 'fisher', 'nonnormal', 'equal', 'identity')
/usr/local/lib64/python2.7/site-packages/scipy/stats/morestats.py:993: UserWarning: 6
  warnings.warn(str(ifault))
Applying extended anova scheme to 3 "crazy" samples: (11.219188123691646, 1.501145547249558e-05, 'fisher', 'nonnormal', 'equal', 'arcsine')

В случае, если нулевая гипотеза (о равенстве средних) отвергается, дальнейшее исследование данных предполагает отыскание тех пар средних, для которых данные различия статистически значимы. Это так называемый post hoc тест для процедуры однофакторного дисперсионного анализа. Кажется нет ничего проще, можно применить тест Стьюдента для всех возможных комбинаций средних, но только все несколько усложняется, если потребовать, чтобы обнаруженные различия проявляли значимость одновременно. Если достаточно только разрозненных фактов попарно-значимых различий, подойдет и тест Стьюдента; если требуется одновременность — используют тест Тьюки (или другие тесты), что будет рассмотрено в последующих заметках.

blog comments powered by Disqus