The model_parameters() function also allows the computation of standard errors, confidence intervals, and p-values based on various covariance matrices: heteroskedasticity-consistent, cluster-robust, bootstrap, etc. This functionality relies on the sandwich and clubSandwich packages. This means that all models supported by either of these packages should work with model_parameters().

The computation of robust standard errors is controlled by two arguments:

1. vcov: accepts 3 types of arguments
• A covariance matrix
• A function which returns a covariance matrix (e.g., stats::vcov())
• A string which indicates the kind of uncertainty estimates to return.
• Heteroskedasticity-consistent: "vcovHC", "HC", "HC0", "HC1", "HC2", "HC3", "HC4", "HC4m", "HC5"
• Cluster-robust: "vcovCR", "CR0", "CR1", "CR1p", "CR1S", "CR2", "CR3"
• Bootstrap: "vcovBS", "xy", "residual", "wild", "mammen", "webb"
• Other sandwich functions: "vcovHAC", "vcovPC", "vcovCL", "vcovPL"
2. vcov_args: list of arguments passed to the sandwich or clubSandwich function used to compute the covariance matrix. See for example ?sandwich::vcovHAC

## Linear Regression Models

### Robust Covariance Matrix Estimation from Model Parameters

There are two arguments (see ?standard_error for further details) that allow for choosing different methods and options of robust estimation: - vcov - vcov_args

Let us start with a simple example, which uses a heteroskedasticity-consistent covariance matrix estimation with estimation-type “HC3” (i.e. sandwich::vcovHC(type = "HC3")).

First let’s create a simple linear regression model, which we know violates homoscedasticity assumption, and thus robust estimation methods are to be considered.

data(cars)
model <- lm(dist ~ speed, data = cars)

library(performance)
check_heteroscedasticity(model)
#> Warning: Heteroscedasticity (non-constant error variance) detected (p = 0.031).

We would extract model parameters both with and without robust estimation to highlight the difference it makes to standard errors, confidence intervals, t-statistic, and p-values. Also, note that the coefficient estimate remains unchanged.

# model parameters, where SE, CI, and p-values are *not* based on robust estimation
model_parameters(model)
#> Parameter   | Coefficient |   SE |          95% CI | t(48) |      p
#> -------------------------------------------------------------------
#> (Intercept) |      -17.58 | 6.76 | [-31.17, -3.99] | -2.60 | 0.012
#> speed       |        3.93 | 0.42 | [  3.10,  4.77] |  9.46 | < .001

# model parameters, where SE, CI, and p-values are based on robust estimation
mp <- model_parameters(model, vcov = "HC3")
mp
#> Parameter   | Coefficient |   SE |          95% CI | t(48) |      p
#> -------------------------------------------------------------------
#> (Intercept) |      -17.58 | 5.93 | [-29.51, -5.65] | -2.96 | 0.005
#> speed       |        3.93 | 0.43 | [  3.07,  4.79] |  9.20 | < .001

# compare standard errors to result from sandwich-package
mp$SE #> [1] 5.93 0.43 unname(sqrt(diag(sandwich::vcovHC(model)))) #> [1] 5.93 0.43 ### Cluster-Robust Covariance Matrix Estimation (sandwich) If a different type of covariance matrix estimation is required, use the vcov-argument. This argument accepts the name of a function from the sandwich or clubSandwich packages as a string, such as "vcovCL" (or just its suffix "CL"). parameters will then call the corresponding function with the content of vcov_args as arguments. The specific estimation type can be controlled by passing a type argument via vcov_args. See ?sandwich::vcovCL for information about the different types of covariance matrices that this function can produce (HC0 to HC3). In the next example, we use a clustered covariance matrix estimation with HC1-estimation type. # let's create a more complicated model data(iris) model <- lm(Petal.Length ~ Sepal.Length * Species + Sepal.Width, data = iris) # change estimation-type mp <- model_parameters( model, vcov = "CL", # type of covariance matrix vcov_args = list(type = "HC1") # type of robust estimation ) mp #> Parameter | Coefficient | SE | 95% CI | t(143) | p #> ------------------------------------------------------------------------------------------ #> (Intercept) | 0.87 | 0.42 | [ 0.03, 1.70] | 2.05 | 0.042 #> Sepal Length | 0.04 | 0.11 | [-0.18, 0.26] | 0.40 | 0.692 #> Species [versicolor] | -0.78 | 0.65 | [-2.07, 0.51] | -1.19 | 0.237 #> Species [virginica] | -0.41 | 0.59 | [-1.57, 0.75] | -0.70 | 0.483 #> Sepal Width | 0.11 | 0.08 | [-0.05, 0.27] | 1.38 | 0.170 #> Sepal Length * Species [versicolor] | 0.61 | 0.12 | [ 0.37, 0.85] | 4.96 | < .001 #> Sepal Length * Species [virginica] | 0.68 | 0.11 | [ 0.46, 0.90] | 6.15 | < .001 # compare standard errors to result from sandwich-package mp$SE
#> [1] 0.422 0.111 0.653 0.587 0.079 0.123 0.111
unname(sqrt(diag(sandwich::vcovCL(model))))
#> [1] 0.422 0.111 0.653 0.587 0.079 0.123 0.111

Usually, clustered covariance matrix estimation is used when there is a cluster-structure in the data. The variable indicating the cluster-structure can be defined in sandwich::vcovCL() with the cluster-argument. In model_parameters(), additional arguments that should be passed down to functions from the sandwich package can be specified in vcov_args:

iris$cluster <- factor(rep(LETTERS[1:8], length.out = nrow(iris))) # change estimation-type, defining additional arguments mp <- model_parameters( model, vcov = "vcovCL", vcov_args = list(type = "HC1", cluster = iris$cluster)
)

mp
#> Parameter                           | Coefficient |   SE |        95% CI | t(143) |      p
#> ------------------------------------------------------------------------------------------
#> (Intercept)                         |        0.87 | 0.34 | [ 0.20, 1.53] |   2.57 | 0.011
#> Sepal Length                        |        0.04 | 0.07 | [-0.10, 0.19] |   0.61 | 0.540
#> Species [versicolor]                |       -0.78 | 0.52 | [-1.80, 0.25] |  -1.49 | 0.137
#> Species [virginica]                 |       -0.41 | 0.26 | [-0.94, 0.11] |  -1.56 | 0.120
#> Sepal Width                         |        0.11 | 0.07 | [-0.03, 0.25] |   1.52 | 0.131
#> Sepal Length * Species [versicolor] |        0.61 | 0.10 | [ 0.42, 0.80] |   6.29 | < .001
#> Sepal Length * Species [virginica]  |        0.68 | 0.05 | [ 0.58, 0.78] |  13.28 | < .001

# compare standard errors to result from sandwich-package
mp$SE #> [1] 0.337 0.072 0.519 0.264 0.072 0.097 0.051 unname(sqrt(diag(sandwich::vcovCL(model, cluster = iris$cluster))))
#> [1] 0.337 0.072 0.519 0.264 0.072 0.097 0.051

### Cluster-Robust Covariance Matrix Estimation (clubSandwich)

using clubSandwich::vcovCR(). Thus, when vcov = "CR", the related function from the clubSandwich package is called. Note that this function requires the specification of the cluster-argument.

# create fake-cluster-variable, to demonstrate cluster robust standard errors
iris$cluster <- factor(rep(LETTERS[1:8], length.out = nrow(iris))) # cluster-robust estimation mp <- model_parameters( model, vcov = "vcovCR", vcov_args = list(type = "CR1", cluster = iris$cluster)
)
mp
#> Parameter                           | Coefficient |   SE |        95% CI | t(143) |      p
#> ------------------------------------------------------------------------------------------
#> (Intercept)                         |        0.87 | 0.33 | [ 0.21, 1.52] |   2.62 | 0.010
#> Sepal Length                        |        0.04 | 0.07 | [-0.10, 0.18] |   0.63 | 0.531
#> Species [versicolor]                |       -0.78 | 0.51 | [-1.78, 0.23] |  -1.53 | 0.129
#> Species [virginica]                 |       -0.41 | 0.26 | [-0.92, 0.10] |  -1.60 | 0.112
#> Sepal Width                         |        0.11 | 0.07 | [-0.03, 0.25] |   1.55 | 0.123
#> Sepal Length * Species [versicolor] |        0.61 | 0.09 | [ 0.42, 0.79] |   6.42 | < .001
#> Sepal Length * Species [virginica]  |        0.68 | 0.05 | [ 0.58, 0.78] |  13.56 | < .001

# compare standard errors to result from clubSsandwich-package
mp$SE #> [1] 0.330 0.070 0.508 0.259 0.071 0.095 0.050 unname(sqrt(diag(clubSandwich::vcovCR(model, type = "CR1", cluster = iris$cluster))))
#> [1] 0.330 0.070 0.508 0.259 0.071 0.095 0.050

### Robust Covariance Matrix Estimation on Standardized Model Parameters

Finally, robust estimation can be combined with standardization. However, robust covariance matrix estimation only works for standardize = "refit".

# model parameters, robust estimation on standardized model
model_parameters(model, standardize = "refit", vcov = "HC3")
#> Parameter                           | Coefficient |   SE |         95% CI | t(143) |      p
#> -------------------------------------------------------------------------------------------
#> (Intercept)                         |       -1.30 | 0.07 | [-1.44, -1.16] | -18.70 | < .001
#> Sepal Length                        |        0.02 | 0.06 | [-0.09,  0.13] |   0.37 | 0.711
#> Species [versicolor]                |        1.57 | 0.09 | [ 1.40,  1.74] |  17.84 | < .001
#> Species [virginica]                 |        2.02 | 0.09 | [ 1.84,  2.20] |  22.49 | < .001
#> Sepal Width                         |        0.03 | 0.02 | [-0.01,  0.07] |   1.32 | 0.190
#> Sepal Length * Species [versicolor] |        0.28 | 0.06 | [ 0.16,  0.41] |   4.65 | < .001
#> Sepal Length * Species [virginica]  |        0.32 | 0.06 | [ 0.21,  0.43] |   5.75 | < .001

## Linear Mixed-Effects Regression Models

### Robust Covariance Matrix Estimation for Mixed Models

For linear mixed-effects models, that by definition have a clustered (hierarchical or multilevel) structure in the data, it is also possible to estimate a cluster-robust covariance matrix. This is possible due to the clubSandwich package, thus we need to define the same arguments as in the above example.

library(lme4)
data(iris)
set.seed(1234)
iris$grp <- as.factor(sample(1:3, nrow(iris), replace = TRUE)) # fit example model model <- lme4::lmer( Sepal.Length ~ Species * Sepal.Width + Petal.Length + (1 | grp), data = iris ) # model parameters without robust estimation model_parameters(model) #> # Fixed Effects #> #> Parameter | Coefficient | SE | 95% CI | t(141) | p #> ------------------------------------------------------------------------------------------ #> (Intercept) | 1.55 | 0.40 | [ 0.76, 2.35] | 3.87 | < .001 #> Species [versicolor] | 0.41 | 0.55 | [-0.67, 1.50] | 0.75 | 0.454 #> Species [virginica] | -0.41 | 0.58 | [-1.56, 0.74] | -0.70 | 0.483 #> Sepal Width | 0.66 | 0.11 | [ 0.44, 0.89] | 5.83 | < .001 #> Petal Length | 0.82 | 0.07 | [ 0.69, 0.95] | 12.52 | < .001 #> Species [versicolor] * Sepal Width | -0.48 | 0.19 | [-0.85, -0.12] | -2.60 | 0.010 #> Species [virginica] * Sepal Width | -0.36 | 0.18 | [-0.71, 0.00] | -1.99 | 0.048 #> #> # Random Effects #> #> Parameter | Coefficient | SE | 95% CI #> ------------------------------------------------------- #> SD (Intercept: grp) | 0.08 | 0.05 | [0.02, 0.29] #> SD (Residual) | 0.30 | 0.02 | [0.27, 0.33] # model parameters with cluster robust estimation model_parameters( model, vcov = "vcovCR", vcov_args = list(type = "CR1", cluster = iris$grp)
)
#> # Fixed Effects
#>
#> Parameter                          | Coefficient |   SE |         95% CI | t(141) |      p
#> ------------------------------------------------------------------------------------------
#> (Intercept)                        |        1.55 | 0.40 | [ 0.76,  2.35] |   3.87 | < .001
#> Species [versicolor]               |        0.41 | 0.80 | [-1.17,  1.99] |   0.51 | 0.608
#> Species [virginica]                |       -0.41 | 0.19 | [-0.78, -0.03] |  -2.15 | 0.033
#> Sepal Width                        |        0.66 | 0.10 | [ 0.46,  0.86] |   6.64 | < .001
#> Petal Length                       |        0.82 | 0.05 | [ 0.72,  0.91] |  17.27 | < .001
#> Species [versicolor] * Sepal Width |       -0.48 | 0.35 | [-1.18,  0.21] |  -1.37 | 0.172
#> Species [virginica] * Sepal Width  |       -0.36 | 0.11 | [-0.57, -0.15] |  -3.39 | < .001
#>
#> # Random Effects
#>
#> Parameter           | Coefficient |   SE |       95% CI
#> -------------------------------------------------------
#> SD (Intercept: grp) |        0.08 | 0.05 | [0.02, 0.29]
#> SD (Residual)       |        0.30 | 0.02 | [0.27, 0.33]

Notice that robust estimation returns different standard errors, confidence intervals, test statistic and p-values compared to the standard estimation. Also, note that the coefficient estimate remains unchanged.

### Robust Covariance Matrix Estimation on Standardized Mixed Model Parameters

Once again, robust estimation can be combined with standardization for linear mixed-effects models as well and works only with standardize = "refit".

# model parameters, cluster robust estimation of standardized mixed model
model_parameters(
model,
standardize = "refit",
vcov = "vcovCR",
vcov_args = list(type = "CR1", cluster = iris\$grp)
)
#> # Fixed Effects
#>
#> Parameter                          | Coefficient |   SE |         95% CI | t(141) |      p
#> ------------------------------------------------------------------------------------------
#> (Intercept)                        |        0.97 | 0.08 | [ 0.82,  1.12] |  12.55 | < .001
#> Species [versicolor]               |       -1.29 | 0.33 | [-1.95, -0.63] |  -3.85 | < .001
#> Species [virginica]                |       -1.81 | 0.23 | [-2.26, -1.37] |  -8.01 | < .001
#> Sepal Width                        |        0.35 | 0.05 | [ 0.24,  0.45] |   6.64 | < .001
#> Petal Length                       |        1.74 | 0.10 | [ 1.54,  1.94] |  17.27 | < .001
#> Species [versicolor] * Sepal Width |       -0.25 | 0.19 | [-0.62,  0.11] |  -1.37 | 0.172
#> Species [virginica] * Sepal Width  |       -0.19 | 0.06 | [-0.30, -0.08] |  -3.39 | < .001
#>
#> # Random Effects
#>
#> Parameter           | Coefficient |   SE |       95% CI
#> -------------------------------------------------------
#> SD (Intercept: grp) |        0.10 | 0.06 | [0.03, 0.35]
#> SD (Residual)       |        0.36 | 0.02 | [0.32, 0.40]

Notice how drastically some of the p-values change between robust-unstandardized model and robust-standardized model.