--- title: "Vignette: Intro to healthiar" output: rmarkdown::html_vignette: toc: true toc_depth: 2 bibliography: ../inst/REFERENCES.bib vignette: > %\VignetteIndexEntry{Vignette: Intro to healthiar} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo=TRUE, eval=TRUE, include=TRUE) ``` ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) options(rmarkdown.html_vignette.check_title = FALSE) options(knitr.kable.max_rows = 10) set.seed(1) ## Avoid using pacman here, as it causes error in installation if it's not installed already library(healthiar) library(tibble) library(dplyr) library(purrr) library(tidyr) library(knitr) ## terra, sf and exactextractr are only suggested by healthiar. The chunks of ## the section on prepare_exposure() are skipped when they are not installed, ## for instance on the CRAN check machines that install no suggested package spatial_packages_available <- requireNamespace("terra", quietly = TRUE) && requireNamespace("sf", quietly = TRUE) && requireNamespace("exactextractr", quietly = TRUE) ``` Hi there! This vignette will tell you about `healthiar` and show you how to use `healthiar` with the help of examples. *Note*: Before using `healthiar`, please read carefully the information provided in the [readme file](https://github.com/SwissTPH/healthiar?tab=readme-ov-file#readme) or the [welcome webpage](https://swisstph.github.io/healthiar/). By using `healthiar`, you agree to the [terms of use and disclaimer](https://github.com/SwissTPH/healthiar?tab=readme-ov-file#readme). *Tip*: Do you want to be notified when a new `healthiar` version is released? See the [Stay updated](https://github.com/SwissTPH/healthiar?tab=readme-ov-file#stay-updated) section of the readme file. ------------------------------------------------------------------------ # About `healthiar` The `healthiar` functions allow you to quantify and monetize the health impacts (or burden of disease) attributable to exposure. The main focus of the EU project that initiated the development of `healthiar` (BEST-COST) has been two environmental exposures: air pollution and noise. However, `healthiar` could be used for other exposures such as green spaces, chemicals, physical activity... See below a an overview of the `healthiar`, which is the first page of the [cheat sheet](https://swisstph.github.io/healthiar/articles/cheatsheet.html). The whole list of functions included in `healthiar` is linked there and available in the [reference](https://swisstph.github.io/healthiar/reference/index.html). ![Figure: Overview of `healthiar`](../man/figures/cheatsheet_healthiar_1st_page.png) # Input & output data ## 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()`. ### Hard coded vs. columns #### Hard coded Depending on the function argument, you will need to enter numeric or character values. ```{r eval=FALSE, include=TRUE} 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 ) ``` #### Columns `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 let's `attribute_health()` with input data from the `healthiar` example data. Note that you can easily provide input data to the function argument using the `$` operator. ```{r, eval=FALSE, include=TRUE} 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 ) ``` ### Tidy data Be aware that `healthiar` functions are easier to use if your data is prepared in a tidy format, i.e.: - Each variable is a column; each column is a variable. - Each observation is a row; each row is an observation. - Each value is a cell; each cell is a single value. To know more about the concept of tidy format, see the article by [@Wickham2014_jss]. For example, in `attribute health()` the length of the input vectors to be entered in the arguments must be either 1 or the result of the combinations of the different values of: - `geo_id_micro` - `exp_...` - `sex` - `age` - (`info` for further sub-group analysis) ## Output #### Structure The output of the `healthiar`function `attribute_health()` and `attribute_lifetable` consists of two lists ("folders"): - `health_main` contains the main results - `health_detailed` contained detailed results and additional info about the assessment. In other `healthiar` functions you can find a similar output structure but using different prefixes. E.g., `social_`in `socialize()` and `monetization_`in `monetitize()`. #### Access A similar structure can be found in other large functions in `helathiar`, e.g., `attribute_lifetable()`, `compare()`, `socialize()` or `monetize()`. In some functions, different elements are available in the output. For instance, `attribute_lifetable()` creates additional output that is specific to life table calculations. 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")` --------------------------------------------------------------------------- # Function examples The descriptions of the [`healthiar` functions](https://swisstph.github.io/healthiar/reference/index.html) provide examples that you can execute (with `healthiar` loaded) by running `example("function_name")`, e.g. `example("attribute_health")`. In the sections below in this vignette, you find additional examples and more detailed explanations. # Relative risk ### Goal E.g., to quantify the COPD cases attributable to PM2.5 (air pollution) exposure in a country. ### Methodology The comparative risk assessment approach [@Murray2003_e] 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 counterfactual scenario. This approach has been extensive documented and applied [e.g., @WHO2003_report; @Steenland2006-e; @Soares2020_report; @Pozzer2023_gh; @GBD2020_tl; @Lehtomaki_2025_eh]. ![Figure: Relative risk approach](../man/figures/bod_rr.png){width="700"} #### Population attributable fraction General integral form for the **population attributable fraction (PAF)**: $$PAF = \frac{\int rr\_at\_exp(x) \times PE(x)dx - 1}{\int rr\_at\_exp(x) \times pop\_exp(x)dx}$$ Where: * $x$ = exposure level * $PE(x)$ = population distribution of exposure * $rr\_at\exp(x)$ = relative risk at exposure level compared to reference ##### Simplified for categorical exposure distribution If exposure is categorical, the integrals are converted to sums: $$PAF = \frac{\sum rr\_at\_exp_i \times PE_i - 1}{\sum rr\_at\_exp_i \times PE_i}$$ Alternatively, an equivalent form is: $$PAF = \frac{\sum PE_i \times (rr\_at\_exp_i - 1)}{\sum PE_i\times (rr\_at\_exp_i - 1) + 1}$$ `healthiar` uses this second form, which is the one published as Formula 3 in [ETC HE Report 2023/11](https://www.eionet.europa.eu/etcs/etc-he/products/etc-he-products/etc-he-reports/etc-he-report-2023-11-environmental-noise-health-risk-assessment-methodology-for-assessing-health-risks-using-data-reported-under-the-environmental-noise-directive). Note that $PE_i$ is the proportion of the **total** population in the exposure category $i$, so the $PE_i$ do **not** have to add up to 1: the part of the population that they do not cover is treated as unexposed, i.e. it gets the relative risk of the reference level. In air pollution assessments everybody is usually exposed and the $PE_i$ add up to 1, whereas in noise assessments they typically add up to less than 1, because exposure is only reported above a given level. Accordingly, the baseline health data entered in `bhd_...` must always refer to the total population. The first form above divides by $\sum rr\_at\_exp_i \times PE_i$ instead and is therefore only equivalent when $\sum PE_i = 1$. ##### Simplified for single exposure value If there is one single single exposure value, corresponding to the population weighted mean concentration, the equation can be simplified as follows: $$PAF = \frac{rr\_at\_exp - 1}{rr\_at\_exp }$$ #### Scaling relative risk How to get this relative risk at exposure level (`rr_at_exp`)? This is normally different to the relative risk published in the epidemiological literature (`rr`) together with the (concentration/dose) `increment` that corresponds to this relative risk. The equations used for scaling relative risk depend on the chosen exposure-response function shapes: - linear [@Lehtomaki_2025_eh] $$RRexp = 1 + \frac{rr - 1}{increment} \times (exp - cutoff)$$ - log-linear [@Lehtomaki_2025_eh] $$RRexp = e^{\frac{\log(rr)}{increment} \times (exp - cutoff)}$$ - log-log [@Lehtomaki_2025_eh] $$RRexp = \left( \frac{exp + 1}{cutoff + 1} \right)^{\frac{\log(rr)}{\log(increment + cutoff + 1) - \log(cutoff + 1)}}$$ - linear-log [@Pozzer2023_gh] $$RRexp = 1 + \frac{\log(rr - 1)}{\log(increment + cutoff + 1) - \log(cutoff + 1)} \times \frac{\log(exp + 1)}{\log(cutoff + 1)}$$ The relative risk at exposure level (`rr_at_exp`) and is part of the output of `attribute_health()` and `attribute_lifetable()`. `rr_at_exp` can also be calculated using `get_risk()`. For conversion of hazard ratios and/or odds ratios to relative risks refer to [@VanderWeele2019_biom] and/or use the conversion tools developed by the Teaching group in EBM in 2022 for hazard ratios (https://ebm-helper.cn/en/Conv/HR_RR.html) and/or odds ratios (https://ebm-helper.cn/en/Conv/OR_RR.html). ### Function call ```{r, eval=TRUE, include=TRUE, echo=TRUE} 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 ```{r, include=TRUE, eval=TRUE, echo=TRUE} results_pm_copd$health_main ``` 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. ```{r} 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) # Absolute risk ### 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 [@WHO2011_report]. ![Figure: Absolute risk approach](../man/figures/bod_ar.png){width="700"} $$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 [Shifting the ERF](#shifting-the-erf). ### Function call ```{r} 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 ```{r echo=FALSE} results_noise_ha |> purrr::pluck("health_main") |> dplyr::select(erf_eq, erf_ci, impact_rounded) |> knitr::kable() # Prints tibble in a minimal layout ``` ### Results per noise exposure band ```{r eval=FALSE, include=TRUE, echo=TRUE} results_noise_ha$health_detailed$results_raw ``` ```{r echo=FALSE, include=TRUE, eval=TRUE} results_noise_ha[["health_detailed"]][["results_raw"]] |> select(exp_category, exp, pop_exp, impact) |> knitr::kable() ``` 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: ```{r eval=FALSE, include=TRUE, echo=TRUE} 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). ```{r} 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" ) ``` ```{r eval=FALSE, include=FALSE, echo=TRUE} results_noise_ha$health_main ``` ```{r eval=TRUE, include=TRUE, echo=FALSE} results_noise_ha[["health_main"]] |> dplyr::select(exp_category, impact) |> knitr::kable() ``` ```{r eval=TRUE, include=FALSE, echo=FALSE} rm(results_noise_ha) ``` # Multiple geographic units ## using relative risk ### Goal E.g., to quantify the disease cases attributable to PM2.5 exposure in multiple cities using one single command. ### Function call - Enter unique ID's as a vector (`numeric` or `character`) to the `geo_id_micro` argument (e.g., municipality names or province abbreviations) - 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). ```{r} 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 ```{r eval=FALSE, include=FALSE, echo=TRUE} results_iteration$health_main ``` ```{r eval=TRUE, include=TRUE, echo=FALSE} results_iteration[["health_main"]] |> dplyr::select(geo_id_macro, impact_rounded, erf_ci, exp_ci, bhd_ci) |> knitr::kable() ``` 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 `r sum(results_iteration[["health_main"]]$impact_rounded)` (using a relative risk of 1.369). ### Detailed results The geo unit-specific information and results are stored under `health_detailed`\>`results_raw` . ```{r eval=FALSE, include=FALSE, echo=TRUE} results_iteration$health_detailed$results_raw ``` ```{r eval=TRUE, include=TRUE, echo=FALSE} results_iteration[["health_detailed"]][["results_raw"]] |> dplyr::select(geo_id_micro, impact_rounded, geo_id_macro) |> knitr::kable() ``` `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). ```{r eval=TRUE, include=FALSE, echo=FALSE} rm(results_iteration) ``` ## using absolute risk ### Goal E.g., to quantify high annoyance cases attributable to noise exposure in rural and urban areas. ### Function call ```{r} data <- exdat_noise |> ## Filter for urban and rural regions dplyr::filter(region == "urban" | region == "rural") ``` ```{r} 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). ```{r eval=FALSE, include=FALSE, echo=TRUE} results_iteration_ar$health_main ``` ```{r eval=TRUE, include=TRUE, echo=FALSE} results_iteration_ar[["health_main"]] |> dplyr::select(geo_id_macro, impact_rounded, erf_ci, exp_ci) |> knitr::kable() ``` ### Detailed results Impact by geo unit, in this case impact in the rural and in the urban area. ```{r eval=FALSE, include=FALSE, echo=TRUE} results_iteration_ar$health_detailed$results_by_geo_id_micro ``` ```{r eval=TRUE, include=TRUE, echo=FALSE} results_iteration_ar[["health_detailed"]][["results_by_geo_id_micro"]] |> dplyr::select(geo_id_micro, geo_id_macro, impact) |> knitr::kable() ``` ```{r eval=TRUE, include=FALSE, echo=FALSE} rm(results_iteration_ar, data) ``` # Uncertainty ## 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. ### Function call ```{r} 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: ```{r eval=FALSE, echo=TRUE, include=FALSE} results_pm_copd$health_detailed$results_raw ``` ```{r echo=FALSE} results_pm_copd$health_detailed$results_raw |> dplyr::select(erf_ci, exp_ci, bhd_ci, impact_rounded) |> dplyr::slice(1:9) |> knitr::kable() # Prints tibble in a minimal layout ``` 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 `r length(1:9)` of the `r length(results_pm_copd$health_detailed$results_raw$impact)` possible combinations are displayed due to space constraints. *Note*: only a selection of columns are shown. ## 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. ### Methodology #### General concepts A Monte Carlo simulation is a statistical method that generates repeated random sampling [@Robert2004_book; @Rubinstein2016_book]. In `healthiar`, you can use the function `summarize_uncertainty()` to simulate values in the arguments with uncertainty and estimate a single confidence interval in the results. For each entered input argument that includes a (95%) confidence interval (i.e. `_lower` and `_upper` bound value) a distribution is fitted (see distributions below). The median of the simulated attributable impacts is reported as the central estimate. The 2.5th and 97.5th percentiles of these simulated impacts define the lower and upper bounds of the 95% summary uncertainty interval. Aggregated estimates (e.g. for a `geo_id_macro` covering several `geo_id_micro`) are obtained by first summing the impacts of all lower level units within each simulation and only then taking the median and the 2.5th and 97.5th percentiles of these sums. Summing instead the central, lower and upper estimates of each lower level unit would assume that the uncertainty is perfectly correlated across the units and would therefore overestimate the width of the confidence interval. #### Distributions used for simulation `summarize_uncertainty()` assumes the following shapes of the distributions in the simulations: - Relative risk: The values are simulated based on an optimized *gamma* distribution, which fits well as relative risks are positive and their distributions are usually right-skewed. The gamma distribution is parametrized such that its mean is equal to the central relative risk estimate (`rate= shape/rr_central`). The shape parameter is then optimized using `stats::optimize()` to match the inputed 95% confidence interval bounds, with `stats::qgamma()` used to evaluate candidate distributions. Finally, `n_sim` relative risk values are simulated using `stats::rgamma()`. - Exposure, cutoff, baseline health data and duration: The values are simulated based on a *normal* distribution using `stats::rnorm()` with `mean = exp_central`, `mean = cutoff_central`, `mean = bhd_central` or `mean = duration_central` and a standard deviation based on corresponding lower and upper 95% exposure confidence interval values. The standard deviation is calculated as $$(upper-lower)/(2*1.96)$$, since for a normal distribution the 95% CI spans approximately two standard deviations on either side of the mean. As these four variables cannot be negative, the distribution is truncated at zero, i.e. negative values are discarded and drawn again until they are positive. Please note that the normal distribution is symmetric. Therefore, the simulated values reproduce the entered confidence interval only if the entered `_lower` and `_upper` values are symmetric around the `_central` value. Only the width of the entered confidence interval is used, not its shape. The more asymmetric the entered confidence interval, the more the simulated values will depart from it. Moreover, if the entered confidence interval is so wide that a relevant share of the distribution falls below zero, the truncation at zero shifts the simulated central estimate above the entered `_central` value. - Disability weights: The values are simulated based on a *beta* distribution, as both the disability weights and the beta distribution are bounded by 0 and 1. The beta distribution best fitting the inputted central disability weight estimate and corresponding lower and upper 95% confidence interval values is fitted using `stats::qbeta()` (the best fitting distribution parameters `shape1` and `shape2` are determined using `stats::optimize()`). For this purpose, we partly adapted the R function `prevalence::beta_expert` with permission of one of the authors [@Devleesschauwer2022_package]. Finally, `n_sim` disability weight values are simulated using `stats::rbeta()`. For stability of the 95% confidence interval, a large number of simulations (e.g., 10,000) is recommended in practice. The example below uses n_sim = 100 for brevity. ### Function call ```{r} results_pm_copd_summarized <- summarize_uncertainty( output_attribute = results_pm_copd, n_sim = 100 ) ``` ### Main results 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 and - `uncertainty_detailed` contains all `n_sim` simulations of the Monte Carlo assessment. ```{r eval=FALSE, include=FALSE, echo=TRUE} results_pm_copd_summarized$uncertainty_main ``` ```{r eval=TRUE, include=TRUE, echo=FALSE} results_pm_copd_summarized$uncertainty_main |> knitr::kable() ``` ### Detailed results The folder `uncertainty_detailed` contains all single simulations. Let's look at the impact of the first 10 simulations. The columns `erf_ci`, `exp_ci`, `bhd_ci`, and `cutoff_ci` indicate the source of uncertainty component used for that simulation (in the first 10 simulations, all use central estimates). ```{r eval=FALSE, include=FALSE, echo=TRUE} results_pm_copd_summarized$uncertainty_detailed$impact_by_sim ``` ```{r eval=TRUE, include=TRUE, echo=FALSE} results_pm_copd_summarized$uncertainty_detailed$impact_by_sim |> knitr::kable() ``` ```{r, echo=FALSE, eval=TRUE, include=FALSE} rm(results_pm_copd_summarized) ``` # User-defined ERF ### 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`. ```{r} 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 ```{r fig.alt="ERF curve", eval=TRUE, include=TRUE, echo=FALSE} 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) spline_func <- stats::splinefun( x = x, y = y, method = "natural") ## Generate finer x-values x_dense <- seq(min(x), max(x), length.out = 200) ## Evaluate the spline function at these points y_dense <- spline_func(x_dense) ## Plot plot(x, y, # main = "User-defined ERF", main = "", xlab = "Exposure", ylab = "Relative risk", col = "blue", pch = 19) lines(x_dense, y_dense, col = "red", lwd = 2) legend("topleft", legend = c("Original points", "Spline curve"), col = c("blue", "red"), pch = c(19, NA), lty = c(NA, 1), lwd = c(NA, 2)) ``` Alternatively, other functions (e.g. `approxfun()`) can be used to create the ERF # Sub-group analysis Every sub-group stratifies the assessment: by age group, by sex, or by anything entered in `info`, such as education level. What differs is what happens to that stratification in the main results. For the sub-groups in this chapter the impacts are **summed**, so the stratification disappears from `health_main` and only the total per geographic unit is shown. It is still available in the `results_by_...` tables of the detailed output, and that is usually what you want: the total across education levels, or across age groups, is a meaningful number. Not every dimension of an assessment is a sub-group in this sense, though. A stratification splits one population into parts that add up to the whole, which is exactly why summing them is meaningful. Different exposure-outcome pairs do not: PM2.5 and mortality and NO2 and asthma cover the *same* whole population, so they are parallel analyses rather than parts of one. Adding them up counts the same people twice. Those dimensions therefore have to be **kept** in the main results, which is what the argument `main_results_by` does. See the chapter *Multiple exposure-outcome pairs*. ## by age group ### Goal E.g., to quantify health impacts attributable to air pollution in a country *by age group*. ### Function call 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`. ```{r} 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 by age group ```{r} results_age_group$health_detailed$results_by_age_group |> dplyr::select(age_group, impact_rounded, exp, bhd) |> knitr::kable() ``` ```{r eval=TRUE, include=FALSE, echo=FALSE} rm(results_age_group) ``` ## by sex ### Goal E.g., to quantify health impacts attributable to air pollution in a country *by sex*. ### Function call The baseline health data (and possibly exposure) must be entered by sex. ```{r eval=TRUE, include=TRUE, echo=TRUE} 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 by sex If the `sex` argument was specified, sex-specific results are available under `health_detailed` in the sub-folder `results_by_sex`. ```{r eval=TRUE, include=TRUE, echo=TRUE} results_sex$health_detailed$results_by_sex |> dplyr::select(sex, impact_rounded, exp, bhd) |> knitr::kable() ``` ```{r eval=TRUE, include=FALSE, echo=FALSE} rm(results_sex) ``` ## by other sub-groups ### Goal E.g., to quantify attributable health impacts *stratified by a sub-group different to age and sex, e.g., education level*. ### Function call 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. Be aware that, if the argument `info`is a data frame or tibble, the columns in the `info` argument are renamed adding the prefix `info_`. If the argument `info` is a vector, the column in the results is called `info`. ```{r eval=TRUE, include=TRUE, echo=TRUE} 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 ) ``` ### Results by other sub-group ```{r eval=TRUE, include=TRUE, echo=TRUE} output_stratified <- output_attribute$health_detailed$results_raw |> dplyr::group_by(info_education) |> dplyr::summarize(mean_impact = mean(impact))|> dplyr::pull(mean_impact) |> print() ``` ```{r eval=TRUE, include=FALSE, echo=FALSE} rm(output_attribute, output_stratified) ``` ## by age, sex and other sub-groups ### Goal E.g., to quantify attributable health impacts *stratified by age, sex and additional sub-group e.g. education level*. ### Function call ```{r eval=TRUE, include=TRUE, echo=TRUE} 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 ) ``` ### Results by all sub-groups ```{r eval=TRUE, include=TRUE, echo=TRUE} output_stratified <- output_attribute$health_detailed$results_raw |> dplyr::group_by(info_education) |> dplyr::summarize(mean_impact = mean(impact))|> dplyr::pull(mean_impact) |> print() ``` ```{r eval=TRUE, include=FALSE, echo=FALSE} rm(output_attribute, output_stratified) ``` # Multiple exposure-outcome pairs ### 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. ```{r eval=TRUE, include=TRUE, echo=TRUE} 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. ```{r eval=TRUE, include=TRUE, echo=TRUE} output_attribute$health_main |> dplyr::select(info_pair, geo_id_micro, erf_ci, impact_rounded) |> print() ``` ### Sensitivity analysis with several exposure-response functions The same mechanism compares one exposure-outcome pair quantified with the relative risks of two different studies. ```{r eval=TRUE, include=TRUE, echo=TRUE} 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() ``` ### 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. ```{r eval=TRUE, include=FALSE, echo=FALSE} rm(output_attribute, output_sensitivity) ``` # YLL & deaths with life table ## All-cause mortality only The life table approach of `attribute_lifetable()` is only applicable to **all-cause mortality**. The baseline deaths (`bhd_central`) are used for two purposes at the same time: to project how the population survives over the years and to obtain the deaths that are attributable to the exposure. The projection of the population needs the deaths of all causes, so `bhd_central` must contain all-cause deaths and the exposure-response function must refer to all-cause mortality. Cause-specific mortality (e.g. lung cancer) cannot be assessed with `attribute_lifetable()`: - If only the cause-specific deaths are entered, the projected population hardly dies of any other cause and reaches unrealistically high ages. The YLL and the attributable deaths are then overestimated. - If all-cause deaths are entered with a cause-specific exposure-response function, the relative risk is applied to all causes of death, which overestimates the impacts as well. For cause-specific mortality, the attributable deaths can be obtained with `attribute_health()`. ## 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+ [@WHO2020_report]. 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. ```{r Examples using prepare_lifetable()} 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) ) ``` ### Fraction of the age interval lived The argument `fraction_lived`, available in both `prepare_lifetable()` and `attribute_lifetable()`, is the average fraction of an age interval that is lived by those who die in it. It corresponds to $_na_x$ in the life table literature [@Chiang1984_book; @Preston2001_book]. The default of 0.5 for all age groups means that deaths are evenly distributed over the age interval, which is the value that AirQ+ assumes [@WHO2020_report]. It can be entered as a single value applying to all age groups or as one value per age group. #### What it determines With $_nm_x$ the death rate of an age interval of width $n$, the probability of dying in that interval is $$_nq_x = \frac{n \times {}_nm_x}{1 + n \times (1 - {}_na_x) \times {}_nm_x}$$ and the person-years lived at a single year of age, i.e. its mid-year population, is one year for each person reaching the next birthday plus $fraction\_lived$ years for each person dying at that age: $$midyear\_population = entry\_population - (1 - fraction\_lived) \times deaths$$ In `prepare_lifetable()` these two relations fix the conversion. The entry population of the first single year of age of an age group follows from the condition that the conversion may neither create nor lose population, i.e. that the single-year mid-year populations have to add up to the mid-year population of the age group that was entered: $$entry\_population_{first} = \frac{population + (1 - fraction\_lived) \times bhd}{\sum_{k = 0}^{n - 1} prob\_survival^{\,k}}$$ For $fraction\_lived = 0.5$ this reproduces the values published in the AirQ+ manual, and unlike the formula given there it holds for any age interval width and any value of `fraction_lived`. Both relations are checked in internal tests in `healthiar`against a published life table, the abridged life table for Austrian males in 1992 [@Preston2001_book], in which the fraction of the age interval lived ranges from 0.48 to 0.63 in the 5-year age groups and is 0.068 at age 0. Applying the probability of dying obtained by `prepare_lifetable()` to the published survivors reproduces the published deaths of that life table for every age group, and the entry population it determines from the published mid-year population and deaths at age 0 reproduces the published number of survivors. #### What it does not determine `fraction_lived` is *not* a way to control how the deaths of an age group are distributed over the single years of age it contains. That distribution follows from the survival probabilities, as in AirQ+, so `fraction_lived` influences it only indirectly and far too weakly to be used for it. Both quantities are distinct in the life table literature and the distribution of the deaths cannot be derived from $_na_x$ [@Preston2001_book]. `prepare_lifetable()` can therefore not reproduce the concentration of deaths at age 0 that is observed in low-mortality settings. If the deaths at age 0 matter for your assessment, obtain single-year population and death data for the youngest ages instead of relying on the conversion. #### How to choose it The sensitivity of the survival probability to `fraction_lived` is $$\frac{\partial \; prob\_survival}{\partial \; fraction\_lived} = - \frac{hazard\_rate^2}{(1 + (1 - fraction\_lived) \times hazard\_rate)^2}$$ i.e. it grows with the square of the hazard rate. `fraction_lived` therefore matters little in the age groups with a low mortality and most in those with the highest mortality, which are the last ones. With the example data of this vignette, entering 0.1 for age 0 changes the YLL by 0.01%, whereas adapting the last age group changes them by 0.25%. Values of $_na_x$ other than 0.5 are usually derived from the mortality level itself. For age 0 and for the ages 1 to 4, @Preston2001_book report rules of thumb going back to Coale and Demeny, e.g. $_1a_0 = 0.045 + 2.684 \times {}_1m_0$ for males with $_1m_0 < 0.107$. For an open-ended last age group, $_na_x$ is the remaining life expectancy $1 / {}_nm_x$, which can exceed the width of the age interval. Whichever `fraction_lived` is used in `prepare_lifetable()` has to be passed on to `attribute_lifetable()`, which needs the same assumption to obtain the survival probabilities. The output column `fraction_lived_for_attribute` is provided for that purpose. ## 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. ### Methodology #### General concept The life table methodology of `attribute_lifetable()` follows that of the WHO tool AirQ+ [@WHO2020_report], which is described in more detail by @Miller2003_jech. The generalized formulas, which do not assume a fixed fraction of life lived (fraction_lived = 0.5), are based on the foundational work by @Chiang1984_book. In short, two scenarios are compared: a) a scenario with the exposure level specified in the function ("exposed scenario") and b) a scenario with no exposure ("unexposed scenario"). First, the entry and mid-year populations of the (first) year of analysis in the unexposed scenario are determined using modified survival probabilities. Second, age-specific population projections using scenario-specific survival probabilities are done for both scenarios. Third, by subtracting the populations in the unexposed scenario from the populations in the exposed scenario the premature deaths/years of life lost attributable to the exposure are determined. An expansive life table case study for is available in a report by @Miller2010_report. Use the two arguments `approach_exposure` and `approach_newborns` to modify the life table calculation: - `approach_exposure` - `"single_year"` (default): Population is exposed for only one year (year of analysis) and attributable health impacts reflect the single-year snapshot. Applicable only for deaths. AirQ+ approach for deaths in 2025. - `"constant"` : Population is exposed every year, i.e. air pollution exposure is sustained across the full projection horizon and attributable health impacts are projected and cumulated across the full lifetime period. Applicable for both deaths and YLL. AirQ+ approach for YLL in 2025. - `approach_newborns` - `"without_newborns"` (default): Population in the year of analysis is followed over time, without considering newborns being born. - `"with_newborns"`: 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. #### Last age group The life table is closed at the last age group, i.e. the survivors of the last age group are not projected into a further age. This follows the usual practice in the life table method, which extends the life table up to an age past which the probability of survival is negligible and then sets the survival above it to zero [@Miller2003_jech]. The data must therefore reach an age at which survival is negligible, i.e. the last age group should be one in which essentially all remaining deaths occur. National life tables usually already close in this way, with a probability of dying of 1 in the last age group. If the last age group still has an appreciable number of survivors, those survivors are dropped from the projection: neither the life years they would still live, nor the difference in those life years between the exposed and the unexposed scenario, are counted. The attributable health impacts are then underestimated. With the example data used in this vignette, which end at age 99 and still have a survival probability of about 0.34 there, the underestimation is around 0.7%. If that is the case, condensate the last age group, i.e. sum the populations and the deaths of the highest ages into it, so that all remaining deaths fall within it. Do *not* add extra age groups beyond the data instead: the added ages would need a population greater than 0, that population would be considered exposed and would generate attributable health impacts of its own, and the results would then grow with the number of added age groups instead of converging. #### Comparison with AirQ+ `healthiar` aligns with AirQ+ in the equations of the life table. The survival probability of AirQ+, $S_i = (2 - h_i) / (2 + h_i)$, is identical to the one used here, and the modified survival probability is obtained in both tools by dividing the hazard rate by the relative risk (which is equivalent to multiplying it by $1 - PAF$) [@WHO2020_report]. `healthiar` reproduces the results of AirQ+. For the same input data, the attributable YLL and premature deaths of both tools agree to the last decimal that AirQ+ prints, for an exposure of one single year as well as for a constant exposure, year by year and accumulated over several years, and so do the life expectancies with and without the exposure. The comparison has to be made by sex, though, because the "All genders" results of AirQ+ are not the sum of its male and female results as in `healthiar`: AirQ+ obtains them from a life table in which the populations and the deaths of both sexes are pooled. With the data of the example below, the pooled life table for all genders in AirQ+ gives 1.6% more YLL than the sum of two sex-specific ones in `healthiar`. `healthiar` returns that pooled value too, when the populations and the deaths of both sexes are entered as one single life table, i.e. passing the same value in `sex` for every age group. Be aware when comparing with AirQ+ that its results accumulated over several years (e.g. "over 10 years") correspond to `approach_exposure = "constant"` and not to `approach_exposure = "single_year"`. #### Determination of populations in the (first) year of analysis ##### Entry population The entry (i.e. start of year) populations in both scenarios (exposed and unexposed) is determined as follows: $$entry\_population_{year_1} = midyear\_population_{year_1} + (1 - fraction\_lived) \times deaths_{year_1}$$ ##### Survival probabilities ###### Exposed scenario The survival probabilities in the exposed scenario from start of year $i$ to start of year $i+1$ are calculated as follows: $$prob\_survival = \frac{midyear\_population_i - fraction\_lived \times deaths_i}{midyear\_population_i + (1 - fraction\_lived) \times deaths_i}$$ Analogously, the probability of survival from start of year $i$ to mid-year $i$: $$prob\_survival\_until\_midyear = 1 - (1 - fraction\_lived) \times (1 - prob\_survival)$$ ###### Unexposed scenario The survival probabilities in the unexposed scenario are calculated as follows: First, the age-group specific hazard rate in the exposed scenario is calculated using the inputted age-specific mid-year populations and deaths. $$hazard\_rate = \frac{deaths}{mid\_year\_population}$$ Second, the hazard rate is multiplied by the modification factor ($= 1 - PAF$) to obtain the age-specific hazard rate in the unexposed scenario. $$hazard\_rate\_mod = hazard\_rate \times modification\_factor$$ Third, the age-specific survival probabilities (from the start until the end in a given age group) in the unexposed scenario are calculated as follows: $$prob\_survival\_mod = \frac{1 - fraction\_lived \times hazard\_rate\_mod}{1 + (1 - fraction\_lived) \times hazard\_rate\_mod}$$ ##### Mid-year population The mid-year population of the (first) year of analysis (year_1) in the unexposed scenario are determined as follows: First, the survival probabilities from start of year $i$ to mid-year $i$ in the unexposed scenario is calculated as: $$prob\_survival\_until\_midyear_{mod} = 1 - (1 - fraction\_lived) \times (1 - prob\_survival\_mod)$$ Second, the mid-year populations of the (first) year of analysis (year_1) in the unexposed scenario is calculated: $$midyear\_population\_unexposed_{year_1} = entry\_population_{year_1} \times prob\_survival\_until\_midyear_{mod}$$ #### Population projection Using the age group-specific and scenario-specific survival probabilities calculated above, future populations of each age-group under each scenario are calculated. ###### Exposed scenario The population projections for the two possible options of `approach_exposure` (`"single_year"` and `"constant"`) for the unexposed scenario are different. In the case of `"single_year"` exposure, the population projection for the years after the year of exposure is the same as in the unexposed scenario. In the case of `"constant"` the population projection is done as follows: First, the entry population of year $i+1$ is calculated (which is the same as the end of year population of year $i$) using the entry population of year $i$. $$entry\_population_{i+1} = entry\_population_i \times prob\_survival$$ Second, the mid-year population of year $i+1$ is calculated. $$midyear\_population_{i+1} = entry\_population_{i+1} \times prob\_survival\_until\_midyear$$ ###### Unexposed scenario The entry and mid-year population projections in the unexposed scenario are done as follows: First, the entry population of year $i+1$ is calculated (which is the same as the end of year population of year $i$) by multiplying the entry population of year $i$ and the modified survival probabilities. $$entry\_population_{i+1} = entry\_population_i \times prob\_survival\_mod$$ Second, the mid-year population of year $i+1$ is calculated. $$midyear\_population_{i+1} = entry\_population_{i+1} \times prob\_survival\_until\_midyear_{mod}$$ ### Function call We can use `attribute_lifetable()` combined with life table input data to determine YLL attributable to an environmental stressor. ```{r} 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). ```{r eval=FALSE, include=FALSE, echo=TRUE} results_pm_yll$health_main ``` ```{r eval=TRUE, include=TRUE, echo=FALSE} results_pm_yll$health_main |> dplyr::slice_head() |> dplyr::select(impact_rounded, erf_ci, exp_ci, bhd_ci) |> knitr::kable() ``` ### 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. ```{r} results_pm_yll$health_detailed$results_raw |> dplyr::summarize( .by = year, impact = sum(impact, na.rm = TRUE) ) ``` ```{r} results_pm_yll$health_detailed$results_raw |> dplyr::summarize( .by = year, impact = sum(impact, na.rm = TRUE)) |> knitr::kable() ``` #### YLL ```{r eval=FALSE, include=FALSE, echo=TRUE} results_pm_yll$health_detailed$intermediate_calculations |> dplyr::filter(sex == "female") |> dplyr::pull(impact_by_age_and_year) |> purrr::pluck(1) ``` ```{r echo=FALSE} results_pm_yll$health_detailed$intermediate_calculations |> dplyr::filter(sex == "female") |> dplyr::pull(impact_by_age_and_year) |> purrr::pluck(1) |> dplyr::select(age_start, age_end, impact_2019) |> dplyr::slice( ( dplyr::n() - 8 ) : dplyr::n() ) |> knitr::kable() ``` #### Population (baseline scenario) Baseline scenario refers to the scenario with exposure (i.e. the scenario specified in the assessment). ```{r eval=FALSE, include=FALSE, echo=TRUE} results_pm_yll$health_detailed$intermediate_calculations|> dplyr::filter(sex == "female") |> dplyr::pull(projection_if_exposed_by_age_and_year) |> purrr::pluck(1) ``` ```{r echo=FALSE} results_pm_yll$health_detailed$intermediate_calculations|> dplyr::filter(sex == "female") |> dplyr::pull(projection_if_exposed_by_age_and_year) |> purrr::pluck(1) |> dplyr::select(age_start, midyear_population_2019, midyear_population_2020, midyear_population_2021, midyear_population_2022) |> dplyr::slice( ( dplyr::n() - 8 ) : dplyr::n() ) |> knitr::kable() ``` #### Population (unexposed scenario) Impacted scenario refers to the scenario without exposure. ```{r eval=FALSE, include=FALSE, echo=TRUE} results_pm_yll$health_detailed$intermediate_calculations|> dplyr::filter(sex == "female") |> dplyr::pull(projection_if_unexposed_by_age_and_year) |> purrr::pluck(1) ``` ```{r echo=FALSE} results_pm_yll$health_detailed$intermediate_calculations|> dplyr::filter(sex == "female") |> dplyr::pull(projection_if_unexposed_by_age_and_year) |> purrr::pluck(1) |> dplyr::select(age_start, midyear_population_2019, midyear_population_2020, midyear_population_2021, midyear_population_2022) |> dplyr::slice( ( dplyr::n() - 8 ) : dplyr::n() ) |> knitr::kable() ``` ## Deaths ### Goal (e.g.) E.g., to determine premature deaths from COPD attributable to PM2.5 exposure during one year. ### Function call See example on YLL for additional info on `attribute_lifetable()` calculations and its output. ```{r} 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). ```{r eval=FALSE, include=FALSE, echo=TRUE} results_pm_deaths$health_main$impact ``` ```{r, eval=TRUE, include=TRUE, echo=FALSE} results_pm_deaths$health_main |> dplyr::slice_head() |> dplyr::select(impact_rounded, erf_ci, exp_ci, bhd_ci) |> knitr::kable() ``` ### 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 inspect below the results for females; male results are also available. *Note*: Because we set the function argument `approach_exposure = "constant"` results are projected and accumulated across the full time_horizon. If you set approach_exposure = "single_year", results for premature deaths are available for one year (the year of analysis). ```{r eval=FALSE, include=FALSE, echo=TRUE} results_pm_deaths$health_detailed$intermediate_calculations|> dplyr::filter(sex == "female") |> dplyr::pull(projection_if_unexposed_by_age_and_year) |> purrr::pluck(1) ``` ```{r eval=TRUE, include=TRUE, echo=FALSE} results_pm_deaths$health_detailed$intermediate_calculations|> dplyr::filter(sex == "female") |> dplyr::pull(projection_if_unexposed_by_age_and_year) |> purrr::pluck(1) |> dplyr::slice( ( dplyr::n() - 8 ) : dplyr::n() ) |> knitr::kable() ``` ```{r eval=TRUE, include=FALSE, echo=FALSE} rm(results_pm_deaths) ``` # YLD ### 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 [@Kim2022_kspm]. * 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 ```{r} 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 ```{r eval=FALSE, include=FALSE, echo=TRUE} results_pm_copd_yld$health_main ``` ```{r eval=TRUE, include=TRUE, echo=FALSE} results_pm_copd_yld$health_main |> dplyr::select(erf_ci, impact) |> knitr::kable() ``` # DALY ### 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 [@GBD2020_tl]. $$DALY = YLL + YLD$$ ### Function call This is possible using the function `daly()`. ```{r} results_daly <- daly( output_attribute_yll = results_pm_yll, output_attribute_yld = results_pm_copd_yld ) ``` ### Main results YLL, YLD & DALY ```{r eval=FALSE, include=FALSE, echo=TRUE} results_daly$health_main |> select(impact_yll_rounded, impact_yld_rounded, impact_rounded) ``` ```{r eval=TRUE, include=TRUE, echo=FALSE} results_daly$health_main |> dplyr::select(impact_yll_rounded, impact_yld_rounded, impact_rounded) |> knitr::kable() ``` ```{r eval=TRUE, include=FALSE, echo=FALSE} rm(results_daly, results_pm_yll, results_pm_copd_yld) ``` # Modification of scenarios ### Goal E.g., to quantify health impacts using `attribute_health`in an scenario B very similar to a previous scenario A. ### Function call ```{r} 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_health`in an existing scenario, e.g. `scenario_A`. ```{r} 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. ```{r} 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) ``` # Comparison of two health scenarios ### 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) [@WHO2014_book] - Population impact fraction (PIF) [@WHO2003_report; @Murray2003_spbm; @Askari2020_ijph]. 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). #### Population Impact Fraction (PIF) The Population Impact Fraction (PIF) is defined as the proportional change in disease or mortality when exposure to a risk factor is changed, for instance due to an intervention. ##### General Integral Form The most general equation describing this mathematically is an integral form [@WHO2003_report; @Murray2003_spbm]: $$PIF = \frac{\int rr\_at\_exp(x)PE(x)dx - \int rr\_at\_exp(x)PE'(x)dx}{\int rr\_at\_exp(x)PE(x)dx}$$ Where: * $x$ = exposure level * $PE(x)$ = population distribution of exposure * $PE'(x)$ = alternative population distribution of exposure * $rr\_at\_exp(x)$ = relative risk at exposure level compared to the reference level ##### Categorical Exposure Form If the population exposure is described as a categorical rather than continuous exposure, the integrals may be converted to sums [@WHO2003_report; @Murray2003_spbm]: $$PIF = \frac{\sum rr\_at\_exp_{i} \times PE_{i} - \sum rr\_at\_exp_{i}PE'_{i}}{\sum rr\_at\_exp_{i}PE_{i}}$$ Where: * $i$ = the exposure category (e.g., in bins of 1 $\mu g/m^3$ $PM_{2.5}$ or 5 dB noise exposure) * $PE_i$ = fraction of population in exposure category $i$ * $PE'_i$ = fraction of population in category $i$ for alternative (ideal) exposure scenario * $rr\_at\_exp_i$ = relative risk for exposure category level $i$ compared to the reference level ##### Population weighted mean concentration form Finally, if the exposure is provided as the population weighted mean concentration, the equation for the PIF is reduced to: $$PIF = \frac{rr\_at\_exp - rr\_at\_exp_{alt}}{rr}$$ Where: * $rr\_at\_exp$ = relative risk at the exposure level * $rr\_at\_exp_{alt}$ = relative risk at the exposure level for the alternative exposure scenario ### Function call 1. Use `attribute_health()` to calculate burden of scenarios A & B. ```{r} 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) ``` ```{r} scenario_B <- attribute_mod( output_attribute = scenario_A, exp_central = 6 ) ``` 2. Use `compare()` to compare scenarios A & B. ```{r} 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 ```{r eval=FALSE, include=FALSE, echo=TRUE} results_comparison$health_main ``` ```{r eval=TRUE, include=TRUE, echo=FALSE} results_comparison[["health_main"]] |> dplyr::select( impact, impact_rounded, impact_scen_1, impact_scen_2, bhd, dplyr::starts_with("exp_"), -dplyr::starts_with("exp_ci"), # remove col "exp_ci" dplyr::starts_with("rr_con")) |> dplyr::slice_head() |> knitr::kable() ``` ### 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) ```{r, echo=FALSE, eval=TRUE, include=FALSE} rm(results_comparison, scenario_A, scenario_B) ``` # Two correlated exposures ### 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 [@Strak2024_report] identified three approaches to add up attributable health impacts from correlated exposures: - Additive approach [@Steenland2006-e]: $$PAF_{additive} = PAF_{exposure1} + PAF_{exposure2}$$ - Multiplicative approach [@Jerrett2013-oup]: $$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 [@Steenland2006-e]: $$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()`. ```{r} 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" ) ``` ```{r, echo=FALSE, eval=TRUE, include=FALSE} rm(results_no2, results_pm) ``` ### Main results ```{r eval=FALSE, include=TRUE, echo=TRUE} results_multiplicative$health_main ``` ```{r eval=TRUE, include=TRUE, echo=FALSE} results_multiplicative$health_main |> dplyr::select(impact_rounded) |> knitr::kable() ``` ```{r eval=TRUE, include=FALSE, echo=FALSE} rm(results_multiplicative) ``` # Standardization ### 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 [@GBD2020_tldemo; @Ahmad2001_report]. 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 ```{r} 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: ```{r} print(results$health_main$impact_per_100k_inhab) ``` Age group-specific impact rate: ```{r} print(results$health_detailed$results_raw$impact_per_100k_inhab) ``` ```{r eval=TRUE, include=FALSE, echo=FALSE} rm(results) ``` # Preparation of exposure data ### Goal E.g., to determine population-weighted mean PM2.5 exposure for several neighborhoods of Brussels (Belgium) ### Methodology The `healthiar`function `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 ```{r eval=spatial_packages_available} # 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`. ```{r echo=FALSE, eval=spatial_packages_available} knitr::kable(as.data.frame(pwm$exposure_main)) ``` ```{r eval=spatial_packages_available, include=FALSE, echo=FALSE} rm(exdat_pwm_1) rm(exdat_pwm_2) rm(pwm) ``` # Cut-off vs. threshold {#cutoff-vs-threshold} ### Goal E.g., to quantify health impacts with an exposure-response function whose effect threshold is lower than the lowest exposure level that the assessment quantifies, 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 [@Engelmann2025_report]. Because the risk is expressed relative to the threshold, the threshold is also the **counterfactual exposure** of the assessment, i.e. the exposure of the hypothetical scenario that the actual exposure is compared with [@Murray2003_e]. * `cutoff_...` is the exposure level below which no health impacts are quantified. It determines which exposures are assessed, not the level that they are compared with. Usually, the cut-off is assumed to be the effect threshold, and then it is the counterfactual exposure as well. This is how the term is used in health risk assessment guidance, e.g., for air pollution [@WHO2025_report], where one single level plays both roles. Accordingly, there are three possible combinations: **a) Cut-off equal to the effect threshold** This is the default situation. Since both threshold and cut-off are usually identical, `healthiar` assumes that the values are the same, if one of them is missing. Therefore, in that case, it is enough to enter only one value. 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, threshold) - threshold) $$ The same result is obtained by entering only `cutoff_...`, only `threshold`, or both with the same value. **b) 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 get the risk at the reference level: $$ 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). **c) 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. ```{r} 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'. ```{r} #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 select(erf_ci, monetized_impact) |> knitr::kable() ``` 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. ```{r} results <- monetize( impact = 1151, valuation = 100 ) ``` ## Cost-benefit analysis ### Goal (e.g.) E.g., to perform an economic evaluation for an intervention by comparing its benefits and costs via Cost-Benefit Analysis (CBA). ### Methodology The CBA is a type of economic evaluation that compares the costs and the benefits of an intervention, considering both measures expressed in monetary terms. To perform a CBA, you can use the function `cba()`. This approach requires monetizing benefits so they can be directly compared with costs. Since interventions typically generate costs and benefits over multi-year time horizons, discounting is a common practice to obtain the present value of future costs and benefits. Depending on the reference guidelines, the discount rate can be specified as the same for costs and benefits or different across them. The outputs of a Cost-Benefit Analysis can be expressed as three main indicators [@Boardman2018_book]: - intervention’s net benefit: the difference between monetized benefits and costs - Cost-Benefit Ratio (CBR): monetized benefits divided by costs and - Return on Investment (ROI): return generated per unit of expenditure by relating net benefits to the intervention’s costs. An intervention is recommended from a Cost-Benefit Analysis perspective, if it yields a positive net benefit or a positive ROI, or equivalently, a CBR greater than one, meaning that the intervention’s monetized benefits exceed its costs. These three outputs are available when running `cba()` and are calculated considering the following formulas. **Net Benefit** $$net\_benefit = benefit - cost$$ **Cost-Benefit Ratio (CBR)** $$cbr = \frac{benefit}{cost}$$ **Return on Investment (ROI)** $$roi = \frac{benefit - cost}{cost} \times 100$$ ### Function call 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 the functions `healthiar` and `cba()` ```{r} 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 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 ```{r} cba$cba_main |> dplyr::select(benefit, cost, net_benefit) |> knitr::kable() ``` ```{r, echo=FALSE, eval=TRUE, include=FALSE} rm(cba) ``` 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! # Social aspects ## Health impact attributable to social indicator ### Goal E.g., to estimate the health impact that is theoretically attributable to the difference in degree of deprivation of the population exposed. ### Methodology Taking into account socio-economic indicators, e.g. a multiple deprivation index [@Mogin2025_ejph], the differences in attributable health impacts across the study areas can be estimated [@Renard2019_bmc; @Otavova_2022_bmc]. Social inequalities are quantified as the difference between the least deprived areas (the last n-quantile) and - the most deprived areas or - the population overall. These differences can be - absolute or - relative. #### Difference most deprived vs. least deprived $$ absolute\_quantile = first - last $$ Where: * $absolute\_quantile$ = Absolute difference between quantiles. * $first$ = Average health impacts in _most_ deprived quantile. * $last$ = Average health impacts in _least_ deprived quantile. $$ relative\_quantile = \frac{absolute\_quantile}{last} $$ #### Difference overall vs. least deprived $$ absolute\_overall = overall - last $$ Where: * $absolute\_overall$ = Absolute difference regarding the overall average. * $overall$ = _Overall_ average health impacts in the study area. * $last$ = Average health impacts in _least_ deprived quantile. If you assume that the least deprived areas are similar to counterfactual cases (no exposure to deprivation), the relative difference regarding the overall average health impact could be interpreted as some kind of relative risk attributable to social inequalities. ### Function call First, quantify health impacts. ```{r eval=TRUE, include=TRUE, echo=TRUE} 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) ``` Second, use the function `socialize()` entering the whole output of `attribute_health()` in the argument `output_attribute`. ```{r eval=TRUE, include=TRUE, echo=TRUE} 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) ``` Alternatively, you can directly enter the health impact in the `socialize()` argument `impact`. ```{r eval=TRUE, include=TRUE, echo=TRUE} 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) ``` ### Main results ```{r echo=FALSE} social$social_main |> dplyr::select(c(parameter, difference_type, difference_compared_with, difference_value, comment)) |> print() ``` ## Multiple deprivation index ### Goal E.g., to estimate the multiple deprivation index (MDI) to use it for the argument `social_indicator` in the function `socialize()`. ### Methodology Socio-economic indicators (e.g., education level, employment status and family structure) can be condensed into a multiple deprivation index (MDI) [@Mogin2025_ejph]. For this purpose, the indicators can be normalized using min-max scaling. The reliability of the MDI can be assessed using Cronbach's alpha [@Cronbach1951_p]. $$ \alpha = \frac{k}{k - 1} \left( 1 - \frac{\sum_{i=1}^{k} \sigma^2_{y_i}}{\sigma^2_x} \right) $$ where: * $k$ is the number of items/variables. * $\sigma^2_{y_i}$ is the variance of the $i$-th item. * $\sum_{i=1}^{k} \sigma^2_{y_i}$ is the sum of the variances of all items. * $\sigma^2_x$ is the total variance of the observed total scores (the sum of all items). To apply this approach, you should ensure that the data set is as complete as possible. Otherwise, you can try to impute missing data using: - Time-Based Imputation: Linear regression based on historical trends if prior years' data is complete. - Indicator-Based Imputation: Multiple linear regression if the missing indicator correlates strongly with others. Imputation models should have an R^2 greater than or equal to 0.7. If R^2 lower than 0.7, consider alternative data sources or methods. ### Function call ```{r eval=TRUE, include=TRUE, echo=TRUE} 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 ```{r eval=FALSE, include=TRUE, echo=TRUE} mdi$mdi_main |> select(geo_id_micro, MDI, MDI_index) ``` ```{r eval=TRUE, include=TRUE, echo=FALSE} mdi$mdi_main |> dplyr::select(geo_id_micro, MDI, MDI_index) |> knitr::kable() ``` The function assesses the reliability of the MDI based on the Cronbach's alpha value as follows: - 0.9 and higher: Excellent reliability - between 0.8 (included) and 0.9: Good reliability - between 0.7 (included) and 0.8: Acceptable reliability - between 0.6 (included) and 0.7: Questionable reliability - lower than 0.6: Poor reliability ### 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 ```{r fig.alt="Boxplot of Normalized Indicators and MDI", eval=TRUE, include=TRUE, echo=TRUE} eval(mdi$mdi_detailed$boxplot) ``` Analogeously, to reproduce the histogram run ```{r fig.alt="Histogram of MDI with normal curve", eval=TRUE, include=TRUE, echo=TRUE} eval(mdi$mdi_detailed$histogram) ``` ```{r, eval=TRUE, include=FALSE, echo=FALSE} rm(mdi) ``` --------------------------------------------------------------------------- # Inside pipes ## Pipe |> `healthiar` can be used inside the _native_ pipes `|>`. See the example below. ```{r eval=FALSE, include=TRUE, echo = TRUE} 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()`. ```{r} 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 })() ``` ## Pipe %>% `healthiar` can also be used inside _magrittr_ pipes `%>%` as follows. ```{r eval=FALSE, include=TRUE, echo = TRUE} 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 and visualize Exporting and visualizing results is out of scope of `healthiar`. To export and visualize, you can make use of existing functions in other packages beyond `healthiar` as indicated below. ## Export results Export as `.csv` file ```{r eval=FALSE, include=FALSE, echo=TRUE} write.csv(x = results_pm_copd$health_main, file = "exported_results/results_pm_copd.csv") ``` Save as `.Rdata` file ```{r eval=FALSE, include=FALSE, echo=TRUE} save(results_pm_copd, file = "exported_results/results_pm_copd.Rdata") ``` Export to Excel (as `.xlsx` file) ```{r eval=FALSE, include=FALSE, echo=TRUE} openxlsx::write.xlsx(x = results_pm_copd$health_main, file = "exported_results/results_pm_copd.xlsx") ``` ## Visualize results Visualization is out of scope of `healthiar`. You can visualize in: - R using base programming or packages such as `ggplot2` [@Wickham2016_ggplot], - Excel (export results first) or - Other tools. ------------------------------------------------------------------------ # Abbreviations BHD/bhd = baseline health data CI = confidence interval CBA/cba = cost-benefit analysis exp = exposure ERF = exposure-response function RR/rr = relative risk WHO = World Health Organization YLL/yll = years of life lost ------------------------------------------------------------------------ # References