Methodological refresher

Figure: Relative risk approach

Function call

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)
)

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.

Refresher - Burden of disease with absolute risk

Figure: Absolute risk approach

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
)

results_noise_ha$health_detailed$results_raw

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"
)

using relative risk

Goal: e.g. To quantify the disease cases attributable to PM2.5 exposure in multiple cities using one single command.

• Enter unique ID’s as a vector (numeric or character) to the geo_id_micro argument (e.g. municipality names or province abbrevations)

• Optional: aggregate unit-specific results by providing higher-level ID’s (e.g. region names or country abbreviations) as a vector (numeric or character) to the geo_id_macro argument

Input to the other function arguments is specified as usual, either as a vector or a single values (which will be recycled to match the length of the other input vectors).

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)
)

using absolute risk

Goal: E.g.To quantify high annoyance cases attributable to noise exposure in rural and urban areas.

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"
)

Confidence interval

Goal: E.g. To quantify the COPD cases attributable to PM2.5 exposure taking into account uncertainty (lower and upper bound of confidence interval) in several input arguments: relative risk, exposure and baseline health data.

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
) 

Monte Carlo simulation

Goal: E.g. To summarize uncertainty of attributable health impacts (i.e. to get a single confidence interval instead of many combinations) by using a Monte Carlo simulation.

You can do this carrying out a Monte Carlo uncertainty analysis via the summarize_uncertainty() function.

The outcome of the Monte Carlo analysis is added to the variable entered as the results argument, which is results_pm_copd in our case.

Two lists (“folders”) are added:

• uncertainty_main contains the central estimate and the corresponding 95% confidence intervals obtained through the Monte Carlo assessment

• uncertainty_detailed contains all n_sim simulations of the Monte Carlo assessment

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

by age group

Goal: e.g. To quantify health impacts attributable to air pollution in a country by age group.

To obtain age-group-specific results, the baseline health data (and possibly exposure) must be available by age group.

If the age argument was specified, age-group-specific results are available under health_detailed in the sub-folder results_by_age_group.

results_age_group <- attribute_health(
        approach_risk = "relative_risk",
        age = c("below_65", "65_plus"),
        exp_central = c(8, 7),
        cutoff_central = c(5, 5),
        bhd_central = c(1000, 5000),
        rr_central = 1.06,
        rr_increment = 10,
        erf_shape = "log_linear"
      )
results_age_group$health_detailed$results_by_age_group |> 
  dplyr::select(age_group, impact_rounded, exp, bhd) |> 
  knitr::kable()

by sex

Goal: e.g. To quantify health impacts attributable to air pollution in a country by sex.

To obtain sex-specific results, the baseline health data (and possibly exposure) must be entered by sex.

If the sex argument was specified, sex-specific results are available under health_detailed in the sub-folder results_by_sex.

results_sex <- attribute_health(
        approach_risk = "relative_risk",
        sex = c("female", "male"),
        exp_central = c(8, 8),
        cutoff_central = c(5, 5),
        bhd_central = c(1000, 1100),
        rr_central = 1.06,
        rr_increment = 10,
        erf_shape = "log_linear"
      )
results_sex$health_detailed$results_by_sex |> 
  dplyr::select(sex, impact_rounded, exp, bhd) |> 
  knitr::kable()

by other sub-groups

Goal: e.g. To quantify health impacts attributable to air pollution in a country by education level.

A (random) ID vector of unique values must be entered to the argument geo_id_micro, to indicate that each row is an observation on its own (and not just exposure categories)

A single vector (or a data frame / tibble with multiple columns) to group the results by can be entered to the info argument. In this example, this will be information about the education level.

In a second step one can group the results based on one or more columns and so summarize the results by the preferred sub-groups

info <- data.frame(
  education = rep(c("secondary", "bachelor", "master"), each = 5) # education level
)
output_attribute <- attribute_health(
  rr_central = 1.063,
  rr_increment = 10,
  erf_shape = "log_linear",
  cutoff_central =  0,
  exp_central = sample(6:10, 15, replace = TRUE),
  bhd_central = sample(100:500, 15, replace = TRUE),
  geo_id_micro = c(1:nrow(info)), # (random) ID must be entered
  info = info
)
output_stratified <- output_attribute$health_detailed$results_raw |>
  dplyr::group_by(info_column_1) |>
  dplyr::summarize(mean_impact = mean(impact)) |>
  print()
#> # A tibble: 3 × 2
#>   info_column_1 mean_impact
#>   <chr>               <dbl>
#> 1 bachelor             18.2
#> 2 master               17.1
#> 3 secondary            11.1

YLL

Goal: E.g. To quantify the years of life lost (YLL) due to deaths from COPD attributable to PM2.5 exposure during one year.

We can use attribute_lifetable() combined with life table input data to determine YLL attributable to an environmental stressor.

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
) 

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()

Deaths

Goal: E.g. To determine premature deaths from COPD attributable to PM2.5 exposure during one year.

See example on YLL for additional info on attribute_lifetable() calculations and its output.

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
)

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

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 <- 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 <- attribute_mod(
  output_attribute = results_pm,
  exp_central = 10.9,
  rr_central = 1.031
)

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

Main results

results_multiplicative$health_main

Monetization

Goal: E.g. To monetize the attributable health impact of a policy that will have health benefits five years from now.

The outcome of the monetization is added to the variable entered to the output_attribute argument, which is results_pm_copd in our case.

Two folders are added:

• monetization_main contains the central monetization estimate and the corresponding 95% confidence intervals obtained through the specified monetization

• monetization_detailed contains the monetized results for each unique combination of the input variable estimates that were provided to the initial attribute_health() call

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

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

Cost-benefit analysis

Goal: E.g. To estimate the net-benefit of a health policy via cost-benefit analysis (CBA).

Let’s imagine we design a policy that would reduce air pollution to 5 \(\mu g/m^3\), which is the concentration specified in the cutoff_central argument in the initial attribute_health() call. So we could avoid all COPD cases attributed to air pollution.

Considering the cost to implement the policy (estimated at 100 million EURO), what would be the monetary net benefit of such a policy, ? We can find out using healthiar’s cba() function.

The outcome of the CBA is contained in two folders, which are added to the existing assessment:

• cba_main contains the central estimate and the corresponding 95% confidence intervals obtained

• cba_detailed contains additional intermediate results for both cost and benefit

  • benefit contains results by_year and raw results health_raw

  • cost contains the costs of the policy at the end of the period specified in the n_years_cost argument

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
)

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

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)

Export results

Export as .csv file

Save as .Rdata file

Export to Excel (as .xlsx file)

Visualize results

Visualization is out of scope of healthiar. You can visualize in

• R, e.g. with the ggplot2 package (online book by the creator)
• Excel (export results first)
• Other tools