ecmerkle / blavaan

An R package for Bayesian structural equation modeling
https://ecmerkle.github.io/blavaan
87 stars 23 forks source link

not happy about missing values #88

Closed jpritikin closed 2 months ago

jpritikin commented 3 months ago

Got the following error:

Error in if (sd(Y1) == 0 || sd(Y2) == 0) { : 
  missing value where TRUE/FALSE needed
Calls: blavaan ... lav_samplestats_step2 -> lav_bvord_cor_twostep_fit -> lav_bvord_init_cache

Same script as before:

library(rstan)
options(mc.cores = parallel::detectCores()/2)
library(blavaan)

future::plan("multicore")

NoYesLevels <- c('No','Not sure','Yes')
EuphoricLevels <- c('Euphoric', 'Mildly euphoric', 'No change', 'Mild discomfort', 'Severe discomfort')
ImpactLevels <- c('Helped','Slightly helped','No impact','Slightly hindered','Hindered')
AbsItemNames <- c('boredom','frustration','anxiety', 'abort','bodyAcute','bodyAfter','sleep')

ImportStudyItems <- function(df) {
  if (any(is.na(match(AbsItemNames, colnames(df))))) stop("Items are mising")

  df[['boredom']] <- ordered(df[['boredom']], levels=NoYesLevels)
  df[['frustration']] <- ordered(df[['frustration']], levels=NoYesLevels)
  df[['anxiety']] <- ordered(df[['anxiety']], levels=NoYesLevels)
  df[['abort']] <- ordered(df[['abort']], levels=NoYesLevels)
  df[['bodyAcute']] <- ordered(df[['bodyAcute']], levels=EuphoricLevels)
  df[['bodyAfter']] <- ordered(df[['bodyAfter']], levels=EuphoricLevels)
  df[['sleep']] <- ordered(df[['sleep']], levels=ImpactLevels)
  df
}

psiData <- ImportStudyItems(read.csv("psi.csv"))
pmData <- ImportStudyItems(read.csv("pm.csv"))
combinedData <- rbind(cbind(psiData, group="psi"),
                      cbind(pmData, group="pm"))
combinedData$group <- factor(combinedData$group)

spec <- '
phenom =~ boredomL*boredom + frustrationL*frustration + anxietyL*anxiety + abortL*abort + bodyAcuteL*bodyAcute + bodyAfterL*bodyAfter + sleepL*sleep

phenom ~ c(psiM, pmM)*1
boredom ~ boredomM*1
frustration ~ frustrationM*1
anxiety ~ anxietyM*1
abort ~ abortM*1
bodyAcute ~ bodyAcuteM*1
bodyAfter ~ bodyAfterM*1
sleep ~ sleepM*1

phenom ~~ 1*phenom
boredom ~~ 1*boredom
frustration ~~ 1*frustration
anxiety ~~ 1*anxiety
abort ~~ 1*abort
bodyAcute ~~ 1*bodyAcute
bodyAfter ~~ 1*bodyAfter
sleep ~~ 1*sleep

boredom | -.431*t1 + .431*t2
frustration | -.431*t1 + .431*t2
anxiety | -.431*t1 + .431*t2
abort | -.431*t1 + .431*t2
bodyAcute | -.842*t1 + -.253*t2 + .253*t3 + .842*t4
bodyAfter | -.842*t1 + -.253*t2 + .253*t3 + .842*t4
sleep | -.842*t1 + -.253*t2 + .253*t3 + .842*t4
'

fit <- blavaan(spec, data=combinedData,
               n.chains = 4,
               ordered = AbsItemNames,
               dp = dpriors(nu = "normal(0,1)", alpha="normal(0,1)", lambda="normal(0,1)"),
               group = "group",
               group.equal = c("loadings", "intercepts", "regressions"),
               std.lv=TRUE, allow.empty.cell = TRUE, cat.wls.w = FALSE, test = "none")
save(fit, file="main2.rda")

stanfit <- blavInspect(fit, 'mcobj')
divergent <- sapply(get_sampler_params(stanfit, inc_warmup=FALSE),
                    function(x) sum(x[,'divergent__']))
fitS <- summary(stanfit, probs=c())$summary
fitS <- fitS[!grepl('^log_lik|^ppp', rownames(fitS)),]
sims <- extract(stanfit, pars=rownames(fitS), permuted = FALSE, inc_warmup=FALSE)
bulkess <- apply(sims, 3, ess_bulk)
tailess <- apply(sims, 3, ess_tail)
print(paste("divergent", sum(divergent)))
print(max(apply(sims, 3, Rhat)))
print(min(c(bulkess, tailess)))

summary(fit)
blavInspect(fit, "hpd", level=.95)['pmM',]

pm.csv psi.csv

ecmerkle commented 3 months ago

Thanks, and it should be fixed. I always forget about na.rm = TRUE.