ТО САМОЕ ДЗ

Задание №1

Выберите свою любую любимую картинку в рамках приличия.

Используя пакет hexSticker, нарисуйте шестиугольный логотип

#в рамках приличия
sticker("./1-126-512.png", package = '',
        p_size=20, s_x=1, s_y=0.95, s_width=0.8, s_height=0.8,
        filename="img/sticker_hex.png", h_fill=NA, h_color="#7B018C")

## Warning: Using ragg device as default. Ignoring `type` and `antialias`
## arguments

#в рамках неприличия
sticker("./clck.png", package = '',
        p_size=20, s_x=1, s_y=1, s_width=0.65, s_height=0.65,
        filename="img/sticker_hex_lol.png", h_fill="white", h_color="#7B018C")

## Warning: Using ragg device as default. Ignoring `type` and `antialias`
## arguments

Задание №2

Воспользуйтесь данными mtcars.

Каждая из 32 машин охарактеризована 11ю показателями.

Скоррелируйте показатели машин друг с другом. Выделите высокоскоррелированные показатели: коэффициент корреляции Спирмена > 0.8 без учета значимости. Оставьте только такие показатели, которые не являются высокоскоррелированными между собой согласно вышеописанным критериям (из группы должен остаться любой один).
Визуализируйте кореляционные матрицы до и после фильтрации так, чтобы были видны в явном виде значения коэффициентов корреляции.

(3*) Выделите на рисунке, построенном по исходной матрице корреляций, значения, которые подлежат фильтрации.

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

correlation_matrix <- cor(mtcars, method = "spearman")
corrplot(correlation_matrix, method = 'number', number.cex = 5, tl.cex = 5, cl.cex = 5,
         col=colorRampPalette(c("#7B018C","white","#BF9315"))(200))

# мои жалкие потуги отобрать что-либо в R, который вообще не умеет обрабатывать таблицы
# correlation_matrix %>% 
#   as.data.frame %>% 
#   rownames_to_column(var = "first") %>%
#   pivot_longer(cols = where(is.numeric), 
#                names_to = "second", values_to = "corr") %>% 
#   filter(abs(corr) > 0.8) %>% 
#   pivot_wider(names_from = second, values_from = corr)
# рандомная функция, которую подсказал ChatGPT
highly_correlated <- findCorrelation(correlation_matrix, cutoff = 0.8)
correlation_matrix[-highly_correlated, -highly_correlated] %>% 
  corrplot(method = 'number', number.cex = 5, tl.cex = 5, cl.cex = 5,
           col=colorRampPalette(c("#7B018C","white","#BF9315"))(200))

fuck_cor <- mtcars %>% 
  select(-highly_correlated) %>%
  scale() %>%
  t() %>%
  psych::corr.test(adjust = 'BH') #%>% пипец конвейер тут не работает

## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
##   # Was:
##   data %>% select(highly_correlated)
## 
##   # Now:
##   data %>% select(all_of(highly_correlated))
## 
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

corrplot(fuck_cor$r,
         p.mat = fuck_cor$p,
         method = "color", insig = "label_sig", pch = "F", pch.cex = 5, number.cex = 5, tl.cex = 5, cl.cex = 5,
         col=colorRampPalette(c("#7B018C","white","#BF9315"))(200))

Задание №3

Для каждого численного параметра из данных про ирисы (iris) (должно получиться 4 параметра) нарисуйте график, состоящий из двух графиков, расположенных рядом: qqplot и гистограмма (на гистограме подписаны текстом значения асимметрии и эксцесса).

В итоге должно получиться 4 графика, каждый из которых составлен из двух графиков.

Подпишите все оси, графики, название параметров и пр.

hist_plot <- function(df, var, name){
df %>% 
    ggplot(mapping = aes(x = var)) +
    geom_histogram(fill = '#598392', color = '#01161e') + 
  annotate("text", x = mean(var), y=20, size=20,
           label = paste("skew: ", round(skew(var), 2))) +
  annotate("text", x = mean(var), y=17, size=20,
           label = paste("kurtosi", round(kurtosi(var), 2))) + #R настолько ужасен, что он даже не позволяет вставить пeренос строки
    labs(x = name,
       y = "Count")
}
qq_plot <- function(df, var, name){
df %>% 
    ggplot(aes(sample = var)) +
    stat_qq() + 
    stat_qq_line()
}

plots <- c(pmap(list(rep(list(iris), ncol(select(iris, where(is.numeric)))),
     map(select(iris, where(is.numeric)), ~ .x), 
     map(select(iris, where(is.numeric)), ~ .x) %>% names() %>% map(~ .x)),
     hist_plot),
     
  pmap(list(rep(list(iris), ncol(select(iris, where(is.numeric)))),
     map(select(iris, where(is.numeric)), ~ .x), 
     map(select(iris, where(is.numeric)), ~ .x) %>% names() %>% map(~ .x)),
     qq_plot))#%>% и тут конвейер не работает( Почему функции в R не дружат друг с другом?
ggarrange(plotlist = plots, ncol = 4, nrow = 2)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Все параметрические тесты предполагают нормальность распределения!!!

Задание №4

Воспользуйтесь данными mice из пакета datarium.

Узнайте, отличается ли средний вес мышей от 25г?

Не забудьте сформулировать нулевую и альтернативную гипотезы, задать уровень значимости, вывести и p-value и сделать вывод.

if (shapiro.test(mice$weight)$p.value > 0.05) {
 t.test(mice$weight, mu = 25) 
} else {
  print("Не нормальное распределение")
  }

## 
##  One Sample t-test
## 
## data:  mice$weight
## t = -8.1045, df = 9, p-value = 1.995e-05
## alternative hypothesis: true mean is not equal to 25
## 95 percent confidence interval:
##  18.78346 21.49654
## sample estimates:
## mean of x 
##     20.14

Нулевая гипотеза: средний вес мышей равен 25г

Альтернатива: средний вес мышей не равен 25г (двусторонняя)

Уровень значимости: 5%

Т.к. p-value=1.995e-05 < 0.05, то на заданном уровне значимости мы отвергаем нулевую гипотезу и принимаем альтернативную. Вывод: на уровне значимости 5% средний вес мышей не равен 25г ### Задание №5

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

if (shapiro.test(mice$weight)$p.value > 0.05) {
 t.test(mice$weight, mu = 25, alternative = "greater")
} else {
  print("Не нормальное распределение")
  }

## 
##  One Sample t-test
## 
## data:  mice$weight
## t = -8.1045, df = 9, p-value = 1
## alternative hypothesis: true mean is greater than 25
## 95 percent confidence interval:
##  19.04074      Inf
## sample estimates:
## mean of x 
##     20.14

Нулевая гипотеза: средний вес мышей меньше или равен 25г

Альтернатива: средний вес мышей больше 25г (односторонняя)

Уровень значимости: 5%

Т.к. p-value=1 > 0.05, то на заданном уровне значимости мы принимаем нулевую гипотезу и отвергаем альтернативную. Вывод (с учётом задания 4): на уровне значимости 5% средний вес мышей меньше 25г

Задание №6

Воспользуйтесь набором данных genderweight из пакета datarium.

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

Объясните выбор теста, сформулируйте необходимые нулевую и альтернативную гипотезы.

Визуализируйте результаты наиболее наглядным образом.

if (shapiro.test(filter(genderweight, group == "F")$weight)$p.value > 0.05 &
    shapiro.test(filter(genderweight, group == "M")$weight)$p.value > 0.05) {
  print(t.test(weight ~ group, genderweight))
  ggerrorplot(genderweight, x = "group", y = "weight", color = "group",
            desc_stat = "mean_sd",
            add = "jitter", add.params = list(color = "group"), size = 0.8
            )
} else {
  print("Не нормальное распределение")
  }

## 
##  Welch Two Sample t-test
## 
## data:  weight by group
## t = -20.791, df = 26.872, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -24.53135 -20.12353
## sample estimates:
## mean in group F mean in group M 
##        63.49867        85.82612

Используем двувыборочный t-тест, т.к. хотим сравнить средние в двух независимых (неспаренных) выборках

Нулевая гипотеза: средний вес мужчин и женщин не различается

Альтернатива: средний вес мужчин и женщин различается (двусторонняя)

Уровень значимости: 5%

Т.к. p-value<2.2e-16 < 0.05, то на заданном уровне значимости мы отвергаем нулевую гипотезу и принимаем альтернативную. Вывод: на уровне значимости 5% средний вес мужчин и женщин различается

Задание №7

Воспользуйтесь набором данных mice2 из пакета datarium.

Изменился ли в среднем вес мышей после воздействия?

Если изменился, то уменьшился или увеличился?

if (shapiro.test(mice2$before)$p.value > 0.05 &
    shapiro.test(mice2$after)$p.value > 0.05) {
  df <- mice2 %>% 
  pivot_longer(cols = c("before", after), names_to = "time", values_to = "weight")
  t.test(weight ~ time, df, alternative = "greater", paired = T)
} else {
  print("Не нормальное распределение")
  }

## 
##  Paired t-test
## 
## data:  weight by time
## t = 25.546, df = 9, p-value = 5.195e-10
## alternative hypothesis: true mean difference is greater than 0
## 95 percent confidence interval:
##  185.166     Inf
## sample estimates:
## mean difference 
##          199.48

Используем двувыборочный парный t-тест (по факту это одновыборочный t-тест для разностей), т.к. хотим сравнить средние в двух спаренных выборках

Нулевая гипотеза: средний вес до и после воздействия не различается

Альтернатива: средний вес после воздействия увеличился (односторонняя)

Уровень значимости: 5%

Т.к. p-value=5,195e-10 < 0.05, то на заданном уровне значимости мы отвергаем нулевую гипотезу и принимаем альтернативную. Вывод: средний вес мышей после воздействия увеличился

Задание №8

Воспользуйтесь набором данных о звездных войнах.

Разделите персонажей на высоких (рост (height) > 180) и невысоких (рост <= 180).

Связаны ли рост и гендер (gender) персонажей?

test_result <- starwars %>% 
  transmute(name=name, gender=gender, height=if_else(height>180, "tall", "short")) %>% 
  drop_na() %>% 
  count(gender, height) %>% 
  pivot_wider(names_from = height, values_from = n) %>% 
  column_to_rownames("gender") %>% 
  chisq.test()
test_result$p.value

## [1] 0.0246992

Используем χ2 - тест на сопряженность

Нулевая гипотеза: рост и гендер персонажей не связаны

Альтернатива: рост и гендер персонажей связаны

Уровень значимости: 5%

Т.к. p-value=0.025 < 0.05, то на заданном уровне значимости мы отвергаем нулевую гипотезу и принимаем альтернативную. Вывод: между ростом и гендором персонажа имеется связь

Задание №9

Воспользуйтесь данными из таблицы pois.tsv (вообще-то там запятые, а не табы)

Различается ли в среднем время (time) между группами poison? Воспользуйтесь ANOVA.

Если различается, то между какими конкретно группами?

Вспомним матстат.

Для честного проведения ANOVA необходимо проверирить соблюдаются ли следующие предположения:

нормальность распределения каждой из сравниваемых выборок
равенство дисперсий сравниваемых выборок
независимость выборок (нельзя формально проверить, принимаем на веру)

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

pois <- read_csv("./pois.tsv")

## Rows: 48 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): treat
## dbl (3): N, time, poison
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

c(shapiro.test(filter(pois, poison == 1)$time)$p.value > 0.05,
  shapiro.test(filter(pois, poison == 2)$time)$p.value > 0.05,
  shapiro.test(filter(pois, poison == 3)$time)$p.value > 0.05,
  var.test(filter(pois, poison == 1)$time, filter(pois, poison == 2)$time)$p.value > 0.05,
  var.test(filter(pois, poison == 1)$time, filter(pois, poison == 3)$time)$p.value > 0.05,
  var.test(filter(pois, poison == 2)$time, filter(pois, poison == 3)$time)$p.value > 0.05)

## [1]  TRUE FALSE  TRUE  TRUE FALSE FALSE

Понимаем, что к нашим данным нельзя применять ANOVA. Плачем. Ещё раз плачем. Забываем матстат (для этого предварительно отключаем мозг). Проводим ANOVA (но при этом помним, что полученный результат туфта с точки зрения математической статистики).

res.aov <- aov(time ~ poison, data = pois)
summary(res.aov)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## poison       1 0.9316  0.9316   20.67 3.96e-05 ***
## Residuals   46 2.0735  0.0451                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

На уровне значимости 5% есть различия времени хотя бы в одной паре групп poison. Найдём с помощью попарного t-теста, какие именно группы различаются.

pairwise.t.test(pois$time, pois$poison)

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  pois$time and pois$poison 
## 
##   1      2     
## 2 0.3284 -     
## 3 1e-04  0.0015
## 
## P value adjustment method: holm

На уровне значимости 5% есть различия времени между группами в парах 1 – 3 и 2 – 3.

Задание №10

Воспользуйтесь данными из таблицы pois.tsv

Различаются ли в среднем времена отклика между каждой группой по столбцу treat? Проведите необходимое количество сравнений.

Для решения воспользуйтесь непараметрическим тестом. Укажите допущения, если необходимо.

Создайте таблицу, где в строках будут указаны пары сравниваемых групп, для каждого сравнения укажите p-value.

С помощью функции p.adj дополните таблицу поправленными p-value, используйте два любых разных способа поправки, принадлежащих к разным классам (FDR, FWER).

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

Сделайте вывод о том, какой из методов “строже” и почему.

FWER <- pairwise.wilcox.test(pois$time, pois$treat, "holm")$p.value

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): не могу
## подсчитать точное p-значение при наличии повторяющихся наблюдений

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): не могу
## подсчитать точное p-значение при наличии повторяющихся наблюдений

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): не могу
## подсчитать точное p-значение при наличии повторяющихся наблюдений

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): не могу
## подсчитать точное p-значение при наличии повторяющихся наблюдений

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): не могу
## подсчитать точное p-значение при наличии повторяющихся наблюдений

FDR <- pairwise.wilcox.test(pois$time, pois$treat, "fdr")$p.value

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): не могу
## подсчитать точное p-значение при наличии повторяющихся наблюдений

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): не могу
## подсчитать точное p-значение при наличии повторяющихся наблюдений

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): не могу
## подсчитать точное p-значение при наличии повторяющихся наблюдений

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): не могу
## подсчитать точное p-значение при наличии повторяющихся наблюдений

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): не могу
## подсчитать точное p-значение при наличии повторяющихся наблюдений

df <- data.frame(FWER = as.vector(FWER),
                 FDR = as.vector(FDR)) %>% drop_na()
df %>% 
  ggplot(aes(x = FWER, y = FDR)) +
  geom_point() +
  geom_abline(slope=1, intecept=0)

## Warning in geom_abline(slope = 1, intecept = 0): Ignoring unknown parameters:
## `intecept`

FWER строже, т.к. сильнее уменьшает p-value и отвергает больше нулевых гипотез. ***

Помните, что контрольная будет ОЧЕНЬ похожа на это ДЗ?

ТО САМОЕ ДЗ

Муравьёв Георгий

2023-12-08

Задание №1

Задание №2

Задание №3

Все параметрические тесты предполагают нормальность распределения!!!

Задание №4

Задание №6

Задание №7

Задание №8

Задание №9

Задание №10