# install.packages("pacman")
pacman::p_load(readxl, ggplot2, dplyr, nortest, car, emmeans, multcomp, multcompView, boot)Пример решения зачета
Загружаем данные и указываем здесь названия переменной и группы, формулу
data <- read_excel("hw4.shelepin.data.xlsx")
var <- "loss"
group <- "species"
formula <- "loss ~ species"Задание 1
- Проведите статистический анализ [числовая переменная] у разных [групп]. Рассчитайте средние и стандартные отклонения [числовой переменной] для всех [групп]. Есть ли отличия между [группами]? Между какими парами есть статистически значимые отличия и в какую сторону? Обоснуйте выбор метода анализа, при использовании статистических критериев в комментарии приведите p-значения и сформулируйте вывод.
analyze_data <- function(data, var, group) {
# Подготовка формул
formula <- as.formula(formula)
# Линейная модель
fit <- lm(formula, data)
res <- resid(fit)
# Проверка на нормальность распределения остатков
shapiro_result <- shapiro.test(res)
print(shapiro_result)
if (shapiro_result$p.value < 0.05) {
cat("Данные распределены ненормально. Используем непараметрические тесты.\n")
# Краскала-Уоллиса тест
kruskal_result <- kruskal.test(formula, data)
print(kruskal_result)
# Парные сравнения Вилкоксона
wilcox_result <- pairwise.wilcox.test(data[[var]], data[[group]])
print(wilcox_result)
} else {
cat("Данные распределены нормально. Проверяем гомогенность дисперсий.\n")
# Левенов тест
library(car)
levene_result <- leveneTest(formula, data)
print(levene_result)
if (levene_result$`Pr(>F)`[1] < 0.05) {
cat("Дисперсии негомогенны. Используем тест Welch и парные t-тесты без предположения о равенстве дисперсий.\n")
# Тест Вэлча
welch_result <- oneway.test(formula, data)
print(welch_result)
# Парные t-тесты без предположения о равенстве дисперсий
pairwise_t_result <- pairwise.t.test(data[[var]], data[[group]], pool.sd = FALSE)
print(pairwise_t_result)
} else {
cat("Дисперсии гомогенны. Проводим ANOVA и парные t-тесты.\n")
# ANOVA
anova_result <- anova(fit)
print(anova_result)
# Парные t-тесты
pairwise_t_homog_result <- pairwise.t.test(data[[var]], data[[group]])
print(pairwise_t_homog_result)
}
}
}
# Пример использования:
# analyze_data(df, "var", "group")analyze_data(data, var, group)##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.98053, p-value = 0.7091
##
## Данные распределены нормально. Проверяем гомогенность дисперсий.## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 1.0666 0.3754
## 36
## Дисперсии гомогенны. Проводим ANOVA и парные t-тесты.
## Analysis of Variance Table
##
## Response: loss
## Df Sum Sq Mean Sq F value Pr(>F)
## species 3 14.517 4.8391 1.8435 0.1567
## Residuals 36 94.496 2.6249
##
## Pairwise comparisons using t tests with pooled SD
##
## data: data[[var]] and data[[group]]
##
## Allolobophora caliginosa Allolobophora longa
## Allolobophora longa 1.00 -
## Eisenia balatonica 1.00 1.00
## Lumbricus rubellus 0.20 0.45
## Eisenia balatonica
## Allolobophora longa -
## Eisenia balatonica -
## Lumbricus rubellus 0.48
##
## P value adjustment method: holmЗадание 2
- Рассчитайте 99.9 %-ный доверительный интервал медианы [числовой переменной] для каждой [группы] с использованием бутстрепа. Используйте метод процентилей.
library(boot)
# Функция для бутстрепа медианы
boot_median <- function(data, indices) {
median(data[indices])
}
# Функция для вычисления доверительного интервала
calculate_bootstrap_ci <- function(data, var, group) {
unique_groups <- unique(data[[group]])
print(unique_groups)
ci_results <- list()
for (g in unique_groups) {
group_data <- data %>%
filter(group == group) %>%
pull(var)
# Бутстреп
boot_result <- boot(group_data, boot_median, R = 10000)
# Доверительный интервал
ci <- boot.ci(boot_result, type = "perc", conf = 0.999)
ci_results[[g]] <- ci$percent[4:5] # Нижний и верхний пределы
}
return(ci_results)
}
# calculate_bootstrap_ci <- function(data, var, group, R = 1000) {
# res <- boot(data[[var]], boot_median, R)
# ci <- boot.ci(res, conf=0.99, type="perc")
# ci_result <- ci$percent[4:5] # Нижний и верхний пределы
# }
# Пример использования:
# ci_results <- calculate_bootstrap_ci(df, "var", "group")
# print(ci_results)ci_results <- calculate_bootstrap_ci(data, var, group)## [1] "Allolobophora longa" "Eisenia balatonica"
## [3] "Lumbricus rubellus" "Allolobophora caliginosa"print(ci_results)## $`Allolobophora longa`
## [1] 7.942497 10.824933
##
## $`Eisenia balatonica`
## [1] 7.942498 10.824851
##
## $`Lumbricus rubellus`
## [1] 7.942497 10.854745
##
## $`Allolobophora caliginosa`
## [1] 7.940506 10.854728Задание 3
- Постройте итоговый график: столбчатая диаграмма, отображающая медианы [числовой переменной] у разных [групп] с доверительными интервалами. [Группы] должны следовать в порядке убывания [числовой переменной]. Оси должны быть подписаны, столбики для разных [групп] должны отличаться цветом. Разместите на графике компактные буквенные обозначения, отражающие значимость отличий между [группами].
Используя ggplot
library(ggplot2)
library(dplyr)
plot_median_ci <- function(data, var, group, ci_results) {
# Вычисление медиан
medians <- data %>%
group_by(!!sym(group)) %>%
summarize(median_value = median(!!sym(var))) %>%
arrange(desc(median_value))
# Подготовка данных для графика
plot_data <- medians %>%
mutate(lower_ci = sapply(medians[[group]], function(g) ci_results[[g]][1]),
upper_ci = sapply(medians[[group]], function(g) ci_results[[g]][2]))
# Проведение множественных сравнений с использованием emmeans
model <- lm(as.formula(formula), data = data)
model_means <- emmeans(model, specs = group)
model_means_cld <- cld(object = model_means, Letters = letters, alpha = 0.05)
# Получение компактных буквенных обозначений
significance_labels <- model_means_cld$.group
# Создание цветов для групп
group_levels <- levels(factor(plot_data[[group]]))
group_colors <- setNames(rainbow(length(group_levels)), group_levels)
# Построение графика
p <- ggplot(plot_data, aes(x = reorder(!!sym(group), -median_value), y = median_value, fill = !!sym(group))) +
geom_bar(stat = "identity", show.legend = FALSE) +
geom_errorbar(aes(ymin = lower_ci, ymax = upper_ci), width = 0.2) +
labs(x = "Группа", y = paste("Медиана", var), title = "Медианы и доверительные интервалы") +
theme_minimal() +
coord_flip() + # Переворачиваем оси для лучшей читаемости
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
geom_text(aes(label = significance_labels, y = median_value + (upper_ci - median_value) * 0.1),
position = position_dodge(width = 0.9), vjust = -0.5, color = "red", size = 3) +
scale_fill_manual(values = group_colors)
print(p)
}
# Пример использования:
# plot_median_ci(df, "var", "group", ci_results)plot_median_ci(data, var, group, ci_results)
### Используя barplot
plot_median_ci <- function(data, var, group, ci_results) {
# Вычисление медиан
medians <- data %>%
group_by(!!sym(group)) %>%
summarize(median_value = median(!!sym(var))) %>%
arrange(desc(median_value))
# Подготовка данных для графика
plot_data <- medians %>%
mutate(lower_ci = sapply(medians[[group]], function(g) ci_results[[g]][1]),
upper_ci = sapply(medians[[group]], function(g) ci_results[[g]][2]))
# Проведение множественных сравнений с использованием emmeans
model <- lm(as.formula(formula), data = data)
model_means <- emmeans(model, specs = group)
model_means_cld <- cld(object = model_means, Letters = letters, alpha = 0.05)
# Получение компактных буквенных обозначений
significance_labels <- model_means_cld$.group
# Создание цветов для групп
group_levels <- levels(factor(plot_data[[group]]))
group_colors <- setNames(rainbow(length(group_levels)), group_levels)
par(mar = c(5, 8, 4, 2)) # Увеличиваем левый отступ
# Построение горизонтальной столбчатой диаграммы с помощью barplot
barplot_heights <- barplot(plot_data$median_value, names.arg = plot_data[[group]], col = group_colors[plot_data[[group]]], las = 1,
xlab = paste("Медиана", var), ylab = "Группа",
main = "Медианы и доверительные интервалы", horiz = TRUE, xlim=c(0, 12), cex.names = 0.7)
# Добавление доверительных интервалов
arrows(x0 = plot_data$lower_ci, y0 = barplot_heights, x1 = plot_data$upper_ci, y1 = barplot_heights, angle = 90, code = 3, length = 0.1)
# Добавление компактных буквенных обозначений значимости отличий
text(x = plot_data$median_value + (plot_data$upper_ci - plot_data$median_value) * 0.1, y = barplot_heights, labels = significance_labels, col = "red", cex = 0.8, pos = 4)
}plot_median_ci(data, var, group, ci_results)