Near-monotone mGP case-crossover for London PM10

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. 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).

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)
#> [1]  8.1 99.6
range(london$tmean, na.rm = TRUE)
#> [1] -1.400507 28.174999

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.

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

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
)
#> N = 4, M = 5 machine precision = 2.22045e-16
#> At X0, 0 variables are exactly at the bounds
#> At iterate     0  f=       5369.7  |proj g|=        10889
#> At iterate     5  f =       5366.4  |proj g|=        127.26
#> At iterate    10  f =       5366.3  |proj g|=        1741.9
#> At iterate    15  f =       5363.6  |proj g|=        524.47
#> At iterate    20  f =       5350.8  |proj g|=        39.187
#> At iterate    25  f =       5350.8  |proj g|=        33.787
#> At iterate    30  f =       5350.8  |proj g|=        621.41
#> At iterate    35  f =       5344.5  |proj g|=        621.55
#> 
#> iterations 38
#> function evaluations 50
#> segments explored during Cauchy searches 38
#> BFGS updates skipped 0
#> active bounds at final generalized Cauchy point 0
#> norm of the final projected gradient 3.1147
#> final function value 5343.78
#> 
#> F = 5343.78
#> final  value 5343.777651 
#> converged
fit_mgp$par_info
#>                       label          mle          se  mle_internal se_internal
#> 1                 rh_linear 6.276465e-04 0.000209456  0.0006276465 0.000209456
#> 2                  pm10_mgp 2.203352e-04          NA -8.4203607005 0.982698620
#> 3 all_dirichlet_multinom_sd 4.984877e-05          NA -9.9065167220 0.000000000
#> 4                 tmean_iwp 1.063009e-02          NA -4.5440663004 0.412293683
#>        init     lower upper parscale   log
#> 1  0.000000      -Inf   Inf        1 FALSE
#> 2 -3.912023 -20.72327   Inf        1  TRUE
#> 3 -3.912023      -Inf   Inf        1  TRUE
#> 4 -6.907755 -20.72327   Inf        1  TRUE

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.

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]))
)
#>         gamma max_abs_diff max_abs_diff_hess
#> 1 pm10_mgp_k5  0.001303743       0.002607486

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.

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))
)
#>   parameter max_abs_diff
#> 1  pm10_mgp  0.002090781

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)

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
)
#> N = 4, M = 5 machine precision = 2.22045e-16
#> At X0, 0 variables are exactly at the bounds
#> At iterate     0  f=       5351.6  |proj g|=        10561
#> At iterate     5  f =       5348.6  |proj g|=        122.43
#> At iterate    10  f =       5348.5  |proj g|=        1673.1
#> At iterate    15  f =       5345.8  |proj g|=        454.54
#> At iterate    20  f =       5340.8  |proj g|=        45.111
#> At iterate    25  f =       5340.8  |proj g|=        188.79
#> At iterate    30  f =       5335.7  |proj g|=        8809.7
#> ys=-1.336e+00  -gs= 6.152e-01, BFGS update SKIPPED
#> At iterate    35  f =       5329.4  |proj g|=        43.464
#> 
#> iterations 38
#> function evaluations 57
#> segments explored during Cauchy searches 38
#> BFGS updates skipped 1
#> active bounds at final generalized Cauchy point 0
#> norm of the final projected gradient 0.0956055
#> final function value 5329.4
#> 
#> F = 5329.4
#> final  value 5329.397015 
#> converged
fit_iwp$par_info
#>                       label         mle           se mle_internal  se_internal
#> 1                 rh_linear 0.000605080 0.0002271874   0.00060508 0.0002271874
#> 2                  pm10_iwp 0.001009738           NA  -6.89806416 0.0000000000
#> 3 all_dirichlet_multinom_sd 0.016985913           NA  -4.07537092 0.1158929861
#> 4                 tmean_iwp 0.010380336           NA  -4.56784204 0.4108226114
#>        init     lower upper parscale   log
#> 1  0.000000      -Inf   Inf        1 FALSE
#> 2 -6.907755 -20.72327   Inf        1  TRUE
#> 3 -3.912023      -Inf   Inf        1  TRUE
#> 4 -6.907755 -20.72327   Inf        1  TRUE

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
)
#> N = 4, M = 5 machine precision = 2.22045e-16
#> At X0, 0 variables are exactly at the bounds
#> At iterate     0  f=       5351.7  |proj g|=        10498
#> At iterate     5  f =       5348.7  |proj g|=         119.2
#> At iterate    10  f =       5348.6  |proj g|=        1635.6
#> At iterate    15  f =       5345.7  |proj g|=        217.86
#> At iterate    20  f =       5340.7  |proj g|=         1.209
#> 
#> iterations 21
#> function evaluations 27
#> segments explored during Cauchy searches 21
#> BFGS updates skipped 0
#> active bounds at final generalized Cauchy point 0
#> norm of the final projected gradient 1.20971
#> final function value 5340.67
#> 
#> F = 5340.67
#> final  value 5340.667903 
#> converged
fit_tiwp$par_info
#>                       label          mle           se  mle_internal
#> 1                 rh_linear 0.0006130639 0.0002088531  0.0006130639
#> 2               sqrt_pm_iwp 0.0010012734           NA -6.9064826600
#> 3 all_dirichlet_multinom_sd 0.0025869677           NA -5.9572688641
#> 4                 tmean_iwp 0.0089191784           NA -4.7195514484
#>    se_internal      init     lower upper parscale   log
#> 1 0.0002088531  0.000000      -Inf   Inf        1 FALSE
#> 2 0.0000000000 -6.907755 -20.72327   Inf        1  TRUE
#> 3 0.0000000000 -3.912023      -Inf   Inf        1  TRUE
#> 4 0.3679125344 -6.907755 -20.72327   Inf        1  TRUE

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.

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$

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)
#> Warning in adlaplace::inner_opt(parameters = x_num, gamma = gamma, control =
#> control_inner, : inner_opt: non-positive or non-finite D; half_log_det not
#> computed
#> Warning: Laplace objective was not finite; retrying from gamma = 0
#> Warning: inner_opt status is Radius of trust region is less than
#> stop.trust.radius; the Laplace approximation used this inner solution anyway
logliks <- vapply(fits, function(f) as.numeric(stats::logLik(f)), numeric(1))
data.frame(a = c(1, 2, Inf), logLik = logliks)
#>     a    logLik
#> 1   1 -5350.276
#> 2   2 -5343.772
#> 3 Inf -5343.886

as.numeric(stats::logLik(fit_mgp))
#> [1] -5343.778

Predictive SD at the MLE

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.