version 0.5-4

cran · Apr 16, 2024 · 0335fd6 · 0335fd6
1 parent 3173e6a
commit 0335fd6
Show file tree

Hide file tree

Showing 26 changed files with 308 additions and 89 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: blavaan
 Title: Bayesian Latent Variable Analysis
-Version: 0.5-3
+Version: 0.5-4
 Authors@R: c(person(given = "Edgar", family = "Merkle",
                     role = c("aut", "cre"),
                     email = "merklee@missouri.edu",
@@ -51,7 +51,7 @@ URL: https://ecmerkle.github.io/blavaan/,
 BugReports: https://github.com/ecmerkle/blavaan/issues
 Additional_repositories: https://mc-stan.org/r-packages/
 Config/Needs/website: brms
-Packaged: 2024-01-19 18:35:45 UTC; merkle
+Packaged: 2024-04-15 02:32:25 UTC; merkle
 Author: Edgar Merkle [aut, cre] (<https://orcid.org/0000-0001-7158-0653>),
   Yves Rosseel [aut],
   Ben Goodrich [aut],
@@ -66,4 +66,4 @@ Author: Edgar Merkle [aut, cre] (<https://orcid.org/0000-0001-7158-0653>),
   Matthew Emery [ctb] (loo moment_match)
 Maintainer: Edgar Merkle <merklee@missouri.edu>
 Repository: CRAN
-Date/Publication: 2024-01-19 22:50:02 UTC
+Date/Publication: 2024-04-15 07:00:02 UTC
diff --git a/MD5 b/MD5
@@ -1,6 +1,6 @@
-461286f25d5434597cbe125092fb9b9a *DESCRIPTION
-f167729dc2e8d93a9a96d4bfd9b8ce20 *NAMESPACE
-93527bada0f4719eb997384579482769 *NEWS.md
+afd24c31d091d76c7a32993334cab6fe *DESCRIPTION
+552d10d059a11a1db946204834266c5f *NAMESPACE
+324987a58dd135cc53b6281fc6e20278 *NEWS.md
 be8bb39b5618c2843ad576dd22765133 *R/00class.R
 00f32ef632f71eb006f1f4737b396334 *R/blav_adapt_quad.R
 58d31abedadd6c0ff3c7ce03c6c2add4 *R/blav_compare.R
@@ -9,22 +9,22 @@ d509e00f26effa24a635e682e31a3708 *R/blav_fit.R
 69048d7cd1b2f0c83c405613189ff7cb *R/blav_fit_measures.R
 850f341db92db5ea4df2e4fcc892e8ed *R/blav_model_loglik.R
 1c07731afacbe31ea0a3bfdfa9c65b39 *R/blav_object_inspect.R
-3d8c1f4ca672fceb64fa0c7e9b2904ac *R/blav_object_methods.R
-1b7771efa9c716d4985a2fef23b17c26 *R/blav_predict.R
+bfad69a019feafe389241a5ab011fe60 *R/blav_object_methods.R
+be1a830301b293518081a9d44edf4864 *R/blav_predict.R
 4bc1824825dce53725bc9cb00341b5bb *R/blav_test.R
 96a5dba3d202785751dbbd80f7daae29 *R/blav_utils.R
-92e79f9f69d001a781de8c7d5a0af203 *R/blavaan.R
+ddccc56649a152ef0f10480fef118e25 *R/blavaan.R
 af5158ba89256b0e98476b1249ab36b1 *R/ctr_bayes_fit.R
 4a8143e8807c05e2acc34291d22a21cf *R/ctr_ppmc.R
 e5d788f75302c12cb449c6aca64a5720 *R/dpriors.R
 bb0152fa9316b1c59a4aed14d85782d8 *R/jags2r.R
 d92d7cf8b732a8fe306c5dd9596077fe *R/lav_export_mcmc.R
 24ccfb7a3d8578e6233aafb5bb4b540d *R/lav_export_stanclassic.R
 660a0c5a7ad3357dfd4da654fc52e653 *R/lav_export_stancond.R
-c6984321d3115edfa05ed102f18a4dc4 *R/lav_export_stanmarg.R
-22b44d8a2f16f27486ce42234a5db874 *R/lvgqs.R
-b9239d9bbc8c302192771cc335f325c4 *R/margloglik.R
-dfc2e5bf745b6378fbbbdbd18bcde5fb *R/postpred.R
+7ecb1ecde7e42a9f136a4cd682e21996 *R/lav_export_stanmarg.R
+a0b1ada069510ce6d3b9e7200f748cc2 *R/lvgqs.R
+da501d5f5f67c59db09635125137dc34 *R/margloglik.R
+6d9c85e5ba7034302e95cdac76b773bc *R/postpred.R
 c90cf5bb58c1b5b50992aabde2deb6a1 *R/set_inits.R
 56a7c46edca960c89bd5a8ccdbf9860d *R/set_partable.R
 96c8a092fd316f412961dffc78111577 *R/set_priors.R
@@ -33,35 +33,36 @@ ee4e848a496cca28960d36aedcf3f9e6 *R/set_stanpars.R
 f0814cd9d865ee24ea57aeb60f4ec048 *R/stanmarg_data.R
 4e143ae1c959b2fa684926623565212f *R/stanmodels.R
 b13e5cca00122f286a8bad36b892c856 *R/zzz.R
-a104d24a545951ff405c9c49be05f480 *README.md
+9d451c9cd7f0ca8723b97a0f264a1876 *README.md
 1cc60bdd8610449d30184e0bb374bf46 *configure
 6bea7bdbd44ffa386677bfdec2d1984b *configure.win
 cbe0b566d817384e99d661ab83ea342c *inst/CITATION
 855c064dbf8deaff96c8dc8a87169109 *inst/include/stan_meta_header.hpp
 0ae76234cc0bacf4d108bff24e7d43a6 *inst/stan/include/license.stan
-3858f354ae3c6c87d53ef3c9cceaffb2 *inst/stan/stanmarg.stan
+e9bca3ea6cf3be0808778c58fe41b320 *inst/stan/stanmarg.stan
 1ae599465625c7c515da1e8ec5c3d7dd *inst/stanfuns/fill_lower.stan
 265b850f409b5196d59b3ef56e719db1 *inst/stanfuns/sem_lv.stan
 412c72bcc727c654bbf78a67bec0312b *inst/stanfuns/sem_lv_missing.stan
 126bc060cc1017c189560452d649fc52 *inst/stanfuns/sem_mean.stan
 b59addbe96ca232c7aa1e09c30c57db3 *inst/stanfuns/sem_mean_eta.stan
-3ed641c01614d9e0e9f7d78f77edb488 *inst/testdata/sysdata.rda
+0aa104c561efb07956c17c425754269b *inst/testdata/sysdata.rda
 28367418aa9e6dd6df963b288893af1c *inst/tinytest/tests.blavaan.R
 5c680bcac5c1d56f05fb7594816d93c9 *inst/tinytest/tests.blavaanobject-methods.R
-46f43bca5303e22d672041bf4f808f64 *man/bcfa.Rd
-aab75d72743709f5c3d288bfa52bf295 *man/bgrowth.Rd
-4bdffd1399d7dbbdc97171e53e38b868 *man/blavCompare.Rd
-aba7d0601546c15cf35a56c8735b7ab5 *man/blavFitIndices.Rd
-fc7ebac76807724da94ace4935f485e9 *man/blavInspect.Rd
-4305ca4f1c246abd13042e5a10f69047 *man/blavPredict.Rd
+dd93575f1e12eaa8e576a8d75dddd5e8 *man/bcfa.Rd
+7505365a89a8ec25a58b783a01d6fc1b *man/bgrowth.Rd
+5e739e724f9faccd4893f9c67a9eef21 *man/blavCompare.Rd
+8df1bb611e3a3174d25ded3391e014f2 *man/blavFitIndices.Rd
+11b149e7ab172b38fdc63aba76299f81 *man/blavInspect.Rd
+229eac402835d3ae27c4c343d4d34f35 *man/blavPredict.Rd
 fe83c6ddce1d9dbae013f417bbc9eb24 *man/blav_internal.Rd
 3a3043db5c4b4d06331d360a69374996 *man/blavaan-class.Rd
-c0617d10963b1a9aea38c2df935b956d *man/blavaan.Rd
-c04399107ea5cd13d3f8513d6bb86067 *man/bsem.Rd
+dafea8149a9cf11d37de9c93156f334d *man/blavaan.Rd
+687f9fbe523258827cac1c03b0c17ee5 *man/bsem.Rd
 05369042b3e04ef3863b7106d831e1ee *man/dpriors.Rd
-f13a1b152b2cb7648dd82f034b6b5068 *man/plot.blavaan.Rd
-3ccbc39b7ccd3d706b9b13ae044b326e *man/ppmc.Rd
-2a6303ebd30f409294e6b126c8c8049d *man/standardizedPosterior.Rd
+f0ea57da2e88f0a943988f50a33008ce *man/plot.blavaan.Rd
+af4913d0e1ef8cebc84f4429b5fdc992 *man/ppmc.Rd
+5e6d263485ade45338eceb1aa25d1b19 *man/sampleData.Rd
+2d3670863bf5d4f2a4a72ef4c4625a34 *man/standardizedPosterior.Rd
 564c40d0e86bb1e58d039ff240e97b9f *src/Makevars
 f3a200f96ee84d7f928a73432f70dc8e *src/RcppExports.cpp
 f0610a919082b019534974e8b1a50b93 *src/stanExports_stanmarg.cc

diff --git a/NAMESPACE b/NAMESPACE
@@ -36,7 +36,7 @@ importFrom("lavaan",
            "lav_model_set_parameters", "lav_model_vcov_se",
            "lav_partable_attributes",
            "modificationIndices", "parTable", "parameterEstimates",
-           "lavPredict", "standardizedSolution")
+           "lavPredict", "standardizedSolution", "lav_data_update")
 
 importFrom("coda",
            "mcmc.list",
@@ -94,7 +94,7 @@ import(rstantools)
 
 export("blavaan", "bcfa", "bsem", "bgrowth", "dpriors", "BF", "blavCompare",
        "blavTech", "blavInspect", "blavFitIndices", "labelfun", "standardizedposterior",
-       "standardizedPosterior", "ppmc", "blavPredict")
+       "standardizedPosterior", "ppmc", "blavPredict", "sampleData")
 
 exportClasses("blavaan", "blavPPMC", "blavFitIndices")
 

diff --git a/NEWS.md b/NEWS.md
@@ -1,4 +1,14 @@
 
+# Version 0.5-4
+## New features
+* New function sampleData() for generating data from a blavaan object.
+
+* Functionality for the newdata argument in blavPredict(), which can generate lv (and other) predictions for new data from a model that has already been estimated (continuous data only, target = "stan" only).
+
+* Refinements to two-level models (handle missingness via listwise deletion; lvs only at one level; improved messaging)
+
+* Bugs from 0.5-3 are fixed.
+
 # Version 0.5-3
 ## New features
 * Functionality to find unrestricted blocks of the model's psi matrix (lv covariance matrix). lkj priors are assigned to these unrestricted blocks, improving the positive definite issue described in the "Opaque priors" paper.
@@ -13,6 +23,8 @@
 
 * Bugs from 0.5-2 are fixed.
 
+## Bugs/glitches discovered after the release:
+* Some models with exogenous covariates, fixed.x=TRUE, and missing data fail to converge and yield implausible parameter values (reported by DeAnne Hunter).
 
 # Version 0.5-2
 ## New features

diff --git a/R/blav_object_methods.R b/R/blav_object_methods.R
@@ -70,7 +70,7 @@ function(object, header       = TRUE,
                  ci           = TRUE,
                  standardized = FALSE,
                  rsquare      = FALSE,
-                 std.nox      = FALSE,
+                 std.nox      = FALSE, #TODO: remove deprecated argument in early 2025
                  psrf         = TRUE,
                  neff         = FALSE,
                  postmedian   = FALSE,
@@ -79,7 +79,7 @@ function(object, header       = TRUE,
                  bf           = FALSE,
                  nd = 3L) {
 
-    if(std.nox) standardized <- TRUE
+    #TODO: remove (deprecated):  if(std.nox) standardized <- TRUE
 
     # print the 'short' summary
     if(header) {
@@ -122,9 +122,10 @@ function(object, header       = TRUE,
 
         if(!("group" %in% names(PE))) PE$group <- 1
         if(!("level" %in% names(PE))) PE$level <- "within"
-        if(standardized && std.nox) {
-            PE$std.all <- PE$std.nox
-        }
+        #TODO: remove deprecated argument in early 2025
+        # if(standardized && std.nox) {
+        #     PE$std.all <- PE$std.nox
+        # }
         PE$group[PE$group == 0] <- 1
 
         if("target" %in% names(object@call)){

diff --git a/R/blav_predict.R b/R/blav_predict.R
@@ -36,7 +36,12 @@ blavPredict <- function(object, newdata = NULL, type = "lv", level = 1L) {
     if(type %in% c("yhat", "ypred", "ymis")) stop("blavaan ERROR: option", type, "is not yet implemented for two-level models.", call. = FALSE)
   }
 
-  if(!is.null(newdata)) stop("blavaan ERROR: posterior predictions for newdata are not currently supported")
+  if(!is.null(newdata)) {
+    if(!stantarget) stop("blavaan ERROR: newdata is currently only available for target='stan'")
+    if(lavInspect(object, "categorical")) stop("blavaan ERROR: newdata is not yet available for ordinal data.")
+    object <- blav_fill_newdata(object, newdata)
+  }
+
 
   ## lv: posterior dist of lvs (use blavInspect functionality); matrix frame
   ## lvmeans: use blavInspect functionality; matrix
@@ -120,3 +125,44 @@ blavPredict <- function(object, newdata = NULL, type = "lv", level = 1L) {
 
   out
 }
+
+## fill blavaan object with newdata, then sample lvs given already-sampled parameters
+blav_fill_newdata <- function(object, newdat) {
+
+  lavd <- getFromNamespace("lavData", "lavaan")
+  olddata <- object@Data
+  OV <- olddata@ov
+  object@Data <- lavd(data = newdat,
+                      group = olddata@group,
+                      ov.names = olddata@ov.names,
+                      ov.names.x = olddata@ov.names.x,
+                      ordered = OV$names[ OV$type == "ordered" ],
+                      lavoptions = object@Options, allow.single.case = TRUE)
+
+  ## Stan-formatted newdata
+  l2s <- lav2stanmarg(object, dp = blavInspect(object, 'options')$dp,
+                      n.chains = blavInspect(object, 'nchains'), inits = "simple")
+  l2slev2 <- lav2stanmarg(object, dp = blavInspect(object, 'options')$dp,
+                          n.chains = blavInspect(object, 'nchains'),
+                          inits = "simple", level = 2, indat = l2s$dat)
+  l2s$dat <- c(l2s$dat, l2slev2$dat)
+  l2s$dat <- l2s$dat[!duplicated(names(l2s$dat))]
+  l2s$free2 <- c(l2s$free2, l2slev2$free2)
+  l2s$lavpartable <- rbind(l2s$lavpartable, l2slev2$lavpartable)
+  l2s$wigpris <- c(l2s$wigpris, l2slev2$wigpris)
+  l2s$init <- lapply(1:length(l2s$init), function(i) c(l2s$init[[i]], l2slev2$init[[i]]))
+  ldargs <- c(l2s$dat, list(lavpartable = l2s$lavpartable, dumlv = l2s$dumlv, dumlv_c = l2slev2$dumlv,
+                            save_lvs = TRUE, do_test = FALSE))
+  smd <- do.call("stanmarg_data", ldargs)
+
+  newlvs <- samp_lvs(object@external$mcmcout, object@Model, object@ParTable, smd, eeta = NULL, categorical = FALSE)
+  lvsumm <- as.matrix(rstan::monitor(newlvs, print=FALSE))
+  cmatch <- match(colnames(object@external$stansumm), colnames(lvsumm))
+  stansumm <- object@external$stansumm
+  lvcols <- grep("^eta", rownames(stansumm))
+  if (length(lvcols) > 0) stansumm <- stansumm[-lvcols, ]
+  object@external$stansumm <- rbind(stansumm, lvsumm[,cmatch])
+  object@external$stanlvs <- newlvs
+
+  object
+}
diff --git a/R/blavaan.R b/R/blavaan.R
@@ -321,7 +321,15 @@ blavaan <- function(...,  # default lavaan arguments
     dotdotdot$do.fit <- TRUE; dotdotdot$warn <- FALSE
     if(LAV@Data@data.type != "moment" && target == "stan"){
         ## if no missing, set missing = "listwise" to avoid meanstructure if possible
-        if(!any(is.na(unlist(lavInspect(LAV, 'data'))))) dotdotdot$missing <- "listwise"
+        if(!any(is.na(unlist(lavInspect(LAV, 'data'))))){
+            dotdotdot$missing <- "listwise"
+        } else {
+            if("cluster" %in% dotNames) {
+                ## set missing = "listwise" for two-level models
+                dotdotdot$missing <- "listwise"
+                cat("blavaan NOTE: listwise deletion is in use! (currently the only missingness option for two-level models)\n\n")
+            }
+        }
     }
 
     # for initial values/some parameter setup:
@@ -345,8 +353,8 @@ blavaan <- function(...,  # default lavaan arguments
     # save.lvs in a model with no lvs
     if(save.lvs){
         clv <- lavInspect(LAV, 'cov.lv')
-        if(is.list(clv)) clv <- clv[[1]]
-        if(nrow(clv) == 0) warning("blavaan WARNING: save.lvs=TRUE, but there are no lvs in the model.", call. = FALSE)
+        if(!is.list(clv)) clv <- list(clv)
+        if(all(sapply(clv, nrow) == 0)) warning("blavaan WARNING: save.lvs=TRUE, but there are no lvs in the model.", call. = FALSE)
     }
 
     # turn warnings back on by default
@@ -402,7 +410,10 @@ blavaan <- function(...,  # default lavaan arguments
                      LAV@ParTable$rhs %in% LAV@ParTable$plabel[cov.pars]) &
                      LAV@ParTable$op == "=="))
     if(con.cov) LAV@Options$auto.cov.lv.x <- FALSE
-
+
+    # ensure group is in the parTable
+    if(!("group" %in% names(LAV@ParTable))) LAV@ParTable$group <- rep(1L, length(LAV@ParTable$lhs))
+
     # if std.lv, truncate the prior of each lv's first loading
     loadpt <- LAV@ParTable$op == "=~"
     lvs <- unique(LAV@ParTable$lhs[loadpt])
@@ -502,6 +513,9 @@ blavaan <- function(...,  # default lavaan arguments
     lavoptions$prisamp   <- prisamp
     lavoptions$target    <- target
     lavoptions$optim.method <- "mcmc"
+    lavoptions$burnin <- burnin
+    lavoptions$sample <- sample
+    lavoptions$n.chains <- n.chains
     if("llnsamp" %in% names(mcmcextra$data)){
         if(length(mcmcextra$data$llnsamp) > 1 ||
            (!inherits(mcmcextra$data$llnsamp, "numeric") &&
@@ -1104,13 +1118,15 @@ blavaan <- function(...,  # default lavaan arguments
         lavInspect(blavaan, "post.check")
     }
 
-    if( "psi" %in% lavpartable$mat &&
-        ( (target == "stan" && !l2s$blkpsi) ||
-          (target != "stan" && with(covres, !(diagpsi | fullpsi))) ) ) {
-      warning("blavaan WARNING: As specified, the psi covariance matrix is neither diagonal nor unrestricted, so the actual prior might differ from the stated prior. See\n https://arxiv.org/abs/2301.08667", call. = FALSE)
-    }
-    if( "theta" %in% lavpartable$mat && with(covres, !(diagthet | fullthet)) ) {
-      warning("blavaan WARNING: As specified, the theta covariance matrix is neither diagonal nor unrestricted, so the actual prior might differ from the stated prior. See\n https://arxiv.org/abs/2301.08667", call. = FALSE)
+    if(!lavoptions$.multilevel) { # because checkcovs() has not been adapted to it
+      if( "psi" %in% lavpartable$mat &&
+          ( (target == "stan" && !l2s$blkpsi) ||
+            (target != "stan" && with(covres, !(diagpsi | fullpsi))) ) ) {
+        warning("blavaan WARNING: As specified, the psi covariance matrix is neither diagonal nor unrestricted, so the actual prior might differ from the stated prior. See\n https://arxiv.org/abs/2301.08667", call. = FALSE)
+      }
+      if( "theta" %in% lavpartable$mat && with(covres, !(diagthet | fullthet)) ) {
+        warning("blavaan WARNING: As specified, the theta covariance matrix is neither diagonal nor unrestricted, so the actual prior might differ from the stated prior. See\n https://arxiv.org/abs/2301.08667", call. = FALSE)
+      }
     }
 
     if(jag.do.fit & lavoptions$warn & !prisamp & !usevb & !grepl("stan", target)){

diff --git a/R/lav_export_stanmarg.R b/R/lav_export_stanmarg.R
@@ -147,22 +147,24 @@ matattr <- function(free, est, constraint, mat, Ng, std.lv, wig, ...) {
       }
     } else if (lammat) {
       for (i in 1:length(free2)) {
-        for (j in 1:NCOL(free2[[i]])) {
-          col <- free[[i]][,j]
-          col2 <- free2[[i]][,j]
-          porg <- col[col != 0L]
-          parnums <- col2[col2 != 0L]
-          if (length(parnums) > 0) {
-            psign <- parnums[which.min(porg)] # sign constrain the first loading from user syntax
-            ## if equality constraint, sign must involve the
-            ## "free" parameter
-            if (wskel[psign,1] == 1L) {
-              psign <- wskel[psign,2]
-            } else {
-              psign <- psign - sum(wskel[1:(psign-1),1] == 1)
+        if (NCOL(free2[[i]]) > 0) {
+          for (j in 1:NCOL(free2[[i]])) {
+            col <- free[[i]][,j]
+            col2 <- free2[[i]][,j]
+            porg <- col[col != 0L]
+            parnums <- col2[col2 != 0L]
+            if (length(parnums) > 0) {
+              psign <- parnums[which.min(porg)] # sign constrain the first loading from user syntax
+              ## if equality constraint, sign must involve the
+              ## "free" parameter
+              if (wskel[psign,1] == 1L) {
+                psign <- wskel[psign,2]
+              } else {
+                psign <- psign - sum(wskel[1:(psign-1),1] == 1)
+              }
+              sign[parnums, 1] <- 1L
+              sign[parnums, 2] <- psign
             }
-            sign[parnums, 1] <- 1L
-            sign[parnums, 2] <- psign
           }
         }
       }

diff --git a/R/lvgqs.R b/R/lvgqs.R
@@ -320,7 +320,9 @@ samp_lvs_2lev <- function(mcobj, lavmodel, lavsamplestats, lavdata, lavpartable,
         standata$Ntot <- sum(standata$nclus[,2])
         standata$Nobs <- with(standata, rep(N_between + N_both, Np))
         standata$Obsvar <- with(standata, matrix(1:standata$Nobs[1], Np, N_between + N_both, byrow = TRUE))
-        tmpmat2[j,,] <- lvgqs(modmat2, standata, eeta[2*(1:standata$Ng)])
+        if (standata$m > 0) {
+          tmpmat2[j,,] <- lvgqs(modmat2, standata, eeta[2*(1:standata$Ng)])
+        }
 
         ## now level 1
         standata <- stanorig
@@ -331,7 +333,9 @@ samp_lvs_2lev <- function(mcobj, lavmodel, lavsamplestats, lavdata, lavpartable,
           if(!("beta" %in% names(modmat1[[g]]))) modmat1[[g]]$beta <- matrix(0, standata$m, standata$m)
         }
 
-        tmpmat[j,,] <- lvgqs(modmat1, standata, eeta[2 * (1:standata$Ng) - 1])
+        if (standata$m > 0) {
+          tmpmat[j,,] <- lvgqs(modmat1, standata, eeta[2 * (1:standata$Ng) - 1])
+        }
       }
       list(tmpmat, tmpmat2)})
   if(!debug) {