Практическая работа: бутстреп

1 Задание 1-2

data(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 ...

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  
##                 
##                 
## 

table(iris$Species)

## 
##     setosa versicolor  virginica 
##         50         50         50

pairs(iris[, 1:4], col = as.numeric(iris$Species), pch = 19,
      main = "Матрица диаграмм рассеяния (iris)")

boxplot(Petal.Length ~ Species, data = iris,
        main = "Длина лепестка по видам", col = c("lightblue", "lightgreen", "lightpink"))

2 Задание 4

library(boot)
set.seed(0)

rsq_function <- function(formula, data, indices) {
  d <- data[indices, ]
  fit <- lm(formula, data = d)
  return(summary(fit)$r.squared)
}

reps <- boot(data = iris, statistic = rsq_function, R = 2000, formula = Petal.Length ~ Petal.Width)
reps

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = iris, statistic = rsq_function, R = 2000, formula = Petal.Length ~ 
##     Petal.Width)
## 
## 
## Bootstrap Statistics :
##      original       bias    std. error
## t1* 0.9271098 0.0003607995 0.009777404

plot(reps)

ci <- boot.ci(reps, type = "bca")
ci

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = reps, type = "bca")
## 
## Intervals : 
## Level       BCa          
## 95%   ( 0.9018,  0.9427 )  
## Calculations and Intervals on Original Scale

plot(iris$Petal.Width, iris$Petal.Length,
     xlab = "Petal.Width", ylab = "Petal.Length", pch = 19, col = "steelblue",
     main = "Регрессия Petal.Length ~ Petal.Width (iris)")
abline(lm(Petal.Length ~ Petal.Width, data = iris), col = "red", lwd = 2)

hist(reps$t, breaks = 50, col = "steelblue", border = "white",
     main = "Бутстреп-распределение R² (2000 повторов)",
     xlab = "R²")
abline(v = ci$bca[4], col = "red", lty = 2, lwd = 2)
abline(v = ci$bca[5], col = "red", lty = 2, lwd = 2)
abline(v = reps$t0, col = "darkgreen", lwd = 2)
legend("topleft", legend = c("Наблюдаемое R²", "Границы 95% BCa CI"),
       col = c("darkgreen", "red"), lwd = 2, lty = c(1, 2))