## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse  = TRUE,
  comment   = "#>",
  fig.width = 7,
  fig.height = 4.5,
  warning   = FALSE,
  message   = FALSE
)

## ----simulate-----------------------------------------------------------------
library(achieveGap)

# Generate data: 400 students in 30 schools, grades K-7
# Gap shape: monotone widening (true gap increases steadily across grades)
sim <- simulate_gap(
  n_students = 400,
  n_schools  = 30,
  gap_shape  = "monotone",
  seed       = 2024
)

# Preview the data
head(sim$data)

## ----true-gap-table-----------------------------------------------------------
# The true gap at each grade
sim$true_gap

## ----fit, cache = TRUE--------------------------------------------------------
# Formula interface (recommended) — same as lme4/nlme style
fit <- achieve_gap(
  score ~ grade,
  group  = "SES_group",
  random = ~ 1 | school/student,
  data   = sim$data,
  k      = 6,     # spline basis dimension (must be < unique grade values)
  n_sim  = 5000   # posterior draws for simultaneous bands
)

## ----print--------------------------------------------------------------------
print(fit)

## ----summary------------------------------------------------------------------
summary(fit)

## ----plot-with-truth, fig.cap = "Estimated gap with both confidence bands and true gap overlaid."----
# Grade labels for x-axis
grade_labs <- c("K", "G1", "G2", "G3", "G4", "G5", "G6", "G7")

# True gap evaluated at the model's grade grid
true_gap_vec <- sim$f1_fun(fit$grade_grid)

plot(
  fit,
  true_gap     = true_gap_vec,
  grade_labels = grade_labs,
  title        = "SES Achievement Gap Trajectory (Simulated Data)"
)

## ----plot-sim-only, fig.cap = "Gap trajectory with simultaneous band only."----
plot(fit, band = "simultaneous", grade_labels = grade_labs)

## ----test---------------------------------------------------------------------
tryCatch(
  test_gap(fit, type = "both"),
  error = function(e) message("test_gap: ", conditionMessage(e))
)

## ----separate, cache = TRUE---------------------------------------------------
sep <- fit_separate(
  data    = sim$data,
  score   = "score",
  grade   = "grade",
  group   = "SES_group",
  school  = "school",
  student = "student"
)

# Compare gap estimates side by side
cat("Joint model gap at grade 4: ", round(fit$gap_hat[fit$grade_grid >= 3.9][1], 3), "\n")
cat("Separate model gap at grade 4:", round(sep$gap_hat[sep$grade_grid >= 3.9][1], 3), "\n")

# Separate model CIs are narrower (anti-conservative)
cat("\nMean CI width - Joint (pointwise):", round(mean(fit$pw_upper - fit$pw_lower), 3))
cat("\nMean CI width - Separate:         ", round(mean(sep$ci_upper - sep$ci_lower), 3), "\n")

## ----nonmono, cache = TRUE, fig.cap = "Non-monotone gap trajectory."----------
sim2 <- simulate_gap(n_students = 400, n_schools = 30,
                     gap_shape = "nonmonotone", seed = 99)

fit2 <- gap_trajectory(
  data    = sim2$data,
  score   = "score",
  grade   = "grade",
  group   = "SES_group",
  school  = "school",
  student = "student",
  n_sim   = 1000,
  verbose = FALSE
)

plot(fit2,
     true_gap     = sim2$f1_fun(fit2$grade_grid),
     grade_labels = grade_labs,
     title        = "Non-Monotone SES Gap: Widens Early, Plateaus, Narrows")

## ----sim-study, eval = FALSE--------------------------------------------------
# results <- run_simulation(
#   n_reps = 5,    # increase to 500 for paper results
#   n_sim  = 1000  # increase to 5000 for final run
# )

## ----sim-summary, eval = FALSE------------------------------------------------
# summarize_simulation(results)

## ----session------------------------------------------------------------------
sessionInfo()