Merge branch 'main' into IK-edits

Author: haziqj

.github/workflows/pkgdown.yaml

@@ -1,12 +1,12 @@
 # Workflow derived from https://github.com/r-lib/actions/tree/v2/examples
 # Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
 on:
-  push:
-    branches: [main, master]
-  pull_request:
-    branches: [main, master]
-  release:
-    types: [published]
+  # push:
+  #   branches: [main, master]
+  # pull_request:
+  #   branches: [main, master]
+  # release:
+  #   types: [published]
   workflow_dispatch:
 
 name: 'pkgdown Quarto manuscript'

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: brlavaan
 Title: Bias-Reduced Latent Variable Analysis
-Version: 0.1.1.9005
+Version: 0.1.1.9007
 Authors@R: c(
     person("Haziq", "Jamil", , "haziq.jamil@gmail.com", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0003-3298-1010")),

R/10-sem_rbm_functions.R

@@ -286,14 +286,26 @@ fit_sem <- function(
     if (!is.null(lavargs$estimator.args$plugin_pen))
       plugin_pen <- lavargs$estimator.args$plugin_pen
   }
-  if ("information" %in% names(lavargs)) information <- lavargs$information
-  else information <- "observed"
-  if ("estimator" %in% names(lavargs)) estimator <- lavargs$estimator
-  else estimator <- "ML"
+  # Default estimator is ML and observed information for se
+  if ("information" %in% names(lavargs))
+    information <- lavargs$information
+  else
+    information <- "observed"
+  if ("se" %in% names(lavargs))
+    se <- lavargs$se
+  else
+    se <- "standard"
+  if ("estimator" %in% names(lavargs))
+    estimator <- lavargs$estimator
+  else
+    estimator <- "ML"
 
   # Validate arguments
   rbm <- validate_rbm(rbm)
-  if (estimator != "ML") cli::cli_abort("Bias reduction methods are currently only available for ML estimation.")
+  if (estimator != "ML")
+    cli::cli_abort("Bias reduction methods are currently only available for ML estimation.")
+  if (!(se %in% c("none", "standard", "robust.huber.white")))
+    cli::cli_abort("Bias reduction methods are currently only for standard and robust (Huber White) standard errors.")
 
   # Which method?
   is_ML     <- rbm == "none"
@@ -469,37 +481,57 @@ fit_sem <- function(
   }
 
   # Standard errors ------------------------------------------------------------
-  j <- information_matrix(
-    x = as.numeric(est),
-    lavmodel = lavmodel,
-    lavsamplestats = lavsamplestats,
-    lavdata = lavdata,
-    lavoptions = lavoptions,
-    kind = information
-  )
-  if (lavmodel@ceq.simple.only) {
-    K <- lavmodel@ceq.simple.K
-    j <- t(K) %*% j %*% K
-  }
-  if (check_mat(j)) {
-    sds <- rep(NA, length(est))
-    jinv <- NULL
-
-    # EXPERIMENTAL
-    if (TRUE) {
-      jinv <- try(solve(Matrix::nearPD(j)$mat), silent = !TRUE)
-      sds <- sqrt(diag(jinv))
-    }
-
+  sds <- rep(NA, length(est))
+  V <- NULL
+  if (se == "none") {
+    # Do nothing
   } else {
-    jinv <- try(solve(j), silent = !TRUE)
-    sds <- sqrt(diag(jinv))
+    j <- information_matrix(
+      x = as.numeric(est),
+      lavmodel = lavmodel,
+      lavsamplestats = lavsamplestats,
+      lavdata = lavdata,
+      lavoptions = lavoptions,
+      kind = information
+    )
+    if (lavmodel@ceq.simple.only) {
+      K <- lavmodel@ceq.simple.K
+      j <- t(K) %*% j %*% K
+    }
+    if (check_mat(j)) {
+      # Case: Information matrix is not positive definite
+      if (FALSE) {  ### EXPERIMENTAL ###
+        V <- try(solve(Matrix::nearPD(j)$mat), silent = !TRUE)
+        sds <- sqrt(diag(V))
+      }
+    } else {
+      jinv <- try(solve(j), silent = !TRUE)
+      if (se == "standard") {
+        V <- jinv
+      } else if (se == "robust.huber.white") {
+        e <- information_matrix(
+          x = as.numeric(est),
+          lavmodel = lavmodel,
+          lavsamplestats = lavsamplestats,
+          lavdata = lavdata,
+          lavoptions = lavoptions,
+          kind = "first.order"
+        )
+        if (lavmodel@ceq.simple.only) {
+          K <- lavmodel@ceq.simple.K
+          e <- t(K) %*% e %*% K
+        }
+        V <- jinv %*% e %*% jinv
+      }
+      sds <- sqrt(diag(V))
+    }
   }
 
   # Unpack estimators and se
   idx <- lavpartable$free[lavpartable$free > 0]
   est <- est[idx]
   sds <- sds[idx]
+  V   <- V[idx, idx, drop = FALSE]
   names(est) <- names(lavaan::coef(fit0))
 
   list(
@@ -509,9 +541,10 @@ fit_sem <- function(
     converged = res$convergence == 0L,
     scaled_grad = scaled_grad,
     optim = res,
-    vcov = jinv,
+    vcov = V,
     Sigma = Sigma,
     information = information,
+    se = se,
     lavfun = lavfun,
     estimator = estimator,
     rbm = rbm,

R/11-brlavaan.R

@@ -210,7 +210,7 @@ create_lav_from_fitsem <- function(
   #   else if (fit$estimator == "IBRMP") "IMP-BR ML"
   #   else if (fit$estimator == "EBRM") "EXP-BR ML"
   # Change vcov slot
-  fit0@vcov$se <- "standard"
+  fit0@vcov$se <- fit$se
   fit0@vcov$vcov <- fit$vcov
 
   # fit0@test <- fit_lav@test

R/12-plugins.R

@@ -13,9 +13,9 @@ NULL
 
 #' @rdname plugin-penalties
 #' @export
-pen_ridge <- function(x, call = FALSE, ...) {
+pen_ridge <- function(x, target = rep(0, length(x)), call = FALSE, ...) {
   if (isTRUE(call)) return("Ridge penalty")
-  sum(x ^ 2)
+  sum((x - target) ^ 2)
 }
 
 #' @rdname plugin-penalties

experiments/brsem/_prepDR.R

@@ -1,3 +1,7 @@
+# Purpose of this script is to prepare the raw simulation results from Dhaene &
+# Rosseel (2022) for the two factor model only (the growth model was re-run by
+# us). Need only run this script once.
+
 truth <-
   expand_grid(
     model = c("twofac", "growth"),
@@ -11,7 +15,6 @@ truth <-
   unnest_wider(out) |>
   unnest(c(param, truth))
 
-# Others
 # dr_file1 <- here::here("experiments/GCM_est_combined_final.RData")
 dr_file2 <- here::here("experiments/brsem/2FSEM_est_combined_final.RData")
 # if (!file.exists(dr_file1))

experiments/brsem/_seeds.R

@@ -4,7 +4,6 @@ here::i_am("experiments/obtain_seeds.R")
 simu_id <-
   expand_grid(
     dist = c("Normal", "Kurtosis", "Non-normal"),
-    # model = "twofac",
     rel = c(0.8, 0.5),
     n = c(15, 20, 50, 100, 1000)
   ) |>

experiments/brsem/_setup.R

@@ -2,7 +2,7 @@ library(brlavaan)
 library(tidyverse)
 library(furrr)
 theme_set(theme_bw())
-load(here::here("experiments/simu_id.RData"))
+load(here::here("experiments/brsem/simu_id.RData"))
 
 ncores <- future::availableCores() - 1
 future::plan(multisession, workers = ncores)
@@ -19,6 +19,8 @@ sim_fun <- function(
     lavsim = FALSE,
     whichsims = c("ML", "eRBM", "iRBM", "Ozenne", "REML"),
     bounds = "standard",
+    information = NULL,
+    se = NULL,
     data_scale = 1,
     seeds = NULL
   ) {
@@ -77,7 +79,9 @@ sim_fun <- function(
       maxgrad = FALSE,
       fn.scale = 1,
       bounds = bounds,
-      start = true_vals
+      start = true_vals,
+      information = information,
+      se = se
     )
 
     fit_list <- list()
@@ -98,9 +102,10 @@ sim_fun <- function(
     # D&R SIMS -----------------------------------------------------------------
     if (model == "growth") {
       fit_lav <- try(growth(mod, dat, start = true_vals, bounds = bounds,
-                            ceq.simple = TRUE))
+                            information = information, se = se, ceq.simple = TRUE))
     } else if (model == "twofac") {
-      fit_lav <- try(sem(mod, dat, bounds = bounds, meanstructure = TRUE))
+      fit_lav <- try(sem(mod, dat, bounds = bounds, meanstructure = TRUE,
+                         information = information, se = se))
     }
     if (!inherits(fit_lav, "try-error")) {
       out <- list()

experiments/brsem/prep_res.R

@@ -28,7 +28,9 @@ data_files <- paste0("experiments/brsem/", c(
   "simu_res_twofac_mp1.RData",
   "simu_res_twofac_mp2.RData",
   "simu_res_twofac_mp3.RData",
-  "simu_res_growth.RData"
+  "simu_res_growth.RData",
+  "simu_res_serobust_twofac.RData",
+  "simu_res_serobust_growth.RData"
 ))
 if (any(!file.exists(here::here(data_files)))) {
   cat("Downloading data files...\n")
@@ -57,6 +59,40 @@ simu_res_twofac[20:30] <- simu_res_twofac1[20:30]
 load(here::here("experiments/brsem/simu_res_growth.RData"))
 simu_res <- c(simu_res_twofac, simu_res_growth)
 
+# simu_res is a list of length 60
+# 2 models x 2 reliability x 3 distributions x 5 sample sizes = 60 simulations
+#
+# > glimpse(simu_res[[1]])
+# Rows: 13,979
+# Columns: 13
+# $ seed      <dbl> 1235, 1235, 1235, 1235, 1235, 1235, 1235, 1236, 1236, 1236, 1236, …
+# $ sim       <int> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 4, …
+# $ dist      <chr> "Normal", "Normal", "Normal", "Normal", "Normal", "Normal", "Norma…
+# $ model     <chr> "twofac", "twofac", "twofac", "twofac", "twofac", "twofac", "twofa…
+# $ rel       <chr> "0.8", "0.8", "0.8", "0.8", "0.8", "0.8", "0.8", "0.8", "0.8", "0.…
+# $ n         <dbl> 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15…
+# $ method    <chr> "ML", "eRBM", "iRBM", "lav", "Ozenne", "JB", "BB", "ML", "eRBM", "…
+# $ est       <named list> <0.74697154, 0.50594263, 1.40644613, 0.55439543, 0.23777413…
+# $ se        <named list> <0.14227840, 0.10401636, 0.33865785, 0.16550320, 0.15639856…
+# $ truth     <list> <0.7000, 0.6000, 0.7000, 0.6000, 0.2500, 0.2500, 0.1225, 0.0900, …
+# $ timing    <dbl> 0.042, 1.211, 16.191, 0.068, 0.225, 0.889, 275.691, 0.057, 2.293, …
+# $ converged <lgl> TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, …
+# $ Sigma_OK  <lgl> TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, TRUE, TRUE, TRUE, TRUE,…
+
+# Add new simulations about robust standard errors
+load(here::here("experiments/brsem/simu_res_serobust_twofac.RData"))
+load(here::here("experiments/brsem/simu_res_serobust_growth.RData"))
+simu_res_serobust <- c(simu_res_serobust_twofac, simu_res_serobust_growth)
+
+simu_res <-
+  map2(simu_res, simu_res_serobust, function(X, Y) {
+    left_join(
+      X,
+      select(Y, seed, sim, dist, model, rel, n, method, serob = se)
+    ) |>
+      select(seed:se, serob, everything())
+  })
+
 ## ----- Convergence statistics ------------------------------------------------
 res_ours_nested <-
   simu_res |>
@@ -73,7 +109,7 @@ res_ours_nested <-
         } else {
          return(TRUE)
         }
-      }, se, model)
+      }, serob, model)
   )
 
 res_conv <-
@@ -118,7 +154,7 @@ res <-
       grepl("[xy][0-9]~1", param) ~ "nu",
       TRUE ~ NA
     ),
-    across(c(est, truth, se), ~if_else(model == "growth" & type != "alpha", .x * 100, .x)),  # rescale back
+    across(c(est, truth, se, serob), ~if_else(model == "growth" & type != "alpha", .x * 100, .x)),  # rescale back
     dist = factor(dist, levels = c("Normal", "Kurtosis", "Non-normal")),
     rel = factor(rel, levels = c("0.8", "0.5"), labels = c("Rel = 0.8", "Rel = 0.5")),
     method = factor(
@@ -128,7 +164,8 @@ res <-
     ),
     bias = est - truth,
     relbias = bias / truth,
-    covered = truth <= est + qnorm(0.975) * se & truth >= est - qnorm(0.975) * se
+    covered = truth <= est + qnorm(0.975) * serob & truth >= est - qnorm(0.975) * serob,
+    covrdse = truth <= est + qnorm(0.975) * se & truth >= est - qnorm(0.975) * se
   )
 
 # Load D&R two factor sims
@@ -139,7 +176,9 @@ if (!file.exists(dr_data_file)) {
   load(dr_data_file)
 }
 
-# For two-factor model BB & JB, get from res_dr
+# For two-factor model BB & JB, get from res_dr (more stable and "nicer" results
+# compared to ours -- but the BB & JB for growth models were rerun by us, and
+# results similar to D&R 2022)
 res <-
   res |>
   filter(!(method %in% c("Bootstrap", "Jackknife") & model == "twofac")) |>
@@ -153,7 +192,8 @@ res <-
 plot_df <-
   res |>
   filter(method %in% c("ML", "eRBM", "iRBM")) |>
-  filter( converged, !is.na(se)) |>
+  filter(converged, !is.na(se)) |>
+  filter(param %in% c(twofacpars, growthpars)) |>
   # for each kind of model, filter bad standard errors
   filter(!(model == "twofac" & abs(se) > 5)) |>
   filter(!(model == "growth" & abs(se) > 500)) |>
@@ -167,7 +207,8 @@ plot_df <-
 plot_df50 <-
   res |>
   filter(method %in% c("ML", "eRBM", "iRBM")) |>
-  filter( converged, !is.na(se)) |>
+  filter(converged, !is.na(se)) |>
+  filter(param %in% c(twofacpars, growthpars)) |>
   # for each kind of model, filter bad standard errors
   filter(!(model == "twofac" & abs(se) > 5)) |>
   filter(!(model == "growth" & abs(se) > 500)) |>
@@ -244,12 +285,18 @@ create_summdf <- function(model, rel) {
 
 create_covrdf <- function(model) {
   i <- res$model == model
-  if (model == "growth") i <- i & res$method != "REML"
+  if (model == "twofac") i <- i & res$method != "REML"
 
   res |>
-    filter(dist != "Kurtosis", param %in% c(twofacpars, growthpars), i) |>
+    filter(i) |>
+    filter(dist != "Kurtosis", param %in% c(twofacpars, growthpars)) |>
     filter(!(model == "twofac" & abs(se) > 5)) |>
     filter(!(model == "growth" & abs(se) > 500)) |>
+    filter(method != "lav") |>
+    mutate(covered = case_when(
+      is.na(covered) ~ covrdse,  # plug in the JB and BB standard errors
+      TRUE ~ covered
+    )) |>
     summarise(
       covr = mean(covered, na.rm = TRUE),
       .by = c(dist:method, param)
@@ -412,3 +459,5 @@ save(twofacpars, growthpars, mycols, simu_id,
      bias_growth_80_df, bias_growth_50_df, covr_growth_df,
      tab_bias, tab_covr,
      file = here::here("experiments/brsem/results.RData"))
+
+save(twofacpars, growthpars, mycols, simu_id, res, file = here::here("experiments/brsem/full_res.RData"))