Running 'basic.R' [2s]
  Running 'check_functions.R' [170s]
  Running 'greg.R' [170s]
  Running 'huber2.R' [163s]
  Running 'ratio.R' [157s]
  Running 'reg.R' [155s]
Running the tests in 'tests/check_functions.R' failed.
Complete output:
  > # functions that are required in the testing suite
  > 
  > # utility function
  > all_equal <- function(target, current, label,
  +     tolerance = sqrt(.Machine$double.eps), scale = NULL,
  +     check.attributes = FALSE)
  + {
  +     if (missing(label))
  +         stop("Argument 'label' is missing\n")
  +     res <- all.equal(target, current, tolerance, scale,
  +         check.attributes = check.attributes)
  +     if (is.character(res))
  +         cat(paste0(label, ": ", res, "\n"))
  + }
  > 
  > # check functions for coef() and SE()
  > check <- function(ref, est, what, characteristic)
  + {
  +     all_equal(as.numeric(coef(ref)), as.numeric(coef(est)),
  +         paste0(what, ": ", characteristic, ": coef"))
  +     all_equal(as.numeric(SE(ref)), as.numeric(SE(est)),
  +         paste0(what,": ", characteristic, ": SE"))
  + }
  > 
  > proc.time()
     user  system elapsed 
     0.17    0.01    0.17 
Running the tests in 'tests/greg.R' failed.
Complete output:
  > # NOTE: The estimators of the population mean for a model without a
  > #       regression intercept are excluded from the test because they
  > #       differ from the implementation based on survey::calibrate
  > suppressPackageStartupMessages(library("survey"))
  > library("robsurvey", quietly = TRUE)
  > source("check_functions.R")
  > 
  > #===============================================================================
  > # 1 MU284strat data (weights are calibrated to N)
  > #===============================================================================
  > dataset <- "MU284strat"
  > dn <- svydesign(ids = ~LABEL, strata = ~Stratum, fpc = ~fpc,
  +          weights = ~weights, data = MU284strat)
  > 
  > f0 <- RMT85 ~ -1 + P85 + S82 + CS82
  > f <- RMT85 ~ P85 + S82 + CS82
  > 
  > N <- 284
  > aux_totals <- c(P85 = 8339, S82 = 13500, CS82 = 2583)
  > 
  > #--------------------------------------
  > # 1.1 total
  > #--------------------------------------
  > 
  > # GREG total with intercept
  > ref <- svytotal(~RMT85, calibrate(dn, f, c("(Intercept)" = N, aux_totals)))
  > est <- svytotal_reg(svyreg(f, dn), aux_totals, N, type = "ADU")
  > check(ref, est, dataset, "total with intercept")
  > # type "projective" and "ADU" coincide
  > est <- svytotal_reg(svyreg(f, dn), aux_totals, N, type = "projective")
  > check(ref, est, dataset, "total with intercept: projective")
  > 
  > # GREG total with intercept + heteroscedasticity
  > ref <- svytotal(~RMT85, calibrate(dn, f, c("(Intercept)" = N, aux_totals),
  +     variance = MU284strat$P85))
  > est <- svytotal_reg(svyreg(f, dn, var = ~ P85), aux_totals, N, type = "ADU")
  > check(ref, est, dataset, "total with intercept: hetero")
  > 
  > # GREG total without intercept
  > ref <- svytotal(~RMT85, calibrate(dn, f0, aux_totals))
  > est <- svytotal_reg(svyreg(f0, dn), aux_totals, type = "ADU")
  > check(ref, est, dataset, "total w/o intercept")
  > 
  > #--------------------------------------
  > # 1.2 mean
  > #--------------------------------------
  > # GREG mean with intercept
  > ref <- svymean(~RMT85, calibrate(dn, f, c("(Intercept)" = N, aux_totals)))
  > est <- svymean_reg(svyreg(f, dn), aux_totals, N, type = "ADU")
  > check(ref, est, dataset, "mean with intercept")
  > 
  > # GREG mean with intercept + heteroscedasticity
  > ref <- svymean(~RMT85, calibrate(dn, f, c("(Intercept)" = N, aux_totals),
  +     variance = MU284strat$P85))
  > est <- svymean_reg(svyreg(f, dn, var = ~ P85), aux_totals, N, type = "ADU")
  > check(ref, est, dataset, "mean with intercept: hetero")
  > 
  > # GREG mean without intercept, see NOTE:
  > if (FALSE) {
  +     ref <- svymean(~RMT85, calibrate(dn, f0, aux_totals))
  +     est <- svymean_reg(svyreg(f0, dn), aux_totals, N, type = "ADU")
  +     check(ref, est, dataset, "mean w/o intercept")
  + }
  > #===============================================================================
  > # 2 workplace data (weights are calibrated to N)
  > #===============================================================================
  > dataset <- "workplace"
  > dn <- svydesign(ids = ~ID, strata = ~strat, fpc = ~fpc, weights = ~weight,
  +     data = workplace)
  > 
  > f0 <- payroll ~ -1 + employment
  > f <- payroll ~ employment
  > 
  > N <- 86514
  > aux_totals <- c(employment = 953555 * 1.05)
  > 
  > #--------------------------------------
  > # 2.1 total
  > #--------------------------------------
  > # GREG total with intercept
  > ref <- svytotal(~payroll, calibrate(dn, f, c("(Intercept)" = N, aux_totals)))
  > est <- svytotal_reg(svyreg(f, dn), aux_totals, N, type = "ADU")
  > check(ref, est, dataset, "total with intercept")
  > # type "projective" and "ADU" coincide
  > est <- svytotal_reg(svyreg(f, dn), aux_totals, N, type = "projective")
  > check(ref, est, dataset, "total with intercept: projective")
  > 
  > # GREG total without intercept
  > ref <- svytotal(~payroll, calibrate(dn, f0, aux_totals))
  > est <- svytotal_reg(svyreg(f0, dn), aux_totals, type = "ADU")
  > check(ref, est, dataset, "total w/o intercept")
  > 
  > # GREG total without intercept + heteroscedasticity (ratio estimator)
  > ref <- svytotal(~payroll, calibrate(dn, f0, aux_totals,
  +     variance = workplace$employment))
  > est <- svytotal_reg(svyreg(f0, dn, var = ~employment), aux_totals,
  +     type = "ADU")
  > check(ref, est, dataset, "total w/o intercept: hetero")
  > 
  > #--------------------------------------
  > # 2.2 mean
  > #--------------------------------------
  > # GREG mean with intercept
  > ref <- svymean(~payroll, calibrate(dn, f, c("(Intercept)" = N, aux_totals)))
  > est <- svymean_reg(svyreg(f, dn), aux_totals, N, type = "ADU")
  > check(ref, est, dataset, "mean with intercept")
  > 
  > # GREG mean without intercept; see NOTE:
  > if (FALSE) {
  +     ref <- svymean(~payroll, calibrate(dn, f0, aux_totals))
  +     est <- svymean_reg(svyreg(f0, dn), aux_totals, N, type = "ADU")
  +     check(ref, est, dataset, "mean w/o intercept")
  + }
  > 
  > # GREG mean without intercept + heteroscedasticity (ratio estimator); see NOTE:
  > if (FALSE) {
  +     ref <- svymean(~payroll, calibrate(dn, f0, aux_totals,
  +         variance = workplace$employment))
  +     est <- svymean_reg(svyreg(f0, dn, var = ~employment), aux_totals, N,
  +         type = "ADU")
  +     check(ref, est, dataset, "mean w/o intercept: hetero")
  + }
  > 
  > #===============================================================================
  > # 3 workplace data (weights are NOT calibrated to N)
  > #===============================================================================
  > dataset <- "workplace"
  > # modify workplace data (mulitply weights of the first unit in each stratum
  > # by 0.25; thus, the weights are not calibrated)
  > wp <- workplace
  > wp[wp$ID == 2, "weight"] <- 0.25 * wp[wp$ID == 2, "weight"]
  > wp[wp$ID == 1, "weight"] <- 0.25 * wp[wp$ID == 1, "weight"]
  > wp[wp$ID == 7, "weight"] <- 0.25 * wp[wp$ID == 7, "weight"]
  > dn <- svydesign(ids = ~ID, strata = ~strat, fpc = ~fpc, weights = ~weight,
  +     data = wp)
  > 
  > #--------------------------------------
  > # 3.1 total
  > #--------------------------------------
  > # GREG total with intercept
  > ref <- svytotal(~payroll, calibrate(dn, f, c("(Intercept)" = N, aux_totals)))
  > est <- svytotal_reg(svyreg(f, dn), aux_totals, N, type = "ADU")
  > check(ref, est, paste0(dataset, " (not calibrated)"), "total with intercept")
  > # type "projective" and "ADU" coincide
  > est <- svytotal_reg(svyreg(f, dn), aux_totals, N, type = "projective")
  > check(ref, est, paste0(dataset, " (not calibrated)"),
  +     "total with intercept: projective")
  > 
  > # GREG total without intercept
  > ref <- svytotal(~payroll, calibrate(dn, f0, aux_totals))
  > est <- svytotal_reg(svyreg(f0, dn), aux_totals, type = "ADU")
  > check(ref, est, paste0(dataset, " (not calibrated)"), "total w/o intercept")
  > 
  > # GREG total without intercept + heteroscedasticity (ratio estimator)
  > ref <- svytotal(~payroll, calibrate(dn, f0, aux_totals,
  +     variance = wp$employment))
  > est <- svytotal_reg(svyreg(f0, dn, var = ~employment), aux_totals,
  +     type = "ADU")
  > check(ref, est, paste0(dataset, " (not calibrated)"),
  +     "total w/o intercept: hetero")
  > 
  > #--------------------------------------
  > # 3.2 mean
  > #--------------------------------------
  > # GREG mean with intercept
  > ref <- svymean(~payroll, calibrate(dn, f, c("(Intercept)" = N, aux_totals)))
  > est <- svymean_reg(svyreg(f, dn), aux_totals, N, type = "ADU")
  > check(ref, est, paste0(dataset, " (not calibrated)"), "mean with intercept")
  > 
  > # GREG mean without intercept; see NOTE:
  > if (FALSE) {
  +     ref <- svymean(~payroll, calibrate(dn, f0, aux_totals))
  +     est <- svymean_reg(svyreg(f0, dn), aux_totals, N, type = "ADU")
  +     check(ref, est, paste0(dataset, " (not calibrated)"), "mean w/o intercept")
  + }
  > 
  > # GREG mean without intercept + heteroscedasticity (ratio estimator); see NOTE:
  > if (FALSE) {
  +     ref <- svymean(~payroll, calibrate(dn, f0, aux_totals,
  +         variance = wp$employment))
  +     est <- svymean_reg(svyreg(f0, dn, var = ~employment), aux_totals, N,
  +         type = "ADU")
  +     check(ref, est, paste0(dataset, " (not calibrated)"),
  +         "mean w/o intercept: hetero")
  + }
  > 
  > proc.time()
     user  system elapsed 
     1.67    0.14    1.79 
Running the tests in 'tests/huber2.R' failed.
Complete output:
  > library("robsurvey", quietly = TRUE)
  > 
  > # test on floating point equality
  > all_equal <- function(target, current, label,
  +     tolerance = sqrt(.Machine$double.eps), scale = NULL,
  +     check.attributes = FALSE)
  + {
  +     if (missing(label))
  +         stop("Argument 'label' is missing\n")
  +     res <- all.equal(target, current, tolerance, scale,
  +         check.attributes = check.attributes)
  +     if (is.character(res))
  +         cat(paste0(label, ": ", res, "\n"))
  + }
  > 
  > # test against MASS::hubers
  > if (requireNamespace("MASS", quietly = TRUE)) {
  +     library(MASS)
  + 
  +     # convergence tolerance
  +     TOLERANCE <- 1e-8
  + 
  +     # make a copy of the function MASS::hubers
  +     hubers_mod <- hubers
  +     # replace the mad by the (scaled) IQR as initial scale estimator
  +     body(hubers_mod)[[7]][[3]][[2]] <- substitute(s0 <- IQR(y, type = 2) *
  +         0.7413)
  + 
  +     # payroll data
  +     attach(workplace)
  +     all_equal(huber2(payroll, rep(1, length(payroll)), tol = TOLERANCE),
  +         hubers_mod(payroll, tol = 1e-8)$mu,
  +         label = "huber2: payroll data: test against MASS::hubers failed")
  + 
  +     # x11 data
  +     x11 <- c(1:10, 1e99); w11 <- rep(1, 11)
  +     all_equal(huber2(x11, w11, tol = TOLERANCE),
  +         hubers_mod(x11, tol = 1e-8)$mu,
  +         label = "huber2: x11 data: test against MASS::hubers failed")
  + 
  +     # x11 data (scaled)
  +     x11 <- c(1:10, 1e99) / 1e20; w11 <- rep(1, 11)
  +     all_equal(huber2(x11, w11, tol = TOLERANCE),
  +         hubers_mod(x11, tol = 1e-8)$mu,
  +         label = "huber2: x11 data (scale): test against MASS::hubers failed")
  + }
  > 
  > proc.time()
     user  system elapsed 
     0.18    0.10    0.29 
Running the tests in 'tests/ratio.R' failed.
Complete output:
  > suppressPackageStartupMessages(library("survey"))
  > library("robsurvey", quietly = TRUE)
  > source("check_functions.R")
  > 
  > #===============================================================================
  > # workplace data
  > #===============================================================================
  > dataset <- "workplace"
  > dn <- svydesign(ids = ~ID, strata = ~strat, fpc = ~fpc, weights = ~weight,
  +     data = workplace)
  > TOT_employment <- 1001233
  > 
  > #--------------------------------------
  > # ratio estimator (Huber)
  > ref <- svyratio(~payroll, ~employment, dn)
  > est <- svyratio_huber(~payroll, ~employment, dn, k = Inf)
  > check(ref, est, dataset, "svyratio_huber")
  > 
  > # ratio estimator of the total (with "hajek" variance estimator)
  > ref_tot <- predict(ref, TOT_employment)
  > est_tot <- svytotal_ratio(est, TOT_employment, variance = "hajek")
  > all_equal(as.numeric(ref_tot$total), coef(est_tot),
  +     paste0(dataset, ": svyratio_huber: svytotal_ratio: coef"))
  > all_equal(as.numeric(ref_tot$se), SE(est_tot),
  +     paste0(dataset, ": svyratio_huber: svytotal_ratio: SE"))
  > 
  > # ratio estimator using svyreg_huberM (with "hajek" variance estimator)
  > est_reg <- svyreg_huberM(payroll ~ -1 + employment, dn, var = ~employment,
  +     k = Inf)
  > check(ref, est_reg, dataset, "svyreg_huberM")
  > 
  > # ratio estimator of the total using svyreg_huberM
  > est_reg_tot <- svytotal_reg(est_reg, TOT_employment, type = "projective")
  > est_rat_tot <- svytotal_ratio(est, TOT_employment, variance = "hajek")
  > check(est_reg_tot, est_rat_tot, dataset, "svyreg_huberM: svytotal_reg")
  > 
  > #--------------------------------------
  > # ratio estimator (Tukey)
  > ref <- svyratio(~payroll, ~employment, dn)
  > est <- svyratio_tukey(~payroll, ~employment, dn, k = Inf)
  > check(est, ref, dataset, "svyratio_tukey")
  > 
  > # ratio estimator of the total (with "hajek* variance estimator)
  > ref_tot <- predict(ref, TOT_employment)
  > est_tot <- svytotal_ratio(est, TOT_employment, variance = "hajek")
  > all_equal(as.numeric(ref_tot$total), coef(est_tot),
  +     paste0(dataset, ": svyratio_tukey: svytotal_ratio: coef"))
  > all_equal(as.numeric(ref_tot$se), SE(est_tot),
  +     paste0(dataset, ": svyratio_tukey: svytotal_ratio: SE"))
  > 
  > # ratio estimator using svyreg_tukeyM (with "hajek" variance estimator)
  > est_reg <- svyreg_tukeyM(payroll ~ -1 + employment, dn, var = ~employment,
  +     k = Inf)
  > check(est, est_reg, dataset, "svyreg_tukeyM")
  > 
  > # ratio estimator of the total using svyreg_tukeyM
  > est_reg_tot <- svytotal_reg(est_reg, TOT_employment, type = "projective")
  > est_rat_tot <- svytotal_ratio(est, TOT_employment, variance = "hajek")
  > check(est_reg_tot, est_rat_tot, dataset, "svyreg_tukeyM: svytotal_reg")
  > 
  > #===============================================================================
  > # MU284strat data
  > #===============================================================================
  > dataset <- "MU284strat"
  > dn <- svydesign(ids = ~LABEL, strata = ~Stratum, fpc = ~fpc,
  +          weights = ~weights, data = MU284strat)
  > TOT_P85 <- 8339
  > 
  > #--------------------------------------
  > # ratio estimator (Huber)
  > ref <- svyratio(~RMT85, ~P85, dn)
  > est <- svyratio_huber(~RMT85, ~P85, dn, k = Inf)
  > check(est, ref, dataset, "svyratio_huber")
  > 
  > # ratio estimator of the total (with "hajek" variance estimator)
  > ref_tot <- predict(ref, TOT_P85)
  > est_tot <- svytotal_ratio(est, TOT_P85, variance = "hajek")
  > all_equal(as.numeric(ref_tot$total), coef(est_tot),
  +     paste0(dataset, ": svyratio_huber: svytotal_ratio: coef"))
  > all_equal(as.numeric(ref_tot$se), SE(est_tot),
  +     paste0(dataset, ": svyratio_huber: svytotal_ratio: SE"))
  > 
  > # ratio estimator using svyreg_huberM (with "hajek" variance estimator)
  > est_reg <- svyreg_huberM(RMT85~ -1 + P85, dn, var = ~P85, k = Inf)
  > check(est, est_reg, dataset, "svyreg_huberM")
  > 
  > # ratio estimator of the total using svyreg_huberM
  > est_reg_tot <- svytotal_reg(est_reg, TOT_P85, type = "projective")
  > est_rat_tot <- svytotal_ratio(est, TOT_P85, variance = "hajek")
  > check(est_reg_tot, est_rat_tot, dataset, "svyreg_huberM: svytotal_reg")
  > 
  > #--------------------------------------
  > # ratio estimator (Tukey)
  > ref <- svyratio(~RMT85, ~P85, dn)
  > est <- svyratio_tukey(~RMT85, ~P85, dn, k = Inf)
  > check(est, ref, dataset, "svyratio_tukey")
  > 
  > # ratio estimator of the total (with "hajek" variance estimator)
  > ref_tot <- predict(ref, TOT_P85)
  > est_tot <- svytotal_ratio(est, TOT_P85, variance = "hajek")
  > all_equal(as.numeric(ref_tot$total), coef(est_tot),
  +     paste0(dataset, ": svyratio_tukey: svytotal_ratio: coef"))
  > all_equal(as.numeric(ref_tot$se), SE(est_tot),
  +     paste0(dataset, ": svyratio_tukey: svytotal_ratio: SE"))
  > 
  > # ratio estimator using svyreg_tukeyM (with "hajek" variance estimator)
  > est_reg <- svyreg_tukeyM(RMT85~ -1 + P85, dn, var = ~P85, k = Inf)
  > check(est, est_reg, dataset, "svyreg_tukeyM")
  > 
  > # ratio estimator of the total using svyreg_tukeyM
  > est_reg_tot <- svytotal_reg(est_reg, TOT_P85, type = "projective")
  > est_rat_tot <- svytotal_ratio(est, TOT_P85, variance = "hajek")
  > check(est_reg_tot, est_rat_tot, dataset, "svyreg_tukeyM: svytotal_reg")
  > 
  > proc.time()
     user  system elapsed 
     1.39    0.18    1.56 
Running the tests in 'tests/reg.R' failed.
Complete output:
  > suppressPackageStartupMessages(library("survey"))
  > library("robsurvey", quietly = TRUE)
  > attach(workplace)
  > 
  > # design object
  > dn <- svydesign(ids = ~ID, strata = ~strat, fpc = ~fpc, weights = ~weight,
  +     data = workplace)
  > 
  > source("check_functions.R")
  > 
  > #===============================================================================
  > # 1 Test against MASS::rlm
  > #===============================================================================
  > if (requireNamespace("MASS", quietly = TRUE)) {
  +     library(MASS)
  +     # make a copy of MASS::rlm
  +     rlm_mod <- MASS:::rlm.default
  +     # replace wmad-function with our weighted mad
  +     body(rlm_mod)[[4]] <- substitute(wmad <- function(x, w) weighted_mad(x, w))
  + 
  +     #---------------------------------------------------------------------------
  +     # Huber regression M-estimator
  +     est <- svyreg_huberM(payroll ~ employment, dn, k = 1.345, tol = 1e-9)
  +     ref <- rlm_mod(x = est$model$x, y = est$model$y, weights = est$model$w,
  +         k = 1.345, scale.est = "MAD", method = "M", wt.method = "case",
  +         acc = 1e-9, maxit = 50, test.vec = "coef")
  +     all_equal(coef(est), coef(ref), tolerance = 1e-7,
  +         label = "Huber regression regression: coefficients")
  +     # model-based covariance matrix
  +     all_equal(vcov(est, "model"), vcov(ref), tolerance = 1e-7,
  +         label = "Huber regression M-est: model-based cov")
  +     # design-based covariance matrix (test against fixed values)
  +     ref <- matrix(c(61704317.084028, -3367311.1407500, -3367311.140750,
  +         870905.6027803), ncol = 2)
  +     colnames(ref) <- rownames(ref) <- c("(Intercept)", "employment")
  +     all_equal(vcov(est), ref,
  +         label = "Huber regression M-est: design-based cov")
  + 
  +     #---------------------------------------------------------------------------
  +     # Tukey regression M-estimator
  +     est <- svyreg_tukeyM(payroll ~ employment, dn, k = 4.6, tol = 1e-9)
  +     ref <- rlm_mod(x = est$model$x, y = est$model$y, weights = est$model$w,
  +         c = 4.6, psi = psi.bisquare,  scale.est = "MAD", method = "M",
  +         wt.method = "case", acc = 1e-9, maxit = 50, test.vec = "coef")
  +     all_equal(coef(est), coef(ref), tolerance = 1e-7,
  +         label = "Tukey regression regression: coefficients")
  +     # model-based covariance matrix
  +     all_equal(vcov(est, "model"), vcov(ref), tolerance = 1e-7,
  +         label = "Tukey regression M-est: model-based cov")
  +     # design-based covariance matrix (test against fixed values)
  +     ref <- matrix(c(43513457.824988, -1421272.808406, -1421272.808406,
  +         589121.583227), ncol = 2)
  +     colnames(ref) <- rownames(ref) <- c("(Intercept)", "employment")
  +     all_equal(vcov(est), ref,
  +         label = "Tukey regression M-est: design-based cov")
  + }
  > 
  > #===============================================================================
  > # 2 Test against survey::svyglm
  > #===============================================================================
  > # survey weighted regression
  > ref <- svyglm(payroll ~ employment, dn)
  > est <- svyreg(payroll ~ employment, dn)
  > 
  > all_equal(coef(est), coef(ref),
  +     label = "Survey weighted regression: coefficients")
  > # design-based covariance matrix (test against survey::svyglm)
  > all_equal(vcov(est, "design"), vcov(ref),
  +     label = "Survey weighted regression: design-based cov")
  > # model-based covariance matrix (test against MASS::rlm)
  > ref <- rlm(payroll ~ employment, weights = weights(dn), data = dn$variables,
  +     method = "M", wt.method = "case", k = 1000)
  > all_equal(vcov(est, "model"), vcov(ref),
  +     label = "Survey weighted regression: model-based cov")
  > 
  > #-------------------------------------------------------------------------------
  > # Mallows regression GM-estimator: Huber psi
  > xwgt <- simpsonWgt(employment / weighted_median(employment, weight), a = 1,
  +     b = 20)
  > est <- svyreg_huberGM(payroll ~ employment, dn, k = 1.345, type = "Mallows",
  +     xwgt = xwgt)
  > all_equal(coef(est),
  +     c("(Intercept)" = 7843.4909, employment = 14425.3658),
  +     tolerance = 1e-4, label = "Huber Mallows regression GM-est")
  > # model-based covariance matrix
  > ref <- matrix(c(35768.4163, -572.71408, -572.71408, 52.40226), ncol = 2)
  > colnames(ref) <- c("(Intercept)", "employment")
  > rownames(ref) <- c("(Intercept)", "employment")
  > all_equal(vcov(est, "model"), ref, tolerance = 1e-4,
  +     label = "Huber Mallows regression GM-est: model-based cov")
  > 
  > #-------------------------------------------------------------------------------
  > # Schweppe regression GM-estimator: Huber psi
  > xwgt <- simpsonWgt(employment / weighted_median(employment, weight), a = 1,
  +     b = 20)
  > est <- svyreg_huberGM(payroll ~ employment, dn, k = 1.345, type = "Schweppe",
  +     xwgt = xwgt)
  > all_equal(coef(est),
  +     c("(Intercept)" = 7769.1899, employment = 14424.6617),
  +     tolerance = 1e-4, label = "Huber Schweppe regression GM-est")
  > # model-based covariance matrix
  > ref <- matrix(c(35697.6180, -620.1554, -620.1554, 57.4449), ncol = 2)
  > colnames(ref) <- c("(Intercept)", "employment")
  > rownames(ref) <- c("(Intercept)", "employment")
  > all_equal(vcov(est, "model"), ref, tolerance = 1e-4,
  +     label = "Huber Schweppe regression GM-est: model-based cov")
  > 
  > proc.time()
     user  system elapsed 
     1.75    0.15    1.90