1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
|
skip_if_not_installed("sandwich")
skip_if_not_installed("clubSandwich")
data(mtcars)
mtcars$am <- as.factor(mtcars$am)
model <- lm(mpg ~ wt * am + cyl + gear, data = mtcars)
test_that("model_parameters, robust CL", {
params1 <- model_parameters(
model,
vcov = "CL",
vcov_args = list(type = "HC1"),
verbose = FALSE
)
robust_se <- unname(sqrt(diag(sandwich::vcovCL(model))))
expect_equal(params1$SE, robust_se, tolerance = 1e-3)
expect_equal(params1$p, c(0, 0.00695, 0.00322, 0.00435, 0.94471, 0.00176), tolerance = 1e-3)
})
test_that("model_parameters, robust", {
params <- model_parameters(model, vcov = "HC", verbose = FALSE)
robust_se <- unname(sqrt(diag(sandwich::vcovHC(model))))
expect_equal(params$SE, robust_se, tolerance = 1e-3)
expect_equal(params$p, c(0, 0.0259, 0.01478, 0.01197, 0.95238, 0.01165), tolerance = 1e-3)
})
test_that("ci, robust", {
params <- ci(model, vcov = "HC", verbose = FALSE)
robust_se <- unname(sqrt(diag(sandwich::vcovHC(model))))
upper_ci <- as.vector(coef(model) + qt(0.975, df.residual(model)) * robust_se)
expect_equal(params$CI_high, upper_ci, tolerance = 1e-3, ignore_attr = TRUE)
})
test_that("model_parameters, robust CL", {
params <- model_parameters(model, vcov = "vcovCL", verbose = FALSE)
robust_se <- unname(sqrt(diag(sandwich::vcovCL(model))))
expect_equal(params$SE, robust_se, tolerance = 1e-3)
expect_equal(params$p, c(0, 0.00695, 0.00322, 0.00435, 0.94471, 0.00176), tolerance = 1e-3)
})
model2 <- lm(mpg ~ wt * am + cyl + gear, data = datawizard::standardize(mtcars))
test_that("model_parameters, robust", {
params <- model_parameters(model, standardize = "refit", vcov = "HC", verbose = FALSE)
robust_se <- unname(sqrt(diag(sandwich::vcovHC(model2))))
expect_equal(params$SE, robust_se, tolerance = 1e-3)
expect_equal(params$p, c(0.28624, 0.0259, 0.43611, 0.01197, 0.95238, 0.01165), tolerance = 1e-3)
})
# cluster-robust standard errors, using clubSandwich
data(iris)
model <- lm(Petal.Length ~ Sepal.Length * Species, data = iris)
iris$cluster <- factor(rep(LETTERS[1:8], length.out = nrow(iris)))
test_that("model_parameters, robust CR", {
params <- model_parameters(
model,
vcov = "CR1",
vcov_args = list(cluster = iris$cluster),
verbose = FALSE
)
robust_se <- unname(sqrt(diag(clubSandwich::vcovCR(model, type = "CR1", cluster = iris$cluster))))
expect_equal(params$SE, robust_se, tolerance = 1e-3)
expect_equal(params$p, c(0.01246, 0.04172, 0.18895, 0.57496, 0, 0), tolerance = 1e-3)
})
test_that("model_parameters, normal", {
params <- model_parameters(model)
expect_equal(params$p, c(0.13267, 0.21557, 0.36757, 0.77012, 3e-05, 0), tolerance = 1e-3)
})
data(mtcars)
mtcars$am <- as.factor(mtcars$am)
model <- lm(mpg ~ wt * am + cyl + gear, data = mtcars)
test_that("model_parameters, robust", {
params <- model_parameters(model, vcov = "HC3")
robust_se <- unname(sqrt(diag(sandwich::vcovHC(model))))
expect_equal(params$SE, robust_se, tolerance = 1e-3)
expect_equal(params$p, c(0, 0.0259, 0.01478, 0.01197, 0.95238, 0.01165), tolerance = 1e-3)
})
test_that("ci, robust", {
params <- ci(model, vcov = "HC3")
robust_se <- unname(sqrt(diag(sandwich::vcovHC(model))))
upper_ci <- as.vector(coef(model) + qt(0.975, df.residual(model)) * robust_se)
expect_equal(params$CI_high, upper_ci, tolerance = 1e-3, ignore_attr = TRUE)
})
test_that("model_parameters, robust CL", {
params <- model_parameters(model, vcov = "vcovCL", vcov_args = list(type = "HC1"))
robust_se <- unname(sqrt(diag(sandwich::vcovCL(model))))
expect_equal(params$SE, robust_se, tolerance = 1e-3)
expect_equal(params$p, c(0, 0.00695, 0.00322, 0.00435, 0.94471, 0.00176), tolerance = 1e-3)
})
d <- datawizard::standardize(mtcars)
model2 <- lm(mpg ~ wt * am + cyl + gear, data = d)
test_that("model_parameters, robust", {
params <- model_parameters(model, standardize = "refit", vcov = "HC3")
robust_se <- unname(sqrt(diag(sandwich::vcovHC(model2))))
expect_equal(params$SE, robust_se, tolerance = 1e-3)
expect_equal(params$p, c(0.28624, 0.0259, 0.43611, 0.01197, 0.95238, 0.01165), tolerance = 1e-3)
})
# cluster-robust standard errors, using clubSandwich
data(iris)
model <- lm(Petal.Length ~ Sepal.Length * Species, data = iris)
iris$cluster <- factor(rep(LETTERS[1:8], length.out = nrow(iris)))
test_that("model_parameters, robust CR", {
params <- model_parameters(model, vcov = "vcovCR", vcov_args = list(type = "CR1", cluster = iris$cluster))
robust_se <- unname(sqrt(diag(clubSandwich::vcovCR(model, type = "CR1", cluster = iris$cluster))))
expect_equal(params$SE, robust_se, tolerance = 1e-3)
expect_equal(params$p, c(0.01246, 0.04172, 0.18895, 0.57496, 0, 0), tolerance = 1e-3)
})
test_that("model_parameters, normal", {
params <- model_parameters(model)
expect_equal(params$p, c(0.13267, 0.21557, 0.36757, 0.77012, 3e-05, 0), tolerance = 1e-3)
})
test_that("ci_ml1, robust", {
skip("TODO: this one actually is not correct.")
skip_if_not(packageVersion("parameters") < "0.16.9.9")
skip_if_not_installed("lme4")
model <- lme4::lmer(Petal.Length ~ Sepal.Length + (1 | Species), data = iris)
params <- ci_ml1(model, vcov = "CR", vcov_args = list(cluster = iris$Species))
robust_se <- unname(sqrt(diag(clubSandwich::vcovCR(model, type = "CR1", cluster = iris$Species))))
upper_ci <- fixef(model) + qt(0.975, dof_ml1(model)) * robust_se
})
|
