performance/tests/testthat/test-test_likelihoodratio.R at b021586e48a4e1f2f37909b7fc7eaa4a77125e33 · easystats/performance · GitHub

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
test_that("test_likelihoodratio - regression models", {
  m1 <- lm(mpg ~ wt + cyl + gear + disp, data = mtcars)
  m2 <- lm(mpg ~ wt + cyl + gear, data = mtcars)
  m3 <- lm(mpg ~ wt + cyl, data = mtcars)

  # stats::anova()
  rez <- test_likelihoodratio(m1, m2, m3, estimator = "OLS")
  ref <- anova(m1, m2, m3, test = "LRT")
  expect_equal(ref$`Pr(>Chi)`, rez$p, tolerance = 1e-03)

  rez <- test_likelihoodratio(m3, m2, m1, estimator = "OLS")
  ref <- anova(m3, m2, m1, test = "LRT")
  expect_equal(ref$`Pr(>Chi)`, rez$p, tolerance = 1e-03)

  # lmtest::lrtest()
  rez <- test_likelihoodratio(m1, m2, m3, estimator = "ML")
  p <- c(NA, 0.5302030, 0.4747344) # lmtest::lrtest(m1, m2, m3)$`Pr(>Chisq)`
  expect_equal(p, rez$p, tolerance = 1e-03)

  rez <- test_likelihoodratio(m3, m2, m1, estimator = "ML")
  p <- c(NA, 0.4747344, 0.5302030) # lmtest::lrtest(m3, m2, m1)$`Pr(>Chisq)`
  expect_equal(p, rez$p, tolerance = 1e-03)
})

test_that("test_likelihoodratio - model names", {
  m1 <- lm(mpg ~ wt + cyl + gear + disp, data = mtcars)
  m2 <- lm(mpg ~ wt + cyl + gear, data = mtcars)
  m3 <- lm(mpg ~ wt + cyl, data = mtcars)

  rez <- test_likelihoodratio(m1, m2, m3)
  expect_equal(rez$Name, c("m1", "m2", "m3"))

  rez <- test_likelihoodratio(list(m1, m2, m3))
  expect_equal(rez$Name, c("Model 1", "Model 2", "Model 3"))

  models <- list(m1, m2, m3)
  rez <- test_likelihoodratio(models)
  expect_equal(rez$Name, c("Model 1", "Model 2", "Model 3"))
})

test_that("test_likelihoodratio - reversed order", {
  m1 <- lm(mpg ~ wt + cyl, data = mtcars)
  m2 <- lm(mpg ~ wt + cyl + gear, data = mtcars)
  m3 <- lm(mpg ~ wt + cyl + gear + disp, data = mtcars)

  rez <- test_likelihoodratio(m1, m2, m3, estimator = "OLS")
  ref <- anova(m1, m2, m3, test = "LRT")
  expect_equal(ref$`Pr(>Chi)`, rez$p, tolerance = 1e-03)

  rez <- test_likelihoodratio(m1, m2, m3, estimator = "ML")
  p <- c(NA, 0.4747344, 0.5302030) # lmtest::lrtest(m1, m2, m3)$`Pr(>Chisq)`
  expect_equal(p, rez$p, tolerance = 1e-03)
})

skip_on_cran()
skip_if_not_installed("lme4")

m1 <- suppressMessages(lme4::lmer(Sepal.Length ~ Petal.Width + (1 | Species), data = iris))
m2 <- suppressMessages(lme4::lmer(Sepal.Length ~ Petal.Width + Petal.Length + (1 | Species), data = iris))
m3 <- suppressMessages(lme4::lmer(Sepal.Length ~ Petal.Width * Petal.Length + (1 | Species), data = iris))

test_that("test_likelihoodratio - lme4 ML", {
  t1 <- test_lrt(m1, m2, m3)
  t2 <- suppressMessages(anova(m1, m2, m3))
  expect_equal(attributes(t1)$estimator, "ml")
  expect_equal(t1$Chi2, c(NA, 85.26365, 0.84141), tolerance = 1e-3)
  expect_equal(t1$p, c(NA, 0, 0.35899), tolerance = 1e-3)
  # close, but not the same
  expect_equal(t1$p, t2$`Pr(>Chisq)`, tolerance = 1e-1)
  expect_equal(t1$Chi2, t2$Chisq, tolerance = 1e-1)
})

test_that("test_likelihoodratio - lme4 OLS", {
  t2 <- test_lrt(m1, m2, m3, estimator = "ols")
  expect_equal(attributes(t2)$estimator, "ols")
  expect_equal(t2$Chi2, c(NA, 105.73844, 1.04346), tolerance = 1e-3)
  expect_equal(t2$p, c(NA, 0, 0.30702), tolerance = 1e-3)
})

test_that("test_likelihoodratio - lme4 REML", {
  expect_warning(t3 <- test_lrt(m1, m2, m3, estimator = "REML"))
  expect_equal(attributes(t3)$estimator, "reml")
  expect_equal(t3$Chi2, c(NA, 89.32933, 2.85635), tolerance = 1e-3)
  expect_equal(t3$p, c(NA, 0, 0.09101), tolerance = 1e-3)
})

m1 <- glm(am ~ mpg, data = mtcars, family = binomial())
m2 <- glm(am ~ mpg + hp, data = mtcars, family = binomial())
m3 <- glm(am ~ mpg + hp + vs, data = mtcars, family = binomial())

test_that("test_likelihoodratio - glm", {
  t1 <- anova(m1, m2, m3, test = "LRT")
  t2 <- test_lrt(m1, m2, m3)
  expect_equal(t1$`Pr(>Chi)`, t2$p, tolerance = 1e-3)
  expect_equal(t1$Deviance, t2$Chi2, tolerance = 1e-3)
  expect_equal(attributes(t2)$estimator, "ml")
})

m1 <- lm(Sepal.Length ~ Petal.Width, data = iris)
m2 <- lm(Sepal.Length ~ Petal.Width + Species, data = iris)
m3 <- lm(Sepal.Length ~ Petal.Width * Species, data = iris)

test_that("test_likelihoodratio - lm", {
  t1 <- anova(m1, m2, m3, test = "LRT")
  t2 <- test_lrt(m1, m2, m3)
  expect_equal(t1$`Pr(>Chi)`, t2$p, tolerance = 1e-3)
  expect_equal(attributes(t2)$estimator, "ols")
})

test_that("test_likelihoodratio - print p-digits", {
  data(mtcars)
  m1 <- lm(mpg ~ hp, data = mtcars)
  m2 <- lm(mpg ~ hp * as.factor(gear), data = mtcars)

  expect_snapshot(test_likelihoodratio(m1, m2))
  expect_snapshot(insight::print_md(test_likelihoodratio(m1, m2), p_digits = 3))
})