6 Неделя 2, День 5

6.1 Стиль кода

После того, как мы научились делать всякие сложные вещи в R, нужно снова вернуться к тому, что иногда кажется несколько занудным и неинтересным: к стилю написания кода и workflow. После того, как путем проб и ошибок нужные результаты посчитаны, а подходящий график нарисован, очень хочется все закрыть и поскорее послать туда, куда это нужно было послать еще неделю назад. При сумбурном образе жизни современного ученого мало кому нравится идея о том, что после написания работающего кода нужно выдохнуть и все перепроверить, переписать код (если что-то можно сделать лучше), закомментировать сложные куски кода. Тем не менее, это не займет много времени, а в дальнейшем много раз окупится, когда придется возвращаться к коду или в нем окажется ошибка (а такое, как Вы уже могли догадаться, бывает часто).

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

Например, можно создать переменные с именем T и F, но делать этого не стоит, потому что это перепишет стандартные значения T и F (TRUE и FALSE). И если хотя где-то в коде Вы используете T вместо TRUE, то ничего не будет работать, придется потратить очень много времени, чтобы понять, в чем была ошибка. По этой же причине, кстати говоря, лучше не писать сокращать TRUE до T, а FALSE до F.

Есть же правила, которые не связаны с тем, что они могут привести к ошибкам, а только с тем, что такой код становится проще читать. Особенно хорошо, если все придерживаются более-менее одного стиля, но, увы, такого в R нет: Вы уже могли заметить, что разные пакеты предлагают свой стиль, немного отличающийся от других.

Поэтому вполне нормально, если Вы выработаете свой стиль с опытом. Главное, по возможности придерживаться его. Довольно универсальный стиль описан Хэдли Уикхемом.

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

  • Для файлов. Не использовать пробелы и знак  - лучше использовать нижнее подчеркивание или дефис. Использовать только маленькие буквы. В начале названия можно использовать числа, чтобы обозначить, в каком порядке нужно запускать скрипты.
  • Для переменных. Использовать короткие, но осмысленные названия переменных. Есть несколько вариантов, как разделять слова в названии переменной: использование больших букв, точка или нижнее подчеркивание.
SomeNumber <- 2 + 2 #можно, но не очень удобно
some.number <- 2 + 2 #удобно, но не рекомендуется
some_number <- 2 + 2 #хорошо

В целом, не стоит использовать вариант с точкой (точка обычно означает метод для generic функции), а вариант с нижнем подчеркиванием - самый популярный. Не используйте T, c и название каких-либо привычных функций для называния переменных. Не только из-за того, что это может привести к ошибкам, но и из-за того, что читающий код может не понять, что происходит.

  • Пробелы. Всегда вокруг арифметических операторов (но не вокруг :!), после, но не до запятой, не ставить пробелов около $ и ::

Кроме того, в RStudio есть удобное сочетание клавиш для форматирования кода: ctrl + shift + A (Windows) или cmd + shift + A (MacOS). Если выделить код и зажать эти клавиши, то вот такой код…

for  (i in 1 : 10)  { if(i%%2== 0)
print(
paste(i , "is even") )
                    }
## [1] "2 is even"
## [1] "4 is even"
## [1] "6 is even"
## [1] "8 is even"
## [1] "10 is even"

…превратится в такой:

for  (i in 1:10)  {
  if (i %% 2 == 0)
  print(paste(i , "is even"))
}
## [1] "2 is even"
## [1] "4 is even"
## [1] "6 is even"
## [1] "8 is even"
## [1] "10 is even"

6.2 Вопроизводимость исследований

6.2.1 Спорные исследовательские практики

Если Вы откроете какой-нибудь более-менее классический учебник по психологии, то увидите там много результатов исследований, которые в дальнейшем не удалось воспроизвести. Эта проблема накрыла академическое психологическое сообщество относительно недавно и стала, пожалуй, самой обсуждаемой темой среди ученых-психологов. Действительно, о чем еще говорить, если половина фактов твоей науки попросту неверна? Об историческом смысле методологического кризиса, конечно же. Получается, у теорий нет никакого подтверждения, а в плане знаний о психологических фактах мы находимся примерно там же, где и психологи, которые закупали метрономы и тахистоскопы почти полтора века назад.

Оказалось, что то, как устроена академическая наука, способствует публикации ложноположительных результатов. Если теоретическая гипотеза подтверждается, то это дает ученому больше плюшек, чем если гипотеза не подтверждается. Да и сами ученые как-то неосознанно хотят оказаться правыми. Ну а перепроверка предыдущих достижений в психологии оказалась не в почете: всем хочется быть первопроходцами и открывать новые неожиданные феномены, а не скрупулезно перепроверять чужие открытия. Все это привело психологию к тому, что в ней (да и не только в ней) закрепилось какое-то количество спорных исследовательских практик (questionable research practices). Спорные практики на то и спорные, что не все они такие уж плохие, многие ученые даже считают, что в них нет ничего плохого. В отличие, например, от фальсификации данных - это, очевидно, совсем плохо.

Вот некоторые распространенные спорные исследовательские практики:

  • Выбор зависимых переменных пост-фактум. По своей сути, это проблема скрытых множественных сравнений: если у нас много зависимых переменных, то можно сделать вид, что измерялось только то, на чем нашли эффект. Поэтому стоит задуматься, если возникает идея измерять все подряд на всякий случай - что конкретно предполагается обнаружить? Если конкретной гипотезы нет, то нужно так и написать, делая соответствующие поправки при анализе.

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

  • Остановка набора испытуемых при достижении уровня значимости (Optional stopping). Представьте себе, что Вы не обнаружили значимых результатов, но p-value болтается около .05. Тогда Вам захочется добрать парочку испытуемых. А потом еще. И еще немного. Проблема с таким подходом в том, что рано или поздно Вы получите статистически значимые результаты. В любом случае, даже если эффекта нет. На самом деле, эта проблема не так сильно влияет на результаты как может показаться, но это все-таки достаточно плохая практика, а ее применение увеличивает количество опубликованных ложно-положительных результатов. Количество необходимых испытуемых для исследования лучше определять заранее (см. 6.2.3).

  • Неправильное округление до .05. Всякий раз, когда видите p = .05 будьте внимательны: вообще-то p не может быть равен именно .05. Либо автор не очень этого понимает, либо просто p больше, а не меньше .05. Например, .054, что потом округляется до .05 и преподносится как статистически значимые результаты.

  • Использование односторонних тестов для сравнения средних. Эта практика похожа на предыдущую: исследователь получил p > .05, но меньше, чем .1. Недобросовестный исследователь, который хочет опубликоваться проводит односторонний т-тест вместо двустороннего, т.е. фактически просто делит р на 2. Совсем недобросовестные даже не пишут, что они использовали односторонний т-тест, хотя по умолчанию все используют двусторонний.

Данный список не претендует на полноту. Этих практик стоит избегать и других от этого отучать. Ну а если Вы думаете, что никто не заметит, если Вы так сделаете, то это не так.

Например, последние две практики можно легко обнаружить с помощью сайта http://statcheck.io. Этот сайт делает магию: нужно кинуть ему статью, он распознает в ней статистические тесты и пересчитывает их. На самом деле, ничего принципиально сложного сайт не делает, он сделан с помощью R и уже знакомых нам инструментов. С помощью специальных пакетов из файлов вытаскивается текст, в тексте с помощью регулярных выражений находятся паттерны вроде “t(18) = -1.91, p = .036” - в большинстве журналов используется очень похожее форматирование статистических результатов. Зная степени свободы и т-статистику, можно пересчитать p-value. В данном случае это можно посчитать вот так:

pt(-1.91, df = 18)*2 #умножаем на 2, потому что двусторонний тест
## [1] 0.07219987

Ну а дальше остается сравнить это с тем, что написали авторы. Например, бывает как в данном случае, что пересчитанный p-value в два раза больше того, что написали авторы. Это означает, скорее всего, что авторы использовали односторонний критерий. Если они об этом нигде не сказали, то это очень плохо.

Самостоятельное задание:

  1. Найдите научную статью, в которой, как Вам кажется, авторы использовали спорные исследовательские практики. Если такой статьи нет на примете, то по ищите в научных базах статьи из журналов с низким импакт-фактором - средним количеством цитирований на статью.

  2. Проверьте статью через http://statcheck.io. Проинтерпретируйте все несовпадения, если такие имеются.

6.2.2 Вопроизводимые исследования в R

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

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

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

Конечно, не все готовы к таким разворотам в исследовательской практике. Но лучше быть готовым, потому что в какой-то момент может оказаться, что новые практики станут обязательными. И тут R будет весьма кстати: можно выкладывать данные c RMarkdown документом, который будет сразу собирать из этого статью и графики. Данные со скриптами можно выкладывать на GitHub - это удобно для коллективной работы над проектом. Другой вариант - выкладывать данные и скрипты для анализа на сайте osf.io. Это сайт был специально создан как платформа для публикации данных исследований и скриптов для них.

Самостоятельное задание:

  1. Найдите статью с кодом по на R по интересуемой Вам теме. Посмотрите код, попытайтесь его понять и запустить самостоятельно.

6.2.3 Статистическая мощность

Чтобы избежать optional stopping, нужно определять размер выборки заранее. Как это сделать? Наиболее корректный способ предполагает использование анализа статистической мощности (statistical power analysis). Для этого понадобится пакет pwr.

install.packages("pwr")

Этот пакет предоставляет семейство функций для расчета мощности, размера эффекта, уровня \(\alpha\) или нужного размера выборки для разных статистических тестов.

Статистическая мощность - вероятность обнаружить статистически значимый эффект, если он действительно есть. Размер эффекта - собственно, размер эффекта в исследовании, обычно в универсальных единицах. Например, размер различия средних обычно измеряется в стандартных отклонениях. Это позволяет сравнивать эффекты в разных исследованиях и даже эффекты в разных областях науки.

Если задать 3 из 4 чисел (1 - статистическая мощность - стандартно это .8; 2 - уровень \(\alpha\) - стандартно .05; размер эффекта; размер выборки), то функция выдаст недостающее число.

Я очень советую поиграться с этим самостоятельно. Например, если размер эффекта в единицах стандартных отклонений \(d_{Cohen's} = 2\), то сколько нужно испытуемых, чтобы обнаружить эффект для двухвыборочного т-теста с вероятностью .8?

library(pwr)
pwr.t.test(d = 2, power = 0.8, type = "two.sample")
## 
##      Two-sample t test power calculation 
## 
##               n = 5.089995
##               d = 2
##       sig.level = 0.05
##           power = 0.8
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

Всеего 6 в каждой группе (округляем в большую сторону)!

pwr.t.test(d = 2, power = 0.8, type = "paired")
## 
##      Paired t test power calculation 
## 
##               n = 4.220726
##               d = 2
##       sig.level = 0.05
##           power = 0.8
##     alternative = two.sided
## 
## NOTE: n is number of *pairs*

Еще меньше, если использовать within-subject дизайн исследования и зависимый т-тест.

А если эффект маленький - всего \(d_{Cohen's} = .2\)?

pwr.t.test(d = .2, power = 0.8, type = "two.sample")
## 
##      Two-sample t test power calculation 
## 
##               n = 393.4057
##               d = 0.2
##       sig.level = 0.05
##           power = 0.8
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

394 испытуемых в каждой группе!

Можно проверить такие расчеты самостоятельно с помощью симуляции данных.

mean(replicate(1000, t.test(rnorm(394, 100, 15), rnorm(394, 103, 15))$p.value < 0.05))
## [1] 0.8

Действительно, если много раз делать случайные выборки из двух нормальных распределений с средними, отличающимися на 0.2 стандартных отклонения, то примерно в 80% случаев вы получите статистически значимые различия!

Самостоятельное задание:

  1. Представьте, что эффекта нет - две выборки взяты из нормального распределения. Посчитайте, как часто t.test() будет выдавать ложноположительный результат при \(\alpha\) = .05.
mean(replicate(10000, t.test(rnorm(15), rnorm(15))$p.value < 0.05))
## [1] 0.0473
  1. Будет ли это зависеть от размера выборки?
mean(replicate(10000, t.test(rnorm(150), rnorm(150))$p.value < 0.05))
## [1] 0.0508
mean(replicate(10000, t.test(rnorm(1500), rnorm(1500))$p.value < 0.05))
## [1] 0.0495
  1. Представьте, что мы работает не с данными нормального распределения, а с rlnorm().
mean(replicate(10000, t.test(rlnorm(15), rlnorm(15))$p.value < 0.05))
## [1] 0.0328
mean(replicate(10000, t.test(rlnorm(150), rlnorm(150))$p.value < 0.05))
## [1] 0.0468