set.seed(6)
arm_cols <- c(Placebo="gray40", Drug="#1369B0",
              A="gray40", B="#1369B0",
              "Current smoker"="#1369B0", "Former smoker"="gray40",
              Succimer="#1369B0")
ok_gg  <- requireNamespace("ggplot2", quietly=TRUE)
ok_d   <- requireNamespace("dplyr", quietly=TRUE)
ok_tid <- requireNamespace("tidyr", quietly=TRUE)
ok_nlme<- requireNamespace("nlme", quietly=TRUE)
ok_emm <- requireNamespace("emmeans", quietly=TRUE)
if (ok_nlme) library(nlme)

# Simulate 133 subjects, irregular times (0,3,6,9,12,15,19 years), two groups (A vs B).
# Deliberately build in CURVATURE so the quadratic machinery has something to find.
N <- 133
times <- c(0,3,6,9,12,15,19)
id <- rep(seq_len(N), each=length(times))
time <- rep(times, N)
group <- rep(sample(c("A","B"), N, TRUE, prob=c(0.5,0.5)), each=length(times))

# Decline with mild curvature so a quadratic has a genuine turning point.
# Group B declines slightly faster and has slightly stronger curvature.
mu0 <- 3.50
slope_A <- -0.060; slope_B_delta <- -0.010   # linear coefficient (per year)
curv_A  <-  0.0018; curv_B_delta <- 0.0006    # positive quadratic = bend back up
slope <- ifelse(group=="B", slope_A + slope_B_delta, slope_A)
curv  <- ifelse(group=="B", curv_A + curv_B_delta, curv_A)
mu <- mu0 + slope*time + curv*time^2

# AR(1)-ish within-subject noise (course-standard rho = 0.8)
rho <- 0.8; sig <- 0.20; u <- rnorm(N,0,0.25); y <- numeric(length(id))
for (i in seq_len(N)) {
  idx <- which(id==i); y[idx][1] <- mu[idx][1] + u[i] + rnorm(1,0,sig)
  for (k in 2:length(idx)) y[idx][k] <- mu[idx][k] + u[i] + rho*(y[idx][k-1]-mu[idx][k-1]-u[i]) + rnorm(1,0,sig*sqrt(1-rho^2))
}
# Group A is the reference level so coefficients read as "groupB"
fev <- data.frame(id=factor(id), time=time,
                  group=factor(group, levels=c("A","B")), y=y)

if (ok_gg) {
  library(ggplot2)
  library(dplyr)
  summ <- fev |> dplyr::group_by(group, time) |> dplyr::summarise(m=mean(y), .groups="drop")
  ggplot(summ, aes(time, m, color=group)) +
    geom_line(linewidth=1.1) + geom_point() +
    scale_color_manual(values=arm_cols) +
    labs(title="Simulated curved means by group", x="Time (years)", y="Mean") +
    theme(legend.position="bottom") -> p_fev
  print(p_fev)
} else cat("Install ggplot2 for plots.\n")

if (ok_nlme) {
  fev$time_c <- fev$time - mean(fev$time)  # centering
  # First fit, shown as a plain gls() call (safe_gls is the same thing):
  fit_lin_fev <- gls(y ~ group * time_c, data = fev,
                     correlation = corCompSymm(form = ~1|id), method = "REML")
  print(summary(fit_lin_fev)$tTable)
} else cat("Install nlme to fit GLS models.\n")

if (ok_nlme) {
  fit_lin_fev_ml  <- gls(y ~ group * time_c, data = fev,
                         correlation = corCompSymm(form = ~1|id), method = "ML")
  fit_quad_fev    <- gls(y ~ group * (time_c + I(time_c^2)), data = fev,
                         correlation = corCompSymm(form = ~1|id), method = "ML")
  print(anova(fit_lin_fev_ml, fit_quad_fev))
  p <- anova(fit_lin_fev_ml, fit_quad_fev)$`p-value`[2]
  cat("\nVerdict:", if (p < 0.05) "reject linear, keep the quadratic term."
                    else "no evidence of curvature, linear suffices.", "\n")
} else cat("Install nlme to fit GLS models.\n")

if (ok_nlme) {
  # fit_quad_fev was fit (ML) on the previous slide; reuse it
  cf   <- coef(fit_quad_fev)
  tbar <- mean(fev$time)

  # Reference group (A): turning point = -b_t / (2 b_t2) + tbar
  b_t  <- cf["time_c"]
  b_t2 <- cf["I(time_c^2)"]
  tp_ref <- -b_t / (2 * b_t2) + tbar

  # Group B adds the interaction coefficients to each term
  b_t_C  <- b_t  + cf["groupB:time_c"]
  b_t2_C <- b_t2 + cf["groupB:I(time_c^2)"]
  tp_cur <- -b_t_C / (2 * b_t2_C) + tbar

  data.frame(group = c("A","B"),
             t_star = round(c(tp_ref, tp_cur), 2))
} else {
  cat("Install nlme to fit the quadratic model and compute turning points.\n")
}

if (ok_emm) {
  library(emmeans)
  emt <- emtrends(fit_lin_fev, ~ group, var = "time_c", data = fev)
  print(emt)
  print(contrast(emt, "pairwise"))
} else {
  cat("Install emmeans to compute slope contrasts.\n")
}

if (ok_nlme) {
  fev$time_s <- as.numeric(scale(fev$time, center=TRUE, scale=TRUE))
  fit_poly <- safe_gls(y ~ group * poly(time_s, 2, raw = FALSE), data=fev,
                       corr = corCompSymm(form = ~1|id), method="ML")
  cat("AIC (orthogonal degree-2):", AIC(fit_poly), "\n")
} else cat("Install nlme to fit GLS models.\n")

smk_url  <- "https://content.sph.harvard.edu/fitzmaur/ala2e/smoking-data.txt"
smk_file <- "../../data/smoking-data.txt"
if (!file.exists(smk_file)) {
  dir.create("../../data", showWarnings = FALSE)
  download.file(smk_url, smk_file)
}
smk <- read.table(smk_file, header = FALSE)
names(smk) <- c("id", "smoke", "time", "fev1")   # smoke: 0 = current, 1 = former
smk$smoke  <- factor(smk$smoke, labels = c("Current smoker", "Former smoker"))
cat("subjects:", length(unique(smk$id)),
    " observations:", nrow(smk),
    " time range:", paste(range(smk$time), collapse = " to "),
    " FEV1 range:", paste(round(range(smk$fev1), 2), collapse = " to "), "\n")

if (ok_nlme) {
  fit_lin_ml  <- gls(fev1 ~ smoke * time, data = smk,
                     correlation = corSymm(form = ~1|id),
                     weights = varIdent(form = ~1|time),
                     method = "ML", na.action = na.omit)
  fit_quad_ml <- gls(fev1 ~ smoke * (time + I(time^2)), data = smk,
                     correlation = corSymm(form = ~1|id),
                     weights = varIdent(form = ~1|time),
                     method = "ML", na.action = na.omit)
  print(anova(fit_lin_ml, fit_quad_ml))
  p_vv <- anova(fit_lin_ml, fit_quad_ml)$`p-value`[2]
  cat("\nVerdict (quadratic vs linear): p =", round(p_vv, 3), "->",
      if (p_vv < 0.05) "curvature needed." else "linear adequate.", "\n")
} else cat("Install nlme to fit GLS models.\n")

# TLC-like: four planned weeks 0,1,4,6; Placebo vs Succimer; sharp early drop then flatten/rebound
n <- 200; weeks <- c(0,1,4,6)
id2 <- rep(seq_len(n), each=length(weeks))
week <- rep(weeks, n)
arm <- rep(sample(c("Placebo","Succimer"), n, TRUE), each=length(weeks))

# Means: Placebo ~ flat decline; Succimer big drop at week1 then rebound
muP <- 26 + c(0, -1.6, -2.2, -2.6)
muS <- 26 + c(0, -12.8, -9.6, -7.2)
mu <- ifelse(arm=="Placebo", muP[match(week, weeks)], muS[match(week, weeks)])

rho <- 0.8; sig <- 6; u2 <- rnorm(n,0,3); y2 <- numeric(length(id2))
for (i in seq_len(n)) {
  idx <- which(id2==i); y2[idx][1] <- mu[idx][1] + u2[i] + rnorm(1,0,sig)
  for (k in 2:length(idx)) y2[idx][k] <- mu[idx][k] + u2[i] + rho*(y2[idx][k-1]-mu[idx][k-1]-u2[i]) + rnorm(1,0,sig*sqrt(1-rho^2))
}
tlc <- data.frame(id=factor(id2), week=week, arm=factor(arm), y=y2)

if (ok_gg) {
  library(ggplot2); library(dplyr)
  s2 <- tlc |> dplyr::group_by(arm, week) |> dplyr::summarise(m=mean(y), .groups="drop")
  ggplot(s2, aes(week, m, color=arm)) + geom_line(linewidth=1.1) + geom_point() +
    scale_color_manual(values=arm_cols) +
    labs(title="TLC means by arm", x="Week", y="Mean lead level") +
    theme(legend.position="bottom") -> p_tlc
  print(p_tlc)
}

if (ok_nlme) {
  tlc$w1 <- pmax(tlc$week - 1, 0)
  fit_lin  <- gls(y ~ arm * week,        data = tlc, correlation = corCompSymm(form = ~1|id), method = "ML")
  fit_spl  <- gls(y ~ arm * (week + w1), data = tlc, correlation = corCompSymm(form = ~1|id), method = "ML")
  print(anova(fit_lin, fit_spl))
  p2 <- anova(fit_lin, fit_spl)$`p-value`[2]
  cat("\nVerdict:", if (p2 < 0.05) "reject linear, keep the knot at week 1."
                    else "no evidence for a knot, single linear trend suffices.", "\n")
  print(summary(fit_spl)$tTable)
} else cat("Install nlme to fit GLS models.\n")

if (ok_nlme) {
  cf <- coef(fit_spl)
  # Pre-knot slope uses 'week'; post-knot adds 'w1'. Arm B adds its interactions.
  pre_ref  <- cf["week"]
  post_ref <- cf["week"] + cf["w1"]
  pre_oth  <- cf["week"] + cf["armSuccimer:week"]
  post_oth <- cf["week"] + cf["w1"] + cf["armSuccimer:week"] + cf["armSuccimer:w1"]
  data.frame(arm     = c("Placebo","Placebo","Succimer","Succimer"),
             segment = c("pre-knot","post-knot","pre-knot","post-knot"),
             slope   = round(c(pre_ref, post_ref, pre_oth, post_oth), 2))
} else cat("Install nlme to fit the spline and show numeric slopes.\n")

if (ok_gg && ok_nlme) {
  tlc$fit <- as.numeric(predict(fit_spl))
  library(dplyr)
  s_fit <- tlc |> dplyr::group_by(arm, week) |> dplyr::summarise(m=mean(y), f=mean(fit), .groups="drop")
  library(ggplot2)
  ggplot(s_fit, aes(week, m, color=arm)) +
    geom_point() + geom_line(linewidth=1.1) +
    geom_line(aes(y=f), linetype="dashed", linewidth=1.1) +
    scale_color_manual(values=arm_cols) +
    labs(title="Observed vs fitted means (spline)", y="Mean", x="Week") -> p_fit
  print(p_fit)
}

if (ok_emm) {
  library(emmeans)
  em6 <- emmeans(fit_spl, ~ arm, at = list(week = 6, w1 = 5), data = tlc)
  print(em6)
  print(contrast(em6, "pairwise"))
} else {
  cat("Install emmeans to estimate mean differences at a specific time.\n")
}

if (ok_nlme && ok_emm) {
  # Profiles model with time as factor
  tlc$time_f <- factor(tlc$week)
  fit_prof <- safe_gls(y ~ arm * time_f, data=tlc, corr = corCompSymm(form = ~1|id), method="ML")
  library(emmeans)
  # Curve-implied means at weeks 1 and 6 from the spline model (fit_spl, fit earlier)
  wsel <- c(1,6)
  neww <- expand.grid(week=wsel, arm=levels(tlc$arm))
  neww$w1 <- pmax(neww$week - 1, 0)
  pred_curve <- aggregate(predict(fit_spl, newdata=neww), by=list(week=neww$week, arm=neww$arm), FUN=mean)
  names(pred_curve)[3] <- "mu_curve"
  # Profiles LS-means
  em_prof <- emmeans(fit_prof, ~ arm | time_f, data = tlc)
  summ_prof <- as.data.frame(em_prof)
  summ_prof$week <- as.numeric(as.character(summ_prof$time_f))
  tab_prof <- subset(summ_prof, week %in% wsel, select=c("week","arm","emmean"))
  names(tab_prof)[3] <- "mu_profile"
  comp <- merge(pred_curve, tab_prof, by=c("week","arm"), all=TRUE)
  # Small table
  comp <- comp[order(comp$week, comp$arm), ]
  print(comp)
} else {
  cat("Install nlme + emmeans to compare curve vs profiles means.\n")
}

if (ok_nlme && ok_gg) {
  # Profile LR over candidate knots for TLC-like data
  ks <- c(0.5, 1, 2, 3)
  base_ml <- safe_gls(y ~ arm * week, data=tlc, corr=corCompSymm(form = ~1|id), method="ML")
  LRs <- sapply(ks, function(k) {
  d <- tlc
  d$wkp <- pmax(d$week - k, 0)
  fitk <- safe_gls(y ~ arm * (week + wkp), data = d, corr = corCompSymm(form = ~1|id), method = "ML")
  2*(as.numeric(logLik(fitk)) - as.numeric(logLik(base_ml)))
})
  df <- data.frame(knot=ks, LR=LRs)
  library(ggplot2)
  ggplot(df, aes(knot, LR)) + geom_line() + geom_point() +
    geom_hline(yintercept=qchisq(0.95, df=2), linetype="dashed") +
    labs(title="Profile LR by knot (spline vs linear)", y="LR stat", x="Candidate knot")
} else {
  cat("Install nlme + ggplot2 to visualize LR by knot.\n")
}

tvals <- c(0, 1, 4, 6)
data.frame(t=tvals, hinge = pmax(tvals - 1, 0))