Разбор программы Statistica

Навигация для удобства чтения:

-Знакомство с программой
-Анализ
-Выводы

Устав от многочисленных вопросов, а чем же конкретно плоха программа Statistica, и почему стоит использовать языки программирования (например, R и Pyhton) для статистических расчетов, я решила написать исчерпывающий пост, обосновывающий проблему. Итак, давайте начнем.

Начнем с того, что так называемая программа "статистика" платная. Цену узнать просто на сайте невозможно, нужно оставлять свои контакты и ждать ответа официального представителя. Ну что же, я попросила коллегу поинтересоваться, вот ответ:

Итак, т.н. программа статистика стоит 290$ в год на одного пользователя. При этом такая однопользовательская версия не поддерживает удаленную работу по протоколу RDP (написано на сайте). То есть условно, подключившись из дома к рабочему компу, я уже не смогу использовать программу для стат расчетов, а для использования на домашнем ПК, нужно покупать снова. Цена не такая небольшая, учитывая что даже дорогущий матлаб (версия для академического использования) стоит 275 долларов на год или 550 навсегда. При этом у матлаба есть возможность купить студенческую лицензию за 29 долларов. Напоминаю, что R и Python бесплатны и доступны для скачивания любым желающим.
Недостаток 1: программа статистика платная и достаточно дорогая.

Естественно, я знаю, что обычно все используют взломанные версии программы, и со мной любезно поделились одним из дистрибутивов, версия STATISTICA 10.

Знакомство с программой

Открываем. При запуске она немного поругалась, видимо на то что кряк, но это решается, если нажать ОК на сообщение об ошибке. Далее появляется интересное предложение:

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

Открывается стартовое окно:

При попытке загрузить тестовый датасет я столкнулась с тем, что программа "статистика" не импортирует файлы в формате .txt и .csv

Тестовик у меня был хороший, только что благополучно загружала в R, выставила здесь настройки импорта, например tab как разделитель и точку как десятичный разделитель, однако импортировать не удалось. ОК, перевела табличку в эксель и загрузила ее. Может быть это не является прямо существенным минусом для всех, поскольку многие работают с анализом данных непосредственно в экселе, а все эти ваши текстовые файлы и тем более comma-separated - это профдеформация, связанная с написанием кода, а не человеческим восприятием. Но все равно отсутствие нормального импорта таб- или запятой-разделенных файлов это тревожный звоночек лично для меня (недостаток 2).

Анализ

Далее, я попыталась запустить разные виды анализа и сравнить результаты, полученные в программе статистика, с результатами вычисления в R и вручную по формулам. Признаться, мне хотелось бы найти очевидные ошибки при вычислении p-value и таким образом исчерпывающе доказать недопустимость использования программы для работы. Однако, все оказалось не так просто, поскольку откровенных ошибок не обнаружилось, что конечно хорошо для текущих пользователей программы. Но у меня есть ряд замечаний по работе тестов, выдаче результатов и отсутствии нормальной документации.

Тест Стьюдента

Начнем с всеми любимого и известного теста Стьюдента (т-тест). Я создала в R датасет с известными параметрами:

``` r

library(dplyr)

set.seed(3)

df_unbalanced <- data.frame(value = c(rnorm(n = 10, mean = 2, sd = 0.5), 

                   rnorm(n = 20, mean = 2.8, sd = 0.5)), 

                   condition = c(rep('a', each = 10), rep('b', each = 20)))

df_unbalanced %>% 

 group_by(condition) %>% 

 summarise(var = var(value), sd = sd(value), mean = mean(value))

#> # A tibble: 2 x 4

#>  condition  var  sd mean

#>  <chr>   <dbl> <dbl> <dbl>

#> 1 a     0.187 0.433 1.97

#> 2 b     0.156 0.396 2.64

```

Если вы запустите у себя этот кусочек кода, то получите тот же самый результат! Для подсчета суммарных статистик понадобится установить пакет dplyr. К сожалению, местный текстовый редактор немного некрасиво отформатировал код, но главное, что тут нужно отметить, это что дисперсии и стандартные отклонения достаточно близкие, а средние различаются. При этом я специально создала выборки разного размера, это будет важно в дальнейшем. Табличку сохраняем в Excel и импортируем в Statistica.
Выбираем t-test independent by groups. Пока что выставляем все параметры по умолчанию.

Следим за руками:

Теперь попробуем запустить t-test в R, также используя параметры по умолчанию:

``` r

t.test(df_unbalanced[df_unbalanced$condition == 'a', 1], df_unbalanced[df_unbalanced$condition == 'b', 1])

```

Итак, мы видим что в R по умолчанию запускается тест Велча (Welch), являющийся модификацией теста Стьюдента для образцов с неравной дисперсией. Об этом сообщает аутпут, что исключает неверную интерпретацию результатов, если читать заголовок конечно, а не сразу копировать p-value. Теперь давайте обратим внимание на значение p-value, полученное в т.н. программе статистика (0.0002) и посчитанное в R (0.0007). В данном случае значение отличается в три раза, а при создании датасета при запуске рандомизатора со значением set.seed(0), различие в p-value оказывалось на порядок (!).
Но погодите радоваться ошибке! Я разобралась, что дело в том, что "статистика" считает по умолчанию обычный тест Стьюдента, то есть с предположением о равенстве дисперсий, а в R по умолчанию тест Велча, как я упоминала выше.
В данном примере я специально подобрала значения так, чтобы дисперсии были близкими (0.187 и 0.156), при этом все равно p-value для теста Стьюдента и теста Велча отличаются в три раза.
Кроме того, я нашла более экстремальный пример, где тест Стьюдента показывает значимые различия (p-value = 0.02996), а тест Велча нет (p-value = 0.07775), ссылка на скрипт будет на гитхабе.

В R посчитать тест Стьюдента вместо теста Велча можно выставив параметр var.equal = T при запуске функции t.test(). Также я нашла, как запустить тест Велча в "статистике", причем конечно не юзеры статистики мне подсказали, а нашла в мануале по SPSS. Скидываю, где нажать волшебную галочку для запуска.

Ура, результаты вычислений в статистике и в R сошлись, но какой ценой? Мне не очень понравилось например, что нет явной подписи, что за тест использовался. Даже в экселе при подсчете т-теста необходимо выставить тип теста: парный, гомоскедастичный (то есть с предположением о равенстве дисперсий, "обычный" тест Стьюдента) или гетероскедастичный (тест Велча). То есть даже при подсчете в экселе нужно хоть чуточку отрефлексировать, что собираешься обсчитывать. В каком-то смысле эксель опережает по осознанности т-теста даже расчеты в R, поскольку нужно также задуматься, односторонний или двусторонний тест, а в R можно не задумываться и запустить все по умолчанию.
Обратите внимание, что подобное значительное расхождение p-value получается при несбалансированном размере выборок, в данном случае в одной выборке у нас 10 образцов, в другой 20, а при одинаковых объемах выборок результаты вычисления отличаются на тысячные и десятитысячные.

Однако по осознанности статистика получает жирный минус.

Это еще не все. Я не нашла возможности сделать односторонний т-тест, в саппорте по SPSS пишут, что это действительно невозможно сделать и предлагают просто поделить p-value, полученное по результатам двустороннего, на 2.
???
Ну ок, это конечно не очень сложно сделать, но все равно это как-то тупо, уж простите.

Мини-вывод: минус односторонний тест и не очень осознанный тест Стьюдента, с неочевидным вызовом теста Велча.

Тест Манна-Уитни

Теперь попробуем запустить тест Манна-Уитни на этих же данных.

Запустили:

Здесь у нас загадка трех p-value. Опять предлагаю попробовать угадать, какое из p-value нужно использовать, и что они значат. Запускаю wilcox.test в R с параметрами по умолчанию.

Можно увидеть, что U (или W) значение подсчитано обеими программами одинаково, кроме того, я самостоятельно тоже посчитала U и получила такой же результат. Вижу, что результат вычисления p-value в R совпал с третьим p-value "2*1 sided exact p", однако не очень ясно, что это значит в т.н. статистике. Также совсем неясно, что значат остальные два p-value и какое нужно использовать. Искала очень долго и мучительно в документации к программе и просто в гугле, но не нашла внятного определения. А загадка решилась в итоге методом научного тыка, чем я не очень довольна. Я предположила, что здесь ситуация может быть похожа на тест Велча и решила проверить параметры по умолчанию в R. Оказалось, что по умолчанию стоит параметр exact = NULL, что превращается в exact = TRUE, если число измерений в обеих выборках (суммарно) меньше 50. Это означает вычисление точного теста, о чем кстати написано в шапке теста, и соответствует 2*1 sided exact p в статистике, и теперь становится понятно, что это значит. Если выставить exact = FALSE, то при расчете p-value используется аппроксимация нормальным распределением (continuity correction), и значение p-value совпадает с первым и вторым p-value в "статистике".

Далее я попыталась разобраться, в чем различие между первым и вторым p-value, оказалось, что p-value могут отличаться, когда в наших выборках встречаются повторяющиеся значения (ties). Однако, при наличии таких повторяющихся значений в R появляется варнинг, что невозможно посчитать точное p-value.

Можно увидеть, что теперь у нас три разных p-value в статистике и четвертый вариант p-value в R. К счастью, они отличаются на десятитысячные, только третье "точное" вычисление p-value отличается на одну сотую. Однако, не очень хорошая новость в том, что никакая манипуляция параметрами в R мне не помогла воспроизвести результаты в статистике. Еще я попробовала вручную посчитать по формуле z-значение, у меня получилось -2.041008, что совпало с расчетом в статистике.

Однако при вычислении p-value без поправки на ties у меня получился вообще пятый вариант p-value = 0.0496, что еще и вовсе на границе значимости.

Это я обнаружила три недели назад и очень расстроилась, искренне пыталась разобраться, но в итоге решила, что я так никогда не опубликую разбор, поэтому решила оставить на следующий раз. В конце концов, отличия не такие фатальные, хотя получается, что "точное" вычисление теста в "статистике" занижает p-value и может привести к перепредсказанию. Стоит отметить, что тест Манн-Уитни может плохо работать с повторяющимися значениями и возможно, в определенных ситуациях есть смысл заменять его бутстреп-тестами.

Итак, я точно не могу сказать, какая из формул более правильная, но в пользу R наличие формулы, то есть открытость кода функции wilcox.test на гитхабе и возможность посмотреть формулы, по которым идет вычисление. То есть в R у нас возможность разобраться с математикой вычислений теста есть, в то время как программа "статистика" не предоставляет такой возможности и работает как черный ящик. Вот единственное, что я нашла по этой теме, здесь есть обоснование, что такое 2*1 sided exact p, но все равно нет формул вычисления p-value. Очень грустно.

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

Пермутационные тесты

Возможности их посчитать нет. И это заметила не только я.

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

Ну и в целом, мне работа с программой показалась намного более неудобной чем написание скриптов. Графический интерфейс с большим количеством кнопок, окон, непонятно как разложенных, оказался очень непростым для освоения. Понимаю, что существуют целые инструкции на ютубе и в конкретной лабе, пользующейся программой, но вопрос в том, зачем? Зачем с трудом осваивать инструмент, который годится для конкретных задач, если можно сразу осваивать язык программирования, который более универсален для сразу широкого класса задач (и сейчас я говорю скорее о Python :). Но для задач биоинформатики R будет более полезен, как минимум потому что ряд инструментов, например для анализа дифференциальной экспрессии (библиотеки DESeq2, edgeR) написаны только на R. Однако это тема для отдельного обсуждения и сейчас не об этом. Основная мысль, что для запуска той же двухфакторной ановы достаточно загрузить табличку в R и написать простейший код, состоящий из двух строчек. Не понимаю, почему многие так боятся "страшной" консоли, ведь страшный графический интерфейс "статистики" намного хуже!

Выводы

1. Программа платная с необходимостью продления лицензии каждый год.
2. Нет возможности нормально импортировать файлы в формате comma-separated и tab-separated.
3. Настройки по умолчанию не прописаны подробно. Например, при вычислении теста Стьюдента и теста Велча нет явного прописывания типа теста. Нет возможности посчитать односторонний т-тест.
4. При вычислении формул сумм квадратов в двухфакторной анове с несбалансированным дизайном используется формула, не описанная нигде, кроме как в инструкции к этой же программе (про это подробнее в другой раз).
5. Нет возможности считать пермутационные тесты, которые являются очень мощным статистическим инструментом
6. (Лично для меня). Крайне неинтуитивный дизайн, мне пришлось потратить какое-то время чтобы запустить даже самый простой анализ.

Подписывайтесь на мой телеграм-канал https://t.me/stats_for_science
возможно, постараюсь писать чаще. В планах у меня много разборов интересных вещей, покруче, чем разборы программы, только бы хватило времени и сил писать.

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

В начало