diff --git a/DESCRIPTION b/DESCRIPTION index 6ae8c981b..aa72d4252 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Type: Package Package: parameters Title: Processing of Model Parameters -Version: 0.21.2.7 +Version: 0.21.2.8 Authors@R: c(person(given = "Daniel", family = "Lüdecke", diff --git a/NEWS.md b/NEWS.md index 610394769..b613a990b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -20,6 +20,8 @@ * Minor fixes for `print_html()` method for `model_parameters()`. +* Robust standard errors (argument `vcov`) now works for `plm` models. + # parameters 0.21.2 ## Changes diff --git a/R/methods_plm.R b/R/methods_plm.R index 4a93da5d3..6704b1445 100644 --- a/R/methods_plm.R +++ b/R/methods_plm.R @@ -15,12 +15,41 @@ degrees_of_freedom.plm <- function(model, method = "wald", ...) { #' @export -standard_error.plm <- function(model, ...) { - se <- stats::coef(summary(model)) +standard_error.plm <- function(model, vcov = NULL, vcov_args = NULL, verbose = TRUE, ...) { + dots <- list(...) + se <- NULL + se_standard <- stats::coef(summary(model)) + + # vcov: matrix + if (is.matrix(vcov)) { + se <- sqrt(diag(vcov)) + } + + # vcov: function which returns a matrix + if (is.function(vcov)) { + args <- c(list(model), vcov_args, dots) + se <- .safe(sqrt(diag(do.call("vcov", args)))) + } + + # vcov: character (with backward compatibility for `robust = TRUE`) + if (is.character(vcov) || isTRUE(dots[["robust"]])) { + .vcov <- insight::get_varcov( + model, + vcov = vcov, + vcov_args = vcov_args, + verbose = verbose, + ... + ) + se <- sqrt(diag(.vcov)) + } + + if (is.null(se)) { + se <- as.vector(se_standard[, 2]) + } .data_frame( - Parameter = .remove_backticks_from_string(rownames(se)), - SE = as.vector(se[, 2]) + Parameter = .remove_backticks_from_string(rownames(se_standard)), + SE = se ) } diff --git a/tests/testthat/_snaps/plm.md b/tests/testthat/_snaps/plm.md new file mode 100644 index 000000000..8490ec38b --- /dev/null +++ b/tests/testthat/_snaps/plm.md @@ -0,0 +1,34 @@ +# vcov standard errors + + Code + print(model_parameters(ran)) + Output + # Fixed Effects + + Parameter | Coefficient | SE | 95% CI | z | df | p + ----------------------------------------------------------------------------- + (Intercept) | 817.10 | 188.26 | [445.84, 1188.36] | 4.34 | 197 | < .001 + capital | -0.58 | 0.15 | [ -0.88, -0.29] | -3.95 | 197 | < .001 + inv | 2.92 | 0.30 | [ 2.33, 3.51] | 9.80 | 197 | < .001 + Message + + Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed + using a Wald z-distribution approximation. + +--- + + Code + print(model_parameters(ran, vcov = "HC1")) + Output + # Fixed Effects + + Parameter | Coefficient | SE | 95% CI | z | df | p + ---------------------------------------------------------------------------- + (Intercept) | 817.10 | 274.75 | [275.28, 1358.92] | 2.97 | 197 | 0.003 + capital | -0.58 | 0.43 | [ -1.43, 0.26] | -1.37 | 197 | 0.173 + inv | 2.92 | 0.89 | [ 1.16, 4.67] | 3.28 | 197 | 0.001 + Message + + Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed + using a Wald z-distribution approximation. + diff --git a/tests/testthat/test-plm.R b/tests/testthat/test-plm.R index 380e0bb29..1dae25c3a 100644 --- a/tests/testthat/test-plm.R +++ b/tests/testthat/test-plm.R @@ -33,6 +33,7 @@ m3 <- suppressWarnings(plm::plm( index = c("ID", "Year") )) + test_that("ci", { expect_equal( ci(m1)$CI_low, @@ -47,6 +48,7 @@ test_that("ci", { expect_equal(ci(m3)$CI_low, -2.60478, tolerance = 1e-3) }) + test_that("se", { expect_equal( standard_error(m1)$SE, @@ -61,6 +63,7 @@ test_that("se", { expect_equal(standard_error(m3)$SE, 0.5166726, tolerance = 1e-3) }) + test_that("p_value", { expect_equal( p_value(m1)$p, @@ -75,6 +78,7 @@ test_that("p_value", { expect_equal(p_value(m3)$p, 0.53696, tolerance = 1e-3) }) + test_that("model_parameters", { expect_equal( model_parameters(m1)$Coefficient, @@ -88,3 +92,22 @@ test_that("model_parameters", { ) expect_equal(model_parameters(m3)$Coefficient, -0.381721, tolerance = 1e-3) }) + + +test_that("vcov standard errors", { + skip_if_not_installed("sandwich") + data("Grunfeld", package = "plm") + ran <- suppressWarnings( + plm::plm(value ~ capital + inv, data = Grunfeld, model = "random", effect = "twoways") + ) + out1 <- standard_error(ran) + out2 <- standard_error(ran, vcov = "HC1") + validate1 <- coef(summary(ran))[, 2] + validate2 <- sqrt(diag(sandwich::vcovHC(ran, type = "HC1"))) + + expect_equal(out1$SE, validate1, tolerance = 1e-3, ignore_attr = TRUE) + expect_equal(out2$SE, validate2, tolerance = 1e-3, ignore_attr = TRUE) + + expect_snapshot(print(model_parameters(ran))) + expect_snapshot(print(model_parameters(ran, vcov = "HC1"))) +})