This set of functions allows you to do dlnm parameter sweep. Testing many different iterations of DLNM options quickly and efficiently. Github page is here

Current method of doing DLNM

From Gasparrini 2014:

library(dlnm) ; library(splines) ; library(foreign) ; library(tsModel)

# LOAD THE DATA
lndn <- read.csv(paste0("https://raw.githubusercontent.com/gasparrini/",
"2014_gasparrini_BMCmrm_Rcodedata/e629be7809853f390d305d8e73ac4b152a54f092/london.csv"))

# KNOTS FOR EXPOSURE-RESPONSE FUNCTION
vk <- equalknots(lndn$tmean, fun="bs", degree = 2, df = 4)

# KNOTS FOR THE LAG-RESPONSE FUNCTION
maxlag <- 25
lk <- logknots(x = maxlag, nk = 3)

# COMPUTE THE CROSS-BASIS
cb <- crossbasis(x = lndn$tmean, 
                 lag=maxlag, 
                 argvar=list(fun="bs", degree = 2, knots = vk), 
                 arglag=list(knots=lk))

# RUN THE MODEL
model <- glm(death~cb+ns(time,10*14)+dow,family=quasipoisson(),lndn)

And create outputs

# PREDICTION AT 99TH PERCENTILES
# COMPUTE 99TH PERCENTILES
perc <- quantile(lndn$tmean,c(0.99))
# CENTERING VALUE (AND PERCENTILE)
cen <- 20
sum(lndn$tmean<cen,na.rm=T)/sum(!is.na(lndn$tmean))

## [1] 0.9360454

pred <- crosspred(cb,model,at=perc,bylag=0.2,cen=cen)

# This is already exponentiated
o <- data.frame(est = pred$allRRfit, lb = pred$allRRlow, ub = pred$allRRhigh)
row.names(o) <- 1
o

##       est       lb       ub
## 1 1.25294 1.188849 1.320487

But what if you wanted to test other spline options? Well now you can:

DLNM Sweep

DLNM Sweep uses a custom function called make_cb_list to expand crossbasis options:

cb2 <- make_cb_list(x = lndn,
                    var_list = c('tmean'), 
                    argvar_list = list(set1 = list(fun="bs", degree = 2, knots = list(vk))) ,
                    lag_list = maxlag,
                    arglag_list = list(set1 = list(knots = list(lk))))

You can see that this creates the same cb object as in the basic DLNM (with some novel attributes removed):

attr(cb2[[1]], 'summary_str') <- NULL
attr(cb2[[1]], 'var_name') <- NULL
identical(cb, cb2[[1]])

## [1] TRUE

Ok so now how to use this to its fullest potential. One application might look like this:

var_list <- c('tmean', 'tmin')

cb_all <- list()

for(var in var_list) {
  
  # make var knots that vary by var
  vk2 = equalknots(lndn[, var], fun="bs", degree = 2, df = 4)
  vk3 = equalknots(lndn[, var], fun="ns", df = 4)
  k50 = quantile(lndn[, var], probs = .5, na.rm = T)
  k10 = quantile(lndn[, var], probs = c(1:9)/10, na.rm = T)
  
  cb_var = make_cb_list(x = lndn,
                         var_list = var, 
                         argvar_list = list(set1 = list(fun="bs", degree = c(2, 3), 
                                                        knots = list(vk2, k50, k10)),
                                            set2 = list(fun="ns",  
                                                        knots = list(vk3, k50, k10))) ,
                         lag_list = maxlag,
                         arglag_list = list(set1 = list(knots = list(lk))))
  
  cb_all <- append(cb_all, cb_var)
}

So now you’ve created the crossbasis matrices, how many versions are there? Should be: 2 variables x (2 bs degrees x 3 bs argvar knots + 1 ns degrees x 3 ns argvar knots)

length(cb_all)

## [1] 18

Next step is to create the model objects (implemented below in parallel):

# setup for parallel processing
library(future)          # parallel processing
library(future.apply)    # parallel processing
library(progressr)       # progress updates during parallel
plan(multisession)
handlers("progress")

get_model_obj <- function(cb_list) {
  cb_seq <- seq(length(cb_list))
  p      <- progressr::progressor(along = cb_seq)
  future_lapply(cb_seq, function(i) {
    cb <- cb_list[[i]]
    p(message = sprintf("%s", attr(cb, "summary_str")))
    # **********
    # model object, same as above
    m  <- glm(death ~ cb + ns(time,10*14) + dow, family = quasipoisson(), lndn)
    # **********
    attr(m, "var_name") <- attr(cb, "var_name")
    attr(m, "summary_str") <- attr(cb, "summary_str")
    m
  })
}

# get objects
with_progress(models <- get_model_obj(cb_all))

Clearly this could be expanded to a two-cross basis example by creating an expanded grid of the length of each cb_list.

Final step is get model predictions, implemented again in parallel:

get_pred <- function(cb_list, model_list) {
  stopifnot(length(cb_list) == length(model_list))
  cb_seq <- seq(length(cb_list))
  p      <- progressr::progressor(along = cb_seq)
  out <- future_lapply(cb_seq, function(i) {
    cb <- cb_list[[i]]
    m  <- model_list[[i]]
    p(message = sprintf("%s", attr(cb, "summary_str")))
    # ********
    # get prediction and centering values
    perc <- quantile(lndn[, attr(cb, "var_name")],c(0.99), na.rm = T)
    cen  <- quantile(lndn[, attr(cb, "var_name")],c(0.94), na.rm = T) # ~ 20c for Tmax
    pred <- crosspred(cb, m, at=perc, bylag=0.2, cen=cen)
    # ********
    data.frame(var_name = attr(cb, "var_name"), 
               summary_str = attr(cb, "summary_str"),
               est = pred$allRRfit[1], 
               lb = pred$allRRlow[1], 
               ub = pred$allRRhigh[1])
  })
  o <- data.frame(do.call(rbind, out))
  row.names(o) <- 1:nrow(o)
  return(o)
}
# get betas
with_progress(RRs <- get_pred(cb_all, models))

And now plot

RRs %>%
  ggplot(., aes(y = reorder(summary_str, est))) +
  geom_point(aes(x=est), shape=15, size=3) +
  geom_linerange(aes(xmin=lb, xmax=ub)) + ylab(NULL)

DLNM parameter sweep

Chad Milando

2023-07-27

Current method of doing DLNM

DLNM Sweep