Question: Interpretation when marginal all-cause survival appears to deviate from observed survival

[connor-ballinger]

Hi,

Thanks for the package, the documentation has been a really useful learning resource.

I am using flexsurv to extrapolate survival using all-cause mortality data in an excess hazard / relative survival framework. When plotting marginal survival from the relative survival model, I can see that it does not fit the observed data very well. In contrast, the standard flexsurv models which do not include excess mortality fit the observed data much better (and yet the literature strongly recommends incorporating population mortality).

Why do these models appear to systematically underestimate observed survival when matched population mortality is incorporated?

It appears this is more likely to happen when there is little difference between survival in the sample and the population (excess hazard is very low). My dataset cannot be shared but I have attached two images: both are Weibull curves and only the second includes background mortality in the estimation. Clearly, the fitted curves have shifted downwards in the relative survival model.

I have tried to provide an illustrative example using the breast cancer dataset from flexsurv. I have crudely adjusted the population hazard by doubling it, such that there is little difference between survival in the population and the sample. This results in the estimated model shifting to a position significantly below the Kaplan-Meier.

Any guidance would be greatly appreciated.

Thanks.

library(tidyverse)
library(flexsurv)
library(ggsurvfit)
set.seed(1)

## Data ------------------------------------------------------------------------
data(bc, package = "flexsurv")

# mostly code from flexsurv vignette
bc$age <- rnorm(dim(bc)[1], mean = 65 - scale(bc$recyrs, scale=F), sd = 5)
## Create age at diagnosis in days - used later for matching to expected rates
bc$agedays <- floor(bc$age * 365.25)
## Create some random diagnosis dates between 01/01/1984 and 31/12/1989
bc$diag <- as.Date(floor(runif(dim(bc)[1], as.Date("01/01/1984", "%d/%m/%Y"), 
                               as.Date("31/12/1989", "%d/%m/%Y"))), 
                   origin="1970-01-01")
## Create sex (assume all are female)
bc$sex <- factor("female")
## 2-level prognostic variable
bc$group2 <- ifelse(bc$group=="Good", "Good", "Medium/Poor")
head(bc)

# Kaplan-Meier
km <- survfit(Surv(recyrs, censrec) ~ group2, data=bc)

## Excess mortality ------------------------------------------------------------
# reshape US lifetable to be a tidy data.frame, and convert rates to per 
# person-year as flexsurv regression is in years
survexp.us.df <- as.data.frame.table(survexp.us, responseName = "exprate") %>%
  mutate(exprate = exprate * 365.25)
survexp.us.df$age <- as.numeric(as.character(survexp.us.df$age))
survexp.us.df$year <- as.numeric(as.character(survexp.us.df$year))

## Obtain attained age and attained calendar year in (whole) years
bc <- bc |> 
  mutate(attained.age.yr = floor(age + recyrs),
         attained.year = lubridate::year(diag + rectime))

## merge in (left join) expected rates at event time
bc <- bc |> 
  left_join(
    survexp.us.df,
    by = c("attained.age.yr"="age", 
           "attained.year"="year", 
           "sex"="sex")
  ) 

## A standard model ------------------------------------------------------------
weib_fit1 <- flexsurvreg(
  Surv(recyrs, censrec) ~ group2, 
  anc = list(shape = ~ group2), 
  data = bc, 
  dist = "weibullPH"
)

### Standardise ----------------------------------------------------------------
weib_stand1 <- standsurv(
  weib_fit1,
  type = "survival",
  at = list(
    list(group2 = "Good"),
    list(group2 = "Medium/Poor")
  ),
  t = seq(0, 25, length = 50),
  rmap = list(
    sex = sex,
    year = diag,
    age = agedays
  ),
  ratetable = survexp.us,
  scale.ratetable = 365.25,
  newdata = bc
)
### Plot -----------------------------------------------------------------------
plot1 <- plot(weib_stand1, expected = TRUE) + 
  geom_step(
    data = ggsurvfit::tidy_survfit(km),
    aes(x = time, y = estimate, colour = strata)
  )

## A relative survival model ---------------------------------------------------
weib_fit2 <- flexsurvreg(
  Surv(recyrs, censrec) ~ group2, 
  anc = list(shape = ~ group2), 
  data = bc, 
  dist = "weibullPH",
  bhazard = exprate # incorporate population hazard
)

### Standardise ----------------------------------------------------------------
weib_stand2 <- standsurv(
  weib_fit2,
  type = "survival",
  at = list(
    list(group2 = "Good"),
    list(group2 = "Medium/Poor")
  ),
  t = seq(0, 25, length = 50),
  rmap = list(
    sex = sex,
    year = diag,
    age = agedays
  ),
  ratetable = survexp.us,
  scale.ratetable = 365.25,
  newdata = bc
)
### Plot -----------------------------------------------------------------------
plot2 <- plot(weib_stand2, expected = TRUE) + 
  geom_step(
    data = ggsurvfit::tidy_survfit(km),
    aes(x = time, y = estimate, colour = strata)
  )

## Increase background mortality -----------------------------------------------
weib_fit3 <- flexsurvreg(
  Surv(recyrs, censrec) ~ group2, 
  anc = list(shape = ~ group2), 
  data = bc, 
  dist = "weibullPH",
  bhazard = exprate * 2 # Double the hazard in the matched population
)

weib_stand3 <- standsurv(
  weib_fit3,
  type = "survival",
  at = list(
    list(group2 = "Good"),
    list(group2 = "Medium/Poor")
  ),
  t = seq(0, 25, length = 50),
  rmap = list(
    sex = sex,
    year = diag,
    age = agedays
  ),
  ratetable = survexp.us,
  scale.ratetable = 365.25 * 2, # Double the hazard in the matched population
  newdata = bc
)

plot3 <- plot(weib_stand3, expected = TRUE) + 
  geom_step(
    data = ggsurvfit::tidy_survfit(km),
    aes(x = time, y = estimate, colour = strata)
  )

[chjackson]

I'm having trouble finding time to look at this in detail, but I can see there's something worth investigating.

@mikesweeting I wonder if you had time - given that this appears to relate to how standsurv constructs the standardised survival curves from a relative survival model, which I'm not very familiar with.

At the same time @connor-ballinger if there is any way the example could be made simpler to demonstrate the problem more directly, that would be appreciated. e.g. no covariates, no "tweaks", so we have the simplest possible relative survival model in which the study population is identical to the background population.

[mikesweeting]

Hi @connor-ballinger , I will try to have a look in more detail at this later. But one thing I just noticed in your reproducible example is that when you doubled the expected rates used to fit the relative survival model, you didn't also double the rates in the lifetable used in standsurv to calculate marginal all-cause survival. I'm not sure if this is causing the issue you see in your actual example though?

[connor-ballinger]

Hi, thank you both for your time. I feel I should disclose that I do not have much experience with survival analysis, so there is certainly a chance that this issue is due to me not comprehending some detail or concept, either in the package or more generally.

For context, I am hoping to publish this analysis and intended to use the background mortality models, as they are encouraged in the literature. However, it is difficult to justify choosing these models over regular models when the visual fits are comparatively so poor.

I have simplified the example code below so that there are no treatment arms or shape coefficients, as @chjackson suggested. I have provided some plots to reflect these changes to my actual model, in case that is helpful.

Regarding the doubling of the population mortality rates @mikesweeting, I have doubled the bhazard in flexsurvreg and the scale.ratetable in standsurv. (Ctrl+F: "# Double the hazard in the matched population" - strange that GitHub does not display line numbers...). I think multiplying the scale.ratetable argument works the same as multiplying the mortality array itself, but perhaps not. I am struggling to check it now, I find the ratetables difficult to work with.

Below are the plots for the breast cancer example. First is the regular survival model, followed by the relative survival figure, followed by the figure with population rates doubled. The marginal curve shifts left in the last case, but the fit against the Kaplan-Meier is perhaps not noticeably worse than the two previous plots. Compare that to my relative survival model (further below) where the estimated survival is entirely below the Kaplan-Meier.

Note in the plot below, marginal survival has shifted left/down compared to the above model, but the fit isn't visibly worse.

(Note I ran this code under the first section I posted, which contained some data wrangling)

## Version 2 -------------------------------------------------------------------
km <- survfit(Surv(recyrs, censrec) ~ 1, data=bc)

### A standard model ----------------------------------------------------------
weib_fit1 <- flexsurvreg(
  Surv(recyrs, censrec) ~ 1, 
  data = bc, 
  dist = "weibullPH"
)

### Standardise ----------------------------------------------------------------
weib_stand1 <- standsurv(
  weib_fit1,
  type = "survival",
  t = seq(0, 25, length = 50),
  rmap = list(
    sex = sex,
    year = diag,
    age = agedays
  ),
  ratetable = survexp.us,
  scale.ratetable = 365.25,
  newdata = bc
)
### Plot -----------------------------------------------------------------------
plot1 <- plot(weib_stand1, expected = TRUE) + 
  geom_step(
    data = ggsurvfit::tidy_survfit(km),
    aes(x = time, y = estimate)
  )

### A relative survival model ---------------------------------------------------
weib_fit2 <- flexsurvreg(
  Surv(recyrs, censrec) ~ 1, 
  data = bc, 
  dist = "weibullPH",
  bhazard = exprate # incorporate population hazard
)

### Standardise ----------------------------------------------------------------
weib_stand2 <- standsurv(
  weib_fit2,
  type = "survival",
  t = seq(0, 25, length = 50),
  rmap = list(
    sex = sex,
    year = diag,
    age = agedays
  ),
  ratetable = survexp.us,
  scale.ratetable = 365.25,
  newdata = bc
)
### Plot -----------------------------------------------------------------------
plot2 <- plot(weib_stand2, expected = TRUE) + 
  geom_step(
    data = ggsurvfit::tidy_survfit(km),
    aes(x = time, y = estimate)
  )

### Increase background mortality -----------------------------------------------
weib_fit3 <- flexsurvreg(
  Surv(recyrs, censrec) ~ 1, 
  data = bc, 
  dist = "weibullPH",
  bhazard = exprate * 2 # Double the hazard in the matched population
)

weib_stand3 <- standsurv(
  weib_fit3,
  type = "survival",
  t = seq(0, 25, length = 50),
  rmap = list(
    sex = sex,
    year = diag,
    age = agedays
  ),
  ratetable = survexp.us,
  scale.ratetable = 365.25 * 2, # Double the hazard in the matched population
  newdata = bc
)

plot3 <- plot(weib_stand3, expected = TRUE) + 
  geom_step(
    data = ggsurvfit::tidy_survfit(km),
    aes(x = time, y = estimate)
  )

My own model, simplified in the same ways (apologies if this is overkill, I figure you can look at as much or as little as you like): First plot, without background mortality: Hazard associated with the above model: Relative survival model, where the fit appears worse than the standard model above: Hazard associated with the relative survival model above.

I completely understand if you don't have time to investigate. I'm sure you get all sorts of anomalies and issues brought to you when maintaining a package.

Thanks.

[irtimmins]

Hi Connor, running your second example with the following code:

I can't replicate the figure with very poor fit. (sorry that may have been a plot corresponding to a different dataset of yours?) So perhaps I have replicated everything above fine. Each plot below looks okay fit.

Might be worth re-installing and updating flexsurv package to version 2.3, survival package to 3.7 and rerunning?

Running the following data wrangling steps as above:

library(tidyverse)
library(flexsurv)
library(ggsurvfit)
set.seed(1)

## Data ------------------------------------------------------------------------
data(bc, package = "flexsurv")

# mostly code from flexsurv vignette
bc$age <- rnorm(dim(bc)[1], mean = 65 - scale(bc$recyrs, scale=F), sd = 5)
## Create age at diagnosis in days - used later for matching to expected rates
bc$agedays <- floor(bc$age * 365.25)
## Create some random diagnosis dates between 01/01/1984 and 31/12/1989
bc$diag <- as.Date(floor(runif(dim(bc)[1], as.Date("01/01/1984", "%d/%m/%Y"), 
                               as.Date("31/12/1989", "%d/%m/%Y"))), 
                   origin="1970-01-01")
## Create sex (assume all are female)
bc$sex <- factor("female")
## 2-level prognostic variable
bc$group2 <- ifelse(bc$group=="Good", "Good", "Medium/Poor")
head(bc)

## Excess mortality ------------------------------------------------------------
# reshape US lifetable to be a tidy data.frame, and convert rates to per 
# person-year as flexsurv regression is in years
survexp.us.df <- as.data.frame.table(survexp.us, responseName = "exprate") %>%
  mutate(exprate = exprate * 365.25)
survexp.us.df$age <- as.numeric(as.character(survexp.us.df$age))
survexp.us.df$year <- as.numeric(as.character(survexp.us.df$year))

## Obtain attained age and attained calendar year in (whole) years
bc <- bc |> 
  mutate(attained.age.yr = floor(age + recyrs),
         attained.year = lubridate::year(diag + rectime))

## merge in (left join) expected rates at event time
bc <- bc |> 
  left_join(
    survexp.us.df,
    by = c("attained.age.yr"="age", 
           "attained.year"="year", 
           "sex"="sex")
  ) 

First model without background mortality:

##############
km <- survfit(Surv(recyrs, censrec) ~ 1, data=bc)

### A standard model ----------------------------------------------------------
weib_fit1 <- flexsurvreg(
  Surv(recyrs, censrec) ~ 1, 
  data = bc, 
  dist = "weibullPH"
)

### Standardise ----------------------------------------------------------------
weib_stand1 <- standsurv(
  weib_fit1,
  type = "survival",
  t = seq(0, 25, length = 50),
  rmap = list(
    sex = sex,
    year = diag,
    age = agedays
  ),
  ratetable = survexp.us,
  scale.ratetable = 365.25,
  newdata = bc
)

### Plot -----------------------------------------------------------------------
plot1 <- plot(weib_stand1, expected = TRUE) + 
  geom_step(
    data = ggsurvfit::tidy_survfit(km),
    aes(x = time, y = estimate)
  )

Here's plot1:

Second model with relative survival:

### A relative survival model ---------------------------------------------------
weib_fit2 <- flexsurvreg(
  Surv(recyrs, censrec) ~ 1, 
  data = bc, 
  dist = "weibullPH",
  bhazard = exprate # incorporate population hazard
)

### Standardise ----------------------------------------------------------------
weib_stand2 <- standsurv(
  weib_fit2,
  type = "survival",
  t = seq(0, 25, length = 50),
  rmap = list(
    sex = sex,
    year = diag,
    age = agedays
  ),
  ratetable = survexp.us,
  scale.ratetable = 365.25,
  newdata = bc
)
### Plot -----------------------------------------------------------------------
plot2 <- plot(weib_stand2, expected = TRUE) + 
  geom_step(
    data = ggsurvfit::tidy_survfit(km),
    aes(x = time, y = estimate)
  )

Here's plot2:

Third model with increased background mortality:

### Increase background mortality -----------------------------------------------
weib_fit3 <- flexsurvreg(
  Surv(recyrs, censrec) ~ 1, 
  data = bc, 
  dist = "weibullPH",
  bhazard = exprate * 2 # Double the hazard in the matched population
)

weib_stand3 <- standsurv(
  weib_fit3,
  type = "survival",
  t = seq(0, 25, length = 50),
  rmap = list(
    sex = sex,
    year = diag,
    age = agedays
  ),
  ratetable = survexp.us,
  scale.ratetable = 365.25 * 2, # Double the hazard in the matched population
  newdata = bc
)

plot3 <- plot(weib_stand3, expected = TRUE) + 
  geom_step(
    data = ggsurvfit::tidy_survfit(km),
    aes(x = time, y = estimate)
  )

Here's plot3: The background mortality rate is higher as expected, and the fit still seems ok.

[irtimmins]

Sorry may have misread your post - I understand that the figure above with the deviated observed and all-cause was from another dataset that couldn't be shared. So I think the figures I showed are the same as you had above.

Sometimes these issues happen when the background/expected rates are on different scales (years vs months vs days), could be worth checking that?

The scale in the above example from bc datasets is all in years. Does your model code have the Surv(years, event) survival model with event time measured in years and not months/days?

[connor-ballinger]

Hi @irtimmins, thanks for your time.

Yes, the figures which show extrapolation to just 15 years are related to my own project, whereas the breast cancer example has been extrapolated to 25 years (that may be the easiest distinction to make given I failed to label them).

Yep, it looks to me like you have replicated my three plots accurately. Perhaps the example is too crude, but I was trying to show how low excess hazard may result in the estimated marginal survival being pushed further and further away from the Kaplan-Meier (which is what I assume is happening with my model, but that idea could be wildly incorrect).

I believe everything in my dataset is measured in years, certainly follow-up time is in years. I sourced background mortality from the Human Mortality Database, using the HMDHFDplus package. I replicated the method described above where US rates were joined to the bc dataset. I supplied a value of 1 to scale.ratetable (compared to 365.25 to the bc dataset above) as I believe that annual mortality rates are provided (so no correction needed). My method of testing this was just to use different arguments for scale.ratetable (divided or multiplied by 365.25). Neither of those scenarios produce anything sensible.

I have been using survival version 3.7-0 from CRAN and flexsurv 2.3 from GitHub. I have now reinstalled them, as suggested.

I can send through my code for extracting background mortality if you think it helpful. I am a bit concerned I may just be sending you down some very time-consuming rabbit holes.

Thanks.

[irtimmins]

Hi @connor-ballinger, thanks for that info. I've redone the bc dataset example using US rates from HMD. I get quite similar results to above. Perhaps worth seeing first how this compares to your code?

I would also recommend the tutorial paper by Mike Sweeting, which has a nice R code example in the supplementary material (https://pubmed.ncbi.nlm.nih.gov/37448102/) for doing relative survival with HMD rates.

Preparing the HMD rates for US:

library(survival)
library(flexsurv)
require(HMDHFDplus)
require(relsurv) 
require(dplyr)

## Download 1x1 Period lifetables from HMD (fltper_1x1 for females, mltper_1x1 for males)
myHMDusername <- "..." # fill in login details
myHMDpassword <- "..." # ....

## create a popmort_files subdirectory if it does not exist
if(!dir.exists("popmort_files"))  dir.create("popmort_files")

## Countries that we will get ratetables 
country <- c("USA")  ## US rates
for(ctry in country){
  print(ctry)
  # Read in popmort file from HMD if the file does not exist
  if(!file.exists(paste0("popmort_files/",ctry,".f.txt"))){
    # Females
    HMD.f <- readHMDweb(CNTRY = ctry, item = "fltper_1x1", username = myHMDusername,
                        password = myHMDpassword)
    write.table(HMD.f, file=paste0("popmort_files/",ctry,".f.txt"))
    # Males
    HMD.m <- readHMDweb(CNTRY = ctry, item = "mltper_1x1", username = myHMDusername,
                        password = myHMDpassword)
    write.table(HMD.m, file=paste0("popmort_files/",ctry,".m.txt"))
  }
  ratetable <- relsurv::transrate.hmd(male = paste0("popmort_files/",ctry,".m.txt"),
                                      female = paste0("popmort_files/",ctry,".f.txt"))

  ## reshape lifetable to be a tidy data.frame, and convert rates to per person-year 
  # (needed if time scale of survival model is in years)

  ratetable.df <- as.data.frame.table(ratetable, responseName = "exprate") %>%
    mutate(exprate = 365.25 * exprate) %>%
    mutate(age = as.numeric(as.character(age)),
           year = as.numeric(as.character(year)))

  assign(paste0("ratetable.",ctry), ratetable)
  assign(paste0("ratetable.",ctry,".df"), ratetable.df)
}

Then for bc example:

## Prep data
set.seed(1)
data(bc, package = "flexsurv")

#
bc$age <- rnorm(dim(bc)[1], mean = 65 - scale(bc$recyrs, scale=F), sd = 5)
## Create age at diagnosis in days - used later for matching to expected rates
bc$agedays <- floor(bc$age * 365.25)
## Create some random diagnosis dates between 01/01/1984 and 31/12/1989
bc$diag <- as.Date(floor(runif(dim(bc)[1], as.Date("01/01/1984", "%d/%m/%Y"), 
                               as.Date("31/12/1989", "%d/%m/%Y"))), 
                   origin="1970-01-01")

## Create sex (assume all are female)
bc$sex <- factor("female")
## 2-level prognostic variable
bc$group2 <- ifelse(bc$group=="Good", "Good", "Medium/Poor")

## Obtain attained age and attained calendar year in (whole) years
bc <- bc |> 
  mutate(attained.age.yr = floor(age + recyrs),
         attained.year = lubridate::year(diag + rectime))

bc <- bc |> 
  left_join(
    ratetable.USA.df,
    by = c("attained.age.yr"="age", 
           "attained.year"="year", 
           "sex"="sex")
  ) 

km <- survfit(Surv(recyrs, censrec) ~ 1, data=bc)

# fit relative survival model
weib_fit2_hmd <- flexsurvreg(
  Surv(recyrs, censrec) ~ 1, 
  data = bc, 
  dist = "weibullPH",
  bhazard = exprate # incorporate population hazard
)

### marginal predictions using HMD US rates.

weib_stand2_hmd  <- standsurv(
  weib_fit2_hmd,
  type = "survival",
  t = seq(0, 25, length = 50),
  rmap = list(
    sex = sex,
    year = diag,
    age = agedays
  ),
  ratetable = ratetable.USA,
  scale.ratetable = 365.25,
  newdata = bc
)

### Plot 
plot2_hmd <- plot(weib_stand2_hmd, expected = TRUE) + 
  geom_step(
    data = ggsurvfit::tidy_survfit(km),
    aes(x = time, y = estimate)
  )

Here's plot2_hmd

[irtimmins]

Hi @connor-ballinger, going back the earlier point about how you double the background rates.

The scale.ratetable argument is telling standsurv how the units of time in the survival model (the "time" in Surv(time,event), which in the bc example is in years) compares with the units of time in the ratetable object (which is rates per day). So we should keep that as 365.25 as the scale.

If we want to double the background rates, we can do this by doubling the rates in the ratetable object:

# create new ratetable with doubled rates.
survexp.us.high.rate <- survexp.us
survexp.us.high.rate[,,] <- 2*survexp.us.high.rate[,,] # doubles all the rates

weib_fit3 <- flexsurvreg(
  Surv(recyrs, censrec) ~ 1, 
  data = bc, 
  dist = "weibullPH",
  bhazard = 2*exprate  # Double the hazard in the matched population
)

weib_stand3 <- standsurv(
  weib_fit3,
  type = "survival",
  t = seq(0, 25, length = 50),
  rmap = list(
    sex = sex,
    year = diag,
    age = agedays
  ),
  ratetable = survexp.us.high.rate, # use the new ratetable
  scale.ratetable = 365.25 , # keep this the same as before
  newdata = bc
)

plot3 <- plot(weib_stand3, expected = TRUE) + 
  geom_step(
    data = ggsurvfit::tidy_survfit(km),
    aes(x = time, y = estimate)
  )

plot3

This gets a bit better visual fit than before. I think in your earlier example, the rates in the bhazard argument no longer matched the rates in the ratetable, and I think this may have caused issues with the model fit.

Thanks, Iain

[connor-ballinger]

Yes, I did my best to follow along with the Sweeting 2023 paper. My other key resource was the standsurv vignette. Admittedly, the code you used to extract mortality rates from the HMD is a bit different to mine - perhaps this could be an issue. I could not find how to construct a ratetable so I may have made an error there. I will try out relsurv::transrate.hmd tomorrow (~12 hours).

Thank you for the direction on the doubling of the ratetable, I couldn't figure out how to do it properly. The code for the bc example is certainly more intuitive now. I agree that the model appears to fit the KM curve better than in the previous iterations. However, is this possibly because background mortality appears to be occurring at a slower rate (or, more simply, the dotted blue line is higher in this last plot compared to previous plots, meaning that marginal survival is not being depressed to the same extent, so marginal survival stays slightly higher and fits the KM reasonably well). I'm sorry, that may sound like egregious nitpicking.

The code below has been adapted for the bc example. It shows my method for extracting background mortality. There is perhaps a better way to construct the ratetable object but I haven't been able to find any examples.

(Skip to the lifetable section)

# bc w/ aus mortality ----------------------------------------------------------

library(tidyverse)
library(flexsurv)
library(ggsurvfit)
set.seed(1)

## Data ------------------------------------------------------------------------
data(bc, package = "flexsurv")

# mostly code from flexsurv vignette
bc$age <- rnorm(dim(bc)[1], mean = 65 - scale(bc$recyrs, scale=F), sd = 5)
## Create age at diagnosis in days - used later for matching to expected rates
bc$agedays <- floor(bc$age * 365.25)
## Create some random diagnosis dates between 01/01/1984 and 31/12/1989
bc$diag <- as.Date(floor(runif(dim(bc)[1], as.Date("01/01/1984", "%d/%m/%Y"), 
                               as.Date("31/12/1989", "%d/%m/%Y"))), 
                   origin="1970-01-01")
## Create sex (assume all are female)
bc$sex <- factor("Female")
bc <- tibble(bc)
# Kaplan-Meier
km <- survfit(Surv(recyrs, censrec) ~ 1, data=bc)

## Excess mortality ------------------------------------------------------------
# download australian mortality from hmd

### Download Lifetable Data ----------------------------------------------------

lifetable <- HMDHFDplus::readHMDweb(
  "AUS",
  item = "Mx_1x1", ### Different to your method ###
  username = ___,
  password = ___
)

df_lifetable_long <- lifetable |> 
  select(-c(Total, OpenInterval)) |> 
  pivot_longer(
    cols = matches("ale$"),
    names_to = "sex", 
    values_to = "exprate"
  ) |> 
  rename_with(tolower)

### Join excess mortality ------------------------------------------------------
## Obtain attained age and attained calendar year in (whole) years
bc <- bc |> 
  mutate(attained.age.yr = floor(age + recyrs),
         attained.year = lubridate::year(diag + rectime))

## merge in (left join) expected rates at event time
bc <- bc |> 
  left_join(
    df_lifetable_long,
    by = c("attained.age.yr"="age", 
           "attained.year"="year", 
           "sex"="sex")
  ) 

### Ratetable ------------------------------------------------------------------
# Need to produce same structure as survexp.us. ### This may be dodgy, I wasn't aware of relsurv::transrate.hmd ###

a_mortality <- array(
  data = NA,
  dim = c(111, 2, 100), # lengths of dimnames
  dimnames = list(
    "age" = unique(df_lifetable_long$age),
    "sex" = unique(df_lifetable_long$sex),
    "year" = unique(df_lifetable_long$year)
  )
)
a_mortality[, "Male", ] <- df_lifetable_long |>
  filter(sex == "Male") |>
  pull(exprate)
a_mortality[, "Female", ] <- df_lifetable_long |>
  filter(sex == "Female") |>
  pull(exprate)

attr(a_mortality, "type") <- c(2, 1, 3)
attr(a_mortality, "cutpoints") <- list(
  as.numeric(0:110),
  NULL,
  seq.Date(as.Date("1921-01-01"), as.Date("2020-01-01"), by = "year")
)
class(a_mortality) <- "ratetable"
is.ratetable(a_mortality, verbose = TRUE)

## A standard model ------------------------------------------------------------
weib_fit1 <- flexsurvreg(
  Surv(recyrs, censrec) ~ 1, 
  data = bc, 
  dist = "weibullPH"
)

## Standardise -----------------------------------------------------------------
weib_stand1 <- standsurv(
  weib_fit1,
  type = "survival",
  t = seq(0, 25, length = 50),
  rmap = list(
    sex = sex,
    year = diag,
    age = agedays
  ),
  ratetable = a_mortality,
  scale.ratetable = 365.25,
  newdata = bc
)
## Plot ------------------------------------------------------------------------
plot1 <- plot(weib_stand1, expected = TRUE) + 
  geom_step(
    data = ggsurvfit::tidy_survfit(km),
    aes(x = time, y = estimate)
  )

## A relative survival model ---------------------------------------------------
weib_fit2 <- flexsurvreg(
  Surv(recyrs, censrec) ~ 1, 
  data = bc, 
  dist = "weibullPH",
  bhazard = exprate # incorporate population hazard
)

## Standardise -----------------------------------------------------------------
weib_stand2 <- standsurv(
  weib_fit2,
  type = "survival",
  t = seq(0, 25, length = 50),
  rmap = list(
    sex = sex,
    year = diag,
    age = agedays
  ),
  ratetable = a_mortality,
  scale.ratetable = 365.25,
  newdata = bc
)
## Plot ------------------------------------------------------------------------
plot2 <- plot(weib_stand2, expected = TRUE) + 
  geom_step(
    data = ggsurvfit::tidy_survfit(km),
    aes(x = time, y = estimate)
  )

## Increase background mortality -----------------------------------------------
a_mortality_v2 <- a_mortality * 2 
# a_mortality_v2 <- a_mortality[, , ] * 2 
### I couldn't get the doubling above to work, I now realise it is because a_mortality is wrong in some way ###

weib_fit3 <- flexsurvreg(
  Surv(recyrs, censrec) ~ 1, 
  data = bc, 
  dist = "weibullPH",
  bhazard = exprate * 2 # Double the hazard in the dataframe
)

weib_stand3 <- standsurv(
  weib_fit3,
  type = "survival",
  t = seq(0, 25, length = 50),
  rmap = list(
    sex = sex,
    year = diag,
    age = agedays
  ),
  ratetable = a_mortality_v2, # Doubled hazard
  scale.ratetable = 365.25, 
  newdata = bc
)

plot3 <- plot(weib_stand3, expected = TRUE) + 
  geom_step(
    data = ggsurvfit::tidy_survfit(km),
    aes(x = time, y = estimate)
  )

[connor-ballinger]

Hi @irtimmins,

This one took me a bit longer to resolve - I didn't realise that the age for the ratetable rmap argument (but not the dataframe join) needed to be specified in days.

Using relsurv::transrate.hmd is clearly a better method than my previous hack, thanks for the tip.

Firstly, for the bc example, I think I have reproduced your plots. Doubling the ratetable results in a fit which looks reasonable. To confuse things, however, my early examples involved changing the scale.ratetable argument, which was clearly a mistake, and appears to have exponentially increased background mortality. Sorry, starting off with an error like that has added some unnecessary confusion to this problem. The bc example shows a much higher disease-specific hazard than my dataset and, for it to be an illustrative example of my own problem, I think applying any sufficiently high rate of background mortality (thereby lowering the disease-specific hazard) will result in a gradual deviation from the KM.

Unfortunately, my problem persists as the survival plot for my dataset stubbornly remains below the KM despite the new method of specifying background mortality. I have tried to replicate my dataset with some simulated data in the second script below - perhaps I should've done this from the start instead of using the bc example. The model incorporating background mortality shows a worse fit to the KM.

Thank you for your patience with this.

# bc example -------------------------------------------------------------------

## use Aus mortality

library(tidyverse)
library(flexsurv)
library(ggsurvfit)
library(relsurv) 
library(HMDHFDplus)

set.seed(1)

## Data ------------------------------------------------------------------------
data(bc, package = "flexsurv")

# mostly code from flexsurv vignette
bc$age <- rnorm(dim(bc)[1], mean = 65 - scale(bc$recyrs, scale=F), sd = 5)
## Create age at diagnosis in days - used later for matching to expected rates
bc$agedays <- floor(bc$age * 365.25)
## Create some random diagnosis dates between 01/01/1984 and 31/12/1989
bc$diag <- as.Date(floor(runif(dim(bc)[1], as.Date("01/01/1984", "%d/%m/%Y"), 
                               as.Date("31/12/1989", "%d/%m/%Y"))), 
                   origin="1970-01-01")
## Create sex (assume all are female)
bc$sex <- factor("female")
bc <- tibble(bc)
# Kaplan-Meier
km <- survfit(Surv(recyrs, censrec) ~ 1, data=bc)

### Alternate lifetable --------------------------------------------------------

# Using Iain's code:

## Download 1x1 Period lifetables from HMD 
## (fltper_1x1 for females, mltper_1x1 for males)
myHMDusername <- username # fill in login details
myHMDpassword <- password 

## create a popmort_files subdirectory if it does not exist
if(!dir.exists("popmort_files"))  dir.create("popmort_files")

## Countries that we will get ratetables 
country <- c("AUS")  # CHANGED
for(ctry in country){
  print(ctry)
  # Read in popmort file from HMD if the file does not exist
  if(!file.exists(paste0("popmort_files/",ctry,".f.txt"))){
    # Females
    HMD.f <- readHMDweb(CNTRY = ctry, item = "fltper_1x1", username = myHMDusername,
                        password = myHMDpassword)
    write.table(HMD.f, file=paste0("popmort_files/",ctry,".f.txt"))
    # Males
    HMD.m <- readHMDweb(CNTRY = ctry, item = "mltper_1x1", username = myHMDusername,
                        password = myHMDpassword)
    write.table(HMD.m, file=paste0("popmort_files/",ctry,".m.txt"))
  }
  ratetable <- relsurv::transrate.hmd(male = paste0("popmort_files/",ctry,".m.txt"),
                                      female = paste0("popmort_files/",ctry,".f.txt"))

  ## reshape lifetable to be a tidy data.frame, and convert rates to per person-year 
  # (needed if time scale of survival model is in years)

  ratetable.df <- as.data.frame.table(ratetable, responseName = "exprate") %>%
    mutate(exprate = 365.25 * exprate) %>%
    mutate(age = as.numeric(as.character(age)),
           year = as.numeric(as.character(year)))

  assign(paste0("ratetable.",ctry), ratetable)
  assign(paste0("ratetable.",ctry,".df"), ratetable.df)
}
ratetable.df <- tibble(ratetable.df) 

### Join excess mortality ------------------------------------------------------
## Obtain attained age and attained calendar year in (whole) years
bc <- bc |> 
  mutate(attained.age.yr = floor(age + recyrs),
         attained.year = lubridate::year(diag + rectime))

## merge in (left join) expected rates at event time
bc <- bc |> 
  left_join(
    ratetable.df,
    by = c("attained.age.yr"="age", 
           "attained.year"="year", 
           "sex"="sex")
  )

## A standard model ------------------------------------------------------------
weib_fit1 <- flexsurvreg(
  Surv(recyrs, censrec) ~ 1, 
  data = bc, 
  dist = "weibullPH"
)

## Standardise -----------------------------------------------------------------
weib_stand1 <- standsurv(
  weib_fit1,
  type = "survival",
  t = seq(0, 25, length = 50),
  rmap = list(
    sex = sex,
    year = diag,
    age = agedays
  ),
  ratetable = ratetable,
  scale.ratetable = 365.25,
  newdata = bc
)
## Plot ------------------------------------------------------------------------
plot1 <- plot(weib_stand1, expected = TRUE) + 
  geom_step(
    data = ggsurvfit::tidy_survfit(km),
    aes(x = time, y = estimate)
  )

## A relative survival model ---------------------------------------------------
weib_fit2 <- flexsurvreg(
  Surv(recyrs, censrec) ~ 1, 
  data = bc, 
  dist = "weibullPH",
  bhazard = exprate # incorporate population hazard
)

## Standardise -----------------------------------------------------------------
weib_stand2 <- standsurv(
  weib_fit2,
  type = "survival",
  t = seq(0, 25, length = 50),
  rmap = list(
    sex = sex,
    year = diag,
    age = agedays
  ),
  ratetable = ratetable,
  scale.ratetable = 365.25,
  newdata = bc
)
## Plot ------------------------------------------------------------------------
plot2 <- plot(weib_stand2, expected = TRUE) + 
  geom_step(
    data = ggsurvfit::tidy_survfit(km),
    aes(x = time, y = estimate)
  )

## Increase background mortality -----------------------------------------------

# create new ratetable with doubled rates.
ratetable_high <- ratetable
ratetable_high[, , ] <- ratetable_high[, , ] * 2

weib_fit3 <- flexsurvreg(
  Surv(recyrs, censrec) ~ 1, 
  data = bc, 
  dist = "weibullPH",
  bhazard = exprate * 2 # Double the hazard in the dataframe
)

weib_stand3 <- standsurv(
  weib_fit3,
  type = "survival",
  t = seq(0, 25, length = 50),
  rmap = list(
    sex = sex,
    year = diag,
    age = agedays
  ),
  ratetable = ratetable_high, # Doubled above
  scale.ratetable = 365.25, 
  newdata = bc
)

plot3 <- plot(weib_stand3, expected = TRUE) + 
  geom_step(
    data = ggsurvfit::tidy_survfit(km),
    aes(x = time, y = estimate)
  )

## My old method ---------------------------------------------------------------

weib_fit4 <- flexsurvreg(
  Surv(recyrs, censrec) ~ 1, 
  data = bc, 
  dist = "weibullPH",
  bhazard = exprate * 2 # Double the hazard in the dataframe
)

weib_stand4 <- standsurv(
  weib_fit4,
  type = "survival",
  t = seq(0, 25, length = 50),
  rmap = list(
    sex = sex,
    year = diag,
    age = agedays
  ),
  ratetable = ratetable,
  scale.ratetable = 365.25 * 2, # doubled scale
  newdata = bc
)

# large deviation, like my early plots
plot4 <- plot(weib_stand4, expected = TRUE) + 
  geom_step(
    data = ggsurvfit::tidy_survfit(km),
    aes(x = time, y = estimate)
  )

# Hazard -----------------------------------------------------------------------

# the mortality rate was doubled in Iain's method. 
# The marginal hazard is nearly but not precisely doubled. 
# My previous method increases the hazard by far more than a factor of 2.

pop_haz2 <- attr(weib_stand2, "expected") # unadjusted rates
pop_haz3 <- attr(weib_stand3, "expected") # doubled with ratetable[, , ] * scalar
pop_haz4 <- attr(weib_stand4, "expected") # doubled at the scale.ratetable arg

ggplot() +
  geom_line(data = pop_haz2, aes(time, exphaz), colour = "red") +
  geom_line(data = pop_haz3, aes(time, exphaz), colour = "blue") +
  geom_line(data = pop_haz4, aes(time, exphaz), colour = "orange")

Simulated data:

# Fake data --------------------------------------------------------------------

library(tidyverse)
library(flexsurv)
library(ggsurvfit)
library(relsurv) 
library(HMDHFDplus)

set.seed(1)

df_fake <- tibble(
  id = 1:300, # 300 patients
  age = round(rnorm(300, 60, 10)), # rounded for ratetable join
  sex = factor( # skewed male
    rbinom(300, 1, 0.8), 
    levels = c(0, 1), 
    labels = c("female", "male")
  ),
  date_baseline = seq.Date( # recruit for two years
    as.Date("2010-01-01"), 
    as.Date("2012-12-31"), 
    length.out = 300
  ),
  dead = rbinom(300, 1, 0.2), # 20% die
  # follow-up limited to 5 years
  follow_years = if_else(dead == 1, runif(300, 0, 5), 5), 
  year_baseline = year(date_baseline),
  agedays = age * 365.25
)

## Use ratetable ---------------------------------------------------------------

# Using Iain's code:

## Download 1x1 Period lifetables from HMD 
## (fltper_1x1 for females, mltper_1x1 for males)
myHMDusername <- user # fill in login details
myHMDpassword <- pwd

## create a popmort_files subdirectory if it does not exist
if(!dir.exists("popmort_files"))  dir.create("popmort_files")

## Countries that we will get ratetables 
country <- c("AUS")  # CHANGED
for(ctry in country){
  print(ctry)
  # Read in popmort file from HMD if the file does not exist
  if(!file.exists(paste0("popmort_files/",ctry,".f.txt"))){
    # Females
    HMD.f <- readHMDweb(CNTRY = ctry, item = "fltper_1x1", username = myHMDusername,
                        password = myHMDpassword)
    write.table(HMD.f, file=paste0("popmort_files/",ctry,".f.txt"))
    # Males
    HMD.m <- readHMDweb(CNTRY = ctry, item = "mltper_1x1", username = myHMDusername,
                        password = myHMDpassword)
    write.table(HMD.m, file=paste0("popmort_files/",ctry,".m.txt"))
  }
  ratetable <- relsurv::transrate.hmd(male = paste0("popmort_files/",ctry,".m.txt"),
                                      female = paste0("popmort_files/",ctry,".f.txt"))

  ## reshape lifetable to be a tidy data.frame, and convert rates to per person-year 
  # (needed if time scale of survival model is in years)

  ratetable.df <- as.data.frame.table(ratetable, responseName = "exprate") %>%
    mutate(exprate = 365.25 * exprate) %>%
    mutate(age = as.numeric(as.character(age)),
           year = as.numeric(as.character(year)))

  assign(paste0("ratetable.",ctry), ratetable)
  assign(paste0("ratetable.",ctry,".df"), ratetable.df)
}
ratetable.df <- tibble(ratetable.df) 

## Join excess mortality ------------------------------------------------------
# Obtain attained age and attained calendar year in (whole) years
df_fake <- df_fake |> 
  mutate(
    attained.age.yr = floor(age + follow_years),
    attained.year = lubridate::year(
      date_decimal(follow_years + decimal_date(date_baseline)) # bit weird...
    )
  )

## merge in (left join) expected rates at event time
df_fake <- df_fake |> 
  left_join(
    ratetable.df,
    by = c("attained.age.yr"="age", 
           "attained.year"="year", 
           "sex"="sex")
  )

km_fake <- survfit2(Surv(follow_years, dead) ~ 1, df_fake)

## Models ----------------------------------------------------------------------

fit1_fake <- flexsurvreg(
  Surv(follow_years, dead) ~ 1,
  data = df_fake,
  dist = "weibullPH"
)
stand1_fake <- standsurv(
  fit1_fake,
  t = seq(0.1, 15, by = 0.1),
  newdata = df_fake,
  rmap = list(age = agedays, 
              year = year_baseline,
              sex = sex),
  ratetable = ratetable,
  scale.ratetable = 365.25
)
p1_fake <- plot(stand1_fake, expected = TRUE) +
  geom_step(data = tidy_survfit(km_fake),
            aes(x = time, y = estimate)) +
  ylim(0.5, 1)

fit2_fake <- flexsurvreg(
  Surv(follow_years, dead) ~ 1,
  data = df_fake,
  dist = "weibullPH",
  bhazard = exprate
)
stand2_fake <- standsurv(
  fit2_fake,
  t = seq(0.1, 15, by = 0.1),
  rmap = list(age = agedays, # at diagnosis, not death
              year = year_baseline,
              sex = sex),
  ratetable = ratetable,
  scale.ratetable = 365.25,
  newdata = df_fake
)
p2_fake <- plot(stand2_fake, expected = TRUE) +
  geom_step(data = tidy_survfit(km_fake),
            aes(x = time, y = estimate)) +
  ylim(0.5, 1)

[irtimmins]

Thanks @connor-ballinger for this reproducible example, really helpful. I think we need to make sure that the "year" passed in the list using rmap argument is in date format.

If we use your code and replace

rmap = list(age = agedays, # at diagnosis, not death
              year = year_baseline,
              sex = sex),

with

rmap = list(age = agedays, # at diagnosis, not death
              year = date_baseline,
              sex = sex)

then I think the fit looks better, and the expected rate in blue looks more correct as well. I reckon that might be the issue?

For p2_fake I then get:

It's fair to say also that the quality of fit also depends on the model choice. I don't think it's necessarily true that higher background rates mean the model deviates from KM curve. If we return to the bc dataset using survexp.us.high.rate life tables with doubled rates, we can get good KM fit using a spline rather than a Weibull model:

spline_fit3 <- flexsurvspline(
  Surv(recyrs, censrec) ~ 1, 
  data = bc, 
  k=4,
  bhazard = 2*exprate  # Double the hazard in the matched population
)

spline_stand3 <- standsurv(
  spline_fit3,
  type = "survival",
  t = seq(0, 25, length = 50),
  rmap = list(
    sex = sex,
    year = diag,
    age = agedays
  ),
  ratetable = survexp.us.high.rate, # use the new ratetable
  scale.ratetable = 365.25 , # keep this the same as before
  newdata = bc
)

plot_spline <- plot(spline_stand3, expected = TRUE) + 
  geom_step(
    data = ggsurvfit::tidy_survfit(km),
    aes(x = time, y = estimate)
  )

plot_spline

[connor-ballinger]

Hi @irtimmins, thank you so much, I think you have solved my problem with the change to the rmap argument. Sorry that it took so long. For whatever reason, I seem to have struggled supplying the correct variables to i) the "by" argument in the left-join of the survival data and the lifetable datatframe, and ii) the "rmap" argument in standsurv, but I think everything is working now.

Below is my initial relative survival model (before troubleshooting) followed by the updated model relying on your corrections, which looks much better, even if the fit is not that great (both are Weibull models).

Thanks for the tip regarding the spline, I am fitting the six standard parametric models plus the three spline options. In choosing the final model, a plot of the hazards suggests the splines may display overfitting - the local maximum is not surprising but the subsequent, smaller, minimum is hard to explain (i.e. is it reasonable for the hazard to decrease after a maximum, only to have a second turning point, a minimum, and start increasing again, before leveling off?). This suggests to me that the generalised gamma or log-normal fits may be better choices.

None of the models produce hazard ratios as nice as in Sweeting 2023, and there is limited evidence of attenuation of the ratios towards 1. This may overstate the long-term treatment effects, particularly given the findings in Jennings 2024 on the overestimation of treatment effects even when marginal treatment effect waning (HR=1) is imposed.

I don't really have the knowledge/experience/competence to solve all these issues, and I can only rely on the generosity of others to a point, so I will have to leave some things unresolved and try to describe them in the limitations of the paper.

This is the first time I have had help with my code like this, so I'm ignorant of any etiquette - let me know if there is some way to repay the favour. flexsurv will be referenced in the paper that reports on this analysis.

Thanks.