diff --git a/DESCRIPTION b/DESCRIPTION index e718af9..d4e5115 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: semlbci Title: Likelihood-Based Confidence Interval in Structural Equation Models -Version: 0.10.4 +Version: 0.10.4.24 Authors@R: c( person(given = "Shu Fai", @@ -28,7 +28,7 @@ License: GPL-3 Encoding: UTF-8 LazyData: true Roxygen: list(markdown = TRUE) -RoxygenNote: 7.2.3 +RoxygenNote: 7.3.1 Suggests: testthat (>= 3.0.0), knitr, @@ -44,8 +44,10 @@ Imports: ggplot2, ggrepel, rlang, - pbapply + pbapply, + callr, + methods VignetteBuilder: knitr Config/testthat/parallel: true Config/testthat/edition: 3 -Config/testthat/start-first: semlbci_wn_mg_* +Config/testthat/start-first: semlbci_wn_mg_*, ur_ci_bound_ur* diff --git a/NAMESPACE b/NAMESPACE index 999e4ef..27df2f3 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -6,10 +6,16 @@ S3method(print,ci_order) S3method(print,cibound) S3method(print,semlbci) export(check_sem_out) +export(ci_bound_ur) +export(ci_bound_ur_i) export(ci_bound_wn_i) export(ci_i_one) export(ci_order) +export(gen_est_i) +export(gen_sem_out_userp) +export(gen_userp) export(get_cibound) +export(get_cibound_status_not_0) export(loglike_compare) export(loglike_point) export(loglike_quad_point) diff --git a/NEWS.md b/NEWS.md index 839e821..74d271b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,110 @@ +# semlbci 0.10.4.24 + +## New Feature + +- Add the method `"ur"` for forming the + LBCI. It uses root finding without + derivatives. It is inefficient before + a model with equality constraint is + fitted in each iteration. However, it + may be able to find the bound for a + parameter when the `"wn"` method + cannot. It also supports robust LBCI + using Satorra (2000) chi-square + difference test. (0.10.4.2) + +- Add `get_cibound_status_not_0()`. It + checks the status of each bound in + a `semlbci` object, and returns as + a list the `cibound` objects of + bounds with status not equal to zero. + For diagnostic purpose. (0.10.4.13) + +## Miscellaneous + +- Suppressed a harmless warning in a + test. (0.10.4.1) +- Added the option to use load balancing + when calling `semlbci()` with + `use_pbapply` set to `TRUE`. Enabled + by default. (0.10.4.3) +- Add error handlers in the output + of `set_constraints()` and + `ci_bound_wn_i()`. (0.10.4.4) +- Revised `ci_bound_wn_i()` and + `ci_bound_ur_i()` to make + sure the bound is a numeric object, + even if `NA`. (0.10.4.5) +- Speed up the examples of + `ci_bound_ur_i()` and `ci_bound_ur()`. + (0.10.4.7) +- Precompute the interval in + `ci_bound_ur_i()` before calling + `ci_bound_ur()`. (0.10.4.8) +- Use `kill()` instead of `close()` + in R sessions. (0.10.4.9) +- Add `bound_start` and `user_est` + to the precomputed interval. + (0.10.4.10) +- Modified several functions related to + `"ur"` method to have the option to + use an existing R session instead of + creating a new one. (0.10.4.11) +- Store the error in the search and + display it in the printout. + (0.10.4.12) +- Fixed some typos in doc. + (0.10.4.14) +- Made the test for the match between + `sem_out` and `semlbci_out` in + `semlbci()` more robust. (0.10.4.15) +- Fixed the order of the output when + `standardized` is `TRUE`. (0.10.4.16) +- Fixed a bug in `print.cibound()` when + the call is of the form `xxx::yyy()`. + (0.10.4.17) +- Added `fit_lb` and `fit_ub` arguments + to `ci_bound_wn_i()` for setting + the bounds. (0.10.4.18, 0.10.4.19) +- Fixed a bug with `std_lav()` for + models with only one factor. + (0.10.4.20) +- Added the `timeout` argument to + `ci_bound_wn_i()`, default to 300 + (300 seconds or 5 minutes). + (0.10.4.21) +- In `print.semlbci()`, added + suggestions on what to do if some + bounds could not be found (0.10.4.22) +- Use whitespace instead of tab to + align the output of `print.cibound()`. + (0.10.4.23) +- Fixed a bug with `semlbci()`. + Intercepts are now automatically + skipped if `standardized` is `TRUE`. + (0.10.4.24) + +## (Possibly) Breaking Changes + +- This may or may not be a breaking + changes. The default of + `try_k_more_times` in `semlbci()` was + changed from 2 to 0. These additional + attempts rarely help but usually + increase processing time + unnecessarily. (0.10.4.6) +- Disabled further attempts in + the `wn` method by default + as they rarely help but they usually + increase + the processing time unnecessarily. + If necessary, one of the set of + attempts, successfully lower the + lower limits of variances can be + enabled again by setting `try_lb` + to `TRUE` when calling `semlbci()` + (0.10.4.6) + # semlbci 0.10.4 ## New Feature diff --git a/R/ci_bound_ur_i.R b/R/ci_bound_ur_i.R new file mode 100644 index 0000000..2f35893 --- /dev/null +++ b/R/ci_bound_ur_i.R @@ -0,0 +1,952 @@ +#' @title Likelihood-Based Confidence +#' Bound By Root Finding +#' +#' @description Using root finding +#' 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. +#' +#' For consistency in the interface, +#' most of the arguments in +#' [ci_bound_wn_i()] are also included +#' in this function, even those not used +#' internally. +#' +#' ## Algorithm +#' +#' This function, unlike +#' [ci_bound_wn_i()], use a simple root +#' finding algorithm. Basically, it tries +#' fixing the target parameter to +#' different values until the likelihood +#' ratio test *p*-value, or the +#' corresponding chi-square difference, +#' is equal to the +#' value corresponding to the desired +#' level of confidence. (Internally, +#' the difference between the *p*-value +#' and the target *p*-value, that for +#' the chi-square difference, is the +#' function value.) +#' +#' For finding the bound, this algorithm +#' can be inefficient compared to the +#' one proposed by Wu and Neale (2012). +#' The difference can be less than +#' one second versus 10 seconds. +#' It is included as a backup algorithm +#' for parameters which are difficult +#' for the method by Wu and Neale. +#' +#' Internally, it uses [uniroot()] to +#' find the root. +#' +#' ## Limitation(s) +#' +#' This function does not handle an +#' estimate close to an attainable bound +#' using the method proposed by Wu and +#' Neale (2012). Use it for such +#' parameters with cautions. +#' +#' @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 a lavaan object, +#' generated by +#' [lavaan::parameterTable()]. +#' +#' @param npar Ignored by this function. +#' Included consistency in the +#' interface. +#' +#' @param sem_out The fit object. +#' Currently supports +#' [lavaan::lavaan-class] objects only. +#' +#' @param f_constr Ignored by this +#' function. Included consistency in the +#' interface. +#' +#' @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 Ignored by +#' this function. Included consistency +#' in the interface. +#' +#' @param lb_var Ignored by this +#' function. Included consistency in the +#' interface. +#' +#' @param wald_ci_start Ignored by this +#' function. Included consistency in the +#' interface. +#' +#' @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 +#' [stats::uniroot()]. 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 Ignored by this function. +#' Included consistency in the +#' interface. +#' +#' @param sf2 Ignored by this function. +#' Included consistency in the +#' interface. +#' +#' @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 Ignored by this +#' function. Included consistency in the +#' interface. +#' +#' @param xtol_rel_factor Ignored by +#' this function. Included consistency +#' in the interface. +#' +#' @param ftol_rel_factor Ignored by +#' this function. Included consistency +#' in the interface. +#' +#' @param lb_prop Ignored by this +#' function. Included consistency in the +#' interface. +#' +#' @param lb_se_k Ignored by this +#' function. Included consistency in the +#' interface. +#' +#' @param d A value used to determine +#' the width of the interval in the +#' initial search. Larger this value, +#' *narrow* the interval. Default is 5. +#' Used by [ci_bound_ur()]. +#' +#' @param ... Optional arguments. Not +#' used. +#' +#' @references +#' +#' 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()]; see +#' [ci_bound_wn_i()] on the version +#' for the method by Wu and Neale +#' (2012). +#' +#' @examples +#' +#' library(lavaan) +#' data(simple_med) +#' dat <- simple_med +#' mod <- +#' " +#' m ~ x +#' y ~ m +#' " +#' fit_med <- sem(mod, simple_med, fixed.x = FALSE) +#' # Remove `opts` in real cases. +#' # The options are added just to speed up the example +#' out1l <- ci_bound_ur_i(i = 1, +#' sem_out = fit_med, +#' which = "lbound", +#' opts = list(use_callr = FALSE, +#' interval = c(0.8277, 0.8278))) +#' out1l +#' @export + +ci_bound_ur_i <- function(i = NULL, + npar = NULL, # Not used by ur + sem_out = NULL, + f_constr = NULL, # Not used by ur + which = NULL, + history = FALSE, + perturbation_factor = .9, # Can be used by sem_out_userp() + lb_var = -Inf, # Not used by ur. Set in sem_out + standardized = FALSE, + wald_ci_start = !standardized, # Not used by ur + opts = list(), + ciperc = .95, + ci_limit_ratio_tol = 1.5, + verbose = FALSE, + sf = 1, # Not used by ur + sf2 = 0, # Not used by ur + p_tol = 5e-4, + std_method = "internal", # Not used by ur + bounds = "none", # Not used by ur + xtol_rel_factor = 1, # Not used by ur + ftol_rel_factor = 1, # Not used by ur + lb_prop = .05, # Not used by ur + lb_se_k = 3, # Not used by ur + d = 5, + ...) { + + # Basic info + + # 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_limit <- switch(which, + lbound = p_est[i, "ci.lower"], + ubound = p_est[i, "ci.upper"]) + } 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_limit <- switch(which, + lbound = p_est[i, "ci.lower"], + ubound = p_est[i, "ci.upper"]) + } + + # Initial interval to search + a <- abs(i_org_ci_limit - i_est) / d + interval0 <- c(i_org_ci_limit - a, i_org_ci_limit + a) + attr(interval0, "bound_start") <- i_org_ci_limit + attr(interval0, "user_est") <- i_est + + npar <- lavaan::lavTech(sem_out, "npar") + + # Generate the user-parameter function + est_i_func <- gen_est_i(i = i, + sem_out = sem_out, + standardized = standardized) + + # Find the root + on.exit(try(rs$kill(), silent = TRUE), add = TRUE) + rs <- callr::r_session$new() + ur_opts <- list(sem_out = sem_out, + func = est_i_func, + level = ciperc, + which = which, + interval = interval0, + uniroot_maxiter = 500, + use_callr = FALSE, + rs = rs) + ur_opts <- utils::modifyList(ur_opts, + opts) + out <- do.call(ci_bound_ur, + ur_opts) + + # Process the results + + bound <- as.numeric(out$bound) + bound_unchecked <- bound + + # Initialize the status code + ## Any values other than 0 denotes a problem + status <- 0 + + # Check convergence + # NOTE: + # - Not necessary to check the code + # because the validity of the bound + # can be determined by the LR test. + # - Find the status code of uniroot + # - Need to be done in the call to uniroot() + # 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 + + 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 + + fit_final <- out$sem_out_bound + + fit_post_check <- lavaan::lavTech(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 + + ciperc_final <- 1 - out$lrt[2, "Pr(>Chisq)"] + 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 + } + + coef0 <- methods::getMethod("coef", + signature = "lavaan", + where = asNamespace("lavaan"))(sem_out) + coef_final <- methods::getMethod("coef", + signature = "lavaan", + where = asNamespace("lavaan"))(fit_final) + xstart <- coef0 + xstart[] <- 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 = unname(xstart), + bound_unchecked = bound_unchecked, + org_values = coef0, + final_values = coef_final, + ciperc = ciperc, + ciperc_final = ciperc_final, + i = i, + i_lor = get_lhs_op_rhs(i, sem_out, more = TRUE), + optim_message = "Nil", # out$message + optim_iterations = out$optimize_out$iter, + optim_status = NA, # out$status + optim_termination_conditions = "Nil", # out$termination_conditions, + method = "ur", + which = which, + standardized = standardized, + check_optimization = check_optimization, + check_post_check = check_post_check, + check_level_of_confidence = check_level_of_confidence + ) + if (verbose) { + out0 <- out + out0$sem_out_bound <- NULL + diag$history <- out0 + diag$fit_final <- fit_final + } else { + diag$history <- NULL + diag$fit_final <- NULL + } + if (status < 0) { + # If convergence status < 0, override verbose + diag$history <- out + } + out <- list(bound = bound, + diag = diag, + call = match.call(), + method_name = "Root Finding") + class(out) <- c("cibound", class(out)) + out + } + +#' @title Find a Likelihood-Based +#' Confidence Bound By Root Finding +#' +#' @description 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()] +#' using [uniroot()]. +#' +#' @details +#' This function is called xby +#' [ci_bound_ur_i()]. This function is +#' exported because it is a +#' stand-alone function that can be used +#' directly for any function that +#' receives a lavaan object and returns +#' a scalar. +#' +#' The function [ci_bound_ur_i()] is a +#' wrapper of this function, with an +#' interface similar to that of +#' [ci_bound_wn_i()] and returns a +#' `cibound`-class object. The +#' user-parameter function is generated +#' internally by [ci_bound_wn_i()]. +#' +#' This function, on the other hand, +#' requires users to supply the function +#' directly through the `func` argument. +#' This provides the flexibility to find +#' the bound for any function of the +#' model parameter, even one that cannot +#' be easily coded in `lavaan` model +#' syntax. +#' +#' @param sem_out The fit object. +#' Currently supports +#' [lavaan::lavaan-class] objects only. +#' +#' @param func A function that receives +#' a lavaan object and returns a scalar. +#' This function is to be used by +#' [gen_userp()] and so there are +#' special requirements on it. +#' Alternatively, it can be the output +#' of [gen_est_i()]. +#' +#' @param ... Optional arguments to be +#' passed to `func`. Usually not used +#' but included in case the function +#' has such arguments. +#' +#' @param level The level of confidence +#' of the confidence interval. Default +#' is .95, or 95%. +#' +#' @param which Whether the lower bound +#' or the upper bound is to be found. +#' Must be `"lbound"` or `"ubound"`. +#' +#' @param interval A numeric vector of +#' two values, which is the initial +#' interval to be searched. If `NULL`, +#' the default, it will be determined +#' internally using Wald or delta +#' method confidence interval, if +#' available. +#' +#' @param progress Whether progress will +#' be reported on screen during the +#' search. Default is +#' `FALSE`. +#' +#' @param method The actual function to +#' be used in the search. which can only +#' be `"uniroot"`, the default, for now. +#' May include other function in the +#' future. +#' +#' @param lrt_method The method used in +#' [lavaan::lavTestLRT()]. Default is +#' `"default"`. It is automatically set +#' to `"satorra.2000"` and cannot be +#' overridden if a scaled test statistic +#' is requested in `sem_out`. +#' +#' @param tol The tolerance used in +#' [uniroot()], default is .005. +#' +#' @param root_target Whether the +#' chi-square difference (`"chisq"`), +#' the default, or its *p*-value +#' (`"pvalue"`) is used as the function +#' value in finding the root. Should have +#' little impact on the results. +#' +#' @param d A value used to determine +#' the width of the interval in the +#' initial search. Larger this value, +#' *narrow* the interval. Default is 5. +#' +#' @param uniroot_extendInt To be passed +#' to the argument `extendInt` of +#' [uniroot()]. Whether the interval +#' should be extended if the root is not +#' found. Default is `"yes"`. Refer to +#' the help page of [uniroot()] for +#' possible values. +#' +#' @param uniroot_trace To be passed to +#' the argument `trace` of [uniroot()]. +#' How much information is printed +#' during the search. Default is 0, and +#' no information is printed during the +#' search. Refer to the help page of +#' [uniroot()] for possible values. +#' +#' @param uniroot_maxiter The maximum +#' number of iteration in the search. +#' Default is 1000. +#' +#' @param use_callr Whether the +#' `callr` package will be used to +#' do the search in a separate R +#' process. Default is `TRUE`. Should +#' not set to `FALSE` if used in an +#' interactive environment unless this is +#' intentional. +#' +#' @param rs Optional. If set to +#' a persistent R process created by +#' `callr`, it will be used instead of +#' starting a new one, and it will not +#' be terminated on exit. +#' +#' @return +#' The function [ci_bound_ur()] returns +#' a list with the following elements: +#' +#' - `bound`: The bound found. +#' +#' - `optimize_out`: THe output of the +#' root finding function, [uniroot()] +#' for now. (Called `optimize_out` +#' because an earlier version of this +#' function also uses [optimize()]). +#' +#' - `sem_out_bound`: The `lavaan` model +#' with the user-defined parameter fixed +#' to the bound. +#' +#' - `lrt`: The output of +#' [lavaan::lavTestLRT()] comparing +#' `sem_out` and `sem_out_bound`. +#' +#' - `bound_start`: The Wald or delta +#' method confidence bound returned when +#' determining the interval internally. +#' +#' - `user_est`: The estimate of the +#' user-defined parameter when +#' determining the interval internally. +#' +#' @examples +#' +#' library(lavaan) +#' data(simple_med) +#' dat <- simple_med +#' mod <- +#' " +#' m ~ x +#' y ~ m +#' " +#' fit_med <- lavaan::sem(mod, simple_med, fixed.x = FALSE) +#' parameterTable(fit_med) +#' # Create a function to get the second parameter +#' est_i <- gen_est_i(i = 2, sem_out = fit_med) +#' # Find the lower bound of the likelihood-based confidence interval +#' # of the second parameter. +#' # user_callr should be TRUE or omitted in read research. +#' # Remove interval in read research. It is added to speed up the example. +#' out1l <- ci_bound_ur(sem_out = fit_med, +#' func = est_i, +#' which = "lbound", +#' use_callr = FALSE, +#' interval = c(.39070, .39075)) +#' out1l +#' +#' @rdname ci_bound_ur +#' +#' @export + +ci_bound_ur <- function(sem_out, + func, + ..., + level = .95, + which = c("lbound", "ubound"), + interval = NULL, + progress = FALSE, + method = "uniroot", + lrt_method = "default", + tol = .0005, # This is the tolerance in x, not in y + root_target = c("chisq", "pvalue"), + d = 5, + uniroot_extendInt = "yes", + uniroot_trace = 0, + uniroot_maxiter = 1000, + use_callr = TRUE, + rs = NULL) { + + # Internal Workflow + # - Define the function that works on the lavaan object. + # - Call add_func() to generate a modified lavaan object. + # - Call sem_out_userp_run() to fix to a value + + # satorra.2000 is used automatically and cannot be overridden + scaled <- any(names(lavaan::lavInspect(sem_out, "test")) %in% + c("satorra.bentler", + "yuan.bentler", + "yuan.bentler.mplus", + "mean.var.adjusted", + "scaled.shifted")) + if (scaled) { + lrt_method <- "satorra.2000" + } + + which <- match.arg(which) + root_target <- match.arg(root_target) + # Only support uniroot() for now + method <- match.arg(method) + lrt_alpha <- 1 - level + # Create the modified model + fit_i <- add_func(func = func, + sem_out = sem_out) + if (is.null(interval)) { + # Set the interval to search + if (progress) { + cat("\nFinding the initial interval ...") + } + interval <- uniroot_interval(sem_out = sem_out, + func = func, + which = which, + d = d) + if (progress) { + cat("\rInitial interval computed. ") + } + } + bound_start <- attr(interval, "bound_start") + user_est <- attr(interval, "user_est") + bound_start <- ifelse(is.null(bound_start), NA, bound_start) + user_est <- ifelse(is.null(user_est), NA, user_est) + if (method == "uniroot") { + # Define the function for which to find the root + lrtp_i <- function(x, + alpha, + root_target = "pvalue", + global_ok = FALSE, + rs = NULL) { + sem_out_x <- sem_out_userp_run(target = x, + object = fit_i, + global_ok = global_ok, + rs = rs) + lrt0 <- lavaan::lavTestLRT(sem_out_x, + sem_out, + method = lrt_method) + cname <- switch(root_target, + pvalue = "Pr(>Chisq)", + chisq = "Chisq diff") + out1 <- lrt0[2, cname] + y_target <- switch(root_target, + pvalue = alpha, + chisq = stats::qchisq(level, df = 1)) + out1 - y_target + } + if (use_callr) { + if (is.null(rs)) { + on.exit(try(rs$kill(), silent = TRUE), add = TRUE) + rs <- callr::r_session$new() + } + # May fix: Cannot pass an R session to another R session + if (progress) { + optimize_out <- rs$call(function(f, + interval, + alpha, + root_target, + global_ok, + extendInt, + tol, + trace, + maxiter) { + stats::uniroot(f = f, + interval = interval, + alpha = alpha, + root_target = root_target, + global_ok = global_ok, + extendInt = extendInt, + tol = tol, + trace = trace, + maxiter = maxiter) + }, + args = list(f = lrtp_i, + interval = interval, + alpha = lrt_alpha, + root_target = root_target, + global_ok = TRUE, + extendInt = uniroot_extendInt, + tol = tol, + trace = uniroot_trace, + maxiter = uniroot_maxiter)) + cat("\nSearch started ...\n") + while (rs$poll_process(500) != "ready") { + rtime <- rs$get_running_time()["total"] + # TODO: + # - Need to handle time longer than 1 mins. + cat("\r", "Time spent: ", round(rtime, 2), " sec.", spe = "") + } + cat("\nSearch finished. Finalize the results ... \n") + optimize_out <- rs$read()$result + } else { + optimize_out <- rs$run(function(f, + interval, + alpha, + root_target, + global_ok, + extendInt, + tol, + trace, + maxiter) { + stats::uniroot(f = f, + interval = interval, + alpha = alpha, + root_target = root_target, + global_ok = global_ok, + extendInt = extendInt, + tol = tol, + trace = trace, + maxiter = maxiter) + }, + args = list(f = lrtp_i, + interval = interval, + alpha = lrt_alpha, + root_target = root_target, + global_ok = TRUE, + extendInt = uniroot_extendInt, + tol = tol, + trace = uniroot_trace, + maxiter = uniroot_maxiter)) + } + } else { + # Run locally. Need global_ok == FALSE. + if (progress) { + cat("\nSearch started ...\n") + } + if (is.null(rs)) { + on.exit(try(rs$kill(), silent = TRUE), add = TRUE) + rs <- callr::r_session$new() + } + optimize_out <- stats::uniroot(lrtp_i, + interval = interval, + alpha = lrt_alpha, + root_target = root_target, + global_ok = FALSE, + rs = rs, + extendInt = uniroot_extendInt, + tol = tol, + trace = uniroot_trace, + maxiter = uniroot_maxiter) + if (progress) { + cat("\nSearch finished. Finalize the results ... \n") + } + } + } + out_root <- optimize_out$root + sem_out_final <- sem_out_userp_run(target = out_root, + object = fit_i, + rs = rs) + lrt_final <- lavaan::lavTestLRT(sem_out_final, + sem_out, + method = lrt_method) + out <- list(bound = out_root, + optimize_out = optimize_out, + sem_out_bound = sem_out_final, + lrt = lrt_final, + bound_start = bound_start, + user_est = user_est) + out + } + +#' @param i The position of the target +#' parameter as appeared in the +#' parameter table of an lavaan object, +#' generated by +#' [lavaan::parameterTable()]. +#' +#' @param standardized If `TRUE`, the +#' standardized estimate is to be +#' retrieved. Default is `FALSE`. +#' Only support `"std.all"` for now. +#' +#' @return +#' The function [gen_est_i()] returns +#' a special function can inspects the +#' `Model` slot (and `implied` slot +#' if necessary) of a modified `lavaan` +#' object and return the parameter +#' estimate. This function is to be used +#' by [ci_bound_ur()] or +#' [gen_sem_out_userp()]. +#' +#' @rdname ci_bound_ur +#' +#' @export + +gen_est_i <- function(i, + sem_out, + standardized = FALSE) { + ptable <- lavaan::parameterTable(sem_out) + # i_user <- (ptable$free > 0) | (ptable$user > 0) + # i_est <- which(which(i_user) == i) + # From the help page of lavaan-class object, + # type = "user" should return *all* parameters + # in the parameter table, including equality constraints. + # Therefore, the row number should be equal to the order + # in the vector of coefficients. + i_user <- i + i_est <- i + is_def <- (ptable$op[i] == ":=") + i_lhs <- (ptable$lhs[i]) + i_rhs <- (ptable$rhs[i]) + i_op <- ptable$op[i] + i_gp <- ptable$group[i] + if (standardized) { + if (is_def) { + out <- function(object) { + # Just in case + force(i_est) + std <- lavaan::standardizedSolution(object, + est = lavaan::lav_model_get_parameters(object@Model, type = "user"), + GLIST = object@Model@GLIST, + se = FALSE, + zstat = FALSE, + pvalue = FALSE, + ci = FALSE, + remove.eq = FALSE, + remove.ineq = FALSE, + remove.def = FALSE, + type = "std.all") + unname(std[i_est, "est.std"]) + } + } else { + out <- function(object) { + # Just in case + force(i_est) + force(i_gp) + force(i_op) + force(i_lhs) + force(i_rhs) + object@implied <- lavaan::lav_model_implied(object@Model) + est <- lavaan::lav_model_get_parameters(object@Model, type = "user")[i_est] + implied <- lavaan::lavInspect(object, "cov.all", drop.list.single.group = FALSE) + implied_sd <- sqrt(diag(implied[[i_gp]])) + est1 <- switch(i_op, + `~~` = est / (implied_sd[i_lhs] * implied_sd[i_rhs]), + `~` = est * implied_sd[i_rhs] / implied_sd[i_lhs], + `=~` = est * implied_sd[i_lhs] / implied_sd[i_rhs]) + unname(est1) + } + } + } else { + out <- function(object) { + force(i_est) + unname(lavaan::lav_model_get_parameters(object@Model, type = "user")[i_est]) + } + } + return(out) + } diff --git a/R/ci_bound_ur_i_helpers.R b/R/ci_bound_ur_i_helpers.R new file mode 100644 index 0000000..7676646 --- /dev/null +++ b/R/ci_bound_ur_i_helpers.R @@ -0,0 +1,267 @@ + +#' @title Initial Interval for the +#' Search by 'ci_bound_ur()' +#' +#' @description Determine the interval to +#' be used by [ci_bound_ur()]. +#' +#' @details +#' Used when the parameter is the output +#' of a function a lavaan object, for +#' which it is not easy to get the Wald +#' or delta method confidence interval. +#' +#' It fits the model to the data in +#' `sem_out` and gets the confidence +#' interval of the function value of +#' `func`. The distance from the target +#' confidence bound (determiend by +#' `which`) to the estimate is computed, +#' divided by `d`. The confidence bound +#' plus and minus this value is the +#' interval to be used. +#' +#' @param sem_out The fit object. +#' Currently supports +#' [lavaan::lavaan-class] objects only. +#' +#' @param func A function that receives +#' a lavaan object and returns a scalar. +#' Usually the output of +#' [gen_sem_out_userp()]. +#' +#' @param which Whether the lower bound +#' or the upper bound is to be found. +#' Must be `"lbound"` or `"ubound"`. +#' +#' @param d A value used to determine +#' the width of the interval in the +#' initial search. Larger this value, +#' *narrow* the interval. Default is 5. +#' +#' @param level The level of confidence +#' of the confidence interval. Default +#' is .95, or 95%. +#' +#' @noRd +# CHECKED: The only change is the level argument +uniroot_interval <- function(sem_out, + func, + which, + d = 5, + level = .95) { + # Set the interval to search + + # NOTE: + # No need to know standardized or not because + # standardization, if any, is done inside func. + + fit_i_free <- add_func(func = func, + sem_out = sem_out, + fix = FALSE) + fit_start <- sem_out_userp_run(target = 0, + object = fit_i_free) + est_start <- lavaan::parameterEstimates(fit_start, + ci = TRUE, + level = level) + i <- which(est_start$rhs == paste0(fit_i_free$userp_name, "()")) + user_est <- est_start[i, "est"] + bound_start <- est_start[i, switch(which, + lbound = "ci.lower", + ubound = "ci.upper")] + bound_start_d <- abs(bound_start - user_est) + i_l <- bound_start - bound_start_d / d + i_u <- bound_start + bound_start_d / d + out <- c(i_l, i_u) + attr(out, "bound_start") <- bound_start + attr(out, "user_est") <- user_est + return(out) + } + + +#' @title Fix a User Parameter To a +#' Value and Refit a Model +#' +#' @description Fix a user parameter in +#' a lavaan model to a target value, +#' refit the model, and return the +#' results. +#' +#' @details +#' The model must be the output of +#' [add_func()]. +#' +#' @param target The value to which the +#' user parameter will be fixed to. +#' +#' @param object The output of +#' [add_func()]. +#' +#' @param verbose Whether diagnostic +#' information will be printed. Default +#' is `FALSE`. +#' +#' @param control To be passed to the +#' argument of the same name in +#' [lavaan::lavaan()]. Default is +#' `list()`. +#' +#' @param seed Numeric. If supplied, it +#' will be used in [set.seed()] to +#' initialize the random number +#' generator. Necessary to reproduce +#' some results because random numbers +#' are used in some steps in `lavaan`. +#' If `NULL`, the default, [set.seed()] +#' will not be called. +#' +#' @param global_ok Logical. Whether +#' is the user function can be stored +#' in the global environment. This +#' option is disabled for now because +#' it is not a good practice to change +#' the global environment. +#' +#' @return +#' A `lavaan` object, with the value +#' of the user-defined parameters +#' fixed to the target value. +#' +#' @noRd +# CHECKED: Identical to the experimental version +sem_out_userp_run <- function(target, + object, + verbose = FALSE, + control = list(), + seed = NULL, + global_ok = FALSE, + rs = NULL) { + + userp <- object$userp + userp_name <- object$userp_name + + # This option is disabled to avoid + # storing things to the global + # environment + global_ok <- FALSE + if (global_ok) { + # Disabled for now + } else { + if (is.null(rs)) { + on.exit(try(rs$kill(), silent = TRUE)) + rs <- callr::r_session$new() + } + out <- rs$run(function(target, + verbose, + control, + seed, + sem_out_userp, + userp, + userp_name) { + assign(userp_name, + value = userp, + envir = parent.frame()) + do.call(sem_out_userp, + list(target = target, + verbose = verbose, + control = control, + seed = seed)) + }, + args = list(target = target, + verbose = verbose, + control = control, + seed = seed, + sem_out_userp = object$sem_out_userp, + userp = userp, + userp_name = userp_name)) + } + out + } + + +#' @title Add a User Function to a 'lavaan' +#' model +#' +#' @description It adds a user function +#' to a `lavaan` model. +#' +#' @details +#' +#' The function is one accepts a +#' `lavaan` object and returns a scalar. +#' An R session will be started to +#' create the model. This is necessary +#' because the function needs to be +#' present in the global environment +#' even if the model is only created +#' but not fitted. +#' +#' @param func A function that accepts +#' a `lavaan` object and returns a +#' scalar. +#' +#' @param sem_out The fit object. +#' Currently supports +#' [lavaan::lavaan-class] objects only. +#' +#' @param userp_name The name of the +#' function (`func`) to be used in the +#' `lavaan` model. Should be changed +#' only if it conflicts with another +#' object in the parent environment, +#' which should not happen if the model +#' is always fitted in a clean R +#' session. +#' +#' @param fix To be passed to +#' [gen_sem_out_userp()]. If `TRUE`, the +#' default, the function generated is +#' used to fix +#' the value of `func` to a target value +#' using an equality constraint. If +#' `FALSE`, then the function simply +#' fits the model to the data. +#' +#' @return +#' A list with the following elements: +#' +#' - `sem_out_userp`: The `lavaan` model +#' with the user function, `func`, added +#' as a user-defined parameter. +#' +#' - `userp`: A copy of `func`. +#' +#' - `userp_name`: The value of +#' `userp_name`, the label of the +#' user-defined parameter in the model. +#' +#' @noRd +# CHECKED: Identical to the experimental version +add_func <- function(func, + sem_out, + userp_name = "semlbciuserp1234", + fix = TRUE, + rs = NULL) { + userp <- gen_userp(func = func, + sem_out = sem_out) + # Create a child process + if (is.null(rs)) { + on.exit(try(rs$kill(), silent = TRUE)) + rs <- callr::r_session$new() + } + fit_i <- rs$run(function(..., + userp, + userp_name) { + assign(userp_name, + value = userp, + parent.frame()) + semlbci::gen_sem_out_userp(...) + }, args = list(userp = userp, + sem_out = sem_out, + userp_name = userp_name, + fix = fix)) + list(sem_out_userp = fit_i, + userp = userp, + userp_name = userp_name) + } + diff --git a/R/ci_bound_wn_i.R b/R/ci_bound_wn_i.R index c850eac..f535e47 100644 --- a/R/ci_bound_wn_i.R +++ b/R/ci_bound_wn_i.R @@ -1,6 +1,7 @@ -#' @title Likelihood-based Confidence Bound For One parameter +#' @title Likelihood-based Confidence Bound By Wu-Neale-2012 #' -#' @description Find the lower or upper bound of the likelihood-based +#' @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()]. #' @@ -114,7 +115,7 @@ #' 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"``, +#' 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 @@ -141,6 +142,29 @@ #' 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 @@ -204,6 +228,10 @@ ci_bound_wn_i <- function(i = NULL, 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, @@ -310,7 +338,7 @@ ci_bound_wn_i <- function(i = NULL, } start0 <- lavaan::parameterTable(sem_out) # The function to be minimized. - lbci_b_f <- function(param, sem_out, debug, lav_warn, sf, sf2) { + 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 @@ -356,13 +384,22 @@ ci_bound_wn_i <- function(i = NULL, } start0 <- lavaan::parameterTable(sem_out) # The function to be minimized. - lbci_b_f <- function(param, sem_out, debug, lav_warn, sf, sf2) { - std0 <- std_lav(param, sem_out) - k * std0[i_std] + 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) { + 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, @@ -380,11 +417,21 @@ ci_bound_wn_i <- function(i = NULL, # 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) { - k * sem_out@Model@def.function(param)[i_name] + 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) { + 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, @@ -400,7 +447,7 @@ ci_bound_wn_i <- function(i = NULL, } 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) { + 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] @@ -408,7 +455,7 @@ ci_bound_wn_i <- function(i = NULL, # 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) { + lbci_b_grad <- function(param, sem_out, debug, lav_warn, sf, sf2, stop_on_error = FALSE) { grad_c } } @@ -425,8 +472,17 @@ ci_bound_wn_i <- function(i = NULL, 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 - fit_lb <- rep(-Inf, npar) + 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, @@ -437,7 +493,17 @@ ci_bound_wn_i <- function(i = NULL, } # Set upper bounds - fit_ub <- rep(+Inf, npar) + 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)) { @@ -470,25 +536,66 @@ ci_bound_wn_i <- function(i = NULL, "xtol_rel" = 1.0e-5 * xtol_rel_factor, "ftol_rel" = 1.0e-5 * ftol_rel_factor, "maxeval" = 500, - "print_level" = 0), + "print_level" = 0, + "maxtime" = timeout), opts) - 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_final, - sem_out = sem_out, - lav_warn = FALSE, - debug = FALSE, - sf = sf, - sf2 = sf2) + try_harder <- as.integer(max(try_harder, 0)) + try_harder_count <- 0 + 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), fit_ub_na) + try_harder_count <- try_harder_count + 1 + } else { + try_harder_count <- Inf + } + } - # Process the results + # Failed after try harder k times - bound <- k * out$objective + 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 @@ -600,7 +707,8 @@ ci_bound_wn_i <- function(i = NULL, standardized = standardized, check_optimization = check_optimization, check_post_check = check_post_check, - check_level_of_confidence = check_level_of_confidence + check_level_of_confidence = check_level_of_confidence, + search_error = search_error ) if (verbose) { diag$history <- out @@ -615,7 +723,8 @@ ci_bound_wn_i <- function(i = NULL, } out <- list(bound = bound, diag = diag, - call = match.call()) + call = match.call(), + method_name = "Wu-Neale-2012") class(out) <- c("cibound", class(out)) out } diff --git a/R/ci_i_one.R b/R/ci_i_one.R index 1dbf45b..a66a5c1 100644 --- a/R/ci_i_one.R +++ b/R/ci_i_one.R @@ -54,16 +54,19 @@ #' [lavaan::lavaan-class] outputs only. #' #' @param method The approach to be used. Default is `"wn"` -#' (Wu-Neale-2012 Method), the only supported method. +#' (Wu-Neale-2012 Method). Another method is "ur", +#' root finding by [stats::uniroot()]. #' #' @param standardized Logical. Whether the bound of the LBCI of the #' standardized solution is to be searched. Default is `FALSE`. #' #' @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. If `"none"``, the +#' [lavaan::lavTestLRT()] is supported for now. If `"none"`, the #' default, then likelihood ratio test based on maximum likelihood -#' estimation will be used. +#' estimation will be used. For "ur", `"satorra.2000"` is +#' automatically used if a scaled test statistic is requested +#' in `sem_out`. #' #' @param sf_full A list with the scaling and shift factors. Ignored #' if `robust` is `"none"`. If `robust` is `"satorra.2000"` and @@ -124,7 +127,7 @@ ci_i_one <- function(i, which = NULL, sem_out, - method = "wn", + method = c("wn", "ur"), standardized = FALSE, robust = "none", sf_full = NA, @@ -132,9 +135,7 @@ ci_i_one <- function(i, sem_out_name = NULL, try_k_more_times = 0, ...) { - if (!(which %in% c("lbound", "ubound"))) { - stop("Must be 'lbound' or 'ubound' for the 'which' argument.") - } + method <- match.arg(method) # It should be the job of the calling function to check whether it is # appropriate to use the robust method. if (tolower(robust) == "satorra.2000") { @@ -165,91 +166,33 @@ ci_i_one <- function(i, } if (method == "wn") { - # Wu-Neale-2012 method. - # The only method supported for now. - ## Attempt 1 - wald_ci_start <- !standardized - std_method_i <- "internal" - b_time <- system.time(b <- try(suppressWarnings(ci_bound_wn_i(i, - sem_out = sem_out, - which = which, - standardized = standardized, - sf = sf, - sf2 = sf2, - std_method = std_method_i, - wald_ci_start = wald_ci_start, - ...)), silent = TRUE)) - ## If "internal" failed, switches to "lavaan". - if (inherits(b, "try-error")) { - std_method_i <- "lavaan" - b_time <- b_time + system.time(b <- suppressWarnings(ci_bound_wn_i(i, - sem_out = sem_out, - which = which, - standardized = standardized, - sf = sf, - sf2 = sf2, - std_method = std_method_i, - ...))) - } - attempt_lb_var <- 0 - # Attempt 2 - if (b$diag$status != 0) { - ## Successively reduce the positive lower bounds for free variances - lb_se_k0 <- 10 - lb_prop0 <- .11 - lb_prop1 <- .01 - lb_propi <- lb_prop0 - while ((lb_propi > lb_prop1) & (b$diag$status != 0)) { - attempt_lb_var <- attempt_lb_var + 1 - lb_propi <- lb_propi - .01 - b_time <- b_time + system.time(b <- suppressWarnings(ci_bound_wn_i(i, - sem_out = sem_out, - which = which, - standardized = standardized, - sf = sf, - sf2 = sf2, - std_method = std_method_i, - lb_prop = lb_propi, - lb_se_k = lb_se_k0, - ...))) - } - } - # Attempt 3 - attempt_more_times <- 0 - if (b$diag$status != 0) { - # Try k more times - # Successively change the tolerance for convergence - ki <- try_k_more_times - fxi <- 1 - fti <- 1 - while ((ki > 0) & (b$diag$status != 0)) { - attempt_more_times <- attempt_more_times + 1 - ki <- ki - 1 - fxi <- fxi * .1 - if (ki > 0) { - fti <- fti * .1 - } else { - # Try hard in the last attempt - fti <- 0 - } - b_time <- b_time + system.time(b <- suppressWarnings(ci_bound_wn_i(i, - sem_out = sem_out, - which = which, - standardized = standardized, - sf = sf, - sf2 = sf2, - std_method = std_method_i, - xtol_rel_factor = fxi, - ftol_rel_factor = fti, - lb_prop = lb_propi, - lb_se_k = lb_se_k0, - ...))) - } - } + b_out <- ci_i_one_wn(i = i, + which = which, + sem_out = sem_out, + standardized = standardized, + sf = sf, + sf2 = sf2, + try_k_more_times = try_k_more_times, + ...) + } + if (method == "ur") { + # satorra.2000 is turned on automatically for ur + b_out <- ci_i_one_ur(i = i, + which = which, + sem_out = sem_out, + standardized = standardized, + try_k_more_times = try_k_more_times, + ...) } if (method == "nm") { stop("The method 'nm' is no longer supported.") } + + b <- b_out$b + b_time <- b_out$b_time + attempt_lb_var <- b_out$attempt_lb_var + attempt_more_times <- b_out$attempt_more_times + # MAY-FIX: # Should not name the elements based on the bound (lower/upper). # But this is not an issue for now. @@ -275,3 +218,148 @@ ci_i_one <- function(i, } out } + +#' @noRd + +ci_i_one_wn <- function(i, + which = c("lbound", "ubound"), + sem_out, + standardized = FALSE, + sf, + sf2, + try_k_more_times = 0, + lb_se_k0 = 10, + lb_prop0 = .11, + lb_prop1 = .01, + try_lb = FALSE, + ...) { + # Wu-Neale-2012 method. + # The only method supported for now. + ## Attempt 1 + wald_ci_start <- !standardized + std_method_i <- "internal" + b_time <- system.time(b <- try(suppressWarnings(ci_bound_wn_i(i, + sem_out = sem_out, + which = which, + standardized = standardized, + sf = sf, + sf2 = sf2, + std_method = std_method_i, + wald_ci_start = wald_ci_start, + ...)), silent = TRUE)) + ## If "internal" failed, switches to "lavaan". + if (inherits(b, "try-error")) { + std_method_i <- "lavaan" + b_time <- b_time + system.time(b <- suppressWarnings(ci_bound_wn_i(i, + sem_out = sem_out, + which = which, + standardized = standardized, + sf = sf, + sf2 = sf2, + std_method = std_method_i, + ...))) + } + attempt_lb_var <- 0 + # Attempt 2 + if ((b$diag$status != 0) && try_lb) { + ## Successively reduce the positive lower bounds for free variances + lb_propi <- lb_prop0 + while ((lb_propi > lb_prop1) & (b$diag$status != 0)) { + attempt_lb_var <- attempt_lb_var + 1 + lb_propi <- lb_propi - .01 + b_time <- b_time + system.time(b <- suppressWarnings(ci_bound_wn_i(i, + sem_out = sem_out, + which = which, + standardized = standardized, + sf = sf, + sf2 = sf2, + std_method = std_method_i, + lb_prop = lb_propi, + lb_se_k = lb_se_k0, + ...))) + } + } + # Attempt 3 + attempt_more_times <- 0 + if (b$diag$status != 0) { + # Try k more times + # Successively change the tolerance for convergence + if (!try_lb) { + lb_propi <- 1e-3 + } + ki <- try_k_more_times + fxi <- 1 + fti <- 1 + while ((ki > 0) && (b$diag$status != 0)) { + attempt_more_times <- attempt_more_times + 1 + ki <- ki - 1 + fxi <- fxi * .1 + if (ki > 0) { + fti <- fti * .1 + } else { + # Try hard in the last attempt + fti <- 0 + } + b_time <- b_time + system.time(b <- suppressWarnings(ci_bound_wn_i(i, + sem_out = sem_out, + which = which, + standardized = standardized, + sf = sf, + sf2 = sf2, + std_method = std_method_i, + xtol_rel_factor = fxi, + ftol_rel_factor = fti, + lb_prop = lb_propi, + lb_se_k = lb_se_k0, + ...))) + } + } + out <- list(b = b, + b_time = b_time, + attempt_lb_var = attempt_lb_var, + attempt_more_times = attempt_more_times) + return(out) + } + +#' @noRd + +ci_i_one_ur <- function(i, + which = c("lbound", "ubound"), + sem_out, + standardized = FALSE, + sf, + sf2, + try_k_more_times = 0, + ...) { + # Root finding by uniroot + ## Rarely need to try more than once. + std_method_i <- "internal" + b_time <- system.time(b <- try(suppressWarnings(ci_bound_ur_i(i, + sem_out = sem_out, + which = which, + standardized = standardized, + sf = sf, + sf2 = sf2, + std_method = std_method_i, + ...)), silent = TRUE)) + ## If "internal" failed, switches to "lavaan". + if (inherits(b, "try-error")) { + std_method_i <- "lavaan" + b_time <- b_time + system.time(b <- suppressWarnings(ci_bound_ur_i(i, + sem_out = sem_out, + which = which, + standardized = standardized, + sf = sf, + sf2 = sf2, + std_method = std_method_i, + ...))) + } + attempt_lb_var <- 0 + attempt_more_times <- 0 + out <- list(b = b, + b_time = b_time, + attempt_lb_var = attempt_lb_var, + attempt_more_times = attempt_more_times) + return(out) + } + diff --git a/R/gen_userp.R b/R/gen_userp.R new file mode 100644 index 0000000..235e783 --- /dev/null +++ b/R/gen_userp.R @@ -0,0 +1,491 @@ +#' @title Create a Wrapper To Be Used +#' in 'lavaan' Models +#' +#' @description Make a function on +#' `lavaan` object usable in a `lavaan` +#' model syntax. +#' +#' @details +#' +#' ## `gen_userp` +#' +#' There are cases in which we want to +#' create a user parameter which is a +#' function of other free parameters, +#' computed by a function. However such +#' a function may work only on a +#' `lavaan` object. +#' +#' If the target function works by +#' extracting parameter estimates stored +#' in the `Model` slot and/or the +#' `implied` slot, then [gen_userp()] +#' can be used to convert it to a +#' function that retrieves the parameter +#' estimates when being called by +#' [lavaan::lavaan()] or its wrappers, +#' modifies the stored +#' `lavaan` object using +#' [lavaan::lav_model_set_parameters()] +#' and [lavaan::lav_model_implied()] to +#' change the estimates, and call the +#' target function. +#' +#' Note that this is an unconventional +#' way to define a user parameter and +#' the generated function should always +#' be checked to see whether it works as +#' expected. +#' +#' As shown in the examples, the +#' parameter computed this may not +#' have standard error nor *p*-value. +#' +#' The main purpose is for the point +#' estimate, for searching the +#' likelihood-based confidence bound +#' using [ci_bound_ur()] and +#' [ci_bound_ur_i()]. +#' +#' Note that the target function +#' specified in `func` should work +#' directly on the parameter estimates +#' stored in the `Model` slot and then +#' get the estimates using +#' [lavaan::lav_model_get_parameters()]. +#' Functions that work on the unmodified +#' output generated by +#' [lavaan::lavaan()] usually do not +#' work. +#' +#' Users are not recommended to use +#' [gen_userp()] and [gen_sem_out_userp()] +#' directly because they require +#' unconventional way to extract +#' parameter estimates from a lavaan +#' model. However, developers may use +#' them to include functions +#' they wrote in a lavaan model. This +#' is the technique used by +#' [ci_bound_ur_i()] to constrain any +#' parameter in a model to an arbitrary +#' value. +#' +#' @param func A function that receives a +#' `lavaan`-object and returns a scalar. +#' See Details on the restriction on this +#' function. +#' +#' @param sem_out A `lavaan`-class object +#' to be modified. +#' +#' @return +#' +#' ## `gen_userp` +#' +#' It returns a function that +#' accepts a numeric vector of length +#' equals to the number of free parameters +#' in `sem_out`, and returns a scalar +#' which is the output of `func`. If this +#' vector is not supplied, it will try to +#' find it in the `parent.frame()`. This +#' is how it works inside a `lavaan` +#' model. +#' +#' @examples +#' +#' library(lavaan) +#' +#' data(simple_med) +#' dat <- simple_med +#' mod <- +#' " +#' m ~ a*x +#' y ~ b*m +#' ab := a*b +#' " +#' fit_med <- sem(mod, simple_med, fixed.x = FALSE) +#' parameterEstimates(fit_med) +#' +#' # A trivial example for verifying the results +#' my_ab <- function(object) { +#' # Need to use lav_model_get_parameters() +#' # because the object is only a modified +#' # lavaan-object, not one directly +#' # generated by lavaan function +#' est <- lavaan::lav_model_get_parameters(object@Model, type = "user") +#' unname(est[1] * est[2]) +#' } +#' +#' # Check the function +#' my_ab(fit_med) +#' coef(fit_med, type = "user")["ab"] +#' +#' # Create the function +#' my_userp <- gen_userp(func = my_ab, +#' sem_out = fit_med) +#' # Try it on the vector of free parameters +#' my_userp(coef(fit_med)) +#' +#' # Generate a modified lavaan model +#' fit_userp <- gen_sem_out_userp(userp = my_userp, +#' userp_name = "my_userp", +#' sem_out = fit_med) +#' +#' # This function can then be used in the model syntax. +#' +#' # Note that the following example only work when called inside the +#' # workspace or inside other functions such as ci_bound_ur()` +#' # and `ci_bound_ur_i()` because `lavaan::sem()` will +#' # search `my_userp()` in the global environment. +#' +#' # Therefore, the following lines are commented out. +#' # They should be run only in a "TRUE" interactive +#' # session. +#' +#' # mod2 <- +#' # " +#' # m ~ x +#' # y ~ m +#' # ab := my_userp() +#' # " +#' # fit_med2 <- sem(mod2, simple_med, fixed.x = FALSE) +#' # parameterEstimates(fit_med2) +#' # +#' # # Fit the model with the output of the function, a*b +#' # # fixed to .50 +#' # +#' # fit_new <- fit_userp(.50) +#' # +#' # # Check if the parameter ab is fixed to .50 +#' # parameterEstimates(fit_new) +#' +#' +#' @rdname gen_userp +#' +#' @export + +gen_userp <- function(func, + sem_out) { + stopifnot(is.function(func)) + stopifnot(inherits(sem_out, "lavaan")) + out <- function(.x., ...) { + if (missing(.x.)) { + .x. <- get(".x.", parent.frame()) + } + # Just in case ... + force(sem_out) + + sem_out@Model <- lavaan::lav_model_set_parameters(sem_out@Model, .x.) + sem_out@implied <- lavaan::lav_model_implied(sem_out@Model) + f0 <- func(sem_out) + if (is.na(f0) || is.nan(f0)) { + f0 <- Inf + } + f0 + } + out + } + +# @title Generate a Function to Fix a +# User-Defined Parameter +# +# @description Add a function to a +# `lavaan` object as a user-defined +# parameter and return a function which +# can fix this user-defined parameter +# to a value. +# +#' @details +#' +#' ## `gen_sem_out_userp` +#' +#' The function [gen_sem_out_userp()] +#' is to be used internally +#' for generating a function for searching +#' a likelihood-based confidence bound. +#' It is exported because it needs to +#' be run in an fresh external R process, +#' usually created by `callr` in other +#' internal functions. +#' +#' @param userp A function that is +#' generated by [gen_userp()]. +#' +#' @param userp_name The name of the +#' function `userp` to be used in the +#' `lavaan` model. It does not have to +#' be the name of the function in +#' `userp`. Should be changed only if it +#' conflicts with another object in the +#' parent environment, which should not +#' happen if the model is always fitted +#' in a clean R session. +#' +#' @param fix If `TRUE`, the default, +#' the function generated is used to fix +#' the value of `userp` to a target +#' value using an equality constraint. +#' If `FALSE`, then the function simply +#' fits the model to the data. +#' +#' @param control_args To be passed to +#' the argument of the same name in +#' [lavaan::lavaan()]. Default is +#' `list()`. Can be used to set the +#' default values of this argument +#' in the generated function. +#' +#' @param iter.max The maximum number of +#' iteration when the generated function +#' fit the model. Default is 10000. +#' +#' @param max_attempts If the initial +#' fit with the equality constraint +#' fails, how many more attempts +#' will be made by the generated +#' function. Default is 5. +#' +#' @return +#' +#' ## `gen_sem_out_userp` +#' +#' If `fix` is `TRUE`, it returns a +#' function with these arguments: +#' +#' - `target`: The value to which the +#' user-defined parameter will be fixed +#' to. +#' +#' - `verbose`: If `TRUE`, additional +#' information will be printed when +#' fitting the model. +#' +#' - `control`: The values to be passed +#' as a list to the argument of the +#' same name in [lavaan::lavaan()]. +#' +#' - `seed`: Numeric. If supplied, it +#' will be used in [set.seed()] to +#' initialize the random number +#' generator. Necessary to reproduce +#' some results because random numbers +#' are used in some steps in `lavaan`. +#' If `NULL`, the default, [set.seed()] +#' will not be called. +#' +#' If `fix` is `FALSE, then it returns a +#' function with optional arguments that +#' will be ignored, Calling it will +#' simply fit the modified model to the +#' data. Useful for getting the value of +#' the user-defined parameter. +#' +#' @rdname gen_userp +#' +#' @export + +gen_sem_out_userp <- function(userp, + sem_out, + userp_name = "semlbciuserp1234", + fix = TRUE, + control_args = list(), + iter.max = 10000, + max_attempts = 5) { + + # Generate a function that: + # - refits the model with the user-parameter fixed to a target value. + # CHECKED: Identical to the experimental version + # For add_fun() using callr. + + iter.max <- max(lavaan::lavInspect(sem_out, "optim")$iterations, + iter.max) + ptable <- lavaan::parameterTable(sem_out) + # Create a unique label: userp_label + labels <- unique(ptable$label) + labels <- labels[nchar(labels) > 0] + userp_label <- make.unique(c(labels, "user")) + userp_label <- userp_label[length(userp_label)] + # Generate the user parameter row using userp() + user_pt1 <- data.frame(lhs = userp_label[length(userp_label)], + op = ":=", + rhs = paste0(userp_name, "()"), + user = 1, + block = 0, + group = 0, + free = 0, + ustart = NA, + label = userp_label, + exo = 0) + user_pt1 <- lavaan::lav_partable_complete(partable = user_pt1, + start = TRUE) + # Add the user parameter to the model + ptable1 <- lavaan::lav_partable_add(partable = ptable, + add = user_pt1) + ptable1$se[which(ptable1$label == userp_label)] <- NA + + # Extract slots to be used + slot_opt <- sem_out@Options + slot_dat <- sem_out@Data + + if (fix) { + # Create the row of equality constraint + user_pt2 <- data.frame(lhs = userp_label, + op = "==", + rhs = NA, + user = 1, + block = 0, + group = 0, + free = 0, + ustart = NA, + exo = 0) + user_pt2 <- lavaan::lav_partable_complete(partable = user_pt2, + start = TRUE) + # Add the equality constraint to the parameter tahle + ptable2 <- lavaan::lav_partable_add(partable = ptable1, + add = user_pt2) + # Store the position of the constraint + i_eq <- which((ptable2$lhs == userp_label) & + (ptable2$op == "==")) + + # Fix the ids of the rows + ptable2$id <- seq_along(ptable2$id) + + # Set the options + slot_opt$do.fit <- TRUE + # The "mplus" version does not take into account the user-parameter + slot_opt$start <- "simple" + # Heywood cases are allowed + slot_opt$check.start <- FALSE + slot_opt$check.post <- FALSE + # Disable some options for efficiency + slot_opt$se <- "none" + # slot_opt$h1 <- FALSE + # slot_opt$baseline <- FALSE + + out <- function(target, + verbose = FALSE, + control = list(), + seed = NULL) { + if (!is.null(seed)) set.seed(seed) + # Just in case ... + force(control_args) + force(slot_opt) + force(slot_dat) + force(ptable2) + force(userp_label) + force(i_eq) + + ptable3 <- ptable2 + ptable3$rhs[i_eq] <- target + control_args <- utils::modifyList(control, + list(iter.max = iter.max)) + slot_opt$control <- control_args + + # Fit once to update the model slot + attempted <- 0 + fit_new_converged <- FALSE + if (verbose) { + cat("Target:", target, "\n") + } + if (verbose) { + cat("First attempt ...\n") + } + # Fail when + # eq.constraints.K is 0x0 + # Succeed when + # fit_new@Model@eq.constraints, or + # length(fit_new@Model@ceq.nonlinear.idx) > 0 + slot_opt_tmp <- slot_opt + slot_opt_tmp$do.fit <- FALSE + eq_not_ok <- TRUE + while (eq_not_ok) { + fit_new <- lavaan::lavaan(model = ptable3, + slotOptions = slot_opt_tmp, + slotData = slot_dat) + eq_not_ok <- !fit_new@Model@eq.constraints && + !(length(fit_new@Model@ceq.nonlinear.idx) > 0) + } + + # Extract the updated slots + slot_pat_new <- fit_new@ParTable + slot_model_new <- fit_new@Model + slot_smp_new <- fit_new@SampleStats + + # Fit the model with the user parameter constrained + # to the target value + fit_new <- tryCatch(lavaan::lavaan(slotParTable = slot_pat_new, + slotModel = slot_model_new, + slotSampleStats = slot_smp_new, + slotOptions = slot_opt, + slotData = slot_dat), + error = function(e) e) + if (inherits(fit_new, "error")) { + stop("Something's with the function.", + "Error occurred when fitting the modified model.") + } + + fit_new_converged <- lavaan::lavInspect(fit_new, "converged") + + # Try max_attempts more times + if (!fit_new_converged) { + if (verbose) { + cat("More attempts ...\n") + } + while (!((attempted > max_attempts) || + fit_new_converged)) { + if (verbose) {cat("Attempt", attempted, "\n")} + ptable_new <- fit_new@ParTable + i_free <- ptable_new$free > 0 + ptable_new$est[i_free] <- jitter(ptable_new$est[i_free], + factor = 50) + slot_opt$start <- as.data.frame(ptable_new) + fit_old <- fit_new + fit_new <- tryCatch(lavaan::lavaan(slotParTable = ptable_new, + slotModel = slot_model_new, + slotSampleStats = slot_smp_new, + slotOptions = slot_opt, + slotData = slot_dat), + error = function(e) e) + if (inherits(fit_new, "error")) { + fit_new <- fit_old + fit_new_converged <- FALSE + } else { + fit_new_converged <- lavaan::lavInspect(fit_new, "converged") + } + attempted <- attempted + 1 + } + } + if (!fit_new_converged) { + if (verbose) { + cat("The modified model failed to converge.") + } + } + if (verbose) { + cat("Done!\n") + } + fit_new + } + } else { + # The arguments are not used in this version + out <- function(target, + verbose = FALSE, + control = list(), + seed = NULL) { + # Generate a model with the user parameter added. + # For debugging + # Just in case ... + force(slot_opt) + force(slot_dat) + force(ptable1) + slot_opt$do.fit <- TRUE + fit_new <- lavaan::lavaan(model = ptable1, + slotOptions = slot_opt, + slotData = slot_dat) + fit_new + } + } + out + } + diff --git a/R/get_cibound.R b/R/get_cibound.R index cbd034a..2073b09 100644 --- a/R/get_cibound.R +++ b/R/get_cibound.R @@ -3,11 +3,27 @@ #' @description Get the `cibound` output of a bound from #' a `semlbci` object, the output of [semlbci()]. #' -#' @details It returns the original output of [ci_bound_wn_i()] -#' for a bound. Usually for diagnosis. +#' @details The function [get_cibound()] +#' returns the original output of +#' [ci_bound_wn_i()] for a bound. +#' Usually for diagnosis. #' -#' @return A `cibound`-class object. See [ci_bound_wn_i()] +#' The function [get_cibound_status_not_0()] +#' checks the status code of each bound, +#' and returns the `cibound` outputs of +#' bounds with status code not equal to +#' zero (i.e., something wrong in the +#' search). Printing it can print the +#' diagnostic information for all bounds +#' that failed in the search. +#' +#' @return +#' [get_cibound()] returns a `cibound`-class object. See [ci_bound_wn_i()] #' for details. +#' [get_cibound_status_not_0()] returns a list of +#' `cibound`-class objects with `status` not equal +#' to zero. If all bounds have `status` equal to +#' zero, it returns an empty list. #' #' @param x The output of [semlbci()]. #' @@ -16,7 +32,7 @@ #' may be omitted in the printout of `x`. #' #' @param which The bound for which the [ci_bound_wn_i()] is -#' to be extracted. Either `"lbound"`` or `"ubound"``. +#' to be extracted. Either `"lbound"` or `"ubound"`. #' #' @author Shu Fai Cheung #' @@ -46,6 +62,7 @@ #' # bound of the LBCI for the indirect effect: #' get_cibound(lbci_med, row_id = 6, which = "ubound") #' +#' @rdname get_cibound #' @export get_cibound <- function(x, @@ -63,3 +80,43 @@ get_cibound <- function(x, } cb_out } + +#' @rdname get_cibound +#' @export + +get_cibound_status_not_0 <- function(x) { + status_not_0_lb <- status_not_0(x, which = "lbound") + status_not_0_ub <- status_not_0(x, which = "ubound") + if (any(status_not_0_lb)) { + out_lb <- lapply(which(status_not_0_lb), + get_cibound, + x = x, + which = "lbound") + } else { + out_lb <- list() + } + if (any(status_not_0_ub)) { + out_ub <- lapply(which(status_not_0_ub), + get_cibound, + x = x, + which = "ubound") + } else { + out_ub <- list() + } + c(out_lb, out_ub) + } + +#' @noRd + +status_not_0 <- function(x, which = c("lbound", "ubound")) { + which <- match.arg(which) + x_diag <- attr(x, switch(which, + lbound = "lb_diag", + ubound = "ub_diag")) + i <- sapply(x_diag, function(x) { + ifelse(length(x) > 1, + x$status != 0, + FALSE) + }) + i + } \ No newline at end of file diff --git a/R/loglike_at_user.R b/R/loglike_at_user.R new file mode 100644 index 0000000..ee47262 --- /dev/null +++ b/R/loglike_at_user.R @@ -0,0 +1,32 @@ +#' @noRd +# TODO: +# - To be exported for plotting loglikelihood + +loglik_user <- function(x, + sem_out_userp, + sem_out, + lrt_method = "default", + ...) { + sem_out_userp_tmp <- sem_out_userp_run(target = x, + object = sem_out_userp, + ...) + lrt_x <- tryCatch(lavaan::lavTestLRT(sem_out_userp_tmp, + sem_out, + method = lrt_method), + error = function(e) e, + warning = function(w) w) + if (inherits(lrt_x, "warning")) { + tmp <- as.character(lrt_x) + if (grepl("scaling factor is negative", + tmp, fixed = TRUE)) { + lrt_x <- tryCatch(lavaan::lavTestLRT(sem_out_userp_tmp, + sem_out, + method = "satorra.2000"), + error = function(e) e, + warning = function(w) w) + } + } + out <- lrt_x[2, "Chisq diff"] / (-2) + attr(out, "sem_out_userp_x") <- sem_out_userp_tmp + out + } \ No newline at end of file diff --git a/R/print_cibound.R b/R/print_cibound.R index 062fda7..3589f2f 100644 --- a/R/print_cibound.R +++ b/R/print_cibound.R @@ -47,7 +47,7 @@ print.cibound <- function(x, digits = 5, ...) { call_org <- x$call out_diag <- x$diag - ci_method <- switch(out_diag$method, wn = "Wu-Neale-2012") + ci_method <- x$method_name if (is.null(out_diag$standardized)) { std <- "No" } else { @@ -56,37 +56,49 @@ print.cibound <- function(x, digits = 5, ...) { lor <- out_diag$i_lor lor2 <- paste0(lor$lhs, " ", lor$op, " ", lor$rhs, " (group = ", lor$group, ", block = ", lor$block, ")") - cat(paste0("Target Parameter:\t", lor2)) - cat(paste0("\nPosition:\t\t", out_diag$i)) - cat(paste0("\nWhich Bound:\t\t", switch(out_diag$which, + cat(paste0( "Target Parameter: ", lor2)) + cat(paste0("\nPosition: ", out_diag$i)) + cat(paste0("\nWhich Bound: ", switch(out_diag$which, lbound = "Lower Bound", ubound = "Upper Bound"))) - cat(paste0("\nMethod:\t\t\t", ci_method)) - cat(paste0("\nConfidence Level:\t", out_diag$ciperc)) - cat(paste0("\nAchieved Level:\t\t", out_diag$ciperc_final)) - cat(paste0("\nStandardized:\t\t", std)) - cat(paste0("\nLikelihood-Based Bound:\t", + cat(paste0("\nMethod: ", ci_method)) + cat(paste0("\nConfidence Level: ", out_diag$ciperc)) + if (!is.null(out_diag$search_error)) { + cat(paste0("\nSearch Error Message:\n", out_diag$search_error, "\n")) + cat(paste0(rep("*", round(getOption("width") * .8)), collapse = "")) + cat("\n") + cat(strwrap(c("The search was terminated due to the error above.", + "The following lines should not be interpreted.")), + sep = "\n") + cat(paste0(rep("*", round(getOption("width") * .8)), collapse = "")) + cat("\n") + } + cat(paste0("\nAchieved Level: ", out_diag$ciperc_final)) + cat(paste0("\nStandardized: ", std)) + cat(paste0("\nLikelihood-Based Bound: ", ifelse(is.na(x$bound), "Not valid", round(x$bound, digits)))) - cat(paste0("\nWald Bound:\t\t", round(out_diag$ci_org_limit, digits))) - cat(paste0("\nPoint Estimate:\t\t", round(out_diag$est_org, digits))) - cat(paste0("\nRatio to Wald Bound:\t", ifelse(is.na(x$bound), "Not valid", + cat(paste0("\nWald Bound: ", round(out_diag$ci_org_limit, digits))) + cat(paste0("\nPoint Estimate: ", round(out_diag$est_org, digits))) + cat(paste0("\nRatio to Wald Bound: ", ifelse(is.na(x$bound), "Not valid", round(out_diag$ci_limit_ratio, digits)))) cat(paste0("\n\n-- Check --")) # Get p_tol in call p_tol_call <- call_org$p_tol if (is.null(p_tol_call)) { ci_method <- as.character(call_org[[1]]) + # ci_method can be of the form xxx::yyy + ci_method <- ci_method[length(ci_method)] p_tol_call <- formals(ci_method)$p_tol } ciperc_diff <- abs(out_diag$ciperc - out_diag$ciperc_final) - cat(paste0("\nLevel achieved?\t\t", + cat(paste0("\nLevel achieved? ", ifelse(ciperc_diff <= p_tol_call, "Yes", "No"), " (Difference: ", formatC(ciperc_diff, digits = digits), ";", " Tolerance: ", p_tol_call, ")")) - cat(paste0("\nSolution admissible?\t", + cat(paste0("\nSolution admissible? ", ifelse(out_diag$fit_post_check, "Yes", "No"))) bound_unchecked <- out_diag$bound_unchecked if (is.na(bound_unchecked)) { @@ -99,13 +111,16 @@ print.cibound <- function(x, digits = 5, ...) { "Yes", "No") ) } - cat(paste0("\nDirection valid?\t", direct_valid)) + cat(paste0("\nDirection valid? ", direct_valid)) cat("\n") cat("\n-- Optimization Information --") - cat(paste0("\nSolver Status:\t\t", out_diag$optim_status)) - cat(paste0("\nConvergence Message:\t", out_diag$optim_message)) - cat(paste0("\nIterations:\t\t", out_diag$optim_iterations)) + cat(paste0("\nSolver Status: ", out_diag$optim_status)) + cat(paste0("\nConvergence Message: ", out_diag$optim_message)) + if (any(grepl("NLOPT_MAXTIME_REACHED: ", out_diag$optim_message))) { + cat("\n - Set 'timeout' to a larger value to increase maximum time.") + } + cat(paste0("\nIterations: ", out_diag$optim_iterations)) t_cond <- out_diag$optim_termination_conditions t_cond2 <- strsplit(t_cond, "\t")[[1]] t_cond2 <- gsub("maxeval: ", "maxeval: ", t_cond2) @@ -122,9 +137,9 @@ print.cibound <- function(x, digits = 5, ...) { cat("\n-- Parameter Estimates --\n") print(coef_c) - cat(paste0("\nBound before check:\t", + cat(paste0("\nBound before check: ", round(out_diag$bound_unchecked, digits))) - cat(paste0("\nStatus Code:\t\t", out_diag$status)) + cat(paste0("\nStatus Code: ", out_diag$status)) cat("\nCall: ") call_print <- call_org if (!is.name(call_print$f_constr)) { diff --git a/R/print_semlbci.R b/R/print_semlbci.R index bd5ce9c..61b4693 100644 --- a/R/print_semlbci.R +++ b/R/print_semlbci.R @@ -195,21 +195,21 @@ print.semlbci <- function(x, out$label <- NULL } class(out) <- "data.frame" - out$lbci_lb <- formatC(out$lbci_lb, digits, format = "f") + out$lbci_lb <- formatC(out$lbci_lb, digits, format = "f", mode = "double") if ("est.std" %in% colnames(x)) { standardized <- TRUE - out$est.std <- formatC(out$est.std, digits, format = "f") + out$est.std <- formatC(out$est.std, digits, format = "f", mode = "double") } else { standardized <- FALSE - out$est <- formatC(out$est, digits, format = "f") + out$est <- formatC(out$est, digits, format = "f", mode = "double") } - out$lbci_ub <- formatC(out$lbci_ub, digits, format = "f") - out$ci_org_lb <- formatC(out$ci_org_lb, digits, format = "f") - out$ci_org_ub <- formatC(out$ci_org_ub, digits, format = "f") - out$ratio_lb <- formatC(out$ratio_lb, digits, format = "f") - out$ratio_ub <- formatC(out$ratio_ub, digits, format = "f") - out$cl_lb <- formatC(out$cl_lb, digits, format = "f") - out$cl_ub <- formatC(out$cl_ub, digits, format = "f") + out$lbci_ub <- formatC(out$lbci_ub, digits, format = "f", mode = "double") + out$ci_org_lb <- formatC(out$ci_org_lb, digits, format = "f", mode = "double") + out$ci_org_ub <- formatC(out$ci_org_ub, digits, format = "f", mode = "double") + out$ratio_lb <- formatC(out$ratio_lb, digits, format = "f", mode = "double") + out$ratio_ub <- formatC(out$ratio_ub, digits, format = "f", mode = "double") + out$cl_lb <- formatC(out$cl_lb, digits, format = "f", mode = "double") + out$cl_ub <- formatC(out$cl_ub, digits, format = "f", mode = "double") call_org <- attr(out, "call") if (drop_no_lbci) { out <- out[!is.na(out$status_lb), ] @@ -298,7 +298,9 @@ print.semlbci <- function(x, msg <- c(msg, paste0("* ok_l, ok_u: Whether the search encountered any problem. ", "If no problem encountered, it is equal to 0. Any value ", - "other than 0 indicates something was wrong in the search.")) + "other than 0 indicates something was wrong in the search. ", + "Try running 'semlbci()' again, setting 'semlbci_out' to this output, ", + "and set 'try_k_more_times' to a positive number, e.g., 2 to 5.")) } if (standardized) { msg <- c(msg, @@ -562,7 +564,9 @@ print_semlbci_text <- function(x, tmp0 <- "- Status: " tmp <- paste0(tmp0, "Whether the search encountered any problem for the two bounds. ", "If no problem encountered, the code is 0 (displayed as 'OK'). Any value ", - "other than 0 indicates something was wrong in the search.") + "other than 0 indicates something was wrong in the search. ", + "Try running 'semlbci()' again, setting 'semlbci_out' to this output, ", + "and set 'try_k_more_times' to a positive number, e.g., 2 to 5.") tmp <- strwrap(tmp, indent = 0, exdent = nchar(tmp0)) @@ -699,7 +703,8 @@ format_ratio <- function(object, where <- match.arg(where) tmp <- formatC(as.numeric(object), digits = digits, - format = "f") + format = "f", + mode = "double") tmp2 <- tmp if (where == "left") { tmp2[which(!flag)] <- paste0(" ", tmp[which(!flag)]) diff --git a/R/ram_to_lav_mod.R b/R/ram_to_lav_mod.R index 24211be..3b902ce 100644 --- a/R/ram_to_lav_mod.R +++ b/R/ram_to_lav_mod.R @@ -20,15 +20,14 @@ ram_to_lav_mod <- function(ram, lav_mod, standardized = FALSE) { ov_names <- rownames(lav_mod$theta) lv_names <- rownames(lav_mod$psi) lv_names <- lv_names[!(lv_names %in% ov_names)] - mA <- ram$A if (!is.null(lav_mod$lambda)) { # A to lambda - lav_mod$lambda[] <- 0 + lav_mod$lambda[is.numeric(lav_mod$lambda)] <- 0 lambda_rnames <- rownames(lav_mod$lambda) lambda_cnames <- colnames(lav_mod$lambda) - lav_mod$lambda <- mA[lambda_rnames, lambda_cnames] + lav_mod$lambda <- mA[lambda_rnames, lambda_cnames, drop = FALSE] lambda1 <- lambda_rnames[lambda_rnames %in% lambda_cnames] lav_mod$lambda[lambda1, ] <- 0 lav_mod$lambda[lambda1, lambda1] <- 1 @@ -36,18 +35,18 @@ ram_to_lav_mod <- function(ram, lav_mod, standardized = FALSE) { if (!is.null(lav_mod$beta)) { # A to beta - lav_mod$beta[] <- 0 + lav_mod$beta[is.numeric(lav_mod$beta)] <- 0 beta_names <- rownames(lav_mod$beta) - lav_mod$beta <- mA[beta_names, beta_names] + lav_mod$beta <- mA[beta_names, beta_names, drop = FALSE] } mS <- ram$S if (!is.null(lav_mod$theta)) { # S to theta - lav_mod$theta[] <- 0 + lav_mod$theta[is.numeric(lav_mod$theta)] <- 0 theta_names <- rownames(lav_mod$theta) - lav_mod$theta <- mS[theta_names, theta_names] + lav_mod$theta <- mS[theta_names, theta_names, drop = FALSE] if (standardized) { tmp <- diag(lav_mod$theta) lav_mod$theta <- stats::cov2cor(lav_mod$theta) @@ -57,9 +56,9 @@ ram_to_lav_mod <- function(ram, lav_mod, standardized = FALSE) { } if (!is.null(lav_mod$psi)) { # S to psi - lav_mod$psi[] <- 0 + lav_mod$psi[is.numeric(lav_mod$psi)] <- 0 psi_names <- rownames(lav_mod$psi) - lav_mod$psi <- mS[psi_names, psi_names] + lav_mod$psi <- mS[psi_names, psi_names, drop = FALSE] if (standardized) { tmp <- diag(lav_mod$psi) lav_mod$psi <- stats::cov2cor(lav_mod$psi) @@ -69,7 +68,7 @@ ram_to_lav_mod <- function(ram, lav_mod, standardized = FALSE) { if (!is.null(lav_mod$nu)) { # M to nu - lav_mod$nu[] <- 0 + lav_mod$nu[is.numeric(lav_mod$nu)] <- 0 nu_names <- rownames(lav_mod$nu) lav_mod$nu[] <- t(ram$M[, nu_names, drop = FALSE]) lav_mod$nu[lambda1] <- 0 @@ -77,9 +76,9 @@ ram_to_lav_mod <- function(ram, lav_mod, standardized = FALSE) { if (!is.null(lav_mod$alpha)) { # M to alpha - lav_mod$alpha[] <- 0 + lav_mod$alpha[is.numeric(lav_mod$alpha)] <- 0 alpha_names <- rownames(lav_mod$alpha) - lav_mod$alpha[] <- t(ram$M[, alpha_names, drop = FALSE]) + lav_mod$alpha[is.numeric(lav_mod$alpha)] <- t(ram$M[, alpha_names, drop = FALSE]) } lav_mod diff --git a/R/semlbci.R b/R/semlbci.R index effae88..f52da08 100644 --- a/R/semlbci.R +++ b/R/semlbci.R @@ -71,8 +71,8 @@ #' @param standardized If `TRUE`, the LBCI is for the standardized estimates. #' #' @param method The method to be used to search for the confidence -#' bounds. Currently only `"wn"` (Wu-Neale-2012), the default, is -#' supported. +#' 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()] @@ -106,6 +106,11 @@ #' 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 #' #' @references @@ -168,15 +173,16 @@ semlbci <- function(sem_out, remove_intercepts = TRUE, ciperc = .95, standardized = FALSE, - method = "wn", + method = c("wn", "ur"), robust = c("none", "satorra.2000"), - try_k_more_times = 2, + try_k_more_times = 0, semlbci_out = NULL, check_fit = TRUE, ..., parallel = FALSE, ncpus = 2, - use_pbapply = TRUE) { + use_pbapply = TRUE, + loadbalancing = TRUE) { if (!inherits(sem_out, "lavaan")) { stop("sem_out is not a supported object.") } @@ -189,6 +195,20 @@ semlbci <- function(sem_out, 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) @@ -208,8 +228,12 @@ semlbci <- function(sem_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("semblci_out and the parameter table of sem_out do not match.") + 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) & @@ -253,7 +277,7 @@ semlbci <- function(sem_out, pars <- remove_variances(pars, sem_out) } if (remove_intercepts) { - pars <- remove_variances(pars, sem_out) + pars <- remove_intercepts(pars, sem_out) } # Remove parameters already with LBCIs in semlbci_out. pars <- pars[!pars %in% i_id[pars_lbci_yes]] @@ -303,6 +327,11 @@ semlbci <- function(sem_out, 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, @@ -414,6 +443,7 @@ semlbci <- function(sem_out, 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")] } diff --git a/R/set_constraint.R b/R/set_constraint.R index d4e499f..20cd658 100644 --- a/R/set_constraint.R +++ b/R/set_constraint.R @@ -81,7 +81,8 @@ set_constraint <- function(sem_out, ciperc = .95) { debug = FALSE, lav_warn = FALSE, sf = 1, - sf2 = 0) { + sf2 = 0, + stop_on_error = FALSE) { target <- fmin + sf * (qcrit - sf2) / (2 * n) if (debug) { cat(ls()) @@ -91,31 +92,49 @@ set_constraint <- function(sem_out, ciperc = .95) { eq_out <- sem_out@Model@ceq.function(param) eq_jac <- sem_out@Model@con.jac if (lav_warn) { - fit2 <- lavaan::lavaan( + fit2 <- tryCatch(lavaan::lavaan( slotOptions = slot_opt3, slotParTable = slot_pat2, slotModel = slot_mod3, slotSampleStats = slot_smp2, - slotData = slot_dat2) + slotData = slot_dat2), + error = function(e) e) } else { - suppressWarnings(fit2 <- lavaan::lavaan( + fit2 <- tryCatch(suppressWarnings(lavaan::lavaan( slotOptions = slot_opt3, slotParTable = slot_pat2, slotModel = slot_mod3, slotSampleStats = slot_smp2, - slotData = slot_dat2)) - } - if (lav_warn) { - fit2_gradient <- rbind(lavaan::lavTech(fit2, "gradient")) - fit2_jacobian <- rbind(eq_jac, - lavaan::lavTech(fit2, "gradient")) - } else { - suppressWarnings(fit2_gradient <- - rbind(lavaan::lavTech(fit2, "gradient"))) - suppressWarnings(fit2_jacobian <- - rbind(eq_jac, - lavaan::lavTech(fit2, "gradient"))) + slotData = slot_dat2)), + error = function(e) e) } + is_error <- inherits(fit2, "error") + if (is_error) { + if (stop_on_error) { + stop("Error in evaluating the constraint(s).") + } + fit2_gradient <- rbind(rep(-Inf, length(param))) + fit2_jacobian <- rbind(eq_jac, + rep(-Inf, length(param))) + out <- list(objective = Inf, + gradient = fit2_gradient, + constraints = Inf, + jacobian = fit2_jacobian, + parameterTable = NA) + return(out) + } else { + if (lav_warn) { + fit2_gradient <- rbind(lavaan::lavTech(fit2, "gradient")) + fit2_jacobian <- rbind(eq_jac, + lavaan::lavTech(fit2, "gradient")) + } else { + fit2_gradient <- tryCatch(suppressWarnings(rbind(lavaan::lavTech(fit2, "gradient"))), + error = function(e) rep(-Inf, length(param))) + suppressWarnings(fit2_jacobian <- + rbind(eq_jac, + fit2_gradient)) + } + } list( objective = lavaan::lavTech(fit2, "optim")$fx, gradient = fit2_gradient, @@ -131,7 +150,8 @@ set_constraint <- function(sem_out, ciperc = .95) { debug = FALSE, lav_warn = FALSE, sf = 1, - sf2 = 0) { + sf2 = 0, + stop_on_error = FALSE) { target <- fmin + sf * (qcrit - sf2) / (2 * n) if (debug) { cat(ls()) @@ -139,29 +159,46 @@ set_constraint <- function(sem_out, ciperc = .95) { } slot_mod3 <- lavaan::lav_model_set_parameters(slot_mod2, param) if (lav_warn) { - fit2 <- lavaan::lavaan( + fit2 <- tryCatch(lavaan::lavaan( slotOptions = slot_opt3, slotParTable = slot_pat2, slotModel = slot_mod3, slotSampleStats = slot_smp2, - slotData = slot_dat2) + slotData = slot_dat2), + error = function(e) e) } else { - suppressWarnings(fit2 <- lavaan::lavaan( + fit2 <- tryCatch(suppressWarnings(lavaan::lavaan( slotOptions = slot_opt3, slotParTable = slot_pat2, slotModel = slot_mod3, slotSampleStats = slot_smp2, - slotData = slot_dat2)) - } - if (lav_warn) { - fit2_gradient <- rbind(lavaan::lavTech(fit2, "gradient")) - fit2_jacobian <- rbind(lavaan::lavTech(fit2, "gradient")) - } else { - suppressWarnings(fit2_gradient <- - rbind(lavaan::lavTech(fit2, "gradient"))) - suppressWarnings(fit2_jacobian <- - rbind(lavaan::lavTech(fit2, "gradient"))) + slotData = slot_dat2)), + error = function(e) e) } + is_error <- inherits(fit2, "error") + if (is_error) { + if (stop_on_error) { + stop("Error in evaluating the constraint(s).") + } + fit2_gradient <- rbind(rep(-Inf, length(param))) + fit2_jacobian <- rbind(rep(-Inf, length(param))) + out <- list(objective = Inf, + gradient = fit2_gradient, + constraints = Inf, + jacobian = fit2_jacobian, + parameterTable = NA) + return(out) + } else { + if (lav_warn) { + fit2_gradient <- rbind(lavaan::lavTech(fit2, "gradient")) + fit2_jacobian <- rbind(lavaan::lavTech(fit2, "gradient")) + } else { + fit2_gradient <- tryCatch(suppressWarnings(rbind(lavaan::lavTech(fit2, "gradient"))), + error = function(e) rep(-Inf, length(param))) + suppressWarnings(fit2_jacobian <- + rbind(fit2_gradient)) + } + } list( objective = lavaan::lavTech(fit2, "optim")$fx, gradient = fit2_gradient, diff --git a/README.md b/README.md index 0ed490a..4833be4 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,7 @@ [![DOI](https://img.shields.io/badge/doi-10.1080/10705511.2023.2183860-blue.svg)](https://doi.org/10.1080/10705511.2023.2183860) -(Version 0.10.4, updated on 2023-10-31, [release history](https://sfcheung.github.io/semlbci/news/index.html)) +(Version 0.10.4.24, updated on 2024-06-01 [release history](https://sfcheung.github.io/semlbci/news/index.html)) # semlbci diff --git a/_pkgdown.yml b/_pkgdown.yml index 87e676c..debf77a 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -66,9 +66,13 @@ reference: desc: Low level functions for advanced users. - contents: - ci_bound_wn_i + - ci_bound_ur_i + - ci_bound_ur - ci_i_one - set_constraint - print.cibound + - gen_sem_out_userp + - gen_userp - title: Datasets desc: Datasets used in examples. - contents: diff --git a/man/ci_bound_ur.Rd b/man/ci_bound_ur.Rd new file mode 100644 index 0000000..0590572 --- /dev/null +++ b/man/ci_bound_ur.Rd @@ -0,0 +1,235 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ci_bound_ur_i.R +\name{ci_bound_ur} +\alias{ci_bound_ur} +\alias{gen_est_i} +\title{Find a Likelihood-Based +Confidence Bound By Root Finding} +\usage{ +ci_bound_ur( + sem_out, + func, + ..., + level = 0.95, + which = c("lbound", "ubound"), + interval = NULL, + progress = FALSE, + method = "uniroot", + lrt_method = "default", + tol = 5e-04, + root_target = c("chisq", "pvalue"), + d = 5, + uniroot_extendInt = "yes", + uniroot_trace = 0, + uniroot_maxiter = 1000, + use_callr = TRUE, + rs = NULL +) + +gen_est_i(i, sem_out, standardized = FALSE) +} +\arguments{ +\item{sem_out}{The fit object. +Currently supports +\link[lavaan:lavaan-class]{lavaan::lavaan} objects only.} + +\item{func}{A function that receives +a lavaan object and returns a scalar. +This function is to be used by +\code{\link[=gen_userp]{gen_userp()}} and so there are +special requirements on it. +Alternatively, it can be the output +of \code{\link[=gen_est_i]{gen_est_i()}}.} + +\item{...}{Optional arguments to be +passed to \code{func}. Usually not used +but included in case the function +has such arguments.} + +\item{level}{The level of confidence +of the confidence interval. Default +is .95, or 95\%.} + +\item{which}{Whether the lower bound +or the upper bound is to be found. +Must be \code{"lbound"} or \code{"ubound"}.} + +\item{interval}{A numeric vector of +two values, which is the initial +interval to be searched. If \code{NULL}, +the default, it will be determined +internally using Wald or delta +method confidence interval, if +available.} + +\item{progress}{Whether progress will +be reported on screen during the +search. Default is +\code{FALSE}.} + +\item{method}{The actual function to +be used in the search. which can only +be \code{"uniroot"}, the default, for now. +May include other function in the +future.} + +\item{lrt_method}{The method used in +\code{\link[lavaan:lavTestLRT]{lavaan::lavTestLRT()}}. Default is +\code{"default"}. It is automatically set +to \code{"satorra.2000"} and cannot be +overridden if a scaled test statistic +is requested in \code{sem_out}.} + +\item{tol}{The tolerance used in +\code{\link[=uniroot]{uniroot()}}, default is .005.} + +\item{root_target}{Whether the +chi-square difference (\code{"chisq"}), +the default, or its \emph{p}-value +(\code{"pvalue"}) is used as the function +value in finding the root. Should have +little impact on the results.} + +\item{d}{A value used to determine +the width of the interval in the +initial search. Larger this value, +\emph{narrow} the interval. Default is 5.} + +\item{uniroot_extendInt}{To be passed +to the argument \code{extendInt} of +\code{\link[=uniroot]{uniroot()}}. Whether the interval +should be extended if the root is not +found. Default is \code{"yes"}. Refer to +the help page of \code{\link[=uniroot]{uniroot()}} for +possible values.} + +\item{uniroot_trace}{To be passed to +the argument \code{trace} of \code{\link[=uniroot]{uniroot()}}. +How much information is printed +during the search. Default is 0, and +no information is printed during the +search. Refer to the help page of +\code{\link[=uniroot]{uniroot()}} for possible values.} + +\item{uniroot_maxiter}{The maximum +number of iteration in the search. +Default is 1000.} + +\item{use_callr}{Whether the +\code{callr} package will be used to +do the search in a separate R +process. Default is \code{TRUE}. Should +not set to \code{FALSE} if used in an +interactive environment unless this is +intentional.} + +\item{rs}{Optional. If set to +a persistent R process created by +\code{callr}, it will be used instead of +starting a new one, and it will not +be terminated on exit.} + +\item{i}{The position of the target +parameter as appeared in the +parameter table of an lavaan object, +generated by +\code{\link[lavaan:lavParTable]{lavaan::parameterTable()}}.} + +\item{standardized}{If \code{TRUE}, the +standardized estimate is to be +retrieved. Default is \code{FALSE}. +Only support \code{"std.all"} for now.} +} +\value{ +The function \code{\link[=ci_bound_ur]{ci_bound_ur()}} returns +a list with the following elements: +\itemize{ +\item \code{bound}: The bound found. +\item \code{optimize_out}: THe output of the +root finding function, \code{\link[=uniroot]{uniroot()}} +for now. (Called \code{optimize_out} +because an earlier version of this +function also uses \code{\link[=optimize]{optimize()}}). +\item \code{sem_out_bound}: The \code{lavaan} model +with the user-defined parameter fixed +to the bound. +\item \code{lrt}: The output of +\code{\link[lavaan:lavTestLRT]{lavaan::lavTestLRT()}} comparing +\code{sem_out} and \code{sem_out_bound}. +\item \code{bound_start}: The Wald or delta +method confidence bound returned when +determining the interval internally. +\item \code{user_est}: The estimate of the +user-defined parameter when +determining the interval internally. +} + +The function \code{\link[=gen_est_i]{gen_est_i()}} returns +a special function can inspects the +\code{Model} slot (and \code{implied} slot +if necessary) of a modified \code{lavaan} +object and return the parameter +estimate. This function is to be used +by \code{\link[=ci_bound_ur]{ci_bound_ur()}} or +\code{\link[=gen_sem_out_userp]{gen_sem_out_userp()}}. +} +\description{ +Find the lower or upper +bound of the likelihood-based +confidence interval (LBCI) for one +parameter in a structural equation +model fitted in \code{\link[lavaan:lavaan]{lavaan::lavaan()}} +using \code{\link[=uniroot]{uniroot()}}. +} +\details{ +This function is called xby +\code{\link[=ci_bound_ur_i]{ci_bound_ur_i()}}. This function is +exported because it is a +stand-alone function that can be used +directly for any function that +receives a lavaan object and returns +a scalar. + +The function \code{\link[=ci_bound_ur_i]{ci_bound_ur_i()}} is a +wrapper of this function, with an +interface similar to that of +\code{\link[=ci_bound_wn_i]{ci_bound_wn_i()}} and returns a +\code{cibound}-class object. The +user-parameter function is generated +internally by \code{\link[=ci_bound_wn_i]{ci_bound_wn_i()}}. + +This function, on the other hand, +requires users to supply the function +directly through the \code{func} argument. +This provides the flexibility to find +the bound for any function of the +model parameter, even one that cannot +be easily coded in \code{lavaan} model +syntax. +} +\examples{ + +library(lavaan) +data(simple_med) +dat <- simple_med +mod <- +" +m ~ x +y ~ m +" +fit_med <- lavaan::sem(mod, simple_med, fixed.x = FALSE) +parameterTable(fit_med) +# Create a function to get the second parameter +est_i <- gen_est_i(i = 2, sem_out = fit_med) +# Find the lower bound of the likelihood-based confidence interval +# of the second parameter. +# user_callr should be TRUE or omitted in read research. +# Remove interval in read research. It is added to speed up the example. +out1l <- ci_bound_ur(sem_out = fit_med, + func = est_i, + which = "lbound", + use_callr = FALSE, + interval = c(.39070, .39075)) +out1l + +} diff --git a/man/ci_bound_ur_i.Rd b/man/ci_bound_ur_i.Rd new file mode 100644 index 0000000..ec98322 --- /dev/null +++ b/man/ci_bound_ur_i.Rd @@ -0,0 +1,289 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ci_bound_ur_i.R +\name{ci_bound_ur_i} +\alias{ci_bound_ur_i} +\title{Likelihood-Based Confidence +Bound By Root Finding} +\usage{ +ci_bound_ur_i( + i = NULL, + npar = NULL, + sem_out = NULL, + f_constr = NULL, + which = NULL, + history = FALSE, + perturbation_factor = 0.9, + lb_var = -Inf, + standardized = FALSE, + wald_ci_start = !standardized, + opts = list(), + ciperc = 0.95, + ci_limit_ratio_tol = 1.5, + verbose = FALSE, + sf = 1, + sf2 = 0, + p_tol = 5e-04, + std_method = "internal", + bounds = "none", + xtol_rel_factor = 1, + ftol_rel_factor = 1, + lb_prop = 0.05, + lb_se_k = 3, + d = 5, + ... +) +} +\arguments{ +\item{i}{The position of the target +parameter as appeared in the +parameter table of a lavaan object, +generated by +\code{\link[lavaan:lavParTable]{lavaan::parameterTable()}}.} + +\item{npar}{Ignored by this function. +Included consistency in the +interface.} + +\item{sem_out}{The fit object. +Currently supports +\link[lavaan:lavaan-class]{lavaan::lavaan} objects only.} + +\item{f_constr}{Ignored by this +function. Included consistency in the +interface.} + +\item{which}{Whether the lower bound +or the upper bound is to be found. +Must be \code{"lbound"} or \code{"ubound"}.} + +\item{history}{Not used. Kept for +backward compatibility.} + +\item{perturbation_factor}{Ignored by +this function. Included consistency +in the interface.} + +\item{lb_var}{Ignored by this +function. Included consistency in the +interface.} + +\item{standardized}{If \code{TRUE}, the +LBCI is for the requested estimate in +the standardized solution. Default is +\code{FALSE}.} + +\item{wald_ci_start}{Ignored by this +function. Included consistency in the +interface.} + +\item{opts}{Options to be passed to +\code{\link[stats:uniroot]{stats::uniroot()}}. Default is +\code{list()}.} + +\item{ciperc}{The intended coverage +probability for the confidence +interval. Default is .95, and the +bound for a 95\% confidence interval +will be sought.} + +\item{ci_limit_ratio_tol}{The +tolerance for the ratio of \code{a} to +\code{b}, where \code{a} is the distance +between an LBCI limit and the point +estimate, and the \code{b} is the distance +between the original confidence limit +(by default the Wald CI in +\code{\link[lavaan:lavaan]{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.} + +\item{verbose}{If \code{TRUE}, the +function will store more diagnostic +information in the attribute \code{diag}. +Default is \code{FALSE}.} + +\item{sf}{Ignored by this function. +Included consistency in the +interface.} + +\item{sf2}{Ignored by this function. +Included consistency in the +interface.} + +\item{p_tol}{Tolerance for checking +the achieved level of confidence. If +the absolute difference between the +achieved level and \code{ciperc} is +greater than this amount, a warning +is set in the status code and the +bound is set to \code{NA}. Default is +5e-4.} + +\item{std_method}{The method used to +find the standardized solution. If +equal to \code{"lavaan"}, +\code{\link[lavaan:standardizedSolution]{lavaan::standardizedSolution()}} will +be used. If equal to \code{"internal"}, an +internal function will be used. The +\code{"lavaan"} method should work in all +situations, but the \code{"internal"} +method is usually much faster. +Default is \code{"internal"}.} + +\item{bounds}{Ignored by this +function. Included consistency in the +interface.} + +\item{xtol_rel_factor}{Ignored by +this function. Included consistency +in the interface.} + +\item{ftol_rel_factor}{Ignored by +this function. Included consistency +in the interface.} + +\item{lb_prop}{Ignored by this +function. Included consistency in the +interface.} + +\item{lb_se_k}{Ignored by this +function. Included consistency in the +interface.} + +\item{d}{A value used to determine +the width of the interval in the +initial search. Larger this value, +\emph{narrow} the interval. Default is 5. +Used by \code{\link[=ci_bound_ur]{ci_bound_ur()}}.} + +\item{...}{Optional arguments. Not +used.} +} +\value{ +A \code{cibound}-class object +which is a list with three elements: +\itemize{ +\item \code{bound}: A single number. The value +of the bound located. \code{NA} is the +search failed for various reasons. +\item \code{diag}: A list of diagnostic +information. +\item \code{call}: The original call. +} + +A detailed and organized output can +be printed by the default print +method (\code{\link[=print.cibound]{print.cibound()}}). +} +\description{ +Using root finding +to find the lower or +upper bound of the likelihood-based +confidence interval (LBCI) for one +parameter in a structural equation +model fitted in \code{\link[lavaan:lavaan]{lavaan::lavaan()}}. +} +\details{ +\subsection{Important Notice}{ + +This function is not supposed to be +used directly by users in typical +scenarios. Its interface is +user-\emph{unfriendly} because it should +be used through \code{\link[=semlbci]{semlbci()}}. It is +exported such that interested users +can examine how a confidence bound is +found, or use it for experiments or +simulations. +} + +\subsection{Usage}{ + +This function is the lowest level +function used by \code{\link[=semlbci]{semlbci()}}. +\code{\link[=semlbci]{semlbci()}} calls this function once +for each bound of each parameter. + +For consistency in the interface, +most of the arguments in +\code{\link[=ci_bound_wn_i]{ci_bound_wn_i()}} are also included +in this function, even those not used +internally. +} + +\subsection{Algorithm}{ + +This function, unlike +\code{\link[=ci_bound_wn_i]{ci_bound_wn_i()}}, use a simple root +finding algorithm. Basically, it tries +fixing the target parameter to +different values until the likelihood +ratio test \emph{p}-value, or the +corresponding chi-square difference, +is equal to the +value corresponding to the desired +level of confidence. (Internally, +the difference between the \emph{p}-value +and the target \emph{p}-value, that for +the chi-square difference, is the +function value.) + +For finding the bound, this algorithm +can be inefficient compared to the +one proposed by Wu and Neale (2012). +The difference can be less than +one second versus 10 seconds. +It is included as a backup algorithm +for parameters which are difficult +for the method by Wu and Neale. + +Internally, it uses \code{\link[=uniroot]{uniroot()}} to +find the root. +} + +\subsection{Limitation(s)}{ + +This function does not handle an +estimate close to an attainable bound +using the method proposed by Wu and +Neale (2012). Use it for such +parameters with cautions. +} +} +\examples{ + +library(lavaan) +data(simple_med) +dat <- simple_med +mod <- +" +m ~ x +y ~ m +" +fit_med <- sem(mod, simple_med, fixed.x = FALSE) +# Remove `opts` in real cases. +# The options are added just to speed up the example +out1l <- ci_bound_ur_i(i = 1, + sem_out = fit_med, + which = "lbound", + opts = list(use_callr = FALSE, + interval = c(0.8277, 0.8278))) +out1l +} +\references{ +Wu, H., & Neale, M. C. (2012). +Adjusted confidence intervals for a +bounded parameter. \emph{Behavior +Genetics, 42}(6), 886-898. +\doi{10.1007/s10519-012-9560-z} +} +\seealso{ +\code{\link[=print.cibound]{print.cibound()}}, +\code{\link[=semlbci]{semlbci()}}, \code{\link[=ci_i_one]{ci_i_one()}}; see +\code{\link[=ci_bound_wn_i]{ci_bound_wn_i()}} on the version +for the method by Wu and Neale +(2012). +} diff --git a/man/ci_bound_wn_i.Rd b/man/ci_bound_wn_i.Rd index 3cb4f9f..e13fca9 100644 --- a/man/ci_bound_wn_i.Rd +++ b/man/ci_bound_wn_i.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/ci_bound_wn_i.R \name{ci_bound_wn_i} \alias{ci_bound_wn_i} -\title{Likelihood-based Confidence Bound For One parameter} +\title{Likelihood-based Confidence Bound By Wu-Neale-2012} \usage{ ci_bound_wn_i( i = NULL, @@ -28,6 +28,10 @@ ci_bound_wn_i( ftol_rel_factor = 1, lb_prop = 0.05, lb_se_k = 3, + try_harder = 0, + fit_lb = -Inf, + fit_ub = +Inf, + timeout = 300, ... ) } @@ -101,7 +105,11 @@ and \code{ciperc} is greater than this amount, a warning is set in the status code and the bound is set to \code{NA}. Default is 5e-4.} \item{std_method}{The method used to find the standardized -solution. If equal to \verb{"lavaan"``, [lavaan::standardizedSolution()] will be used. If equal to }"internal"\verb{, an internal function will be used. The }"lavaan"\verb{method should work in all situations, but the}"internal"\verb{method is usually much faster. Default is}"internal"`.} +solution. If equal to \code{"lavaan"}, +\code{\link[lavaan:standardizedSolution]{lavaan::standardizedSolution()}} will be used. If equal to +\code{"internal"}, an internal function will be used. The \code{"lavaan"} +method should work in all situations, but the \code{"internal"} method +is usually much faster. Default is \code{"internal"}.} \item{bounds}{Default is \code{""} and this function will set the lower bounds to \code{lb_var} for variances. Other valid values are those @@ -124,6 +132,29 @@ negative.} for free variances. Default is 3, the estimate minus 3 standard error. If negative, the lower bound is set using \code{lb_prop}.} +\item{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.} + +\item{fit_lb}{The vector of lower +bounds of parameters. Default is +\code{-Inf}, setting the lower bounds to +\code{-Inf} for all parameters except for +free variances which are controlled +by \code{lb_var}.} + +\item{fit_ub}{The vector of upper +bounds of parameters. Default is +\code{+Inf}, setting the lower bounds to +\code{+Inf} for all parameters.} + +\item{timeout}{The approximate +maximum time for the search, in +second. Default is 300 seconds +(5 minutes).} + \item{...}{Optional arguments. Not used.} } \value{ @@ -139,7 +170,8 @@ A detailed and organized output can be printed by the default print method (\code{\link[=print.cibound]{print.cibound()}}). } \description{ -Find the lower or upper bound of the likelihood-based +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 \code{\link[lavaan:lavaan]{lavaan::lavaan()}}. } diff --git a/man/ci_i_one.Rd b/man/ci_i_one.Rd index 5f044f6..bc84a91 100644 --- a/man/ci_i_one.Rd +++ b/man/ci_i_one.Rd @@ -8,7 +8,7 @@ ci_i_one( i, which = NULL, sem_out, - method = "wn", + method = c("wn", "ur"), standardized = FALSE, robust = "none", sf_full = NA, @@ -30,16 +30,19 @@ found. Must be \code{"lbound"} or \code{"ubound"}.} \link[lavaan:lavaan-class]{lavaan::lavaan} outputs only.} \item{method}{The approach to be used. Default is \code{"wn"} -(Wu-Neale-2012 Method), the only supported method.} +(Wu-Neale-2012 Method). Another method is "ur", +root finding by \code{\link[stats:uniroot]{stats::uniroot()}}.} \item{standardized}{Logical. Whether the bound of the LBCI of the standardized solution is to be searched. Default is \code{FALSE}.} \item{robust}{Whether the LBCI based on robust likelihood ratio test is to be found. Only \code{"satorra.2000"} in -\code{\link[lavaan:lavTestLRT]{lavaan::lavTestLRT()}} is supported for now. If `"none"``, the +\code{\link[lavaan:lavTestLRT]{lavaan::lavTestLRT()}} is supported for now. If \code{"none"}, the default, then likelihood ratio test based on maximum likelihood -estimation will be used.} +estimation will be used. For "ur", \code{"satorra.2000"} is +automatically used if a scaled test statistic is requested +in \code{sem_out}.} \item{sf_full}{A list with the scaling and shift factors. Ignored if \code{robust} is \code{"none"}. If \code{robust} is \code{"satorra.2000"} and diff --git a/man/gen_userp.Rd b/man/gen_userp.Rd new file mode 100644 index 0000000..3c07ce9 --- /dev/null +++ b/man/gen_userp.Rd @@ -0,0 +1,265 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/gen_userp.R +\name{gen_userp} +\alias{gen_userp} +\alias{gen_sem_out_userp} +\title{Create a Wrapper To Be Used +in 'lavaan' Models} +\usage{ +gen_userp(func, sem_out) + +gen_sem_out_userp( + userp, + sem_out, + userp_name = "semlbciuserp1234", + fix = TRUE, + control_args = list(), + iter.max = 10000, + max_attempts = 5 +) +} +\arguments{ +\item{func}{A function that receives a +\code{lavaan}-object and returns a scalar. +See Details on the restriction on this +function.} + +\item{sem_out}{A \code{lavaan}-class object +to be modified.} + +\item{userp}{A function that is +generated by \code{\link[=gen_userp]{gen_userp()}}.} + +\item{userp_name}{The name of the +function \code{userp} to be used in the +\code{lavaan} model. It does not have to +be the name of the function in +\code{userp}. Should be changed only if it +conflicts with another object in the +parent environment, which should not +happen if the model is always fitted +in a clean R session.} + +\item{fix}{If \code{TRUE}, the default, +the function generated is used to fix +the value of \code{userp} to a target +value using an equality constraint. +If \code{FALSE}, then the function simply +fits the model to the data.} + +\item{control_args}{To be passed to +the argument of the same name in +\code{\link[lavaan:lavaan]{lavaan::lavaan()}}. Default is +\code{list()}. Can be used to set the +default values of this argument +in the generated function.} + +\item{iter.max}{The maximum number of +iteration when the generated function +fit the model. Default is 10000.} + +\item{max_attempts}{If the initial +fit with the equality constraint +fails, how many more attempts +will be made by the generated +function. Default is 5.} +} +\value{ +\subsection{\code{gen_userp}}{ + +It returns a function that +accepts a numeric vector of length +equals to the number of free parameters +in \code{sem_out}, and returns a scalar +which is the output of \code{func}. If this +vector is not supplied, it will try to +find it in the \code{parent.frame()}. This +is how it works inside a \code{lavaan} +model. +} + +\subsection{\code{gen_sem_out_userp}}{ + +If \code{fix} is \code{TRUE}, it returns a +function with these arguments: +\itemize{ +\item \code{target}: The value to which the +user-defined parameter will be fixed +to. +\item \code{verbose}: If \code{TRUE}, additional +information will be printed when +fitting the model. +\item \code{control}: The values to be passed +as a list to the argument of the +same name in \code{\link[lavaan:lavaan]{lavaan::lavaan()}}. +\item \code{seed}: Numeric. If supplied, it +will be used in \code{\link[=set.seed]{set.seed()}} to +initialize the random number +generator. Necessary to reproduce +some results because random numbers +are used in some steps in \code{lavaan}. +If \code{NULL}, the default, \code{\link[=set.seed]{set.seed()}} +will not be called. +} + +If \code{fix} is `FALSE, then it returns a +function with optional arguments that +will be ignored, Calling it will +simply fit the modified model to the +data. Useful for getting the value of +the user-defined parameter. +} +} +\description{ +Make a function on +\code{lavaan} object usable in a \code{lavaan} +model syntax. +} +\details{ +\subsection{\code{gen_userp}}{ + +There are cases in which we want to +create a user parameter which is a +function of other free parameters, +computed by a function. However such +a function may work only on a +\code{lavaan} object. + +If the target function works by +extracting parameter estimates stored +in the \code{Model} slot and/or the +\code{implied} slot, then \code{\link[=gen_userp]{gen_userp()}} +can be used to convert it to a +function that retrieves the parameter +estimates when being called by +\code{\link[lavaan:lavaan]{lavaan::lavaan()}} or its wrappers, +modifies the stored +\code{lavaan} object using +\code{\link[lavaan:lav_model]{lavaan::lav_model_set_parameters()}} +and \code{\link[lavaan:lav_model]{lavaan::lav_model_implied()}} to +change the estimates, and call the +target function. + +Note that this is an unconventional +way to define a user parameter and +the generated function should always +be checked to see whether it works as +expected. + +As shown in the examples, the +parameter computed this may not +have standard error nor \emph{p}-value. + +The main purpose is for the point +estimate, for searching the +likelihood-based confidence bound +using \code{\link[=ci_bound_ur]{ci_bound_ur()}} and +\code{\link[=ci_bound_ur_i]{ci_bound_ur_i()}}. + +Note that the target function +specified in \code{func} should work +directly on the parameter estimates +stored in the \code{Model} slot and then +get the estimates using +\code{\link[lavaan:lav_model]{lavaan::lav_model_get_parameters()}}. +Functions that work on the unmodified +output generated by +\code{\link[lavaan:lavaan]{lavaan::lavaan()}} usually do not +work. + +Users are not recommended to use +\code{\link[=gen_userp]{gen_userp()}} and \code{\link[=gen_sem_out_userp]{gen_sem_out_userp()}} +directly because they require +unconventional way to extract +parameter estimates from a lavaan +model. However, developers may use +them to include functions +they wrote in a lavaan model. This +is the technique used by +\code{\link[=ci_bound_ur_i]{ci_bound_ur_i()}} to constrain any +parameter in a model to an arbitrary +value. +} + +\subsection{\code{gen_sem_out_userp}}{ + +The function \code{\link[=gen_sem_out_userp]{gen_sem_out_userp()}} +is to be used internally +for generating a function for searching +a likelihood-based confidence bound. +It is exported because it needs to +be run in an fresh external R process, +usually created by \code{callr} in other +internal functions. +} +} +\examples{ + +library(lavaan) + +data(simple_med) +dat <- simple_med +mod <- +" +m ~ a*x +y ~ b*m +ab := a*b +" +fit_med <- sem(mod, simple_med, fixed.x = FALSE) +parameterEstimates(fit_med) + +# A trivial example for verifying the results +my_ab <- function(object) { + # Need to use lav_model_get_parameters() + # because the object is only a modified + # lavaan-object, not one directly + # generated by lavaan function + est <- lavaan::lav_model_get_parameters(object@Model, type = "user") + unname(est[1] * est[2]) + } + +# Check the function +my_ab(fit_med) +coef(fit_med, type = "user")["ab"] + +# Create the function +my_userp <- gen_userp(func = my_ab, + sem_out = fit_med) +# Try it on the vector of free parameters +my_userp(coef(fit_med)) + +# Generate a modified lavaan model +fit_userp <- gen_sem_out_userp(userp = my_userp, + userp_name = "my_userp", + sem_out = fit_med) + +# This function can then be used in the model syntax. + +# Note that the following example only work when called inside the +# workspace or inside other functions such as ci_bound_ur()` +# and `ci_bound_ur_i()` because `lavaan::sem()` will +# search `my_userp()` in the global environment. + +# Therefore, the following lines are commented out. +# They should be run only in a "TRUE" interactive +# session. + +# mod2 <- +# " +# m ~ x +# y ~ m +# ab := my_userp() +# " +# fit_med2 <- sem(mod2, simple_med, fixed.x = FALSE) +# parameterEstimates(fit_med2) +# +# # Fit the model with the output of the function, a*b +# # fixed to .50 +# +# fit_new <- fit_userp(.50) +# +# # Check if the parameter ab is fixed to .50 +# parameterEstimates(fit_new) + + +} diff --git a/man/get_cibound.Rd b/man/get_cibound.Rd index 2298ed8..a9f20cd 100644 --- a/man/get_cibound.Rd +++ b/man/get_cibound.Rd @@ -2,9 +2,12 @@ % Please edit documentation in R/get_cibound.R \name{get_cibound} \alias{get_cibound} +\alias{get_cibound_status_not_0} \title{A 'cibound' Output From a 'semlbci' Object} \usage{ get_cibound(x, row_id, which = c("lbound", "ubound")) + +get_cibound_status_not_0(x) } \arguments{ \item{x}{The output of \code{\link[=semlbci]{semlbci()}}.} @@ -14,19 +17,34 @@ on the left, not the actual row number, because some rows may be omitted in the printout of \code{x}.} \item{which}{The bound for which the \code{\link[=ci_bound_wn_i]{ci_bound_wn_i()}} is -to be extracted. Either \verb{"lbound"`` or }"ubound"``.} +to be extracted. Either \code{"lbound"} or \code{"ubound"}.} } \value{ -A \code{cibound}-class object. See \code{\link[=ci_bound_wn_i]{ci_bound_wn_i()}} +\code{\link[=get_cibound]{get_cibound()}} returns a \code{cibound}-class object. See \code{\link[=ci_bound_wn_i]{ci_bound_wn_i()}} for details. +\code{\link[=get_cibound_status_not_0]{get_cibound_status_not_0()}} returns a list of +\code{cibound}-class objects with \code{status} not equal +to zero. If all bounds have \code{status} equal to +zero, it returns an empty list. } \description{ Get the \code{cibound} output of a bound from a \code{semlbci} object, the output of \code{\link[=semlbci]{semlbci()}}. } \details{ -It returns the original output of \code{\link[=ci_bound_wn_i]{ci_bound_wn_i()}} -for a bound. Usually for diagnosis. +The function \code{\link[=get_cibound]{get_cibound()}} +returns the original output of +\code{\link[=ci_bound_wn_i]{ci_bound_wn_i()}} for a bound. +Usually for diagnosis. + +The function \code{\link[=get_cibound_status_not_0]{get_cibound_status_not_0()}} +checks the status code of each bound, +and returns the \code{cibound} outputs of +bounds with status code not equal to +zero (i.e., something wrong in the +search). Printing it can print the +diagnostic information for all bounds +that failed in the search. } \examples{ diff --git a/man/semlbci-package.Rd b/man/semlbci-package.Rd index 561a3c3..0cfc260 100644 --- a/man/semlbci-package.Rd +++ b/man/semlbci-package.Rd @@ -3,7 +3,6 @@ \docType{package} \name{semlbci-package} \alias{semlbci-package} -\alias{_PACKAGE} \title{semlbci: Likelihood-Based Confidence Interval in Structural Equation Models} \description{ \if{html}{\figure{logo.png}{options: style='float: right' alt='logo' width='120'}} diff --git a/man/semlbci.Rd b/man/semlbci.Rd index 7888da6..adffe3b 100644 --- a/man/semlbci.Rd +++ b/man/semlbci.Rd @@ -12,15 +12,16 @@ semlbci( remove_intercepts = TRUE, ciperc = 0.95, standardized = FALSE, - method = "wn", + method = c("wn", "ur"), robust = c("none", "satorra.2000"), - try_k_more_times = 2, + try_k_more_times = 0, semlbci_out = NULL, check_fit = TRUE, ..., parallel = FALSE, ncpus = 2, - use_pbapply = TRUE + use_pbapply = TRUE, + loadbalancing = TRUE ) } \arguments{ @@ -56,8 +57,8 @@ interval.} \item{standardized}{If \code{TRUE}, the LBCI is for the standardized estimates.} \item{method}{The method to be used to search for the confidence -bounds. Currently only \code{"wn"} (Wu-Neale-2012), the default, is -supported.} +bounds. Supported methods are\code{"wn"} (Wu-Neale-2012), the default, +and \code{"ur"} (root finding by \code{\link[stats:uniroot]{stats::uniroot()}}).} \item{robust}{Whether the LBCI based on robust likelihood ratio test is to be found. Only \code{"satorra.2000"} in \code{\link[lavaan:lavTestLRT]{lavaan::lavTestLRT()}} @@ -90,6 +91,11 @@ cores.} is installed, \code{\link[pbapply:pbapply]{pbapply::pbapply()}} will be used to display a progress bar when finding the intervals. Default is \code{TRUE}. Ignored if \code{pbapply} is not installed.} + +\item{loadbalancing}{Whether load +balancing is used when \code{parallel} +is \code{TRUE} and \code{use_pbapply} is +\code{TRUE}.} } \value{ A \code{semlbci}-class object similar to the parameter table diff --git a/tests/testthat/test-check_sem_out.R b/tests/testthat/test-check_sem_out.R index 617480b..87f9beb 100644 --- a/tests/testthat/test-check_sem_out.R +++ b/tests/testthat/test-check_sem_out.R @@ -26,7 +26,8 @@ fit02 <- lavaan::sem(mod, dat, estimator = "MLR") # fit03 <- lavaan::sem(mod, dat, estimator = "ML", se = "robust") # (out_03 <- check_sem_out(fit03)) -fit04<- lavaan::sem(mod, dat, estimator = "DWLS") +# The warning can be ignored because this problem is intentional +suppressWarnings(fit04 <- lavaan::sem(mod, dat, estimator = "DWLS")) (out_04 <- check_sem_out(fit04)) suppressWarnings(fit05 <- lavaan::sem(mod, dat, group = "gp")) diff --git a/tests/testthat/test-ci_bound_wn_i_z_print_maxtime.R b/tests/testthat/test-ci_bound_wn_i_z_print_maxtime.R new file mode 100644 index 0000000..fa36c72 --- /dev/null +++ b/tests/testthat/test-ci_bound_wn_i_z_print_maxtime.R @@ -0,0 +1,33 @@ +skip_on_cran() + +library(testthat) +library(semlbci) + +# Fit the model + +library(lavaan) + +data(simple_med_mg) +dat <- simple_med_mg +mod <- +" +m ~ c(a1, a2)*x +y ~ c(b1, b2)*m +ab := a2*b2 +a1 == a2 +" +fit <- lavaan::sem(mod, simple_med_mg, fixed.x = FALSE, group = "gp") + +# Find the LBCIs + +ciperc <- .96 + +fn_constr0 <- set_constraint(fit, ciperc = ciperc) + +opts0 <- list() +opts0 <- list(ftol_rel = 1e-7, + maxtime = .01 + ) +time1l <- system.time(out1l <- ci_bound_wn_i(17, 16, sem_out = fit, f_constr = fn_constr0, which = "lbound", standardized = TRUE, ciperc = ciperc, opts = opts0, verbose = TRUE, wald_ci_start = FALSE, std_method = "internal")) + +expect_output(print(out1l), "timeout", fixed = TRUE) diff --git a/tests/testthat/test-ci_i_one_std_pa_ur.R b/tests/testthat/test-ci_i_one_std_pa_ur.R new file mode 100644 index 0000000..90d1371 --- /dev/null +++ b/tests/testthat/test-ci_i_one_std_pa_ur.R @@ -0,0 +1,32 @@ +skip_on_cran() +library(testthat) +library(semlbci) + +# Fit the model + +library(lavaan) + +data(simple_med) +dat <- simple_med +mod <- +" +m ~ x +y ~ m +" +fit <- lavaan::sem(mod, simple_med, fixed.x = FALSE) + +# Find the LBCIs + +ciperc <- .96 + +out1l <- ci_i_one(1, npar = 5, which = "lbound", sem_out = fit, method = "ur", + ciperc = ciperc, + standardized = TRUE, + opts = list(use_callr = FALSE)) + +# Check with known results + +test_that("Check with know results", { + expect_equal(round(unname(out1l$bounds["lbound"]), 2), .13) + }) + diff --git a/tests/testthat/test-ci_i_one_ustd_pa_ur.R b/tests/testthat/test-ci_i_one_ustd_pa_ur.R new file mode 100644 index 0000000..03d9251 --- /dev/null +++ b/tests/testthat/test-ci_i_one_ustd_pa_ur.R @@ -0,0 +1,31 @@ +skip_on_cran() +library(testthat) +library(semlbci) + +# Fit the model + +library(lavaan) + +data(simple_med) +dat <- simple_med +mod <- +" +m ~ x +y ~ m +" +fit <- lavaan::sem(mod, simple_med, fixed.x = FALSE) + +# Find the LBCIs + +ciperc <- .96 + +out1l <- ci_i_one(1, npar = 5, which = "lbound", sem_out = fit, method = "ur", + ciperc = ciperc, + opts = list(use_callr = FALSE)) + +# Check with known results + +test_that("Check with know results", { + expect_equal(round(unname(out1l$bounds["lbound"]), 2), .79) + }) + diff --git a/tests/testthat/test-semlbci_lavaan_sem_ur.R b/tests/testthat/test-semlbci_lavaan_sem_ur.R new file mode 100644 index 0000000..5490f63 --- /dev/null +++ b/tests/testthat/test-semlbci_lavaan_sem_ur.R @@ -0,0 +1,40 @@ +skip("Test parallel processing: Test in interactive sections") + +library(testthat) +library(semlbci) + +# lavaan example: sem() + +library(lavaan) +model <- ' + # latent variable definitions + ind60 =~ x1 + x2 + x3 + dem60 =~ y1 + a*y2 + b*y3 + c*y4 + dem65 =~ y5 + a*y6 + b*y7 + c*y8 + + # regressions + dem60 ~ ind60 + dem65 ~ ind60 + dem60 + + # residual correlations + y1 ~~ y5 + y2 ~~ y4 + y6 + y3 ~~ y7 + y4 ~~ y8 + y6 ~~ y8 +' + +fit <- sem(model, data = PoliticalDemocracy) +# summary(fit, fit.measures = TRUE) + +p_table <- parameterTable(fit) +i_free <- which(p_table$free > 0) + +fit_lbci_parallel <- semlbci(fit, pars = i_free[9:11], parallel = TRUE, ncpus = 3, method = "ur") +fit_lbci_no_parallel <- semlbci(fit, pars = i_free[9:11], parallel = FALSE, ncpus = 3, method = "ur") + +test_that("Compare parallel and non-parallel results", { + expect_true(all.equal(data.frame(fit_lbci_parallel[12:14, -c(18, 19)]), + data.frame(fit_lbci_no_parallel[12:14, -c(18, 19)]), + check_attributes = FALSE)) +}) diff --git a/tests/testthat/test-semlbci_lavaan_sem_ur_rb.R b/tests/testthat/test-semlbci_lavaan_sem_ur_rb.R new file mode 100644 index 0000000..80ca702 --- /dev/null +++ b/tests/testthat/test-semlbci_lavaan_sem_ur_rb.R @@ -0,0 +1,40 @@ +skip("Test parallel processing: Test in interactive sections") + +library(testthat) +library(semlbci) + +# lavaan example: sem() + +library(lavaan) +model <- ' + # latent variable definitions + ind60 =~ x1 + x2 + x3 + dem60 =~ y1 + a*y2 + b*y3 + c*y4 + dem65 =~ y5 + a*y6 + b*y7 + c*y8 + + # regressions + dem60 ~ ind60 + dem65 ~ ind60 + dem60 + + # residual correlations + y1 ~~ y5 + y2 ~~ y4 + y6 + y3 ~~ y7 + y4 ~~ y8 + y6 ~~ y8 +' + +fit <- sem(model, data = PoliticalDemocracy, test = "satorra.bentler") +# summary(fit, fit.measures = TRUE) + +p_table <- parameterTable(fit) +i_free <- which(p_table$free > 0) + +fit_lbci_parallel <- semlbci(fit, pars = i_free[9:10], parallel = TRUE, ncpus = 3, method = "ur") +fit_lbci_no_parallel <- semlbci(fit, pars = i_free[9:10], parallel = FALSE, ncpus = 3, method = "ur") + +test_that("Compare parallel and non-parallel results", { + expect_true(all.equal(data.frame(fit_lbci_parallel[12:14, -c(18, 19)]), + data.frame(fit_lbci_no_parallel[12:14, -c(18, 19)]), + check_attributes = FALSE)) +}) diff --git a/tests/testthat/test-std_lav_one_factor.R b/tests/testthat/test-std_lav_one_factor.R new file mode 100644 index 0000000..a18023a --- /dev/null +++ b/tests/testthat/test-std_lav_one_factor.R @@ -0,0 +1,51 @@ +skip_on_cran() + +library(testthat) +library(semlbci) +library(lavaan) + +dat_cov <- matrix(c(1, .30, .40, + .30, 1, .20, + .40, .20, 1), 3, 3) +colnames(dat_cov) <- rownames(dat_cov) <- c("x1", "x2", "x3") +dat_mean <- c(1, 2, 3) +names(dat_mean) <- colnames(dat_cov) + +mod <- +" +f1 =~ x1 + x2 + x3 +" +fit <- lavaan::sem(mod, + sample.cov = dat_cov, + sample.nobs = 50) + +expect_no_error(std_lav(lavaan::coef(fit), fit)) + +mod <- +" +x1 ~ x2 + x3 +" +fit <- lavaan::sem(mod, + sample.cov = dat_cov, + sample.nobs = 50) +expect_no_error(std_lav(lavaan::coef(fit), fit)) + +mod <- +" +x1 ~ x2 + x3 +" +fit <- lavaan::sem(mod, + sample.cov = dat_cov, + sample.nobs = 50) +expect_no_error(std_lav(lavaan::coef(fit), fit)) + +mod <- +" +x1 ~ x2 + x3 +" +fit <- lavaan::sem(mod, + sample.cov = dat_cov, + sample.mean = dat_mean, + sample.nobs = 50, + meanstructure = TRUE) +expect_no_error(std_lav(lavaan::coef(fit), fit)) diff --git a/tests/testthat/test-ur_add_func.R b/tests/testthat/test-ur_add_func.R new file mode 100644 index 0000000..ae81630 --- /dev/null +++ b/tests/testthat/test-ur_add_func.R @@ -0,0 +1,39 @@ +# Ready + +library(testthat) +library(lavaan) + +# cfa() example +HS.model <- ' visual =~ x1 + x2 + x3 + textual =~ x4 + x5 + x6 + speed =~ x7 + x8 + x9 + textual ~ a*visual + speed ~ b*textual + ab := a*b' +fit <- cfa(HS.model, data = HolzingerSwineford1939) + +# Get one parameter + +test_that("add_func", { + ind <- function(object) { + # Need to use lav_model_get_parameters() + # because coef() may not work. + est <- lavaan::lav_model_get_parameters(object@Model, type = "user") + unname(est[10] * est[11]) + } + fit_i_free <- add_func(func = ind, + sem_out = fit, + fix = FALSE) + fit_i <- add_func(func = ind, + sem_out = fit) + # Need to see whether fit_i() can run without ind + rm(ind) + + fit_i_tmp <- sem_out_userp_run(target = .400, + object = fit_i) + expect_equal(lavTestLRT(fit_i_tmp, fit)[2, "Df diff"], + 1) + expect_equal(unname(coef(fit_i_tmp, type = "user")["user"]), + .400, + tolerance = 1e-4) +}) diff --git a/tests/testthat/test-ur_ci_bound_ur_i_mg_std_1.R b/tests/testthat/test-ur_ci_bound_ur_i_mg_std_1.R new file mode 100644 index 0000000..334dd20 --- /dev/null +++ b/tests/testthat/test-ur_ci_bound_ur_i_mg_std_1.R @@ -0,0 +1,44 @@ +skip_on_cran() +skip("To be run in an interactive session") + +# Ready + +library(testthat) +library(lavaan) + +HS.model_gp <- ' visual =~ x1 + x2 + x3 + textual =~ x4 + x5 + x6 + visual ~~ c(a1, a2)*textual + adiff := a1^2' +fit_gp <- sem(HS.model_gp, data = HolzingerSwineford1939, + group = "school", + group.equal = "loadings") + +## Two groups + +test_that("ci_bound_ur: Two groups, standardized", { + +# Take more than one minute to run +# Test one bound is enough + +system.time( +ci_lb <- ci_bound_ur_i(i = 29, + sem_out = fit_gp, + which = "lbound", + standardized = TRUE) +) +# system.time( +# ci_ub <- ci_bound_ur_i(i = 29, +# sem_out = fit_gp, +# which = "ubound", +# standardized = TRUE) +# ) + +expect_equal(round(1 - ci_lb$diag$ciperc_final, 2), + .05, + tolerance = 1e-3) +# expect_equal(round(1 - ci_ub$diag$ciperc_final, 2), +# .05, +# tolerance = 1e-3) + +}) diff --git a/tests/testthat/test-ur_ci_bound_ur_i_mg_ustd_1.R b/tests/testthat/test-ur_ci_bound_ur_i_mg_ustd_1.R new file mode 100644 index 0000000..a0ad7e2 --- /dev/null +++ b/tests/testthat/test-ur_ci_bound_ur_i_mg_ustd_1.R @@ -0,0 +1,44 @@ +skip_on_cran() + +# Ready + +library(testthat) +library(lavaan) + +HS.model <- ' visual =~ x1 + x2 + x3 + textual =~ x4 + x5 + x6 + speed =~ x7 + x8 + x9 + visual ~~ a*textual + visual ~~ b*speed + textual ~~ d*speed + ab := a*b' +fit_gp <- sem(HS.model, data = HolzingerSwineford1939, + group = "school", + group.equal = "loadings") + +## Two groups + +test_that("ci_bound_ur: Two groups, unstandardized", { + +# Test one bound is enough + +est_i <- gen_est_i(i = 46, sem_out = fit_gp) +system.time( +ci_lb <- ci_bound_ur_i(i = 46, + sem_out = fit_gp, + which = "ubound") +) +# system.time( +# ci_ub <- ci_bound_ur_i(i = 46, +# sem_out = fit_gp, +# which = "ubound") +# ) + +expect_equal(round(1 - ci_lb$diag$ciperc_final, 2), + .05, + tolerance = 1e-3) +# expect_equal(round(1 - ci_ub$diag$ciperc_final, 2), +# .05, +# tolerance = 1e-3) + +}) diff --git a/tests/testthat/test-ur_ci_bound_ur_i_rb_ustd_1.R b/tests/testthat/test-ur_ci_bound_ur_i_rb_ustd_1.R new file mode 100644 index 0000000..87ff67b --- /dev/null +++ b/tests/testthat/test-ur_ci_bound_ur_i_rb_ustd_1.R @@ -0,0 +1,48 @@ +skip_on_cran() + +library(testthat) +library(semlbci) + +# Fit the model + +library(lavaan) + +data(simple_med) +dat <- simple_med +mod <- +" +m ~ a*x +y ~ b*m +ab := a*b +" +fit <- lavaan::sem(mod, simple_med, fixed.x = FALSE, test = "satorra.bentler") + +## One group + +test_that("ci_bound_ur: One group, unstandardized, satorra.2000", { + +# Test one bound is enough + +system.time( +ci_lb <- ci_bound_ur_i(i = 2, + sem_out = fit, + which = "lbound", + verbose = TRUE) +) +# system.time( +# ci_ub <- ci_bound_ur_i(i = 9, +# sem_out = fit, +# which = "ubound", +# verbose = TRUE) +# ) +expect_true(grepl("satorra.2000", + attr(ci_lb$diag$history$lrt, "head"), + fixed = TRUE)) +expect_equal(round(1 - ci_lb$diag$ciperc_final, 2), + .05, + tolerance = 1e-3) +# expect_equal(round(1 - ci_ub$diag$ciperc_final, 2), +# .05, +# tolerance = 1e-3) + +}) \ No newline at end of file diff --git a/tests/testthat/test-ur_ci_bound_ur_i_std_1.R b/tests/testthat/test-ur_ci_bound_ur_i_std_1.R new file mode 100644 index 0000000..40ce736 --- /dev/null +++ b/tests/testthat/test-ur_ci_bound_ur_i_std_1.R @@ -0,0 +1,43 @@ +skip_on_cran() + +# Ready + +library(testthat) +library(lavaan) + +# Standardized + +HS.model <- ' visual =~ x1 + x2 + d1*x3 + textual =~ x4 + x5 + d2*x6 + visual ~~ textual + d12 := d1^2' +fit <- sem(HS.model, data = HolzingerSwineford1939) + +## One group + +test_that("ci_bound_ur: One group, standardized", { + +# Take more than one minute to run +# Test one bound is enough + +system.time( +ci_lb <- ci_bound_ur_i(i = 7, + sem_out = fit, + which = "ubound", + standardized = TRUE) +) +# system.time( +# ci_ub <- ci_bound_ur_i(i = 7, +# sem_out = fit, +# which = "ubound", +# standardized = TRUE) +# ) + +expect_equal(round(1 - ci_lb$diag$ciperc_final, 2), + .05, + tolerance = 1e-3) +# expect_equal(round(1 - ci_ub$diag$ciperc_final, 2), +# .05, +# tolerance = 1e-3) + +}) diff --git a/tests/testthat/test-ur_ci_bound_ur_i_ustd_1.R b/tests/testthat/test-ur_ci_bound_ur_i_ustd_1.R new file mode 100644 index 0000000..1023b4e --- /dev/null +++ b/tests/testthat/test-ur_ci_bound_ur_i_ustd_1.R @@ -0,0 +1,61 @@ +skip_on_cran() + +# Ready + +library(testthat) +library(lavaan) + +HS.model <- ' visual =~ x1 + x2 + x3 + textual =~ x4 + x5 + x6 + speed =~ x7 + x8 + x9 + visual ~~ a*textual + visual ~~ b*speed + textual ~~ d*speed + ab := a*b' +fit <- sem(HS.model, data = HolzingerSwineford1939) + +# Unstandardized + +## One group + +test_that("ci_bound_ur: One group, unstandardized", { + +# Test one bound is enough + +# system.time( +# ci_ub <- ci_bound_ur_i(i = 9, +# sem_out = fit, +# which = "lbound") +# ) +system.time( +ci_ub <- ci_bound_ur_i(i = 9, + sem_out = fit, + which = "ubound") +) +# expect_equal(round(1 - ci_lb$diag$ciperc_final, 2), +# .05, +# tolerance = 1e-3) +expect_equal(round(1 - ci_ub$diag$ciperc_final, 2), + .05, + tolerance = 1e-3) + + +# system.time( +# ci_ub <- ci_bound_ur_i(i = 25, +# sem_out = fit, +# which = "lbound") +# ) +system.time( +ci_ub <- ci_bound_ur_i(i = 25, + sem_out = fit, + which = "ubound") +) +# expect_equal(round(1 - ci_lb$diag$ciperc_final, 2), +# .05, +# tolerance = 1e-3) +expect_equal(round(1 - ci_ub$diag$ciperc_final, 2), + .05, + tolerance = 1e-3) + + +}) diff --git a/tests/testthat/test-ur_ci_bound_ur_mg_std_1.R b/tests/testthat/test-ur_ci_bound_ur_mg_std_1.R new file mode 100644 index 0000000..dce239d --- /dev/null +++ b/tests/testthat/test-ur_ci_bound_ur_mg_std_1.R @@ -0,0 +1,55 @@ +skip_on_cran() +skip("To be run in an interactive session") + +# Ready + +library(testthat) +library(lavaan) + +HS.model_gp <- ' visual =~ x1 + x2 + x3 + textual =~ x4 + x5 + x6 + visual ~~ c(a1, a2)*textual + adiff := a1^2' +fit_gp <- sem(HS.model_gp, data = HolzingerSwineford1939, + group = "school", + group.equal = "loadings") + +## Two groups + +test_that("ci_bound_ur: Two groups, standardized", { + +# Take more than one minute to run +# Test one bound is enough + +est_i <- gen_est_i(i = 29, sem_out = fit_gp, standardized = TRUE) +system.time( +ci_lb <- ci_bound_ur(sem_out = fit_gp, + func = est_i, + which = "lbound", + progress = FALSE) +) +# system.time( +# ci_ub <- ci_bound_ur(sem_out = fit_gp, +# func = est_i, +# which = "ubound", +# progress = FALSE) +# ) + +expect_equal(round(ci_lb$lrt[2, "Pr(>Chisq)"], 2), + .05, + tolerance = 1e-3) +# expect_equal(round(ci_ub$lrt[2, "Pr(>Chisq)"], 2), +# .05, +# tolerance = 1e-3) + +ci_lbci <- semlbci(sem_out = fit_gp, pars = "textual =~ 2*x6", standardized = TRUE, use_pbapply = FALSE) + +expect_equal(unname(round(c(ci_lb$bound), 3)), + unname(round(unlist(confint(ci_lbci))[1], 3)), + tolerance = 1e-2) + +# expect_equal(unname(round(c(ci_lb$bound, ci_ub$bound), 3)), +# unname(round(unlist(confint(ci_lbci)), 3)), +# tolerance = 1e-2) + +}) diff --git a/tests/testthat/test-ur_ci_bound_ur_mg_ustd_1.R b/tests/testthat/test-ur_ci_bound_ur_mg_ustd_1.R new file mode 100644 index 0000000..59708f6 --- /dev/null +++ b/tests/testthat/test-ur_ci_bound_ur_mg_ustd_1.R @@ -0,0 +1,48 @@ +skip_on_cran() + +# Ready + +library(testthat) +library(lavaan) + +HS.model <- ' visual =~ x1 + x2 + x3 + textual =~ x4 + x5 + x6 + speed =~ x7 + x8 + x9 + visual ~~ a*textual + visual ~~ b*speed + textual ~~ d*speed + ab := a*b' +fit_gp <- sem(HS.model, data = HolzingerSwineford1939, + group = "school", + group.equal = "loadings") + +## Two groups + +test_that("ci_bound_ur: Two groups, unstandardized", { + +# Test one bound is enough + +est_i <- gen_est_i(i = 46, sem_out = fit_gp) +# system.time( +# ci_lb <- ci_bound_ur(sem_out = fit_gp, +# func = est_i, +# which = "lbound", +# progress = FALSE, +# use_callr = FALSE) +# ) +system.time( +ci_ub <- ci_bound_ur(sem_out = fit_gp, + func = est_i, + which = "ubound", + progress = FALSE, + use_callr = TRUE) +) + +# expect_equal(round(ci_lb$lrt[2, "Pr(>Chisq)"], 2), +# .05, +# tolerance = 1e-3) +expect_equal(round(ci_ub$lrt[2, "Pr(>Chisq)"], 2), + .05, + tolerance = 1e-3) + +}) diff --git a/tests/testthat/test-ur_ci_bound_ur_rb_mg_ustd_1.R b/tests/testthat/test-ur_ci_bound_ur_rb_mg_ustd_1.R new file mode 100644 index 0000000..30f3e52 --- /dev/null +++ b/tests/testthat/test-ur_ci_bound_ur_rb_mg_ustd_1.R @@ -0,0 +1,51 @@ +skip_on_cran() + +# Ready + +library(testthat) +library(lavaan) + +HS.model <- ' visual =~ x1 + x2 + x3 + textual =~ x4 + x5 + x6 + speed =~ x7 + x8 + x9 + visual ~~ a*textual + visual ~~ b*speed + textual ~~ d*speed + ab := a*b' +fit_gp <- sem(HS.model, data = HolzingerSwineford1939, + group = "school", + group.equal = "loadings", + test = "satorra.bentler") + +## Two groups + +test_that("ci_bound_ur: Two groups, unstandardized, satorra.2000", { + +# Test one bound is enough + +est_i <- gen_est_i(i = 46, sem_out = fit_gp) +# system.time( +# ci_lb <- ci_bound_ur(sem_out = fit_gp, +# func = est_i, +# which = "lbound", +# progress = FALSE, +# use_callr = FALSE) +# ) +system.time( +ci_ub <- ci_bound_ur(sem_out = fit_gp, + func = est_i, + which = "ubound", + progress = FALSE, + use_callr = FALSE) +) +expect_true(grepl("satorra.2000", + attr(ci_ub$lrt, "head"), + fixed = TRUE)) +# expect_equal(round(ci_lb$lrt[2, "Pr(>Chisq)"], 2), +# .05, +# tolerance = 1e-3) +expect_equal(round(ci_ub$lrt[2, "Pr(>Chisq)"], 2), + .05, + tolerance = 1e-3) + +}) diff --git a/tests/testthat/test-ur_ci_bound_ur_rb_ustd_1.R b/tests/testthat/test-ur_ci_bound_ur_rb_ustd_1.R new file mode 100644 index 0000000..56570c3 --- /dev/null +++ b/tests/testthat/test-ur_ci_bound_ur_rb_ustd_1.R @@ -0,0 +1,52 @@ +skip_on_cran() + +library(testthat) +library(semlbci) + +# Fit the model + +library(lavaan) + +data(simple_med) +dat <- simple_med +mod <- +" +m ~ a*x +y ~ b*m +ab := a*b +" +fit <- lavaan::sem(mod, simple_med, fixed.x = FALSE, test = "satorra.bentler") + +## One group + +test_that("ci_bound_ur: One group, unstandardized, satorra.2000", { + +# Test one bound is enough + +est_i <- gen_est_i(i = 2, sem_out = fit) +system.time( +ci_lb <- ci_bound_ur(sem_out = fit, + func = est_i, + which = "lbound", + progress = FALSE, + user_callr = FALSE) +) +# system.time( +# ci_ub <- ci_bound_ur(sem_out = fit, +# func = est_i, +# which = "ubound", +# progress = FALSE, +# user_callr = FALSE) +# ) + +expect_true(grepl("satorra.2000", + attr(ci_lb$lrt, "head"), + fixed = TRUE)) +expect_equal(round(ci_lb$lrt[2, "Pr(>Chisq)"], 2), + .05, + tolerance = 1e-3) +# expect_equal(round(ci_ub$lrt[2, "Pr(>Chisq)"], 2), +# .05, +# tolerance = 1e-3) + +}) \ No newline at end of file diff --git a/tests/testthat/test-ur_ci_bound_ur_std_1.R b/tests/testthat/test-ur_ci_bound_ur_std_1.R new file mode 100644 index 0000000..e90abf2 --- /dev/null +++ b/tests/testthat/test-ur_ci_bound_ur_std_1.R @@ -0,0 +1,77 @@ +skip_on_cran() + +# Ready + +library(testthat) +library(lavaan) + +# Standardized + +HS.model <- ' visual =~ x1 + x2 + d1*x3 + textual =~ x4 + x5 + d2*x6 + visual ~~ textual + d12 := d1^2' +fit <- sem(HS.model, data = HolzingerSwineford1939) + +## One group + +test_that("ci_bound_ur: One group, standardized", { + +# Take more than one minute to run +# Test one bound is enough + +est_i <- gen_est_i(i = 7, sem_out = fit, standardized = TRUE) +system.time( +ci_lb <- ci_bound_ur(sem_out = fit, + func = est_i, + which = "lbound", + progress = FALSE) +) +# system.time( +# ci_ub <- ci_bound_ur(sem_out = fit, +# func = est_i, +# which = "ubound", +# progress = FALSE) +# ) + +expect_equal(round(ci_lb$lrt[2, "Pr(>Chisq)"], 2), + .05, + tolerance = 1e-3) +# expect_equal(round(ci_ub$lrt[2, "Pr(>Chisq)"], 2), +# .05, +# tolerance = 1e-3) + +ci_lbci <- semlbci(sem_out = fit, pars = "visual ~~ textual", standardized = TRUE, use_pbapply = FALSE) + +expect_equal(unname(round(c(ci_lb$bound), 3)), + unname(round(unlist(confint(ci_lbci))[1], 3)), + tolerance = 1e-2) + +# est_i <- gen_est_i(i = 3, sem_out = fit, standardized = TRUE) +# system.time( +# ci_lb <- ci_bound_ur(sem_out = fit, +# func = est_i, +# which = "lbound", +# progress = FALSE) +# ) +# system.time( +# ci_ub <- ci_bound_ur(sem_out = fit, +# func = est_i, +# which = "ubound", +# progress = FALSE) +# ) + +# expect_equal(round(ci_lb$lrt[2, "Pr(>Chisq)"], 2), +# .05, +# tolerance = 1e-3) +# expect_equal(round(ci_ub$lrt[2, "Pr(>Chisq)"], 2), +# .05, +# tolerance = 1e-3) + +# ci_lbci <- semlbci(sem_out = fit, pars = "visual =~ x3", standardized = TRUE, use_pbapply = FALSE) + +# expect_equal(unname(round(c(ci_lb$bound), 3)), +# unname(round(unlist(confint(ci_lbci))[1], 3)), +# tolerance = 1e-2) + +}) diff --git a/tests/testthat/test-ur_ci_bound_ur_ustd_1.R b/tests/testthat/test-ur_ci_bound_ur_ustd_1.R new file mode 100644 index 0000000..1a8fef4 --- /dev/null +++ b/tests/testthat/test-ur_ci_bound_ur_ustd_1.R @@ -0,0 +1,73 @@ +skip_on_cran() + +# Ready + +library(testthat) +library(lavaan) + +HS.model <- ' visual =~ x1 + x2 + x3 + textual =~ x4 + x5 + x6 + speed =~ x7 + x8 + x9 + visual ~~ a*textual + visual ~~ b*speed + textual ~~ d*speed + ab := a*b' +fit <- sem(HS.model, data = HolzingerSwineford1939) + +# Unstandardized + +## One group + +test_that("ci_bound_ur: One group, unstandardized", { + +# Test one bound is enough + +est_i <- gen_est_i(i = 9, sem_out = fit) +# system.time( +# ci_lb <- ci_bound_ur(sem_out = fit, +# func = est_i, +# which = "lbound", +# progress = FALSE, +# user_callr = FALSE) +# ) +system.time( +ci_ub <- ci_bound_ur(sem_out = fit, + func = est_i, + which = "ubound", + progress = FALSE, + user_callr = TRUE) +) +# expect_equal(round(ci_lb$lrt[2, "Pr(>Chisq)"], 2), +# .05, +# tolerance = 1e-3) +expect_equal(round(ci_ub$lrt[2, "Pr(>Chisq)"], 2), + .05, + tolerance = 1e-3) + +est_i <- gen_est_i(i = 25, sem_out = fit) +# system.time( +# ci_lb <- ci_bound_ur(sem_out = fit, +# func = est_i, +# which = "lbound", +# progress = FALSE, +# use_callr = TRUE) +# ) +system.time( +ci_ub <- ci_bound_ur(sem_out = fit, + func = est_i, + which = "ubound", + progress = FALSE, + user_callr = FALSE) +) +# expect_equal(round(ci_lb$lrt[2, "Pr(>Chisq)"], 2), +# .05, +# tolerance = 1e-3) +expect_equal(round(ci_ub$lrt[2, "Pr(>Chisq)"], 2), + .05, + tolerance = 1e-3) + +ci_lbci <- semlbci(sem_out = fit, pars = "ab :=", use_pbapply = FALSE) + +expect_equal(unname(round(c(ci_ub$bound), 3)), + unname(round(unlist(confint(ci_lbci))[2], 3))) +}) diff --git a/tests/testthat/test-ur_gen_est_i.R b/tests/testthat/test-ur_gen_est_i.R new file mode 100644 index 0000000..c6ad342 --- /dev/null +++ b/tests/testthat/test-ur_gen_est_i.R @@ -0,0 +1,135 @@ +skip_on_cran() + +# Ready + +library(testthat) +library(lavaan) + +# cfa() example +HS.model <- ' visual =~ x1 + x2 + x3 + textual =~ x4 + x5 + x6 + speed =~ x7 + x8 + x9 + textual ~ a*visual + speed ~ b*textual + ab := a*b' +fit <- sem(HS.model, data = HolzingerSwineford1939) +fit_gp <- sem(HS.model, data = HolzingerSwineford1939, + group = "school", + group.equal = "intercepts") + +# Unstandardized + +test_that("Unstandardized", { + est_i <- gen_est_i(i = 9, sem_out = fit) + expect_equal(est_i(fit), + unname(coef(fit, type = "user")["speed=~x9"])) + est_i <- gen_est_i(i = 24, sem_out = fit) + expect_equal(est_i(fit), + unname(coef(fit, type = "user")["ab"])) + est_i <- gen_est_i(i = 46, sem_out = fit_gp) + expect_equal(est_i(fit_gp), + unname(coef(fit_gp, type = "user")["b"])) + est_i <- gen_est_i(i = 71, sem_out = fit_gp) + expect_equal(est_i(fit_gp), + unname(coef(fit_gp, type = "user")["ab"])) +}) + +# Standardized + +test_that("Standardized", { + std_all <- standardizedSolution(fit) + std_all_gp <- standardizedSolution(fit_gp) + + est_i <- gen_est_i(i = 9, sem_out = fit, standardized = TRUE) + expect_equal(est_i(fit), + unname(std_all[9, "est.std"])) + est_i <- gen_est_i(i = 24, sem_out = fit, standardized = TRUE) + expect_equal(est_i(fit), + unname(std_all[24, "est.std"])) + est_i <- gen_est_i(i = 46, sem_out = fit_gp, standardized = TRUE) + expect_equal(est_i(fit_gp), + unname(std_all_gp[46, "est.std"])) + est_i <- gen_est_i(i = 71, sem_out = fit_gp, standardized = TRUE) + expect_equal(est_i(fit_gp), + unname(std_all_gp[71, "est.std"])) +}) + +# gne_userp + +test_that("Unstandardized", { + +est_i <- gen_est_i(i = 9, sem_out = fit) +userp <- gen_userp(func = est_i, sem_out = fit) +expect_equal(userp(coef(fit)), + coef(fit, type = "user")[9], + ignore_attr = TRUE) +tmp <- coef(fit) +tmp[6] <- .30 +expect_equal(userp(tmp), + .30) + +est_i <- gen_est_i(i = 24, sem_out = fit) +userp <- gen_userp(func = est_i, sem_out = fit) +expect_equal(userp(coef(fit)), + coef(fit, type = "user")[24], + ignore_attr = TRUE) +tmp <- coef(fit) +tmp[7] <- .30 +tmp[8] <- .40 +expect_equal(userp(tmp), + .30 * .40) + +est_i <- gen_est_i(i = 46, sem_out = fit_gp) +userp <- gen_userp(func = est_i, sem_out = fit_gp) +expect_equal(userp(coef(fit_gp)), + coef(fit_gp, type = "user")[46], + ignore_attr = TRUE) +tmp <- coef(fit_gp) +tmp[37] <- .30 +expect_equal(userp(tmp), + .30) + +est_i <- gen_est_i(i = 71, sem_out = fit_gp) +userp <- gen_userp(func = est_i, sem_out = fit_gp) +expect_equal(userp(coef(fit_gp)), + coef(fit_gp, type = "user")[71], + ignore_attr = TRUE) +tmp <- coef(fit_gp) +tmp[7] <- .30 +tmp[8] <- .40 +tmp[37] <- .30 +tmp[38] <- .40 +expect_equal(userp(tmp), + .30 * .40) +}) + + +test_that("Standardized", { + + std_all <- standardizedSolution(fit) + std_all_gp <- standardizedSolution(fit_gp) + + est_i <- gen_est_i(i = 9, sem_out = fit, standardized = TRUE) + userp <- gen_userp(func = est_i, sem_out = fit) + expect_equal(userp(coef(fit)), + std_all[9, "est.std"], + ignore_attr = TRUE) + + est_i <- gen_est_i(i = 24, sem_out = fit, standardized = TRUE) + userp <- gen_userp(func = est_i, sem_out = fit) + expect_equal(userp(coef(fit)), + std_all[24, "est.std"], + ignore_attr = TRUE) + + est_i <- gen_est_i(i = 46, sem_out = fit_gp, standardized = TRUE) + userp <- gen_userp(func = est_i, sem_out = fit_gp) + expect_equal(userp(coef(fit_gp)), + std_all_gp[46, "est.std"], + ignore_attr = TRUE) + + est_i <- gen_est_i(i = 71, sem_out = fit_gp, standardized = TRUE) + userp <- gen_userp(func = est_i, sem_out = fit_gp) + expect_equal(userp(coef(fit_gp)), + std_all_gp[71, "est.std"], + ignore_attr = TRUE) +}) diff --git a/tests/testthat/test-ur_gen_userp.R b/tests/testthat/test-ur_gen_userp.R new file mode 100644 index 0000000..61a3f4e --- /dev/null +++ b/tests/testthat/test-ur_gen_userp.R @@ -0,0 +1,48 @@ +# Ready + +library(testthat) +library(lavaan) + +# cfa() example +HS.model <- ' visual =~ x1 + x2 + x3 + textual =~ x4 + x5 + x6 + speed =~ x7 + x8 + x9 + textual ~ a*visual + speed ~ b*textual + ab := a*b' +fit <- cfa(HS.model, data = HolzingerSwineford1939) + +# Get one parameter + +ind <- function(object) { + # Need to use lav_model_get_parameters() + # because coef() may not work. + est <- lavaan::lav_model_get_parameters(object@Model, type = "user") + unname(est[10] * est[11]) + } + +test_that("Indirect effect", { + userp <- gen_userp(func = ind, sem_out = fit) + expect_equal(userp(coef(fit)), + coef(fit, type = "user")[24], + ignore_attr = TRUE) + tmp <- coef(fit) + tmp[7] <- .30 + tmp[8] <- .40 + expect_equal(userp(tmp), + .30 * .40) +}) + +# gen_sem_out_userp + +skip("To be run in an interactive session") + +userp1234 <- gen_userp(func = ind, sem_out = fit) +fit_userp <- gen_sem_out_userp(userp = userp1234, + userp_name = "userp1234", + sem_out = fit) +fit_test <- fit_userp(.12) +expect_equal(coef(fit_test, type = "user")["user"], + .12, + ignore_attr = TRUE, + tolerance = 1e-4) \ No newline at end of file