Nothing
#' @title Likelihood-based Confidence Bound By Wu-Neale-2012
#'
#' @description User the method proposed by Wu and Neale
#' (2012) to find the lower or upper bound of the likelihood-based
#' confidence interval (LBCI) for one parameter in a structural
#' equation model fitted in [lavaan::lavaan()].
#'
#' @details
#'
#' ## Important Notice
#'
#' This function is not supposed to be used directly by users in
#' typical scenarios. Its interface is user-*unfriendly* because it
#' should be used through [semlbci()]. It is exported such that
#' interested users can examine how a confidence bound is found, or
#' use it for experiments or simulations.
#'
#' ## Usage
#'
#' This function is the lowest level function used by [semlbci()].
#' [semlbci()] calls this function once for each bound of each
#' parameter. To use it, [set_constraint()] needs to be called first
#' to create the equality constraint required by the algorithm
#' proposed by Wu and Neale (2012).
#'
#' ## Algorithm
#'
#' This function implements the algorithm presented in Wu and Neale
#' (2012; see also Pek & Wu, 2015, Equation 12) that estimates all
#' free parameters in the optimization.
#'
#' ## Limitation(s)
#'
#' This function does not yet implement the method by Wu and Neale
#' (2012) for an estimate close to an attainable bound.
#'
#' @return A `cibound`-class object which is a list with three elements:
#'
#' - `bound`: A single number. The value of the bound located. `NA` is
#' the search failed for various reasons.
#'
#' - `diag`: A list of diagnostic information.
#'
#' - `call`: The original call.
#'
#' A detailed and organized output can be printed by the default print
#' method ([print.cibound()]).
#'
#' @param i The position of the target parameter as appeared in the
#' parameter table of an lavaan object, generated by
#' [lavaan::parameterTable()].
#'
#' @param npar The number of free parameters, including those
#' constrained to be equal.
#'
#' @param sem_out The fit object. Currently supports
#' [lavaan::lavaan-class] objects only.
#'
#' @param f_constr The constraint function generated by
#' [set_constraint()].
#'
#' @param which Whether the lower bound or the upper bound is to be
#' found. Must be `"lbound"` or `"ubound"`.
#'
#' @param history Not used. Kept for backward compatibility.
#'
#' @param perturbation_factor A number multiplied to the parameter
#' estimates in `sem_out`. Using the parameter estimates as starting
#' values may lead to errors in the first few iterations. Default is
#' .90. This argument is ignored if `wald_ci_start` is `TRUE.
#'
#' @param lb_var The lower bound for free parameters that are
#' variances. If equal to `-Inf`, the default, `lb_prop` and
#' `lb_se_k` will be used to set the lower bounds for free variances.
#' If it is a number, it will be used to set the lower bounds for all
#' free variances.
#'
#' @param wald_ci_start If `TRUE`, there are no equality
#' constraints in the model, and the target parameter is not a
#' user-defined parameter, the Wald confidence bounds will be used as
#' the starting value.
#'
#' @param standardized If `TRUE`, the LBCI is for the requested
#' estimate in the standardized solution. Default is `FALSE`.
#'
#' @param opts Options to be passed to [nloptr::nloptr()], the current
#' optimizer. Default is `list()`.
#'
#' @param ciperc The intended coverage probability for the confidence
#' interval. Default is .95, and the bound for a 95% confidence
#' interval will be sought.
#'
#' @param ci_limit_ratio_tol The tolerance for the ratio of `a` to
#' `b`, where `a` is the distance between an LBCI limit and the point
#' estimate, and the `b` is the distance between the original
#' confidence limit (by default the Wald CI in [lavaan::lavaan()])
#' and the point estimate. If the ratio is larger than this value or
#' smaller than the reciprocal of this value, a warning is set in the
#' status code. Default is 1.5.
#'
#' @param verbose If `TRUE`, the function will store more diagnostic
#' information in the attribute `diag`. Default is `FALSE`.
#'
#' @param sf A scaling factor. Used for robust confidence bounds.
#' Default is 1. Computed by an internal function called by
#' [semlbci()] when `robust = "satorra.2000"`.
#'
#' @param sf2 A shift factor. Used for robust confidence bounds.
#' Default is 0. Computed by an internal function called by
#' [semlbci()] when `robust = "satorra.2000"`.
#'
#' @param p_tol Tolerance for checking the achieved level of
#' confidence. If the absolute difference between the achieved level
#' and `ciperc` is greater than this amount, a warning is set in the
#' status code and the bound is set to `NA`. Default is 5e-4.
#'
#' @param std_method The method used to find the standardized
#' solution. If equal to `"lavaan"`,
#' [lavaan::standardizedSolution()] will be used. If equal to
#' `"internal"`, an internal function will be used. The `"lavaan"`
#' method should work in all situations, but the `"internal"` method
#' is usually much faster. Default is `"internal"`.
#'
#' @param bounds Default is `""` and this function will set the lower
#' bounds to `lb_var` for variances. Other valid values are those
#' accepted by [lavaan::lavaan()]. Ignored for now.
#'
#' @param xtol_rel_factor Multiply the default `xtol_rel` by a number,
#' usually a positive number equal to or less than 1, to change the
#' default termination criterion. Default is 1.
#'
#' @param ftol_rel_factor Multiply the default `ftol_rel` by a number,
#' usually a positive number equal to or less than 1, to change the
#' default termination criterion. Default is 1.
#'
#' @param lb_prop Used by an internal function to set the lower bound
#' for free variances. Default is .05, setting the lower bound to
#' .05 * estimate. Used only if the lower bound set by `lb_se_k` is
#' negative.
#'
#' @param lb_se_k Used by an internal function to set the lower bound
#' for free variances. Default is 3, the estimate minus 3 standard
#' error. If negative, the lower bound is set using `lb_prop`.
#'
#' @param try_harder If error occurred
#' in the optimization, how many more
#' times to try. In each new attempt,
#' the starting values will be randomly
#' jittered. Default is 0.
#'
#' @param fit_lb The vector of lower
#' bounds of parameters. Default is
#' `-Inf`, setting the lower bounds to
#' `-Inf` for all parameters except for
#' free variances which are controlled
#' by `lb_var`.
#'
#' @param fit_ub The vector of upper
#' bounds of parameters. Default is
#' `+Inf`, setting the lower bounds to
#' `+Inf` for all parameters.
#'
#' @param timeout The approximate
#' maximum time for the search, in
#' second. Default is 300 seconds
#' (5 minutes).
#'
#' @param ... Optional arguments. Not used.
#'
#' @references
#'
#' 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}
#'
#'
#' @seealso [print.cibound()], [semlbci()], [ci_i_one()]
#'
#' @examples
#'
#' data(simple_med)
#' dat <- simple_med
#'
#' mod <-
#' "
#' m ~ x
#' y ~ m
#' "
#'
#' fit_med <- lavaan::sem(mod, simple_med, fixed.x = FALSE)
#'
#' fn_constr0 <- set_constraint(fit_med)
#'
#' out1l <- ci_bound_wn_i(i = 1,
#' npar = 5,
#' sem_out = fit_med,
#' f_constr = fn_constr0,
#' which = "lbound")
#' out1l
#'
#' @export
ci_bound_wn_i <- function(i = NULL,
npar = NULL,
sem_out = NULL,
f_constr = NULL,
which = NULL,
history = FALSE,
perturbation_factor = .9,
lb_var = -Inf,
standardized = FALSE,
wald_ci_start = !standardized,
opts = list(),
ciperc = .95,
ci_limit_ratio_tol = 1.5,
verbose = FALSE,
sf = 1,
sf2 = 0,
p_tol = 5e-4,
std_method = "internal",
bounds = "none",
xtol_rel_factor = 1,
ftol_rel_factor = 1,
lb_prop = .05,
lb_se_k = 3,
try_harder = 0,
fit_lb = -Inf,
fit_ub = +Inf,
timeout = 300,
...) {
k <- switch(which,
lbound = 1,
ubound = -1)
# Check if the parameter is a user-defined parameter
p_table <- lavaan::parameterTable(sem_out)
i_op <- p_table[i, "op"]
# The id in the vector of free parameters
i_in_free <- p_table[i, "free"]
# Get original point estimate and CI
if (standardized) {
## Standardized solution
p_est <- lavaan::standardizedSolution(sem_out,
type = "std.all",
se = TRUE,
zstat = FALSE,
pvalue = FALSE,
ci = TRUE,
level = ciperc,
remove.eq = FALSE,
remove.ineq = FALSE,
remove.def = FALSE,
output = "data.frame")
i_est <- p_est[i, "est.std"]
i_org_ci_lower <- p_est[i, "ci.lower"]
i_org_ci_upper <- p_est[i, "ci.upper"]
if (std_method == "internal") {
# Test whether the internal function can reproduce the
# solution from lavaan::standardizedSolution().
p_est_int <- std_lav(lavaan::coef(sem_out),
sem_out)
p_est_test <- lavaan::parameterEstimates(
sem_out,
standardized = TRUE,
se = FALSE,
zstat = FALSE,
pvalue = FALSE,
ci = FALSE,
cov.std = FALSE,
remove.system.eq = FALSE,
remove.eq = FALSE,
remove.ineq = FALSE,
remove.def = FALSE,
remove.nonfree = FALSE)
if (lavaan::lavTech(sem_out, "ngroups") > 1) {
p_est_test <- p_est_test[, c("lhs", "op", "rhs",
"group",
"std.all")]
} else {
p_est_test <- p_est_test[, c("lhs", "op", "rhs",
"std.all")]
}
p_est_int_tmp <- p_est_test
p_est_int_tmp$std.all <- p_est_int
p_est_test <- p_est_test[!(p_est_test$op %in% c("~1", "==")), ]
p_est_int_tmp <- p_est_int_tmp[!(p_est_int_tmp$op %in% c("~1", "==")), ]
if (all(p_est_test$std.all != p_est_int_tmp$std.all)) {
stop(paste0("The internal method cannot reproduce std.all.\n",
"Please use std_method = ", dQuote("lavaan"), "."))
}
}
} else {
## Unstandardized solution
p_est <- lavaan::parameterEstimates(sem_out,
se = TRUE,
ci = TRUE,
level = ciperc,
zstat = FALSE,
fmi = FALSE,
rsquare = TRUE,
output = "data.frame")
i_est <- p_est[i, "est"]
i_org_ci_lower <- p_est[i, "ci.lower"]
i_org_ci_upper <- p_est[i, "ci.upper"]
}
# Create the objective function and the gradient function
## Standardized solution
if (standardized && (i_op %in% c("=~", "~", "~~", ":="))) {
if (!(std_method %in% c("lavaan", "internal"))) {
stop(paste0("std_method must be either ",
dQuote("lavaan"), " or ",
dQuote("internal")))
}
if (std_method == "lavaan") {
p_std <- lavaan::standardizedSolution(sem_out,
type = "std.all",
se = FALSE,
zstat = FALSE,
pvalue = FALSE,
ci = FALSE,
remove.eq = FALSE,
remove.ineq = FALSE,
remove.def = FALSE,
output = "data.frame")
p_std$id <- seq_len(nrow(p_std))
if (lavaan::lavTech(sem_out, "ngroups") > 1) {
i_lor <- get_lhs_op_rhs(i, sem_out, more = TRUE)
i_std <- merge(p_std, i_lor, by = c("lhs", "op", "rhs", "group"))$id
} else {
i_lor <- get_lhs_op_rhs(i, sem_out)
i_std <- merge(p_std, i_lor, by = c("lhs", "op", "rhs"))$id
}
start0 <- lavaan::parameterTable(sem_out)
# The function to be minimized.
lbci_b_f <- function(param, sem_out, debug, lav_warn, sf, sf2, stop_on_error = FALSE) {
start1 <- start0
start1[start1$free > 0, "est"] <- param
sem_out2 <- sem_out
sem_out2@ParTable <- as.list(start1)
sem_model <- sem_out2@Model
sem_model <- lavaan::lav_model_set_parameters(sem_model,
start1[start1$free > 0, "est"])
sem_out2@Model <- sem_model
std0 <- lavaan::standardizedSolution(sem_out2,
type = "std.all",
se = FALSE,
zstat = FALSE,
pvalue = FALSE,
ci = FALSE,
remove.eq = FALSE,
remove.ineq = FALSE,
remove.def = FALSE,
output = "data.frame")
k * std0[i_std, "est.std"]
}
}
if (std_method == "internal") {
p_std <- lavaan::parameterEstimates(
sem_out,
standardized = TRUE,
se = FALSE,
zstat = FALSE,
pvalue = FALSE,
ci = FALSE,
cov.std = FALSE,
remove.system.eq = FALSE,
remove.eq = FALSE,
remove.ineq = FALSE,
remove.def = FALSE,
remove.nonfree = FALSE)
p_std$id <- seq_len(nrow(p_std))
if (lavaan::lavTech(sem_out, "ngroups") > 1) {
i_lor <- get_lhs_op_rhs(i, sem_out, more = TRUE)
i_std <- merge(p_std, i_lor, by = c("lhs", "op", "rhs", "group"))$id
} else {
i_lor <- get_lhs_op_rhs(i, sem_out)
i_std <- merge(p_std, i_lor, by = c("lhs", "op", "rhs"))$id
}
start0 <- lavaan::parameterTable(sem_out)
# The function to be minimized.
lbci_b_f <- function(param, sem_out, debug, lav_warn, sf, sf2, stop_on_error = FALSE) {
std0 <- tryCatch(std_lav(param, sem_out),
error = function(e) e)
if (inherits(std0, "error")) {
if (stop_on_error) {
stop("Error in std_lav().")
} else {
return(Inf)
}
} else {
return(k * std0[i_std])
}
}
}
# The gradient of the function to be minimized
lbci_b_grad <- function(param, sem_out, debug, lav_warn, sf, sf2, stop_on_error = FALSE) {
lavaan::lav_func_gradient_complex(
lbci_b_f, param, sem_out = sem_out,
debug = debug, lav_warn = lav_warn,
sf = sf,
sf2 = sf2
)
}
}
## Unstandardized solution
if (i_op == ":=") {
## User-defined parameter
if (standardized) {
# Not used
} else {
# Get the name of the defined parameter
i_name <- p_table[i, "label"]
# The function to be minimized.
lbci_b_f <- function(param, sem_out, debug, lav_warn, sf, sf2, stop_on_error = FALSE) {
out <- tryCatch(k * sem_out@Model@def.function(param)[i_name],
error = function(e) e)
if (inherits(out, "error")) {
if (stop_on_error) {
stop("Error in cmoputing user-parameter(s).")
} else {
return(Inf)
}
} else {
return(out)
}
}
# The gradient of the function to be minimized
lbci_b_grad <- function(param, sem_out, debug, lav_warn, sf, sf2, stop_on_error = FALSE) {
lavaan::lav_func_gradient_complex(
lbci_b_f, param, sem_out = sem_out,
debug = debug, lav_warn = lav_warn,
sf = sf,
sf2 = sf2
)
}
}
} else {
## Free parameter
if (standardized) {
# Not used
} else {
# The function to be minimized.
# force() may not be needed but does not harm to keep them.
lbci_b_f <- function(param, sem_out, debug, lav_warn, sf, sf2, stop_on_error = FALSE) {
force(k)
force(i_in_free)
k * param[i_in_free]
}
# The gradient of the function to be minimized
grad_c <- rep(0, npar)
grad_c[i_in_free] <- k
lbci_b_grad <- function(param, sem_out, debug, lav_warn, sf, sf2, stop_on_error = FALSE) {
grad_c
}
}
}
# Set the starting values of the parameters
if (wald_ci_start) {
xstart <- set_start_wn(i, sem_out,
which = which,
standardized = standardized,
ciperc = ciperc)
xstart <- xstart[xstart$free > 0, "est"]
} else {
xstart <- perturbation_factor * lavaan::coef(sem_out)
}
# Get SEs. For bounds
p_table <- lavaan::parameterTable(sem_out)
se_free <- p_table[p_table$free > 0, "se"]
est_free <- p_table[p_table$free > 0, "est"]
# Set lower bounds
if (fit_lb == -Inf) {
fit_lb <- rep(-Inf, npar)
} else if (length(fit_lb) == 1) {
fit_lb <- est_free - se_free * fit_lb
}
if (isTRUE(identical(lb_var, -Inf))) {
# Override lb_var = -Inf for free variances
fit_lb[find_variance_in_free(sem_out)] <- find_variance_in_free_lb(sem_out,
prop = lb_prop,
se_k = lb_se_k)
} else {
fit_lb[find_variance_in_free(sem_out)] <- lb_var
}
# Set upper bounds
if (fit_ub == -Inf) {
fit_ub <- rep(-Inf, npar)
} else if (length(fit_ub) == 1) {
fit_ub <- est_free + se_free * fit_ub
}
# For pmax() and pmin()
fit_lb_na <- fit_lb
fit_ub_na <- fit_ub
fit_lb_na[fit_lb_na == -Inf] <- NA
fit_ub_na[fit_ub_na == +Inf] <- NA
# Need further testing on using lavaan bounds
# if (!is.null(p_table$lower)) {
# fit_lb <- p_table$lower[p_table$free > 0]
# } else {
# if (bounds == "") {
# fit_lb <- rep(-Inf, npar)
# fit_lb[find_variance_in_free(sem_out)] <- lb_var
# } else {
# p_bounds <- set_bounds(sem_out = sem_out, bounds = bounds)
# fit_lb <- p_bounds$lower[p_bounds$free > 0]
# }
# }
# if (!is.null(p_table$upper)) {
# fit_ub <- p_table$upper[p_table$free > 0]
# } else {
# if (bounds == "") {
# fit_ub <- rep(+Inf, npar)
# } else {
# p_bounds <- set_bounds(sem_out = sem_out, bounds = bounds)
# fit_ub <- p_bounds$upper[p_bounds$free > 0]
# }
# }
# Optimization
# Default options in nloptr
opts_final <- utils::modifyList(list("algorithm" = "NLOPT_LD_SLSQP",
# "xtol_rel" = 1.0e-10,
"xtol_rel" = 1.0e-5 * xtol_rel_factor,
"ftol_rel" = 1.0e-5 * ftol_rel_factor,
"maxeval" = 500,
"print_level" = 0,
"maxtime" = timeout),
opts)
try_harder <- as.integer(max(try_harder, 0))
try_harder_count <- 0
# Make sure the starting values are within the bounds
xstart <- pmin(pmax(xstart,
fit_lb_na,
na.rm = TRUE),
fit_ub_na,
na.rm = TRUE)
xstart_i <- xstart
while(try_harder_count <= try_harder) {
out <- tryCatch(nloptr::nloptr(
x0 = xstart_i,
eval_f = lbci_b_f,
lb = fit_lb,
ub = fit_ub,
eval_grad_f = lbci_b_grad,
eval_g_eq = f_constr,
opts = opts_final,
sem_out = sem_out,
lav_warn = FALSE,
debug = FALSE,
sf = sf,
sf2 = sf2,
stop_on_error = TRUE),
error = function(e) e)
if (inherits(out, "error")) {
xstart_i <- stats::runif(length(xstart_i),
min = -2,
max = 2) * xstart_i
xstart_i <- pmin(pmax(xstart_i,
fit_lb_na,
na.rm = TRUE),
fit_ub_na,
na.rm = TRUE)
try_harder_count <- try_harder_count + 1
} else {
try_harder_count <- Inf
}
}
# Failed after try harder k times
if (inherits(out, "error")) {
search_error <- as.character(out)
opts_error <- utils::modifyList(opts_final,
list("maxeval" = 1))
out <- nloptr::nloptr(x0 = xstart,
eval_f = lbci_b_f,
lb = fit_lb,
ub = fit_ub,
eval_grad_f = lbci_b_grad,
eval_g_eq = f_constr,
opts = opts_error,
sem_out = sem_out,
lav_warn = FALSE,
debug = FALSE,
sf = sf,
sf2 = sf2,
stop_on_error = TRUE)
} else {
search_error <- NULL
}
# Process the results
## Need as.numeric() to handle NA
bound <- as.numeric(k * out$objective)
bound_unchecked <- bound
# Initialize the status code
## Any values other than 0 denotes a problem
status <- 0
# Check convergence
if (out$status < 0) {
# Verify the range of nloptr that denotes an error
status <- 1
check_optimization <- FALSE
# bound <- NA
} else {
check_optimization <- TRUE
}
# Check the limit
i_org_ci_limit <- switch(which,
lbound = i_org_ci_lower,
ubound = i_org_ci_upper)
if (!is.na(bound)) {
ci_limit_ratio <- abs((bound - i_est) / (i_org_ci_limit - i_est))
ci_limit_ratio2 <- 1 / ci_limit_ratio
if (any(ci_limit_ratio > ci_limit_ratio_tol,
ci_limit_ratio2 > ci_limit_ratio_tol)) {
# status <- 1
# Do not set the bound to NA because the limit may still be valid.
}
} else {
ci_limit_ratio <- NA
}
# Check whether admissible
start0 <- lavaan::parameterTable(sem_out)
i_free <- find_free(sem_out)
start0[i_free, "est"] <- out$solution
start0[i_free, "start"] <- out$solution
tmp <- sem_out@Options
tmp$check.start <- TRUE
tmp$check.post <- TRUE
tmp$check.vcov <- TRUE
tmp$warn <- FALSE
tmp$do.fit <- FALSE
tmp$start <- start0
fit_final <- lavaan::lavaan(model = start0,
slotOptions = tmp,
# slotParTable = tmp2,
# slotModel = tmp3,
slotSampleStats = sem_out@SampleStats,
slotData = sem_out@Data)
# fit_final <- lavaan::update(sem_out, start = start0, do.fit = FALSE,
# check.start = TRUE,
# check.post = TRUE,
# check.vcov = TRUE)
fit_post_check <- lavaan::lavInspect(fit_final, "post.check")
if (!fit_post_check) {
status <- 1
check_post_check <- FALSE
# bound <- NA
# The warning should be raised by the calling function, not this one.
} else {
check_post_check <- TRUE
}
# Achieved level of confidence
if (lavaan::lavTech(sem_out, "ngroups") > 1) {
ntotal <- lavaan::lavTech(sem_out, "ntotal")
} else {
ntotal <- lavaan::lavTech(sem_out, "nobs")
}
fmin_org <- lavaan::lavTech(sem_out, "optim")$fx
fmin_final <- lavaan::lavTech(fit_final, "optim")$fx
chisq_diff <- (fmin_final - fmin_org) * 2 * ntotal / sf + sf2
ciperc_final <- stats::pchisq(chisq_diff, 1)
if (abs(ciperc_final - ciperc) > p_tol) {
status <- 1
check_level_of_confidence <- FALSE
# bound <- NA
} else {
check_level_of_confidence <- TRUE
}
if (!all(check_optimization,
check_post_check,
check_level_of_confidence)) {
bound <- NA
}
## Collect diagnostic information
diag <- list(status = status,
est_org = i_est,
ci_org_limit = i_org_ci_limit,
ci_limit_ratio = ci_limit_ratio,
fit_post_check = fit_post_check,
start_values = xstart,
bound_unchecked = bound_unchecked,
org_values = lavaan::coef(sem_out),
final_values = lavaan::coef(fit_final),
ciperc = ciperc,
ciperc_final = ciperc_final,
i = i,
i_lor = get_lhs_op_rhs(i, sem_out, more = TRUE),
optim_message = out$message,
optim_iterations = out$iterations,
optim_status = out$status,
optim_termination_conditions = out$termination_conditions,
method = "wn",
which = which,
standardized = standardized,
check_optimization = check_optimization,
check_post_check = check_post_check,
check_level_of_confidence = check_level_of_confidence,
search_error = search_error
)
if (verbose) {
diag$history <- out
diag$fit_final <- fit_final
} else {
diag$history <- NULL
diag$fit_final <- NULL
}
if (out$status < 0) {
# If convergence status < 0, override verbose
diag$history <- out
}
out <- list(bound = bound,
diag = diag,
call = match.call(),
method_name = "Wu-Neale-2012")
class(out) <- c("cibound", class(out))
out
}
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.