Input

You can enter data in healthiar functions using: - hard coded values or - columns inside pre-loaded data frames or tibbles.

Let’s see some examples calling the most important function in healthiar: attribute_health().

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

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

Goal

E.g., to quantify the COPD cases attributable to PM2.5 (air pollution) exposure in a country.

Methodology

The comparative risk assessment approach (C. J. Murray et al. 2003) is applied obtaining the population attributable fraction (percent of cases that are attributable to the exposure) based on the relative risk. The exposure scenario is compared with a counter-factual scenario.

This approach has been extensive documented and applied (e.g., WHO 2003; Steenland and Armstrong 2006; Soares et al. 2022; Pozzer et al. 2023; GBD 2019 Risk Factors Collaborators 2020; Lehtomäki et al. 2025).

Figure: Relative risk approach

Function call

results_pm_copd <- attribute_health(
  approach_risk = "relative_risk", # If you do not call this argument, "relative_risk" will be assigned by default.
  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
)

Main results

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 of 3 rows and 23 columns. Be aware that this main output contains input data, some intermediate steps and the final results in different formats.

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)

Goal

E.g., to quantify the number incidence cases of high annoyance attributable to (road traffic) noise exposure.

Methodology

In the absolute risk calculation pathway, estimates are based on the size and distribution of the exposed population, rather than on baseline health data, as is the case in the relative risk pathway (WHO 2011).

Figure: Absolute risk approach

\[N = \sum AR_i \times PE_i\]

Where:

• \(N\) = attributed cases

• \(AR_i\) = absolute risk at category \(i\)

• \(PE_i\) = absolute population exposed at category \(i\)

healthiar consistently handles cut-off values across both relative and absolute risk approaches. When a cutoff_ argument is specified alongside an exposure–response function (erf_eq_), the function evaluates exposure as excess exposure relative to the cut-off (c as \(\text{exp} - \text{cutoff}\)), effectively shifting the function by the cut-off. However, many absolute risk curves published in the literature are parameterized on raw exposure levels that already incorporate the cut-off directly into the function definition. Therefore, healthiar issues an informative warning if both erf_eq_ and a cutoff_ are specified when using attribute_health(). If your absolute risk function is parameterized on raw exposure levels rather than on the exposure above the cut-off, see chapter Shifted vs. unshifted exposure-response functions.

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 on user-defined ERF).

Main results

Results per noise exposure band

results_noise_ha$health_detailed$results_raw

Remember, that if the equation of the exposure-response function (erf_eq_...) requires taking a maximum in a vectorised context, pmax() must be used instead of max(). pmax() should be used whenever an element-wise maximum is required (the output will be a vector), while max() returns a single global maximum for the entire vector. For example:

erf_eq_central <- 
  "exp(0.2969*log((pmax(0,c-2.4)/1.9+1))/(1+exp(-(pmax(0,c-2.4)-12)/40.2)))"  

One exposure category

Alternatively, it’s also possible to only assess the absolute risk impacts for one exposure category (e.g., 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"
)

using relative risk

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

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

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

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

Goal

E.g., to quantify COPD cases attributable to air pollution exposure by applying a user-defined exposure-response function (ERF), such as the MR-BRT curves from Global Burden of Disease study.

Function call

In this case, the function arguments erf_eq_... require a function as input, so we use an auxiliary function (splinefun()) to transform the points on the ERF into type function.

results_pm_copd_mr_brt <- attribute_health(
  exp_central = 8.85,
  bhd_central = 30747,
  cutoff_central = 0,
  # Specify the function based on x-y point pairs that lie on the ERF
  erf_eq_central = splinefun(
    x = c(0, 5, 10, 15, 20, 25, 30, 50, 70, 90, 110),
    y = c(1.00, 1.04, 1.08, 1.12, 1.16, 1.20, 1.23, 1.35, 1.45, 1.53, 1.60),
    method = "natural")
)

The ERF curve created looks as follows

Alternatively, other functions (e.g. approxfun()) can be used to create the ERF

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

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

output_attribute <- healthiar::attribute_health(
    rr_central = 1.063,
    rr_increment = 10,
    erf_shape = "log_linear",
    cutoff_central =  0,
    exp_central = c(6, 7, 8,
                    7, 8, 9,
                    8, 9, 10,
                    9, 10, 11),
    bhd_central = c(600, 700, 800,
                    700, 800, 900,
                    800, 900, 1000,
                    900, 1000, 1100),
    geo_id_micro = rep(c("a", "b", "c", "d"), each = 3),
    info = data.frame(
      education = rep(c("secondary", "bachelor", "master"), times = 4)) # education level
  )

output_stratified <- output_attribute$health_detailed$results_raw |>
      dplyr::group_by(info_education) |>
      dplyr::summarize(mean_impact = mean(impact))|>
      dplyr::pull(mean_impact) |>
      print()
#> [1] 43.72087 54.26844 34.30332

by age, sex and other sub-groups

output_attribute <- healthiar::attribute_health(
    rr_central = 1.063,
    rr_increment = 10,
    erf_shape = "log_linear",
    cutoff_central =  0,
    age_group = rep(c("50_and_younger", "50_plus"), each = 4, times= 2),
    sex = rep(c("female", "male"), each = 2, times = 4),
    exp_central = c(6, 7, 8, 7, 8, 9, 8, 9,
                    10, 9, 10, 11, 10, 11, 12, 13),
    bhd_central = c(600, 700, 800, 700, 800, 900, 800, 900,
                    1000, 900, 1000, 1100, 1000, 1100, 1200, 1000),
    geo_id_micro = rep(c("a", "b"), each = 8),
    info = data.frame(
      education = rep(c("without_master", "with_master"), times = 8)) # education level
  )

output_stratified <- output_attribute$health_detailed$results_raw |>
      dplyr::group_by(info_education) |>
      dplyr::summarize(mean_impact = mean(impact))|>
      dplyr::pull(mean_impact) |>
      print()
#> [1] 52.80090 49.83826

Goal

E.g., to quantify in one single call the impacts of several exposure-outcome pairs, such as PM2.5 and COPD together with NO2 and asthma, or to compare the same pair quantified with the relative risks of two different studies.

Parallel analyses, not sub-groups

As explained in the chapter Sub-group analysis, the impacts of the sub-groups are summed by default, so only the total per geographic unit is shown in health_main.

That is not what you want here. Premature deaths attributable to PM2.5 and asthma cases attributable to NO2 are not two parts of one population: they are two analyses of the same population, so adding them up counts the same people twice. The same applies to a sensitivity analysis in which one pair is quantified with the relative risks of two different studies: the two results are alternatives to each other, not summands.

Enter therefore in main_results_by the names of the dimensions that the main results have to be reported by. They stay as separate rows in health_main and in all results_by_... tables. Note that main_results_by does not create those tables, which are available in the detailed output in any case: it determines which dimensions survive in them and in the main results.

Function call

The pairs are identified in a column of info, and that column is named in main_results_by. All other arguments are entered as parallel vectors, exactly as for multiple geographic units. Note that the pairs can have different exposure-response functions, here log-linear for PM2.5 and linear for NO2.

output_attribute <- healthiar::attribute_health(
    info = data.frame(
      pair = rep(c("pm2.5_copd", "no2_asthma"), each = 2)),
    main_results_by = "pair",
    exp_central    = c(8.85, 9.20, 22.1, 24.5),
    cutoff_central = c(5, 5, 10, 10),
    rr_central     = c(1.369, 1.369, 1.041, 1.041),
    rr_increment   = 10,
    erf_shape      = c("log_linear", "log_linear", "linear", "linear"),
    bhd_central    = c(30747, 31500, 12000, 12500),
    geo_id_micro   = rep(c("a", "b"), times = 2)
  )

Main results

One row per exposure-outcome pair and geographic unit. Without main_results_by the four impacts would have been added up into one single number per geographic unit.

output_attribute$health_main |>
  dplyr::select(info_pair, geo_id_micro, erf_ci, impact_rounded) |>
  print()
#> # A tibble: 4 × 4
#>   info_pair  geo_id_micro erf_ci  impact_rounded
#>   <chr>      <chr>        <chr>            <dbl>
#> 1 pm2.5_copd a            central           3502
#> 2 pm2.5_copd b            central           3893
#> 3 no2_asthma a            central            567
#> 4 no2_asthma b            central            701

Sensitivity analysis with several exposure-response functions

The same mechanism compares one exposure-outcome pair quantified with the relative risks of two different studies.

output_sensitivity <- healthiar::attribute_health(
    info = data.frame(rr_source = c("Liu 2020", "WHO 2021")),
    main_results_by = "rr_source",
    exp_central = c(8.85, 8.85),
    cutoff_central = 5,
    rr_central = c(1.369, 1.080),
    rr_increment = 10,
    erf_shape = "log_linear",
    bhd_central = c(30747, 30747)
  )

output_sensitivity$health_main |>
  dplyr::select(info_rr_source, erf_ci, impact_rounded) |>
  print()
#> # A tibble: 2 × 3
#>   info_rr_source erf_ci  impact_rounded
#>   <chr>          <chr>            <dbl>
#> 1 Liu 2020       central           3502
#> 2 WHO 2021       central            898

Which dimensions can be entered

main_results_by is available in both attribute_health() and attribute_lifetable(). It accepts:

• the columns of info, using the names that you gave them (e.g. "pair"), or "info" itself if you entered a vector instead of a data frame;
• the arguments sex, age_group, exp_category, geo_id_micro and geo_id_macro;
• in attribute_lifetable() also year, i.e. the years of the time horizon, which are otherwise added up.

For example, main_results_by = "sex" shows one row per sex in health_main instead of the total across sexes. This changes nothing in the calculation, only which rows are kept apart in the results, and the rows still add up to the total that you would get without the argument.

A name that is not available in the assessment produces an error listing the options, so a typo in a column name of info cannot pass unnoticed.

Limitations

All pairs in one call must use the same way of defining the exposure-response function, i.e. either rr_... or erf_eq_... for all of them, and the same approach_risk.

standardize() keeps the pairs apart as well, so each pair gets its own age-standardized impacts. Nothing has to be done for that: its results simply carry the columns of main_results_by.

socialize() and summarize_uncertainty() analyse one pair at a time and stop with an error otherwise, so for them attribute_health() has to be called once per pair. The reason is different in each case. socialize() groups the impacts only by geographic unit, age group and social quantile, and its arguments have no way of identifying the pairs when the input is entered directly instead of as an attribute_health() output. summarize_uncertainty() simulates one value per variable and geographic unit, assuming that e.g. the relative risk is the same everywhere, which is not the case when each pair has its own.

Data preparation

The life table approach to obtain YLL and deaths requires population and baseline mortality data to be stratified by one year age groups. However, in some cases these data are only available for larger age groups (e.g., 5-year data: 0-4 years old, 5-9 years old, …). What to do?

• If your population and mortality data are not available by one-year age group, data must be prepared by interpolating values. The healthiar function prepare_lifetable() makes this conversion using the same approach as the WHO tool AirQ+ (WHO 2020). In standard AirQ+ life table disaggregations, deaths are assumed to be distributed uniformly across age intervals with a fraction of the age interval lived of 0.5. The argument fraction_lived allows other assumptions, see the section on the fraction of the age interval lived below.

• If your population and death data are stratified by one-year age group, you are lucky, you can ignore this initial step. prepare_lifetable() passes single-year data on unchanged, so it can also be used for a life table that mixes already available single-year data with converted larger age groups.

age_groups <- c(0, 5, 10, 15)
pop <- c(438200, 445100, 439800, 421500)
bhd_counts <- c(1420, 45, 50, 125)

# Standard uniform disaggregation (fraction_lived = 0.5)
prepared_data <- healthiar::prepare_lifetable(
  age_group = age_groups,
  population = pop,
  bhd = bhd_counts
  # fraction_lived = 0.5 (default)
)

# Age-group specific fraction of the age interval lived (first fraction_lived = 0.1)
prepared_data_with_age_specific_fraction_lived <- healthiar::prepare_lifetable(
  age_group = age_groups,
  population = pop,
  bhd = bhd_counts,
  fraction_lived = c(0.1, 0.5, 0.5, 0.5)
)

YLL

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

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
)

Goal

E.g., to quantify the years lived with disability (YLD) attributable to air pollution exposure using disability weights.

Methodology

To quantify the YLDs, you can use a prevalence-based or an incidence-based approach (Kim et al. 2022).

• Prevalence-based : Enter 1 (year) in the argument(s) dw_... and prevalence cases in bhd_....

• Incidence-based: Enter a value above 1 in dw_... and incidence cases in bhd_....

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

Goal (e.g.)

E.g., to obtain the disability-adjusted life years (DALY) as the sum of YLLs and YLDs.

Methodology

To obtain the attributable DALY, its two components, i.e. years of life lost (YLL) and years lived with disability (YLD), must be summed (GBD 2019 Risk Factors Collaborators 2020).

\[DALY = YLL + YLD\]

Function call

This is possible using the function daly().

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

Main results

YLL, YLD & DALY

Goal

E.g., to quantify health impacts using attribute_healthin an scenario B very similar to a previous scenario A.

Function call

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)

The function attribute_mod() can be used to modify one or multiple arguments of attribute_healthin an existing scenario, e.g. scenario_A.

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

This is equivalent to building the whole scenario again (see below), but more time and code efficient.

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)

Goal

E.g., to compare the health impacts in the scenario “before intervention” vs. “after intervention”.

Methodology

Two approaches can be used for the comparison of scenarios:

• Delta: Subtraction of health impact in scenario 1 minus in scenarios 2 (i.e. two PAF) (WHO Regional Office for Europe 2014)

• Population impact fraction (PIF) (WHO 2003; C. J. L. Murray et al. 2003; Askari and Namayandeh 2020).

Note that the PIF comparison approach assumes same baseline health data for scenario 1 and 2 (e.g., comparison of two scenarios at the same time point), while the delta comparison approach, the difference between two scenarios is obtained by subtraction. Therefore, the delta approach is suited for comparison of scenarios in different time points.

IMPORTANT If your aim is to quantify health impacts from a policy intervention, be aware that you should use the same year of analysis and therefore same health baseline data in both scenarios. The only variable that should change in the second scenario is the exposure (change as a result of the intervention).

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_mod(
  output_attribute = scenario_A, 
  exp_central = 6
)

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

results_comparison <- healthiar::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)

Goal

E.g., to quantify the total health impact attributable to PM2.5 and NO2.

Methodology

A methodological report of the EU project BEST-COST (Strak et al. 2024) identified three approaches to add up attributable health impacts from correlated exposures:

• Additive approach (Steenland and Armstrong 2006):

\[PAF_{additive} = PAF_{exposure1} + PAF_{exposure2}\]

• Multiplicative approach (Jerrett et al. 2013):

\[PAF_{multiplicative} = \frac{\sum PE \times (rr\_at\_exp_{multiplicative} - 1)}{\sum PE \times (rr\_at\_exp_{multiplicative}-1) + 1}\]

\[rr\_at\_exp_{multiplicative} = rr\_at\_exp_{exposure1} * rr\_at\_exp_{exposure2}\]

• Combined approach (Steenland and Armstrong 2006):

\[PAF_{combined} = 1-[(1-PAF_{exposure1}) \times (1-PAF_{exposure2})]\]

Attention: To apply any of these approaches, the relative risks for one exposure must be adjusted for the second exposure and the way round.

Function call

For this purpose, you can use the function multiexpose().

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

Goal

E.g., to obtain the age-standardized attributable health impacts of two age groups.

Methodology

Age standardization involves adjusting the observed rates of a particular outcome to a standard population with a specific age structure. This is a technique used to allow the comparison of populations with different age structures (GBD 2019 Demographics Collaborators 2020; Ahmad et al. 2001). In healthiar, the function standardize() applies the direct method, where the age-specific rates observed in a study population are applied to a standard (reference) population distribution.

The standardized health impact rate is computed as \[ impact\_per\_100k\_inhab_{std} = \sum_{i=1}^{k} (impact\_per\_100k\_inhab_i \times ref\_prop\_pop_i) \]

where:

• \(impact\_per\_100k\_inhab_{std}\) is the age-standardized health impact rate.
• \(impact\_per\_100k\_inhab_i\) is the impact rate observed in age group \(i\) (e.g., impact per 100,000 inhabitants).
• \(ref\_prop\_pop_i\) is the proportion of the reference population in age group \(i\) .
• \(k\) is the number of age groups.

Function call

output_attribute <- attribute_health(
  rr_central = 1.063,
  rr_increment = 10,
  erf_shape = "log_linear",
  cutoff_central =  0,
  age_group = c("below_40", "above_40"),
  exp_central = c(8.1, 10.9),
  bhd_central = c(1000, 4000),
  population = c(100000, 500000)
  )

results <- standardize(
  output_attribute = output_attribute,
  age_group = c("below_40", "above_40"),
  ref_prop_pop = c(0.5, 0.5)
  )

Main results

Age-standardized impact rate:

print(results$health_main$impact_per_100k_inhab)  
#> [1] 49.91113

Age group-specific impact rate:

print(results$health_detailed$results_raw$impact_per_100k_inhab)  
#> [1] 48.28250 51.53977

Goal

E.g., to determine population-weighted mean PM2.5 exposure for several neighborhoods of Brussels (Belgium)

Methodology

The healthiarfunction prepare_exposure() helps users that do not have the exposure data (needed for healthiar functions), but only spatial concentration and population data. The function calculates an average concentration value in each geographic unit, weighted with population at each location.

\[ exp = \frac{\sum_{i=1}^{n} (C_i \times population_i)}{\sum_{i=1}^{n} population_i} \]

where:

• \(exp\) = population-weighted mean exposure for the geographic unit.
• \(C_i\) = pollutant concentration in grid cell \(i\).
• \(population_i\) = population count in grid cell \(i\).
• \(n\) = total number of grid cells contained by the geographic unit.

In case population is entered as count by geographic sub-unit, the function calculates the mean concentration in each sub-unit and aggregates it to higher-level geographic units. If no population data is entered, the function calculates a simple spatial mean concentration as exposure value.

The output of prepare_exposure() can be entered in the argument exp_mean, exp_lower and/or exp_upper in healthiar functions such as attribute_health().

Function call

# exdat_pwm_1 = Pollution grid data
exdat_pwm_1 <- terra::rast(system.file("extdata", "exdat_pwm_1.tif", package = "healthiar"))

# exdat_pwm_2 = Data with the geo units and population data. This is pre-loaded in healthiar.
# If your raw data are in .gpkg format, you can use e.g.  sf::st_read() 

pwm <- healthiar::prepare_exposure(
  poll_grid = exdat_pwm_1, # Formal class SpatRaster,
  geo_units = exdat_pwm_2, # sf of the geographic sub-units
  population = sf::st_drop_geometry(exdat_pwm_2$population), # population per geographic sub-unit
  geo_id_macro = sf::st_drop_geometry(exdat_pwm_2$region)) # higher-level IDs to aggregate

Main results

Within the function output, the list exposure_main contains the population-weighted mean exposures for the (higher-level) geographic units in the column exposure_mean and the total population in each unit in column population_total.

Goal

E.g., to quantify health impacts with an exposure-response function whose effect threshold differs from the counterfactual exposure of the assessment, such as noise assessments in which exposure data are only available above 55 dB although health effects already occur at lower exposures.

Methodology

healthiar distinguishes two exposure levels:

• threshold is the effect threshold, i.e. the exposure level from which the exposure-response function starts to show an effect. It is the anchor of the curve and is therefore subtracted from the exposure. The term is used in this sense in health risk assessment guidance, e.g., for environmental noise (Engelmann et al. 2025).
• cutoff_... is the exposure level below which no health impacts are quantified. As long as it is identical to the effect threshold, it is also the counterfactual exposure of the assessment. If it is higher, it only determines which exposures are assessed, while the effect threshold remains the counterfactual exposure. The term is used in this sense in health risk assessment guidance, e.g., for air pollution (WHO 2025). Be aware that a cut-off recommended in guidance documents may reflect the exposure range covered by the evidence rather than a demonstrated no-effect level, i.e. it does not necessarily indicate that there are no health effects below it.

Both levels are usually identical. Therefore it is enough to enter only one of them, and the other one takes the same value in the background. Accordingly, there are three possible combinations:

1. Cut-off equal to the effect threshold

This is the default situation. The exposure-response function is anchored at that level, which is therefore subtracted from the exposure, and the exposures below it get the risk at the reference level:

\[ rr_{at\_exp} = f(\max(exp, cutoff) - cutoff) \]

The same result is obtained by entering only cutoff_..., only threshold, or both with the same value.

2. Cut-off higher than the effect threshold

The exposure-response function is truncated at the cut-off. The exposures above the cut-off are assessed with the risk referring to the effect threshold, while the exposures below the cut-off are treated as unexposed:

\[ rr_{at\_exp} = \begin{cases} 1, & exp < cutoff \\ f(exp - threshold), & exp \ge cutoff \end{cases} \]

In the absolute risk approach, the exposures below the cut-off get an absolute risk of 0 instead of a relative risk of 1, i.e. in both approaches no health impacts are attributed to them.

This is the case, e.g., in noise assessments under the Environmental Noise Directive, in which exposure data are only reported above a given exposure level (cut-off) although the exposure-response functions show effects already at lower exposures (effect threshold). Be aware that this results in a step in the exposure-response function at the cut-off and that the people exposed between the effect threshold and the cut-off get no attributable health impacts (i.e. a conservative estimate).

3. Cut-off lower than the effect threshold

The cut-off has no effect on the results and healthiar shows a warning. The reason is that the exposure-response function is anchored at the effect threshold, so the exposures below the effect threshold already get the risk at the reference level. A cut-off below the effect threshold cannot remove any additional exposure and the results are the same as if only threshold had been entered.

Function call

The example below refers to case 2 above: health impacts in the exposure group 55 dB+ (cut-off) that are affected by an exposure above the effect threshold of 45 dB.

results_cutoff_above_threshold <- healthiar::attribute_health(
  approach_risk = "relative_risk",
  erf_shape = "log_linear",
  rr_central = 1.055,
  rr_increment = 10,
  prop_pop_exp = c(300000,200000,150000,120000,100000,70000,60000)/10000000,
  exp_central = c(47,52,57,62,67,72,77),
  # Effect threshold, i.e. the anchor of the exposure-response function
  threshold = 45,
  # Cut-off, i.e. the exposures below it are not quantified
  cutoff_central = 55,
  bhd_central = 50000)$health_main$impact_rounded

Alternatively, the same result can be obtained by entering a user-defined exposure-response function in erf_eq_... which already includes the case distinction. This is only needed if the shape of the curve above the cut-off is not covered by the erf_shape options. The function arguments erf_eq_... require a function as input. Instead of using a splinefun() this can also be fulfilled by using a ‘function(c)’ which is of type ‘function’.

#setting up function parameters
threshold <- 45
RR <- 1.055
cutoff <- 55
rr_increment <- 10

# define categorical function, the ifelse condition enables the case distinction
erf_function <- function(c){
  output <- ifelse(c<cutoff, 1, exp((log(RR)/rr_increment)*(c-threshold)))
  return(output)
}
# attribute_health
results_catERF_different_cutoff <- healthiar::attribute_health(
  approach_risk = "relative_risk",
  erf_eq_central = erf_function,
  prop_pop_exp = c(300000,200000,150000,120000,100000,70000,60000)/10000000,
  exp_central = c(47,52,57,62,67,72,77),
  cutoff_central=0,
  bhd_central=50000)$health_main$impact_rounded

The categorical ERF curve created looks as follows. The step at the cut-off is clearly visible.

Goal

E.g., to use an exposure-response function published in the literature which is parameterized on raw exposure levels instead of on the exposure above a cut-off.

Methodology

healthiar evaluates the exposure-response function entered in erf_eq_... at \(c = exp - threshold\), and the threshold takes the value of the cut-off if it is not entered (see chapter Cut-off vs. threshold). Therefore, the cut-off moves the exposure-response function horizontally:

1. Shifted exposure-response function. If a cut-off is entered, the exposure-response function is moved by the cut-off, i.e. \(c = exp - cutoff\). The function must therefore be parameterized on the exposure above the cut-off.

2. Unshifted exposure-response function. If no cut-off is entered, nothing is subtracted from the exposure, i.e. \(c = exp\). The function must therefore be parameterized on raw exposure levels.

Be aware that healthiar shows a warning if both erf_eq_... and cutoff_... are entered, because many exposure-response functions published in the literature already incorporate the cut-off in the function definition.

Function call

# Case A: Shifted exposure-response function
# The function is parameterized on the exposure above the cut-off
shifted_erf <- healthiar::attribute_health(
  approach_risk = "absolute_risk",
  exp_central = c(50, 60, 70),
  pop_exp = c(300000,200000,150000),
  cutoff_central = 53,
  erf_eq_central = "78.9270 - 3.1162 * (c + 53) + 0.0342 * (c + 53)^2"
)
#> Warning: You entered a value for: cutoff_central alongside absolute risk.
#> Be aware that healthiar shifts the exposure in 'erf_eq' as c = (exp - cutoff).

# Case B: Unshifted exposure-response function
# The function is parameterized on raw exposure levels
unshifted_erf <- healthiar::attribute_health(
  approach_risk = "absolute_risk",
  exp_central = c(50, 60, 70),
  pop_exp = c(300000,200000,150000),
  erf_eq_central = "78.9270 - 3.1162 * c + 0.0342 * c^2"
)

If you want the exposure-response function to stay unshifted but the exposures below a given level not to be quantified, enter a cutoff_... value higher than the threshold (see chapter Cut-off vs. threshold).

Monetization

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

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

Health impact attributable to social indicator

 health_impact <- healthiar::attribute_health(
   age_group = exdat_socialize$age_group,
   exp_central = exdat_socialize$pm25_mean,
   cutoff_central = 0,
   rr_central = exdat_socialize$rr,
   erf_shape = "log_linear",
   rr_increment = 10,
   bhd_central = exdat_socialize$mortality,
   population = exdat_socialize$population,
   geo_id_micro = exdat_socialize$geo_unit)

social_t <- healthiar::socialize(
  output_attribute = health_impact,
  age_group = exdat_socialize$age_group, # They have to be the same in socialize() and in attribute_health()
  ref_prop_pop = exdat_socialize$ref_prop_pop, # Population already provided in output_attribute
  geo_id_micro = exdat_socialize$geo_unit,
  social_indicator = exdat_socialize$score,
  n_quantile = 10,
  increasing_deprivation = TRUE)
#> Warning: 194 geographic unit(s) have no value in social_indicator. They get no
#> social_quantile and are therefore not included in the comparison between the
#> quantiles (first and last), but they are still part of the overall values.

social <- healthiar::socialize(
  impact = health_impact$health_detailed$results_by_age_group$impact,
  age_group = exdat_socialize$age_group, # They have to be the same in socialize() and in attribute_health()
  ref_prop_pop = exdat_socialize$ref_prop_pop,
  geo_id_micro = exdat_socialize$geo_unit,
  social_indicator = exdat_socialize$score,
  population = exdat_socialize$population, # Population has to be provided because no output_attribute
  n_quantile = 10,
  increasing_deprivation = TRUE)
#> Warning: 194 geographic unit(s) have no value in social_indicator. They get no
#> social_quantile and are therefore not included in the comparison between the
#> quantiles (first and last), but they are still part of the overall values.

#> # A tibble: 4 × 5
#>   parameter      difference_type difference_compared_…¹ difference_value comment
#>   <chr>          <chr>           <chr>                             <dbl> <chr>  
#> 1 impact_rate_s… absolute        last_quantile                    26.5   <NA>   
#> 2 impact_rate_s… relative        last_quantile                     0.597 <NA>   
#> 3 impact_rate_s… absolute        overall                          14.2   It can…
#> 4 impact_rate_s… relative        overall                           0.242 It can…
#> # ℹ abbreviated name: ¹​difference_compared_with

Multiple deprivation index

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
)

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

eval(mdi$mdi_detailed$boxplot)

eval(mdi$mdi_detailed$histogram)

Pipe |>

healthiar can be used inside the native pipes |>. See the example below.

exdat_noise |>
  (\(df) {
    healthiar::attribute_health(
      approach_risk = df$risk_estimate_type,
      exp_central = df$exposure_mean,
      pop_exp = df$exposed,
      erf_eq_central = df$erf
      )$health_main$impact_rounded
    })()

Shorter making used of the base R function with().

exdat_noise |>
      (\(df) {
        with(df, healthiar::attribute_health(
         approach_risk = risk_estimate_type,
         exp_central = exposure_mean,
         pop_exp = exposed,
         erf_eq_central = erf
         ))$health_main$impact_rounded
        })()
#> [1] 348464

Pipe %>%

healthiar can also be used inside magrittr pipes %>% as follows.

exdat_noise %>%
  {
    healthiar::attribute_health(
      approach_risk = .$risk_estimate_type,
      exp_central = .$exposure_mean,
      pop_exp = .$exposed,
      erf_eq_central = .$erf
    )$health_main$impact_rounded
  }

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 using base programming or packages such as ggplot2 (Wickham 2016),
• Excel (export results first) or
• Other tools.