clean up blavPredict and blavInspect for two-level lvs

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: blavaan
 Title: Bayesian Latent Variable Analysis
-Version: 0.4-9.1145
+Version: 0.4-9.1147
 Authors@R: c(person(given = "Edgar", family = "Merkle",
                     role = c("aut", "cre"),
                     email = "merklee@missouri.edu",

R/blav_object_inspect.R

@@ -16,6 +16,8 @@ blavInspect <- function(blavobject, what, ...) {
     dotNames <- names(dotdotdot)
     add.labels <- TRUE
     if(any(dotNames == "add.labels")) add.labels <- dotdotdot$add.labels
+    level <- 1L
+    if(any(dotNames == "level")) level <- dotdotdot$level
 
     jagtarget <- lavInspect(blavobject, "options")$target == "jags"
 
@@ -92,6 +94,7 @@ blavInspect <- function(blavobject, what, ...) {
             if(what == "hpd"){
                 pct <- .95
                 if("level" %in% dotNames) pct <- dotdotdot$level
+                if("prob" %in% dotNames) pct <- dotdotdot$prob
                 draws <- mcmc(do.call("rbind", draws))
                 draws <- HPDinterval(draws, pct)
                 if(add.labels) rownames(draws) <- labs
@@ -111,36 +114,36 @@ blavInspect <- function(blavobject, what, ...) {
             if(!inherits(lvmn, "list")){
                 lvmn <- list(lvmn)
             }
-            nlv <- length(lvmn[[1]])
-            nlv2 <- 0
-            if(length(lvmn) > 1) nlv2 <- length(lvmn[[2]])
+            if(level == 1L){
+              nlv <- length(lvmn[[1]])
+              nsamp <- sum(lavInspect(blavobject, "nobs"))
+            } else {
+              nlv <- length(lvmn[[2]])
+              nsamp <- unlist(lavInspect(blavobject, "nclusters"))
+            }
+
+            if(nlv == 0) stop("blavaan ERROR: no latent variables are at this level of the model")
+            if(!etas) stop("blavaan ERROR: factor scores not saved; set save.lvs=TRUE")            
 
-            if(nlv == 0) stop("blavaan ERROR: no latent variables are in the model")
-            if(!etas) stop("blavaan ERROR: factor scores not saved; set save.lvs=TRUE")
-
-            nsamp <- sum(lavInspect(blavobject, "nobs"))
-            nclus <- lavInspect(blavobject, "nclusters")
+            if(level == 2L & all(nsamp == 1)) stop("blavaan ERROR: level 2 was requested but the model is not multilevel")
 
             draws <- make_mcmc(blavobject@external$mcmcout, blavobject@external$stanlvs)
-            drawcols <- grep("^eta\\[", colnames(draws[[1]]))
-            drawcols2 <- grep("^eta_b", colnames(draws[[1]]))
 
             if(jagtarget){
                 ## remove phantoms
                 drawcols <- drawcols[1:(nlv * nsamp)]
             } else {
-                nfound <- length(drawcols)/nsamp
-                drawcols <- drawcols[as.numeric(matrix(1:length(drawcols),
-                                                       nsamp, nfound,
-                                                       byrow=TRUE)[,1:nlv])]
-
-                if(any(nclus > 1)){
-                  nfound2 <- length(drawcols2)/sum(nclus)
-                  drawcols2 <- drawcols2[as.numeric(matrix(1:length(drawcols2),
-                                                           sum(nclus), nfound2,
-                                                           byrow=TRUE)[,1:nlv2])]
-                  drawcols <- c(drawcols, drawcols2)
-                }
+              if(level == 1L){
+                drawcols <- grep("^eta\\[", colnames(draws[[1]]))
+              } else {
+                drawcols <- grep("^eta_b", colnames(draws[[1]]))
+                nsamp <- sum(nsamp)
+              }
+              nfound <- length(drawcols)/nsamp
+
+              drawcols <- drawcols[as.numeric(matrix(1:length(drawcols),
+                                                     nsamp, nfound,
+                                                     byrow=TRUE)[,1:nlv])]
             }
             draws <- lapply(draws, function(x) mcmc(x[,drawcols]))
 
@@ -183,21 +186,15 @@ blavInspect <- function(blavobject, what, ...) {
                     summ <- blavobject@external$stansumm
                     summname <- "mean"
                 }
-                mnrows <- grep("^eta\\[", rownames(summ))
-                mnrows2 <- grep("^eta_b", rownames(summ))
+                if(level == 1L){
+                  mnrows <- grep("^eta\\[", rownames(summ))
+                } else {
+                  mnrows <- grep("^eta_b", rownames(summ))
+                }
 
                 draws <- matrix(summ[mnrows,summname], nsamp,
                                 length(mnrows)/nsamp, byrow=br)[,1:nlv,drop=FALSE]
-                ## FIXME multiple groups?
-                colnames(draws) <- names(lvmn[[1]])
-
-                if(any(nclus > 1)){
-                  draws2 <- matrix(summ[mnrows2,summname], sum(nclus),
-                                   length(mnrows2)/sum(nclus), byrow=br)[,1:nlv2,drop=FALSE]
-                  colnames(draws2) <- names(lvmn[[2]])
-
-                  draws <- list(draws, draws2)
-                }
+                colnames(draws) <- names(lvmn[[level]])
 
                 if(blavobject@Options$target == "stan" & mis){
                     draws[rank(rorig),] <- draws

R/blav_predict.R

@@ -5,7 +5,7 @@ function(object, newdata = NULL) {
     blavPredict(object, newdata = newdata)
 })
 
-blavPredict <- function(object, newdata = NULL, type = "lv") {
+blavPredict <- function(object, newdata = NULL, type = "lv", level = 1L) {
 
   stopifnot(inherits(object, "blavaan"))
   blavmodel <- object@Model
@@ -30,6 +30,11 @@ blavPredict <- function(object, newdata = NULL, type = "lv") {
   stantarget <- lavopt$target == "stan"
 
   if(lavInspect(object, "categorical") & type == "ymis") stop("blavaan ERROR: ymis is not yet implemented for ordinal models.", call. = FALSE)
+
+  if(level == 2L){
+    if(all(unlist(lavInspect(object, "nclusters")) == 1)) stop("blavaan ERROR: level 2 was requested but this does not appear to be a 2-level model.", call. = FALSE)
+    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")
 
@@ -39,14 +44,25 @@ blavPredict <- function(object, newdata = NULL, type = "lv") {
   ## ypred: posterior predictive distribution of ovs conditioned on lv samples; mcmc list
   ## ymis: posterior predictive distribution of missing values conditioned on observed values; matrix
   if(type == "lv") {
-    FS <- do.call("rbind", blavInspect(object, 'lvs'))
+    FS <- do.call("rbind", blavInspect(object, 'lvs', level = level))
 
     ## N and latent variable names, to set dimensions
-    N <- sum(lavInspect(object, "ntotal"))
-    etas <- lavNames(object, "lv")
+    lvmn <- lavInspect(object, "mean.lv")
+    if(!inherits(lvmn, "list")){
+      lvmn <- list(lvmn)
+    }
+    if(level == 1L){
+      nlv <- length(lvmn[[1]])
+      N <- sum(lavInspect(object, "ntotal"))
+      etas <- names(lvmn[[1]])
+    } else {
+      nlv <- length(lvmn[[2]])
+      N <- sum(unlist(lavInspect(object, "nclusters")))
+      etas <- names(lvmn[[2]])
+    }
 
     out <- lapply(1:NROW(FS), function(i) {
-      rowmat <- matrix(FS[i,], N, length(etas))
+      rowmat <- matrix(FS[i,], N, nlv)
       colnames(rowmat) <- etas
       rowmat } )
   } else if(type == "lvmeans") {

man/blavInspect.Rd

@@ -14,7 +14,7 @@ blavTech(blavobject, what, ...)
 \arguments{
 \item{blavobject}{An object of class blavaan.}
 \item{what}{Character. What needs to be inspected/extracted? See Details for Bayes-specific options, and see \code{\link{lavaan}}'s \code{lavInspect()} for additional options. Note: the \code{what} argument is not case-sensitive (everything is converted to lower case.)}
-\item{...}{Default lavaan arguments supplied to \code{lavInspect()}; see \code{\link{lavaan}}.}
+\item{...}{lavaan arguments supplied to \code{lavInspect()}; see \code{\link{lavaan}}.}
 }
 \details{
 Below is a list of Bayesian-specific values for the \code{what}
@@ -37,11 +37,9 @@ documentation.
   \item{\code{"postmode"}:}{Estimated posterior mode of each free parameter.}
   \item{\code{"postmean"}:}{Estimated posterior mean of each free parameter.}
   \item{\code{"postmedian"}:}{Estimated posterior median of each free parameter.}
-  \item{\code{"lvs"}:}{An object of class \code{mcmc} containing latent
-  variable (factor score) draws.}
-  \item{\code{"lvmeans"}:}{A matrix of mean factor scores (rows are
-  observations, columns are variables).}
-  \item{\code{"hpd"}:}{HPD interval of each free parameter. In this case, an additional argument \code{level} can be supplied to specify a number in (0,1) reflecting the percentage of the interval.}
+  \item{\code{"lvs"}:}{An object of class \code{mcmc} containing latent variable (factor score) draws. In two-level models, use \code{level = 1} or \code{level = 2} to specify which factor scores you want.}
+  \item{\code{"lvmeans"}:}{A matrix of mean factor scores (rows are observations, columns are variables). Use the additional \code{level} argument in the same way.}
+  \item{\code{"hpd"}:}{HPD interval of each free parameter. In this case, the \code{prob} argument can be used to specify a number in (0,1) reflecting the desired percentage of the interval.}
 }
 }
 \seealso{
@@ -59,6 +57,6 @@ fit <- bcfa(HS.model, data = HolzingerSwineford1939,
 
 # extract information
 blavInspect(fit, "psrf")
-blavInspect(fit, "hpd", level = .9)
+blavInspect(fit, "hpd", prob = .9)
 }
 }

man/blavPredict.Rd

@@ -7,7 +7,7 @@
   types of model predictions, conditioned on observed data. This differs
   somewhat from \code{lavPredict()} in lavaan.}
 \usage{
-blavPredict(object, newdata = NULL, type = "lv")
+blavPredict(object, newdata = NULL, type = "lv", level = 1L)
 }
 \arguments{
 \item{object}{An object of class \code{\linkS4class{blavaan}}.}
@@ -20,6 +20,8 @@ the observed variables in the model are computed. If
 observed variables (including residual noise) are computed. If
 \code{"ymis"} or \code{"ovmis"}, model predicted values ("imputations")
 for the missing data are computed. See details for further information.}
+\item{level}{For \code{type = "lv"}, used to specify whether one desires
+the level 1 latent variables or level 2 latent variables.}
 }
 \details{
 The \code{predict()} function calls the \code{blavPredict()} function