06_Bootstrap.R

# Мы рассмотрели классические основы статистики и сегодня коснёмся 
# одного из более современных способов оценки параметров и проверки гипотез

# Автор: Алексей Замалутдинов

#-------------------------------------------------------------------------
library(dplyr)
library(ggplot2)
library(boot)
main_dir <- dirname(rstudioapi::getSourceEditorContext()$path) 
setwd(main_dir)
library(devtools)
devtools::install_github("ProcessMiner/nlcor")
library(nlcor)
#-------------------------------------------------------------------------

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

# На самом деле, нам необязательно пытаться подстроить наши данные под то или иное типовое распределение,
# ведь наши данные сами являются распределением. 
# Построим график кумулятивной вероятности (от 0 до 1)

x <- rnorm(50000)
x1 <- rnorm(20)
x4 <- rnorm(200)
x <- as.data.frame(x)
x1 <- as.data.frame(x1)
x4 <- as.data.frame(x4)
ggplot() + 
  stat_ecdf(data = x, aes(x)) + 
  stat_ecdf(data = x1, aes(x1),colour = "red") + 
  stat_ecdf(data = x4, aes(x4),colour = "blue") 

stat_fun <- function(x, idx) mean(x[idx])
boot_obj <- boot(x4$x4, R=10000, statistic=stat_fun)
boot_obj
mean(x4$x4)
sd(x4$x4)
mean(boot_obj$t)
# Что мы здесь видим? Чем больше наша выборка, тем ближе она приближается к кривой вероятности, 
# созданной исходным распределением (нормальным в данном случае). Таким образом, нам не нужно 
# определять исходное распределение и определять его параметры. Такая функция распределения называется
# эмпирической функцией распределения.
# 
# Если наша функция распределения довольно хорошо приближается к 
# оригинальному распределению то, мы можем использовать её для нашего анализа.
# 
# Вернёмся к вопросу о медиане, как оценить её доверительный интервал? 
# Было бы круто, если бы у нас было много таких выборок, 
# мы бы посчитали медиану для каждой и рассчитали бы среднее и стандратное отклонение.
# Но у нас всего одна, что делать?
#
# Сгенерировать новые! На основе нашей выборки создадим 1000 новых выборок с повторением.
#
(a <- sample(filter(iris, Species == "setosa")$Sepal.Length, size = 50, replace=T))
median(a)

stat_fun <- function(x, idx) median(x[idx])
boot_obj <- boot(filter(iris, Species == "setosa")$Sepal.Length, R=1000, statistic=stat_fun)
boot_obj
median(filter(iris, Species == "setosa")$Sepal.Length) # медиана выборки равна 5

# Полученная 1000 медиан не очень нормально распределена
t <- as.data.frame(boot_obj$t)
ggplot(t, aes(V1)) +
  geom_density()
shapiro.test(boot_obj$t)


boot.ci(boot_obj, R=1000, statistic=stat_fun)
# Здесь для нас сразу расчитывается несколько версий доверительных интервалов, рассмотрим 2 из них
# Если наши оценки медиан распределены нормально, мы можем расчитать 
# доверительный интервал с критическими значениями нормального распределения 
# В данном случае делается ещё поправка на смещение величины
# sample_median - bias +- se * qnorm(0.975)

# Верхняя граница тогда будет такой, как и расчитано функцией в типе Normal
5 - 0.01225 + 0.04743087 * qnorm(0.975)
# Нижняя
5 - 0.01225 - 0.04743087 * qnorm(0.975)


# Другой интересный способ оценки, это перцентиль. 
# В случае, если оценка получилась ненормально распределённой, то мы можем задать доверительный интервал, 
# выбросив самые экстремальные значения
# В нашем случае это 2.5% оценок с одной стороны и 2.5% с другой
stat_fun <- function(x, idx) median(x[idx])
boot_diamonds <- boot(slice_sample(diamonds, n = 1000)$price, R=1000, statistic=stat_fun)
boot_diamonds
t <- as.data.frame(boot_diamonds$t)
ggplot(t, aes(V1)) +
  geom_density()
boot.ci(boot_diamonds, R=1000, statistic=stat_fun)

# А ещё мы можем считать доверительные интервалы не только для параметров зависящих от одной переменной, 
# но и более, например, корреляцию. Да и сразу несколько, чтобы получилось быстрее:)

foo <- function(data, indices){
  dt<-data[indices,]
  c(
    cor(dt[,1], dt[,2], method='s'),
    median(dt[,1])
  )
}
boot_2_metrics <- boot(filter(iris, Species == "setosa"), R=1000, statistic=foo)
boot_2_metrics
t <- as.data.frame(boot_2_metrics$t)
ggplot(t, aes(V1)) +
  geom_density()
# Как раз случай, когда метрика распределена не очень нормально

boot.ci(boot_2_metrics, R=1000, statistic=foo, index=1)
#
#-------------------------------------------------------------------------
# Более того, бустреп позволяет нам проверять гипотезы. Если у нас есть доверительные интервалы,
# а значит и есть критические значения. Всё, что за них выходит является экстремальными значениями.
# Таким образом, у нас есть основания принять или отвергнуть наши предположения (гипотезы)

# Пример:

# Когда-то мы уже проверяли гипотезу о равенстве средних ширин лепестков двух популяций ирисов
# Давайте проверим это снова:
setosa <- filter(iris, Species == "setosa")$Sepal.Width
virginica <- filter(iris, Species == "virginica")$Sepal.Width
t.test(setosa,virginica) # различаются по т тесту

# Теперь сделаем проверку бутстрепом
mean_difference <- function(data, indices){
  dt<-data[indices,]
  mean_1 <- mean(dt[,1])
  mean_2 <- mean(dt[,2])
  mean_1 - mean_2
}
data_iris <- data.frame(setosa,virginica)
boot_test <- boot(data_iris, R=1000, statistic=mean_difference)
boot_test
t <- as.data.frame(boot_test$t)
ggplot(t, aes(V1)) +
  geom_density()

boot.ci(boot_test, R=1000, statistic=mean_difference)
# Так, что мы видим
# Получается, что 95% доверительный интервал для разницы средних у нас от 0.3136 до 0.5983
# Мы же в нулевой гипотезе предполагаем, что разницы между средними нет, равна 0
# То есть 0 не находится в нашем доверительном интервале, а значит мы можем отвергнуть нулевую гипотезу
pnorm(0,0.454+-0.001712,0.07002185)*2 # умножаем на 2, чтобы учесть двухстороннесть оценки
# 
# Прелесть бутстрапа в том, что мы можем сравнивать разные параметры (средние, медианы, корреляции и тд),
# не особо заботясь о распределении наших данных. 

# P.S.
# Важно понимать, что происходит внутри тулов

#-------------------------------------------------------------------------

# Но это не единственный "неклассический" способ проверки гипотез, мы можем использовать ещё пермутацию
# Идея довольно логичная. Если мы предполагаем, что наборы данных не отличаются, то и в случае
# их перемешивания между собой, они тоже не будут отличаться.
# Алгоритм:
# 1) Определить нулевую гипотезу
# 2) Расчитать параметр для исходных данных (например, разницу между средними)
# 3) Объединить данные и случайно набрать группы заново (без повторений)
# 4) Расчитать параметр для новых наборов
# 5) Повторить п. 3-4 много раз (10000)
# 6) Сранить реальное значение с полученным распределением
# Вернёмся к нашим ирисам и попробуем на них
setosa <- filter(iris, Species == "setosa")$Sepal.Width
virginica <- filter(iris, Species == "virginica")$Sepal.Width
# Определим разницу
real_diff <- mean(setosa) - mean(virginica)
real_diff
# Объединим данные
setosa_virginica <- c(setosa,virginica)
# Зададим функцию, которая будет перемешивать наши данные и считать разницу
perm_fun <- function(x, nA, nB) {
  n <- nA+nB
  idx_b <- sample(1:n, nB)
  idx_a <- setdiff(1:n, idx_b)
  mean_diff <- mean(x[idx_b]) - mean(x[idx_a]) 
  return(mean_diff)
}
# Проверим её
perm_fun(setosa_virginica,50,50)
# Теперь повторим 10000
i = 1
result <- c()
for (i in seq(1:10000000)) {
  result[i] <- perm_fun(setosa_virginica,50,50)
}
result <- as.data.frame(result)
ggplot(result, aes(result)) + 
  geom_density()
# А это наше p value. Вероятность получить группы более разные чем у нас есть.
mean(result > real_diff)
# Как и всегда, есть более автоматизированный вариант анализа
iris_2 <- filter(iris, Species != "versicolor")
library(coin)
oneway_test(Sepal.Width ~ Species, data = iris_2,distribution=approximate(nresample=10000))
#
#-------------------------------------------------------------------------
#
# Ещё примеры и описания
# https://habr.com/ru/companies/X5Tech/articles/679842/
# https://www.datacamp.com/tutorial/bootstrap-r
# https://www.r-bloggers.com/2019/09/understanding-bootstrap-confidence-interval-output-from-the-r-boot-package/
# https://www.dataanalysisclassroom.com/lesson98/
#-------------------------------------------------------------------------
# Поговорим ещё немного о корреляциях
# Мы обсуждали в самом начале корреляцию Пирсона
# Однако, она хорошо приспособлена для работы с нормально распределёнными данными.
# А что если это не так? У нас порядковые величины или ненормальное распределение
# Можно применить непараметрические методы расчёта корреляции, например, Спирмена или Кендалла
# Идейно они очень похожи и используют ранги для расчёта коэффициента
# Возьмём для примера не очень нормальные данные mtcars
data <- mtcars
ggplot(data, aes(disp))+
  geom_density()
ggplot(data, aes(hp))+
  geom_density()
# Рассчитаем классическую корреляцию Пирсона
cor(data$disp,data$hp)
# А теперь корреляцию Спирмена
# Строго говоря, она рассчитывается также как и Пирсона, но вместо значений мы берём ранги 
# Cov(X,Y) / (std(x)*std(y))
# Где Cov(X,Y)=sum((x-mean(x))*(y - mean(y)))/(n-1)
rank_disp <- rank(data$disp)
rank_hp <- rank(data$hp)

sum((rank_disp-mean(rank_disp))*(rank_hp - mean(rank_hp)))/(length(rank_disp)-1)/(sd(rank_disp)*sd(rank_hp))
# Проверим расчёты
cor(rank_disp,rank_hp)
cor(data$disp,data$hp,method = "spearman")
# Практически везде фигурирует такая формула
# 1 - 6*sum(d**2)/(n*(n**2-1))
# d = rank_1 - rank_2
# Она возможна при допущении, что все ранги целые числа (не всегда верно)
# Можно посмотреть вывод тут
# https://en.wikipedia.org/wiki/Spearman%27s_rank_correlation_coefficient
1 - 6*sum((rank_disp - rank_hp)**2)/(length(rank_disp)*(length(rank_disp)**2-1))
# Но результат не совсем верный, так как у нас есть дробные ранги
#
#-------------------------------------------------------------------------
# Тем не менее, часто наши данные имеют нелинейную зависимость
library(nlcor)
# Суть работы данной функции заключается в локальном вычислении корреляции 
# https://github.com/ProcessMiner/nlcor?ysclid=lnbqzupxge497670557
x <- rnorm(100)
x2 <- x**2
plot(x,x2)
cor(x,x2) # выглядит не очень
result <- nlcor(x, x2,plt = T)
result$cor.estimate # теперь намного лучше
plot(result$cor.plot)
#