version 0.2-1

cran · Sep 20, 2016 · 14fc7a8 · 14fc7a8
1 parent 8154f08
commit 14fc7a8
Show file tree

Hide file tree

Showing 19 changed files with 1,190 additions and 1,391 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: blavaan
 Title: Bayesian Latent Variable Analysis
-Version: 0.1-4
+Version: 0.2-1
 Authors@R: c(person(given = "Edgar", family = "Merkle",
                     role = c("aut", "cre"),
                     email = "merklee@missouri.edu"),
@@ -10,14 +10,15 @@ Authors@R: c(person(given = "Edgar", family = "Merkle",
             )
 Description: Fit a variety of Bayesian latent variable models, including confirmatory
    factor analysis, structural equation models, and latent growth curve models.
-Depends: R(>= 3.2.0), methods, runjags, lavaan(>= 0.5-20)
-Imports: stats, utils, graphics, MASS, MCMCpack, mnormt, nonnest2, loo
-Suggests: modeest, rjags, semTools
+Depends: R(>= 3.2.0), methods, runjags, lavaan(>= 0.5-21)
+Imports: stats, utils, graphics, MASS, MCMCpack, coda, mnormt,
+        nonnest2(>= 0.4-1), loo
+Suggests: modeest, rjags, semTools, parallel
 License: GPL (>= 2)
 NeedsCompilation: no
-Packaged: 2016-06-16 18:18:28 UTC; merkle
+Packaged: 2016-09-20 20:53:14 UTC; merkle
 Author: Edgar Merkle [aut, cre],
   Yves Rosseel [aut]
 Maintainer: Edgar Merkle <merklee@missouri.edu>
 Repository: CRAN
-Date/Publication: 2016-06-17 08:18:33
+Date/Publication: 2016-09-20 23:31:45
diff --git a/MD5 b/MD5
@@ -1,27 +1,28 @@
-8c9975ccb89f475b615fdbd23f4d7806 *DESCRIPTION
-ef4fcbecec4da9d82e4b40891475d3ec *NAMESPACE
-dace3b12f1eccef83ef5aaf1c089307a *NEWS
+7153a1d9df7b7ce88e412d811ffd7d4a *DESCRIPTION
+ac250dac20b3b310a1f747e71f6e92c0 *NAMESPACE
+81383663d5d5eb03a49838c04cb687e0 *NEWS
 be8bb39b5618c2843ad576dd22765133 *R/00class.R
 617fe666a4fef38b109e1f5d4a458330 *R/blav_fit.R
-b7c3f87a045f19fd75545be19034ae7c *R/blav_fit_measures.R
-30e1518cf04a4cea77ba5267948e4f94 *R/blav_object_methods.R
+e3dafa4dffdf2b13f6c8b12d3c560973 *R/blav_fit_measures.R
+39a4d1f861f4ba5b2f1631ce291dbe62 *R/blav_object_inspect.R
+3a255ce14aab087a682a23b8ee4ea6a0 *R/blav_object_methods.R
 2a75e183762a72bdae46f04ace1b6df3 *R/blav_test.R
-047ef73d006bf2eef9b1cd069d941dad *R/blav_utils.R
-e3bdc1abba0c93c752d2a314603fe294 *R/blavaan.R
+45453bb7fd344fc1065e7670b73ee9af *R/blav_utils.R
+3486d4c474be633aad72195f066b7293 *R/blavaan.R
 4d108f3aca9fc62aebd5d6382700ed40 *R/dpriors.R
 cfd89b2dcbe1cfb05c97fc18281cf08c *R/jags2r.R
-c14c36bd5e40148932cd553270d378f4 *R/lav_export_jags.R
-24dbc57700367f1acee296a1544cc6ff *R/margloglik.R
-01b5ad74843f03fd492ec4fa3ec0030a *R/postpred.R
-e3afb5a2ba0e78b879ff7653589737dc *R/set_inits.R
-78ccaafd2f60b8c6f6077a3e0e183ff8 *R/set_partable.R
-404a2cad6e31a78d6187eeb3a1d812ab *R/set_priors.R
+dcbdc1d438e18b9724d60fbe6e0cf31d *R/lav_export_jags.R
+d7e280f00cd92b5d60e688d885d49a3a *R/margloglik.R
+ce10d6846481213b7c27cd0f3c21b2e0 *R/postpred.R
+e616613f7b1d079d71b3df4ff95be06d *R/set_inits.R
+d99f910bc6a9d41c784dccb09c5c4171 *R/set_partable.R
+50a7a20ee6e854ef367b8d21f339f916 *R/set_priors.R
 ff6d505d1b13ce7fd4bd92949882f7ad *R/zzz.R
 76551b9ac5c7c9ff3c644a7b105b5742 *README
 da57c30fa52ef9a3b8563116a55b98dc *man/BF.Rd
-1a12761cb965497b78522609d837f850 *man/bcfa.Rd
-aec4d2ba5b4bbc04f65982e3a886c859 *man/bgrowth.Rd
+b504db7713efbe19cdbe15adad29d187 *man/bcfa.Rd
+9bb0e03fe7257d8056914df0f3eda1ba *man/bgrowth.Rd
 ca196de7665b078604b09bac468f523a *man/blav_internal.Rd
-93ae21624faf53a3c28cf3f6cad7e948 *man/blavaan.Rd
-1eeb38b23594547f11af815cffd40891 *man/bsem.Rd
+9634ebeb7252c0b82c2257c0ee67d728 *man/blavaan.Rd
+ac2aa5c341b95d6cd89dec943dd4abe2 *man/bsem.Rd
 c4c5c587307e68850481f82d60c2502c *man/dpriors.Rd
diff --git a/NAMESPACE b/NAMESPACE
@@ -1,5 +1,5 @@
 importFrom("methods", 
-           "is", "new")
+           "is", "new", "slot")
 
 importFrom("utils", 
            # "sessionInfo",
@@ -23,7 +23,7 @@ importFrom("lavaan",
            "lavaan", "logLik",
            "fitMeasures", "fitmeasures",
            "inspect", "lavInspect", "lavTech", "lavNames", 
-           "lavParseModelString", 
+           "lavParseModelString", "lavMatrixRepresentation",
            "lav_func_jacobian_complex", "lav_func_jacobian_simple", 
            "lav_partable_labels", "lavaanify", 
            "lav_model_get_parameters", "lav_model_implied", 
@@ -37,6 +37,11 @@ importFrom("runjags",
 importFrom("MCMCpack",
            "dwish")
 
+importFrom("coda",
+           "mcmc.list",
+           "mcmc",
+           "HPDinterval")
+
 importFrom("mnormt",
            "dmnorm")
 

diff --git a/NEWS b/NEWS
@@ -1,3 +1,18 @@
+Changes in Version 0.2-1
+    o Major update to internals: Model matrices/parameters now correspond to the Lisrel representation used in lavaan.
+
+    o General parameter equality constraints using '==' are now available (with one parameter on the lhs).
+
+    o New function blavInspect() for extracting various pieces of the MCMC run, including HPDs using an optional 'level' argument.
+
+    o JAGS syntax now uses the original observed variable names. It also assigns all prior/constraints to a single parameter vector, then defines model matrices based on this parameter vector.
+
+    o A list of user-defined initial values can be supplied via the inits argument.
+
+    o Posterior predictive computations are parallelized, if package parallel is installed.
+
+    o Improved timings for various parts of the model estimation.
+
 Changes in Version 0.1-4
     o New convergence="auto" option to run chains until convergence.
 

diff --git a/R/blav_fit_measures.R b/R/blav_fit_measures.R
@@ -16,7 +16,7 @@ blav_fit_measures <- function(object, fit.measures = "all",
                               baseline.model = NULL) {
 
     # has the model converged?
-    if(object@Fit@npar > 0L && !object@Fit@converged) {
+    if(object@Fit@npar > 0L && !object@optim$converged) {
         warning("blavaan WARNING: the chains may not have converged.")
     }
 

diff --git a/R/blav_object_inspect.R b/R/blav_object_inspect.R
@@ -0,0 +1,105 @@
+## inspect blavaan object (wrapper around lavInspect with
+## some additions)
+blavTech <- function(blavobject,
+                     what                   = "free",
+                     add.labels             = FALSE,
+                     add.class              = FALSE,
+                     list.by.group          = FALSE,  
+                     drop.list.single.group = FALSE,
+                     ...) {
+
+    blavInspect(blavobject = blavobject, what = what,
+                add.labels = add.labels, add.class = add.class,
+                list.by.group = list.by.group,
+                drop.list.single.group =  drop.list.single.group,
+                ...)
+}
+
+## use lavInspect everywhere we can:
+blavInspect <- function(blavobject,
+                        what                   = "free",
+                        add.labels             = TRUE,
+                        add.class              = TRUE,
+                        list.by.group          = TRUE,
+                        drop.list.single.group = TRUE,
+                        ...) {
+
+    stopifnot(inherits(blavobject, "blavaan"))
+    dotdotdot <- list(...)
+
+    ## only a single argument
+    if(length(what) > 1) {
+        stop("`what' arguments contains multiple arguments; only one is allowed")
+    }
+
+    ## be case insensitive
+    what <- tolower(what)
+
+    ## whats unique to blavaan
+    blavwhats <- c("start", "starting.values", "inits", "psrf",
+                   "ac.10", "neff", "mcmc", "draws", "samples",
+                   "n.chains", "cp", "dp", "postmode", "postmean",
+                   "postmedian", "hpd")
+
+    ## whats that are not handled (or modified handling)
+    nowhats <- c("mi", "modindices", "modification.indices",
+                 "wls.est", "wls.obs", "wls.v")
+
+    if(what %in% blavwhats){
+        idx <- blavobject@ParTable$jagpnum
+        idx <- idx[!is.na(idx)]
+        labs <- lav_partable_labels(blavobject@ParTable, type = "free")
+        if(what %in% c("start", "starting.values", "inits")){
+            blavobject@external$runjags$inits
+        } else if(what %in% c("psrf", "ac.10", "neff")){
+            if(what == "psrf"){
+                OUT <- blavobject@external$runjags$summaries[idx,'psrf']
+                ## blavobject@ParTable$psrf[!is.na(blavobject@ParTable$psrf)]
+            }else if(what == "ac.10"){
+                OUT <- blavobject@external$runjags$summaries[idx,'AC.10']
+            } else{
+                OUT <- blavobject@external$runjags$summaries[idx,'SSeff']
+            }
+            if(add.labels) names(OUT) <- labs
+            OUT
+        } else if(what %in% c("mcmc", "draws", "samples", "hpd")){
+            draws <- blavobject@external$runjags$mcmc
+            draws <- lapply(draws, function(x) x[,idx])
+            draws <- mcmc.list(draws)
+            if(what == "hpd"){
+                pct <- .95
+                if("level" %in% names(dotdotdot)) pct <- dotdotdot$level
+                draws <- mcmc(do.call("rbind", draws))
+                draws <- HPDinterval(draws, pct)
+                if(add.labels) rownames(draws) <- labs
+            }
+            draws
+        } else if(what == "n.chains"){
+            length(blavobject@external$runjags$mcmc)
+        } else if(what == "cp"){
+            blavobject@Options$cp
+        } else if(what == "dp"){
+            blavobject@Options$dp
+        } else if(what %in% c("postmode", "postmean", "postmedian")){
+            if(what == "postmean"){
+                OUT <- blavobject@external$runjags$summaries[idx,'Mean']
+            }else if(what == "postmedian"){
+                OUT <- blavobject@external$runjags$summaries[idx,'Median']
+            } else{
+                OUT <- blavobject@external$runjags$summaries[idx,'Mode']
+            }
+            if(add.labels) names(OUT) <- labs
+            OUT
+        }
+    } else if(what %in% nowhats){
+        stop(paste("blavaan ERROR: argument", what,
+                   "not available for Bayesian models."))
+    } else {
+        ## we can use lavInspect
+        lavInspect(lavobject = blavobject,
+                   what = what, add.labels = add.labels,
+                   add.class = add.class,
+                   list.by.group = list.by.group,
+                   drop.list.single.group = drop.list.single.group)
+    }
+}
diff --git a/R/blav_object_methods.R b/R/blav_object_methods.R
@@ -165,7 +165,6 @@ function(object, header       = TRUE,
                  standardized = FALSE,
                  rsquare      = FALSE,
                  std.nox      = FALSE,
-                 modindices   = FALSE,
                  psrf         = TRUE,
                  neff         = FALSE,
                  postmedian   = FALSE,
@@ -209,28 +208,16 @@ function(object, header       = TRUE,
         ## 95% HPD; FIXME display blanks for equality-constrained parameters
         ##                (like Std.Err column)
 
-        ## remove rho parameters from PE
-        rhos <- grep("rho", object@ParTable$jlabel[object@ParTable$op != "=="])
-        if(length(rhos) > 0) PE <- PE[-rhos,]
-        ## move rho priors to covariance rows
-        rhos <- grep("rho", object@external$runjags$origpt$jlabel)
-        covrhos <- grep("@rho", object@external$runjags$origpt$plabel)
-        object@ParTable$prior[covrhos] <- object@external$runjags$origpt$prior[rhos]
-        ## remove equality constraints from ParTable (rhos removed in blavaan())
-        eqc <- which(object@ParTable$op == "==")
-        if(length(eqc) > 0){
-            newpt <- lapply(object@ParTable, function(x) x[-eqc])
-        } else {
-            newpt <- object@ParTable
-        }
+        ## TODO put parameter priors in partable
+
+        newpt <- object@ParTable
 
         ## match jags names to partable, then partable to PE
-        ptentry <- match(rownames(object@external$runjags$HPD), newpt$jlabel) #object@ParTable$jlabel)
-        pte2 <- ptentry[!is.na(ptentry)]
+        pte2 <- which(!is.na(newpt$jagpnum))
         peentry <- match(with(newpt, paste(lhs[pte2], op[pte2], rhs[pte2], group[pte2], sep="")),
                          paste(PE$lhs, PE$op, PE$rhs, PE$group, sep=""))
-        PE$ci.lower[peentry] <- object@external$runjags$HPD[!is.na(ptentry),'Lower95']
-        PE$ci.upper[peentry] <- object@external$runjags$HPD[!is.na(ptentry),'Upper95']
+        PE$ci.lower[peentry] <- object@external$runjags$HPD[newpt$jagpnum[pte2],'Lower95']
+        PE$ci.upper[peentry] <- object@external$runjags$HPD[newpt$jagpnum[pte2],'Upper95']
 
         ## NB This is done so that we can remove fixed parameter hpd intervals without
         ##    making changes to lavaan's print.lavaan.parameterEstimates(). But maybe
@@ -245,13 +232,13 @@ function(object, header       = TRUE,
 
         ## FIXME defined parameters never get psrf + others;
         ## see line 200 of lav_print.R
-        if(psrf & class(object@external$runjags$psrf) != "character"){
+        if(psrf){
           PE$psrf <- rep(NA, nrow(PE))
-          PE$psrf[peentry] <- object@external$runjags$psrf$psrf[!is.na(ptentry),'Point est.']
+          PE$psrf[peentry] <- newpt$psrf[pte2]
         }
         if(neff){
           PE$neff <- rep(NA, nrow(PE))
-          PE$neff[peentry] <- object@external$runjags$summaries[!is.na(ptentry),'SSeff']
+          PE$neff[peentry] <- object@external$runjags$summaries[newpt$jagpnum[pte2],'SSeff']
         }
         if(priors){
           PE$prior <- rep(NA, nrow(PE))
@@ -260,11 +247,11 @@ function(object, header       = TRUE,
         }
         if(postmedian){
           PE$Post.Med <- rep(NA, nrow(PE))
-          PE$Post.Med[peentry] <- object@external$runjags$summaries[!is.na(ptentry),'Median']
+          PE$Post.Med[peentry] <- object@external$runjags$summaries[newpt$jagpnum[pte2],'Median']
         }
         if(postmode){
           PE$Post.Mode <- rep(NA, nrow(PE))
-          PE$Post.Mode[peentry] <- object@external$runjags$summaries[!is.na(ptentry),'Mode']
+          PE$Post.Mode[peentry] <- object@external$runjags$summaries[newpt$jagpnum[pte2],'Mode']
           if(all(is.na(PE$Post.Mode))) warning("blavaan WARNING: Posterior modes require installation of the modeest package.")
         }
         if(bf){
@@ -281,24 +268,17 @@ function(object, header       = TRUE,
         }
         ## alternative names because this is not ML
         penames <- names(PE)
-        ## FIXME we need an est column for print.lavaan if we have constraints
-        names(PE)[penames == "est"] <- "Post.Mean"
+        ## This could be called "Post.Mean" except constraints
+        ## require "est"
+        #names(PE)[penames == "est"] <- "Post.Mean"
+        #PE$est <- PE$Post.Mean
         names(PE)[penames == "se"] <- "Post.SD"
         names(PE)[penames == "ci.lower"] <- "HPD.025"
         names(PE)[penames == "ci.upper"] <- "HPD.975"
         names(PE)[penames == "psrf"] <- "PSRF"
         print(PE, nd = nd)
 
     } # parameter estimates
-
-    # modification indices?
-    if(modindices) {
-        cat("Modification Indices:\n\n")
-        object@Options$estimator <- "ML"
-        object@Fit@test[[2]] <- NULL
-        print( modificationIndices(object, standardized=TRUE) )
-    }
-
 })
 
 
@@ -382,7 +362,7 @@ function(object, header       = TRUE,
 
 plot.blavaan <- function(x, pars, plot.type="trace", ...){
     # NB: arguments go to plot.runjags()
-    parnames <- rownames(x@external$runjags$summaries)[pars]
+    parnames <- x@ParTable$pxnames[match(pars, x@ParTable$free)]
     plot(x@external$runjags, plot.type=plot.type, vars=parnames, ...)
 }