Bootstrap estimates, standardized solution, R squares, HTMT, group differences, consistent p-values, and other coefficients with a single bootstrap run
Christian Arnold
I occasionally read questions about various problems and challenges with the bootstrap procedure. I have taken the liberty of developing a small application that may be able to contribute to eliminating a problem or two. It is important for me to clarify that this is only (only!) a prototype. There is no guarantee that the application will produce correct results. Nevertheless, I am happy to receive any feedback. Criticism. Praise. Ideas. Everything is very welcome. I have attached the application. If you want to test 0it, you need to save it to your working directory and include it in your R script. The procedure is as follows:
(1) Bootstrap with x.boot
x.boot accepts all arguments that bootstrapLavaan processes (except parallel - sorry I'm a WIN user and couldn't test this). Unlike bootstrapLavaan, all bootstrapped lavaan objects are saved. This probably has some advantages (see below), but is memory intensive. x.boot additionally accepts the argument x.autofill = TRUE (default is FALSE to mimic the lavaan behavior). x.autofill = TRUE bootstraps until all non-converged bootstrap draws are replaced by converged draws. x.boot also accepts the argument x.disc. Here you can specify a path to the local file system. Please do not use the working directory here and use this feature only if you request a lot of bootstrap draws and the RAM of your machine is not sufficient. "c:/temp/x.boot" would store all bootstrap draws on disk and in the specified folder. Note: this feature is completely experimental. Please report bugs if you wish to use this feature. x.boot can also add more bootstrap draws (bootstrap [... do something] add more bootstrap draws [... do something].
(2) Request coefficients with x.boot.out
x.boot.out accepts the following arguments:
ci = The confidence level.
ci.type = Type of the interval. Specifically, "norm", "basic", "perc" (aka "percentile"), "bc" (aka "bca.simple").
what = For which coefficients should the confidence interval be calculated? Accepted are "est" (estimates) "std.lv", "std.all", "std.nox" (see lavaan parameterEstimates) and all coefficients that lavInspect provides and make sense to bootstrap (examples: "rsquared", "cov.lv", "cor.ov", "resid").
myFUN = Custom Function.
myArgs = Arguments passed to the Custom Function as a list
myFUN.post = Manipulate results and custom output after calculating CI.
p.adjust = Request p-values that are consistent with the confidence intervals. This feature is entirely experimental. I cannot even name any reference for this. Please report any inconsistency!
Here are some examples (Note: Every line of code is at your disposal. You can customize everything and use or change any Custom Function - see below. I stress again: Everything is experimental and I do not guarantee the correctness of the output!
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)
Shu Fai Cheung
HG Wells
Christian Arnold
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
Ibrahim Nasser
This should be converted into a package