Merge pull request #113 from sfcheung/devel

Update to 0.1.6.2
sfcheung · Mar 15, 2024 · 0d96f55 · 0d96f55
2 parents bf1e481 + 512733f
commit 0d96f55
Show file tree

Hide file tree

Showing 17 changed files with 441 additions and 21 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: semfindr
 Title: Influential Cases in Structural Equation Modeling
-Version: 0.1.6.1
+Version: 0.1.6.2
 Authors@R: c(
     person(given = "Shu Fai",
            family = "Cheung",
@@ -25,13 +25,14 @@ License: GPL-3
 Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.1
 Suggests:
     testthat (>= 3.0.0),
     parallel,
     knitr,
     rmarkdown,
-    modi
+    modi,
+    MASS
 Imports:
     lavaan,
     ggplot2,

diff --git a/NEWS.md b/NEWS.md
@@ -1,10 +1,17 @@
-# semfindr 0.1.6.1
+# semfindr 0.1.6.2
+
+## Others
 
 - Updated the two tests for nonconvergence.
   Will check against results generated
   by directly calling `lavaan` functions.
   (0.1.6.1)
 
+## Bug Fixes
+
+- Fixed a bug with listwise deletion
+  in `lavaan_rerun()`. (0.1.6.2)
+
 # semfindr 0.1.6
 
 ## New Features

diff --git a/R/lavaan_rerun.R b/R/lavaan_rerun.R
@@ -212,6 +212,7 @@ lavaan_rerun <- function(fit,
       n <- nrow(lavaan::lavInspect(fit, "data"))
       n_j <- n
     }
+  n_orig <- sum(lavaan::lavInspect(fit, "norig"))
 
   if (is.null(case_id)) {
       case_ids <- lavaan::lavInspect(fit, "case.idx",
@@ -220,10 +221,11 @@ lavaan_rerun <- function(fit,
     } else {
       case_ids <- lavaan::lavInspect(fit, "case.idx",
                                     drop.list.single.group = FALSE)
-      if (length(case_id) != n) {
+      case_ids <- sort(unlist(case_ids, use.names = FALSE))
+      if (length(case_id) != n_orig) {
           stop("The length of case_id is not equal to the number of cases.")
         } else {
-          case_ids <- case_id
+          case_ids <- case_id[case_ids]
         }
     }
 
@@ -249,12 +251,17 @@ lavaan_rerun <- function(fit,
               stop("Some elements in to_rerun is not in the case_id vectors.")
             }
         } else {
-          if (!all(to_rerun %in% seq_len(n))) {
+          if (!all(to_rerun %in% seq_len(n_orig))) {
               stop("Some elements in to_rerun is not valid row numbers.")
             }
+          if (!all(to_rerun %in% case_ids)) {
+              stop("Some cases in to_rerun is not used in lavaan output. Probably due to listwise deletion.")
+            }
+          to_reun_org <- to_rerun
+          to_rerun <- match(to_rerun, case_ids)
         }
     } else {
-      to_rerun <- sort(unlist(case_ids, use.names = FALSE))
+      to_rerun <- order(unlist(case_ids, use.names = FALSE))
     }
 
   if (!missing(md_top)) {
@@ -291,15 +298,26 @@ lavaan_rerun <- function(fit,
       fit_resid_md <- unlist(fit_resid_md, use.names = FALSE)
       tmp1 <- lavaan::lavInspect(fit, "case.idx",
                                  drop.list.single.group = FALSE)
-      tmp2 <- unlist(tmp1, use.names = FALSE)
+      tmp2 <- sort(unlist(tmp1, use.names = FALSE))
+      if (ngroups > 1) {
+          tmp <- order(unlist(tmp1, use.names = FALSE))
+          fit_resid_md <- fit_resid_md[tmp]
+        }
       names(fit_resid_md) <- tmp2
       fit_resid_md_ordered <- order(fit_resid_md, decreasing = TRUE, na.last = NA)
       fit_resid_md_ordered <- fit_resid_md_ordered[!is.na(fit_resid_md_ordered)]
       fit_resid_md_selected <- fit_resid_md_ordered[seq_len(resid_md_top)]
       fit_resid_md_selected <- fit_resid_md_selected[!is.na(fit_resid_md_selected)]
-      to_rerun <- case_ids[fit_resid_md_selected]
+      if (!is.null(case_id)) {
+          to_rerun <- case_ids[fit_resid_md_selected]
+        } else {
+          to_rerun <- fit_resid_md_selected
+        }
     }
-
+  # listwise:
+  #  to_rerun:
+  #    no case_id: The positions in the *listwise* dataset
+  #    case_id: The case id to rerun
   if (!is.null(case_id)) {
       case_ids <- to_rerun
       id_to_rerun <- match(to_rerun, case_id)
@@ -311,12 +329,15 @@ lavaan_rerun <- function(fit,
       id_to_rerun <- tmp[to_rerun]
     }
   fit_total_time <- lavaan::lavInspect(fit, "timing")$total
+  lav_case_idx <- sort(unlist(lavaan::lavInspect(fit, "case.idx",
+                    drop.list.single.group = FALSE),
+                    use.names = FALSE))
   if (rerun_method == "lavaan") {
-      rerun_i <- gen_fct_use_lavaan(fit)
+      rerun_i <- gen_fct_use_lavaan(fit, lav_case_idx = lav_case_idx)
     }
   if (rerun_method == "update") {
       environment(gen_fct_use_update) <- parent.frame()
-      rerun_i <- gen_fct_use_update(fit)
+      rerun_i <- gen_fct_use_update(fit, lav_case_idx = lav_case_idx)
     }
   rerun_test <- suppressWarnings(rerun_i(NULL))
   if (!isTRUE(all.equal(unclass(coef(fit)),
@@ -393,7 +414,8 @@ lavaan_rerun <- function(fit,
   out
 }
 
-gen_fct_use_lavaan <- function(fit) {
+gen_fct_use_lavaan <- function(fit,
+                               lav_case_idx) {
   slot_opt <- fit@Options
   slot_pat <- data.frame(fit@ParTable)
   slot_pat$est <- NULL
@@ -403,16 +425,21 @@ gen_fct_use_lavaan <- function(fit) {
   ngroups <- lavaan::lavInspect(fit, "ngroups")
   if (ngroups > 1) {
       gp_var <- lavaan::lavInspect(fit, "group")
+      gp_label <- lavaan::lavInspect(fit, "group.label")
+      slot_opt$group.label <- gp_label
       out <- function(i = NULL) {
           if (is.null(i)) {
               return(lavaan::lavaan(data = data_full,
                                     model = slot_pat,
                                     group = gp_var,
+                                    group.label = gp_label,
                                     slotOptions = slot_opt))
             } else {
-              return(lavaan::lavaan(data = data_full[-i, ],
+              i1 <- match(i, lav_case_idx)
+              return(lavaan::lavaan(data = data_full[-i1, ],
                                     model = slot_pat,
                                     group = gp_var,
+                                    group.label = gp_label,
                                     slotOptions = slot_opt))
             }
         }
@@ -423,7 +450,8 @@ gen_fct_use_lavaan <- function(fit) {
                                     model = slot_pat,
                                     slotOptions = slot_opt))
             } else {
-              return(lavaan::lavaan(data = data_full[-i, ],
+              i1 <- match(i, lav_case_idx)
+              return(lavaan::lavaan(data = data_full[-i1, ],
                                     model = slot_pat,
                                     slotOptions = slot_opt))
             }

diff --git a/README.md b/README.md
@@ -9,7 +9,7 @@
 [![R-CMD-check](https://github.com/sfcheung/semfindr/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/sfcheung/semfindr/actions/workflows/R-CMD-check.yaml)
 <!-- badges: end -->
 
-(Version 0.1.6.1, updated on 2024-01-03, [release history](https://sfcheung.github.io/semfindr/news/index.html))
+(Version 0.1.6.2, updated on 2024-03-15, [release history](https://sfcheung.github.io/semfindr/news/index.html))
 
 # semfindr: Finding influential cases in SEM <img src="man/figures/logo.png" align="right" height="150" />
 

diff --git a/tests/testthat/test-lavaan_rerun.R b/tests/testthat/test-lavaan_rerun.R
@@ -4,7 +4,7 @@ library(semfindr)
 
 #context("Test lavaan_rerun")
 
-mod <- 
+mod <-
 '
 iv1 ~~ iv2
 m1 ~ iv1 + iv2
@@ -26,3 +26,16 @@ test_that("Compare parameter estimates of omitting an arbitrary case", {
       )
   })
 
+datm <- dat[1:20, ]
+datm[1, 2] <- datm[2, 3] <- datm[3, 4] <- datm[4, ] <- NA
+fitm <- lavaan::sem(mod, datm)
+
+test_that("Works for missing data", {
+    expect_no_error(rerunm_out <- lavaan_rerun(fitm))
+    fitm_10 <- lavaan::sem(mod, datm[-10, ])
+    rerunm_10 <- rerunm_out$rerun[["10"]]
+    expect_equal(ignore_attr = TRUE,
+        parameterEstimates(fitm_10), parameterEstimates(rerunm_10)
+      )
+  })
+
diff --git a/tests/testthat/test-lavaan_rerun_multi.R b/tests/testthat/test-lavaan_rerun_multi.R
@@ -33,3 +33,15 @@ test_that("Compare parameter estimates of omitting an arbitrary case", {
       )
   })
 
+datm <- dat0[1:25, ]
+datm[1, 2] <- datm[2, 3] <- datm[3, 4] <- datm[4, 1:4] <- NA
+fitm <- lavaan::sem(mod, datm, group = "gp")
+
+test_that("Works for missing data", {
+    expect_no_error(rerunm_out <- lavaan_rerun(fitm))
+    fitm_10 <- lavaan::sem(mod, datm[-10, ], group = "gp")
+    rerunm_10 <- rerunm_out$rerun[["10"]]
+    expect_equal(ignore_attr = TRUE,
+        parameterEstimates(fitm_10), parameterEstimates(rerunm_10)
+      )
+  })
diff --git a/tests/testthat/test-lavaan_rerun_multi_heywood.R b/tests/testthat/test-lavaan_rerun_multi_heywood.R
@@ -1,3 +1,4 @@
+skip_if_not_installed("MASS")
 library(testthat)
 library(lavaan)
 library(semfindr)
@@ -40,3 +41,55 @@ test_that("Warnings", {
                             what = "simpleWarning")),
                  0)
   })
+
+# With Listwise
+
+n <- 50
+sigma <- matrix(.3, 3, 3)
+diag(sigma) <- 1
+sigma <- sigma * (n - 1) / n
+set.seed(12345)
+dat0 <- MASS::mvrnorm(n, rep(0, 3), sigma, empirical = TRUE)
+dat0 <- as.data.frame(dat0)
+colnames(dat0) <- paste0("x", 1:3)
+cov(dat0)
+
+dat1 <- dat0
+dat1[1, 2] <- dat1[2, 3] <- dat1[4, ] <- NA
+dat1[6, 1] <- -3
+dat1[6, 2] <-  3
+set.seed(856041)
+dat1$gp <- sample(c("gp2", "gp1"), size = nrow(dat1), replace = TRUE)
+cov(dat1[dat1$gp == "gp1", -4], use = "complete.obs")
+cov(dat1[dat1$gp == "gp2", -4], use = "complete.obs")
+
+cfa_dat_heywood <- dat1
+
+mod <-
+"
+f1 =~ x1 + x2 + x3
+"
+suppressWarnings(fit <- lavaan::cfa(mod, cfa_dat_heywood, group = "gp"))
+
+attr(lavaan_rerun_check(fit), "info")
+
+test_that("Reject inadmissible solution", {
+    expect_error(lavaan_rerun(fit))
+  })
+
+tmp <- c(3, 6, 5)
+fit_rerun <- lavaan_rerun(fit, to_rerun = tmp, allow_inadmissible = TRUE)
+suppressWarnings(fit1 <- lavaan::cfa(mod, cfa_dat_heywood[-3, ], group = "gp", se = "none"))
+suppressWarnings(fit2 <- lavaan::cfa(mod, cfa_dat_heywood[-6, ], group = "gp", se = "none"))
+suppressWarnings(fit3 <- lavaan::cfa(mod, cfa_dat_heywood[-5, ], group = "gp", se = "none"))
+chk <- sum(sapply(list(fit1, fit2, fit3), function(x) {
+                    inherits(tryCatch(lavInspect(x, "post.check"),
+                            warning = function(w) w),
+                            "simpleWarning")
+                  }))
+
+test_that("Warnings", {
+    expect_equal(sum(sapply(fit_rerun$post_check, inherits,
+                            what = "simpleWarning")),
+                 chk)
+  })
diff --git a/tests/testthat/test-lavaan_rerun_multi_nonconvergence.R b/tests/testthat/test-lavaan_rerun_multi_nonconvergence.R
@@ -1,5 +1,6 @@
 skip_on_cran()
 # "Essential but may be machine dependent"
+skip_if_not_installed("MASS")
 
 library(testthat)
 library(lavaan)
@@ -75,3 +76,56 @@ test_that("fit_measures_change", {
     expect_equal(as.vector(!is.na(tmp)),
                  nc_check0)
   })
+
+
+# With Listwise
+
+n <- 50
+sigma <- matrix(.3, 3, 3)
+diag(sigma) <- 1
+sigma <- sigma * (n - 1) / n
+set.seed(12345)
+dat0 <- MASS::mvrnorm(n, rep(0, 3), sigma, empirical = TRUE)
+dat0 <- as.data.frame(dat0)
+colnames(dat0) <- paste0("x", 1:3)
+cov(dat0)
+
+dat1 <- dat0
+dat1[1, 2] <- dat1[2, 3] <- dat1[4, ] <- NA
+dat1[6, 1] <- -10
+dat1[6, 2] <-  10
+set.seed(856041)
+dat1$gp <- sample(c("gp2", "gp1"), size = nrow(dat1), replace = TRUE)
+cor(dat1[dat1$gp == "gp1", -4], use = "complete.obs")
+cor(dat1[dat1$gp == "gp2", -4], use = "complete.obs")
+
+cfa_dat_heywood <- dat1
+
+mod <-
+"
+f1 =~ x1 + x2 + x3
+"
+suppressWarnings(fit <- lavaan::cfa(mod, cfa_dat_heywood, group = "gp", control = list(iter.max = 10)))
+
+attr(lavaan_rerun_check(fit), "info")
+
+test_that("Reject inadmissible solution", {
+    expect_error(lavaan_rerun(fit))
+  })
+
+j <- c(3, 7, 5)
+fit_rerun <- lavaan_rerun(fit, to_rerun = j, skip_all_checks = TRUE)
+tmp <- sapply(j, function(x) {suppressWarnings(lavaan::cfa(mod, dat[-x, ], control = list(iter.max = 10), group = "gp"))})
+nc_check0 <- sapply(tmp, lavInspect, what = "converged")
+nc_check1 <- sum(nc_check0)
+pc_check0 <- suppressWarnings(sapply(tmp, lavInspect, what = "post.check"))
+pc_check1 <- sum(!pc_check0)
+
+test_that("Convergence", {
+    expect_equal(sum(sapply(fit_rerun$rerun, lavInspect,
+                            what = "converged")),
+                 nc_check1)
+    expect_equal(fit_rerun$converged,
+                 nc_check0,
+                 ignore_attr = TRUE)
+  })
diff --git a/tests/testthat/test-lavaan_rerun_multi_select_by_id.R b/tests/testthat/test-lavaan_rerun_multi_select_by_id.R
@@ -80,3 +80,28 @@ test_that("Check selected", {
     expect_equal(case_id_test[rerun_out$selected], case_id_to_rerun)
   })
 
+# Listwise
+
+dat0 <- dat[1:60, ]
+set.seed(856041)
+dat0$gp <- sample(c("gp2", "gp1"), size = nrow(dat0), replace = TRUE)
+
+dat0[1, 2] <- dat0[2, 3] <- dat0[3, 4] <- dat0[5, 1:4] <- NA
+fit0 <- lavaan::sem(mod, dat0, group = "gp")
+
+set.seed(4345)
+case_id_test <- paste0(sample(letters, nrow(dat0), replace = TRUE),
+                       sample(letters, nrow(dat0), replace = TRUE))
+case_id_to_rerun <- case_id_test[c(6, 4, 7)]
+rerun_out <- lavaan_rerun(fit0, case_id = case_id_test,
+                                to_rerun = case_id_to_rerun, parallel = FALSE)
+id_test <- which(case_id_test %in% case_id_to_rerun)[3]
+fit0_test <- lavaan::sem(mod, dat0[-id_test, ], group = "gp")
+
+rerun_test <- rerun_out$rerun[[case_id_test[id_test]]]
+
+test_that("Compare parameter estimates of omitting an arbitrary case", {
+    expect_equal(ignore_attr = TRUE,
+        parameterEstimates(fit0_test), parameterEstimates(rerun_test)
+      )
+  })
diff --git a/tests/testthat/test-lavaan_rerun_multi_select_by_md.R b/tests/testthat/test-lavaan_rerun_multi_select_by_md.R
@@ -44,6 +44,13 @@ test_that("Check selected", {
     expect_equal(rerun_out$selected, md_selected)
   })
 
+# For listwise
+
+fit0_case_ids <- unlist(lavInspect(fit0, "case.idx"), use.names = FALSE)
+test_that("Check case ids", {
+    expect_equal(sort(as.numeric(rownames(md_fit0))),
+                 sort(fit0_case_ids))
+  })
 
 # With Case ID