how to get relative risk from a GAM with a distributed lag model

Question

I would like to obtain a relative risk/ risk ratio from a GAM with a distributed lag model. I have a GAM (implemented in mgcv in R) predicting daily deaths from time series data consisting of daily temperature, humidity and rainfall. The GAM includes a distributed lag model because deaths may occur over several days following a high heat day.

What I'd like to do is compute and plot the relative risk (accumulated across all lags) against temperature for a given reference temperature, e.g. the temperature at which the risk is lowest, with corresponding confidence intervals.(Additionally, I'd also like to plot the relative risk against temperature for different lags separately). I'm not sure how to go about this and would really appreciate some pointers. (I am aware of the predict.gam function but am not sure if and how it should be used in this case).

I found a post by Simon Wood on a related issue and also a paper (see supplement/ additional file for code) which extracts an RR from a GAM, but I understand not all of it, and I'm also not sure how I would apply this to a GAM with time series data and distributed lag.

Below is a reproducible example of my GAM:

library(mgcv)
set.seed(3) # make reproducible example
simdat <- gamSim(1,400) # simulate data
g <- exp(simdat$f/5)
simdat$y <- rnbinom(g,size=3,mu=g)  # negative binomial response var
simdat$time <- 1:400  # create time series
names(simdat) <- c("deaths", "temp", "humidity", "rain", "x3", "f",
"f0", "f1", "f2", "f3", "time")

# lag function based on Simon Wood (book 2017, p.349 and gamair
# package documentation p.54
# https://cran.rstudio.com/web/packages/gamair/gamair.pdf)
lagard <- function(x,n.lag=7) {
n <- length(x); X <- matrix(NA,n,n.lag)
for (i in 1:n.lag) X[i:n,i] <- x[i:n-i+1]
X
}

# set up lag, temp, rain and humidity as 7-column matrices
# to create lagged variables - based on Simon Wood's example
dat <- list(lag=matrix(0:6,nrow(simdat),7,byrow=TRUE),
deaths=simdat$deaths, time = simdat$time)
dat$temp <- lagard(simdat$temp)
dat$rain <- lagard(simdat$rain)
dat$humidity <- lagard(simdat$humidity)

mod <- gam(deaths~s(time, k=70) + te(temp, lag, k=c(12, 4)) +
te(humidity, lag, k=c(12, 4)) + te(rain, lag, k=c(12, 4)), data = dat,
family = nb, method = 'REML', select = TRUE)

summary(mod)
plot(mod, scheme = 1)
plot(mod, scheme = 2)

I'd be very grateful for any suggestions you may have!

UPDATE: There is a way to get overall relative risk, by using predict.gam() with the type = terms argument, and then for each prediction temperature summing over the lag periods. However, this does not (I think?) allow calculation of confidence intervals around the RRs. Via the r-help mailing list Simon Wood pointed out that using type = lpmatrix allows to get CIs. My problem is that I don't know how to work this for a model with distributed lag.