---
title: "Near-monotone mGP case-crossover for London PM10"
author: "Patrick Brown"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{Near-monotone mGP case-crossover for London PM10}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

## Background

Daily mortality and PM10 in London (2002–2006) often show a near-monotone
exposure–response: risk rises with pollution, with a base shape close to
$\sqrt{\mathrm{PM}_{10}}$. The near-monotone Gaussian process (mGP) of
[Zhang et al.](https://github.com/AgueroZZ/near_monotone_smoothing) places a
spline penalty on the departure from that base model (Ramsay, 1998), while
keeping $f(x_0)=f'(x_0)=0$ so case-crossover strata can absorb the intercept.

This vignette fits a time-stratified case-crossover model with
`mgp(pm10, a = 2)` and compares it to `iwp(pm10, p = 2)`.

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5
)
library(adlaplace)
library(adlaplaceFem)
if (file.exists("vignette_full.R")) {
  source("vignette_full.R", local = TRUE)
} else {
  vignette_full <- TRUE
  vignette_num_threads <- 8L
}
fit_num_shards <- if (vignette_full) 200L else 4L
fit_maxit <- if (vignette_full) 200L else 10L
```


## Data

```{r data}
data("london", package = "adlaplace")
london <- london[!is.na(london$pm10), ]
if (!vignette_full) {
  london <- london[london$year == 2003, ]
}
london$sqrt_pm <- sqrt(london$pm10)
range(london$pm10)
range(london$tmean, na.rm = TRUE)
```

PM10 ranges from about 8.1 to 99.6 $\mu\mathrm{g}/\mathrm{m}^3$.

## Model

Strata absorb the intercept (`0 +` on the RHS). Bare pre-built term objects
are not formula-safe (they become `linear()`), so constructors are written
inline. Copies of the terms are kept for later `psd()` plots.

```{r model-spec}
pm10_iwp <- iwp(
  "pm10",
  p = 2,
  knots = seq(5, 105, by = 20),
  ref_value = 25,
  init = 0.001
)
pm10_mgp <- mgp("pm10",
  a = 2, knots = pm10_iwp[[1]]@knots,
  ref_value = pm10_iwp[[1]]@ref_value
)
pm10_tiwp <- iwp("sqrt_pm",
  p = 2,
  knots = seq(1, 11, by = 2),
  ref_value = 5, init = 0.001
)

cc_base <- collect_terms(
  dirichlet_multinom(
    all,
    by = c("year", "month", "dow"),
    init = 0.02
  ) ~ 0 + rh +
    iwp(tmean, p = 3, knots = seq(-2, 30, by = 2), ref_value = 20, init = 0.001)
)
```




## Fit

```{r fit}
fit_mgp <- adlaplace(
  formula = c(pm10_mgp, cc_base),
  data = london,
  config = list(
    num_threads = vignette_num_threads,
    num_shards = fit_num_shards
  ),
  control = list(maxit = fit_maxit, trace = 3, REPORT = 5),
  verbose = FALSE,
  hessian = vignette_full
)
fit_mgp$par_info
```

The inner solver can stop because the trust radius falls below
`stop.trust.radius` and still hand that $\gamma$ to the Laplace
approximation. That solver differentiates the joint log density, so compare
its derivative to first differences along the random-effect coordinate with
the largest gradient at the fitted point.

```{r joint-deriv, fig.height = 4, fig.width = 9, eval = vignette_full}
dens_ptr <- clone_ad_pack_ptr(fit_mgp$ad_pack@ptr)
x_here <- as.numeric(fit_mgp$details$full_parameters)
g_here <- grad(dens_ptr, x_here, negative = FALSE)
n_beta <- as.integer(fit_mgp$ad_pack@sizes[["beta"]])
n_gamma <- as.integer(fit_mgp$ad_pack@sizes[["gamma"]])
gamma_idx <- seq.int(n_beta + 1L, length.out = n_gamma)
Dpar <- gamma_idx[which.max(abs(g_here[gamma_idx]))]
gamma_label <- names(fit_mgp$gamma)[Dpar - n_beta]

Dpar <- length(x_here)

scan_j <- deriv_scan_joint(
  dens_ptr,
  x = x_here,
  index = Dpar,
  width = 0.1,
  n = 7L,
  negative = FALSE
)

grad_at_mid <- (
  scan_j$grad[-1L, Dpar] + scan_j$grad[-nrow(scan_j$grad), Dpar]
) / 2
hess_at_mid <- (
  scan_j$hess_row[-1L, Dpar] + scan_j$hess_row[-nrow(scan_j$hess_row), Dpar]
) / 2
data.frame(
  gamma = gamma_label,
  max_abs_diff = max(abs(grad_at_mid - scan_j$fd_value)),
  max_abs_diff_hess = max(abs(hess_at_mid - scan_j$fd_grad[, Dpar]))
)

oldpar <- par(no.readonly = TRUE)
par(mfrow = c(1, 3), mar = c(3, 3, 2, 0.5), mgp = c(1.6, 0.5, 0))
plot(
  scan_j$grid, scan_j$value,
  type = "b",
  xlab = gamma_label, ylab = "joint log density"
)
plot(
  scan_j$grid, scan_j$grad[, Dpar],
  type = "l",
  xlab = gamma_label, ylab = "derivative"
)
points(scan_j$mid, scan_j$fd_value, pch = 16)
legend(
  "topright",
  legend = c("AD", "first difference"),
  lty = c(1, NA), pch = c(NA, 16), bty = "n"
)
plot(
  scan_j$grid, scan_j$hess_row[, Dpar],
  type = "l",
  xlab = gamma_label, ylab = "second derivative",
  ylim = range(
    c(scan_j$hess_row[, Dpar], scan_j$fd_grad[, Dpar]),
    na.rm = TRUE
  )
)
points(scan_j$mid, scan_j$fd_grad[, Dpar], pch = 16)
par(oldpar)
```

The outer objective is the Laplace approximation to the log likelihood.
The same first-difference check, along the PM10 mGP standard deviation,
compares that profile derivative to the log likelihood itself.

```{r loglik-deriv, fig.height = 4, eval = vignette_full}
x_outer <- fit_mgp$par_info$mle_internal
Dpar_ll <- match("pm10_mgp", fit_mgp$par_info$label)
scan_label <- fit_mgp$par_info$label[Dpar_ll]

scan_ll <- deriv_scan_laplace(
  fit_mgp$ad_pack,
  x = x_outer,
  index = Dpar_ll,
  width = 0.2,
  n = 5L,
  config = fit_mgp$config,
  gamma = as.numeric(fit_mgp$gamma),
  control = fit_mgp$control_inner
)

grad_ll_mid <- (
  scan_ll$grad[-1L, Dpar_ll] + scan_ll$grad[-nrow(scan_ll$grad), Dpar_ll]
) / 2
data.frame(
  parameter = scan_label,
  max_abs_diff = max(abs(grad_ll_mid - scan_ll$fd_value))
)

oldpar <- par(no.readonly = TRUE)
par(mfrow = c(1, 2), mar = c(3, 3, 2, 0.5), mgp = c(1.6, 0.5, 0))
plot(
  scan_ll$grid, scan_ll$value,
  type = "b",
  xlab = scan_label, ylab = "negative log likelihood"
)
plot(
  scan_ll$grid, scan_ll$grad[, Dpar_ll],
  type = "l",
  xlab = scan_label, ylab = "derivative"
)
points(scan_ll$mid, scan_ll$fd_value, pch = 16)
legend(
  "topright",
  legend = c("AD", "first difference"),
  lty = c(1, NA), pch = c(NA, 16), bty = "n"
)
par(oldpar)
```

```{r fit-iwp, eval = vignette_full}
fit_iwp <- adlaplace(
  formula = c(pm10_iwp, cc_base),
  data = london,
  config = list(
    num_threads = vignette_num_threads,
    num_shards = fit_num_shards
  ),
  control = list(maxit = fit_maxit, trace = 3, REPORT = 5),
  verbose = FALSE,
  hessian = vignette_full
)
fit_iwp$par_info
```

```{r fit-tiwp, eval = vignette_full}
fit_tiwp <- adlaplace(
  formula = c(pm10_tiwp, cc_base),
  data = london,
  config = list(
    num_threads = vignette_num_threads,
    num_shards = fit_num_shards
  ),
  control = list(maxit = fit_maxit, trace = 3, REPORT = 5),
  verbose = FALSE,
  hessian = vignette_full
)
fit_tiwp$par_info
```

## Relative risk vs PM10

`sim_random()` needs `fit$gamma` and `fit$details$hessian$chol_inner`.
Those are always filled by the post-fit Laplace evaluation, but the curves
below are full-build only (abbreviated `R CMD check` stops after `par_info`).
The grey histogram is the distribution of PM10 in the case-crossover data.

```{r plot-rr, eval = vignette_full}
num_sim <- 200L
# Pass new_x explicitly: pretty() can step outside the knot range for some n.
pm10_grid <- data.frame(pm10 = seq(8, 100, length.out = 81L))
set.seed(1)
sim_mgp <- sim_random(
  "pm10",
  fit = fit_mgp,
  new_x = pm10_grid,
  num_sim = num_sim,
  transform = exp
)
set.seed(1)
sim_iwp <- sim_random(
  "pm10",
  fit = fit_iwp,
  new_x = pm10_grid,
  num_sim = num_sim,
  transform = exp
)
sim_tiwp <- sim_random(
  "sqrt_pm",
  fit = fit_tiwp,
  new_x = data.frame(sqrt_pm = sqrt(pm10_grid$pm10)),
  num_sim = num_sim,
  transform = exp
)
rr_ylim <- range(sim_mgp$pointwise, sim_iwp$pointwise)
oldpar <- par(no.readonly = TRUE)
par(mar = c(3, 3, 1, 1), bty = "l", mgp = c(1.8, 0.5, 0))
plot(
  NA, NA,
  xlim = range(pm10_grid$pm10), ylim = rr_ylim,
  xlab = expression(PM[10]), ylab = "RR"
)
hx <- london$pm10
hx <- hx[is.finite(hx) & hx >= min(pm10_grid$pm10) & hx <= max(pm10_grid$pm10)]
if (length(hx) > 10L) {
  hb <- hist(hx, breaks = 30, plot = FALSE)
  dens <- hb$counts / max(hb$counts)
  y0 <- rr_ylim[1]
  y1 <- y0 + 0.12 * diff(rr_ylim)
  rect(
    hb$breaks[-length(hb$breaks)], y0,
    hb$breaks[-1L], y0 + dens * (y1 - y0),
    col = adjustcolor("grey50", alpha.f = 0.25),
    border = NA
  )
}
lines(sim_mgp$x, sim_mgp$pointwise[, "0.5"], lwd = 2, col = "#1b4f72")
matlines(sim_mgp$x, sim_mgp$pointwise[, c("0.1", "0.9")], lty = 2, col = "#1b4f72")
lines(sim_iwp$x, sim_iwp$pointwise[, "0.5"], lwd = 2, col = "#b03a2e")
matlines(sim_iwp$x, sim_iwp$pointwise[, c("0.1", "0.9")], lty = 2, col = "#b03a2e")
matlines(sim_tiwp$x^2, sim_tiwp$pointwise[, c("0.1", "0.5", "0.9")],
  lty = c(2, 1, 2), col = "green"
)
legend(
  "topleft",
  legend = c("mGP (a = 2)", "IWP", "TIWP"),
  col = c("#1b4f72", "#b03a2e", "green"),
  lwd = 2,
  bty = "n"
)
par(oldpar)
```

## Sensitivity in $a$

```{r sensitivity, eval = vignette_full}
fit_a <- function(a) {
  adlaplace(
    c(
      mgp("pm10",
        a = a, knots = pm10_mgp[[1]]@knots,
        ref_value = pm10_mgp[[1]]@ref_value
      ),
      cc_base
    ),
    data = london,
    config = list(num_threads = vignette_num_threads, num_shards = fit_num_shards),
    control = list(maxit = fit_maxit),
    verbose = FALSE,
    hessian = FALSE
  )
}
fits <- lapply(c(1, 2, Inf), fit_a)
logliks <- vapply(fits, function(f) as.numeric(stats::logLik(f)), numeric(1))
data.frame(a = c(1, 2, Inf), logLik = logliks)

as.numeric(stats::logLik(fit_mgp))
```


## Predictive SD at the MLE

```{r psd-plots, eval = vignette_full}
sd_col <- if ("mle" %in% names(fit_mgp$par_info)) "mle" else "init"
sd_mgp <- fit_mgp$par_info[[sd_col]][fit_mgp$par_info$label == "pm10_mgp"]
sd_iwp <- fit_iwp$par_info[[sd_col]][fit_iwp$par_info$label == "pm10_iwp"]

h <- 5
x_grid <- seq(8, 90, length.out = 40)
psd_mgp_x <- psd(pm10_mgp$pm10_mgp, h = h, x = x_grid, sd = sd_mgp)
psd_iwp_x <- rep(psd(pm10_iwp$pm10_iwp, h = h, sd = sd_iwp), length(x_grid))

oldpar <- par(no.readonly = TRUE)
par(mar = c(3, 3, 1.5, 1), bty = "l", mgp = c(1.8, 0.5, 0))
plot(
  x_grid, psd_mgp_x,
  type = "l", lwd = 2, col = "#1b4f72",
  xlab = expression(PM[10]), ylab = expression(sigma[x](5)),
  ylim = range(psd_mgp_x, psd_iwp_x)
)
lines(x_grid, psd_iwp_x, lwd = 2, col = "#b03a2e")
legend("topleft", c("mGP", "IWP"), col = c("#1b4f72", "#b03a2e"), lwd = 2, bty = "n")
par(oldpar)
```

The mGP predictive SD depends on $x$ (nonstationary), while the IWP predictive
SD is constant.