- Загрузите необходимые данные из файла hw2.shelepin.data.xlsx
library(readxl)
hw2_data <- read_excel("hw2.shelepin.data.xlsx",
sheet = "activity")
summary(hw2_data)## type rate
## Length:60 Min. : 43.67
## Class :character 1st Qu.: 52.19
## Mode :character Median : 58.24
## Mean : 60.41
## 3rd Qu.: 66.50
## Max. :100.00- Рассчитайте описательные статистики и заполните таблицу
library(dplyr)##
## Присоединяю пакет: 'dplyr'## Следующие объекты скрыты от 'package:stats':
##
## filter, lag## Следующие объекты скрыты от 'package:base':
##
## intersect, setdiff, setequal, unionresult <- hw2_data %>%
group_by(type) %>%
summarise(
Среднее = mean(rate),
Дисперсия = var(rate),
Стандартное_отклонение = sd(rate),
Объем_выборки = n(),
Минимум = min(rate),
Максимум = max(rate)
)
# Вывод результата
print(result)## # A tibble: 4 × 7
## type Среднее Дисперсия Стандартное_отклонение Объем_выборки Минимум Максимум
## <chr> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
## 1 белко… 64.3 105. 10.3 16 46.6 83
## 2 разде… 67.7 206. 14.4 14 46.4 100
## 3 среди… 56.7 97.8 9.89 14 43.7 74
## 4 фрукт… 53.4 65.6 8.10 16 44.0 71- Преобразуйте фактор таким образом, чтобы категории следовали в порядке роста средней активности амилазы. Переименуйте имена категорий для более компактного отображения на графиках (сократите имена, можно использовать аббревиатуры)
# Рассчитаем средние значения для каждой ассоциации
means <- aggregate(rate ~ type, data = hw2_data, FUN = mean)
means <- means[order(means$rate), ] # сортируем по возрастанию средних
# Создаем словарь сокращений
abbreviations <- c(
"белково-бургерная" = "ББ",
"раздельное питание" = "РП",
"средиземноморская" = "СМ",
"фруктово-ягодно-овощная" = "ФЯО"
)
# Создаем новый фактор с правильным порядком и сокращениями
hw2_data$type_n <- factor(hw2_data$type,
levels = means$type,
labels = abbreviations[means$type])
# Проверяем результат
levels(hw2_data$type_n) # покажет новые уровни фактора в правильном порядке## [1] "ФЯО" "СМ" "ББ" "РП"- Новая таблица
result <- hw2_data %>%
group_by(type_n) %>%
summarise(
Среднее = mean(rate),
Дисперсия = var(rate),
Стандартное_отклонение = sd(rate),
Объем_выборки = n(),
Минимум = min(rate),
Максимум = max(rate)
)
# Вывод результата
print(result)## # A tibble: 4 × 7
## type_n Среднее Дисперсия Стандартное_отклонение Объем_выборки Минимум Максимум
## <fct> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
## 1 ФЯО 53.4 65.6 8.10 16 44.0 71
## 2 СМ 56.7 97.8 9.89 14 43.7 74
## 3 ББ 64.3 105. 10.3 16 46.6 83
## 4 РП 67.7 206. 14.4 14 46.4 100- Сформируйте модель дисперсионного анализа исследуемой зависимости
# Создаем модель дисперсионного анализа
fit <- lm(rate ~ type_n, data = hw2_data)
# Выводим результаты
summary(fit)##
## Call:
## lm(formula = rate ~ type_n, data = hw2_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.331 -7.665 -1.259 6.314 32.314
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 53.381 2.697 19.795 < 2e-16 ***
## type_nСМ 3.277 3.948 0.830 0.409955
## type_nББ 10.958 3.814 2.873 0.005730 **
## type_nРП 14.305 3.948 3.624 0.000627 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.79 on 56 degrees of freedom
## Multiple R-squared: 0.2327, Adjusted R-squared: 0.1916
## F-statistic: 5.66 on 3 and 56 DF, p-value: 0.001851# Получаем остатки из модели
residuals <- residuals(fit)
# 1. Критерий Шапиро-Уилка
shapiro_test <- shapiro.test(residuals)
# 2. Критерий Колмогорова-Смирнова
# Сначала стандартизируем остатки
standardized_residuals <- scale(residuals)
ks_test <- ks.test(standardized_residuals, "pnorm")## Warning in ks.test.default(standardized_residuals, "pnorm"): в тесте согласия
## Колмогорова не должно быть повторяющихся значений# 3. Критерий Андерсона-Дарлинга
library(nortest)
ad_test <- ad.test(residuals)# Выводим результаты
cat("Тест Шапиро-Уилка: p =", shapiro_test$p.value, "\n")## Тест Шапиро-Уилка: p = 0.01024024# Выводим результаты
cat("Тест Колмогорова-Смирнова: p =", ks_test$p.value, "\n")## Тест Колмогорова-Смирнова: p = 0.1802803# Выводим результаты
cat("Тест Андерсона-Дарлинга: p =", ad_test$p.value, "\n")## Тест Андерсона-Дарлинга: p = 0.001824153# Гистограмма остатков
# Визуальная проверка
# Q-Q график
qqnorm(residuals)
# Гистограмма остатков
hist(residuals, freq = FALSE)
Тест Шапиро-Уилка: p = 0.01024024 Это значение меньше стандартного уровня значимости 0.05. Следовательно, мы отвергаем нулевую гипотезу о нормальности распределения. Данные, вероятно, не распределены нормально. Тест Колмогорова-Смирнова: p = 0.1802803 Это значение больше 0.05, что означает, что мы не можем отвергнуть нулевую гипотезу. Согласно этому тесту, данные могут быть нормально распределены. Тест Андерсона-Дарлинга: p = 0.001824153 Это значение значительно меньше 0.05, что снова приводит к отвержению нулевой гипотезы о нормальности. Вывод: Два из трех тестов (Шапиро-Уилка и Андерсона-Дарлинга) указывают на то, что данные, вероятно, не распределены нормально. Тест Колмогорова-Смирнова не отвергает гипотезу о нормальности, но он менее мощный для обнаружения отклонений от нормальности по сравнению с другими двумя тестами.
kruskal.test(rate ~ type_n, data = hw2_data)##
## Kruskal-Wallis rank sum test
##
## data: rate by type_n
## Kruskal-Wallis chi-squared = 15.278, df = 3, p-value = 0.001594- P-значение (0.001594) значительно меньше стандартного уровня значимости 0.05. Это означает, что мы отвергаем нулевую гипотезу.
- Нулевая гипотеза в тесте Краскела-Уоллиса предполагает, что все группы имеют одинаковое распределение (или одинаковые медианы).
- Отвержение нулевой гипотезы указывает на то, что существуют статистически значимые различия между группами.
Вывод: Существуют статистически значимые различия в распределении переменной “rate” между четырьмя группами, определенными переменной “type_n”. Иными словами, тип (type_n) оказывает существенное влияние на показатель rate.
pairwise.wilcox.test(hw2_data$rate, hw2_data$type_n)## 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-значение при наличии повторяющихся наблюдений
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): не могу
## подсчитать точное p-значение при наличии повторяющихся наблюдений##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: hw2_data$rate and hw2_data$type_n
##
## ФЯО СМ ББ
## СМ 0.8832 - -
## ББ 0.0098 0.1592 -
## РП 0.0098 0.0809 0.8832
##
## P value adjustment method: holmlibrary(emmeans)## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'library(multcomp)## Загрузка требуемого пакета: mvtnorm## Загрузка требуемого пакета: survival## Загрузка требуемого пакета: TH.data## Загрузка требуемого пакета: MASS##
## Присоединяю пакет: 'MASS'## Следующий объект скрыт от 'package:dplyr':
##
## select##
## Присоединяю пакет: 'TH.data'## Следующий объект скрыт от 'package:MASS':
##
## geyserlibrary(multcompView)
# set up model
model <- lm(rate ~ type_n, data = hw2_data)
# get (adjusted) weight means per group
model_means <- emmeans(object = model,
specs = "type_n")
# add letters to each mean
model_means_cld <- cld(object = model_means,
Letters = letters,
alpha = 0.05)
# show output
model_means_cld## type_n emmean SE df lower.CL upper.CL .group
## ФЯО 53.4 2.70 56 48.0 58.8 a
## СМ 56.7 2.88 56 50.9 62.4 ab
## ББ 64.3 2.70 56 58.9 69.7 bc
## РП 67.7 2.88 56 61.9 73.5 c
##
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.library(tidyverse) # ggplot & helper functions## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ readr 2.1.5
## ✔ ggplot2 3.5.1 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ MASS::select() masks dplyr::select()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(scales) # more helper functions##
## Присоединяю пакет: 'scales'
##
## Следующий объект скрыт от 'package:purrr':
##
## discard
##
## Следующий объект скрыт от 'package:readr':
##
## col_factor# optional: sort factor levels of groups column according to highest mean
# ...in means table
model_means_cld <- model_means_cld %>%
mutate(type_n = fct_reorder(type_n, emmean))
# ...in data table
hw2_data <- hw2_data %>%
mutate(type_n = fct_relevel(type_n, levels(model_means_cld$type_n)))
# base plot setup
ggplot() +
# y-axis
scale_y_continuous(
name = "Rate",
limits = c(0, NA),
breaks = pretty_breaks(),
expand = expansion(mult = c(0,0.1))
) +
# x-axis
scale_x_discrete(
name = "Type"
) +
# general layout
theme_classic() +
# colored boxplots (вместо черного)
geom_boxplot(
data = hw2_data,
aes(y = rate, x = type_n, fill = type_n),
width = 0.5
) +
# добавляем цветовую палитру
scale_fill_brewer(palette = "Set3") +
# black data points
geom_point(
data = hw2_data,
aes(y = rate, x = type_n),
shape = 16,
alpha = 0.5,
position = position_jitter(width = 0.1)
) +
# red letters
geom_text(
data = model_means_cld,
aes(
y = emmean,
x = type_n,
label = str_trim(.group)
),
position = position_nudge(x = 0, y = -24),
hjust = 0.5,
color = "red"
) +
# добавляем стрелку сверху
annotate("segment",
x = 1,
xend = length(unique(hw2_data$type_n)),
y = max(hw2_data$rate) + 2,
yend = max(hw2_data$rate) + 2,
arrow = arrow(length = unit(0.5, "cm"))) +
annotate("text",
x = length(unique(hw2_data$type_n))/2,
y = max(hw2_data$rate) + 6,
label = "Увеличение средней активности") +
# убираем легенду
theme(legend.position = "none") +
# caption
labs(
caption = str_wrap("Black dots represent raw data.Means not sharing any letter are significantly different by the Tukey-test at the 5% level of significance.", width = 70)
)