## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5
)
options(hcinfer.use_emoji = FALSE)

## -----------------------------------------------------------------------------
library(hcinfer)

schools <- PublicSchools
schools$income_scaled <- schools$income / 10000
schools$income_scaled_sq <- schools$income_scaled^2

fit <- lm(expenditure ~ income_scaled + income_scaled_sq, data = schools)
fit

## -----------------------------------------------------------------------------
hc_methods()

## -----------------------------------------------------------------------------
result <- hcinfer(fit)
result

## -----------------------------------------------------------------------------
summary(result)

## -----------------------------------------------------------------------------
result_hc3 <- hcinfer(fit, type = "hc3")
result_hc3

## -----------------------------------------------------------------------------
tests(result)

## -----------------------------------------------------------------------------
tests(result, parm = "income_scaled_sq")
tests(result, parm = 3)

## -----------------------------------------------------------------------------
tests(result, alpha = 0.10)

## -----------------------------------------------------------------------------
confint(result)

## -----------------------------------------------------------------------------
confint(result, parm = "income_scaled_sq")
confint(result, parm = "income_scaled_sq", level = 0.90)

## -----------------------------------------------------------------------------
coef(result)

## -----------------------------------------------------------------------------
robust_vcov <- vcov(result)
robust_vcov

## -----------------------------------------------------------------------------
sqrt(diag(robust_vcov))

## -----------------------------------------------------------------------------
cov_hcbeta <- vcov_hc(fit)
cov_hcbeta

## -----------------------------------------------------------------------------
summary(cov_hcbeta)
vcov(cov_hcbeta)

## -----------------------------------------------------------------------------
cov_hc5 <- vcov_hc(fit, type = "hc5", k = 0.7)
cov_hc5

## ----fig.alt = "Robust confidence intervals for the public-schools regression coefficients."----
plot(result)

## ----fig.alt = "Robust confidence interval for the quadratic income coefficient."----
plot(result, parm = "income_scaled_sq")

## ----fig.alt = "HCbeta adjustment factors plotted against leverage values for the public-schools regression."----
plot(cov_hcbeta)

## ----fig.alt = "HC3 adjustment factors plotted against leverage values with the two largest weights labeled."----
plot(vcov_hc(fit, type = "hc3"), label_top = 2)

## -----------------------------------------------------------------------------
test_hcbeta <- tests(result, parm = "income_scaled_sq")
test_hc3 <- tests(result_hc3, parm = "income_scaled_sq")

comparison <- data.frame(
  estimator = c("HCbeta", "HC3"),
  estimate = c(test_hcbeta$estimate, test_hc3$estimate),
  robust_se = c(test_hcbeta$std_error, test_hc3$std_error),
  p_value = c(test_hcbeta$p_value, test_hc3$p_value)
)

comparison

## ----eval = FALSE-------------------------------------------------------------
# fit <- lm(y ~ x1 + x2, data = data)
# result <- hcinfer(fit, type = "hcbeta")
# 
# summary(result)
# tests(result)
# confint(result)
# plot(result)