Drop dependency on assertthat and refactor calcTheil

Author: johanneskoch94

DESCRIPTION

@@ -46,7 +46,6 @@ Depends:
     R (>= 2.10.0)
 Imports:
     assertr,
-    assertthat,
     broom,
     car,
     countrycode,

NAMESPACE

@@ -199,7 +199,6 @@ importFrom(stats,na.exclude)
 importFrom(stats,na.omit)
 importFrom(stats,na.pass)
 importFrom(stats,nls)
-importFrom(stats,qnorm)
 importFrom(stats,sd)
 importFrom(stats,weighted.mean)
 importFrom(tibble,as_tibble)

R/calcTheil.R

@@ -1,101 +1,69 @@
 #' Calculate regional Theil-T index
 #'
-
-#' To calculate the regional Theil-T index (= correction to welfare function for
-#' a lognormal income distribution) we do the following:
-#' (1) convert country-level Gini coefficients to Theil (2) calculate contribution
-#' to Theil-T index that includes both between-countries and within-country inequality
-#' (see e.g. https://en.wikipedia.org/wiki/Theil_index). The latter can then be
-#' aggregated with calcOutput().
+#' To calculate the regional Theil-T index (= correction to welfare function for a lognormal income distribution) we do
+#' the following: (1) convert country-level Gini coefficients to Theil (2) calculate contribution to Theil-T index that
+#' includes both between-countries and within-country inequality (see e.g. https://en.wikipedia.org/wiki/Theil_index).
+#' The latter can then be aggregated with calcOutput().
 #'
 #' NB 1: the aggregation depends on the region mapping. It is implemented such
 #' that the regionmapping specified in getConfig()$regionmapping is used.
 #'
 #' NB 2: the result of calcOutput('Theil', aggregate = FALSE), is NOT the country
 #' Theil-T, but the unweighted contribution from a given country to the regional value.
 #'
-#' @return magpie objects of unweighted contribution to Theil,
-#' weights (= country shares of regional GDP), docstring
+#' @return magpie objects of unweighted contribution to Theil, weights (= country shares of regional GDP)
 #' @author Bjoern Soergel
 #' @seealso \code{\link{calcOutput}} \code{\link{convertGini},\link{readGini}}
 #' @examples
 #' \dontrun{
-#' a <- calcOutput("Theil")
+#'   calcOutput("Theil")
 #' }
 #'
-#' @importFrom stats qnorm
-
 calcTheil <- function() {
-
-  ##  helper functions.
-  TheilT.from.sigma <- function(sigma) {
-    # Theil T coefficient for lognormal distribution
-    TheilT <- sigma^2 / 2.
-    return(TheilT)
-  }
-
-  sigma.from.Gini <- function(G) {
-    # assuming lognormal distribution: convert Gini to sigmas
-    sigma <- sqrt(2) * qnorm((G + 1) / 2)
-    return(sigma)
-  }
-
-  # Gini and Theil
+  # Read Gini
   Gini <- readSource("Gini")
-  years <- getYears(Gini)
-  TheilT <- TheilT.from.sigma(sigma.from.Gini(Gini))
 
-  # population (in 1e6)
-  pop <- calcOutput(type = "Population", aggregate = FALSE)
-  sspnames <- c(paste0("SSP", 1:5), "SDP", "SSP2EU")
-  pop <- pop[, years, paste0("pop_", sspnames)]
-  getNames(pop) <- sspnames
-  getSets(pop) <- c("iso3c", "year", "scenario")
+  # Convert Gini to sigmas assuming lognormal distribution
+  sigma <- sqrt(2) * stats::qnorm((Gini + 1) / 2)
+  # Theil T coefficient for lognormal distribution
+  TheilT <- sigma^2 / 2
 
-  # GDP (in million $ PPP 2005)
-  GDPppp <- calcOutput(type = "GDP", aggregate = FALSE)
-  GDPnames <- paste0("gdp_", sspnames)
-  GDPppp <- GDPppp[, years, GDPnames]
-  getNames(GDPppp) <- sspnames
-  getSets(GDPppp) <- c("iso3c", "year", "scenario")
+  # We need the GDP and GDP per capita scenarios, for the scenarios and years of Gini.
+  # We set extension2150 = "constant" because the Gini coefficients are also extended in the same way.
+  # As a regionmapping we use the one set in the config (which is the default behavior). We call it explicitly as it is
+  # used in the calculations of the Theil contribution and weights.
+  s <- getNames(Gini)
+  y <- getYears(Gini)
+  gdp      <- calcOutput("GDP",   naming = "scenario", extension2150 = "constant", years = y, aggregate = FALSE)[, , s]
+  gdpReg   <- calcOutput("GDP",   naming = "scenario", extension2150 = "constant", years = y)[, , s]
+  gdppc    <- calcOutput("GDPpc", naming = "scenario", extension2150 = "constant", years = y, aggregate = FALSE)[, , s]
+  gdppcReg <- calcOutput("GDPpc", naming = "scenario", extension2150 = "constant", years = y)[, , s]
 
-  # allocate empty objects for storing Theil contribution and weights
-  contribTheilT <- pop
+  # Allocate empty objects for storing Theil contribution and weights
+  contribTheilT <- TheilT
   contribTheilT[, , ] <- NA
-  s_i <- pop
-  s_i[, , ] <- NA
-
-  # contribution to Theil index depends on region mapping. We always use the one specified in getConfig().
-  regionmapping <- read.csv(
-    toolGetMapping(
-      type = "regional", name = getConfig()$regionmapping,
-      returnPathOnly = TRUE, where = "mappingfolder"
-    ),
-    sep = ";", colClasses = "character"
-  )
+  weight <- TheilT
+  weight[, , ] <- NA
 
-  # GDP per capita
-  xbar_i <- GDPppp / pop
-  for (rr in unique(regionmapping$RegionCode)) {
+  # Compute Theil contribution and weights
+  regionmapping <- toolGetMapping(getConfig("regionmapping"), type = "regional", where = "mappingfolder")
+  for (rr in getRegions(gdppcReg)) {
     rrCountries <- regionmapping$CountryCode[regionmapping$RegionCode == rr]
-    # regional GDP per capita
-    GDPppp_rr <- dimSums(GDPppp[rrCountries, , ], dim = 1)
-    Ntot_rr <- dimSums(pop[rrCountries, , ], dim = 1)
-    xbar_rr <- GDPppp_rr / Ntot_rr
-    # contribution to Theil index (unweighted)
-    contribTheilT[rrCountries, , ] <- TheilT[rrCountries, , ] + log(xbar_i[rrCountries, , ] / xbar_rr)
-    # weights = income share of country i:
-    # s_i = N_i/N * xbar_i/xbar = GDP_i/GDP_rr # nolint
-    s_i[rrCountries, , ] <- GDPppp[rrCountries, , ] / GDPppp_rr
-    # sanity check: ensure that weights for a region sum to one (within floating point precision)
-    assertthat::assert_that(max(abs(dimSums(s_i[rrCountries, , ], dim = 1) - 1)) < 1e-10)
+    # Contribution to Theil index (unweighted)
+    contribTheilT[rrCountries, , ] <- TheilT[rrCountries, , ] + log(gdppc[rrCountries, , ] / gdppcReg[rr, , ])
+    # Weights = country shares of regional GDP
+    weight[rrCountries, , ] <- gdp[rrCountries, , ] / gdpReg[rr, , ]
+    # Sanity check: ensure that weights for a region sum to one (within floating point precision)
+    stopifnot(max(abs(dimSums(weight[rrCountries, , ], dim = 1) - 1)) < 1e-10)
   }
 
-  # for easier REMIND integration use same names as GDP scenarios for Theil
-  # change this if we later want to test effect of per capita income growth vs. inequality
-  getNames(contribTheilT) <- GDPnames
-  getNames(s_i) <- GDPnames
+  # For easier REMIND integration use same names as GDP scenarios for Theil
+  # Change this if we later want to test effect of per capita income growth vs. inequality
+  getNames(contribTheilT) <- paste0("gdp_", getNames(contribTheilT))
+  getNames(weight)        <- paste0("gdp_", getNames(weight))
 
-  return(list(x = contribTheilT, weight = s_i, unit = "-",
-              description = "aggregated: Theil-T index, not-aggregated: unweighted contribution to Theil-T"))
+  list(x = contribTheilT,
+       weight = weight,
+       unit = "-",
+       description = "Aggregated: Theil-T index. Not-aggregated: unweighted contribution to Theil-T")
 }