Bootstrap estimation of student loneliness scores and links to screen time and close friends

# Contoh kode R: Bootstrap Confidence Intervals
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(boot)
})
set.seed(42)
# ---- Gaya grafik ----
pal <- c("#1f6feb", "#e8590c", "#2f9e44", "#ae3ec9", "#f59f00")
warna <- list(bg = "#f4f6fa", panel = "#ffffff", teks = "#1f2933", redup = "#5f6b78", grid = "#e3e8ee", bar = "#1d3557")
tata <- 1
theme_set(theme_minimal(base_size = 11) + theme(legend.position = "bottom"))

# ---- Data sintetis: Y menceng kanan, dipengaruhi X1 dan X2 ----
n <- 200
d <- data.frame(daily_screen_time = round(pmax(rnorm(n, 50, 10), 0), 3),
                number_of_close_friends = round(pmax(rnorm(n, 30, 6), 0), 3))
z1 <- as.numeric(scale(d$daily_screen_time))
z2 <- as.numeric(scale(d$number_of_close_friends))
lat <- runif(1, 0.45, 0.8) * z1 + runif(1, 0.15, 0.4) * sample(c(-1, 1), 1) * z2 + rnorm(n, 0, runif(1, 0.6, 0.85))
w <- exp(runif(1, 0.45, 0.7) * lat)                                  # transformasi eksponensial -> menceng kanan
d$loneliness_score <- round(pmax(100 + 15 * (w - mean(w)) / sd(w), 0.05 * abs(100)), 3)
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Statistik yang di-bootstrap (median, korelasi, koefisien regresi) ----
statistik <- function(data, i) {
  s <- data[i, ]
  c(median = median(s$loneliness_score),
    korelasi = cor(s$loneliness_score, s$daily_screen_time),
    koef = unname(coef(lm(loneliness_score ~ daily_screen_time + number_of_close_friends, data = s))[2]))
}
print(summary(d))
cat("Skewness Y:", round(mean((d$loneliness_score - mean(d$loneliness_score))^3) / sd(d$loneliness_score)^3, 3), "\n")

# ---- Bootstrap nonparametrik R = 2000 ----
b <- boot::boot(d, statistik, R = 2000)
print(b)
ci <- lapply(1:3, function(j) boot::boot.ci(b, conf = 0.95, type = c("perc", "bca"), index = j))
print(ci)

# ---- Interval kepercayaan klasik sebagai pembanding ----
ys <- sort(d$loneliness_score)
kb <- qbinom(c(0.025, 0.975), n, 0.5)                                # CI median berbasis statistik tataan (binomial)
ci_med <- c(ys[max(kb[1], 1)], ys[min(kb[2] + 1, n)])
ct <- cor.test(d$loneliness_score, d$daily_screen_time)                                     # CI Fisher z untuk korelasi
fit <- lm(loneliness_score ~ daily_screen_time + number_of_close_friends, data = d)
ci_koef <- confint(fit)[2, ]                                          # CI t (asumsi normal) untuk koefisien
print(summary(fit))
tabel <- do.call(rbind, lapply(1:3, function(j) {
  klasik <- list(ci_med, as.numeric(ct$conf.int), as.numeric(ci_koef))[[j]]
  data.frame(stat = c("Median", "Correlation", "Coefficient")[j], estimasi = unname(b$t0[j]),
             metode = c("Classic", "Percentile", "BCa"),
             bawah = c(klasik[1], ci[[j]]$percent[4], ci[[j]]$bca[4]),
             atas = c(klasik[2], ci[[j]]$percent[5], ci[[j]]$bca[5]))
}))
tabel$lebar <- tabel$atas - tabel$bawah
print(tabel)
se_boot <- apply(b$t, 2, sd)


# ---- Grafik ----
nm_stat <- c("Median", "Correlation", "Coefficient")
db <- data.frame(stat = factor(rep(nm_stat, each = nrow(b$t)), levels = nm_stat), nilai = as.vector(b$t))
tabel$stat <- factor(tabel$stat, levels = nm_stat)
tabel$metode <- factor(tabel$metode, levels = c("BCa", "Percentile", "Classic"))
garis <- rbind(transform(tabel[tabel$metode != "Classic", ], x = bawah),
               transform(tabel[tabel$metode != "Classic", ], x = atas))
p1 <- ggplot(db, aes(nilai)) +
  geom_histogram(bins = 35, fill = pal[1], colour = warna$panel, alpha = 0.75) +
  geom_vline(data = garis, aes(xintercept = x, colour = metode, linetype = metode), linewidth = 0.8) +
  geom_vline(data = tabel[tabel$metode == "BCa", ], aes(xintercept = estimasi), colour = warna$teks, linewidth = 0.9) +
  facet_wrap(~stat, scales = "free") +
  scale_colour_manual(values = c(BCa = pal[2], Percentile = pal[3]), name = "95% CI") +
  scale_linetype_manual(values = c(BCa = "solid", Percentile = "dashed"), name = "95% CI") +
  scale_x_continuous(n.breaks = 4) +
  labs(title = "Bootstrap distributions (R = 2000)", x = "Bootstrap estimate", y = "Count") +
  theme(legend.position = "bottom")
p2 <- ggplot(tabel, aes(estimasi, metode, colour = metode)) +
  geom_errorbarh(aes(xmin = bawah, xmax = atas), height = 0.3, linewidth = 1) +
  geom_point(size = 2.4) +
  facet_wrap(~stat, scales = "free_x", ncol = 1) +
  scale_colour_manual(values = c(BCa = pal[2], Percentile = pal[3], Classic = warna$redup), guide = "none") +
  scale_x_continuous(n.breaks = 5) +
  labs(title = "95% CI: classic vs percentile vs BCa", x = "Estimate", y = NULL)
idx <- replicate(60, sample(n, replace = TRUE), simplify = FALSE)
gb <- do.call(rbind, lapply(seq_along(idx), function(k) {
  cf <- coef(lm(loneliness_score ~ daily_screen_time, data = d[idx[[k]], ]))
  data.frame(rep = k, a = cf[1], b = cf[2])
}))
p3 <- ggplot(d, aes(daily_screen_time, loneliness_score)) +
  geom_point(colour = pal[4], alpha = 0.45, size = 1.3) +
  geom_abline(data = gb, aes(intercept = a, slope = b), colour = pal[2], alpha = 0.15, linewidth = 0.6) +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE, colour = warna$teks, linewidth = 1) +
  labs(title = "Bootstrap regression lines (60 resamples)", x = "daily_screen_time", y = "loneliness_score")
p_hasil <- if (tata == 1) p1 / (p3 | p2) + plot_layout(heights = c(0.9, 1.2)) else p1 / (p2 | p3) + plot_layout(heights = c(0.9, 1.2))

print(p_hasil)