Case-crossover models with adlaplace

The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

Patrick Brown

2026-10-09

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

## Warning: inner_opt status is Radius of trust region is less than
## stop.trust.radius; the Laplace approximation used this inner solution anyway

## 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.0492283849
## 3                  pm10_iwp 0.0010411985           NA -6.8673827804
## 4                 tmean_iwp 0.0090517336           NA -4.7047989815
##    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.3720226591 -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.3351
## At iterate    30  f =       5330.5  |proj g|=        35.953
## At iterate    35  f =       5330.4  |proj g|=        3.2288
## At iterate    40  f =       5330.1  |proj g|=      0.026001
## 
## iterations 43
## function evaluations 71
## 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.293903"

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

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.