Nothing
## inspect blavaan object (wrapper around lavInspect with
## some additions)
blavTech <- function(blavobject, what, ...) {
blavInspect(blavobject, what, ...)
}
## use lavInspect everywhere we can:
blavInspect <- function(blavobject, what, ...) {
stopifnot(inherits(blavobject, "blavaan"))
what <- tolower(what)
dotdotdot <- list(...)
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"
## 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", "jagnames", "stannames",
"fscores", "lvs", "fsmeans", "lvmeans", "mcobj",
"rhat", "n_eff", "nchain", "nchains")
## blavwhats that don't require do.fit
blavnofit <- c("start", "starting.values", "inits", "n.chains", "cp", "dp",
"jagnames", "stannames", "nchain", "nchains")
## whats that are not handled
nowhats <- c("mi", "modindices", "modification.indices",
"wls.est", "wls.obs", "wls.v")
if(what %in% blavwhats){
if(!(what %in% blavnofit) & !blavobject@Options$do.fit){
stop(paste0("blavaan ERROR: ", what, " does not exist when do.fit = FALSE"))
}
if(jagtarget){
idx <- blavobject@ParTable$jagpnum
idx <- idx[!is.na(idx)]
} else {
idx <- blavobject@ParTable$stansumnum
if("pxnames" %in% names(blavobject@ParTable)){
drows <- grepl("^def", blavobject@ParTable$pxnames)
} else {
drows <- grepl("def", blavobject@ParTable$mat)
}
idx <- idx[blavobject@ParTable$free > 0 | drows]
}
labs <- lav_partable_labels(blavobject@ParTable, type = "free")
if(what %in% c("start", "starting.values", "inits")){
blavobject@external$inits
} else if(what %in% c("psrf", "ac.10", "neff", "rhat", "n_eff")){
if(jagtarget){
mcmcsumm <- blavobject@external$mcmcout$summaries
} else {
mcmcsumm <- rstan::summary(blavobject@external$mcmcout)$summary
}
if(what %in% c("psrf", "rhat")){
if(jagtarget){
OUT <- mcmcsumm[idx,'psrf']
} else {
OUT <- mcmcsumm[idx,'Rhat']
}
}else if(what == "ac.10"){
if(jagtarget){
OUT <- mcmcsumm[idx,'AC.10']
} else {
stop("blavaan ERROR: autocorrelation stat currently unavailable for Stan.")
}
} else {
if(jagtarget){
OUT <- mcmcsumm[idx,'SSeff']
} else {
OUT <- mcmcsumm[idx,'n_eff']
}
}
if(add.labels) names(OUT) <- labs
OUT
} else if(what %in% c("mcmc", "draws", "samples", "hpd")){
## add defined parameters to labels
pt <- blavobject@ParTable
pt$free[pt$op == ":="] <- max(pt$free, na.rm = TRUE) + 1:sum(pt$op == ":=")
labs <- lav_partable_labels(pt, type = "free")
draws <- make_mcmc(blavobject@external$mcmcout)
draws <- lapply(draws, function(x) mcmc(x[,idx]))
draws <- mcmc.list(draws)
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
}
draws
} else if(what == "mcobj"){
blavobject@external$mcmcout
} else if(what %in% c("fscores","lvs","fsmeans","lvmeans")){
if(jagtarget){
etas <- any(blavobject@external$mcmcout$monitor == "eta")
} else {
etas <- any(grepl("^eta", rownames(blavobject@external$stansumm)))
}
## how many lvs, excluding phantoms
lvmn <- lavInspect(blavobject, "mean.lv")
if(!inherits(lvmn, "list")){
lvmn <- list(lvmn)
}
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(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)
if(jagtarget){
drawcols <- grep("^eta\\[", colnames(draws[[1]]))
## remove phantoms
drawcols <- drawcols[1:(nlv * nsamp)]
} else {
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]))
## for target="stan" + missing, use @Data@Mp to reorder rows to correspond
## to original data
mis <- any(is.na(unlist(blavobject@Data@X)))
Mp <- blavobject@Data@Mp
if(blavobject@Options$target == "stan" & mis){
rorig <- sapply(Mp, function(x) unlist(x$case.idx))
empties <- sapply(Mp, function(x) x$empty.idx)
if(inherits(rorig, "list")){
## multiple groups
for(ii in length(rorig)){
rorig[[ii]] <- blavobject@Data@case.idx[[ii]][rorig[[ii]]]
}
rorig <- unlist(rorig)
}
cids <- Mp2dataidx(Mp, blavobject@Data@case.idx)
## reordering for lvs:
nfit <- sum(lavInspect(blavobject, 'nobs'))
rsamps <- rep(NA, nlv*nfit)
for(j in 1:nlv){
rsamps[((j-1)*nfit + 1):(j*nfit)] <- (j-1)*nfit + cids
}
for(j in 1:length(draws)){
draws[[j]][,rsamps] <- draws[[j]]
}
}
draws <- mcmc.list(draws)
if((what %in% c("lvmeans", "fsmeans")) | ("means" %in% dotdotdot)){
br <- TRUE
if(jagtarget){
summ <- blavobject@external$mcmcout$summaries
summname <- "Mean"
br <- FALSE
} else {
summ <- blavobject@external$stansumm
summname <- "mean"
}
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]
colnames(draws) <- names(lvmn[[level]])
if(blavobject@Options$target == "stan" & mis){
draws[rank(rorig),] <- draws
}
}
draws
} else if(what %in% c("n.chains", "nchain", "nchains")){
draws <- make_mcmc(blavobject@external$mcmcout)
length(draws)
} else if(what == "cp"){
blavobject@Options$cp
} else if(what == "dp"){
blavobject@Options$dp
} else if(what %in% c("postmode", "postmean", "postmedian")){
if(jagtarget){
mcmcsumm <- blavobject@external$mcmcout$summaries
} else {
mcmcsumm <- rstan::summary(blavobject@external$mcmcout)$summary
}
if(what == "postmean"){
if(jagtarget){
OUT <- mcmcsumm[idx,'Mean']
} else {
OUT <- mcmcsumm[idx,'mean']
}
}else if(what == "postmedian"){
if(jagtarget){
OUT <- mcmcsumm[idx,'Median']
} else {
OUT <- mcmcsumm[idx,'50%']
}
} else {
if(jagtarget){
OUT <- mcmcsumm[idx,'Mode']
} else {
stop("blavaan ERROR: Modes unavailable for Stan.")
}
}
if(add.labels) names(OUT) <- labs
OUT
} else if(what == "jagnames"){
if(!jagtarget) stop("blavaan ERROR: JAGS was not used for model estimation.")
OUT <- blavobject@ParTable$pxnames[blavobject@ParTable$free > 0]
OUT <- OUT[order(blavobject@ParTable$free[blavobject@ParTable$free > 0])]
if(add.labels) names(OUT) <- labs
OUT
} else if(what == "stannames"){
if(jagtarget) stop("blavaan ERROR: Stan was not used for model estimation.")
mcmcsumm <- rstan::summary(blavobject@external$mcmcout)$summary
OUT <- rownames(mcmcsumm)[idx]
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
lavargs <- c(dotdotdot, list(object = blavobject, what = what))
do.call("lavInspect", lavargs)
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.