Nothing
#' @title Likelihood-Based Confidence Interval
#'
#' @description Find the likelihood-based confidence intervals (LBCIs) for
#' selected free parameters in an SEM output.
#'
#' @details [semlbci()] finds the positions of the selected parameters
#' in the parameter table and then calls [ci_i_one()] once for each
#' of them. For the technical details, please see [ci_i_one()] and
#' the functions it calls to find a confidence bound, currently
#' [ci_bound_wn_i()]. [ci_bound_wn_i()] uses the approach proposed by
#' Wu and Neale (2012) and illustrated by Pek and Wu (2015).
#'
#' It supports updating an output of [semlbci()] by setting
#' `semlbci_out`. This allows forming LBCIs for some parameters after
#' those for some others have been formed.
#'
#' If possible, parallel processing should be used (see `parallel` and
#' `ncpus`), especially for a model with many parameters.
#'
#' If the search for some of the confidence bounds failed, with `NA` for the
#' bounds, try increasing `try_k_more_times`.
#'
#' The SEM output will first be checked by [check_sem_out()] to see
#' whether the model and the estimation method are supported. To skip this
#' test (e.g., for testing or experimenting with some models and estimators),
#' set `check_fit` to `FALSE`.
#'
#' Examples and technical details can be found at Cheung
#' and Pesigan (2023), the website of the `semlbci`
#' package (https://sfcheung.github.io/semlbci/),
#' and the technical appendices at
#' (https://sfcheung.github.io/semlbci/articles/).
#'
#' It currently supports [lavaan::lavaan-class] outputs only.
#'
#' @return A `semlbci`-class object similar to the parameter table
#' generated by [lavaan::parameterEstimates()], with the LBCIs for
#' selected parameters added. Diagnostic information, if requested,
#' will be included in the attributes. See [print.semlbci()] for options
#' available.
#'
#' @param sem_out The SEM output. Currently supports
#' [lavaan::lavaan-class] outputs only.
#'
#' @param pars The positions of the parameters for which the LBCIs are
#' to be searched. Use the position as appeared on the parameter
#' tables of the `sem_out`. If `NULL`, the default, then LBCIs for
#' all free parameters will be searched. Can also be a vector of
#' strings to indicate the parameters on the parameter table. The
#' parameters should be specified in [lavaan::lavaan()] syntax. The
#' vector of strings will be converted by [syntax_to_i()] to
#' parameter positions. See [syntax_to_i()] on how to specify the parameters.
#'
#' @param include_user_pars Logical. Whether all user-defined parameters
#' are automatically included when `pars` is not set. Default is `TRUE`.
#' If `pars` is explicitly set, this argument will be ignored.
#'
#' @param remove_variances Logical. Whether variances and error variances
#' will be removed. Default is `TRUE`, removing all variances and error
#' variances even if specified in `pars`.
#'
#' @param remove_intercepts Logical. Whether intercepts will be removed.
#' Default is `TRUE`, removing all intercepts (parameters with operator `~1`).
#' Intercepts are not yet supported in standardized solution and so will
#' always be removed if `standardized = TRUE`.
#'
#' @param ciperc The proportion of coverage for the confidence
#' interval. Default is .95, requesting a 95 percent confidence
#' interval.
#'
#' @param standardized If `TRUE`, the LBCI is for the standardized estimates.
#'
#' @param method The method to be used to search for the confidence
#' bounds. Supported methods are`"wn"` (Wu-Neale-2012), the default,
#' and `"ur"` (root finding by [stats::uniroot()]).
#'
#' @param robust Whether the LBCI based on robust likelihood ratio
#' test is to be found. Only `"satorra.2000"` in [lavaan::lavTestLRT()]
#' is supported for now, implemented by the method proposed by Falk
#' (2018). If `"none"`, the default, then likelihood ratio test based
#' on maximum likelihood estimation will be used.
#'
#' @param try_k_more_times How many more times to try if failed.
#' Default is 2.
#'
#' @param semlbci_out An `semlbci-class` object. If provided, parameters already
#' with LBCIs formed will be excluded from `pars`.
#'
#' @param check_fit If `TRUE` (default), the input (`sem_out`) will
#' be checked by [check_sem_out()]. If not
#' supported, an error will be raised. If `FALSE`,
#' the check will be skipped and the LBCIs will be
#' searched even for a model or parameter not
#' supported. Set to `TRUE` only for testing.
#'
#' @param ... Arguments to be passed to [ci_bound_wn_i()].
#'
#' @param parallel If `TRUE`, will use parallel processing to do the search.
#'
#' @param ncpus The number of workers, if `parallel` is `TRUE`.
#' Default is 2. This number should not be larger than the number CPU
#' cores.
#'
#' @param use_pbapply If `TRUE` and `pbapply`
#' is installed, [pbapply::pbapply()] will be used to display a
#' progress bar when finding the intervals. Default is `TRUE`.
#' Ignored if `pbapply` is not installed.
#'
#' @param loadbalancing Whether load
#' balancing is used when `parallel`
#' is `TRUE` and `use_pbapply` is
#' `TRUE`.
#'
#' @author Shu Fai Cheung <https://orcid.org/0000-0002-9871-9448>
#'
#' @references
#'
#' Cheung, S. F., & Pesigan, I. J. A. (2023). *semlbci*:
#' An R package for forming likelihood-based confidence
#' intervals for parameter estimates, correlations,
#' indirect effects, and other derived parameters.
#' *Structural Equation Modeling: A Multidisciplinary Journal*,
#' *30*(6), 985--999.
#' \doi{10.1080/10705511.2023.2183860}
#'
#' Falk, C. F. (2018). Are robust standard errors the best approach
#' for interval estimation with nonnormal data in structural equation
#' modeling? *Structural Equation Modeling: A Multidisciplinary
#' Journal, 25*(2), 244-266.
#' \doi{10.1080/10705511.2017.1367254}
#'
#' Pek, J., & Wu, H. (2015). Profile likelihood-based confidence
#' intervals and regions for structural equation models.
#' *Psychometrika, 80*(4), 1123-1145.
#' \doi{10.1007/s11336-015-9461-1}
#'
#' Wu, H., & Neale, M. C. (2012). Adjusted confidence intervals for a
#' bounded parameter. *Behavior Genetics, 42*(6), 886-898.
#' \doi{10.1007/s10519-012-9560-z}
#'
#'
#' Pritikin, J. N., Rappaport, L. M., & Neale, M. C. (2017). Likelihood-based
#' confidence intervals for a parameter with an upper or lower bound.
#' *Structural Equation Modeling: A Multidisciplinary Journal, 24*(3), 395-401.
#' \doi{10.1080/10705511.2016.1275969}
#'
#' @seealso [print.semlbci()], [confint.semlbci()], [ci_i_one()], [ci_bound_wn_i()]
#'
#' @examples
#'
#' library(lavaan)
#' mod <-
#' "
#' m ~ a*x
#' y ~ b*m
#' ab := a * b
#' "
#' fit_med <- sem(mod, simple_med, fixed.x = FALSE)
#' p_table <- parameterTable(fit_med)
#' p_table
#' lbci_med <- semlbci(fit_med,
#' pars = c("m ~ x",
#' "y ~ m",
#' "ab :="))
#' lbci_med
#'
#' @export
semlbci <- function(sem_out,
pars = NULL,
include_user_pars = TRUE,
remove_variances = TRUE,
remove_intercepts = TRUE,
ciperc = .95,
standardized = FALSE,
method = c("wn", "ur"),
robust = c("none", "satorra.2000"),
try_k_more_times = 2,
semlbci_out = NULL,
check_fit = TRUE,
...,
parallel = FALSE,
ncpus = 2,
use_pbapply = TRUE,
loadbalancing = TRUE) {
if (!inherits(sem_out, "lavaan")) {
stop("sem_out is not a supported object.")
}
if (!is.null(semlbci_out)) {
if (!inherits(semlbci_out, "semlbci")) {
stop("semlbci_out is not an output of semlbci().")
}
}
method <- match.arg(method)
robust <- match.arg(robust)
sem_out_name <- deparse(substitute(sem_out))
# Use satorra.2000 automatically if
# - a scaled test is used and
# - the method is "ur"
scaled <- any(names(lavaan::lavInspect(sem_out, "test")) %in%
c("satorra.bentler",
"yuan.bentler",
"yuan.bentler.mplus",
"mean.var.adjusted",
"scaled.shifted"))
if (scaled && (method == "ur")) {
robust <- "satorra.2000"
}
# Check sem_out
if (check_fit) {
sem_out_check <- check_sem_out(sem_out = sem_out, robust = robust)
if (sem_out_check > 0) {
msg <- paste0(paste(attr(sem_out_check, "info"), collapse = "\n"), "\n")
warning(msg, immediate. = TRUE)
}
if (sem_out_check < 0) {
msg <- paste0(paste(attr(sem_out_check, "info"), collapse = "\n"), "\n")
stop(msg)
}
}
ptable <- as.data.frame(lavaan::parameterTable(sem_out))
pars_lbci_yes <- rep(FALSE, nrow(ptable))
if (!is.null(semlbci_out)) {
# Check whether semlbci_out and sem_out match
tmp0 <- ptable[, c("id", "lhs", "op", "rhs", "group", "label")]
tmp1 <- as.data.frame(semlbci_out)[, c("id", "lhs", "op", "rhs", "group", "label")]
tmp0 <- tmp0[order(tmp0$id), ]
tmp1 <- tmp1[order(tmp1$id), ]
rownames(tmp0) <- NULL
rownames(tmp1) <- NULL
if (!identical(tmp0, tmp1)) {
stop("semlbci_out and the parameter table of sem_out do not match.")
}
# Find pars with both lbci_lb and lbci_ub
pars_lbci_yes <- !sapply(semlbci_out$lbci_lb, is.na) &
!sapply(semlbci_out$lbci_ub, is.na)
}
i <- ptable$free > 0
i_id <- ptable$id
i_id_free <- i_id[i]
i_id_user <- i_id[(ptable$free == 0) & (ptable$op == ":=")]
# pars must be the row numbers as in the lavaan parameterTable.
# Support operators: "~", "~~", "=~", ":="
pars_tmp <- fix_pars(pars = pars,
sem_out = sem_out,
standardized = standardized,
include_user_pars = include_user_pars,
remove_variances = remove_variances,
remove_intercepts = remove_intercepts,
pars_lbci_yes = pars_lbci_yes,
i_id = i_id,
i_id_free = i_id_free,
i_id_user = i_id_user)
pars <- pars_tmp$pars
i_selected <- pars_tmp$i_selected
if (length(pars) == 0) {
if (is.null(semlbci_out)) {
stop("The number of parameters selected is zero.")
} else {
warning("No parameter selected. semlbci_out returned.")
return(semlbci_out)
}
}
npar <- sum(i)
# KEEP for reference: environment(set_constraint) <- parent.frame()
if (method == "wn") {
f_constr <- eval(set_constraint(sem_out = sem_out, ciperc = ciperc),
envir = parent.frame())
} else {
f_constr <- NULL
}
# Compute the scaling factors
if (robust == "satorra.2000") {
sf_full_list <- lapply(pars,
scaling_factor3,
sem_out = sem_out,
standardized = standardized,
std_method = "internal",
sem_out_name = sem_out_name)
} else {
sf_full_list <- rep(NA, length(pars))
}
if (parallel) {
# Parallel
cl <- parallel::makeCluster(ncpus)
# Rebuild the environment
pkgs <- .packages()
pkgs <- rev(pkgs)
parallel::clusterExport(cl, "pkgs", envir = environment())
parallel::clusterEvalQ(cl, {sapply(pkgs,
function(x) library(x, character.only = TRUE))
})
parallel::clusterExport(cl, ls(envir = parent.frame()),
envir = environment())
if (requireNamespace("pbapply", quietly = TRUE) &&
use_pbapply) {
if (loadbalancing) {
pboptions_old <- pbapply::pboptions(use_lb = TRUE)
# Restore pboptions on exit
on.exit(pbapply::pboptions(pboptions_old))
}
# Use pbapply
args_final <- utils::modifyList(list(...),
list(npar = npar,
sem_out = sem_out,
standardized = standardized,
debug = FALSE,
f_constr = f_constr,
method = method,
ciperc = ciperc,
robust = robust,
try_k_more_times = try_k_more_times,
sem_out_name = sem_out_name))
pars2 <- rep(pars, each = 2)
sf_full_list2 <- rep(sf_full_list, each = 2)
which2 <- rep(c("lbound", "ubound"), times = length(pars))
plist <- mapply(function(x, y, z) list(i = x, sf_full = y, which = z),
x = pars2, y = sf_full_list2, z = which2,
SIMPLIFY = FALSE)
tmpfct <- function(x) {
force(args_final)
args_final1 <- utils::modifyList(args_final, x)
do.call(semlbci::ci_i_one, args_final1)
}
out_raw2 <- pbapply::pblapply(plist,
tmpfct,
cl = cl)
} else {
# No progress bar
args_final <- utils::modifyList(list(...),
list(npar = npar,
sem_out = sem_out,
standardized = standardized,
debug = FALSE,
f_constr = f_constr,
method = method,
ciperc = ciperc,
robust = robust,
try_k_more_times = try_k_more_times,
sem_out_name = sem_out_name))
pars2 <- rep(pars, each = 2)
sf_full_list2 <- rep(sf_full_list, each = 2)
which2 <- rep(c("lbound", "ubound"), times = length(pars))
out_raw2 <- parallel::clusterMap(cl,
semlbci::ci_i_one,
i = pars2,
sf_full = sf_full_list2,
which = which2,
MoreArgs = args_final,
SIMPLIFY = FALSE,
.scheduling = ifelse(loadbalancing,
yes = "dynamic",
no = "static"))
}
parallel::stopCluster(cl)
} else {
# Serial
# Progress bar not support when ran in serial.
args_final <- utils::modifyList(list(...),
list(npar = npar,
sem_out = sem_out,
standardized = standardized,
debug = FALSE,
f_constr = f_constr,
method = method,
ciperc = ciperc,
robust = robust,
try_k_more_times = try_k_more_times,
sem_out_name = sem_out_name))
pars2 <- rep(pars, each = 2)
sf_full_list2 <- rep(sf_full_list, each = 2)
which2 <- rep(c("lbound", "ubound"), times = length(pars))
if (requireNamespace("pbapply", quietly = TRUE) &&
use_pbapply) {
out_raw2 <- pbapply::pbmapply(
ci_i_one,
i = pars2,
sf_full = sf_full_list2,
which = which2,
MoreArgs = args_final,
SIMPLIFY = FALSE)
} else {
out_raw2 <- mapply(
ci_i_one,
i = pars2,
sf_full = sf_full_list2,
which = which2,
MoreArgs = args_final,
SIMPLIFY = FALSE)
}
}
tmpfct2 <- function(x, y) {
out <- mapply(c, x, y, SIMPLIFY = FALSE)
out$method <- out$method[[1]]
out$sf_full <- out$sf_full[[1]]
out
}
tmp1 <- seq(1, 2 * length(pars), 2)
tmp2 <- tmp1 + 1
out_raw <- mapply(tmpfct2,
out_raw2[tmp1],
out_raw2[tmp2],
SIMPLIFY = FALSE)
out <- do.call(rbind, lapply(out_raw, function(x) x$bounds))
if (is.null(semlbci_out)) {
# Build the output
out_p <- ptable[, c("id", "lhs", "op", "rhs", "group", "label")]
if (standardized) {
pstd <- lavaan::standardizedSolution(sem_out)
if (lavaan::lavTech(sem_out, "ngroups") == 1) {
pstd$group <- 1
}
pstd$group[pstd$op == ":="] <- 0
out_p <- merge(out_p, pstd[, c("lhs", "op", "rhs", "group", "est.std")],
by = c("lhs", "op", "rhs", "group"), all.x = TRUE, sort = FALSE)
out_p <- out_p[order(out_p$id), ]
} else {
out_p$est <- ptable[, c("est")]
}
out_p$lbci_lb <- NA
out_p$lbci_ub <- NA
} else {
# Use existing output
out_p <- semlbci_out
}
out_p[i_selected, "lbci_lb"] <- out[, 1]
out_p[i_selected, "lbci_ub"] <- out[, 2]
# Collect diagnostic info
if (is.null(semlbci_out)) {
# Create the holder
q <- length(i_id)
lb_diag <- as.list(rep(NA, q))
ub_diag <- as.list(rep(NA, q))
lb_time <- as.numeric(rep(NA, q))
ub_time <- as.numeric(rep(NA, q))
lb_out <- as.list(rep(NA, q))
ub_out <- as.list(rep(NA, q))
ci_method <- as.character(rep(NA, q))
scaling_factor <- as.list(rep(NA, q))
} else {
# Retrieve from existing output
q <- length(i_id)
lb_diag <- attr(semlbci_out, "lb_diag")
ub_diag <- attr(semlbci_out, "ub_diag")
lb_time <- attr(semlbci_out, "lb_time")
ub_time <- attr(semlbci_out, "ub_time")
lb_out <- attr(semlbci_out, "lb_out")
ub_out <- attr(semlbci_out, "ub_out")
ci_method <- attr(semlbci_out, "ci_method")
scaling_factor <- attr(semlbci_out, "scaling_factor")
}
# Update the output
lb_diag[pars] <- lapply(out_raw, function(x) x$diags$lb_diag)
ub_diag[pars] <- lapply(out_raw, function(x) x$diags$ub_diag)
lb_time[pars] <- sapply(out_raw, function(x) x$times$lb_time)
ub_time[pars] <- sapply(out_raw, function(x) x$times$ub_time)
lb_out[pars] <- lapply(out_raw, function(x) x$ci_bound_i_out$lb_out)
ub_out[pars] <- lapply(out_raw, function(x) x$ci_bound_i_out$ub_out)
ci_method[pars] <- sapply(out_raw, function(x) x$method)
scaling_factor[pars] <- lapply(out_raw, function(x) x$sf_full)
p_names <- sapply(pars, i_to_name, sem_out = sem_out)
names(lb_diag)[pars] <- p_names
names(ub_diag)[pars] <- p_names
names(lb_time)[pars] <- p_names
names(ub_time)[pars] <- p_names
names(lb_out)[pars] <- p_names
names(ub_out)[pars] <- p_names
names(ci_method)[pars] <- p_names
names(scaling_factor)[pars] <- p_names
attr(out_p, "lb_diag") <- lb_diag
attr(out_p, "ub_diag") <- ub_diag
attr(out_p, "lb_time") <- lb_time
attr(out_p, "ub_time") <- ub_time
attr(out_p, "lb_out") <- lb_out
attr(out_p, "ub_out") <- ub_out
attr(out_p, "ci_method") <- ci_method
attr(out_p, "scaling_factor") <- scaling_factor
attr(out_p, "call") <- match.call()
# Append diagnostic info
out_p[pars, "status_lb"] <- sapply(lb_diag[pars], function(x) x$status)
out_p[pars, "status_ub"] <- sapply(ub_diag[pars], function(x) x$status)
out_p[pars, "ci_org_lb"] <- sapply(lb_diag[pars], function(x) x$ci_org_limit)
out_p[pars, "ci_org_ub"] <- sapply(ub_diag[pars], function(x) x$ci_org_limit)
out_p[pars, "ratio_lb"] <- sapply(lb_diag[pars], function(x) x$ci_limit_ratio)
out_p[pars, "ratio_ub"] <- sapply(ub_diag[pars], function(x) x$ci_limit_ratio)
out_p[pars, "post_check_lb"] <-
sapply(lb_diag[pars], function(x) x$fit_post_check)
out_p[pars, "post_check_ub"] <-
sapply(ub_diag[pars], function(x) x$fit_post_check)
out_p[pars, "time_lb"] <- as.numeric(lb_time[pars])
out_p[pars, "time_ub"] <- as.numeric(ub_time[pars])
out_p[pars, "cl_lb"] <- sapply(lb_diag[pars], function(x) x$ciperc_final)
out_p[pars, "cl_ub"] <- sapply(ub_diag[pars], function(x) x$ciperc_final)
out_p[pars, "method"] <- ci_method[pars]
if (robust != "none") {
out_p[pars, "robust"] <- robust
}
class(out_p) <- c("semlbci", class(out_p))
out_p
}
#' @noRd
# Process the `pars` argument
fix_pars <- function(pars,
sem_out,
standardized,
include_user_pars,
remove_variances,
remove_intercepts,
pars_lbci_yes,
i_id,
i_id_free,
i_id_user) {
if (!is.null(pars)) {
# Parameters supplied
if (is.character(pars)) {
pars <- pars_op(pars, sem_out = sem_out)
# Convert syntax to row numbers.
pars <- syntax_to_i(pars, sem_out)
}
if (standardized) {
# Always exclude variances fixed in the standardized solution
pars <- remove_v1(pars, sem_out)
}
# Remove parameters already with LBCIs in semlbci_out.
pars <- pars[!pars %in% i_id[pars_lbci_yes]]
# Ids in the vector of free parameters
i_selected <- i_id[pars]
} else {
# Start with all free parameters
pars <- i_id_free
if (standardized) {
# Intercepts and means not supported for standardized solution
remove_intercepts <- TRUE
# Always exclude variances fixed in the standardized solution
pars <- remove_v1(pars, sem_out)
# Include parameters free in the standardized solution,
# e.g., fixed loadings.
pars <- sort(unique(c(pars, free_in_std(i_id, sem_out))))
}
if (include_user_pars && length(i_id_user) > 0) {
pars <- c(pars, i_id_user)
}
if (remove_variances) {
pars <- remove_variances(pars, sem_out)
}
if (remove_intercepts) {
pars <- remove_intercepts(pars, sem_out)
}
# Remove parameters already with LBCIs in semlbci_out.
pars <- pars[!pars %in% i_id[pars_lbci_yes]]
# Ids in the vector of free parameters
i_selected <- i_id[pars]
}
list(pars = pars,
i_selected = i_selected)
}
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.