…Given estimates and denominators

The following data are given in the R data frame ihdengland supplied in the package.

• inc_num,inc_denom: incidence of IHD by age

• mort_num,mort_denom: estimates of IHD-specific mortality by age

• prev_num,prev_denom: estimates of prevalence by age

and we wish to estimate case fatality by age, given these inputs. A selection of rows from the full data frame are shown here. We select only the data from one area and gender (Bristol, male) for this illustration. Remission rates are assumed to be zero, hence remission data are not included.

library(disbayes)
library(dplyr,quietly = TRUE)
ihdbristol <- ihdengland %>% filter(area=="Bristol", gender=="Male")
ihdbristol %>% filter(between(age, 50, 55))

## # A tibble: 6 × 9
##     age gender area    inc_num inc_denom prev_num prev_denom mort_num mort_denom
##   <dbl> <chr>  <fct>     <dbl>     <dbl>    <dbl>      <dbl>    <dbl>      <dbl>
## 1    50 Male   Bristol      28      8326      142       5810        4       7817
## 2    51 Male   Bristol      32      8294      140       5041        4       7771
## 3    52 Male   Bristol      32      7934      140       4433        4       7664
## 4    53 Male   Bristol      30      7245      144       3985        5       7498
## 5    54 Male   Bristol      28      6223      150       3699        5       7274
## 6    55 Male   Bristol      23      4874      158       3574        5       6990

These data were obtained from the Global Burden of Disease, 2017, and transformed in three ways to obtain the data seen here.

Firstly, the published point estimates and confidence intervals for each annual risk probability were converted to an implicit numerator and denominator. This is done by assuming that

• the incidence has been estimated as \(r_a^{(inc)}/n_a^{(inc)}\), where out of \(n_a^{(inc)}\) people who are alive at age \(a\), \(r_a^{(inc)}\) of these get IHD before age \(a+1\).

• the mortality has been estimated as \(r_a^{(mort)}/n_a^{(mort)}\), where out of \(n_a^{(mort)}\) people who are alive at age \(a\), \(r_a^{(mort)}\) of these die of IHD before age \(a+1\). These ratios are estimates of the true, underlying disease-specific mortalities \(dm_a\).

• the prevalence has been estimated as \(r_a^{(prev)}/n_a^{(prev)}\), where out of \(n_a^{(prev)}\) people who are alive at age \(a\), \(r_a^{(prev)}\) of them have IHD.

The published estimates and uncertainty interval are then used to derive the implicit denominator and numerator, using the method described in the next subsection.

Secondly, estimates for five-year age groups were converted to one-year groups. Most simply this could be done by assuming the counts are equal for each year of age, but a more sophisticated method was used based on “temporal disaggregation” using one of the methods implemented by Sax and Steiner in the tempdisagg R package (Sax and Steiner (2013)). This ensures that the disaggregated counts vary smoothly with age, while preserving the aggregate totals.

Thirdly, counts from smaller geographical areas were simply added up to obtain the corresponding counts for the larger areas (English city regions) required in this analysis.

ihdbristol[ihdbristol$age %in% 50:55, ]

## # A tibble: 6 × 9
##     age gender area    inc_num inc_denom prev_num prev_denom mort_num mort_denom
##   <dbl> <chr>  <fct>     <dbl>     <dbl>    <dbl>      <dbl>    <dbl>      <dbl>
## 1    50 Male   Bristol      28      8326      142       5810        4       7817
## 2    51 Male   Bristol      32      8294      140       5041        4       7771
## 3    52 Male   Bristol      32      7934      140       4433        4       7664
## 4    53 Male   Bristol      30      7245      144       3985        5       7498
## 5    54 Male   Bristol      28      6223      150       3699        5       7274
## 6    55 Male   Bristol      23      4874      158       3574        5       6990

Smoothly disaggregating data from age groups to years of age

As an aside, here we show how the temporal disaggregation procedure is done in R. We have a disease count outcome that is measured for 5-year age groups, as in the following data frame.

dat <- data.frame(agefrom=seq(0,20,5), ageto=seq(4,24,5),
                  measure=c(15,20,24,35,29))
dat

##   agefrom ageto measure
## 1       0     4      15
## 2       5     9      20
## 3      10    14      24
## 4      15    19      35
## 5      20    24      29

The following code shows how to produce a smoothly disaggregated version of this count, and compares the smooth disaggregation to a crude disaggregation where the 5-year counts are simply divided by 5 to estimate one-year counts. We can check that the smoothly-disaggregated measures add up to the original data that were provided for the 5-year age groups.

if (requireNamespace("tempdisagg")) { 
   disagg_crude <- rep(dat$measure/5, each=5)
   disagg_smooth <- predict(tempdisagg::td(dat$measure ~ 1, to=5, method="fast"))
   ageyr <- 0:24
   plot(ageyr, disagg_crude, type="l", xlab="Age", ylab="Measure", ylim=c(0,8))
   lines(ageyr, disagg_smooth, col="blue")
   agegroup <- cut(ageyr, seq(0,25,5), right = FALSE)
   tapply(disagg_smooth, agegroup, sum)
}

## Loading required namespace: tempdisagg

##   [0,5)  [5,10) [10,15) [15,20) [20,25) 
##      15      20      24      35      29

See the documentation for the tempdisagg package, including the vignettes and the help page for the td function in this package, for more information about smooth temporal disaggregation. It is unclear how far the method generalises, e.g. to intervals of varying lengths.

Determining numerators and denominators from point and/or interval estimates

Sometimes a point estimate \(\hat{p}\) for a quantity such as incidence is published alongside an (e.g. 95%) interval estimate \((\hat{p}^{(lower)},\hat{p}^{(upper)})\). The interval estimate can be assumed to express the uncertainty associated with the estimate. Such information can be converted to an implicit numerator \(r\) and denominator \(n\) as follows. We assume the point and interval estimate are summaries of a Beta posterior distribution which has been obtained by combining a vague prior with an observation of \(r\) events occurring out of a sample of \(n\) individuals. If a Beta(0,0) prior is used, then the posterior is known to be Beta(r, n-r). We can then search for the best-fitting Beta(r,n-r)$ distribution which has median \(\hat{p}\) and (2.5,97.5) quantiles \((\hat{p}^{(lower)},\hat{p}^{(upper)})\), and set \(r=a, n=a+b\). A utility to perform this search is provided in the SHELF package for expert elicitation (O’Hagan et al. (2006)).

In other applications, the denominators might be directly available, e.g. we might know that the prevalence was estimated from a survey of \(r_a^{(prev)}\) people, or we might assume that the incidence and mortality were obtained from a comprehensive registry of the population, in which case the denominator would equal the population size (within the specific stratum, e.g. as defined by age and gender, that the estimate refers to).

The uncertainty inherent in the information supplied about each of incidence, prevalence and mortality is measured by the denominator. In the example above, this is equated to the size of the population used to produce the estimate. However, if we also suspected that one of the data sources may be biased, but were unsure about the direction of bias, we could downweight that data source by multiplying both the numerator and denominator by the same amount, e.g. 0.5 if we wanted to give a data source half its original weight. Note that if counts are scaled in this way, they should then be rounded to the nearest integer.

Bayesian modelling process

The four steps of the Bayesian modelling process are then implemented as follows:

1. write down the theoretical disease model, as given above

2. write down the statistical model for the data. All count data are assumed to arise from a Binomial distribution with the corresponding denominator, and a probability parameter which is a function of the parameters in the theoretical disease model.

   • (incidence) \(r_a^{(inc)} \sim ~ Binomial(n_a, 1 - P_{a11})\), where \(P_{a11}\) is the annual probability of remaining free of the disease.

   • (mortality) \(r_a^{(mort)} \sim ~ Binomial(n_a, dm_a)\), where the disease-specific mortality \(dm_a\) is a deterministic function of the incidences and case fatalities \(\{i_j,f_j\}\) for ages \(j\) up to \(a\), as described in the theoretical disease model.

   • (prevalence) \(r_a^{(prev)} \sim ~ Binomial(n_a^{(prev)}, pr_a)\), where \(pr_a\) is the true prevalence, defined as a deterministic function of the incidences and case fatalities.

   • (remission, if assumed possible) \(r_a^{(rem)} \sim ~ Binomial(n_a^{(rem)}, P_{a21})\), where \(P_{a21}\) is the probability that someone with the disease will have recovered from the disease one year from now, defined as a deterministic function of the remission rates, case fatalities and incidences.

3. define prior distributions for the unknown parameters, as explained in the next section.

4. compute the posterior distribution \(p(\theta | \mathbf{y})\) for parameters \(\theta = \{i_a,f_a\}\) given data \(\mathbf{y} = \{i_a, r_a, n_a\}\)

Alternative models / prior assumptions

In example such as this one, where latent quantities are estimated from indirect data, it is important to consider what substantive prior information is available about the unknown quantities, and what influence the assumed prior has on the resulting estimates.

In these examples, the case fatality is only informed indirectly. The disbayes package implements alternative models for how case fatality depends on age. The same methods can be used for how incidence depends on age.

1. the case fatality rate \(f_a\) for each year of age \(a\) is assumed to have an independent vague prior distribution.

2. the case fatality rate is assumed to be a smooth function \(g()\) of age, \(f_a = g(a)\). This smooth function is defined by a spline basis \(log(f_a) = \sum_k \beta_k g_k(a)\), where \(g_k()\) are basis functions. A “thin plate” spline is used, following the default in the mgcv R package (Wood), and the amount of smoothness is determined automatically through a hierarchical prior on the coefficients.

3. the case fatality rate is a smoothly increasing function of age.

4. the case fatality rate is constant with age.

In addition, for all ages below a given cut-off age \(a_{base}\), case fatalities are assumed to be equal to a constant \(f_{base}\). This cut-off age \(a_{base}\) needs to be supplied when calling the disbayes function. \(f_{base}\) does not need to be supplied however, and is estimated from the data under the assumption that the dependence on age is a smooth function.

Model (2) is more realistic. However it may not give sensible estimates if the information provided by the data is too weak. The results of fitting Model (1) can help to diagnose where the indirect information on case fatality provided by the rest of the data is weaker or stronger. If the information on \(f_a\) for a particular age is too weak, then the posterior distribution will be identical to the prior. In those cases, substantive prior information about case fatality at that age is necessary for the estimates to be meaningful.

This information might come from nearby ages that are better-informed, through extrapolation from the smoothing model. However extrapolating in this way is only feasible for a limited time span (perhaps around 10-20 years) before stronger assumptions are necessary, e.g. that case fatality is constant below a certain age.

The increasing model (3) may be required in cases where the data alone are not sufficient to identify an increasing shape with age. A constant model (4) may be necessary if the data on age variations are very weak, or simply if this model is thought to be plausible.

Alternative models of how the rates depend on age are also supported for incidence and remission rates.

Running disbayes: optimisation method

The following call to the disbayes function is then made to fit the model to the ihdbristol data. The argument eqage=40 defines an assumption that case fatality is constant for all ages below 40.

dbres <- disbayes(data = ihdbristol, age = "age",
                 inc_num = "inc_num", inc_denom = "inc_denom", 
                 prev_num = "prev_num", prev_denom = "prev_denom",
                 mort_num = "mort_num", mort_denom = "mort_denom",
                 eqage = 40)

For comparison, we also fit an “unsmoothed” version of this model where the age-specific case fatality rates have independent priors, rather than being related through a spline function. This is specified with cf_model="indep".

dbresu <- disbayes(data = ihdbristol, age = "age",
                 inc_num = "inc_num", inc_denom = "inc_denom", 
                 prev_num = "prev_num", prev_denom = "prev_denom",
                 mort_num = "mort_num", mort_denom = "mort_denom",
                 cf_model = "indep", 
                 eqage = 40)

The optimisation method produces accurate point estimates (defined by the posterior mode), and approximate credible intervals.

Occasionally the optimisation may produce an error message that mentions chol(H), but still produce a point estimate. This indicates that the level of uncertainty cannot be determined, and casts doubt on the validity of the point estimate. In these cases it may work if the model is simplified, e.g. by fixing the level of smoothness for the age-dependence of the rates by using the hp_fixed option to disbayes

Running disbayes: MCMC method

If more accurate credible intervals are desired, then MCMC sampling can be used instead by setting method="mcmc", as follows. This is more computationally intensive, typically taking a few minutes per fit for the standard disbayes model.

options(mc.cores = parallel::detectCores())
dbresm <- disbayes(data = ihdbristol, age = "age",
                 inc_num = "inc_num", inc_denom = "inc_denom", 
                 prev_num = "prev_num", prev_denom = "prev_denom",
                 mort_num = "mort_num", mort_denom = "mort_denom",
                 method="mcmc", chains=2, iter=1000,
                 eqage = 40)

The sampling by rstan can be controlled, e.g. by setting the number of chains to run (perhaps on parallel cores) and the number of iterations - see help(disbayes) and the Stan documentation for full details. The number of iterations iter should be high enough that the results are not affected by Monte Carlo error.

After using MCMC, we should check that the sampler has converged, by examining the “trace plot” of the simulation progress for selected parameters. The simulated chains should mix together and look like white noise. Also if the disbayes call resulted in many reports of “divergent transitions”, then the chains are unlikely to have converged to a valid posterior. If convergence hasn’t occurred, then the model probably needs to be simplified (e.g. by increasing eqage, or fixing the degree of smoothness using hp_fixed), or the data are insufficient.

rstan::traceplot(dbres$fit, pars=paste0("cf[", 60:65, "]"))

Other options to disbayes

Full documentation is given in help(disbayes). Some more useful options include

• cf_model. Model for how case fatality varies with age. As well as the default spline smooth model and unsmoothed models illustrated above, cf_model="increasing" fits a model where case fatality is restricted to increase with age, and cf_model="const" specifies a constant case fatality with age.

• inc_model. By default, incidence rates are smoothed as a function of age using a spline, as with case fatality rates. inc_model="indep" estimates independent incidence rates for each age.

• rem_model. If no remission data are supplied, then the disease is assumed to have zero remission rate. If remission data are supplied, then remission rates are estimated. By default, these are assumed to be a constant function of age (rem_model="const"). To estimate independent remission rates for each age, supply rem_model="indep". To model remission rates as a smooth spline function of age, use rem_model="smooth".

• hp_fixed. A list of “hyperparameters” to fix at constant values. If using smooth spline functions, this option may be required if the appropriate level of smoothness for the spline cannot be determined from the data. When using the default optimisation method, this can lead to an error message mentioning "Error in chol.default(-H)". In these cases, the smoothness can be fixed in advance, e.g. hp_fixed = list(scf=1) to fix the case fatality smoothness parameter to a typical value of 1. Lower values give more smoothness, higher values give wigglier curves. You may need to experiment with different values to determine what gives a plausible degree of smoothness.

Hierarchical models

disbayes_hier is used to fit a joint model for data from multiple areas.

Hierarchical models with additive area and gender effects

disbayes_hier can also be used to fit a joint model for data from multiple areas and different genders. It is assumed that the relative case fatality between women and men is the same for each area. While this is labelled gender it can be used for any single binary predictor, though for chronic diseases, we will usually want to consider the effect of gender.

Trends through time (non-hierarchical models only)

disbayes can include information about trends through calendar time, which can be age dependent. The ratio of incidence between previous years and the current year (the one represented by the data) can be supplied as a matrix, with one row per year of age and one column per calendar year. See the disbayes, arguments inc_trend and cf_trend.