set.seed(667)
# Defensive package loads: a missing package does not hard-fail the render
pkgs <- c("nlme", "dplyr", "ggplot2")
for (p in pkgs) {
  if (requireNamespace(p, quietly = TRUE)) suppressMessages(library(p, character.only = TRUE))
}

# Six Cities data: local copy first, download fallback
fev_url  <- "https://content.sph.harvard.edu/fitzmaur/ala2e/fev1.txt"
fev_file <- "../../data/fev1.txt"
if (!file.exists(fev_file)) download.file(fev_url, fev_file)
fev <- read.table(fev_file, header = TRUE)   # id, ht, age, baseht, baseage, logfev1

# FLW Ch. 9 outcome: log(FEV1 / height^2)
fev$y <- fev$logfev1 - 2 * log(fev$ht)

# Subject-mean and within-centered age (for the decomposition later)
fev <- fev |>
  group_by(id) |>
  mutate(age_bar = mean(age), age_w = age - age_bar) |>
  ungroup()

c(n_subjects = length(unique(fev$id)), n_obs = nrow(fev))

set.seed(667)
samp_ids <- sample(unique(fev$id), 25)
ggplot(fev |> filter(id %in% samp_ids),
       aes(x = age, y = y, group = id)) +
  geom_line(alpha = 0.5, colour = "steelblue") +
  geom_point(alpha = 0.4, size = 0.8, colour = "steelblue") +
  labs(title = "Six Cities: log(FEV1 / height^2) vs age (25 girls)",
       x = "Age (years)", y = "y = log(FEV1) - 2 log(height)") +
  theme_minimal()

# RE = the Chapter 8 random-intercept mixed model (PROC MIXED in FLW)
re <- lme(y ~ age, random = ~ 1 | id, data = fev, method = "REML")

# FE = OLS with a subject indicator (dummy) per subject; within estimator.
# Regressing on within-centered age + subject dummies gives the within slope
# (PROC GLM with CLASS/ABSORB id in FLW Section 9.7).
fe <- lm(y ~ age_w + factor(id), data = fev)

data.frame(
  model       = c("RE (mixed model, lme)", "FE (OLS + subject dummies)"),
  age_slope   = c(fixef(re)["age"], coef(fe)["age_w"])
)

# FLW eq. 9.5: between (age_bar) + within (age_w) effects in one mixed model
dec <- lme(y ~ age_w + age_bar, random = ~ 1 | id, data = fev, method = "REML")

data.frame(
  effect      = c("within  beta^(L)  (age_w)", "between beta^(C)  (age_bar)"),
  estimate    = c(fixef(dec)["age_w"], fixef(dec)["age_bar"])
)

# Equivalent fit: x and x-bar. The age_bar coefficient is the CONTEXTUAL effect
# beta^(C) - beta^(L); testing it = 0 tests H0: beta^(C) = beta^(L).
dec2 <- lme(y ~ age + age_bar, random = ~ 1 | id, data = fev, method = "ML")
round(summary(dec2)$tTable["age_bar", ], 4)

have_plm <- requireNamespace("plm", quietly = TRUE)
if (have_plm) suppressMessages(library(plm))
if (requireNamespace("lmtest", quietly = TRUE)) suppressMessages(library(lmtest))
if (requireNamespace("sandwich", quietly = TRUE)) suppressMessages(library(sandwich))

if (have_plm) {
  fe_p <- plm(y ~ age, data = fev, index = "id", model = "within")
  re_p <- plm(y ~ age, data = fev, index = "id",
              model = "random", random.method = "swar")  # Swamy-Arora
  print(data.frame(
    model     = c("FE (plm within)", "RE (plm Swamy-Arora)"),
    age_slope = c(coef(fe_p)["age"], coef(re_p)["age"])
  ))
} else {
  message("plm not installed; skipping econometrics supplement.")
}

if (have_plm) {
  print(phtest(fe_p, re_p))
}

if (have_plm && requireNamespace("sandwich", quietly = TRUE)) {
  # Robust (clustered) Hausman, for heteroskedasticity / serial correlation
  print(phtest(fe_p, re_p,
        vcov = function(z) vcovHC(z, method = "arellano", type = "HC0", cluster = "group")))

  # Two-way FE: also absorb common time shocks (here, rounded age groups)
  fev$age_grp <- factor(round(fev$age))
  fe_tw <- plm(y ~ age, data = fev, index = c("id", "age_grp"),
               model = "within", effect = "twoways")
  cat("Two-way FE age slope:", round(coef(fe_tw)["age"], 4), "\n")
}

gen_panel <- function(n = 100, T = 6, beta = 1, rho_xa = 0) {
  alpha <- rnorm(n, 0, 1)
  data.frame(id = rep(1:n, each = T)) |>
    mutate(a = rep(alpha, each = T),
           x = rnorm(n * T, 0, 1) + rho_xa * a,   # rho_xa > 0: x correlated with alpha
           y = 2 + beta * x + a + rnorm(n * T, 0, 1))
}
set.seed(667)
fit_both <- function(rho) {
  d <- gen_panel(rho_xa = rho)
  d <- d |> group_by(id) |> mutate(xb = mean(x), xw = x - xb) |> ungroup()
  c(FE_within = coef(lm(y ~ xw + factor(id), data = d))["xw"],
    RE_mixed  = fixef(lme(y ~ x, random = ~ 1 | id, data = d))["x"])
}
rbind(`Scenario A: rho_xa = 0 (congruent)` = round(fit_both(0), 3),
      `Scenario B: rho_xa = 0.5 (RE biased)` = round(fit_both(0.5), 3))