epsP and epsH included as arguments to recalib(); tests updated

statistikat · Jan 7, 2020 · 3def4390b8e87f92ee726242ed4e1057cad88ae3 · 3def439 · lintr-bot · Jan 7, 2020
1 parent 3fee6ed
commit 3def4390b8e87f92ee726242ed4e1057cad88ae3
diff --git a/R/ipf.r b/R/ipf.r
@@ -267,10 +267,10 @@ calibH <- function(i, dat, error, valueH, hColNames, bound, verbose, calIter,
           !is.na(fVariableForCalibrationIPF) &
             (abs(1 / fVariableForCalibrationIPF - 1) > epsHcur),
           list(maxFac = max(abs(1 / fVariableForCalibrationIPF - 1)), .N,
-               head(epsHcur, 1),
+               epsH = head(epsHcur, 1),
                sumCalibWeight = sum(get(variableKeepingTheCalibWeight) *
                                       representativeHouseholdForCalibration),
-               head(value, 1)),
+               PopMargin = head(value, 1)),
           by = eval(hColNames[[i]])]
         print(tmp[order(maxFac, decreasing = TRUE), ])
 
@@ -549,7 +549,7 @@ ipf <- function(
 
   OriginalSortingVariable <- V1 <- epsvalue <-
     f <- temporary_hvar <-
-    value <- wValue <- representativeHouseholdForCalibration <- NULL
+    value <- wValue <- representativeHouseholdForCalibration <- ..hid <- NULL
   dat_original <- dat
   dat <- copy(dat)
   ## originalsorting is fucked up without this

diff --git a/R/recalib.R b/R/recalib.R
@@ -46,6 +46,8 @@
 #' @param conH.var character vector containig household-specific variables to
 #'   which weights should be calibrated or a list of such character vectors. Contingency tables for the population
 #'   are calculated per `period` using `weights`.
+#' @param epsP numeric value specifying the convergence limit for `conP.var` or `conP`, see [ipf()].
+#' @param epsH numeric value specifying the convergence limit for `conH.var` or `conH`, see [ipf()].
 #' @param ... additional arguments passed on to function [ipf()] from this
 #'   package.
 #'
@@ -71,19 +73,20 @@
 #'                            strata = "region", period = "year")
 #'
 #' # calibrate weight for bootstrap replicates
-#' dat_boot_calib <- recalib(dat_boot, conP.var = "gender", conH.var = "region")
+#' dat_boot_calib <- recalib(dat_boot, conP.var = "gender", conH.var = "region",
+#'                           verbose = TRUE)
 #'
 #'
 #' # calibrate on other variables
 #' dat_boot_calib <- recalib(dat_boot, conP.var = c("gender", "age"),
-#'                           conH.var = c("region", "hsize"))
+#'                           conH.var = c("region", "hsize"), verbose = TRUE)
 #' 
 #' # supply contingency tables directly                           
 #' conP <- xtabs(pWeight ~ age + gender + year, data = eusilc)
 #' conH <- xtabs(pWeight ~ hsize + region + year, data = eusilc[!duplicated(paste(db030,year))])
 #' dat_boot_calib <- recalib(dat_boot, conP.var = NULL,
 #'                           conH.var = NULL, conP = list(conP),
-#'                           conH = list(conH))                           
+#'                           conH = list(conH), verbose = TRUE)                           
 #' }
 #'
 #' @export recalib
@@ -92,10 +95,11 @@
 recalib <- function(
   dat, hid = attr(dat, "hid"), weights = attr(dat, "weights"), b.rep =
     attr(dat, "b.rep"), period = attr(dat, "period"), conP.var = NULL,
-  conH.var = NULL, ...) {
+  conH.var = NULL, epsP = 1e-2, epsH = 2e-2, ...) {
 
-  verbose <- epsP <- epsH <- bound <- maxIter <- meanHH <- check_hh_vars <-
-    hidfactor <- calibWeight <- FirstPersonInHousehold_ <- NULL
+  hidfactor <- calibWeight <- FirstPersonInHousehold_ <- verbose <-
+    bound <- maxiter <- meanHH <- check_hh_vars <- allPthenH <-
+    returnNA <- conversion_messages <- maxIter <- NULL
 
   ##########################################################
   # INPUT CHECKING
@@ -160,33 +164,27 @@ recalib <- function(
   ##########################################################
 
   # define default values for ipf
+  ipfDefaults <- formals(ipf)
+  ipfDefaults <- ipfDefaults[!names(ipfDefaults)%in%names(formals(recalib))]
   ellipsis <- list(...)
+  # set these to FALSE by default
+  ellipsis[["check_hh_vars"]] <- getEllipsis2("check_hh_vars", 
+                                             FALSE, ellipsis)
+  ellipsis[["conversion_messages"]] <- getEllipsis2("conversion_messages", 
+                                                   FALSE, ellipsis)
   ellipsis[["verbose"]] <- getEllipsis2("verbose", TRUE, ellipsis)
-  ellipsis[["epsP"]] <- getEllipsis2("epsP", 1e-2, ellipsis)
-  ellipsis[["epsH"]] <- getEllipsis2("epsH", 5e-2, ellipsis)
-  ellipsis[["conP"]] <- getEllipsis2("conP", NULL, ellipsis)
-  ellipsis[["conH"]] <- getEllipsis2("conH", NULL, ellipsis)
-  ellipsis[["bound"]] <- getEllipsis2("bound", 4, ellipsis)
-  ellipsis[["maxIter"]] <- getEllipsis2("maxIter", 100, ellipsis)
-  ellipsis[["meanHH"]] <- getEllipsis2("meanHH", TRUE, ellipsis)
-  ellipsis[["check_hh_vars"]] <- getEllipsis2("check_hh_vars", FALSE, ellipsis)
-  ellipsis[["conversion_messages"]] <- getEllipsis2("conversion_messages", FALSE,
-                                                   ellipsis)
-
+  ellipsis <- lapply(names(ipfDefaults),function(z){
+    getEllipsis2(z,ipfDefaults[[z]],ellipsis)
+  })
+  names(ellipsis) <- names(ipfDefaults)
+
   ellipsisNames <- names(ellipsis)
   ellipsisContent <- paste0("ellipsis[['",ellipsisNames,"']]")
   eval(parse(text = paste(
     ellipsisNames,
     ellipsisContent, sep = "<-"
   )))
-
-  if(length(ellipsis[["conH"]])==0){
-    conH <- NULL
-  }
-  if(length(ellipsis[["conP"]])==0){
-    conP <- NULL
-  }
-
+
   # check conP and conH
   conPnames <- lapply(conP,function(z){
     z <- names(dimnames(z))
@@ -228,7 +226,7 @@ recalib <- function(
   # recode household and person variables to factor
   # improves runtime for ipf
   #
-  vars <- c(period, unique(unlist(c(conP.var,conH.var))))
+  vars <- c(period, unique(unlist(c(conP.var,conH.var,conPnames,conHnames))))
   vars.class <- unlist(lapply(dat[, mget(vars)], function(z) {
     z.class <- class(z)
     if (z.class[1] == "labelled"){
@@ -246,26 +244,23 @@ recalib <- function(
 
   # calculate contingency tables
   for(p in seq_along(conP.var)){
-    existTab <- sapply(conPnames,setequal,y=conP.var[p])
+    existTab <- sapply(conPnames,setequal,y=conP.var[[p]])
     if(all(!existTab)){
-      conP <- c(conP,
-                makeCalibTable(dat,weights=weights,period=period,
-                               vars=conP.var[[p]]))
+      formTab <- paste(weights,"~",paste(c(period,conP.var[[p]]),collapse="+"))
+      conP <- c(conP,list(xtabs(formTab,data=dat)))
     }else{
       stop("contingency table for ",paste(conP.var[p],collapse = ", ")," was supplied through paramter conP AND conP.var")
     } 
   }
 
   dat[,FirstPersonInHousehold_:=c(1L,rep(0,.N-1)),by=c(hid,period)]
   for(h in seq_along(conH.var)){
-    existTab <- sapply(conHnames,setequal,y=conH.var[h])
+    existTab <- sapply(conHnames,setequal,y=conH.var[[h]])
     if(all(!existTab)){
-      conH <- c(conH,
-                makeCalibTable(dat[FirstPersonInHousehold_==1],
-                               weights=weights,period=period,
-                               vars=conH.var[[h]]))
+      formTab <- paste(weights,"~",paste(c(period,conH.var[[h]]),collapse="+"))
+      conH <- c(conH,list(xtabs(formTab,data=dat[FirstPersonInHousehold_==1])))
     }else{
-      stop("contingency table for ",paste(conH.var[h],collapse = ", ")," was supplied through paramter conH AND conH.var")
+      stop("contingency table for ",paste(conH.var[[h]],collapse = ", ")," was supplied through paramter conH AND conH.var")
     } 
   }
 
@@ -276,7 +271,7 @@ recalib <- function(
   # calibrate weights to conP and conH
   select.var <- unique(c("hidfactor", weights, period, unlist(c(conP.var,conH.var,conPnames,conHnames))))
   calib.fail <- c()
-
+ 
   for (g in b.rep) {
     set(dat, j = g, value = dt.eval("dat[,", g, "*", weights, "]"))
 
@@ -292,7 +287,7 @@ recalib <- function(
       nrow(check.z) > 0
     })
 
-    if (!is.null(conP.var) | !is.null(conH.var)) {
+    if (!is.null(conP) | !is.null(conH)) {
       if (any(unlist(c(check.conH, check.conP)))) {
         calib.fail <- c(calib.fail, g)
         set(dat, j = g, value = NA_real_)
@@ -301,7 +296,9 @@ recalib <- function(
           dat = copy(dat[, mget(c(g, select.var))]), conP = conP,
           conH = conH, verbose = verbose, epsP = epsP, epsH = epsH,
           w = g, bound = bound, maxIter = maxIter, meanHH = meanHH,
-          hid = "hidfactor", check_hh_vars = check_hh_vars
+          hid = "hidfactor", check_hh_vars = check_hh_vars,
+          allPthenH = allPthenH, returnNA = returnNA,
+          conversion_messages = conversion_messages
         )[, calibWeight])
         if (dat[, any(is.na(get(g)))]) {
           calib.fail <- c(calib.fail, g)

diff --git a/man/recalib.Rd b/man/recalib.Rd
diff --git a/src/symbols.rds b/src/symbols.rds
diff --git a/tests/testthat/test_calcstError.R b/tests/testthat/test_calcstError.R
@@ -13,7 +13,8 @@ eusilc <- eusilc[!db040 %in% c("Vienna", "Lower Austria", "Upper Austria")]
 
 eusilc <- draw.bootstrap(eusilc, REP = 2, hid = "db030", weights = "db090",
                          period = "year", strata = "db040")
-eusilc <- recalib(eusilc, conP.var = "rb090", conH.var = "db040")
+eusilc <- recalib(eusilc, conP.var = c("rb090", "age"),
+                  conH.var = c("db040", "hsize"))
 
 # test input parameter
 test_that("test para - data", {

diff --git a/tests/testthat/test_calcstError_long.R b/tests/testthat/test_calcstError_long.R
@@ -15,7 +15,8 @@ if (Sys.getenv("SURVEYSD_ADDITIONAL_TEST") == "TRUE") {
                            period = "year", strata = "db040")
   eusilc <- recalib(
     eusilc, hid = "db030", weights = "db090", b.rep = paste0("w", 1:2),
-    period = "year", conP.var = "rb090", conH.var = "db040")
+    period = "year", conP.var = c("rb090", "age"),
+    conH.var = c("db040", "hsize"))
 
   # these are some longer tests
 

diff --git a/tests/testthat/test_plot.R b/tests/testthat/test_plot.R
@@ -13,7 +13,8 @@ eusilc <- eusilc[!db040 %in% c("Vienna", "Lower Austria", "Upper Austria")]
 
 eusilc <- draw.bootstrap(eusilc, REP = 2, hid = "db030", weights = "db090",
                          period = "year", strata = "db040")
-eusilc <- recalib(eusilc, conP.var = "rb090", conH.var = "db040")
+eusilc <- recalib(eusilc, conP.var = c("rb090", "age"),
+                  conH.var = c("db040", "hsize"))
 res <- calc.stError(eusilc, var = "povmd60", group =
                       list("rb090", "db040", c("rb090", "db040"), "hsize"))
 

diff --git a/tests/testthat/test_recalib.R b/tests/testthat/test_recalib.R
@@ -7,7 +7,7 @@ library(surveysd)
 library(laeken)
 library(data.table)
 
-eusilc <- surveysd:::demo.eusilc()
+eusilc <- surveysd:::demo.eusilc(n=4)
 eusilc <- draw.bootstrap(eusilc, REP = 2, hid = "db030", weights = "db090",
                          period = "year", strata = "db040")
 
@@ -82,7 +82,7 @@ test_that("test para - conP conH", {
   conP2 <- xtabs(db090 ~ rb090 + year, data = eusilc)
   # conP2 <- xtabs(db090 ~ rb090 + db040 + year, data = eusilc)
 
-  conH1 <- xtabs(db090 ~ hsize + db040 + year, data = eusilc[!duplicated(paste(db030,year))])
+  conH1 <- xtabs(db090 ~ hsize + year, data = eusilc[!duplicated(paste(db030,year))])
   conH2 <- xtabs(db090 ~ db040 + year, data = eusilc[!duplicated(paste(db030,year))])
 
   expect_error(
@@ -103,7 +103,7 @@ test_that("test para - conP conH", {
 test_that("test return", {
   dat.calib <- recalib(
     eusilc, hid = "db030", weights = "db090", b.rep = paste0("w", 1:2), period =
-      "year", conP.var = "rb090", conH.var = "db040")
+      "year", conP.var = c("rb090", "age"), conH.var = c("db040", "hsize"))
   rb090.compare <- dat.calib[, lapply(.SD, sum), by = list(year, rb090),
                              .SDcols = c("db090", paste0("w", 1:2))]
   expect_true(all(!is.na(rb090.compare[, .SD, .SDcols = paste0("w", 1:2)])))
@@ -113,7 +113,7 @@ test_that("test return", {
       abs(db090 - z) / db090
     }
     ), .SDcols = paste0("w", 1:2)]
-  expect_true(all(rb090.compare < 0.01))
+  expect_true(all(rb090.compare < formals(recalib)$epsP))
 
   db040.compare <- dat.calib[, lapply(
     .SD,
@@ -127,5 +127,5 @@ test_that("test return", {
       abs(db090 - z) / db090
     }
   ), .SDcols = paste0("w", 1:2)]
-  expect_true(all(db040.compare < 0.05))
+  expect_true(all(db040.compare < formals(recalib)$epsH))
 })
diff --git a/vignettes/error_estimation.Rmd b/vignettes/error_estimation.Rmd
@@ -34,7 +34,8 @@ set.seed(1234)
 eusilc <- demo.eusilc(prettyNames = TRUE)
 dat_boot <- draw.bootstrap(eusilc, REP = 10, hid = "hid", weights = "pWeight",
                            strata = "region", period = "year")
-dat_boot_calib <- recalib(dat_boot, conP.var = "gender", conH.var = "region")
+dat_boot_calib <- recalib(dat_boot, conP.var = "gender", conH.var = "region",
+                          verbose = TRUE)
 dat_boot_calib[, onePerson := nrow(.SD) == 1, by = .(year, hid)]
 
 ## print part of the dataset

diff --git a/vignettes/surveysd.Rmd b/vignettes/surveysd.Rmd
@@ -47,7 +47,8 @@ Calibrate each sample according to the distribution of `gender` (on a personal l
 (on a household level). 
 
 ```{r}
-dat_boot_calib <- recalib(dat_boot, conP.var = "gender", conH.var = "region")
+dat_boot_calib <- recalib(dat_boot, conP.var = "gender", conH.var = "region",
+                          verbose = TRUE)
 dat_boot_calib[1:5, .(year, povertyRisk, gender, pWeight, w1, w2, w3, w4)]
 ```