План занятий

(2 пары по 1:30):

Welcome to RStudio!

Интерфейс, проекты и git (вкратце)


Распределения данных

Нормальное распределение

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

Один из классических примером подбора выборки – random sampling. Если мы действительно рандомно выбираем из генеральной совокупности, то, скорее всего, наши данные распределены нормально.

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

Какими свойствами обладает нормальное распределение? Что отличает его от других?

Spoiler warning

Нормальное распределение:

  • Унимодально (одна мода, один пик)
  • Симметрично, не скошено вправо или влево, пик ровно в середине распределения (skewness = 0)
  • Средней “сплющенности”, не вытянуто или сплющено по оси y (kurtosis = 0)
  • Описывается определенной функцией вероятности: \[f(x) = \frac{1}{(\sigma\sqrt{2 \pi})} e^{-(\frac{(x - \mu)^2}{2 \sigma^2})}\]


Cтандартное нормальное распределение

Отдельный вид нормального распределения – стандартное нормальное распределение, у которого помимо всех свойств нормального распределения M = 0 и sd = 0. Именно к нему мы приводим наше нормальное распределение, когда делаем Z-преобразование. Такая трансформация дает возможность сравнивать совершенно разные распределения друг с другом. Такое распределение в основном теоретическое, применяемое для удосбвта сравнения и почти не встречается в природе.

t-распределение

Близкое к нормальному распределению, но не совсем его повторяющее – t-распределение, или распределение Стьюдента. t-распределение описывается немного другой формулой:\[f(x) = Gamma((n+1)/2) / (sqrt(n pi) Gamma(n/2)) (1 + x^2/n)^-((n+1)/2)\] Пусть формала не страшит, ее все равно считает за нас R. Важно понимать, что t-распределение очень похоже на нормальное, но не равно ему, и чем меньше выборка, тем сильнее заметна эта разница. При больших выбороках (условно при n>30) разница почти незаметна.

Спасибо Толе Карпову за шайни https://gallery.shinyapps.io/dist_calc/

Почему нам так важно, чтобы данные описывались каким-либо распределением?

Это сильно упрощает работу с этими данными и позволяет сравнивать данные между собой. Буквально дает основания для сравнения. Мы уже знаем, что данные, распределенные нормально, а также данные, распределенные в соответствии с t-распределением, описываются математическими законами – и это дает нам отличные опорные точки для сравнения. В данных, описываемых стандартных отклонением, параметров два: выборочное среднее (M, mean) и выборочное стандартное отклонение (sd, standart deviation). И именно они используются для сравнения выборок между собой.

Хорошие новости – распределений много. Плохая – очень много. http://www.math.wm.edu/~leemis/chart/UDR/UDR.html

Генерация данных из нормального распределения и визуализация

Вы уже должны были видеть, как выглядит нормальное распределение – а теперь давайте попробуем нарисовать его сами.

Предположим, что мы придумали супер-пупер инновационную методику повышения интеллекта. Например, снабдили студентов определенного университета новой методикой повышения осознанности, и теперь они решают задачки со всем обстоятельным подходом к делу. Мы хотим проверить, насколько наша методика работает, поэтому нам нужны будут данные тестов IQ у этих студентов и контрольной группы (все совпадения случайны, а исследовать интеллект в вакууме очень тяжело, но зато он стал популярным примером из-за когда-то кем-то на каком-то тесте полученных данных, что IQ распределен нормально со средним 100 и стандартным отклонением 15)

rnorm()

n = 100
M = 100
sd = 15
IQdataFromRnorm <- round(rnorm(n, mean = M, sd = sd))
IQdataFromRnorm
##   [1]  70 135 102 100  97 109  92  97  89  59 101  99 105  93  64 109 102 106
##  [19] 108 110 114  94 100  85  95  99 124 110  83 107  79  86 129  88  89 115
##  [37] 114  83 117 129 103 102 113 110 105  87 119  94 105 128 102  71  95 110
##  [55]  98 102  94  96  81  87  76  96 114  92 101 114  86 134 125  97  79 125
##  [73]  86 110  78 107 114  95  90  88 104 110  91 100 124  84  60 100 104  88
##  [91] 100  98 119  98 115  83  74 116  97  94
plot(IQdataFromRnorm, 51:150)

hist(IQdataFromRnorm, probability = T)
lines(x = density(IQdataFromRnorm), col="red")
abline(v = median(IQdataFromRnorm), col = 'green')
abline(v = mean(IQdataFromRnorm), col = 'blue')

psych::describe(IQdataFromRnorm)
##    vars   n  mean   sd median trimmed   mad min max range  skew kurtosis   se
## X1    1 100 99.55 15.6    100   99.66 14.83  59 135    76 -0.12     0.01 1.56

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

IQtable <- data.frame(
  subject = seq(1:100),
  gender = rep(c('f','m'),50),
  age = round(rnorm(100,22,0.8)),
  iq = IQdataFromRnorm
)

head(IQtable)
##   subject gender age  iq
## 1       1      f  22  70
## 2       2      m  24 135
## 3       3      f  22 102
## 4       4      m  23 100
## 5       5      f  23  97
## 6       6      m  21 109

dnorm()

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

data <- seq(50,150,1)
data
##   [1]  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65  66  67
##  [19]  68  69  70  71  72  73  74  75  76  77  78  79  80  81  82  83  84  85
##  [37]  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100 101 102 103
##  [55] 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121
##  [73] 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139
##  [91] 140 141 142 143 144 145 146 147 148 149 150
IQdataFromDnorm <- dnorm(data, mean = M, sd = sd)
IQdataFromDnorm
##   [1] 0.0001028186 0.0001281199 0.0001589392 0.0001962978 0.0002413624
##   [6] 0.0002954566 0.0003600704 0.0004368688 0.0005276968 0.0006345818
##  [11] 0.0007597324 0.0009055313 0.0010745239 0.0012694000 0.0014929687
##  [16] 0.0017481259 0.0020378140 0.0023649729 0.0027324837 0.0031431044
##  [21] 0.0035993978 0.0041036534 0.0046578051 0.0052633439 0.0059212307
##  [26] 0.0066318093 0.0073947223 0.0082088348 0.0090721655 0.0099818310
##  [31] 0.0109340050 0.0119238944 0.0129457370 0.0139928197 0.0150575218
##  [36] 0.0161313816 0.0172051884 0.0182690978 0.0193127702 0.0203255285
##  [41] 0.0212965337 0.0222149735 0.0230702595 0.0238522286 0.0245513427
##  [46] 0.0251588818 0.0256671250 0.0260695129 0.0263607894 0.0265371151
##  [51] 0.0265961520 0.0265371151 0.0263607894 0.0260695129 0.0256671250
##  [56] 0.0251588818 0.0245513427 0.0238522286 0.0230702595 0.0222149735
##  [61] 0.0212965337 0.0203255285 0.0193127702 0.0182690978 0.0172051884
##  [66] 0.0161313816 0.0150575218 0.0139928197 0.0129457370 0.0119238944
##  [71] 0.0109340050 0.0099818310 0.0090721655 0.0082088348 0.0073947223
##  [76] 0.0066318093 0.0059212307 0.0052633439 0.0046578051 0.0041036534
##  [81] 0.0035993978 0.0031431044 0.0027324837 0.0023649729 0.0020378140
##  [86] 0.0017481259 0.0014929687 0.0012694000 0.0010745239 0.0009055313
##  [91] 0.0007597324 0.0006345818 0.0005276968 0.0004368688 0.0003600704
##  [96] 0.0002954566 0.0002413624 0.0001962978 0.0001589392 0.0001281199
## [101] 0.0001028186
plot(data, IQdataFromDnorm)
lines(x = data, y = IQdataFromDnorm, col="red")
abline(v = median(data), col = 'green')
abline(v = mean(data), col = 'blue')

Внимательно посмотрите на данные, которые выдаются по вызову dataFromDnorm. Что в них особенного? Является ли распределение dataFromDnorm стандартным нормальным распределением?

Возьмем пока данные, сгенерированные функцией rnorm IQdataFromRnorm, и проверим, насколько они нормальны.

Проверка данных на нормальность

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

Возьмем наш уже известный график распредедения IQdataFromRnorm

hist(IQdataFromRnorm, probability = T)
lines(x = density(IQdataFromRnorm), col="red")

Визуально мы уже посмотрели на гистограмму и функцию плотности. Выглядит неплохо!

На что еще можно посмотреть?

Вторым действием после рисования гистограммы может быть рисование QQ plot: график, показывающий, как распределены данные по квартилям.

qqnorm(IQdataFromRnorm, pch = 1, frame = FALSE)
qqline(IQdataFromRnorm, col = "red", lwd = 2)

В целом, этого скорее достаточно, чтобы оценить нормальность данных. Если функциая плотности выглядит как равномерный колоколо и QQplot расположился вдоль линии y=x – примите поздравления! Ваша выборка близка к нормальной. И можно применять статистические тесты, которые требует нормального распределения данных. Например, t-test, который мы будем здесь рассматривать.

Если вы очень хотите, можно провести тест на нормальност данных, Shapiro-Wilk или Колмогоров-Смирнова. Но учтите, что они очень чувствительны и свидетельствуют о ненормальности даже там, где отклонение от нормального распределение настолько не существенно, что им обычно пренебрегают.

shapiro.test(IQdataFromRnorm)
## 
##  Shapiro-Wilk normality test
## 
## data:  IQdataFromRnorm
## W = 0.99114, p-value = 0.756

А теперь давайте посмотрим на то, как данные могут выглядеть в реальности:

 library(SimMultiCorrData)
IQNotNormalDataFrame <- nonnormvar1(method = c("Fleishman", "Polynomial"), means = M, vars = sd,
  skews = 1, skurts = 3, fifths = 0, sixths = 0, Six = NULL,
  cstart = NULL, n = 25, seed = 1234)
## Constants: Distribution  1  
## 
## Constants calculation time: 0.025 minutes 
## Total Simulation time: 0.025 minutes
IQNotNormalData <- IQNotNormalDataFrame$continuous_variable$V1
IQNotNormalData
##  [1]  96.46461 101.50508 105.70789  92.54114 102.17805 102.53764  98.39658
##  [8]  98.48625  98.42933  97.42091  98.71034  97.09356  97.76785 100.63301
## [15] 104.94617  99.97228  98.60074  97.35681  97.58158 117.32967 100.90947
## [22]  98.66650  98.83018 102.31891  98.02269
psych::describe(IQNotNormalData)
##    vars  n  mean   sd median trimmed  mad   min    max range skew kurtosis   se
## X1    1 25 100.1 4.55  98.67   99.54 1.94 92.54 117.33 24.79 2.06     5.76 0.91
hist(IQNotNormalData, probability = T)
lines(x = density(IQNotNormalData), col="red")

qqnorm(IQNotNormalData, pch = 1, frame = FALSE)
qqline(IQNotNormalData, col = "red", lwd = 2)

shapiro.test(IQNotNormalData)
## 
##  Shapiro-Wilk normality test
## 
## data:  IQNotNormalData
## W = 0.77555, p-value = 9.132e-05

Сравнение двух средних

Перейдем к тому, чтобы проверить эффективность нашей супер-пупер методики и статистически сравнить экспериментальную выборку, в которой применялась методика, с контрольной.

Возьмем уже знакомые нам данные IQdataFromRnorm и приблизим их к реальным данным экспериментов: уменьшим выборку до 30 человек. Пусть это будет наша контрольная группа.

nreal <- 30
IQdata_real1 <- round(rnorm(nreal, mean = M, sd = sd))
IQdata_real1
##  [1]  78 109  85 100  86 117  93  89  92  76  82  67  80  96  93 122  84  87  96
## [20]  85  85  83  81  92  93  73  91  83  85  98
hist(IQdata_real1, probability = T)
lines(x = density(IQdata_real1), col="red")

Теперь сгенерируем экспериментальную группу: пусть мы обработали данные и получили среднее M=115 во стандартным отклонением sd=13.

M2 <- 115
sd2 <- 13
IQdata_real2 <- round(rnorm(nreal, mean = M2, sd = sd2))
IQdata_real2
##  [1] 122 136 105 136 100 124 148 115 106 115 138 100 133 132 119 115 109 110 123
## [20] 142 113  97 106 118 111 113 113  97 113 126
hist(IQdata_real2, probability = T)
lines(x = density(IQdata_real2), col="red")

Null hypothesis testing

В классической частотной парадигме обработки данных сравнения происходят на основании null hypothesis testing, принципах тестирования нулевой гипотезы:

  • Нулевая гипотеза, H0 – всегда о том что нет различий.
  • Альтернативная гипотеза, H1 – о том, что сравниваемые параметры статистически значимо отличаются друг от друга.
  • Уровень статистической значимости, alpha, определяется заранее. В наших экспериментах чаще всего alpha=0.05. Этот процент возможных ложноположительных выводов: если мы провели эксперимент 100 раз, при alpha=0.05 мы допускаем, что в 5% случаев мы можем получить значимый результат просто случайно.
  • p value – вероятность получить эффект, когда он на самом деле отсутствет, в наших конкретных условиях. Полученное на практике p value определяют большое число различных факторов, таких как величина выборки, количество сравниваемых за один раз гипотез, появилась ли гипотеза до того, как мы посмотрели на данные, или после, принятие решений о прекращении данных с достижением первого значимого результата и тд.
  • Если наше p value < alpha – считаем, что вероятность получить ложноположительный результат мала, и отвергаем H0 и делаем вывод об истинности H1.
  • Если мы получили, что p value > alpha – увы, мы ничего не можем сказать про данные. Не факт, что нулевая гипотеза верна, возможно, просто наш тест оказался слабоват.

Assumptions

Чтобы применять тот или иной статистический критерий, нужно понимать, удовлетворяют ли данных условиям его применимости. Это также одна из причин симуляции данных: проверить, на что будут похожи наши экспериментальные данные, как мы потом будем их обрабатывать? Предположим, что у нас есть гипотеза, что данные будут распределены нормально, и мы хотим попробовать применить t-test. Для этого нужно проверить, удовлетворяют ли наши данные условиям его применимости, assumptions, для t-test это, в основном:

  • Количественные данные
  • Нормальность распределения данных
qqnorm(IQdata_real1, pch = 1, frame = FALSE)
qqline(IQdata_real1, col = "red", lwd = 2)

qqnorm(IQdata_real2, pch = 1, frame = FALSE)
qqline(IQdata_real2, col = "red", lwd = 2)

Будем считать, что эти данные выглядят похоже на нормальные.

Применение t-test

Теперь применим сам тест:

t.test(IQdata_real1, IQdata_real2)
## 
##  Welch Two Sample t-test
## 
## data:  IQdata_real1 and IQdata_real2
## t = -8.646, df = 56.823, p-value = 6.022e-12
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -35.06018 -21.87316
## sample estimates:
## mean of x mean of y 
##  89.36667 117.83333
t.test(IQdata_real1, IQdata_real2)$statistic
##         t 
## -8.645981
t.test(IQdata_real1, IQdata_real2)$p.value
## [1] 6.022186e-12

Как проинтерпретировать получившийся результат?

Нас здесь волнуют две цифры: t-значение и p-value. t-значение показывает, где в какую часть t-распределения разницы между средними наших выборок попала наша разница в средних, то есть насколько M2-M1 далеко от нуля. p-value говорит нам, можем ли мы сейчас, на именно сейчас получившихся данных отбросить нулевую гипотезу H0 об отсутствии различий или нет.

https://rpsychologist.com/d3/nhst/

Еще один важный параметр, без которого в 2020 нельзя интерпретировать данные – это размер эффекта. Для t-testa обычно в качестве метрики размера эффекта используется Cohen’s d: она показывает, насколько средние в наших выборках далеко отстоят друг от друга.

CohenD <- lsr::cohensD(IQdata_real1, IQdata_real2 )
CohenD
## [1] 2.232383

https://rpsychologist.com/d3/cohend/

А теперь давайте представим, что у нас есть возможность провести такой эксперимент 100 раз, где каждая попытка независима от других (симуляция же!). Больше всего нас интересует p-value.

# replicate(20, t.test(IQdata_real1, IQdata_real2, var.equal = F)$p.value)
replication <- replicate(100, t.test(round(rnorm(nreal, mean = M, sd = sd)), round(rnorm(nreal, mean = M2, sd = sd2)), var.equal = F)$p.value)
replication
##   [1] 4.899060e-05 1.330380e-03 3.748846e-03 7.162242e-05 1.751456e-05
##   [6] 4.471011e-05 1.031421e-04 2.060662e-03 1.189335e-06 4.388035e-04
##  [11] 9.637346e-04 2.036633e-02 5.292333e-04 1.861156e-03 4.888665e-05
##  [16] 4.223028e-03 1.268740e-06 8.364222e-06 9.034219e-06 8.619191e-04
##  [21] 1.610832e-04 5.348003e-03 5.945277e-05 2.483617e-07 5.234175e-05
##  [26] 3.107215e-04 3.883168e-04 1.884023e-07 7.566260e-03 1.648249e-03
##  [31] 1.523608e-03 1.965842e-04 3.945341e-06 1.203193e-03 4.344030e-05
##  [36] 2.208230e-03 2.104044e-03 1.383979e-06 4.189102e-06 3.773417e-05
##  [41] 1.685932e-07 1.205009e-05 7.656423e-06 7.764309e-06 1.651622e-06
##  [46] 6.136179e-04 3.573856e-03 8.707174e-06 9.589222e-06 3.357796e-05
##  [51] 6.903567e-04 3.383171e-03 2.324175e-05 6.784578e-05 9.037145e-04
##  [56] 1.528185e-03 8.701251e-07 2.059495e-04 4.411936e-03 5.916645e-03
##  [61] 2.159921e-04 3.277256e-04 2.444768e-05 1.384962e-04 3.344884e-06
##  [66] 7.500219e-05 3.596988e-05 1.181423e-04 4.427287e-03 4.027477e-04
##  [71] 1.213587e-03 5.278856e-04 3.988449e-03 2.408692e-06 7.561533e-07
##  [76] 1.035566e-03 1.190255e-03 7.004663e-05 2.304293e-04 8.864144e-06
##  [81] 1.613631e-03 2.631144e-04 1.812245e-04 1.323637e-03 2.481068e-03
##  [86] 1.647994e-04 7.375072e-04 1.311151e-03 1.211557e-05 8.217431e-04
##  [91] 2.459120e-05 9.994130e-03 1.528191e-03 2.632219e-03 9.038199e-07
##  [96] 4.784812e-08 3.195907e-05 3.844426e-03 2.241332e-08 5.060940e-07
hist(replication, breaks = 10)

Рассчитаем мощность, основываясь на полученном размере эффекта:

pwr::pwr.t.test(nreal, CohenD, sig.level = 0.05)
## 
##      Two-sample t test power calculation 
## 
##               n = 30
##               d = 2.232383
##       sig.level = 0.05
##           power = 1
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

Итог: как анализировать данные

Первый шаг – понять, что у вас за данные:

  1. Симулировать данные, чтобы посмотреть, как будут выгледть данные, и чем их можно будет обработать, или взять старые похожие данные (например, из открытых источников) и потренироваться на них.
  2. Визуализировать распределение данных
  3. Посчитать описательные статистики (среднее, медиана, стандартное отклонение от среднего и тд)
  4. Сравнить с нормальным: визуально, с помощью QQ-plot и, возможно, с помощью статистических тестов (например, Колмагоров-Смирнов)

Второй шаг – выбрать статистический тест. Для выбора теста надо уже понимать, что у вас за данные, и что вы хотите получить в результате теста. Подробнее об этом можно почитать, например, здесь https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3116565/

Третий шаг:

  1. Посчитать размер эффекта для полученных статистик
  2. Проинтерпретировать данные