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

formula <- Petal.Length ~ Sepal.Length

set.seed(42)
reps <- boot(data = iris, statistic = rsq_function, R = 2000, formula = formula)

cat("Bootstrap statistics:\n")
## Bootstrap statistics:
print(reps)
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = iris, statistic = rsq_function, R = 2000, formula = formula)
## 
## 
## Bootstrap Statistics :
##      original        bias    std. error
## t1* 0.7599546 -0.0009143115  0.02967252
ci_rsq <- boot.ci(reps, type = "bca")
cat("\n95% BCa Confidence Interval for R²:\n")
## 
## 95% BCa Confidence Interval for R²:
print(ci_rsq)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = reps, type = "bca")
## 
## Intervals : 
## Level       BCa          
## 95%   ( 0.6951,  0.8134 )  
## Calculations and Intervals on Original Scale
par(mfrow = c(1,2))

hist(reps$t, breaks = 50, main = "Histogram of Bootstrap R²", xlab = "R²", col = "lightblue")
abline(v = ci_rsq$bca[4:5], col = c("green", "orange"), lty = 3, lwd = 2)
legend("topright", legend = "95% BCa CI", col = c("green", "orange"), lty = 3)

qqnorm(reps$t, main = "Q-Q Plot of Bootstrap R²")
qqline(reps$t, col = "red")

dev.copy(png, filename = "bootstrap_rsq_r.png", width = 800, height = 400)
## png 
##   3
dev.off()
## png 
##   2