Анализ6.4
Семенова Кристина
2026-02-15
4. Бутстрап на языке R
Загрузка и разведочный анализ данных
# Загрузка данных iris
data(iris)
cat("Датасет iris загружен.\n")## Датасет iris загружен.cat("Размерность:", dim(iris), "\n")## Размерность: 150 5cat("Структура данных:\n")## Структура данных: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 строк:print(head(iris))## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa# Основные статистики
cat("\nСводная статистика:\n")##
## Сводная статистика:summary(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100
## 1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300
## Median :5.800 Median :3.000 Median :4.350 Median :1.300
## Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199
## 3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800
## Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500
## Species
## setosa :50
## versicolor:50
## virginica :50
##
##
## Визуализация исходных данных
# График зависимости Petal.Length от Sepal.Length
plot(iris$Sepal.Length, iris$Petal.Length,
main = "Зависимость Petal.Length от Sepal.Length",
xlab = "Sepal.Length (см)",
ylab = "Petal.Length (см)",
col = iris$Species, pch = 19, cex = 1.2)
# Добавляем линию регрессии
abline(lm(Petal.Length ~ Sepal.Length, data = iris),
col = "red", lwd = 2)
# Легенда
legend("topleft",
legend = levels(iris$Species),
col = 1:3, pch = 19, cex = 0.8)
# Добавляем текст с уравнением
model <- lm(Petal.Length ~ Sepal.Length, data = iris)
text(5, 6,
paste("y =", round(coef(model)[1], 2), "+",
round(coef(model)[2], 2), "x"),
cex = 0.8)
Функции для бутстрапа
# Функция для бутстрапа R-квадрат
rsq_function <- function(formula, data, indices) {
d <- data[indices, ] # выборка с возвращением
fit <- lm(formula, data = d) # регрессионная модель
return(summary(fit)$r.square) # возвращаем R-квадрат
}
# Функция для бутстрапа коэффициента
coef_function <- function(formula, data, indices) {
d <- data[indices, ]
fit <- lm(formula, data = d)
return(coef(fit)[2]) # коэффициент при предикторе
}Запуск бутстрапа
set.seed(123) # для воспроизводимости
# Бутстрап для R-квадрат (2000 итераций)
cat("Запуск бутстрапа для R-квадрат...\n")## Запуск бутстрапа для R-квадрат...reps_rsq <- boot(data = iris,
statistic = rsq_function,
R = 2000,
formula = Petal.Length ~ Sepal.Length)
# Бутстрап для коэффициента (2000 итераций)
cat("Запуск бутстрапа для коэффициента...\n")## Запуск бутстрапа для коэффициента...reps_coef <- boot(data = iris,
statistic = coef_function,
R = 2000,
formula = Petal.Length ~ Sepal.Length)Результаты бутстрапа
cat("\n=== РЕЗУЛЬТАТЫ БУТСТРАПА ДЛЯ R-КВАДРАТ ===\n")##
## === РЕЗУЛЬТАТЫ БУТСТРАПА ДЛЯ R-КВАДРАТ ===print(reps_rsq)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = iris, statistic = rsq_function, R = 2000, formula = Petal.Length ~
## Sepal.Length)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.7599546 -0.0008324978 0.02971748cat("\n=== РЕЗУЛЬТАТЫ БУТСТРАПА ДЛЯ КОЭФФИЦИЕНТА ===\n")##
## === РЕЗУЛЬТАТЫ БУТСТРАПА ДЛЯ КОЭФФИЦИЕНТА ===print(reps_coef)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = iris, statistic = coef_function, R = 2000, formula = Petal.Length ~
## Sepal.Length)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 1.858433 0.001390632 0.06947054Доверительные интервалы
# Доверительные интервалы для R-квадрат
ci_rsq <- boot.ci(reps_rsq, type = c("norm", "basic", "perc", "bca"))
# Доверительные интервалы для коэффициента
ci_coef <- boot.ci(reps_coef, type = c("norm", "basic", "perc", "bca"))
cat("\n=== ДОВЕРИТЕЛЬНЫЕ ИНТЕРВАЛЫ ДЛЯ R-КВАДРАТ ===\n")##
## === ДОВЕРИТЕЛЬНЫЕ ИНТЕРВАЛЫ ДЛЯ R-КВАДРАТ ===print(ci_rsq)## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = reps_rsq, type = c("norm", "basic", "perc",
## "bca"))
##
## Intervals :
## Level Normal Basic
## 95% ( 0.7025, 0.8190 ) ( 0.7052, 0.8197 )
##
## Level Percentile BCa
## 95% ( 0.7002, 0.8148 ) ( 0.6969, 0.8114 )
## Calculations and Intervals on Original Scalecat("\n=== ДОВЕРИТЕЛЬНЫЕ ИНТЕРВАЛЫ ДЛЯ КОЭФФИЦИЕНТА ===\n")##
## === ДОВЕРИТЕЛЬНЫЕ ИНТЕРВАЛЫ ДЛЯ КОЭФФИЦИЕНТА ===print(ci_coef)## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = reps_coef, type = c("norm", "basic", "perc",
## "bca"))
##
## Intervals :
## Level Normal Basic
## 95% ( 1.721, 1.993 ) ( 1.716, 1.993 )
##
## Level Percentile BCa
## 95% ( 1.724, 2.001 ) ( 1.710, 1.985 )
## Calculations and Intervals on Original ScaleВизуализация результатов бутстрапа
# Установка параметров графиков
par(mfrow = c(2, 3))
# 1. Исходные данные с линией регрессии
plot(iris$Sepal.Length, iris$Petal.Length,
main = "Исходные данные",
xlab = "Sepal.Length", ylab = "Petal.Length",
col = iris$Species, pch = 19)
abline(lm(Petal.Length ~ Sepal.Length, data = iris),
col = "red", lwd = 2)
# 2. Гистограмма распределения R-квадрат
hist(reps_rsq$t, breaks = 30,
main = "Распределение R-квадрат",
xlab = "R-квадрат",
col = "lightblue",
border = "black")
abline(v = reps_rsq$t0, col = "red", lwd = 2, lty = 2)
abline(v = ci_rsq$percent[4:5], col = "blue", lwd = 2, lty = 2)
legend("topright",
legend = c("Точечная оценка", "95% ДИ"),
col = c("red", "blue"), lty = 2, lwd = 2, cex = 0.7)
# 3. QQ-plot для R-квадрат
qqnorm(reps_rsq$t, main = "QQ-plot: R-квадрат")
qqline(reps_rsq$t, col = "red", lwd = 2)
# 4. Гистограмма распределения коэффициента
hist(reps_coef$t, breaks = 30,
main = "Распределение коэффициента",
xlab = "Коэффициент",
col = "lightgreen",
border = "black")
abline(v = reps_coef$t0, col = "red", lwd = 2, lty = 2)
abline(v = ci_coef$percent[4:5], col = "blue", lwd = 2, lty = 2)
# 5. QQ-plot для коэффициента
qqnorm(reps_coef$t, main = "QQ-plot: Коэффициент")
qqline(reps_coef$t, col = "red", lwd = 2)
# 6. Boxplot бутстрап-статистик
boxplot(list("R-квадрат" = reps_rsq$t,
"Коэффициент" = reps_coef$t),
main = "Boxplot бутстрап-статистик",
col = c("lightblue", "lightgreen"),
ylab = "Значение")
Визуализация результатов бутстрапа
# Сброс параметров
par(mfrow = c(1, 1))Дополнительная визуализация с ggplot2
# Создаем data frame для ggplot
df_rsq <- data.frame(
value = reps_rsq$t,
type = "R-квадрат"
)
df_coef <- data.frame(
value = reps_coef$t,
type = "Коэффициент"
)
df_all <- rbind(df_rsq, df_coef)
# График плотностей
ggplot(df_all, aes(x = value, fill = type)) +
geom_density(alpha = 0.5) +
geom_vline(xintercept = reps_rsq$t0,
color = "red", linetype = "dashed", size = 1) +
geom_vline(xintercept = reps_coef$t0,
color = "blue", linetype = "dashed", size = 1) +
labs(title = "Плотности распределения бутстрап-статистик",
x = "Значение",
y = "Плотность") +
theme_minimal() +
scale_fill_manual(values = c("lightblue", "lightgreen")) +
theme(legend.title = element_blank())
Визуализация плотностей распределения
Сравнение доверительных интервалов
cat("\n=== СРАВНЕНИЕ ТИПОВ ДОВЕРИТЕЛЬНЫХ ИНТЕРВАЛОВ ===\n")##
## === СРАВНЕНИЕ ТИПОВ ДОВЕРИТЕЛЬНЫХ ИНТЕРВАЛОВ ===# Таблица для R-квадрат
ci_table_rsq <- data.frame(
Метод = c("Normal", "Basic", "Percentile", "BCa"),
Нижняя = c(ci_rsq$normal[2], ci_rsq$basic[4], ci_rsq$percent[4], ci_rsq$bca[4]),
Верхняя = c(ci_rsq$normal[3], ci_rsq$basic[5], ci_rsq$percent[5], ci_rsq$bca[5]),
Ширина = c(ci_rsq$normal[3] - ci_rsq$normal[2],
ci_rsq$basic[5] - ci_rsq$basic[4],
ci_rsq$percent[5] - ci_rsq$percent[4],
ci_rsq$bca[5] - ci_rsq$bca[4])
)
cat("\nДоверительные интервалы для R-квадрат:\n")##
## Доверительные интервалы для R-квадрат:# Используем kable для красивого вывода
kable(ci_table_rsq, digits = 4,
caption = "Доверительные интервалы для R-квадрат")| Метод | Нижняя | Верхняя | Ширина |
|---|---|---|---|
| Normal | 0.7025 | 0.8190 | 0.1165 |
| Basic | 0.7052 | 0.8197 | 0.1145 |
| Percentile | 0.7002 | 0.8148 | 0.1145 |
| BCa | 0.6969 | 0.8114 | 0.1145 |
# Таблица для коэффициента
ci_table_coef <- data.frame(
Метод = c("Normal", "Basic", "Percentile", "BCa"),
Нижняя = c(ci_coef$normal[2], ci_coef$basic[4], ci_coef$percent[4], ci_coef$bca[4]),
Верхняя = c(ci_coef$normal[3], ci_coef$basic[5], ci_coef$percent[5], ci_coef$bca[5]),
Ширина = c(ci_coef$normal[3] - ci_coef$normal[2],
ci_coef$basic[5] - ci_coef$basic[4],
ci_coef$percent[5] - ci_coef$percent[4],
ci_coef$bca[5] - ci_coef$bca[4])
)
cat("\nДоверительные интервалы для коэффициента:\n")##
## Доверительные интервалы для коэффициента:kable(ci_table_coef, digits = 4,
caption = "Доверительные интервалы для коэффициента")| Метод | Нижняя | Верхняя | Ширина |
|---|---|---|---|
| Normal | 1.7209 | 1.9932 | 0.2723 |
| Basic | 1.7161 | 1.9930 | 0.2769 |
| Percentile | 1.7238 | 2.0007 | 0.2769 |
| BCa | 1.7100 | 1.9847 | 0.2747 |
Анализ результатов
cat("\n=== АНАЛИЗ РЕЗУЛЬТАТОВ ===\n")##
## === АНАЛИЗ РЕЗУЛЬТАТОВ ===# Основные статистики для R-квадрат
cat("\n1. R-квадрат:\n")##
## 1. R-квадрат:cat(" - Точечная оценка:", round(reps_rsq$t0, 4), "\n")## - Точечная оценка: 0.76cat(" - Среднее бутстрап-оценок:", round(mean(reps_rsq$t), 4), "\n")## - Среднее бутстрап-оценок: 0.7591cat(" - Медиана:", round(median(reps_rsq$t), 4), "\n")## - Медиана: 0.76cat(" - Стандартная ошибка:", round(sd(reps_rsq$t), 4), "\n")## - Стандартная ошибка: 0.0297cat(" - Смещение:", round(mean(reps_rsq$t) - reps_rsq$t0, 4), "\n")## - Смещение: -8e-04# Основные статистики для коэффициента
cat("\n2. Коэффициент регрессии:\n")##
## 2. Коэффициент регрессии:cat(" - Точечная оценка:", round(reps_coef$t0, 4), "\n")## - Точечная оценка: 1.8584cat(" - Среднее бутстрап-оценок:", round(mean(reps_coef$t), 4), "\n")## - Среднее бутстрап-оценок: 1.8598cat(" - Медиана:", round(median(reps_coef$t), 4), "\n")## - Медиана: 1.8611cat(" - Стандартная ошибка:", round(sd(reps_coef$t), 4), "\n")## - Стандартная ошибка: 0.0695cat(" - Смещение:", round(mean(reps_coef$t) - reps_coef$t0, 4), "\n")## - Смещение: 0.0014График сходимости
# График сходимости для R-квадрат
plot(reps_rsq, main = "Сходимость бутстрапа для R-квадрат")
График сходимости бутстрапа
# График сходимости для коэффициента
plot(reps_coef, main = "Сходимость бутстрапа для коэффициента")
График сходимости бутстрапа
Выводы
cat("\n=== ВЫВОДЫ ===\n")##
## === ВЫВОДЫ ===cat("1. Модель линейной регрессии Petal.Length ~ Sepal.Length:\n")## 1. Модель линейной регрессии Petal.Length ~ Sepal.Length:cat(" - Уравнение: Petal.Length =", round(coef(model)[1], 2),
"+", round(coef(model)[2], 2), "* Sepal.Length\n")## - Уравнение: Petal.Length = -7.1 + 1.86 * Sepal.Lengthcat(" - R-квадрат =", round(reps_rsq$t0, 4),
"(модель объясняет", round(reps_rsq$t0 * 100, 1),
"% дисперсии)\n\n")## - R-квадрат = 0.76 (модель объясняет 76 % дисперсии)cat("2. Бутстрап-оценки:\n")## 2. Бутстрап-оценки:cat(" - R-квадрат: 95% ДИ [", round(ci_rsq$percent[4], 4),
",", round(ci_rsq$percent[5], 4), "]\n")## - R-квадрат: 95% ДИ [ 0.7002 , 0.8148 ]cat(" - Коэффициент: 95% ДИ [", round(ci_coef$percent[4], 4),
",", round(ci_coef$percent[5], 4), "]\n\n")## - Коэффициент: 95% ДИ [ 1.7238 , 2.0007 ]cat("3. Статистическая значимость:\n")## 3. Статистическая значимость:cat(" - Все доверительные интервалы не включают 0\n")## - Все доверительные интервалы не включают 0cat(" - Связь между Sepal.Length и Petal.Length статистически значима\n\n")## - Связь между Sepal.Length и Petal.Length статистически значимаcat("4. Качество бутстрапа:\n")## 4. Качество бутстрапа:cat(" - Смещение минимально (порядка",
round(abs(mean(reps_rsq$t) - reps_rsq$t0), 4), ")\n")## - Смещение минимально (порядка 8e-04 )cat(" - Распределения близки к нормальным\n")## - Распределения близки к нормальнымcat(" - 2000 итераций достаточно для сходимости\n")## - 2000 итераций достаточно для сходимости