Hard coded vs. columns

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
)

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 (Wickham 2014).

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)

Structure

The output of the healthiarfunction 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")

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

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

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

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

Main results

The main output contains aggregated results

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

Detailed results

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

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

Goal

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

Function call

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

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

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

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

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

Main results

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

Detailed results

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

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

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

Detailed results

Let’s inspect the detailed results:

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

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

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

• …

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

Note: only a selection of columns are shown.

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

Function call

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.

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

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.

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

results_age_group$health_detailed$results_by_age_group |> 
  dplyr::select(age_group, impact_rounded, exp, bhd) |> 
  knitr::kable()

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.

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.

results_sex$health_detailed$results_by_sex |> 
  dplyr::select(sex, impact_rounded, exp, bhd) |> 
  knitr::kable()

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.

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

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

Goal

E.g., to quantify attributable health impacts stratified by age, sex and additional sub-group e.g. education level.

Function call

output_attribute <- healthiar::attribute_health(
    rr_central = 1.063,
    rr_increment = 10,
    erf_shape = "log_linear",
    cutoff_central =  0,
    age_group = base::rep(c("50_and_younger", "50_plus"), each = 4, times= 2),
    sex = base::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 = base::rep(c("a", "b"), each = 8),
    info = base::data.frame(
      education = base::rep(c("without_master", "with_master"), times = 8)) # education level
  )

Results by all sub-groups

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

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

Function call

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

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

Main results

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

Detailed results

Attributable YLL results

• per year

• per age (group)

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

are available.

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

Results per year

Note: only a selection of years is shown.

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

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

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.

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

Main results

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

Detailed results

Attributable premature deaths results

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

• per age (group)

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

are available.

Note: We 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).

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.

Goal

E.g., to monetize the attributable health impact of a policy that will have health benefits five years from now.

Methodology

In health economic evaluations and economic burden of disease assessments, health impacts may need to be converted into monetary values. For this purpose, you can use monetize().

Several valuation metrics are available, depending on how outcomes are quantified in natural or health units (e.g. cases reduced, deaths prevented, reductions in mortality risk, life-years gained, QALYs gained, DALYs averted). Common metrics include the Value of a Statistical Life (VSL) (OECD 2025) , the Value of a Life-Year (VOLY) (Hammitt 2007) and the Value of a Quality-Adjusted Life-Year (VAQALY) (Bobinac et al. 2010).

Discounting is the practice of converting future costs (or health impacts as previous step to valuating them) into their present value. The underlying rationale is that the value placed on outcomes declines as they occur further in the future. Therefore, future costs and effects are converted into present-value terms to make them comparable over time (Attema et al. 2018).

Discounting is implemented by selecting a discount rate, which is used to compute a discount factor for each time period. This factor is then multiplied by the corresponding future cost (or effect) to express it in present-value terms.

If you need the discounted values of a cost or health outcome, you can call the healthiar function discount(). If you just need the discount factor, you can alternatively call get_discount_factor() (entering is_deflation = TRUE). If you just need the inflation factor, you can get_inflation_factor().

Different functional forms can be used to apply discounting. The most common is exponential discounting, also referred to as constant discounting, since outcomes are discounted proportionally as time increases. An alternative is hyperbolic discounting, which tends to better capture human behavior by discounting the near present more heavily than outcomes further in the future (Lipman and Attema 2024)

See below the equations that are used behind these functions.

\[monetized\_impact = impact \times valuation \times discount\_factor \times deflator\_factor \times real\_growth\_factor\] The arguments discount_factor, deflator_factor and real_growth_factor are only used if a value is entered in the arguments discount_rate, ingflation_rate and real_growth_rate respectively
(otherwise ignored).

Function call

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

Main results

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

Two folders are added:

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

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

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

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

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

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 (Boardman et al. 2018): - 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()

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

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

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

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 (Mogin et al. 2025), the differences in attributable health impacts across the study areas can be estimated (Renard et al. 2019; Otavova et al. 2022).

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.

Function call

First, quantify health impacts.

 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.

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.

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

#> # A tibble: 4 × 5
#>   parameter      difference_type difference_compared_…¹ difference_value comment
#>   <chr>          <chr>           <chr>                             <dbl> <chr>  
#> 1 impact_rate_s… absolute        last_quantile                   11.5    <NA>   
#> 2 impact_rate_s… relative        last_quantile                    0.193  <NA>   
#> 3 impact_rate_s… absolute        overall                         -0.834  It can…
#> 4 impact_rate_s… relative        overall                         -0.0143 It can…
#> # ℹ abbreviated name: ¹​difference_compared_with

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) (Mogin et al. 2025). For this purpose, the indicators can be normalized using min-max scaling.

The reliability of the MDI can be assessed using Cronbach’s alpha (Cronbach 1951).

\[ \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

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

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

Main results

Function output includes:

• mdi_main, a tibble containing the BEST-COST MDI

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

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

eval(mdi$mdi_detailed$boxplot)

Analogeously, to reproduce the histogram run

eval(mdi$mdi_detailed$histogram)