R/weighted_mean_huber.R

In robsurvey: Robust Survey Statistics Estimation

# Huber M-estimator of the weighted mean
weighted_mean_huber <- function(x, w, k, type = "rwm", asym = FALSE,
    info = FALSE, na.rm = FALSE, verbose = TRUE, ...)
{
    type <- match.arg(type, c("rht", "rhj", "rwm"))
    if (type == "rhj")
        type <- "rwm"

    string <- paste0("Huber M-estimator (type = ", type, ifelse(asym,
        "; asym. psi", ""), ")")
    psi <- ifelse(asym, 1, 0)
    dat <- .check_data_weights(x, w, na.rm)
    # empty data
    if (is.null(dat))
        return(NA)
    # only one observation
    if (dat$n == 1) {
        if (info)
            return(list(characteristic = "mean", estimator = list(string =
                string, type = type, psi = psi, psi_fun = ifelse(asym,
                "asymHuber", "Huber"), k = k), estimate = dat$x, scale = NA,
                residuals = 0, model = list(y = dat$x, w = dat$w),
                call = match.call()))
        else
            return(dat$x)
    }
    # otherwise
    if (type == "rwm") {
        res <- robsvyreg(rep(1, dat$n), dat$x, dat$w, k, psi, 0, NULL, NULL,
            verbose, ...)
    }
    if (type == "rht") {
        xvar <- mean(dat$w) / dat$w
        res <- robsvyreg(xvar, dat$x, dat$w, k, psi, 0, NULL, xvar, verbose,
            ...)
    }
    # check for failure of convergence
    if (!res$optim$converged) {
        res$estimate <- NA
        res$scale <- NA
    }
    # return
    if (info) {
        res$model[c("n", "p")] <- NULL
        if (type == "rwm")
            res$model$x <- NULL
        res$characteristic <- "mean"
        res$estimator$string <- string
        res$estimator$type <- type
        res$call <- match.call()
        return(res)
    } else {
        return(res$estimate)
    }
}
# Huber M-estimator of the weighted total
weighted_total_huber <- function(x, w, k, type = "rwm", asym = FALSE,
    info = FALSE, na.rm = FALSE, verbose = TRUE, ...)
{
    res <- weighted_mean_huber(x, w, k, type, asym, info, na.rm, verbose, ...)
    if (length(res) == 1) {
        res <- res * sum(w)
    } else {
        res$characteristic <- "total"
        res$estimate <- res$estimate * sum(w)
        res$call <- match.call()
    }
    return(res)
}