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

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

schools <- PublicSchools2
schools$income_scaled <- schools$income / 10000

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

## -----------------------------------------------------------------------------
gls_fit <- gls_mult(fit)
gls_fit

## -----------------------------------------------------------------------------
summary(gls_fit)

## -----------------------------------------------------------------------------
coef(gls_fit)                        # mean coefficients
coef(gls_fit, model = "dispersion")  # log-variance coefficients

## -----------------------------------------------------------------------------
vcov(gls_fit)                        # mean covariance
vcov(gls_fit, model = "dispersion")  # dispersion covariance

## -----------------------------------------------------------------------------
tests(gls_fit)
confint(gls_fit)

## -----------------------------------------------------------------------------
head(fitted(gls_fit))
head(gls_fit$fitted_variances)   # estimated conditional variances
head(gls_fit$weights)            # GLS weights exp(-eta)

## -----------------------------------------------------------------------------
two_step <- gls_mult(fit, estimator = "two_step")
coef(two_step)

## -----------------------------------------------------------------------------
data.frame(
  term = names(coef(gls_fit)),
  ml = coef(gls_fit),
  two_step = coef(two_step)
)

## -----------------------------------------------------------------------------
gls_income <- gls_mult(fit, variance = ~ income_scaled)
coef(gls_income, model = "dispersion")

## -----------------------------------------------------------------------------
full <- gls_mult(fit)                       # variance ~ income + south
income_only <- gls_mult(fit, variance = ~ income_scaled)
homoskedastic <- gls_mult(fit, variance = ~ 1)
AIC(full, income_only, homoskedastic)
BIC(full, income_only, homoskedastic)

## -----------------------------------------------------------------------------
data.frame(
  term = names(coef(fit)),
  ols = sqrt(diag(vcov(fit))),
  fgls_ml = sqrt(diag(vcov(gls_fit))),
  hcbeta = sqrt(diag(vcov(hcinfer(fit, type = "hcbeta")))),
  hc3 = sqrt(diag(vcov(hcinfer(fit, type = "hc3"))))
)