permutation_tests

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.4     ✓ dplyr   1.0.2
## ✓ tidyr   1.1.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(modelr)
library(broom)
## 
## Attaching package: 'broom'
## The following object is masked from 'package:modelr':
## 
##     bootstrap
library(coin)
## Loading required package: survival
load("student_mat_full.RData")

Задачи на видео:

student_mat$absences %>% summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   4.000   5.709   8.000  75.000
student_mat %>% 
  group_by(help) %>% 
  summarise(m_absences = mean(absences))
## `summarise()` ungrouping output (override with `.groups` argument)
## # A tibble: 2 x 2
##   help         m_absences
##   <fct>             <dbl>
## 1 need_support       6.19
## 2 no_support         5.28

Quick plots to find a difference in

student_mat %>% 
  ggplot(aes(x = absences, fill= help))+
  geom_density(alpha = 0.5)

Вроде бы отличаются, но очень слабо

student_mat %>% 
  ggplot(aes(x = G1, fill= help))+
  geom_density(alpha = 0.5)

студенты из 2 групп отличаются сильно, хотя есть и те кто попал на пересечение

student_mat %>% 
  ggplot(aes(x = goout, fill= help))+
  geom_density(alpha = 0.5)

Разница уже очень слабо видна и не понятно, случайно ли она появилась

Давайте посмотрим на номинальные переменные и их доли

student_mat %>% 
  ggplot(aes(x = paid, fill= help))+
  geom_bar(position = "fill")

практически равное распределение

student_mat %>% 
  ggplot(aes(x = higher, fill= help))+
  geom_bar(position = "fill")

и кажется есть перекос в данных

Немного про деревья

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

Если мы наращиваем наше дерево и добавляем туда правило, которое мы выучили на данных, которые случайно немного отличаются

То наше правило не будет хорошо использоваться для предсказания на новых данных

И мы получим переобученную модель

Test

Статистические тесты

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

Последовательность шагов

  1. Определить нулевую гипотезу
  2. Посчитать статистику
  3. Вычислить p-value (ту самую вероятность получить эффект случайно)
  4. Интерпретировать результаты

Если p-value низок, то считается, что эффект статистически значим, т.е. маловероятно, что мы получили его случайно. Вывод: эффект, скорее всего, существует и в популяции, не только в выборке. (Важно: 1) все равно есть шанс случайности, 2) если наши данные со смещениями, то и результат будет со смещениями)

Этот процесс по принципу действия похож на доказательство от обратного. Сначала предполагаем, что эффекта (взаимосвязи) нет – это нулевая гипотеза. На следующем шаге мы считаем вероятность получить тот эффект (или сильнее), что мы наблюдаем в данных (p-value), случайно в предположении, что нулевая гипотеза верна. Если p-value низок (<0.05 как общепринятое значение), мы отвергаем нулевую гипотезу.

Hypothesis

H0: оценки G1 связаны с тем, успешно ли студент закончит год (help) H1: есть связь между оценкой и тем, успешно ли студент закончит год

number_of_permutations = 1000

pb = txtProgressBar(max = length(number_of_permutations), title = "Loop")

diff_random  <- rep(NA_real_, number_of_permutations) # вектор с неопределенными 

for (i in 1 : number_of_permutations) {
  
    # перемешиваем значения
    shuffled = student_mat %>% 
      select(G1) %>% 
      slice_sample(n =nrow(student_mat)) %>% 
      bind_cols(student_mat %>% select(help)) %>% 
      group_by(help) %>% 
      summarise(m = mean(G1), .groups = "drop")
    
    # diff_random[i] = shuffled %>% 
    #   summarise(diff = m-lag(m)) %>% 
    #   na.omit() %>% pull()
    
    # записываем разницу между теми кого случайно отнесли в одну или другую группу
    diff_random[i] = shuffled$m[1] - shuffled$m[2]
    
    # progress bar
    setTxtProgressBar(pb,i)
  }
## ================================================================================
close(pb)

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

df_r <- tibble(diff.random = diff_random)

df_r %>% 
ggplot() +
  geom_histogram(aes(diff.random), fill="cyan4")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

добавим туда разницу, что мы наблюдаем в настоящих данных

group_diff <- function(data, group, target){
  diff_real = data %>% 
  group_by({{group}}) %>% 
  summarise(m = mean({{target}}), .groups = "drop") %>% 
  summarise(diff = m - lag(m)) %>% 
  na.omit() %>% 
  pull()
  
  return(diff_real)
}


diff_real = student_mat %>% 
  group_diff(help, G1)

df_r %>% 
ggplot() +
  geom_histogram(aes(diff.random), fill="cyan4")+
  geom_vline(xintercept = diff_real)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

и доверительные интервалы, т.е. границы, куда попали 99.8% наших наблюдений

df_r %>% 
ggplot() +
  geom_histogram(aes(diff.random), fill="cyan4", binwidth = 0.1)+
  geom_vline(xintercept = diff_real)+
  geom_vline(xintercept = quantile(df_r$diff.random, 0.001), color = "red")+
  geom_vline(xintercept = quantile(df_r$diff.random, 0.999), color = "red")

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

# number_of_permutations = 1000

df_r %>% 
  mutate(abs_diff = abs(diff.random)) %>% 
  mutate(more = abs_diff >= abs(diff_real)) %>% 
  count(more) %>% 
  filter(more == TRUE) %>% 
  summarise(pvalue = n / number_of_permutations)
## # A tibble: 0 x 1
## # … with 1 variable: pvalue <dbl>

Broom mean from t-test

сделаем такой же тест но чуть более tidy

Используем permute, чтобы library(modelr) сделал за нас 100 перемешанных версий нашего датафрейма, а затем для каждого из них посчитаем среднее

perms = student_mat %>% 
  select(absences, help) %>% 
  modelr::permute(n = 500, help)

models <-  map(perms$perm, ~ t.test(absences ~ help, data = .)) 
glanced <- map_df(models, broom::glance, .id = "id")
# hist(glanced$statistic)

и нарисуем, чтобы сравнить с нашими настоящими

observed = broom::glance(t.test(absences ~ help, data = student_mat))$estimate
observed_2 = student_mat %>% group_diff(help, absences)

observed_2 + observed
## [1] 0
glanced %>% 
ggplot(aes(x = estimate)) +
  geom_histogram(fill = "cyan4") +
  geom_vline(xintercept = observed, color = "black") +
  geom_vline(aes(xintercept = quantile(estimate, 0.005)), color = "red")+
  geom_vline(aes(xintercept = quantile(estimate, 0.995)), color = "red")+
  theme_minimal()+
  ggtitle("absences and grades", "perm_test")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

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

tidy but in a different way

soure: https://mikedecr.github.io/post/randomization-inference-purrr/#fn3

observed_estimates <- student_mat %>%
  summarize(
    mean_control = mean(G1[help == "no_support"]),
    mean_treatment = mean(G1[help == "need_support"]),
    diff_means = mean_treatment - mean_control
  ) %>% print()
## # A tibble: 1 x 3
##   mean_control mean_treatment diff_means
##          <dbl>          <dbl>      <dbl>
## 1         13.3           8.26      -5.01

с помощью crossing мы переберём сделаем 3 датасета, а с помощью nest запишим их в ячейки нашего тибла

# ?crossing
student_mat %>%
  select(G1, help) %>% 
  mutate(id = row_number()) %>% 
  crossing(m = 1:3) %>%
  group_by(m) %>%
  nest()
## # A tibble: 3 x 2
## # Groups:   m [3]
##       m data              
##   <int> <list>            
## 1     1 <tibble [395 × 3]>
## 2     2 <tibble [395 × 3]>
## 3     3 <tibble [395 × 3]>

сперва сделаем много датасетов

# total number of iterations
n_sims <- 2000

# nested data frame of all M iterations
sim_table <- student_mat %>%
  select(G1, help) %>% 
  mutate(id = row_number()) %>% 
  crossing(m = 1:n_sims) %>%
  group_by(m) %>%
  nest()  %>% print()
## # A tibble: 2,000 x 2
## # Groups:   m [2,000]
##        m data              
##    <int> <list>            
##  1     1 <tibble [395 × 3]>
##  2     2 <tibble [395 × 3]>
##  3     3 <tibble [395 × 3]>
##  4     4 <tibble [395 × 3]>
##  5     5 <tibble [395 × 3]>
##  6     6 <tibble [395 × 3]>
##  7     7 <tibble [395 × 3]>
##  8     8 <tibble [395 × 3]>
##  9     9 <tibble [395 × 3]>
## 10    10 <tibble [395 × 3]>
## # … with 1,990 more rows

затем создадим новую колонку из датасетов, где перемешаем значения в колонке help

perm_table <- sim_table %>%
  mutate(
    permuted_data = map(
      .x = data, 
      .f = mutate, 
      help = sample(help)
    )
  ) %>% print()
## # A tibble: 2,000 x 3
## # Groups:   m [2,000]
##        m data               permuted_data     
##    <int> <list>             <list>            
##  1     1 <tibble [395 × 3]> <tibble [395 × 3]>
##  2     2 <tibble [395 × 3]> <tibble [395 × 3]>
##  3     3 <tibble [395 × 3]> <tibble [395 × 3]>
##  4     4 <tibble [395 × 3]> <tibble [395 × 3]>
##  5     5 <tibble [395 × 3]> <tibble [395 × 3]>
##  6     6 <tibble [395 × 3]> <tibble [395 × 3]>
##  7     7 <tibble [395 × 3]> <tibble [395 × 3]>
##  8     8 <tibble [395 × 3]> <tibble [395 × 3]>
##  9     9 <tibble [395 × 3]> <tibble [395 × 3]>
## 10    10 <tibble [395 × 3]> <tibble [395 × 3]>
## # … with 1,990 more rows

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

perm_table %>%
  ungroup() %>%
  sample_n(5) %>%
  unnest(permuted_data) %>%
  select(m, id, help) %>%
  pivot_wider(
    names_from = m, 
    values_from = help,
    names_prefix = "help, m = "
  )
## # A tibble: 395 x 6
##       id `help, m = 130` `help, m = 97` `help, m = 170` `help, m = 215`
##    <int> <fct>           <fct>          <fct>           <fct>          
##  1   249 no_support      need_support   no_support      need_support   
##  2   138 no_support      no_support     no_support      no_support     
##  3     1 no_support      need_support   no_support      no_support     
##  4     2 no_support      need_support   need_support    no_support     
##  5    80 need_support    no_support     no_support      need_support   
##  6   145 no_support      no_support     need_support    need_support   
##  7   154 no_support      need_support   no_support      no_support     
##  8   162 no_support      no_support     no_support      no_support     
##  9   165 no_support      need_support   need_support    no_support     
## 10     5 no_support      need_support   need_support    no_support     
## # … with 385 more rows, and 1 more variable: `help, m = 720` <fct>

ну и посчитаем для каждого из наших переменных датасетов разницу в среднем среди тех кто нуждается в помощи и не нуждается

est_table <- perm_table %>%
  mutate(
    estimates = map(
      .x = permuted_data,
      .f = summarize, 
      mean_treatment_m = mean(G1[help == "need_support"]),
      mean_control_m = mean(G1[help == "no_support"]),
      diff_means_m = mean_treatment_m - mean_control_m
    )
  ) %>% print()
## # A tibble: 2,000 x 4
## # Groups:   m [2,000]
##        m data               permuted_data      estimates       
##    <int> <list>             <list>             <list>          
##  1     1 <tibble [395 × 3]> <tibble [395 × 3]> <tibble [1 × 3]>
##  2     2 <tibble [395 × 3]> <tibble [395 × 3]> <tibble [1 × 3]>
##  3     3 <tibble [395 × 3]> <tibble [395 × 3]> <tibble [1 × 3]>
##  4     4 <tibble [395 × 3]> <tibble [395 × 3]> <tibble [1 × 3]>
##  5     5 <tibble [395 × 3]> <tibble [395 × 3]> <tibble [1 × 3]>
##  6     6 <tibble [395 × 3]> <tibble [395 × 3]> <tibble [1 × 3]>
##  7     7 <tibble [395 × 3]> <tibble [395 × 3]> <tibble [1 × 3]>
##  8     8 <tibble [395 × 3]> <tibble [395 × 3]> <tibble [1 × 3]>
##  9     9 <tibble [395 × 3]> <tibble [395 × 3]> <tibble [1 × 3]>
## 10    10 <tibble [395 × 3]> <tibble [395 × 3]> <tibble [1 × 3]>
## # … with 1,990 more rows

вытащим датасеты из ячеек и переделаем в привычный вид с колонками

est_summary <- est_table %>%
  select(m, estimates) %>%
  unnest(cols = estimates) %>%
  ungroup() %>% print()
## # A tibble: 2,000 x 4
##        m mean_treatment_m mean_control_m diff_means_m
##    <int>            <dbl>          <dbl>        <dbl>
##  1     1             11.0           10.8        0.182
##  2     2             10.8           11         -0.194
##  3     3             11.0           10.8        0.203
##  4     4             11.2           10.7        0.528
##  5     5             10.7           11.1       -0.448
##  6     6             10.6           11.2       -0.539
##  7     7             10.8           11.0       -0.224
##  8     8             10.7           11.1       -0.356
##  9     9             10.8           11.0       -0.163
## 10    10             11.2           10.7        0.518
## # … with 1,990 more rows

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

est_summary %>% 
  ggplot(aes(x = diff_means_m))+
  geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

и посчитаем p-value

pval <- est_summary %>%
  summarize(
    p.value = mean(abs(diff_means_m) > abs(observed_estimates$diff_means))
  ) %>%
  pull(p.value) %>%
  print()
## [1] 0

coin

И посмотрим как этот тест реализован в пакете library(coin)

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

а затем Y – переменная, которой мы эту связь пытаемся объяснить

independence_test(G1 ~ help, data = student_mat)
## 
##  Asymptotic General Independence Test
## 
## data:  G1 by help (need_support, no_support)
## Z = -14.974, p-value < 2.2e-16
## alternative hypothesis: two.sided

Z statistics – это мерика , которая показывает, насколько реальные данные далеко от случайных

в нашем случае мы на 14.9 стандартных отклонений удалены от среднего случайного значения.

p-value эта наша вероятность ошибки

independence_test(absences ~ help, data = student_mat)
## 
##  Asymptotic General Independence Test
## 
## data:  absences by help (need_support, no_support)
## Z = 1.1229, p-value = 0.2615
## alternative hypothesis: two.sided
independence_test(higher ~ help, data = student_mat)
## 
##  Asymptotic General Independence Test
## 
## data:  higher by help (need_support, no_support)
## Z = 3.0225, p-value = 0.002507
## alternative hypothesis: two.sided
independence_test(paid ~ help, data = student_mat)
## 
##  Asymptotic General Independence Test
## 
## data:  paid by help (need_support, no_support)
## Z = 0.85477, p-value = 0.3927
## alternative hypothesis: two.sided

end

Задачи на видео:

  • Мы не хотим строить модель на основании зависимости, которая могла появиться случайно
  • в тесте перестановок мы сравниваем настоящее распроделение переменной со множеством случайных распределений
  • independence_test(x ~ y), где x это переменная которую мы пытаемся объяснить предсказать, а у эта та переменная, которая должна с помощью которой мы пытаемся объяснить