Merge pull request #57 from SoerenPannier/dev

Dev into master

.Rbuildignore

@@ -16,3 +16,4 @@ revdep
 
 CONTRIBUTING.md
 
+^CRAN-SUBMISSION$

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: emdi
 Title: Estimating and Mapping Disaggregated Indicators
-Version: 2.1.2
-Date: 2022-08-19
+Version: 2.2.0
+Date: 2023-06-13
 Authors@R: c(person("Sylvia", "Harmening", role="aut", email = "sylvia.harmening@fu-berlin.de"), 
     person("Ann-Kristin", "Kreutzmann", role="aut", email="ann-kristin.kreutzmann@fu-berlin.de"),
     person("Soeren", "Pannier", role=c("aut", "cre"), email="soeren.pannier@fu-berlin.de"),
@@ -28,13 +28,13 @@ Description: Functions that support estimating, assessing and mapping regional
     by Kreutzmann et al. (2019) <doi:10.18637/jss.v091.i07> and the second package vignette 
     "A Framework for Producing Small Area Estimates Based on Area-Level Models in R".
 Depends:
-    R (>= 3.5.0)
+    R (>= 4.2.0)
 License: GPL-2
 URL: https://github.com/SoerenPannier/emdi
 LazyData: true
 Encoding: UTF-8
 Copyright: inst/COPYRIGHTS
-RoxygenNote: 7.2.1
+RoxygenNote: 7.2.3
 Imports:
     nlme,
     moments,
@@ -53,13 +53,13 @@ Imports:
     readODS,
     formula.tools,
     saeRobust,
+    rlang, 
     spdep
 Suggests:
     testthat,
     R.rsp,
     simFrame,
     laeken,
-    rgeos,
-    maptools,
+    sf,
     graphics
 VignetteBuilder: R.rsp

NAMESPACE

@@ -143,7 +143,6 @@ importFrom(ggplot2,coord_fixed)
 importFrom(ggplot2,coord_flip)
 importFrom(ggplot2,element_blank)
 importFrom(ggplot2,facet_wrap)
-importFrom(ggplot2,fortify)
 importFrom(ggplot2,geom_abline)
 importFrom(ggplot2,geom_boxplot)
 importFrom(ggplot2,geom_density)
@@ -152,16 +151,15 @@ importFrom(ggplot2,geom_point)
 importFrom(ggplot2,geom_polygon)
 importFrom(ggplot2,geom_qq)
 importFrom(ggplot2,geom_segment)
+importFrom(ggplot2,geom_sf)
 importFrom(ggplot2,geom_smooth)
 importFrom(ggplot2,geom_text)
 importFrom(ggplot2,geom_vline)
 importFrom(ggplot2,ggplot)
 importFrom(ggplot2,ggtitle)
-importFrom(ggplot2,guides)
 importFrom(ggplot2,labs)
 importFrom(ggplot2,qplot)
 importFrom(ggplot2,scale_color_manual)
-importFrom(ggplot2,scale_colour_gradient)
 importFrom(ggplot2,scale_fill_gradient)
 importFrom(ggplot2,scale_fill_manual)
 importFrom(ggplot2,scale_linetype_discrete)
@@ -206,6 +204,7 @@ importFrom(parallelMap,parallelLibrary)
 importFrom(parallelMap,parallelStop)
 importFrom(readODS,write_ods)
 importFrom(reshape2,melt)
+importFrom(rlang,.data)
 importFrom(spdep,geary.test)
 importFrom(spdep,mat2listw)
 importFrom(spdep,moran.test)

NEWS.md

@@ -1,3 +1,14 @@
+# emdi 2.2.0
+* Extension of the ebp function to allow for population weights
+* Extension of the ebp function to allow the aggregation of the estimates to different levels  
+* a more flexible use of the custom_indicator agrument within the ebp function
+
+# emdi 2.1.3
+* Improved summary with clearer notation of R2
+* Updated Area Level (FH) vignette
+* Minor improvements in checking T/F in if clauses
+* Increased dependency to R (>= 4.2.0) corresponding to the imported package MuMin
+
 # emdi 2.1.2
 * Improved messages
 * Copyrights updated
@@ -44,21 +55,21 @@
 * Updated R version dependency
 
 # emdi 1.1.7
- 
+
 * Minor typos corrected
 * Unit-Tests adjusted to the forthcoming R-Version 4.0
 
 # emdi 1.1.6
- 
+
 * New and updated references.
 
 # emdi 1.1.5
- 
+
 * Tests updated to deal with new random number generation in R
 * Some spelling improved
 
 # emdi 1.1.4
- 
+
 * Fixed Bug in summary: R2 calculation with MuMIn is now fully working
 * Added feature: formula used in fixed is now preserved, even if passed to ebp as a variable
 
@@ -87,12 +98,12 @@ reproducibility even when the function is run in parallel mode.
 * Updated Vignette
 
 # emdi 1.1.0
-  
+
 * A new function `direct` is made available, which provides direct estimation for small areas.
 * The function `ebp` now allows for a user-defined threshold.
 * The function `ebp` is now able to perform a semi-parametric wild bootstrap for MSE estimation.
 * The function `ebp` has new default value for parallelization that automatically adopts for the operating system.
 * The two data sets `eusilcA_smp` and `eusilcA_pop` have been updated.
 * For the function `map_plot` additional customization is now applicable.
 * All methods for `emdi model` except plot are extended for `emdi direct`.
-* `subset` and `as.data.frame` have been added as methods for class `emdi.estimators`
+* `subset` and `as.data.frame` have been added as methods for class `emdi.estimators`

R/FH.R

@@ -354,7 +354,7 @@ fh <- function(fixed, vardir, combined_data, domains = NULL, method = "reml",
     B <- c(B, 0)
   }
 
-  if (!(method == "reblup" | method == "reblupbc")) {
+  if (!(method == "reblup" || method == "reblupbc")) {
     # Estimate sigma u ---------------------------------------------------------
     sigmau2 <- wrapper_estsigmau2(
       framework = framework, method = method,
@@ -385,7 +385,7 @@ fh <- function(fixed, vardir, combined_data, domains = NULL, method = "reml",
           combined_data = framework$combined_data
         )
       }
-      if ((method == "ml" | method == "reml") & correlation == "spatial") {
+      if ((method == "ml" || method == "reml") && correlation == "spatial") {
         # Spatial EBLUP --------------------------------------------------------
         eblup <- eblup_SFH(
           framework = framework, sigmau2 = sigmau2,
@@ -436,9 +436,9 @@ fh <- function(fixed, vardir, combined_data, domains = NULL, method = "reml",
         )
         MSE <- mse_data$mse_data
         MSE_method <- mse_data$MSE_method
-        if (mse_type == "spatialnonparboot" |
-          mse_type == "spatialnonparbootbc" |
-          mse_type == "spatialparboot" |
+        if (mse_type == "spatialnonparboot" ||
+          mse_type == "spatialnonparbootbc" ||
+          mse_type == "spatialparboot" ||
           mse_type == "spatialparbootbc") {
           successful_bootstraps <- mse_data$successful_bootstraps
         }
@@ -685,7 +685,7 @@ fh <- function(fixed, vardir, combined_data, domains = NULL, method = "reml",
         successful_bootstraps = NULL
       )
     }
-  } else if (method == "reblup" | method == "reblupbc") {
+  } else if (method == "reblup" || method == "reblupbc") {
 
     # Standard EBLUP -----------------------------------------------------------
     eblup <- eblup_robust(

R/back_transformation.R

@@ -198,7 +198,7 @@ arcsin_mse <- function(sigmau2 = sigmau2, combined_data = combined_data,
     )
     mse_backtransformed <- mse_backtransformed[[2]]$MSE
     mse_method <- "bootstrap"
-  } else if (mse_type == "jackknife" | mse_type == "weighted_jackknife") {
+  } else if (mse_type == "jackknife" || mse_type == "weighted_jackknife") {
     transformation <- "arcsin"
     jack_mse <- wrapper_MSE(
       framework = framework, combined_data = combined_data,

R/check_direct_arguments.R

@@ -164,15 +164,19 @@ direct_check <- function(y,
                      not be included in the indicator. For help, see Example 3
                      in help(direct)."))
       } else if (inherits(custom_indicator[[i]], "function") &&
-        !all(names(formals(custom_indicator[[i]])) ==
+        !all(names(formals(custom_indicator[[i]])) %in%
           c("y", "weights", "threshold"))) {
         stop(strwrap(prefix = " ", initial = "",
-                     "Functions for custom indicators need to have exactly the
-                     following three arguments: y, weights threshold; even
-                     though weights might not be needed and a threshold might
-                     not be included in the indicator. For help, see Example 3
-                     in help(direct)."))
-      }
+                     "Functions for custom indicators need to have the argument
+                     y and optional the argument weights and threshold. For
+                     help, see Example 3 in help(direct)."))
+      } else if (inherits(custom_indicator[[i]], "function") &&
+        !("y" %in% names(formals(custom_indicator[[i]])))) {
+        stop(strwrap(prefix = " ", initial = "",
+                     "Functions for custom indicators need to have the argument
+                     y and optional the argument weights and threshold. For
+                     help, see Example 3 in help(direct)."))
+     }
     }
   }
 }

R/check_ebp_arguments.R

@@ -47,7 +47,8 @@ ebp_check1 <- function(fixed, pop_data, pop_domains, smp_data, smp_domains, L) {
 }
 
 ebp_check2 <- function(threshold, transformation, interval, MSE, boot_type, B,
-                       custom_indicator, cpus, seed, na.rm, weights) {
+                       custom_indicator, cpus, seed, na.rm, weights,
+                       pop_weights) {
   if (!is.null(threshold) && !(is.numeric(threshold) &&
     length(threshold) == 1) && !inherits(threshold, "function")) {
     stop(strwrap(prefix = " ", initial = "",
@@ -124,21 +125,29 @@ ebp_check2 <- function(threshold, transformation, interval, MSE, boot_type, B,
 
     N_custom <- length(custom_indicator)
     for (i in seq_len(N_custom)) {
+      if(!is.null(pop_weights)) {
+        if(!all(c("y", "pop_weights") %in%
+                names(formals(custom_indicator[[i]])))) {
+          stop(strwrap(prefix = " ", initial = "",
+                     "Please provide the argument pop_weights to the your
+                     custom_indicator. All other indicators will be
+                     calculated using population weights."))
+        }
+      }
       if (!inherits(custom_indicator[[i]], "function")) {
         stop(strwrap(prefix = " ", initial = "",
                      "The elements of the list need to be functions. These
-                     functions for custom indicators need to have exactly the
-                     following two arguments: y, threshold; even though a
-                     threshold might not included in the indicator. For help
-                     see Example 2 in help(ebp)."))
+                     functions for custom indicators need to have the
+                     argument y and optional the agruments pop_weights and
+                     threshold. For help see Example 2 in help(ebp)."))
       } else if (inherits(custom_indicator[[i]], "function") &&
-        !all(names(formals(custom_indicator[[i]])) ==
-          c("y", "threshold"))) {
+        !all(names(formals(custom_indicator[[i]])) %in%
+          c("y", "pop_weights", "threshold"))) {
         stop(strwrap(prefix = " ", initial = "",
                      "Functions for custom indicators need to have exactly the
-                     following two arguments: y, threshold; even though a
-                     threshold might not included in the indicator. For help
-                     see Example 2 in help(ebp)."))
+                     following argument y and optional the arguments
+                     pop_weights, threshold. For help see Example 2 in
+                     help(ebp)."))
       }
     }
   }
@@ -165,12 +174,21 @@ ebp_check2 <- function(threshold, transformation, interval, MSE, boot_type, B,
                  "The weighted version of ebp is only available with the
                  ''parametric'' bootstrap."))
   }
+  if (is.character(pop_weights) && length(pop_weights) != 1 ||
+      !is.character(pop_weights) && !is.null(pop_weights)) {
+    stop(strwrap(prefix = " ", initial = "",
+                 "Pop_weights must be a vector of length 1 and of class
+                 character specifying the variable name of a numeric variable
+                 indicating weights in the population data. See also
+                 help(ebp)."))
+  }
 }
 
 
 # Functions called in notation
-fw_check1 <- function(pop_data, mod_vars, pop_domains, smp_data,
-                      fixed, smp_domains, threshold, weights) {
+fw_check1 <- function(pop_data, mod_vars, pop_domains, smp_data, fixed,
+                      smp_domains, aggregate_to, threshold, weights,
+                      pop_weights) {
   if (!all(mod_vars %in% colnames(pop_data))) {
     stop(strwrap(prefix = " ", initial = "",
                  paste0("Variable ",
@@ -235,12 +253,39 @@ fw_check1 <- function(pop_data, mod_vars, pop_domains, smp_data,
     }
   }
 
-  if (dim(pop_data)[1] < dim(smp_data)[1]) {
-    stop(strwrap(prefix = " ", initial = "",
+  if(is.null(aggregate_to) != TRUE){
+    if (!(aggregate_to %in% colnames(pop_data))) {
+      stop(paste0("The domain variable ", aggregate_to, " is not contained in
+                  pop_data. Please provide valid variable name for the
+                  aggregation."))
+    }
+  }
+
+  if (is.character(pop_weights)) {
+    if (!is.numeric(pop_data[[pop_weights]])) {
+      stop(strwrap(prefix = " ", initial = "",
+                   paste0("The variable ", pop_weights, " must be the name of a
+                          variable that is a numeric vector.")))
+    }
+  }
+  if (is.character(pop_weights)) {
+    if (!all(pop_data[[pop_weights]] >= 1)) {
+      stop(strwrap(prefix = " ", initial = "",
+                   paste0("Negative or zero weights are included in ",
+                          pop_weights, " Please remove obersvations with weight
+                          values smaller than 1.")))
+    }
+  }
+
+  if (is.null(pop_weights)) {
+    if (dim(pop_data)[1] < dim(smp_data)[1]) {
+      stop(strwrap(prefix = " ", initial = "",
                  "The population data set cannot have less observations than
                  the sample data set."))
+    }
   }
 
+
   if (inherits(threshold, "function") &&
     (!is.numeric(threshold(smp_data[[paste(fixed[2])]])) ||
       length(threshold(smp_data[[paste(fixed[2])]])) != 1)) {
@@ -254,7 +299,7 @@ fw_check1 <- function(pop_data, mod_vars, pop_domains, smp_data,
 
 
 fw_check2 <- function(pop_domains, pop_domains_vec, smp_domains,
-                      smp_domains_vec) {
+                      smp_domains_vec, aggregate_to, aggregate_to_vec) {
   if (!(is.numeric(pop_domains_vec) ||
     any(inherits(pop_domains_vec, "factor")))) {
     stop(strwrap(prefix = " ", initial = "",
@@ -267,6 +312,13 @@ fw_check2 <- function(pop_domains, pop_domains_vec, smp_domains,
                  paste0(smp_domains, " needs to be the name of a variable that
                         is numeric or a (ordered) factor.")))
   }
+  if(is.null(aggregate_to) != TRUE){
+    if (!(is.numeric(aggregate_to_vec) ||
+          any(inherits(aggregate_to_vec, "factor")))) {
+      stop(paste0(aggregate_to, " needs to be the name of a variable that is
+                  numeric or a (ordered) factor."))
+    }
+  }
   if ((is.numeric(pop_domains_vec) &&
     any(inherits(smp_domains_vec, "factor"))) ||
     (is.numeric(smp_domains_vec) &&

R/compare_plot.R

@@ -205,7 +205,7 @@ compare_plots <- function(object, type, selected_indicators, MSE, CV, label,
   area <- NULL
 
 
-  if (MSE == FALSE & CV == FALSE) {
+  if (MSE == FALSE && CV == FALSE) {
     plotList <- vector(mode = "list", length = length(selected_indicators) * 2)
     names(plotList) <- paste(rep(
       c("scatter", "line"),
@@ -214,7 +214,7 @@ compare_plots <- function(object, type, selected_indicators, MSE, CV, label,
     rep(selected_indicators, each = 2),
     sep = "_"
     )
-  } else if ((MSE == TRUE | CV == TRUE) & !(MSE == TRUE & CV == TRUE)) {
+  } else if ((MSE == TRUE || CV == TRUE) && !(MSE == TRUE && CV == TRUE)) {
     plotList <- vector(mode = "list", length = length(selected_indicators) * 4)
     names(plotList) <- paste(rep(
       c("scatter", "line"),
@@ -223,7 +223,7 @@ compare_plots <- function(object, type, selected_indicators, MSE, CV, label,
     rep(selected_indicators, each = 4),
     sep = "_"
     )
-  } else if (MSE == TRUE & CV == TRUE) {
+  } else if (MSE == TRUE && CV == TRUE) {
     plotList <- vector(mode = "list", length = length(selected_indicators) * 6)
     names(plotList) <- paste(rep(
       c("scatter", "line"),

R/compare_plot.ebp.R

@@ -22,12 +22,12 @@ compare_plot.ebp <- function(model = NULL, direct = NULL, indicator = "all",
                         )))
   }
 
-  if ((inherits(model, "ebp") & is.null(direct)) |
-    (inherits(direct, "ebp") & is.null(model))) {
+  if ((inherits(model, "ebp") && is.null(direct)) |
+    (inherits(direct, "ebp") && is.null(model))) {
     stop(strwrap(prefix = " ", initial = "",
                  paste0("If the model is of type 'ebp', the input argument
                         direct is required.")))
-  } else if (inherits(model, "ebp") & inherits(direct, "direct")) {
+  } else if (inherits(model, "ebp") && inherits(direct, "direct")) {
     compare_plot_ebp(
       model = model, direct = direct, indicator = indicator,
       MSE = MSE, CV = CV,