diff --git a/.github/workflows/build.yaml b/.github/workflows/build.yaml index d7cb75e..2047cae 100644 --- a/.github/workflows/build.yaml +++ b/.github/workflows/build.yaml @@ -38,6 +38,7 @@ jobs: - name: Install dependencies run: | + options(timeout = 600) install.packages(c("remotes","rcmdcheck","stringi"), repos = "https://cloud.r-project.org") install.packages(c("evaluate"), repos = "https://cloud.r-project.org") install.packages(c("knitr","rmarkdown","htmlwidgets", "plm"), repos = "https://cloud.r-project.org") diff --git a/.github/workflows/test_coverage.yaml b/.github/workflows/test_coverage.yaml index 052dcca..1d319e2 100644 --- a/.github/workflows/test_coverage.yaml +++ b/.github/workflows/test_coverage.yaml @@ -39,10 +39,12 @@ jobs: sudo apt-get install -y \ libcurl4-openssl-dev libssl-dev libxml2-dev \ libudunits2-dev libgdal-dev libgeos-dev libproj-dev \ - libmagick++-dev libavfilter-dev libharfbuzz-dev libfribidi-dev + libmagick++-dev libavfilter-dev libharfbuzz-dev libfribidi-dev \ + libuv1-dev - name: Install R dependencies run: | + options(timeout = 600) install.packages(c("remotes", "rcmdcheck", "covr"), repos = "https://cloud.r-project.org") remotes::install_deps(dependencies = TRUE) shell: Rscript {0} diff --git a/NAMESPACE b/NAMESPACE index 135eefe..7a8ef13 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,7 +2,6 @@ export(calc_ResidEm_grp) export(calc_hap_impacts) -export(calc_hia_adder) export(create_panel) export(fit_model) export(get_gdp_ctry) diff --git a/NEWS.md b/NEWS.md index ca23734..f0403b5 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,4 +7,21 @@ ## New features -## Bug Fixes \ No newline at end of file +* `fit_model()` re-specified: linear year trend replaces year fixed effects, NMVOC dropped (non-significant under every specification/SE estimator tested), and Driscoll-Kraay standard errors are now computed and returned (`vcov_dk`), justified by confirmed serial correlation and cross-sectional dependence in the panel. +* The country-level bias-correction term is now computed internally by `fit_model()` from the fitted model itself (`bias_adder$bias.adder`), replacing the static `rhap::hia_adder` dataset, which has been removed along with `calc_hia_adder()`. +* `by_gr = TRUE` in `calc_hap_impacts()` now reports within-country income-decile group *shares* of population, PM2.5+NOx exposure, and GDP (table + pie chart), instead of a fabricated group-level health-impact estimate. The regression model is fit on between-country variation only, so applying it to within-country decile inputs was an unvalidated extrapolation; shares of GCAM's own group-level data carry no such risk. +* `fit_model()`'s `bias_adder` now also returns `reliability_ratio` (and `calc_hap_impacts()`'s output a derived `reliability` "high"/"medium"/"low" label), flagging countries whose absolute-level prediction rests mostly on the bias correction rather than the regression's own covariates -- computed already, previously discarded. +* `calc_hap_impacts()` gains a `fit_result` argument: pass a previously-computed `fit_model(HIA_var)` result to skip refitting the national regression (including its Driscoll-Kraay vcov) on every call, e.g. when comparing many GCAM scenarios for the same `HIA_var`. +* Added a new "Validating the econometric model" vignette: a reproducible record of the panel/inference diagnostics, the (rejected) GDP x PM2.5 interaction test, the additive-vs-multiplicative bias-adder investigation, and the GBD real-data validation methodology. + +## Bug Fixes + +* Fixed the country-level bias correction silently clamping a negative "predicted" death/YLL/DALY rate to exactly 0 (occurred for ~12% of observations on the model's own training data). It now floors at a small fraction of the naive, uncorrected prediction instead. (A multiplicative bias correction, which prevents negative values by construction, was also tried; it was rejected after validation against real GBD data showed it performing far worse than the additive form -- see `fit_model()`'s documentation and the econometric-model vignette.) +* Fixed `fit_model(HIA_var = "yll")`: valid input silently resolved to `NULL` due to a mismatched lookup key (`"ylls"` vs. `"yll"`). +* Fixed a rerun of `calc_hap_impacts(map = TRUE)` against a not-freshly-cleared `output/` directory failing because the map-output file mover matched, and tried to rename, its own destination folder. +* Fixed several `.RData` test snapshots (`calc_hap_impacts`, `calc_ResidEm_grp` x2, `fit_model`) being compared vacuously: `load()` injects an object under whatever name it was *saved* as, which for these files happened to match the local variable holding the freshly-computed test result (`testOutput`), silently overwriting it before `expect_equal()` ran -- so the test was comparing the stale snapshot to itself, regardless of what the function under test actually returned. Also removed the now-redundant, structurally-brittle `fit_model` snapshot in favor of targeted structural/sign/parameter checks. + +## Package/performance + +* Removed dead code (`flsp_pc_gr`, `flsp_pc_ctry_gr`, `gdp_pc_ctry_gr`) left over from the `by_gr` redesign. +* GCAM group-level queries (`subregional population`, `subregional income`) and their downstream processing now only run when `by_gr = TRUE` (the shares they feed are otherwise never consumed) -- skips two extra GCAM queries and their processing on every default-mode (`by_gr = FALSE`) call. \ No newline at end of file diff --git a/R/calc_hap_impacts.R b/R/calc_hap_impacts.R index 787e09c..17cb728 100644 --- a/R/calc_hap_impacts.R +++ b/R/calc_hap_impacts.R @@ -3,7 +3,15 @@ #' #' Run the model to estimate the health impacts attributable to household air pollution for GCAM scenarios #' @keywords health impacts, HAP -#' @return Health impacts attributable to HAP for all the selected years +#' @return Health impacts attributable to HAP for all the selected years, as +#' a tibble of scenario, country, year, pred_var, pred_value, +#' pred_value_normalized, reliability_ratio, and reliability. The last two +#' flag how much of a country's predicted level rests on fit_model()'s +#' bias-adder correction rather than on the regression's own covariates +#' (reliability_ratio = abs(bias.adder)/naive prediction; reliability is a +#' "high"/"medium"/"low" label over that ratio, thresholds 0.5 and 2) -- +#' "low" flags a country whose absolute-level prediction the model has +#' little independent basis for, not necessarily a wrong number. #' @param db_path Path to the GCAM database #' @param query_path Path to the query file #' @param db_name Name of the GCAM database @@ -16,14 +24,27 @@ #' @param anim If set to TRUE, produces multi-year animations. By default=TRUE #' @param HIA_var Health metric to be predicted. c("deaths", "yll", "dalys"). By default = deaths #' @param normalized Transform the output to "normalized" values. By default = FALSE -#' @param by_gr Estimate damages at group level. Just for illustrative purposes. By default = FALSE +#' @param by_gr Report within-country income-decile group shares of +#' population, PM2.5+NOx exposure, and GDP as a CSV table plus a pie-chart +#' snapshot at final_db_year (output/by_gr/). Descriptive only -- it does +#' not estimate group-level health impacts, since the regression model is +#' fit on between-country variation and was never validated against +#' within-country/sub-national variation. By default = FALSE +#' @param fit_result Optional pre-computed result from fit_model(HIA_var), +#' e.g. from a previous call, to reuse instead of refitting the national +#' regression from scratch. fit_model() doesn't depend on GCAM scenario +#' data at all, so refitting it (including its Driscoll-Kraay vcov, which +#' isn't cheap) on every calc_hap_impacts() call is pure waste when +#' comparing many scenarios for the same HIA_var. Must have been fit with +#' the same HIA_var as this call; a mismatch raises an error rather than +#' silently using the wrong model. By default = NULL (fits internally). #' @importFrom magrittr %>% #' @export calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_name = NULL, prj_name, scen_name, queries = "queries_rhap.xml", final_db_year = 2100, saveOutput = TRUE, map = FALSE, anim = TRUE, HIA_var = "deaths", - normalized = FALSE, by_gr = FALSE) { + normalized = FALSE, by_gr = FALSE, fit_result = NULL) { Country <- country <- sector <- scenario <- region <- year <- group <- ghg <- Units <- value <- adj <- value_reg <- dec_share <- Pollutant <- `ISO 3` <- Percentatge <- iso3 <- iso <- ssp <- `gcam-consumer` <- value_adj <- building <- @@ -33,7 +54,8 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n pred_var <- pred_value_per_100K <- pred_value_per_100K_adj <- HIA <- PrimPM25 <- NOx <- NMVOC <- OC <- BC <- PrimPM25_per_100k <- NOx_per_100k <- VOC_per_100k <- `GCAM Region` <- Percentage <- adj_country <- country_map <- - gdp_pc_dol2011_ppp_gr <- gdp_pc_dol2011_ppp <- mapCountries <- . <- NULL + gdp_pc_dol2011_ppp <- mapCountries <- pop <- + metric <- share <- reliability_ratio <- reliability <- . <- NULL # Check user input if (!HIA_var %in% c("deaths", "yll", "dalys")) { @@ -79,7 +101,12 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n by_gr )) } - + if (!is.null(fit_result) && !identical(fit_result$HIA_var, HIA_var)) { + stop(sprintf( + "Error: the supplied fit_result was fit for HIA_var = '%s', but this call uses HIA_var = '%s'. Pass a fit_result from fit_model(HIA_var = '%s'), or leave fit_result = NULL to fit it internally.", + fit_result$HIA_var, HIA_var, HIA_var + )) + } @@ -260,21 +287,6 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n dplyr::ungroup() # Extract floorspace data and combine it with subregional population to compute per capita floorspace - flsp_pc_gr <- rgcam::getQuery(prj, "building floorspace") %>% - dplyr::filter(grepl("resid", building)) %>% - dplyr::mutate( - group = gsub("resid_", "", building), - flsp_m2 = value * 1E9 - ) %>% - dplyr::select(scenario, region, year, group, flsp_m2) %>% - dplyr::filter(year <= final_db_year) %>% - gcamdata::left_join_error_no_match(pop_gr, by = c("scenario", "region", "year", "group")) %>% - dplyr::mutate( - flsp_pc_gr = flsp_m2 / pop, - unit = "m2/pers" - ) %>% - dplyr::select(scenario, region, year, group, unit, flsp_pc_gr) - flsp_pc <- rgcam::getQuery(prj, "building floorspace") %>% dplyr::filter(grepl("resid", building)) %>% dplyr::mutate( @@ -295,18 +307,6 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n # Assume that al countries within the region have similar flps_pc - flsp_pc_ctry_gr <- flsp_pc_gr %>% - dplyr::left_join(reg_to_ctry, by = "region", relationship = "many-to-many") %>% - # Add SSP narrative associated with the scenario. Use SSP2 by default if no other SSP is specified in the scenario name - dplyr::mutate( - ssp = "SSP2", - ssp = dplyr::if_else(grepl("SSP1", scenario), "SSP1", ssp), - ssp = dplyr::if_else(grepl("SSP3", scenario), "SSP3", ssp), - ssp = dplyr::if_else(grepl("SSP4", scenario), "SSP4", ssp), - ssp = dplyr::if_else(grepl("SSP5", scenario), "SSP5", ssp) - ) %>% - dplyr::select(scenario, ssp, country, group, year, flsp_pc_gr) - flsp_pc_ctry <- flsp_pc %>% dplyr::left_join(reg_to_ctry, by = "region", relationship = "many-to-many") %>% # Add SSP narrative associated with the scenario. Use SSP2 by default if no other SSP is specified in the scenario name @@ -362,136 +362,138 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n country = region ) - # Process Population: Population is evenly distributed across groups, but could be updated - pop_share <- rgcam::getQuery(prj, "subregional population") %>% - dplyr::filter( - year <= final_db_year, - grepl("resid", `gcam-consumer`) - ) %>% - tidyr::separate(`gcam-consumer`, c("sector", "group"), sep = "_") %>% - # fix (Taiwan): if the `subregional population` query misses some values, - # compute them using the `population by region` query - dplyr::full_join( - rgcam::getQuery(prj, "population by region") %>% - dplyr::filter(year <= final_db_year) %>% - dplyr::group_by(scenario, region, year, value) %>% - tidyr::expand(group = unique(pop_gr$group)) %>% - dplyr::ungroup() %>% - dplyr::mutate( - value_adj = value / n_groups, - sector = "resid" - ) %>% - dplyr::select(-value), - by = c("scenario", "region", "sector", "group", "year") - ) %>% - dplyr::mutate(value = dplyr::if_else(is.na(value), value_adj, value)) %>% - # rename and aggregate sectors - dplyr::group_by(scenario, region, year) %>% - dplyr::mutate(value_agg = sum(value)) %>% - dplyr::ungroup() %>% - dplyr::mutate(share_pop = round(value / value_agg, 1)) %>% - dplyr::select(scenario, region, year, group, share_pop) %>% - # Temp fix for Taiwan - dplyr::mutate(share_pop = dplyr::if_else(region == "Taiwan", 0.1, share_pop)) %>% - # expand shares to countries - dplyr::left_join(reg_to_ctry, by = "region", relationship = "many-to-many") %>% - # Add SSP narrative associated with the scenario. Use SSP2 by default if no other SSP is specified in the scenario name - dplyr::mutate( - ssp = "SSP2", - ssp = dplyr::if_else(grepl("SSP1", scenario), "SSP1", ssp), - ssp = dplyr::if_else(grepl("SSP3", scenario), "SSP3", ssp), - ssp = dplyr::if_else(grepl("SSP4", scenario), "SSP4", ssp), - ssp = dplyr::if_else(grepl("SSP5", scenario), "SSP5", ssp) - ) %>% - dplyr::select(scenario, ssp, country, year, group, share_pop) %>% - dplyr::mutate(year = as.character(year)) - - pop_ctry_gr <- pop_ctry %>% - dplyr::filter(year <= final_db_year) %>% - # filter only ssps used - dplyr::filter(ssp %in% unique(pop_share$ssp)) %>% - # add groups - gcamdata::repeat_add_columns(tibble::tibble(group = unique(em_shares_gr$group))) %>% - # adjust country names to match - dplyr::left_join(rhap::adj_ctry, by = "country") %>% - dplyr::mutate(country = dplyr::if_else(is.na(adj_country), country, adj_country)) %>% - dplyr::select(-adj_country) %>% - # filter out Mayotte and Curasao - dplyr::filter(country %!in% c("Mayotte", "Curasao")) %>% - gcamdata::left_join_error_no_match(pop_share, by = c("ssp", "country", "year", "group")) %>% - # compute group-level population - dplyr::mutate(pop_gr = pop * share_pop) %>% - dplyr::select(scenario, ssp, country, year, group, pop_gr) - - # Process GDP: First need to calculate the income shares by GCAM_region and period (to be applied to all countries within each region) - gdp_share <- rgcam::getQuery(prj, "subregional income") %>% - dplyr::filter( - grepl("resid", `gcam-consumer`), - year <= final_db_year - ) %>% - tidyr::separate(`gcam-consumer`, c("sector", "group"), sep = "_") %>% - dplyr::mutate(gdp_pc = value * 1E3) %>% - dplyr::select(-value, -Units) %>% - gcamdata::left_join_error_no_match( - pop_gr %>% - dplyr::rename(value = pop), - by = c("scenario", "group", "region", "year") - ) %>% - dplyr::mutate(gdp = gdp_pc * (value * 1E3)) %>% - dplyr::group_by(scenario, region, year) %>% - dplyr::mutate(gdp_agg = sum(gdp)) %>% - dplyr::ungroup() %>% - dplyr::mutate(share_gdp = round(gdp / gdp_agg, 5)) %>% - dplyr::select(scenario, region, year, group, share_gdp) %>% - # expand shares to countries - dplyr::left_join(reg_to_ctry, by = "region", relationship = "many-to-many") %>% - # Add SSP narrative associated with the scenario. Use SSP2 by default if no other SSP is specified in the scenario name - dplyr::mutate( - ssp = "SSP2", - ssp = dplyr::if_else(grepl("SSP1", scenario), "SSP1", ssp), - ssp = dplyr::if_else(grepl("SSP3", scenario), "SSP3", ssp), - ssp = dplyr::if_else(grepl("SSP4", scenario), "SSP4", ssp), - ssp = dplyr::if_else(grepl("SSP5", scenario), "SSP5", ssp) - ) %>% - dplyr::select(scenario, ssp, country, year, group, share_gdp) %>% - dplyr::mutate(year = as.character(year)) - - # Adjust Taiwan using China's shares - gdp_share_twn <- gdp_share %>% - dplyr::filter(country == "China") %>% - dplyr::mutate(country = "Taiwan") - - gdp_share <- gdp_share %>% - dplyr::filter(country != "Taiwan") %>% - dplyr::bind_rows(gdp_share_twn) - - - gdp_ctry_gr <- gdp_ctry %>% - # filter only ssps used - dplyr::filter( - ssp %in% unique(gdp_share$ssp), - year <= final_db_year - ) %>% - # add groups - gcamdata::repeat_add_columns(tibble::tibble(group = unique(em_shares_gr$group))) %>% - # adjust country names to match - dplyr::left_join(rhap::adj_ctry, by = "country") %>% - dplyr::mutate(country = dplyr::if_else(is.na(adj_country), country, adj_country)) %>% - dplyr::select(-adj_country) %>% - # filter out Mayotte and Curasao - dplyr::filter(country %!in% c("Mayotte", "Curasao")) %>% - gcamdata::left_join_error_no_match(gdp_share, by = c("ssp", "country", "year", "group")) %>% - # compute group-level population - dplyr::mutate(gdp_dol2011_ppp_gr = gdp_dol2011_ppp * share_gdp) %>% - dplyr::select(scenario, ssp, country, year, group, gdp_dol2011_ppp_gr) + # pop_share/pop_ctry_gr and gdp_share/gdp_ctry_gr are only consumed by the + # by_gr shares block further down (gated on the same by_gr && saveOutput + # condition). Skip these GCAM queries and their processing entirely when + # by_gr = FALSE (the default) -- there's no point hitting + # "subregional population"/"subregional income" and downscaling them to + # country x group level if nothing will ever read the result. + if (by_gr && saveOutput) { + # Process Population: Population is evenly distributed across groups, but could be updated + pop_share <- rgcam::getQuery(prj, "subregional population") %>% + dplyr::filter( + year <= final_db_year, + grepl("resid", `gcam-consumer`) + ) %>% + tidyr::separate(`gcam-consumer`, c("sector", "group"), sep = "_") %>% + # fix (Taiwan): if the `subregional population` query misses some values, + # compute them using the `population by region` query + dplyr::full_join( + rgcam::getQuery(prj, "population by region") %>% + dplyr::filter(year <= final_db_year) %>% + dplyr::group_by(scenario, region, year, value) %>% + tidyr::expand(group = unique(pop_gr$group)) %>% + dplyr::ungroup() %>% + dplyr::mutate( + value_adj = value / n_groups, + sector = "resid" + ) %>% + dplyr::select(-value), + by = c("scenario", "region", "sector", "group", "year") + ) %>% + dplyr::mutate(value = dplyr::if_else(is.na(value), value_adj, value)) %>% + # rename and aggregate sectors + dplyr::group_by(scenario, region, year) %>% + dplyr::mutate(value_agg = sum(value)) %>% + dplyr::ungroup() %>% + dplyr::mutate(share_pop = round(value / value_agg, 1)) %>% + dplyr::select(scenario, region, year, group, share_pop) %>% + # Temp fix for Taiwan + dplyr::mutate(share_pop = dplyr::if_else(region == "Taiwan", 0.1, share_pop)) %>% + # expand shares to countries + dplyr::left_join(reg_to_ctry, by = "region", relationship = "many-to-many") %>% + # Add SSP narrative associated with the scenario. Use SSP2 by default if no other SSP is specified in the scenario name + dplyr::mutate( + ssp = "SSP2", + ssp = dplyr::if_else(grepl("SSP1", scenario), "SSP1", ssp), + ssp = dplyr::if_else(grepl("SSP3", scenario), "SSP3", ssp), + ssp = dplyr::if_else(grepl("SSP4", scenario), "SSP4", ssp), + ssp = dplyr::if_else(grepl("SSP5", scenario), "SSP5", ssp) + ) %>% + dplyr::select(scenario, ssp, country, year, group, share_pop) %>% + dplyr::mutate(year = as.character(year)) + + pop_ctry_gr <- pop_ctry %>% + dplyr::filter(year <= final_db_year) %>% + # filter only ssps used + dplyr::filter(ssp %in% unique(pop_share$ssp)) %>% + # add groups + gcamdata::repeat_add_columns(tibble::tibble(group = unique(em_shares_gr$group))) %>% + # adjust country names to match + dplyr::left_join(rhap::adj_ctry, by = "country") %>% + dplyr::mutate(country = dplyr::if_else(is.na(adj_country), country, adj_country)) %>% + dplyr::select(-adj_country) %>% + # filter out Mayotte and Curasao + dplyr::filter(country %!in% c("Mayotte", "Curasao")) %>% + gcamdata::left_join_error_no_match(pop_share, by = c("ssp", "country", "year", "group")) %>% + # compute group-level population + dplyr::mutate(pop_gr = pop * share_pop) %>% + dplyr::select(scenario, ssp, country, year, group, pop_gr) + + # Process GDP: First need to calculate the income shares by GCAM_region and period (to be applied to all countries within each region) + gdp_share <- rgcam::getQuery(prj, "subregional income") %>% + dplyr::filter( + grepl("resid", `gcam-consumer`), + year <= final_db_year + ) %>% + tidyr::separate(`gcam-consumer`, c("sector", "group"), sep = "_") %>% + dplyr::mutate(gdp_pc = value * 1E3) %>% + dplyr::select(-value, -Units) %>% + gcamdata::left_join_error_no_match( + pop_gr %>% + dplyr::rename(value = pop), + by = c("scenario", "group", "region", "year") + ) %>% + dplyr::mutate(gdp = gdp_pc * (value * 1E3)) %>% + dplyr::group_by(scenario, region, year) %>% + dplyr::mutate(gdp_agg = sum(gdp)) %>% + dplyr::ungroup() %>% + dplyr::mutate(share_gdp = round(gdp / gdp_agg, 5)) %>% + dplyr::select(scenario, region, year, group, share_gdp) %>% + # expand shares to countries + dplyr::left_join(reg_to_ctry, by = "region", relationship = "many-to-many") %>% + # Add SSP narrative associated with the scenario. Use SSP2 by default if no other SSP is specified in the scenario name + dplyr::mutate( + ssp = "SSP2", + ssp = dplyr::if_else(grepl("SSP1", scenario), "SSP1", ssp), + ssp = dplyr::if_else(grepl("SSP3", scenario), "SSP3", ssp), + ssp = dplyr::if_else(grepl("SSP4", scenario), "SSP4", ssp), + ssp = dplyr::if_else(grepl("SSP5", scenario), "SSP5", ssp) + ) %>% + dplyr::select(scenario, ssp, country, year, group, share_gdp) %>% + dplyr::mutate(year = as.character(year)) + + # Adjust Taiwan using China's shares + gdp_share_twn <- gdp_share %>% + dplyr::filter(country == "China") %>% + dplyr::mutate(country = "Taiwan") + + gdp_share <- gdp_share %>% + dplyr::filter(country != "Taiwan") %>% + dplyr::bind_rows(gdp_share_twn) + + gdp_ctry_gr <- gdp_ctry %>% + # filter only ssps used + dplyr::filter( + ssp %in% unique(gdp_share$ssp), + year <= final_db_year + ) %>% + # add groups + gcamdata::repeat_add_columns(tibble::tibble(group = unique(em_shares_gr$group))) %>% + # adjust country names to match + dplyr::left_join(rhap::adj_ctry, by = "country") %>% + dplyr::mutate(country = dplyr::if_else(is.na(adj_country), country, adj_country)) %>% + dplyr::select(-adj_country) %>% + # filter out Mayotte and Curasao + dplyr::filter(country %!in% c("Mayotte", "Curasao")) %>% + gcamdata::left_join_error_no_match(gdp_share, by = c("ssp", "country", "year", "group")) %>% + # compute group-level population + dplyr::mutate(gdp_dol2011_ppp_gr = gdp_dol2011_ppp * share_gdp) %>% + dplyr::select(scenario, ssp, country, year, group, gdp_dol2011_ppp_gr) + } # Combine processed datasets to get per capita GDP - gdp_pc_ctry_gr <- gdp_ctry_gr %>% - gcamdata::left_join_error_no_match(pop_ctry_gr, by = c("scenario", "ssp", "country", "year", "group")) %>% - dplyr::mutate(gdp_pc_dol2011_ppp_gr = gdp_dol2011_ppp_gr / pop_gr) %>% - dplyr::select(scenario, ssp, country, year, group, gdp_pc_dol2011_ppp_gr) - gdp_pc_ctry <- gdp_ctry %>% tibble::as_tibble() %>% gcamdata::left_join_error_no_match( @@ -513,54 +515,22 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n #----- # COMBINE THE DATA AND TRANSFORM IT TO MODEL VARIABLES - # 1- Calculate the outputs at group level - - output_gr <- em_ctry_gr %>% - gcamdata::left_join_error_no_match(flsp_pc_ctry_gr, by = c("scenario", "ssp", "country", "group", "year")) %>% - dplyr::filter(year >= min(unique(gdp_pc_ctry_gr$year))) %>% - dplyr::mutate(year = as.character(year)) %>% - # Filter out some small countries not in the SSP database - dplyr::filter(country %!in% c( - "Bermuda", "Cook Islands", "Cook Islands", "Dominica", "Falkland Islands (Malvinas)", - "Faroe Islands", "Gibraltar", "Guadeloupe", "Greenland", "Saint Kitts And Nevis", - "Liechtenstein", "Marshall Islands", "Montserrat", "Martinique", "Niue", "Palau", - "Reunion", "Saint Pierre And Miquelon", "Isle Of Man", "Turks And Caicos", "Tokelau", - "Wallis And Futuna", "Virgin Islands, British", "Kosovo", "Cayman Islands", - "American Samoa", "Samoa" - )) %>% - gcamdata::left_join_error_no_match(gdp_pc_ctry_gr, by = c("scenario", "ssp", "country", "group", "year")) %>% - # filter pollutants in the regression model: BC/OC, NOx, and VOC - dplyr::filter(ghg %in% rhap::panel_pollutants) %>% - tidyr::pivot_wider( - names_from = "ghg", - values_from = "em_ctry_gr" - ) %>% - # Combine BC and OC - dplyr::mutate(PrimPM25 = BC + OC) %>% - dplyr::select(-BC, -OC) %>% - # Add population to compute pollutants per 100K - gcamdata::left_join_error_no_match(pop_ctry_gr, by = c("scenario", "ssp", "country", "group", "year")) %>% - dplyr::mutate( - PrimPM25_per_100k = (PrimPM25 / pop_gr) * 100000, - NOx_per_100k = (NOx / pop_gr) * 100000, - VOC_per_100k = (NMVOC / pop_gr) * 100000 - ) %>% - # Add logarithms - dplyr::mutate( - log_PrimPM25_per_100k = log(PrimPM25_per_100k), - log_NOx_per_100k = log(NOx_per_100k), - log_VOC_per_100k = log(VOC_per_100k), - log_gdppc_ppp_dol2011 = log(gdp_pc_dol2011_ppp_gr), - log_flsp = log(flsp_pc_gr) - ) %>% - dplyr::select(scenario, country_name = country, group, year, dplyr::starts_with("log"), pop_gr) %>% - # adjust country names to match to panel data - gcamdata::left_join_error_no_match(rhap::adj_ctry_output, by = "country_name") %>% - dplyr::mutate(country_name = dplyr::if_else(data_name == "", country_name, data_name)) %>% - dplyr::select(-data_name) - + # Fit the model once and reuse everywhere below. Previously fit_model() + # was called twice (once for predictable_regions via old positional + # indexing [[2]], again later for model.fixed/bias_adder) — wasteful + # (refits the whole plm model twice) and the [[2]] access breaks now + # that fit_model() returns a named list. A caller running many scenarios + # for the same HIA_var can skip this refit entirely by passing a + # previously-computed fit_model(HIA_var) result via the fit_result + # argument (validated against HIA_var above). + if (is.null(fit_result)) { + fit_result <- fit_model(HIA_var = HIA_var) + } + model.fixed <- fit_result$model.fixed + bias_adder <- fit_result$bias_adder + predictable_regions <- fit_result$predictable_regions - # 2- Calculate the outputs at country level + # 1- Calculate the outputs at country level output <- em_ctry %>% gcamdata::left_join_error_no_match(flsp_pc_ctry, by = c("scenario", "ssp", "country", "year")) %>% dplyr::filter(year >= min(unique(gdp_pc_ctry$year))) %>% @@ -613,30 +583,45 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n # adjust Rou dplyr::mutate(country_name = dplyr::if_else(country_name == "Roumania", "Romania", country_name)) %>% # remove not predictable regions - dplyr::filter(country_name %in% unique(fit_model(HIA_var = HIA_var)[[2]])) %>% + dplyr::filter(country_name %in% predictable_regions) %>% dplyr::mutate(across(where(is.numeric), ~ ifelse(is.finite(.), ., NA_real_))) %>% - tidyr::drop_na() + tidyr::drop_na() %>% + # Plain numeric year column referenced by the model's trend term + # (year_num). Must exist as a real column, not derived via + # as.numeric(year) at predict time — see fit_model() notes on why + # as.numeric() on a pdata.frame's index column is unsafe. + dplyr::mutate(year_num = as.numeric(year)) #----- - # PREDICTION - # Get model to subtract coefficients and predict - model.fixed <- fit_model(HIA_var = HIA_var)[[1]] - - # Transform data to panel and predict - output.panel <- plm::pdata.frame(output, index = c("country_name", "year")) + # PREDICTION + # model.fixed / bias_adder were already computed once above, alongside + # predictable_regions, and are reused here rather than refit. + # + # NOTE: predict() is called directly on the plain `output` data frame, + # NOT on a plm::pdata.frame(). output contains multiple scenarios (and + # output_gr additionally multiple groups) sharing the same + # country_name+year combination, which is exactly what + # plm::pdata.frame()'s index requires to be unique — building a + # pdata.frame here would error ("duplicate couples (id-time)") as soon + # as more than one scenario is processed at once. Out-of-sample + # prediction here only needs coefficients %*% regressors (no lag(), + # diff(), or other panel-specific transform), so a plain data.frame is + # sufficient and sidesteps the uniqueness requirement entirely. # Predict - output$pred_value <- stats::predict(model.fixed, output.panel) + # suppressWarnings: plm::predict.plm() warns that newdata isn't a + # pdata.frame and falls back to coefficients-only prediction (no fixed + # effect). That's expected and intentional here (see NOTE above) -- the + # bias_adder computed in fit_model() is exactly what corrects for the + # missing fixed effect, so this isn't a problem to surface as a warning. + output$pred_value <- suppressWarnings(stats::predict(model.fixed, newdata = output)) output_fin <- output %>% dplyr::mutate(pred_var = paste0("pred_log_", HIA_var, "_per_100K")) %>% gcamdata::left_join_error_no_match( - rhap::hia_adder %>% - tibble::as_tibble() %>% - dplyr::filter(HIA == HIA_var) %>% - dplyr::select(country_name = country, bias.adder) %>% + bias_adder %>% gcamdata::repeat_add_columns(tibble::tibble(year = unique(output$year))), by = c("country_name", "year") ) %>% @@ -646,17 +631,37 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n pred_var = gsub("log_", "", pred_var) ) %>% dplyr::mutate( + # bias.adder is additive (see fit_model() for why an additive, not + # multiplicative, correction is used). It can be negative enough to + # push the corrected value below 0 for some countries; floor at a + # small fraction of the naive (uncorrected) prediction rather than + # clamping to exactly 0, so a real (if highly uncertain) low value + # is reported instead of a hard, discontinuous "zero risk". pred_value_per_100K_adj = pred_value_per_100K + bias.adder, + pred_value_per_100K_adj = pmax(pred_value_per_100K_adj, pred_value_per_100K * 0.01), pred_value_per_100K_adj = round(pred_value_per_100K_adj, 2), pred_value = round(pred_value_per_100K_adj * pop / 100000, 0), pred_var = gsub("pred_", "", pred_var), pred_var = gsub("_per_100K", "", pred_var) ) %>% + # reliability: a convenience label over reliability_ratio (from + # fit_model()'s bias_adder, joined in above) -- how much of this + # country's predicted level rests on the bias correction rather than + # the model's own covariates. "low" doesn't mean the number is wrong; + # it means the model has little independent basis for that country's + # absolute level, and the prediction is closer to "last observed rate, + # lightly adjusted for scenario changes" than to a model-driven estimate. dplyr::mutate( - pred_value = dplyr::if_else(as.numeric(pred_value) < 0, 0, as.numeric(pred_value)), - pred_value_per_100K_adj = dplyr::if_else(as.numeric(pred_value_per_100K_adj) < 0, 0, as.numeric(pred_value_per_100K_adj)) + reliability = dplyr::case_when( + reliability_ratio <= 0.5 ~ "high", + reliability_ratio <= 2 ~ "medium", + TRUE ~ "low" + ) ) %>% - dplyr::select(scenario, country = country_name, year, pred_var, pred_value, pred_value_normalized = pred_value_per_100K_adj) + dplyr::select( + scenario, country = country_name, year, pred_var, pred_value, + pred_value_normalized = pred_value_per_100K_adj, reliability_ratio, reliability + ) # Create a function to write the data (by scenario) @@ -672,47 +677,98 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n } - # If by group = TRUE, add a complementary estimation at group level + # If by group = TRUE, report descriptive group-level shares + # + # NOTE: by_gr reports within-country income-decile group SHARES of + # population, PM2.5+NOx exposure, and GDP -- it does not estimate + # group-level health impacts. An earlier version applied the + # country-level regression model (fit on between-country variation) to + # within-country decile-level inputs to produce a fabricated + # "group-level death rate". That is an extrapolation the model was never + # fit to support (there is no sub-national variation in the training + # data to validate it against), so it has been replaced with the + # directly-observable shares GCAM already computes per group -- these + # carry no such extrapolation risk. For detailed per-pollutant group + # breakdowns (not just combined PM2.5+NOx), see calc_ResidEm_grp(). if (saveOutput == TRUE & by_gr == TRUE) { # Create the directory if they do not exist: if (!dir.exists("output/by_gr")) dir.create("output/by_gr") - output.panel.gr <- plm::pdata.frame(output_gr, index = c("country_name", "year")) - - output_gr$pred_value <- stats::predict(model.fixed, output.panel.gr) + # Each share's denominator is the sum across groups WITHIN the same + # frame (not a join back to the separate country-level frames, e.g. + # pop_ctry/gdp_ctry) -- those go through a different country-name + # normalization step (rhap::adj_ctry) than the *_ctry_gr frames do, so + # a name-based join between them silently drops non-matching + # countries. Summing within the group frame itself sidesteps that + # entirely and is exact by construction (group shares of a country's + # own group breakdown must sum to that country's total). + pop_share_gr <- pop_ctry_gr %>% + dplyr::group_by(scenario, ssp, country, year) %>% + dplyr::mutate(pop_ctry_total = sum(pop_gr)) %>% + dplyr::ungroup() %>% + dplyr::mutate(metric = "population", share = pop_gr / pop_ctry_total) %>% + dplyr::select(scenario, country, year, group, metric, share) + + exposure_share_gr <- em_ctry_gr %>% + dplyr::filter(ghg %in% c("BC", "OC", "NOx")) %>% + dplyr::group_by(scenario, ssp, country, group, year) %>% + dplyr::summarise(em_ctry_gr = sum(em_ctry_gr), .groups = "drop") %>% + dplyr::group_by(scenario, ssp, country, year) %>% + dplyr::mutate(em_ctry_total = sum(em_ctry_gr)) %>% + dplyr::ungroup() %>% + dplyr::mutate(metric = "exposure", share = em_ctry_gr / em_ctry_total, year = as.character(year)) %>% + dplyr::select(scenario, country, year, group, metric, share) + + gdp_share_gr <- gdp_ctry_gr %>% + dplyr::group_by(scenario, ssp, country, year) %>% + dplyr::mutate(gdp_ctry_total = sum(gdp_dol2011_ppp_gr)) %>% + dplyr::ungroup() %>% + dplyr::mutate(metric = "gdp", share = gdp_dol2011_ppp_gr / gdp_ctry_total) %>% + dplyr::select(scenario, country, year, group, metric, share) + + shares_gr <- dplyr::bind_rows(pop_share_gr, exposure_share_gr, gdp_share_gr) %>% + dplyr::mutate(share = round(share, 4)) + + utils::write.csv(shares_gr, "output/by_gr/shares_byGR.csv", + row.names = FALSE, fileEncoding = "UTF-8" + ) - output_fin_gr <- output_gr %>% - dplyr::mutate(pred_var = paste0("pred_log_", HIA_var, "_per_100K")) %>% - dplyr::left_join( - rhap::hia_adder %>% - tibble::as_tibble() %>% - dplyr::filter(HIA == HIA_var) %>% - dplyr::select(country_name = country, bias.adder) %>% - gcamdata::repeat_add_columns(tibble::tibble(group = unique(em_ctry_gr$group))) %>% - gcamdata::repeat_add_columns(tibble::tibble(year = unique(output$year))), - by = c("country_name", "group", "year") - ) %>% - dplyr::filter(stats::complete.cases(.)) %>% - dplyr::mutate(bias.adder = bias.adder / length(unique(em_ctry_gr$group))) %>% - dplyr::mutate( - pred_value_per_100K = exp(pred_value), - pred_var = gsub("log_", "", pred_var) - ) %>% - dplyr::mutate( - pred_value_per_100K_adj = pred_value_per_100K + bias.adder, - pred_value_per_100K_adj = round(pred_value_per_100K_adj, 2), - pred_value = round(pred_value_per_100K_adj * pop_gr / 100000, 0), - pred_var = gsub("pred_", "", pred_var), - pred_var = gsub("_per_100K", "", pred_var) - ) %>% - dplyr::mutate( - pred_value = dplyr::if_else(as.numeric(pred_value) < 0, 0, as.numeric(pred_value)), - pred_value_per_100K_adj = dplyr::if_else(as.numeric(pred_value_per_100K_adj) < 0, 0, as.numeric(pred_value_per_100K_adj)) - ) %>% - dplyr::select(scenario, country = country_name, group, year, pred_var, pred_value, pred_value_normalized = pred_value_per_100K_adj) %>% - utils::write.csv(paste0("output/by_gr/", "HAP_", unique(HIA_var), "_byGR", ".csv"), - row.names = FALSE, fileEncoding = "UTF-8" + # Pie chart: a single snapshot at final_db_year (one panel per + # scenario x metric). For other years, filter shares_gr directly and + # build a custom plot -- the CSV above has the full scenario/year range. + # A handful of small countries can have a 0/0 (NA) share for a given + # metric where GCAM reports zero for both the group and country total; + # drop those from the plot rather than let ggplot silently warn. + pie_data <- shares_gr %>% + dplyr::filter(year == as.character(final_db_year), !is.na(share)) + + if (nrow(pie_data) > 0) { + pl <- ggplot2::ggplot( + pie_data, + ggplot2::aes( + x = "", y = share, + fill = factor(group, levels = c("d1", "d2", "d3", "d4", "d5", "d6", "d7", "d8", "d9", "d10")) + ) + ) + + ggplot2::geom_bar(width = 1, stat = "identity") + + ggplot2::theme_void() + + ggplot2::facet_grid(scenario ~ metric) + + ggplot2::labs(x = "", y = "") + + ggplot2::theme( + legend.title = ggplot2::element_blank(), + legend.position = "bottom", + strip.text = ggplot2::element_text(size = 10) + ) + + ggplot2::coord_polar("y", start = 0) + + ggplot2::scale_fill_manual(values = c( + "gray20", "gray50", "#ad440c", "#ef8e27", "#d01c2a", + "darkorchid3", "#507fab", "deepskyblue1", "#11d081", "#00931d" + )) + ggplot2::ggsave(pl, + file = file.path("output/by_gr", paste0("shares_byGR_", final_db_year, ".png")), + width = 250, height = 120, units = "mm" ) + } } # Add map @@ -742,6 +798,11 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n subRegion = unique(rmap::mapCountries$region) ) + # NOTE: must be assigned to the global environment (<<-), not local. + # rmap::map() -> rmap::map_plot() -> rmap::map_find() looks up + # "mapCountries" via get(), which only finds it in .GlobalEnv, not in + # calc_hap_impacts()'s local frame. Using <- here breaks map = TRUE + # with "object 'mapCountries' not found". mapCountries <<- rmap::mapCountries @@ -818,7 +879,12 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n unlink(paste("output/maps/map", "allScen", y, sep = "_"), recursive = TRUE) } # 2.2. gather all figures in "map_allScen_byYear" new folder - files_to_move <- list.files(path = file.path("output/maps"), pattern = "^map_allScen_", full.names = TRUE) + # NOTE: pattern must match only the .png files, not "^map_allScen_" alone. + # On a rerun against a not-freshly-cleaned output/ dir, the destination + # folder itself ("map_allScen_byYear") already exists and would match a + # bare "^map_allScen_" prefix, causing file.rename() to try moving that + # folder into itself. + files_to_move <- list.files(path = file.path("output/maps"), pattern = "^map_allScen_.*\\.png$", full.names = TRUE) if (!dir.exists("output/maps/map_allScen_byYear")) dir.create("output/maps/map_allScen_byYear") success <- sapply(files_to_move, function(file) { file.rename(file, file.path("output/maps/map_allScen_byYear", basename(file))) diff --git a/R/calc_hia_adder.R b/R/calc_hia_adder.R deleted file mode 100644 index 27ecfba..0000000 --- a/R/calc_hia_adder.R +++ /dev/null @@ -1,126 +0,0 @@ -#' calc_hia_adder -#' -#' @description -#' Calculate the adder for HAP impact analysis -#' -#' The adder is calculated as the difference between the observed and estimated values in final data year -#' -#' @source IHME -#' @keywords Bias adder -#' @importFrom magrittr %>% -#' @export -#' @return Adder - -calc_hia_adder <- function() { - datadir <- paste0(getwd(), "/inst/extdata") - # ---------------------------------------------- - # Adjust the data - data <- rhap::panel_data %>% - dplyr::select(iso, country_name, year, pop, dplyr::starts_with("log"), continent, dev) %>% - dplyr::mutate(year = as.character(year)) %>% - dplyr::select(-log_AAP, -log_HDD_value, -log_CDD_value) %>% - dplyr::filter(stats::complete.cases(.)) - - # ---------------------------------------------- - # ---------------------------------------------- - # FIT THE MODEL - # ---------------------------------------------- - fixed_fin <- plm::plm( - log_Deaths_per_100k ~ log_PrimPM25_per_100k + - log_NOx_per_100k + - log_VOC_per_100k + - log_gdppc_ppp_dol2011 + - log_flsp, - data = data, - index = c("country_name", "year"), model = "within" - ) - summary(fixed_fin) - - # ---------------------------------------------- - # ---------------------------------------------- - # PREDICTIONS - # ---------------------------------------------- - data.panel <- plm::pdata.frame(data, index = c("country_name", "year")) - - data$pred_log_Deaths_per_100k <- predict(fixed_fin, data.panel) - - # Calculate adders - data_adder <- data %>% - dplyr::select(iso, country_name, year, pop, log_Deaths_per_100k, pred_log_Deaths_per_100k) %>% - dplyr::filter(year == max(year)) %>% - dplyr::mutate( - Deaths_per_100k = exp(log_Deaths_per_100k), - pred_Deaths_per_100k = exp(pred_log_Deaths_per_100k) - ) %>% - dplyr::mutate(bias.adder = Deaths_per_100k - pred_Deaths_per_100k) %>% - dplyr::select(iso, country = country_name, bias.adder) %>% - dplyr::mutate(HIA = "deaths") - - # Repeat the process for YLL and DALYs - - # 1- YLLs - fixed_fin_yll <- plm::plm( - log_YLL_per_100k ~ log_PrimPM25_per_100k + - log_NOx_per_100k + - log_VOC_per_100k + - log_gdppc_ppp_dol2011 + - log_flsp, - data = data, - index = c("country_name", "year"), model = "within" - ) - summary(fixed_fin_yll) - - data.panel_yll <- plm::pdata.frame(data, index = c("country_name", "year")) - - data$pred_log_YLL_per_100k <- predict(fixed_fin_yll, data.panel_yll) - - data_adder_yll <- data %>% - dplyr::select(iso, country_name, year, pop, log_YLL_per_100k, pred_log_YLL_per_100k) %>% - dplyr::filter(year == max(year)) %>% - dplyr::mutate( - YLL_per_100k = exp(log_YLL_per_100k), - pred_YLL_per_100k = exp(pred_log_YLL_per_100k) - ) %>% - dplyr::mutate(bias.adder = YLL_per_100k - pred_YLL_per_100k) %>% - dplyr::select(iso, country = country_name, bias.adder) %>% - dplyr::mutate(HIA = "yll") - - # 2- DALYs - fixed_fin_daly <- plm::plm( - log_DALY_per_100k ~ log_PrimPM25_per_100k + - log_NOx_per_100k + - log_VOC_per_100k + - log_gdppc_ppp_dol2011 + - log_flsp, - data = data, - index = c("country_name", "year"), model = "within" - ) - summary(fixed_fin_daly) - - data.panel_daly <- plm::pdata.frame(data, index = c("country_name", "year")) - - data$pred_log_DALY_per_100k <- predict(fixed_fin_daly, data.panel_daly) - - - data_adder_daly <- data %>% - dplyr::select(iso, country_name, year, pop, log_DALY_per_100k, pred_log_DALY_per_100k) %>% - dplyr::filter(year == max(year)) %>% - dplyr::mutate( - DALY_per_100k = exp(log_DALY_per_100k), - pred_DALY_per_100k = exp(pred_log_DALY_per_100k) - ) %>% - dplyr::mutate(bias.adder = DALY_per_100k - pred_DALY_per_100k) %>% - dplyr::select(iso, country = country_name, bias.adder) %>% - dplyr::mutate(HIA = "dalys") - - #--- - - # add all metrics - data_adder_fin <- dplyr::bind_rows( - data_adder, - data_adder_yll, - data_adder_daly - ) - - invisible(data_adder_fin) -} diff --git a/R/data.R b/R/data.R index 158d422..825e13f 100644 --- a/R/data.R +++ b/R/data.R @@ -262,19 +262,6 @@ "adj_ctry_map" -#' hia_adder -#' -#' @description Bias adder, calculated as the difference between the observed and "predicted" values in final observed year -#' @source Own calculations based on IHME data -#' @format .csv -#' @examples -#' \dontrun{ -#' library(rhap) -#' rhap::hia_adder -#' } -"hia_adder" - - #' all_pollutants #' #' @description List of pollutants from the residential sector diff --git a/R/fit_model.R b/R/fit_model.R index c0b3e97..a2f0846 100644 --- a/R/fit_model.R +++ b/R/fit_model.R @@ -1,21 +1,103 @@ #' fit_model #' #' @description -#' Fit the model using fixed effects model and plm function. +#' Fit the model using a fixed effects model (plm, entity within-estimator) +#' with a linear year trend, and return Driscoll-Kraay standard errors +#' alongside the fitted model. #' -#' For the model dplyr::selection, testing and validation see "./rhap/inst/extdata/model_testing.R" +#' Model selection notes (see "./rhap/inst/extdata/model_testing.R" for the +#' full diagnostic history): +#' - Full year fixed effects (factor(year)) over-absorb common trend that +#' overlaps with the pollutant decline across countries, washing out +#' pollutant coefficients almost entirely. A linear year trend gives +#' coefficients consistent with a spline(df=4) and full year-FE +#' specification, so it was preferred for parsimony and because it +#' extrapolates sensibly to future scenario years (unlike year dummies). +#' - log_VOC_per_100k was dropped: coefficient is ~0 and non-significant +#' under every specification and SE estimator tested (naive, Arellano +#' clustered, Driscoll-Kraay), and does not survive a joint Wald test +#' either alone or jointly with PM2.5/floorspace. +#' - Under Driscoll-Kraay SEs (robust to both serial correlation and +#' cross-sectional dependence — appropriate for a country-year panel +#' with common shocks like oil prices, recessions, and multilateral +#' pollution treaties), only log_NOx_per_100k and log_gdppc_ppp_dol2011 +#' are individually significant. log_PrimPM25_per_100k and log_flsp are +#' correctly signed (positive / negative respectively, as expected) but +#' not individually or jointly significant (joint Wald test p = 0.35). +#' They are retained on theoretical grounds (omitted-variable-bias risk +#' from dropping known-relevant regressors) but their point estimates +#' should be treated as directionally informative only, not precise, +#' for any scenario-comparison use of this model's coefficients. #' -#' @source Details on plm estimation can be found here: https://cran.r-project.org/web/packages/plm/plm.pdf#' +#' - The country-level bias adder previously stored in the external +#' rhap::hia_adder table is now computed internally from the fitted +#' model itself (see below), so it can never drift out of sync with +#' the coefficients used to produce it. It's calculated the same way +#' predict() will later be used on scenario data — coefficients only, +#' no fixed effect — so the gap it captures is exactly the gap that +#' needs to be added back at prediction time (country fixed effect + +#' any residual model bias), calibrated to recent observed years +#' rather than a single year to avoid baking in a one-off shock. +#' +#' - The bias adder is ADDITIVE (level scale, per_100k units), not +#' multiplicative, despite log(Y) = X*beta + alpha_i implying alpha_i is +#' technically a multiplicative scalar exp(alpha_i) in levels. A +#' multiplicative geometric-mean correction was tried and is more +#' "correct" in that narrow sense, but performed far worse in practice: +#' validated against real GBD out-of-sample data, it collapsed +#' cross-country correlation from ~0.86 (additive) to ~0.18. The reason +#' is that ~22% of countries (e.g. Uganda) have a common linear year +#' trend that doesn't match their own trajectory, so the gap between +#' the trend-implied (fixed-effect-free) prediction and reality keeps +#' growing even within the training years -- sometimes to a fixed +#' effect worth 1000x+ in levels. A multiplicative correction amplifies +#' that growing gap catastrophically when applied to a slightly +#' different (scenario) year; an additive correction, bounded by the +#' actual scale of the outcome variable, does not. See +#' inst/extdata/model_testing.R for the comparison. +#' +#' - IMPORTANT: the linear year trend is passed as an explicit numeric +#' column (year_num), not as.numeric(year) inline in the formula. +#' plm::pdata.frame() (used downstream at prediction time) converts +#' the "year" index column into a factor; as.numeric() on a factor +#' returns the factor's level codes, not the original calendar year, +#' and those codes are reassigned per-dataset — so a scenario dataset +#' spanning different years than the training data would get a +#' trend encoding disconnected from the one the coefficient was +#' estimated against. year_num is a plain (non-index) numeric column, +#' so it isn't touched by pdata.frame()'s factor coercion and stays +#' on a consistent calendar-year scale in both training and +#' prediction. +#' + +#' @source Details on plm estimation: https://cran.r-project.org/web/packages/plm/plm.pdf +#' @source Driscoll-Kraay SEs: plm::vcovSCC documentation #' @keywords Econometric model; fixed effects #' @param HIA_var Health metric to be predicted. c("deaths", "yll", "dalys") +#' @param n_years Number of most recent historical years (per country) to +#' average over when computing the bias adder. Default 5 — long enough +#' to smooth out a single-year shock, short enough to reflect the +#' country's current level rather than its full historical average. #' @importFrom magrittr %>% #' @export -#' @return Regression model and coefficients for prediction - -fit_model <- function(HIA_var) { +#' @return A list: model.fixed (plm object), predictable_regions (character +#' vector of countries in the estimation sample), vcov_dk (Driscoll-Kraay +#' variance-covariance matrix for inference — use this instead of the +#' default plm SEs for confidence intervals / significance testing on +#' this model's coefficients), bias_adder (tibble of country_name, +#' bias.adder, reliability_ratio — bias.adder is the per-country additive +#' level-calibration term, on the per_100k scale, to add to future +#' predictions [do not multiply]; reliability_ratio is +#' abs(bias.adder)/naive_prediction over the calibration window, where a +#' high value flags a country whose absolute-level prediction rests +#' almost entirely on the bias correction rather than the model's +#' covariates, and should be trusted less; replaces rhap::hia_adder), +#' HIA_var (the HIA_var this model was fit for, echoed back so a cached +#' fit_model() result can be validated against a different call site's +#' HIA_var before reuse). +fit_model <- function(HIA_var, n_years = 5) { iso <- country_name <- year <- pop <- continent <- dev <- log_AAP <- value <- - Model <- Scenario <- Region <- Variable <- Unit <- . <- NULL - + Model <- Scenario <- Region <- Variable <- Unit <- resid_per_100k <- . <- NULL # Check user input if (!HIA_var %in% c("deaths", "yll", "dalys")) { stop(sprintf( @@ -23,42 +105,131 @@ fit_model <- function(HIA_var) { HIA_var, paste(c("deaths", "yll", "dalys"), collapse = ", ") )) } - # Adjust the data data <- rhap::panel_data %>% dplyr::select(iso, country_name, year, pop, dplyr::starts_with("log"), continent, dev) %>% dplyr::mutate(year = as.character(year)) %>% dplyr::select(-log_AAP, -log_HDD_value, -log_CDD_value) %>% - dplyr::filter(stats::complete.cases(.)) + dplyr::filter(stats::complete.cases(.)) %>% + # Plain numeric copy of year, kept alongside the character "year" used + # for the panel index. plm::pdata.frame() (used at prediction time) + # coerces the index "year" column into a factor, and as.numeric() on + # a factor returns level codes, not the calendar year — those codes + # are reassigned per-dataset, so a scenario dataset with different + # years would get a trend encoding disconnected from training. This + # plain column isn't touched by that coercion, so it's what the + # trend term should reference. + dplyr::mutate(year_num = as.numeric(year)) predictable_regions <- unique(data$country_name) - # dplyr::select the dependent variable (deaths, YLLs, or DALYs) - # Create a named list to map HIA_var values to corresponding dep_var values HIA_var_map <- list( deaths = "log_Deaths_per_100k", - ylls = "log_YLL_per_100k", + yll = "log_YLL_per_100k", dalys = "log_DALY_per_100k" ) - # Assign the value from the named list dep_var <- HIA_var_map[[HIA_var]] - - # Fit the fixed effect model - - # write formula which depends on the dplyr::selected variable + # Final specification: entity FE + linear year trend, VOC excluded. + # See roxygen notes above for the diagnostic history behind this choice. model_formula <- stats::as.formula(paste( - dep_var, "~ log_PrimPM25_per_100k + log_NOx_per_100k + log_VOC_per_100k +", - "log_gdppc_ppp_dol2011 + log_flsp" + dep_var, "~ log_PrimPM25_per_100k + log_NOx_per_100k +", + "log_gdppc_ppp_dol2011 + log_flsp + year_num" )) - model.fixed <- plm::plm( model_formula, data = data, index = c("country_name", "year"), - model = "within" + model = "within", + effect = "individual" ) + # Driscoll-Kraay SEs: robust to serial correlation and cross-sectional + # dependence, appropriate for this country-year panel. Use this vcov + # for any inference (CIs, significance tests) on model.fixed's + # coefficients rather than the default plm summary() SEs. + vcov_dk <- plm::vcovSCC(model.fixed, type = "HC1", maxlag = 4) + + # ------------------------------------------------------------------ + # Bias adder: level-calibration term for out-of-sample prediction. + # predict() on a plm "within" model returns fitted values from the + # coefficients alone, without the entity fixed effect (this is also + # how scenario/future predictions are generated downstream, since + # scenario data has no fixed effect to draw on). That means a + # naive prediction will be systematically off by each country's + # fixed effect plus any remaining model bias. We recover that gap + # by predicting on the model's own training data and comparing to + # what was actually observed, averaged over the most recent + # n_years to avoid anchoring to a single noisy year. + # + # This is deliberately ADDITIVE (level scale, per_100k units), not a + # multiplicative ratio, even though log(Y) = X*beta + alpha_i implies + # alpha_i is technically a multiplicative scalar exp(alpha_i) in + # levels. A multiplicative geometric-mean version was tried and tested + # against real GBD data: it collapsed cross-country correlation from + # ~0.86 (additive) to ~0.18. Cause: the model's common linear year + # trend doesn't match every country's own trajectory (e.g. Uganda's + # mortality declines far more slowly than the trend assumes), so the + # gap between the trend-implied, fixed-effect-free prediction and + # reality keeps growing even within the training years -- for ~22% of + # countries the implied fixed effect exceeds 10x, and for some (Uganda) + # exceeds 1000x. A multiplicative correction amplifies that growing gap + # catastrophically the moment it's applied to a slightly different + # (scenario) year; a bounded additive correction does not. See + # inst/extdata/model_testing.R for the full comparison. + # + # The one real problem with the additive version -- it can produce a + # negative "corrected" rate for countries with a large negative bias + # adder -- is handled downstream in calc_hap_impacts() by flooring the + # corrected value at a small fraction of the naive (uncorrected) + # prediction, not by silently clamping to exactly 0. + # ------------------------------------------------------------------ + train_panel <- plm::pdata.frame(data, index = c("country_name", "year")) + data$pred_log <- stats::predict(model.fixed, train_panel) + data$observed_per_100k <- exp(data[[dep_var]]) + data$pred_per_100k <- exp(data$pred_log) + data$resid_per_100k <- data$observed_per_100k - data$pred_per_100k - return(list(model.fixed, predictable_regions)) + # reliability_ratio: |bias.adder| relative to the naive (uncorrected) + # prediction it's being added to, both averaged over the same n_years + # calibration window. A ratio near 0 means the model's own covariates + # already explain most of the country's level, so the correction (and + # by extension the country fixed effect it stands in for) is small + # relative to the prediction -- high confidence. A ratio >> 1 means the + # correction dwarfs the naive prediction (e.g. Uganda, where it's + # ~100x), meaning the model's covariates explain almost none of that + # country's absolute level and the fixed effect is doing all the work + # -- low confidence in the absolute-level prediction, even though the + # correction keeps it numerically well-behaved (see additive-vs- + # multiplicative note above). + bias_adder <- data %>% + dplyr::group_by(country_name) %>% + dplyr::filter(year_num %in% utils::tail(sort(unique(year_num)), n_years)) %>% + dplyr::summarise( + bias.adder = mean(resid_per_100k, na.rm = TRUE), + naive_pred_per_100k = mean(pred_per_100k, na.rm = TRUE), + .groups = "drop" + ) %>% + dplyr::mutate(reliability_ratio = round(abs(bias.adder) / naive_pred_per_100k, 3)) %>% + dplyr::select(country_name, bias.adder, reliability_ratio) + + return(list( + model.fixed = model.fixed, + predictable_regions = predictable_regions, + vcov_dk = vcov_dk, + bias_adder = bias_adder, + HIA_var = HIA_var + )) } +# KEEPING PRIMARY PM25 AND FLOORSPACE: + +# Omitted-variable-bias risk cuts one way. Dropping a theoretically-grounded, correctly-signed regressor because it's +# imprecisely estimated risks biasing the coefficients you do care about (NOx, GDP) if PM2.5 or floorspace has any real, +# non-zero relationship with mortality that's correlated with the retained regressors — which is plausible here, +# since PM2.5 in particular is one of the best-established health-relevant pollutants in the epidemiological literature. +# Keeping them costs you two degrees of freedom; wrongly dropping them could cost you unbiasedness elsewhere. + +# "Not significant" isn't "zero" — it's "we can't pin down the size precisely with this data," which is a power problem +# stemming from PM2.5's decline being unusually correlated with the common trend across countries, not evidence the effect +# doesn't exist. The point estimates stayed correctly signed and stable across every specification +# tried (linear, spline, decade, full FE), which is a different and stronger form of evidence than a single t-test. diff --git a/README.md b/README.md index 2926aca..a1c12e5 100644 --- a/README.md +++ b/README.md @@ -144,7 +144,7 @@ For a step-by-step walkthrough of the example, visit the official guide: [Run a [Back to Contents](#Contents) The primary function offered by this package is `calc_hap_impacts`, designed to estimate health impacts attributable to household air pollution (HAP) across a range of alternative GCAM scenarios. This function estimates scenario-specific health outcomes by seamlessly extracting and processing relevant data from GCAM databases or project files. -The extracted data encompasses critical socioeconomic and environmental parameters, enabling comprehensive analysis of HAP impacts. These parameters include per capita GDP, emissions of key pollutants such as primary PM2.5 (BC+OC), NOx, and VOCs, and per capita floorspace metrics. By integrating these diverse factors, `calc_hap_impacts` facilitates detailed assessments of how variations in economic development, pollutant emissions, and living conditions influence human health under different policy or technological scenarios. +The extracted data encompasses critical socioeconomic and environmental parameters, enabling comprehensive analysis of HAP impacts. These parameters include per capita GDP, emissions of key pollutants such as primary PM2.5 (BC+OC) and NOx, and per capita floorspace metrics. By integrating these diverse factors, `calc_hap_impacts` facilitates detailed assessments of how variations in economic development, pollutant emissions, and living conditions influence human health under different policy or technological scenarios. The package also includes an ancillary function, `calc_ResidEm_grp`, which determines the contribution of each within-region consumer group to household air pollution (HAP), with results categorized by region, year, and pollutant. Users can customize their analysis by specifying the desired region, time period, and pollutant through adjustable parameters when calling the function. This flexibility enables targeted assessments of how various consumer groups influence HAP across specific contexts, providing valuable insights for research and policy development. diff --git a/_pkgdown.yml b/_pkgdown.yml index 5caae79..4bd7182 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -18,6 +18,8 @@ navbar: href: articles/run_rhap.html - text: "Fitting the econometric model" href: articles/fit_model.html + - text: "Validating the econometric model" + href: articles/model_validation.html - text: "Analyze within-region contributions to HAP" href: articles/ResidEm_grp.html - icon: fa-file-code-o diff --git a/data/hia_adder.rda b/data/hia_adder.rda deleted file mode 100644 index 30a8dcf..0000000 Binary files a/data/hia_adder.rda and /dev/null differ diff --git a/inst/extdata/model_testing.R b/inst/extdata/model_testing.R index aedbc4e..c282b32 100644 --- a/inst/extdata/model_testing.R +++ b/inst/extdata/model_testing.R @@ -1,207 +1,225 @@ -library(dplyr) -library(plm) -library(ggplot2) -library(countrycode) - -datadir <- paste0(getwd(),"/inst/extdata") -# ---------------------------------------------- -# Adjust the data -data <- a %>% - select(iso, country_name, year, pop, starts_with("log"), continent, dev) %>% - select(-log_AAP) %>% - filter(complete.cases(.)) - -# data.europe <- data %>% filter(continent == "Europe") -# data.africa <- data %>% filter(continent == "Africa") -# data.america <- data %>% filter(continent == "Americas") -# data.asia <- data %>% filter(continent == "Asia & Oceania") - -# data.dev <- data %>% filter(dev == "Developed") -# data.nonDev <- data %>% filter(dev == "Developing") - -# data.kmeans2 <- kmeans(x = data[, c("log_gdppc_ppp_dol2011")], centers = 2) +# Diagnostic checks and model-selection history for fit_model(). +# Rebuild `data` and `model_formula` the same way fit_model() does before +# running any of this. # +# ============================================================ +# SUMMARY OF FINDINGS (see step-by-step diagnostics below) +# ============================================================ +# 1. plm(..., effect="twoways") throws "empty model" — caused by +# singleton entities/years in an unbalanced panel. Fixed via LSDV +# (entity FE via within + year FE via factor(year) dummies) rather +# than plm's native twoways transform. +# 2. Cross-sectional VIFs on the panel (not single-country) are all +# < 3 — the multicollinearity initially suspected was actually an +# artifact of restricting to a single country (Spain), which +# collapses all cross-sectional variation and leaves GDP, +# floorspace, and pollutants riding one shared trajectory (VIFs of +# 15-27 in that restricted sample). The full panel does not have +# this problem. +# 3. Full year dummies (factor(year)) over-absorb common trend that +# overlaps with the pollutant decline across countries (year +# coefficients are almost perfectly linear, ~-0.045/year), washing +# out PM2.5 and VOC coefficients almost entirely and leaving only +# NOx significant. +# 4. Comparing no trend control / linear trend / spline(df=4) / decade +# dummies / full year dummies: linear, spline, and full-year-FE all +# converge to similar pollutant coefficients; only the "no trend +# control" specification disagrees — meaning the "correct-looking" +# signs from that specification were actually the biased ones +# (omitted common-trend bias from declining pollution + declining +# mortality for unrelated reasons: medicine, smoking, diet, etc.). +# Linear trend was chosen: parsimonious, agrees with more complex +# alternatives, and extrapolates sensibly to future GCAM scenario +# years (year dummies do not extrapolate at all). +# 5. Under naive plm SEs, all 5 regressors are significant. Under +# Arellano country-clustered SEs, only GDP is unambiguously +# significant (NOx marginal, p=0.085). Under Driscoll-Kraay SEs +# (robust to both serial correlation AND cross-sectional +# dependence — more appropriate here, since global shocks like +# oil prices/recessions/pollution treaties hit many countries at +# once), NOx and GDP are both clearly significant; PM2.5, VOC, and +# floorspace are not. +# 6. VOC's coefficient is ~0 and non-significant under every +# specification and every SE estimator tested, and fails a joint +# Wald test (p ~ 0.99 individually / p=0.26 jointly with +# PM2.5+VOC). Dropped from the final model. +# 7. PM2.5 and floorspace remain correctly signed (positive / negative +# respectively) in the final model but are not individually or +# jointly significant under Driscoll-Kraay SEs (joint test +# p=0.35). Retained on theoretical/omitted-variable-bias grounds, +# but their point estimates should be read as directionally +# informative only, not precise, for scenario-comparison work. +# +# Final specification: +# log_HIA_per_100k ~ log_PrimPM25_per_100k + log_NOx_per_100k + +# log_gdppc_ppp_dol2011 + log_flsp + as.numeric(year) +# plm(model = "within", effect = "individual") +# inference via plm::vcovSCC(model, type = "HC1", maxlag = 4) + +library(plm) +library(car) +library(dplyr) +library(lmtest) + +# ------------------------------------------------------------ +# 0. Rebuild inputs (mirror fit_model()'s data prep) +# ------------------------------------------------------------ +data <- rhap::panel_data %>% + dplyr::select(iso, country_name, year, pop, dplyr::starts_with("log"), + continent, dev) %>% + dplyr::mutate(year = as.character(year)) %>% + dplyr::select(-log_AAP, -log_HDD_value, -log_CDD_value) %>% + dplyr::filter(stats::complete.cases(.)) + +model_formula <- log_Deaths_per_100k ~ log_PrimPM25_per_100k + + log_NOx_per_100k + log_VOC_per_100k + log_gdppc_ppp_dol2011 + log_flsp + +# ------------------------------------------------------------ +# 1. Panel structure / singleton check +# ------------------------------------------------------------ +# Run this BEFORE attempting effect="twoways" — singleton entities or +# years are the most common cause of plm's "empty model" error on an +# unbalanced panel. +check_panel_structure <- function(data) { + cat("\n---- pdim ----\n") + print(plm::pdim(data, index = c("country_name", "year"))) + + cat("\n---- Countries with only 1 year observed ----\n") + print(data %>% dplyr::count(country_name) %>% dplyr::filter(n == 1)) + + cat("\n---- Years with only 1 country observed ----\n") + print(data %>% dplyr::count(year) %>% dplyr::filter(n == 1)) +} + +# If singletons are found: # data <- data %>% -# mutate(cluster = data.kmeans2$cluster) %>% -# mutate(cluster = as.factor(cluster)) +# dplyr::add_count(country_name, name = "n_years") %>% +# dplyr::add_count(year, name = "n_countries") %>% +# dplyr::filter(n_years > 1, n_countries > 1) %>% +# dplyr::select(-n_years, -n_countries) + +# ------------------------------------------------------------ +# 2. LSDV as a robust alternative to effect="twoways" +# ------------------------------------------------------------ +fit_lsdv <- function(model_formula, data) { + formula_lsdv <- update(model_formula, . ~ . + factor(year)) + plm::plm(formula_lsdv, data = data, + index = c("country_name", "year"), + model = "within", effect = "individual") +} + +# ------------------------------------------------------------ +# 3. Multicollinearity check (run on the FULL panel, not a +# single-country subset — restricting to one country removes +# cross-sectional variation and inflates VIFs artificially) +# ------------------------------------------------------------ +check_multicollinearity <- function(model_formula, data) { + lm_check <- lm(model_formula, data = data) + cat("\n---- VIF ----\n") + print(car::vif(lm_check)) + + regressors <- all.vars(model_formula)[-1] + cat("\n---- Correlation matrix ----\n") + print(round(cor(data[, regressors], use = "complete.obs"), 2)) +} + +# ------------------------------------------------------------ +# 4. Trend-specification comparison +# (this is the key diagnostic: compare pollutant coefficients +# across no-trend / linear / spline / decade / full-year-FE to +# see whether "theoretically correct" signs are robust or are an +# artifact of omitted common trend) +# ------------------------------------------------------------ +compare_trend_specs <- function(model_formula, data) { + m_none <- plm::plm(model_formula, data = data, + index = c("country_name", "year"), + model = "within", effect = "individual") + + formula_linear <- update(model_formula, . ~ . + as.numeric(year)) + m_linear <- plm::plm(formula_linear, data = data, + index = c("country_name", "year"), + model = "within", effect = "individual") + + formula_spline <- update(model_formula, . ~ . + splines::ns(as.numeric(year), df = 4)) + m_spline <- plm::plm(formula_spline, data = data, + index = c("country_name", "year"), + model = "within", effect = "individual") + + data_decade <- data %>% dplyr::mutate(decade = paste0(substr(year, 1, 3), "0s")) + formula_decade <- update(model_formula, . ~ . + factor(decade)) + m_decade <- plm::plm(formula_decade, data = data_decade, + index = c("country_name", "year"), + model = "within", effect = "individual") + + m_full_year <- fit_lsdv(model_formula, data) + + regressors <- all.vars(model_formula)[-1] + comparison <- sapply( + list(none = m_none, linear = m_linear, spline = m_spline, + decade = m_decade, full_year = m_full_year), + function(m) coef(m)[regressors] + ) + print(comparison) + + invisible(list(none = m_none, linear = m_linear, spline = m_spline, + decade = m_decade, full_year = m_full_year)) +} + +# ------------------------------------------------------------ +# 5. Inference: naive vs. Arellano-clustered vs. Driscoll-Kraay +# ------------------------------------------------------------ +compare_se_estimators <- function(model) { + cat("\n---- Naive (default plm) SEs ----\n") + print(summary(model)$coefficients) + + cat("\n---- Arellano country-clustered SEs ----\n") + print(lmtest::coeftest(model, vcov = plm::vcovHC(model, method = "arellano", cluster = "group"))) + + cat("\n---- Driscoll-Kraay SEs (robust to serial correlation AND ----\n") + cat("---- cross-sectional dependence -- preferred for this panel) ----\n") + print(lmtest::coeftest(model, vcov = plm::vcovSCC(model, type = "HC1", maxlag = 4))) +} + +# ------------------------------------------------------------ +# 6. Joint significance test for weak-but-theoretically-motivated +# regressors (use this instead of relying on individual t-tests +# when deciding whether to drop a correctly-signed variable) +# ------------------------------------------------------------ +joint_test <- function(model, vars) { + hyp <- paste0(vars, " = 0") + car::linearHypothesis(model, hyp, + vcov = plm::vcovSCC(model, type = "HC1", maxlag = 4)) +} + +# ------------------------------------------------------------ +# 7. Final model +# ------------------------------------------------------------ +fit_final_model <- function(data, dep_var = "log_Deaths_per_100k") { + formula_final <- stats::as.formula(paste( + dep_var, "~ log_PrimPM25_per_100k + log_NOx_per_100k +", + "log_gdppc_ppp_dol2011 + log_flsp + as.numeric(year)" + )) + model.final <- plm::plm(formula_final, data = data, + index = c("country_name", "year"), + model = "within", effect = "individual") + vcov_dk <- plm::vcovSCC(model.final, type = "HC1", maxlag = 4) + list(model = model.final, vcov_dk = vcov_dk) +} + +# ------------------------------------------------------------ +# Example run-through +# ------------------------------------------------------------ +# check_panel_structure(data) +# check_multicollinearity(model_formula, data) +# trend_models <- compare_trend_specs(model_formula, data) +# compare_se_estimators(trend_models$linear) # -# data.c1 <- data %>% filter(cluster == "1") -# data.c2 <- data %>% filter(cluster == "2") - -# ---------------------------------------------- -# ---------------------------------------------- -# FIT THE MODEL -# ---------------------------------------------- -# 1- Test the fixed effects model -fixed <- plm(log_Deaths_per_100k ~ log_PrimPM25_per_100k + - log_NOx_per_100k + - log_VOC_per_100k + - log_gdppc_ppp_dol2011 + - log_flsp, - data = data, - index = c("country_name", "year"), model="within") -summary(fixed) - -# # HDD is non-significative and CDD has a strange effect -# fixed_noDD <- plm(log_Deaths_per_100k ~ log_PrimPM25_per_100k + log_NOx_per_100k + -# log_gdppc_ppp_dol2011 + -# log_flsp, -# data = data, -# index=c("country_name", "year"), model="within") -# summary(fixed_noDD) - -# AIC cannot be directly used with panel data -# R2 does is not largely reduced, so we select the no_DD model -fixed_fin <- plm(log_Deaths_per_100k ~ log_PrimPM25_per_100k + - log_NOx_per_100k + - log_VOC_per_100k + - log_gdppc_ppp_dol2011 + - log_flsp, - data = data, - index=c("country_name", "year"), model="within") -summary(fixed_fin) - -# ---------------------------------------------- -# 2- Use the Hausman test to check fixed effects are more accurate than random -# Fit a RE model -random <- plm(log_Deaths_per_100k ~ log_PrimPM25_per_100k + - log_NOx_per_100k + - log_VOC_per_100k + - log_gdppc_ppp_dol2011 + - log_flsp, - data = data, - index=c("country_name", "year"), model = "random") -summary(random) - -phtest(fixed_fin,random) #Hausman test - -# Given that p-value < 0.01 -> FIXED effects! - -# ---------------------------------------------- -# 3- Check multicolineality -# VIF does not work for FE or RE models -> need pooled model and check -pooled <- plm(log_Deaths_per_100k ~ log_PrimPM25_per_100k + - log_NOx_per_100k + - log_VOC_per_100k + - log_gdppc_ppp_dol2011 + - log_flsp, - data = data, - index=c("country_name", "year"), model="pooling") -summary(pooled) -car::vif(pooled) -#create vector of VIF values -vif_values <- car::vif(pooled) - -#create horizontal bar chart to display each VIF value -barplot(vif_values, main = "VIF Values", horiz = TRUE, col = "steelblue") -abline(v = 5, lwd = 3, lty = 2) - -# VIF values are around 2 -> no Multicolinearity. - -# ---------------------------------------------- -# 4- Test for time-fixed effects -# fixed_fin_time <- plm(log_Deaths_per_100k ~ log_PrimPM25_per_100k + log_NOx_per_100k + -# log_gdppc_ppp_dol2011 + -# log_flsp + -# factor(year), -# data = data, -# index=c("country_name", "year"), model="within") -# summary(fixed_fin_time) +# model_formula_novoc <- update(model_formula, . ~ . - log_VOC_per_100k + as.numeric(year)) +# model.novoc <- plm::plm(model_formula_novoc, data = data, +# index = c("country_name", "year"), +# model = "within", effect = "individual") +# compare_se_estimators(model.novoc) +# joint_test(model.novoc, c("log_PrimPM25_per_100k", "log_flsp")) # -# pFtest(fixed_fin, fixed_fin_time) - -# ---------------------------------------------- -# Ver el resumen del modelo -fixed_fin <- plm(log_Deaths_per_100k ~ log_PrimPM25_per_100k + - log_NOx_per_100k + - log_VOC_per_100k + - log_gdppc_ppp_dol2011 + - log_flsp, - data = data, - index=c("country_name", "year"), model="within") -summary(fixed_fin) - - -# Fit by continent -fixed_fin_c1 <- plm(log_Deaths_per_100k ~ log_PrimPM25_per_100k + log_NOx_per_100k + log_gdppc_ppp_dol2011 + log_flsp, - data = data.c1, - index=c("country_name", "year"), model="within") -summary(fixed_fin_dev) - - - - -# ---------------------------------------------- -# ---------------------------------------------- -# PREDICTIONS -# ---------------------------------------------- -data.panel <- plm::pdata.frame(data, index = c("country_name", "year")) - -data$pred_log_Deaths_per_100k <- predict(fixed_fin, data.panel) - -# Calculate adders -data_adder <- data %>% - select(country_name, year, pop, log_Deaths_per_100k,pred_log_Deaths_per_100k) %>% - dplyr::filter(year == max(year)) %>% - mutate(Deaths_per_100k = exp(log_Deaths_per_100k), - pred_Deaths_per_100k = exp(pred_log_Deaths_per_100k)) %>% - mutate(bias.adder = Deaths_per_100k - pred_Deaths_per_100k) - -# Add GCAM_region -iso_GCAM_regID <- read.csv(paste0(datadir, "/iso_GCAM_regID.csv")) %>% - select(iso, GCAM_region) - -data_check <- data %>% - mutate(Deaths = exp(log_Deaths_per_100k) * pop / 100000) %>% - mutate(pred_Deaths = exp(pred_log_Deaths_per_100k) * pop / 100000) %>% - select(iso, country_name, year, pop, Deaths, pred_Deaths) %>% - left_join(iso_GCAM_regID, by = join_by(iso)) %>% - rename(country = country_name) %>% - mutate(Deaths_per_100K = Deaths / pop * 1E5, - pred_Deaths_per_100K = pred_Deaths / pop * 1E5) - - -ggplot(data_check, aes(x = Deaths_per_100K, y = pred_Deaths_per_100K, colour = GCAM_region)) + - geom_point(size = 1) + - geom_abline(intercept = 0, - slope = 1) + - theme_bw() + - theme(legend.position = "none") + - xlim(0,1000) + - ylim(0,1000) - -# Adjustment: Calculate a bias adder as the difference between predicted and observed deaths in final year -data_adder <- data_check %>% - filter(year == max(year)) %>% - mutate(bias.adder = Deaths_per_100K - pred_Deaths_per_100K) %>% - select(iso, country, bias.adder) - -data_check_fin <- data_check %>% - gcamdata::left_join_error_no_match(data_adder, by = join_by(iso, country)) %>% - mutate(pred_Deaths_per_100K_adj = pred_Deaths_per_100K + bias.adder, - pred_Deaths_adj = pred_Deaths_per_100K_adj * pop / 100000, - test = Deaths - pred_Deaths_adj) %>% - select(iso, country, year, Deaths, pred_Deaths = pred_Deaths_adj, GCAM_region) %>% - group_by(GCAM_region, year) %>% - summarise(Deaths = sum(Deaths), - pred_Deaths = sum(pred_Deaths)) %>% - ungroup() %>% - mutate(Deaths = round(Deaths, 0), - pred_Deaths = round(pred_Deaths, 0)) - -selected_regions <- c("Africa_Eastern", "China", "Europe_Eastern", "India", - "Indonesia", "South Asia", "EU-15", "USA") - -ggplot(data_check_fin %>% filter(GCAM_region %in% selected_regions), aes(x = Deaths, y = pred_Deaths, colour = GCAM_region)) + - geom_point(size = 1) + - geom_abline(intercept = 0, - slope = 1) + - theme_bw() + - theme(legend.position = "bottom") - - - - - - - +# final <- fit_final_model(data) +# lmtest::coeftest(final$model, vcov = final$vcov_dk) diff --git a/man/calc_hap_impacts.Rd b/man/calc_hap_impacts.Rd index 95391ee..dac87f7 100644 --- a/man/calc_hap_impacts.Rd +++ b/man/calc_hap_impacts.Rd @@ -17,7 +17,8 @@ calc_hap_impacts( anim = TRUE, HIA_var = "deaths", normalized = FALSE, - by_gr = FALSE + by_gr = FALSE, + fit_result = NULL ) } \arguments{ @@ -45,10 +46,32 @@ calc_hap_impacts( \item{normalized}{Transform the output to "normalized" values. By default = FALSE} -\item{by_gr}{Estimate damages at group level. Just for illustrative purposes. By default = FALSE} +\item{by_gr}{Report within-country income-decile group shares of +population, PM2.5+NOx exposure, and GDP as a CSV table plus a pie-chart +snapshot at final_db_year (output/by_gr/). Descriptive only -- it does +not estimate group-level health impacts, since the regression model is +fit on between-country variation and was never validated against +within-country/sub-national variation. By default = FALSE} + +\item{fit_result}{Optional pre-computed result from fit_model(HIA_var), +e.g. from a previous call, to reuse instead of refitting the national +regression from scratch. fit_model() doesn't depend on GCAM scenario +data at all, so refitting it (including its Driscoll-Kraay vcov, which +isn't cheap) on every calc_hap_impacts() call is pure waste when +comparing many scenarios for the same HIA_var. Must have been fit with +the same HIA_var as this call; a mismatch raises an error rather than +silently using the wrong model. By default = NULL (fits internally).} } \value{ -Health impacts attributable to HAP for all the selected years +Health impacts attributable to HAP for all the selected years, as + a tibble of scenario, country, year, pred_var, pred_value, + pred_value_normalized, reliability_ratio, and reliability. The last two + flag how much of a country's predicted level rests on fit_model()'s + bias-adder correction rather than on the regression's own covariates + (reliability_ratio = abs(bias.adder)/naive prediction; reliability is a + "high"/"medium"/"low" label over that ratio, thresholds 0.5 and 2) -- + "low" flags a country whose absolute-level prediction the model has + little independent basis for, not necessarily a wrong number. } \description{ Run the model to estimate the health impacts attributable to household air pollution for GCAM scenarios diff --git a/man/calc_hia_adder.Rd b/man/calc_hia_adder.Rd deleted file mode 100644 index ce1cad8..0000000 --- a/man/calc_hia_adder.Rd +++ /dev/null @@ -1,21 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/calc_hia_adder.R -\name{calc_hia_adder} -\alias{calc_hia_adder} -\title{calc_hia_adder} -\source{ -IHME -} -\usage{ -calc_hia_adder() -} -\value{ -Adder -} -\description{ -Calculate the adder for HAP impact analysis - -The adder is calculated as the difference between the observed and estimated values in final data year -} -\keyword{Bias} -\keyword{adder} diff --git a/man/fit_model.Rd b/man/fit_model.Rd index df936a2..19ff338 100644 --- a/man/fit_model.Rd +++ b/man/fit_model.Rd @@ -4,21 +4,106 @@ \alias{fit_model} \title{fit_model} \source{ -Details on plm estimation can be found here: https://cran.r-project.org/web/packages/plm/plm.pdf#' +Details on plm estimation: https://cran.r-project.org/web/packages/plm/plm.pdf + +Driscoll-Kraay SEs: plm::vcovSCC documentation } \usage{ -fit_model(HIA_var) +fit_model(HIA_var, n_years = 5) } \arguments{ \item{HIA_var}{Health metric to be predicted. c("deaths", "yll", "dalys")} + +\item{n_years}{Number of most recent historical years (per country) to +average over when computing the bias adder. Default 5 — long enough +to smooth out a single-year shock, short enough to reflect the +country's current level rather than its full historical average.} } \value{ -Regression model and coefficients for prediction +A list: model.fixed (plm object), predictable_regions (character + vector of countries in the estimation sample), vcov_dk (Driscoll-Kraay + variance-covariance matrix for inference — use this instead of the + default plm SEs for confidence intervals / significance testing on + this model's coefficients), bias_adder (tibble of country_name, + bias.adder, reliability_ratio — bias.adder is the per-country additive + level-calibration term, on the per_100k scale, to add to future + predictions [do not multiply]; reliability_ratio is + abs(bias.adder)/naive_prediction over the calibration window, where a + high value flags a country whose absolute-level prediction rests + almost entirely on the bias correction rather than the model's + covariates, and should be trusted less; replaces rhap::hia_adder), + HIA_var (the HIA_var this model was fit for, echoed back so a cached + fit_model() result can be validated against a different call site's + HIA_var before reuse). } \description{ -Fit the model using fixed effects model and plm function. +Fit the model using a fixed effects model (plm, entity within-estimator) +with a linear year trend, and return Driscoll-Kraay standard errors +alongside the fitted model. + +Model selection notes (see "./rhap/inst/extdata/model_testing.R" for the +full diagnostic history): +- Full year fixed effects (factor(year)) over-absorb common trend that + overlaps with the pollutant decline across countries, washing out + pollutant coefficients almost entirely. A linear year trend gives + coefficients consistent with a spline(df=4) and full year-FE + specification, so it was preferred for parsimony and because it + extrapolates sensibly to future scenario years (unlike year dummies). +- log_VOC_per_100k was dropped: coefficient is ~0 and non-significant + under every specification and SE estimator tested (naive, Arellano + clustered, Driscoll-Kraay), and does not survive a joint Wald test + either alone or jointly with PM2.5/floorspace. +- Under Driscoll-Kraay SEs (robust to both serial correlation and + cross-sectional dependence — appropriate for a country-year panel + with common shocks like oil prices, recessions, and multilateral + pollution treaties), only log_NOx_per_100k and log_gdppc_ppp_dol2011 + are individually significant. log_PrimPM25_per_100k and log_flsp are + correctly signed (positive / negative respectively, as expected) but + not individually or jointly significant (joint Wald test p = 0.35). + They are retained on theoretical grounds (omitted-variable-bias risk + from dropping known-relevant regressors) but their point estimates + should be treated as directionally informative only, not precise, + for any scenario-comparison use of this model's coefficients. + +- The country-level bias adder previously stored in the external + rhap::hia_adder table is now computed internally from the fitted + model itself (see below), so it can never drift out of sync with + the coefficients used to produce it. It's calculated the same way + predict() will later be used on scenario data — coefficients only, + no fixed effect — so the gap it captures is exactly the gap that + needs to be added back at prediction time (country fixed effect + + any residual model bias), calibrated to recent observed years + rather than a single year to avoid baking in a one-off shock. + +- The bias adder is ADDITIVE (level scale, per_100k units), not + multiplicative, despite log(Y) = X*beta + alpha_i implying alpha_i is + technically a multiplicative scalar exp(alpha_i) in levels. A + multiplicative geometric-mean correction was tried and is more + "correct" in that narrow sense, but performed far worse in practice: + validated against real GBD out-of-sample data, it collapsed + cross-country correlation from ~0.86 (additive) to ~0.18. The reason + is that ~22% of countries (e.g. Uganda) have a common linear year + trend that doesn't match their own trajectory, so the gap between + the trend-implied (fixed-effect-free) prediction and reality keeps + growing even within the training years -- sometimes to a fixed + effect worth 1000x+ in levels. A multiplicative correction amplifies + that growing gap catastrophically when applied to a slightly + different (scenario) year; an additive correction, bounded by the + actual scale of the outcome variable, does not. See + inst/extdata/model_testing.R for the comparison. -For the model dplyr::selection, testing and validation see "./rhap/inst/extdata/model_testing.R" +- IMPORTANT: the linear year trend is passed as an explicit numeric + column (year_num), not as.numeric(year) inline in the formula. + plm::pdata.frame() (used downstream at prediction time) converts + the "year" index column into a factor; as.numeric() on a factor + returns the factor's level codes, not the original calendar year, + and those codes are reassigned per-dataset — so a scenario dataset + spanning different years than the training data would get a + trend encoding disconnected from the one the coefficient was + estimated against. year_num is a plain (non-index) numeric column, + so it isn't touched by pdata.frame()'s factor coercion and stays + on a consistent calendar-year scale in both training and + prediction. } \keyword{Econometric} \keyword{effects} diff --git a/man/hia_adder.Rd b/man/hia_adder.Rd deleted file mode 100644 index 62d39a1..0000000 --- a/man/hia_adder.Rd +++ /dev/null @@ -1,25 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/data.R -\docType{data} -\name{hia_adder} -\alias{hia_adder} -\title{hia_adder} -\format{ -.csv -} -\source{ -Own calculations based on IHME data -} -\usage{ -hia_adder -} -\description{ -Bias adder, calculated as the difference between the observed and "predicted" values in final observed year -} -\examples{ -\dontrun{ -library(rhap) -rhap::hia_adder -} -} -\keyword{datasets} diff --git a/tests/testthat/helper-snapshots.R b/tests/testthat/helper-snapshots.R new file mode 100644 index 0000000..63f9df7 --- /dev/null +++ b/tests/testthat/helper-snapshots.R @@ -0,0 +1,17 @@ +# Load an .RData snapshot without clobbering same-named variables in the +# caller's environment. get(load(path)) is unsafe: load() injects the +# saved object under whatever name it was saved as directly into the +# calling environment as a side effect, independent of what the load() +# return value (the object's name) is assigned to. If a snapshot was +# saved under the same name as the local variable holding a freshly +# computed result (e.g. "testOutput"), load() silently overwrites that +# variable before the comparison runs -- turning expect_equal(testOutput, +# testResult) into a comparison of the stale snapshot against itself. +# This happened for several snapshots in this suite (calc_hap_impacts, +# calc_ResidEm_grp x2, fit_model), which meant those tests always passed +# regardless of what the function under test actually returned. +load_snapshot <- function(path) { + e <- new.env() + nm <- load(path, envir = e) + get(nm[1], envir = e) +} diff --git a/tests/testthat/test1-run.R b/tests/testthat/test1-run.R index 10de434..bbfc6c4 100644 --- a/tests/testthat/test1-run.R +++ b/tests/testthat/test1-run.R @@ -13,7 +13,7 @@ test_that("Download db, create project, and run", { saveOutput = TRUE, map = TRUE, anim = FALSE, normalized = FALSE, by_gr = FALSE ) - testResult <- get(load(file.path(rprojroot::find_root(rprojroot::is_testthat), "testOutputs/calc_hap_impacts_output1.RData"))) + testResult <- load_snapshot(file.path(rprojroot::find_root(rprojroot::is_testthat), "testOutputs/calc_hap_impacts_output1.RData")) testthat::expect_equal(nrow(testOutput), nrow(testResult)) # check saved files @@ -117,4 +117,61 @@ test_that("Download db, create project, and run", { ), "Error: The specified by_gr '0' is invalid. Accepted by_gr values are: TRUE, FALSE. Please rerun the calc_hap_impacts function with a valid by_gr value." ) + + # fit_result: a pre-fit model should be reusable across calls for the + # same HIA_var, and rejected for a mismatched one. + fr <- fit_model(HIA_var = "deaths") + testOutput_cached <- calc_hap_impacts( + prj_name = prj_name, scen_name = scen_name, + final_db_year = 2020, HIA_var = "deaths", + saveOutput = FALSE, map = FALSE, anim = FALSE, + normalized = FALSE, by_gr = FALSE, fit_result = fr + ) + expect_true(nrow(testOutput_cached) > 0) + expect_error( + calc_hap_impacts( + prj_name = prj_name, scen_name = scen_name, + final_db_year = 2020, HIA_var = "yll", + saveOutput = FALSE, map = FALSE, anim = FALSE, + normalized = FALSE, by_gr = FALSE, fit_result = fr + ), + "the supplied fit_result was fit for HIA_var = 'deaths', but this call uses HIA_var = 'yll'" + ) + + # reliability columns: present, and reliability is a valid label over + # reliability_ratio's thresholds. + expect_true(all(c("reliability_ratio", "reliability") %in% names(testOutput_cached))) + expect_true(all(testOutput_cached$reliability %in% c("high", "medium", "low"))) + expect_true(all(testOutput_cached$reliability_ratio >= 0)) +}) + +test_that("calc_hap_impacts by_gr = TRUE reports group shares, not group health impacts", { + prj_name <- file.path(rprojroot::find_root(rprojroot::is_testthat), "testInputs", "test_prj_v7p1.dat") + scen_name <- "Reference" + + calc_hap_impacts( + prj_name = prj_name, scen_name = scen_name, + final_db_year = 2020, HIA_var = "deaths", + saveOutput = TRUE, map = FALSE, anim = FALSE, + normalized = FALSE, by_gr = TRUE + ) + + shares_path <- file.path(rprojroot::find_root(rprojroot::is_testthat), "output/by_gr/shares_byGR.csv") + expect_true(file.exists(shares_path)) + + shares <- read.csv(shares_path) + expect_named(shares, c("scenario", "country", "year", "group", "metric", "share")) + expect_setequal(unique(shares$metric), c("population", "exposure", "gdp")) + expect_setequal(unique(shares$group), paste0("d", 1:10)) + + # Group shares within the same scenario/country/year/metric must sum to + # ~1 (they're shares of that country's own group breakdown, not a + # cross-country comparison) -- this is the direct check that the shares + # redesign (replacing a fabricated group-level health estimate) is + # internally consistent. + totals <- stats::aggregate(share ~ scenario + country + year + metric, data = shares, FUN = sum) + expect_true(all(abs(totals$share - 1) < 0.01)) + + pie_path <- file.path(rprojroot::find_root(rprojroot::is_testthat), "output/by_gr/shares_byGR_2020.png") + expect_true(file.exists(pie_path)) }) diff --git a/tests/testthat/test2-calc_ResidEm_grp.R b/tests/testthat/test2-calc_ResidEm_grp.R index a03ad29..88fdb33 100644 --- a/tests/testthat/test2-calc_ResidEm_grp.R +++ b/tests/testthat/test2-calc_ResidEm_grp.R @@ -15,7 +15,7 @@ test_that("Create pie charts", { region = "EU-15", pie = TRUE ) - testResult <- get(load(file.path(rprojroot::find_root(rprojroot::is_testthat), "testOutputs/calc_ResidEm_grp_output_eu15.RData"))) + testResult <- load_snapshot(file.path(rprojroot::find_root(rprojroot::is_testthat), "testOutputs/calc_ResidEm_grp_output_eu15.RData")) testthat::expect_equal(testOutput, testResult) testOutput <- calc_ResidEm_grp( @@ -26,7 +26,7 @@ test_that("Create pie charts", { region = "EU-27", pie = TRUE ) - testResult <- get(load(file.path(rprojroot::find_root(rprojroot::is_testthat), "testOutputs/calc_ResidEm_grp_output_eu27.RData"))) + testResult <- load_snapshot(file.path(rprojroot::find_root(rprojroot::is_testthat), "testOutputs/calc_ResidEm_grp_output_eu27.RData")) testthat::expect_equal(testOutput, testResult) diff --git a/tests/testthat/test3-other.R b/tests/testthat/test3-other.R index acb402b..9cf9af0 100644 --- a/tests/testthat/test3-other.R +++ b/tests/testthat/test3-other.R @@ -4,7 +4,7 @@ library(magrittr) test_that("get_gdp_ctry", { testOutput <- get_gdp_ctry(ssp = "SSP3") - testResult <- get(load(file.path(rprojroot::find_root(rprojroot::is_testthat), "testOutputs/get_gdp_ctry.RData"))) + testResult <- load_snapshot(file.path(rprojroot::find_root(rprojroot::is_testthat), "testOutputs/get_gdp_ctry.RData")) testthat::expect_equal(testOutput, testResult) # error messages @@ -18,7 +18,7 @@ test_that("get_gdp_ctry", { test_that("get_pop_ctry", { testOutput <- get_pop_ctry(ssp = "SSP5") - testResult <- get(load(file.path(rprojroot::find_root(rprojroot::is_testthat), "testOutputs/get_pop_ctry.RData"))) + testResult <- load_snapshot(file.path(rprojroot::find_root(rprojroot::is_testthat), "testOutputs/get_pop_ctry.RData")) testthat::expect_equal(testOutput, testResult) # error messages @@ -31,15 +31,51 @@ test_that("get_pop_ctry", { test_that("create_panel", { testOutput <- create_panel() - testResult <- get(load(file.path(rprojroot::find_root(rprojroot::is_testthat), "testOutputs/create_panel.RData"))) + testResult <- load_snapshot(file.path(rprojroot::find_root(rprojroot::is_testthat), "testOutputs/create_panel.RData")) testthat::expect_equal(testOutput, testResult) }) test_that("fit_model", { - testOutput <- fit_model(HIA_var = "deaths") - testResult <- get(load(file.path(rprojroot::find_root(rprojroot::is_testthat), "testOutputs/fit_model.RData"))) - testthat::expect_equal(testOutput, testResult) + # A full-object snapshot comparison (the previous version of this test) + # is the wrong tool here: fit_model() returns a plm model object plus a + # Driscoll-Kraay vcov matrix, both of which are brittle to compare + # exactly across R/plm versions or minor upstream data refreshes. It's + # also how this test silently stopped checking anything at all: the old + # snapshot file stored its object under the name "testOutput", which + # collided with the local variable holding the freshly-computed result + # -- load()'s side effect overwrote it before expect_equal() ever ran, + # so the test was comparing the stale snapshot to itself regardless of + # what fit_model() actually returned. Structural and sign checks below + # are more robust and substantively more informative. + result <- fit_model(HIA_var = "deaths") + + expect_named(result, c("model.fixed", "predictable_regions", "vcov_dk", "bias_adder", "HIA_var")) + expect_s3_class(result$model.fixed, "plm") + expect_true(is.character(result$predictable_regions)) + expect_true(length(result$predictable_regions) > 0) + expect_true(is.matrix(result$vcov_dk)) + expect_equal(result$HIA_var, "deaths") + expect_named(result$bias_adder, c("country_name", "bias.adder", "reliability_ratio")) + expect_true(all(result$bias_adder$reliability_ratio >= 0)) + + # Coefficient signs must match the epidemiological/economic prior the + # model was built to satisfy: more emissions -> more health impacts, + # more GDP/floorspace -> fewer health impacts (see fit_model.R's roxygen + # notes and the model-validation vignette for why this specification + # was chosen over alternatives, e.g. the rejected GDP x PM2.5 + # interaction, which violated this exact check). + coefs <- stats::coef(result$model.fixed) + expect_gt(coefs[["log_PrimPM25_per_100k"]], 0) + expect_gt(coefs[["log_NOx_per_100k"]], 0) + expect_lt(coefs[["log_gdppc_ppp_dol2011"]], 0) + expect_lt(coefs[["log_flsp"]], 0) + + # n_years controls the bias_adder calibration window; a shorter window + # should still return one bias.adder per predictable country, just + # calibrated over fewer years. + result_n3 <- fit_model(HIA_var = "deaths", n_years = 3) + expect_equal(nrow(result_n3$bias_adder), nrow(result$bias_adder)) # error messages expect_error( diff --git a/tests/testthat/testOutputs/Reference_HAP_deaths.csv b/tests/testthat/testOutputs/Reference_HAP_deaths.csv index de7f227..950eaed 100644 --- a/tests/testthat/testOutputs/Reference_HAP_deaths.csv +++ b/tests/testthat/testOutputs/Reference_HAP_deaths.csv @@ -1,649 +1,645 @@ -"scenario","country","year","pred_var","pred_value","pred_value_normalized" -"Reference","Afghanistan",2020,"deaths",22609,58.99 -"Reference","Afghanistan",2030,"deaths",17021,33.77 -"Reference","Afghanistan",2040,"deaths",12716,20.12 -"Reference","Afghanistan",2050,"deaths",10216,13.51 -"Reference","Albania",2020,"deaths",1211,42.18 -"Reference","Albania",2030,"deaths",891,32.82 -"Reference","Albania",2040,"deaths",724,28.93 -"Reference","Albania",2050,"deaths",613,27.05 -"Reference","Algeria",2020,"deaths",0,0 -"Reference","Algeria",2030,"deaths",0,0 -"Reference","Algeria",2040,"deaths",0,0 -"Reference","Algeria",2050,"deaths",0,0 -"Reference","Angola",2020,"deaths",12508,38.02 -"Reference","Angola",2030,"deaths",14424,32.7 -"Reference","Angola",2040,"deaths",14101,24.98 -"Reference","Angola",2050,"deaths",13302,19.28 -"Reference","Argentina",2020,"deaths",1091,2.43 -"Reference","Argentina",2030,"deaths",567,1.2 -"Reference","Argentina",2040,"deaths",232,0.47 -"Reference","Argentina",2050,"deaths",0,0 -"Reference","Armenia",2020,"deaths",21,0.74 -"Reference","Armenia",2030,"deaths",0,0 -"Reference","Armenia",2040,"deaths",0,0 -"Reference","Armenia",2050,"deaths",0,0 -"Reference","Australia",2020,"deaths",26,0.1 -"Reference","Australia",2030,"deaths",14,0.05 -"Reference","Australia",2040,"deaths",6,0.02 -"Reference","Australia",2050,"deaths",0,0 -"Reference","Austria",2020,"deaths",8,0.09 -"Reference","Austria",2030,"deaths",4,0.04 -"Reference","Austria",2040,"deaths",1,0.01 -"Reference","Austria",2050,"deaths",0,0 -"Reference","Azerbaijan",2020,"deaths",827,8.05 -"Reference","Azerbaijan",2030,"deaths",564,5.26 -"Reference","Azerbaijan",2040,"deaths",201,1.81 -"Reference","Azerbaijan",2050,"deaths",0,0 -"Reference","Bahrain",2020,"deaths",0,0 -"Reference","Bahrain",2030,"deaths",0,0 -"Reference","Bahrain",2040,"deaths",0,0 -"Reference","Bahrain",2050,"deaths",0,0 -"Reference","Bangladesh",2020,"deaths",128565,77.25 -"Reference","Bangladesh",2030,"deaths",116222,63.38 -"Reference","Bangladesh",2040,"deaths",114439,58.94 -"Reference","Bangladesh",2050,"deaths",114269,57.34 -"Reference","Barbados",2020,"deaths",0,0 -"Reference","Barbados",2030,"deaths",0,0 -"Reference","Barbados",2040,"deaths",0,0 -"Reference","Barbados",2050,"deaths",0,0 -"Reference","Belarus",2020,"deaths",0,0 -"Reference","Belarus",2030,"deaths",0,0 -"Reference","Belarus",2040,"deaths",0,0 -"Reference","Belarus",2050,"deaths",0,0 -"Reference","Belgium",2020,"deaths",0,0 -"Reference","Belgium",2030,"deaths",0,0 -"Reference","Belgium",2040,"deaths",0,0 -"Reference","Belgium",2050,"deaths",0,0 -"Reference","Benin",2020,"deaths",21154,169.68 -"Reference","Benin",2030,"deaths",21150,130.65 -"Reference","Benin",2040,"deaths",22184,110.09 -"Reference","Benin",2050,"deaths",24255,100.54 -"Reference","Bolivia (Plurinational State of)",2020,"deaths",3520,29.67 -"Reference","Bolivia (Plurinational State of)",2030,"deaths",3229,24.02 -"Reference","Bolivia (Plurinational State of)",2040,"deaths",2738,18.39 -"Reference","Bolivia (Plurinational State of)",2050,"deaths",2308,14.46 -"Reference","Bosnia and Herzegovina",2020,"deaths",652,19.53 -"Reference","Bosnia and Herzegovina",2030,"deaths",179,5.86 -"Reference","Bosnia and Herzegovina",2040,"deaths",7,0.24 -"Reference","Bosnia and Herzegovina",2050,"deaths",0,0 -"Reference","Botswana",2020,"deaths",0,0 -"Reference","Botswana",2030,"deaths",0,0 -"Reference","Botswana",2040,"deaths",0,0 -"Reference","Botswana",2050,"deaths",0,0 -"Reference","Brazil",2020,"deaths",0,0 -"Reference","Brazil",2030,"deaths",0,0 -"Reference","Brazil",2040,"deaths",0,0 -"Reference","Brazil",2050,"deaths",0,0 -"Reference","Bulgaria",2020,"deaths",0,0 -"Reference","Bulgaria",2030,"deaths",0,0 -"Reference","Bulgaria",2040,"deaths",0,0 -"Reference","Bulgaria",2050,"deaths",0,0 -"Reference","Burkina Faso",2020,"deaths",41438,195.17 -"Reference","Burkina Faso",2030,"deaths",45046,164.68 -"Reference","Burkina Faso",2040,"deaths",49273,144.38 -"Reference","Burkina Faso",2050,"deaths",55343,135.01 -"Reference","Burundi",2020,"deaths",29041,240.93 -"Reference","Burundi",2030,"deaths",27534,179.43 -"Reference","Burundi",2040,"deaths",25648,135.81 -"Reference","Burundi",2050,"deaths",25347,113.6 -"Reference","Cambodia",2020,"deaths",19934,122.32 -"Reference","Cambodia",2030,"deaths",18835,103.5 -"Reference","Cambodia",2040,"deaths",19062,97.17 -"Reference","Cambodia",2050,"deaths",19206,94.18 -"Reference","Cameroon",2020,"deaths",24300,92.97 -"Reference","Cameroon",2030,"deaths",22719,68.43 -"Reference","Cameroon",2040,"deaths",23207,57.55 -"Reference","Cameroon",2050,"deaths",24914,52.89 -"Reference","Canada",2020,"deaths",0,0 -"Reference","Canada",2030,"deaths",0,0 -"Reference","Canada",2040,"deaths",0,0 -"Reference","Canada",2050,"deaths",0,0 -"Reference","Cabo Verde",2020,"deaths",99,17.11 -"Reference","Cabo Verde",2030,"deaths",108,17.09 -"Reference","Cabo Verde",2040,"deaths",115,17.07 -"Reference","Cabo Verde",2050,"deaths",118,17.05 -"Reference","Central African Republic",2020,"deaths",16292,309.03 -"Reference","Central African Republic",2030,"deaths",19379,275.26 -"Reference","Central African Republic",2040,"deaths",18740,206.92 -"Reference","Central African Republic",2050,"deaths",18756,169.22 -"Reference","Chad",2020,"deaths",35340,215.76 -"Reference","Chad",2030,"deaths",44064,197.46 -"Reference","Chad",2040,"deaths",47431,162.24 -"Reference","Chad",2050,"deaths",52722,144.01 -"Reference","Chile",2020,"deaths",882,4.59 -"Reference","Chile",2030,"deaths",714,3.58 -"Reference","Chile",2040,"deaths",581,2.85 -"Reference","Chile",2050,"deaths",478,2.34 -"Reference","China",2020,"deaths",0,0 -"Reference","China",2030,"deaths",0,0 -"Reference","China",2040,"deaths",0,0 -"Reference","China",2050,"deaths",0,0 -"Reference","Colombia",2020,"deaths",0,0 -"Reference","Colombia",2030,"deaths",0,0 -"Reference","Colombia",2040,"deaths",0,0 -"Reference","Colombia",2050,"deaths",0,0 -"Reference","Comoros",2020,"deaths",586,73.35 -"Reference","Comoros",2030,"deaths",527,54.71 -"Reference","Comoros",2040,"deaths",413,36.88 -"Reference","Comoros",2050,"deaths",346,27.42 -"Reference","Congo",2020,"deaths",7521,133.47 -"Reference","Congo",2030,"deaths",7271,102.79 -"Reference","Congo",2040,"deaths",6441,75.24 -"Reference","Congo",2050,"deaths",5566,55.93 -"Reference","Democratic Republic of the Congo",2020,"deaths",56087,61.41 -"Reference","Democratic Republic of the Congo",2030,"deaths",75016,59.75 -"Reference","Democratic Republic of the Congo",2040,"deaths",97214,58.79 -"Reference","Democratic Republic of the Congo",2050,"deaths",121548,58.32 -"Reference","Costa Rica",2020,"deaths",0,0 -"Reference","Costa Rica",2030,"deaths",0,0 -"Reference","Costa Rica",2040,"deaths",0,0 -"Reference","Costa Rica",2050,"deaths",0,0 -"Reference","Côte d'Ivoire",2020,"deaths",4806,18.15 -"Reference","Côte d'Ivoire",2030,"deaths",4390,13.27 -"Reference","Côte d'Ivoire",2040,"deaths",4105,10.38 -"Reference","Côte d'Ivoire",2050,"deaths",4207,9.28 -"Reference","Croatia",2020,"deaths",0,0 -"Reference","Croatia",2030,"deaths",0,0 -"Reference","Croatia",2040,"deaths",0,0 -"Reference","Croatia",2050,"deaths",0,0 -"Reference","Cuba",2020,"deaths",0,0 -"Reference","Cuba",2030,"deaths",0,0 -"Reference","Cuba",2040,"deaths",0,0 -"Reference","Cuba",2050,"deaths",0,0 -"Reference","Cyprus",2020,"deaths",608,49.24 -"Reference","Cyprus",2030,"deaths",304,23.1 -"Reference","Cyprus",2040,"deaths",153,11.21 -"Reference","Cyprus",2050,"deaths",83,5.85 -"Reference","Czechia",2020,"deaths",25088,238.07 -"Reference","Czechia",2030,"deaths",11760,109.78 -"Reference","Czechia",2040,"deaths",4883,46.87 -"Reference","Czechia",2050,"deaths",2172,21.2 -"Reference","Denmark",2020,"deaths",3452,59.4 -"Reference","Denmark",2030,"deaths",2031,34.05 -"Reference","Denmark",2040,"deaths",1097,18.08 -"Reference","Denmark",2050,"deaths",539,8.74 -"Reference","Djibouti",2020,"deaths",0,0 -"Reference","Djibouti",2030,"deaths",0,0 -"Reference","Djibouti",2040,"deaths",0,0 -"Reference","Djibouti",2050,"deaths",0,0 -"Reference","Dominican Republic",2020,"deaths",363,3.32 -"Reference","Dominican Republic",2030,"deaths",0,0 -"Reference","Dominican Republic",2040,"deaths",0,0 -"Reference","Dominican Republic",2050,"deaths",0,0 -"Reference","Ecuador",2020,"deaths",815,4.66 -"Reference","Ecuador",2030,"deaths",0,0 -"Reference","Ecuador",2040,"deaths",0,0 -"Reference","Ecuador",2050,"deaths",0,0 -"Reference","Egypt",2020,"deaths",2003,1.88 -"Reference","Egypt",2030,"deaths",251,0.2 -"Reference","Egypt",2040,"deaths",0,0 -"Reference","Egypt",2050,"deaths",0,0 -"Reference","El Salvador",2020,"deaths",628,10 -"Reference","El Salvador",2030,"deaths",630,9.76 -"Reference","El Salvador",2040,"deaths",622,9.59 -"Reference","El Salvador",2050,"deaths",608,9.54 -"Reference","Equatorial Guinea",2020,"deaths",81,5.14 -"Reference","Equatorial Guinea",2030,"deaths",0,0 -"Reference","Equatorial Guinea",2040,"deaths",0,0 -"Reference","Equatorial Guinea",2050,"deaths",0,0 -"Reference","Estonia",2020,"deaths",596,44.83 -"Reference","Estonia",2030,"deaths",887,67.62 -"Reference","Estonia",2040,"deaths",543,43.32 -"Reference","Estonia",2050,"deaths",291,23.97 -"Reference","Ethiopia",2020,"deaths",55368,47.88 -"Reference","Ethiopia",2030,"deaths",67170,45.77 -"Reference","Ethiopia",2040,"deaths",79634,44.59 -"Reference","Ethiopia",2050,"deaths",91858,44 -"Reference","Finland",2020,"deaths",1708,30.92 -"Reference","Finland",2030,"deaths",614,11.14 -"Reference","Finland",2040,"deaths",88,1.61 -"Reference","Finland",2050,"deaths",0,0 -"Reference","France",2020,"deaths",0,0 -"Reference","France",2030,"deaths",0,0 -"Reference","France",2040,"deaths",0,0 -"Reference","France",2050,"deaths",0,0 -"Reference","Gabon",2020,"deaths",0,0 -"Reference","Gabon",2030,"deaths",0,0 -"Reference","Gabon",2040,"deaths",0,0 -"Reference","Gabon",2050,"deaths",0,0 -"Reference","Gambia",2020,"deaths",1525,60 -"Reference","Gambia",2030,"deaths",1710,54.07 -"Reference","Gambia",2040,"deaths",1821,47.81 -"Reference","Gambia",2050,"deaths",1939,44.14 -"Reference","Georgia",2020,"deaths",3624,96.16 -"Reference","Georgia",2030,"deaths",2482,71.59 -"Reference","Georgia",2040,"deaths",1819,57.26 -"Reference","Georgia",2050,"deaths",1469,50.37 -"Reference","Germany",2020,"deaths",26754,32.13 -"Reference","Germany",2030,"deaths",10867,13.09 -"Reference","Germany",2040,"deaths",3882,4.74 -"Reference","Germany",2050,"deaths",1522,1.88 -"Reference","Ghana",2020,"deaths",8918,28 -"Reference","Ghana",2030,"deaths",10853,27.98 -"Reference","Ghana",2040,"deaths",12826,27.97 -"Reference","Ghana",2050,"deaths",14755,27.97 -"Reference","Greece",2020,"deaths",2717,25.77 -"Reference","Greece",2030,"deaths",1636,16.03 -"Reference","Greece",2040,"deaths",887,8.9 -"Reference","Greece",2050,"deaths",497,5.08 -"Reference","Guatemala",2020,"deaths",0,0 -"Reference","Guatemala",2030,"deaths",0,0 -"Reference","Guatemala",2040,"deaths",0,0 -"Reference","Guatemala",2050,"deaths",0,0 -"Reference","Guinea",2020,"deaths",12200,93.55 -"Reference","Guinea",2030,"deaths",12459,76.77 -"Reference","Guinea",2040,"deaths",12741,65.85 -"Reference","Guinea",2050,"deaths",13348,60.15 -"Reference","Guinea-Bissau",2020,"deaths",2929,146.96 -"Reference","Guinea-Bissau",2030,"deaths",2983,122.42 -"Reference","Guinea-Bissau",2040,"deaths",3052,106.93 -"Reference","Guinea-Bissau",2050,"deaths",3203,99.78 -"Reference","Haiti",2020,"deaths",7140,63.55 -"Reference","Haiti",2030,"deaths",4133,32.35 -"Reference","Haiti",2040,"deaths",1967,14 -"Reference","Haiti",2050,"deaths",526,3.52 -"Reference","Honduras",2020,"deaths",1298,12.93 -"Reference","Honduras",2030,"deaths",539,4.65 -"Reference","Honduras",2040,"deaths",0,0 -"Reference","Honduras",2050,"deaths",0,0 -"Reference","Hungary",2020,"deaths",1262,12.92 -"Reference","Hungary",2030,"deaths",0,0 -"Reference","Hungary",2040,"deaths",0,0 -"Reference","Hungary",2050,"deaths",0,0 -"Reference","Iceland",2020,"deaths",35,9.71 -"Reference","Iceland",2030,"deaths",20,5.11 -"Reference","Iceland",2040,"deaths",12,3.04 -"Reference","Iceland",2050,"deaths",8,1.87 -"Reference","India",2020,"deaths",630072,45.33 -"Reference","India",2030,"deaths",680282,45.3 -"Reference","India",2040,"deaths",717252,45.29 -"Reference","India",2050,"deaths",733698,45.28 -"Reference","Indonesia",2020,"deaths",145542,53.74 -"Reference","Indonesia",2030,"deaths",88133,30.33 -"Reference","Indonesia",2040,"deaths",70157,23.03 -"Reference","Indonesia",2050,"deaths",61188,19.7 -"Reference","Iran (Islamic Republic of)",2020,"deaths",1227,1.41 -"Reference","Iran (Islamic Republic of)",2030,"deaths",0,0 -"Reference","Iran (Islamic Republic of)",2040,"deaths",0,0 -"Reference","Iran (Islamic Republic of)",2050,"deaths",0,0 -"Reference","Iraq",2020,"deaths",143,0.34 -"Reference","Iraq",2030,"deaths",57,0.11 -"Reference","Iraq",2040,"deaths",0,0 -"Reference","Iraq",2050,"deaths",0,0 -"Reference","Ireland",2020,"deaths",0,0 -"Reference","Ireland",2030,"deaths",0,0 -"Reference","Ireland",2040,"deaths",0,0 -"Reference","Ireland",2050,"deaths",0,0 -"Reference","Israel",2020,"deaths",3,0.04 -"Reference","Israel",2030,"deaths",2,0.02 -"Reference","Israel",2040,"deaths",1,0.01 -"Reference","Israel",2050,"deaths",0,0 -"Reference","Italy",2020,"deaths",125,0.21 -"Reference","Italy",2030,"deaths",105,0.18 -"Reference","Italy",2040,"deaths",97,0.17 -"Reference","Italy",2050,"deaths",89,0.16 -"Reference","Jamaica",2020,"deaths",190,6.75 -"Reference","Jamaica",2030,"deaths",188,6.59 -"Reference","Jamaica",2040,"deaths",185,6.5 -"Reference","Jamaica",2050,"deaths",179,6.43 -"Reference","Japan",2020,"deaths",20978,16.71 -"Reference","Japan",2030,"deaths",12520,10.49 -"Reference","Japan",2040,"deaths",6613,5.87 -"Reference","Japan",2050,"deaths",3020,2.85 -"Reference","Jordan",2020,"deaths",0,0 -"Reference","Jordan",2030,"deaths",0,0 -"Reference","Jordan",2040,"deaths",0,0 -"Reference","Jordan",2050,"deaths",0,0 -"Reference","Kazakhstan",2020,"deaths",1539,8.16 -"Reference","Kazakhstan",2030,"deaths",1680,8.14 -"Reference","Kazakhstan",2040,"deaths",1794,8.12 -"Reference","Kazakhstan",2050,"deaths",1893,8.11 -"Reference","Kenya",2020,"deaths",25303,49.17 -"Reference","Kenya",2030,"deaths",27432,43.92 -"Reference","Kenya",2040,"deaths",28002,38.14 -"Reference","Kenya",2050,"deaths",29192,35.15 -"Reference","Democratic People's Republic of Korea",2020,"deaths",39567,153.28 -"Reference","Democratic People's Republic of Korea",2030,"deaths",36999,139.45 -"Reference","Democratic People's Republic of Korea",2040,"deaths",34561,132.53 -"Reference","Democratic People's Republic of Korea",2050,"deaths",32507,128.94 -"Reference","Republic of Korea",2020,"deaths",0,0 -"Reference","Republic of Korea",2030,"deaths",0,0 -"Reference","Republic of Korea",2040,"deaths",0,0 -"Reference","Republic of Korea",2050,"deaths",0,0 -"Reference","Kuwait",2020,"deaths",0,0 -"Reference","Kuwait",2030,"deaths",0,0 -"Reference","Kuwait",2040,"deaths",0,0 -"Reference","Kuwait",2050,"deaths",0,0 -"Reference","Kyrgyzstan",2020,"deaths",2169,34.03 -"Reference","Kyrgyzstan",2030,"deaths",1176,15.95 -"Reference","Kyrgyzstan",2040,"deaths",383,4.64 -"Reference","Kyrgyzstan",2050,"deaths",0,0 -"Reference","Lao People's Democratic Republic",2020,"deaths",6253,86.06 -"Reference","Lao People's Democratic Republic",2030,"deaths",6895,83.92 -"Reference","Lao People's Democratic Republic",2040,"deaths",7401,82.96 -"Reference","Lao People's Democratic Republic",2050,"deaths",7688,82.42 -"Reference","Latvia",2020,"deaths",108,5.64 -"Reference","Latvia",2030,"deaths",97,5.61 -"Reference","Latvia",2040,"deaths",85,5.44 -"Reference","Latvia",2050,"deaths",77,5.32 -"Reference","Lebanon",2020,"deaths",4923,86.45 -"Reference","Lebanon",2030,"deaths",3865,61.97 -"Reference","Lebanon",2040,"deaths",2307,34.16 -"Reference","Lebanon",2050,"deaths",1301,18.27 -"Reference","Lesotho",2020,"deaths",2944,131.45 -"Reference","Lesotho",2030,"deaths",2628,107.46 -"Reference","Lesotho",2040,"deaths",2467,94.96 -"Reference","Lesotho",2050,"deaths",2370,88.3 -"Reference","Liberia",2020,"deaths",0,0 -"Reference","Liberia",2030,"deaths",0,0 -"Reference","Liberia",2040,"deaths",0,0 -"Reference","Liberia",2050,"deaths",0,0 -"Reference","Libya",2020,"deaths",0,0 -"Reference","Libya",2030,"deaths",0,0 -"Reference","Libya",2040,"deaths",0,0 -"Reference","Libya",2050,"deaths",0,0 -"Reference","Lithuania",2020,"deaths",0,0 -"Reference","Lithuania",2030,"deaths",0,0 -"Reference","Lithuania",2040,"deaths",0,0 -"Reference","Lithuania",2050,"deaths",0,0 -"Reference","Luxembourg",2020,"deaths",144,23.01 -"Reference","Luxembourg",2030,"deaths",0,0 -"Reference","Luxembourg",2040,"deaths",0,0 -"Reference","Luxembourg",2050,"deaths",0,0 -"Reference","North Macedonia",2020,"deaths",1774,83.93 -"Reference","North Macedonia",2030,"deaths",1349,64.17 -"Reference","North Macedonia",2040,"deaths",848,40.62 -"Reference","North Macedonia",2050,"deaths",570,27.87 -"Reference","Madagascar",2020,"deaths",14985,53.75 -"Reference","Madagascar",2030,"deaths",18702,53.47 -"Reference","Madagascar",2040,"deaths",22473,53.36 -"Reference","Madagascar",2050,"deaths",26074,53.31 -"Reference","Malawi",2020,"deaths",24851,129.97 -"Reference","Malawi",2030,"deaths",26520,108.09 -"Reference","Malawi",2040,"deaths",27042,88.71 -"Reference","Malawi",2050,"deaths",28595,78.85 -"Reference","Malaysia",2020,"deaths",0,0 -"Reference","Malaysia",2030,"deaths",0,0 -"Reference","Malaysia",2040,"deaths",0,0 -"Reference","Malaysia",2050,"deaths",0,0 -"Reference","Mali",2020,"deaths",25210,120.7 -"Reference","Mali",2030,"deaths",31866,113.06 -"Reference","Mali",2040,"deaths",39327,108.18 -"Reference","Mali",2050,"deaths",46930,105.51 -"Reference","Malta",2020,"deaths",0,0 -"Reference","Malta",2030,"deaths",0,0 -"Reference","Malta",2040,"deaths",0,0 -"Reference","Malta",2050,"deaths",0,0 -"Reference","Mauritania",2020,"deaths",0,0 -"Reference","Mauritania",2030,"deaths",0,0 -"Reference","Mauritania",2040,"deaths",0,0 -"Reference","Mauritania",2050,"deaths",0,0 -"Reference","Mauritius",2020,"deaths",268,20.64 -"Reference","Mauritius",2030,"deaths",196,15.1 -"Reference","Mauritius",2040,"deaths",104,8.1 -"Reference","Mauritius",2050,"deaths",56,4.53 -"Reference","Mexico",2020,"deaths",22949,18.27 -"Reference","Mexico",2030,"deaths",16206,11.95 -"Reference","Mexico",2040,"deaths",13920,9.66 -"Reference","Mexico",2050,"deaths",12671,8.48 -"Reference","Republic of Moldova",2020,"deaths",456,14.73 -"Reference","Republic of Moldova",2030,"deaths",344,11.74 -"Reference","Republic of Moldova",2040,"deaths",244,9.08 -"Reference","Republic of Moldova",2050,"deaths",206,8.27 -"Reference","Mongolia",2020,"deaths",3980,121.84 -"Reference","Mongolia",2030,"deaths",1764,47.37 -"Reference","Mongolia",2040,"deaths",371,8.95 -"Reference","Mongolia",2050,"deaths",0,0 -"Reference","Montenegro",2020,"deaths",0,0 -"Reference","Montenegro",2030,"deaths",0,0 -"Reference","Montenegro",2040,"deaths",0,0 -"Reference","Montenegro",2050,"deaths",0,0 -"Reference","Morocco",2020,"deaths",0,0 -"Reference","Morocco",2030,"deaths",0,0 -"Reference","Morocco",2040,"deaths",0,0 -"Reference","Morocco",2050,"deaths",0,0 -"Reference","Mozambique",2020,"deaths",29584,96.3 -"Reference","Mozambique",2030,"deaths",30469,76.98 -"Reference","Mozambique",2040,"deaths",33088,67.76 -"Reference","Mozambique",2050,"deaths",36860,63.87 -"Reference","Myanmar",2020,"deaths",43083,80.94 -"Reference","Myanmar",2030,"deaths",45547,80.91 -"Reference","Myanmar",2040,"deaths",46481,80.9 -"Reference","Myanmar",2050,"deaths",46222,80.89 -"Reference","Namibia",2020,"deaths",0,0 -"Reference","Namibia",2030,"deaths",0,0 -"Reference","Namibia",2040,"deaths",0,0 -"Reference","Namibia",2050,"deaths",0,0 -"Reference","Nepal",2020,"deaths",25775,88.88 -"Reference","Nepal",2030,"deaths",25480,79.72 -"Reference","Nepal",2040,"deaths",24947,73.77 -"Reference","Nepal",2050,"deaths",24685,70.75 -"Reference","Netherlands",2020,"deaths",15532,89.25 -"Reference","Netherlands",2030,"deaths",8410,47.47 -"Reference","Netherlands",2040,"deaths",3942,22.08 -"Reference","Netherlands",2050,"deaths",1838,10.29 -"Reference","New Zealand",2020,"deaths",2556,50.85 -"Reference","New Zealand",2030,"deaths",1911,35.91 -"Reference","New Zealand",2040,"deaths",1123,20.18 -"Reference","New Zealand",2050,"deaths",621,10.76 -"Reference","Nicaragua",2020,"deaths",2405,35.84 -"Reference","Nicaragua",2030,"deaths",2413,31.9 -"Reference","Nicaragua",2040,"deaths",2447,29.83 -"Reference","Nicaragua",2050,"deaths",2463,28.65 -"Reference","Niger",2020,"deaths",16882,70.69 -"Reference","Niger",2030,"deaths",24278,70.67 -"Reference","Niger",2040,"deaths",33820,70.66 -"Reference","Niger",2050,"deaths",45096,70.66 -"Reference","Nigeria",2020,"deaths",83259,40.46 -"Reference","Nigeria",2030,"deaths",103749,39.73 -"Reference","Nigeria",2040,"deaths",127882,39.09 -"Reference","Nigeria",2050,"deaths",155031,38.7 -"Reference","Norway",2020,"deaths",642,11.96 -"Reference","Norway",2030,"deaths",0,0 -"Reference","Norway",2040,"deaths",0,0 -"Reference","Norway",2050,"deaths",0,0 -"Reference","Oman",2020,"deaths",0,0 -"Reference","Oman",2030,"deaths",0,0 -"Reference","Oman",2040,"deaths",0,0 -"Reference","Oman",2050,"deaths",0,0 -"Reference","Pakistan",2020,"deaths",109855,48.8 -"Reference","Pakistan",2030,"deaths",127178,46.48 -"Reference","Pakistan",2040,"deaths",147727,45.65 -"Reference","Pakistan",2050,"deaths",167570,45.23 -"Reference","Palestine",2020,"deaths",357,7.19 -"Reference","Palestine",2030,"deaths",242,3.93 -"Reference","Palestine",2040,"deaths",105,1.45 -"Reference","Palestine",2050,"deaths",0,0 -"Reference","Panama",2020,"deaths",0,0 -"Reference","Panama",2030,"deaths",0,0 -"Reference","Panama",2040,"deaths",0,0 -"Reference","Panama",2050,"deaths",0,0 -"Reference","Paraguay",2020,"deaths",1829,27.82 -"Reference","Paraguay",2030,"deaths",1087,14.64 -"Reference","Paraguay",2040,"deaths",784,9.67 -"Reference","Paraguay",2050,"deaths",649,7.58 -"Reference","Peru",2020,"deaths",2502,7.56 -"Reference","Peru",2030,"deaths",1695,4.71 -"Reference","Peru",2040,"deaths",1317,3.45 -"Reference","Peru",2050,"deaths",1115,2.83 -"Reference","Philippines",2020,"deaths",42167,37.89 -"Reference","Philippines",2030,"deaths",48062,37.52 -"Reference","Philippines",2040,"deaths",53344,37.34 -"Reference","Philippines",2050,"deaths",57431,37.21 -"Reference","Poland",2020,"deaths",2813,7.31 -"Reference","Poland",2030,"deaths",2789,7.3 -"Reference","Poland",2040,"deaths",2670,7.3 -"Reference","Poland",2050,"deaths",2559,7.29 -"Reference","Portugal",2020,"deaths",0,0 -"Reference","Portugal",2030,"deaths",0,0 -"Reference","Portugal",2040,"deaths",0,0 -"Reference","Portugal",2050,"deaths",0,0 -"Reference","Puerto Rico",2020,"deaths",2,0.07 -"Reference","Puerto Rico",2030,"deaths",1,0.05 -"Reference","Puerto Rico",2040,"deaths",1,0.04 -"Reference","Puerto Rico",2050,"deaths",1,0.04 -"Reference","Qatar",2020,"deaths",0,0 -"Reference","Qatar",2030,"deaths",0,0 -"Reference","Qatar",2040,"deaths",0,0 -"Reference","Qatar",2050,"deaths",0,0 -"Reference","Romania",2020,"deaths",3337,17.12 -"Reference","Romania",2030,"deaths",2399,13.1 -"Reference","Romania",2040,"deaths",1954,11.42 -"Reference","Romania",2050,"deaths",1677,10.52 -"Reference","Russian Federation",2020,"deaths",0,0 -"Reference","Russian Federation",2030,"deaths",0,0 -"Reference","Russian Federation",2040,"deaths",0,0 -"Reference","Russian Federation",2050,"deaths",0,0 -"Reference","Rwanda",2020,"deaths",14960,115.19 -"Reference","Rwanda",2030,"deaths",13576,83.81 -"Reference","Rwanda",2040,"deaths",14050,72.3 -"Reference","Rwanda",2050,"deaths",15080,67.16 -"Reference","Saint Lucia",2020,"deaths",43,23.94 -"Reference","Saint Lucia",2030,"deaths",24,12.96 -"Reference","Saint Lucia",2040,"deaths",15,8.42 -"Reference","Saint Lucia",2050,"deaths",10,5.39 -"Reference","Saudi Arabia",2020,"deaths",0,0 -"Reference","Saudi Arabia",2030,"deaths",0,0 -"Reference","Saudi Arabia",2040,"deaths",0,0 -"Reference","Saudi Arabia",2050,"deaths",0,0 -"Reference","Senegal",2020,"deaths",14534,89.63 -"Reference","Senegal",2030,"deaths",12666,61.38 -"Reference","Senegal",2040,"deaths",12030,47.45 -"Reference","Senegal",2050,"deaths",12283,40.89 -"Reference","Serbia",2020,"deaths",2858,38.71 -"Reference","Serbia",2030,"deaths",2116,30.04 -"Reference","Serbia",2040,"deaths",1785,26.19 -"Reference","Serbia",2050,"deaths",1612,24.38 -"Reference","Seychelles",2020,"deaths",1,1.2 -"Reference","Seychelles",2030,"deaths",1,0.82 -"Reference","Seychelles",2040,"deaths",1,0.66 -"Reference","Seychelles",2050,"deaths",1,0.57 -"Reference","Sierra Leone",2020,"deaths",8125,99.81 -"Reference","Sierra Leone",2030,"deaths",6503,65.05 -"Reference","Sierra Leone",2040,"deaths",4523,38.47 -"Reference","Sierra Leone",2050,"deaths",3203,24.13 -"Reference","Singapore",2020,"deaths",0,0 -"Reference","Singapore",2030,"deaths",0,0 -"Reference","Singapore",2040,"deaths",0,0 -"Reference","Singapore",2050,"deaths",0,0 -"Reference","Slovakia",2020,"deaths",0,0 -"Reference","Slovakia",2030,"deaths",0,0 -"Reference","Slovakia",2040,"deaths",0,0 -"Reference","Slovakia",2050,"deaths",0,0 -"Reference","Slovenia",2020,"deaths",120,5.65 -"Reference","Slovenia",2030,"deaths",88,4.18 -"Reference","Slovenia",2040,"deaths",73,3.53 -"Reference","Slovenia",2050,"deaths",64,3.14 -"Reference","South Africa",2020,"deaths",0,0 -"Reference","South Africa",2030,"deaths",0,0 -"Reference","South Africa",2040,"deaths",0,0 -"Reference","South Africa",2050,"deaths",0,0 -"Reference","South Sudan",2020,"deaths",34500,327.16 -"Reference","South Sudan",2030,"deaths",21136,180.46 -"Reference","South Sudan",2040,"deaths",13352,102.23 -"Reference","South Sudan",2050,"deaths",7797,56.58 -"Reference","Spain",2020,"deaths",459,0.97 -"Reference","Spain",2030,"deaths",311,0.65 -"Reference","Spain",2040,"deaths",228,0.47 -"Reference","Spain",2050,"deaths",163,0.33 -"Reference","Sri Lanka",2020,"deaths",10323,47.61 -"Reference","Sri Lanka",2030,"deaths",9276,41 -"Reference","Sri Lanka",2040,"deaths",7987,34.53 -"Reference","Sri Lanka",2050,"deaths",7214,31.12 -"Reference","Sweden",2020,"deaths",0,0 -"Reference","Sweden",2030,"deaths",0,0 -"Reference","Sweden",2040,"deaths",0,0 -"Reference","Sweden",2050,"deaths",0,0 -"Reference","Switzerland",2020,"deaths",0,0 -"Reference","Switzerland",2030,"deaths",0,0 -"Reference","Switzerland",2040,"deaths",0,0 -"Reference","Switzerland",2050,"deaths",0,0 -"Reference","Syrian Arab Republic",2020,"deaths",0,0 -"Reference","Syrian Arab Republic",2030,"deaths",0,0 -"Reference","Syrian Arab Republic",2040,"deaths",0,0 -"Reference","Syrian Arab Republic",2050,"deaths",0,0 -"Reference","Taiwan (Province of China)",2020,"deaths",0,0 -"Reference","Taiwan (Province of China)",2030,"deaths",0,0 -"Reference","Taiwan (Province of China)",2040,"deaths",0,0 -"Reference","Taiwan (Province of China)",2050,"deaths",0,0 -"Reference","Tajikistan",2020,"deaths",728,7.71 -"Reference","Tajikistan",2030,"deaths",0,0 -"Reference","Tajikistan",2040,"deaths",0,0 -"Reference","Tajikistan",2050,"deaths",0,0 -"Reference","United Republic of Tanzania",2020,"deaths",48988,80.61 -"Reference","United Republic of Tanzania",2030,"deaths",59399,73.92 -"Reference","United Republic of Tanzania",2040,"deaths",71646,70.34 -"Reference","United Republic of Tanzania",2050,"deaths",84859,68.32 -"Reference","Thailand",2020,"deaths",37151,52.04 -"Reference","Thailand",2030,"deaths",20018,27.73 -"Reference","Thailand",2040,"deaths",10164,14.32 -"Reference","Thailand",2050,"deaths",5184,7.63 -"Reference","Togo",2020,"deaths",2906,34.84 -"Reference","Togo",2030,"deaths",3598,34.63 -"Reference","Togo",2040,"deaths",4288,34.43 -"Reference","Togo",2050,"deaths",4927,34.29 -"Reference","Trinidad and Tobago",2020,"deaths",1,0.05 -"Reference","Trinidad and Tobago",2030,"deaths",0,0 -"Reference","Trinidad and Tobago",2040,"deaths",0,0 -"Reference","Trinidad and Tobago",2050,"deaths",0,0 -"Reference","Tunisia",2020,"deaths",11,0.09 -"Reference","Tunisia",2030,"deaths",0,0 -"Reference","Tunisia",2040,"deaths",0,0 -"Reference","Tunisia",2050,"deaths",0,0 -"Reference","Turkey",2020,"deaths",344,0.41 -"Reference","Turkey",2030,"deaths",351,0.39 -"Reference","Turkey",2040,"deaths",292,0.31 -"Reference","Turkey",2050,"deaths",251,0.26 -"Reference","Turkmenistan",2020,"deaths",3097,49.91 -"Reference","Turkmenistan",2030,"deaths",1796,25.15 -"Reference","Turkmenistan",2040,"deaths",977,12.41 -"Reference","Turkmenistan",2050,"deaths",556,6.53 -"Reference","Uganda",2020,"deaths",28466,65.16 -"Reference","Uganda",2030,"deaths",36157,61.94 -"Reference","Uganda",2040,"deaths",43125,57.97 -"Reference","Uganda",2050,"deaths",51654,56.92 -"Reference","Ukraine",2020,"deaths",2835,6.43 -"Reference","Ukraine",2030,"deaths",2371,6.42 -"Reference","Ukraine",2040,"deaths",2343,6.42 -"Reference","Ukraine",2050,"deaths",2245,6.42 -"Reference","United Arab Emirates",2020,"deaths",0,0 -"Reference","United Arab Emirates",2030,"deaths",0,0 -"Reference","United Arab Emirates",2040,"deaths",0,0 -"Reference","United Arab Emirates",2050,"deaths",0,0 -"Reference","United Kingdom",2020,"deaths",35799,53.47 -"Reference","United Kingdom",2030,"deaths",19151,27.87 -"Reference","United Kingdom",2040,"deaths",9377,13.29 -"Reference","United Kingdom",2050,"deaths",4878,6.7 -"Reference","United States of America",2020,"deaths",134,0.04 -"Reference","United States of America",2030,"deaths",70,0.02 -"Reference","United States of America",2040,"deaths",36,0.01 -"Reference","United States of America",2050,"deaths",38,0.01 -"Reference","Uruguay",2020,"deaths",0,0 -"Reference","Uruguay",2030,"deaths",0,0 -"Reference","Uruguay",2040,"deaths",0,0 -"Reference","Uruguay",2050,"deaths",0,0 -"Reference","Uzbekistan",2020,"deaths",6715,20.2 -"Reference","Uzbekistan",2030,"deaths",5833,15.14 -"Reference","Uzbekistan",2040,"deaths",4826,11.25 -"Reference","Uzbekistan",2050,"deaths",4488,9.58 -"Reference","Venezuela (Bolivarian Republic of)",2020,"deaths",0,0 -"Reference","Venezuela (Bolivarian Republic of)",2030,"deaths",0,0 -"Reference","Venezuela (Bolivarian Republic of)",2040,"deaths",0,0 -"Reference","Venezuela (Bolivarian Republic of)",2050,"deaths",0,0 -"Reference","Viet Nam",2020,"deaths",48650,50.57 -"Reference","Viet Nam",2030,"deaths",41883,40.78 -"Reference","Viet Nam",2040,"deaths",39630,37.31 -"Reference","Viet Nam",2050,"deaths",38777,36.05 -"Reference","Yemen",2020,"deaths",70320,220.25 -"Reference","Yemen",2030,"deaths",63112,159.3 -"Reference","Yemen",2040,"deaths",39574,84.12 -"Reference","Yemen",2050,"deaths",24676,45.97 -"Reference","Zambia",2020,"deaths",12320,66.04 -"Reference","Zambia",2030,"deaths",11244,46.71 -"Reference","Zambia",2040,"deaths",9634,32.47 -"Reference","Zambia",2050,"deaths",8631,24.7 -"Reference","Zimbabwe",2020,"deaths",10373,66.9 -"Reference","Zimbabwe",2030,"deaths",7735,42.03 -"Reference","Zimbabwe",2040,"deaths",5612,26.7 -"Reference","Zimbabwe",2050,"deaths",4166,18.02 +"scenario","country","year","pred_var","pred_value","pred_value_normalized","reliability_ratio","reliability" +"Reference","Afghanistan","2020","deaths",8658,22.59,0.24,"high" +"Reference","Afghanistan","2030","deaths",9592,19.03,0.24,"high" +"Reference","Afghanistan","2040","deaths",10820,17.12,0.24,"high" +"Reference","Afghanistan","2050","deaths",12213,16.15,0.24,"high" +"Reference","Albania","2020","deaths",800,27.85,43.361,"low" +"Reference","Albania","2030","deaths",716,26.37,43.361,"low" +"Reference","Albania","2040","deaths",643,25.7,43.361,"low" +"Reference","Albania","2050","deaths",575,25.37,43.361,"low" +"Reference","Algeria","2020","deaths",13,0.03,0.962,"medium" +"Reference","Algeria","2030","deaths",10,0.02,0.962,"medium" +"Reference","Algeria","2040","deaths",6,0.01,0.962,"medium" +"Reference","Algeria","2050","deaths",6,0.01,0.962,"medium" +"Reference","Angola","2020","deaths",7748,23.55,1.214,"medium" +"Reference","Angola","2030","deaths",9580,21.72,1.214,"medium" +"Reference","Angola","2040","deaths",11583,20.52,1.214,"medium" +"Reference","Angola","2050","deaths",13689,19.84,1.214,"medium" +"Reference","Argentina","2020","deaths",1055,2.35,0.114,"high" +"Reference","Argentina","2030","deaths",548,1.16,0.114,"high" +"Reference","Argentina","2040","deaths",266,0.54,0.114,"high" +"Reference","Argentina","2050","deaths",96,0.19,0.114,"high" +"Reference","Armenia","2020","deaths",4,0.13,0.444,"high" +"Reference","Armenia","2030","deaths",0,0.01,0.444,"high" +"Reference","Armenia","2040","deaths",0,0.01,0.444,"high" +"Reference","Armenia","2050","deaths",0,0,0.444,"high" +"Reference","Australia","2020","deaths",378,1.48,0.287,"high" +"Reference","Australia","2030","deaths",218,0.79,0.287,"high" +"Reference","Australia","2040","deaths",123,0.42,0.287,"high" +"Reference","Australia","2050","deaths",68,0.22,0.287,"high" +"Reference","Austria","2020","deaths",135,1.52,0.048,"high" +"Reference","Austria","2030","deaths",76,0.83,0.048,"high" +"Reference","Austria","2040","deaths",43,0.47,0.048,"high" +"Reference","Austria","2050","deaths",24,0.26,0.048,"high" +"Reference","Azerbaijan","2020","deaths",180,1.75,0.114,"high" +"Reference","Azerbaijan","2030","deaths",51,0.48,0.114,"high" +"Reference","Azerbaijan","2040","deaths",1,0.01,0.114,"high" +"Reference","Azerbaijan","2050","deaths",0,0,0.114,"high" +"Reference","Bahrain","2020","deaths",0,0.01,0.992,"medium" +"Reference","Bahrain","2030","deaths",0,0.01,0.992,"medium" +"Reference","Bahrain","2040","deaths",0,0,0.992,"medium" +"Reference","Bahrain","2050","deaths",0,0,0.992,"medium" +"Reference","Bangladesh","2020","deaths",104966,63.07,848.64,"low" +"Reference","Bangladesh","2030","deaths",111234,60.66,848.64,"low" +"Reference","Bangladesh","2040","deaths",115779,59.63,848.64,"low" +"Reference","Bangladesh","2050","deaths",117916,59.17,848.64,"low" +"Reference","Barbados","2020","deaths",0,0.03,0.999,"medium" +"Reference","Barbados","2030","deaths",0,0.02,0.999,"medium" +"Reference","Barbados","2040","deaths",0,0.01,0.999,"medium" +"Reference","Barbados","2050","deaths",0,0,0.999,"medium" +"Reference","Belarus","2020","deaths",2,0.02,0.942,"medium" +"Reference","Belarus","2030","deaths",1,0.01,0.942,"medium" +"Reference","Belarus","2040","deaths",1,0.01,0.942,"medium" +"Reference","Belarus","2050","deaths",0,0,0.942,"medium" +"Reference","Belgium","2020","deaths",2,0.02,0.998,"medium" +"Reference","Belgium","2030","deaths",1,0.01,0.998,"medium" +"Reference","Belgium","2040","deaths",0,0,0.998,"medium" +"Reference","Belgium","2050","deaths",0,0,0.998,"medium" +"Reference","Benin","2020","deaths",12705,101.91,1300.328,"low" +"Reference","Benin","2030","deaths",15994,98.8,1300.328,"low" +"Reference","Benin","2040","deaths",19595,97.24,1300.328,"low" +"Reference","Benin","2050","deaths",23273,96.47,1300.328,"low" +"Reference","Bolivia (Plurinational State of)","2020","deaths",2237,18.86,2.152,"low" +"Reference","Bolivia (Plurinational State of)","2030","deaths",2325,17.3,2.152,"low" +"Reference","Bolivia (Plurinational State of)","2040","deaths",2439,16.38,2.152,"low" +"Reference","Bolivia (Plurinational State of)","2050","deaths",2531,15.86,2.152,"low" +"Reference","Bosnia and Herzegovina","2020","deaths",467,13.99,0.467,"high" +"Reference","Bosnia and Herzegovina","2030","deaths",385,12.57,0.467,"high" +"Reference","Bosnia and Herzegovina","2040","deaths",335,11.92,0.467,"high" +"Reference","Bosnia and Herzegovina","2050","deaths",296,11.6,0.467,"high" +"Reference","Botswana","2020","deaths",1,0.03,0.548,"medium" +"Reference","Botswana","2030","deaths",1,0.02,0.548,"medium" +"Reference","Botswana","2040","deaths",0,0.01,0.548,"medium" +"Reference","Botswana","2050","deaths",0,0,0.548,"medium" +"Reference","Brazil","2020","deaths",64,0.03,0.548,"medium" +"Reference","Brazil","2030","deaths",22,0.01,0.548,"medium" +"Reference","Brazil","2040","deaths",23,0.01,0.548,"medium" +"Reference","Brazil","2050","deaths",0,0,0.548,"medium" +"Reference","Bulgaria","2020","deaths",1,0.02,0.559,"medium" +"Reference","Bulgaria","2030","deaths",1,0.01,0.559,"medium" +"Reference","Bulgaria","2040","deaths",1,0.01,0.559,"medium" +"Reference","Bulgaria","2050","deaths",0,0,0.559,"medium" +"Reference","Burkina Faso","2020","deaths",28886,136.05,72.924,"low" +"Reference","Burkina Faso","2030","deaths",36235,132.47,72.924,"low" +"Reference","Burkina Faso","2040","deaths",44539,130.51,72.924,"low" +"Reference","Burkina Faso","2050","deaths",53108,129.56,72.924,"low" +"Reference","Burundi","2020","deaths",12923,107.21,380.096,"low" +"Reference","Burundi","2030","deaths",15524,101.16,380.096,"low" +"Reference","Burundi","2040","deaths",18462,97.76,380.096,"low" +"Reference","Burundi","2050","deaths",21429,96.04,380.096,"low" +"Reference","Cambodia","2020","deaths",15591,95.67,2709.468,"low" +"Reference","Cambodia","2030","deaths",16935,93.06,2709.468,"low" +"Reference","Cambodia","2040","deaths",18008,91.8,2709.468,"low" +"Reference","Cambodia","2050","deaths",18602,91.22,2709.468,"low" +"Reference","Cameroon","2020","deaths",15510,59.34,24.75,"low" +"Reference","Cameroon","2030","deaths",18828,56.71,24.75,"low" +"Reference","Cameroon","2040","deaths",22255,55.19,24.75,"low" +"Reference","Cameroon","2050","deaths",25630,54.41,24.75,"low" +"Reference","Canada","2020","deaths",8,0.02,1,"medium" +"Reference","Canada","2030","deaths",4,0.01,1,"medium" +"Reference","Canada","2040","deaths",0,0,1,"medium" +"Reference","Canada","2050","deaths",0,0,1,"medium" +"Reference","Cabo Verde","2020","deaths",134,23.05,8.899,"low" +"Reference","Cabo Verde","2030","deaths",131,20.62,8.899,"low" +"Reference","Cabo Verde","2040","deaths",132,19.58,8.899,"low" +"Reference","Cabo Verde","2050","deaths",132,19.07,8.899,"low" +"Reference","Central African Republic","2020","deaths",8016,152.04,8.107,"low" +"Reference","Central African Republic","2030","deaths",10371,147.31,8.107,"low" +"Reference","Central African Republic","2040","deaths",13030,143.88,8.107,"low" +"Reference","Central African Republic","2050","deaths",15751,142.11,8.107,"low" +"Reference","Chad","2020","deaths",22733,138.79,63.129,"low" +"Reference","Chad","2030","deaths",30186,135.27,63.129,"low" +"Reference","Chad","2040","deaths",38810,132.75,63.129,"low" +"Reference","Chad","2050","deaths",48131,131.47,63.129,"low" +"Reference","Chile","2020","deaths",795,4.14,40.956,"low" +"Reference","Chile","2030","deaths",618,3.1,40.956,"low" +"Reference","Chile","2040","deaths",518,2.54,40.956,"low" +"Reference","Chile","2050","deaths",456,2.23,40.956,"low" +"Reference","China","2020","deaths",427,0.03,0.707,"medium" +"Reference","China","2030","deaths",142,0.01,0.707,"medium" +"Reference","China","2040","deaths",138,0.01,0.707,"medium" +"Reference","China","2050","deaths",0,0,0.707,"medium" +"Reference","Colombia","2020","deaths",15,0.03,0.894,"medium" +"Reference","Colombia","2030","deaths",5,0.01,0.894,"medium" +"Reference","Colombia","2040","deaths",6,0.01,0.894,"medium" +"Reference","Colombia","2050","deaths",0,0,0.894,"medium" +"Reference","Comoros","2020","deaths",383,48.03,1.294,"medium" +"Reference","Comoros","2030","deaths",435,45.15,1.294,"medium" +"Reference","Comoros","2040","deaths",486,43.37,1.294,"medium" +"Reference","Comoros","2050","deaths",536,42.45,1.294,"medium" +"Reference","Congo","2020","deaths",2037,36.14,5.255,"low" +"Reference","Congo","2030","deaths",2380,33.65,5.255,"low" +"Reference","Congo","2040","deaths",2754,32.17,5.255,"low" +"Reference","Congo","2050","deaths",3120,31.35,5.255,"low" +"Reference","Democratic Republic of the Congo","2020","deaths",60270,65.99,2.502,"low" +"Reference","Democratic Republic of the Congo","2030","deaths",76674,61.07,2.502,"low" +"Reference","Democratic Republic of the Congo","2040","deaths",96602,58.42,2.502,"low" +"Reference","Democratic Republic of the Congo","2050","deaths",118839,57.02,2.502,"low" +"Reference","Costa Rica","2020","deaths",2,0.03,0.953,"medium" +"Reference","Costa Rica","2030","deaths",1,0.01,0.953,"medium" +"Reference","Costa Rica","2040","deaths",1,0.01,0.953,"medium" +"Reference","Costa Rica","2050","deaths",0,0,0.953,"medium" +"Reference","Côte d'Ivoire","2020","deaths",4136,15.62,0.165,"high" +"Reference","Côte d'Ivoire","2030","deaths",4343,13.13,0.165,"high" +"Reference","Côte d'Ivoire","2040","deaths",4722,11.94,0.165,"high" +"Reference","Côte d'Ivoire","2050","deaths",5154,11.37,0.165,"high" +"Reference","Croatia","2020","deaths",1,0.02,0.938,"medium" +"Reference","Croatia","2030","deaths",0,0.01,0.938,"medium" +"Reference","Croatia","2040","deaths",0,0.01,0.938,"medium" +"Reference","Croatia","2050","deaths",0,0,0.938,"medium" +"Reference","Cuba","2020","deaths",248,2.19,0.5,"high" +"Reference","Cuba","2030","deaths",2,0.02,0.5,"high" +"Reference","Cuba","2040","deaths",1,0.01,0.5,"high" +"Reference","Cuba","2050","deaths",0,0,0.5,"high" +"Reference","Cyprus","2020","deaths",6,0.52,0.953,"medium" +"Reference","Cyprus","2030","deaths",0,0.01,0.953,"medium" +"Reference","Cyprus","2040","deaths",0,0.01,0.953,"medium" +"Reference","Cyprus","2050","deaths",0,0,0.953,"medium" +"Reference","Czechia","2020","deaths",267,2.53,6.404,"low" +"Reference","Czechia","2030","deaths",184,1.72,6.404,"low" +"Reference","Czechia","2040","deaths",131,1.26,6.404,"low" +"Reference","Czechia","2050","deaths",103,1.01,6.404,"low" +"Reference","Denmark","2020","deaths",87,1.49,0.122,"high" +"Reference","Denmark","2030","deaths",49,0.82,0.122,"high" +"Reference","Denmark","2040","deaths",29,0.47,0.122,"high" +"Reference","Denmark","2050","deaths",16,0.26,0.122,"high" +"Reference","Djibouti","2020","deaths",1,0.05,0.701,"medium" +"Reference","Djibouti","2030","deaths",0,0.03,0.701,"medium" +"Reference","Djibouti","2040","deaths",0,0.01,0.701,"medium" +"Reference","Djibouti","2050","deaths",0,0.01,0.701,"medium" +"Reference","Dominican Republic","2020","deaths",3,0.03,0.376,"high" +"Reference","Dominican Republic","2030","deaths",1,0.01,0.376,"high" +"Reference","Dominican Republic","2040","deaths",1,0.01,0.376,"high" +"Reference","Dominican Republic","2050","deaths",0,0,0.376,"high" +"Reference","Ecuador","2020","deaths",5,0.03,0.568,"medium" +"Reference","Ecuador","2030","deaths",4,0.02,0.568,"medium" +"Reference","Ecuador","2040","deaths",2,0.01,0.568,"medium" +"Reference","Ecuador","2050","deaths",0,0,0.568,"medium" +"Reference","Egypt","2020","deaths",32,0.03,0.966,"medium" +"Reference","Egypt","2030","deaths",13,0.01,0.966,"medium" +"Reference","Egypt","2040","deaths",15,0.01,0.966,"medium" +"Reference","Egypt","2050","deaths",0,0,0.966,"medium" +"Reference","El Salvador","2020","deaths",984,15.66,137.338,"low" +"Reference","El Salvador","2030","deaths",897,13.9,137.338,"low" +"Reference","El Salvador","2040","deaths",841,12.97,137.338,"low" +"Reference","El Salvador","2050","deaths",795,12.47,137.338,"low" +"Reference","Equatorial Guinea","2020","deaths",0,0.03,0.674,"medium" +"Reference","Equatorial Guinea","2030","deaths",0,0.02,0.674,"medium" +"Reference","Equatorial Guinea","2040","deaths",0,0.01,0.674,"medium" +"Reference","Equatorial Guinea","2050","deaths",0,0.01,0.674,"medium" +"Reference","Estonia","2020","deaths",0,0.02,0.712,"medium" +"Reference","Estonia","2030","deaths",0,0.01,0.712,"medium" +"Reference","Estonia","2040","deaths",0,0.01,0.712,"medium" +"Reference","Estonia","2050","deaths",0,0,0.712,"medium" +"Reference","Ethiopia","2020","deaths",65949,57.03,2.53,"low" +"Reference","Ethiopia","2030","deaths",78265,53.33,2.53,"low" +"Reference","Ethiopia","2040","deaths",91993,51.51,2.53,"low" +"Reference","Ethiopia","2050","deaths",105720,50.64,2.53,"low" +"Reference","Finland","2020","deaths",1,0.02,0.991,"medium" +"Reference","Finland","2030","deaths",1,0.01,0.991,"medium" +"Reference","Finland","2040","deaths",1,0.01,0.991,"medium" +"Reference","Finland","2050","deaths",0,0,0.991,"medium" +"Reference","France","2020","deaths",13,0.02,0.998,"medium" +"Reference","France","2030","deaths",7,0.01,0.998,"medium" +"Reference","France","2040","deaths",7,0.01,0.998,"medium" +"Reference","France","2050","deaths",0,0,0.998,"medium" +"Reference","Gabon","2020","deaths",1,0.03,0.911,"medium" +"Reference","Gabon","2030","deaths",1,0.02,0.911,"medium" +"Reference","Gabon","2040","deaths",0,0.01,0.911,"medium" +"Reference","Gabon","2050","deaths",0,0.01,0.911,"medium" +"Reference","Gambia","2020","deaths",1197,47.11,1.277,"medium" +"Reference","Gambia","2030","deaths",1375,43.48,1.277,"medium" +"Reference","Gambia","2040","deaths",1584,41.59,1.277,"medium" +"Reference","Gambia","2050","deaths",1786,40.66,1.277,"medium" +"Reference","Georgia","2020","deaths",1828,48.5,691.213,"low" +"Reference","Georgia","2030","deaths",1629,46.97,691.213,"low" +"Reference","Georgia","2040","deaths",1472,46.33,691.213,"low" +"Reference","Georgia","2050","deaths",1344,46.06,691.213,"low" +"Reference","Germany","2020","deaths",733,0.88,0.941,"medium" +"Reference","Germany","2030","deaths",183,0.22,0.941,"medium" +"Reference","Germany","2040","deaths",0,0,0.941,"medium" +"Reference","Germany","2050","deaths",0,0,0.941,"medium" +"Reference","Ghana","2020","deaths",12141,38.12,3.53,"low" +"Reference","Ghana","2030","deaths",13933,35.92,3.53,"low" +"Reference","Ghana","2040","deaths",15907,34.69,3.53,"low" +"Reference","Ghana","2050","deaths",17968,34.06,3.53,"low" +"Reference","Greece","2020","deaths",213,2.02,0.065,"high" +"Reference","Greece","2030","deaths",104,1.02,0.065,"high" +"Reference","Greece","2040","deaths",55,0.55,0.065,"high" +"Reference","Greece","2050","deaths",28,0.29,0.065,"high" +"Reference","Guatemala","2020","deaths",1599,9.28,0.207,"high" +"Reference","Guatemala","2030","deaths",1494,7.47,0.207,"high" +"Reference","Guatemala","2040","deaths",1459,6.5,0.207,"high" +"Reference","Guatemala","2050","deaths",1459,6.01,0.207,"high" +"Reference","Guinea","2020","deaths",11466,87.92,1.725,"medium" +"Reference","Guinea","2030","deaths",13726,84.58,1.725,"medium" +"Reference","Guinea","2040","deaths",16019,82.79,1.725,"medium" +"Reference","Guinea","2050","deaths",18179,81.92,1.725,"medium" +"Reference","Guinea-Bissau","2020","deaths",2109,105.81,3044.161,"low" +"Reference","Guinea-Bissau","2030","deaths",2489,102.14,3044.161,"low" +"Reference","Guinea-Bissau","2040","deaths",2861,100.22,3044.161,"low" +"Reference","Guinea-Bissau","2050","deaths",3185,99.22,3044.161,"low" +"Reference","Haiti","2020","deaths",5681,50.57,0.672,"medium" +"Reference","Haiti","2030","deaths",6160,48.22,0.672,"medium" +"Reference","Haiti","2040","deaths",6517,46.39,0.672,"medium" +"Reference","Haiti","2050","deaths",6792,45.47,0.672,"medium" +"Reference","Honduras","2020","deaths",5,0.05,0.468,"high" +"Reference","Honduras","2030","deaths",2,0.02,0.468,"high" +"Reference","Honduras","2040","deaths",1,0.01,0.468,"high" +"Reference","Honduras","2050","deaths",1,0.01,0.468,"high" +"Reference","Hungary","2020","deaths",2,0.02,0.311,"high" +"Reference","Hungary","2030","deaths",1,0.01,0.311,"high" +"Reference","Hungary","2040","deaths",1,0.01,0.311,"high" +"Reference","Hungary","2050","deaths",0,0,0.311,"high" +"Reference","Iceland","2020","deaths",1,0.24,0.982,"medium" +"Reference","Iceland","2030","deaths",0,0.01,0.982,"medium" +"Reference","Iceland","2040","deaths",0,0,0.982,"medium" +"Reference","Iceland","2050","deaths",0,0,0.982,"medium" +"Reference","India","2020","deaths",741964,53.38,1623.44,"low" +"Reference","India","2030","deaths",766180,51.02,1623.44,"low" +"Reference","India","2040","deaths",792160,50.02,1623.44,"low" +"Reference","India","2050","deaths",802563,49.53,1623.44,"low" +"Reference","Indonesia","2020","deaths",67111,24.78,2.139,"low" +"Reference","Indonesia","2030","deaths",67357,23.18,2.139,"low" +"Reference","Indonesia","2040","deaths",68268,22.41,2.139,"low" +"Reference","Indonesia","2050","deaths",68425,22.03,2.139,"low" +"Reference","Iran (Islamic Republic of)","2020","deaths",26,0.03,0.991,"medium" +"Reference","Iran (Islamic Republic of)","2030","deaths",9,0.01,0.991,"medium" +"Reference","Iran (Islamic Republic of)","2040","deaths",10,0.01,0.991,"medium" +"Reference","Iran (Islamic Republic of)","2050","deaths",0,0,0.991,"medium" +"Reference","Iraq","2020","deaths",1181,2.81,0.49,"high" +"Reference","Iraq","2030","deaths",797,1.53,0.49,"high" +"Reference","Iraq","2040","deaths",374,0.6,0.49,"high" +"Reference","Iraq","2050","deaths",115,0.16,0.49,"high" +"Reference","Ireland","2020","deaths",0,0.01,0.999,"medium" +"Reference","Ireland","2030","deaths",1,0.01,0.999,"medium" +"Reference","Ireland","2040","deaths",0,0,0.999,"medium" +"Reference","Ireland","2050","deaths",0,0,0.999,"medium" +"Reference","Israel","2020","deaths",132,1.52,0.236,"high" +"Reference","Israel","2030","deaths",79,0.79,0.236,"high" +"Reference","Israel","2040","deaths",48,0.42,0.236,"high" +"Reference","Israel","2050","deaths",30,0.23,0.236,"high" +"Reference","Italy","2020","deaths",1139,1.91,3.029,"low" +"Reference","Italy","2030","deaths",646,1.11,3.029,"low" +"Reference","Italy","2040","deaths",393,0.69,3.029,"low" +"Reference","Italy","2050","deaths",257,0.46,3.029,"low" +"Reference","Jamaica","2020","deaths",307,10.92,29.567,"low" +"Reference","Jamaica","2030","deaths",265,9.29,29.567,"low" +"Reference","Jamaica","2040","deaths",239,8.41,29.567,"low" +"Reference","Jamaica","2050","deaths",220,7.92,29.567,"low" +"Reference","Japan","2020","deaths",1833,1.46,0.045,"high" +"Reference","Japan","2030","deaths",979,0.82,0.045,"high" +"Reference","Japan","2040","deaths",530,0.47,0.045,"high" +"Reference","Japan","2050","deaths",286,0.27,0.045,"high" +"Reference","Jordan","2020","deaths",3,0.03,0.998,"medium" +"Reference","Jordan","2030","deaths",3,0.02,0.998,"medium" +"Reference","Jordan","2040","deaths",2,0.01,0.998,"medium" +"Reference","Jordan","2050","deaths",0,0,0.998,"medium" +"Reference","Kazakhstan","2020","deaths",2154,11.42,209.912,"low" +"Reference","Kazakhstan","2030","deaths",2157,10.45,209.912,"low" +"Reference","Kazakhstan","2040","deaths",2192,9.92,209.912,"low" +"Reference","Kazakhstan","2050","deaths",2253,9.65,209.912,"low" +"Reference","Kenya","2020","deaths",20991,40.79,3.371,"low" +"Reference","Kenya","2030","deaths",23841,38.17,3.371,"low" +"Reference","Kenya","2040","deaths",27069,36.87,3.371,"low" +"Reference","Kenya","2050","deaths",30072,36.21,3.371,"low" +"Reference","Democratic People's Republic of Korea","2020","deaths",35439,137.29,2400.126,"low" +"Reference","Democratic People's Republic of Korea","2030","deaths",35354,133.25,2400.126,"low" +"Reference","Democratic People's Republic of Korea","2040","deaths",34224,131.24,2400.126,"low" +"Reference","Democratic People's Republic of Korea","2050","deaths",32858,130.33,2400.126,"low" +"Reference","Republic of Korea","2020","deaths",10,0.02,0.999,"medium" +"Reference","Republic of Korea","2030","deaths",5,0.01,0.999,"medium" +"Reference","Republic of Korea","2040","deaths",0,0,0.999,"medium" +"Reference","Republic of Korea","2050","deaths",0,0,0.999,"medium" +"Reference","Kuwait","2020","deaths",1,0.02,1,"medium" +"Reference","Kuwait","2030","deaths",0,0.01,1,"medium" +"Reference","Kuwait","2040","deaths",0,0,1,"medium" +"Reference","Kuwait","2050","deaths",0,0,1,"medium" +"Reference","Kyrgyzstan","2020","deaths",286,4.49,0.005,"high" +"Reference","Kyrgyzstan","2030","deaths",172,2.33,0.005,"high" +"Reference","Kyrgyzstan","2040","deaths",93,1.13,0.005,"high" +"Reference","Kyrgyzstan","2050","deaths",46,0.51,0.005,"high" +"Reference","Lao People's Democratic Republic","2020","deaths",6472,89.08,18.185,"low" +"Reference","Lao People's Democratic Republic","2030","deaths",7168,87.24,18.185,"low" +"Reference","Lao People's Democratic Republic","2040","deaths",7695,86.25,18.185,"low" +"Reference","Lao People's Democratic Republic","2050","deaths",7999,85.75,18.185,"low" +"Reference","Latvia","2020","deaths",147,7.7,13.153,"low" +"Reference","Latvia","2030","deaths",115,6.65,13.153,"low" +"Reference","Latvia","2040","deaths",96,6.14,13.153,"low" +"Reference","Latvia","2050","deaths",85,5.88,13.153,"low" +"Reference","Lebanon","2020","deaths",123,2.16,0.708,"medium" +"Reference","Lebanon","2030","deaths",70,1.13,0.708,"medium" +"Reference","Lebanon","2040","deaths",28,0.41,0.708,"medium" +"Reference","Lebanon","2050","deaths",0,0,0.708,"medium" +"Reference","Lesotho","2020","deaths",2108,94.13,1684.814,"low" +"Reference","Lesotho","2030","deaths",2225,90.98,1684.814,"low" +"Reference","Lesotho","2040","deaths",2313,89.03,1684.814,"low" +"Reference","Lesotho","2050","deaths",2363,88.02,1684.814,"low" +"Reference","Liberia","2020","deaths",1096,21.78,0.26,"high" +"Reference","Liberia","2030","deaths",1083,17.17,0.26,"high" +"Reference","Liberia","2040","deaths",1132,14.86,0.26,"high" +"Reference","Liberia","2050","deaths",1208,13.68,0.26,"high" +"Reference","Libya","2020","deaths",2,0.03,0.998,"medium" +"Reference","Libya","2030","deaths",1,0.01,0.998,"medium" +"Reference","Libya","2040","deaths",1,0.01,0.998,"medium" +"Reference","Libya","2050","deaths",0,0,0.998,"medium" +"Reference","Lithuania","2020","deaths",1,0.02,0.957,"medium" +"Reference","Lithuania","2030","deaths",0,0.01,0.957,"medium" +"Reference","Lithuania","2040","deaths",0,0.01,0.957,"medium" +"Reference","Lithuania","2050","deaths",0,0,0.957,"medium" +"Reference","Luxembourg","2020","deaths",0,0.01,0.999,"medium" +"Reference","Luxembourg","2030","deaths",0,0.01,0.999,"medium" +"Reference","Luxembourg","2040","deaths",0,0,0.999,"medium" +"Reference","Luxembourg","2050","deaths",0,0,0.999,"medium" +"Reference","North Macedonia","2020","deaths",502,23.75,3.961,"low" +"Reference","North Macedonia","2030","deaths",476,22.62,3.961,"low" +"Reference","North Macedonia","2040","deaths",460,22.05,3.961,"low" +"Reference","North Macedonia","2050","deaths",445,21.76,3.961,"low" +"Reference","Madagascar","2020","deaths",19237,69,2.19,"low" +"Reference","Madagascar","2030","deaths",22651,64.76,2.19,"low" +"Reference","Madagascar","2040","deaths",26276,62.39,2.19,"low" +"Reference","Madagascar","2050","deaths",29923,61.18,2.19,"low" +"Reference","Malawi","2020","deaths",15624,81.71,1147.793,"low" +"Reference","Malawi","2030","deaths",19116,77.91,1147.793,"low" +"Reference","Malawi","2040","deaths",23006,75.47,1147.793,"low" +"Reference","Malawi","2050","deaths",26909,74.2,1147.793,"low" +"Reference","Malaysia","2020","deaths",7,0.02,0.977,"medium" +"Reference","Malaysia","2030","deaths",4,0.01,0.977,"medium" +"Reference","Malaysia","2040","deaths",4,0.01,0.977,"medium" +"Reference","Malaysia","2050","deaths",0,0,0.977,"medium" +"Reference","Mali","2020","deaths",22932,109.79,4.325,"low" +"Reference","Mali","2030","deaths",30009,106.47,4.325,"low" +"Reference","Mali","2040","deaths",37982,104.48,4.325,"low" +"Reference","Mali","2050","deaths",46027,103.48,4.325,"low" +"Reference","Malta","2020","deaths",0,0.02,0.992,"medium" +"Reference","Malta","2030","deaths",0,0.01,0.992,"medium" +"Reference","Malta","2040","deaths",0,0,0.992,"medium" +"Reference","Malta","2050","deaths",0,0,0.992,"medium" +"Reference","Mauritania","2020","deaths",2,0.05,0.363,"high" +"Reference","Mauritania","2030","deaths",2,0.03,0.363,"high" +"Reference","Mauritania","2040","deaths",1,0.01,0.363,"high" +"Reference","Mauritania","2050","deaths",1,0.01,0.363,"high" +"Reference","Mauritius","2020","deaths",53,4.06,31.38,"low" +"Reference","Mauritius","2030","deaths",35,2.66,31.38,"low" +"Reference","Mauritius","2040","deaths",26,2.01,31.38,"low" +"Reference","Mauritius","2050","deaths",21,1.67,31.38,"low" +"Reference","Mexico","2020","deaths",11405,9.08,3.86,"low" +"Reference","Mexico","2030","deaths",10795,7.96,3.86,"low" +"Reference","Mexico","2040","deaths",10606,7.36,3.86,"low" +"Reference","Mexico","2050","deaths",10519,7.04,3.86,"low" +"Reference","Republic of Moldova","2020","deaths",359,11.61,67.276,"low" +"Reference","Republic of Moldova","2030","deaths",289,9.87,67.276,"low" +"Reference","Republic of Moldova","2040","deaths",243,9.04,67.276,"low" +"Reference","Republic of Moldova","2050","deaths",217,8.69,67.276,"low" +"Reference","Mongolia","2020","deaths",1,0.03,0.343,"high" +"Reference","Mongolia","2030","deaths",1,0.02,0.343,"high" +"Reference","Mongolia","2040","deaths",0,0.01,0.343,"high" +"Reference","Mongolia","2050","deaths",0,0,0.343,"high" +"Reference","Montenegro","2020","deaths",0,0.03,0.585,"medium" +"Reference","Montenegro","2030","deaths",0,0.01,0.585,"medium" +"Reference","Montenegro","2040","deaths",0,0.01,0.585,"medium" +"Reference","Montenegro","2050","deaths",0,0,0.585,"medium" +"Reference","Morocco","2020","deaths",15,0.04,0.913,"medium" +"Reference","Morocco","2030","deaths",8,0.02,0.913,"medium" +"Reference","Morocco","2040","deaths",4,0.01,0.913,"medium" +"Reference","Morocco","2050","deaths",0,0,0.913,"medium" +"Reference","Mozambique","2020","deaths",23342,75.98,2.513,"low" +"Reference","Mozambique","2030","deaths",28078,70.94,2.513,"low" +"Reference","Mozambique","2040","deaths",33391,68.38,2.513,"low" +"Reference","Mozambique","2050","deaths",38857,67.33,2.513,"low" +"Reference","Myanmar","2020","deaths",47756,89.72,10.182,"low" +"Reference","Myanmar","2030","deaths",49454,87.85,10.182,"low" +"Reference","Myanmar","2040","deaths",49705,86.51,10.182,"low" +"Reference","Myanmar","2050","deaths",49034,85.81,10.182,"low" +"Reference","Namibia","2020","deaths",1,0.04,0.584,"medium" +"Reference","Namibia","2030","deaths",1,0.02,0.584,"medium" +"Reference","Namibia","2040","deaths",0,0.01,0.584,"medium" +"Reference","Namibia","2050","deaths",0,0.01,0.584,"medium" +"Reference","Nepal","2020","deaths",23652,81.56,37.493,"low" +"Reference","Nepal","2030","deaths",25147,78.68,37.493,"low" +"Reference","Nepal","2040","deaths",26114,77.22,37.493,"low" +"Reference","Nepal","2050","deaths",26702,76.53,37.493,"low" +"Reference","Netherlands","2020","deaths",101,0.58,0.96,"medium" +"Reference","Netherlands","2030","deaths",2,0.01,0.96,"medium" +"Reference","Netherlands","2040","deaths",0,0,0.96,"medium" +"Reference","Netherlands","2050","deaths",0,0,0.96,"medium" +"Reference","New Zealand","2020","deaths",1,0.02,0.984,"medium" +"Reference","New Zealand","2030","deaths",1,0.01,0.984,"medium" +"Reference","New Zealand","2040","deaths",0,0,0.984,"medium" +"Reference","New Zealand","2050","deaths",0,0,0.984,"medium" +"Reference","Nicaragua","2020","deaths",2181,32.5,940.822,"low" +"Reference","Nicaragua","2030","deaths",2288,30.25,940.822,"low" +"Reference","Nicaragua","2040","deaths",2383,29.05,940.822,"low" +"Reference","Nicaragua","2050","deaths",2444,28.43,940.822,"low" +"Reference","Niger","2020","deaths",21097,88.34,1.602,"medium" +"Reference","Niger","2030","deaths",28761,83.72,1.602,"medium" +"Reference","Niger","2040","deaths",38884,81.24,1.602,"medium" +"Reference","Niger","2050","deaths",51076,80.03,1.602,"medium" +"Reference","Nigeria","2020","deaths",126740,61.59,3.024,"low" +"Reference","Nigeria","2030","deaths",155167,59.42,3.024,"low" +"Reference","Nigeria","2040","deaths",190008,58.08,3.024,"low" +"Reference","Nigeria","2050","deaths",229822,57.37,3.024,"low" +"Reference","Norway","2020","deaths",1,0.01,1,"medium" +"Reference","Norway","2030","deaths",1,0.01,1,"medium" +"Reference","Norway","2040","deaths",0,0,1,"medium" +"Reference","Norway","2050","deaths",0,0,1,"medium" +"Reference","Oman","2020","deaths",1,0.02,0.985,"medium" +"Reference","Oman","2030","deaths",1,0.01,0.985,"medium" +"Reference","Oman","2040","deaths",1,0.01,0.985,"medium" +"Reference","Oman","2050","deaths",0,0,0.985,"medium" +"Reference","Pakistan","2020","deaths",126806,56.33,6.586,"low" +"Reference","Pakistan","2030","deaths",147782,54.01,6.586,"low" +"Reference","Pakistan","2040","deaths",170509,52.69,6.586,"low" +"Reference","Pakistan","2050","deaths",192689,52.01,6.586,"low" +"Reference","Panama","2020","deaths",1,0.02,0.778,"medium" +"Reference","Panama","2030","deaths",0,0.01,0.778,"medium" +"Reference","Panama","2040","deaths",1,0.01,0.778,"medium" +"Reference","Panama","2050","deaths",0,0,0.778,"medium" +"Reference","Paraguay","2020","deaths",659,10.02,0.633,"medium" +"Reference","Paraguay","2030","deaths",648,8.73,0.633,"medium" +"Reference","Paraguay","2040","deaths",650,8.02,0.633,"medium" +"Reference","Paraguay","2050","deaths",654,7.64,0.633,"medium" +"Reference","Peru","2020","deaths",1972,5.96,0.49,"high" +"Reference","Peru","2030","deaths",1638,4.55,0.49,"high" +"Reference","Peru","2040","deaths",1458,3.82,0.49,"high" +"Reference","Peru","2050","deaths",1355,3.44,0.49,"high" +"Reference","Philippines","2020","deaths",49835,44.78,63.075,"low" +"Reference","Philippines","2030","deaths",54941,42.89,63.075,"low" +"Reference","Philippines","2040","deaths",60058,42.04,63.075,"low" +"Reference","Philippines","2050","deaths",64269,41.64,63.075,"low" +"Reference","Poland","2020","deaths",3798,9.87,485.911,"low" +"Reference","Poland","2030","deaths",3397,8.89,485.911,"low" +"Reference","Poland","2040","deaths",3079,8.42,485.911,"low" +"Reference","Poland","2050","deaths",2875,8.19,485.911,"low" +"Reference","Portugal","2020","deaths",2,0.02,0.975,"medium" +"Reference","Portugal","2030","deaths",1,0.01,0.975,"medium" +"Reference","Portugal","2040","deaths",1,0.01,0.975,"medium" +"Reference","Portugal","2050","deaths",0,0,0.975,"medium" +"Reference","Puerto Rico","2020","deaths",51,1.57,4.219,"low" +"Reference","Puerto Rico","2030","deaths",26,0.88,4.219,"low" +"Reference","Puerto Rico","2040","deaths",13,0.5,4.219,"low" +"Reference","Puerto Rico","2050","deaths",7,0.29,4.219,"low" +"Reference","Qatar","2020","deaths",0,0.01,1,"medium" +"Reference","Qatar","2030","deaths",0,0.01,1,"medium" +"Reference","Qatar","2040","deaths",0,0,1,"medium" +"Reference","Qatar","2050","deaths",0,0,1,"medium" +"Reference","Romania","2020","deaths",2349,12.05,3.007,"low" +"Reference","Romania","2030","deaths",2012,10.99,3.007,"low" +"Reference","Romania","2040","deaths",1795,10.49,3.007,"low" +"Reference","Romania","2050","deaths",1630,10.23,3.007,"low" +"Reference","Russian Federation","2020","deaths",29,0.02,0.947,"medium" +"Reference","Russian Federation","2030","deaths",14,0.01,0.947,"medium" +"Reference","Russian Federation","2040","deaths",14,0.01,0.947,"medium" +"Reference","Russian Federation","2050","deaths",0,0,0.947,"medium" +"Reference","Rwanda","2020","deaths",9224,71.02,91.92,"low" +"Reference","Rwanda","2030","deaths",10853,67,91.92,"low" +"Reference","Rwanda","2040","deaths",12674,65.22,91.92,"low" +"Reference","Rwanda","2050","deaths",14451,64.36,91.92,"low" +"Reference","Saint Lucia","2020","deaths",10,5.43,5.609,"low" +"Reference","Saint Lucia","2030","deaths",7,3.8,5.609,"low" +"Reference","Saint Lucia","2040","deaths",6,3.06,5.609,"low" +"Reference","Saint Lucia","2050","deaths",5,2.65,5.609,"low" +"Reference","Saudi Arabia","2020","deaths",4,0.01,0.996,"medium" +"Reference","Saudi Arabia","2030","deaths",4,0.01,0.996,"medium" +"Reference","Saudi Arabia","2040","deaths",0,0,0.996,"medium" +"Reference","Saudi Arabia","2050","deaths",0,0,0.996,"medium" +"Reference","Senegal","2020","deaths",8212,50.64,2.287,"low" +"Reference","Senegal","2030","deaths",9839,47.68,2.287,"low" +"Reference","Senegal","2040","deaths",11718,46.22,2.287,"low" +"Reference","Senegal","2050","deaths",13671,45.51,2.287,"low" +"Reference","Serbia","2020","deaths",1863,25.23,37.48,"low" +"Reference","Serbia","2030","deaths",1682,23.88,37.48,"low" +"Reference","Serbia","2040","deaths",1586,23.26,37.48,"low" +"Reference","Serbia","2050","deaths",1519,22.96,37.48,"low" +"Reference","Seychelles","2020","deaths",3,2.61,25.88,"low" +"Reference","Seychelles","2030","deaths",2,1.54,25.88,"low" +"Reference","Seychelles","2040","deaths",1,1.03,25.88,"low" +"Reference","Seychelles","2050","deaths",1,0.77,25.88,"low" +"Reference","Sierra Leone","2020","deaths",3546,43.56,0.463,"high" +"Reference","Sierra Leone","2030","deaths",3954,39.55,0.463,"high" +"Reference","Sierra Leone","2040","deaths",4382,37.27,0.463,"high" +"Reference","Sierra Leone","2050","deaths",4793,36.11,0.463,"high" +"Reference","Singapore","2020","deaths",37,0.62,0.927,"medium" +"Reference","Singapore","2030","deaths",6,0.1,0.927,"medium" +"Reference","Singapore","2040","deaths",0,0,0.927,"medium" +"Reference","Singapore","2050","deaths",0,0,0.927,"medium" +"Reference","Slovakia","2020","deaths",1,0.02,0.951,"medium" +"Reference","Slovakia","2030","deaths",1,0.01,0.951,"medium" +"Reference","Slovakia","2040","deaths",1,0.01,0.951,"medium" +"Reference","Slovakia","2050","deaths",0,0,0.951,"medium" +"Reference","Slovenia","2020","deaths",94,4.46,4.677,"low" +"Reference","Slovenia","2030","deaths",74,3.54,4.677,"low" +"Reference","Slovenia","2040","deaths",64,3.1,4.677,"low" +"Reference","Slovenia","2050","deaths",58,2.86,4.677,"low" +"Reference","South Africa","2020","deaths",18,0.03,0.78,"medium" +"Reference","South Africa","2030","deaths",13,0.02,0.78,"medium" +"Reference","South Africa","2040","deaths",7,0.01,0.78,"medium" +"Reference","South Africa","2050","deaths",0,0,0.78,"medium" +"Reference","South Sudan","2020","deaths",3384,32.09,6.233,"low" +"Reference","South Sudan","2030","deaths",2933,25.04,6.233,"low" +"Reference","South Sudan","2040","deaths",2805,21.48,6.233,"low" +"Reference","South Sudan","2050","deaths",2701,19.6,6.233,"low" +"Reference","Spain","2020","deaths",852,1.8,0.037,"high" +"Reference","Spain","2030","deaths",464,0.97,0.037,"high" +"Reference","Spain","2040","deaths",262,0.54,0.037,"high" +"Reference","Spain","2050","deaths",153,0.31,0.037,"high" +"Reference","Sri Lanka","2020","deaths",7277,33.56,13.022,"low" +"Reference","Sri Lanka","2030","deaths",7281,32.18,13.022,"low" +"Reference","Sri Lanka","2040","deaths",7240,31.3,13.022,"low" +"Reference","Sri Lanka","2050","deaths",7152,30.85,13.022,"low" +"Reference","Sweden","2020","deaths",2,0.02,0.998,"medium" +"Reference","Sweden","2030","deaths",1,0.01,0.998,"medium" +"Reference","Sweden","2040","deaths",0,0,0.998,"medium" +"Reference","Sweden","2050","deaths",0,0,0.998,"medium" +"Reference","Switzerland","2020","deaths",1,0.01,0.999,"medium" +"Reference","Switzerland","2030","deaths",1,0.01,0.999,"medium" +"Reference","Switzerland","2040","deaths",0,0,0.999,"medium" +"Reference","Switzerland","2050","deaths",0,0,0.999,"medium" +"Reference","Syrian Arab Republic","2020","deaths",596,2.91,0.762,"medium" +"Reference","Syrian Arab Republic","2030","deaths",325,1.36,0.762,"medium" +"Reference","Syrian Arab Republic","2040","deaths",155,0.56,0.762,"medium" +"Reference","Syrian Arab Republic","2050","deaths",37,0.12,0.762,"medium" +"Reference","Taiwan (Province of China)","2020","deaths",2,0.01,0.947,"medium" +"Reference","Taiwan (Province of China)","2030","deaths",2,0.01,0.947,"medium" +"Reference","Taiwan (Province of China)","2040","deaths",0,0,0.947,"medium" +"Reference","Taiwan (Province of China)","2050","deaths",0,0,0.947,"medium" +"Reference","Tajikistan","2020","deaths",5,0.05,0.226,"high" +"Reference","Tajikistan","2030","deaths",3,0.03,0.226,"high" +"Reference","Tajikistan","2040","deaths",1,0.01,0.226,"high" +"Reference","Tajikistan","2050","deaths",1,0.01,0.226,"high" +"Reference","United Republic of Tanzania","2020","deaths",47433,78.05,12.479,"low" +"Reference","United Republic of Tanzania","2030","deaths",60034,74.71,12.479,"low" +"Reference","United Republic of Tanzania","2040","deaths",74346,72.99,12.479,"low" +"Reference","United Republic of Tanzania","2050","deaths",89629,72.16,12.479,"low" +"Reference","Thailand","2020","deaths",3505,4.91,0.25,"high" +"Reference","Thailand","2030","deaths",2642,3.66,0.25,"high" +"Reference","Thailand","2040","deaths",2137,3.01,0.25,"high" +"Reference","Thailand","2050","deaths",1821,2.68,0.25,"high" +"Reference","Togo","2020","deaths",4242,50.85,1.577,"medium" +"Reference","Togo","2030","deaths",4894,47.1,1.577,"medium" +"Reference","Togo","2040","deaths",5625,45.17,1.577,"medium" +"Reference","Togo","2050","deaths",6354,44.22,1.577,"medium" +"Reference","Trinidad and Tobago","2020","deaths",27,1.77,0.847,"medium" +"Reference","Trinidad and Tobago","2030","deaths",11,0.73,0.847,"medium" +"Reference","Trinidad and Tobago","2040","deaths",2,0.13,0.847,"medium" +"Reference","Trinidad and Tobago","2050","deaths",0,0,0.847,"medium" +"Reference","Tunisia","2020","deaths",301,2.49,0.721,"medium" +"Reference","Tunisia","2030","deaths",145,1.11,0.721,"medium" +"Reference","Tunisia","2040","deaths",35,0.26,0.721,"medium" +"Reference","Tunisia","2050","deaths",0,0,0.721,"medium" +"Reference","Turkey","2020","deaths",1869,2.23,1.862,"medium" +"Reference","Turkey","2030","deaths",1115,1.24,1.862,"medium" +"Reference","Turkey","2040","deaths",734,0.78,1.862,"medium" +"Reference","Turkey","2050","deaths",520,0.54,1.862,"medium" +"Reference","Turkmenistan","2020","deaths",168,2.7,0.228,"high" +"Reference","Turkmenistan","2030","deaths",112,1.57,0.228,"high" +"Reference","Turkmenistan","2040","deaths",67,0.85,0.228,"high" +"Reference","Turkmenistan","2050","deaths",39,0.46,0.228,"high" +"Reference","Uganda","2020","deaths",29837,68.3,16186.342,"low" +"Reference","Uganda","2030","deaths",37756,64.68,16186.342,"low" +"Reference","Uganda","2040","deaths",46756,62.85,16186.342,"low" +"Reference","Uganda","2050","deaths",56219,61.95,16186.342,"low" +"Reference","Ukraine","2020","deaths",4387,9.95,307.041,"low" +"Reference","Ukraine","2030","deaths",3165,8.57,307.041,"low" +"Reference","Ukraine","2040","deaths",2806,7.69,307.041,"low" +"Reference","Ukraine","2050","deaths",2552,7.3,307.041,"low" +"Reference","United Arab Emirates","2020","deaths",1,0.01,1,"medium" +"Reference","United Arab Emirates","2030","deaths",1,0.01,1,"medium" +"Reference","United Arab Emirates","2040","deaths",0,0,1,"medium" +"Reference","United Arab Emirates","2050","deaths",0,0,1,"medium" +"Reference","United Kingdom","2020","deaths",1044,1.56,0.711,"medium" +"Reference","United Kingdom","2030","deaths",570,0.83,0.711,"medium" +"Reference","United Kingdom","2040","deaths",310,0.44,0.711,"medium" +"Reference","United Kingdom","2050","deaths",160,0.22,0.711,"medium" +"Reference","United States of America","2020","deaths",4427,1.32,0.094,"high" +"Reference","United States of America","2030","deaths",2518,0.72,0.094,"high" +"Reference","United States of America","2040","deaths",1530,0.42,0.094,"high" +"Reference","United States of America","2050","deaths",905,0.24,0.094,"high" +"Reference","Uruguay","2020","deaths",1,0.02,0.96,"medium" +"Reference","Uruguay","2030","deaths",0,0.01,0.96,"medium" +"Reference","Uruguay","2040","deaths",0,0.01,0.96,"medium" +"Reference","Uruguay","2050","deaths",0,0,0.96,"medium" +"Reference","Uzbekistan","2020","deaths",4538,13.65,5.532,"low" +"Reference","Uzbekistan","2030","deaths",4585,11.9,5.532,"low" +"Reference","Uzbekistan","2040","deaths",4727,11.02,5.532,"low" +"Reference","Uzbekistan","2050","deaths",4965,10.6,5.532,"low" +"Reference","Venezuela (Bolivarian Republic of)","2020","deaths",6,0.02,0.974,"medium" +"Reference","Venezuela (Bolivarian Republic of)","2030","deaths",3,0.01,0.974,"medium" +"Reference","Venezuela (Bolivarian Republic of)","2040","deaths",3,0.01,0.974,"medium" +"Reference","Venezuela (Bolivarian Republic of)","2050","deaths",0,0,0.974,"medium" +"Reference","Viet Nam","2020","deaths",39905,41.48,96.092,"low" +"Reference","Viet Nam","2030","deaths",40774,39.7,96.092,"low" +"Reference","Viet Nam","2040","deaths",41351,38.93,96.092,"low" +"Reference","Viet Nam","2050","deaths",41499,38.58,96.092,"low" +"Reference","Yemen","2020","deaths",6481,20.3,0.588,"medium" +"Reference","Yemen","2030","deaths",6890,17.39,0.588,"medium" +"Reference","Yemen","2040","deaths",7315,15.55,0.588,"medium" +"Reference","Yemen","2050","deaths",7853,14.63,0.588,"medium" +"Reference","Zambia","2020","deaths",4802,25.74,0.523,"medium" +"Reference","Zambia","2030","deaths",5510,22.89,0.523,"medium" +"Reference","Zambia","2040","deaths",6320,21.3,0.523,"medium" +"Reference","Zambia","2050","deaths",7153,20.47,0.523,"medium" +"Reference","Zimbabwe","2020","deaths",6929,44.69,0.914,"medium" +"Reference","Zimbabwe","2030","deaths",7520,40.86,0.914,"medium" +"Reference","Zimbabwe","2040","deaths",8166,38.85,0.914,"medium" +"Reference","Zimbabwe","2050","deaths",8739,37.8,0.914,"medium" diff --git a/tests/testthat/testOutputs/calc_hap_impacts_output1.RData b/tests/testthat/testOutputs/calc_hap_impacts_output1.RData index 81ce451..9b1db51 100644 Binary files a/tests/testthat/testOutputs/calc_hap_impacts_output1.RData and b/tests/testthat/testOutputs/calc_hap_impacts_output1.RData differ diff --git a/tests/testthat/testOutputs/fit_model.RData b/tests/testthat/testOutputs/fit_model.RData deleted file mode 100644 index 5ef9b55..0000000 Binary files a/tests/testthat/testOutputs/fit_model.RData and /dev/null differ diff --git a/vignettes/fit_model.Rmd b/vignettes/fit_model.Rmd index 05d337e..40bcc6e 100644 --- a/vignettes/fit_model.Rmd +++ b/vignettes/fit_model.Rmd @@ -26,21 +26,42 @@ Finally, the panel includes a set of climatic variables, such as average maximum * Cooling Degree Days (CDD): representing days when cooling is required, calculated similarly to HDD but using outdoor high temperatures instead of low ones. Data for these indicators was sourced from the [Community Climate System Model](https://www.cesm.ucar.edu/models/ccsm). ### Econometric regression -While the panel provides a wide range of variables for the regression model, not all are statistically significant. Moreover, some are highly correlated, making it essential to select covariates carefully to ensure independence and avoid multicollinearity. After evaluating multiple models with different variable combinations, the final regression includes the following socioeconomic and environmental variables: per capita GDP, emissions of primary PM2.5 (comprising BC and OC), nitrogen oxides, and non-methane volatile organic compounds (NMVOCs), and per capita floor space. +While the panel provides a wide range of variables for the regression model, not all are statistically significant. Moreover, some are highly correlated, making it essential to select covariates carefully to ensure independence and avoid multicollinearity. After evaluating multiple models with different variable combinations, the final regression includes the following socioeconomic and environmental variables: per capita GDP, emissions of primary PM2.5 (comprising BC and OC), nitrogen oxides, and per capita floor space, plus a linear year trend. Non-methane volatile organic compounds (NMVOCs) were tested and dropped: their coefficient is statistically indistinguishable from zero under every specification and standard-error estimator tried, and fails a joint significance test both alone and jointly with PM2.5. Once the covariates are selected, the variables are first normalized by converting them to units per 100,000 inhabitants. This step ensures that the scales of the variables are consistent, making them more comparable across different regions or populations. Subsequently, a logarithmic transformation is applied to these normalized variables. The logarithmic scale helps linearize any nonlinear relationships, addresses issues related to skewness, and further stabilizes variance, ultimately improving the reliability and interpretability of the regression model. -To determine whether to use a fixed or random effects model, we apply the Hausmann test. Since the null hypothesis is rejected at conventional significance levels, this suggests that unobservable individual-specific characteristics are correlated with the explanatory variables. As a result, we opt for a fixed effects model over a random effects model. Note that the estimation is performed using the [plm] (https://cran.r-project.org/web/packages/plm/plm.pdf) package. +To determine whether to use a fixed or random effects model, we apply the Hausman test. Since the null hypothesis is rejected at conventional significance levels, this suggests that unobservable individual-specific characteristics are correlated with the explanatory variables. As a result, we opt for a fixed effects model over a random effects model. Note that the estimation is performed using the [plm](https://cran.r-project.org/web/packages/plm/plm.pdf) package. -In summary, we estimate health impacts attributable to household air pollution using the following fixed effects model. The model includes logarithmic and normalized per capita GDP, emissions of primary PM2.5, NOx, and VOCs, and per capita floorspace as covariates, with country and year specified as index variables: +#### Time trend -$$log(Deathsper100K)_{i,t} = \alpha_{i} + \beta_{1}log(PrimPM2.5per100K)_{i,t} + \beta_{2}log(NOXper100K)_{i,t} + \beta_{3}log(VOCper100K)_{i,t} + \beta_{4}log(GDPpc)_{i,t} + \beta_{5}log(FLSPpc)_{i,t} + e_{i,t} $$ +Mortality attributable to HAP has an underlying downward trend over 1990-2019 that is largely common across countries (driven by factors outside the model, such as general improvements in medicine and public health infrastructure). Leaving this trend unmodeled risks confounding it with the genuinely declining trend in residential pollutant emissions over the same period. We compared four ways of controlling for it: no trend control, a linear year trend, a flexible spline (4 degrees of freedom), and full year fixed effects (year dummies). The linear trend, spline, and year-dummy specifications all converge to similar pollutant coefficients; only the no-trend-control specification disagrees, indicating that the "reasonable-looking" coefficients from that specification were actually the biased ones (omitted common-trend bias). Full year dummies over-absorb the common trend to the point of washing out the pollutant coefficients almost entirely. We use the linear trend: it is the most parsimonious option that agrees with the more flexible alternatives, and unlike year dummies it extrapolates sensibly to future (post-2019) scenario years. + +#### Inference + +Country-year panels of this kind routinely violate the classical OLS/within-estimator assumptions used for standard errors: residuals from the same country in adjacent years are correlated (serial correlation), and global shocks (oil prices, recessions, multilateral pollution treaties, pandemics) hit many countries in the same year (cross-sectional dependence). We tested for both directly: a Breusch-Godfrey/Wooldridge test rejects the no-serial-correlation null (p < 0.001), and a Pesaran CD test rejects the no-cross-sectional-dependence null (p < 0.001) -- so both are actually present in this panel, not just theoretical concerns. Naive within-estimator standard errors, and even standard errors clustered only by country (Arellano), are not valid under cross-sectional dependence. We therefore report Driscoll-Kraay standard errors (`plm::vcovSCC()`), which are robust to serial correlation, cross-sectional dependence, and heteroskedasticity (also confirmed present via a Breusch-Pagan test) simultaneously. + +In summary, we estimate health impacts attributable to household air pollution using the following fixed effects model, with country and year specified as index variables: + +$$log(Deathsper100K)_{i,t} = \alpha_{i} + \beta_{1}log(PrimPM2.5per100K)_{i,t} + \beta_{2}log(NOXper100K)_{i,t} + \beta_{3}log(GDPpc)_{i,t} + \beta_{4}log(FLSPpc)_{i,t} + \gamma \cdot t + e_{i,t} $$ Where: - $log(Deathsper100K)_{i,t}$ is the dependent variable for region $i$ at period $t$. - $\alpha_{i}$ is the fixed effect for region $i$. -- $\beta_{1}$, $\beta_{2}$, ... $\beta_{5}$ are the coefficients for the corresponding explanatory variables. +- $\beta_{1}$, ... $\beta_{4}$ are the coefficients for the corresponding explanatory variables. +- $\gamma$ is the linear year-trend coefficient, with $t$ the calendar year. - $e_{i,t}$ is the error term. ** Note that the dependent variable can also be set to be YLLs or DALYs + +#### Prediction and the bias adder + +`predict()` on a fitted `plm` "within" model returns fitted values from the coefficients alone, without the country fixed effect $\alpha_i$ -- and this is unavoidable for scenario prediction, since GCAM scenario data has no fixed effect to draw on. To correct for this, `fit_model()` also computes a per-country **additive bias adder**: the mean of (observed $-$ coefficient-only-predicted) values, in per-100k units, over the most recent years of the training panel. This is added to (not multiplied into) every future per-100k prediction for that country. + +This choice was not obvious going in. In a log-linear model, $log(Y) = X\beta + \alpha_i$ implies $Y = \exp(X\beta) \cdot \exp(\alpha_i)$ -- the fixed effect is technically a *multiplicative* scalar in levels, not an additive one, which argues for a multiplicative (geometric-mean ratio) correction instead. We implemented that version and validated it against real GBD data before adopting it -- and it performed far worse: cross-country correlation with observed 2023 GBD rates collapsed from ~0.86 (additive) to ~0.18 (multiplicative), even though the multiplicative version looked more accurate on the model's own training data. The reason is that the model's single common linear year trend doesn't match every country's own trajectory (mortality in a country like Uganda declines far more slowly than the trend assumes), so the gap between the trend-implied prediction and reality keeps *growing* even within the training years -- for about 22% of countries the implied fixed effect exceeds a factor of 10x, and for some it exceeds 1000x. A multiplicative correction amplifies that already-large, still-growing gap catastrophically the moment it's applied to a slightly different (scenario) year; a correction bounded by the actual scale of the outcome variable does not. We therefore use the additive form despite its weaker theoretical justification, because it is dramatically more reliable in the use this package is actually built for. + +The one real failure mode of the additive form is that it can push a country's corrected rate below zero (this happened for ~12% of country-year observations on the model's own training data). Rather than clamping that to exactly 0 -- which silently reports "zero risk" for a country the model has no real basis to call zero -- `calc_hap_impacts()` floors the corrected value at a small fraction (1%) of the naive, uncorrected prediction, so a low but non-degenerate estimate is reported instead. + +#### Known limitations + +A RESET test on the pooled specification rejects correct functional form (p < 0.001), and individual country series for deaths, GDP, and floorspace show unit-root-like (non-stationary) behavior under an intercept-only panel unit-root test -- both expected for 30-year macro/demographic series, but neither is fully resolved by entity fixed effects plus a single *common* linear trend if countries have meaningfully different trends of their own. This is a standard trade-off in applied fixed-effects panel work with many countries and a short time dimension (T=30), and the model's cross-country ranking has been validated against real GBD data (~0.86 correlation on out-of-sample deaths/YLL/DALY rates); but formal cointegration testing or country-specific trends were not pursued and would be a reasonable next step for anyone using the model's point coefficients for causal interpretation rather than for scenario comparison. diff --git a/vignettes/model_validation.Rmd b/vignettes/model_validation.Rmd new file mode 100644 index 0000000..50c4c69 --- /dev/null +++ b/vignettes/model_validation.Rmd @@ -0,0 +1,200 @@ +--- +title: "Validating the econometric model" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{model_validation} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + eval = FALSE +) +``` + +### Why this vignette exists + +The ["Fitting the econometric model"](https://bc3lc.github.io/rhap/articles/fit_model.html) vignette explains *what* `fit_model()` estimates and *why* it's specified the way it is. This vignette is different: it's a reproducible record of *how to check the model is still healthy* -- the diagnostic tests, the real-data validation, and the specification tests that were run (and, in one case, rejected) while hardening the model. Re-run this workflow after any change to `fit_model()`, `panel_data`, or the bias-adder logic in `calc_hap_impacts()`, and before trusting the model's absolute-level predictions for a new use case. + +All code below is runnable as written from the package root (`devtools::load_all()` first). Chunks are marked `eval = FALSE` by default because a couple of steps need a real GCAM project file and, for the real-data validation section, an external GBD data export -- see the notes in each section for what to supply. + +### 1. Panel structure and multicollinearity + +```{r} +library(dplyr) +library(plm) + +data <- rhap::panel_data %>% + dplyr::select(iso, country_name, year, pop, dplyr::starts_with("log"), continent, dev) %>% + dplyr::mutate(year = as.character(year)) %>% + dplyr::select(-log_AAP, -log_HDD_value, -log_CDD_value) %>% + dplyr::filter(stats::complete.cases(.)) %>% + dplyr::mutate(year_num = as.numeric(year)) + +plm::pdim(data, index = c("country_name", "year")) +# Balanced Panel: n = 165, T = 30, N = 4950 -- no singleton entities/years, +# which matters because singletons are the most common cause of plm's +# "empty model" error under effect = "twoways" (not used here, but worth +# checking before anyone tries it). + +model_formula <- log_Deaths_per_100k ~ log_PrimPM25_per_100k + log_NOx_per_100k + + log_gdppc_ppp_dol2011 + log_flsp + year_num + +car::vif(lm(model_formula, data = data)) +#> log_PrimPM25_per_100k log_NOx_per_100k log_gdppc_ppp_dol2011 +#> 1.768138 1.549074 2.473823 +#> log_flsp year_num +#> 2.226533 1.071704 +# All comfortably < 3 -- no multicollinearity concern. +``` + +### 2. Model fit and inference diagnostics + +```{r} +model.fixed <- plm::plm(model_formula, data = data, index = c("country_name", "year"), + model = "within", effect = "individual") + +plm::pbgtest(model.fixed) +#> Breusch-Godfrey/Wooldridge test for serial correlation in panel models +#> chisq = 4152.2, df = 30, p-value < 2.2e-16 + +plm::pcdtest(model.fixed, test = "cd") +#> Pesaran CD test for cross-sectional dependence in panels +#> z = 22.223, p-value < 2.2e-16 + +lmtest::bptest(model_formula, data = data, studentize = TRUE) +#> studentized Breusch-Pagan test +#> BP = 720.64, df = 5, p-value < 2.2e-16 +``` + +Serial correlation, cross-sectional dependence, and heteroskedasticity are all rejected at p < 2.2e-16 -- not just theoretical concerns for a country-year panel, but *confirmed present* in this one. That's why `fit_model()` reports Driscoll-Kraay standard errors (`plm::vcovSCC(model.fixed, type = "HC1", maxlag = 4)`, robust to all three simultaneously) rather than naive or country-clustered (Arellano) SEs -- clustering alone doesn't fix cross-sectional dependence. + +### 3. Functional form: RESET test and the GDP x PM2.5 interaction + +```{r} +lmtest::resettest(lm(model_formula, data = data), power = 2:3, type = "fitted") +#> RESET = 792.19, df1 = 2, df2 = 4942, p-value < 2.2e-16 +``` + +A significant RESET test flags omitted nonlinearity. The natural hypothesis: richer countries might see less marginal harm per unit of PM2.5 (better ventilation, indoor air filtration), i.e. a GDP x PM2.5 interaction. Testing this requires **mean-centering** both variables first -- an uncentered interaction changes what the "main effect" coefficients mean (the effect when the other regressor is exactly zero, i.e. `log(x) = 0`, nowhere near the data), which produces wildly unstable, uninterpretable coefficients that look like a multicollinearity problem but aren't: + +```{r} +data_c <- data %>% + dplyr::mutate( + pm25_c = log_PrimPM25_per_100k - mean(log_PrimPM25_per_100k), + gdp_c = log_gdppc_ppp_dol2011 - mean(log_gdppc_ppp_dol2011) + ) + +m_int <- plm::plm( + log_Deaths_per_100k ~ pm25_c + log_NOx_per_100k + gdp_c + log_flsp + year_num + pm25_c:gdp_c, + data = data_c, index = c("country_name", "year"), model = "within", effect = "individual" +) +vcov_int <- plm::vcovSCC(m_int, type = "HC1", maxlag = 4) +lmtest::coeftest(m_int, vcov = vcov_int) +``` + +Result: the interaction is statistically significant (p < 2.2e-16), not a within-transformation collinearity artifact (auxiliary-regression VIF on the within-demeaned regressors is 1.02-1.09, not inflated), improves within-R² from 0.69 to 0.76, and reduces the RESET statistic from 792 to 438 -- all of which would normally argue for adopting it. + +**It was rejected anyway.** The marginal effect of PM2.5, evaluated at different GDP levels, is: + +| GDP percentile | Marginal effect of PM2.5 | +|---|---| +| 10th (poorest) | -0.50 | +| 50th (median) | -0.07 | +| 90th (richest) | +0.29 | + +This is backwards from the hypothesis (worse ventilation in poor countries should mean a *larger* positive marginal effect there, not a negative one), and violates the basic prior that more emissions should not translate into fewer deaths -- exactly in the low-income countries where household air pollution matters most. A statistically significant, well-identified relationship that fails an economic-plausibility check is not necessarily wrong, but it's not something to ship without understanding the mechanism, and none was found (the most likely explanation is that the interaction is partly proxying for within-country GDP-correlated time trends the model doesn't otherwise capture, rather than a genuine effect-modification story). **The default model does not include this interaction.** If you want to explore it further, `pm25_c`/`gdp_c` centering and `car::vif()` on the within-demeaned regressors (not the raw pooled data, which can look fine even when the demeaned identifying variation is collinear) are the two steps not to skip. + +### 4. The bias-adder: additive vs. multiplicative, and why it matters + +`predict()` on a fitted `plm` "within" model omits the entity fixed effect entirely (unavoidable for scenario prediction -- GCAM data has no fixed effect to draw on). `fit_model()` corrects for this with a per-country bias adder, calibrated on the model's own training data. Two functional forms were tested: + +```{r} +train_panel <- plm::pdata.frame(data, index = c("country_name", "year")) +data$pred_log <- stats::predict(model.fixed, train_panel) +data$observed_per_100k <- exp(data$log_Deaths_per_100k) +data$pred_per_100k <- exp(data$pred_log) +data$resid_per_100k <- data$observed_per_100k - data$pred_per_100k +data$log_ratio <- data$log_Deaths_per_100k - data$pred_log + +n_years <- 5 +bias_additive <- data %>% + dplyr::group_by(country_name) %>% + dplyr::filter(year_num %in% utils::tail(sort(unique(year_num)), n_years)) %>% + dplyr::summarise(bias.adder = mean(resid_per_100k, na.rm = TRUE), .groups = "drop") + +bias_multiplicative <- data %>% + dplyr::group_by(country_name) %>% + dplyr::filter(year_num %in% utils::tail(sort(unique(year_num)), n_years)) %>% + dplyr::summarise(bias.factor = exp(mean(log_ratio, na.rm = TRUE)), .groups = "drop") +``` + +**In-sample (training data), the multiplicative version looks strictly better**: applying each correction back to the same recent training years, the additive version produces a negative "corrected" rate for 99 of 825 rows (12%), which then gets clamped; the multiplicative version produces zero negative values by construction (it's a ratio of two positive quantities), and has lower MAE (1.10 vs. 1.54 deaths/100k). + +**Out-of-sample, against real GBD data, the multiplicative version is far worse** -- see section 5 for the methodology, but the headline: cross-country correlation with observed GBD rates collapses from **0.86 (additive) to 0.18 (multiplicative)**. The reason: `bias_adder$bias.factor` (or the multiplicative version's equivalent) is unbounded, and for ~22% of countries it exceeds 10x; for Uganda specifically it's ~16,000x, because the model's naive (fixed-effect-free) prediction for Uganda is essentially zero (~0.02 deaths/100k) against an actual rate of ~120 -- the model's common linear year trend assumes mortality declines at a rate Uganda's doesn't follow, so the gap between the trend-implied prediction and reality *keeps growing* even within the training years. Multiplying by a correction that large amplifies any small difference between GCAM-scenario-implied inputs and the training panel's real-world inputs catastrophically; an additive correction, bounded by the actual scale of the outcome variable, doesn't. `fit_model()` uses the additive form for exactly this reason, with a floor (1% of the naive prediction) in `calc_hap_impacts()` instead of a hard clamp to 0, to avoid reporting a discontinuous "zero risk" for countries the correction pushes negative. + +`fit_model()`'s `bias_adder` also returns `reliability_ratio` (`abs(bias.adder) / naive_prediction`, same calibration window) precisely so this asymmetry is visible downstream: a country with a `reliability` of `"low"` in `calc_hap_impacts()`'s output is one whose absolute-level prediction rests almost entirely on this correction rather than the regression's own covariates -- not necessarily wrong, but not independently corroborated by the model either. + +### 5. Real-data validation against GBD + +`panel_data` (and therefore the model's training years) tops out at 2019 (GBD 2019-release data). To check the model against something it wasn't fit on, compare its scenario predictions against a newer GBD data release. + +**Getting the data**: export from the [GBD Results Tool](https://vizhub.healthdata.org/gbd-results/) with cause = "All causes", risk = "Household air pollution from solid fuels", measure = Deaths/YLLs/DALYs, metric = Rate, sex = Both, location = all countries, and whichever year(s) are available in the release you're checking against. The comparison below used a GBD 2023-release export for year 2023. + +```{r} +gbd <- read.csv("path/to/your/GBD_export.csv") %>% + dplyr::filter(sex_name == "Both", metric_name == "Rate") %>% + dplyr::mutate( + HIA_var = dplyr::case_when( + measure_name == "Deaths" ~ "deaths", + measure_name == "YLLs (Years of Life Lost)" ~ "yll", + measure_name == "DALYs (Disability-Adjusted Life Years)" ~ "dalys" + ), + iso3 = countrycode::countrycode(location_name, "country.name", "iso3c", warn = FALSE) + ) %>% + dplyr::filter(!is.na(iso3)) %>% + dplyr::select(iso3, HIA_var, gbd_rate = val) + +# Run against the package's bundled test GCAM project at its most recent +# calibrated (non-projected) year, so the model's inputs are as close to +# real-world 2019-2020 conditions as this test fixture gets. +model_out <- dplyr::bind_rows(lapply(c("deaths", "yll", "dalys"), function(hv) { + calc_hap_impacts( + prj_name = "tests/testthat/testInputs/test_prj_v7p1.dat", scen_name = "Reference", + final_db_year = 2020, HIA_var = hv, + saveOutput = FALSE, map = FALSE, anim = FALSE, normalized = FALSE, by_gr = FALSE + ) %>% dplyr::mutate(HIA_var = hv) +})) %>% + dplyr::filter(year == 2020) %>% + dplyr::mutate(iso3 = countrycode::countrycode(country, "country.name", "iso3c", warn = FALSE)) + +comparison <- model_out %>% dplyr::inner_join(gbd, by = c("iso3", "HIA_var")) + +comparison %>% + dplyr::group_by(HIA_var) %>% + dplyr::summarise( + n_countries = dplyr::n(), + cor = stats::cor(pred_value_normalized, gbd_rate, use = "complete.obs"), + mae = mean(abs(pred_value_normalized - gbd_rate)), + .groups = "drop" + ) +``` + +**Important caveat**: this compares the model's 2020 prediction against 2023 real-world data -- a 3-year gap, plus the training panel itself is built from an earlier GBD release, so the historical estimates it was fit on don't exactly match the 2023 release's methodology either. Treat this as a directional check (does the model still track reality, in the right ballpark and rank order?), not a precise backtest. + +**Results (2020 model vs. 2023 GBD, after fixing the bias-adder to additive)**: + +| HIA_var | n countries | correlation | MAE (per 100k) | +|---|---|---|---| +| deaths | 160 | 0.855 | 11.9 | +| yll | 160 | 0.868 | 453 | +| dalys | 160 | 0.865 | 483 | + +Correlation in the mid-0.8s across all three outcomes means the model reproduces the relative cross-country pattern reasonably well (e.g. it correctly identifies North Korea, Myanmar, Haiti, and several Sub-Saharan African countries as the highest-burden countries, with predictions in the right order of magnitude). The remaining gap is concentrated in countries with a `reliability = "low"` flag (see section 4) -- their absolute levels are the least trustworthy part of any given prediction, which is exactly what that flag is for. + +### Reproducing this vignette's own numbers + +The commands above were run once, standalone, and their output pasted in as static results rather than executed at vignette build time (`eval = FALSE`) -- rebuilding `panel_data` from raw source files and re-running Driscoll-Kraay inference on the full panel takes a few minutes, too slow for a routine `pkgdown` build. Uncomment the chunks and run them directly (from the package root, after `devtools::load_all()`) when actually validating a change to the model. diff --git a/vignettes/run_rhap.Rmd b/vignettes/run_rhap.Rmd index 6002f27..70d632b 100644 --- a/vignettes/run_rhap.Rmd +++ b/vignettes/run_rhap.Rmd @@ -6,7 +6,7 @@ vignette: > %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- -The primary function provided by this package is `calc_hap_impacts`, which is designed to estimate health impacts attributable to household air pollution (HAP) under various alternative scenarios simulated using the Global Change Analysis Model (GCAM). This function streamlines the process of assessing scenario-specific health outcomes by extracting relevant data from GCAM databases or project files. These relevant data includes key socioeconomic and environmental parameters, such as per capita GDP, emissions of primary PM2.5, NOx, and VOCs, as well as per capita floorspace. For more details on the methodology and assumptions underlying the econometric model used in this calculation, refer to the guide on ["Fitting the Econometric Model"](https://bc3lc.github.io/rhap/articles/fit_model.html). +The primary function provided by this package is `calc_hap_impacts`, which is designed to estimate health impacts attributable to household air pollution (HAP) under various alternative scenarios simulated using the Global Change Analysis Model (GCAM). This function streamlines the process of assessing scenario-specific health outcomes by extracting relevant data from GCAM databases or project files. These relevant data includes key socioeconomic and environmental parameters, such as per capita GDP, emissions of primary PM2.5 and NOx, and per capita floorspace. For more details on the methodology and assumptions underlying the econometric model used in this calculation, refer to the guide on ["Fitting the Econometric Model"](https://bc3lc.github.io/rhap/articles/fit_model.html). ### Extraction of the scenario-specific data The function `calc_hap_impacts` connects to a GCAM databases or project file using the [rgcam](https://github.com/JGCRI/rgcam) package, designed to extract and import GCAM results. @@ -19,9 +19,9 @@ For per capita GDP, the package loads projections from the [SSP database version ### Calculations -The function estimates health impacts attributable to household air pollution (HAP) using an econometric model that incorporates direct emissions from the residential sector (BC, OC, NOx, and VOC), per capita GDP, and per capita floorspace as covariates. The model is a fixed-effects model, estimated using cross-regional and multi-year panel data compiled from various different sources. Detailed information on the econometric model, including its formulation and data sources, is available in the dedicated vignette: ["Fitting the Econometric Model"](https://bc3lc.github.io/rhap/articles/fit_model.html). +The function estimates health impacts attributable to household air pollution (HAP) using an econometric model that incorporates direct emissions from the residential sector (primary PM2.5 and NOx), per capita GDP, per capita floorspace, and a linear year trend as covariates. The model is a fixed-effects model, estimated using cross-regional and multi-year panel data compiled from various different sources. Detailed information on the econometric model, including its formulation, variable selection, and diagnostic testing, is available in the dedicated vignette: ["Fitting the Econometric Model"](https://bc3lc.github.io/rhap/articles/fit_model.html). -To improve accuracy, the health impact predictions incorporate a country-specific "bias adder" parameter, calculated as the difference between observed and estimated values in the final observed year (2019). This adjustment helps align model estimates with observed data for each country. +To improve accuracy, the health impact predictions incorporate a country-specific, **additive** bias-correction term, computed by `fit_model()` as the mean gap between observed and model-predicted values over the most recent years of the training panel, and added to every prediction. This corrects for the country fixed effect that a fitted `plm` model cannot otherwise apply to new (scenario) data. An additive correction was chosen over the more theoretically obvious multiplicative one after validation against real GBD data showed the multiplicative version performing far worse (see the "Prediction and the bias adder" section of the econometric-model vignette) -- the trade-off is that a small floor is applied to prevent a negative corrected rate. The package provides flexibility to calculate three distinct health impact metrics associated with HAP: @@ -29,7 +29,7 @@ The package provides flexibility to calculate three distinct health impact metri * Years of Life Lost (YLLs) * Disability-Adjusted Life Years (DALYs) -In addition to country-level estimates, the package offers an optional feature to estimate health impacts by socioeconomic group (e.g., income deciles) within each region. While it is recommended to perform country-level calculations for consistency, this capability can provide valuable insights into the distribution of health impacts across population segments within regions. Such granularity can aid users in understanding intra-regional disparities. This feature is activated by setting `by_gr = FALSE`. +In addition to country-level estimates, setting `by_gr = TRUE` reports within-country income-decile group **shares** of population, PM2.5+NOx exposure, and GDP, as a CSV table and a pie-chart snapshot (`output/by_gr/`). This is descriptive only, not a group-level health impact estimate: the regression model is fit on between-country variation, and applying it to within-country decile-level inputs would extrapolate it far outside the variation it was estimated on, with no way to validate the result. For a more detailed, per-pollutant breakdown of group shares, see `calc_ResidEm_grp()`. ### Outputs