Mortality adjustments

Notation

We denote functions of survival, hazard and cumulative hazard as \(S(t), h(t)\) and \(H(t)\). There are ‘unadjusted’ values for relevant endpoints (e.g. PFS) dependent on the model structure (PSM, STM-CF or STM-CR). We assume that there is a ‘general background mortality’, against which any extrapolations and derivations of RMD should be constrained, denoted with subscript ‘gen’. Consequently there are ‘adjusted’ values, denoted with superscript ‘adj’.

Background mortality survival and hazard functions

Our lifetable gives us decreasing \(l_x\) values at points between \(t=0, ..., t_{max}\), with \(l_{t_{max}}=0\). Let us define the background mortality survival function as follows.

\[ S_{gen}(t) = 1 - R \left(1 - \frac{l_{x_0 + t}}{l_{x_0}} \right) \]

where \(x_0\) is the mean age of the study population at baseline, \(R\) is a Standardized Mortality Ratio, and \(l_x\) is an appropriate lifetable.

We can also derive the average background mortality hazard for \(t\) in the range \([t_1, t_2)\).

\[ h_{gen}(t) = \frac{\log[S_{gen}(t_1)] - \log[S_{gen}(t_2)]}{t_2 - t_1} \]

Approach 1 - Constraining only the survival function

Under this approach, the background mortality is assumed to limit the survival function as follows. There is no direct constraint to the hazard.

\[ S^{adj}(t) = \min \left[ S(t), S_{gen}(t) \right] \]

This is a simple approach, and has been used at least one NICE technology appraisal (National Institute for Health and Care Excellence 2022):

“Patient longevity is always the lesser of values generated from the disease-specific survival curve (after adjustment for treatment and functional status) and the survival curve for the general population according to age and sex.”

However, this approach does not ensure that the hazard is at least as great as the hazard of background mortality.

Approach 2 - Constraining the hazard function

A more rigorous adjustment would be that the hazard of background mortality acts as a constraint to the unadjusted hazard (Sweeting et al. 2023).

\[ h^{adj}(t) = \max \left[ h(t), h_{gen}(t) \right] \] \[ \implies S^{adj}(t) = \exp \left[ - \int_0^t h^{adj}(u) du \right] \]

In either case, we must assume that the original dataset was NOT subject to background mortality. Although unlikely to be true, this is a common, pragmatic assumption in cost-effectiveness models as long as background mortality is relatively insignificant during trial follow-up.

Approach 3 - Modeling excess hazard in the original dataset

The above approaches work by adjusting extrapolations made from models fitted to data assumed to be not subject to background mortality. If - as is likely - the population in the dataset were in fact subject to background mortality, then it would be better to model that dataset using excess hazard methods rather than to seek to make adjustments to extrapolations after the fact. Further discussion of survival extrapolation incorporating general population mortality is provided by Sweeting et al. (2023).

Set-up

First we load the packages we need - all of which are suggested for or imported to psm3mkv. If you haven’t installed psm3mkv yet, please see the installation instructions to install it with its dependencies. Of these, we will making further use of Wickham et al. (2023), Riffe (2015) and Müller and Wickham (2023).

library("dplyr")
library("psm3mkv")
library("tibble")

Creating the life table

In order to apply constraints to background mortality, we need some background mortality data, in the form of a lifetable. We can take this from the Human Mortality Database using the HMDHFDplus package mentioned above (with thanks to Robert Hettle for the recommendation). The lifetable will need to start from an assumed age at baseline. Lifetables are constructed by time in years.

# HMD HFD package
library("HMDHFDplus")

# Assumed population age at baseline (time=0)
baseage <- 51.0

# Mortality data - England & Wales, Female, 2019
mort <- readHMDweb(
  CNTRY = "GBRTENW",
  item = "fltper_1x1",
  username = "",
  password = ""
) |>
  filter(Year == 2019) |>
  select(Age, lx, dx) |>
  mutate(Timey = ceiling(Age - baseage)) |>
  filter(Timey >= 0)

# The table needs to end with lx=0
mort <- add_row(mort, Age = 111, lx = 0, Timey = 60)

You will see that the above code cannot run without a login for the Human Mortality Database. Alternatively, we could just make up a mortality table.

mort <- tibble(
  Timey = 0:30,
  lx = 10000 * exp(-0.03 * Timey^1.1),
  dx = lx - lead(lx),
  qx = dx / lx
)

However we do it, once we have the mortality data, we can apply the SMR (\(R\)) and derive a lifetable from time zero as required.

# Assumed Standardized Mortality Ratio
SMR <- 2

# Recalculate the lifetable with the SMR applied
mort$adjlx <- mort$lx
for (i in 2:length(mort$lx)) {
  mort$adjlx[i] <- mort$adjlx[i - 1] * (1 - SMR * mort$qx[i - 1])
}

# Ensure lx>=0
mort$adjlx[mort$adjlx < 0] <- 0

# Create and view lifetable
ltable <- tibble(lttime = mort$Timey, lx = mort$adjlx)
head(ltable)
#> # A tibble: 6 × 2
#>   lttime     lx
#>    <int>  <dbl>
#> 1      0 10000 
#> 2      1  9409.
#> 3      2  8774.
#> 4      3  8151.
#> 5      4  7553.
#> 6      5  6985.

Obtaining and fitting distributions to patient-level data

Next we get some patient-level data and fit the PSM and STMs. So far, we are closely following vignette("example"). Time is assumed to be recorded in the patient-level data in weeks, with approximately 52.18 weeks per year.

# Get some data
bosonc <- create_dummydata("flexbosms")

# We'll make the durations a lot longer
# so that they will be definitely constrained by the lifetable
bosonc <- bosonc |>
  mutate(
    pfs.durn = 20 * pfs.durn,
    os.durn = 20 * os.durn,
    ttp.durn = 20 * ttp.durn
  )

# Fit parametric models to each endpoint
allfits <- fit_ends_mods_par(bosonc)

# Bring together preferred fits for each endpoint
params <- list(
  ppd = find_bestfit(allfits$ppd, "aic")$fit,
  ttp = find_bestfit(allfits$ttp, "aic")$fit,
  pfs = find_bestfit(allfits$pfs, "aic")$fit,
  os = find_bestfit(allfits$os, "aic")$fit,
  pps_cf = find_bestfit(allfits$pps_cf, "aic")$fit,
  pps_cr = find_bestfit(allfits$pps_cr, "aic")$fit
)

Making the projections

We are skipping over any internal or external validation for the purposes of this vignette, and jump instead straight to making the survival projections. We will assume a 10 year time horizon.

First we derive a projection using integral formulae, without lifetable constraints.

# Set time horizon
thoz <- 10

# Run the calculations
proj1 <- calc_allrmds(bosonc, dpam = params, Ty = thoz)

# Present the results
res1 <- proj1$results |> mutate(method = "int", lxadj = "no", psmmeth = "simple")
res1 |> mutate(
  across(c(pf, pd, os), ~ num(.x, digits = 1))
)
#> # A tibble: 3 × 7
#>          pf        pd        os model  method lxadj psmmeth
#>   <num:.1!> <num:.1!> <num:.1!> <chr>  <chr>  <chr> <chr>  
#> 1     274.4     118.3     392.7 PSM    int    no    simple 
#> 2     270.2     119.7     389.9 STM-CF int    no    simple 
#> 3     270.2     123.8     393.9 STM-CR int    no    simple

Next we derive a projection using a discretization approximation, still without any lifetable constraints.

# Run the calculation
proj2 <- calc_allrmds(bosonc, dpam = params, Ty = thoz, rmdmethod = "disc")

# Present the results
res2 <- proj2$results |>
  mutate(method = "disc", lxadj = "no", psmmeth = "simple")
res2 |> mutate(
  across(c(pf, pd, os), ~ num(.x, digits = 1))
)
#> # A tibble: 3 × 7
#>          pf        pd        os model  method lxadj psmmeth
#>   <num:.1!> <num:.1!> <num:.1!> <chr>  <chr>  <chr> <chr>  
#> 1     274.3     118.3     392.6 PSM    disc   no    simple 
#> 2     270.1     119.7     389.8 STM-CF disc   no    simple 
#> 3     270.1     123.7     393.8 STM-CR disc   no    simple

Discretization has proven fairly accurate. The STM-CR estimate of mean time alive has reduced by just 0.1 weeks for example, from 393.9 to 393.8 weeks.

Finally, we use a discretization approximation and apply the lifetable constraints.

# Run the calculations
proj3 <- calc_allrmds(bosonc, dpam = params, Ty = thoz, rmdmethod = "disc", lifetable = ltable)

# Present the results
res3 <- proj3$results |>
  mutate(method = "disc", lxadj = "yes", psmmeth = "simple")
res3 |> mutate(
  across(c(pf, pd, os), ~ num(.x, digits = 1))
)
#> # A tibble: 3 × 7
#>          pf        pd        os model  method lxadj psmmeth
#>   <num:.1!> <num:.1!> <num:.1!> <chr>  <chr>  <chr> <chr>  
#> 1     235.4     105.5     340.9 PSM    disc   yes   simple 
#> 2     230.1     107.6     337.7 STM-CF disc   yes   simple 
#> 3     230.1     111.6     341.6 STM-CR disc   yes   simple

# Proportion of time alive spent in PF
proppf3 <- res3$pf / res3$os

There is now some difference between the STM-CR compared to the STM-CF and PSM structures. All structures estimate fairly similar RMD estimates of being alive (range: 337.7 to 341.6 weeks). The estimate of mean time alive in the STM-CR structure has reduced by 52.3 weeks from 393.9 to 341.6 weeks.

The proportion of that time in the PF rather than PD state varies somewhat (PSM: 69%, STM-CF: 68.1%, STM-CR: 67.3%). Given that the hazard of mortality would typically be lower before rather than after progression, the pre-progression mortality is likely to be more affected by any mortality constraint that post-progression mortality. The STM-CR model of post-progression survival depends on time from progression and is therefore rather different from the STM-CF and PSM structures, which model mortality based on time from baseline.

We can check whether using a ‘complex’ type of PSM structure rather than the default ‘simple’ affects RMD estimates from the PSM structure.

# Run the calculations
proj4 <- calc_allrmds(bosonc, dpam = params, Ty = thoz, psmtype = "complex", rmdmethod = "disc", lifetable = ltable)

# Present the results
res4 <- proj4$results |>
  mutate(method = "disc", lxadj = "yes", psmmeth = "complex")
res4 |> mutate(
  across(c(pf, pd, os), ~ num(.x, digits = 1))
)
#> # A tibble: 3 × 7
#>          pf        pd        os model  method lxadj psmmeth
#>   <num:.1!> <num:.1!> <num:.1!> <chr>  <chr>  <chr> <chr>  
#> 1     231.0     103.7     334.7 PSM    disc   yes   complex
#> 2     230.1     107.6     337.7 STM-CF disc   yes   complex
#> 3     230.1     111.6     341.6 STM-CR disc   yes   complex

# Proportion of time alive spent in PF
proppf4 <- res4$pf / res4$os

All structures estimate fairly similar RMD estimates of being alive (range: 334.7 to 341.6 weeks). The proportion of that time in the PF rather than PD state varied a little by model structure (PSM: 69%, STM-CF: 68.1%, STM-CR: 67.3%).

Comparing the results

A summary of the STM-CF estimates of mean time in PF state and mean time alive is given in the table below.

In this fictional case, the application of a lifetable constraint had fairly sizeable effects on the results - although this was by design of this demonstration.

A summary of the STM-CR estimates of mean time in PF state and mean time alive is given in the table below.

The application of a lifetable constraint had similar and fairly sizeable effects on the results here also.

A summary of the PSM estimates of mean time in PF state and mean time alive is given in the table below.

The change in the RMDs for the ‘complex’ PSM structure compared to the ‘simple’ PSM structure were -4.3, -1.9 and -6.2 weeks respectively for time in PF state, PD state and time alive.

Possible extensions

Not applied here, the package also allows application of discounting through the discrate optional call to calc_allrnds(). This work could also be easily extended through calculating RMDs or Life Years that are Quality Adjusted (QALYs).