Октябрь 2016

Том 31 номер 10

Тесты - ANOVA на C#

Джеймс Маккафри

Джеймс МаккафриДисперсионный анализ (analysis of variance, ANOVA) — классический статистический метод, применяемый для заключения о том, являются ли средние трех или более групп равными, в ситуациях, где в наличии только одна выборка данных. Например, в университете есть три разных вводных курса в компьютерную науку. Каждый курс ведет один и тот же преподаватель, но использует другой учебник и другую методику преподавания. Вы хотите узнать, одинакова ли академическая успеваемость студентов.

Вы проводите экзамен, в ходе которого оцениваете подготовку в области компьютерной науки по шкале от 1 до 15, но, поскольку этот экзамен очень дорогостоящий и занимающий много времени, вы можете проэкзаменовать лишь шесть случайным образом выбранных студентов с каждого курса. Вы ведете экзамен и выполняете ANOVA на выборках, чтобы сделать заключение, являются ли средние оценки по всем трем курсам одинаковыми.

Если вы новичок в ANOVA, название этого метода может слегка сбить с толку, так как его цель — анализ средних значений наборов данных. ANOVA назван так потому, что закулисно он анализирует дисперсии (variances) для логических заключений о средних.

Хороший способ получить представление о том, что такое ANOVA и куда я клоню в этой статье, — взглянуть на демонстрационную программу на рис. 1. В демонстрации подготавливаются «зашитые» в код оценки для трех групп. Заметьте, что в Group1 только четыре оценки, а в Group3 — пять. Размеры выборок весьма часто неравные, поскольку испытуемые объекты могут выпадать либо данные могут быть утрачены или повреждены.

Демонстрация ANOVA на C#
Рис. 1. Демонстрация ANOVA на C#

В ANOVA две главные стадии. На первой вычисляются F-статистическое значение (критерий Фишера) и пара значений, называемых степенями свободы (degrees of freedom, df), используя данные выборки. На второй стадии на основе значений F и df определяется вероятность того, что все популяционные средние (population means) одинаковы (p-значение). Первая стадия сравнительно несложна, а вторая очень трудна.

В демонстрации значение F равно 15.884. В целом, чем больше F, тем менее вероятно, что все популяционные средние одинаковы. Вскоре я поясню, почему df = (2, 12). На основе F и df вычисленное p-значение равно 0.000425. Оно очень мало, поэтому вы заключили бы, что популяционные средние скорее всего не одинаковы. В этот момент вы могли бы выполнить дополнительные статистические проверки, чтобы определить, какие популяционные средние отличаются друг от друга. В случае демонстрационных данных оказывается, что Group1 (среднее выборки = 4.50) хуже Group2 (среднее = 9.67) и Group3 (среднее = 10.60).

Демонстрационная программа

Чтобы создать демонстрационную программу, я запустил Visual Studio, открыл File | New | Project и выбрал шаблон C# Console Application. В этой программе нет значимых зависимостей от .NET Framework, поэтому подойдет любая версия Visual Studio. После загрузки кода шаблона в окно редактора я переименовал в окне Solution Explorer файл Program.cs в более описательный AnovaProgram.cs, и Visual Studio автоматически переименовала класс Program за меня.

В начале кода я удалил все лишние выражения using, оставив только ссылку на пространство имен верхнего уровня System. Общая структура программы показана на рис. 2. Демонстрационная программа слишком длинна, чтобы представить ее здесь во всей полноте, но весь исходный код вы найдете в пакете, сопутствующем этой статье.

Рис. 2. Структура демонстрационной программы

using System;
namespace Anova
{
  class AnovaProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("Begin ANOVA using C# demo");
      // Подготавливаем данные выборки.
      // Используем данные для вычисления F и df.
      // Используем F и df для вычисления p-значения
      Console.WriteLine("End ANOVA demo");
    }

    static double Fstat(double[][] data, out int[] df) { . . }
    static double LogGamma(double z) { . . }
    static double BetaIncCf(double a, double b, double x)
      { . . }
    static double BetaInc(double a, double b, double x) { . . }
    static double PF(double a, double b, double x) { . . }
    static double QF(double a, double b, double x) { . . }
    static void ShowData(double[][] data, string[] colNames)
     { . . }
  }
}

Статический метод Fstat вычисляет и возвращает F-статистику на основе данных, хранящихся в объекте «массив массивов». Этот метод также вычисляет и возвращает два значения df в выходном параметре-массиве. ShowData — это просто небольшая вспомогательная функция для отображения средних по выборке.

Остальные пять методов используются для вычисления p-значения. QF — основной метод. Он вызывает метод PF, который в свою очередь обращается к методу BetaInc, а тот вызывает методы BetaIncCf и LogGamma.

После некоторых предваряющих сообщений через WriteLine метод Main подготавливает и отображает данные выборок:

double[][] data = new double[3][]; // 3 группы
data[0] = new double[] { 3, 4, 6, 5 };
data[1] = new double[] { 8, 12, 9, 11, 10, 8 };
data[2] = new double[] { 13, 9, 11, 8, 12 };
string[] colNames = new string[]
  {  "Group1", "Group2", "Group3" };
ShowData(data, colNames);

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

F-статистика и df вычисляются так:

int[] df = null;
double F = Fstat(data, out df);

В случае ANOVA df для набора данных представляет собой пару значений. Первое является K–1, где K — количество групп, а второе — N–K, где N — общее количество значений в выборке. Поэтому для демонстрационных данных df = (K–1, N–K) = (3–1, 15–3) = (2, 12).

Затем p-значение вычисляется и отображается следующим образом:

double pValue = QF(df[0], df[1], F);
Console.Write("p-value = ");

Итак, при выполнении ANOVA выражения вызовов очень просты. Но за кулисами осуществляется большой объем работы.

Вычисление F-статистики

Вычисление значения F-статистики (критерия Фишера) происходит в несколько этапов. Допустим, значения данных в выборке такие же, как в демонстрации:

Group1: 3.00, 4.00, 6.00, 5.00
Group2: 8.00, 12.00, 9.00, 11.00, 10.00, 8.00
Group3: 13.00, 9.00, 11.00, 8.00, 12.00

Первый этап — вычисление средних для каждой группы и совокупного среднего по всем значениям выборки. В случае демонстрационных данных имеем следующее:

means[0] = (3.0 + 4.0 + 6.0 + 5.0) / 4 = 4.50
means[1] = (8.0 + 12.0 + 9.0 + 11.0 + 10.0 + 8.0) / 6 = 9.67
means[2] = (13.0 + 9.0 + 11.0 + 8.0 + 12.0) / 5 = 10.60
gMean = (3.0 + 4.0 + . . . + 12.0) / 15 = 8.60

Определение метода Fstat начинается с:

static double Fstat(double[][] data, out int[] df)
{
  int K = data.Length;  // количество групп
  int[] n = new int[K]; // число элементов в каждой группе
  int N = 0;            // общее количество точек данных
  for (int i = 0; i < K; ++i) {
    n[i] = data[i].Length;
    N += data[i].Length;
  }
...

К этому моменту локальный массив n содержит количество значений в каждой группе, K — число групп, а N — общее количество значений во всех группах. Затем вычисляются групповые средние (group means) и помещаются в массив means, а вычисленное среднее по совокупности (grand mean) записывается в переменную gMean:

double[] means = new double[K];
double gMean = 0.0;
for (int i = 0; i < K; ++i) {
  for (int j = 0; j < data[i].Length; ++j) {
    means[i] += data[i][j];
    gMean += data[i][j];
  }
  means[i] /= n[i];
}
gMean /= N;

Следующий этап — вычисление суммы квадратов между группами (SSb) и среднего квадрата между группами (MSb). SSb — это взвешенная сумма квадратов разностей между средним каждой группы и совокупным средним. MSb = SSb / (K–1), где K — количество групп. В случае демонстрационных данных это выглядит так:

SSb = (4 * (4.50 - 8.60)^2) + (6 * (9.67 - 8.60)^2) +
  (5 * (10.60 - 8.60)^2) = 94.07
MSb = 94.07 / (3-1) = 47.03

Вот код, вычисляющий SSb и MSb:

double SSb = 0.0;
for (int i = 0; i < K; ++i)
  SSb += n[i] * (means[i] - gMean) * (means[i] - gMean);
double MSb = SSb / (K - 1);

Очередной этап — вычисление суммы квадратов внутри групп (SSw) и среднего квадрата внутри групп (MSw). SSw — это сумма квадратов разностей между каждым значением выборки и ее групповым средним. MSw = SSw / (N–K). В случае демонстрационных данных это выглядит так:

SSw = (3.0 - 4.50)^2 + . . + (8.0 - 9.67)^2 + . . +
  (12.0 - 10.60)^2 = 35.53
MSw = 35.53 / (15-3) = 2.96

Вот код, вычисляющий SSw и MSw:

double SSw = 0.0;
for (int i = 0; i < K; ++i)
  for (int j = 0; j < data[i].Length; ++j)
    SSw += (data[i][j] - means[i]) * (data[i][j] - means[i]);
double MSw = SSw / (N - K);

Последний этап — вычисление двух значений df и F-статистики. Двумя значениями df являются K–1 и N–K. F = MSb / MSw. В случае демонстрационных данных это выглядит так:

df = (K-1, N-K) = (3-1, 15-3) = (2, 12)
F = 47.03 / 2.96 = 15.88.

Демонстрационный код, вычисляющий df и F, представляет собой следующее:

...
  df = new int[2];
  df[0] = K - 1;
  df[1] = N - K;
  double F = MSb / MSw;
  return F;
} // Fstat

Полагаю, вы согласитесь, что при знании математических уравнений вычислений значений F-статистики и df по набору данных является механистичным и сравнительно простым.

Вычисление p-значения

Преобразование значений F-статистики и df в p-значение, сообщающее вам вероятность того, что все популяционные средние равны на основе данных выборки, которая дает F и df, несложное в принципе, но крайне трудное на практике. Я постараюсь пояснить это как можно короче, не вдаваясь в колоссальное количество деталей, которые потребовали бы огромного объема отдельных разъяснений. Взгляните на график на рис. 3.

Вычисление p-значения по значениям F-статистики и df
Рис. 3. Вычисление p-значения по значениям F-статистики и df

 

F(x) F(x)
F-Distribution, df1 = 4, df2 = 12 F-распределение, df1 = 4, df2 = 12
The area under curve from calculated F statistic to +infinity is probability that all means are equal Область под кривой от вычисленного значения F-статистики до +бесконечности отражает вероятность того, что все средние равны
F-statistic = 2.67 F-статистика = 2.67

 

Каждая возможная пара значений df определяет график, называемый F-распределением (или распределением отношения дисперсии). Форма F-распределения может меняться в широких пределах в зависимости от значений df. График на рис. 3 показывает F-распределение для df = (4, 12). Я использовал df = (4, 12) вместо df = (2, 12) из демонстрационных данных, потому что форма F-распределения при df = (2, 12) очень атипична..

Общая площадь под любым F-распределением составляет ровно 1.0. Если вам известно значение F-статистики, то p-значение — это площадь под F-распределением от F до +бесконечности. Немного сбивает с толку, что площадь под F-распределением от 0 до F-статистики часто называют PF, а площадь под F-распределением от F-статистики до +бесконечности (представляющей p-значение) — QF. Поскольку общая площадь под распределением равна 1, значит, PF + QF = 1. Оказывается, вычислить PF немного легче, чем QF, поэтому, чтобы найти p-значение (QF), как правило, вычисляют PF, а затем вычитают его из единицы и получают QF.

Вычислить PF чертовски трудно, но, к счастью, уже в течение десятилетий известны волшебные оценочные уравнения. Эти математические уравнения и сотни других вы найдете в знаменитом справочнике Абрамовица (Abramowitz) и Стигана (Stegun) «Handbook of Mathematical Functions». Среди программистов в научных областях эту книгу часто называют просто «A&S». Каждое уравнение в A&S имеет идентификационный номер.

В демонстрации метод PF на самом деле является не более чем оболочкой метода BetaInc:

static double PF(double a, double b, double x)
{
  double z = (a * x) / (a * x + b);
  return BetaInc(a / 2, b / 2, z);
}

Имя метода BetaInc образовано от сокращения «incomplete Beta» (неполная бета-функция). Метод BetaInc использует уравнения A&S под номерами 6.6.2 и 26.5.8. Из этих уравнений вызываются функции LogGamma и BetaIncCf. Функция LogGamma крайне трудна в объяснении и реализации. Если кратко, то математическая гамма-функция расширяет концепцию факториала до вещественных чисел. Как и факториалы, значения гамма-функции могут достигать астрономических величин, поэтому для их обработки обычно вычисляют логарифм гамма-функции, чтобы значения были меньше.

Вычисление LogGamma очень хитроумно, и вы можете использовать несколько алгоритмов. В демонстрационной программе используется алгоритм аппроксимации Ланцоша (Lanczos approximation) с (g=5, n=7). В справочнике A&S даны другие алгоритмы, способные вычислять LogGamma, но аппроксимация Ланцоша, не известная на момент публикации A&S, дает более точные результаты.

Имя метода BetaIncCf образовано от сокращения «incomplete Beta computed by continued fraction» (неполная бета-функция, вычисленная по непрерывной дроби). В демонстрационной программе для метода BetaIncCf применяется уравнение из A&S под номером 26.5.8.

Заключение

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

Существует несколько способов проверки этих допущений, но интерпретация их результатов является трудной задачей. Проблема заключается в следующем. Крайне маловероятно, что реальные данные имеют одинаковое нормальное распределение и в точности равные дисперсии, хотя ANOVA все равно работает, когда данные имеют несколько отличные дисперсии и небольшие отклонения от нормального распределения. Вывод состоит в том, что чрезвычайно сложно доказать допущения ANOVA, поэтому следует проявлять максимальную консервативность при интерпретации результатов.

ANOVA тесно связана с t-критерием (t-test). Последний определяет, действительно ли популяционные средние ровно двух групп одинаковы в ситуациях, где вы располагаете только данными выборки. Так что, если у вас есть три группы, как в этой демонстрации, то вместо использования ANOVA вы могли бы три проверки по t-критериям, сравнивая группы 1 и 2, группы 1 и 3 и группы 2 и 3. Однако этот подход не рекомендуется, поскольку он вводит так называемую ошибку Type 1 (ложноположительный результат).

Разновидность ANOVA, рассмотренная в этой статье, называется однонаправленной (или однофакторной) ANOVA. Другой метод, двунаправленная ANOVA, используется при наличии двух факторов.

ANOVA основана на вычислении значения F-статистики по набору данных. Существуют другие статистические проверки, в которых применяется F-статистика (критерий Фишера). Например, вы можете использовать F-статистику для логического определения, одинаковы ли дисперсии двух групп данных.


Джеймс Маккафри (Dr. James McCaffrey) — работает на Microsoft Research в Редмонде (штат Вашингтон). Принимал участие в создании нескольких продуктов Microsoft, в том числе Internet Explorer и Bing. С ним можно связаться по адресу jammc@microsoft.com.

Выражаю благодарность за рецензирование статьи экспертам Microsoft Крису Ли (Chris Lee) и Кирку Олинику (Kirk Olynyk).