library(knitr) options(knitr.kable.NA = "") options(digits = 2) knitr::opts_chunk$set( echo = TRUE, collapse = TRUE, warning = FALSE, message = FALSE, comment = "#>", out.width = "100%" ) if (!requireNamespace("poorman", quietly = TRUE) || !requireNamespace("clubSandwich", quietly = TRUE) || !requireNamespace("sandwich", quietly = TRUE) || !requireNamespace("lme4", quietly = TRUE)) { knitr::opts_chunk$set(eval = FALSE) } else { library(parameters) library(poorman) } set.seed(333)
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:
vcov
: accepts 3 types of argumentsstats::vcov()
)"vcovHC"
, "HC"
, "HC0"
, "HC1"
, "HC2"
, "HC3"
, "HC4"
, "HC4m"
, "HC5"
"vcovCR"
, "CR0"
, "CR1"
, "CR1p"
, "CR1S"
, "CR2"
, "CR3"
"vcovBS"
, "xy"
, "residual"
, "wild"
, "mammen"
, "webb"
sandwich
functions: "vcovHAC"
, "vcovPC"
, "vcovCL"
, "vcovPL"
vcov_args
: list of arguments passed to the sandwich or clubSandwich function used to compute the covariance matrix. See for example ?sandwich::vcovHAC
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)
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) # model parameters, where SE, CI, and p-values are based on robust estimation mp <- model_parameters(model, vcov = "HC3") mp # compare standard errors to result from sandwich-package mp$SE unname(sqrt(diag(sandwich::vcovHC(model))))
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 # compare standard errors to result from sandwich-package mp$SE unname(sqrt(diag(sandwich::vcovCL(model))))
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 # compare standard errors to result from sandwich-package mp$SE unname(sqrt(diag(sandwich::vcovCL(model, cluster = iris$cluster))))
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 # compare standard errors to result from clubSsandwich-package mp$SE unname(sqrt(diag(clubSandwich::vcovCR(model, type = "CR1", cluster = iris$cluster))))
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")
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) # model parameters with cluster robust estimation model_parameters( model, vcov = "vcovCR", vcov_args = list(type = "CR1", cluster = iris$grp) )
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.
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) )
Notice how drastically some of the p-values change between robust-unstandardized model and robust-standardized model.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.