Case-crossover models with adlaplace

Introduction

This vignette fits a time-stratified case-crossover model to daily mortality and PM10 in London (2002–2006), using the example dataset from Tobias et al. (2024).

Strata are defined by day-of-week within month and year via dirichlet_multinom(..., by = "year_month_dow"). PM10 enters through an integrated Wiener process smooth (iwp(pm10, ...)), with relative humidity as a linear fixed effect. Fitting is done with adlaplace().

Data

data("london", package = "adlaplace")
if (!vignette_full) {
  london <- london[london$year == 2003, ]
}
head(london)

##     city       date year month day dow all       tmean      rh     pm10
## 1 London 2002-01-01 2002   Jan   1 Tue 199 -0.22500001 75.6783 71.70000
## 2 London 2002-01-02 2002   Jan   2 Wed 231  0.08749995 77.5250 40.20000
## 3 London 2002-01-03 2002   Jan   3 Thu 210  0.85000002 81.3250 41.80000
## 4 London 2002-01-04 2002   Jan   4 Fri 203  0.53750002 85.4500 50.40000
## 5 London 2002-01-05 2002   Jan   5 Sat 224  4.25000000 93.5250 49.36429
## 6 London 2002-01-06 2002   Jan   6 Sun 198  7.06874990 96.4000 31.10000

quantile(london$pm10, na.rm = TRUE)

##       0%      25%      50%      75%     100% 
##  8.10000 30.70000 39.72857 50.50000 99.60000

quantile(london$tmean, na.rm = TRUE)

##        0%       25%       50%       75%      100% 
## -1.400507  7.506250 11.472748 16.198438 28.174999

Model specification

Case-crossover strata absorb the intercept, so the formula uses 0 + on the right-hand side.

london$temp_div10 <- london$tmean / 10
cc_formula <- dirichlet_multinom(
  all,
  by = c("year", "month", "dow"),
  init = 0.02
) ~ 0 + rh + iwp(
  pm10,
  p = 2,
  knots = seq(0, 100, by = 10),
  ref_value = 40,
  init = 0.001
) + iwp(
  tmean,
  p = 3, knots = seq(-2, 30, by = 2),
  init = 0.001,
  ref_value = 20
)

Fit

fit <- adlaplace(
  formula = cc_formula,
  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=       5352.3  |proj g|=        10604
## At iterate     5  f =       5349.2  |proj g|=        121.32
## At iterate    10  f =       5349.1  |proj g|=          1659
## At iterate    15  f =       5346.5  |proj g|=        438.17
## At iterate    20  f =       5341.5  |proj g|=        8.1316
## 
## iterations 21
## function evaluations 30
## 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.00993
## final function value 5341.48
## 
## F = 5341.48
## final  value 5341.481082 
## converged

fit$par_info

##                       label          mle           se  mle_internal
## 1                 rh_linear 0.0006178345 0.0002083658  0.0006178345
## 2 all_dirichlet_multinom_sd 0.0023596821           NA -6.0492283834
## 3                  pm10_iwp 0.0010411985           NA -6.8673827805
## 4                 tmean_iwp 0.0090517336           NA -4.7047989832
##    se_internal      init     lower upper parscale   log
## 1 0.0002083658  0.000000      -Inf   Inf        1 FALSE
## 2 0.0000000000 -3.912023      -Inf   Inf        1  TRUE
## 3 0.0000000000 -6.907755 -20.72327   Inf        1  TRUE
## 4 0.3720226588 -6.907755 -20.72327   Inf        1  TRUE

The same model with overdispersion optimized on the natural scale (log = FALSE) should recover the same MLE and log-likelihood.

cc_formula_nat <- dirichlet_multinom(
  all,
  by = c("year", "month", "dow"),
  init = 0.015,
  lower = 0,
  upper = 0.5,
  parscale = 1e-3,
  log = FALSE
) ~ 0 + rh + iwp(
  pm10,
  p = 2,
  knots = seq(0, 100, by = 10),
  ref_value = 40,
  init = 0.0001
) + iwp(
  tmean,
  p = 3, knots = seq(-2, 30, by = 2),
  init = 0.002,
  ref_value = 20
)

fit_nat <- adlaplace(
  formula = cc_formula_nat,
  data = london,
  config = list(
    num_threads = vignette_num_threads,
    num_shards = 100L
  ),
  control = list(
    maxit = fit_maxit,
    trace = 3, REPORT = 5,
    factr = 1e3,
    pgtol = 1e-8
  ),
  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=       5345.1  |proj g|=        12225
## At iterate     5  f =       5341.4  |proj g|=        105.58
## At iterate    10  f =       5341.3  |proj g|=        1492.9
## At iterate    15  f =       5335.5  |proj g|=         11033
## At iterate    20  f =       5330.5  |proj g|=        2.0286
## At iterate    25  f =       5330.5  |proj g|=        4.3353
## At iterate    30  f =       5330.5  |proj g|=        35.951
## At iterate    35  f =       5330.4  |proj g|=        3.3148
## At iterate    40  f =       5330.1  |proj g|=      0.027709
## 
## iterations 43
## function evaluations 60
## segments explored during Cauchy searches 43
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.000377107
## final function value 5330.09
## 
## F = 5330.09
## final  value 5330.091906 
## converged

all.equal(coef(fit), coef(fit_nat), tolerance = 1e-5)

## [1] "Mean relative difference: 1.293902"

all.equal(as.numeric(logLik(fit)), as.numeric(logLik(fit_nat)), tolerance = 1e-5)

## [1] "Mean relative difference: 0.002132213"

num_sim <- 200L
num_grid <- 101L
gamma_sims <- rmvnldl(n = num_sim, fit = fit)

smooth_vars <- unique(fit$model_data$term_data$info$gamma$term)

oldpar <- par(no.readonly = TRUE)
par(mfrow = c(1, 2), mar = c(2.25, 2.25, 0, 0), bty = "l", mgp = c(1, 0.5, 0))
for (D in smooth_vars) {
  sim_here <- sim_random(
    D,
    fit = fit,
    gamma_sims = gamma_sims,
    num_grid = num_grid,
    transform = exp
  )
  q <- if (!is.null(sim_here$global)) sim_here$global else sim_here$pointwise

  matplot(
    sim_here$x, sim_here$sim,
    type = "l", lty = 1, col = "#00000010",
    xlab = D, ylab = "RR", xaxs = "i", yaxs = "i",
    ylim = quantile(sim_here$sim, prob = c(0.005, 0.95))
  )
  matlines(sim_here$x, q, lty = c(2, 1, 2), col = "red")
}

par(oldpar)

AD development check

Build model_data, shard map, and ad_pack without running outer optimization, then maximize the Laplace profile likelihood with optim directly.

md <- model_data(cc_formula, data = london, verbose = FALSE)
config_dev <- list(
  transform_theta = TRUE,
  num_threads = vignette_num_threads,
  verbose = FALSE
)
config_dev$obs_groups <- obs_groups(
  A = md$term_data$A,
  elgm_matrix = md$term_data$elgm_matrix,
  num_shards = 50L,
  min_shards = min(50L, 2L * 4L)
)
par_info <- md$term_data$info$parameters
cache <- new.env(parent = emptyenv())
cache$gamma <- rep(0, nrow(md$term_data$info$gamma))
control_dev <- list(
  maxit = 10, trace = 3, REPORT = 5,
  parscale = par_info$parscale
)
control_inner <- list(report.level = 0)
af <- ad_pack(md, config_dev, num_threads = vignette_num_threads)
slotNames(af)

##  [1] "ptr"             "group_sparsity"  "outer"           "inner"          
##  [5] "map_outer"       "map_inner"       "parallel_map"    "chol_inner"     
##  [9] "chol_inner_list" "sizes"           "info"

dim(config_dev$obs_groups)

## [1] 420  50

par_info

##    term              model                     label      init     lower upper
## 1    rh             linear                 rh_linear  0.000000      -Inf   Inf
## 2   all dirichlet_multinom all_dirichlet_multinom_sd -3.912023      -Inf   Inf
## 3  pm10                iwp                  pm10_iwp -6.907755 -20.72327   Inf
## 4 tmean                iwp                 tmean_iwp -6.907755 -20.72327   Inf
##   parscale   log
## 1        1 FALSE
## 2        1  TRUE
## 3        1  TRUE
## 4        1  TRUE

bob <- stats::optim(
  par = par_info$init,
  fn = outer_fn,
  gr = outer_gr,
  method = "L-BFGS-B",
  control = control_dev,
  lower = par_info$lower,
  upper = par_info$upper,
  config = config_dev,
  ad_pack = af,
  cache = cache,
  control_inner = control_inner,
  hessian = TRUE
)

## N = 4, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=       5352.3  |proj g|=        10604
## At iterate     5  f =       5349.2  |proj g|=        121.32
## At iterate    10  f =       5349.1  |proj g|=          1659
## final  value 5348.990265 
## stopped after 11 iterations

Joint log density

At starting values, evaluate the full log-density and its gradient over all parameters (beta, inner gamma, and theta).

sz <- af@sizes
gamma0 <- rep(0, nrow(md$term_data$info$gamma))
x_full <- c(
  par_info$init[seq_len(sz["beta"])],
  gamma0,
  par_info$init[seq(sz["beta"] + 1L, length.out = sz["theta"])]
)

# Serial dens/grad/hess on a dens-safe clone (multi-thread affinity cleared).
dens_ptr <- clone_ad_pack_ptr(af@ptr)
shards <- seq.int(from = 0L, length.out = adlaplace:::n_groups(dens_ptr))
by_shard <- vapply(
  shards,
  function(s) joint_log_dens(dens_ptr, x_full, ad_shards = s),
  numeric(1)
)
by_shard[1:5]

## [1] 4386.024 2790.146 4277.597 6248.223 3486.712

by_shard[seq(to = length(by_shard), length.out = 5)]

## [1]    9445.66083    8740.04252     -59.88817     -95.82107 -363562.89841

sum(by_shard)

## [1] 5271.516

joint_log_dens(dens_ptr, x_full)

## [1] 5271.516

grad(dens_ptr, x_full)

##  [1]  6.734860e+02  6.169219e+02  1.701712e+02  1.065871e+01  3.620924e-02
##  [6] -2.058877e+03 -1.421416e+03 -8.792575e+02 -4.298762e+02 -1.374054e+02
## [11] -7.191833e+00 -3.433567e+03 -1.181943e+03 -1.960381e+03 -2.089588e+03
## [16] -1.849754e+03 -1.501943e+03 -1.151819e+03 -7.451234e+02 -3.623431e+02
## [21] -1.257438e+02 -1.793431e+01 -2.509312e-01 -1.504503e+03 -5.623891e+02
## [26] -1.529105e+02 -1.805906e+01 -8.759916e-03 -3.446074e+03 -5.494638e+03
## [31] -3.552390e+01  1.000000e+01  1.600000e+01

bob <- hessian(dens_ptr, x_full)

Likelihood

With outer parameters fixed at their starting values, log_lik_laplace() optimizes random effects and returns the log-likelihood.

x_outer <- par_info$init

lik <- log_lik_laplace(
  x = x_outer,
  ad_pack = af,
  config = modifyList(config_dev, list(verbose = TRUE)),
  gamma = gamma0,
  control = list(report.level = 2, report.freq = 1),
  deriv = TRUE
)

## inner_opt: threads = 2, shards = 53, params = 33 (beta = 1, gamma = 29, theta = 3)
## inner_opt: constructing AD_Func_Opt (hess maps)...
## inner_opt: AD_Func_Opt done nnz=755
## inner_opt: entering CppadParallelScope...
## CppadParallelScope: drain shard eval buffers (entry)...
## CppadParallelScope: setup num_threads=2
## CppadParallelScope: setup done
## inner_opt get_fdf: before OpenMP threads=2 Nparams=33 shards=53
## inner_opt get_fdf: after OpenMP
## inner_opt get_fdfh: before OpenMP threads=2 Nparams=33 hess_size=755 shards=53
## inner_opt get_fdfh parallel_team=2
## inner_opt get_fdfh: after OpenMP api_err_shard=-1
## inner_opt get_fdfh: scatter Hessian...
## inner_opt get_fdfh: scatter done
## inner_opt: opt.run() starting (first get_fdfh)...
## Beginning optimization
## inner_opt get_f: before OpenMP threads=2 shards=53
## inner_opt get_f: after OpenMP
## inner_opt get_fdf: before OpenMP threads=2 Nparams=33 shards=53
## inner_opt get_fdf: after OpenMP
## 
## iter           f          nrm_gr                     status
##   1   5178.319260    254.593240                 Continuing
## inner_opt get_fdfh: before OpenMP threads=2 Nparams=33 hess_size=755 shards=53
## inner_opt get_fdfh parallel_team=2
## inner_opt get_fdfh: after OpenMP api_err_shard=-1
## inner_opt get_fdfh: scatter Hessian...
## inner_opt get_fdfh: scatter done
## inner_opt get_f: before OpenMP threads=2 shards=53
## inner_opt get_f: after OpenMP
## inner_opt get_fdf: before OpenMP threads=2 Nparams=33 shards=53
## inner_opt get_fdf: after OpenMP
##   2   5178.248523      0.158533                 Continuing
## inner_opt get_fdfh: before OpenMP threads=2 Nparams=33 hess_size=755 shards=53
## inner_opt get_fdfh parallel_team=2
## inner_opt get_fdfh: after OpenMP api_err_shard=-1
## inner_opt get_fdfh: scatter Hessian...
## inner_opt get_fdfh: scatter done
## inner_opt get_f: before OpenMP threads=2 shards=53
## inner_opt get_f: after OpenMP
## inner_opt get_fdf: before OpenMP threads=2 Nparams=33 shards=53
## inner_opt get_fdf: after OpenMP
##   3   5178.248523      0.000001                 Continuing
## 
## Iteration has terminated
##   3   5178.248523      0.000001                    Success
## 
## inner_opt get_fdf: before OpenMP threads=2 Nparams=33 shards=53
## inner_opt get_fdf: after OpenMP
## inner_opt get_fdfh: before OpenMP threads=2 Nparams=33 hess_size=755 shards=53
## inner_opt get_fdfh parallel_team=2
## inner_opt get_fdfh: after OpenMP api_err_shard=-1
## inner_opt get_fdfh: scatter Hessian...
## inner_opt get_fdfh: scatter done
## inner_opt: opt.run() done iters=3 fval=5178.248523
## inner_opt: chol_update starting...
## inner_opt: chol_update done log_det=401.365487
## inner_opt: trust region done; outer get_fdfh next (2 threads)
## outer get_fdfh: before OpenMP threads=2 Nparams=33 hess_size=929 shards=53
## outer get_fdfh parallel_team=2
## outer get_fdfh: after OpenMP api_err_shard=-1
## outer get_fdfh: scatter Hessian...
## outer get_fdfh: scatter done
## inner_opt: outer get_fdfh done
## inner_opt: trace_hinv_t next (2 threads)
## trace_hinv_t: threads = 2, shards = 53, params = 33
## trace_hinv_t: finished
## inner_opt: leaving CppadParallelScope (teardown next)...
## CppadParallelScope: drain shard eval buffers (exit)...
## CppadParallelScope: teardown (free_available + free_all)...
## CppadParallelScope: teardown done
## inner_opt: parallel block ended
## inner_opt: finished

lik$neg_log_lik

## [1] 5352.282

lik$deriv$d_neg_log_lik

## [1] -1.060442e+04  8.159608e+00 -1.261264e-01 -9.724369e+00

head(lik$inner_opt$solution)

## [1]  6.175304e-05 -1.200242e-05 -2.041847e-07 -4.323359e-08  1.271120e-04
## [6]  2.751351e-04

neg_log_lik is the negative profile log-likelihood after inner optimization over random effects.