diff --git a/.github/workflows/build.yaml b/.github/workflows/build.yaml index d7cb75e..f2477e6 100644 --- a/.github/workflows/build.yaml +++ b/.github/workflows/build.yaml @@ -1,49 +1,74 @@ name: build - on: push: branches: [main, k_dev] pull_request: branches: [main, k_dev] - jobs: R-CMD: runs-on: ${{ matrix.config.os }} name: ${{ matrix.config.os }} (${{ matrix.config.r }}) - strategy: fail-fast: false matrix: config: - { os: ubuntu-latest, r: 'devel', cran: "https://demo.rstudiopm.com/all/__linux__/xenial/latest"} - env: R_REMOTES_NO_ERRORS_FROM_WARNINGS: true GITHUB_PAT: ${{ secrets.GH_RHAP_TOKEN}} - steps: - - uses: actions/checkout@v1 - - uses: r-lib/actions/setup-r@v2-branch + - uses: actions/checkout@v4 + + - uses: r-lib/actions/setup-r@v2 with: - r-version: '4.1' - - uses: r-lib/actions/setup-pandoc@v2-branch - - uses: r-lib/actions/setup-tinytex@v2-branch + r-version: ${{ matrix.config.r }} + + - uses: r-lib/actions/setup-pandoc@v2 - - name: Install libraries + - uses: r-lib/actions/setup-tinytex@v2 + + - name: Install system libraries if: startsWith(matrix.config.os, 'ubuntu') run: | - sudo apt-add-repository ppa:ubuntugis/ubuntugis-unstable sudo apt-get update - 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 + 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 \ + libfontconfig1-dev libfreetype6-dev libpng-dev libtiff5-dev libjpeg-dev \ + libuv1-dev - - name: Install dependencies + - name: Install R dependencies run: | - 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") - remotes::install_deps(dependencies = TRUE) + options(Ncpus = parallel::detectCores()) + options(repos = c(CRAN = "https://cloud.r-project.org")) + install.packages(c("remotes", "rcmdcheck", "covr", "fs", "usethis", "testthat") + remotes::install_deps(dependencies = TRUE, upgrade = "never") + shell: Rscript {0} + + - name: Install the package + run: remotes::install_local(".") shell: Rscript {0} + # - name: Install dependencies + # run: | + # install.packages(c("remotes", "rcmdcheck", "stringi", "devtools"), 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") + # remotes::install_deps(dependencies = TRUE) + # shell: Rscript {0} + # + # - name: Install GitHub packages + # run: | + # remotes::install_github("JGCRI/rmap") + # remotes::install_github("JGCRI/gcamdata") + # shell: Rscript {0} + - name: Check run: rcmdcheck::rcmdcheck(args = c("--no-manual", "--ignore-vignettes", "--no-build-vignettes"), build_args = c("--no-manual", "--ignore-vignettes", "--no-build-vignettes"), error_on = "error") shell: Rscript {0} + + - name: Run tests + run: testthat::test() + shell: Rscript {0} diff --git a/.github/workflows/test_coverage.yaml b/.github/workflows/test_coverage.yaml index 052dcca..3af57f0 100644 --- a/.github/workflows/test_coverage.yaml +++ b/.github/workflows/test_coverage.yaml @@ -1,26 +1,21 @@ name: test_coverage - on: push: branches: [main] pull_request: branches: [main] - jobs: R-CMD: runs-on: ${{ matrix.config.os }} name: ${{ matrix.config.os }} (${{ matrix.config.r }}) - strategy: fail-fast: false matrix: config: - - { os: ubuntu-latest, r: 'devel', cran: "https://demo.rstudiopm.com/all/__linux__/xenial/latest" } - + - { os: ubuntu-latest, r: 'release' } env: R_REMOTES_NO_ERRORS_FROM_WARNINGS: true GITHUB_PAT: ${{ secrets.GH_RHAP_TOKEN }} - steps: - uses: actions/checkout@v4 @@ -29,27 +24,31 @@ jobs: r-version: ${{ matrix.config.r }} - uses: r-lib/actions/setup-pandoc@v2 + - uses: r-lib/actions/setup-tinytex@v2 - name: Install system libraries if: startsWith(matrix.config.os, 'ubuntu') run: | - sudo apt-add-repository ppa:ubuntugis/ubuntugis-unstable sudo apt-get update 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 \ + libfontconfig1-dev libfreetype6-dev libpng-dev libtiff5-dev libjpeg-dev \ + libuv1-dev - name: Install R dependencies run: | - install.packages(c("remotes", "rcmdcheck", "covr"), repos = "https://cloud.r-project.org") - remotes::install_deps(dependencies = TRUE) + options(Ncpus = parallel::detectCores()) + options(repos = c(CRAN = "https://cloud.r-project.org")) + install.packages(c("remotes", "rcmdcheck", "covr", "fs", "usethis", "testthat") + remotes::install_deps(dependencies = TRUE, upgrade = "never") shell: Rscript {0} - name: Install the package - run: | - remotes::install_local(".") + run: remotes::install_local(".") shell: Rscript {0} - name: Run tests and generate coverage @@ -60,7 +59,6 @@ jobs: - uses: codecov/codecov-action@v4.0.1 with: - token: ${{ secrets.CODECOV_TOKEN }} # Optional if repo is public + token: ${{ secrets.CODECOV_TOKEN }} files: ./coverage.xml verbose: true - diff --git a/.gitignore b/.gitignore index eaec7e4..d359fa1 100644 --- a/.gitignore +++ b/.gitignore @@ -8,3 +8,7 @@ R/*.txt /tests/testthat/testInputs/database_basexdb_ref/* /tests/testthat/output inst/doc +/output_esp +plot* +*pptx +*check* diff --git a/DESCRIPTION b/DESCRIPTION index 7d956ad..c72165f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -19,7 +19,7 @@ Depends: R (>= 3.5) License: MIT + file LICENSE Encoding: UTF-8 LazyData: true -RoxygenNote: 7.3.2 +RoxygenNote: 7.3.3 Imports: usethis, ggplot2, diff --git a/NAMESPACE b/NAMESPACE index 135eefe..40cbf5e 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,9 +2,12 @@ export(calc_ResidEm_grp) export(calc_hap_impacts) +export(calc_hap_impacts_predict) +export(calc_hap_impacts_preprocess) export(calc_hia_adder) export(create_panel) export(fit_model) export(get_gdp_ctry) export(get_pop_ctry) +export(listYears) importFrom(magrittr,"%>%") diff --git a/R/calc_ResidEm_grp.R b/R/calc_ResidEm_grp.R index f50b1b5..157b99b 100644 --- a/R/calc_ResidEm_grp.R +++ b/R/calc_ResidEm_grp.R @@ -34,26 +34,6 @@ calc_ResidEm_grp <- function(db_path = NULL, query_path = "./inst/extdata", db_n if (length(pollutant) != 1) stop("Error: Please provide a single pollutant as input to the `calc_ResidEm_grp` function.") if (length(region) != 1) stop("Error: Please provide a single region as input to the `calc_ResidEm_grp` function.") - year <- as.numeric(as.character(year)) - if (!year %in% seq(2020, 2100, 5)) { - stop(sprintf( - "Error: The specified year '%s' is invalid. Accepted years are: %s. Please rerun the `calc_ResidEm_grp` function with a valid year.", - year, paste(seq(2020, 2100, 5), collapse = ", ") - )) - } - if (!pollutant %in% rhap::all_pollutants) { - stop(sprintf( - "Error: The specified pollutant '%s' is invalid. Accepted pollutants are: %s. Please rerun the `calc_ResidEm_grp` function with a valid pollutant.", - pollutant, paste(rhap::all_pollutants, collapse = ", ") - )) - } - if (!region %in% rhap::gcam_regions) { - stop(sprintf( - "Error: The specified region '%s' is invalid. Accepted regions are: %s. Please rerun the `calc_ResidEm_grp` function with a valid region.", - region, paste(rhap::gcam_regions, collapse = ", ") - )) - } - # Then, load the rgcam project if prj not passed as a parameter: if (!is.null(db_path) & !is.null(db_name)) { @@ -79,11 +59,39 @@ calc_ResidEm_grp <- function(db_path = NULL, query_path = "./inst/extdata", db_n final_db_year, max(rgcam::getQuery(prj, "nonCO2 emissions by sector (excluding resource production)")$year) ) + # Find last historical year + years_in_prj <- listYears(prj) + years_in_prj <- years_in_prj[!is.na(years_in_prj)] + years_in_prj <- setdiff(years_in_prj, NA) + base_year <- dplyr::if_else(2021 %in% years_in_prj, 2021, 2015) + base_year_p <- dplyr::if_else(base_year == 2021, base_year, 2020) + + # More checks + year <- as.numeric(as.character(year)) + if (!year %in% c(base_year_p, seq(2025, 2100, 5))) { + stop(sprintf( + "Error: The specified year '%s' is invalid. Accepted years are: %s. Please rerun the `calc_ResidEm_grp` function with a valid year.", + year, paste(c(base_year_p, seq(2025, 2100, 5)), collapse = ", ") + )) + } + if (!pollutant %in% rhap::all_pollutants) { + stop(sprintf( + "Error: The specified pollutant '%s' is invalid. Accepted pollutants are: %s. Please rerun the `calc_ResidEm_grp` function with a valid pollutant.", + pollutant, paste(rhap::all_pollutants, collapse = ", ") + )) + } + if (!region %in% rhap::gcam_regions) { + stop(sprintf( + "Error: The specified region '%s' is invalid. Accepted regions are: %s. Please rerun the `calc_ResidEm_grp` function with a valid region.", + region, paste(rhap::gcam_regions, collapse = ", ") + )) + } + if (year > final_db_year) { stop(sprintf( "Error: The specified year '%s' is invalid. The database only contains data up to %s. Accepted years are: %s. Please rerun the `calc_ResidEm_grp` function with a valid year.", - year, final_db_year, paste(seq(2020, final_db_year, 5), collapse = ", ") + year, final_db_year, paste(c(base_year_p, seq(2025, final_db_year, 5)), collapse = ", ") )) } diff --git a/R/calc_hap_impacts.R b/R/calc_hap_impacts.R index 787e09c..3fa83f9 100644 --- a/R/calc_hap_impacts.R +++ b/R/calc_hap_impacts.R @@ -1,3 +1,5 @@ + + ##' calc_hap_impacts #' #' @@ -15,6 +17,8 @@ #' @param map Produce the maps. By default = FALSE #' @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 countries Target countries over which the model will be fit (defaults to \code{"All"}). +#' For a complete list of valid country names, check \code{unique(rhap::panel_data$country_name)}. #' @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 #' @importFrom magrittr %>% @@ -22,7 +26,7 @@ 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", + saveOutput = TRUE, map = FALSE, anim = TRUE, HIA_var = "deaths", countries = 'All', normalized = FALSE, by_gr = FALSE) { Country <- country <- sector <- scenario <- region <- year <- group <- ghg <- Units <- value <- adj <- value_reg <- dec_share <- Pollutant <- `ISO 3` <- @@ -35,6 +39,54 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n VOC_per_100k <- `GCAM Region` <- Percentage <- adj_country <- country_map <- gdp_pc_dol2011_ppp_gr <- gdp_pc_dol2011_ppp <- mapCountries <- . <- NULL + + output <- calc_hap_impacts_preprocess(db_path, query_path, db_name, prj_name, + scen_name, queries, final_db_year, + HIA_var, countries) + + output_fin <- calc_hap_impacts_predict(data = output, saveOutput, map, anim, + HIA_var, countries, normalized, by_gr) + + return(invisible(output_fin)) + +} + + + +##' calc_hap_impacts_preprocess +#' +#' +#' Prepare the data from GCAM to latter estimate the health impacts attributable to household air pollution for GCAM scenarios. +#' Run `calc_hap_impacts_predict` to compute the health estimates +#' @keywords health impacts, HAP +#' @return Health impacts attributable to HAP for all the selected years +#' @param db_path Path to the GCAM database +#' @param query_path Path to the query file +#' @param db_name Name of the GCAM database +#' @param prj_name Name of the rgcam project. This can be an existing project, or, if not, this will be the name +#' @param scen_name Vector names of the GCAM scenarios to be processed +#' @param queries Name of the GCAM query file. The file by default includes the queries required to run rfasst +#' @param final_db_year Final year in the GCAM database (this allows to process databases with user-defined "stop periods") +#' @param HIA_var Health metric to be predicted. c("deaths", "yll", "dalys"). By default = deaths +#' @param countries Target countries over which the model will be fit (defaults to \code{"All"}). +#' For a complete list of valid country names, check \code{unique(rhap::panel_data$country_name)}. +#' @importFrom magrittr %>% +#' @export + +calc_hap_impacts_preprocess <- function(db_path = NULL, query_path = "./inst/extdata", db_name = NULL, prj_name, + scen_name, queries = "queries_rhap.xml", final_db_year = 2100, + HIA_var = "deaths", countries = 'All') { + 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 <- + unit <- flsp_m2 <- Model <- Scenario <- Region <- Variable <- Unit <- gdp <- + value_agg <- share_pop <- gdp_pc <- gdp_agg <- share_gdp <- gdp_dol2011_ppp <- + gdp_dol2011_ppp_gr <- data_name <- country_name <- bias.adder <- pred_value <- + 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 <- + gdp_pc_dol2011_ppp_gr <- gdp_pc_dol2011_ppp <- . <- NULL + # Check user input if (!HIA_var %in% c("deaths", "yll", "dalys")) { stop(sprintf( @@ -49,62 +101,40 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n final_db_year, paste(seq(2015, 2100, 5), collapse = ", ") )) } - if (!is.logical(saveOutput)) { - stop(sprintf( - "Error: The specified saveOutput '%s' is invalid. Accepted saveOutput values are: TRUE, FALSE. Please rerun the calc_hap_impacts function with a valid saveOutput value.", - saveOutput - )) - } - if (!is.logical(map)) { - stop(sprintf( - "Error: The specified map '%s' is invalid. Accepted map values are: TRUE, FALSE. Please rerun the calc_hap_impacts function with a valid map value.", - map - )) - } - if (!is.logical(anim)) { - stop(sprintf( - "Error: The specified anim '%s' is invalid. Accepted anim values are: TRUE, FALSE. Please rerun the calc_hap_impacts function with a valid anim value.", - anim - )) - } - if (!is.logical(normalized)) { - stop(sprintf( - "Error: The specified normalized '%s' is invalid. Accepted normalized values are: TRUE, FALSE. Please rerun the calc_hap_impacts function with a valid normalized value.", - normalized - )) - } - if (!is.logical(by_gr)) { - stop(sprintf( - "Error: The specified by_gr '%s' is invalid. Accepted by_gr values are: TRUE, FALSE. Please rerun the calc_hap_impacts function with a valid by_gr value.", - by_gr - )) + if ((length(countries) == 1 && countries != "All") || (length(countries) != 1 && !any(countries == "All"))) { + wrong_countries <- setdiff(countries, unique(rhap::panel_data$country_name)) + if (length(wrong_countries) == 1) { + stop(sprintf( + "Error: The specified country '%s' is invalid. Run `unique(rhap::panel_data$country_name)` to see the accepted country names are. Please rerun the `fit_model` function with valid `coutries`.", + wrong_countries)) + } else if (length(wrong_countries) > 1) { + stop(sprintf( + "Error: The specified countries '%s' are invalid. Run `unique(rhap::panel_data$country_name)` to see the accepted country names are. Please rerun the `fit_model` function with valid `coutries`.", + paste(wrong_countries, collapse = ", "))) + } } - # Ancillary Functions `%!in%` <- Negate(`%in%`) # Add converter: Teragram to kt CONV_Tg_kt <- 1E3 - # Create the directory if they do not exist: - if (!dir.exists("output")) dir.create("output") - # Then, load the rgcam project if prj not passed as a parameter: if (!is.null(db_path) & !is.null(db_name)) { rlang::inform("Creating project ...") conn <- rgcam::localDBConn(db_path, - db_name, - migabble = FALSE + db_name, + migabble = FALSE ) prj <- suppressWarnings( rgcam::addScenario(conn, - prj_name, - scen_name, - paste0(query_path, "/", queries), - saveProj = TRUE + prj_name, + scen_name, + paste0(query_path, "/", queries), + saveProj = TRUE ) ) } else { @@ -117,6 +147,13 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n max(rgcam::getQuery(prj, "nonCO2 emissions by sector (excluding resource production)")$year) ) + # Find last historical year + years_in_prj <- listYears(prj) + years_in_prj <- years_in_prj[!is.na(years_in_prj)] + years_in_prj <- setdiff(years_in_prj, NA) + base_year <- dplyr::if_else(2021 %in% years_in_prj, 2021, 2015) + base_year_p <- dplyr::if_else(base_year == 2021, base_year, 2020) + #----- # EXTRACT DATA FROM GCAM SCENARIO OUTPUTS @@ -127,7 +164,11 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n reg_to_ctry <- rhap::Percen %>% dplyr::select(region = `GCAM Region`, country = Country) %>% dplyr::distinct() %>% - dplyr::mutate(country = stringr::str_to_title(country)) + dplyr::mutate(country = stringr::str_to_title(country), + # update regions if necessary + region = dplyr::if_else(base_year == 2021 & region == 'Europe_Eastern', 'Ukraine', region) + ) + # 1- Emissions em <- rgcam::getQuery(prj, "nonCO2 emissions by sector (excluding resource production)") %>% @@ -182,7 +223,10 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n tibble::as_tibble() %>% dplyr::mutate( Pollutant = dplyr::if_else(Pollutant == "POM", "OC", Pollutant), - year = as.numeric(as.character(year)) + year = as.numeric(as.character(year)), + # update base year and regions if necessary + year = dplyr::if_else(base_year == 2021 & year == 2020, 2021, year), + `GCAM Region` = dplyr::if_else(base_year == 2021 & `GCAM Region` == 'Europe_Eastern', 'Ukraine', `GCAM Region`) ) %>% dplyr::rename( region = `GCAM Region`, @@ -333,12 +377,25 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n dplyr::rename( ssp = scenario, country = region + ) %>% + dplyr::mutate( + year = as.numeric(year), + # update base year if necessary + year = dplyr::if_else(base_year == 2021 & year == 2020, 2021, year), + year = as.character(year) ) + # Add missing countries gdp_adj <- rhap::ssp_gdp_adj %>% tidyr::gather(year, value, -Model, -Scenario, -Region, -Variable, -Unit) %>% - dplyr::mutate(year = gsub("X", "", year)) %>% + dplyr::mutate( + year = gsub("X", "", year), + # update base year if necessary + year = as.numeric(year), + year = dplyr::if_else(base_year == 2021 & year == 2020, 2021, year), + year = as.character(year) + ) %>% dplyr::filter( year %in% unique(gdp_ctry$year), Scenario != "Historical Reference" @@ -360,6 +417,12 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n dplyr::rename( ssp = scenario, country = region + ) %>% + dplyr::mutate( + year = as.numeric(year), + # update base year if necessary + year = dplyr::if_else(base_year == 2021 & year == 2020, 2021, year), + year = as.character(year) ) # Process Population: Population is evenly distributed across groups, but could be updated @@ -454,7 +517,12 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n ssp = dplyr::if_else(grepl("SSP5", scenario), "SSP5", ssp) ) %>% dplyr::select(scenario, ssp, country, year, group, share_gdp) %>% - dplyr::mutate(year = as.character(year)) + dplyr::mutate( + year = as.numeric(year), + # update base year if necessary + year = dplyr::if_else(base_year == 2021 & year == 2020, 2021, year), + year = as.character(year) + ) # Adjust Taiwan using China's shares gdp_share_twn <- gdp_share %>% @@ -613,22 +681,178 @@ 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% unique(fit_model(HIA_var = HIA_var, countries = countries)[[2]])) %>% dplyr::mutate(across(where(is.numeric), ~ ifelse(is.finite(.), ., NA_real_))) %>% tidyr::drop_na() + # 3 - Return output + return(invisible(output)) +} + + + +##' calc_hap_impacts_predict +#' +#' +#' Run the model to estimate the health impacts attributable to household air pollution for GCAM scenarios. +#' Notice that the required dataset can be obtained through GCAM or uploaded manually through the `data` parameter +#' @keywords health impacts, HAP +#' @return Health impacts attributable to HAP for all the selected years +#' @param data Dataset with ghg, gdppc and flsp estimates to compute health impacts. Structure defined +#' as the output of `calc_hap_impacts_preprocess` +#' @param saveOutput Writes the emission files. By default=TRUE +#' @param map Produce the maps. By default = FALSE +#' @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 countries Target countries over which the model will be fit (defaults to \code{"All"}). +#' For a complete list of valid country names, check \code{unique(rhap::panel_data$country_name)}. +#' @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 +#' @importFrom magrittr %>% +#' @export + +calc_hap_impacts_predict <- function(data, saveOutput = TRUE, map = FALSE, anim = TRUE, + HIA_var = "deaths", countries = 'All', + normalized = FALSE, by_gr = FALSE) { + 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 <- + unit <- flsp_m2 <- Model <- Scenario <- Region <- Variable <- Unit <- gdp <- + value_agg <- share_pop <- gdp_pc <- gdp_agg <- share_gdp <- gdp_dol2011_ppp <- + gdp_dol2011_ppp_gr <- data_name <- country_name <- bias.adder <- pred_value <- + 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 + + # Check user input + if (!HIA_var %in% c("deaths", "yll", "dalys")) { + stop(sprintf( + "Error: The specified HIA_var '%s' is invalid. Accepted HIA_var are: %s. Please rerun the calc_hap_impacts function with a valid HIA_var value.", + HIA_var, paste(c("deaths", "yll", "dalys"), collapse = ", ") + )) + } + if (!is.logical(saveOutput)) { + stop(sprintf( + "Error: The specified saveOutput '%s' is invalid. Accepted saveOutput values are: TRUE, FALSE. Please rerun the calc_hap_impacts function with a valid saveOutput value.", + saveOutput + )) + } + if (!is.logical(map)) { + stop(sprintf( + "Error: The specified map '%s' is invalid. Accepted map values are: TRUE, FALSE. Please rerun the calc_hap_impacts function with a valid map value.", + map + )) + } + if (!is.logical(anim)) { + stop(sprintf( + "Error: The specified anim '%s' is invalid. Accepted anim values are: TRUE, FALSE. Please rerun the calc_hap_impacts function with a valid anim value.", + anim + )) + } + if (!is.logical(normalized)) { + stop(sprintf( + "Error: The specified normalized '%s' is invalid. Accepted normalized values are: TRUE, FALSE. Please rerun the calc_hap_impacts function with a valid normalized value.", + normalized + )) + } + if (!is.logical(by_gr)) { + stop(sprintf( + "Error: The specified by_gr '%s' is invalid. Accepted by_gr values are: TRUE, FALSE. Please rerun the calc_hap_impacts function with a valid by_gr value.", + by_gr + )) + } + if(!identical(colnames(data), c('scenario','country_name','year','log_PrimPM25_per_100k','log_NOx_per_100k','log_VOC_per_100k','log_gdppc_ppp_dol2011','log_flsp','pop'))){ + stop(sprintf( + "Error: The specified data is invalid. Run `calc_hap_impacts_preproces` to obtain a valid dataset. Please rerun the calc_hap_impacts_predict function with a valid data input." + )) + } + if ((length(countries) == 1 && countries != "All") || (length(countries) != 1 && !any(countries == "All"))) { + wrong_countries <- setdiff(countries, unique(rhap::panel_data$country_name)) + if (length(wrong_countries) == 1) { + stop(sprintf( + "Error: The specified country '%s' is invalid. Run `unique(rhap::panel_data$country_name)` to see the accepted country names are. Please rerun the `fit_model` function with valid `coutries`.", + wrong_countries)) + } else if (length(wrong_countries) > 1) { + stop(sprintf( + "Error: The specified countries '%s' are invalid. Run `unique(rhap::panel_data$country_name)` to see the accepted country names are. Please rerun the `fit_model` function with valid `coutries`.", + paste(wrong_countries, collapse = ", "))) + } + } + + + + # Ancillary Functions + `%!in%` <- Negate(`%in%`) + + # Create the directory if they do not exist: + if (countries == 'All') { + output_dir <- 'output' + } else { + list_iso <- rhap::panel_data %>% + dplyr::filter(country_name %in% countries) %>% + dplyr::pull(iso) %>% + unique() + output_dir <- paste0('output_', paste(list_iso, collapse = '_')) + } + if (!dir.exists(output_dir)) dir.create(output_dir) + + #----- # PREDICTION + output <- data # Get model to subtract coefficients and predict - model.fixed <- fit_model(HIA_var = HIA_var)[[1]] + model.fixed <- fit_model(HIA_var = HIA_var, countries = countries)[[1]] + data_train <- fit_model(HIA_var = HIA_var, countries = countries)[[3]] + + # Restrict data to selected countries + if (countries != 'All') output <- output %>% dplyr::filter(country_name %in% countries) # Transform data to panel and predict output.panel <- plm::pdata.frame(output, index = c("country_name", "year")) - # Predict - output$pred_value <- stats::predict(model.fixed, output.panel) + # Predict - manual implementation of stats::predict(model.fixed, output.panel) + predict_plm_within <- function(model, train_data, newdata) { + # train_data = the original plain data frame used to fit (data) + # newdata = plain data frame to predict on (output) + + coefs <- coef(model) + rhs_vars <- names(coefs) + + # 1. Compute Xβ on TRAINING data + X_train <- as.matrix(train_data[, rhs_vars]) + xb_train <- as.numeric(X_train %*% coefs) + + # 2. Recover fixed effect per country: + # FE_i = mean(y_i) - mean(Xβ_i) [the within estimator identity] + train_data$.xb <- xb_train + fe <- tapply( + train_data$log_Deaths_per_100k - train_data$.xb, + train_data$country_name, + mean + ) + + # 3. Compute Xβ on newdata + missing_vars <- setdiff(rhs_vars, names(newdata)) + if (length(missing_vars) > 0) + stop("Missing in newdata: ", paste(missing_vars, collapse = ", ")) + + X_new <- as.matrix(newdata[, rhs_vars]) + xb_new <- as.numeric(X_new %*% coefs) + + # 4. Look up country FE for each row of newdata + country_fe <- fe[as.character(newdata$country_name)] + + if (any(is.na(country_fe))) { + warning("Unseen countries — using mean fixed effect as fallback") + country_fe[is.na(country_fe)] <- mean(fe) + } + + xb_new + as.numeric(country_fe) + } + output$pred_value <- predict_plm_within(model.fixed, data_train, output) output_fin <- output %>% dplyr::mutate(pred_var = paste0("pred_log_", HIA_var, "_per_100K")) %>% @@ -662,8 +886,8 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n # Create a function to write the data (by scenario) output.write <- function(df) { df <- as.data.frame(df) - utils::write.csv(df, paste0("output/", unique(df$scenario), "_HAP_", unique(HIA_var), ".csv"), - row.names = FALSE, fileEncoding = "UTF-8" + utils::write.csv(df, paste0(output_dir, "/", unique(df$scenario), "_HAP_", unique(HIA_var), ".csv"), + row.names = FALSE, fileEncoding = "UTF-8" ) } @@ -675,7 +899,9 @@ 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 (saveOutput == TRUE & by_gr == TRUE) { # Create the directory if they do not exist: - if (!dir.exists("output/by_gr")) dir.create("output/by_gr") + if (!dir.exists(file.path(output_dir, "by_gr"))) dir.create(file.path(output_dir, "by_gr")) + # Restrict data to selected countries + if (countries != 'All') output_gr <- output_gr %>% dplyr::filter(country_name %in% countries) output.panel.gr <- plm::pdata.frame(output_gr, index = c("country_name", "year")) @@ -710,15 +936,15 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n 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" + utils::write.csv(paste0(file.path(output_dir, "by_gr"), "/HAP_", unique(HIA_var), "_byGR", ".csv"), + row.names = FALSE, fileEncoding = "UTF-8" ) } # Add map if (map == TRUE) { # Create the directory if they do not exist: - if (!dir.exists("output/maps")) dir.create("output/maps") + if (!dir.exists(file.path(output_dir, "maps"))) dir.create(file.path(output_dir, "maps")) # The variable to be plotted depends on if the user selects or no to use normalized values if (normalized == TRUE) { @@ -739,7 +965,7 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n value = var_to_plot ) %>% tidyr::complete(tidyr::nesting(scenario, year, pred_var), - subRegion = unique(rmap::mapCountries$region) + subRegion = unique(rmap::mapCountries$region) ) mapCountries <<- rmap::mapCountries @@ -752,7 +978,7 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n # 1. Plot rmap::map( data = output_fin_map %>% dplyr::filter(scenario == sc), - folder = paste("output/maps/map", sc, "allYears", sep = "_"), + folder = paste(file.path(output_dir, "maps", "map"), sc, "allYears", sep = "_"), legendType = "pretty", background = TRUE, animate = anim, @@ -764,32 +990,32 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n # 2. Reorder folders and rename figures # 2.1. move the allYears figure file.rename( - from = file.path(paste("output/maps/map", sc, "allYears", sep = "_"), "map_param_PRETTY.png"), - to = paste0("output/maps/map_", sc, "_allYears", ".png") + from = file.path(paste(file.path(output_dir, "maps", "map"), sc, "allYears", sep = "_"), "map_param_PRETTY.png"), + to = paste0(file.path(output_dir, "maps", "map"), "_", sc, "_allYears", ".png") ) # 2.2. move all annual figures - files_to_move <- list.files(file.path(paste("output/maps/map", sc, "allYears", sep = "_"), "byYear"), full.names = TRUE) + files_to_move <- list.files(file.path(paste(file.path(output_dir, "maps", "map"), sc, "allYears", sep = "_"), "byYear"), full.names = TRUE) success <- sapply(files_to_move, function(file) { - file.rename(file, file.path(paste("output/maps/map", sc, "allYears", sep = "_"), basename(file))) + file.rename(file, file.path(paste(file.path(output_dir, "maps", "map"), sc, "allYears", sep = "_"), basename(file))) }) # 2.3. remove unnecessary directories and files if (all(success)) { - unlink(file.path(paste("output/maps/map", sc, "allYears", sep = "_"), "byYear"), recursive = TRUE) + unlink(file.path(paste(file.path(output_dir, "maps", "map"), sc, "allYears", sep = "_"), "byYear"), recursive = TRUE) } else { message("Some files could not be moved. The source folder was not deleted.") } - unlink(file.path(paste("output/maps/map", sc, "allYears", sep = "_"), "map_param_MEAN_PRETTY.png"), recursive = TRUE) - unlink(file.path(paste("output/maps/map", sc, "allYears", sep = "_"), "map_param.csv"), recursive = TRUE) + unlink(file.path(paste(file.path(output_dir, "maps", "map"), sc, "allYears", sep = "_"), "map_param_MEAN_PRETTY.png"), recursive = TRUE) + unlink(file.path(paste(file.path(output_dir, "maps", "map"), sc, "allYears", sep = "_"), "map_param.csv"), recursive = TRUE) # 2.4. rename folder - if (dir.exists(file.path(paste("output/maps/map", sc, "byYear", sep = "_")))) { - unlink(file.path(paste("output/maps/map", sc, "byYear", sep = "_")), recursive = TRUE) + if (dir.exists(file.path(paste(file.path(output_dir, "maps", "map"), sc, "byYear", sep = "_")))) { + unlink(file.path(paste(file.path(output_dir, "maps", "map"), sc, "byYear", sep = "_")), recursive = TRUE) } file.rename( - file.path(paste("output/maps/map", sc, "allYears", sep = "_")), - file.path(paste("output/maps/map", sc, "byYear", sep = "_")) + file.path(paste(file.path(output_dir, "maps", "map"), sc, "allYears", sep = "_")), + file.path(paste(file.path(output_dir, "maps", "map"), sc, "byYear", sep = "_")) ) } @@ -800,7 +1026,7 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n rmap::map( data = output_fin_map %>% dplyr::filter(year == y) %>% dplyr::rename(class = scenario), - folder = paste("output/maps/map", "allScen", y, sep = "_"), + folder = paste(file.path(output_dir, "maps", "map"), "allScen", y, sep = "_"), legendType = "pretty", background = TRUE, animate = anim, @@ -812,16 +1038,16 @@ calc_hap_impacts <- function(db_path = NULL, query_path = "./inst/extdata", db_n # 2. Reorder folders and rename figures # 2.1. remove an intermediate folder file.rename( - from = file.path(paste("output/maps/map", "allScen", y, sep = "_"), "map_param_PRETTY.png"), - to = paste0("output/maps/map_", "allScen_", y, ".png") + from = file.path(paste(file.path(output_dir, "maps", "map"), "allScen", y, sep = "_"), "map_param_PRETTY.png"), + to = paste0(file.path(output_dir, "maps"), "/map_", "allScen_", y, ".png") ) - unlink(paste("output/maps/map", "allScen", y, sep = "_"), recursive = TRUE) + unlink(paste(file.path(output_dir, "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) - if (!dir.exists("output/maps/map_allScen_byYear")) dir.create("output/maps/map_allScen_byYear") + files_to_move <- list.files(path = file.path(file.path(output_dir, "maps")), pattern = "^map_allScen_", full.names = TRUE) + if (!dir.exists(file.path(output_dir, "maps", "map_allScen_byYear"))) dir.create(paste0(file.path(output_dir, "maps", "map_allScen_byYear"))) success <- sapply(files_to_move, function(file) { - file.rename(file, file.path("output/maps/map_allScen_byYear", basename(file))) + file.rename(file, file.path(file.path(output_dir, "maps", "map_allScen_byYear"), basename(file))) }) } diff --git a/R/fit_model.R b/R/fit_model.R index c0b3e97..e0c2f39 100644 --- a/R/fit_model.R +++ b/R/fit_model.R @@ -8,11 +8,13 @@ #' @source Details on plm estimation can be found here: https://cran.r-project.org/web/packages/plm/plm.pdf#' #' @keywords Econometric model; fixed effects #' @param HIA_var Health metric to be predicted. c("deaths", "yll", "dalys") +#' @param countries Target countries over which the model will be fit (defaults to \code{"All"}). +#' For a complete list of valid country names, check \code{unique(rhap::panel_data$country_name)}. #' @importFrom magrittr %>% #' @export #' @return Regression model and coefficients for prediction -fit_model <- function(HIA_var) { +fit_model <- function(HIA_var, countries = 'All') { iso <- country_name <- year <- pop <- continent <- dev <- log_AAP <- value <- Model <- Scenario <- Region <- Variable <- Unit <- . <- NULL @@ -23,6 +25,18 @@ fit_model <- function(HIA_var) { HIA_var, paste(c("deaths", "yll", "dalys"), collapse = ", ") )) } + if ((length(countries) == 1 && countries != "All") || (length(countries) != 1)) { + wrong_countries <- as.character(setdiff(countries, unique(rhap::panel_data$country_name))) + if (length(wrong_countries) == 1) { + stop(sprintf( + "Error: The specified country '%s' is invalid. Run `unique(rhap::panel_data$country_name)` to see the accepted country names are. Please rerun the `fit_model` function with valid `coutries`.", + wrong_countries)) + } else if (length(wrong_countries) > 1) { + stop(sprintf( + "Error: The specified countries '%s' are invalid. Run `unique(rhap::panel_data$country_name)` to see the accepted country names are. Please rerun the `fit_model` function with valid `coutries`.", + paste(wrong_countries, collapse = ", "))) + } + } # Adjust the data data <- rhap::panel_data %>% @@ -30,7 +44,6 @@ fit_model <- function(HIA_var) { dplyr::mutate(year = as.character(year)) %>% dplyr::select(-log_AAP, -log_HDD_value, -log_CDD_value) %>% dplyr::filter(stats::complete.cases(.)) - predictable_regions <- unique(data$country_name) # dplyr::select the dependent variable (deaths, YLLs, or DALYs) @@ -53,12 +66,28 @@ fit_model <- function(HIA_var) { "log_gdppc_ppp_dol2011 + log_flsp" )) - model.fixed <- plm::plm( - model_formula, - data = data, - index = c("country_name", "year"), - model = "within" - ) + if (countries == 'All') { + model.fixed <- plm::plm( + model_formula, + data = data, + index = c("country_name", "year"), + model = "within" + ) + + predictable_regions <- unique(data$country_name) + data_train <- data + + } else { + data_ctry_full <- data %>% dplyr::filter(country_name %in% countries) + model.fixed <- plm::plm( + model_formula, + data = data_ctry_full, + index = c("country_name", "year"), + model = "within") + + predictable_regions <- unique(data_ctry_full$country_name) + data_train <- data_ctry_full + } - return(list(model.fixed, predictable_regions)) + return(list(model.fixed, predictable_regions, data_train)) } diff --git a/R/utils.R b/R/utils.R new file mode 100644 index 0000000..c312dd5 --- /dev/null +++ b/R/utils.R @@ -0,0 +1,58 @@ +#' listYears +#' +#' Return the years of the queries available for a scenario in a project data set. +#' This function requires the data set to have been previously loaded, so it cannot take a file name. +#' Source: gcamreport +#' +#' @param projData The data set to report on. +#' @param scenarios The name(s) of the scenario(s) to report on. If NULL, report on all of them. +#' @param queries The name(s) of the queries(s) to report on. If NULL, report on all of them. +#' @param anyscen If TRUE, then list queries that are in any scenario. If FALSE, list queries that are in all scenarios. +#' @return list of years reported in the project/scenario/queries. +#' @export +listYears <- function (projData, scenarios = NULL, queries = NULL, anyscen = TRUE) { + if (is.character(projData)) { + stop("listYears() function requires the data set to have been already loaded.") + } + if (is.null(scenarios)) { + scenarios <- rgcam::listScenarios(projData) + } + if (is.null(queries)) { + queries <- rgcam::listQueries(projData) + } + sqlist <- lapply(scenarios, function(scen) { + lapply(queries, function(quer) { + if ("year" %in% names(projData[[scen]][[quer]])) { + yy = unique(projData[[scen]][[quer]][['year']]) + if (length(yy) > 100) { + NULL + } else { + yy + } + } else { + NULL + } + }) + }) + + combine <- if (anyscen) union else intersect + + if (identical(combine, union)) { + # Union case: count appearances and keep values appearing >10 times + # (avoid problems with 2020 and 2021) + all_years <- unlist(sqlist) + all_years <- all_years[!is.na(all_years)] + year_counts <- table(all_years) + if (length(queries) == 1) { + result <- sort(as.numeric(names(year_counts))) + } else { + result <- sort(as.numeric(names(year_counts[year_counts > (length(queries)/2 + 1)]))) + } + } else { + # Intersect case: just intersect all elements + result <- Reduce(intersect, Reduce(intersect, sqlist)) + } + + return(result) +} + diff --git a/man/calc_hap_impacts.Rd b/man/calc_hap_impacts.Rd index 95391ee..27918e3 100644 --- a/man/calc_hap_impacts.Rd +++ b/man/calc_hap_impacts.Rd @@ -16,6 +16,7 @@ calc_hap_impacts( map = FALSE, anim = TRUE, HIA_var = "deaths", + countries = "All", normalized = FALSE, by_gr = FALSE ) @@ -43,6 +44,9 @@ calc_hap_impacts( \item{HIA_var}{Health metric to be predicted. c("deaths", "yll", "dalys"). By default = deaths} +\item{countries}{Target countries over which the model will be fit (defaults to \code{"All"}). +For a complete list of valid country names, check \code{unique(rhap::panel_data$country_name)}.} + \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} diff --git a/man/calc_hap_impacts_predict.Rd b/man/calc_hap_impacts_predict.Rd new file mode 100644 index 0000000..98c79fa --- /dev/null +++ b/man/calc_hap_impacts_predict.Rd @@ -0,0 +1,46 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/calc_hap_impacts.R +\name{calc_hap_impacts_predict} +\alias{calc_hap_impacts_predict} +\title{calc_hap_impacts_predict} +\usage{ +calc_hap_impacts_predict( + data, + saveOutput = TRUE, + map = FALSE, + anim = TRUE, + HIA_var = "deaths", + countries = "All", + normalized = FALSE, + by_gr = FALSE +) +} +\arguments{ +\item{data}{Dataset with ghg, gdppc and flsp estimates to compute health impacts. Structure defined +as the output of `calc_hap_impacts_preprocess`} + +\item{saveOutput}{Writes the emission files. By default=TRUE} + +\item{map}{Produce the maps. By default = FALSE} + +\item{anim}{If set to TRUE, produces multi-year animations. By default=TRUE} + +\item{HIA_var}{Health metric to be predicted. c("deaths", "yll", "dalys"). By default = deaths} + +\item{countries}{Target countries over which the model will be fit (defaults to \code{"All"}). +For a complete list of valid country names, check \code{unique(rhap::panel_data$country_name)}.} + +\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} +} +\value{ +Health impacts attributable to HAP for all the selected years +} +\description{ +Run the model to estimate the health impacts attributable to household air pollution for GCAM scenarios. +Notice that the required dataset can be obtained through GCAM or uploaded manually through the `data` parameter +} +\keyword{HAP} +\keyword{health} +\keyword{impacts,} diff --git a/man/calc_hap_impacts_preprocess.Rd b/man/calc_hap_impacts_preprocess.Rd new file mode 100644 index 0000000..091ba03 --- /dev/null +++ b/man/calc_hap_impacts_preprocess.Rd @@ -0,0 +1,48 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/calc_hap_impacts.R +\name{calc_hap_impacts_preprocess} +\alias{calc_hap_impacts_preprocess} +\title{calc_hap_impacts_preprocess} +\usage{ +calc_hap_impacts_preprocess( + db_path = NULL, + query_path = "./inst/extdata", + db_name = NULL, + prj_name, + scen_name, + queries = "queries_rhap.xml", + final_db_year = 2100, + HIA_var = "deaths", + countries = "All" +) +} +\arguments{ +\item{db_path}{Path to the GCAM database} + +\item{query_path}{Path to the query file} + +\item{db_name}{Name of the GCAM database} + +\item{prj_name}{Name of the rgcam project. This can be an existing project, or, if not, this will be the name} + +\item{scen_name}{Vector names of the GCAM scenarios to be processed} + +\item{queries}{Name of the GCAM query file. The file by default includes the queries required to run rfasst} + +\item{final_db_year}{Final year in the GCAM database (this allows to process databases with user-defined "stop periods")} + +\item{HIA_var}{Health metric to be predicted. c("deaths", "yll", "dalys"). By default = deaths} + +\item{countries}{Target countries over which the model will be fit (defaults to \code{"All"}). +For a complete list of valid country names, check \code{unique(rhap::panel_data$country_name)}.} +} +\value{ +Health impacts attributable to HAP for all the selected years +} +\description{ +Prepare the data from GCAM to latter estimate the health impacts attributable to household air pollution for GCAM scenarios. +Run `calc_hap_impacts_predict` to compute the health estimates +} +\keyword{HAP} +\keyword{health} +\keyword{impacts,} diff --git a/man/fit_model.Rd b/man/fit_model.Rd index df936a2..a02efcb 100644 --- a/man/fit_model.Rd +++ b/man/fit_model.Rd @@ -7,10 +7,13 @@ Details on plm estimation can be found here: https://cran.r-project.org/web/packages/plm/plm.pdf#' } \usage{ -fit_model(HIA_var) +fit_model(HIA_var, countries = "All") } \arguments{ \item{HIA_var}{Health metric to be predicted. c("deaths", "yll", "dalys")} + +\item{countries}{Target countries over which the model will be fit (defaults to \code{"All"}). +For a complete list of valid country names, check \code{unique(rhap::panel_data$country_name)}.} } \value{ Regression model and coefficients for prediction diff --git a/man/listYears.Rd b/man/listYears.Rd new file mode 100644 index 0000000..e22cee5 --- /dev/null +++ b/man/listYears.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils.R +\name{listYears} +\alias{listYears} +\title{listYears} +\usage{ +listYears(projData, scenarios = NULL, queries = NULL, anyscen = TRUE) +} +\arguments{ +\item{projData}{The data set to report on.} + +\item{scenarios}{The name(s) of the scenario(s) to report on. If NULL, report on all of them.} + +\item{queries}{The name(s) of the queries(s) to report on. If NULL, report on all of them.} + +\item{anyscen}{If TRUE, then list queries that are in any scenario. If FALSE, list queries that are in all scenarios.} +} +\value{ +list of years reported in the project/scenario/queries. +} +\description{ +Return the years of the queries available for a scenario in a project data set. +This function requires the data set to have been previously loaded, so it cannot take a file name. +Source: gcamreport +} diff --git a/tests/testthat/test1-run.R b/tests/testthat/test1-run.R index 10de434..b618eba 100644 --- a/tests/testthat/test1-run.R +++ b/tests/testthat/test1-run.R @@ -117,4 +117,26 @@ 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." ) + expect_error( + calc_hap_impacts( + prj_name = prj_name, scen_name = scen_name, + final_db_year = 2050, HIA_var = "deaths", + countries = "dummy", + saveOutput = TRUE, map = FALSE, anim = FALSE, + normalized = TRUE, by_gr = TRUE + ), + "Error: The specified country 'dummy' is invalid. Run `unique(rhap::panel_data$country_name)` to see the accepted country names are. Please rerun the `fit_model` function with valid `coutries`.", + fixed = T + ) + expect_error( + calc_hap_impacts( + prj_name = prj_name, scen_name = scen_name, + final_db_year = 2050, HIA_var = "deaths", + countries = c("dummy1", "dummy2", "Portugal"), + saveOutput = TRUE, map = FALSE, anim = FALSE, + normalized = TRUE, by_gr = TRUE + ), + "Error: The specified countries 'dummy1, dummy2' are invalid. Run `unique(rhap::panel_data$country_name)` to see the accepted country names are. Please rerun the `fit_model` function with valid `coutries`.", + fixed = T + ) }) diff --git a/tests/testthat/test3-other.R b/tests/testthat/test3-other.R index acb402b..17a3730 100644 --- a/tests/testthat/test3-other.R +++ b/tests/testthat/test3-other.R @@ -41,9 +41,37 @@ test_that("fit_model", { testResult <- get(load(file.path(rprojroot::find_root(rprojroot::is_testthat), "testOutputs/fit_model.RData"))) testthat::expect_equal(testOutput, testResult) + + testOutput <- fit_model(HIA_var = "deaths", countries = "Spain") + testResult <- get(load(file.path(rprojroot::find_root(rprojroot::is_testthat), "testOutputs/fit_model_esp.RData"))) + testthat::expect_equal(testOutput, testResult) + # error messages expect_error( fit_model(HIA_var = "Death"), "Error: The specified HIA_var 'Death' is invalid. Accepted HIA_var are: deaths, yll, dalys. Please rerun the `fit_model` function with a valid HIA_var." ) + + expect_error( + fit_model(HIA_var = "deaths", countries = "dummy"), + "Error: The specified country 'dummy' is invalid. Run `unique(rhap::panel_data$country_name)` to see the accepted country names are. Please rerun the `fit_model` function with valid `coutries`.", + fixed = TRUE + ) + + expect_error( + fit_model(HIA_var = "deaths", countries = c("dummy","Zambia","Uzbekistan2")), + "Error: The specified countries 'dummy, Uzbekistan2' are invalid. Run `unique(rhap::panel_data$country_name)` to see the accepted country names are. Please rerun the `fit_model` function with valid `coutries`.", + fixed = TRUE + ) + }) + + +test_that("listYears", { + prj_name <- file.path(rprojroot::find_root(rprojroot::is_testthat), "testInputs", "test_prj_v7p1.dat") + prj <- rgcam::loadProject(prj_name) + + prj_years <- listYears(prj) + testthat::expect_equal(prj_years, c(1975, 1990, 2005, 2010, 2015, 2020, 2025, 2030, 2035, 2040, 2045, 2050, 2055)) +}) + diff --git a/tests/testthat/testOutputs/fit_model_esp.RData b/tests/testthat/testOutputs/fit_model_esp.RData new file mode 100644 index 0000000..82b4d1a Binary files /dev/null and b/tests/testthat/testOutputs/fit_model_esp.RData differ