## ----setup, include=FALSE-----------------------------------
library(NHANES)       # install.packages("NHANES") if needed

data(NHANES)

set.seed(42)

## NHANES contains 10,000 rows pooled over two survey cycles; a few
## participants appear twice, so keep one record per ID.
nh <- NHANES[!duplicated(NHANES$ID), ]

## Analysis sample: adults with a recorded BMI
adults <- subset(nh, Age >= 18 & !is.na(BMI))

bmi <- adults$BMI            # continuous response, right-skewed

plot(density(bmi))


## ----bootstrap-fn, echo=TRUE, size="scriptsize"-------------
bootstrap_ci <- function(x, n_boot = 2000, seed = 42) {

  set.seed(seed)
  n          <- length(x)
  obs_median <- median(x)
  boot_meds  <- numeric(n_boot)      # pre-allocate storage

  for (b in 1:n_boot) {
    idx          <- sample(1:n,      # resample row indices
                           size    = n,
                           replace = TRUE)
    x_star       <- x[idx]          # bootstrap sample via indexing
    boot_meds[b] <- median(x_star)  # store one bootstrap replicate
  }

  ci <- quantile(boot_meds, probs = c(0.025, 0.975))

  list(
    estimate  = obs_median,
    ci_lower  = ci[[1]],
    ci_upper  = ci[[2]],
    boot_dist = boot_meds           # keep for diagnostic plots
  )
}





## ----boxplot-gender, echo=F, fig.width=7.5, fig.height=2.9, out.width="95%"----
par(mar = c(4, 4, 2, 1))
boxplot(BMI ~ Gender, data = adults, notch = TRUE,
        col = c("lightsteelblue", "lightsalmon"), border = "grey30",
        main = "Adult BMI by Gender (NHANES)",
        xlab = "Gender", ylab = "BMI (kg/m^2)", cex.main = 0.9,
        outpch = 20, outcex = 0.5, outcol = adjustcolor("grey40", 0.5))
abline(h = median(bmi), col = "firebrick", lwd = 2, lty = 2)
legend("topright", legend = "Overall median", bty = "n", cex = 0.8,
       col = "firebrick", lwd = 2, lty = 2)


## ----run-bootstrap, echo=TRUE-------------------------------
result <- bootstrap_ci(bmi, n_boot = 2000, seed = 42)

cat("Observed median :", round(result$estimate, 2), "\n")
cat("95% Bootstrap CI: [",
    round(result$ci_lower, 2), ",",
    round(result$ci_upper, 2), "]\n")


## ----plot-boot, echo=F, fig.width=8, fig.height=3.0, out.width="100%"----
par(mar = c(4, 4, 2, 1))
hist(result$boot_dist, breaks = 10,
     col = "steelblue", border = "white",
     main = "Bootstrap Distribution of the Median BMI",
     xlab = "Bootstrap Median BMI (kg/m^2)", ylab = "Frequency",
     cex.main = 0.9)
abline(v = result$ci_lower, col = "firebrick", lwd = 2, lty = 2)
abline(v = result$ci_upper, col = "firebrick", lwd = 2, lty = 2)
abline(v = result$estimate,  col = "black",    lwd = 2)
legend("topright", bty = "n", cex = 0.8,
       legend = c("Observed median", "95% CI bounds"),
       col = c("black", "firebrick"), lwd = 2, lty = c(1, 2))


## ----qqplot, echo=F, fig.width=6.5, fig.height=3.0, out.width="90%"----
par(mar = c(4, 4, 2, 1))
qqnorm(result$boot_dist,
       main = "Normal Q-Q Plot: Bootstrap Medians",
       pch = 20, col = adjustcolor("steelblue", 0.4),
       cex.main = 0.9)
qqline(result$boot_dist, col = "firebrick", lwd = 2)


## ----sensitivity, echo=F------------------------------------
B_values  <- c(200, 500, 1000, 2000, 5000)
res_table <- data.frame(B     = B_values,
                         lower = NA_real_,
                         upper = NA_real_,
                         width = NA_real_)

for (i in 1:length(B_values)) {
  r <- bootstrap_ci(bmi, n_boot = B_values[i], seed = 42)
  res_table$lower[i] <- round(r$ci_lower, 2)
  res_table$upper[i] <- round(r$ci_upper, 2)
  res_table$width[i] <- round(r$ci_upper - r$ci_lower, 2)
}

knitr::kable(res_table,
             col.names = c("B", "Lower", "Upper", "Width"),
             caption = "95% Bootstrap CI by number of resamples")


## ----session, echo=FALSE, size="tiny"-----------------------
sessionInfo()