Зачет

Шелепин Никита 3724М1БЛнб

Описание данных

Набор данных iris содержит данные о длине и ширине чашелистиков (переменные Sepal.Length и Sepal.Width), а также о длине и ширине лепестков (переменные Petal.Length и Petal.Width) в цветках трех видов ирисов: Iris setosa, I. versicolor и I. virginica (переменная Species: принимает значения setosa, versicolor и virginica соответственно).

Загрузите данные по ирисам. Указание: воспользуйтесь функцией data().

data(iris)
data <- iris
var <- "Sepal.Width"
group <- "Species"
formula <- "Sepal.Width ~ Species"

Проведите статистический анализ ширины чашелистика у разных видов ириса. Рассчитайте средние и стандартные отклонения ширины чашелистика для всех видов. Есть ли отличия между видами? Между какими парами видов есть статистически значимые отличия и в какую сторону? Обоснуйте выбор метода анализа, при использовании статистических критериев в комментарии приведите p-значения и сформулируйте вывод.

# Подготовка формул
formula <- as.formula(formula)

# Линейная модель
fit <- lm(formula, data)
res <- resid(fit)

# Проверка на нормальность распределения остатков
shapiro_result <- shapiro.test(res)
print(shapiro_result)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.98948, p-value = 0.323

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)
  }
}

## Данные распределены нормально. Проверяем гомогенность дисперсий.
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   2  0.5902 0.5555
##       147               
## Дисперсии гомогенны. Проводим ANOVA и парные t-тесты.
## Analysis of Variance Table
## 
## Response: Sepal.Width
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## Species     2 11.345  5.6725   49.16 < 2.2e-16 ***
## Residuals 147 16.962  0.1154                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  data[[var]] and data[[group]] 
## 
##            setosa  versicolor
## versicolor < 2e-16 -         
## virginica  9.1e-10 0.0031    
## 
## P value adjustment method: holm

т. к. все p-value в попарных сравнениях < 0.05, все виды статистически отличаются

Рассчитайте 99.9 %-ный доверительный интервал медианы ширины чашелистика для каждой породы ирисов с использованием бутстрепа. Используйте метод процентилей.

# Функция для бутстрепа медианы
boot_median <- function(data, indices) {
  median(data[indices])
}

unique_groups <- unique(data[[group]])
print(unique_groups)

## [1] setosa     versicolor virginica 
## Levels: setosa versicolor virginica

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]  # Нижний и верхний пределы
}

print(ci_results)

## $setosa
## [1] 2.9 3.2
## 
## $versicolor
## [1] 2.900027 3.200000
## 
## $virginica
## [1] 2.9 3.2

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

# Вычисление медиан
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(8, 4, 4, 2))
barplot_heights <- barplot(plot_data$median_value, names.arg = plot_data[[group]], 
                          col = group_colors[plot_data[[group]]], las = 1,
                          ylab = paste("Медиана", var), xlab = "Группа",
                          main = "Медианы и доверительные интервалы", 
                          ylim=c(0, 4), cex.names = 0.7)

# Добавление доверительных интервалов
arrows(x0 = barplot_heights, y0 = plot_data$lower_ci, 
       x1 = barplot_heights, y1 = plot_data$upper_ci, 
       angle = 90, code = 3, length = 0.1)

# Добавление компактных буквенных обозначений значимости отличий
text(x = barplot_heights * 0.95, 
     y = plot_data$median_value + (plot_data$upper_ci - plot_data$median_value) * 1.5, 
     labels = significance_labels, cex = 1)