-
Notifications
You must be signed in to change notification settings - Fork 0
/
test-mahalanobis_multi_missing.R
49 lines (41 loc) · 1.35 KB
/
test-mahalanobis_multi_missing.R
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
skip_if_not_installed("modi")
skip_if_not_installed("norm2")
library(testthat)
library(lavaan)
# Can handle missing data
mod <-
'
m1 ~ iv1 + c(a1, a2) * iv2
dv ~ c(b, b) * m1
a2b := a2 * b
'
dat <- pa_dat
dat0 <- dat[1:40, ]
set.seed(856041)
dat0$gp <- sample(c("gp2", "gp1"), size = nrow(dat0), replace = TRUE)
dat0[1, 2] <- dat0[2, 3] <- dat0[3, 4] <- dat0[4, 1:4] <- NA
head(dat0)
suppressWarnings(fit0 <- lavaan::sem(mod, dat0, missing = "fiml.x",
group = "gp"))
fit0_data <- lav_data_used(fit0)
head(fit0_data)
fit0_data_g1 <- fit0_data[fit0_data$gp == "gp2", -5]
fit0_data_g2 <- fit0_data[fit0_data$gp == "gp1", -5]
head(fit0_data_g1)
head(fit0_data_g2)
em_out_g1 <- norm2::emNorm(fit0_data_g1, estimate.worst = FALSE, criterion = 1e-6)
md_check_g1 <- modi::MDmiss(fit0_data_g1,
em_out_g1$param$beta,
em_out_g1$param$sigma)
em_out_g2 <- norm2::emNorm(fit0_data_g2, estimate.worst = FALSE, criterion = 1e-6)
md_check_g2 <- modi::MDmiss(fit0_data_g2,
em_out_g2$param$beta,
em_out_g2$param$sigma)
md_check <- c(md_check_g1, md_check_g2)
md_rerun <- mahalanobis_rerun(fit0)
test_that("Compare Mahalanobis distances", {
expect_equal(ignore_attr = TRUE,
sort(as.vector(md_rerun)),
sort(md_check)
)
})