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")
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.2,
    parscale = 5
) ~ 0 + rh + iwp(
    pm10,
    p = 2,
    knots = seq(0, 100, by = 10),
    ref_value = 40,
    init = 0.001, parscale = 100
) + iwp(
    tmean,
    p = 3, knots = seq(-2, 30, by = 2),
    init = 0.01, parscale = 100,
    ref_value = 20
)

Fit

fit <- adlaplace(
    formula = cc_formula,
    data = london,
    num_threads = 2L,
    num_groups = 50L,
    config = list(transform_theta = TRUE),
    control = list(
        maxit = 200, trace = 0, REPORT = 1,
        factr = 1e3,
        pgtol = 1e-8
    ),
    verbose = FALSE
)

fit

## Laplace-approximate maximum likelihood fit (adlaplace)
## 
## Call:
## adlaplace(formula = cc_formula, data = london, config = list(transform_theta = TRUE), 
##     num_threads = 2L, num_groups = 50L, control = list(maxit = 200, 
##         trace = 0, REPORT = 1, factr = 1000, pgtol = 1e-08), 
##     verbose = FALSE)
## 
## Coefficients (natural scale):
##                 rh_linear all_dirichlet_multinom_sd                  pm10_iwp 
##                    0.0006                    0.0169                    0.0001 
##                 tmean_iwp 
##                    0.0018 
## 
## Log-likelihood: -5334   df: 4   n: 1825   random effects: 29

summary(fit)

## Laplace-approximate maximum likelihood fit (adlaplace)
## 
## Call:
## adlaplace(formula = cc_formula, data = london, config = list(transform_theta = TRUE), 
##     num_threads = 2L, num_groups = 50L, control = list(maxit = 200, 
##         trace = 0, REPORT = 1, factr = 1000, pgtol = 1e-08), 
##     verbose = FALSE)
## 
## Log-likelihood: -5334   df: 4   n: 1825 
## Convergence code: 0   max|gradient|: 0.0026 
## 
## Coefficients (natural scale, 95% CI):
##                           Estimate Std. Error  lower  upper
## rh_linear                   0.0006      2e-04 0.0002 0.0011
## all_dirichlet_multinom_sd   0.0169         NA 0.0135 0.0213
## pm10_iwp                    0.0001         NA 0.0000 0.0059
## tmean_iwp                   0.0018         NA 0.0008 0.0041
## 
## Random effects: 29 coefficients

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)


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
    )
    rr <- exp(sim_here$sim)
    if (requireNamespace("GET", quietly = TRUE)) {
        q <- GET::central_region(
            GET::create_curve_set(list(obs = rr)),
            coverage = 0.8
        )[, c("lo", "central", "hi")]
    } else {
        q <- t(apply(rr, 1, quantile, probs = c(0.1, 0.5, 0.9)))
    }

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

    matlines(sim_here$x, q, lty = c(2, 1, 2), col = "red")
}

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 = 2L,
    verbose = FALSE
)
config_dev$obs_groups <- obs_groups(
    A = md$term_data$A,
    elgm_matrix = md$term_data$elgm_matrix,
    num_groups = 50L,
    min_groups = 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 = 200, trace = 0, REPORT = 1,
    factr = 1e3,
    pgtol = 1e-8,
    parscale = par_info$parscale
)
control_inner <- list(report.level = 0)
af <- ad_pack(md, config_dev, num_threads = 2L)
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 -1.609438      -Inf   Inf
## 3  pm10                iwp                  pm10_iwp -6.907755 -20.72327   Inf
## 4 tmean                iwp                 tmean_iwp -4.605170 -20.72327   Inf
##   parscale   log
## 1        1 FALSE
## 2        5  TRUE
## 3      100  TRUE
## 4      100  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
)

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] -2785.855 -1890.515 -2484.733 -4613.652 -2210.951

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

## [1]  -6907.58198  -5915.20166    -71.40109    -64.52488 266320.89838

sum(by_shard)

## [1] 6733.595

joint_log_dens(dens_ptr, x_full)

## [1] 6733.595

grad(dens_ptr, x_full)

##  [1]  6.016247e+01  3.017066e+03  7.925580e+02  4.325034e+01  1.816206e-01
##  [6] -9.138977e+03 -5.996974e+03 -3.562910e+03 -1.699437e+03 -5.317709e+02
## [11] -2.426920e+01 -1.589603e+03 -3.087635e+02 -6.337118e+02 -7.072894e+02
## [16] -6.355527e+02 -5.199163e+02 -3.996765e+02 -2.562626e+02 -1.222508e+02
## [21] -4.154798e+01 -5.795445e+00 -8.536926e-02 -5.010426e+02 -1.784458e+02
## [26] -4.580742e+01 -5.072920e+00 -2.398693e-03 -3.239870e+02 -7.960846e+02
## [31]  1.334465e+03  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)
## Beginning optimization
## 
## iter           f           nrm_gr                     status
##   1   6729.011705    1298.802480                 Continuing
##   2   6728.943320      18.431171                 Continuing
##   3   6728.943305       0.004503                 Continuing
##   4   6728.943305       0.000000                 Continuing
## 
## Iteration has terminated
##   4   6728.943305       0.000000                    Success
## 
## trace_hinv_t: threads = 2, shards = 53, params = 33
## trace_hinv_t: finished
## inner_opt: finished

lik$neg_log_lik

## [1] 6888.477

lik$deriv$d_neg_log_lik

## [1] -480.018802 1337.537165    1.464126    3.223225

head(lik$inner_opt$solution)

## [1] -1.697941e-06 -6.855554e-06  6.603516e-07 -1.536926e-08 -1.994934e-05
## [6]  1.777748e-05

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