source("x.boot.r") #Include application
HS.model <- "
visual =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed =~ x7 + x8 + x9
x3 ~~ x5
x4 ~~ x7
x7 ~~ x8
"
# Create x.boot.object
set.seed(1)
fit <- sem(HS.model, HolzingerSwineford1939)
HS.boot <- x.boot(fit, R = 400, verbose = TRUE)
summary(HS.boot)
# Inspect some bootstrap draws
inspect(HS.boot, boot.num = c(1, 9), FUN = lavInspect, args = list("what" = "resid"))
# 95% BC CI (estimates)
summary(
x.boot.out(HS.boot, what = "est", ci.type = "bca.simple")
)
# 95% BC CI (standardized solution)
summary(
x.boot.out(HS.boot, what = "std.all", ci.type = "bca.simple")
)
# Bootstrap Multi-Group Model and autofill non-converged bootstrap draws
set.seed(1)
fit.groups <- sem(HS.model, HolzingerSwineford1939, group = "school")
HS.boot.groups <- x.boot(fit.groups, R = 300, verbose = TRUE, x.autofill = TRUE)
HS.boot.groups <- x.boot(HS.boot.groups, R = 100, verbose = TRUE, x.autofill = TRUE)
summary(HS.boot)
# 95% BC CI (estimates) with CI consistent p-value
summary(
x.boot.out(HS.boot.groups, what = "est", ci.type = "bca.simple", p.adjust = TRUE)
)
# 95% percentile CI for reliability coefficients (semTools)
summary(
x.boot.out(HS.boot.groups, myFUN = reliability, ci.type = "perc")
)
# R-squared differences for 2 group models
myRSquareDiff <- function(object) {
if(lavInspect(object, "ngroups") != 2) stop("Only for multi-group models with 2 groups")
res <- lavInspect(object, "rsquare")
c(res[[2]] - res[[1]])
}
summary(
x.boot.out(HS.boot.groups, myFUN = myRSquareDiff, ci.type = "bc", p.adjust = FALSE)
)
# HTMT for one or more groups
myHTMT <- function(object) {
res <- list()
groups = lavInspect(object, "ngroups")
for(g in 1 : groups) {
cov <- unclass(lavInspect(object, "sampstat"))[[g]]
if(typeof(cov) == "list") cov <- cov$cov
pt <- parameterTable(object)
model <- paste(apply(pt, 1, function(x) {
if(x["op"] == "=~" && as.numeric(x["group"]) == g) {
paste(x["lhs"], x["op"], if(as.numeric(x["free"]) == 0) paste(x["ustart"], "*", sep =""), x["rhs"], sep = "")
}
else ""
}), collapse ="\n")
res[[if(groups > 1) lavInspect(object, "group.label")[g] else g]] <- htmt(model, sample.cov = cov)
}
return(res)
}
summary(
x.boot.out(HS.boot, myFUN = myHTMT)
)
# The differences of all estimated parameters (could be further developed to mimic the Cheung/Lau invariance test without messy parametrization: https://journals.sagepub.com/doi/abs/10.1177/1094428111421987)
myBootDifEstimatedCoef <- function(object) {
if(lavInspect(object, "ngroups") != 2) stop("Only for multi-group models with 2 groups")
res <- list()
pe <- lavInspect(object, "x")
free <- lavInspect(object, "free")
for(m in names(pe[[1]])) {
free.g1 <- replace(replace(free[[1]][[m]], which(free[[1]][[m]] > 0), 1), which(free[[1]][[m]] == 0), NA)
free.g2 <- replace(replace(free[[2]][[m]], which(free[[2]][[m]] > 0), 1), which(free[[2]][[m]] == 0), NA)
g1 <- pe[[1]][[m]] * free.g1
g2 <- pe[[2]][[m]] * free.g2
res[[m]] <- g2 - g1
}
return(res)
}
summary(
x.boot.out(HS.boot.groups, myFUN = myBootDifEstimatedCoef, remove.na = TRUE)
)
# Bollen Stine p-value and an alternative calculation of the p-value using the adjusted p-value approach
set.seed(1)
fit.bs <- sem(HS.model, HolzingerSwineford1939)
HS.boot.bs <- x.boot(fit.bs, R = 400, verbose = TRUE, type = "bollen.stine")
myChisqDif <- function(object, chisq) {
fitMeasures(object, "chisq") - chisq
}
myBollenStine <- function(object) {
if(object$args$type != "bollen.stine") stop("Bootstrap type is not suitable.")
hist(object$c[object$t.names,])
get.ci <- function(t, alpha) {
res <- c(
as.numeric(quantile(t, alpha[1], type = 6, na.rm = TRUE)),
as.numeric(quantile(t, alpha[2], type = 6, na.rm = TRUE))
)
}
t <- object$c[object$t.names,]
p.bs <- sum(sapply(t, function(x) if(x > 0) 1 else 0)) / length(t)
p.adj <- if(p.bs < 1) p.bs else 2
if(!p.bs %in% c(0, 1)) {
p.adj = 0.0001
repeat {
p.adj <- p.adj + 0.0001
alpha <- 0.5 * c(p.adj, 2 - p.adj)
ci.adj <- get.ci(t, alpha)
if(sign(ci.adj[2]) == -1) break
}
}
out <- paste(
paste("bias ", round(object$c["bias", 1], 3), sep = ""),
paste("se ", round(object$c["se", 1], 3), sep = ""),
paste("ci ", "[", round(object$c["ci.lower", 1], 3), ", ", round(object$c["ci.upper", 1], 3), "]", sep = ""),
paste("p(adj) ", format(round(p.adj / 2, 3), nsmall = 3),sep = ""),
sep = "\n"
)
object$print <- paste(out, "\np(Bollen-Stine)", format(round(p.bs, 3), nsmall = 3), sep = " ")
return(object)
}
res <- x.boot.out(
myFUN = myChisqDif,
myArgs = list(chisq = fitMeasures(fit.bs, "chisq")),
myFUN.post = myBollenStine
)
summary(res)
lavaan community,
here is an additional example. It is about obtaining confidence intervals for indirect and total effects without having to create self-defined parameters (ab := a*b). This example is completely experimental. Please do not trust the results. I already know some bugs (see below). The source code (get.effects.r) is bad style and urgently needs to be reworked. I still want to share the example and am happy to get feedback. In particular, I'm interested if this solution could add value for you.
To test, you
need to save the attached file (get.effects.r) in your working directory. Here
are a few examples:
source("x.boot.r") # Include x.boot
source("get.effects.r") # Include get.effects
# Simulate a crazy model with two groups
pop.model.g1
<- "
f2 ~ 0.3 * f1
f3 ~ 0.5 * f1 + 0.2 * f2
f4 ~ 0.4 * f1 + 0.6 * f2 + 0.1 * f3
f5 ~ 0.1 * f2 + -0.2 * f4
f6 ~ 0.3 * f3 + 0.3 * f4 + -0.3 * f5
"
pop.model.g2 <- "
f2 ~ 0.1 * f1
f3 ~ 0.5 * f1 + 0.5 * f2
f4 ~ 0.3 * f1 + 0.2 * f2 + 0.4 * f3
f5 ~ 0.1 * f2 + -0.2 * f4
f6 ~ 0.3 * f3 + 0.3 * f4 + -0.7 * f5
"
s <- 1000
g <- c(rep("Group1", s), rep("Group2", s))
data <- rbind(
simulateData(pop.model.g1, sample.nobs = 1000),
simulateData(pop.model.g1, sample.nobs = 1000)
)
data <- cbind(data, g = c(rep("Group1", s), rep("Group2", s)))
# x.boot custom function
myEffects <- function(object, effects) {
sapply(effects, function(x, object) {
r <- unname(eval(parse(text = x)))
}, object)
}
# A simple mediation example (single group)
model <- "
f2 ~ f1
f3 ~ f1 + f2
"
set.seed(1)
fit <- sem(model, data[data$g == "Group1",])
boot <- x.boot(fit, R = 400, verbose = TRUE)
# x.boot results
x.boot.out(boot, myFUN = myEffects, myArgs = list(effects = get.effects(fit)), ci.type = "bca.simple")
# 2 groups mediation model
set.seed(1)
fit.groups <- sem(model, data, group = "g")
boot.groups <- x.boot(fit.groups, R = 400, verbose = TRUE)
# x.boot results
x.boot.out(boot.groups, myFUN = myEffects, myArgs = list(effects = get.effects(fit.groups)), ci.type = "bca.simple")
# Check against "native" lavaan solution
model.lavaan <- "
f2 ~ c(a1, a2) * f1
f3 ~ c(c1, c2) * f1 + c(b1, b2) * f2
G1.IND := a1 * b1
G1.TOT := a1 * b1 + c1
G2.IND := a2 * b2
G2.TOT := a2 * b2 + c2
"
fit.groups.native <- sem(model.lavaan, data, group = "g")
pe <- parameterEstimates(fit.groups.native)
pe[pe$op == ":=",]
# The crazy model (this could make the application interesting)
model.crazy <- "
f2 ~ f1
f3 ~ f1 + f2
f4 ~ f1 + f2 + f3
f5 ~ f2 + f4
f6 ~ f3 + f4 + f5"
set.seed(1)
fit.groups.crazy <- sem(model.crazy, data[data$g == "Group1",])
boot.groups.crazy <- x.boot(fit.groups.crazy, R = 400, verbose = TRUE)
#
x.boot results
x.boot.out(boot.groups.crazy, myFUN = myEffects, myArgs = list(effects = get.effects(fit.groups.crazy)), ci.type = "bca.simple")
The last example contains a lot of indirect effects. I did not check if the application could correctly identify all of them. With a simple mediation model it is easy to define new parameters. With the model above I perceive it more difficult and that is why I created this application.
Note: The application does not
work if cross group constraints are used. Strictly speaking, the algorithm already fails if any parameter labels are defined. This does not have to be the case. That would be solvable. Group
differences are not calculated. This could also be easily implemented. The meanstructure is not taken into account.
Before I completely rework the application
(if I have the time), I would appreciate your feedback.
Regards
Christian