Как ответить?

Кто больше времени уделяет спорту:

женщины или мужчины?

Как ответить?

Способ 1: выяснить у всех женщин и мужчин, насколько активный образ жизни они ведут

Как ответить?

Способ 2: набрать репрезентативную выборку и выяснить их спортивную активность

  • Сколько людей набрать?

  • Кого брать?

  • Как вообще выбирать людей?

  • Как оценивать спортивную активность?

Генеральная совокупность и выборка

Хотим получить знание о генеральной совокупности, т.е. о всех возможных участниках, которых касается наш опрос или эксперимент

Исследовать генеральную совокупность долго, дорого и часто просто невозможно

Решение: набрать такую выборку, чтобы по результатам исследования этой выборки можно было бы судить о всей генеральной совокупности

Генеральная совокупность

При планировании исследования важно понимать, какие именно объекты составляют для этого исследования генеральную совокупность.

От этого зависит стратегия составления репрезентативной выборки

Способы создания репрезентативной выборки

  • Простая случайная выборка

  • Кластерная случайная выборка

  • Стратифицированная случайная выборка

  • Систематическая выборка

  • …

Простая случайная выборка

Каждый элемент генеральной совокупности имеет равный шанс попасть в выборку

Кластерная случайная выборка

  • Разбиваем генеральную совокупность на кластеры

  • Случайным образом отбираем несколько кластеров

  • Включаем ВСЕХ участников этого кластера

  • Больше подходит для однородной генеральной совокупности

Стратифицированная случайная выборка

  • Разбиваем генеральную совокупность на кластеры

  • Случайным образом из каждого кластера отбираем случайных участников

  • Больше подходит для неоднородной генеральной совокупности

Систематическая случайная выборка

  • Упорядочиваем всех членов генеральной совокупности в некотором порядке (по алфавиту, по возрасту, …)

  • Выбираем случайным образом начальную точку

  • Отбираем каждого n-ного участника

Невероятностные способы набора выборок

  • “Удобная” выборка

Стоим у дверей факультета и опрашиваем прохожих

  • Метод добровольного опроса

Просим подписчиков заполнить анкету

  • Снежный ком

Отбираем несколько человек для начала исследования и просим их поспособствовать дальнейшему набору

  • Целевой отбор

Решаем, кто нам нужен, и набираем таких людей

Систематическая ошибка отбора

  • Исследовать только интервал, укладывающийся в теорию

  • Ошибка выжившего

  • Некорректная “удобная” выборка

  • Ошибка меткого стрелка

Систематическая ошибка неответа

Участники, попавшие в выборку, значительно отличаются от непопавших.

В итоге выборка не репрезентативна.

Пример: В кофейне решили выяснить, кто больше предпочитает новый сорт кофе: “совы” или “жаворонки”. Опрос устроили в 8 утра.

И еще много других ошибок…

Переменные

Влияет ли занятия спортом на настроение?

  • Независимая переменная: наличие занятий спортом (количество раз в неделю)

  • Зависимая переменная: оценка настроения (по какой-то шкале)

  • Конфаундеры: влияют на значения зависимой переменной, но не контролируются в эксперименте (время года, возраст и пол участников, масса тела, сфера деятельности, ….)

Конфаундеры

При планировании исследования нужно тщательно продумать, что может влиять на результат

По возможности устранить или снизить эти влияния

Ослепление

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

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

  • Слепое исследование: испытуемый не знает, в какой он группе

  • Двойное слепое исследование: ни испытуемый, ни экспериментатор (врач) не знают, кто в какой группе

  • Тройное слепое исследование: даже тот, кто обрабатывает данные, не знает, какая группа получала плацебо, а какая лекарство

Типы данных

Визуализация категориальных данных

library(tidyverse)
iris %>%
  sample_n(60) %>%
  ggplot(aes(x = Species, 
             fill = Species)) +
  geom_bar() + theme_bw() +
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=14,face="bold"))

Для порядковых категориальных данных важен порядок столбцов

Визуализация непрерывных данных

iris %>%
  ggplot(aes(x = Sepal.Length)) + 
  geom_histogram(fill = 'orange', color = 'brown') + theme_bw() +
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=14,face="bold"))

Гистограмма

Многое зависит от ширины бина

iris %>%
  ggplot(aes(x = Sepal.Length)) + 
  geom_histogram(fill = 'orange', color = 'brown', binwidth = 1) + theme_bw() +
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=14,face="bold"))

Статистика

  • Описательная статистика

  • Индуктивная статистика

Описательная статистика

Любой набор данных нужно сначала изучить

  • список переменных

  • количество наблюдений

  • наличие пропущенных значений

  • распределения численных переменных

  • композиция категориальных переменных

  • наличие выбросов

  • визуализация

Описательные статистики

Функция от выборки

Характеризуем выборку одним числом

  • медиана

  • среднее

  • мода

  • дисперсия

  • разброс

  • …

Выборочное среднее

\[\overline{X} = {x_{1} + x_{2} + ... + x_{n} \over n}\]

set.seed(123)
x <- rpois(6, 3)
x <- sort(x)
x
[1] 0 2 2 4 5 6
mean(x)
[1] 3.166667

Медиана

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

x = c(1, 4, 8, 9, 11)
x
[1]  1  4  8  9 11
median(x)
[1] 8
x = c(1, 4, 8, 9, 11, 13)
x
[1]  1  4  8  9 11 13
median(x)
[1] 8.5

Процентили

N-ный перцентиль - это такое число, что N% элементов выборки не больше N

a <- rpois(30, 2)
quantile(a)
  0%  25%  50%  75% 100% 
0.00 2.00 2.00 3.75 6.00

  • 25-ый перцентиль - 1ый квартиль

  • 75-ый перцентиль - 3ий квартиль

  • 50-ый перцентиль - медиана

Для описания большинства данных можно использовать 95-ый перцентиль

Мода

Наиболее частое значение

set.seed(123)
a <- rpois(20, 2)
sort(a)
 [1] 0 0 0 1 1 1 2 2 2 2 2 2 3 3 4 4 4 4 5 5

Мода

library(tidyverse)
ggplot(mapping = aes(x = a)) + geom_dotplot() + theme_bw()

Мод может быть несколько

Разброс

Минимум и максимум

sort(a)
 [1] 0 0 0 1 1 1 2 2 2 2 2 2 3 3 4 4 4 4 5 5
range(a)
[1] 0 5

Боксплот - ящик с усами

Представление числовых данных через квартили

Боксплот

quantile(iris$Sepal.Length)
  0%  25%  50%  75% 100% 
 4.3  5.1  5.8  6.4  7.9 
ggplot(iris) +
  geom_boxplot(aes(Sepal.Length)) + theme_bw()

Для одной выборки не очень подходит

Боксплот

ggplot(iris) +
  geom_boxplot(aes(Sepal.Length, fill = Species)) + theme_bw()

Интерквартильный размах

ggplot(iris) +
  geom_boxplot(aes(Sepal.Length, fill = Species)) + theme_bw()

summary(iris[iris$Species == 'setosa',]$Sepal.Length)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  4.300   4.800   5.000   5.006   5.200   5.800 
IQR(iris[iris$Species == 'setosa',]$Sepal.Length)
[1] 0.4

Выбросы

Есть много методов поиска выбросов

Выброс должен быть объясним с точки зрения структуры эксперимента или исследования

Выброс - это ошибка или интересное наблюдение?

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

Данные про зубы

library(ggpubr)
data("ToothGrowth")
head(ToothGrowth)
   len supp dose
1  4.2   VC  0.5
2 11.5   VC  0.5
3  7.3   VC  0.5
4  5.8   VC  0.5
5  6.4   VC  0.5
6 10.0   VC  0.5

Визуализируем среднее

ggerrorplot(ToothGrowth, x = "dose", y = "len", 
            desc_stat = "mean_sd",
            error.plot = "errorbar",            
            add = "mean"                        
            )

Визуализируем среднее

ggerrorplot(ToothGrowth, x = "dose", y = "len", 
            desc_stat = "mean_sd", color = "black",
            add = "jitter", add.params = list(color = "darkgray")
            )

Визуализируем среднее

ggerrorplot(ToothGrowth, x = "dose", y = "len", 
            desc_stat = "mean_sd", 
            color = "dose", palette = "jco")

Визуализируем среднее

ggerrorplot(ToothGrowth, x = "dose", y = "len", 
            desc_stat = "mean_sd", 
            color = "supp", palette = "jco",
            position = position_dodge(0.3))

Визуализируем среднее

ggbarplot(ToothGrowth, x = "dose", y = "len", 
       add = c("mean_se", "jitter"))

Визуализируем среднее

ggbarplot(ToothGrowth, x = "dose", y = "len", 
          add = c("mean_se", "jitter"),
          color = "supp", palette = "jco",
          position = position_dodge(0.8))

Barplot

Обязательно указывайте, что именно отображает ваш барплот:

  • количество наблюдений в группе

  • количество каких-то конкретных значений в группе

  • среднее значение параметра в группе

И обязательно указывайте, что именно отображают “усы”

Еще раз про NA

Есть ли хотя бы одно пропущенное значение?

anyNA(airquality)
[1] TRUE

Сколько всего NA?

colSums(is.na(airquality))
  Ozone Solar.R    Wind    Temp   Month     Day 
     37       7       0       0       0       0

В каких столбцах и в каких колонках находятся пропущенные значения?

Как вы это узнали?

А если столбцов 1000?

Если вы что-то оцениваете методом глаза, то вы делаете это неправильно

Задание №1

Смоделируйте скошенное распределение любым способом (положите в вектор)

Визуализируйте распределение, выбрав наиболее оптимальный тип графика

Решение 1.1

set.seed(123)
a <- c(rnorm(10000), rnorm(500, 4, 3), rnorm(500, 6, 5))

Решение 1.2

ggplot(mapping = aes(x = a)) + 
  geom_histogram(fill = '#598392', color = '#01161e') + 
  theme_bw()

Усеченное среднее

mean(a)
[1] 0.4604595
mean(a, trim = .2)
[1] 0.1117887

Что такое среднее с параметром trim = 0.5?

Асимметрия и эксцесс

skew(a)
[1] 3.24231
kurtosi(a)
[1] 15.81482

Асимметрия и эксцесс: примеры

Сводное описание

summary(a)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-7.9186 -0.6191  0.0928  0.4605  0.8825 20.0108 
describe(a) # library(psych)
   vars     n mean   sd median trimmed  mad   min   max range skew kurtosis
X1    1 11000 0.46 2.18   0.09    0.13 1.11 -7.92 20.01 27.93 3.24    15.81
     se
X1 0.02

describe + group_by/summarise

iris %>%
  group_by(Species) %>%
  summarise(describe(Petal.Width))
# A tibble: 3 × 14
  Species  vars     n  mean    sd median trimmed   mad   min   max range    skew
  <fct>   <dbl> <dbl> <dbl> <dbl>  <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl>
1 setosa      1    50 0.246 0.105    0.2   0.238 0       0.1   0.6   0.5  1.18  
2 versic…     1    50 1.33  0.198    1.3   1.32  0.222   1     1.8   0.8 -0.0293
3 virgin…     1    50 2.03  0.275    2     2.03  0.297   1.4   2.5   1.1 -0.122 
# ℹ 2 more variables: kurtosis <dbl>, se <dbl>

library(skimr)

Data summary
Name iris
Number of rows 150
Number of columns 5
_______________________
Column type frequency:
factor 1
numeric 4
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
Species 0 1 FALSE 3 set: 50, ver: 50, vir: 50

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
Sepal.Length 0 1 5.84 0.83 4.3 5.1 5.80 6.4 7.9 ▆▇▇▅▂
Sepal.Width 0 1 3.06 0.44 2.0 2.8 3.00 3.3 4.4 ▁▆▇▂▁
Petal.Length 0 1 3.76 1.77 1.0 1.6 4.35 5.1 6.9 ▇▁▆▇▂
Petal.Width 0 1 1.20 0.76 0.1 0.3 1.30 1.8 2.5 ▇▁▇▅▃

qqplot

set.seed(1)
my_df <- data.frame(y=rexp(200, rate=3))
ggplot(my_df, mapping = aes(x = y)) + geom_histogram(fill = '#598392', color = '#01161e') + theme_bw()

qqplot

ggplot(my_df, aes(sample=y)) +
  stat_qq() + 
  stat_qq_line() + theme_bw()

qqplot

Задание 2

Отличается ли в среднем длина лепестков у ирисов от 4 см?

Сначала рисуйте!

Нарисуем!

ggplot(iris) +
  geom_histogram(aes(x = Petal.Length, fill = Species)) + theme_bw()

Какие идеи?

Подготовим данные

df <- iris %>%
  filter(Species %in% c('versicolor', 'virginica'))

Решение 2.1

t.test(df$Petal.Length)

    One Sample t-test

data:  df$Petal.Length
t = 59.425, df = 99, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 4.742187 5.069813
sample estimates:
mean of x 
    4.906

Что не так?

Какая нулевая гипотеза?

Решение 2.2

t.test(df$Petal.Length, mu = 4)

    One Sample t-test

data:  df$Petal.Length
t = 10.974, df = 99, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 4
95 percent confidence interval:
 4.742187 5.069813
sample estimates:
mean of x 
    4.906

Без интерпретации результата - ответа нет!

Неверная интерпретация - ответа нет!

Вывод результата теста

t.test(df$Petal.Length, mu = 4) -> x
x$p.value
[1] 8.307922e-19

Задание 3

Больше ли в среднем длина лепестков у ирисов, чем 3 см?

Какая нулевая гипотеза?

Решение 3.1

t.test(df$Petal.Length, mu = 3, alternative = 'greater')

    One Sample t-test

data:  df$Petal.Length
t = 23.087, df = 99, p-value < 2.2e-16
alternative hypothesis: true mean is greater than 3
95 percent confidence interval:
 4.768922      Inf
sample estimates:
mean of x 
    4.906

Задание 4

Различаются ли в среднем длины лепестков у ирисов видов versicolor и virginica?

Какая нулевая гипотеза?

Решение 4.1

ggplot(df) + geom_boxplot(aes(Sepal.Length, fill = Species)) + theme_bw()

Теперь визуализируйте правильно

Решение 4.2

ggerrorplot(df, x = "Species", y = "Petal.Length", 
            desc_stat = "mean_sd", color = "black",
            add = "jitter", add.params = list(color = "darkgray"), size = 0.8
            )

Решение 4.3

a <- filter(df, Species == 'versicolor')
b <- filter(df, Species == 'virginica')

t.test(a$Petal.Length, b$Petal.Length)

    Welch Two Sample t-test

data:  a$Petal.Length and b$Petal.Length
t = -12.604, df = 95.57, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.49549 -1.08851
sample estimates:
mean of x mean of y 
    4.260     5.552

Решение 4.4

Лучше

t.test(Petal.Length ~ Species, df)

    Welch Two Sample t-test

data:  Petal.Length by Species
t = -12.604, df = 95.57, p-value < 2.2e-16
alternative hypothesis: true difference in means between group versicolor and group virginica is not equal to 0
95 percent confidence interval:
 -1.49549 -1.08851
sample estimates:
mean in group versicolor  mean in group virginica 
                   4.260                    5.552

Решение 4.5

Теперь визуализируем еще лучше

ggbarplot(df, x = "Species", y = "Petal.Length", 
          fill = ("Species"), 
          add = c('mean', 'jitter'),
          palette = c("#9e2a2b", "#62929e")) +
  stat_compare_means(method = "t.test")

Задание 5

Пусть виды ирисов будут “до” обработки супер-средством и “после”.

Влияет ли супер-средство на длину лепестков ирисов?

Решение 5.1

t.test(Pair(Petal.Length, Species) ~ 1, df)

    Paired t-test

data:  Pair(Petal.Length, Species)
t = 45.341, df = 99, p-value < 2.2e-16
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 2.300709 2.511291
sample estimates:
mean difference 
          2.406

Решение 5.2

Визуализируем

ggpaired(df, x = "Species", y = "Petal.Length", 
          fill = ("Species"), 
          add = c('mean', 'jitter'),
          palette = c("#9e2a2b", "#62929e"),
         line.color = "gray", line.size = 0.4) +
  stat_compare_means(method = "t.test")

t.test(Pair(Petal.Length, Species) ~ 1, df)$p.value -> a

ggpaired(df, x = "Species", y = "Petal.Length", 
          fill = ("Species"), 
          add = c('mean', 'jitter'),
          palette = c("#9e2a2b", "#62929e"),
         line.color = "gray", line.size = 0.4) + 
annotate("text", x=1.3, y=7.5, 
         label = paste('p.value (t.test): ', a, sep = ''), size = 4)

immer

Урожайность ячменя в 1931(Y1) и 1932(Y2) годах.

Задание 6

Различается ли урожайность?

Смотрим

library(MASS)
head(immer)
  Loc Var    Y1    Y2
1  UF   M  81.0  80.7
2  UF   S 105.4  82.3
3  UF   V 119.7  80.4
4  UF   T 109.7  87.2
5  UF   P  98.3  84.2
6   W   M 146.6 100.4

Поправим данные

im <- immer %>%
  pivot_longer(Y1:Y2, names_to = 'year', values_to = 'value')

Рисуем

Рисуем

im %>% 
  ggplot() + geom_boxplot(aes(value, fill = year)) + theme_bw()

Еще рисуем

im %>%
  ggplot(aes(x = year, y = value)) + 
  geom_bar(position='dodge', stat='summary', fun='mean', 
                      fill = 'lightblue', color = 'brown') + 
  geom_jitter(width = 0.25) + theme_bw()

Решение 6.1

im_wid <- im %>% 
  pivot_wider(names_from = year)
wilcox.test(Pair(Y1, Y2) ~ 1, im_wid)

    Wilcoxon signed rank test with continuity correction

data:  Pair(Y1, Y2)
V = 368.5, p-value = 0.005318
alternative hypothesis: true location shift is not equal to 0

Задание 7

Есть ли зависимость между длиной лепестка и длиной чашелистика у ирисов?

Рисуем

Рисуем

ggplot(data = iris) + 
  geom_point(aes(x = Sepal.Length, y = Petal.Length)) + theme_bw()
cor(iris$Sepal.Length, iris$Petal.Length)
[1] 0.8717538

Пирсона или Спирмена?

cor(iris$Sepal.Length, iris$Petal.Length)
[1] 0.8717538
cor(iris$Sepal.Length, iris$Petal.Length, method = 'spearman')
[1] 0.8818981

А значимость?

cor.test(iris$Sepal.Length, iris$Petal.Length)

    Pearson's product-moment correlation

data:  iris$Sepal.Length and iris$Petal.Length
t = 21.646, df = 148, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.8270363 0.9055080
sample estimates:
      cor 
0.8717538 
cor.test(iris$Sepal.Length, iris$Petal.Length, method = 'spearman')

    Spearman's rank correlation rho

data:  iris$Sepal.Length and iris$Petal.Length
S = 66429, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.8818981

Конец!