Case-crossover models with hpolcc

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.

Data

london <- fread(system.file("extdata", "london.csv", package = "hpolcc"))
london$date <- as.Date(ISOdate(london$year, london$month, london$day))
london$year_month_dow <- format(london$date, "%Y-%m-%a")
head(london)

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

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

london$temp_div10 <- london$tmean / 10
cc_formula <- dirichlet_multinom(
    all,
    by = "year_month_dow",
    init = 0.2,
    parscale = 5
) ~ 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 <- hnlm(
    formula = cc_formula,
    data = london,
    config = list(
        transform_theta = TRUE,
        num_threads = 2L,
        num_shards = 50L,
        verbose = FALSE
    ),
    control = list(
        maxit = 200, trace = 0, REPORT = 1,
        factr = 1e3,
        pgtol = 1e-8
    )
)

fit

## hnlm fit
##   formula: cc_formula 
##   n = 420 
##   log-lik: -5334.4 
##   converged: yes 
##   beta:
##      label          mle
##  rh_linear 0.0006173671
##   theta:
##                      label          mle
##  all_dirichlet_multinom_sd 0.0169254661
##                   pm10_iwp 0.0001281147
##                  tmean_iwp 0.0018065549

summary(fit)

## 
## Call:
## hnlm(formula = cc_formula, data = london, config = list(transform_theta = TRUE, 
##     num_threads = 2L, num_shards = 50L, verbose = FALSE), control = list(maxit = 200, 
##     trace = 0, REPORT = 1, factr = 1000, pgtol = 1e-08))
## 
## Log-likelihood: -5334 
## Observations: 420 
## Converged: yes 
## 
## Parameters:
##                      label  term              model       mle        se
##                  rh_linear    rh             linear 0.0006174 0.0002286
##  all_dirichlet_multinom_sd   all dirichlet_multinom 0.0169300 0.0019800
##                   pm10_iwp  pm10                iwp 0.0001281 0.0002505
##                  tmean_iwp tmean                iwp 0.0018070 0.0007461

par(mfrow = c(1, 2), mar = c(2.25, 2.25, 0, 0), bty = "l", mgp = c(1, 0.5, 0))
for (D in names(fit$sample$sim)) {
    the_xlim <- range(london[[D]], na.rm = TRUE)
    sample_in_xlim <- fit$sample$sim[[D]][
        fit$sample$x[[D]] > min(the_xlim) &
            fit$sample$x[[D]] < max(the_xlim),
    ]
    matplot(fit$sample$x[[D]], fit$sample$sim[[D]],
        type = "l", lty = 1, col = "#00000010",
        xlab = D, ylab = "RR", xaxs = "i", yaxs = "i",
        xlim = the_xlim,
        ylim = quantile(sample_in_xlim, prob = c(0.005, 0.99))
    )
    matlines(fit$sample$x[[D]], fit$sample$envelope$common[[D]], lty = c(2, 1, 2), col = "red")
}

AD development check

for_dev = TRUE builds model_data, the combined ad_fun, and config without running outer optimization.

config_dev <- list(
    transform_theta = TRUE,
    num_threads = 2L,
    num_shards = 50L,
    verbose = FALSE
)
forres <- hnlm(
    cc_formula,
    data = london,
    config = config_dev,
    for_dev = TRUE
)
slotNames(forres$ad_fun)

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

dim(forres$config$shards)

## [1] 420  50

forres$model_data$data$info$parameters

##    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 transform
## 1        1     FALSE
## 2        5      TRUE
## 3      100      TRUE
## 4      100      TRUE

bob <- stats::optim(
    par = forres$config$opt$init,
    fn = adlaplace::outer_fn,
    gr = adlaplace::outer_gr,
    method = "L-BFGS-B",
    control = forres$control,
    lower = forres$config$opt$lower,
    upper = forres$config$opt$upper,
    config = forres$config,
    ad_fun = forres$ad_fun,
    cache = forres$cache,
    control_inner = forres$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=       6888.5  |proj g|=       6687.7
## At iterate     1  f =       6569.7  |proj g|=        6200.3
## At iterate     2  f =       5443.1  |proj g|=         45112
## At iterate     3  f =       5360.5  |proj g|=        7030.4
## At iterate     4  f =       5347.5  |proj g|=          5890
## At iterate     5  f =       5343.7  |proj g|=        2944.3
## At iterate     6  f =       5341.6  |proj g|=        394.61
## At iterate     7  f =       5341.2  |proj g|=        263.86
## At iterate     8  f =       5336.8  |proj g|=        53.006
## At iterate     9  f =       5334.8  |proj g|=        1704.6
## At iterate    10  f =       5334.7  |proj g|=        179.67
## At iterate    11  f =       5334.7  |proj g|=        57.641
## At iterate    12  f =       5334.7  |proj g|=        47.978
## At iterate    13  f =       5334.7  |proj g|=        65.761
## At iterate    14  f =       5334.6  |proj g|=        139.58
## At iterate    15  f =       5334.5  |proj g|=        164.76
## At iterate    16  f =       5334.5  |proj g|=        280.28
## At iterate    17  f =       5334.4  |proj g|=        218.44
## At iterate    18  f =       5334.4  |proj g|=        33.008
## At iterate    19  f =       5334.4  |proj g|=        2.6841
## At iterate    20  f =       5334.4  |proj g|=        9.7068
## At iterate    21  f =       5334.4  |proj g|=        3.0484
## At iterate    22  f =       5334.4  |proj g|=        10.847
## 
## iterations 22
## function evaluations 32
## segments explored during Cauchy searches 24
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 10.8475
## final function value 5334.39
## 
## F = 5334.39
## final  value 5334.389277 
## converged

Joint log density

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

sz <- forres$ad_fun@sizes
x_full <- c(
    forres$config$opt$init[seq_len(sz["beta"])],
    forres$cache$gamma,
    forres$config$opt$init[seq(sz["beta"] + 1L, length.out = sz["theta"])]
)

shards <- seq.int(from = 0L, length.out = adlaplace:::n_groups(forres$ad_fun@ptr))
by_shard <- vapply(
    shards,
    function(s) joint_log_dens(forres$ad_fun, x_full, shards = s),
    numeric(1)
)
by_shard[1:5]

## [1]  -2785.857  -2459.259  -4613.676 -10358.298  -6973.934

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

## [1]  -4845.69226  -4576.26136    -71.39340    -64.47563 266320.89838

sum(by_shard)

## [1] 6729.06

joint_log_dens(forres$ad_fun, x_full)

## [1] 6729.06

grad(forres$ad_fun, x_full)

##  [1] -5.892584e+02  5.264055e+01  7.987757e+01  6.527927e-01  1.307021e-01
##  [6]  1.049319e+03  5.248790e+02  5.808055e+01 -2.476329e+01 -3.767879e+01
## [11] -1.888044e-01  1.055043e+02  1.216595e+02  1.647132e+02  1.587220e+02
## [16]  1.348805e+02  1.001427e+02  6.204113e+01  3.397821e+01  1.524148e+01
## [21]  3.630802e+00  3.203617e-01 -1.713720e-03  7.487946e+01  2.425810e+01
## [26]  5.833047e+00  6.533312e-01  2.126975e-04  5.135125e+01  1.958144e+02
## [31]  1.344413e+03  9.984609e+00  1.590149e+01

bob <- hessian(forres$ad_fun, x_full)

Likelihood

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

x_outer <- forres$model_data$data$info$parameters$init

lik <- log_lik_laplace(
    x = x_outer,
    ad_fun = forres$ad_fun,
    config = modifyList(forres$config, list(verbose = TRUE)),
    gamma = forres$cache$gamma,
    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   6728.943312     13.858325                 Continuing
##   2   6728.943305      0.001874                 Continuing
##   3   6728.943305      0.000000                 Continuing
## 
## Iteration has terminated
##   3   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$grad

## [1] -480.018802 1337.537165    1.464126    3.223225

head(lik$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.

Conditional simulations of IWP smooths use adlaplaceHgp::cond_sim_iwp() with the Laplace output and model-data bundle:

sim <- adlaplaceHgp::cond_sim_iwp(
    fit = lik,
    model_data = forres$model_data,
    n = 25
)
names(sim)

## [1] "quantiles" "envelope"  "x"         "gamma"     "sim"

# inspect assembled inputs:
# str(adlaplaceHgp::cond_sim_iwp_inputs(lik, forres$model_data))