The flexsurv package in R facilitates the use of background mortality estimates to model excess hazard. In the examples I can find (the standsurv vignette and the 2023 paper by Sweeting et al.), the source of this background mortality is the survival::survexp.us object. I would like to use more relevant sources of background mortality but cannot figure out how to do this.

I tried using the HMDHFDplus to obtain mortality rates (recommended by flexsurv). This produces a data.frame which looks helpful, though it is not a ratetable object or an array like survexp.us. Below, I use the survival::lung dataset to fit a basic Weibull model and attempt to incorporate population mortality.

df <- lung |> 
  select(time, status, age, sex) |>
  mutate(
    sex = if_else(sex == 1, "Male", "Female"), # same as lifetable below
    status = status - 1, # coded as 1 and 2 otherwise
    year = 2020 # just as an example
  ) |> 
  as_tibble()

lifetable <- HMDHFDplus::readHMDweb(
  "AUS", # Australia
  item = "Mx_1x1", # Male and Female annual mortality rates
  username = _username_, 
  password = _password_ 
) |> 
  filter(Year == 2020) |> 
  pivot_longer(names_to = "sex", values_to = "rate", cols = contains("male")) |> 
  mutate(rate = rate / 365.25) |> # from annual to daily rate
  select(year = Year, age = Age, sex, rate)

df <- df |> 
  left_join(lifetable, by = c("age", "sex", "year"))

weibull <- flexsurvreg(
  Surv(time, status) ~ sex + age, 
  data = df, 
  dist = "weibullPH",
  bhazard = rate
)

standsurv(weibull,
          t = max(df$time),
          at = list(list(sex = "Male"),
                    list(sex = "Female")),
          rmap = list(age = age,
                      year = year,
                      sex = sex),
          ratetable = ratetable(lifetable),
          newdata = df)

Marginal all-cause survival will be calculated
Calculating marginal expected survival and hazard
Error in survexp(formula = t.temp ~ 1, rmap = list(age = age, year = year, :
Invalid rate table