greg.R
In robsurvey: Robust Survey Statistics Estimation

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "",
    prompt = TRUE
)

## ---- eval=FALSE--------------------------------------------------------------
#  library("robsurvey", quietly = TRUE)
#  library("survey")
#  data("MU284pps")

## ---- echo = FALSE------------------------------------------------------------
library(robsurvey, quietly = TRUE)
suppressPackageStartupMessages(library(survey))

## -----------------------------------------------------------------------------
dn <- svydesign(ids = ~LABEL, fpc = ~pi, data = MU284pps, pps = "brewer")

## -----------------------------------------------------------------------------
dn

## ---- fig.show="hold", out.width="80%"----------------------------------------
svyplot(RMT85 ~ P85, dn, xlab = "P85", ylab = "RMT85", inches = 0.1)

## -----------------------------------------------------------------------------
rat <- svyratio_huber(~RMT85, ~P85, dn, k = Inf)
rat

## -----------------------------------------------------------------------------
tot <- svytotal_ratio(rat, total = 8339)
tot

## -----------------------------------------------------------------------------
mse(tot) / 1e6

## -----------------------------------------------------------------------------
rat_rob <- svyratio_huber(~RMT85, ~P85, dn, k = 20)
tot_rob <- svytotal_ratio(rat_rob, total = 8339)
tot_rob

## -----------------------------------------------------------------------------
mse(tot_rob) / 1e6

## -----------------------------------------------------------------------------
wls <- svyreg(RMT85 ~ P85 + SS82, dn)
wls

## -----------------------------------------------------------------------------
summary(wls)

## ---- fig.show="hold", out.width="50%"----------------------------------------
plot(wls)

## -----------------------------------------------------------------------------
tot <- svytotal_reg(wls, totals = c(P85 = 8339, SS82 = 6301), N = 284, type = "ADU")
tot

## -----------------------------------------------------------------------------
mse(tot) / 1e6

## -----------------------------------------------------------------------------
rob <- svyreg_tukeyM(RMT85 ~ P85 + SS82, dn, k = 15)
rob

## ---- fig.show="hold", out.width="50%"----------------------------------------
plot(rob)

## -----------------------------------------------------------------------------
tot <- svytotal_reg(rob, totals = c(P85 = 8339, SS82 = 6301), N = 284, type = "huber", k = 50)
tot

## -----------------------------------------------------------------------------
mse(tot) / 1e6