---
title: "Practical example: (Cluster) Robust Standard Errors"
author: "Daniel Lüdecke"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Practical example: (Cluster) Robust Standard Errors}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r set-options, echo = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", dev = "png", fig.width = 7, fig.height = 5, message = FALSE, warning = FALSE)
if (!requireNamespace("sandwich", quietly = TRUE) ||
    !requireNamespace("clubSandwich", quietly = TRUE)) {
  knitr::opts_chunk$set(eval = FALSE)
}
```

This vignette demonstrate how to compute confidence intervals based on (cluster) robust variance-covariance matrices for standard errors.

First, we load the required packages and create a sample data set with a binomial and continuous variable as predictor as well as a group factor. To avoid convergence warnings, the continuous variable is standardized.

```{r}
library(ggeffects)
set.seed(123)

# example taken from "?clubSandwich::vcovCR"
m <- 8
cluster <- factor(rep(LETTERS[1:m], 3 + rpois(m, 5)))
n <- length(cluster)
X <- matrix(rnorm(3 * n), n, 3)
nu <- rnorm(m)[cluster]
e <- rnorm(n)
y <- X %*% c(.4, .3, -.3) + nu + e
dat <- data.frame(y, X, cluster, row = 1:n)

# fit linear model
model <- lm(y ~ X1 + X2 + X3, data = dat)
```

## Predictions with normal standard errors

In this example, we use the normal standard errors, as returned by `predict()`, to compute confidence intervals.

```{r message = TRUE}
ggpredict(model, "X1")
```

```{r message = FALSE}
me <- ggpredict(model, "X1")
plot(me)
```

## Predictions with HC-estimated standard errors

Now, we use `sandwich::vcovHC()` to estimate heteroskedasticity-consistent standard errors. To do so, first the function name, `vcovHC()`, must be supplied to the `vcov.fun`-argument. `sandwich::vcovHC()`, in turn, has different types of estimation. This must be specified in `vcov.type`.

```{r message = TRUE}
ggpredict(model, "X1", vcov.fun = "vcovHC", vcov.type = "HC0")
```

```{r message = FALSE}
me <- ggpredict(model, "X1", vcov.fun = "vcovHC", vcov.type = "HC0")
plot(me)
```

## Predictions with cluster-robust standard errors

The last example shows how to define cluster-robust standard errors. These are based on `clubSandwich::vcovCR()`. Thus, `vcov.fun = "vcovCR"` is always required when estimating cluster robust standard errors. `clubSandwich::vcovCR()` has also different estimation types, which must be specified in `vcov.type`. Furthermore, `clubSandwich::vcovCR()` *requires* the `cluster`-argument, which must be specified in `vcov.args`:

```{r message = TRUE}
ggpredict(
  model, "X1", vcov.fun = "vcovCR", vcov.type = "CR0", 
  vcov.args = list(cluster = dat$cluster)
)
```

```{r message = FALSE}
me <- ggpredict(
  model, "X1", vcov.fun = "vcovCR", vcov.type = "CR0", 
  vcov.args = list(cluster = dat$cluster)
)
plot(me)
```

