x <- c(65, 80, 67, 78, 69, 67, 60, 70, 72, 67, 80, 78, 74, 70, 64)
y <- c(70, 73, 70, 80, 80, 83, 78, 76, 78, 74, 78, 80, 65, 85, 80)

model <- lm(y ~ x)
summary(model)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.615  -3.123   1.353   3.377   8.322 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 77.78107   17.23109   4.514 0.000582 ***
## x           -0.01575    0.24275  -0.065 0.949240    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.583 on 13 degrees of freedom
## Multiple R-squared:  0.0003239,  Adjusted R-squared:  -0.07657 
## F-statistic: 0.004212 on 1 and 13 DF,  p-value: 0.9492

y_pred <- fitted(model)
residual <- resid(model)

B <- 1000

intercept_boot <- numeric(B)
slope_boot <- numeric(B)

set.seed(123)

for(i in 1:B){
  residual_boot <- sample(residual,
                          size = length(residual),
                          replace = TRUE)
  y_boot <- y_pred + residual_boot
  model_boot <- lm(y_boot ~ x)
  
  intercept_boot[i] <- coef(model_boot)[1]
  slope_boot[i] <- coef(model_boot)[2]
}

mean(intercept_boot)

## [1] 77.21981

mean(slope_boot)

## [1] -0.00878834

cat("Estimasi bootstrap intercept =", mean(intercept_boot), "\n")

## Estimasi bootstrap intercept = 77.21981

cat("Estimasi bootstrap slope =", mean(slope_boot), "\n")

## Estimasi bootstrap slope = -0.00878834

quantile(intercept_boot, c(0.025, 0.975))

##      2.5%     97.5% 
##  46.95323 108.10543

quantile(slope_boot, c(0.025, 0.975))

##       2.5%      97.5% 
## -0.4491062  0.4241199

hist(slope_boot,
     main = "Histogram Bootstrap Slope",
     xlab = "Nilai Slope")