1. Выбор датасета и разведочный анализ данных (EDA)
# Загрузка данных
data(mtcars)
# Структура данных
str(mtcars)
## 'data.frame': 32 obs. of 11 variables:
## $ mpg : num 21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
## $ cyl : num 6 6 4 6 8 6 8 4 4 6 ...
## $ disp: num 160 160 108 258 360 ...
## $ hp : num 110 110 93 110 175 105 245 62 95 123 ...
## $ drat: num 3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
## $ wt : num 2.62 2.88 2.32 3.21 3.44 ...
## $ qsec: num 16.5 17 18.6 19.4 17 ...
## $ vs : num 0 0 1 1 0 1 0 1 1 1 ...
## $ am : num 1 1 1 0 0 0 0 0 0 0 ...
## $ gear: num 4 4 4 3 3 3 3 4 4 4 ...
## $ carb: num 4 4 1 1 2 1 4 2 2 4 ...
# Первые несколько строк
head(mtcars)
## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4
## Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4
## Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
## Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
## Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
## Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1
# Основные статистики
summary(mtcars)
## mpg cyl disp hp
## Min. :10.40 Min. :4.000 Min. : 71.1 Min. : 52.0
## 1st Qu.:15.43 1st Qu.:4.000 1st Qu.:120.8 1st Qu.: 96.5
## Median :19.20 Median :6.000 Median :196.3 Median :123.0
## Mean :20.09 Mean :6.188 Mean :230.7 Mean :146.7
## 3rd Qu.:22.80 3rd Qu.:8.000 3rd Qu.:326.0 3rd Qu.:180.0
## Max. :33.90 Max. :8.000 Max. :472.0 Max. :335.0
## drat wt qsec vs
## Min. :2.760 Min. :1.513 Min. :14.50 Min. :0.0000
## 1st Qu.:3.080 1st Qu.:2.581 1st Qu.:16.89 1st Qu.:0.0000
## Median :3.695 Median :3.325 Median :17.71 Median :0.0000
## Mean :3.597 Mean :3.217 Mean :17.85 Mean :0.4375
## 3rd Qu.:3.920 3rd Qu.:3.610 3rd Qu.:18.90 3rd Qu.:1.0000
## Max. :4.930 Max. :5.424 Max. :22.90 Max. :1.0000
## am gear carb
## Min. :0.0000 Min. :3.000 Min. :1.000
## 1st Qu.:0.0000 1st Qu.:3.000 1st Qu.:2.000
## Median :0.0000 Median :4.000 Median :2.000
## Mean :0.4062 Mean :3.688 Mean :2.812
## 3rd Qu.:1.0000 3rd Qu.:4.000 3rd Qu.:4.000
## Max. :1.0000 Max. :5.000 Max. :8.000
# Матрица корреляций
cor_matrix <- cor(mtcars)
print(round(cor_matrix, 3))
## mpg cyl disp hp drat wt qsec vs am gear
## mpg 1.000 -0.852 -0.848 -0.776 0.681 -0.868 0.419 0.664 0.600 0.480
## cyl -0.852 1.000 0.902 0.832 -0.700 0.782 -0.591 -0.811 -0.523 -0.493
## disp -0.848 0.902 1.000 0.791 -0.710 0.888 -0.434 -0.710 -0.591 -0.556
## hp -0.776 0.832 0.791 1.000 -0.449 0.659 -0.708 -0.723 -0.243 -0.126
## drat 0.681 -0.700 -0.710 -0.449 1.000 -0.712 0.091 0.440 0.713 0.700
## wt -0.868 0.782 0.888 0.659 -0.712 1.000 -0.175 -0.555 -0.692 -0.583
## qsec 0.419 -0.591 -0.434 -0.708 0.091 -0.175 1.000 0.745 -0.230 -0.213
## vs 0.664 -0.811 -0.710 -0.723 0.440 -0.555 0.745 1.000 0.168 0.206
## am 0.600 -0.523 -0.591 -0.243 0.713 -0.692 -0.230 0.168 1.000 0.794
## gear 0.480 -0.493 -0.556 -0.126 0.700 -0.583 -0.213 0.206 0.794 1.000
## carb -0.551 0.527 0.395 0.750 -0.091 0.428 -0.656 -0.570 0.058 0.274
## carb
## mpg -0.551
## cyl 0.527
## disp 0.395
## hp 0.750
## drat -0.091
## wt 0.428
## qsec -0.656
## vs -0.570
## am 0.058
## gear 0.274
## carb 1.000
# Визуализация попарных зависимостей
pairs(mtcars[, c("mpg", "disp", "hp", "wt")],
main = "Парные графики для основных переменных mtcars",
col = "steelblue", pch = 19)
2. Адаптация программы бутстрапа на R
set.seed(42)
# Функция для вычисления R-квадрат на бутстрап выборках
rsq_function <- function(formula, data, indices) {
d <- data[indices, ]
fit <- lm(formula, data = d)
return(summary(fit)$r.square)
}
# Запуск бутстрапа с 2000 повторений
reps <- boot(data = mtcars, statistic = rsq_function, R = 2000,
formula = mpg ~ disp)
# Вывод результатов
print(reps)
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = mtcars, statistic = rsq_function, R = 2000, formula = mpg ~
## disp)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.7183433 0.002527745 0.065857
# Доверительный интервал
ci <- boot.ci(reps, type = c("norm", "basic", "perc", "bca"))
print(ci)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = reps, type = c("norm", "basic", "perc", "bca"))
##
## Intervals :
## Level Normal Basic
## 95% ( 0.5867, 0.8449 ) ( 0.6029, 0.8662 )
##
## Level Percentile BCa
## 95% ( 0.5705, 0.8338 ) ( 0.5426, 0.8161 )
## Calculations and Intervals on Original Scale
3. Визуализация результатов бутстрапа
# Подготовка данных для визуализации
boot_stats <- data.frame(r_squared = reps$t)
original_r2 <- round(reps$t0, 4)
ci_bca <- boot.ci(reps, type = "bca")
ci_lower <- round(ci_bca$bca[4], 4)
ci_upper <- round(ci_bca$bca[5], 4)
# Гистограмма распределения R^2
p1 <- ggplot(boot_stats, aes(x = r_squared)) +
geom_histogram(aes(y = after_stat(density)), bins = 30,
fill = "steelblue", color = "black", alpha = 0.7) +
geom_density(color = "darkred", linewidth = 1) +
geom_vline(xintercept = original_r2, color = "darkgreen",
linetype = "dashed", linewidth = 1) +
geom_vline(xintercept = ci_lower, color = "orange",
linetype = "dotted", linewidth = 1) +
geom_vline(xintercept = ci_upper, color = "orange",
linetype = "dotted", linewidth = 1) +
labs(title = "Бутстрап распределение R²",
subtitle = paste0("Модель: mpg ~ disp (исходный R² = ", original_r2, ")"),
x = expression(R^2), y = "Плотность") +
annotate("text", x = ci_lower - 0.03, y = 4,
label = paste0("95% BCa CI: [", ci_lower, ", ", ci_upper, "]"),
color = "orange", size = 3) +
theme_minimal()
# Q-Q plot для проверки нормальности
p2 <- ggplot(boot_stats, aes(sample = r_squared)) +
stat_qq() +
stat_qq_line(color = "red") +
labs(title = "Q-Q график",
x = "Теоретические квантили",
y = "Выборочные квантили") +
theme_minimal()
# График зависимости с доверительной полосой
p3 <- ggplot(mtcars, aes(x = disp, y = mpg)) +
geom_point(size = 2, color = "steelblue") +
geom_smooth(method = "lm", se = TRUE, color = "darkred", fill = "lightpink") +
labs(title = "Зависимость mpg от disp с 95% доверительной полосой",
x = "Объем двигателя (куб. дюймы)",
y = "Расход топлива (мили на галлон)") +
theme_minimal()
# Вывод графиков
grid.arrange(p1, p2, p3, ncol = 2, nrow = 2)
4. Выводы
## ========================================
## РЕЗУЛЬТАТЫ БУТСТРАП АНАЛИЗА
## ========================================
## Исходная модель: mpg ~ disp
## R² на исходных данных: 0.7183
## Стандартная ошибка R²: 0.0659
## Смещение (bias): 0.0025
## 95% доверительный интервал для R²: [ 0.5426 , 0.8161 ]
## ВЫВОДЫ:
## 1. Модель объясняет 71.8 % дисперсии расхода топлива (mpg) с помощью объема двигателя (disp).
## 2. Стандартная ошибка R² составляет 0.0659 , что указывает на умеренную стабильность оценки.
## 3. Доверительный интервал [ 0.5426 , 0.8161 ] не включает ноль, что подтверждает статистическую значимость модели.
## 4. Распределение бутстрап-статистик близко к нормальному, что подтверждается Q-Q графиком.
## 5. Метод бутстрапа позволил оценить неопределенность параметров модели
## без параметрических предположений о распределении ошибок.