Refresher - Burden of disease with relative risk

Function call - Hard coded

results_pm_copd <- attribute_health(
  erf_shape = "log_linear", # shape of the exposure-response function (ERF)
  rr_central = 1.369, # relative risk (RR) central estimate
  rr_increment = 10,  # increment for which relative risk is valid (in \\mu g / m^3)
  exp_central = 8.85, # PM2.5 exposure (in \\mu g / m^3) (here: population-weighted)
  cutoff_central = 5, # cutoff (in \\mu g / m^3) below which no health effects occur 
  bhd_central = 30747 # baseline health data (BHD; here: COPD incidence)
)

For alternative ERF shapes see the function documentation of attribute_health().

Function call - Pre-loaded data

healthiar comes with some example data that start with exdat_ that allow you to test functions. Some of these example data will be used in some examples in this vignette.

Now we call attribute_health with input data from the healthiar example data. Note that we can easily provide input data to the function argument using the $ operator.

results_pm_copd <- attribute_health(
  erf_shape = "log_linear",
  rr_central = exdat_pm$relative_risk, 
  rr_increment = 10, 
  exp_central = exdat_pm$mean_concentration,
  cutoff_central = exdat_pm$cut_off_value,
  bhd_central = exdat_pm$incidence
)

Output structure

Every attribute_health() output consists of two lists (“folders”)

• health_main contains the main results

• health_detailed detailed results and additional info about the assessment

  • results_raw tibble that contains the detailed results

  • input_args contains a complete list of all arguments used in the assessment (including some internal arguments)

  • input_table tibble that contains all input data entered by the user or derived for the calculations

  • results_by_...

    • geo_id_micro: contains results by the inputted or default micro (e.g. high resolution) geo IDs

    • sex: contains results by sex (only present if function argument sex specified)

    • age_group: contains results by age groups (only present if function argument age_group sepcified)

NOTE: attribute_lifetable() creates additional output that is specific to life table calculations

Main results

Let’s inspect the main results

There exist different, equivalent ways of accessing the output

• With $ operator: results_pm_copd$health_main$impact_rounded (as in the example above)

• By mouse: go to the Environment tab in RStudio and click on the variable you want to inspect, and then open the health_main results table

• With [[]] operator results_pm_copd[["health_main"]]

• With pluck() & pull(): use the purrr::pluck function to select a list and then the dplyr::pull function extract values from a specified column, e.g. results_pm_copd |> purrr::pluck("health_main") |> dplyr::pull("impact_rounded")

results_pm_copd$health_main
#> # A tibble: 1 × 24
#>   geo_id_micro erf_ci  exp_ci  bhd_ci  cutoff_ci exp_category sex   age_group
#>   <chr>        <chr>   <chr>   <chr>   <chr>            <int> <chr> <chr>    
#> 1 a            central central central central              1 all   all      
#> # ℹ 16 more variables: impact <dbl>, impact_rounded <dbl>, approach_risk <chr>,
#> #   rr_increment <dbl>, erf_shape <chr>, prop_pop_exp <dbl>, exp_length <int>,
#> #   exp_type <chr>, exp <dbl>, bhd <dbl>, cutoff <dbl>, rr <dbl>,
#> #   is_lifetable <lgl>, pop_fraction_type <chr>, rr_at_exp <dbl>,
#> #   pop_fraction <dbl>

It is a table of the format tibble (very similar to a data.frame) of 3 rows and 23 columns. Let’s zoom in on some relevant aspects

results_pm_copd$health_main |> 
  dplyr::select(exp, bhd, rr, erf_ci, pop_fraction, impact_rounded) |> 
  knitr::kable() # For formatting reasons only: prints tibble in nice layout

Interpretation: this table shows us that exposure was 8.85 \(\mu g/m^3\), the baseline health data (bhd_central) was 30747 (COPD incidence in this instance). The 1st row further shows that the impact attributable to this exposure using the central relative risk (rr_central) estimate of 1.369 is 3502 COPD cases, or ~11% of all baseline cases.

Some of the most results columns include:

• impact_rounded rounded attributable health impact/burden
• impact raw impact/burden
• pop_fraction population attributable fraction (PAF) or population impact fraction (PIF)

NOTE: the main output contains more columns that provide additional information about the assessment, such as intermediate results and assessment specifications, …

Function call

results_pm_copd <- attribute_health(
    erf_shape = "log_linear",
    rr_central = 1.369, 
    rr_lower = 1.124, # lower 95% confidence interval (CI) bound of RR
    rr_upper = 1.664, # upper 95% CI bound of RR
    rr_increment = 10, 
    exp_central = 8.85, 
    exp_lower = 8, # lower 95% CI bound of exposure
    exp_upper = 10, # upper 95% CI bound of exposure
    cutoff_central = 5,
    bhd_central = 30747, 
    bhd_lower = 28000, # lower 95% confidence interval estimate of BHD
    bhd_upper = 32000 # upper 95% confidence interval estimate of BHD
) 

Detailed results

Let’s inspect the detailed results:

Each row contains the estimated attributable cases (impact_rounded) obtained by the input data specified in the columns ending in “_ci” and the other calculation pathway specifications in that row (not shown).

• The 1st contains the estimated attributable impact when using the central estimates of relative risk, exposure and baseline health data

• The 2nd row shows the impact when using the central estimates of the relative risk, exposure in combination with the lower estimate of the baseline health data

• …

NOTE: only 9 of the 27 possible combinations are displayed due to space constraints

NOTE: only a selection of columns are shown

Refresher - Burden of disease with absolute risk

Function call

results_noise_ha <- attribute_health(
  approach_risk = "absolute_risk", # default is "relative_risk"
  exp_central = c(57.5, 62.5, 67.5, 72.5, 77.5), # mean of the exposure categories
  pop_exp = c(387500, 286000, 191800, 72200, 7700), # population exposed per exposure category
  erf_eq_central = "78.9270-3.1162*c+0.0342*c^2" # exposure-response function
)

The erf_eq_central argument can digest other types of functions (see section [Example: attribute_health() with user-defined ERF])

Main results

Results per noise exposure band

results_noise_ha$health_detailed$results_raw

Alternatively, it’s also possible to only assess the impacts for a single noise exposure band

results_noise_ha <- attribute_health(
  approach_risk = "absolute_risk",
  exp_central = 57.5,
  pop_exp = 387500,
  erf_eq_central = "78.9270-3.1162*c+0.0342*c^2"
)

Function call

results_iteration <- attribute_health(
    # Names of Swiss cantons
    geo_id_micro = c("Zurich", "Basel", "Geneva", "Ticino", "Jura"),
    # Names of languages spoken in the selected Swiss cantons
    geo_id_macro = c("German","German","French","Italian","French"),
    rr_central = 1.369,
    rr_increment = 10, 
    cutoff_central = 5,
    erf_shape = "log_linear",
    exp_central = c(11, 11, 10, 8, 7),
    bhd_central = c(4000, 2500, 3000, 1500, 500)
)

In this example we want to aggregate the lower-level geographic units (municipalities) by the higher-level language region ("German", "French", "Italian")

Main results

The main output contains aggregated results

In this case health_main contains the cumulative / summed number of stroke cases attributable to PM2.5 exposure in the 5 geo units, which is 1717 (using a relative risk of 1.369).

Detailed results

The geo unit-specific information and results are stored under health_detailed>results_raw

health_detailed also contains impacts obtained through all combinations of input data central, lower and upper estimates (as usual), besides the results per geo unit (not shown above)

Function call

data <- exdat_noise |> 
  ## Filter for urban and rural regions
  dplyr::filter(region == "urban" | region == "rural")

results_iteration_ar <- attribute_health( 
    # Both the rural and urban areas belong to the higher-level "total" region
    geo_id_macro = "total",
    geo_id_micro = data$region,
    approach_risk = "absolute_risk",
    exp_central = data$exposure_mean,
    pop_exp = data$exposed,
    erf_eq_central = "78.9270-3.1162*c+0.0342*c^2"
)

NOTE: the length of the input vectors fed to geo_id_micro, exp_central, pop_exp must match and must be

(number of geo units) x (number of exposure categories) = 2 x 5 = 10,

because we have 2 geo units ("rural" and "urban") and 5 exposure categories.

Main results

health_main contains the aggregated results (i.e. sum of impacts in rural and urban areas)

Detailed results

Impact by geo unit, in this case impact in the rural and in the urban area

Function call

1. Use attribute_health() to calculate burden of scenarios A & B

scenario_A <- attribute_health(
    exp_central = 8.85,   # EXPOSURE 1
    cutoff_central = 5, 
    bhd_central = 25000,
    approach_risk = "relative_risk",
    erf_shape = "log_linear",
    rr_central = 1.118,
    rr_increment = 10)

scenario_B <- attribute_health(
    exp_central = 6,     # EXPOSURE 2
    cutoff_central = 5, 
    bhd_central = 25000,
    approach_risk = "relative_risk",
    erf_shape = "log_linear",
    rr_central = 1.118,
    rr_increment = 10)

Alternatively, the function attribute_mod() can be used to modify an existing scenario, e.g. scenario_A

scenario_B <- attribute_mod(
  output_attribute = scenario_A, 
  exp_central = 6
)

2. Use compare() to compare scenarios A & B

results_comparison <- compare(
  approach_comparison = "delta", # or "pif" (population impact fraction)
  output_attribute_scen_1 = scenario_A,
  output_attribute_scen_2 = scenario_B
)

The default value for the argument approach_comparison is "delta". The alterntive is "pif" (population impact fraction). See the function documentation of compare() for more details.

Main results

Detailed results

The compare() results contain two additional outputs in addition to those we have already seen

• health_detailed

  • scen_1 contains results of scenario 1 (scenario A in our case)

  • scen_2 contains results of scenario 2 (scenario B in our case)

Function call

results_pm_copd_summarized <- 
  summarize_uncertainty(
    output_attribute = results_pm_copd,
    n_sim = 100
)

Main results

Detailed results

The folder uncertainty_detailed contains all single simulations. Let’s look at the impact of the first 10 simulations.

Function call

results_pm_copd <- monetize(
    output_attribute = results_pm_copd,
    discount_shape = "exponential",
    discount_rate = 0.03,
    n_years = 5,
    valuation = 50000 # E.g. EURO
)

Main results

We see that the monetized impact (discounted) is more than 160 million EURO.

Alternatively, you can also monetize (attributable) health impacts from a non-healthiar source.

results <- monetize(
  impact = 1151,
  valuation = 100
)

Function call

cba <- cba(
    output_attribute = results_pm_copd,
    valuation = 50000,
    cost = 100000000,
    discount_shape = "exponential",
    discount_rate_benefit = 0.03,
    discount_rate_cost = 0.03,
    n_years_benefit = 5,
    n_years_cost = 5
)

Main results

cba$cba_main |>  
  dplyr::select(benefit, cost, net_benefit) |> 
  knitr::kable()

We see that the central and upper 95% confidence interval estimates of avoided attributable COPD cases result in a net monetary benefit of the policy, while the lower 95% confidence interval estimate results in a net cost!

Function call

results_pm_yll <- attribute_lifetable(
  year_of_analysis = 2019, 
  health_outcome = "yll",
  rr_central =  1.118, 
  rr_increment = 10,
  erf_shape = "log_linear",
  exp_central = 8.85,
  cutoff_central = 5,
  min_age = 20, # age from which population is affected by the exposure
  # Life table information
  age_group = exdat_lifetable$age_group,
  sex = exdat_lifetable$sex,
  population = exdat_lifetable$midyear_population,
  # In the life table case, BHD refers to deaths
  bhd_central = exdat_lifetable$deaths
) 

Main results

Total YLL attributable to exposure (sum of sex-specific impacts)

Use the two arguments approach_exposure and approach_newborns to modify the life table calculation

• approach_exposure

  • "single_year" (default) population is exposed for a single year and the impacts of that exposure are calculated

  • "constant" population is exposed every year

• approach_newborns

  • "without_newborns" (default) assumes that the population in the year of analysis is followed over time, without considering newborns being born

  • "with_newborns" assumes that for each year after the year of analysis n babies are born, with n being equal to the (male and female) population aged 0 that is provided in the argument population

Detailed results

Attributable YLL results

• per year

• per age (group)

• per sex (if sex-specific life table data entered)

are available

Note: we will inspect the results for females; male results are also available

Results per year

NOTE: only a selection of years is shown

results_pm_yll$health_detailed$results_raw |>
  dplyr::summarize(
    .by = year, 
    impact = sum(impact, na.rm = TRUE)
  )
#> # A tibble: 100 × 2
#>    year  impact
#>    <chr>  <dbl>
#>  1 2019   1300.
#>  2 2020   2422.
#>  3 2021   2221.
#>  4 2022   2033.
#>  5 2023   1858.
#>  6 2024   1695.
#>  7 2025   1545.
#>  8 2026   1409.
#>  9 2027   1284.
#> 10 2028   1171.
#> # ℹ 90 more rows

results_pm_yll$health_detailed$results_raw |>
  dplyr::summarize(
    .by = year, 
    impact = sum(impact, na.rm = TRUE)) |>
  knitr::kable()

Results by age

Results by year and age

NOTE: only a selection of years and ages is shown

Function call

results_pm_deaths <- attribute_lifetable(
  health_outcome = "deaths",
  year_of_analysis = 2019,
  rr_central =  1.118, 
  rr_increment = 10,
  erf_shape = "log_linear",
  exp_central = 8.85,
  cutoff_central = 5,
  min_age = 20, # age from which population is affected by the exposure   
  # Life table information
  age_group = exdat_lifetable$age_group,   
  sex = exdat_lifetable$sex,
  population = exdat_lifetable$midyear_population, 
  bhd_central = exdat_lifetable$deaths
)

Main results

Total premature deaths attributable to exposure (sum of sex-specific impacts)

Detailed results

Attributable premature deaths results

• per year (if argument approach_exposure = "constant")

• per age (group)

• per sex (if sex-specific life table data entered)

are available

Note: we will inspect the results for females; male results are also available

Note: because we set the function argument approach_exposure = "constant" in the function call results are available for one year (the year of analysis)

Function call

results_pm_copd_yld  <- attribute_health(
  rr_central = 1.1, 
  rr_increment = 10, 
  erf_shape = "log_linear",  
  exp_central = 8.85,
  cutoff_central = 5,
  bhd_central = 1000,
  duration_central = 10,
  dw_central = 0.2
)

Main results

Function call

results_daly <- daly(
     output_attribute_yll = results_pm_yll,
     output_attribute_yld = results_pm_copd_yld
)

Main results

YLL, YLD & DALY

Function call

results_pm_copd <- attribute_health(
  erf_shape = "log_linear",
  rr_central = 1.369, 
  rr_increment = 10,
  exp_central = 8.85,
  cutoff_central = 5,
  bhd_central = 30747
) 

results_no2_copd <- attribute_mod(
  output_attribute = results_pm_copd,
  exp_central = 10.9,
  rr_central = 1.031
)

results_multiplicative <- multiexpose(
  output_attribute_exp_1 = results_pm_copd,
  output_attribute_exp_2 = results_no2_copd,
  exp_name_1 = "pm2.5",
  exp_name_2 = "no2",
  approach_multiexposure = "multiplicative"
)

Main results

results_multiplicative$health_main

Function call

mdi <- prepare_mdi(
  geo_id_micro = exdat_prepare_mdi$id,
  edu = exdat_prepare_mdi$edu,
  unemployed = exdat_prepare_mdi$unemployed,
  single_parent = exdat_prepare_mdi$single_parent,
  pop_change = exdat_prepare_mdi$pop_change,
  no_heating = exdat_prepare_mdi$no_heating,
  n_quantile = 10,
  verbose = FALSE
)

Note: verbose = FALSE suppresses any output to the console (default: verbose = TRUE, i.e. with printing turned on).

Main results

Function output includes

• mdi_main, a tibble containing the BEST-COST MDI

mdi$mdi_main |> 
  select(geo_id_micro, MDI, MDI_index)

Detailed results

• mdi_detailed

  • DESCRIPTIVE STATISTICS

  • PEARSON’S CORRELATION COEFFICIENTS

  • CRONBACH’S α, including the reliability rating this value indicates

  • Code for boxplots of the single indicators

  • Code for histogram of the MDI’s for the geo units with a normal distribution curve

To reproduce the boxlots run

eval(mdi$mdi_detailed$boxplot)

Analogeously, to reproduce the histogram run

eval(mdi$mdi_detailed$histogram)