Полный бутстрап анализ данных Iris

1. Подготовка и описание данных

data(iris)
cat("Структура данных iris:\n")
## Структура данных iris:
str(iris)
## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
cat("\nПервые 6 наблюдений:\n")
## 
## Первые 6 наблюдений:
head(iris) %>% kable()
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
5.1 3.5 1.4 0.2 setosa
4.9 3.0 1.4 0.2 setosa
4.7 3.2 1.3 0.2 setosa
4.6 3.1 1.5 0.2 setosa
5.0 3.6 1.4 0.2 setosa
5.4 3.9 1.7 0.4 setosa

Набор данных содержит измерения 150 цветков ириса по 4 параметрам:

Для трех видов: * setosa * versicolor * virginica

2. Разведочный анализ

2.1 Основные статистики

summary(iris) %>% kable(caption = "Сводная статистика")
Сводная статистика
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100 setosa :50
1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300 versicolor:50
Median :5.800 Median :3.000 Median :4.350 Median :1.300 virginica :50
Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199 NA
3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800 NA
Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500 NA

2.2 Визуализация распределений

# График зависимости признаков
p1 <- ggplot(iris, aes(Sepal.Width, Sepal.Length, color = Species)) +
  geom_point(size = 3, alpha = 0.7) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Зависимость длины от ширины чашелистика",
       x = "Ширина чашелистика", y = "Длина чашелистика") +
  theme_minimal(base_size = 12)

# Распределение по видам
p2 <- ggplot(iris, aes(Species, Sepal.Length, fill = Species)) +
  geom_boxplot(alpha = 0.7) +
  geom_jitter(width = 0.1, alpha = 0.5) +
  labs(title = "Распределение длины чашелистика",
       y = "Длина чашелистика") +
  theme_minimal(base_size = 12)

gridExtra::grid.arrange(p1, p2, ncol = 2)

3. Бутстрап анализ R-квадрата

3.1 Определение функции для бутстрапа

rsq_function <- function(formula, data, indices) {
  d <- data[indices, ]  # Бутстрап выборка
  fit <- lm(formula, data = d)  # Линейная регрессия
  return(summary(fit)$r.squared)  # Возвращаем R²
}

3.2 Запуск бутстрапа

set.seed(123)
bootstrap_results <- boot(
  data = iris,
  statistic = rsq_function,
  R = 2000,
  formula = Sepal.Length ~ Sepal.Width
)

cat("Результаты бутстрап анализа:\n")
## Результаты бутстрап анализа:
print(bootstrap_results)
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = iris, statistic = rsq_function, R = 2000, formula = Sepal.Length ~ 
##     Sepal.Width)
## 
## 
## Bootstrap Statistics :
##       original      bias    std. error
## t1* 0.01382265 0.005852842  0.01898089

3.3 Доверительные интервалы

ci_perc <- boot.ci(bootstrap_results, type = "perc")
ci_bca <- boot.ci(bootstrap_results, type = "bca")

intervals <- data.frame(
  Тип = c("Перцентильный", "BCa"),
  Нижняя_граница = c(ci_perc$percent[4], ci_bca$bca[4]),
  Верхняя_граница = c(ci_perc$percent[5], ci_bca$bca[5])
)

kable(intervals, caption = "Доверительные интервалы для R²")
Доверительные интервалы для R²
Тип Нижняя_граница Верхняя_граница
Перцентильный 9.93e-05 0.0692268
BCa 6.16e-05 0.0632579

4. Визуализация результатов

4.1 Распределение бутстрап оценок

ggplot(data.frame(R_squared = bootstrap_results$t), aes(R_squared)) +
  geom_histogram(bins = 30, fill = "lightblue", color = "black") +
  geom_vline(xintercept = bootstrap_results$t0, 
             color = "red", linetype = "dashed", size = 1) +
  labs(title = "Распределение R² из бутстрап выборок",
       x = "R-squared", y = "Частота") +
  theme_minimal(base_size = 14)

4.2 Диагностические графики

par(mfrow = c(2, 2))
plot(bootstrap_results)
title(main = "Диагностические графики бутстрап анализа", outer = TRUE, line = -1)

5. Анализ результатов

5.1 Основные показатели

results <- data.frame(
  Показатель = c("Оригинальное R²", "Среднее R² (бутстрап)", 
                "Смещение", "Стандартная ошибка"),
  Значение = c(
    round(bootstrap_results$t0, 4),
    round(mean(bootstrap_results$t), 4),
    round(mean(bootstrap_results$t) - bootstrap_results$t0, 4),
    round(sd(bootstrap_results$t), 4)
  )
)

kable(results, caption = "Основные результаты анализа")
Основные результаты анализа
Показатель Значение
Оригинальное R² 0.0138
Среднее R² (бутстрап) 0.0197
Смещение 0.0059
Стандартная ошибка 0.0190